一種高階線性時(shí)不變系統(tǒng)的辨識算法仿真

吳相甫

（航空工業(yè)中國飛機(jī)強(qiáng)度研究所，西安710065）

在現(xiàn)代控制理論系統(tǒng)辨識的研究領(lǐng)域里，線性系統(tǒng)建模是控制學(xué)科的重要課題。在古典控制中，常常通過測試被控對象的單位階躍響應(yīng)來求其傳遞函數(shù)[1-3]。利用階躍響應(yīng)曲線來確定典型工業(yè)過程傳遞函數(shù)的方法很多，常用的有近似法、半對數(shù)法、切線法、兩點(diǎn)法和面積法等。當(dāng)階躍響應(yīng)曲線比較規(guī)則時(shí)，近似法、切線法就、半對數(shù)法和兩點(diǎn)法都能比較有效地導(dǎo)出傳遞函數(shù)[4-6]。而在現(xiàn)代控制中，伴隨著自適應(yīng)控制算法的完善和發(fā)展，以系統(tǒng)階躍響應(yīng)或脈沖響應(yīng)為依據(jù)而進(jìn)行建模的算法也應(yīng)運(yùn)而生，如動(dòng)態(tài)矩陣控制算法，它的模型是建立在對象的單位階躍響應(yīng)的基礎(chǔ)之上。目前，工程上用得最多的是用階躍擾動(dòng)和方波擾動(dòng)來測取階躍響應(yīng)。

從某種意義上講，階躍響應(yīng)建模法提供了過程數(shù)學(xué)模型建立的一種解決方案。階躍響應(yīng)曲線形象地表示了對象的動(dòng)態(tài)特性，但對于進(jìn)一步地分析和研究系統(tǒng)卻很不方便，因此需要從階躍響應(yīng)曲線求出對象的傳遞函數(shù)，過去常用的方法是圖表法和積分法。圖表法包括切線法和兩點(diǎn)法，但都有它的局限性對于有自平衡能力的對象，要在階躍響應(yīng)曲線上找出拐點(diǎn)，在其拐點(diǎn)處作切線，由于拐點(diǎn)位置不易選準(zhǔn)，且切線的方向斜率也難于確定，以致傳遞函數(shù)的特征參量不能準(zhǔn)確確定。再者，根據(jù)曲線所得的數(shù)據(jù)使傳遞函數(shù)不便于運(yùn)算和模擬，給傳遞函數(shù)的確定帶來了一定的誤差，同樣，數(shù)值積分法雖是一種從階躍響應(yīng)曲線求出傳遞函數(shù)的通用數(shù)學(xué)方法，可以適用于各種形式的傳遞函數(shù)，但是由于數(shù)值計(jì)算的誤差，實(shí)際上只能有效地定出少量的常數(shù)[7-8]。

本文利用系統(tǒng)階躍響應(yīng)的曲線特征來確定傳遞函數(shù)，思路是利用傳遞函數(shù)可以分解為一二階典型環(huán)節(jié)的并聯(lián)原理，假設(shè)模型的階數(shù)后，調(diào)節(jié)典型環(huán)節(jié)并聯(lián)模型中的參數(shù)值，用擬合曲線逼近實(shí)際的階躍響應(yīng)曲線，從而得到系統(tǒng)的傳遞函數(shù)[9]。研究時(shí)借助Matlab 平臺開發(fā)辨識軟件，通過曲線擬合確定典型環(huán)節(jié)并聯(lián)模型的參數(shù)值，可以在一定精度內(nèi)獲得系統(tǒng)的傳遞函數(shù)。

1 高階線性時(shí)不變系統(tǒng)模型分析

一般情況下，高階線性時(shí)不變系統(tǒng)的傳遞函數(shù)分子分母都是s 的多項(xiàng)式，可以寫為

式（1）可表示成如下形式：

式中：m＜n，q+2r=n

由式（2）可知，任何高階線性時(shí)不變系統(tǒng)都可以分解為若干個(gè)一階系統(tǒng)和二階系統(tǒng)的典型環(huán)節(jié)之和。

另外，式（1）可表示為如下因式的乘積形式：

在輸入為單位階躍函數(shù)時(shí)，輸出量可表示為

式中：q 為實(shí)數(shù)極點(diǎn)的個(gè)數(shù)；r 為共軛復(fù)數(shù)極點(diǎn)的對數(shù)。設(shè)0<ζk<1，將式（5）展成部分分式

設(shè)初始條件全部為零，將式（6）進(jìn)行拉氏反變換，可得高階系統(tǒng)的單位階躍響應(yīng)：

式（7）表明，高階系統(tǒng)的時(shí)間響應(yīng)是由一階系統(tǒng)和二階系統(tǒng)的時(shí)間響應(yīng)疊加組成的。

綜上所述，任何線性高階系統(tǒng)都可以分解為若干個(gè)一階和二階典型環(huán)節(jié)的并聯(lián)，且其響應(yīng)也可以分解成若干個(gè)一階和二階典型環(huán)節(jié)的階躍響應(yīng)之和。更進(jìn)一步，通過一個(gè)一階系統(tǒng)和一個(gè)二階系統(tǒng)的典型環(huán)節(jié)的并聯(lián)模型的階躍響應(yīng)曲線逼近實(shí)際系統(tǒng)的階躍響應(yīng)曲線，通過調(diào)節(jié)參數(shù)

改變曲線形狀，不斷地逼近高階系統(tǒng)的階躍響應(yīng)曲線，能夠反映原系統(tǒng)的整體特性[10]，即可用一個(gè)一階系統(tǒng)和一個(gè)二階系統(tǒng)典型環(huán)節(jié)的并聯(lián)模型去等效高階線性系統(tǒng)的傳遞函數(shù)。

2 階躍響應(yīng)曲線特征分析

傳遞函數(shù)的參數(shù)值對系統(tǒng)的動(dòng)態(tài)性能和穩(wěn)態(tài)性能有著直接的影響，各參數(shù)的取值情況直接決定系統(tǒng)的階躍響應(yīng)曲線特征。一階系統(tǒng)與二階系統(tǒng)典型環(huán)節(jié)并聯(lián)模型可表示為

各參數(shù)對階躍響應(yīng)曲線形狀影響的規(guī)律如下：

增益值的影響系統(tǒng)階躍響應(yīng)穩(wěn)態(tài)值為增益值K1、K2之和，K1、K2對上升時(shí)間、峰值時(shí)間、調(diào)節(jié)時(shí)間無影響；K1、K2之和不變時(shí)，K2越大，峰值越大，即超調(diào)量越大。

時(shí)間常數(shù)的影響T1、T2由小變大時(shí)，峰值時(shí)間和調(diào)節(jié)時(shí)間變大，這是因?yàn)橐浑A環(huán)節(jié)響應(yīng)速度變慢；T2相對于T1會(huì)使響應(yīng)的速度變慢的更快。按照T1、T2的比例關(guān)系將該結(jié)構(gòu)的三階系統(tǒng)分為3 類：T1≥5T2、T2＜T1＜5T2、T1＜T2，時(shí)間常數(shù)對其階躍響應(yīng)曲線形狀的影響如下[11]：

（1）T1≥5T2： 階躍響應(yīng)曲線的第一個(gè)波峰值要小于穩(wěn)態(tài)值。二階系統(tǒng)在階躍響應(yīng)曲線前期產(chǎn)生作用，一階系統(tǒng)在后期產(chǎn)生作用。

（2）T2＜T1＜5T2： 階躍響應(yīng)曲線的第一個(gè)波峰值大于穩(wěn)態(tài)值。動(dòng)態(tài)部分上升速度慢，超調(diào)小甚至沒有超調(diào)，但曲線上有拐點(diǎn)。

（3）T1＜T2： 階躍響應(yīng)曲線的第一個(gè)波峰值大于穩(wěn)態(tài)值。一階環(huán)節(jié)發(fā)生作用比二階環(huán)節(jié)早，一階環(huán)節(jié)在上升段初期發(fā)生作用，二階環(huán)節(jié)產(chǎn)生影響的時(shí)間靠后。

阻尼比的影響阻尼比ζ 增大時(shí)，超調(diào)減小，振蕩次數(shù)減少。峰值時(shí)間tp由其中的二階系統(tǒng)決定，根據(jù)二階系統(tǒng)峰值時(shí)間的計(jì)算方法，T2與ζ 滿足關(guān)系式

該式在K1、K2之和不變時(shí)，K1、K2的變化不影響該結(jié)論。

3 辨識方法

根據(jù)上一節(jié)的分析，階躍響應(yīng)穩(wěn)態(tài)值由K1、K2之和決定，T1、T2決定曲線的響應(yīng)速度，阻尼比決定曲線的超調(diào)。按照T1、T2的比例關(guān)系對階躍響應(yīng)曲線特征的影響將分3 類：T1≥5T2、T2＜T1＜5T2、T1＜T2，分析辨識方法。

T1≥5T2這種情況下，一階環(huán)節(jié)的時(shí)間常數(shù)比二階環(huán)節(jié)時(shí)間常數(shù)大得多，因此一階環(huán)節(jié)發(fā)生作用比二階環(huán)節(jié)晚。此時(shí)，階躍響應(yīng)曲線的第一個(gè)波峰值要小于穩(wěn)態(tài)值，該特征可用來判斷此系統(tǒng)屬于T1≥5T2這一類。

具體參數(shù)辨識算法如下：

（1）一階時(shí)間常數(shù)T1：由于一階環(huán)節(jié)發(fā)生作用要比二階環(huán)節(jié)要晚，曲線會(huì)先發(fā)生振蕩，然后曲線會(huì)平緩上升直到達(dá)到穩(wěn)態(tài)值，因此可利用一階系統(tǒng)求時(shí)間常數(shù)的方法求出T1，即在階躍響應(yīng)曲線上找到輸出量變化至終值95%時(shí)的坐標(biāo)點(diǎn)，它所對應(yīng)的時(shí)刻為調(diào)節(jié)時(shí)間ts，則T1=ts/3。

（2）ζ、T2、K1、K2：這4 個(gè)參數(shù)要由循環(huán)搜索逼近的方法確定。

算法包括兩個(gè)循環(huán)，外層循環(huán)調(diào)節(jié)ζ、T2，內(nèi)層循環(huán)調(diào)節(jié)K1、K2：首先，求出峰值時(shí)間tp，令阻尼比ζ=0.05，由式（10）可以算出T2，并令K1=K2=K/2，在該組數(shù)值下求出階躍響應(yīng)數(shù)據(jù)，算出第一個(gè)波峰差p=y1-y1′（實(shí)際階躍曲線值減去擬合曲線值），若p>0，則令K2=K2+0.1，K1=K1-0.1（若p<0，則K2=K2-0.1，K1=K1+0.1），此為內(nèi)循環(huán)（用來調(diào)節(jié)K1、K2）；當(dāng)p 的絕對值小于0.1 時(shí)，在當(dāng)前參數(shù)值之下算出第一個(gè)波谷差bg=y2-y2′（實(shí)際階躍曲線值減去擬合曲線值），若bg>0.1，則說明阻尼比ζ 的值太?。ㄕ袷幪螅?，此時(shí)令ζ=ζ+0.05，重新計(jì)算T2，此為外循環(huán)（用來調(diào)節(jié)ζ、T2），并令K1=K2=K/2，再根據(jù)波峰差p 調(diào)節(jié)K1、K2。參數(shù)辨識流程見圖1。

圖1 參數(shù)辨識流程（T1≥5T2）Fig.1 Flow chart of parameter identification（T1≥5T2）

T1＜T2這種情況下，一階環(huán)節(jié)的時(shí)間常數(shù)比二階環(huán)節(jié)時(shí)間常數(shù)小，因此一階環(huán)節(jié)發(fā)生作用比二階環(huán)節(jié)早。此時(shí)，階躍響應(yīng)曲線的第一個(gè)波峰值大于穩(wěn)態(tài)值，該特征可用來判斷此系統(tǒng)不屬于T1≥5T2類。

在T1＜T2以及T2＜T1＜5T2時(shí)都首先假設(shè)該系統(tǒng)為二階系統(tǒng)，求出阻尼比ζ、時(shí)間常數(shù)T2、增益K，在這3 個(gè)參數(shù)下求出這個(gè)擬合的階躍響應(yīng)擬合值y′，用實(shí)際階躍響應(yīng)曲線值y 減去y′，求出動(dòng)態(tài)誤差d，這樣得到的動(dòng)態(tài)誤差d 的值會(huì)很大。經(jīng)過多組實(shí)驗(yàn)仿真，當(dāng)d>0.4 時(shí)，將待辨識系統(tǒng)歸為T1＜T2類；當(dāng)d<0.4 時(shí)，將待辨識系統(tǒng)歸為T2＜T1＜5T2類。

對 于T1＜T2類，ζ、T1、T2、K1、K2這5個(gè)參數(shù)由3個(gè)循環(huán)搜索逼近的方法確定，具體參數(shù)辨識算法如下：

（1）ζ、T2： 這2 個(gè)參數(shù)由最外層循環(huán)搜索確定。由于含1 個(gè)一階環(huán)節(jié)，該環(huán)節(jié)會(huì)減緩系統(tǒng)的振蕩，因此由假設(shè)的二階系統(tǒng)求出的阻尼比ζ 是偏大的，因此要向下調(diào)整，T2則由式（10）計(jì)算得出。

（2）T1、K1、K2：這3 個(gè)參數(shù)由中層和最內(nèi)層循環(huán)搜索確定。第一步確定了阻尼比ζ 和T2的值后，令T1=T2，K1=K2=K/2。由于T1＜T2，T1應(yīng)向下調(diào)整，每次減小0.1（此值可調(diào)整），約束條件為T1>0.1。每賦T1一個(gè)值時(shí)，K1、K2在最內(nèi)層循環(huán)從K/2 開始根據(jù)峰值差p 調(diào)整，算出峰值差p=y1-y1′（實(shí)際階躍曲線值減去擬合曲線值），若p>0，則令K2=K2+0.1，K1=K1-0.1；若p<0，則令K2=K2-0.1，K1=K1+0.1。

（3）當(dāng)T1減小至接近0.1 時(shí)，若誤差仍然有d>0.05，則另阻尼比ζ 減小0.05，重新計(jì)算，然后再進(jìn)入第二步，如此循環(huán)搜索直至d<0.05。

參數(shù)辨識流程見圖2。

圖2 參數(shù)辨識流程圖（T1＜T2）Fig.2 Flow chart of parameter identification（T1＜T2）

T2＜T1＜5T2這種情況與T1＜T2相同，階躍響應(yīng)曲線的第一個(gè)波峰值大于穩(wěn)態(tài)值，該特征可用來判斷此系統(tǒng)不屬于T1≥5T2類。建模思想同T1＜T2相同，首先認(rèn)為待辨識系統(tǒng)為二階系統(tǒng)，并算出動(dòng)態(tài)誤差d，當(dāng)d<0.4 時(shí)，將待辨識系統(tǒng)歸為T2＜T1＜5T2類。

對于T2＜T1＜5T2類，ζ、T1、T2、K1、K2這5 個(gè) 參 數(shù)同樣由3 個(gè)循環(huán)搜索逼近的方法確定，具體參數(shù)辨識算法如下：

（1）ζ、T2：同T1＜T2類辨識算法第一步；

（2）T1、K1、K2：這3 個(gè)參數(shù)由中層和最內(nèi)層循環(huán)搜索確定。第一步確定了阻尼比ζ 和T2的值后，令T1=T2，K1=K2=K/2。由于T2＜T1＜5T2，T1應(yīng)向上調(diào)整，每次增加0.1（此值可調(diào)整）。每賦T1一個(gè)值時(shí)，K1、K2在最內(nèi)層循環(huán)從K/2 開始根據(jù)峰值差p 調(diào)整，算出峰值差p=y1-y1′（實(shí)際階躍曲線值減去擬合曲線值），若p>0，則令K2=K2+0.1，K1=K1-0.1；若p<0，則令K2=K2-0.1，K1=K1+0.1。當(dāng)K1、K2的調(diào)整使擬合模型的峰值滿足要求時(shí)，再計(jì)算已知階躍響應(yīng)和擬合模型的第一個(gè)波谷差bg，判斷當(dāng)前阻尼比是否正確。

（3）當(dāng)T1調(diào)整到使估計(jì)模型的峰值與穩(wěn)態(tài)值大約相等時(shí)，說明T1已足夠大，此時(shí)若仍有動(dòng)態(tài)誤差d>0.05、bg 的絕對值大于0.15 則進(jìn)入第一步重新開始搜索。

參數(shù)辨識流程見圖3。

圖3 參數(shù)辨識流程圖（T2＜T1＜5T2）Fig.3 Flow chart of parameter Identification（T2＜T1＜5T2）

4 辨識仿真

為了驗(yàn)證辨識算法，利用Matlab 圖形用戶界面GUI 來設(shè)計(jì)辨識軟件，將辨識算法嵌入到軟件操作界面中。基于某實(shí)驗(yàn)平臺搭建模擬電路，如圖4 所示，輸入方波信號，用示波器觀察響應(yīng)曲線，當(dāng)該電路的階躍響應(yīng)進(jìn)入穩(wěn)態(tài)值后，點(diǎn)擊示波器的停止按鈕，保存并導(dǎo)出階躍響應(yīng)數(shù)據(jù)，保存為Excel 格式文件，并導(dǎo)入到辨識軟件中，計(jì)算出模型中的參數(shù)。

圖4 模擬電路圖Fig.4 Analog circuit

該電路由1 個(gè)一階系統(tǒng)與1 個(gè)二階系統(tǒng)并聯(lián)而成，當(dāng)電阻、電容參數(shù)值不相同時(shí)，這類結(jié)構(gòu)的三階系統(tǒng)可分為3 種情況。

T1≥5T2模擬電路各參數(shù)取值：R1=500 kΩ，C1=2 μF，R01=250 kΩ；R2=50 kΩ，C2=1 μF；R3=200 kΩ，C3=1 μF，R02=200 kΩ。

一階環(huán)節(jié)的傳遞函數(shù)為

式中：時(shí)間常數(shù)T1為1，K1為2。

二階環(huán)節(jié)的開環(huán)傳遞函數(shù)為

那么該二階環(huán)節(jié)的閉環(huán)傳遞函數(shù)為

式中：K2為1，T2為0.1，阻尼比ζ 為0.25。

故該三階系統(tǒng)的傳遞函數(shù)為

辨識結(jié)果如圖5 所示，根據(jù)辨識算法計(jì)算出的參數(shù)值與原電路模型參數(shù)值基本接近。

T1＜T2模擬電路各參數(shù)取值：R1=100 kΩ，C1=1 μF，R01=50 kΩ；R2=500 kΩ，C2=1 μF；R3=500 kΩ，C3=1 μF，R02=500 kΩ。

一階環(huán)節(jié)的傳遞函數(shù)為

圖5 辨識結(jié)果（T1≥5T2）Fig.5 Identification results（T1≥5T2）

式中：時(shí)間常數(shù)T1為0.1，K1為2。

二階環(huán)節(jié)的開環(huán)傳遞函數(shù)為

那么該二階環(huán)節(jié)的閉環(huán)傳遞函數(shù)為

式中：K2為1，T2為0.5，阻尼比ζ 為0.5。

故該三階系統(tǒng)的傳遞函數(shù)為

辨識結(jié)果如圖6 所示，根據(jù)辨識算法計(jì)算出的參數(shù)值與原電路模型參數(shù)值基本接近。

圖6 辨識結(jié)果（T1＜T2）Fig.6 Dentification results（T1＜T2）

T2＜T1＜5T2模擬電路各參數(shù)取值：R1=200 kΩ，C1=5 μF，R01=200 kΩ；R2=50 kΩ，C2=10 μF；R3=500 kΩ，C3=1 μF，R02=500 kΩ。

一階環(huán)節(jié)的傳遞函數(shù)為

式中：時(shí)間常數(shù)T1為1，K1為1。

二階環(huán)節(jié)的開環(huán)傳遞函數(shù)為

那么該二階環(huán)節(jié)的閉環(huán)傳遞函數(shù)為

式中：K2為2，T2為0.5，阻尼比ζ 為0.5。

故該三階系統(tǒng)的傳遞函數(shù)為

辨識結(jié)果如圖7 所示，根據(jù)辨識算法計(jì)算出的參數(shù)值與原電路模型參數(shù)值基本接近。

圖7 辨識結(jié)果（T2＜T1＜5T2）Fig.7 Dentification results（T2＜T1＜5T2）

5 結(jié)語

由辨識仿真結(jié)果看出，根據(jù)辨識算法計(jì)算出的參數(shù)值與原電路模型參數(shù)值基本接近，并且擬合曲線和原曲線十分接近，動(dòng)、靜態(tài)誤差也滿足辨識要求。辨識結(jié)果存在誤差主要有兩個(gè)原因：

（1）實(shí)驗(yàn)箱電阻、電容的標(biāo)識值與實(shí)際值可能有偏差：階躍響應(yīng)數(shù)據(jù)是基于實(shí)驗(yàn)箱模擬電路得到的，若電阻、電容的標(biāo)識值與實(shí)際值有偏差，那么從實(shí)驗(yàn)箱獲取的階躍響應(yīng)數(shù)據(jù)也就和事先設(shè)計(jì)的系統(tǒng)參數(shù)不同，根據(jù)這樣的階躍響應(yīng)數(shù)據(jù)辨識出來的參數(shù)也就必然與事先設(shè)計(jì)的系統(tǒng)不同。

（2）算法需深入研究：搜索逼近參數(shù)時(shí)設(shè)定的步長在算法中是固定的，會(huì)影響辨識精度和軟件搜索時(shí)間；另外，用于判斷終止搜索條件的動(dòng)態(tài)誤差計(jì)算方法也會(huì)影響辨識的結(jié)果。