系统辨识最小二乘Matlab算法

%辅助变量法--递推算法(包括Tally法和纯滞后) clc,clear; L=2^9-1; lambda=0.3; a1=-1.5;a2=0.7;b1=1;b2=0.5;c1=1;c2=0.41; %给z,e,h,hNew,theta赋初值 z(1)=0.2;z(2)=0.3;e(1)=0.2;e(2)=0.3; p0=10^6*eye(4); h1(:,1)=[0.001;0.001;0.001;0.001]; h1(:,2)=[0.001;0.001;0.001;0.001]; h1(:,3)=[0.001;0.001;0.001;0.001]; h1(:,4)=[0.001;0.001;0.001;0.001]; hNew1(:,1)=[0.001;0.001;0.001;0.001]; hNew1(:,2)=[0.001;0.001;0.001;0.001]; hNew1(:,3)=[0.001;0.001;0.001;0.001]; hNew1(:,4)=[0.001;0.001;0.001;0.001]; theta1(:,1)=[0.001;0.001;0.001;0.001]; theta1(:,2)=[0.001;0.001;0.001;0.001]; theta1(:,3)=[0.001;0.001;0.001;0.001]; theta1(:,4)=[0.001;0.001;0.001;0.001]; p00=10^6*eye(4); h2(:,1)=[0.001;0.001;0.001;0.001]; h2(:,2)=[0.001;0.001;0.001;0.001]; h2(:,3)=[0.001;0.001;0.001;0.001]; h2(:,4)=[0.001;0.001;0.001;0.001]; hNew2(:,1)=[0.001;0.001;0.001;0.001]; hNew2(:,2)=[0.001;0.001;0.001;0.001]; hNew2(:,3)=[0.001;0.001;0.001;0.001]; hNew2(:,4)=[0.001;0.001;0.001;0.001]; theta2(:,1)=[0.001;0.001;0.001;0.001]; theta2(:,2)=[0.001;0.001;0.001;0.001]; theta2(:,3)=[0.001;0.001;0.001;0.001]; theta2(:,4)=[0.001;0.001;0.001;0.001]; %%产生幅值为+1，-1的M序列，作为输入信号u M=[0 1 0 1 0 1 0 1 1]; for i=10:L M(i)=xor(M(i-4),M(i-9)); end u=2*M-1; %产生服从标准正态分布的噪声v v=randn(1,L); %辨识系统模型搭建 for i=3:L e(i)=v(i)+c1*v(i-1)+c2*v(i-2);%有色噪声 z(i)=-a1*z(i-1)-a2*z(i-2)+b1*u(i-1)+b2*u(i-2)+e(i); end for i=5:L h1(:,i)=[-z(i-1);-z(i-2);u(i-1);u(i-2)]; hNew1(:,i)=[-z(i-2-1);-z(i-2-2);u(i-1);u(i-2)];%辅助变量(Tally) k1(:,i)=p0*hNew1(:,i)*inv(h1(:,i)'*p0*hNew1(:,i)+1); p1=(eye(4)-k1(:,i)*h1(:,i)')*p0; p0=p1; theta1(:,i)=theta1(:,i-1)+k1(:,i)*(z(i)-h1(:,i)'*theta1(:,i-1)); end for i=5:L h2(:,i)=[-z(i-1);-z(i-2);u(i-1);u(i-2)]; hNew2(:,i)=[-u(i-2-1);-u(i-2-2);u(i-1);u(i-2)];%辅助变量（纯滞后） k2(:,i)=p00*hNew2(:,i)*inv(h2(:,i)'*p00*hNew2(:,i)+1); p2=(eye(4)-k2(:,i)*h2(:,i)')*p00; p00=p2; theta2(:,i)=theta2(:,i-1)+k2(:,i)*(z(i)-h2(:,i)'*theta2(:,i-1)); end figure(1) i=1:L; plot(i,theta1(1,:),i,theta1(2,:),i,theta1(3,:),i,theta1(4,:));legend('a1','a2','b1','b2'); axis([0 512 -2 1.5]); title('辅助变量法(Tally)的递推算法'); figure(2) i=1:L; plot(i,theta2(1,:),i,theta2(2,:),i,theta2(3,:),i,theta2(4,:));legend('a1','a2','b1','b2'); axis([0 512 -2 1.5]); title('辅助变量法(纯滞后)的递推算法');