程序1.m

我做的关于广义预测控制的小程序

clear all
close all
ysp=1;
lanbuda=0.00055;
sysc1=tf([2.45],[30.5 1],'inputdelay',35.8);
sysd1=c2d(sysc1,1);
[num1,den1]=tfdata(sysd1,'v');
a_1=-den1(2);
b0_1=num1(2);
b1_1=num1(3);
sysc2=tf([2.45],[35.2 1],'inputdelay',35.8);
sysd2=c2d(sysc2,1);
[num2,den2]=tfdata(sysd2,'v');
a_2=-den2(2);
b0_2=num2(2);
b1_2=num2(3);
sysc3=tf([2.45],[35.2 1],'inputdelay',35.1);
sysd3=c2d(sysc3,1);
[num3,den3]=tfdata(sysd3,'v');
a_3=-den3(2);
b0_3=num3(2);
b1_3=num3(3);
theta=[0.96774 0.016013 0.063012]';
a_esti=theta(1);
b0_esti=theta(2);
b1_esti=theta(3);
a=0.82;
lw=b0_esti/(b0_esti^2+lanbuda);
lw1=lw*(1-a^2);
ly1=-lw*(1+a_esti-a^2);
ly2=lw*a_esti;
lu=-lw*b1_esti;
P=(10^6)*eye(3);
for k=1:1:200
     switch k
      case 1,
         yout(k)=0;
         deltu(k)=lw1*ysp;
         u(k)=0+deltu(k);
       case 2,
          yout(k)=0;
          yout_esti(k+35)=b0_1*u(1);
          deltu(k)=lw1*ysp+ly1*yout_esti(k+35)+lu*deltu(k-1);
          u(k)=u(k-1)+deltu(k);
      case {3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36}
       yout(k)=0;
       yout_esti(k+35)=a_1*yout_esti(k+34)+b0_1*u(k-1)+b1_1*u(k-2);    
       deltu(k)=lw1*ysp+ly1*yout_esti(k+35)+ly2*yout_esti(k+34)+lu*deltu(k-1);
       u(k)=u(k-1)+deltu(k);
      case 37,
       yout(k)=b0_1*u(k-36);
       yout_esti(k+35)=a_1*yout_esti(k+34)+b0_1*u(k-1)+b1_1*u(k-2);       deltu(k)=lw1*ysp+ly1*yout_esti(k+35)+ly2*yout_esti(k+34)+lu*deltu(k-1);
       u(k)=u(k-1)+deltu(k);
       X=[yout(37) u(2) u(1)]';
     otherwise 
         if(38<=k)&(k<=50)
             yout(k)=a_1*yout(k-1)+b0_1*u(k-36)+b1_1*u(k-37);
         end
          if(51<=k)&(k<=100)
              yout(k)=a_2*yout(k-1)+b0_2*u(k-36)+b1_2*u(k-37);
          end
           if(101<=k)&(k<=150)
            yout(k)=a_3*yout(k-1)+b0_3*u(k-36)+b1_3*u(k-37);
           end
           if(151<=k)&(k<=200)
            yout(k)=a_1*yout(k-1)+b0_1*u(k-36)+b1_1*u(k-37);
           end
           K=P*X/(0.7+X'*P*X);
P=(P-K*X'*P)/0.7;
theta=theta+K*(yout(k)-X'*theta);
X=[yout(k) u(k-35) u(k-36)]';
a_esti=theta(1);
b0_esti=theta(2);
b1_esti=theta(3);
lw=b0_esti/(b0_esti^2+lanbuda);
lw1=lw*(1-a^2);
ly1=-lw*(1+a_esti-a^2);
ly2=lw*a_esti;
lu=-lw*b1_esti;
yout_esti(k+35)=a_esti*yout_esti(k+34)+b0_esti*u(k-1)+b1_esti*u(k-2);
deltu(k)=lw1*ysp+ly1*yout_esti(k+35)+ly2*yout_esti(k+34)+lu*deltu(k-1);
u(k)=u(k-1)+deltu(k);
     end
end
figure(1)
plot(u)
title('u')
xlabel('t/s')
ylabel('u')
axis([0,200,0,2])
grid
figure(2)
plot(yout)
title('Y')
xlabel('t/s')
ylabel('Y')
axis([0,200,0,2])
grid