对于双输入双输出系统的模型预测控制DMC的MATLAB实现

对于双输入双输出系统的模型预测控制（DMC）的MATLAB实现，能直接运行得到结论，是学习DMC的一个较好的资源。

clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%基础参数设置
delt=1;%采样周期
N=150; %模型长度
p=2;m=2; %系统维数-2输入2输出
t=0:delt:delt*N;
P=10;%优化时域 
M=5; %控制时域
nn=N/delt+1; %在线计算步数
%输出误差加权系数-单位阵 
ramena_y=1;q=ramena_y*eye（P）;
%控制加权系数-对角线元素为0.02的对角阵
ramena_u=0.02;r=ramena_u*eye（M）;
rd=0;%扰动

%传递函数建立
sys11=tf（0.35[500 45 1]‘inputdelay‘1）;
sys12=tf（0.50[150 25 1]‘inputdelay‘2.5）;
sys21=tf（0.85[2000 90 1]‘inputdelay‘2）;
sys22=tf（-0.90[1200 70 1]‘inputdelay‘1）;

%建立阶跃响应模型
[y11t11x11]=step（sys11t）;model{11}=y11（1:N）;
[y12t12x12]=step（sys12t）;model{12}=y12（1:N）;
[y21t21x21]=step（sys21t）;model{21}=y21（1:N）;
[y22t22x22]=step（sys22t）;model{22}=y22（1:N）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%模型预测
%建立动态响应矩阵 A
for i=1:p
  for j=1:m  %对 p*m 个优化变量赋值
    for i1=1:M
      for i2=1:P
        if i2           B{ij}（i2i1）=0; %超出控制时域时，控制增量无变化
        else
           B{ij}（i2i1）=model{ij}（i2-i1+1）;%未超出时赋值为对应的阶跃响应系数
        end
      end
    end
  end
  C{i}=cat（2B{i1:m}）; %联结为动态响应矩阵的行数
end
A=cat（1C{1:p}）;
%按列连接为动态响应矩阵 由阶跃响应系数组成的 P*M阵
%{}中下标 i 表示所预测未来输出的个数，下标 j 表示控制量变化的次数
for i=1:p
  aa{i}=cat（2model{i1:m}）;
end
A0=cat（1aa{1:p}）;
ym=zeros（p*N1）;
y=zeros（nnp）;
du=zeros（m1）;
uu=zeros（m1）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%滚动优化
%求控制向量 dT（由 Q R P M 四个参数确定）
for i=1:p 
    qq{i}=q;
end
Q=blkdiag（qq{1:p}）;
for i=1:m
    rr{i}=r;
end
R=blkdiag（rr{1:m}）;
L=zeros（mm*M）;  

for i=1:m
      L（iM*（i-1）+1）=1;
 end
D=L*inv（A‘*Q*A+R）*A‘*Q; %控制向量计算公式

sp1=1;sp2=1;%预测输出设定值 可调整分别不一样
w{1}=sp1*ones（P1）;
w{2}=sp2*ones（P1）;
W=cat（1w{1:p}）;
%输出期望值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%反馈校正
%建立矩阵HhstSS0
hst=ones（N1）;
for i=1:p
    hh{i}=hst;
end
  H=blkdiag（hh{1:p}）;
  ss=zeros（pp*N）;
  
for i=1:m
    ss（iN*（i-1）+1）=1;
end
  S=zeros（N）;
  
for i=1:N
    for j=1:N
       if j==i+1
          S（ij）=1;
       end
    end
end
S（NN）=1;

for i=1:p
    s{i}=S;
end
S0=blkdiag（s{1:p}）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%在线计算初始状态设置
u1_1=0;u1_2=0;
u2_1=0;u2_2=0;
y1_1=0;y1_2=0;y1_3=0;y1_4=0;
y2_1=0;y2_2=0;
y3_1=0;y3_2=0;
y4_1=0;y4_2=0;
u1=0;u2=0;
z11=c2d（sys110.5*delt‘zoh‘）;
[num1den1]=tfdata（z11‘v‘）;
z12=c2d（sys120.5*delt‘zoh‘）;
[num2den2]=tfdata（z12‘v‘）;
z21=c2d（sys210.5*delt‘zoh‘）;
[num3den3]=tfdata（z21‘v‘）;
z22=c2d（sys220.5*delt‘zoh‘）;
[num4den4]=tfdata（z22‘v‘）;

%在线计算
for k=1:nn
y111=num1（1）*u1+num1（2）*u1_

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件    2134489  2019-06-21 16:14  DMC_ttf_huitudou.pdf

     文件       4599  2019-06-21 20:34  DMC.m

----------- ---------  ---------- -----  ----

              2139088                    2