基于预测控制模型mpc的matlab仿真.zip

clear disp('单变量系统的隐式广义预测控制算法的研究') disp('广义预测控制算法初始值') nn=input('时域长度nn='); n=input('预测长度n='); m=input('控制长度m='); t0=input('控制加权系数t0='); a=input('柔化系数a='); disp('最小二乘公式初始值') t1=1; %input('遗忘因子=')； d1=input('n+1阶方阵P的形式：0-对角阵，1-方阵：'); d2=input('n+1阶方阵P的初始值：0-自动赋值1e+5，1-键盘输入：'); if(d1==1) if(d2==1) P=input('在方括号【】中，输入（n+1）阶方阵P的值：'); else p=(1e+5)*ones(n+1); end else if(d2==1) PP=input('在方括号【】中，输入（n+1）阶对角阵P对角线上的值：'); P=diag(PP); else P=(1e+5)*eye(n+1); end end %参数初始值 uuu=0;yyy=0; uu=zeros(n,1);u=zeros(m,1); yy=zeros(n,1);y1=zeros(n,1); Q=zeros(n+1,1);Q(1,1)=1;Q(n+1,1)=1; %产生周期为100，时间为T，幅值为1的方波信号的给定值 T=300;[yr0,t]=gensig('square',100,T,1); d3=input('输出曲线是否去掉前100步：0-不，1-去掉：'); nm=length(t); for ij=2:nm yr=yr0(ij)+1; %产生周期为100，时间为T，幅值在1和2之间变化的方波信号的给定值 %根据系统模型，计算k时刻的输出值y（k） y=1.496585*yy(n,1)-0.496585*yy(n-1,1)+0.5*uu(n-1,1); %(nn=6;n=6;m=2;t0=0.8;a=0.3;t1=1) %example9-1仿真模型 %y=2.001676*yy(n,1)-1.24339*yy(n-1,1)+0.24171*yy(n-2,1)+0.23589*uu(n,1); %example9-2仿真模型 %产生均匀分布的白噪声 a9=0; for i=1:1 a9=a9+rand; end a8=0.01*(a9-6); %保存k时刻及以前的n个输出值y(k),y(k-1),...,y(k-n),以供模型运算 for i=1:n-1 yy(i,1)=yy(i+1,1); end yy(n,1)=y; yyy=[yyy;y]; %保存各k时刻nm个输出以便绘图 %根据最小二乘公式，由y(k)计算G阵的各元素值g0,g1,...gn for i=1:n X(1,i)=uu(i,1); end X(1,n+1)=1; K=P*X'*inv(t1+X*P*X'); P=(eye(n+1)-K*X)*P/t1; Q=Q+K*(y-X*Q); %根据元素值g0,g1,...,gn,求G阵 for j=1:m for i=n:-1:j i1=n-i+j; G(i1,j)=Q(i,1); end end %求nn维y0向量（y1为上一时刻的y0向量） y0(1:nn-1,1)=y1(2:nn,1);y0(nn,1)=y1(nn,1); y0=y0+(y-y1(1,1)); for i=1:n y1(i,1)=y0(i,1)+u(1,1)*Q(n+1-i,1); end for i=n+1:nn y1(i,1)=y1(n,1); end %根据y0,求n维f向量f（k+1）,...f(k+n) f(1:n,1)=y0(1:n,1); %由当前k时刻的输出值y(k)和给定值yr，求k时刻以后的n个参考轨迹w(k+1),...w(k+n) w=a*y+(1-a)*yr; for i=2:n w=[w;a^i*y+(1-a^i)*yr]; end %计算k时刻及以后的m个控制增量Du（k），,,,Du(k+m) u=inv(G'*G+t0*eye(m))*G'*(w-f); %保存k时刻及以后的n个控制增量，以供模型运算 for i=1:n-1 uu(i,1)=uu(i+1,1); end uu(n,1)=u(1,1); uuu=[uuu;u(1,1)]; %保存各k时刻的nm个控制增量以便绘图 %控制量限幅 if(u(1,1)>1) u(1,1)=1; end if(u(1,1)<-1) u(1,1)=-1; end end %绘制给定值、输出值和控制增量曲线 if(d3==1) %绘制去掉前100步的给定值、输出值和控制增量曲线 yyy1(1:(T-100),1)=yyy(101:T,1); uuu1(1:(T-100),1)=uuu(101:T,1); t1(1:(T-100),1)=t(101:T,1)-100; yr01(1:(T-100),1)=yr0(101:T,1); subplot(2,1,1):plot(t1,(yr01+1),t1,yyy1); axis([0,nm-100,0,2.5]); xlabel('t');ylabel('yr,y') subplot(2,1,2):plot(t1,uuu1); axis([0,nm-100,-1.5,1.5]); xlabel('t');ylabel('u') else %绘制完整的给定值，输出值和控制增量曲线 subplot(2,1,1);plot(t,(yr0+1),t,yyy); axis([0,nm,0,2.5]); xlabel('t');ylabel('yr,y') subplot(2,1,2);plot(t,uuu); axis([0,nm,-1.5,1.5]); xlabel('t');ylabel('u') end