系統(tǒng)辨識(shí)與參數(shù)估計(jì)大作業(yè)教案

1、系統(tǒng)辨識(shí)與參數(shù)估計(jì)大作業(yè)第一題遞推最小二乘估計(jì)參數(shù)考慮如上圖所示的仿真對(duì)象，選擇模型結(jié)構(gòu)為： z(k) a1z(k -1) a2z(k -2) =bu(k -1) b2u(k -2) v(k),其中v(k)是服從N(0,1)正態(tài)分布的不相關(guān)隨機(jī)噪聲； 輸入信號(hào)u(k)采用4階逆M序列, 特征多項(xiàng)式取F(s) =13s$ s4,幅度為1,循環(huán)周期為Np =62bit ;控制九值，使數(shù)據(jù) 的噪信比分別為 10%, 73%, 100%三種情況。加權(quán)因子 A(k)=1；數(shù)據(jù)長度L=500;初始 條件取 P (0) =1061 , ?(0) =0.001 ,利用遞推最小二乘算法在線估計(jì)參數(shù)，利用模型階次

2、辨識(shí)方法(AIC準(zhǔn)則)，確定模型的階次。估計(jì)噪聲v(k)的方差和模型靜態(tài)增益 K作出參數(shù)估計(jì)值隨時(shí)間的變化圖答：設(shè)過程的輸入輸出關(guān)系可以描述成z(k) = hT (k)H + n(k)z(k)是輸出量，h(k)是可觀測(cè)的數(shù)據(jù)向量，n(k)是均值為0的隨機(jī)噪聲,Jh(k) - I -z(k -1), -z(k -2), u(k -1),u(k -2)3 - 1al ,a2,b1, b2 T選取的模型為結(jié)構(gòu)是z(k) = -a1Z(k 一1) 一a2z(k 2) bu1 (k -1) bu2(k -2) ai 二1.5,a2 =0.7,“ =1.0,b2 =0.5加權(quán)最小二乘參數(shù)估計(jì)遞推算法RWL

3、S的公式如下，1K(k) = p(k -1)h(k) hT(k)p(k-1)h(k) TOC o 1-5 h z IL(k)9(k) =0(k 1) + K(k) z(k) hT(k聲(k 1)p(k)-|I -K(k)hT(k) p(k -1)為了把p(k)的對(duì)稱性，可以把 p(k)寫成T T1p(k) = p(k-1)-K(k)KT(k) hT(k)p(k -1)h(k)IL”(k)如果把A(k)設(shè)成1的時(shí)候，加權(quán)最小二乘法就退化成最小二乘法。用AIC準(zhǔn)則定階法來定階，所用公式乙VnZn = lz(1),z(2), z(3),., z(L)T-jT*n =回，-a2,-ana回以。z(0)

4、z(-1).z(1 - n)u(0)u(-1).u(1-n)z(1)z(0).z(2 -n)u(1)u(0).u(2 - n)Hn= ._z(L-1) z(L-2) . z(L-n) u(L-1) u(L-2) . u(L-n)_2其中模型參數(shù) 我和 噪聲V(k)方差的極大似然估計(jì)值為19 ML ，仃v=(HnTHn)HnTZn21T(Z = (7 - H ) (Z - H )v l (乙n n n ml )(乙 n n ml )AIC的定階公式寫成2AIC (n) = L log 4n取n=1,2,3,4;分別計(jì)算AIC(n),找到使AIC(n)最小的那個(gè)n作為模型的階次。一般說來,這樣得到

5、的模型階次都能比較接近實(shí)際過程的真實(shí)階次。 信噪比為10%時(shí)：060100150 200250300350400 450500參數(shù)ala2b1b2噪聲力差靜態(tài)增益模型階次真值-1.50.710.51估計(jì)值-1.5190.722591.03140.509231.09517.56612信噪比為 參數(shù)73% 時(shí)：a1a2b1b2噪聲力差靜態(tài)增益模型階次真值-1.50.710.51傳計(jì)值-1.5190.722591.03140.509231.09517.566124111fli111nl050100150200250300350400450500信噪比為 參數(shù)100% 時(shí)：ala2b1b2噪聲力差靜態(tài)

6、增益模型階次真值-1.50.710.51傳計(jì)值-1.5190.722591.03140.50923情計(jì)值7.56612源程序:%function a1 a2 bl b2 na nb fangcha Kk=rwls(L,syn,Np) %na,nb 為模型階次，fangcha 為噪聲 方差，Kk為靜態(tài)增益a1=0;a2=0;b1=0;b2=0;na=0;nb=0;fangcha=0;Kk=0;L=500;Np=62;syn=1;x(1:4)=1 0 1 0;for i=1:Nptemp=xor(x(1),x(4);M(i)=x(4);for j=4:-1:2x(j)=x(j-1);endx(1)

7、=temp;endS=ones(1,Np);%先產(chǎn)生一個(gè)全是 1的序列Sif mod(Np,2)=0% 判斷Np是奇數(shù)還是偶數(shù)p=Np/2;else p=(Np-1)/2;endfor j=1:pS(2*j)=0;%將S序列的偶數(shù)位值均置為0,從而使S序列是0或1的方波序列endIM=xor(M,S); %使用M序列與方波序列 S復(fù)合生成逆 M序列IMu=IM*2-1;for i=(Np+1):Lu(i)=u(i-Np);endrandn(seed,2);v=randn(1,L);syms c;y(1)=0;y(2)=u(1);e(1)=c*v(1);e(2)=1.5*e(1)+c*v(2);

8、for i=3:Ly(i)=1.5*y(i-1)-0.7*y(i-2)+u(i-1)+0.5*u(i-2);e(i)=1.5*e(i-1)-0.7*e(i-2)+c*v(i);endm=sum(e.A2);n=sum(y.A2);%c=solve(m/n-syn*syn,c);c=solve(m/n-syn*syn);c1=abs(double(c(1);z(1)=c1*v(1);z(2)=u(1)+1.5*z(1)+c1*v(2);for i=3:Lz(i)=1.5*z(i-1)-0.7*z(i-2)+u(i-1)+0.5*u(i-2)+c1*v(i);endcta=zeros(4,L);c

9、ta(:,1)=0.001 0.001 0.001 0.001;P=diag(10A6 10A6 10A6 10A6);%k=2 時(shí)h=-z(1) 0 u(1) 0;K= P*h*inv(h*P*h+1);P=P-K*K*(h*P*h+1);cta(:,2)=cta(:,1)+K*(z(2)-h*cta(:,1);for k=3:Lh=-z(k-1) -z(k-2) u(k-1) u(k-2);K=P*h*inv(h*P*h+1);P=P-K*K*(h*P*h+1);cta(:,k)=cta(:,(k-1)+K*(z(k)-h*cta(:,(k-1);end %以上為參數(shù)估計(jì)值z(mì)=z;for

10、na=1:4for nb=1:4A=zeros(L,na);B=zeros(L,nb);for i=1:Lfor j=1:naif ijA(i,j)=z(i-j);endendendfor i=1:Lfor j=1:nbif ijB(i,j)=u(i-j);endendendH=A B;cta1=inv(H*H)*H*z;%cta1為模型參數(shù)極大似然估計(jì)值cgm(na,nb)=(z-H*cta1)*(z-H*cta1)/L;%cgm為噪聲方差極大似然估計(jì)值A(chǔ)IC(na,nb)=L*log(cgm(na,nb)+2*(na+nb);endendna nb=find(AIC=min(min(AIC

11、);fangcha=cgm(na,nb)/(c1A2);a1=cta(1,500);a2=cta(2,500);b1=cta(3,500);b2=cta(4,500);Kk=(cta(3,L)+cta(4,L)/(1+cta(1,L)+cta(2,L);m=1:L;plot(m,cta(1,:),b-,m,cta(2,:),k-,m,cta(3,:),y-,m,cta(4,:),r-) 第二大題卡爾曼濾波一個(gè)系統(tǒng)模型為x1(k 1) = x1(k) x2(k) w(k), k = 0,1 x2 (k 1) = x2 (k) w(k)同時(shí)有下列條件:初始條件已知且有。x(0)=0,0Tw(k)是

12、一個(gè)標(biāo)量零均值白高斯序列，且自相關(guān)函數(shù)已知為 Ew(j)w(k) =3jk ,乙(k 1) = x1(k 1) v1(k 1), k =0,1另外，我們有下列觀測(cè)模型，即,，且z2 (k 1) = x2(k 1) v2(k 1)有下列條件：v1(k+1)和v2(k+1)是獨(dú)立的零均值白高斯序列，且有Ev1(jW(k)=染，Ev2(j)v2(k)=25jk, k = 0,1,2.對(duì)于所有的j和k , w(k)與觀測(cè)噪聲過程 *+1)和丫2%+1)是不相關(guān)的，即Ew(j)v1(k) =0, Ew(jM(k)=0我們希望得到由觀測(cè)矢量z(k+1)=乙(k+1),z2(k+1)T估計(jì)狀態(tài)矢量x(k +

13、1) =x1(k+1),x2(k +1)T的卡爾曼濾波的公式表示，并求解以下問題：求出卡爾曼增益矩陣，并得出最優(yōu)估計(jì)x(k +1)和觀測(cè)矢量z(1), z(2),., z(k +1)之間的遞推關(guān)系。用模擬數(shù)據(jù)確定狀態(tài)矢量x(k)的估計(jì)值 父(k|k), k =0,1,2.10 ,并畫出當(dāng)k=0,1,2.10 時(shí)，R(k|k), x2(k|k)的圖通常，狀態(tài)矢量的真實(shí)值是得不到的，但是為了用作圖來說明問題，表 1和表2給出了狀態(tài)矢量元素的真實(shí)值。對(duì)于k = 0,1,2.10 ,在同一幅圖中畫出真實(shí)值和在(b)中確定的x1(k)的估計(jì)值。對(duì)x2(k)重復(fù)這一過程。當(dāng) k從1變到10時(shí)，對(duì)每個(gè)元素i

14、 =1,2 ,計(jì)算并畫出各自的誤差圖，即xi(k)-?i(k|k)o當(dāng)k從i變到10時(shí)，通過用由卡爾曼濾波器決定的狀態(tài)誤差協(xié)方差矩陣畫出E 2 (k | k)和 Ezl (k | k),而備(k | k)=為(k) ? (k | k),i = 1,2討論一下(C)中計(jì)算的誤差與(d)中方差之間的關(guān)系。表1觀測(cè)值時(shí)間卜標(biāo)k觀測(cè)值z(mì)1 (k)觀測(cè)值Z2(k)13. 296919692. 1013429423. 387365150. 4754079737. 028306413. 1768889849. 712125212. 49811140511, 420183152. 91992424615. 9

15、79805836. 069342855. 42519274828, 302127813. 05365741930. 446838315. 980511411038. 758755954. 51016361表2由模擬得到的實(shí)際狀態(tài)值時(shí)間卜標(biāo)k實(shí)際狀態(tài)值x1(k)實(shí)際狀態(tài)值x2(k)00.000000000.0000000011.654287141.6542871423.503007021.8487198835.997852922.4755222249.150407403508739103.35833170616.921925944.413186

16、84721.344833524.42290758825.893351444.54851792931.541353305.648001861036.936056705.39447034答:(a)卡爾曼增益矩陣：Hk=Pk*CTM(CMPkMCT+R)，估計(jì)值與觀測(cè)值之間的遞歸關(guān)系為：Xk=Ak Xk-i Hk (Yk-Ck Ak Xk-i)(b)狀態(tài)矢量估計(jì)值的計(jì)算框圖：(c) xi (k/k)和 x2 (k/k)的圖:(d)真實(shí)值與估計(jì)值的比較圖:(e)通過用卡爾曼濾波器的狀態(tài)誤差協(xié)方差矩陣畫出的E與2(k/k)和E42(k/k):狀態(tài)矢量X每個(gè)單獨(dú)分量估計(jì)的方差(f)分析：(e)中的方差是(

17、d)中的誤差平方后取均值，是均方誤差。誤差直 接由真實(shí)值減去估計(jì)值，有正有負(fù)，而均方誤差沒有這個(gè)缺陷，更能綜合的表示 濾波的效果。源程序：%P K X Z均為2*2矩陣，用三維矩陣表示他們，記錄每一個(gè)k值下的情況，如 P(:,:,k)%k從1開始,k=1時(shí)取初值clear all;P(:,:,1)=10 0;0 10;%P(0|0)初值，P 表示 P(k|k)X(:,:,1)=0;0 ; %X(0)初值fy= 1 1;0 1; %()tao=1;1;% rH= 1 0;0 1;Q=1;R=1 0;0 2;Z(:,:,2)=3.29691969;2.10134294;Z(:,:,3)=3.387

18、36515;0.47540797;Z(:,:,4)=7.02830641;3.17688898;Z(:,:,5)=9.71212521;2.49811140;Z(:,:,6)=11.42018315;2.91992424;Z(:,:,7)=15.97870583;6.17307616;Z(:,:,8)=22.06934285;5.42519274;Z(:,:,9)=28.30212781;3.05365741;Z(:,:,10)=30.44683831;5.98051141;Z(:,:,11)=38.75875595;4.51016361;for k=1:10P1(:,:,k尸fy*P(:,:

19、,k)*fy+tao*Q*tao; %P1(:,:,k)表示 P(k+1|k)K(:,:,k+1)=P1(:,:,k)*H*inv(H*P1(:,:,k)*H+R);X(:,:,k+1)=fy*X(:,:,k)+K(:,:,k+1)*(Z(:,:,k+1)-H*fy*X(:,:,k);P(:,:,k+1)=(eye(2,2)-K(:,:,k+1)*H)*P1(:,:,k);endgx1(1:11)=X(1,1,1:11); %gx1,gx2 分別表示 X 的估計(jì)值gx2(1:11)=X(2,1,1:11);plot(0:10,gx1,0:10,gx2);legend(x1 估計(jì)值,x2 估計(jì)值,location,northwest);title(狀態(tài)矢量 X的估計(jì)值x1(k)和x2(k);xlabel(k);ylabel(X);x1=0,1.65428714,3.50300702,5.997852924,9.15040740,12.50873910,.16.92192594,21.344