matlab_MIMO_实现集中式MPC和分散式MPC的控制算法

clc; clear; Ts=4;%采样周期 N=100;%建模时域 t=0:Ts:N*Ts;%设置一个时间矢量 M=3;%控制时域 P=15;%优化时域 %% 求对象的阶跃响应模型 sys11=tf(4.05,[50,1],'inputdelay',27)%系统模型 sysd11=c2d(sys11,Ts);%将连续传递函数转换为离散传递函数模型（采样时间为Ts） [yp11,tp11,xp11]=step(sysd11,t);%计算阶跃响应，yp为输出 a11=yp11(2:N+1); sys12=tf(1.77,[60,1],'inputdelay',28) sysd12=c2d(sys12,Ts); [yp12,tp12,xp12]=step(sysd12,t); a12=yp12(2:N+1); sys13=tf(5.88,[50,1],'inputdelay',27) sysd13=c2d(sys13,Ts); [yp13,tp13,xp13]=step(sysd13,t); a13=yp13(2:N+1); sys21=tf(5.39,[50,1],'inputdelay',18) sysd21=c2d(sys21,Ts); [yp21,tp21,xp21]=step(sysd21,t); a21=yp21(2:N+1); sys22=tf(5.72,[60,1],'inputdelay',14) sysd22=c2d(sys22,Ts); [yp22,tp22,xp22]=step(sysd22,t) a22=yp22(2:N+1); sys23=tf(6.9,[40,1],'inputdelay',15) sysd23=c2d(sys23,Ts); [yp23,tp23,xp23]=step(sysd23,t); a23=yp23(2:N+1); sys31=tf(4.38,[33,1],'inputdelay',20) sysd31=c2d(sys31,Ts); [yp31,tp31,xp31]=step(sysd31,t); a31=yp31(2:N+1); sys32=tf(4.42,[44,1],'inputdelay',22) sysd32=c2d(sys32,Ts); [yp32,tp32,xp32]=step(sysd32,t); a32=yp32(2:N+1); sys33=tf(7.2,[19,1]) sysd33=c2d(sys33,Ts); [yp33,tp33,xp33]=step(sysd33,t); a33=yp33(2:N+1); %% 参数设置 A=decentralization_matrixA(a11,a22,a33,P,M);%构造动态矩阵A Q = 1; Q1=Q*eye(P);%误差权矩阵 Q2=Q*eye(P); Q3=Q*eye(P); Q=mdiag(Q1,Q2,Q3); R = 1; R1=R*eye(M);%控制权矩阵 R2=R*eye(M); R3=R*eye(M); R=mdiag(R1,R2,R3); H1=ones(N,1);%校正向量 H2=ones(N,1); H3=ones(N,1); H=mdiag(H1,H2,H3); S=zeros(N);% 求移位矩阵S for i=1:N for j=1:N if(j==i+1) S(i,j)=1; end end end S(N,N)=1; ss=eye(N); %% DMC程序初始化,设置初值 yp01=zeros(N,1);yp02=zeros(N,1);yp03=zeros(N,1);%设置预测初值yn0 ypm1=zeros(N,1);ypm2=zeros(N,1);ypm3=zeros(N,1);%设置预测值yn1 y1=zeros(N,1);y2=zeros(N,1);y3=zeros(N,1);%实际预测输出值 u=zeros(N,3);%控制量 delat_u=zeros(N,3);%控制增量 uu=zeros(1,3);%控制量初值 du=zeros(1,3);%控制增量初值 %% 在线计算优化 for i=1:N y1(i)=ss(i,1:N)*ypm1;y2(i)=ss(i,1:N)*ypm2;y3(i)=ss(i,1:N)*ypm3;%当前时刻真实值 yp01=S*ypm1;yp02=S*ypm2;yp03=S*ypm3;%移位设置该时刻的预测初值 [du(1),du(2),du(3)]=decentralization_getdu(A,P,M,Q,R,[yp01(1:P);yp02(1:P);yp03(1:P)],u(i,:),Ts);%计算控制增量 du=[du(1),du(2),du(3)]; uu=uu+du;%计算控制量 u(i,:)=uu; delat_u(i,:)=du; ypm1=yp01+a11*du(1);%该时刻预测值向量 ypm2=yp02+a22*du(2); ypm3=yp03+a33*du(3); end %% 画图输出 t=0:N-1; figure; subplot(311);plot(t,y1);title('y1');xlabel('Time(s)');ylabel('y');grid on; subplot(312);plot(t,y2);title('y2');xlabel('Time(s)');ylabel('y');grid on; subplot(313);plot(t,y3);title('y3');xlabel('Time(s)');ylabel('y');grid on; figure; subplot(311);stairs(t,u(:,1));title('u1');xlabel('Time(s)');ylabel('u');grid on; subplot(312);stairs(t,u(:,2));title('u2');xlabel('Time(s)');ylabel('u');grid on; subplot(313);stairs(t,u(:,3));title('u3');xlabel('Time(s)');ylabel('u');grid on; figure; subplot(311);stairs(t,delat_u(:,1));title('delat_u1');xlabel('Time(s)');ylabel('du');grid on; subplot(312);stairs(t,delat_u(:,2));title('delat_u2');xlabel('Time(s)');ylabel('du');grid on; subplot(313);stairs(t,delat_u(:,3));title('delat_u3');xlabel('Time(s)');ylabel('du');grid on;