北航《系统辨识》课程作业

姓名：xx学号：SYxx2013/11/301《系统辨识》课程作业一、题目1已知系统的差分方程为：其中是均值为0，,并服从正态分布的不相关随机噪声，采用4级移位寄存器产生的幅度为1的M序列，1ms。数据长度取N=240，请绘出Newton-Raphson方法求参数的极大似然估计程序流程图并附上Matlab程序。(一)流程图(1),(1),(2),...(n),(1)(n),...,,(n1),(n2),...,(n)N(j)(k)22(j)(k)()J(j)(k)J22(j)(j)(k)(k)(j1)(j)[]()JJ(j1)(k)(k)图1.1流程框图(二)程序见附录1。(三)运行结果参见表1.1和图1.2、图1.3。(四)结果分析：收敛较快；参数a1,a2,b1,b2（图1.3中蓝色曲线）的估计准确，由于噪声的随机性，d1,d2（图2.1中红色曲线）的估计误差较大；系统输出误差随噪声变量均方根增大而增大，即使当噪声水平较高时也()1.5(1)0.7(2)(1)0.5(2)()(1)0.2(2)ykykykukukkkk()k7.2()uk姓名：xx学号：SYxx2013/11/302可以获得较好的估计；但是必须进行N次观测，才可以递推一次，不适用于实时应用。表1.1模型参数辨识结果被估参数a1a2b1b2d1d2真实值-1.50.71.00.5-10.2估计值-1.46800.67620.97500.6143-0.94600.1688图1.2参数估计收敛值姓名：xx学号：SYxx2013/11/303图1.3参数估计误差二、题目2动态系统模型为其中是均值为零，方差为1并服从正态分布的不相关随机噪声；是由四位移位寄存器产生的幅度为1的M序列,利用递推极大似然算法对系统参数进行辨识，递推停止条件为N=1500。(要求仿真程序附注释)。(一)程序见附录2。(二)运行结果见表2.1和图2.1、图2.2。(三)结果分析：收敛较快；参数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()0.5(1)0.2(2)1.2(1)0.3(2)()(1)0.8(2)ykykykukukkkk()k()uk姓名：xx学号：SYxx2013/11/304图2.1参数估计误差图2.2系统输出误差姓名：xx学号：SYxx2013/11/305三、附录1.题目1程序clccloseallsigma=7.2;%均方差epsilon=0.001;%迭代终止条件total=16;N=240;V=sigma*randn(total,1);%噪声%M序列产生y1=1;y2=1;y3=1;y4=0;fori=1:15x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5,u(i)=-1;elseu(i)=1;endy1=x1;y2=x2;y3=x3;y4=x4;end%最小二乘一般算法,产生初始估计值a1,a2,b1,b2;z=zeros(1,total);fork=3:totalz(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)-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);-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]';%参数估计初值姓名：xx学号：SYxx2013/11/306v(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)=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);%执行算法whilesum(abs(theta-bef))epsilon%采集输入输出fori=1:N+3x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5,u(i)=-1;elseu(i)=1;endy1=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;fork=3:N+3y(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);fork=3:N+3v(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-2);姓名：xx学号：SYxx2013/11/307v_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=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_db1(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)holdonplot(error2)holdonplot(error3)holdonplot(error4)holdonplot(error5,'r')holdonplot(error6,'r')title('参数估计误差')姓名：xx学号：SYxx2013/11/308xlabel('迭代次数')ylabel('误差')holdofffigure(2);plot(-1.5-error1,'b');holdonplot(0.7-error2,'c')holdonplot(1-error3,'g')holdonplot(0.5-error4,'y')holdonplot(-1-error5,'m')holdonplot(0.2-error6,'r')legend('a1','a2','b1','b2','d1','d2',-1)title('参数估计值变化')xlabel('迭代次数')ylabel('参数')holdoff2.题目2程序clcclosealltotal=1500;sigma=1.0;%M序列输入A1=1;A2=1;A3=1;A4=0;fori=1:1:totalx1=xor(A3,A4);x2=A1;x3=A2;x4=A3;OUT(i)=A4;ifOUT(i)0.5u(i)=1;elseu(i)=-1;endA1=x1;A2=x2;A3=x3;A4=x4;endfigure(1)stem(u,'r')gridon姓名：xx学号：SYxx2013/11/309%产生输出y(1)=0;y(2)=0;v=sigma*randn(total,1);%噪声y(1)=1;y(2)=0.01;fork=3:totaly(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);b1=theta0(3);b2=theta0(4);d1=theta0(5);d2=theta0(6);P0=eye(6,6);fai0=0.1*[-1-11111]';e(1)=1;e(2)=1;%递推计算fori=3:totalkusai=[-y(i-1)-y(i-2)u(i-1)u(i-2)e(i-1)e(i-2)]';Q=[-y(i-1)0u(i-1)0e(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,3)=1;A(6,5)=1;fai=A*fai0+Q;K=P0*fai*inv(fai'*P0*fai