最小二乘估计及模型阶次辨

1、实验二最小二乘估计及模型阶次辨识、实验目的 通过实验,掌握最小二乘参数辨识方法 通过实验,掌握模型阶次辨识方法、实验内容1、仿真模型实验所用的仿真模型如下：框图表小模型表示z(k)1,5z(k1)0.7z(k2)u(k1)0.5u(k2)v(k)其中u(k)和z(k)分别为模型的输入和输出变量；v(k)为零均值、方差为1、服从正态分布的白噪声；为噪声的标准差(实验时,可取0.0、0.1、0.5、1.0);输入变量u(k)采用M序列,其特征多项式取F(s)s4s1,幅度取1.0.2、辨识模型辨识模型的形式取1_1A(z)z(k)B(z)u(k)e(k)为方便起见,取naAbn,即A(z)1az1

2、a2z2anznB(z)biz1dz2bnzn根据仿真模型生成的数据u(k),k1,L和z(k),k1,L,确定模型n,并辨识模型的参数a1,a2,自和6也,bn;3、辨识算法模型阶次辨识根据行列式比确定模型的阶次令：HnZnUn,其中,z(n)z(n1)zz(n1)z(n)z(2)Znz(nL1)z(nL2)z(L)u(n)u(n1)u(1)u(n1)u(n)u(2)Unu(nL1)u(nL2)u(L)定义判别式子(行列式比)：DR加detH(i?1)其中:八1八1H(?)-HrrHn.,H(?1)-Hr?1H?1假设DR(的较DR(i?1)有显著增加时,那么取阶次估计值为no?模型参数辨识

3、:最小二乘一次完成算法:设输入信号的取值是从k=1到k=16的M序列,那么待辨识参数靠为I=(HLHL)1HlZ.其中,被辨识参数/、观测矩阵zl、HL的表达式为?Ls31a2bib2z(3)zlZ,z(16)Hlz(2)z(1)u(2)u(1)z(3)z(2)u(3)u(2)2(192(14)u(15)u(1.1LS(HlHl)HlZl三、实验步骤(1)掌握最小二乘一次完成算法和根据行列式比确定模型的阶次的基本原理.(2)设计实验方案.(3)编制实验程序.(4)调试程序,研究实验问题,记录数据.(5)分析实验结果,完成实验报告.四、实验结果(1)实验对象及参数实验模型如下列图所示：z(k)1

4、,5z(k1)0.7z(k2)u(k1)0.5u(k2)v(k)(2)程序代码：a.主函数functionleastSquaresMainFuca1=1.5;a2=-0.7;b1=1;b2=0.5;DR=estModelOrder(a1,a2,b1,b2);display(DR);estimate=leastSquares(a1,a2,b1,b2);display(estimate);recursiveLeastSquares(a1,a2,b1,b2)endb.模型阶次辨识函数functionDR=estModelOrder(a1,a2,b1,b2)x=010110111;%initialva

5、luen=1003;%n为脉冲数目,L=1000,且存在k-2,故取1003M=;%存放M序列fori=1:ntemp=xor(x(4),x(9);M(i)=x(9);forj=9:-1:2x(j)=x(j-1);endx(1)=temp;end%产生高斯白噪声v=randn(1,1003);z=;z(1)=-1;z(2)=0;L=1000;fori=3:L+3z(i)=a1*z(i-1)+a2*z(i-2)+b1*M(i-1)+b2*M(i-2)+v(i);end%n=1fori=1:LH1(i,1)=z(i);H1(i,2)=M(i);endA=H1'*H1/L;%n=2fori=

6、1:LH2(i,1)=z(i+1);H2(i,2)=z(i);H2(i,3)=M(i+1);H2(i,4)=M(i);endB=H2'*H2/L;%n=3fori=1:LH3(i,1)=z(i+2);H3(i,2)=z(i+1);H3(i,3)=z(i);H3(i,4)=M(i+2);H3(i,5)=M(i+1);H3(i,6)=M(i);endC=H3'*H3/L;%n=4fori=1:LH4(i,1)=z(i+3);H4(i,2)=z(i+2);H4(i,3)=z(i+1);H4(i,4)=z(i);H4(i,5)=M(i+3);H4(i,6)=M(i+2);H4(i,7)

7、=M(i+1);H4(i,8)=M(i);endD=H4'*H4/L;DR(1)=det(A)/det(B);DR(2)=det(B)/det(C);DR(3)=det(C)/det(D);i=1:3;figure(1)stem(i,DR);%display(DR)title('利用行列式比估计模型阶次')xlabel('阶次')ylabel('行歹！J式比')endc.批量最小二乘估计functionestimate=leastSquares(a1,a2,b1,b2)x=010110111;%initialvaluen=403;%n为脉

8、冲数目M=;%存放M序列fori=1:ntemp=xor(x(4),x(9);M(i)=x(9);forj=9:-1:2x(j)=x(j-1);endx(1)=temp;end%产生均值为0,方差为1的高斯白噪声v=randn(1,400);z=；z(1)=-1；z(2)=0;fori=3:402z(i)=a1*z(i-1)+a2*z(i-2)+b1*M(i-1)+b2*M(i-2)+v(i-2);endH=zeros(400,4);fori=1:400H(i,1)=-z(i+1);H(i,2)=-z(i);H(i,3)=M(i+1);endH(i,4)=M(i);estimate=inv(H

9、'*H)*H'*(z(3:402)'endd.最小二乘的递推算法的参数估计functionrecursiveLeastSquares(a1,a2,b1,b2)x=010110111;%initialvaluen=403;%n为脉冲数目M=;%存放M序列fori=1:ntemp=xor(x(4),x(9);M(i)=x(9);forj=9:-1:2x(j)=x(j-1);endx(1)=temp;end%=产生均值为0,方差为1的高斯白噪声v=randn(1,400);%=产生观测序列z=z=zeros(402,1);z(1)=-1；z(2)=0;fori=3:402z(

10、i)=a1*z(i-1)+a2*z(i-2)+b1*M(i-1)+b2*M(i-2)+v(i-2);end%递推求解P=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=zeros(4,400);%增益矩阵K=10;10;10;10;fori=3:402h=-z(i-1);-z(i-2);M(i-1);M(i-2);K=P*h*inv(h'*P*h+1);Theta(:,i

11、-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);endi=1:401;figure(2);plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)t田e('待估参数过渡过程,);flgure(3);plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title('估计方差变化过程');end(3)