clear
disp（‘多變量系統的隱式廣義預測控制算法的研究（pt70）‘）
disp（‘廣義預測控制算法初始值‘）
nn=input（‘時域長度nn=‘）;
n=input（‘預測長度n=‘）;
m=input（‘控制長度m=‘）;
t0=input（‘控制加權系數t0=‘）;
a=input（‘柔化系數a=‘）;
disp（‘最小二乘公式初始值‘）
t1=1;%input（‘遺忘因子t1=‘）;
t2=zeros（2001）;
t3=zeros（3001）;
d1=input（‘（2*n+1）階方陣P的形式:0-對角陣1-方陣:‘）;
d2=input（‘（2*n+1）階方陣P的初始值:0-自動賦值1e+51-鍵盤輸入:‘）;
if（d1==1）
????if（d2==1）
????????P=input（‘在方括號[]中，輸入（2*n+1）階方陣P的值:‘）;
????else
????????P=（1e+5）*ones（2*n+1）;
????end
else
????if（d2==1）
????????PP=input（‘在方括號[]中，輸入（2*n+1）階對角方陣P對角線上的值:‘）;
????????P=diag（PP）;
????else
????????P=（1e+5）*eye（2*n+1）;
????end
end
%參數初始值
uuu1=0;yyy1=0;uuu2=0;yyy2=0;
uu1=zeros（n1）;u1=zeros（m1）;uu2=zeros（n1）;u2=zeros（m1）;
yy1=zeros（n1）;y11=zeros（n1）;yy2=zeros（n1）;y12=zeros（n1）;
Q1=zeros（2*n+11）;Q1（11）=1;Q1（n+11）=1;Q1（2*n+11）=1;
Q2=Q1;
%產生時間為T，幅值為1的階躍信號的給定值
T=300;
t=1:T;
yr0（t1）=heaviside（t）;
d3=input（‘輸出曲線是否去掉前100步:0-不1-去掉:‘）;
nm=length（t）;%確定循環次數
for?ij=2:nm
????yr1=yr0（ij）;
????yr2=yr1;
????%根據系統模型，計算K時刻的輸出值y1（k）和y2（k）
????y1=1.629*yy1（n1）-0.6787*yy1（n-11）+0.0202*yy1（n-21）-0.0012*uu1（n1）+0.2167*uu1（n-11）-0.4678*uu1（n-21）+0.2543*uu1（n-31）+0.0003*uu2（n-11）+0.0396*uu2（n-21）+0.0301*uu2（n-31）-0.0025*uu2（n-41）;
????y2=2.208*yy2（n1）-1.5457*yy2（n-11）+0.3339*yy2（n-21）-0.0008*uu1（n-11）-0.0061*uu1（n-21）+0.0091*uu1（n-31）-0.0027*uu1（n-41）+0.0052*uu2（n-11）-0.0007*uu2（n-21）-0.0033*uu2（n-31）;
????%（nn=n=6;m=2;t=0.8;a=0.3;t1=1）
????%保存k時刻及以前的n個輸出值y（k）y（k-1）…y（k-n）以供模型運算
????for?i=1:n-1
????????yy1（i1）=yy1（i+11）;yy2（i1）=yy2（i+11）;
????end
????yy1（n1）=y1;yy2（n1）=y2;
????yyy1=[yyy1;（y1-1.2970）*210];yyy2=[yyy2;（y2+0.7014）*18];%保存各k時刻的nm個輸出量以便繪圖
????%根據最小二乘公式，由y1（k）和y2（k）計算G11，G12，G21和G22陣的各元素值g0g1…gn
????for?i=1:n
????????X（1i）=uu1（i1）;X（1i+n）=uu2（i1）;
????end
????X（12*n+1）=1;
????K=P*X‘*inv（t1+X*P*X‘）;
????P=（eye（2*n+1）-K*X）*P/t1;
????Q1=Q1+K*（y1-X*Q1）;Q2=Q2+K*（y2-X*Q2）;
????%根據元素值g0g1…gn求G11G12G21和G22陣
????for?j=1:m
????????for?i=n:-1:j
????????????i1=n-i+j;
????????????G11（i1j）=Q1（i1）;G12（i1j）=Q1（i+n1）;
????????????G21（i1j）=Q2（i1）;G22（i1j）=Q2（i+n1）;
????????end
????end
????%for?i=1:nG12（i2）=0;G21（i2）=0;end
????%求nn維y01y02向量（y11和y12為上一時刻的y01和y02向量）
????e1=y1-y11（11）;e2