1、實驗二 最小二乘估計及模型階次辨識一、實驗目的 通過實驗，掌握最小二乘參數辨識方法 通過實驗，掌握模型階次辨識方法二、實驗內容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、辨識算法 模型階次辨識 根據行列式比確定模型的階次

2、令：，其中，定義判別式子（行列式比）：其中：若較有顯著增加時，則取階次估計值為。 模型參數辨識：最小二乘一次完成算法：設輸入信號的取值是從k =1到k =16的M序列，則待辨識參數為=。其中，被辨識參數、觀測矩陣z L、H L的表達式為 , ， 三、實驗步驟(1) 掌握最小二乘一次完成算法和根據行列式比確定模型的階次的基本原理。(2) 設計實驗方案。(3) 編制實驗程序。(4) 調試程序，研究實驗問題，記錄數據。(5) 分析實驗結果，完成實驗報告。四、實驗結果(1) 實驗對象及參數 實驗模型如下圖所示：(2) 程序代碼：a.主函數function leastSquaresMainFuca1 =

3、 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序

4、列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(

5、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

6、(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

7、 為脈沖數目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)

8、=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);%=產生觀

9、測序列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);

10、-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）實驗結果及