(完整)系統(tǒng)辨識(shí)大作業(yè)匯總,推薦文檔

1、一、 問(wèn)題描述考慮仿真對(duì)象：z(k ) + 0.9z(k -1) + 0.15z(k - 2) + 0.02z(k - 3) = 0.7u(k -1) - 0.15u(k - 2) + e(k )e(k ) +1.0e(k -1) + 0.41e(k - 2) = av(k ), v n (0,1)式中，u(k)和 z(k)是輸入輸出數(shù)據(jù)，v(k)是零均值、方差為 1 的不相關(guān)的隨機(jī)噪聲；u(k)采用與 e(k)不相關(guān)的隨機(jī)序列。1. 設(shè)計(jì)實(shí)驗(yàn)，產(chǎn)生輸入輸出數(shù)據(jù)；2. 使用基本最小二乘法估計(jì)參數(shù)；3. 考慮其他適用于有色噪聲的辨識(shí)方法估計(jì)參數(shù)；4. 模型驗(yàn)證。二、 問(wèn)題分析對(duì)于單輸入單輸出系統(tǒng)

2、（single input single output, siso），如圖 1所示，將待辨識(shí)的系統(tǒng)看作“灰箱”，它只考慮系統(tǒng)的輸入輸出特性，而不強(qiáng)調(diào)系統(tǒng)的內(nèi)部機(jī)理。圖1 中，輸入 u(k)和輸出 z(k)是可以測(cè)量的， g(z-1) 是系統(tǒng)模型，用來(lái)描述系統(tǒng)的輸入輸出特性，y(k)是系統(tǒng)的實(shí)際輸出。 n (z-1) 是噪聲模型，v(k)是均值為零的不相關(guān)隨機(jī)噪聲，e(k)是有色噪聲。圖 1 siso 系統(tǒng)的“灰箱”結(jié)構(gòu)對(duì)于 siso 隨機(jī)系統(tǒng)，被辨識(shí)模型g(z) 為：y(z)b z-1 + b z-2 +l+ b z-ng(z) = u(z)=121+-1a z+ a z-2n+l+ a z-

3、n其相應(yīng)的差分方程為12nn ny(k ) = - ai y(k - i) + biu(k - i)i=1i=1若考慮被辨識(shí)系統(tǒng)或觀測(cè)信息中含有噪聲，被辨識(shí)模型可改寫為nnz(k ) = - ai y(k - i) + biu(k - i) + v(k )i=1i=1式中，z(k)為系統(tǒng)輸出量的第 k 次觀測(cè)值；y(k)為系統(tǒng)輸出量的第 k 次真值，y(k-1)為系統(tǒng)輸出量的第 k-1 次真值，以此類推；u(k)為系統(tǒng)的第 k 個(gè)輸入值，u(k-1)為系統(tǒng)的第 k-1 個(gè)輸入值；v(k)為均值為 0 的不相關(guān)隨機(jī)噪聲。1. 數(shù)據(jù)生成本部分需要生成系統(tǒng)的輸入輸出數(shù)據(jù)以及噪聲數(shù)據(jù)。1） 白噪聲的生

4、成辨識(shí)數(shù)據(jù)通常包含有噪聲，如果該噪聲相關(guān)性較弱或者強(qiáng)度很小，可近似看作白噪聲。本次實(shí)驗(yàn)問(wèn)題中v(t) : n (0,1) ，即服從標(biāo)準(zhǔn)正態(tài)分布，可以將噪聲看作為服從正態(tài)分布的白噪聲過(guò)程，在 matlab 中可以由 randn 函數(shù)生成。2） 輸入數(shù)據(jù)的生成偽隨機(jī)二進(jìn)制序列（pseudo random binary sequence, prbs）是廣泛應(yīng)用的一種偽隨機(jī)序列，所謂“二進(jìn)制”是指序列中每個(gè)隨機(jī)變量只有“0”或“1”兩種邏輯狀態(tài)。偽隨機(jī)二進(jìn)制序列可由多級(jí)線性反饋移位寄存器組成的隨機(jī)信號(hào)發(fā)生器產(chǎn)生，其中具有最長(zhǎng)循環(huán)周期的線性移位寄存器序列是偽隨機(jī)二進(jìn)制序列最常見(jiàn)的一種形式，簡(jiǎn)稱 m 序列

5、（maximal length sequence）。m 序列由于具有近似白噪聲的性質(zhì)，而且工程上易于實(shí)現(xiàn)，能夠保證較好的系統(tǒng)辨識(shí)精度，是普遍采用的一種辨識(shí)用輸入信號(hào)。圖 2 線性反饋移位寄存器產(chǎn)生偽隨機(jī)二進(jìn)制序列結(jié)構(gòu)圖以一個(gè) 4 級(jí)線性反饋移位寄存器產(chǎn)生偽隨機(jī)二進(jìn)制序列為例，如圖2所示。假設(shè) 4 個(gè)移位寄存器a0 , a1, a2 , a3 輸出的初態(tài)非全零，移位寄存器的工作原理是：一個(gè)移位脈沖來(lái)到后，每級(jí)移位寄存器的輸出移到下一級(jí)移位寄存器作為輸入，最末一級(jí)移位寄存器的輸出即為偽隨機(jī)二進(jìn)制序列。3）輸出數(shù)據(jù)的生成根據(jù)給定的 siso 系統(tǒng)，可以求出 z(k)的表達(dá)式：z(k ) = -0.9

6、z(k -1) - 0.15z(k - 2) - 0.02z(k - 3) + 0.7u(k -1)-0.15u(k - 2) + av(k ) -1.0e(k -1) - 0.41e(k - 2)其理想系數(shù)值為a1 = 0.9, a2 = 0.15, a3 = 0.02, b1 = 0.7, b2 = -1.5, c1 = 1.0, c2 = 0.41.可以根據(jù)生成的白噪聲序列和輸入序列，以及必要的 0 初始值，帶入表達(dá)式即可得到采樣輸出數(shù)據(jù)。2. 差分模型階檢驗(yàn)在實(shí)際場(chǎng)景中，辨識(shí)模型的階數(shù)和純時(shí)延往往是未知的，在很多情況下僅僅依靠猜測(cè)。在模型的階數(shù)和純時(shí)延不確定時(shí)，設(shè)系統(tǒng)模型為nny(t)

7、 = - ai y(t - i) + bj y(t - i) +a(t)i=1j =1其中 n 為模型的階數(shù)，a(t) = c(z-1 )e(t)模型的階估計(jì)可以采用多種方法，本實(shí)驗(yàn)采用比較簡(jiǎn)單易行的損失函數(shù)檢驗(yàn)法。定義預(yù)報(bào)誤差（噪聲方差的估計(jì)值）的平方和為損失函數(shù)，即n j= 1 n (t)ni=1當(dāng) n 從小增大時(shí)， jn 應(yīng)隨之減小，當(dāng) n 增大到某一值時(shí)， jn 應(yīng)近似白噪聲過(guò)程。采用以下的檢驗(yàn)原則：在 n-1 這一點(diǎn)， jn 最后一次出現(xiàn)陡峭的下降，此后就近似地保持不變或只有微小的下降，則取n= n 。3. 參數(shù)辨識(shí)模型在系統(tǒng)辨識(shí)和參數(shù)估計(jì)領(lǐng)域中，最廣泛的估計(jì)方法時(shí)最小二乘法和極大似

8、 然估計(jì)法。最小二乘法作為一種最基本的估計(jì)方法應(yīng)用極為廣泛，其他的大多數(shù)估計(jì)算法都與最小二乘法有關(guān)。它既可用于動(dòng)態(tài)系統(tǒng)，也可用于靜態(tài)系統(tǒng)； 可用于線性系統(tǒng)，也可用于非線性系統(tǒng)；可用于離線估計(jì)，也可用于在線估計(jì)。在隨機(jī)環(huán)境下利用最小二乘法時(shí)，無(wú)須知道觀測(cè)數(shù)據(jù)的概率統(tǒng)計(jì)信息，而這種方法獲得的結(jié)果，卻有相當(dāng)好的統(tǒng)計(jì)性質(zhì)。最小二乘參數(shù)估計(jì)方法來(lái)源于數(shù)理統(tǒng)計(jì)的回歸分析，它能提供一個(gè)在最小 方差意義上與實(shí)驗(yàn)數(shù)據(jù)最好的你和的模型。該估計(jì)在一定條件下有最佳的統(tǒng)計(jì)特性，即它們是一致的，無(wú)偏的和有效的。最小二乘法時(shí)一個(gè)經(jīng)典的方法，概念簡(jiǎn)明，適應(yīng)范圍廣，在一些情況下，可得到與極大似然法一樣好的統(tǒng)計(jì)效果， 它能很方便

9、地與其他辨識(shí)算法建立關(guān)系。（1） 遞推最小二乘算法當(dāng)獲得一批數(shù)據(jù)后，可一次求得相應(yīng)的參數(shù)估計(jì)值，這樣處理問(wèn)題的方法 稱為一次完成算法或批處理算法。它在理論研究方面有許多方便之處，但當(dāng)矩陣的維數(shù)增加時(shí)，矩陣求逆運(yùn)算的計(jì)算量會(huì)急劇增加，將給計(jì)算機(jī)的計(jì)算速度和存儲(chǔ)量帶來(lái)負(fù)擔(dān)，而且不適合在線辨識(shí)，無(wú)法跟蹤參數(shù)隨時(shí)間變化的情況。為了減少計(jì)算量，減少數(shù)據(jù)在計(jì)算機(jī)中所占的存儲(chǔ)量，也為了實(shí)時(shí)地辨識(shí)出動(dòng)態(tài)系統(tǒng)的特性，在用最小二乘法進(jìn)行參數(shù)估計(jì)時(shí)，把它轉(zhuǎn)化成參數(shù)遞推的估計(jì)。參數(shù)遞推估計(jì)是指被辨識(shí)的系統(tǒng)，每取得一次新的測(cè)量數(shù)據(jù)后，就在前一次估計(jì)結(jié)果的基礎(chǔ)上，利用新引入的測(cè)量數(shù)據(jù)對(duì)前一次估計(jì)的結(jié)果進(jìn)行修正， 從而遞推

10、地得出新的參數(shù)估計(jì)值。這樣，隨著新測(cè)量數(shù)據(jù)的引入，一次接一次地進(jìn)行參數(shù)估計(jì)，直到估計(jì)值達(dá)到滿意的精確程度為止。最小二乘遞推算法的基本思想可以概括為：當(dāng)前的估計(jì)值a(k ) =上次估計(jì)值(k -1) +修正項(xiàng)即新的估計(jì)值a(k ) 是在舊的估計(jì)值(k -1) 的基礎(chǔ)上，利用新的觀測(cè)數(shù)據(jù)對(duì)舊的估計(jì)值進(jìn)行修正而成的?？紤]如下模型： a(z-1)z(k ) = b(z-1)u(k ) + v(k ) ，其中u(k ) ， z(k ) 分別是系統(tǒng)的輸入和輸出； v(k ) 是均值為零，方差為一的不相關(guān)白噪聲。且滿足： a(z-1) = 1+ a z-1 + a z-2 + a z-n b(z-1) =

11、b z-11 + b z-2 2 +. + b zn -m12mh(k ) = -z(k -1),., -z(k - n), u(k -1),.u(k - m)t令a= a , a ,.a , b , b ,.b 12n12m則使用加權(quán)最小二乘參數(shù)估計(jì)遞推算法（recursiveweightedleast squares, rwls）有：a(k ) =(k -1) + k (k )z(k ) - ht (k )(k -1)1 k (k ) = p(k -1)h(k )ht (k ) p(k -1)h(k ) +-1l(k )p(k ) = i - k (k )ht (k )p(k -1)h(k

12、 ) = -z(k -1),. - z(k - n), u(k -1),., u(k - m)t圖 3 最小二乘遞推辨識(shí)參數(shù)估計(jì)過(guò)程中信息的變換可以看出，取l(k ) = 1的時(shí)候，加權(quán)最小二乘估計(jì)就退化成了最小二乘參數(shù)估計(jì)的遞推算法（recursiveleastsquares,rls）。加權(quán)參數(shù) 1 可以在l(0,1 范圍內(nèi)選擇，如果選 1l(k )= 1 ，所有的采樣數(shù)據(jù)都是等同加權(quán)的，如果1l(k ) 0, v(k ) = z(k ) - ht (k )a(k -1) ，即可得遞推公式如下：a(k ) =(k -1) + k (k )v(k )k (k ) = p(k -1)h(k )h

13、t (k )p(k -1)h(k ) +1-1p(k ) = i - k (k )ht (k )p(k -1)a(k ) = z(k ) - ht (k )(k -1)h(k ) = -z(k -1),. - z(k - n), -u(k -1),. - u(k - m), -v(k -1),. - v(k - v)t一般說(shuō)來(lái)選取初值p(0) = ia(0) = 0增廣最小二乘遞推算法擴(kuò)充了最小二乘法的參數(shù)向量a和數(shù)據(jù)向量h(k ) 的維數(shù)，把噪聲模型的辨識(shí)同時(shí)考慮進(jìn)去。最小二乘法的許多結(jié)論對(duì)它都是適用的，但最小二乘法只能獲得模型的參數(shù)估計(jì)，對(duì)于有色噪聲，也就是噪聲模型必須用c(z-1)v(k

14、 ) 表示時(shí)，只能采用增廣最小二乘遞推算法，方可獲得無(wú)偏估計(jì)。這是最小二乘參數(shù)估計(jì)的遞推算法不可代替的。增廣最小二乘算法的流程如圖 5 所示。圖 4 遞推最小二乘參數(shù)辨識(shí)流程圖圖 5 增廣最小二乘算法流程圖三、matlab 仿真及運(yùn)行結(jié)果在本次實(shí)驗(yàn)中， l 取值為 1，則系統(tǒng)辨識(shí)對(duì)象為:z(k ) = -0.9z(k -1) - 0.15z(k - 2) - 0.02z(k - 3) + 0.7u(k -1) -1.5u(k - 2) +av(k ) -1.0e(k -1) - 0.41e(k - 2)1. 差分模型的階檢驗(yàn)在實(shí)際應(yīng)用場(chǎng)景中，模型的階數(shù)和純時(shí)延往往是未知的，在此種情況下， 一般

15、采用如下方法來(lái)確定模型的階數(shù)：固定一個(gè)階數(shù) n，當(dāng) n 從小增大時(shí)，定義損失函數(shù) jn 應(yīng)隨之減小，當(dāng) n 增大到某一值時(shí)， jn 應(yīng)近似白噪聲過(guò)程。在n-1 這一點(diǎn)， jn 最后一次出現(xiàn)陡峭的下降，此后就近似地保持不變或只有微小的下降，則取n= n 。使用 matlab 進(jìn)行模型階數(shù)確定代碼如下：%噪聲方差的估計(jì)值和損失函數(shù)法定階%clear close all; l=100;y1=1;y2=1;y3=1;y4=0;for i=1:l+5;x1=xor(y3,y4);%第一個(gè)移位寄存器的輸入信號(hào)x2=y1;%第二個(gè)移位寄存器的輸入信號(hào)x3=y2;%第三個(gè)移位寄存器的輸入信號(hào)x4=y3;%第四

16、個(gè)移位寄存器的輸入信號(hào)y(i)=y4;%第四個(gè)移位寄存器的輸出，即 m 序列，幅值為“0”或“1” if y(i)0.5;u(i)=-1;%m 序列的值為 1 時(shí)，辨識(shí)的輸入取-1elseu(i)=1;%m 序列為 0 時(shí)，輸入取 1endy1=x1;y2=x2;y3=x3;y4=x4;%移位寄存器的輸出序列endv=randn(l+5,1);%- ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)%ksi(2) = 0;ksi(1) = 0;%設(shè) ksi 的前兩個(gè)初始值為 0lamda = 1; for k=3:l+5;ksi(k) = lamda*v(k) -

17、 ksi(k-1) - 0.41*ksi(k-2);endz(3)=0;z(2)=0;z(1)=0;%取 z 的前三個(gè)初始值為零for k=4:l+5;z(k) = -0.9*z(k-1) - 0.15*z(k-2) - 0.02*z(k-3) + 0.7*u(k-1) - 1.5*u(k-2) + ksi(k);%理想辨識(shí)輸出采樣信號(hào)end% 模型階次 n=1for i=1:lh1(i,1)=z(i);h1(i,2)=u(i);end estimate1=inv(h1*h1)*h1*z(2:l+1);d1=(z(2:l+1)-h1*estimate1)*(z(2:l+1)-h1*estima

18、te1)/l; %噪聲方差的估計(jì)值aic1=l*log(d1)+4*1;% 模型階次 n=2for i=1:lh2(i,1)=z(i+1);h2(i,2)=z(i);h2(i,3)=u(i+1);h2(i,4)=u(i);end estimate2=inv(h2*h2)*h2*z(3:l+2);d2=(z(3:l+2)-h2*estimate2)*(z(3:l+2)-h2*estimate2)/l; %噪聲方差的估計(jì)值aic2=l*log(d2)+4*2;% 模型階次 n=3for i=1:lh3(i,1)=z(i+2);h3(i,2)=z(i+1);h3(i,3)=z(i);h3(i,4)=

19、u(i+2);h3(i,5)=u(i+1);h3(i,6)=u(i);end estimate3=inv(h3*h3)*h3*z(4:l+3);d3=(z(4:l+3)-h3*estimate3)*(z(4:l+3)-h3*estimate3)/l; %噪聲方差的估計(jì)值aic3=l*log(d3)+4*3;% 模型階次 n=4for i=1:lh4(i,1)=z(i+3);h4(i,2)=z(i+2);h4(i,3)=z(i+1);h4(i,4)=z(i);h4(i,5)=u(i+3);h4(i,6)=u(i+2);h4(i,7)=u(i+1);h4(i,8)=u(i);end estimat

20、e4=inv(h4*h4)*h4*z(5:l+4);d4=(z(5:l+4)-h4*estimate4)*(z(5:l+4)-h4*estimate4)/l; %噪聲方差的估計(jì)值aic4=l*log(d4)+4*4;% 模型階次 n=5 for i=1:lh5(i,1)=z(i+4);h5(i,2)=z(i+3);h5(i,3)=z(i+2);h5(i,4)=z(i+1);h5(i,5)=z(i);h5(i,6)=u(i+4);h5(i,7)=u(i+3);h5(i,8)=u(i+2);h5(i,9)=u(i+1);h5(i,10)=u(i);end estimate5=inv(h5*h5)*

21、h5*z(6:l+5);d5=(z(6:l+5)-h5*estimate5)*(z(6:l+5)-h5*estimate5)/l; %噪聲方差的估計(jì)值aic5=l*log(d5)+4*5;plot(1:5,aic1 aic2 aic3 aic4 aic5)grid ontitle(損失函數(shù)法判斷模型階次) xlabel(階次)ylabel(損失函數(shù))程序運(yùn)行結(jié)果如圖 6 所示：圖 6 模型階次的判定由圖 6 可知，當(dāng)系統(tǒng)階次從 1 到 5 變化時(shí)，在 n=3 及以后，損失函數(shù)只有微小的下降，在 n=4 以后又出現(xiàn)上升的趨勢(shì)，所以取 n=3 作為模型的階次。2. 參數(shù)辨識(shí)（1）遞推最小二乘法根據(jù)

22、前文所述遞推最小二乘法算法流程，編寫遞推最小二乘法程序如下：% 遞推最小二乘辨識(shí)clcclear all%采用 m 序列方法產(chǎn)生輸入數(shù)據(jù)l=90;%90 個(gè)數(shù)據(jù)namenda=1; y1=1;y2=1;y3=1;y4=0;u_f=2; for i=1:lx1=xor(y3,y4);x2=y1; x3=y2; x4=y3; y(i)=y4;if y(i)0.5u(i)=-u_f; elseu(i)=u_f;end y1=x1;y2=x2;y3=x3;y4=x4;end v=1*randn(l,1); eb(2)=0;eb(1)=0;for k=3:leb(k)=-eb(k-1)-0.41*eb(

23、k-2)+namenda*v(k);endz(3)=0;z(2)=0;z(1)=0;for k=4:l;z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+eb(k);end figure(1) i=1:30;subplot(2,1,1) stem(u(i);title(輸入 m 序列)grid on subplot(2,1,2) stem(z(i);title(輸出序列)grid on;%c0=0.001 0.001 0.001 0.001 0.001;p0=103*eye(5,5);e= 0.0000005;c=c0

24、,zeros(5,l-1); e=zeros(5,l);%for k=4:lh1=-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2);k1=p0*h1*inv(h1*p0*h1+1); c1=c0+k1*(z(k)-h1*c0); e1=c1-c0;e2=e1./c0; e(:,k)=e2; c0=c1; c(:,k)=c1;p1(:,:,k)=p0-k1*k1*h1*p0*h1+1; p0=p1(:,:,k);if abs(e2)0.5;u(i)=-1;%m 序列的值為 1 時(shí)，辨識(shí)的輸入為-1 elseu(i)=1;%m 序列的值為 0 時(shí)，辨識(shí)的輸入為 1endy

25、1=x1;y2=x2;y3=x3;y4=x4;%移位寄存器的輸出序列,為下一次的輸入信號(hào)做準(zhǔn)備end figure(1);subplot(3,1,1);stem(u),grid on%畫出 m 序列輸入信號(hào)徑線圖并給圖形加上網(wǎng)格title(輸入信號(hào)（m 序列）)%-隨機(jī)噪聲的產(chǎn)生%v=randn(1,l);%產(chǎn)生一組 l=60 個(gè)正態(tài)分布的噪聲subplot(3,1,2);plot(v),grid on%畫出隨機(jī)噪聲信號(hào)title(隨機(jī)噪聲信號(hào) e(k)%-實(shí)際噪聲的產(chǎn)生%- ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)%ksi(2) = 0;ksi(1)

26、= 0;%設(shè) ksi 的前兩個(gè)初始值為 0lamda = 1; for k=3:l;ksi(k) = lamda*v(k) - ksi(k-1) - 0.41*ksi(k-2);end subplot(3,1,3)plot(ksi), grid ontitle(實(shí)際噪聲信號(hào)xi(k)%-辨識(shí)過(guò)程%z(3)=0;z(2)=0;z(1)=0;zm(3)=0;zm(2)=0;zm(1)=0;%輸出采樣，不考慮噪聲時(shí)系統(tǒng)輸出，不考慮噪聲時(shí)模型輸出，模型輸出的初值c0=0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001; %直接給出辨識(shí)被辨識(shí)參數(shù)的初始值，即一

27、個(gè)充分小的實(shí)向量p0=106*eye(8,8);%直接給出初始向量 p0，即一個(gè)充分大的實(shí)數(shù)單位矩陣c=c0,zeros(8,l-1);%被辨識(shí)參數(shù)參數(shù)矩陣的初始值及大小e=zeros(8,l);%相對(duì)誤差的初始值及大小for k=4:l;%開始求 kz(k) = -0.9*z(k-1) - 0.15*z(k-2) - 0.02*z(k-3) + 0.7*u(k-1) - 1.5*u(k-2) + ksi(k);%系統(tǒng)在 m 序列下的輸出采樣信號(hào)h1=-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),v(k),ksi(k-1),ksi(k-2);%為求 k（k）做準(zhǔn)備

28、x=h1*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%k d1=z(k)- h1*c0;c1=c0+k1*d1;%辨識(shí)參數(shù) ce1=c1-c0;e2=e1./c0;%求參數(shù)誤差的相對(duì)變化e(:,k)=e2;c0=c1;%給下一次用c(:,k)=c1;%把遞推出的辨識(shí)參數(shù) c 的列向量加入辨識(shí)參數(shù)矩陣p1=p0-k1*k1*h1*p0*h1+1;% find p(k)p0=p1;%給下次的計(jì)算使用end%整個(gè)循環(huán)結(jié)束c,e,%顯示被辨識(shí)參數(shù)及參數(shù)收斂情況% 分離變量a1=c(1,:);a2=c(2,:);a3=c(3,:);b1=c(4,:);b2=c(5,:); %分離出 a

29、1,a2,a3,b1,b2d1=c(6,:);d2=c(7,:);d3=c(8,:);%分離出的 d1,d2,d3% 分離變量的收斂情況ea1=e(1,:);ea2=e(2,:);ea3=e(3,:);eb1=e(4,:);eb2=e(5,:); %分離出 a1,a2,b1,b2 的收斂情況ed1=e(6,:);ed2=e(7,:);ed3=e(8,:);%分離出 d1,d2,d3 的收斂情況figure(2);i=1:l;plot(i,a1,r,i,a2,r:,i,a3,r+,i,b1,b,i,b2,b:,i,d1,g,i,d2,g:,i,d3,g+) grid onlegend(a1,a2

30、,a3,b1,b2,c1,c2,c3)title(辨識(shí)曲線)figure(3); i=1:l;plot(i,ea1,r,i,ea2,r:,i,ea3,r+,i,eb1,b,i,eb2,b:,i,ed1,g,i,ed2,g:,i,ed3,g+)%畫出各個(gè)參數(shù)的收斂情況grid on legend(ea1,ea2,ea3,eb1,eb2,ec1,ec2,ec3) title(辨識(shí)誤差曲線)figure(4); i=1:l;plot(i,z(i),g), grid on;title(辨識(shí)系統(tǒng)的輸出響應(yīng))%畫出被辨識(shí)系統(tǒng)的輸出響應(yīng)程序運(yùn)行結(jié)果如圖 10 至圖 13 所示。圖 10 輸入數(shù)據(jù)及噪聲圖 1

31、1 辨識(shí)曲線圖 12 辨識(shí)誤差曲線圖 13 辨識(shí)的輸出響應(yīng)使用增廣最小二乘算法進(jìn)行系統(tǒng)辨識(shí)的結(jié)果如表 2 所示。表 2 增廣最小二乘算法的辨識(shí)結(jié)果參數(shù) 真值0.90.150.020.7-1.51.0-1.0-0.4估計(jì)值 0.90.150.020.7-1.51.0-1.0-0.4從上表中可以看出，對(duì)于含有有色噪聲的系統(tǒng)，采用增廣最小二乘法進(jìn)行系統(tǒng)辨識(shí)可以取得無(wú)偏估計(jì)。四、模型與誤差分析對(duì)系統(tǒng)進(jìn)行辨識(shí)后，需要對(duì)辨識(shí)所得的系統(tǒng)進(jìn)行進(jìn)一步的模型驗(yàn)證，以觀察所得的模型是否是對(duì)當(dāng)前觀測(cè)數(shù)據(jù)的最佳選擇。對(duì)模型進(jìn)行驗(yàn)證，可以采用事先生成的輸入數(shù)據(jù)，分別輸入到已知模型和經(jīng)過(guò)系統(tǒng)辨識(shí)所得的模型中，檢驗(yàn)所得模型是

32、否存在問(wèn)題。如果檢驗(yàn)失敗，可能存在的問(wèn)題是：辨識(shí)所用的一組數(shù)據(jù)包含的信息部足或所選的模型類不合適。對(duì)本實(shí)驗(yàn)中的已知對(duì)象及經(jīng)過(guò)系統(tǒng)辨識(shí)所得的模型，選取一組輸入信號(hào)為：u = 1, -1,1, -1,1,1, -1,1,1, -1,1,1, -1, -1, -1,1比較原參考對(duì)象的輸出以及辨識(shí)模型的輸出。具體代碼如下：%模型檢驗(yàn)%clear all close allu=1,-1,1,-1,1,1,-1,1,1,-1,1,1,-1,-1,-1,1;%系統(tǒng)輸入信號(hào)l=20z=zeros(1,l);x=zeros(1,l);y=zeros(1,l);v=randn(1,l);%-實(shí)際噪聲的產(chǎn)生%- k

33、si(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)%ksi(2) = 0;ksi(1) = 0;%設(shè) ksi 的前兩個(gè)初始值為 0lamda = 1; for k=3:l;ksi(k) = lamda*v(k) - ksi(k-1) - 0.41*ksi(k-2);endfor k=4:16z(k)=-0.9*z(k-1)-0.15*z(k-2)-0.02*z(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k)%真實(shí)模型的輸出y(k)=-1.2317*y(k-1)-0.4735*y(k-2)-0.0600*y(k-3)+0.8534*u(k-1)-1

34、.2373*u(k-2)%遞推最小二乘辨識(shí)模型輸出x(k)=-0.9*x(k-1)-0.15*x(k-2)-0.02*x(k-3)+0.7*u(k-1)-1.5*u(k-2)+ksi(k)%增廣最小二乘辨識(shí)模型的輸出end figure(1);subplot(3,1,1) stem(u),grid on title(系統(tǒng)輸入信號(hào))subplot(3,1,2) i=1:1:l;plot(i,z,r,i,y,b),grid onlegend(真實(shí)模型輸出,遞推最小二乘辨識(shí)模型) title(真實(shí)模型輸出和遞推最小二乘辨識(shí)模型輸出值)subplot(3,1,3) i=1:1:l;plot(i,z,r

35、,i,x,g),grid onlegend(真實(shí)模型輸出,增廣最小二乘辨識(shí)模型) title(真實(shí)模型輸出和增廣最小二乘辨識(shí)模型輸出值) figure(2);subplot(2,1,1)i=1:1:l;plot(i,z-y), grid on;title(真實(shí)模型與遞推最小二乘辨識(shí)模型的差值)subplot(2,1,2) i=1:1:lplot(i,z-x), grid on;title(真實(shí)模型與增廣最小二乘辨識(shí)模型的差值)程序運(yùn)行結(jié)果如圖 14 和圖 15 所示。圖 14 不同辨識(shí)模型輸出和真實(shí)系統(tǒng)輸出對(duì)比圖 15 真實(shí)模型與辨識(shí)模型輸出信號(hào)的差值由上圖并結(jié)合表 1 和表 2，分析可得表 1 遞推最小二乘算法的辨識(shí)結(jié)果參數(shù)真值0.90.150.020.7-1.5估計(jì)值0.93670.1591-0.01380.6