1、北航系统辨识课程作业系统辨识课程作业一、 题目1已知系统的差分方程为：其中 是均值为0， ,并服从正态分布的不相关随机噪声， 采用4级移位寄存器产生的幅度为1的M序列， 。数据长度取N=240，请绘出Newton-Raphson方法求参数的极大似然估计程序流程图并附上Matlab程序。 (一) 流程图图1.1 流程框图(二) 程序见附录1。(三) 运行结果参见表1.1和图1.2、图1.3。(四) 结果分析：收敛较快；参数a1,a2,b1,b2（图1.3中蓝色曲线）的估计准确，由于噪声的随机性，d1,d2（图2.1中红色曲线）的估计误差较大；系统输出误差随噪声变量均方根增大而增大，即使当噪声水平
2、较高时也可以获得较好的估计；但是必须进行N次观测，才可以递推一次，不适用于实时应用。表1.1 模型参数辨识结果被估参数a1a2b1b2d1d2真实值-1.50.71.00.5-10.2估计值-1.46800.67620.97500.6143-0.94600.1688图1.2 参数估计收敛值图1.3 参数估计误差二、 题目2动态系统模型为其中 是均值为零，方差为1并服从正态分布的不相关随机噪声； 是由四位移位寄存器产生的幅度为1的M序列,利用递推极大似然算法对系统参数进行辨识，递推停止条件为N=1500。(要求仿真程序附注释)。(一) 程序见附录2。(二) 运行结果见表2.1和图2.1、图2.2
3、。(三) 结果分析：收敛较快；参数a1,a2,b1,b2（图2.1中蓝色曲线）的估计准确，由于噪声的随机性，d1,d2（图2.1中红色曲线）的估计误差较大；系统输出误差随噪声变量均方根增大而增大，即其决定估计的准确程度。表2.1 动态模型参数辨识结果被估参数a1a2b1b2d1d2真实值0.5-0.21.20.3-10.8估计值0.4806-0.19331.19100.2862-0.93490.7357图2.1 参数估计误差图2.2 系统输出误差三、 附录1. 题目1程序clcclose allsigma=7.2;%均方差epsilon=0.001;%迭代终止条件total=16;N=240;
4、V=sigma*randn(total,1); %噪声%M序列产生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%最小二乘一般算法,产生初始估计值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%给样本系数矩阵H=-z(2
5、) -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(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)
6、; -z(15) -z(14) u(15) u(14);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=inv(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;%参数估计初值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_d
7、d1(1)=0;v_dd2(1)=0;v_da1(2)=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);%执行算法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
8、+3,1); %噪声y(1)=1;y(2)=0.01; 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
9、-1)-d1*v_da1(k-1)-d2*v_da1(k-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_d
10、b2(k),v_dd1(k),v_dd2(k); J_d=J_d+v(k)*d_theta; JJ_d=JJ_d+d_theta*d_theta; endbef=theta;theta=theta-inv(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_d
11、b1(2)=v_db1(N+2);v_db2(2)=v_db2(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%输出估计参数%作图figure(1);plot(error1)hold onplot(error2
12、)hold onplot(error3)hold onplot(error4)hold onplot(error5,r)hold onplot(error6,r)title(参数估计误差)xlabel(迭代次数)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
13、,-1)title(参数估计值变化)xlabel(迭代次数)ylabel(参数)hold off2. 题目2程序clcclose alltotal=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%产生输出y(1)=0;y(2)=0;v=sigma*rand
14、n(total,1); %噪声y(1)=1;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=the
15、ta0(1);a2=theta0(2);b1=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;%递推计算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,
16、6)=-d2; A(2,1)=1; A(4,3)=1; A(6,5)=1; fai=A*fai0+Q; K=P0*fai*inv(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);
17、error3(i)=1.2-theta(3); error4(i)=0.3-theta(4); error5(i)=-1-theta(5); error6(i)=0.8-theta(6);endfigure(2);plot(error1)hold onplot(error2)hold onplot(error3)hold onplot(error4)hold onplot(error5,r)hold onfigure(2);plot(error6,r)title(参数估计误差)xlabel(迭代次数)ylabel(误差变化)hold offfigure(3);plot(e)title(输出误差)xlabel(迭代次数)ylabel(误差变化)%输出辨识参数theta
