系統(tǒng)辨識(shí)實(shí)驗(yàn)報(bào)告xjtu

1、系統(tǒng)辨識(shí)實(shí)驗(yàn)報(bào)告高海南碩037一、實(shí)驗(yàn)問題單輸入單輸出系統(tǒng)：A(zt)y(t) = B(zT)(t) + C(zT)e(t) 其中：A(z-i) = 1 - 2.851z-i + 2.717z-2 -0.865z-3B( z-1) = z-1 + z -2 + z -3C (z-1) = 1 + 0.7 z-1 + 0.22 z -2e(t)N(0,1)產(chǎn)生必要隨機(jī)數(shù)；設(shè)計(jì)實(shí)驗(yàn)，產(chǎn)生輸出數(shù)據(jù)；用最小二乘辨識(shí)系統(tǒng)模型；任用一種適合有色噪聲辨識(shí)算法辨識(shí)系統(tǒng)模型。二、實(shí)驗(yàn)原理遞推最小二乘遞推最小二乘辨識(shí)算法可減少運(yùn)算量和數(shù)據(jù)在計(jì)算機(jī)中所占的存儲(chǔ)量，同時(shí) 也能實(shí)時(shí)辨識(shí)出系統(tǒng)的特性，將最小二乘轉(zhuǎn)化為遞推

2、估計(jì)。最小二乘遞推算法RLS的 基本思想是一人- 一人一新的估計(jì)值(Q=老的估計(jì)值(kD+修正項(xiàng)算法迭代方程為0(k) =0(k -1) + K (k) z (k) - hr (k )0 (k -1) K(k) = P(k - 1)h(k)hr(k)P(k - 1)h(k) +1-1(1)P(k) = P(k -1) - K (k )hr (k) P(k -1)增廣最小二乘增廣最小二乘遞推算法，擴(kuò)充了最小二乘法的參數(shù)向量0和數(shù)據(jù)向量h(k)的 維數(shù)，把噪聲模型的辨識(shí)同時(shí)考慮進(jìn)去，因此被稱為增廣最小二乘法。最小二乘 法的許多結(jié)論對(duì)它都是適用的，但最小二乘法只能獲得模型的參數(shù)估計(jì)。如果噪 聲模型必

3、須用D(z-1)v(k)表示時(shí)，只能用RELS算法，方可獲得無偏估計(jì)，這是 RLS算法所不能代替的?？紤]SISO的動(dòng)態(tài)系統(tǒng)，輸入u(k)和輸出z(k)是可以觀測(cè)的；G(z-1)是系統(tǒng)模型，用來描述系統(tǒng)的輸入輸出特性，y模型，用來描述系統(tǒng)的輸入輸出特性，y(k)是系統(tǒng)的實(shí)際輸出。N(z-1)是噪聲模型，v(k)是均值為零的不相關(guān)隨機(jī)噪聲。通常G(z-i) = B1), N(z-i) = D(Z一 1)A(z-i)C(z -i)式中A(z-i) = i + a z-i + a z-2 hf a z-na TOC o 1-5 h z 12氣B(z-i) = b z-i + b z-2 + b z-n

4、bi2nbC(z -i) = i + c z -i + c z -2 + c z-nc12氣D(z-i) = i + d z-i + d z-2 + d z-ndi2nd 若SISO系統(tǒng)采用平均滑動(dòng)模型，即C(z-i) = A(z-i)A(z-i)z(k) = B(z-i)u (k) + D( z -i)v (k)若假定模型階次氣、nb、nd已經(jīng)確定，則這類問題的辨識(shí)可用增廣最小二 乘法，以便獲得滿意的結(jié)果。令0 = a ,a ，a ,b ,b，,b ,d ,d ，d T i 2 匕 i 2nbi 2nd(5)I h(k) = -z(k 一 i),，-z(k 一 n ), u(k 一 i),u

5、(k 一 n ), v(k 一 i),v(k 一 n )tabd將模型(3)化為最小二乘格式z(k) = h t (k )0 + v(k)由于v(k)是白噪聲，所以利用最小二乘法即可獲得參數(shù)0的無偏估計(jì)。但是數(shù)據(jù)向量h(k)中包含著不可測(cè)的噪聲量v(k -i),v(k -七),它可用相應(yīng)的估計(jì)值替 代。置 TOC o 1-5 h z h(k) = -z(k -1),，一z(k -n ),u(k -1),u(k -n ), v(k -1), v(k -n )t (7) abd式中，v(k) = 0,k 0 時(shí)，v(k) = z(k) 一 hT (k)0 (k -1)(8)或v(k) = z(k)

6、 一 hT (k)0 (k)(9)增廣最小二乘遞推算法RELSB(k) =9(k -1) + K (k) z (k) - hr (k )9 (k -1)K(k) = P(k - 1)h(k)hr(k)P(k - 1)h(k) +1/A(k)-i(10)P (k) = P(k -1) - K (k) K (k )t hr (k) P( k - 1)h(k) +1 如果1 = 1,即所有采樣數(shù)據(jù)都是等同加權(quán)時(shí)，增廣最小二乘遞推算法RELS可以 A寫為9 (k) =9(k -1) + K (k) z (k) - hr (k )9 (k -1)K(k) = P(k - 1)h(k)hr(k)P(k -

7、 1)h(k) +1-1(11)P(k) = I - K (k )hr (k )P(k -1)三、實(shí)驗(yàn)過程產(chǎn)生白噪聲和m序列白噪聲序列e(t)可由乘同余法自己編寫函數(shù)實(shí)現(xiàn)，也可以由matlab自帶的 randn函數(shù)產(chǎn)生，本實(shí)驗(yàn)采用randn函數(shù)產(chǎn)生；系統(tǒng)的輸入u(t)選擇m序列，其 級(jí)數(shù)為4，初相為1 1 1 0,本原多項(xiàng)式為f(x) = x4 + x3 +1。產(chǎn)生輸出采樣信號(hào)考慮噪聲對(duì)輸出的影響，則系統(tǒng)輸出為z (k)= 2.851z (k-1)-2.717 z (k - 2 )+ 0.865 z (k - 3)+ u (k -1)+u (k - 2)+ u (k - 3)+ v (k)+

8、0.7v (k -1)+ 0.22v (k - 2)遞推最小二乘辨識(shí)給被辨識(shí)參數(shù)9和P賦初值。根據(jù)(1)式計(jì)算K(k)、9 (k)、P(k)，當(dāng)參數(shù) 收斂滿足要求時(shí)停止迭代計(jì)算，輸出辨識(shí)結(jié)果。增廣最小二乘辨識(shí)給被辨識(shí)參數(shù)9和P賦初值。根據(jù)(11)式計(jì)算K(k)、9 (k)、P(k)，當(dāng)參數(shù) 收斂滿足要求時(shí)停止迭代計(jì)算，輸出辨識(shí)結(jié)果。實(shí)驗(yàn)結(jié)果與分析系統(tǒng)白噪聲和系統(tǒng)輸入山序列如圖1所示。圖1系統(tǒng)噪聲及m序列噪聲序列圖1系統(tǒng)噪聲及m序列噪聲序列F 呷：pp件卜輸入M序列J甲 ：甲孑 PB-2010203040 SO EO 708090100米用遞推最小二乘辨識(shí)算法實(shí)驗(yàn)時(shí)，給定辨識(shí)精度為E=, m序列

9、信號(hào)峰峰大小為10,經(jīng)過157次迭代， 達(dá)到給定經(jīng)度。此時(shí)辨識(shí)結(jié)果為a =, a=, a =b, b2=, b =在=實(shí)驗(yàn)時(shí)，對(duì)于m序列，其峰峰值越大，辨識(shí)的參數(shù)與真值誤差越小。遞推最小二乘的仿真結(jié)果如圖2和圖3所示。圖2遞推最小二乘參數(shù)結(jié)果遂推是小二推法識(shí)黠杲-500 -100D-1500-2Q0D L圖2遞推最小二乘參數(shù)結(jié)果遂推是小二推法識(shí)黠杲-500 -100D-1500-2Q0D L氛敏遮差收皴情況圖3參數(shù)誤差收斂情況增廣最小二乘辨識(shí)實(shí)驗(yàn)時(shí)，給定辨識(shí)精度為E=，m序列峰峰值為2,經(jīng)過40次迭代達(dá)到給定經(jīng) 度，此時(shí)辨識(shí)結(jié)果為a =, a =, a =；1 2 3，b 1，b2=1，b =

10、1； c 1，c2=，c =增廣最小二乘的仿真結(jié)果如圖4和圖5所示。、二乘圖5參數(shù)誤差收斂情況、二乘參數(shù)的真值和估計(jì)值如下表所示：參數(shù)真值遞推最小二乘辨識(shí)增廣最小二乘辨識(shí)aaabbb3-cc從上面兩表對(duì)比可以看到，增廣最小二乘考慮了噪聲模型，與遞推最小 算法相比，速度快、辨識(shí)結(jié)果精確，而且可以得到噪聲模型參數(shù)。四.附錄附錄1遞推最小二乘辨識(shí)matlab仿真程序clearL=200;y1=1;y2=1;y3=1;y4=0;for i=1:Lx1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;if y(i)u(i)=-5;elseu(i)=5;endy1=x1;y2=x

11、2;y3=x3;y4=x4;endfigure(1)i=1:20;stem(u(i);grid on;title(,EaEeMDdADAC 2035OUU,)v=1*randn(L,1);z(3)=0;z(2)=0;z(1)=0;for k=4:L;z(k)=*z(kT)*z(k-2)+*z(k-3)+u(kT)+u(k-2)+u(k-3)+v(k)+*v(kT)+*v(k-2);endc0=;p0=103*eye(6,6);E=;c=c0,zeros(6,LT);e=zeros(6,L);for k=4:Lh1=-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(

12、k-3);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;c(:,k)=c1;p1(:,:,k)=p0-k1*k1*h1*p0*h1+1;p0=p1(:,:,k);if abs(e2)u(i)=-1;elseu(i)=1;endy1=x1;y2=x2;y3=x3;y4=x4;endv=1*randn(1,L);z(3)=0;z(2)=0;z(1)=0;for k=4:L;z(k)=*z(kT)*z(k-2)+*z(k-3)+u(kT)+u(k-2)+u(k-3)+v(k)+*v(

13、kT)+*v(k-2);endfigure(1);subplot(2,1,2);stem(u),grid on ;title(，EaEeMBoAD，);subplot(2,1,1);plot(v),grid on ;title(，OeEuDoAB，);c0=，;p0=106*eye(9,9);E=5e-10;c=c0,zeros(9,LT);e=zeros(9,L);N=L+1;for k=4:L;h1=-z(kT),-z(k-2),z(k-3),u(kT),u(k-2),u(k-3),v(k),v(kT),v(k-2);k1=p0*h1*inv(h1*p0*h1+1);c1=c0+k1*(z

14、(k)-h1*c0);e1=c1-c0;e2=e1./c0;e(:,k)=e2;c0=c1;c(:,k)=c1;p1(:,:,k)=p0-k1*k1*h1*p0*h1+1;p0=p1(:,:,k);if abs(e2)=EN=k;break;endendif N=L+1a1=c(1,1:N); a2=c(2,1:N); a3=c(3,1:N);b1=c(4,1:N); b2=c(5,1:N);b3=c(6,1:N);d1=c(7,1:N); d2=c(8,1:N); d3=c(9,1:N);ea1=e(1,1:N); ea2=e(2,1:N); ea3=e(3,1:N);eb1=e(4,1:N); eb2=e(5,1:N);eb3=e(6,1:N);ed1=e(7,1:N); ed2=e(8,1:N); ed3=e(9,1:N);figure(2);i=1:N;plot(i,a1,r,i,a2,r:,i,a3,r+,i,b1,b,i,b2,，b:，,i