系統(tǒng)辨識之經典辨識法

1、 系統(tǒng)辨識作業(yè)一 學 院 信息科學與工程學院 專 業(yè) 控制科學與工程 班 級 控制二班 姓 名 學 號 2018 年 11 月 系統(tǒng)辨識 所謂辨識就是通過測取研究對象在認為輸入作用的輸出響應，或正常運行時的輸入輸出數(shù)據(jù)記錄，加以必要的數(shù)據(jù)處理和數(shù)學計算，估計出對象的數(shù)學模型。 辨識的內容主要包括四個方面： 實驗設計； 模型結構辨識； 模型參數(shù)辨識； 模型檢驗。 辨識的一般步驟：根據(jù)辨識目的，利用先驗知識，初步確定模型結構；采集數(shù)據(jù)；然后進行模型參數(shù)和結構辨識；最終驗證獲得的最終模型。 根據(jù)辨識方法所涉及的模型形式來說，辨識方法可以分為兩類：一類是非參數(shù)模型辨識方法，另一類是參數(shù)模型辨識方法。

2、其中，非參數(shù)模型辨識方法又稱為經典的辨識方法，它主要獲得的是模型是非參數(shù)模型。在假定過程是線性的前提下，不必事先確定模型的具體結構，廣泛適用于一些復雜的過程。經典辨識方法有很多，其中包括階躍響應法、脈沖響應法、相關分析法和普分析法等等，本次實驗所采用的辨識方法為階躍響應法和脈沖響應法。 1階躍響應法 階躍響應法是一種常用非參數(shù)模型辨識方法。常用的方法有近似法、半對數(shù)法、切線法、兩點法和面積法等。本次作業(yè)采用面積法求傳遞函數(shù)。 1.1面積法 當系統(tǒng)的傳遞函數(shù)無零點時，即系統(tǒng)傳遞函數(shù)如下： G(S) = 𝑎𝑛𝑠𝑛+𝑎&

3、#119899;1𝑠𝑛11+𝑎1𝑠+1 (1-1) 系統(tǒng)的傳遞函數(shù)與微分方程存在著一一對應的關系，因此，可以通過求取微分方程的系數(shù)來辨識系統(tǒng)的傳遞函數(shù)。在求得系統(tǒng)的放大倍數(shù)K后，要得到無因次階躍響應y(t)(設=0)，其中y(t)用下式描述： 𝑑𝑛𝑦(𝑡)𝑛1(𝑡)𝑎𝑛𝑑𝑡𝑑𝑡𝑑𝑡 (1-2) 面積法原則上可以

4、求出n為任意階的個系數(shù)。以n為3為例。有： 𝑑3𝑦(𝑡)𝑑2𝑦(𝑡)𝑑𝑦(𝑡)𝑑𝑡|𝑡 =𝑑𝑡|𝑡 =𝑑𝑡|𝑡 = 0 (1-3) 𝑦(𝑡)|𝑡 = 1將式（1）中的y(t)移至右邊，在0，t上積分，得 𝑑2𝑦(𝑡)

5、𝑎3𝑑𝑡𝑑𝑡 (1-4) 定義： 𝐹1(𝑡) = 0𝑡1 𝑦(𝑡)𝑑𝑡 (1-5) 由式（1-3）條件可知，當t時， 𝑎𝑑𝑡 (1-6) 同理，定義 𝐹2𝑑𝑡 (1-7) 由式（1-,3）條件可知，當t時， 𝑎𝑑𝑡 (1-8) 因此，可得 𝐹

6、19899;(𝑡) = 0𝑡𝐹𝑛1(𝑡) 𝑎𝑛1𝑦(𝑡) dt (1-9) 𝑎𝑛 = 𝐹𝑛() (1-10) 當系統(tǒng)的傳遞函數(shù)存在零點時，傳遞函數(shù)如下： G（s）=k×b smmn +ba smn-1-1smn-1-1 +LL +a sbs11 +1+1，（n ³ m） （1-11） a sn+其中，K h= ( )¥ / U0 定義 D1G(s)=K P(s

7、)其中， P(s) = b sa sn mn +ba smn-1-1smn-1-1 +LL +a sbs11 +1+1 = +1 åi¥=1 C sii （1-12） m根據(jù)1h*(t)的Laplace變換，求出一階面積A1，確定Lh（*1 t ），并定義二階面積A2 ，以此類推，得到i 階面積Ai 。進一步利用est 拉氏變換，得到L1h*（t ）=M sii ，進而得到Ai 的值： i=0Ai = ò0¥éë1h*(t)ûù (i 1)!t)i1dt +åtj=20 Ai j 1ò0

8、5;éë1h*(t)ùû (j!t) j dt （1-13） (根據(jù)A Ci = i ，可得：𝑎𝑛𝑠𝑛 + 𝑎𝑛1𝑠𝑛1 + + 𝑎1𝑠 + 1 = (𝑏𝑚𝑠𝑚 + 𝑏𝑚1𝑠𝑚1 + + 𝑏1𝑠 + 1)(1 + 𝑖=1 &

9、#119860;𝑖𝑠𝑖)。比較上式兩邊s的各次冪，便可得到a, b, A之間的關系，如下： é ùb1é AnAn1 L An m +1ù é1 An+1 ùê úê úb2êê An+1An L An m +2ú êú êAn+2 úú =ê úMê LLLL ú ê Múê úë

10、; ûbmêëAn m+ 1 An m+ 2 LAn ú êû ëAn m+ úûé ùb1é ù éa1110 LL 0 00 0ùúê úê úbM2 +éêê úAA12úù （1-14） ê ú êê ú êa2 = A1úê úê

11、 ú êMLL L L L úê ú ê úMê ú êúê úbmê úAnûë û ëanAn1 An2 L A1 1ûê úë û0ë由此可知，根據(jù)式(1-12)、(1-13)、(1-14)便可得到辨識傳遞函數(shù)的參數(shù)a, b。 1.2實驗過程 1.2.1無零點模型系統(tǒng) 假設系統(tǒng)的傳遞函數(shù)模型為G(s) =2 1 ,為無零點的模型，利用10&#

12、119904; +6.5𝑠+1Matlab 編程，分別在沒有噪聲和有噪聲兩種情況下進行辨識，比較辨識結果。 1.沒有噪聲時，程序如下： clear; %=獲得原傳遞函數(shù)方程=% num=1; den=10 6.5 1; %=產生階躍采樣序列=% T=0.2; %采樣周期 t=0:T:30; %采樣時間 L=length(t); %數(shù)據(jù)長度 h=step(num,den,t); %原傳遞函數(shù)的階躍響應 K=h(L) %系統(tǒng)增益 %=面積法求解參數(shù)=% s1=0; for i=1:L s1=s1+(1-h(i)*T; F(i)=s1; end a1=s1; s2=0; for i=

13、1:L s2=s2+(F(i)-a1*h(i)*T; end a2=s2; num1=1; den1=a2 a1 1; disp('原傳遞函數(shù)為:') G1=tf(num,den) disp('通過辨識得到的傳遞函數(shù)為：') G2=tf(num1,den1) %=原傳遞函數(shù)和辨識函數(shù)的階躍響應對比圖=% step(G1,'b-',G2,'r-.') title('原系統(tǒng)與辨識后所得到系統(tǒng)階躍響應對比') legend('原響應曲線','辨識響應曲線') （1） 當采樣周期T=0.2秒

14、，采樣時間t=30s時，行程序后得到原傳遞函數(shù)G1和辨識得到的傳遞函數(shù)G2如圖1.1： 圖1.1 原系統(tǒng)和辨識后系統(tǒng)的階躍響應對比圖如下： 圖1.2 （2） 當采樣周期T=0.2秒，采樣時間t=50s時，行程序后得到原傳遞函數(shù)G1和辨識得到的傳遞函數(shù)G2如下： 圖1.3 原系統(tǒng)和辨識后系統(tǒng)的階躍響應對比圖如下： （3） 當采樣周期T=0.02秒，采樣時間t=50s時，行程序后得到原傳遞函數(shù)G1 和辨識得到的傳遞函數(shù)G2如下： 圖1.5 原系統(tǒng)和辨識后系統(tǒng)的階躍響應對比圖如下： 2.有噪聲的情況下，系統(tǒng)程序如下：主程序還是用面積法，在程序中添加以下代碼： %產生期望為0，方差為0.01的噪聲 f

15、igure(1) w=randn(1,L); % 建立服從正態(tài)分布的隨機矩陣。 w=w/std(w); w=w-mean(w); qw=0; fc=0.01; w=qw+sqrt(fc)*w; %=階躍采樣序列中加入白噪聲=% h=h+w; plot(t,w); （1） 加入的噪聲如下圖所示： 圖1.7 （2） 當采樣周期T=0.02s,采樣時間t=50s時，辨識結果如下： 圖1.8 原系統(tǒng)與辨識系統(tǒng)階躍響應如圖所示： 結合上述無測量噪聲和有測量噪聲兩種情況下的辨識結果，列出如下所示的表格：表1-1 噪聲情況 條件 增益 a1 a2 參數(shù) 采樣時間 數(shù)據(jù)長度 1.0 10 6.5 真值 無測量

16、噪聲 0.2 30 0.9985 11.43 6.594 估計值 0.2 50 1.0000 11.31 6.6 0.02 50 1.0000 10.13 6.51 有測量噪 (方差為 0.01) 0.02 50 1 7.784 6.709 分析：通過對比不同的采樣周期和不同的采樣時間在無測量噪聲情況的辨識結果可知，在相同的采樣周期下，適當?shù)脑黾硬蓸訒r間，可以提高辨識精度，尤其是對增益的提高有很大影響；而在相同的采樣時間下，適當?shù)臏p小采樣時間，對于系統(tǒng)參數(shù)的辨識精度有很大的提高。因此，可以發(fā)現(xiàn)合理采樣時間和數(shù)據(jù)長度，可以提高辨識的精度，令辨識后的傳遞函數(shù)系數(shù)與原傳遞函數(shù)系數(shù)更接近，差距小，從而

17、得到滿意的辨識結果。 通過對比無測量噪聲和有測量噪聲兩種情況下的辨識結果，我們可以發(fā)現(xiàn)在白噪聲的情況下，曲線擬合較無噪聲情況下要差，說明白噪聲對于面積法辨識系統(tǒng)存在較大的干擾，會對辨識結果產生一定的影響。 1.2.2有零點模型系統(tǒng) 17.5𝑆2+7.5𝑆+1 假設系統(tǒng)的傳遞函數(shù)為G(s) = 4𝑆3+5𝑆2+8𝑆+1,為有零點的模型，其中n=3,m=2, 用面積法需要求解𝐴1𝐴5，利用 Matlab 編程，分別在沒有噪聲和有噪聲兩種情況下進行辨識，比較辨識結果。 （1） 沒有噪聲時，

18、程序如下： clear; %=獲得原傳遞函數(shù)方程=% num=17.5 7.5 1; den=4 5 8 1; %=產生階躍采樣序列=% T=0.02; %采樣周期 t=0:T:100; %采樣時間 L=length(t); %數(shù)據(jù)長度 y=step(num,den,t); k=y(L) %系統(tǒng)增益 %=面積法求解參數(shù)=% sum1=0; for i=1:L-1; sum1=sum1+(1-(y(i)+y(i+1)/2)*T; A(i)=sum1; end A1=sum1 sum2=0; for i=1:L-1; sum2=sum2+(A(i)-A1*(y(i)+y(i+1)/2)*T; B(

19、i)=sum2; end A2=sum2 sum3=0; for i=1:L-1; sum3=sum3+(B(i)-A2*(y(i)+y(i+1)/2)*T; C(i)=sum3; end A3=sum3 sum4=0; for i=1:L-1; sum4=sum4+(C(i)-A3*(y(i)+y(i+1)/2)*T; D(i)=sum4; end A4=sum4 sum5=0; for i=1:L-1; sum5=sum5+(D(i)-A4*(y(i)+y(i+1)/2)*T; end A5=sum5 %=根據(jù)所得A(i)，利用公式求取a、b=% M=(-1)*(inv(A3,A2;A4,

20、A3)*A4;A5; b1=M(1,1); b2=M(2,1); N=1 0 0;A1 1 0;A2 A1 1*b1;b2;0+A1;A2;A3; a1=N(1,1); a2=N(2,1); a3=N(3,1); %=根據(jù)所求a、b，得到辨識后傳遞函數(shù)=% num1=b2 b1 1; den1=a3 a2 a1 1; disp('原傳遞函數(shù)為:') G1=tf(num,den) disp('通過辨識得到的傳遞函數(shù)為：') G2=tf(num1,den1) %=原傳遞函數(shù)和辨識函數(shù)的階躍響應對比圖=% step(G1,'b-',G2,'r-

21、.') title('原系統(tǒng)與辨識后所得到系統(tǒng)階躍響應對比') legend('原響應曲線','辨識響應曲線') 當采樣時間取0.02秒，數(shù)據(jù)長度為100時，辨識結果如下： 圖1.10 原系統(tǒng)與辨識后的系統(tǒng)階躍響應對比圖： 當采樣時間為0.02，數(shù)據(jù)長度為400時，系統(tǒng)辨識結果如下： 圖1-12 原系統(tǒng)與辨識后的系統(tǒng)階躍響應對比圖： 圖1-13 當采樣時間為0.2秒，數(shù)據(jù)長度為400時，系統(tǒng)辨識結果如下： 圖1-14 原系統(tǒng)與辨識后的系統(tǒng)階躍響應對比圖： 綜上所述，結果如表 表1-2 條件 增益 a 3 a 2 a 1 b 2 b 1 參

22、數(shù) 采樣時間 數(shù)據(jù)長度 1.0000 4 5 8 17.5 7.5 真值 0.02 100 1.0000 3.559 5.224 8.051 17.71 7.515 估計值0.02 400 1.0000 4.016 4.914 7.979 17.42 7.479 0. 2 400 1.0000 4.129 4.168 7.792 16.61 7.278 分析：通過對比不同的采樣周期和不同的采樣時間在無測量噪聲情況的辨識結果可知，對于存在有零點的系統(tǒng)來說，通過面積法辨識系統(tǒng)必須合理的選擇分子分母的階次，否則不能得出正確的辨識結果。從表格中也可以發(fā)現(xiàn)，在相同的采樣時間下，由于系統(tǒng)收斂過程較長從，增

23、加數(shù)據(jù)長度對系統(tǒng)的辨識精度有了很大的提高；同樣在相同的數(shù)據(jù)長度下，合理的減小采樣時間，也可以提高系統(tǒng)的辨識精度。 2.脈沖響應法 脈沖響應法（ impulse response method ）是指線性系統(tǒng)在零初始條件輸入脈沖信號，信號后輸出的瞬態(tài)響應，即輸出響應叫脈沖響應。建立系統(tǒng)的非參數(shù)模來觀測系統(tǒng)脈沖響應，以求得系統(tǒng)數(shù)學模型的待定參數(shù)，進而實現(xiàn)系統(tǒng)辨識的目的。 2.1 Hankel矩陣法 一個 n 階過程的脈沖傳遞函數(shù)為 G(z )=11 1+b za z2222+LL+b zna znnn (2-1) b z 11+a z1+將傳遞函數(shù)轉化成為狀態(tài)方程后，進一步推導，可知 G(z )=

24、g(1)z11+g(2)z2 +g(3)z3 +L 根據(jù)上述公式，可得 bz1 1 +b z2 2 +L +b zn n=éëg(1)z1 +g(2)z2 +g(3)z3 +L ùéûë1+a z1 1 +a z2 2 +L +a zn n ûù= g(1)z1 +L +g(2)+ag1 (1)z2 +L +êég(n)+ån1 g(i)an i ùúzn (2-2) ëi=1ûn2n1+êég(n+ +1) åg(

25、i)an+i ùúz +(n 1) +L +êég(2n)+åg(i)a2n i ùúz2n ëi=1ûëi=1û 脈沖相應法參數(shù)計算公式如下： ég(1)g(2) Lg n( ) ùé an ùég n( +1)ùêg(2)g(3) L g n( +1) úúêêan1úú=êêg n( +2)úú (2-3) &#

26、234;ê MM OM úê Múê M úêúêúêúëg n( ) g n( +!) L g n(2 1)ûë a1 ûë g n(2 ) ûé ù éb110 L 00ùég(1)ùê ú êê ú êb2 = a11 L 00úúêêg(2)ú

27、;úê ú êMM M O M Múê Múê ú êúêúë û ëbnan1 an2 L a11ûëg n( )û (2-4) Hankel矩陣的定義： 𝑔(𝑘)𝑔(𝑘 + 1)𝑔(𝑘 + 𝑙 1) 𝑔(𝑘 + 1)𝑔(𝑘 + 2

28、)𝑔(𝑘 + 𝑙)H(l, k) = () (2-5) 𝑔(𝑘 + 𝑙 1)𝑔(𝑘 + 𝑙)𝑔(𝑘 + 2𝑙 2)公式（2-3）左邊的矩陣是一種特定的Hankel矩陣。記作H(n,l),并且它是可逆的。因此，只需要將獲得的脈沖響應值g(k)，k=1,2,3,2n填入式（3-3）、（3-4）便可以求出脈沖傳遞函數(shù)的估計值，即傳遞函數(shù)系數(shù)：𝑎1, 𝑎2 , 𝑎&

29、#119899; , 𝑏1, 𝑏2 , 𝑏𝑛。 2.2.實驗過程 假設系統(tǒng)模型的傳遞函數(shù)為G(s) = 𝑠3+155𝑠2𝑠+61𝑠+100, 利用Matlab編程，分別在沒有噪聲和有噪聲兩種情況下進行辨識，比較辨識結果。 (1)無噪聲的情況下，程序如下： clear %=獲得原傳遞函數(shù)方程=% num=5 1; den=1 15 50 100; %=產生脈沖采樣序列=% T=0.02; %采樣周期 t=0:T:30; %采樣時間 gk=impulse(num,den,

30、t); %原傳遞函數(shù)的脈沖響應 %=Hamkel矩陣求解參數(shù)=% H=gk(1) gk(2) gk(3);gk(2) gk(3) gk(4);gk(3) gk(4) gk(5);%構造Hankel矩陣求a、b系數(shù)矩陣 H1=H-1; %H求逆 G1=-gk(4);gk(5);gk(6); a=H1*G1; %這里a=a3;a2;a1 a=flipud(a).' %a中的值倒轉，a=a1,a2,a3 A=1 0 0;a(1) 1 0;a(2) a(1) 1; %構造下三角陣 G2=gk(1);gk(2);gk(3); b=A*G2; a=1 a; disp('原傳遞函數(shù)及其Z函數(shù)

31、') G1=tf(num,den) dsys1=c2d(G1,T) disp('辨識得到的Z函數(shù)和傳遞函數(shù)') dsys2=filt(b.',a,0.02) %求得辨識得到的Z函數(shù) G2=T*d2c(dsys2,'tustin',0.02) %雙線性變換 %=原傳遞函數(shù)和辨識函數(shù)的脈沖響應對比圖=% impulse(G1,'b-',G2,'r-.') %繪出原系統(tǒng)與辨識得到的系統(tǒng)的脈沖響應曲線 title('原模型與辨識后所得到系統(tǒng)脈沖響應對比') legend('原始響應曲線','辨識曲線') 程序運行結果如下： 原系統(tǒng)的傳遞函數(shù)和脈沖階躍函數(shù)如下： 圖2.1 辨識之后系統(tǒng)的傳遞函數(shù)和脈沖函數(shù)如下： 圖2.2 在沒有噪聲的情況下，原系統(tǒng)與辨識系統(tǒng)脈沖響應對比圖： 圖2.3 （2）在噪聲期望為 0，方差為 0.01 的情況下，對系統(tǒng)進行參數(shù)辨識，在 Mtlab中編程實現(xiàn)，程序如下： clear %=獲得原傳遞函數(shù)方程=% num=5 1; den=1 15 50 100;