【文章內(nèi)容簡介】 c1=c0+h1*d1。 %前面三句是使用遞推公式計算待估計的參數(shù)：c(K+1)=c(N)+K(N+1)*P(N)*x(N+1)*[y(N+1)x39。(N+1)*c(N)] e1=c1c0。 %求估計參數(shù)的誤差 e2=e1./c0。 e(:,k)=e2。 %完善誤差矩陣e c0=c1。 %為下次循環(huán)做準(zhǔn)備 c(:,k)=c1。 %每次循環(huán)都使待估計參數(shù)矩陣增加一列，記錄每一次的估計值 p1=p0h1*x139。*p0。 %利用遞推公式：P(N+1)=P(N)K(N+1)*P(N)*x(N+1)*x39。(N+1)P(N)來計算P(N+1) p0=p1。 %為下次循環(huán)做準(zhǔn)備 if e2=E %判斷估計參數(shù)與實際參數(shù)之間的誤差是否滿足要求 m=k %若滿足要求，則跳出循環(huán)，并且輸出停止估計的次數(shù) break。 endenda1=c(1,1:k)。a2=c(2,1:k)。b1=c(3,1:k)。b2=c(4,1:k)。ea1=e(1,1:k)。ea2=e(2,1:k)。eb1=e(3,1:k)。 %輸出最終估計參數(shù)的矩陣，包含所有估計值eb2=e(4,1:k)。 %（由于矩陣初始值為0，若不改變矩陣值，則仍為0，所以，在循環(huán)跳出時，記錄循環(huán)次數(shù)，不輸出后續(xù)的0）c,e %輸出辨識參數(shù)矩陣和辨識誤差矩陣figure(2)。i=1:k。plot(i,a1,39。r39。,i,a2,39。:39。,i,b1,39。g39。,i,b2,39。y39。) %第二個圖形，用于繪制所估計的參數(shù)全部估計過程gtext(39。\leftarrow39。) %點擊鼠標(biāo)，標(biāo)注左箭頭gtext(39。a139。) %點擊鼠標(biāo)，標(biāo)注a1曲線gtext(39。\leftarrow39。) %點擊鼠標(biāo)，標(biāo)注左箭頭gtext(39。a239。) %點擊鼠標(biāo)，標(biāo)注a2曲線gtext(39。\leftarrow39。)gtext(39。b139。)gtext(39。\leftarrow39。)gtext(39。b239。)figure(3)。i=1:k。plot(i,ea1,39。r39。,i,ea2,39。:39。,i,eb1,39。g39。,i,eb2,39。y39。) %第三個圖形，繪制辨識的誤差title(39。Identification Precision39。) %為圖形注釋標(biāo)題gtext(39。\leftarrow39。) %點擊鼠標(biāo)，標(biāo)注左箭頭gtext(39。ea139。) %點擊鼠標(biāo)，標(biāo)注ea1曲線gtext(39。\leftarrow39。)gtext(39。ea239。)gtext(39。\leftarrow39。)gtext(39。eb139。)gtext(39。\leftarrow39。)gtext(39。eb239。)增廣最小二乘法clc。u=[,,,,,,,,]。y=zeros(1,30)。 %給輸出矩陣賦初值ys=zeros(1,30)。 %ys為不考慮噪聲時的輸出矩陣，為其賦初值ym=zeros(1,30)。 %ym為模型輸出矩陣(不考慮擇噪聲)ymd=zeros(1,30)。 %ymd為模型輸出矩陣(考慮擇噪聲)e0=sqrt()+randn(1,30) %產(chǎn)生隨機序列e0for k=3:30 y(k)=*y(k1)*y(k2)+*u(k1)+*u(k2)+e0(k)+*e0(k1)+*e0(k2)。 %產(chǎn)生輸出矩陣(帶噪聲) ys(k)=*y(k1)*y(k2)+*u(k1)+*u(k2)。 %產(chǎn)生輸出矩陣(不帶噪聲)endfigure(1)。 %第一個圖形，繪制出輸入、輸出的根徑圖subplot(3,1,1)stem(u)subplot(3,1,2)i=1:30。plot(i,y)subplot(313)stem(y)grid on。u,y,ys %輸出輸入、輸出矩陣c0=[ ]39。 %c0為待估計參數(shù)的初始值，辨識參數(shù)為7個，所以為7行p0=10^6*eye(7,7)。 %p0為遞推公式中的變量，設(shè)置其初始值為單位陣E=。 %設(shè)定誤差允許值c=[c0,zeros(7,29)]。 %c為辨識參數(shù)矩陣，設(shè)置其初始值e=zeros(7,30)。 %e為辨識誤差矩陣，該句用于設(shè)置其初值for k=3:30 x1=[y(k1) y(k2) u(k1) u(k2) e0(k) e0(k1) e0(k2)]39。 h=1+x139。*p0*x1。k1=inv(h)。 %該兩句為利用遞推公式K(N+1)=1/[1+x39。(N+1)*P(N)*x(N+1)] h1=k1*p0*x1。 d1=y(k)x139。*c0。 c1=c0+h1*d1。 %前面三句是使用遞推公式計算待估計的參數(shù)：c(K+1)=c(N)+K(N+1)*P(N)*x(N+1)*[y(N+1)x39。(N+1)*c(N)] ym(k)=[y(k1) y(k2) u(k1) u(k2)]*[c1(1)。c1(2)。c1(3)。c1(4)]。 %產(chǎn)生不考慮噪聲時的模型輸出 ymd(k)=x139。*c1。 %產(chǎn)生考慮噪聲時的模型輸出 e1=c1c0。 %求估計參數(shù)的誤差 e2=e1./c0。 e(:,k)=e2。 %完善誤差矩陣e c0=c1。 %為下次循環(huán)做準(zhǔn)備 c(:,k)=c1。 %每次循環(huán)都使待估計參數(shù)矩陣增加一列，記錄每一次的估計值 p1=p0h1*x139。*p0。 %利用遞推公式：P(N+1)=P(N)K(N+1)*P(N)*x(N+1)*x39。(N+1)P(N)來計算P(N+1) p0=p1。 %為下次循環(huán)做準(zhǔn)備 if e2=E %判斷估計參數(shù)與實際參數(shù)之間的誤差是否滿足要求 m=k。 %若滿足要求，則跳出循環(huán)，并且輸出停止估計的次數(shù) break。 endendc,e