dmc_modelmatched.m

一个动态矩阵预测控制程序

clear all;%以下2--5语句与[a,t]=step(tf(num,den,'inputdelay',3))等价
close all;
num=2.994;
den=[82^2 164 1];
;%基本算法对象2
Ts=5;
dt=Ts;
t0=1:dt:1045;
Delta_u=0;y=0;
[aa,t]=step(num,den,t0);
a=aa(2:209);
%也对aaa1=0;aaa2=0;
%for i=2:161
%    aaa1(i)=exp(-Ts/82)*aaa1(i-1)+2.994*(1-exp(-Ts/82))*1;%零阶保持器下的差分方程；
%     aaa2(i)=exp(-Ts/82)*aaa2(i-1)+(1-exp(-Ts/82))*aaa1(i);
%end
%a=aaa2(2:161);
P=30;
for i=1:P                                                                                                                                               
    Q(i)=1;  
end %误差加权系数

Tr=3;
ArfaR=exp(-Ts/Tr);
M=30;

for i=1:M
    R(i)=0.3;
end %控制加权系数
Q=diag(Q);%误差加权矩阵
R=diag(R);%控制加权矩阵

N=length(a);
%构造矩阵A，可以从Workspace复制，也可以由line5直接得出;
 for i=1:P
     for j=1:M  
         if i==j  
              A(i,j)=a(1);
         end 
         if i>j
             A(i,j)=a(i-j+1);
         end
         if i<j
             A(i,j)=0;
         end
     end 
 end
 %构造矩阵A0b;
 for i=1:P
     for j=1:i+1
         A0b(i,j)=a(N);
     end
     for j=i+2:N
         A0b(i,j)=a(N-j+i+1);
     end
 end
 %构造矩阵A0;
 for i=1:P
     for j=2:N-1
         A0(i,j)=A0b(i,j)-A0b(i,j+1);
     end 
     A0(i,N)=A0b(i,N);
 end
 A0(:,1)=[];%删除矩阵的第一列；
  
 h=1*ones(1,P);
 h=h';
 
 SimulationStep=100;
 D=inv(A'*Q*A+R)*A'*Q;%inv为矩阵求逆；需列计算式；
 d1=D(1,:);
 
 U(1:N-1)=0; U=U';DeltaU(1:N)=0; DeltaU=DeltaU';y=[0,0];u=0;w=1;y1=0;ym=0;y1m=0;TT=0;yN=[0,0];
  
for i=2:SimulationStep
    

     u(i)=u(i-1)+Delta_u; 
     y1(i)=exp(-Ts/82)*y1(i-1)+2.994*(1-exp(-Ts/82))*u(i);%零阶保持器下的差分方程；
     y(i)=exp(-Ts/82)*y(i-1)+(1-exp(-Ts/82))*y1(i);
     yN(i)=y(i);
     %if i>=140
     %yN(i)=y(i)+0.1;end
     
     y1m(i)=exp(-Ts/82)*y1m(i-1)+2.994*(1-exp(-Ts/82))*u(i);%零阶保持器下的差分方程；
     ym(i)=exp(-Ts/82)*ym(i-1)+(1-exp(-Ts/82))*y1m(i);
     %ym=y(i-1)+a(1)*Delta_u;


    %构造参考轨迹输出yr；
    for j=1:P
        yr(j)=ArfaR^j*y(i)+(1-ArfaR^j)*w;
    end
    
    e(i)=yN(i)-ym(i);%计算实际输出误差；
      
     Delta_u=d1*(yr'-A0*U-h*e(i));%控制增量；
     
     for j=1:N-2
         U(j)=U(j+1);
     end
     U(N-1)=u(i);
        
     for j=1:N-1
         DeltaU(j)=DeltaU(j+1);
     end
     DeltaU(N)=Delta_u;
     TT(i)=i;
     
 end     

  TT=TT*Ts;
figure;  
% subplot(1,2,1);
plot(TT,yN)
grid
% subplot(1,2,2);
% plot(TT,u)
% grid
% ylabel('output');