應(yīng)用最小二乘一次完成法和遞推最小二乘法算法地系統(tǒng)辨識

1、躺二請分別用最小二乘一次完咸算法和遞推最小二乘法對以下數(shù)學(xué)模型進(jìn)行參j(Ar) + 烏鞏上一 1)-+一2)二既疋(上-1)+-2)+ r(k)式中，理想的系數(shù)值：=1.5)尿=07，= 1 _0.= 0.5玖此為服從 的隨機(jī)噪聲.輸入處)采用4階陽序列幅值為5,權(quán) 陣 W = I °1最小二乘法的理論基礎(chǔ)1.1最小二乘法設(shè)單輸入單輸出線性定長系統(tǒng)的差分方程表示為：y ka1yk 1a2y k2 Lanyknb)u kb|Uk 1L bnuk n k其中3(k)為服從N(0,1)的隨機(jī)噪聲，現(xiàn)分別測出n+N個輸出輸入值 y(1),y(2),y(n+N),u(1),u(2),u(n+N

2、),則可寫出N個方程，寫成向量一矩陣形式y(tǒng) nLy1u n 1Lu 1y n 2y n 1Ly2u n 2Lu 2MMMMy n Ny n N 1LyNu n NLu Na1Mn 1ann 2b0MMnNbnanbon 1n 2Mn N(4.1.1)bny N unN L u N貝U式可寫為y( )式中：y為N維輸出向量； 訪N為維噪聲向量；B為(2n+1 )維參數(shù)向量；為N x(2n+1 )測量矩陣。因此，式()是一個含有(2n+1 )個未知參數(shù)，由 N個方程組成的聯(lián)立方程組。1 1y在給定輸出向量y和測量矩陣的條件下求參數(shù)B的估計，這就是系統(tǒng)辨識問題。設(shè) ?表示的估計值，?表示的最優(yōu)估計，

3、則有y(4.1.3 )式中：? n1M? n2?M,bo1? nNMbn設(shè) e(k)=y(k)-?(k), e(k)稱為殘差，則有 e=y- ?=y- B最小二乘估計要求殘差的平方和最小，即按照指數(shù)函數(shù)J eT ey求j對的偏導(dǎo)數(shù)并令其等于 0可得:由式（4.1.5 ）可得的最小二乘估計:(414 )(4.1.5 )(4.1.6 )j為極小值的充分條件是:2J即矩陣T為正定矩陣，或者說是非奇異的。1.1.1最小二乘法估計中的輸入信號當(dāng)矩陣T的逆陣存在是，式（1.1.6 ）才有解。一般地，如果u（k）是隨機(jī)序列或偽隨機(jī)二位式序列，則矩陣T是非奇異的，即（T）一1存在，式（1.1.6 ）有解?，F(xiàn)在

4、從T必須正定出發(fā)，討論對u（k）的要求。n N 1Tyyyuk nuyuu(4.1.7 )當(dāng)N足夠大時有W.P.1RyRuyRyuRu(4.1.8 )如果矩陣T正定，則Ru是是對稱矩陣，并且各階主子式的行列式為正。當(dāng)N足夠大時，矩陣 Ru才是是對稱的。由此引出矩陣 T正定的必要條件是 u(k)為持續(xù)激勵信號。如果序列u(k)的n+1階方陣Ru是正定的，則稱u(k)的n+1階持續(xù)激勵信號。下列隨機(jī)信號都能滿足Ru正定的要求1. 有色隨機(jī)信號2. 偽隨機(jī)信號3. 白噪聲序列最小二乘估計的概率性質(zhì)最小二乘估計的概率性質(zhì)1 )無偏性由于輸出值y是隨機(jī)的，所以)是隨機(jī)的，但注意B不是隨機(jī)值。如果E ?

5、E，則稱？是 無偏估計2) 一致性估計誤差)的方差矩陣為2彳1Var?E丄TNN()2? 1lim Var lim 一 NN N上式表明當(dāng)N is時，)以概率1趨近于0o當(dāng)gk)為不相關(guān)隨機(jī)序列時，最小二乘估計具有無偏性和一致性3) 有效性如果式()中的E的均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計值為有效值，參數(shù)估計偏差的方差達(dá)到 Cramer-Rao 不等式的下界，即Var2E T 1)式中M為Fisher矩陣，且In p y/M Eln p y /(4.1.11 )4）漸近正態(tài)性如果式（4.1.2 ）中的E是均值為0且服從正態(tài)分布的白噪聲向量,則最小二乘參數(shù)估計值服從正態(tài)分布

6、,?: N2E T(4.1.2 )1.2遞推最小二乘法為了實現(xiàn)實時控制，必須采用遞推算法，這種辨識方法主要用于在線辨識設(shè)已獲得的觀測數(shù)據(jù)長度為N，則Yn(421)Nyn估計方差矩陣為式中Var2PnPNP NNY N如果再獲得1組新的觀測值，則又增加1個方程yN得新的參數(shù)估計值nyn1yN式中PNPn1應(yīng)用矩陣求逆引理，可得Pn i和Pn的遞推關(guān)系式矩陣求逆引理設(shè)A為n Xn矩陣,B和C為n Xm矩陣，并且A, A + BCT 和 I+CTA-1B都是非奇異矩陣，則有矩陣恒等式A BCtA 1 A 1B I CtA 1B 1CtA 1(422 )得到遞推關(guān)系式P NPNN 1 1N PN1 R

7、i由于 n Pn n 1是標(biāo)量，因而上式可以寫成最后，得最小二乘法辨識公式廣?N 1?NKnT ?1yN 1N 1 N(4.2.3 )YK N 1PnN 11N PN N 1IPN 1PnPnT 1 TN 1 1N 1 FNlN 1N 1PNpn 1P nRj n 1 1N PNN 1TN1 PN(4.2.4 )有2種給出初值的辦法(1 )設(shè) NO ( N0>n )為N的初始值，可算出FN oT1NO NO,oPn on oYn o(4.2.5 )(2)假定? O,F02c I ,C是充分大的常數(shù)，I 為(2n+1)X(2n+1)單位矩陣，則經(jīng)過若干次遞推之后能得到較好的參數(shù)估計。2兩種

8、算法的實現(xiàn)方案2.1最小二乘法一次完成算法實現(xiàn)如果把式（）中的？取好足夠的輸入一輸出數(shù)據(jù)以后一次計算出來，那么這種算法就式最小二乘法的一 次完成法。最小二乘一次完成算法程序框圖圖3.2最小二乘一次完成算法程序框圖一次完成法程序具體程序參見附錄4一次完成算法程序運行結(jié)果ans = al =1.50000.70001.00000.50001.5000 a2 =0.7000 b1 =1.0000辨識數(shù)據(jù)比較?aab?b21.50000.70001.0000.5000參數(shù)真值1.5000.7001.000.50誤差00002.1.5程序運行曲線1 k 入0 輸-15101520253035404550

9、0岀0輸-50 10 2030405060岀輸散離2.2遞推最小二乘法的實現(xiàn)遞推算法實現(xiàn)步驟1 ）輸入M序列的產(chǎn)生程序，可以參見3.2當(dāng)中M序列產(chǎn)生方法（具體程序見附錄。）2 ）初始化。2一種簡單的方法是選擇? 0和Pn Cl，其中C為盡量大的數(shù)，然后以它們?yōu)槠鹗贾颠M(jìn)行計算（參照式?？梢宰C明，當(dāng) Cis時，按照遞推公式算得的將與上面用批處理算法求得的結(jié)果相同，當(dāng)N很大時,C的選擇便無關(guān)緊要了。3 ）遞推算法的停機(jī)標(biāo)準(zhǔn)：?. k ?1 k 1 |max 和 ，£是適當(dāng)小的數(shù)。4）最小二乘估計的遞推算法系統(tǒng)參數(shù)的步驟如下：(1) 產(chǎn)生系統(tǒng)輸入信號 M序列，輸入系統(tǒng)階次，計算輸入和輸出數(shù)據(jù)

10、u (i), y (i), i=1,2,n;(2) 置 N=n , ? 0, PN 10101 (I 單位矩陣)；(3) 形成行向量N1 y(n) L y(N 1 n) u(N) L u(N 1 nTT1(4) 計算Kn 1PnN1 N1PnN11；(5) 輸入新測量數(shù)據(jù) y(N+1 )和u (N+1 );(6) 計算N 1NKn 1y(N 1)N1£(7) 計算FN 1IKN 1N 1 PN;(8) 輸出? 1和Pn 1 ；(9) N - N+1,轉(zhuǎn)(3 );程序編制思路：Q(1 )產(chǎn)生M序列，通過計算，依據(jù)算出輸出值，設(shè)定遞推初始值，采用424方法2設(shè)定辨識參數(shù)初值，?為一個很小

11、的量，P0為一個很大值的4 X4矩陣，設(shè)定誤差量，作為停機(jī)標(biāo)準(zhǔn)。(2 )依據(jù)設(shè)定，計算1y(2)y(1) u(1) u(2) T,T1 Po11可以算出為一標(biāo)量，依據(jù)式425算出K1,K1為4 X1矩陣。(3 )依據(jù)式425計算出?，和P，加入估計參數(shù)矩陣；(4)計算誤差大小，求出相對變化率，加入誤差矩陣第一列。(5 )再以?和P為已知參數(shù)，重復(fù)(2)-( 3)計算出參數(shù) ?，并加入估計參數(shù)矩陣。(6)如此循環(huán)，計算出參數(shù)估計?。(7 )判斷誤差變化率滿足要求否，如果滿足，則退出循壞，不再進(jìn)行計算。(8 )分離估計參數(shù)，繪出圖形。遞推最小二乘法程序框圖如下所示:具體程序參見附錄圖5最小二乘遞推

12、算法辨識的Malab程序流程圖Parameter Ide ntificati on with Recursive Least Squares Method224程序運行曲線0.50.40.30.20.10-0.1-0.2-0.3-0.4-0.50102030405060708090100圖6M序列信號的波形1.61.41.210.80.60.40.20-0.2-0.40102030405060708090100Identification Precision1400120010008006004002000-200-4000102030405060708090100測試結(jié)果見附錄6地退數(shù)據(jù)表遞

13、推次數(shù)a1誤差印?a2a2誤;差02?bb2誤差b?誤;差bzb211.500.0011.4990.70.0010.6991.00.0010.9990.50.0010.499921.501.50.700.71.001.00.500.531.50.0011.4990.70.0010.6991.00.750.250.50.750.2541.5-1.500.70.0010.6991.00.750.250.50.750.2551.5-1.500.70.701.00.750.250.50.750.2561.5-1.500.70.701.01.000.50.5071.51.500.70.701.01.00

14、0.50.50從以上數(shù)據(jù)可以看出，隨著遞推次數(shù)的增加，參數(shù)估計的誤差逐漸減少，最終趨于穩(wěn)定。下面是辨識誤 差的分布趨勢。(誤差=新估計值 舊估計值)【5舊估計值6結(jié)論最小二乘法是應(yīng)用廣泛的辨識方法之一?？捎糜趧討B(tài)系統(tǒng)，也可用于靜態(tài)系統(tǒng)；可用于線性系統(tǒng)，也可用 于非線性系統(tǒng)。通過本課程設(shè)計，了解和掌握了基于最小二乘法原理的一次性完成算法和遞推算法的實現(xiàn)方 法，完成了 Matlab下的仿真計算，能夠精確的估計辨識參數(shù)。最小二乘一次性完成算法是離線算法，需要 采集大量數(shù)據(jù)，一次性完成計算，因此，數(shù)據(jù)計算量大，當(dāng)數(shù)據(jù)量很大時，數(shù)據(jù)輸入不方便，但在本課程設(shè) 計過程當(dāng)中，考慮到了此問題，運用相應(yīng)的辦法，解

15、決了矩陣輸入的問題。遞推算法適合于在線算法，利用 原有參數(shù)估計進(jìn)行下一步估計，可以做到運算量小，實時進(jìn)行估計，根據(jù)仿真結(jié)果圖示，可以看出計算一定 量的數(shù)據(jù)后，參數(shù)估計趨于穩(wěn)定，只要滿足誤差要求，不必計算完所有數(shù)據(jù)。兩種算法不足之處，在考慮噪 聲干擾時，看成了白噪聲，具有不相關(guān)性。如果噪聲引起的方程誤差是相關(guān)的，而不是白噪聲，可以通過下 面兩種方法進(jìn)行估計先估計未受噪音污染的系統(tǒng)輸出，然后再估計模型的參數(shù)，就是輔助變量法；使用 一種濾波器，將相關(guān)的方程的誤差轉(zhuǎn)換為白噪聲再進(jìn)行估計，就是增廣矩陣法。7參考文獻(xiàn)1 言俊科，系統(tǒng)辨識理論及應(yīng)用，國防工業(yè)，2006.7 .2 宏才,系統(tǒng)辨識與參數(shù)估計，冶

16、金工業(yè)，1999.3 黃道平，MATLAB與控制系統(tǒng)的數(shù)字仿真及CAD化學(xué)工業(yè)，2004.104 侯媛彬汪梅王立奇系統(tǒng)辨識及其MATLAB仿真科學(xué)2004.25 蔡季冰，系統(tǒng)辨識，理工大學(xué)，1991.3附錄附錄1 M序列產(chǎn)生程序function Out=Mserise(le ngth,ple ngth，a)%plength=4;length=15;%M序列周期=2Aplength-1，% 置M序列總長度和級數(shù) ，a是幅度for n=1:plength%移位寄存器輸入 X(i)初T態(tài)(1111 )X( n)=1;endfor i=1:le ngthfor j=ple ngth:-1:2精彩文檔X

17、(j)=X(j-1);end%異或運算X(1)=xor(X(ple ngth-1),X(ple ngth);if X(ple ngth)=OOut(i)=-1;elseOut(i)=X(ple ngth);endend%save ('mdata','Out');%繪圖i1=ik=1:1:le ngth;plot(k,Out,k,Out,'rx')xlabel('k')ylabel('M 序列')title('移位寄存器產(chǎn)生的M序列')end附錄4程序運行結(jié)果如下圖所示移位寄存器產(chǎn)生的ri序列圖7移位

18、寄存器產(chǎn)生的M序列附錄2最小二乘一次完成法%Oncefinish%model is y(k)+1.5y(k-1)+0.7(k-2)=u(k-1)+0.5(k-2)+v(k);L=50;A=1;n=2;%output nu mber;u=mseries(A,L);%u=Mserise(L,M,A); %系統(tǒng)辨識的輸入信號為一個周期的M序列，信號長度L;幅度A ；z=zeros(1,16);%定義輸出觀測值的長度1*16dd=zeros(1,4);gg=zeros(14,1);for k=n+1:L+1z(k)=-1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %

19、用理想輸出值作為觀測值endsubplot(3,1,1) %畫三行一列圖形窗口中的第一個圖形stem(u) %畫出輸入信號u的經(jīng)線圖形ylabel('輸入 u(k)')subplot(3,1,2)i=1:1:L+1;plot(i,z)ylabel('輸出 z(k)')subplot(3,1,3)%畫三行一列圖形窗口中的第三個圖形stem(z),grid on%畫出輸出觀測值 z的經(jīng)線圖形，并顯示坐標(biāo)網(wǎng)格ylabel('離散輸出 z(k)')%u,z;%顯示輸入信號和輸出觀測信號%L=L-1%數(shù)據(jù)長度for i=n +1:L+1h=-z(i-1)

20、-z(i-2) u(i-1) u(i-2);dd(i-2,:)=h;%dd=(L-1)*2 nend%HL=-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12)

21、u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14) %給樣本矩陣HL賦值for j=n +1:L+1zz=z( j); gg(j-2,:)=z z;end%c=i nv(dd'*dd)*dd'*gg;給樣本矩陣zL賦值%ZL=z (3) ;z (4) ;z (5) ;z (6);z 7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z (佝%calculat ing parameters%計算參數(shù)c1=dd'

22、*dd; c2=in v(c1); c3=dd'*gg; c=c2*c3;c'%"= dd'*dd; c2=i nv(c1); c3=dd'*ZL; c=c2*c3 %計算并顯示%DISPLAY PARAMETERS從中分離出并顯示 al、a2、 bl、 b2a1=c(1), a2=c (2), b1=c(3), b2=c(4); %End附錄3遞推最小二乘算法程序clear;clc;%清理工作間變量L=100;% M 序列的周期y1=1;y2=1;y3=1;y4=0;%四個移位積存器的輸出初始值for i=1:L;%開始循環(huán)，長度為 Lx仁xor(y

23、3,y4);%第一個移位積存器的輸入是第3個與第4個移位積存器的輸出的“或”x2=y1;%第二個移位積存器的輸入是第3個移位積存器的輸出x3=y2;%第三個移位積存器的輸入是第2個移位積存器的輸出x4=y3;%第四個移位積存器的輸入是第3個移位積存器的輸出y(i)=y4;%取出第四個移位積存器幅值為"0"和"1"的輸出信號，if y(i)>0.5,u(i)=-0.5;% 如果M序列的值為"1"時,辨識的輸入信號取“ -0.03 ”else u(i)=0.5;% 當(dāng)M序列的值為"0"時,辨識的輸入信號取“ 0.

24、03 ”end%小循環(huán)結(jié)束y1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號做準(zhǔn)備end%大循環(huán)結(jié)束，產(chǎn)生輸入信號ufigure(1);%第1個圖形stem(u),grid on%以徑的形式顯示出輸入信號并給圖形加上網(wǎng)格z(2)=0;z(1)=0;% 取z的前兩個初始值為零for k=3:L;%循環(huán)變量從3到25z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%給出理想的辨識輸出采樣信號end%RLS遞推最小二乘辨識c0=0.001 0.001 0.001 0.001'%直接給出被辨識參數(shù)的初始值，即一個充分小的實向量p0=10A6*eye(4,4);%直接給出初始狀態(tài)P0，即一個充分大的實數(shù)單位矩陣E=0.000000005;% 相對誤差 E=0.000000005c=c0,zeros(4,99);%被辨識參數(shù)矩陣的初始值及大小e=zeros(4,100);%相對誤差的初始值及大小for k=3:L; % 開始求 Kh仁-z(k-1),-z(k-2),u(k-1),u(k-2)' x=h1'*p0*h1+1; x1=inv(x); %開始求 K(k)k1=p0*