系統(tǒng)辨識與自適應控制--大作業(yè)

1、1 辨識的對象模型假設有一理想數(shù)學模型，它的離散化方程如下式所示： 式中，是服從正態(tài)分布的白噪聲，為系統(tǒng)輸入，為系統(tǒng)輸出。現(xiàn)在輸入信號采用4階序列，其幅值為1。假設系統(tǒng)的模型階次是已知的，即。下面采用遞推最小二乘參數(shù)辨識。 2 遞推最小二乘參數(shù)辨識方法 簡單的最小二乘參數(shù)辨識一次性方法計算復雜，不能夠進行在線辨識，而且所需要的計算存儲空間很大，而很多計算都是重復的計算。為了解決這個問題，并實現(xiàn)在線的實時辨識，引入遞推的最小二乘參數(shù)辨識。遞推最小二乘參數(shù)辨識的整體思想是，最新辨識出來的參數(shù)是建立在上次辨識的參數(shù)基礎上，根據(jù)最新得到的辨識數(shù)據(jù)，對辨識的參數(shù)添加了一個參數(shù)增量。下面利用數(shù)學語言對遞推

2、最小二乘參數(shù)辨識方法進行描述。根據(jù)最小二乘原理，用次觀測數(shù)據(jù)，得出參數(shù)向量的最小二乘估計 (1)其中，表示根據(jù)次觀測數(shù)據(jù)所得到的最小二乘值計量，下表表示該符號代表次觀測數(shù)據(jù)構成的矩陣。 (2)= （3）然后令，且是一方陣，它的維數(shù)取決于未知數(shù)的個數(shù)，而與觀測次數(shù)無關。則 （4）式中表示第次觀測數(shù)據(jù)。利用矩陣反演公式計算（4）式 （5）式中是一個標量，其求逆只是一個簡單的除法。令，于是有，則最小二乘估計量的表達式可改寫為，于是第次觀測數(shù)據(jù)所得到的最小二乘值計量為由（4）式可知 （6）于是 （7） （8）將上述的推導過程總結如下：3 辨識過程整體上的辨識過程流程圖如下圖所示：圖-1 遞推最小二乘參

3、數(shù)辨識流程示意圖詳細步驟如下所示：1 在此選擇4階的M序列。這樣，M序列的最大循環(huán)周期長度為15（bit）。首先利用MATLAB程序生成4階幅值為1的M序列和服從正態(tài)分布的白噪聲序列。2 根據(jù)上述的輸入序列，帶入待辨識的系統(tǒng)模型中，得到一系列對應的輸出信號序列。然后根據(jù)輸入和輸出序列，組建辨識所需的觀測矩陣。本系統(tǒng)的觀測矩陣形式如下：3選取待辨識的參數(shù)向量為模極小的數(shù)值，選取迭代矩陣P的初值為極大的矩陣，如4 計算增益矩陣G其中 5 將上述的初值和增益矩陣G(k)代入上述的遞推最小二乘參數(shù)辨識方程組中進行循環(huán)計算。直到計算的結果精度符合要求或者迭代的次數(shù)達到最大為止。4 辨識結果與分析利用ma

4、tlab對上述系統(tǒng)進行辨識方針，輸入的M序列如下圖所示： 圖-2 M序列信號圖辨識的結果如下圖所示：圖-3 遞推最小二乘參數(shù)辨識結果圖辨識結果與真實的系統(tǒng)參數(shù)對比：參數(shù)a1a2b1b2真值1.80.31.21估計值1.79950.2991.19920.9992對結果的分析：1 利用遞推最小二乘參數(shù)辨識方法，能夠很好的辨識出系統(tǒng)的參數(shù)，并且最后收斂穩(wěn)定。其實，對于穩(wěn)定的系統(tǒng)，遞推最小二乘算法隨著遞推次數(shù)的增加，所得到的辨識結果將會越來越接近真實值。但是對于其他系統(tǒng)，則會存在數(shù)據(jù)飽和的現(xiàn)象，隨著遞推次數(shù)的增加，辨識結果不但不接近真實數(shù)值，反而出現(xiàn)發(fā)散的情況。對于這種情況，需要在遞推最小二乘參數(shù)辨識

5、法的基礎上進一步進行修改，得到限定記憶法和加權最小二乘參數(shù)辨識法。2 本次仿真所采用的M序列階次很小，明顯導致了參數(shù)仿真得到的遞推估計曲線在開始的5秒鐘內(nèi)變化非常劇烈，而且非常不光滑。5 仿真源程序%利用遞推最小二乘法進行參數(shù)辨識%產(chǎn)生M序列L=15;%M序列的周圍位數(shù)y1=1;y2=1;y3=1;y4=0;%四個移位積存器的輸出初始值for i=1:L;%開始循環(huán)，長度為Lx1=xor(y3,y4);%第一個移位積存器的輸入是第3個與第4個移位積存器的輸出的“或”x2=y1;%第二個移位積存器的輸入是第3個移位積存器的輸出x3=y2;%第三個移位積存器的輸入是第2個移位積存器的輸出x4=y3

6、;%第四個移位積存器的輸入是第3個移位積存器的輸出y(i)=y4;%取出第四個移位積存器幅值為"0"和"1"的輸出信號，if y(i)>0.5,u(i)=-1;%如果M序列的值為"1"時,辨識的輸入信號取“-0.03”else u(i)=1;%當M序列的值為"0"時,辨識的輸入信號取“0.03”end%小循環(huán)結束y1=x1;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號做準備end%大循環(huán)結束，產(chǎn)生輸入信號u figure(1); stem(u);grid on;title('輸入信號

7、9;)% 產(chǎn)生白噪聲A=6;x0=1;M=255;for k=1:10000 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*2; x0=x1; v0=v1；endnum=v;k1=k;% 遞推最小二乘辨識程序lamt=1;z(2)=0;z(1)=0;%取z的前兩個初始值為零for k=3:15;%循環(huán)變量從3到15 z(k)=-1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+0*num(1,k);%給出理想辨識輸出采樣信號 endc0=0.001 0.001 0.001 0.001'%直接給出被辨識

8、參數(shù)的初始值,即一個充分小的實向量p0=103*eye(4,4);%直接給出協(xié)方差陣初始狀態(tài)P0，即一個充分大的實數(shù)單位矩陣E=0.000000005;%相對誤差E=0.000000005c=c0,zeros(4,14);%被辨識參數(shù)矩陣的初始值及大小e=zeros(4,15);%相對誤差的初始值及大小for k=3:15; h1=-z(k-1),-z(k-2),u(k-1),u(k-2)' x=h1'*p0*h1+1*lamt; x1=inv(x); %開始求K(k) k1=p0*h1*x1;%求出K的值 d1=z(k)-h1'*c0; c1=c0+k1*d1;%求被

9、辨識參數(shù)c p1=1/lamt*(eye(4)-k1*h1')*p0; e1=c1-c0;%求參數(shù)當前值與上一次的值的差值 e2=e1./c0;%求參數(shù)的相對變化 e(:,k)=e2; %當前相對變化列向量加入誤差矩陣的最后一列 c(:,k)=c1;%把辨識參數(shù)c列向量加入辨識參數(shù)矩陣 c0=c1;%新獲得的參數(shù)作為下一次遞推的舊參數(shù) p0=p1;%給下次用 if e2<=E break;%若參數(shù)收斂滿足要求，終止計算 endend% %得到辨識結果，并作圖a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); %得到辨識結果ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); figure(2);i=1:15;plot(i,a1,'r:',i,a2,'b-.',i,b1,'k-',i,b2,'g',&#