搜索
最小二乘递推算法的MATLAB仿真针对辨识模型，有z(k)-+a1*z(k-1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k)模型结构，对其进行最小二乘递推算法的MATLAB仿真，对比真值与估计值。更改a1、a2、b1、b2参数，观察结果。仿真对象：z(k)-1.5*z(k-1)+0.7*z(k-2)=u(k-1)+0.5*u(k-2)+v(k)程序如下：L=15;y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的初始值fori=1:L;%移位循环x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;%取出作为输出信号，即M序列ify(i)0.5,u(i)=-0.03;%输入信号elseu(i)=0.03;endy1=x1;y2=x2;y3=x3;y4=x4;endfigure(1);stem(u),gridonz(2)=0;z(1)=0;fork=3:15;z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%输出采样信号endc0=[0.0010.0010.0010.001]';%直接给出被识别参数的初始值p0=10^6*eye(4,4);%直接给出初始状态P0E=0.000000005;c=[c0,zeros(4,14)];e=zeros(4,15);fork=3:15;%开始求kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%开始求k的值d1=z(k)-h1'*c0;c1=c0+k1*d1;e1=c1-c0;e2=e1./c0;%求参数的相对变化e(:,k)=e2;c0=c1;c(:,k)=c1;p1=p0-k1*k1'*[h1'*p0*h1+1];%求出P(k)的值p0=p1;ife2=Ebreak;endendc,e%显示被辨识参数及其误差情况a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);figure(2);i=1:15;plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':')title('ParameterIdentificationwithRecursiveLeastSquaresMethod')figure(3);i=1:15;plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')title('IdentificationPrecision')程序运行结果：p0=1000000000010000000000100000000001000000c=Columns1through90.001000.0010-0.4984-1.2325-1.4951-1.4962-1.4991-1.49980.000100.00010.0001-0.23580.69120.69410.69900.69980.001000.25091.24971.06651.00171.00201.00020.99990.00100-0.24890.75000.56680.50200.50160.50080.5002Columns10through15-1.4999-1.5000-1.5000-1.5000-1.4999-1.49990.69990.70000.70000.70000.70000.70000.99980.99990.99990.99990.99990.99990.50020.50000.50000.50000.50000.5000e=1.0e+003*Columns1through9000-0.49940.00150.00020.00000.00000.00000000-2.3592-0.00390.00000.00000.0000000.24990.0040-0.0001-0.00010.0000-0.0000-0.000000-0.2499-0.0040-0.0002-0.0001-0.0000-0.0000-0.0000Columns10through150.00000.00000.0000-0.0000-0.00000.00000.00000.0000-0.00000.00000.00000.0000-0.00000.00000.00000.00000.0000-0.0000-0.0000-0.00000.0000-0.00000.0000-0.0000程序运行曲线：图1.输入信号图2.a1,a2,b1,b2辨识仿真结果图3.a1,a2,b1,b2各次辨识结果收敛情况分析：由运行结果可看出，输出观测值没有任何噪声成分时，辨识结果最大相对误差达到3位数。更改仿真对象参数:z(k)=z(k-1)-1.5*z(k-2)+u(k-1)+1.5*u(k-2);程序如下:L=15;y1=1;y2=1;y3=1;y4=0;fori=1:L;x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;ify(i)0.5,u(i)=-0.03;elseu(i)=0.03;endy1=x1;y2=x2;y3=x3;y4=x4;endfigure(1);stem(u),gridonz(2)=0;z(1)=0;fork=3:15;z(k)=z(k-1)-1.5*z(k-2)+u(k-1)+1.5*u(k-2);endc0=[0.001.00010.0010.001]';p0=10^6*eye(4,4)E=0.000000005;c=[c0,zeros(4,14)];e=zeros(4,15);fork=3:15;h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';x=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1;d1=z(k)-h1'*c0;c1=c0+k1*d1;e1=c1-c0;e2=e1./c0;e(:,k)=e2;c0=c1;c(:,k)=c1;p1=p0-k1*k1'*[h1'*p0*h1+1];p0=p1;ife2=Ebreak;endendc,ea1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);figure(2);i=1:15;plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':')title('ParameterIdentificationwithRecursiveLeastSquaresMethod')figure(3);i=1:15;plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')title('IdentificationPrecision')程序运行结果:p0=1000000000010000000000100000000001000000c=Columns1through90.001000.00100.4447-1.2248-0.9996-0.9999-1.0001-1.00010.000100.00010.00010.37571.49871.49971.50001.50000.00100-0.24890.63851.05601.00001.00010.99990.99990.001000.25091.13821.55571.49971.49951.49951.4996Columns10through15-1.0001-1.0000-1.0000-1.0000-1.0000-1.00001.50001.50001.50001.50001.50001.50000.99990.99990.99990.99990.99990.99991.49961.49961.49991.49991.49991.4999e=1.0e+003*Columns1through90000.4437-0.0038-0.00020.00000.0000-0.000000003.75640.00300.00000.0000-0.000000-0.2499-0.00360.0007-0.00010.0000-0.00000.0000000.24990.00350.0004-0.0000-0.00000.00000.0000Columns10through150.0000-0.0000-0.00000.0000-0.00000.00000.0000-0.0000-0.00000.0000-0.0000-0.0000-0.00000.0000-0.0000-0.0000-0.00000.0000-0.0000-0.00000.00000.00000.0000-0.0000程序运行曲线:图4.输入信号图5.a1,a2,b1,b2辨识仿真结果图6.a1,a2,b1,b2各次辨识结果的收敛情况结果总结及问题分析:辨识结果与程序运行曲线表明,大约递推到五步时,参数辨识的结果基本达到稳定状态。在参数更改前后运行结果可以看出，输出观测值在没有任何噪声成分下，估计值与真值最大相对误差达到3位数。在编写程序时，发现很多语句都不怎么明白，对程序原理也不是完全了解。在通过对程序的逐句编译后，逐渐了解了每一步的原理与作用，不仅理解了最小二乘递推算法的原理，也学习程序的编写及更多MATLAB语句，比如eye(),zeros(),inv()等函数的应用。