系統(tǒng)辨識(shí)之經(jīng)典辨識(shí)法

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

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

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

4、求出n為任意階的個(gè)系數(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）條件可知，當(dāng)t時(shí)， 𝑎𝑑𝑡 (1-6) 同理，定義 𝐹2𝑑𝑡 (1-7) 由式（1-,3）條件可知，當(dāng)t時(shí)， 𝑎𝑑𝑡 (1-8) 因此，可得 𝐹

6、19899;(𝑡) = 0𝑡𝐹𝑛1(𝑡) 𝑎𝑛1𝑦(𝑡) dt (1-9) 𝑎𝑛 = 𝐹𝑛() (1-10) 當(dāng)系統(tǒng)的傳遞函數(shù)存在零點(diǎn)時(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 ，以此類(lèi)推，得到i 階面積Ai 。進(jìn)一步利用est 拉氏變換，得到L1h*（t ）=M sii ，進(jìn)而得到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之間的關(guān)系，如下： é ù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ù)的參數(shù)a, b。 1.2實(shí)驗(yàn)過(guò)程 1.2.1無(wú)零點(diǎn)模型系統(tǒng) 假設(shè)系統(tǒng)的傳遞函數(shù)模型為G(s) =2 1 ,為無(wú)零點(diǎn)的模型，利用10&#

12、119904; +6.5𝑠+1Matlab 編程，分別在沒(méi)有噪聲和有噪聲兩種情況下進(jìn)行辨識(shí)，比較辨識(shí)結(jié)果。 1.沒(méi)有噪聲時(shí)，程序如下： clear; %=獲得原傳遞函數(shù)方程=% num=1; den=10 6.5 1; %=產(chǎn)生階躍采樣序列=% T=0.2; %采樣周期 t=0:T:30; %采樣時(shí)間 L=length(t); %數(shù)據(jù)長(zhǎng)度 h=step(num,den,t); %原傳遞函數(shù)的階躍響應(yīng) 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('通過(guò)辨識(shí)得到的傳遞函數(shù)為：') G2=tf(num1,den1) %=原傳遞函數(shù)和辨識(shí)函數(shù)的階躍響應(yīng)對(duì)比圖=% step(G1,'b-',G2,'r-.') title('原系統(tǒng)與辨識(shí)后所得到系統(tǒng)階躍響應(yīng)對(duì)比') legend('原響應(yīng)曲線(xiàn)','辨識(shí)響應(yīng)曲線(xiàn)') （1） 當(dāng)采樣周期T=0.2秒

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

15、igure(1) w=randn(1,L); % 建立服從正態(tài)分布的隨機(jī)矩陣。 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） 當(dāng)采樣周期T=0.02s,采樣時(shí)間t=50s時(shí)，辨識(shí)結(jié)果如下： 圖1.8 原系統(tǒng)與辨識(shí)系統(tǒng)階躍響應(yīng)如圖所示： 結(jié)合上述無(wú)測(cè)量噪聲和有測(cè)量噪聲兩種情況下的辨識(shí)結(jié)果，列出如下所示的表格：表1-1 噪聲情況 條件 增益 a1 a2 參數(shù) 采樣時(shí)間 數(shù)據(jù)長(zhǎng)度 1.0 10 6.5 真值 無(wú)測(cè)量

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

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

18、程序如下： clear; %=獲得原傳遞函數(shù)方程=% num=17.5 7.5 1; den=4 5 8 1; %=產(chǎn)生階躍采樣序列=% T=0.02; %采樣周期 t=0:T:100; %采樣時(shí)間 L=length(t); %數(shù)據(jù)長(zhǎng)度 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í)后傳遞函數(shù)=% num1=b2 b1 1; den1=a3 a2 a1 1; disp('原傳遞函數(shù)為:') G1=tf(num,den) disp('通過(guò)辨識(shí)得到的傳遞函數(shù)為：') G2=tf(num1,den1) %=原傳遞函數(shù)和辨識(shí)函數(shù)的階躍響應(yīng)對(duì)比圖=% step(G1,'b-',G2,'r-

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

22、數(shù) 采樣時(shí)間 數(shù)據(jù)長(zhǎng)度 1.0000 4 5 8 17.5 7.5 真值 0.02 100 1.0000 3.559 5.224 8.051 17.71 7.515 估計(jì)值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 分析：通過(guò)對(duì)比不同的采樣周期和不同的采樣時(shí)間在無(wú)測(cè)量噪聲情況的辨識(shí)結(jié)果可知，對(duì)于存在有零點(diǎn)的系統(tǒng)來(lái)說(shuō)，通過(guò)面積法辨識(shí)系統(tǒng)必須合理的選擇分子分母的階次，否則不能得出正確的辨識(shí)結(jié)果。從表格中也可以發(fā)現(xiàn)，在相同的采樣時(shí)間下，由于系統(tǒng)收斂過(guò)程較長(zhǎng)從，增

23、加數(shù)據(jù)長(zhǎng)度對(duì)系統(tǒng)的辨識(shí)精度有了很大的提高；同樣在相同的數(shù)據(jù)長(zhǎng)度下，合理的減小采樣時(shí)間，也可以提高系統(tǒng)的辨識(shí)精度。 2.脈沖響應(yīng)法 脈沖響應(yīng)法（ impulse response method ）是指線(xiàn)性系統(tǒng)在零初始條件輸入脈沖信號(hào)，信號(hào)后輸出的瞬態(tài)響應(yīng)，即輸出響應(yīng)叫脈沖響應(yīng)。建立系統(tǒng)的非參數(shù)模來(lái)觀(guān)測(cè)系統(tǒng)脈沖響應(yīng)，以求得系統(tǒng)數(shù)學(xué)模型的待定參數(shù)，進(jìn)而實(shí)現(xiàn)系統(tǒng)辨識(shí)的目的。 2.1 Hankel矩陣法 一個(gè) n 階過(guò)程的脈沖傳遞函數(shù)為 G(z )=11 1+b za z2222+LL+b zna znnn (2-1) b z 11+a z1+將傳遞函數(shù)轉(zhuǎn)化成為狀態(tài)方程后，進(jìn)一步推導(dǎo)，可知 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û 脈沖相應(yīng)法參數(shù)計(jì)算公式如下： é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),并且它是可逆的。因此，只需要將獲得的脈沖響應(yīng)值g(k)，k=1,2,3,2n填入式（3-3）、（3-4）便可以求出脈沖傳遞函數(shù)的估計(jì)值，即傳遞函數(shù)系數(shù)：𝑎1, 𝑎2 , 𝑎&

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

30、t); %原傳遞函數(shù)的脈沖響應(yīng) %=Hamkel矩陣求解參數(shù)=% H=gk(1) gk(2) gk(3);gk(2) gk(3) gk(4);gk(3) gk(4) gk(5);%構(gòu)造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中的值倒轉(zhuǎn)，a=a1,a2,a3 A=1 0 0;a(1) 1 0;a(2) a(1) 1; %構(gòu)造下三角陣 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('辨識(shí)得到的Z函數(shù)和傳遞函數(shù)') dsys2=filt(b.',a,0.02) %求得辨識(shí)得到的Z函數(shù) G2=T*d2c(dsys2,'tustin',0.02) %雙線(xiàn)性變換 %=原傳遞函數(shù)和辨識(shí)函數(shù)的脈沖響應(yīng)對(duì)比圖=% impulse(G1,'b-',G2,'r-.') %繪出原系統(tǒng)與辨識(shí)得到的系統(tǒng)的脈沖響應(yīng)曲線(xiàn) title('原模型與辨識(shí)后所得到系統(tǒng)脈沖響應(yīng)對(duì)比') legend('原始響應(yīng)曲線(xiàn)','辨識(shí)曲線(xiàn)') 程序運(yùn)行結(jié)果如下： 原系統(tǒng)的傳遞函數(shù)和脈沖階躍函數(shù)如下： 圖2.1 辨識(shí)之后系統(tǒng)的傳遞函數(shù)和脈沖函數(shù)如下： 圖2.2 在沒(méi)有噪聲的情況下，原系統(tǒng)與辨識(shí)系統(tǒng)脈沖響應(yīng)對(duì)比圖： 圖2.3 （2）在噪聲期望為 0，方差為 0.01 的情況下，對(duì)系統(tǒng)進(jìn)行參數(shù)辨識(shí)，在 Mtlab中編程實(shí)現(xiàn)，程序如下： clear %=獲得原傳遞函數(shù)方程=% num=5 1; den=1 15 50 100;