预测控制DMC算法matlab

%DMC控制算法clc;clearall;G=input('输入传递函数G=');%输入传递函数%判断是否为稳定系统，若是可以控制，若不是，则无法用DMC算法进行控制den=G.den{1};%取传函的分母p=real(roots(den))%求传函的极点的实部fori=1:length(p)r=p(i);ifr0%若有某一个极点的实部的实部大于零，则为不稳定系统，DMC无法控制p,G%在命令窗口显示极点和传函Error=('您要控制的对象为不稳定系统，DMC算法只适用于稳定系统!')returnendend%设置DMC参数Ts=input('采样周期Ts=');%采样时间P=input('预测时域P=');%预测步长M=input('控制时域M=');%控制步长N=80;%截断步长%设定参考值yr=10;%系统期望输出%建立系统阶跃响应模型[y0,t0]=step(G,0:5:500);%初始化DMCA=zeros(P,M);%动态矩阵a=zeros(N,1);fori=1:Na(i)=y0(i);endfori=1:Pforj=1:Mifi-j+10A(i,j)=a(i-j+1);%构造矩阵Aendendend%初始化向量ys，y,u,e和矩阵A0ys=ones(N,1);y=zeros(N,1);u=zeros(N,1);e=zeros(N,1);A0=zeros(P,N-1);fori=1:Pforj=N-2:-1:1ifN-j+1+i-1=NA0(i,j)=a(N-j+1+i-1)-a(N-j+i-1);%构造矩阵A0elseA0(i,j)=0;endendA0(i,N-1)=a(i+1);end%DMC程序fork=2:NUk_1=zeros(N-1,1);fori=1:N-1ifk-N+i=0Uk_1(i)=0;elseUk_1(i)=u(k-N+i);endendY0=A0*Uk_1;e(k)=y(k-1)-Y0(1);Yr=zeros(P,1);fori=1:PYr(i)=yr;endEk=zeros(P,1);fori=1:PEk(i)=e(k);enddelta_u=inv(A'*A+eye(M))*A'*(Yr-Y0-Ek);%控制增量的计算fori=1:Mifk+i-1=Nu(k+i-1)=u(k+i-1-1)+delta_u(i);%控制律的计算endendtemp=0;%设置在k-j-1时刻以前的控制律forj=1:N-1ifk-j=0temp;elseifk-j-1=0temp=temp+a(j)*u(k-j);elsetemp=temp+a(j)*(u(k-j)-u(k-j-1));endendendifk-N=0y(k)=temp+e(N);elsey(k)=temp+a(N)*u(k-N)+e(N);endend%画图显示结果t=10*(1:N);subplot(211);plot(t,y);title('DMC控制输出曲线');xlabel('t')ylabel('y')gridonsubplot(212);plot(t,u,'r');title('控制作用');xlabel('t')ylabel('u')gridon