clear?all;
clc;
s=tf（‘s‘）;
m=3;n=3;%m輸出，n輸入
G=cell（mn）;%系統總矩陣輸入矩陣的時候，在G里面先第一列用完，再第二列...
G{11}=4.05*exp（-28*s）/（1+50*s）;%實際對象，下面的tau值里面的分子也要對應改
G{21}=1.77*exp（-28*s）/（1+60*s）;
G{12}=5.39*exp（-21*s）/（1+50*s）;
G{22}=5.72*exp（-14*s）/（1+60*s）;
G{13}=3.02*exp（-21*s）/（1+50*s）;
G{23}=4.09*exp（-28*s）/（1+60*s）;
G{31}=3.71*exp（-21*s）/（1+50*s）;
G{32}=5.18*exp（-14*s）/（1+60*s）;
G{33}=4.78*exp（-21*s）/（1+50*s）;
Ts=7;%采樣周期
M=1;%控制時域
P=30*ones（n1）;%優化時域
%P=zeros（n1）;P（1）=20;P（2）=20;P（3）=20;%優化時域簡化處理，優化時域取一樣的值，也可以分開取不同的值
N=50;%模型長度
alpha=0.6*ones（m1）;%誤差校正系數，簡化處理取一樣的值，也可以取不同的值
%alpha=zeros（m1）;alpha（1）=0.6;alpha（2）=0.6;alpha（3）=0.6;%誤差校正參數
tau=zeros（mn）;%相對于采樣時間的滯后時間參數
tau（11）=28/Ts;
tau（21）=28/Ts;
tau（12）=21/Ts;
tau（22）=14/Ts;
tau（13）=21/Ts;
tau（23）=28/Ts;
tau（31）=21/Ts;
tau（32）=14/Ts;
tau（33）=21/Ts;
L=200;%仿真步長
aef=0*ones（m1）;%參考軌跡的柔化系數，DMC取參考軌跡為恒值c，所以取為0，也可以取別的值
%aef=zeros（m1）;aef（1）=0;aef（2）=0;aef（3）=0;%參考軌跡的柔化系數
c=zeros（m1）;%設定值
c（1）=-1;c（2）=1;c（3）=1;
plant=cell（mn）;
for?i=1:m*n
????plant{i}=c2d（G{i}Ts）;
end
stepresp=cell（mn）;
for?i=1:m*n
????stepresp{i}=step（plant{i}[Ts:Ts:Ts*N]）;
end
a=cell（mn）;
for?i=1:m*n
????a{i}=stepresp{i}（1:end）;
end
num=cell（mn）;den=cell（mn）;
for?i=1:m*n
????num{i}=plant{i}.num{1};
????den{i}=plant{i}.den{1};
end
nn=zeros（mn）;
nd=zeros（mn）;
for?i=1:m*n
????nn（i）=length（num{i}）;
????nd（i）=length（den{i}）;
end
AA=cell（mn）;
for?i=1:m*n
????AA{i}（:1）=a{i}（1:P）;
end
for?l=1:m*n
????for?i=2:P
????????for?j=2:M
????????????if?i>=j
???????????????AA{l}（ij）=AA{l}（i-1j-1）;
????????????end
????????end
????end
end
for?i=1:m
????Aa{i1}=cat（2AA{i:}）;
end
????A=cat（1Aa{1:m}）;
AA0=cell（mn）;
for?i=1:m*n
????AA0{i}=a{i}（1:N）;
end
for?i=1:m
????Aa0{i1}=cat（2AA0{i:}）;
end
????A0=cat（1Aa0{1:m}）;
QQ=cell（m1）;
for?i=1:m
????QQ{i}=diag（ones（P1）0）;
end
Q=blkdiag