dmc7arfa.m

此算法能实现预测控制理论中不同衰减因子的影响及鲁棒性验证

clear all;%以下2--5语句与[a,t]=step(tf(num,den,'inputdelay',3))等价
num=2.994;
den=[82^2 164 1];
      %简化算法1对象2
s=tf('s');
sysA=2.994/(82^2*s^2+82*2*s+1);
 
Ts=5;
Delta_u=0;y=0;

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=150;
for i=1:P
    Q(i)= 1;  
end %误差加权系数

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

for i=1:M
    R(i)=0.1;
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';
 
%简化算法；
 arfa=0.9;
 for i=1:M
     Aarfa(i)=arfa^(i-1);
 end
 Aarfa=Aarfa';
 Au=A*Aarfa;
 
 
 SimulationStep=200;
 D=inv(Au'*Q*Au+Aarfa'*R*Aarfa)*Au'*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;y2=0;y3=0;TT=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);
     
     
     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)=y(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;
%简化算法；
 arfa=0.5;
  for i=1:M
     Aarfa1(i)=arfa^(i-1);
 end
 Aarfa1=Aarfa1';
 Au1=A*Aarfa1;
 
 SimulationStep=200;
 D1=inv(Au1'*Q*Au1+Aarfa1'*R*Aarfa1)*Au1'*Q;%inv为矩阵求逆；需列计算式；
 d11=D1(1,:);


 
 U1(1:N-1)=0; U1=U1';DeltaU1(1:N)=0; DeltaU1=DeltaU1';z=[0,0];u1=0;w=1;z1=0;zm=0;z1m=0;z2=0;z3=0;TT=0;Delta_u1=0;

 %仿真；
for i=2:SimulationStep
    

     u1(i)=u1(i-1)+Delta_u1;  
     z1(i)=exp(-Ts/82)*z1(i-1)+2.994*(1-exp(-Ts/82))*u1(i);%零阶保持器下的差分方程，对象；
     z(i)=exp(-Ts/82)*z(i-1)+(1-exp(-Ts/82))*z1(i);
     
     
     z1m(i)=exp(-Ts/82)*z1m(i-1)+2.994*(1-exp(-Ts/82))*u1(i);%零阶保持器下的差分方程，模型；
     zm(i)=exp(-Ts/82)*zm(i-1)+(1-exp(-Ts/82))*z1m(i);
     %ym=y(i-1)+a(1)*Delta_u;
     
     
     %构造参考轨迹输出yr；
    for j=1:P
        zr(j)=ArfaR^j*z(i)+(1-ArfaR^j)*w;
    end
    
    e(i)=z(i)-zm(i);%计算实际输出误差；
      
     Delta_u1=d11*(zr'-A0*U1-h*e(i));%控制增量；

     
     for j=1:N-2
         U1(j)=U1(j+1);
     end
     U1(N-1)=u1(i);
        
     for j=1:N-1
         DeltaU1(j)=DeltaU1(j+1);
     end
     DeltaU1(N)=Delta_u1;
     
     TT(i)=i;
 end 
 TT=TT*Ts;
figure;  
 hold on
 
plot(TT,y,'k-')
plot(TT,z,'k-')
grid