離散模型的參數估計及階次辨識

1、學號： 10011314189、10011314188 北京工商大學實驗報告離散模型的參數估計及階次辨識院 （系） 計算機與信息工程學院 專 業(yè) 控制工程 學生姓名 2015年1月離散模型的參數估計及階次辨識一、 實驗目的通過實驗掌握幾種常用的模型參數估計算法和階次辨識理論，并在實驗的基礎上總結體會不同辨識方法的適用范圍和優(yōu)缺點。二、 實驗內容過程辨識書P538實驗2 離散模型的參數估計及階次辨識，選擇模擬的是第一個過程，采用遞推的隨機逼近算法對被辨識系統(tǒng)進行辨識。1. 系統(tǒng)模擬圖如圖1所示 + + 圖 1 系統(tǒng)模擬圖其中 (1)輸入信號采用幅值為1的M序列，其特征多項式由實驗者根據具體需要確

2、定；是均值為零，方差為服從正態(tài)分布的不相關噪聲，方差的大小由信噪比確定。信噪比定義為： (2) 2. 實驗取=0.1和1兩種情況。數據長度L分別取300、500、1000或2000。3. 估計噪聲的方差和模擬靜態(tài)增益。4. 計算參數估計偏差指標 (3)其中，5. 靜態(tài)增益相對偏差 (4)三、 實驗原理隨機逼近法，是一種應用廣泛的參數估計方法。它是在有隨機誤差干擾的情況下，用逐步逼近的方式估計某一特定值的數理統(tǒng)計方法。尋找?guī)д`差的量測到的未知回歸函數的零點或極值，是系統(tǒng)辨識，適應控制、模式識別、適應濾波和神經元網絡等領域中都要遇到的問題。隨機逼近提供了解決這一問題的遞推方法。隨機逼近法的原理就是

3、利用變量x1,x2,，及對應的隨機變量y(x1),y(x2)，通過迭代計算，逐步逼近方程式的解。常用的迭代算法有Robbins-Monro(RM)算法和Keifer-Wolfowitz(KW)算法等?？紤]的參數辨識問題，設準則函數。式中，為某一標量函數；表示k 時刻之前的輸入輸出數據集合。那么，準則函數的一階負梯度為 (5)模型參數辨識問題歸結為求方程的解。根據隨機逼近原理，有 (6)如果準則函數取為 (7)則參數估計的迭代方程可寫為 (8)其中為收斂因子，它必須滿足收斂性條件，一般地，隨著k的增加要有足夠的下降速度，但又不能下降得太快，否則被處理的數據總量太少。四、實驗過程1. 實驗模型由(

4、1)式知， ，再根據題目給出的系統(tǒng)結構圖，可以得到如下的系統(tǒng)模型： (9)即： (10) 由(10)式即可得到系統(tǒng)的待辨識模型為： (11)將(11)式化為最小二乘形式為： (12)其中，為待辨識參數。四、 實驗結果4.1遞推的隨機逼近辨識結果數據長度L=1000，加入=1白噪聲，辨識結果如圖3和表1所示。(a)運行結果：(b)圖3 數據長度L=1000，=1的白噪聲表1 辨識結果(L=1000，=1)a1a2b1b2辨識值-1.45420.68021.00480.5130理論值-1.50000.700010.5000數據長度L=1000，加入=0.1白噪聲，辨識結果如圖4和表2所示。(a)運

5、行結果：(b)圖4 數據長度L=1000，=0.1的白噪聲表2 辨識結果(L=1000，=0.1)a1a2b1b2辨識值-1.49850.69410.97800.5081理論值-1.50000.700010.5000（二）AIC模型定階結果AIC模型定階結果如圖5所示。圖5 AIC模型定階結果四、源程序參數辨識的源程序：success.m%二階系統(tǒng)的遞推的隨機逼近算法%二階系統(tǒng)的模型% z(k)=a1*z(k-1)+a2*z(k-2)+b1*u(k-1)+b2*u(k-2)+v(k)%辨識出參數 theta=a1,a2,b1,b2%清理工作間變量clear allclose allclc%生成

6、輸入的M序列M=1 0 1 0 1 1 0 1 1 1; %寄存器的初始狀態(tài)N = length(M);T = 2N-1;% M序列的周期for i=1:T; m=xor(M(9),M(4); for j=N:-1:2; M(j)=M(j-1); end M(1)=m; u(i) = M(1); if u(i)=0 u(i)=-1; else u(i)=1; endend%產生白噪聲yita=input('請輸入yita的值='); wn=normrnd(0,1,T,1);v=yita*wn;%構造輸出數據L=input('請輸入序列長度L='); z=zero

7、s(1,L); N=2;z(1)=0;z(2)=0;for k=N+1:L; z(k)=1.5*z(k-1)-0.7*z(k-2)+1*u(k-1)+0.5*u(k-2)+v(k);end %初始值置零theta0=10(-4)*diag(eye(2*N,2*N);%置參數theta0的初始值(很小的實向量)cc=theta0,zeros(4,L-1);%遞推算法for k=N+1:Lrouk=5/(k+800);h1=-z(k-1),-z(k-2),u(k-1),u(k-2)'aa=z(k)-h1'*theta0;theta1=theta0+rouk*h1*aa;cc(:,k

8、)=theta1; theta0=theta1;theta=theta1;end thet=-1.5,0.7,1,0.5'theta=theta,theta1=cc(1,:); a2=cc(2,:); b1=cc(3,:); b2=cc(4,:); %畫出參數估計的沒步迭代結果圖i=1:L;figure(1);plot(i,a1,'r',i,a2,'-.',i,b1,'g',i,b2,':')legend('a1','a2','b1','b2')xlabel(

9、'k'); ylabel('Theta'); title('遞推的隨機逼近辨識結果');模型定階的源程序：clear allclose all%生成輸入的M序列M = 1,0,1,0,1,0,0,1,1,1; %設置寄存器的初始狀態(tài)N = length(M);T =2N-1;%反饋邏輯運算模2運算算規(guī)則u = zeros(1,T); for i = 1:T; m = xor(M(10),M(5);% n=10 10 到2 倒著走 for j = N:-1:2 %n=10 10 µ½2 µ¹×

10、97;×ß M(j) = M(j-1);%移位運算 end M(1) = m; u(i) = M(1); if u(i)=0 u(i)=-2; else u(i)=2; endend % z的前三個初始值置零z(1)=0;z(2)=0;v=randn(1,T);%返回一個m*n的隨機項矩陣e(1)=0.6;r=0.9;%控制信噪比for i=2:T e(i)=-0.1*e(i-1)+r*v(i);end%構造輸出數據(或含噪聲數據)v = normrnd(0,1,T,1);%均值為0,方差為1,t行一列for k=3:T; z(k)=1.5*z(k-1)-0.7*z(k-2

11、)+1*u(k-1)+0.5*u(k-2)+v(k);end%準備階次及參數辨識a=3;b=3;length=400;L=max(size(u);aic=zeros(a,b);for n=1:a for m=1:b theta=10(-9)*ones(2*n+m,1); P=106*eye(2*n+m); aa=max(m,n); zf=zeros(1,L); uf=zeros(1,L); vf(1:aa)=0; for k=(aa+1):L h=-z(k-1):-1:(k-n) u(k-1):-1:(k-n) v(k-1):-1:(k-m); v(k)=z(k)-h*theta; vf(k)

12、=v(k)-vf(k-1):-1:(k-m)*theta(2*n+1):(2*n+m); zf(k)=z(k)-zf(k-1):-1:(k-m)*theta(2*n+1):(2*n+m); uf(k)=u(k)-uf(k-1):-1:(k-m)*theta(2*n+1):(2*n+m); hf=-zf(k-1):-1:(k-n) uf(k-1):-1:(k-n) vf(k-1):-1:(k-m); K=P*hf'*inv(hf*P*hf'+1); P=(eye(2*n+m)-K*hf)*P; theta=theta+K*(z(k)-h*theta);, end v=zeros(1,L); for k=(aa+1):L v(k)=z(k)+z(k-1):-1:(k-n) -u(k-1):-1:(k-n) -v(k-1):-1:(k-m)*theta; end sigma=1/L*v(length+1:L)*v(lengt