一種高階線性時不變系統的辨識算法仿真

吳相甫

（航空工業中國飛機強度研究所，西安710065）

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

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

本文利用系統階躍響應的曲線特征來確定傳遞函數，思路是利用傳遞函數可以分解為一二階典型環節的并聯原理，假設模型的階數后，調節典型環節并聯模型中的參數值，用擬合曲線逼近實際的階躍響應曲線，從而得到系統的傳遞函數[9]。研究時借助Matlab 平臺開發辨識軟件，通過曲線擬合確定典型環節并聯模型的參數值，可以在一定精度內獲得系統的傳遞函數。

1 高階線性時不變系統模型分析

一般情況下，高階線性時不變系統的傳遞函數分子分母都是s 的多項式，可以寫為

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

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

由式（2）可知，任何高階線性時不變系統都可以分解為若干個一階系統和二階系統的典型環節之和。

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

在輸入為單位階躍函數時，輸出量可表示為

式中：q 為實數極點的個數；r 為共軛復數極點的對數。設0<ζk<1，將式（5）展成部分分式

設初始條件全部為零，將式（6）進行拉氏反變換，可得高階系統的單位階躍響應：

式（7）表明，高階系統的時間響應是由一階系統和二階系統的時間響應疊加組成的。

綜上所述，任何線性高階系統都可以分解為若干個一階和二階典型環節的并聯，且其響應也可以分解成若干個一階和二階典型環節的階躍響應之和。更進一步，通過一個一階系統和一個二階系統的典型環節的并聯模型的階躍響應曲線逼近實際系統的階躍響應曲線，通過調節參數

改變曲線形狀，不斷地逼近高階系統的階躍響應曲線，能夠反映原系統的整體特性[10]，即可用一個一階系統和一個二階系統典型環節的并聯模型去等效高階線性系統的傳遞函數。

2 階躍響應曲線特征分析

傳遞函數的參數值對系統的動態性能和穩態性能有著直接的影響，各參數的取值情況直接決定系統的階躍響應曲線特征。一階系統與二階系統典型環節并聯模型可表示為

各參數對階躍響應曲線形狀影響的規律如下：

增益值的影響系統階躍響應穩態值為增益值K1、K2之和，K1、K2對上升時間、峰值時間、調節時間無影響；K1、K2之和不變時，K2越大，峰值越大，即超調量越大。

時間常數的影響T1、T2由小變大時，峰值時間和調節時間變大，這是因為一階環節響應速度變慢；T2相對于T1會使響應的速度變慢的更快。按照T1、T2的比例關系將該結構的三階系統分為3 類：T1≥5T2、T2＜T1＜5T2、T1＜T2，時間常數對其階躍響應曲線形狀的影響如下[11]：

（1）T1≥5T2： 階躍響應曲線的第一個波峰值要小于穩態值。二階系統在階躍響應曲線前期產生作用，一階系統在后期產生作用。

（2）T2＜T1＜5T2： 階躍響應曲線的第一個波峰值大于穩態值。動態部分上升速度慢，超調小甚至沒有超調，但曲線上有拐點。

（3）T1＜T2： 階躍響應曲線的第一個波峰值大于穩態值。一階環節發生作用比二階環節早，一階環節在上升段初期發生作用，二階環節產生影響的時間靠后。

阻尼比的影響阻尼比ζ 增大時，超調減小，振蕩次數減少。峰值時間tp由其中的二階系統決定，根據二階系統峰值時間的計算方法，T2與ζ 滿足關系式

該式在K1、K2之和不變時，K1、K2的變化不影響該結論。

3 辨識方法

根據上一節的分析，階躍響應穩態值由K1、K2之和決定，T1、T2決定曲線的響應速度，阻尼比決定曲線的超調。按照T1、T2的比例關系對階躍響應曲線特征的影響將分3 類：T1≥5T2、T2＜T1＜5T2、T1＜T2，分析辨識方法。

T1≥5T2這種情況下，一階環節的時間常數比二階環節時間常數大得多，因此一階環節發生作用比二階環節晚。此時，階躍響應曲線的第一個波峰值要小于穩態值，該特征可用來判斷此系統屬于T1≥5T2這一類。

具體參數辨識算法如下：

（1）一階時間常數T1：由于一階環節發生作用要比二階環節要晚，曲線會先發生振蕩，然后曲線會平緩上升直到達到穩態值，因此可利用一階系統求時間常數的方法求出T1，即在階躍響應曲線上找到輸出量變化至終值95%時的坐標點，它所對應的時刻為調節時間ts，則T1=ts/3。

（2）ζ、T2、K1、K2：這4 個參數要由循環搜索逼近的方法確定。

算法包括兩個循環，外層循環調節ζ、T2，內層循環調節K1、K2：首先，求出峰值時間tp，令阻尼比ζ=0.05，由式（10）可以算出T2，并令K1=K2=K/2，在該組數值下求出階躍響應數據，算出第一個波峰差p=y1-y1′（實際階躍曲線值減去擬合曲線值），若p>0，則令K2=K2+0.1，K1=K1-0.1（若p<0，則K2=K2-0.1，K1=K1+0.1），此為內循環（用來調節K1、K2）；當p 的絕對值小于0.1 時，在當前參數值之下算出第一個波谷差bg=y2-y2′（實際階躍曲線值減去擬合曲線值），若bg>0.1，則說明阻尼比ζ 的值太小（振蕩太大），此時令ζ=ζ+0.05，重新計算T2，此為外循環（用來調節ζ、T2），并令K1=K2=K/2，再根據波峰差p 調節K1、K2。參數辨識流程見圖1。

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

T1＜T2這種情況下，一階環節的時間常數比二階環節時間常數小，因此一階環節發生作用比二階環節早。此時，階躍響應曲線的第一個波峰值大于穩態值，該特征可用來判斷此系統不屬于T1≥5T2類。

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

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

（1）ζ、T2： 這2 個參數由最外層循環搜索確定。由于含1 個一階環節，該環節會減緩系統的振蕩，因此由假設的二階系統求出的阻尼比ζ 是偏大的，因此要向下調整，T2則由式（10）計算得出。

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

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

參數辨識流程見圖2。

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

T2＜T1＜5T2這種情況與T1＜T2相同，階躍響應曲線的第一個波峰值大于穩態值，該特征可用來判斷此系統不屬于T1≥5T2類。建模思想同T1＜T2相同，首先認為待辨識系統為二階系統，并算出動態誤差d，當d<0.4 時，將待辨識系統歸為T2＜T1＜5T2類。

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

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

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

（3）當T1調整到使估計模型的峰值與穩態值大約相等時，說明T1已足夠大，此時若仍有動態誤差d>0.05、bg 的絕對值大于0.15 則進入第一步重新開始搜索。

參數辨識流程見圖3。

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

4 辨識仿真

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

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

該電路由1 個一階系統與1 個二階系統并聯而成，當電阻、電容參數值不相同時，這類結構的三階系統可分為3 種情況。

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

一階環節的傳遞函數為

式中：時間常數T1為1，K1為2。

二階環節的開環傳遞函數為

那么該二階環節的閉環傳遞函數為

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

故該三階系統的傳遞函數為

辨識結果如圖5 所示，根據辨識算法計算出的參數值與原電路模型參數值基本接近。

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

一階環節的傳遞函數為

圖5 辨識結果（T1≥5T2）Fig.5 Identification results（T1≥5T2）

式中：時間常數T1為0.1，K1為2。

二階環節的開環傳遞函數為

那么該二階環節的閉環傳遞函數為

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

故該三階系統的傳遞函數為

辨識結果如圖6 所示，根據辨識算法計算出的參數值與原電路模型參數值基本接近。

圖6 辨識結果（T1＜T2）Fig.6 Dentification results（T1＜T2）

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

一階環節的傳遞函數為

式中：時間常數T1為1，K1為1。

二階環節的開環傳遞函數為

那么該二階環節的閉環傳遞函數為

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

故該三階系統的傳遞函數為

辨識結果如圖7 所示，根據辨識算法計算出的參數值與原電路模型參數值基本接近。

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

5 結語

由辨識仿真結果看出，根據辨識算法計算出的參數值與原電路模型參數值基本接近，并且擬合曲線和原曲線十分接近，動、靜態誤差也滿足辨識要求。辨識結果存在誤差主要有兩個原因：

（1）實驗箱電阻、電容的標識值與實際值可能有偏差：階躍響應數據是基于實驗箱模擬電路得到的，若電阻、電容的標識值與實際值有偏差，那么從實驗箱獲取的階躍響應數據也就和事先設計的系統參數不同，根據這樣的階躍響應數據辨識出來的參數也就必然與事先設計的系統不同。

（2）算法需深入研究：搜索逼近參數時設定的步長在算法中是固定的，會影響辨識精度和軟件搜索時間；另外，用于判斷終止搜索條件的動態誤差計算方法也會影響辨識的結果。