预测函数控制一阶二阶matlab代码

clear all; close all; clc P1=4; P2=8; P3=15; P4=20; Tr=10; Ts=1; beta=exp(-Ts/Tr); numm=[1.5]; denm=[7 1]; num=[1.5]; den=[7 1]; [A,B,C,D]=tf2ss(num,den); [Gd,Hd,Cd,Dd]=c2dm(A,B,C,D,Ts); [Am,Bm,Cm,Dm]=tf2ss(numm,denm); [Gmd,Hmd,Cmd,Dmd]=c2dm(Am,Bm,Cm,Dm,Ts); u_1=0; X_1=0; Xm_1=0; for j=1:P1 gk1(j)=Cmd*Gmd^(j-1)*Hmd; end gk1=sum(gk1); for k=1:P2 gk2(k)=Cmd*Gmd^(k-1)*Hmd; end gk2=sum(gk2); for m=1:P3 gk3(m)=Cmd*Gmd^(m-1)*Hmd; end gk3=sum(gk3); for n=1:P4 gk4(n)=Cmd*Gmd^(n-1)*Hmd; end gk4=sum(gk4); gk=[gk1 gk2 gk3 gk4]; M=inv(gk*gk')*gk; k0=M*[1-beta^P1 1-beta^P2 1-beta^P3 1-beta^P4]'; km=(-M)*[Cmd*(Gmd^P1-1) Cmd*(Gmd^P2-1) Cmd*(Gmd^P3-1) Cmd*(Gmd^P4-1) ]'; for k=1:1:1200 time(k)=k*Ts; c(k)=1.0+0.005*time(k)-0.0005*time(k)^2+0.0000005*time(k)^3; X(k)=Gd*X_1+Hd*u_1; y(k)=Cd*X(k); Xm(k)=Gmd*Xm_1+Hmd*u_1; ym(k)=Cmd*Xm(k); u(k)=k0*[c(k)-y(k)]+km*Xm(k); u_1=u(k); X_1=X(k); Xm_1=Xm(k); end plot(time,u,'b--',time,y,'r-.',time,c,'k'); xlabel('t/s'); ylabel('c,y,u'); legend('控制量','过程输出','设定值')