%產生周期為200時間為T，幅值為l的方波信號的給定值
T=600;
[yr0t]=gensig（‘square‘200T1）;
%plot（tyr0）

nn=input（‘時域長度?nn=‘）;
n=input（‘預測長度?n=‘）;
m=input（‘控制長度?m=‘）;
lmd=input（‘控制加權系數?lmd=‘）;
a=input（‘柔化系數?a=‘）;
P=（1e+5）*eye（n+1）;
%參數初始值
uuu=0;yyy=0;
uu=zeros（n1）;u=zeros（m1）;
yy=zeros（n1）;y1=zeros（n1）;
Q=zeros（n+11）;
Q（11）=1;
Q（n+11）=1;
t1=1;
nm=length（t）;
for?ij=2:nm
????yr=yr0（ij）;
%根據系統模型，計算k時刻的輸出值y（k）
%?????A=yy（n-2）*yy（n-1）;
%?????B=1.2*yy（n-2）^2+0.8*yy（n-1）^2+1;
%?????y=A./B+0.5*uu（n1）;
y=0.1*（2.571*yy（n-11）-2.806*yy（n-21）+2.04*yy（n-31）-1.095*yy（n-41）+...
???0.2818*yy（n-51）+0.001283+0.0243*uu（n1）+0.0473*uu（n-11）+...
????0.01576*uu（n-21）+0.005498*uu（n-31））;
????%產生均勻分布的白噪聲
????%保存k時刻及以前的n個輸出值y（k）y（k-1）...y（k-n）
????for?i=1:n-1
????????yy（i1）=yy（i+11）;
????end
????yy（n1）=y;
????yyy=[yyy;y];?%保存各k時刻的n個輸出量以便繪圖
????%根據最小二乘公式．由y（k）計算G陣的各元素值g0g0...gn
????for?i=1:n
????????X（1i）=uu（i1）;
????end
????%?????a9=0;
%?????for?i=1:1
%?????????a9=a9+rand;
%?????end
%?????a8=0.01*（a9-6）;

????X（1n+1）=1;
????K=P*X‘*inv（t1+X*P*X‘）;
????P=（eye（n+1