MPC.rar_MPC matlab_control_matlab mpc_model mpc_predictive contr

程序一：跟踪阶跃函数 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:100 time(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-beta^H)+ym(k)*(1-alph^H))/(Km*(1-alph^H)); %参数更新 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:40 time(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-beta^H)+ym(k)*(1-alph^H))/(Km*(1-alph^H)); %参数更新 u_1=u(k); y_1=y(k); ym_1=ym(k); end for k=41:1:44 time(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-beta^H)+ym(k)*(1-alph^H))/(Km*(1-alph^H)); %参数更新 u_1=u(k); y_1=y(k); ym_1=ym(k); end for k=45:1:140 time(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-beta^H)+ym(k)*(1-alph^H))/(Km*(1-alph^H)); %参数更新 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-alph^H1); G1(H2)=Km*(1-alph^H2); for i=1:1:H1-1 s1=s1+i*alph^(H1-1-i); end; for j=1:1:H2-1 s2=s2+i*alph^(H2-1-i); end; G2(H1)=Km*(1-alph)*s1; G2(H2)=Km*(1-alph)*s2; for k=1:1:200 time(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)-beta^H1*(c(k)-ypav(k))其中的e(k+H1)应该为c(k+H1) %即yr(k+H1)=c(k+H1)-beta^H1*(c(k)-ypav(k)) %X1=yr(k+H1)-alph^H1*ym(k)-e(k+H1) yr(k+H1)=(k+H1)*Ts-beta^H1*((k+0)*Ts-y(k)); yr(k+H2)=(k+H2)*Ts-beta^H2*((k+0)*Ts-y(k)); e(k+H1)=y(k)-ym(k); e(k+H2)=y(k)-ym(k); X1=yr(k+H1)-alph^H1*ym(k)-e(k+H1); X2=yr(k+H2)-alph^H2*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-alph^H1); G1(H2)=Km*(1-alph^H2); for i=1:1:H1-1 s1=s1+i*alph^(H1-1-i); end; for j=1:1:H2-1 s2=s2+i*alph^(H2-1-i); end; G2(H1)=Km*(1-alph)*s1; G2(H2)=Km*(1-alph)*s2; for k=1:1:80 time(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)-beta^H1*(c(k)-ypav(k))其中的e(k+H1)应该为c(k+H1) %即yr(k+H1)=c(k+H1)-beta^H1*(c(k)-ypav(k)) %X1=yr(k+H1)-alph^H1*ym(k)-e(k+H1) yr(k+H1)=(k+H1)*Ts-beta^H1*((k+0)*Ts-y(k)); yr(k+H2)=(k+H2)*Ts-beta^H2*((k+0)*Ts-y(k)); e(k+H1)=y(k)-ym(k); e(k+H2)=y(k)-ym(k); X1=yr(k+H1)-alph^H1*ym(k)-e(k+H1); X2=yr(k+H2)-alph^H2*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:160 time(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-beta^H1*(40-y(k)); yr(k+H2)=40-beta^H2*(40-y(k)); e(k+H1)=y(k)-ym(k); e(k+H2)=y(k)-ym(k); X1=yr(k+H1)-alph^H1*ym(k)-e(k+H1); X2=yr(k+H2)-alph^H2*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:200 time(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-beta^H1*(c(k)-y(k)); yr(k+H2)=-1.5*(k+H2)*Ts+160-beta^H2*(c(k)-y(k)); e(k+H1)=y(k)-ym(k); e(k+H2)=y(k)-ym(k); X1=yr(k+H1)-alph^H1*ym(k)-e(k+H1); X2=yr(k+H2)-alph^H2*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:300 time(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-beta^H1*(c(k)-y(k)); yr(k+H2)=10-beta^H2*(c(k)-y(k)); e(k+H1)=y(k)-ym(k); e(k+H2)=y(k)-ym(k); X1=yr(k+H1)-alph^H1*ym(k)-e(k+H1); X2=yr(k+H2)-alph^H2*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 150 0 45]) xlabel('time(s)');ylabel('c,y');