程序2.m

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

clc
clear all;
close all;
lanbuda=0.1;
ysp=20;
sysc1=tf([36.96],[40 14 1],'inputdelay',1.5);
sysd1=c2d(sysc1,1.5);
[num1,den1]=tfdata(sysd1,'v');
a0=-den1(2);
a1=den1(3);
b0=num1(2);
b1=num1(3);
lw=b0/(b0^2+lanbuda);
ly1=-lw*(1+a0);
ly2=lw*(a0+a1);
ly3=-lw*a1;
lu=-lw*b1;
for k=1:1:20
    switch k
        case 1,
            yout(k)=0;
            yout_esti(k+1)=0;
            deltu(k)=lw*ysp;
            u(k)=deltu(k);
        case 2,
            yout(k)=a0*yout(k-1);
            yout_esti(k+1)=a0*yout_esti(k)+b0*u(k-1);            deltu(k)=lw*ysp+ly1*yout_esti(k+1)+ly2*yout_esti(k)+lu*deltu(k-1);
            u(k)=u(k-1)+deltu(k);
        case 3,
yout(k)=a0*yout(k-1)-a1*yout(k-2)+b0*u(k-2);            yout_esti(k+1)=a0*yout_esti(k)-a1*yout_esti(k-1)+b0*u(k-1)+b1*u(k-2);            deltu(k)=lw*ysp+ly1*yout_esti(k+1)+ly2*yout_esti(k)+ly3*yout_esti(k-1)+lu*deltu(k-1);
u(k)=u(k-1)+deltu(k); 
    otherwise
yout(k)=a0*yout(k-1)-a1*yout(k-2)+b0*u(k-2)+b1*u(k-3);            yout_esti(k+1)=a0*yout_esti(k)-a1*yout_esti(k-1)+b0*u(k-1)+b1*u(k-2);            deltu(k)=lw*ysp+ly1*yout_esti(k+1)+ly2*yout_esti(k)+ly3*yout_esti(k-1)+lu*deltu(k-1);
u(k)=u(k-1)+deltu(k); 
    end
end
subplot(2,1,1);
plot(yout);
title('Y')
axis([0,20,0,30]);
subplot(2,1,2);
plot(u*60);
title('U')
axis([0,20,-2000,2000]);