%DMC控制算法
clc;
clear?all;
G=input（‘輸入傳遞函數G=‘）;%輸入傳遞函數
%判斷是否為穩定系統，若是可以控制，若不是，則無法用DMC算法進行控制
den=G.den{1};%取傳函的分母
p=real（roots（den））%求傳函的極點的實部
for?i=1:length（p）
???????r=p（i）;
???if?r>0???%若有某一個極點的實部的實部大于零，則為不穩定系統，DMC無法控制
???????pG???%在命令窗口顯示極點和傳函
Error=（‘您要控制的對象為不穩定系統，DMC算法只適用于穩定系統!‘）
????return
??end
end
%設置DMC參數
Ts=input（‘采樣周期?Ts=?‘）;%采樣時間
P=input（‘預測時域?P=?‘）;%預測步長
M=input（‘控制時域?M=?‘）;%控制步長
N=80;%截斷步長
%設定參考值
yr=10;?%系統期望輸出
%建立系統階躍響應模型
[y0t0]=step（G0:5:500）;
%初始化DMC
A=zeros（PM）;%動態矩陣
a=zeros（N1）;
for?i=1:N
????a（i）=y0（i）;
end
for?i=1:P
????for?j=1:M
????????if?i-j+1>0
????????????A（ij）=a（i-j+1）;?%構造矩陣A
????????end
????end
end
%初始化向量ys，yue和矩陣A0
ys=ones（N1）;
y=zeros（N1）;
u=zeros（N1）;
e=zeros（N1）;
A0=zeros（PN-1）;
for?i=1:P
????for?j=N-2:-