clear?all;
close?all;
clc;

ts=1;
sys1=tf（1[4.7871?1]）;????????%比例系數？？
sys1=sys1*exp（-10）;
dsys1=c2d（sys1ts‘z‘）;
[num1den1]=tfdata（dsys1‘v‘）;?%上水箱

u_2=0;
u1_1=0;
y2_1=0;
y1_1=0.1;
ei=0;
x=[000]‘;
error_1=0;
kp=80;ki=2;kd=0;?
%?p=13;???????%??p=17m=3.%外環
%?m=2;
ywt=[];
uwt=1;
N=500;

%?yr=50;

for?s=1:1:100?

????time（s）=s*ts;?
????%Linear?model
????G=tf（[0.1426][536?46.3024?1]）;
????G=G*y1_1;
????P=50;M=20;
????A=zeros（PM）;%動態矩陣
????a=zeros（N1）;

????[y0t0]=step（G0:1:500）;

????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:-1:1
????????????if?N-j+1+i-1<=N
????????????????A0（ij）=a（N-j+1+i-1）-a（N-j+i-1）;%構造矩陣A0
????????????else
????????????????A0（ij）=0;
????????????end
????????end
????????A0（iN-1）=a（i+1）;
????end
????%DMC程序
????for?k=2:N
????????if?k>=N/2&&k<=N/2+3
????????????yr=60;
????????else
????????????yr=50;
????????end
????????Uk_1=zeros（N-11）;
????????for?i=1:N-1
????????????if?k-N+i<=0
????????????????Uk_1（i）=0;
????????????else
????????????????Uk_1（i）=u（k-N+i）;
????????????end
????????end
????????Y0=A0*Uk_1;
????????e（k）=y（k-1）-Y0（1）;
????????Yr=zeros（P1）;
????????for?i=1:P
????????????Yr（i）=yr;
????????end
????????Ek=zeros（P1）;
????????for?i=1:P
????????????Ek（i）=e（k）;
????????end
????????delta_u=inv（A‘*A+eye（M））*A‘*（Yr-Y0-Ek）;?%控制增量的計算
????????for?i=1:M
????????????if?k+i-1<=N
????????????????u（k+i-1）=u（k+i-1-1）+delta_u（i）;??%控制律的計算
????????????end
????????end
????????temp=0;%設置在k-j-1時刻以前的控制律
????????for?j=1:N-1
????????????if?k-j<=0
????????????????temp;
????????????else
????????????????if?k-j-1<=0
????????????????????temp=temp+a（j）*u（k-j）;
????????????????else
????????????????????temp=temp+a（j）*（u（k-j）-u（k-j-1））;
????????????????end
????????????end
????????end
????????if?k-N<=0
????????????y（k）=temp+e（N）;
????????else
????????????y（k）=temp+a（N）*u（k-N）+e（N）;
????????end
????end
?
????error（s）=u（s）-y1_1;
????x（1）=error（s）;?????????????????%Calculating?P
????x（2）=（error（s）-error_1）/ts;??%Calculating?D
????x（3）=x（3）+error（s）*ts;????????%Calculating?I

????u1（s）=kp*x（1）+kd*x（2）+ki*x（3）;?
????y1（s）=-den1（2）*y1_1+num1（2）*u1（s）;??%內環
????%?
????error_1=error（s）;
????%?u1_1=u1（k）;
????y1_1=y1（s）;
????%?u_2=u（k）;
????y2_1=y（s）;

end
?t=（1:N）;
%?subplot（211）;
plot（ty‘r‘）
hold?on;
%?title（‘PID串級PID和預測PID輸出‘）;
xlabel（‘t‘）
ylabel（‘y‘）
%grid?on
%?subplot（212）;
%?plot（tu‘r‘）;
%?title（‘控制作用‘）;
%?xlabel（‘t‘）
%?ylabel（‘u‘）
%grid?on
ts=1;
sys1=tf（0.1426[4.78711]）;

dsys1=c2d（sys1ts‘z‘）;
[num1den1]=tfdata（dsys1‘v‘）;

sys2=tf（1[23.04941]）;
%?sys2=sys2*exp（-10）;
dsys2=c2d（sys2ts‘z‘）;
[num2den2]=tfdata（dsy