最小二乘递推算法的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; %四个移位寄存器的初始值for i=1:L; %移位循环x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4; %取出作为输出信号，即M序列if y(i)>0.5,u(i)=-0.03; %输入信号else u(i)=0.03;endy1=x1;y2=x2;y3=x3;y4=x4;endfigure(1);stem(u),grid onz(2)=0;z(1)=0;for k=3:15;z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %输出采样信号endc0=[0.001 0.001 0.001 0.001]'; %直接给出被识别参数的初始值p0=10^6*eye(4,4); %直接给出初始状态P0E=0.000000005;c=[c0,zeros(4,14)];e=zeros(4,15);for k=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;if e2<=E break;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('Parameter Identification with Recursive Least Squares Method')figure(3);i=1:15;plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')title('Identification Precision')程序运行结果：p0 =1000000 0 0 00 1000000 0 00 0 1000000 00 0 0 1000000c =Columns 1 through 90.0010 0 0.0010 -0.4984 -1.2325 -1.4951 -1.4962 -1.4991 -1.49980.0001 0 0.0001 0.0001 -0.2358 0.6912 0.6941 0.6990 0.69980.0010 0 0.2509 1.2497 1.0665 1.0017 1.0020 1.0002 0.99990.0010 0 -0.2489 0.7500 0.5668 0.5020 0.5016 0.5008 0.5002Columns 10 through 15-1.4999 -1.5000 -1.5000 -1.5000 -1.4999 -1.49990.6999 0.7000 0.7000 0.7000 0.7000 0.70000.9998 0.9999 0.9999 0.9999 0.9999 0.99990.5002 0.5000 0.5000 0.5000 0.5000 0.5000e =1.0e+003 *Columns 1 through 90 0 0 -0.4994 0.0015 0.0002 0.0000 0.0000 0.00000 0 0 0 -2.3592 -0.0039 0.0000 0.0000 0.00000 0 0.2499 0.0040 -0.0001 -0.0001 0.0000 -0.0000 -0.00000 0 -0.2499 -0.0040 -0.0002 -0.0001 -0.0000 -0.0000 -0.0000Columns 10 through 150.0000 0.0000 0.0000 -0.0000 -0.0000 0.00000.0000 0.0000 -0.0000 0.0000 0.0000 0.0000-0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000-0.0000 -0.0000 0.0000 -0.0000 0.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;for i=1:L;x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;if y(i)>0.5,u(i)=-0.03;else u(i)=0.03;endy1=x1;y2=x2;y3=x3;y4=x4;endfigure(1);stem(u),grid onz(2)=0;z(1)=0;for k=3:15;z(k)=z(k-1)-1.5*z(k-2)+u(k-1)+1.5*u(k-2);endc0=[0.001 .0001 0.001 0.001]';p0=10^6*eye(4,4)E=0.000000005;c=[c0,zeros(4,14)];e=zeros(4,15);for k=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;if e2<=E break;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('Parameter Identification with Recursive Least Squares Method')figure(3);i=1:15;plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:')title('Identification Precision')程序运行结果:p0 =1000000 0 0 00 1000000 0 00 0 1000000 00 0 0 1000000c =Columns 1 through 90.0010 0 0.0010 0.4447 -1.2248 -0.9996 -0.9999 -1.0001 -1.00010.0001 0 0.0001 0.0001 0.3757 1.4987 1.4997 1.50001.50000.0010 0 -0.2489 0.6385 1.0560 1.0000 1.0001 0.99990.99990.0010 0 0.2509 1.1382 1.5557 1.4997 1.4995 1.49951.4996Columns 10 through 15-1.0001 -1.0000 -1.0000 -1.0000 -1.0000 -1.00001.5000 1.5000 1.5000 1.5000 1.5000 1.50000.9999 0.9999 0.9999 0.9999 0.9999 0.99991.4996 1.4996 1.4999 1.4999 1.4999 1.4999e =1.0e+003 *Columns 1 through 90 0 0 0.4437 -0.0038 -0.0002 0.0000 0.0000 -0.00000 0 0 0 3.7564 0.0030 0.0000 0.0000 -0.00000 0 -0.2499 -0.0036 0.0007 -0.0001 0.0000 -0.0000 0.00000 0 0.2499 0.0035 0.0004 -0.0000 -0.0000 0.0000 0.0000Columns 10 through 150.0000 -0.0000 -0.0000 0.0000 -0.0000 0.00000.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000-0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000-0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000程序运行曲线:图4. 输入信号图5. a1,a2,b1,b2辨识仿真结果图6. a1,a2,b1,b2各次辨识结果的收敛情况结果总结及问题分析:辨识结果与程序运行曲线表明,大约递推到五步时,参数辨识的结果基本达到稳定状态。

在参数更改前后运行结果可以看出，输出观测值在没有任何噪声成分下，估计值与真值最大相对误差达到3位数。