先进控制大作业DMC多变量

班级：控制工程17硕学号：6171913006姓名：岳慧慧DMC控制算法求解两输入两输出例题1、MATLAB.m文件程序%DMC控制算法clc;clearall;G11=input('输入传递函数G11=');%输入传递函数%判断是否为稳定系统，若是可以控制，若不是，则无法用DMC算法进行控制den=G11.den{1};%取传函的分母p=real(roots(den));%求传函的极点的实部fori=1:length(p)r=p(i);ifr0%若有某一个极点的实部的实部大于零，则为不稳定系统，DMC无法控制p,G11%在命令窗口显示极点和传函Error=('您要控制的对象为不稳定系统，DMC算法只适用于稳定系统!')returnendendG12=input('输入传递函数G12=');%输入传递函数%判断是否为稳定系统，若是可以控制，若不是，则无法用DMC算法进行控制den=G12.den{1};%取传函的分母p=real(roots(den));%求传函的极点的实部fori=1:length(p)r=p(i);ifr0%若有某一个极点的实部的实部大于零，则为不稳定系统，DMC无法控制p,G12%在命令窗口显示极点和传函Error=('您要控制的对象为不稳定系统，DMC算法只适用于稳定系统!')returnendendG21=input('输入传递函数G21=');%输入传递函数%判断是否为稳定系统，若是可以控制，若不是，则无法用DMC算法进行控制den=G21.den{1};%取传函的分母p=real(roots(den));%求传函的极点的实部fori=1:length(p)r=p(i);ifr0%若有某一个极点的实部的实部大于零，则为不稳定系统，DMC无法控制p,G21%在命令窗口显示极点和传函Error=('您要控制的对象为不稳定系统，DMC算法只适用于稳定系统!')returnendendG22=input('输入传递函数G22=');%输入传递函数%判断是否为稳定系统，若是可以控制，若不是，则无法用DMC算法进行控制den=G22.den{1};%取传函的分母p=real(roots(den));%求传函的极点的实部fori=1:length(p)r=p(i);ifr0%若有某一个极点的实部的实部大于零，则为不稳定系统，DMC无法控制p,G22%在命令窗口显示极点和传函Error=('您要控制的对象为不稳定系统，DMC算法只适用于稳定系统!')returnendend%设置DMC参数Ts=input('采样周期Ts=');%采样时间P=input('预测时域P=');%预测步长M=input('控制时域M=');%控制步长N=50;%截断步长%设定参考值yr=10;%建立系统阶跃响应模型G1=G11+G12;[y01,t01]=step(G1,0:5:500);G2=G21;[y02,t02]=step(G2,0:5:500);%初始化DMCA1=zeros(P,M);%动态矩阵A2=zeros(P,M);%动态矩阵a1=zeros(N,1);a2=zeros(N,1);fori=1:Na1(i)=y01(i);a2(i)=y02(i);endfori=1:Pforj=1:Mifi-j+10A1(i,j)=a1(i-j+1);%构造矩阵AA2(i,j)=a2(i-j+1);%构造矩阵Aendendend%初始化向量ys1，ys2,y11,u11,y12,u12,e1,e2和矩阵A01,A02ys1=ones(N,1);ys2=ones(N,1);y11=zeros(N,1);u11=zeros(N,1);y12=zeros(N,1);u12=zeros(N,1);e1=zeros(N,1);e2=zeros(N,1);A01=zeros(P,N-1);A02=zeros(P,N-1);fori=1:Pforj=N-2:-1:1ifN-j+1+i-1=NA01(i,j)=a1(N-j+1+i-1)-a1(N-j+i-1);%构造矩阵A01A02(i,j)=a2(N-j+1+i-1)-a2(N-j+i-1);%构造矩阵A02elseA01(i,j)=0;A02(i,j)=0;endendA01(i,N-1)=a1(i+1);A02(i,N-1)=a2(i+1);end%DMC程序fork=2:NUk_1=zeros(N-1,1);Uk_2=zeros(N-1,1);fori=1:N-1ifk-N+i=0Uk_1(i)=0;Uk_2(i)=0;elseUk_1(i)=u1(k-N+i);Uk_2(i)=u2(k-N+i);endendY01=A01*Uk_1;Y02=A02*Uk_2;e1(k)=y11(k-1)-Y01(1);e2(k)=y12(k-1)-Y02(1);Yr=zeros(P,1);fori=1:PYr(i)=yr;endEk1=zeros(P,1);Ek2=zeros(P,1);fori=1:PEk1(i)=e1(k);Ek2(i)=e2(k);enddelta_u1=inv(A1'*A1+eye(M))*A1'*(Yr-Y01-Ek1);%控制增量的计算delta_u2=inv(A2'*A2+eye(M))*A2'*(Yr-Y02-Ek2);%控制增量的计算fori=1:Mifk+i-1=Nu11(k+i-1)=u11(k+i-1-1)+delta_u1(i);%控制律的计算u12(k+i-1)=u12(k+i-1-1)+delta_u2(i);%控制律的计算endendtemp1=0;%设置在k-j-1时刻以前的控制律temp2=0;%设置在k-j-1时刻以前的控制律forj=1:N-1ifk-j=0temp1;temp2;elseifk-j-1=0temp1=temp1+a1(j)*u11(k-j);temp2=temp2+a2(j)*u12(k-j);elsetemp1=temp1+a1(j)*(u11(k-j)-u11(k-j-1));temp2=temp2+a2(j)*(u12(k-j)-u12(k-j-1));endendendifk-N=0y11(k)=temp1+e1(N);y12(k)=temp2+e2(N);elsey11(k)=temp1+a1(N)*u1(k-N)+e1(N);y12(k)=temp2+a2(N)*u2(k-N)+e2(N);endend%画图显示结果t=10*(1:N);subplot(211);plot(t,y11);plot(t,y12);title('DMC控制输出曲线');xlabel('t')ylabel('y')gridonsubplot(212);plot(t,u11,'r');plot(t,u12,'r');title('控制作用');xlabel('t')ylabel('u')gridon2、调试输入命令：G11=tf(0.5,[21],'inputdelay',2);G12=tf(0.8,[1.51],'inputdelay',3.2);G21=tf(1.8,[3.61],'inputdelay',1);G22=tf(0,[1]);Ts=30P=40M=5闭环系统的输出y11、y12与控制量u11、u12变化曲线如下：