系統(tǒng)辨識大作業(yè)匯總

1、問題描述考慮仿真對象: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)=Hv(k),vN(0,1)式中，u(k)和z(k)是輸入輸出數(shù)據(jù)，v(k)是零均值、方差為1的不相關(guān)的隨機(jī)噪聲；u(k)采用與e(k)不相關(guān)的隨機(jī)序列。1 .設(shè)計實驗，產(chǎn)生輸入輸出數(shù)據(jù)；2 .使用基本最小二乘法估計參數(shù)；3 .考慮其他適用于有色噪聲的辨識方法估計參數(shù)；4 .模型驗證。問題分析對于單輸入單輸出系統(tǒng)(SingleInputSingleOutput,SISO)，如圖1所示，將待辨識的系統(tǒng)看作“灰箱”，

2、它只考慮系統(tǒng)的輸入輸出特性，而不強(qiáng)調(diào)系統(tǒng)的內(nèi)部機(jī)理。圖1中，輸入u(k)和輸出z(k)是可以測量的，G(z/)是系統(tǒng)模型，用來描述系統(tǒng)的輸入輸出特性，y(k)是系統(tǒng)的實際輸出。N(z)是噪聲模型，v(k)是均值為零的不相關(guān)隨機(jī)噪聲，e(k)是有色噪聲。N(Z,)H燈圖1SISO系統(tǒng)的“灰箱”結(jié)構(gòu)對于SISO隨機(jī)系統(tǒng)，被辨識模型G(z)為:G()_y(z)_叱+|+bnz-Zu(z)1a1za2z，HIanz其相應(yīng)的差分方程為nny(k)=-aiy(k-i)xbu(k-i)i1i1若考慮被辨識系統(tǒng)或觀測信息中含有噪聲，被辨識模型可改寫為nnz(k)-，aiy(k-i)、bu(k-i)v(k)i1

3、i4式中，z(k)為系統(tǒng)輸出量的第k次觀測值；y(k)為系統(tǒng)輸出量的第k次真值，y(k-1)為系統(tǒng)輸出量的第k-1次真值，以此類推；u(k)為系統(tǒng)的第k個輸入值，u(k-1)為系統(tǒng)的第k-1個輸入值；v(k)為均值為0的不相關(guān)隨機(jī)噪聲。1.數(shù)據(jù)生成本部分需要生成系統(tǒng)的輸入輸出數(shù)據(jù)以及噪聲數(shù)據(jù)。1)白噪聲的生成辨識數(shù)據(jù)通常包含有噪聲，如果該噪聲相關(guān)性較弱或者強(qiáng)度很小，可近似看作白噪聲。本次實驗問題中v(t)|_|N(0,1),即服從標(biāo)準(zhǔn)正態(tài)分布，可以將噪聲看作為服從正態(tài)分布的白噪聲過程，在Matlab中可以由randn函數(shù)生成。2)輸入數(shù)據(jù)的生成偽隨機(jī)二進(jìn)制序列(PseudoRandomBina

4、rySequence,PRB0是廣泛應(yīng)用的一種偽隨機(jī)序列，所謂“二進(jìn)制”是指序列中每個隨機(jī)變量只有“0”或“1”兩種邏輯狀態(tài)。偽隨機(jī)二進(jìn)制序列可由多級線性反饋移位寄存器組成的隨機(jī)信號發(fā)生器產(chǎn)生，其中具有最長循環(huán)周期的線性移位寄存器序列是偽隨機(jī)二進(jìn)制序列最常見的一種形式，簡稱M序列(MaximalLengthSequence。M序列由于具有近似白噪聲的性質(zhì)，而且工程上易于實現(xiàn)，能夠保證較好的系統(tǒng)辨識精度，是普遍采用的一種辨識用輸入信號。圖2線性反饋移位寄存器產(chǎn)生偽隨機(jī)二進(jìn)制序列結(jié)構(gòu)圖以一個4級線性反饋移位寄存器產(chǎn)生偽隨機(jī)二進(jìn)制序列為例，如圖2所示。假設(shè)4個移位寄存器a0,a1,a2,a3輸出的初

5、態(tài)非全零，移位寄存器的工作原理是：一個移位脈沖來到后，每級移位寄存器的輸出移到下一級移位寄存器作為輸入，最末一級移位寄存器的輸出即為偽隨機(jī)二進(jìn)制序列。3)輸出數(shù)據(jù)的生成根據(jù)給定的SISO系統(tǒng)，可以求出z(k)的表達(dá)式：z(k)二一0.9z(k-1)-0.15z(k-2)-0.02z(k-3)0.7u(k-1)-0.15u(k-2)v(k)-1.0e(k-1)-0.41e(k-2)其理想系數(shù)值為&-0.9,a2-0.15,a3=0.02,b=0.7,b2=-1.5,g-1.0,02=0.41.可以根據(jù)生成的白噪聲序列和輸入序列，以及必要的0初始值，帶入表達(dá)式即可得到采樣輸出數(shù)據(jù)。2 .差分模型階

6、檢驗在實際場景中，辨識模型的階數(shù)和純時延往往是未知的，在很多情況下僅僅依靠猜測。在模型的階數(shù)和純時延不確定時，設(shè)系統(tǒng)模型為nny(t)-aiy(t-i)八bjy(t-i)(t)i1j1其中n為模型的階數(shù)，口=C(z)e(t)模型的階估計可以采用多種方法，本實驗采用比較簡單易行的損失函數(shù)檢驗法。定義預(yù)報誤差(噪聲方差的估計值)的平方和為損失函數(shù)，即1、N、JN?(t)N7當(dāng)n從小增大時，Jn應(yīng)隨之減小，當(dāng)n增大到某一值時，Jn應(yīng)近似白噪聲過程。采用以下的檢驗原則：在n-1這一點(diǎn)，Jn最后一次出現(xiàn)陡峭的下降，此后就近似地保持不變或只有微小的下降，則取i?=n。3 .參數(shù)辨識模型在系統(tǒng)辨識和參數(shù)估計

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

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

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

10、n),u(k-1),.u(k-m)T?-a1,a2,.an,n,b2,.bm則使用加權(quán)最小二乘參數(shù)估計遞推算法(RecursiveWeightedLeastSquares,RWLS有：k)=九一1)K(k)z(k)hT(k)，(k1)T11K(k)=P(k-1)h(k)hT(k)p(k-1)h(k)+-上(k)P(k)=I-K(k)hT(k)P(k-1)Th(k)=-z(k-1),.-z(k-n),u(k-1),.,u(k-m)新觀測傷中)咐)展小二乘參&估計遞推算法設(shè)定值NQ圖3最小二乘遞推辨識參數(shù)估計過程中信息的變換可以看出，取A(k)=1的時候，加權(quán)最小二乘估計就退化成了最小二乘參數(shù)估計

11、的遞推算法(RecursiveLeastSquares,RLS。加權(quán)參數(shù)?？梢栽?0,1范A圍內(nèi)選擇，如果選，=1,所有的采樣數(shù)據(jù)都是等同加權(quán)的，如果,0.5;u(i)=-1;%M序列的值為1時，辨識的輸入取-1elseu(i)=1;%M序列為0時，輸入取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的前兩個初始值為0lamda=1;fork=3:L+5;ksi(k)=lamda*v(k)-ksi(

12、k-1)-0.41*ksi(k-2);endz(3)=0;z(2)=0;z(1)=0;%取z的前三個初始值為零fork=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);%!I想辨識輸出采樣信號end%模型階次n=1fori=1:LH1(i,1)=z(i);H1(i,2)=u(i);endestimate1=inv(H1*H1)*H1*z(2:L+1);D1=(z(2:L+1)-H1*estimate1)*(z(2:L+1)-H1*estimate1)/L;%噪聲方差的估計值A(chǔ)IC1=L*lo

13、g(D1)+4*1;%模型階次n=2fori=1:LH2(i,1)=z(i+1);H2(i,2)=z(i);H2(i,3)=u(i+1);H2(i,4)=u(i);endestimate2=inv(H2*H2)*H2*z(3:L+2);D2=(z(3:L+2)-H2*estimate2)*(z(3:L+2)-H2*estimate2)/L;%噪聲方差的估計值A(chǔ)IC2=L*log(D2)+4*2;%模型階次n=3fori=1:LH3(i,1)=z(i+2);H3(i,2)=z(i+1);H3(i,3)=z(i);H3(i,4)=u(i+2);H3(i,5)=u(i+1);endestimate3

14、=inv(H3*H3)*H3*z(4:L+3);噪聲方差的估計值D3=(z(4:L+3)-H3*estimate3)*(z(4:L+3)-H3*estimate3)/L;%AIC3=L*log(D3)+4*3;%模型階次n=4fori=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);endestimate4=inv(H4*H4)*H4*z(5:L+4),;D4=(z(5:L+4)-H4*estimate4)

15、*(z(5:L+4)-H4*estimate4)/L;%噪聲方差的估計值A(chǔ)IC4=L*log(D4)+4*4;%模型階次n=5fori=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);endestimate5=inv(H5*H5)*H5*z(6:L+5);D5=(z(6:L+5)-H5*estimate5)*(z(6:L+5)-H5*estima

16、te5)/L;%噪聲方差的估計值A(chǔ)IC5=L*log(D5)+4*5;plot(1:5,AIC1AIC2AIC3AIC4AIC5)gridontitle(損失函數(shù)法判斷模型階次)xlabel(階次)ylabel(損失函數(shù))-IIX程序運(yùn)行結(jié)果如圖6所示：Figure1ilEditI9ols以心t即VindcHylp它白為要一國|口口|口損失函數(shù)法判斷模型階就iao00-004208111160圖6模型階次的判定由圖6可知，當(dāng)系統(tǒng)階次從1到5變化時，在n=3及以后，損失函數(shù)只有微小的下降，在n=4以后又出現(xiàn)上升的趨勢，所以取n=3作為模型的階次2.參數(shù)辨識(1)遞推最小二乘法根據(jù)前文所述遞推最小

17、二乘法算法流程，編寫遞推最小二乘法程序如下:%遞推最小二乘辨識clcclearall標(biāo)用M序列方法產(chǎn)生輸入數(shù)據(jù)L=90;%90個數(shù)據(jù)namenda=1;y1=1;y2=1;y3=1;y4=0;u_f=2;fori=1:Lx1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5u(i)=-u_f;elseu(i)=u_f;endy1=x1;y2=x2;y3=x3;y4=x4;endv=1*randn(L,1);eb(2)=0;eb(1)=0;fork=3:Leb(k)=-eb(k-1)-0.41*eb(k-2)+namenda*v(k);endz(3)=

18、0;z(2)=0;z(1)=0;fork=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);endfigure(1)i=1:30;subplot(2,1,1)stem(u(i);title(輸入M序列)gridonsubplot(2,1,2)stem(z(i);title(輸出序列)gridon;%c0=0.0010.0010.0010.0010.001;p0=10A3*eye(5,5);E=0.0000005;c=c0,zeros(5,L-1);e=zeros(5,L);%fork=4:Lh1=-

19、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;p1(:,:,k)=p0-k1*k1*h1*p0*h1+1;p0=p1(:,:,k);ifabs(e2)0.5;u(i)=-1;%M序列的值為1時，辨識的輸入為-1elseu(i)=1;%M序列的值為0時，辨識的輸入為1endy1=x1;y2=x2;y3=x3;y4=x4;%移位寄存器的輸出序列，為下一次的輸入信號做準(zhǔn)備endfigure(1);subplot

20、(3,1,1);stem(u),gridon%畫出M序列輸入信號徑線圖并給圖形加上網(wǎng)格title(,輸入信號(M序列)1)%隨機(jī)噪聲的產(chǎn)生%v=randn(1,L);%產(chǎn)生一組L=60個正態(tài)分布的噪聲subplot(3,1,2);plot(v),gridon%畫出隨機(jī)噪聲信號實際噪聲的產(chǎn)生%ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)%ksi(2)=0;ksi(1)=0;%設(shè)ksi的前兩個初始值為0lamda=1;fork=3:L;ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2);endsubplot(3,1,3)plot(k

21、si),gridontitle(實際噪聲信號xi(k)%辨識過程%z(3)=0;z(2)=0;z(1)=0;zm(3)=0;zm(2)=0;zm(1)=0;%輸出采樣，不考慮噪聲時系統(tǒng)輸出，不考慮噪聲時模型輸出，模型輸出的初值c0=0.0010.0010.0010.0010.0010.0010.0010.001;%直接給出辨識被辨識參數(shù)的初始值，即一個充分小的實向量p0=10A6*eye(8,8);%c=c0,zeros(8,L-1);%e=zeros(8,L);直接給出初始向量P0,即一個充分大的實數(shù)單位矩陣被辨識參數(shù)參數(shù)矩陣的初始值及大小相對誤差的初始值及大小fork=4:L;%開始求Kz

22、(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序列下的輸出采樣信號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)備x=h1*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%Kd1=z(k)-h1*c0;c1=c0+k1*d1;%辨識參數(shù)ce1=c1-c0;e2=e1./c0;%jt參數(shù)誤差的相對變化e(:,k)=e2;c0=c1;%合下一次用c(:,k)=c1;咐巴遞推出的辨識參數(shù)c

23、的列向量加入辨識參數(shù)矩陣p1=p0-k1*k1*h1*p0*h1+1;%findp(k)p0=p1;%給下次的計算使用end%整個循環(huán)結(jié)束c,e,媼示被辨識參數(shù)及參數(shù)收斂情況%分離變量分離出a1,a2,a3,b1,b2分離出的d1,d2,d3分離出a1,a2,b1,b2的收斂情況分離出d1,d2,d3的收斂情況a1=c(1,:);a2=c(2,:);a3=c(3,:);b1=c(4,:);b2=c(5,:);%d1=c(6,:);d2=c(7,:);d3=c(8,:);%分離變量的收斂情況ea1=e(1,:);ea2=e(2,:);ea3=e(3,:);eb1=e(4,:);eb2=e(5,:

24、);%ed1=e(6,:);ed2=e(7,:);ed3=e(8,:);%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+)gridonIegend(a1,a2,a3,b1,b2,c1,c2,c3)title(辨識曲線)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+)%畫出各個參數(shù)的收斂情況gridonIegend(ea1,ea2,ea3,eb1,eb

25、2,ec1,ec2,ec3)title(辨識誤差曲線)figure(4);i=1:L;plot(i,z(i),g),gridon;title(辨識系統(tǒng)的輸出響應(yīng),)%畫出被辨識系統(tǒng)的輸出響應(yīng)程序運(yùn)行結(jié)果如圖10至圖13所示圖10輸入數(shù)據(jù)及噪聲圖11辨識曲線圖12辨識誤差曲線圖13辨識的輸出響應(yīng)使用增廣最小二乘算法進(jìn)行系統(tǒng)辨識的結(jié)果如表2所示表2增廣最小二乘算法的辨識結(jié)果4比真值0.90.150.020.7-1.51.0-1.0-0.4傳計值0.90.150.020.7-1.51.0-1.0-0.4從上表中可以看出，對于含有有色噪聲的系統(tǒng)，采用增廣最小二乘法進(jìn)行系統(tǒng)辨識可以取得無偏估計。四、模型

26、與誤差分析對系統(tǒng)進(jìn)行辨識后，需要對辨識所得的系統(tǒng)進(jìn)行進(jìn)一步的模型驗證，以觀察所得的模型是否是對當(dāng)前觀測數(shù)據(jù)的最佳選擇。對模型進(jìn)行驗證，可以采用事先生成的輸入數(shù)據(jù)，分別輸入到已知模型和經(jīng)過系統(tǒng)辨識所得的模型中，檢驗所得模型是否存在問題。如果檢驗失敗，可能存在的問題是：辨識所用的一組數(shù)據(jù)包含的信息部足或所選的模型類不合適。對本實驗中的已知對象及經(jīng)過系統(tǒng)辨識所得的模型，選取一組輸入信號為：u=1,-1,1,-1,1,1,-1,1,1,-1,1,1,-1,-1,-1,1比較原參考對象的輸出以及辨識模型的輸出。具體代碼如下：%O%M%clearallcloseallu=1,-1,1,-1,1,1,-1,

27、1,1,-1,1,1,-1,-1,-1,1;%系統(tǒng)輸入信號L=20z=zeros(1,L);x=zeros(1,L);y=zeros(1,L);v=randn(1,L);%實際噪聲的產(chǎn)生%ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2)%ksi(2)=0;ksi(1)=0;%設(shè)ksi的前兩個初始值為0lamda=1;fork=3:L;ksi(k)=lamda*v(k)-ksi(k-1)-0.41*ksi(k-2);endfork=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)%1實模型的輸出y(k)=-1.2317*y(k-1)-0.4735*y(k-2)-0.0600*y(k-3)+0.8534*u(k-1)-1.2373*u(k-2)%遞推最小二乘辨識模型輸出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)%曾廣最小二乘辨識模型的輸出endfigure(1);subplot(3,1,1)stem(u),gridontitle(系統(tǒng)輸入信號)subplot(3,1,