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

实验二 最小二乘估计及模型阶次辨识一、实验目的 通过实验，掌握最小二乘参数辨识方法 通过实验，掌握模型阶次辨识方法二、实验内容1、仿真模型实验所用的仿真模型如下：框图表示e(k) +v(k)u(k)z(k)y(k)模型表示其中u(k)和z(k)分别为模型的输入和输出变量；v(k)为零均值、方差为1、服从正态分布的白噪声；为噪声的标准差（实验时，可取0.0、0.1、0.5、1.0）；输入变量u(k)采用M序列，其特征多项式取,幅度取1.0。2、辨识模型辨识模型的形式取为方便起见，取，即根据仿真模型生成的数据和，确定模型n，并辨识模型的参数； 3、辨识算法 模型阶次辨识 根据行列式比确定模型的阶次令：，其中，定义判别式子（行列式比）：其中：若较有显著增加时，则取阶次估计值为。 模型参数辨识：最小二乘一次完成算法：设输入信号的取值是从k =1到k =16的M序列，则待辨识参数为=。其中，被辨识参数、观测矩阵z L、H L的表达式为 , ， 三、实验步骤(1) 掌握最小二乘一次完成算法和根据行列式比确定模型的阶次的基本原理。(2) 设计实验方案。(3) 编制实验程序。(4) 调试程序，研究实验问题，记录数据。(5) 分析实验结果，完成实验报告。四、实验结果(1) 实验对象及参数 实验模型如下图所示：(2) 程序代码：a.主函数function leastSquaresMainFuca1 = 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.模型阶次辨识函数function DR = estModelOrder(a1,a2,b1,b2)x=0 1 0 1 1 0 1 1 1; %initial valuen=1003; %n为脉冲数目,L = 1000,且存在k-2,故取1003M=; %存放M序列for i=1:n temp=xor(x(4),x(9); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp;end%产生高斯白噪声v=randn(1,1003);z=;z(1)=-1;z(2)=0;L=1000;for i=3:L+3 z(i)=a1*z(i-1)+a2*z(i-2)+b1*M(i-1)+b2*M(i-2)+v(i);end% n=1for i=1:L H1(i,1)=z(i); H1(i,2)=M(i);endA=H1*H1/L;% n=2for i=1:L H2(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=3for i=1:L H3(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=4for i=1:L H4(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)=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(行列式比)endc.批量最小二乘估计function estimate = leastSquares(a1,a2,b1,b2)x=0 1 0 1 1 0 1 1 1; %initial valuen=403; %n 为脉冲数目M=; %存放M 序列for i=1:n temp=xor(x(4),x(9); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp;end%产生均值为0，方差为1 的高斯白噪声v=randn(1,400);z=;z(1)=-1;z(2)=0;for i=3:402 z(i)=a1*z(i-1)+a2*z(i-2)+b1*M(i-1)+b2*M(i-2)+v(i-2);endH=zeros(400,4);for i=1:400 H(i,1)=-z(i+1); H(i,2)=-z(i); H(i,3)=M(i+1);endH(i,4)=M(i);estimate = inv(H*H)*H*(z(3:402);endd.最小二乘的递推算法的参数估计function recursiveLeastSquares(a1,a2,b1,b2)x=0 1 0 1 1 0 1 1 1; %initial valuen=403; %n 为脉冲数目M=; %存放M 序列for i=1:n temp=xor(x(4),x(9); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp;end%=产生均值为0，方差为1 的高斯白噪声=v=randn(1,400);%=产生观测序列z=z=zeros(402,1);z(1)=-1;z(2)=0;for i=3:402 z(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;for i=3:402 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);endi=1:401;figure(2);plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)title(待估参数过渡过程);figure(3);plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title(估计方差变化过程);end（3）实验结果及其分析 1