预测控制的DMC代码

单输入单输出的预测控制算法clear%离线准备工作Ts=1;%采样周期tfinal=70;%采样结束时间N=tfinal/Ts;%建模时域t=[0:1:70];P=6;%优化时域M=2;%控制时域sys=tf(1,[50,15,1],'inputdelay',1);%测得系统模型sysr=tf(1,[60,18,1.2],'inputdelay',1);%实际系统模型[yp,tp,xp]=step(sys,t);[ypr,tpr,xpr]=step(sys,t);model=[yp(2:end)];%因为step的第一个值是从0时刻开始的，而模型阶跃系数是从s1第一个时刻开始的modelr=[ypr(2:end)];A=zeros(P,M);fori1=1:Pfori2=1:Mif(i1i2)A(i1,i2)=0;elseA(i1,i2)=model(i1-i2+1);endendendQ=eye(P);%误差权矩阵R=0*eye(M);%控制权矩阵cT=[1,zeros(1,M-1)];%M维行向量，取首元素K=inv(A'*Q*A+R)*A'*Q;dT=cT*K;%控制向量,P维行向量h=ones(N,1);%选择校正系数%下面完成在线计算%首先初始化模块ym=zeros(N,1);%预测初值即y0tend=50;%运行时间nn=tend/Ts+1;sp=2;%设定值w=ones(P,1)*sp;S=zeros(N);%求移位矩阵Sfori=1:Nforj=1:Nif(i+1==j)S(i,j)=1;endendendS(N,N)=1;ss=[1,zeros(1,N-1)];%取输出的第一个元素du=0;uu=0;rd=-0.2;%扰动yn=zeros(N,1);%预测输出yn0=zeros(N,1);%初始预测值yr=zeros(N,1);%系统未来N个时刻实际的输出值y=zeros(nn,1);%实际值%在线计算%检测实际值fork=1:nnyr=yr+modelr*du;y(k)=ss*yr+rd;%系统实际的输出yr=S*yr;%移位e=y(k)-ss*yn;ycor=yn+h*e;%预测值校正yn0=S*ycor;%初始预测值du=dT*(w-yn0(1:P));uu=uu+du;if(uu30)uu=30;endif(uu-30)uu=-30;endyn=yn0+model*du;%预测输出u(k)=uu;endt=0:Ts:tend;subplot(211);plot(t,y);gridon;subplot(212);stairs(t,u);gridon;