1、 系統(tǒng)辨識(shí)課程作業(yè) 題目1已知系統(tǒng)的差分方程為：其中 是均值為0， ,并服從正態(tài)分布的不相關(guān)隨機(jī)噪聲， 采用4級(jí)移位寄存器產(chǎn)生的幅度為1的M序列， 。數(shù)據(jù)長(zhǎng)度取N=240，請(qǐng)繪出Newton-Raphson方法求參數(shù)的極大似然估計(jì)程序流程圖并附上Matlab程序。 1.流程圖圖1.1 流程框圖2.程序見附錄1。3.運(yùn)行結(jié)果參見表1.1和圖1.2、圖1.3。4.結(jié)果：收斂較快；參數(shù)a1,a2,b1,b2（圖1.3中藍(lán)色曲線）的估計(jì)準(zhǔn)確，由于噪聲的隨機(jī)性，d1,d2（圖2.1中紅色曲線）的估計(jì)誤差較大；系統(tǒng)輸出誤差隨噪聲變量均方根增大而增大，即使當(dāng)噪聲水平較高時(shí)也可以獲得較好的估計(jì)；但是必須進(jìn)行N

2、次觀測(cè)，才可以遞推一次，不適用于實(shí)時(shí)應(yīng)用。表1.1 模型參數(shù)辨識(shí)結(jié)果被估參數(shù)a1a2b1b2d1d2真實(shí)值-1.50.71.00.5-10.2估計(jì)值-1.48840.69270.97460.6057-0.94690.1865圖1.2 參數(shù)估計(jì)收斂值圖1.3 參數(shù)估計(jì)誤差題目2動(dòng)態(tài)系統(tǒng)模型為 其中 是均值為零，方差為1并服從正態(tài)分布的不相關(guān)隨機(jī)噪聲； 是由四位移位寄存器產(chǎn)生的幅度為1的M序列,利用遞推極大似然算法對(duì)系統(tǒng)參數(shù)進(jìn)行辨識(shí)，遞推停止條件為N=1500。(要求仿真程序附注釋)。1程序見附錄2。2運(yùn)行結(jié)果見表2.1和圖2.1、圖2.2。3結(jié)果分析：收斂較快；參數(shù)a1,a2,b1,b2（圖2.

3、1中藍(lán)色曲線）的估計(jì)準(zhǔn)確，由于噪聲的隨機(jī)性，d1,d2（圖2.1中紅色曲線）的估計(jì)誤差較大；系統(tǒng)輸出誤差隨噪聲變量均方根增大而增大，即其決定估計(jì)的準(zhǔn)確程度。表2.1 動(dòng)態(tài)模型參數(shù)辨識(shí)結(jié)果被估參數(shù)a1a2b1b2d1d2真實(shí)值0.5-0.21.20.3-10.8估計(jì)值0.5064-0.18521.18270.3301-0.98290.7742圖2.1 參數(shù)估計(jì)誤差圖2.2 系統(tǒng)輸出誤差附錄題目1程序clear all;close all;clc;sigma=7.2;%均方差epsilon=0.001;%迭代終止條件total=16;N=240;V=sigma*randn(total,1); %噪

4、聲%M序列產(chǎn)生y1=1;y2=1;y3=1;y4=0;for i=1:15 x1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4; if y(i)0.5,u(i)=-1; else u(i)=1; end y1=x1;y2=x2;y3=x3;y4=x4;end%最小二乘一般算法,產(chǎn)生初始估計(jì)值a1,a2,b1,b2;z=zeros(1,total);for k=3:total z(k)=1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2)+V(k);end%給樣本系數(shù)矩陣H=-z(2) -z(1) u(2) u(1);-z(3) -z

5、(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) 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

6、);Z=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);c=(H*H)(H*Z);a1=c(1);a2=c(2);b1=c(3);b2=c(4);d1=0.1;d2=0.1;theta=a1,a2,b1,b2,d1,d2;%參數(shù)估計(jì)初值v(1)=0;v(2)=0;d_theta1=zeros(6,1);d_theta2=zeros(6,1);v_da1(1)=0;v_da2(1)=0;v_db1(1)=0;v_db2(1)=0;v_dd1(1)=0;v_dd2(1)=0;v_da1(2)=

7、0;v_da2(2)=0;v_db1(2)=0;v_db2(2)=0;v_dd1(2)=0;v_dd2(2)=0;j=1;bef=zeros(6,1);%執(zhí)行算法while sum(abs(theta-bef)epsilon %采集輸入輸出 for i=1:N+3 x1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4; if y(i)0.5,u(i)=-1; else u(i)=1; end y1=x1;y2=x2;y3=x3;y4=x4; endy(1)=0;y(2)=0;V=sigma*randn(N+3,1); %噪聲y(1)=1;y(2)=0.01;

8、for k=3:N+3 y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-1)+0.5*u(k-2)+V(k)-V(k-1)+0.2*V(k-2); endJ_d=0;JJ_d=0;a1=theta(1);a2=theta(2);b1=theta(3);b2=theta(4);d1=theta(5);d2=theta(6); for k=3:N+3 v(k)=y(k)+a1*y(k-1)+a2*y(k-2)-b1*u(k-1)-b2*u(k-2)-d1*v(k-1)-d2*v(k-2);%求取v(k) v_da1(k)=y(k-1)-d1*v_da1(k-1)-d2*v_da1(k

9、-2); v_da2(k)=y(k-2)-d1*v_da2(k-1)-d2*v_da2(k-2); v_db1(k)=-u(k-1)-d1*v_db1(k-1)-d2*v_db1(k-2); v_db2(k)=-u(k-2)-d1*v_db2(k-1)-d2*v_db2(k-2); v_dd1(k)=-v(k-1)-d1*v_dd1(k-1)-d2*v_dd1(k-2); v_dd2(k)=-v(k-2)-d1*v_dd2(k-1)-d2*v_dd2(k-2); d_theta=v_da1(k),v_da2(k),v_db1(k),v_db2(k),v_dd1(k),v_dd2(k); J_d

10、=J_d+v(k)*d_theta; JJ_d=JJ_d+d_theta*d_theta; endbef=theta;theta=theta-(JJ_d)J_d;v(1)=v(N+1);v(2)=v(N+2);v_da1(1)=v_da1(N+1);v_da2(1)=v_da2(N+1);v_db1(1)=v_db1(N+1);v_db2(1)=v_db2(N+1);v_dd1(1)=v_dd1(N+1);v_dd2(1)=v_dd2(N+1);v_da1(2)=v_da1(N+2);v_da2(2)=v_da2(N+2);v_db1(2)=v_db1(N+2);v_db2(2)=v_db2(

11、N+2);v_dd1(2)=v_dd1(N+2);v_dd2(2)=v_dd2(N+2);%求取誤差error1(j)=-1.5-theta(1);error2(j)=0.7-theta(2);error3(j)=1-theta(3);error4(j)=0.5-theta(4);error5(j)=-1-theta(5);error6(j)=0.2-theta(6);v_error(j)=v(N+2);j=j+1;endtheta%輸出估計(jì)參數(shù)%作圖figure(1);plot(error1)hold onplot(error2)hold onplot(error3)hold onplot(

12、error4)hold onplot(error5,r)hold onplot(error6,r)title(參數(shù)估計(jì)誤差)xlabel(迭代次數(shù))ylabel(誤差)hold offfigure(2);plot(-1.5-error1,b);hold onplot(0.7-error2,c)hold onplot(1-error3,g)hold onplot(0.5-error4,y)hold onplot(-1-error5,m)hold onplot(0.2-error6,r)legend(a1,a2,b1,b2,d1,d2,-1)title(參數(shù)估計(jì)值變化)xlabel(迭代次數(shù))yl

13、abel(參數(shù))hold off題目2程序clear all;close all;clc;total=1500;sigma=1.0;%M序列輸入A1=1;A2=1;A3=1;A4=0;for i=1:1:total x1=xor(A3,A4); x2=A1; x3=A2; x4=A3; OUT(i)=A4; if OUT(i)0.5 u(i)=1; else u(i)=-1; end A1=x1;A2=x2;A3=x3;A4=x4;endfigure(1);stem(u,r)grid on;%產(chǎn)生輸出y(1)=0;y(2)=0;v=sigma*randn(total,1); %噪聲y(1)=1

14、;y(2)=0.01;for k=3:total y(k)=-0.5*y(k-1)+0.2*y(k-2)+1.2*u(k-1)+0.3*u(k-2)+v(k)-v(k-1)+0.8*v(k-2);end%初始化theta0=0.0001*ones(6,1);error1(1)=0.5-theta0(1);error2(1)=-0.2-theta0(2);error3(1)=1.2-theta0(3);error4(1)=0.3-theta0(4);error5(1)=-1-theta0(5);error6(1)=0.8-theta0(6);a1=theta0(1);a2=theta0(2);b

15、1=theta0(3);b2=theta0(4);d1=theta0(5);d2=theta0(6);P0=eye(6,6);fai0=0.1*-1 -1 1 1 1 1;e(1)=1;e(2)=1;%遞推計(jì)算for i=3:total kusai=-y(i-1) -y(i-2) u(i-1) u(i-2) e(i-1) e(i-2); Q=-y(i-1) 0 u(i-1) 0 e(i-1) 0; A=zeros(6,6); A(1,1)=-d1; A(3,3)=-d1; A(5,5)=-d1; A(1,2)=-d2; A(3,4)=-d2; A(5,6)=-d2; A(2,1)=1; A(4

16、,3)=1; A(6,5)=1; fai=A*fai0+Q; K=P0*fai/(fai*P0*fai+1); P=eye(6,6)-K*fai*P0; e(i)=y(i)-kusai*theta0; theta=theta0+K*e(i); P0=P; theta0=theta; fai0=fai; a1=theta(1);a2=theta(2); b1=theta(3);b2=theta(4); d1=theta(5);d2=theta(6); error1(i)=0.5-theta(1); error2(i)=-0.2-theta(2); error3(i)=1.2-theta(3); error4(i)=0.3-theta(4); error5(i)=-1-theta(5); error6(i)=0