过程辨识最小二乘算法Matlab程序

%多级最小二乘法 clc clear L=400; a=1; disp('辨识对象： z(k)+0.9*z(k-1)+0.15*z(k-2)+0.02*z(k-3)=0.7*u(k-1)-1.5*u(k-2)+e(k);') disp('e(k)+1.0*e(k-1)+0.41*e(k-2)=λ*v(k)') na=3;nb=2;nc=2; %~~~~~~~~~~~~~~~~输入采用4阶M序列，幅值为1 x=[0 1 1 0]; for i=5:L x(i)=xor(x(i-1),x(i-4)); end for i=1:L u(i)=a*(1-2*x(i)); end %~~~~~~~~~~~~~~~白噪声 A=179; M=2^35; x(1)=11; for i=1:L x(i+1)=mod(A*x(i),M); v(i)=x(i)/M; end % v=rand(1,L); lamd=0.92; %~~~~~~~~~~~~~ e=zeros(1,L); for i=3:L e(i)=-1.0*e(i-1)-0.41*e(i-2)+lamd*v(i); end %~~~~~~~~~~~~~~输出z z=zeros(1,L); for i=4:L z(i)=-0.9*z(i-1)-0.15*z(i-2)-0.02*(i-3)+0.7*u(i-1)-1.5*u(i-2)+e(i); end %~~~~~~~~~~~~~~~~辅助模型参数辨识 H1=zeros(L-5,9); for i=6:L H1(i-5,:)=[-z(i-1),-z(i-2),-z(i-3),-z(i-4),-z(i-5),u(i-1),u(i-2),u(i-3),u(i-4)];%H1 end H11=H1'; z1=z(6:L)'; w=inv(H11*H1); theta1=w*H11*z1 %~~~~~~~~~~~~~~~~过程模型参数辨识 e1(1:5)=theta1(1:5)'; f=theta1(6:9)'; o=zeros(1,na); z2=[f,o]'; for i=1:na+nb+nc%~~~~~~~~~~构造H2 for j=1:na if i-j>4 h21(i,j)=0; elseif i>j h21(i,j)=-f(i-j); end end end for i=1:na+nb+nc for j=1:nb if i-j>5 h22(i,j)=0; elseif i==j h22(i,j)=1; elseif i>j h22(i,j)=e1(i-j); end end end H2=[h21,h22]; H22=H2'; w2=inv(H22*H2); theta2=w2*H22*z2 %~~~~~~~~~~~~~~~噪声模型参数辨识 a=theta2(1:3); b=theta2(4:5); for i=1:na z31(i)=e1(i)-a(i); end for i=1:nc z32(i)=e1(na+i); end for i=1:nb-1 z33(i)=f(1+i)-b(1+i); end for i=1:nc z34(i)=f(nb+i); end z3=[z31,z32,z33,z34]'; for i=1:na+nc %~~~~~~~~~~~H3 for j=1:nc if i-j>3 h31(i,j)=0; elseif i==j h31(i,j)=1; elseif i>j h31(i,j)=a(i-j); end end end for i=1:nb+nc-1 for j=1:nc if i-j>=2 h32(i,j)=0; elseif i>=j h32(i,j)=b(i-j+1); end end end H3=[h31;h32]; H33=H3'; w3=inv(H33*H3); theta3=w3*H33*z3 %辨识结果 % theta2 = % % 0.9321 % 0.2013 % -0.0022 % 0.7448 % -1.4670 % % % theta3 = % % -0.1636 % -0.7986 % 结果a3 及theta3存在很大误差 ，未能解决。