预测函数控制预测函数控制 MATLABMATLAB 程序程序程序一：跟踪阶跃函数clear all; close all;Ts=0.5;%采样周期 0.5s Tr=1;H=10;%参考轨迹时间常数为 1,预测时域优化长度为 10 Km=3;Tm=3;%预测模型的参数 beta=exp(-Ts/Tr); alph=exp(-Ts/Tm);%对象离散化 sys=tf(3,3,1); dsys=c2d(sys,Ts,z); num,den=tfdata(dsys,v);%预测模型离散化 sysm=tf(Km,Tm,1); dsysm=c2d(sysm,Ts,z); numm,denm=tfdata(dsysm,v);u_1=0; y_1=0; ym_1=0; c=1;%设定值 for k=1:1:100time(k)=k*Ts;%经 z 变换后的离散化对象y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;%修正后的过程输出值 ypav(k)=y(k)ypav(k)=y(k)%控制器输出u(k)=(c-ypav(k)*(1-betaH)+ym(k)*(1-alphH)/(Km*(1-alphH);%参数更新u_1=u(k); y_1=y(k);ym_1=ym(k);end figure(1); plot(time,c,b,time,y,r); axis(0 50 0 1.2) xlabel(time(s);ylabel(c,y);程序二：跟踪阶跃加扰动后响应clear all; close all;Ts=0.5;%采样周期 0.5s Tr=1;H=10;%参考轨迹时间常数为 1,预测时域优化长度为 10 Km=3;Tm=3;%预测模型的参数 beta=exp(-Ts/Tr); alph=exp(-Ts/Tm);%对象离散化 sys=tf(3,3,1); dsys=c2d(sys,Ts,z); num,den=tfdata(dsys,v);%预测模型离散化 sysm=tf(Km,Tm,1); dsysm=c2d(sysm,Ts,z); numm,denm=tfdata(dsysm,v);u_1=0; y_1=0; ym_1=0; c=1;%设定值 for k=1:1:40time(k)=k*Ts;%经 z 变换后的离散化对象y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;%控制器输出u(k)=(c-y(k)*(1-betaH)+ym(k)*(1-alphH)/(Km*(1-alphH);%参数更新u_1=u(k); y_1=y(k);ym_1=ym(k);end for k=41:1:44time(k)=k*Ts;c=1.2;%经 z 变换后的离散化对象y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;%控制器输出u(k)=(c-y(k)*(1-betaH)+ym(k)*(1-alphH)/(Km*(1-alphH);%参数更新u_1=u(k); y_1=y(k);ym_1=ym(k);end for k=45:1:140time(k)=k*Ts;c=1;%经 z 变换后的离散化对象y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;%控制器输出u(k)=(c-y(k)*(1-betaH)+ym(k)*(1-alphH)/(Km*(1-alphH);%参数更新u_1=u(k); y_1=y(k);ym_1=ym(k);end figure(1); plot(time,y,r); axis(0 70 0 1.2); xlabel(time(s);ylabel(c,y);程序三：跟踪斜坡函数clear all; close all;Ts=0.5;%采样周期 0.5s Tr=1;H1=4;H2=5;%参考轨迹时间常数为 1,预测步长 H1 为 4，H2 为 5 Km=15;Tm=15;%预测模型的参数 beta=exp(-Ts/Tr); alph=exp(-Ts/Tm);%对象离散化 sys=tf(5,5,1); dsys=c2d(sys,Ts,z); num,den=tfdata(dsys,v);%预测模型离散化 sysm=tf(Km,Tm,1); dsysm=c2d(sysm,Ts,z); numm,denm=tfdata(dsysm,v);u_1=0; y_1=0; ym_1=0; s1=0; s2=0;G1(H1)=Km*(1-alphH1); G1(H2)=Km*(1-alphH2); for i=1:1:H1-1s1=s1+i*alph(H1-1-i); end; for j=1:1:H2-1s2=s2+i*alph(H2-1-i); end; G2(H1)=Km*(1-alph)*s1;G2(H2)=Km*(1-alph)*s2;for k=1:1:200time(k)=k*Ts;c(k)=time(k);%设定值%c(k+H1)=(k+H1)*Ts;%c(k+H2)=(k+H2)*Ts;%经 z 变换后的离散化对象y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;%修正后的过程输出值 ypav(k)=y(k)%ypav(k)=y(k);%yr(k+H1)=e(k+H1)-betaH1*(c(k)-ypav(k)其中的 e(k+H1)应该为 c(k+H1)%即 yr(k+H1)=c(k+H1)-betaH1*(c(k)-ypav(k)%X1=yr(k+H1)-alphH1*ym(k)-e(k+H1)yr(k+H1)=(k+H1)*Ts-betaH1*(k+0)*Ts-y(k);yr(k+H2)=(k+H2)*Ts-betaH2*(k+0)*Ts-y(k);e(k+H1)=y(k)-ym(k);e(k+H2)=y(k)-ym(k);X1=yr(k+H1)-alphH1*ym(k)-e(k+H1);X2=yr(k+H2)-alphH2*ym(k)-e(k+H2);%控制器输出u(k)=(G2(H2)*X1-G2(H1)*X2)/(G1(H1)*G2(H2)-G1(H2)*G2(H1);%参数更新u_1=u(k); y_1=y(k);ym_1=ym(k);end figure(1); plot(time,c,b-,time,y,r-); %axis(0 50 0 1.2) xlabel(time(s);ylabel(c,y);程序四：跟踪设定曲线clear all; close all;Ts=0.5;%采样周期 0.5s Tr=1;H1=4;H2=5;%参考轨迹时间常数为 1,预测步长 H1 为 4，H2 为 5 Km=15;Tm=15;%预测模型的参数 beta=exp(-Ts/Tr); alph=exp(-Ts/Tm);%对象离散化 sys=tf(5,5,1); dsys=c2d(sys,Ts,z); num,den=tfdata(dsys,v);%预测模型离散化 sysm=tf(Km,Tm,1); dsysm=c2d(sysm,Ts,z); numm,denm=tfdata(dsysm,v);u_1=0; y_1=0; ym_1=0; s1=0; s2=0;G1(H1)=Km*(1-alphH1); G1(H2)=Km*(1-alphH2); for i=1:1:H1-1s1=s1+i*alph(H1-1-i); end; for j=1:1:H2-1s2=s2+i*alph(H2-1-i); end; G2(H1)=Km*(1-alph)*s1; G2(H2)=Km*(1-alph)*s2;for k=1:1:80time(k)=k*Ts;c(k)=time(k);%设定值%c(k+H1)=(k+H1)*Ts;%c(k+H2)=(k+H2)*Ts;%经 z 变换后的离散化对象y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;%修正后的过程输出值 ypav(k)=y(k)%ypav(k)=y(k);%yr(k+H1)=e(k+H1)-betaH1*(c(k)-ypav(k)其中的 e(k+H1)应该为 c(k+H1)%即 yr(k+H1)=c(k+H1)-betaH1*(c(k)-ypav(k)%X1=yr(k+H1)-alphH1*ym(k)-e(k+H1)yr(k+H1)=(k+H1)*Ts-betaH1*(k+0)*Ts-y(k);yr(k+H2)=(k+H2)*Ts-betaH2*(k+0)*Ts-y(k);e(k+H1)=y(k)-ym(k);e(k+H2)=y(k)-ym(k);X1=yr(k+H1)-alphH1*ym(k)-e(k+H1);X2=yr(k+H2)-alphH2*ym(k)-e(k+H2);%控制器输出u(k)=(G2(H2)*X1-G2(H1)*X2)/(G1(H1)*G2(H2)-G1(H2)*G2(H1);%参数更新u_1=u(k); y_1=y(k);ym_1=ym(k);end for k=81:1:160time(k)=k*Ts;c(k)=40;%设定值y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;yr(k+H1)=40-betaH1*(40-y(k);yr(k+H2)=40-betaH2*(40-y(k);e(k+H1)=y(k)-ym(k);e(k+H2)=y(k)-ym(k);X1=yr(k+H1)-alphH1*ym(k)-e(k+H1);X2=yr(k+H2)-alphH2*ym(k)-e(k+H2);u(k)=(G2(H2)*X1-G2(H1)*X2)/(G1(H1)*G2(H2)-G1(H2)*G2(H1);u_1=u(k); y_1=y(k);ym_1=ym(k);end for k=161:1:200time(k)=k*Ts;c(k)=-1.5*k*Ts+160;y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;yr(k+H1)=-1.5*(k+H1)*Ts+160-betaH1*(c(k)-y(k);yr(k+H2)=-1.5*(k+H2)*Ts+160-betaH2*(c(k)-y(k);e(k+H1)=y(k)-ym(k);e(k+H2)=y(k)-ym(k);X1=yr(k+H1)-alphH1*ym(k)-e(k+H1);X2=yr(k+H2)-alphH2*ym(k)-e(k+H2);u(k)=(G2(H2)*X1-G2(H1)*X2)/(G1(H1)*G2(H2)-G1(H2)*G2(H1);u_1=u(k); y_1=y(k);ym_1=ym(k);end for k=200:1:300time(k)=k*Ts;c(k)=10;y(k)=-den(2)*y_1+num(2)*u_1;ym(k)=-denm(2)*ym_1+numm(2)*u_1;yr(k+H1)=10-betaH1*(c(k)-y(k);yr(k+