一、 最小二乘法（LS ）辨识系统Z(K+2)=1.5*Z(K+1)-0.7*Z(k)+u(K+1)+0.5*u(k)+v(k) 辨识参数 L T LLTL LS y XX X 1)(-Λ=θ其中MAT程序>> x=[0 1 0 1 1 0 1 1 1];>> n=403; >> M=[]; >> for i=1:ntemp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); endx(1)=temp; end>> v=randn(1,400); >> z=[]; >> z(1)=-1; >> z(2)=0; >> for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); end>> H=zeros(400,4); >> for i=1:400 H(i,1)=-z(i+1); H(i,2)=-z(i); H(i,3)=M(i+1); H(i,4)=M(i); end>> Estimate=inv(H'*H)*H'*(z(3:402))' 辨识参数为： Estimate =-1.49161.03640.4268>>二、最小二乘递推法（RLS）辨识Z(K+2)=1.5*Z(K+1)-0.7*Z(k)+u(K+1)+0.5*u(k)+v(k) 递推公式：其中：MATLAB程序：>> x=[0 1 0 1 1 0 1 1 1];n=403;M=[];for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;endv=randn(1,400);z=[];z(1)=-1;z(2)=0;for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); endP=100*eye(4);Pstore=zeros(4,401);>> Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];>> Theta=zeros(4,401);Theta(:,1)=[3;3;3;3];>> K=[10;10;10;10];h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end>> i=1:401;>> figure(1)>> plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) >> title('待估参数过渡过程')>> figure(2)>> plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) >> title('估计方差变化过程')结果图：三、增广最小二乘法（ELS）辨识Z(K+2)=1.5*Z(K+1)-0.7*Z(k)+u(K+1)+0.5*u(k)-v(k+1)+0.2*v(k) 公式：其中：x=[0 1 0 1 1 0 1 1 1];n=403;M=[];for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;endv=randn(1,402);z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:402z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)-v(i-1)+0.2*v(i-2); endP=100*eye(6);Pstore=zeros(6,401);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];>> Theta=zeros(6,401);>> Theta(:,1)=[3;3;3;3;3;3];>> K=[10;10;10;10;10;10];>> for i=3:402h=[-z(i-1);-z(i-2);M(i-1);M(i-2);v(i-1);v(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(6)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];end>> disp('参数a1、a2、b1、b2、d1、d2估计结果：')参数a1、a2、b1、b2、d1、d2估计结果：>> Theta(:,401)ans =-1.50000.70001.00010.5002-0.99990.2000>> i=1:401;>> figure(1)>> plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:),i,Theta(6,:)); >> title('待估参数过渡过程')>> figure(2)>> plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:),i,Pstore(6,:)); >> title('估计方差变化过程')四、广义最小二乘法（GLS）辨识系统：Z(K+2)=1.5*Z(K+1)-0.7*Z(k)+u(K+1)+0.5*u(k)+e(k) e(k+2)+2.1*e(K+1)-2.5*e(k)=v(k+2)公式：>> x=[0 1 0 1 1 0 1 1 1];n=403;M=[];for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end>> v=randn(1,400);>> e=[];>> e(1)=v(1);>> e(2)=v(2);>> for i=3:400e(i)=0*e(i-1)+0*e(i-2)+v(i);end>> z=zeros(400,1);>> z(1)=-1;z(2)=0;>> for i=3:400z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+e(i);end>> zf=[];>> zf(1)=-1;>> zf(2)=0;>> for i=3:400zf(i)=z(i)-0*z(i-1)-0*z(i-2);end>> uf=[];>> uf(1)=M(1);>> uf(2)=M(2);for i=3:400uf(i)=M(i)-0*M(i-1)-0*M(i-2);end>> P=100*eye(4);>> Pstore=zeros(4,400);>> Pstore(:,2)=[P(1,1),P(2,2),P(3,3),P(4,4)];>> Theta=zeros(4,400);>> Theta(:,2)=[3;3;3;3];>> K=[10;10;10;10];PE=10*eye(2);ThetaE=zeros(2,400);>> ThetaE(:,2)=[0.5,0.3];>> KE=[10;10];>> for i=3:400h=[-zf(i-1);-zf(i-2);uf(i-1);uf(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];he=[-e(i-1);-e(i-2)];KE=PE*he*inv(1+he'*PE*he);ThetaE(:,i)=ThetaE(:,i-1)+KE*(e(i)-he'*ThetaE(:,i-1)); PE=(eye(2)-KE*he')*PE;end>> disp('参数a1、a2、b1、b2估计结果：')参数a1、a2、b1、b2估计结果：>> Theta(:,400)ans =-1.48650.66990.97590.4721>> disp('噪声传递系数c1、c2估计结果：')噪声传递系数c1、c2估计结果：>> ThetaE(:,400)ans =-0.0404-0.0228>> i=1:400;>> figure(1)>> plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) >> title('待估参数过渡过程')>> figure(2)>> plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) >> title('估计方差变化过程')>> figure(3)>> plot(i,ThetaE(1,:),i,ThetaE(2,:));结果图方差变化过程五、辅助变量最小二乘法辨识：Z(K+2)=1.5*Z(K+1)-0.7*Z(k)+u(K+1)+0.5*u(k)+e(k) e(k)=v(k)+0.5*v(k-1)+0.2*v(k-2)递推公式：MATLAB程序：>> x=[0 1 0 1 1 0 1 1 1];n=403;M=[];for i=1:ntemp=xor(x(4),x(9));M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end>> v=randn(1,400);e=[];>> e(1)=0.3;>> e(2)=0.5;>> for i=3:400e(i)=v(i)-0.5*v(i-1)+0.2*v(i-2);end>> z=zeros(402,1);>> z(1)=-1;z(2)=0;for i=3:400z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+e(i); end>> P=100*eye(4);>> Pstore=zeros(4,400);>> Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];>> Theta=zeros(4,400);>> Theta(:,1)=[3;3;3;3];>> Theta(:,2)=[3;3;3;3];>> Theta(:,3)=[3;3;3;3];>> Theta(:,4)=[3;3;3;3];>> K=[10;10;10;10];>> for i=5:400h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];hstar=[-z(i-2-1);-z(i-2-2);M(i-1);M(i-2)];K=P*hstar*inv(h'*P*hstar+1);Theta(:,i)=Theta(:,i-1)+K*(z(i)-h'*Theta(:,i-1));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end>> disp('参数a1、a2、b1、b2估计结果：')参数a1、a2、b1、b2估计结果：>> Theta(:,400)ans =-1.57560.74240.86780.3569>> i=1:400;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程')>> figure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程')。