廣義遞推最小二乘辨識(shí)

1、北京工商大學(xué)系統(tǒng)辨識(shí)上機(jī)實(shí)驗(yàn)報(bào)告廣義遞推最小二乘辨識(shí)一、實(shí)驗(yàn)?zāi)康? 通過實(shí)驗(yàn)掌握廣義最小二乘辨識(shí)算法； 2 運(yùn)用MATLAB編程，掌握算法實(shí)現(xiàn)方法。二、實(shí)驗(yàn)原理廣義最小二乘法的基本思想是基于對(duì)數(shù)據(jù)先進(jìn)行一次濾波預(yù)處理，然后利用普通最小二乘法對(duì)濾波后的數(shù)據(jù)進(jìn)行辨識(shí)。如果濾波模型選擇得合適，對(duì)數(shù)據(jù)進(jìn)行了較好的白色化處理，那么直接利用普通最小二乘法就能獲得無偏一致估計(jì)。廣義最小二乘法所用的濾波模型實(shí)際上就是一種動(dòng)態(tài)模型，在整個(gè)迭代過程中不斷靠偏差信息來調(diào)整這個(gè)濾波模型，使它逐漸逼近于一個(gè)較好的濾波模型，以便對(duì)數(shù)據(jù)進(jìn)行較好的白色化處理，使模型參數(shù)估計(jì)稱為無偏一致估計(jì)。理論上說，廣義最小二乘法所用的動(dòng)態(tài)

2、模型經(jīng)過幾次迭代調(diào)整后，便可對(duì)數(shù)據(jù)進(jìn)行較好的白化處理，但是，當(dāng)過程的輸出噪信比比較大或模型參數(shù)比較多時(shí)，這種數(shù)據(jù)白色化處理的可靠性就會(huì)下降。此時(shí)，準(zhǔn)則函數(shù)可能出現(xiàn)多個(gè)局部收斂點(diǎn)，因而辨識(shí)結(jié)果可能使準(zhǔn)則函數(shù)收斂于局部極小點(diǎn)上而不是全局極小點(diǎn)上。這樣，最終的辨識(shí)結(jié)果往往也會(huì)是有偏的。其收斂速度比較慢，需要經(jīng)過多次迭代計(jì)算，才能得到較準(zhǔn)確的參數(shù)估計(jì)值。一般情況下，經(jīng)過多次迭代后，估計(jì)值便會(huì)收斂到穩(wěn)態(tài)值。但在某些情況下（如噪聲比較低時(shí)）存在局部極小值，估計(jì)值不一定收斂到準(zhǔn)則函數(shù)的全局極小值上。為了防止參數(shù)估計(jì)值收斂到局部極小值，最好選定初值接近最優(yōu)解，一般可以用最小二乘法的批處理估計(jì)值作為初值。如果系

3、統(tǒng)是時(shí)變的，或?yàn)榱丝朔?shù)據(jù)飽和現(xiàn)象，可以在兩次RLS算法中分別引進(jìn)遺忘因子。三、實(shí)驗(yàn)內(nèi)容 數(shù)據(jù)獲取實(shí)驗(yàn)數(shù)據(jù)按照表10-1，為二階線性離散系統(tǒng)的輸入輸出數(shù)據(jù) 數(shù)據(jù)處理為了提高辨識(shí)精度，實(shí)驗(yàn)者必須對(duì)原始數(shù)據(jù)進(jìn)行剔除壞數(shù)據(jù)、零均值化、工頻濾波等處理。實(shí)驗(yàn)進(jìn)行了白化濾波處理。 辨識(shí)算法利用處理過的數(shù)據(jù)（取適當(dāng)?shù)臄?shù)據(jù)長度），選擇某種辨識(shí)方法（如RLS遞推最小二乘法、RELS、RIV或RML等參數(shù)估計(jì)算法及F-檢驗(yàn)或AIC定階法），估計(jì)出模型參數(shù)和階次，同時(shí)分析辨識(shí)結(jié)果。本實(shí)驗(yàn)采用廣義遞推最小二乘法進(jìn)行系統(tǒng)辨識(shí)。三、廣義遞推最小二乘法（RLS）原理廣義最小二乘法是用迭代的松弛算法對(duì)最小二乘估計(jì)的一種改進(jìn)，

4、它的基本思想是引入一個(gè)白化濾波器，把相關(guān)噪聲轉(zhuǎn)換為白噪聲，基于對(duì)觀測(cè)數(shù)據(jù)先進(jìn)行一次濾波處理然后利用普通最小二乘法對(duì)濾波后的數(shù)據(jù)進(jìn)行辨識(shí)。廣義最小二乘法的計(jì)算步驟如下：1 給定初始條件：包括給定的輸入輸出數(shù)據(jù)或者產(chǎn)生的數(shù)據(jù)序列，初始狀態(tài)矩陣P0，被辨識(shí)參數(shù)的初始值（取一個(gè)充分小的實(shí)向量），濾波器參數(shù)與矩陣初值。2 利用式 計(jì)算濾波后的輸入輸出序列。3 對(duì)于二階離散系統(tǒng)，利用式構(gòu)造。4 利用 三個(gè)式子遞推計(jì)算辨識(shí)矩陣5 利用式計(jì)算，并根據(jù)構(gòu)造。6 利用三個(gè)式子遞推計(jì)算。7 返回第二步進(jìn)行迭代計(jì)算，直至獲得滿意的辨識(shí)結(jié)果。四、實(shí)驗(yàn)步驟 輸入輸出數(shù)據(jù)：u=1.147,0.201 ,-0.787 ,-1

5、.584 -1.052, 0.866 ,1.152 ,1.573 ,0.626, 0.433 . -0.958,0.810, -0.044, 0.947, -1.474, -0.719, -0.086, 1.099, 1.450, 1.151. 0.485,1.633, 0.043 ,1.326, 1.706, -0.340, 0.890,0.433, -1.177, -0.390. -0.982,1.435 ,-0.119 ,-0.769, -0.899, 0.882,-1.008, -0.844, 0.628,-0.679 . 1.541,1.375, -0.984, -0.582, 1.

6、609,0.090, -0.813, -0.428 ,-0.848,-0.410. 0.048,-1.099 ,-1.108 ,0.259,-1.627, -0.528, 0.203, 1.204,1.691, -1.235. -1.228,-1.267, 0.309,0.043, 0.043, 1.461, 1.585, 0.552,-0.601, -0.319. 0.744 0.829,-1.626, -0.127, -1.578, -0.822, 1.469, -0.379,-0.212, 0.178. 0.493 -0.056, -0.1294, 1.228, -1.606, -0.3

7、82, -0.229, 0.313,-0.161, -0.810. -0.277 0.983, -0.288, 0.846, 1.325, 0.723, 0.713, 0.643 0.463,0.786. 1.161,0.850, -1.349, -0.596, 1.512, 0.795, -0.713, 0.453,-1.604, 0.889. -0.938, 0.056, 0.829,-0.981, -1.232, 1.327, -0.681,0.114,-1.135, 1.284. -1.201 0.758, 0.590, -1.007, 0.390, 0.836,-1.52, -1.0

8、53, -0.083, 0.619. 0.840 -1.258, -0.354, 0.629, -0.242,1.680, -1.236, -0.803,0.537, -1.100. 1.417,-1.024, 0.671, 0.688,-0.123, -0.952, 0.232, -0.793,-1.138, 1.154. 0.206,1.196, 1.013,1.518, -0.553, -0.987, 0.167, -1.445,0.630, 1.255. 0.311,-1.726,0.975, 1.718, 1.360, 1.667, 1.111, 1.018, 0.078,-1.66

9、5. -0.760, 1.184, -0.614, 0.994, -0.089, 0.947, 1.706, -0.395, 1.222, -1.351. 0.231,1.425, 0.114, -0.689, -0.704, 1.070, 0.262, 1.610 ,1.489,-1.602. 0.020, -0.601, -0.020, -0.601, -0.235, 1.245, 1.226, -0.204, 0.926, -1.297;figure(1);stem(u)grid ontitle(圖1 輸入信號(hào))y=1.381 3.794 2.481 -0.280 -2.742 -1.5

10、54 2.129 2.691 3.427 2.199 . 1.679 -1.249 1.371 0.637 3.131 -0.819 0.235 1.262 2.849 3.374. 2.346 0.664 3.015 0.561 2.271 3.650 0.625 2.305 0.364 1.857. -0.912 -1.547 1.940 0.262 -0.379 -0.176 3.720 0.058 -0.752 1.983. -0.923 3.361 4.240 -0.074 -0.481 3.780 2.137 0.086 0.638 -0.971. -0.929 0.679 -0.

11、664 -0.433 1.570 -2.785 -1.153 0.819 3.484 4.091. -2.375 -2.561 -2.778 2.911 1.362 0.735 3.118 3.770 2.381 -0.812. -1.635 0.589 1.550 -3.410 -1.249 -3.692 -2.358 2.552 -0.228 0.554. 2.178 2.471 0.743 -0.004 2.504 -3.204 -1.800 -1.284 0.159 0.426. 0.059 0.395 2.371 -0.157 2.248 3.297 2.329 2.780 2.37

12、5 1.873. 2.411 3.928 2.846 -2.215 -1.104 3.460 2.883 0.245 -0.231 -2.963. 2.072 -0.845 -0.074 1.037 2.468 -3.679 2.149 -0.081 1.639 -1.291. 2.548 -1.681 2.307 2.227 -1.558 0.008 2.055 -1.102 -1.427 0.350. 2.736 2.965 -2.346 -1.510 0.809 -0.592 2.706 -1.941 2.275 2.802. -1.337 2.091 -2.585 0.013 1.21

13、7 0.691 -0.491 2.114 0.333 -0.482. 3.388 2.082 3.797 4.079 5.036 1.250 -1.019 -0.160 -3.201 1.161. 3.926 1.789 -2.703 2.125 5.054 4.678 5.236 -0.241 2.152 0.356 . -3.519 2.213 1.527 -1.206 2.151 0.264 1.595 2.864 -0.539 1.982. -3.104 -0.264 2.433 0.009 -1.360 -0.521 3.319 1.445 3.105 3.783. -1.973 -

14、0.138 -0.452 -0.586 -4.045 -1.743 2.577 3.849 0.367 1.324;初始值設(shè)置，包括被辨識(shí)參數(shù)的初始值、誤差序列以及濾波器參數(shù)初值；迭代循環(huán)，辨識(shí)參數(shù)更新，根據(jù)誤差調(diào)整濾波器參數(shù)，迭代計(jì)算被辨識(shí)參數(shù)，直至參數(shù)符合條件或迭代次數(shù)到。計(jì)算數(shù)據(jù)與圖形顯示，包括辨識(shí)參數(shù)辨識(shí)過程以及誤差的收斂情況。五、運(yùn)行結(jié)果顯示1 輸入數(shù)據(jù)：2 被辨識(shí)參數(shù)：辨識(shí)結(jié)果：ans = -0.4251 0.0116 1.8869 -0.54643 被辨識(shí)參數(shù)誤差的收斂情況：六、實(shí)驗(yàn)源程序%遞推的廣義最小二乘法進(jìn)行參數(shù)估計(jì)clear;close all;display(廣義遞推最

15、小二乘辨識(shí));u=1.147,0.201 ,-0.787 ,-1.584 -1.052, 0.866 ,1.152 ,1.573 ,0.626, 0.433 . -0.958,0.810, -0.044, 0.947, -1.474, -0.719, -0.086, 1.099, 1.450, 1.151. 0.485,1.633, 0.043 ,1.326, 1.706, -0.340, 0.890,0.433, -1.177, -0.390. -0.982,1.435 ,-0.119 ,-0.769, -0.899, 0.882,-1.008, -0.844, 0.628,-0.679 .

16、 1.541,1.375, -0.984, -0.582, 1.609,0.090, -0.813, -0.428 ,-0.848,-0.410. 0.048,-1.099 ,-1.108 ,0.259,-1.627, -0.528, 0.203, 1.204,1.691, -1.235. -1.228,-1.267, 0.309,0.043, 0.043, 1.461, 1.585, 0.552,-0.601, -0.319. 0.744 0.829,-1.626, -0.127, -1.578, -0.822, 1.469, -0.379,-0.212, 0.178. 0.493 -0.0

17、56, -0.1294, 1.228, -1.606, -0.382, -0.229, 0.313,-0.161, -0.810. -0.277 0.983, -0.288, 0.846, 1.325, 0.723, 0.713, 0.643 0.463,0.786. 1.161,0.850, -1.349, -0.596, 1.512, 0.795, -0.713, 0.453,-1.604, 0.889. -0.938, 0.056, 0.829,-0.981, -1.232, 1.327, -0.681,0.114,-1.135, 1.284. -1.201 0.758, 0.590,

18、-1.007, 0.390, 0.836,-1.52, -1.053, -0.083, 0.619. 0.840 -1.258, -0.354, 0.629, -0.242,1.680, -1.236, -0.803,0.537, -1.100. 1.417,-1.024, 0.671, 0.688,-0.123, -0.952, 0.232, -0.793,-1.138, 1.154. 0.206,1.196, 1.013,1.518, -0.553, -0.987, 0.167, -1.445,0.630, 1.255. 0.311,-1.726,0.975, 1.718, 1.360,

19、1.667, 1.111, 1.018, 0.078,-1.665. -0.760, 1.184, -0.614, 0.994, -0.089, 0.947, 1.706, -0.395, 1.222, -1.351. 0.231,1.425, 0.114, -0.689, -0.704, 1.070, 0.262, 1.610 ,1.489,-1.602. 0.020, -0.601, -0.020, -0.601, -0.235, 1.245, 1.226, -0.204, 0.926, -1.297;figure(1);stem(u)grid ontitle(圖1 輸入信號(hào))y=1.38

20、1 3.794 2.481 -0.280 -2.742 -1.554 2.129 2.691 3.427 2.199 . 1.679 -1.249 1.371 0.637 3.131 -0.819 0.235 1.262 2.849 3.374. 2.346 0.664 3.015 0.561 2.271 3.650 0.625 2.305 0.364 1.857. -0.912 -1.547 1.940 0.262 -0.379 -0.176 3.720 0.058 -0.752 1.983. -0.923 3.361 4.240 -0.074 -0.481 3.780 2.137 0.08

21、6 0.638 -0.971. -0.929 0.679 -0.664 -0.433 1.570 -2.785 -1.153 0.819 3.484 4.091. -2.375 -2.561 -2.778 2.911 1.362 0.735 3.118 3.770 2.381 -0.812. -1.635 0.589 1.550 -3.410 -1.249 -3.692 -2.358 2.552 -0.228 0.554. 2.178 2.471 0.743 -0.004 2.504 -3.204 -1.800 -1.284 0.159 0.426. 0.059 0.395 2.371 -0.

22、157 2.248 3.297 2.329 2.780 2.375 1.873. 2.411 3.928 2.846 -2.215 -1.104 3.460 2.883 0.245 -0.231 -2.963. 2.072 -0.845 -0.074 1.037 2.468 -3.679 2.149 -0.081 1.639 -1.291. 2.548 -1.681 2.307 2.227 -1.558 0.008 2.055 -1.102 -1.427 0.350. 2.736 2.965 -2.346 -1.510 0.809 -0.592 2.706 -1.941 2.275 2.802

23、. -1.337 2.091 -2.585 0.013 1.217 0.691 -0.491 2.114 0.333 -0.482. 3.388 2.082 3.797 4.079 5.036 1.250 -1.019 -0.160 -3.201 1.161. 3.926 1.789 -2.703 2.125 5.054 4.678 5.236 -0.241 2.152 0.356 . -3.519 2.213 1.527 -1.206 2.151 0.264 1.595 2.864 -0.539 1.982. -3.104 -0.264 2.433 0.009 -1.360 -0.521 3

24、.319 1.445 3.105 3.783. -1.973 -0.138 -0.452 -0.586 -4.045 -1.743 2.577 3.849 0.367 1.324;%用最小二乘遞推算法辨識(shí)參數(shù)：a1,a2,b1,b2c0=0.001 0.001 0.001 0.001; %被辨識(shí)參數(shù)的初始值直接給取，取一個(gè)充分小的實(shí)向量ce0=0.001;0.001; %濾波器C(q-1)的參數(shù)初值p0=106*eye(4,4); %初始狀態(tài)P0也采用直接取方式，直接給出初始狀態(tài)P0，取一個(gè)充分大的實(shí)數(shù)單位矩陣pe0=106*eye(2,2); %計(jì)算濾波器的P矩陣初值yf=zeros(1,2

25、00);yf(1)=0;yf(2)=0;uf=zeros(1,200);uf(1)=0;uf(2)=0;e=zeros(1,200); %殘差序列e(1)=0;e(2)=0; c=c0,zeros(4,199); %被辨識(shí)參數(shù)矩陣的初始值及大小e=zeros(4,200); %相對(duì)誤差-殘差的初始值及大小ce=zeros(2,200); %濾波器參數(shù)每次辨識(shí)的結(jié)果n=0; %用于統(tǒng)計(jì)遞推次數(shù)%對(duì)200組數(shù)據(jù)進(jìn)行參數(shù)估計(jì)for k=3:200; %開始遞推運(yùn)算，求K yf(k)=y(k)+ce0(1)*y(k-1)+ce0(2)*y(k-2); %濾波后的輸入序列 uf(k)=u(k)+ce0(

26、1)*u(k-1)+ce0(2)*u(k-2); %濾波后的輸出序列 h1=-yf(k-1),-yf(k-2),uf(k-1),uf(k-2); %求h(k)%濾波后的輸出輸出組成的辨識(shí)數(shù)據(jù) x=h1*p0*h1+1; x1=inv(x); k1=p0*h1*x1;%求出K的值 d1=yf(k)-h1*c0; c1=c0+k1*d1;%求被辨識(shí)參數(shù)c e1=c1-c0;%求參數(shù)當(dāng)前值與上一次的值的差值 e2=e1./c0;%求參數(shù)的相對(duì)變化 e(:,k)=e2; %把當(dāng)前相對(duì)變化的列向量加入誤差矩陣的最后一列 c0=c1; %新獲得的參數(shù)作為下一次遞推的舊參數(shù) c(:,k)=c1; %把當(dāng)前所

27、辨識(shí)參數(shù)的c1列向量加入辨識(shí)參數(shù)矩陣的最后一列 p1=p0-k1*k1*h1*p0*h1+1; %求p(k)值 p0=p1; %把當(dāng)前的p(k)值給下次用 %遞推算法的實(shí)質(zhì)就是c和P在原來的基礎(chǔ)上計(jì)算，因此要為下時(shí)刻的這兩個(gè)值做準(zhǔn)備 e(k)=y(k)-y(k-1),-y(k-2),u(k-1),u(k-2)*c1; %計(jì)算殘差e(k),注意這里的h的組成 he1=-e(k-1);-e(k-2); xe=1+he1*pe0*he1; xe1=inv(xe); ke1=pe0*he1*xe1; %求當(dāng)前的K矩陣 de1=e(k)-he1*ce0; ce1=ce0+ke1*de1; ce(:,k)

28、=ce1; %保存濾波器參數(shù) ee1=ce1-ce0; ee2=ee1./ce0; %模型參數(shù)的相對(duì)誤差 Ee(:,k)=ee1; ce0=ce1; %為下一時(shí)刻做準(zhǔn)備 pe1=pe0-ke1*ke1*he1*pe0*he1+1; %為下一時(shí)刻賦值 pe0=pe1; %為下一時(shí)刻P矩陣賦值 n=n+1; %完成一次遞推，統(tǒng)計(jì)值加1end %大循環(huán)結(jié)束c %顯示被辨識(shí)參數(shù)e %顯示辨識(shí)結(jié)果的收斂情況n %顯示收斂條件滿足時(shí)算法最終遞推次數(shù)a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);%分離參數(shù)A1=a1(200)A2=a2(200)B1=b1(200)B2=b2(200)figure(2);i=1:200;plot(i,a1,r,i,a2,:,i,b1,g,i,b2,:)%直角坐標(biāo)下的線性刻度曲線title(圖2 被辨識(shí)參數(shù))figure(3);i=1:200;plot(i,ea1,r,i,ea2,g,i,eb1,b,i,eb2,r:)%a