系統(tǒng)辨識算法

1、傳統(tǒng)系統(tǒng)辨識算法1. 引言迄今為止，已經(jīng)有許多不同的辨識方法。這些辨識方法就其所涉及的模型的形式來說可以分為兩類。一類是非參數(shù)模型的辨識方法，一類是參數(shù)模型的辨識方法。非參數(shù)模型的辨識方法（亦稱經(jīng)典的辨識放法）獲得的模型是非參數(shù)模型。它在假定過程是線性的前提下，不必事先確定模型的具體結(jié)構(gòu)，因而這類方法可適用于任意復(fù)雜的過程，工程上至今仍經(jīng)常采用。參數(shù)模型的辨識方法（亦稱現(xiàn)代的辨識方法）必須假定一種模型結(jié)構(gòu)，通過極小化模型與過程之間的誤差準(zhǔn)則函數(shù)來確定模型的參數(shù)。如果模型的結(jié)構(gòu)無法實(shí)現(xiàn)確定，則必須利用結(jié)構(gòu)辨識方法先確定模型的結(jié)構(gòu)參數(shù)（比如階次、純延遲等），再進(jìn)一步確定模型參數(shù)。參數(shù)模型的辨識方法

2、又可以分為：最小二乘法辨識、梯度校正法辨識以及極大似然法辨識。根據(jù)計(jì)算機(jī)與過程之間的不同聯(lián)接方式，辨識又可以分為離線辨識和在線辨識。離線辨識首先將采集到的數(shù)據(jù)儲存在磁盤或磁帶中，然后將數(shù)據(jù)成批輸入計(jì)算機(jī)進(jìn)行辨識計(jì)算。這種辨識方式多采用成批處理的算法，或稱一次完成算法，其缺點(diǎn)是占用內(nèi)存較大。在線辨識通常要在正常運(yùn)行工況下進(jìn)行，它一般采用實(shí)際處理算法，即每采樣一組數(shù)據(jù)就進(jìn)行一次辨識計(jì)算。這種辨識方式占用內(nèi)存量比較小；尤其對時(shí)變過程的辨識或自適應(yīng)控制問題來說，它比離線辨識方式具有更多的優(yōu)勢。本次作業(yè)使用經(jīng)典的辨識中的一些方法對系統(tǒng)進(jìn)行辨識。在經(jīng)典的控制理論中，線性過程的動態(tài)特性通常用： 傳遞函數(shù) G

3、(s) 頻率響應(yīng) G(jw) 脈沖響應(yīng) g(t) 階躍響應(yīng) h(t)來表示。后三種為非參數(shù)模型，其表現(xiàn)形式是以時(shí)間或頻率為自變量的實(shí)驗(yàn)曲線。對過程施加特定的實(shí)驗(yàn)信號，同時(shí)測定過程的輸出，可以求得這些非參數(shù)模型。經(jīng)過適當(dāng)?shù)臄?shù)學(xué)處理，它們又可以轉(zhuǎn)變成參數(shù)模型傳遞函數(shù)的形式。獲取上述非參數(shù)模型，并把它們轉(zhuǎn)化成傳遞函數(shù)的主要方法有： 階躍響應(yīng)法 脈沖響應(yīng)法 頻率響應(yīng)法 相關(guān)分析法 譜分析法這些辨識方法在工程上有廣泛的應(yīng)用，至今仍受到普遍重視。在本次作業(yè)中，我主要使用了階躍響應(yīng)法中的面積法和脈沖響應(yīng)法來對系統(tǒng)進(jìn)行辨識，并且對系統(tǒng)施加一定的噪聲干擾，比較在有誤噪聲情況下辨識結(jié)果的不同。下面將對其分別進(jìn)行介

4、紹。2. 面積法2.1 面積法的基本原理面積法為階躍響應(yīng)法之中的一種，階躍響應(yīng)法還包括兩點(diǎn)法、切線法、近似法等。在系統(tǒng)的階躍響應(yīng)曲線呈現(xiàn)不規(guī)則的形狀時(shí)，使用面積法比較方便。在使用面積法對系統(tǒng)進(jìn)行辨識，首先就是要測得過程的階躍響應(yīng)，測試路線如圖2-1所示。圖2-1 過程階躍響應(yīng)的測取路線經(jīng)過試驗(yàn)獲取的階躍響應(yīng)要將其轉(zhuǎn)化成無因此的形式，如圖2-2所示：圖2-2 把階躍響應(yīng)轉(zhuǎn)化成無因次的形式設(shè)過程的傳遞函數(shù)為： （2.1.1）其中階次m，n需要根據(jù)先驗(yàn)知識和精度考慮實(shí)現(xiàn)選定。顯然有 （2.1.2）定義： （2.1.3）其中： （2.1.4）那么的Laplace變換為： （2.1.5）則一階面積A1（

5、圖2-2中陰影部分）為： （2.1.6）再令 （2.1.7）定義二階面積為： （2.1.8）同理，令 （2.1.9）可得三階面積為： （2.1.10）以此類推i階面積Ai為： （2.1.11）其中： （2.1.12）進(jìn)一步利用： （2.1.13）可得： （2.1.14）其中： （2.1.15）由式（2.1.5）、（2.1.14）以及可得： （2.1.16）顯然，式（2.1.16）左邊s各次冪項(xiàng)的系數(shù)均為零，固有： （2.1.17）即： （2.1.18）又根據(jù)式（2.1.14）以及可得： （2.1.19）比較上式兩邊s各次冪的系數(shù)，有： （2.1.20）有上式便可以求出傳遞函數(shù)的系數(shù)。應(yīng)當(dāng)指出，

6、為了求（n+m)個(gè)未知數(shù)，需要（n+m）個(gè)方程，因此需要利用式（2.1.18）求階的面積。注意到以及，則式（2.1.20）可以寫成矩陣形式： （2.1.21）由此，便可以利用面積法對系統(tǒng)進(jìn)行辨識。但是，我們可以看出，利用面積法進(jìn)行辨識計(jì)算量是比較大的，但是，計(jì)算機(jī)有著強(qiáng)大的計(jì)算功能，所以，將面積法（主要是式（2.1.18）以及（2.1.21）編制成計(jì)算機(jī)程序之后，再進(jìn)行仿真計(jì)算，就會變得非常簡單。在下一小節(jié)我們將進(jìn)行詳細(xì)的仿真介紹。2.2 MATLAB仿真對于面積法，應(yīng)用MATLAB進(jìn)行仿真時(shí)最主要的幾個(gè)公式在上一節(jié)已經(jīng)提到，就是公式（2.1.8）以及（2.1.21），公式（2.1.8）使用來

7、求取對應(yīng)的面積，公式（2.1.21）使用來求取函數(shù)分子分母多項(xiàng)式的系數(shù)。下面我們對課本中的例題4.1進(jìn)行仿真，運(yùn)行程序如下所示：clear allclcdt=0.01;dt1=0.01;tmax=40;t=0:dt:tmax;M=1+tmax/dt;num=1;den=10,6.5,1;sys=tf(num,den)y=step(sys,t);yn=y/y(M);%把階躍響應(yīng)轉(zhuǎn)化成無因次型A=0,0;%計(jì)算面積A1；A(1)=0;for i=1:M A(1)=A(1)+(1-yn(i)*dt;endA(1);%計(jì)算面積A2A21=0;for i=1:M A21=A21+(1-yn(i)*(-1

8、)*t(i)*dt;endA22=0;for i=1:M A22=A22+(1-yn(i)*dt;endA22=A22*A(1);A(2)=A21+A22;A%bnum1=1;%求aC2=1,0; A(1),1;D2=0;0;den1=;AA=A(1);A(2);den11=(C2*D2+AA)'den1=den11(2),den11(1),1;sysc=tf(num1,den1)yc=step(sysc,t);%沒有噪聲時(shí)的辨識系統(tǒng)yc=yc/yc(M); %加入噪聲之后再進(jìn)行辨識yz=yn+0.5*yn.*rand(size(yn);%yz=yz/yz(M);figure(2)pl

9、ot(t,yz)xlabel('Time Second');ylabel('Step Response');title('加入噪聲之后系統(tǒng)的單位階躍響應(yīng)');Az=0,0;%計(jì)算面積Az1；Az(1)=0;for i=1:M Az(1)=Az(1)+(1-yz(i)*dt;endAz(1);%計(jì)算面積Az2Az21=0;for i=1:M Az21=Az21+(1-yz(i)*(-1)*t(i)*dt;endAz22=0;for i=1:M Az22=Az22+(1-yz(i)*dt;endAz22=Az22*Az(1);Az(2)=Az21+A

10、z22;Az%bznumz1=1;%azCz2=1,0; Az(1),1;Dz2=0;0;denz1=;AAz=Az(1);Az(2);denz11=(Cz2*Dz2+AAz)'denz1=denz11(2),denz11(1),1;syszc=tf(numz1,denz1)yzc=step(syszc,t);%加入噪聲之后的辨識系統(tǒng)yzc=yzc/yzc(M);figure(3)plot(t,y,'r-',t,yc,'b-.',t,yzc,'k')xlabel('Time Second');ylabel('Ste

11、p Response');title('有無噪聲辨識的結(jié)果');legend('原系統(tǒng)階躍響應(yīng)','無噪聲情況下的辨識','加入噪聲之后的辨識')仿真結(jié)果：原系統(tǒng)的傳遞函數(shù)：sys = 1 - 10 s2 + 6.5 s + 1 沒有噪聲情況下辨識系統(tǒng)的傳遞函數(shù)：sysc = 1 - 10.12 s2 + 6.5 s + 1有噪聲情況下辨識系統(tǒng)的傳遞函數(shù)：syszc = 1 - 19.17 s2 + 5.659 s + 1 系統(tǒng)的單位階躍響應(yīng)曲線如圖2-2-1所示：圖2-2-1各系統(tǒng)階躍響應(yīng)曲線仿真分析：我們可以看到，在沒

12、有噪聲的情況下辨識系統(tǒng)和原系統(tǒng)的響應(yīng)曲線完全重合，有噪聲的情況下辨識的結(jié)果就會出先一定的偏差，反復(fù)檢查修改程序之后結(jié)果還是沒有變化，而且這個(gè)程序存在一個(gè)很大的缺點(diǎn)，那就是只能辨識分母為二階，分子為零階的系統(tǒng)，因?yàn)槊娣eA1，A2是分別編寫程序求解的，求解面積的程序并不是一個(gè)通用的程序，想要求解A3，A4或者更多的面積就得重新再加入求解A3，A4的程序，求解系數(shù)矩陣a,b時(shí)所用到的矩陣也是根據(jù)公式人為直接輸入的，并不是通過程序?qū),m的判斷來自動構(gòu)建的。所以，鑒于以上的種種缺點(diǎn)，對程序進(jìn)行了改進(jìn)，改進(jìn)后的程序如下所示：clearclcAs=0;As1=0;Azs=0;Azs1=0;dt=0.005

13、;tmax=40;t=0:dt:tmax;M=1+tmax/dt;num=1;den=10,6.5,1;n=2;%可以更改m=0;%可以更改sys=tf(num,den)%原系統(tǒng)的傳遞函數(shù)y=step(sys,t);yn=y/y(M);%把階躍響應(yīng)轉(zhuǎn)化成無因次型A=;%無噪聲情況下的面積%無噪聲%對求解面積進(jìn)行了改進(jìn)for i=1:n+m Af=0; for j1=1:M Af=Af+(1-yn(j1)*(-t(j1)(i-1)*dt/(factorial(i-1); end for j2=0:i-2 As1=0; As=0; for j3=1:M As1=As1+(1-yn(j3)*(-t(

14、j3)j2)*dt/(factorial(j2); end As=A(i-j2-1)*As1+As; end A(i)=As+Af;endA;%求解b%對求解a,b時(shí)用到的矩陣的構(gòu)建方法進(jìn)行了改進(jìn)Ab=zeros(m,m);for i=1:m; Ab(i,:)=A(n+i-1):-1:(n-m+i);endAb;Anb=A(n+1:n+m);b=-inv(Ab)*Anb'b=b'num1=b(m:-1:1),1;%求解aAa=ones(n,n);for i=1:n for j=1:n if(j>i) Aa(i,j)=0; end endendfor i=1:n for j

15、=1:i-1 Aa(i,j)=A(i-j); endendAa; for i=1:n if(i<=m) Aaa(i,1)=b(i); else Aaa(i,1)=0; end endAaa;Aaaa=A(1:n);a=Aa*Aaa+Aaaa'a=a'den1=a(n:-1:1),1;sysc=tf(num1,den1)%沒有噪聲時(shí)辨識系統(tǒng)的傳遞函數(shù)yc=step(sysc,t);yc=yc/yc(M);%加入噪聲之后再進(jìn)行辨識yz=yn+0.5*yn.*rand(size(yn);figure(1)plot(t,yz)xlabel('Time Second'

16、;);ylabel('Step Response');title('加入噪聲之后系統(tǒng)的單位階躍響應(yīng)');Az=;%有噪聲情況下的面積for i=1:n+m Azf=0; for j1=1:M Azf=Azf+(1-yz(j1)*(-t(j1)(i-1)*dt/(factorial(i-1); end for j2=0:i-2 Azs1=0; Azs=0; for j3=1:M Azs1=As1+(1-yz(j3)*(-t(j3)j2)*dt/(factorial(j2); end Azs=Az(i-j2-1)*Azs1+Azs; end Az(i)=Azs+Az

17、f;endAz;%求解bzAzb=zeros(m,m);for i=1:m; Azb(i,:)=Az(n+i-1):-1:(n-m+i);endAzb;Aznb=Az(n+1:n+m);bz=-inv(Azb)*Aznb'bz=bz'numz1=bz(m:-1:1),1;%求解azAza=ones(n,n);for i=1:n for j=1:n if(j>i) Aza(i,j)=0; end endendfor i=1:n for j=1:i-1 Aza(i,j)=Az(i-j); endendAza; for i=1:n if(i<=m) Azaa(i,1)=b

18、z(i); else Azaa(i,1)=0; end endAzaa;Azaaa=Az(1:n);az=Aza*Azaa+Azaaa'az=az'denz1=az(n:-1:1),1;syszc=tf(numz1,denz1)%加入噪聲之后辨識系統(tǒng)的傳遞函數(shù)yzc=step(syszc,t);yzc=yzc/yzc(M);figure(2)plot(t,y,'r-',t,yc,'b-.',t,yzc,'k')xlabel('Time Second');ylabel('Step Response')

19、;title('有無噪聲辨識的結(jié)果');legend('原系統(tǒng)階躍響應(yīng)','無噪聲情況下的辨識','加入噪聲之后的辨識')改進(jìn)之后程序的仿真結(jié)果為：原系統(tǒng)的傳遞函數(shù)：sys = 1 - 10 s2 + 6.5 s + 1 無噪聲時(shí)辨識系統(tǒng)的傳遞函數(shù)：sysc = 1 - 10.09 s2 + 6.498 s + 1 有噪聲時(shí)辨識系統(tǒng)的傳遞函數(shù)：syszc = 1 - 23.89 s2 + 5.658 s + 1原系統(tǒng)以及各辨識系統(tǒng)單位階躍響應(yīng)曲線如圖2-2-2所示：圖2-2-2各系統(tǒng)單位階躍響應(yīng)曲線仿真分析：程序改進(jìn)之后對原來的系

20、統(tǒng)進(jìn)行仿真辨識，發(fā)現(xiàn)結(jié)果和原系統(tǒng)的幾乎沒有區(qū)別，說明辨識結(jié)果的可信度還是有的。因?yàn)楦倪M(jìn)的程序可以選擇傳遞函數(shù)的階次，所以我們改變原來的系統(tǒng)再次進(jìn)行仿真，將系統(tǒng)的傳遞函數(shù)改成 此時(shí)，令n=4,m=2,得出仿真結(jié)果：原系統(tǒng)的傳遞函數(shù)：sys = 2 s2 + 3 s + 1 - 2 s4 + 12 s3 + 10 s2 + 6.5 s + 1 無噪聲時(shí)辨識系統(tǒng)的傳遞函數(shù)：sysc = 10.37 s2 + 6.049 s + 1 - 3.857 s4 + 41.07 s3 + 29.17 s2 + 9.544 s + 1有噪聲時(shí)辨識系統(tǒng)的傳遞函數(shù)：syszc = -6.301 s2 + 6.308

21、 s + 1 - 4924 s4 + 2.788e04 s3 + 3854 s2 + 8.893 s + 1 圖2-2-3各系統(tǒng)單位階躍響應(yīng)曲線我們可以看到在沒有噪聲的情況下辨識結(jié)果還是不錯(cuò)的，但是在噪聲的影響下辨識效果就變得很差，說明面積法不能很好的消除噪聲的影響。并且，n,m的選擇也很關(guān)鍵，上述仿真是在令n=4,m=2的情況下進(jìn)行的仿真，這是最理想的結(jié)果，但是事實(shí)往往是我們要對m,n的值進(jìn)行推測，如果我們令n=4,m=1，效果便非常的不理想，辨識之后系統(tǒng)的階躍響應(yīng)曲線如圖2-2-4所示，我們可以看出，辨識系統(tǒng)和原系統(tǒng)之間差距非常大，所以在利用面積法求傳遞函數(shù)時(shí)必須正確選擇傳遞函數(shù)的階次。圖

22、2-2-4各系統(tǒng)單位階躍響應(yīng)曲線其次在仿真的過程之中還發(fā)現(xiàn)了一個(gè)問題，當(dāng)原系統(tǒng)系統(tǒng)分子分母之間的階次之差大于2時(shí)，即使m,n選擇正確，也無法正確的對系統(tǒng)進(jìn)行仿真，例如設(shè)系統(tǒng)傳遞函數(shù)為：，我們令n=3,m=0,仿真結(jié)果如圖2-2-5所示，造成這結(jié)果的原因現(xiàn)在還不清楚。圖2-2-5各系統(tǒng)單位階躍響應(yīng)曲線3. 脈沖響應(yīng)法3.1 脈沖響應(yīng)法的基本原理脈沖響應(yīng)是在理想脈沖輸入作用下過程的輸出響應(yīng)?？紤]到工程上實(shí)際輸入理想脈沖是不可能的，因此通常采用矩形脈沖輸入，如圖3-1-11所示：圖3-1-1 過程的脈沖響應(yīng)當(dāng)矩形脈沖的寬度比過程的過渡時(shí)間小得多、且矩形脈沖的面積等于1時(shí)，過程的輸出可以近似為脈沖輸出

23、。由脈沖響應(yīng)確定傳遞函數(shù)的方法有很多，下面介紹幾種常用的方法。1.一階過程如果過程能用一階傳遞函數(shù)描述 ，則傳遞函數(shù)的參數(shù)T和K可以直接在脈沖響應(yīng)曲線上確定，如圖3-1-2所示： （3.1.1）圖3-1-2一階慣性環(huán)節(jié)脈沖響應(yīng)與傳遞函數(shù)參數(shù)的關(guān)系2. 二階過程如果過程能用二階傳遞函數(shù)描述 ，則傳遞函數(shù)的參數(shù)也可以直接由脈沖響應(yīng)曲線確定，如圖3-1-3所示： （3.1.2）圖3-1-3二階慣性環(huán)節(jié)脈沖響應(yīng)與傳遞函數(shù)參數(shù)的關(guān)系3. 差分方程法設(shè)過程的傳遞函數(shù)為： （3.1.3）當(dāng)特征方程具有n個(gè)單根s1,s2,sn時(shí)，傳遞函數(shù)可以寫成 （3.1.4）對應(yīng)的脈沖響應(yīng)為： （3.1.5） 當(dāng)特征方程有

24、n-r個(gè)單根s1,s2,sn-r，r 階重根s0時(shí)，傳遞函數(shù)可寫成 （3.1.6）對應(yīng)的脈沖響應(yīng)為：（3.1.7）為了確定ci和si，從獲取的過程輸出脈沖響應(yīng)g(t)中，選取前n+1個(gè)坐標(biāo)點(diǎn)，每個(gè)坐標(biāo)點(diǎn)間隔相同的采樣時(shí)間T0，如圖所示。各坐標(biāo)點(diǎn)上的脈沖響應(yīng)記為g(k)，g(k+1)，g(k+n)，組成一個(gè)自回歸模型（AR）： （3.1.8）其中為待定系數(shù)。如果特征方程 （3.1.9）有一個(gè)單根 ，則必是AR模型的解，它們的線性組合也是AR模型的解.當(dāng)其特正方程有n個(gè)單根時(shí)，自回歸模型的解為： （3.1.10）當(dāng)其特征方程有n-r個(gè)單根，r階重根時(shí)，自回歸模型的解為： (3.1.11)對照上述公

25、式，不論單根還是重根，都有： （3.1.12）3.2 MATLAB仿真對于脈沖響應(yīng)法一階和二階過程比較簡單，在這里就針對第三種差分方程法進(jìn)行仿真。在第一次進(jìn)行仿真時(shí)和面積法的情況一樣，編寫的程序只能對固定階數(shù)的系統(tǒng)進(jìn)行仿真，在這實(shí)際情況中是不允許的，所以經(jīng)過改進(jìn)的程序可以對任意階的系統(tǒng)進(jìn)行仿真，但是，后來發(fā)現(xiàn)程序有一個(gè)問題，就是當(dāng)特征方程具有重根時(shí)，雖然能夠求得和，但是在求解傳遞函數(shù)時(shí)還沒有編寫出一個(gè)通用的程序，只能是根據(jù)已得到的參數(shù)，然后根據(jù)式（3.1.6）進(jìn)行認(rèn)為的求解傳遞函數(shù)，這是本程序的一個(gè)不足之處。本次的MATLAB仿真是根據(jù)式（3.1.4）、（3.1.8）、（3.1.9）、（3.1

26、.10）、（3.1.12）進(jìn)行編寫的，程序如下所示：clear allclcnum=1.92;den=1,3.8,3.04,1.92;n=3;%可以手動修改階次進(jìn)行辨識直到最佳。disp('原系統(tǒng)傳遞函數(shù)')sys=tf(num,den)t=0:0.01:15;y=impulse(sys,t);%脈沖響應(yīng)；%在T0=1的情況下提取一組數(shù)據(jù),最少需要提取2n個(gè)數(shù)據(jù)a=y(1:100:(2*n-1)*100+1);%取k=0,1，2，.，n-1；A=;for i=1:n A(i,:)=a(i+1:1:i+n);endA;for i=1:n B(i,:)=-a(i);endB;X=A

27、B;%求解方程組的解X=X'den1=;den1=X(n:-1:1),1;P=roots(den1);%求解對應(yīng)的特特征方程的解%求解系統(tǒng)的極點(diǎn)；for i=1:n s(i)=log(P(i);ends;%求取留數(shù)for i=1:n A1(i,:)=(P(1:n).(i-1);endA1;for i=1:n B1(i,:)=a(i);endB1;Y=;Y=A1B1;%求取系數(shù)之后得到傳遞函數(shù)den1=;num1=;for i=1:n zj=1,-s(i); den1(i,:)=zj; num1(i)=Y(i);endsyscom=tf(num1(1),den1(1,:);for i=2

28、:n syscom=parallel(syscom,tf(num1(i),den1(i,:);enddisp('無噪聲時(shí)辨識系統(tǒng)的傳遞函數(shù)')syscomy1=impulse(syscom,t); %加入噪聲之后再辨識zs=y+y.*rand(size(y)*0.5;%加入噪聲figure(1)plot(t,zs)xlabel('Time Second')ylabel('Impulse Response')title('加入噪聲后系統(tǒng)的脈沖響應(yīng)圖')%在T0=1的情況下提取一組數(shù)據(jù)az=zs(1:100:(2*n-1)*100+1);Az=;for i=1:n Az(i,:)=az(i+1:1:i+n);endfor i=1:n Bz(i,:)=-az(i);endBz;Xz=AzBz;%求解方程組的解Xz=Xz'den1z=Xz(n:-1:1),1;Pz=roots(den1z);%求解對應(yīng)的特特征方程的解%求解系統(tǒng)的極點(diǎn)；for i=1:n sz(i)=log(Pz(i);end%求取留數(shù)for i=1:n Az(