【正文】 plot(i,ymd) %繪制考慮噪聲時的模型輸出矩陣圖形title(39。)subplot(413)i=1:30。帶噪聲時的系統(tǒng)輸出39。)gtext(39。)gtext(39。)gtext(39。)gtext(39。)gtext(39。)gtext(39。) %點(diǎn)擊鼠標(biāo)，標(biāo)注左箭頭gtext(39。) %第三個圖形，繪制辨識的誤差title(39。,i,ed1,39。,i,eb2,39。,i,ea2,39。d239。d139。d039。b239。b139。a239。a139。辨識參數(shù)圖形39。,i,d2,39。,i,d0,39。,i,b1,39。plot(i,a1,39。ed1=e(6,1:k)。ea2=e(2,1:k)。d0=c(5,1:k)。 endendc,e %輸出辨識矩陣和辨識誤差矩陣ym,ymd %輸出帶噪聲的輸出矩陣和不帶噪聲的輸出矩陣a1=c(1,1:k)。 %利用遞推公式：P(N+1)=P(N)K(N+1)*P(N)*x(N+1)*x39。 %完善誤差矩陣e c0=c1。*c1。c1(2)。*c0。k1=inv(h)。 %c為辨識參數(shù)矩陣，設(shè)置其初始值e=zeros(7,30)。u,y,ys %輸出輸入、輸出矩陣c0=[ ]39。 %產(chǎn)生輸出矩陣(帶噪聲) ys(k)=*y(k1)*y(k2)+*u(k1)+*u(k2)。 %給輸出矩陣賦初值ys=zeros(1,30)。eb239。eb139。ea239。ea139。Identification Precision39。g39。r39。b239。b139。a239。a139。y39。:39。i=1:k。ea2=e(2,1:k)。a2=c(2,1:k)。 %利用遞推公式：P(N+1)=P(N)K(N+1)*P(N)*x(N+1)*x39。 %完善誤差矩陣e c0=c1。 %前面三句是使用遞推公式計算待估計的參數(shù)：c(K+1)=c(N)+K(N+1)*P(N)*x(N+1)*[y(N+1)x39。(N+1)*P(N)*x(N+1)] h1=k1*p0*x1。 h=1+x139。 %為p0賦初值，為遞推做準(zhǔn)備E=。 %第一個圖形，繪制出輸入、輸出的根徑圖subplot(3,1,1)stem(u)subplot(3,1,2)i=1:30。 %給輸出矩陣賦初值e0=randn(1,30)。*yL。 XL(k,:)=x1 %每次循環(huán)都增加XL的一行 yL(k,:)=y1 %每次循環(huán)都增加yL的一行endc1=XL39。 %生成輸出序列endsubplot(3,1,1)stem(u) %繪制輸入矩陣的根徑圖subplot(3,1,2)i=1:30。七、附錄（程序）最小二乘法一次實(shí)現(xiàn)：%最小二乘法一次實(shí)現(xiàn)u=[,,]。本實(shí)驗(yàn)完成過程中使用以下三種方法：最小二乘法一次完成、遞推最小二乘法、增廣矩陣法。（4） 設(shè)δk+f1δk=ε(k)，用廣義最小二乘法估計參數(shù)θ。二、 實(shí)驗(yàn)內(nèi)容設(shè)單輸入單輸出系統(tǒng)的差分方程為y（k）=a1yk1a2yk2+b1uk1+b2uk2+δ（k） δk=εk+a1εk1+a2ε（k2）取真實(shí)值θT=[a1 a2 b1 b2]=[ ] ，輸入數(shù)據(jù)如下表所列。k u(k)k u(k) k u(k)1 11 21 2 12 22 3 13 23 4 14 24 5 15 25 6 16 26 7 17 27 8 18 28 9 19 29 10 20 30 用θ的真實(shí)值利用查分方程求出yk作為測量值，εk為均值為0。（5） 用增廣矩陣法估計參數(shù)θ詳細(xì)分析和比較所獲得的參數(shù)辨識結(jié)果，并說明上述參數(shù)辨識方法的優(yōu)點(diǎn)。1. 在辨識方法上的比較：最小二乘法的一次實(shí)現(xiàn)是離線辨識，需要采集的數(shù)據(jù)量較大，一次完成計算；遞推最小二乘法為在線計算，在某個初始條件下啟動，按照遞推算法，隨著信息的不斷貨的而不斷修正模型參數(shù)的估計值；增光矩陣法為在線算法，可以處理噪聲不相關(guān)時的系統(tǒng)辨識問題。y=zeros(1,30)。plot(i,y) %繪制輸出矩陣subplot(3,1,3)stem(y) %繪制輸出矩陣的根徑圖grid on。*XL。c=c2*c3a1=c(1),a2=c(2),b1=c(3),b2=c(4) %輸出辨識參數(shù)遞推最小二乘法：%遞推最小二乘法實(shí)現(xiàn)參數(shù)估計clc。 %產(chǎn)生隨機(jī)誤差for k=3:16 %產(chǎn)生隨機(jī)噪聲 v(k)=e0(k)+*e0(k1)+*e0(k2)。plot(i,y)subplot(3,1,3)stem(y)grid on。 %規(guī)定誤差允許范圍c=[c0,zeros(4,29)]。*p0*x1。 d1=y(k)x139。(N+1)*c(N)] e1=c1c0。 %為下次循環(huán)做準(zhǔn)備 c(:,k)=c1。(N+1)P(N)來計算P(N+1) p0=p1。b1=c(3,1:k)。eb1=e(3,1:k)。plot(i,a1,39。,i,b1,39。) %第二個圖形，用于繪制所估計的參數(shù)全部估計過程gtext(39。) %點(diǎn)擊鼠標(biāo)，標(biāo)注a1曲線gtext(39。) %點(diǎn)擊鼠標(biāo)，標(biāo)注a2曲線gtext(39。)gtext(39。)figure(3)。,i,ea2,39。,i,eb2,39。) %為圖形注釋標(biāo)題gtext(39。) %點(diǎn)擊鼠標(biāo)，標(biāo)注ea1曲線gtext(39。)gtext(39。)gtext(39。)增廣最小二乘法clc。 %ys為不考慮噪聲時的輸出矩陣，為其賦初值ym=zeros(1,30)。 %產(chǎn)生輸出矩陣(不帶噪聲)endfigure(1)。 %c0為待估計參數(shù)的初始值，辨識參數(shù)為7個，所以為7行p0=10^6*eye(7,7)。 %e為辨識誤差矩陣，該句用于設(shè)置其初值for k=3:30 x1=[y(k1) y(k2) u(k1) u(k2) e0(k) e0(k1) e0(k2)]39。 %該兩句為利用遞推公式K(N+1)=1/[1+x39。 c1=c0+h1*d1。c1(3)。