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

1、躺二請(qǐng)分別用最小二乘一次完咸算法和遞推最小二乘法對(duì)以下數(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)分別測(cè)出n+N個(gè)輸出輸入值 y(1),y(2),y(n+N),u(1),u(2),u(n+N

2、),則可寫出N個(gè)方程，寫成向量一矩陣形式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 )測(cè)量矩陣。因此，式()是一個(gè)含有(2n+1 )個(gè)未知參數(shù)，由 N個(gè)方程組成的聯(lián)立方程組。1 1y在給定輸出向量y和測(cè)量矩陣的條件下求參數(shù)B的估計(jì)，這就是系統(tǒng)辨識(shí)問題。設(shè) ?表示的估計(jì)值，?表示的最優(yōu)估計(jì)，

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

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

5、E，則稱？是 無偏估計(jì)2) 一致性估計(jì)誤差)的方差矩陣為2彳1Var?E丄TNN()2? 1lim Var lim 一 NN N上式表明當(dāng)N is時(shí)，)以概率1趨近于0o當(dāng)gk)為不相關(guān)隨機(jī)序列時(shí)，最小二乘估計(jì)具有無偏性和一致性3) 有效性如果式()中的E的均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計(jì)值為有效值，參數(shù)估計(jì)偏差的方差達(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ù)估計(jì)值服從正態(tài)分布

6、,?: N2E T(4.1.2 )1.2遞推最小二乘法為了實(shí)現(xiàn)實(shí)時(shí)控制，必須采用遞推算法，這種辨識(shí)方法主要用于在線辨識(shí)設(shè)已獲得的觀測(cè)數(shù)據(jù)長度為N，則Yn(421)Nyn估計(jì)方差矩陣為式中Var2PnPNP NNY N如果再獲得1組新的觀測(cè)值，則又增加1個(gè)方程yN得新的參數(shù)估計(jì)值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)量，因而上式可以寫成最后，得最小二乘法辨識(shí)公式廣?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ù)估計(jì)。2兩種

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

9、0岀0輸-50 10 2030405060岀輸散離2.2遞推最小二乘法的實(shí)現(xiàn)遞推算法實(shí)現(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)行計(jì)算（參照式?？梢宰C明，當(dāng) Cis時(shí)，按照遞推公式算得的將與上面用批處理算法求得的結(jié)果相同，當(dāng)N很大時(shí),C的選擇便無關(guān)緊要了。3 ）遞推算法的停機(jī)標(biāo)準(zhǔn)：?. k ?1 k 1 |max 和 ，£是適當(dāng)小的數(shù)。4）最小二乘估計(jì)的遞推算法系統(tǒng)參數(shù)的步驟如下：(1) 產(chǎn)生系統(tǒng)輸入信號(hào) M序列，輸入系統(tǒng)階次，計(jì)算輸入和輸出數(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) 計(jì)算Kn 1PnN1 N1PnN11；(5) 輸入新測(cè)量數(shù)據(jù) y(N+1 )和u (N+1 );(6) 計(jì)算N 1NKn 1y(N 1)N1£(7) 計(jì)算FN 1IKN 1N 1 PN;(8) 輸出? 1和Pn 1 ；(9) N - N+1,轉(zhuǎn)(3 );程序編制思路：Q(1 )產(chǎn)生M序列，通過計(jì)算，依據(jù)算出輸出值，設(shè)定遞推初始值，采用424方法2設(shè)定辨識(shí)參數(shù)初值，?為一個(gè)很小

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

12、算法辨識(shí)的Malab程序流程圖Parameter Ide ntificati on with Recursive Least Squares Method224程序運(yùn)行曲線0.50.40.30.20.10-0.1-0.2-0.3-0.4-0.50102030405060708090100圖6M序列信號(hào)的波形1.61.41.210.80.60.40.20-0.2-0.40102030405060708090100Identification Precision1400120010008006004002000-200-4000102030405060708090100測(cè)試結(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ù)估計(jì)的誤差逐漸減少，最終趨于穩(wěn)定。下面是辨識(shí)誤 差的分布趨勢(shì)。(誤差=新估計(jì)值 舊估計(jì)值)【5舊估計(jì)值6結(jié)論最小二乘法是應(yīng)用廣泛的辨識(shí)方法之一。可用于動(dòng)態(tài)系統(tǒng)，也可用于靜態(tài)系統(tǒng)；可用于線性系統(tǒng)，也可用 于非線性系統(tǒng)。通過本課程設(shè)計(jì)，了解和掌握了基于最小二乘法原理的一次性完成算法和遞推算法的實(shí)現(xiàn)方 法，完成了 Matlab下的仿真計(jì)算，能夠精確的估計(jì)辨識(shí)參數(shù)。最小二乘一次性完成算法是離線算法，需要 采集大量數(shù)據(jù)，一次性完成計(jì)算，因此，數(shù)據(jù)計(jì)算量大，當(dāng)數(shù)據(jù)量很大時(shí)，數(shù)據(jù)輸入不方便，但在本課程設(shè) 計(jì)過程當(dāng)中，考慮到了此問題，運(yùn)用相應(yīng)的辦法，解

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

16、金工業(yè)，1999.3 黃道平，MATLAB與控制系統(tǒng)的數(shù)字仿真及CAD化學(xué)工業(yè)，2004.104 侯媛彬汪梅王立奇系統(tǒng)辨識(shí)及其MATLAB仿真科學(xué)2004.25 蔡季冰，系統(tǒng)辨識(shí)，理工大學(xué)，1991.3附錄附錄1 M序列產(chǎn)生程序function Out=Mserise(le ngth,ple ngth，a)%plength=4;length=15;%M序列周期=2Aplength-1，% 置M序列總長度和級(jí)數(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%異或運(yùn)算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程序運(yùn)行結(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)辨識(shí)的輸入信號(hào)為一個(gè)周期的M序列，信號(hào)長度L;幅度A ；z=zeros(1,16);%定義輸出觀測(cè)值的長度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、用理想輸出值作為觀測(cè)值endsubplot(3,1,1) %畫三行一列圖形窗口中的第一個(gè)圖形stem(u) %畫出輸入信號(hào)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)%畫三行一列圖形窗口中的第三個(gè)圖形stem(z),grid on%畫出輸出觀測(cè)值 z的經(jīng)線圖形，并顯示坐標(biāo)網(wǎng)格ylabel('離散輸出 z(k)')%u,z;%顯示輸入信號(hào)和輸出觀測(cè)信號(hào)%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%計(jì)算參數(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 %計(jì)算并顯示%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;%四個(gè)移位積存器的輸出初始值for i=1:L;%開始循環(huán)，長度為 Lx仁xor(y

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

24、03 ”end%小循環(huán)結(jié)束y1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號(hào)做準(zhǔn)備end%大循環(huán)結(jié)束，產(chǎn)生輸入信號(hào)ufigure(1);%第1個(gè)圖形stem(u),grid on%以徑的形式顯示出輸入信號(hào)并給圖形加上網(wǎng)格z(2)=0;z(1)=0;% 取z的前兩個(gè)初始值為零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);%給出理想的辨識(shí)輸出采樣信號(hào)end%RLS遞推最小二乘辨識(shí)c0=0.001 0.001 0.001 0.001'%直接給出被辨識(shí)參數(shù)的初始值，即一個(gè)充分小的實(shí)向量p0=10A6*eye(4,4);%直接給出初始狀態(tài)P0，即一個(gè)充分大的實(shí)數(shù)單位矩陣E=0.000000005;% 相對(duì)誤差 E=0.000000005c=c0,zeros(4,99);%被辨識(shí)參數(shù)矩陣的初始值及大小e=zeros(4,100);%相對(duì)誤差的初始值及大小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*