遞推最小二乘估計及模型階次辨識

1、實驗二 遞推最小二乘估計(RLS)及模型階次辨識(F-Test)1實驗方案設(shè)計i.i生成輸入數(shù)據(jù)和噪聲用M序列作為辨識的輸入信號，噪聲采用標(biāo)準(zhǔn)正態(tài)分布的白噪聲。生成白噪聲時，首先利用乘同余法生成 U0,1均勻分布的隨機數(shù)，再利用U0,1均勻分布的隨機數(shù)生成標(biāo)準(zhǔn)正態(tài)分布的白噪聲。1.2過程仿真辨識模型的形式取A(z)z(k)二B(z)u(k) - e(k)，為方便起見，取na = n n 即A(z) =1 az' a2a* anzB(z) =1 dz' b2a，dz用M序列作為辨識的輸入信號。1.3遞推遺忘因子法方(0) =0.001_1 0 0 01數(shù)據(jù)長度L取534,初值0

2、10 0P(0)=0 0 10kJ0 0 0 1 一1.4計算損失函數(shù)、噪聲標(biāo)準(zhǔn)差損失函數(shù)J(k)二J J(k -1)+【z(k)hk)致k1)2hk)P(k1)h(k) + »噪聲標(biāo)準(zhǔn)差？=J(L)L -dim -1.6 F-Test定階法計算模型階次統(tǒng)計量tt(n,n 1)F(2, L - 2n - 2)J(n) - J(n 1) L -2n -2J(n 1)2其中，J()為相應(yīng)階次下的損失函數(shù)值，L為所用的數(shù)據(jù)長度，n為模型的估計階次。若t(n,n 1) ta，拒絕 H0: n n0，若t(n, n 1) ： ta，接受 H 0: n n0，其中L為風(fēng)險水平:下的閥值。這時模型

3、的階次估計值可取n 1 o1.6計算噪信比和性能指標(biāo)參數(shù)估計平方相對偏差、1八i -彳£ (目)2參數(shù)估計平方根偏差6=堺何呼協(xié)j22編程說明M序列中，M序列循環(huán)周期取 叫=24-1=15，時鐘節(jié)拍t=1Sec,幅度a=1，特征多項式為F (s)二S6二S5二1。白噪聲循環(huán)周期為215 = 32768。G(s)采樣時間To設(shè)為 1Sec，ai 一 -1.5, a2 =0.7, 0=1, b2 = 0.53源程序清單3.1正態(tài)分布白噪聲生成函數(shù)fun ctio nv=no ise(N)%生成正態(tài)分布N(0,sigma)%生成N個0 1均勻分布隨機數(shù)A=179; x0=11; M=2A1

4、5;for k=1:Nx2=A*xO;x仁m od(x2,M);v1=x1/(M+1);v(:,k)=v1;x0=x1;endaipi=v;sigma=1;%標(biāo)準(zhǔn)差for k=1:le ngth(aipi)ksai=0;for i=1:12temp=mod(i+k,le ngth(aipi)+1; ksai=ksai+aipi(temp);endv(k)=sigma*(ksai-6);endend3.2 M序列生成函數(shù)fun ctio nNp r M=createM( n,a)%生成長度為n的M序列，周期為Np，周期數(shù)為r x=1 1 1 1;%初始化初態(tài)for i=1: ny=x;x(2:4

5、)=y(1:3);x(1)=xor(y(1),y (4);U(i)=2* y(4)-a;endM=U*a;len x=le ngth(x);Np=2A|e nx-1;r=n/Np;end3.3加權(quán)最小二乘遞推算法函數(shù)fun ctio nAes,Bes,Error=RLS( na, nb,Z,U,f)%Aes、Bes為參數(shù)估計值，na、nb為模型階次，Z、U為輸出輸入數(shù)據(jù)，f為加權(quán)因子N=na+nb;n_ max=le ngth(Z);X=0.001.*o nes(N,1);P=10A5.*eye(N);e=0.0001; stop=1;n=N;Error=zeros(n_max,1);whil

6、e (stop=1 &&n<=n_ max)H=；for i=1: naH=H;-Z( n-i);endfor j=1:nbH=H;U( n-j);endK=P*H*i nv(H'*P*H+f);X_past=X;X=X+K*(Z( n)-H'*X);P=P-K*K'*(H'*P*H+f);temp=abs(X-X_past)./X_past); stop=sum(temp)>=e;Error( n)=Z( n)-H'*X;n=n+1;endAes=X(1: na)'Bes=X( na+1:N)'%初始估計值%

7、初始P%誤差要求，循環(huán)停止信號%新的數(shù)據(jù)向量%計算增益矩陣%計算新的估計值%計算下次遞推用到的 P%相對誤差%判斷精度3.4主函數(shù)clear%清理工作間變量L=534;%M序列的周期，四級移位寄存器生成 M序列，作為輸入信號u(k)ex=60;%在圖像中展示的數(shù)據(jù)個數(shù)a=1;aa仁-1.5; aa2=0.7; bb仁1; bb2=0.5;%提前規(guī)定的 a,b,c,dNp r u=createM(L,a);%生成 M序列figure(1);%畫第 1 個圖形：u(k)stem(u(1:ex),grid;%以徑的形式顯示出部分輸入信號并給圖形加上網(wǎng)格xlabel( 'k' ) %標(biāo)

8、注橫軸變量ylabel('輸入信號')%標(biāo)注縱軸變量title( '四級移位寄存器生成 M序列輸入信號(前，int2str(ex),'位)')%圖形標(biāo)題z(2)=0;z(1)=0;%取z的前兩個初始值為零y=z;v=no ise(L); lamat=0.1;for k=3:L;%生成白噪聲%循環(huán)變量從3到Ly(k)=-aa1*z(k-1)-aa2*z(k-2)+bb1*u(k-1)+bb2*u(k-2);z(k)=y(k)+lamat*v(k);%給出辨識輸出采樣信號endov=fa ngcha(v); oy=fan gcha(y); yita=sqrt

9、(oyov);%計算噪聲方差%計算信號方差%計算噪信比%用最小二乘遞推算法辨識參數(shù)：a,b,c,deO=O.OO1 0.001 0.001 0.001'%被辨識參數(shù)的初始值采用直接取方式，取一個充分小的實向量p0=10A7*eye(4,4);%初始狀態(tài)P0也采用直接取方式，取一個充分大的實數(shù)單位矩陣E=1e-10;%相對誤差E參考值取0.000000005e=e0,zeros(4,L-1);%被辨識參數(shù)矩陣的初始值及大小eee=zeros(4,L);%相對誤差的初始值及大小n=0;%用于統(tǒng)計遞推次數(shù)for k=3:L;%開始遞推運算hk=-z(k-1),-z(k-2),u(k-1),u

10、(k-2)'%求 h(k)K=p0*hk*i nv(hk'*p0*hk+1);%求 K(k)e1=e0+K*z(k)-hk'*e0;%求 0 (k)ee=(K*(z(k)-hk'*e0);%求相對誤差eee(:,k)=ee;%把當(dāng)前相對變化的列向量加入誤差矩陣的最后一列e0=e1;%新獲得的參數(shù)作為下一次遞推的舊參數(shù)e(:,k)=e1;%把當(dāng)前所辨識參數(shù)的cl列向量加入辨識參數(shù)矩陣的最后一列pk=p0-K*K'*hk'*p0*hk+1;%求 p(k) 值p0=pk;%把當(dāng)前的p(k)值給下次用n=n +1;%完成一次遞推，統(tǒng)計值加1if ee=E

11、 break ;%若參數(shù)收斂滿足要求，終止計算end %小循環(huán)結(jié)束end %大循環(huán)結(jié)束 grid on;a1=e(1,:);a2=e(2,:);b1=e(3,:);ae2=eee(2,:); be1=eee(3,:); be2=eee(4,:);figure(2);%畫第2個圖形i=1:e x;%橫坐標(biāo)從1到Lplot(i,a1(i),'r',i,a2(i),'m' ,i,b1(i),d的各次辨識結(jié)果b2=e(4,:);ae1=eee(1,:);%分離參數(shù)'c',i,b2(i),'g' )%畫出 a , b , c,xlabel(

12、 'k')ylabel('辨識參數(shù)')%標(biāo)注縱軸變量title('最小二乘各次遞推參數(shù)估計值')%圖形標(biāo)題str仁',理論值為：'lege nd( 'a1' ,str1, nu m2str(aa1),'a2' ,str1, nu m2str(aa2),r1, nu m2str(bb1),'b2',str1, nu m2str(bb2);figure(3);%畫第3個圖形i=1:e x;%橫坐標(biāo)從1到Lplot(i,ae1(i),'r',i,ae2(i),'g&

13、#39;,i,be1(i),'b',i,be2(i),a,b,c,d的各次辨識結(jié)果的收斂情況grid on;xlabel('k' ) %標(biāo)注橫軸變量ylabel('參數(shù)誤差')%標(biāo)注縱軸變量legend( 'a1', 'a2', 'b1', 'b2');'b1' ,st'r:')%畫出title('參數(shù)的誤差收斂情況')%圖形標(biāo)題d=1;Us=u(11:L);Zs=v(11:L);arfa=0.01;if (arfa=0.05)T_ar

14、fa=5.3;elseif (arfa=0.05)T_arfa=4.6;elseif (arfa=0.025)T_arfa=3.69;elseif (arfa=0.05)T_arfa=2.99;elseif (arfa=0.1)T_arfa=2.3;else T_arfa=0.1;%舍棄前10個數(shù)據(jù)%置信度%由自由度確定百分點%最大階數(shù)%各階殘差方差 %起始階數(shù) %終止信號endno rder_max=10;J=zeros(1, no rder_max); no rder=d;stop=0;N=0;PEs=zeros(norder_max,norder_max);%存儲 RLS估計結(jié)果%RLS

15、參數(shù)估計while (stop=0&&no rder<=no rder_max)Aes,Bes,Error=RLS( norder, norder,Zs,Us,0.9); theta=Aes,Bes'thetale ngth=le ngth(theta);PEs(1:thetale ngth ,no rder)=theta;J(n order)=Error'*Error;%殘差方差if no rder>1t=(J( norder-1)-J( norder)/J( norder)*(L-2* no rder)/2);if t<=T_arfaN1= no rder-1;N2=N1-d;endendnorder=no rder+1;,num2str(N1), ',The n isenddisp( 'The estimated order by F_test is ',num2str(yita);4曲線打印四級移位寄存器生成財序列諭入信號前和位）圖1驅(qū)動序列：M序列最小二乘各次遞推參數(shù)估計值捫.理論值為:-1.5 a 理論值為? 0 711 理