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

1、1最小二乘法的理論基礎(chǔ)1.1最小二乘法設(shè)單輸入單輸出線性定長系統(tǒng)的差分方程表示為：其中(k)為服從N(0,1)的隨機(jī)噪聲,現(xiàn)分別測出n+N個輸出輸入值y(1),y(2),y(n+N),u(1),u(2),u(n+N),則可寫出N個方程，寫成向量矩陣形式(4.1.1)則式(1.1.1)可寫為 （4.1.2）式中：y為N維輸出向量；為N為維噪聲向量；為（2n+1）維參數(shù)向量；為N×（2n+1）測量矩陣。因此，式（4.1.1）是一個含有（2n+1）個未知參數(shù)，由N個方程組成的聯(lián)立方程組。在給定輸出向量y和測量矩陣的條件下求參數(shù)的估計，這就是系統(tǒng)辨識問題。設(shè) 表示 的估計值，表示y的最優(yōu)估計

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

3、）如果矩陣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. 白噪聲序列1.1.2最小二乘估計的概率性質(zhì)最小二乘估計的概率性質(zhì)1） 無偏性由于輸出值y是隨機(jī)的，所以是隨機(jī)的，但注意不是隨機(jī)值。如果，則稱是無偏估計2）一致性估計誤差的方差矩陣為（4.1.9）上式表明當(dāng)N時，以概率1趨近于。當(dāng)(k)為不相關(guān)隨機(jī)序列時，最小二乘估計具有無偏性和一致性3)有

4、效性如果式（1.1.2）中的的均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計值 為有效值，參數(shù)估計偏差的方差達(dá)到Cramer-Rao不等式的下界，即（4.1.10）式中M為Fisher矩陣，且（4.1.11）4）漸近正態(tài)性如果式（4.1.2）中的是均值為0且服從正態(tài)分布的白噪聲向量，則最小二乘參數(shù)估計值服從正態(tài)分布，即（4.1.2）1.2遞推最小二乘法為了實現(xiàn)實時控制，必須采用遞推算法，這種辨識方法主要用于在線辨識（4.2.1）設(shè)已獲得的觀測數(shù)據(jù)長度為N，則估計方差矩陣為式中于是如果再獲得1組新的觀測值，則又增加1個方程得新的參數(shù)估計值式中應(yīng)用矩陣求逆引理，可得和的遞推關(guān)系式矩陣求逆引

5、理 設(shè)A為n×n矩陣，B和C為n×m矩陣，并且A，ABCT和I+CTA-1B都是非奇異矩陣，則有矩陣恒等式（4.2.2）得到遞推關(guān)系式由于是標(biāo)量，因而上式可以寫成（4.2.3）最后，得最小二乘法辨識公式（4.2.4）有2種給出初值的辦法（1）設(shè)N0（N0>n）為N的初始值，可算出（4.2.5）（2）假定C是充分大的常數(shù)，I為(2n+1) ×(2n+1)單位矩陣，則經(jīng)過若干次遞推之后能得到較好的參數(shù)估計。1242兩種算法的實現(xiàn)方案2.1最小二乘法一次完成算法實現(xiàn)如果把式（1.1.6）中的取好足夠的輸入輸出數(shù)據(jù)以后一次計算出來，那么這種算法就式最小二乘法的一次完

6、成法。2.1.1最小二乘一次完成算法程序框圖賦輸入信號初值u定義輸出觀測值的長度并計算系統(tǒng)的輸出值畫出輸入和輸出觀測值的圖形給樣本矩陣HL和zL賦值計算參數(shù)從中分離出并顯示出被辨識參數(shù)a1, a2, b1, b2停機(jī)圖3.2 最小二乘一次完成算法程序框圖2.1.2一次完成法程序具體程序參見附錄42.1.3一次完成算法程序運(yùn)行結(jié)果ans = 1.5000 0.7000 1.0000 0.5000a1 = 1.5000a2 = 0.7000b1 =1.00002.1.4辨識數(shù)據(jù)比較1.50000.70001.0000.5000參數(shù)真值1.5000.7001.000.50誤差00002.1.5程序運(yùn)

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

8、新測量數(shù)據(jù)y ( N+1 )和u（N+1）；（6）計算（7）計算（8）輸出和；（9）NN+1,轉(zhuǎn)（3）；2.2.2程序編制思路：（1）產(chǎn)生M序列，通過計算，依據(jù)算出輸出值，設(shè)定遞推初始值，采用4.2.4方法2設(shè)定辨識參數(shù)初值，為一個很小的量，P0為一個很大值的4×4矩陣，設(shè)定誤差量，作為停機(jī)標(biāo)準(zhǔn)。（2）依據(jù)設(shè)定，計算，可以算出為一標(biāo)量，依據(jù)式4.2.5算出K1,K1為4×1矩陣。（3）依據(jù)式4.2.5計算出，和，加入估計參數(shù)矩陣；（4）計算誤差大小，求出相對變化率，加入誤差矩陣第一列。（5）再以和為已知參數(shù)，重復(fù)（2）（3）計算出參數(shù)，并加入估計參數(shù)矩陣。（6）如此循環(huán)，計

9、算出參數(shù)估計。（7）判斷誤差變化率滿足要求否，如果滿足，則退出循壞，不再進(jìn)行計算。（8）分離估計參數(shù)，繪出圖形。2.2.3遞推最小二乘法程序框圖如下所示：具體程序參見附錄Y工作間清零產(chǎn)生輸出采樣信號給被辨識參數(shù)和P賦初值按照式(4.2.4)的第三式計算P(k)計算被辨識參數(shù)的相對變化量參數(shù)收斂滿足要求？停機(jī)按照式(4.2.4)的第二式計算K(k)按照式(4.2.4)的第一式計算(k)第四個移位寄存器的輸出取反，并將幅值變?yōu)?.03得到辨識系統(tǒng)的輸入信號樣本值給M序列的長度L和移位寄存器的輸入賦初始值畫出被辨識參數(shù)的各次遞推估計值圖形分離參數(shù)畫出被辨識參數(shù)的相對誤差的圖形畫出辨識的輸入信號徑線圖

10、形圖5 最小二乘遞推算法辨識的Malab程序流程圖2.2.4程序運(yùn)行曲線圖6M序列信號的波形2.2.5測試結(jié)果見附錄62.2.6地退數(shù)據(jù)表遞推次數(shù)誤差誤差誤差誤差11.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

11、-1.500.70.701.01.000.50.5071.51.500.70.701.01.000.50.50從以上數(shù)據(jù)可以看出，隨著遞推次數(shù)的增加，參數(shù)估計的誤差逐漸減少，最終趨于穩(wěn)定。下面是辨識誤差的分布趨勢。（誤差=）56結(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ù)輸入不方便，但

12、在本課程設(shè)計過程當(dāng)中，考慮到了此問題，運(yùn)用相應(yīng)的辦法，解決了矩陣輸入的問題。遞推算法適合于在線算法，利用原有參數(shù)估計進(jìn)行下一步估計，可以做到運(yù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è)

13、出版社,2006.7 .2劉宏才,系統(tǒng)辨識與參數(shù)估計,北京,冶金工業(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(length,plength,a)%plength=4;length=15; %M序列周期=2plength-1 ,%置M序列總長度和級數(shù) ，a是幅度 for n=1:plength %移位寄存器輸入X(i)初T態(tài)（1111）

14、X(n)=1;endfor i=1:length for j=plength:-1:2 X(j)=X(j-1); end X(1)=xor(X(plength-1),X(plength); %異或運(yùn)算 if X(plength)=0 Out(i)=-1; else Out(i)=X(plength); end end%save ('mdata','Out');%繪圖i1=ik=1:1:length;plot(k,Out,k,Out,'rx')xlabel('k')ylabel('M序列')title('移位

15、寄存器產(chǎn)生的M序列')end附錄4程序運(yùn)行結(jié)果如下圖所示圖7移位寄存器產(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 number; 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+1 z(k)=-1.5

16、*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想輸出值作為觀測值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; %顯示輸入信號和

17、輸出觀測信號%L=L-1%數(shù)據(jù)長度for i=n+1:L+1 h=-z(i-1) -z(i-2) u(i-1) u(i-2); dd(i-2,:)=h; %dd=(L-1)*2nend%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(1

18、1) -z(10) u(11) u(10);-z(12) -z(11) u(12) 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+1 zz=z(j); gg(j-2,:)=zz;end %c=inv(dd'*dd)*dd'*gg; %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(16)% 給樣本矩陣zL賦值%c

19、alculating parameters%計算參數(shù)c1=dd'*dd; c2=inv(c1); c3=dd'*gg; c=c2*c3;c'%c1=dd'*dd; c2=inv(c1); c3=dd'*ZL; c=c2*c3 %計算并顯示 %DISPLAY PARAMETERSa1=c(1), a2=c(2), b1=c(3), b2=c(4); %從 中分離出并顯示a1 、a2、 b1、 b2%End附錄3遞推最小二乘算法程序clear;clc;%清理工作間變量L=100;% M序列的周期y1=1;y2=1;y3=1;y4=0;%四個移位積存器的輸出初

20、始值for i=1:L;%開始循環(huán)，長度為L x1=xor(y3,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)

21、M序列的值為"0"時,辨識的輸入信號取“0.03” end%小循環(huán)結(jié)束 y1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號做準(zhǔn)備end%大循環(huán)結(jié)束，產(chǎn)生輸入信號u figure(1);%第1個圖形stem(u),grid on%以徑的形式顯示出輸入信號并給圖形加上網(wǎng)格z(2)=0;z(1)=0;%取z的前兩個初始值為零for k=3:L;%循環(huán)變量從3到25 z(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.00

22、1'%直接給出被辨識參數(shù)的初始值,即一個充分小的實向量p0=106*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; %開始求K h1=-z(k-1),-z(k-2),u(k-1),u(k-2)' x=h1'*p0*h1+1; x1=inv(x); %開始求K(k) k1=p0*h1*x1;%求出K的值 d1=z(k)-h1'*c0; c1=c0+k1*d1;%求被辨識參數(shù)c e1=c1-c0;%求參數(shù)當(dāng)前值與上一次的值的差值