（上海理工大學(xué)光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093）

0 引言

工業(yè)生產(chǎn)過(guò)程中，要想良好的控制被控對(duì)象，首先要確定被控對(duì)象模型，但工業(yè)過(guò)程如伺服電機(jī)控制系統(tǒng)、爐內(nèi)脫硫系統(tǒng)等往往會(huì)產(chǎn)生一些非線性、復(fù)雜、高階的［1-3］時(shí)滯對(duì)象。而工業(yè)控制中最重要的控制器就是PID 控制器［4］，但對(duì)于高階時(shí)滯系統(tǒng)的PID 參數(shù)設(shè)計(jì)非常繁瑣。為了方便控制器設(shè)計(jì)，需以低階簡(jiǎn)易的模型進(jìn)行擬合。系統(tǒng)辨識(shí)是工業(yè)建模中的常用方法，一般先選取合適的傳遞函數(shù)模型，然后加上激勵(lì)信號(hào)，通過(guò)分析系統(tǒng)輸入輸出特性，采用一定的辨識(shí)方法得到對(duì)象的等價(jià)模型［5］。

傳統(tǒng)的辨識(shí)方法主要有極大似然法、頻域響應(yīng)法、最小二乘法。工業(yè)生產(chǎn)對(duì)于辨識(shí)高階系統(tǒng)所采用的模型，主要以基于最小二乘法的一階時(shí)滯模型［6-8］和基于1/1Pade逼近的二階模型［9］為主。由于對(duì)象的復(fù)雜性，實(shí)際應(yīng)用這些方法會(huì)與原系統(tǒng)存在較大誤差，并且在擬合精度上較差。為了提高模型擬合精度，Wang 等［10］提出了一種通過(guò)對(duì)信號(hào)進(jìn)行分解，選取10 個(gè)頻率點(diǎn)來(lái)擬合二階純滯后模型的閉環(huán)辨識(shí)方法，該方法對(duì)于一些強(qiáng)耦合、大時(shí)滯的復(fù)雜對(duì)象具有很好的辨識(shí)效果；文獻(xiàn)［11］采用二階時(shí)滯模型對(duì)含積分環(huán)節(jié)的對(duì)象進(jìn)行辨識(shí)；文獻(xiàn)［12］利用頻域分析及選點(diǎn)，通過(guò)二階時(shí)滯模型對(duì)含微分環(huán)節(jié)的對(duì)象進(jìn)行辨識(shí)；文獻(xiàn)［13］通過(guò)二階時(shí)滯模型對(duì)具有振蕩特性或者非振蕩特性的高階時(shí)滯模型進(jìn)行辨識(shí)。這些學(xué)者的實(shí)驗(yàn)均驗(yàn)證了二階時(shí)滯模型的優(yōu)越性，但他們的研究均是在重要頻率段內(nèi)隨機(jī)選取一定數(shù)量的頻率響應(yīng)點(diǎn)，由于選點(diǎn)的隨機(jī)性及選點(diǎn)個(gè)數(shù)的局限性，模型參數(shù)的精度仍有很大提升空間。

本文提出采用蟻群算法實(shí)現(xiàn)高階時(shí)滯系統(tǒng)辨識(shí)優(yōu)化，是在二階時(shí)滯模型基礎(chǔ)上，充分利用蟻群算法的特點(diǎn)，通過(guò)尋優(yōu)適應(yīng)度函數(shù)極小值，找到在重要頻率段內(nèi)的最佳頻率響應(yīng)點(diǎn)，即通過(guò)較少的頻率響應(yīng)點(diǎn)，求出二階時(shí)滯模型的各參數(shù)，利用MATLAB 仿真并與傳統(tǒng)工業(yè)方法作比較。通過(guò)結(jié)果可知，本文方法得出的幅頻特性曲線與原曲線擬合程度更高。

1 對(duì)象的頻域建模方法

1.1 頻率響應(yīng)特性

控制系統(tǒng)中的信號(hào)經(jīng)常出現(xiàn)從一個(gè)穩(wěn)定狀態(tài)變化到另一個(gè)穩(wěn)定狀態(tài)的情況，兩個(gè)穩(wěn)定狀態(tài)變化過(guò)程中包含了極為重要的頻域特性［14］。一般的系統(tǒng)響應(yīng)信號(hào)f（t）是由瞬態(tài)部分Δf(t)和穩(wěn)態(tài)部分fs(t)組成，如式（1）所示。

由于兩個(gè)穩(wěn)定狀態(tài)之間的變化信號(hào)通常不滿足絕對(duì)可積條件，考慮對(duì)等式兩邊進(jìn)行Laplace 變換，得到式（2）。

若系統(tǒng)輸出信號(hào)在t=Tf時(shí)刻達(dá)到穩(wěn)定狀態(tài)，此時(shí)瞬態(tài)響應(yīng)Δf(t)的值為0，由式（2）可得式（3）。

控制系統(tǒng)的輸入信號(hào)和輸出信號(hào)可由式（3）近似得出當(dāng)前工作點(diǎn)的傳遞函數(shù)。

若系統(tǒng)對(duì)象輸入為階躍信號(hào)時(shí)，?u(t)值為0，輸出則為瞬態(tài)和穩(wěn)態(tài)兩個(gè)部分。令s=jω，由式（4）可得式（5）。

將式（5）的積分部分分成N 個(gè)長(zhǎng)度為Δtk的子區(qū)間，展開(kāi)成積分累加的形式，如式（6）所示。

其中，Nf為數(shù)據(jù)的采樣個(gè)數(shù)。

在式（6）的基礎(chǔ)上，進(jìn)一步推導(dǎo)整理，得到對(duì)象各重要頻率點(diǎn)的頻率響應(yīng)特性。

這樣對(duì)于任意給定頻率ω，都可根據(jù)式（7）求出該頻率點(diǎn)時(shí)對(duì)象的幅值和相位。

1.2 重要頻率段確定

假設(shè)對(duì)象模型收斂，由于系統(tǒng)穩(wěn)定性主要由中低頻段決定，而高頻段幾乎不影響系統(tǒng)的穩(wěn)定，所以重要頻率段位于Nyquist 圖的第二、三象限，即從零開(kāi)始到對(duì)象的臨界頻率點(diǎn)是其重要的頻域范圍。構(gòu)建弦截迭代公式如式（8）所示［15］。

相角與頻率間的函數(shù)關(guān)系如式（9）所示。

初值ω0和φ0均為零，ω1取一個(gè)極小的數(shù)，例如0.000 01；其中G(jωn)可由式（7）得到，φ1通過(guò)式（9）計(jì)算得到。經(jīng)過(guò)數(shù)次的反復(fù)迭代后，可以求出對(duì)象的臨界頻率ωc（對(duì)應(yīng)相角為-180°），則（0，ωc）為確定的重要頻率段。

1.3 傳遞函數(shù)模型

實(shí)際生產(chǎn)過(guò)程中的工業(yè)控制非常復(fù)雜，文中采用具有任意極點(diǎn)的二階時(shí)滯模型去擬合高階時(shí)滯對(duì)象，傳遞函數(shù)模型如式（10）所示。

其參數(shù)是在重要頻率段內(nèi)選取多個(gè)頻率響應(yīng)點(diǎn)，本文利用蟻群算法，在重要頻率段內(nèi)通過(guò)選取5 個(gè)較優(yōu)點(diǎn)求取，再通過(guò)幅值和相角關(guān)系擬合此模型。幅值關(guān)系如式（11）所示。

可用矩陣表示如式（12）所示。

其中：

根據(jù)線性最小二乘法式（13）：

從而得到式（14）。

可求得參數(shù)a、b、c，若對(duì)象的輸出隨著系統(tǒng)輸入增加而增加，則參數(shù)值取正，反之取負(fù)。相角關(guān)系如式（15）所示。

同樣采用最小二乘法求出參數(shù)L 如式（16）所示。

2 蟻群算法

2.1 算法介紹

蟻群算法最早由Dorigo 等［16］于1991 年提出，他們?cè)谘芯啃滦退惴ǖ倪^(guò)程中，發(fā)現(xiàn)蟻群在尋找食物時(shí)，通過(guò)分泌一種稱為信息素的生物激素交流覓食信息從而能快速找到目標(biāo)，據(jù)此提出了基于信息正反饋原理的蟻群算法。

蟻群算法主要應(yīng)用于求解01 背包問(wèn)題、最優(yōu)解問(wèn)題、TSP 問(wèn)題等，目前已逐漸應(yīng)用于其它領(lǐng)域。作為啟發(fā)式優(yōu)化算法，其以較佳的全局搜索能力、較快的收斂速度而廣泛應(yīng)用于各類尋優(yōu)案例中。其優(yōu)點(diǎn)：①采用正反饋機(jī)制，使得搜索過(guò)程不斷收斂，最終逼近最優(yōu)解；②每個(gè)個(gè)體可以通過(guò)釋放信息素來(lái)改變周圍環(huán)境，且每個(gè)個(gè)體能夠感知周圍環(huán)境的實(shí)時(shí)變化，個(gè)體間通過(guò)環(huán)境進(jìn)行間接通訊；③搜索過(guò)程采用分布式計(jì)算方式，多個(gè)個(gè)體同時(shí)進(jìn)行并行計(jì)算，極大提高了算法計(jì)算能力和運(yùn)行效率；④啟發(fā)式的概率搜索方式不容易陷入局部最優(yōu)，易于尋找到全局最優(yōu)解。

其算法步驟如下：

步驟1：初始化參數(shù)，包含蟻群規(guī)模N*D，信息素?fù)]發(fā)程度因子Rho，信息素常數(shù)Q，轉(zhuǎn)移概率常數(shù)p0，最大迭代次數(shù)iter_max，信息素因子Alpha，啟發(fā)函數(shù)因子Belta。

步驟2：構(gòu)建解空間，將各螞蟻隨機(jī)置于不同位置，對(duì)每個(gè)螞蟻，按照轉(zhuǎn)移概率計(jì)算公式，確定其下一個(gè)位置。

步驟3：更新信息素，計(jì)算各螞蟻所在位置的信息素含量，根據(jù)信息素迭代公式對(duì)各位置上的信息素濃度進(jìn)行更新，同時(shí)記錄當(dāng)前迭代次數(shù)的最優(yōu)解。

步驟四：判斷是否終止，若達(dá)到最大迭代次數(shù)，則終止計(jì)算，輸出最優(yōu)解，否則，返回步驟2。

利用蟻群算法進(jìn)行優(yōu)化的參數(shù)為重要頻率段中的D個(gè)頻率點(diǎn)，即算法隨機(jī)生成N*D 維種群，代入到設(shè)定的適應(yīng)度函數(shù)中求解出適應(yīng)度值，再通過(guò)不斷更新種群位置，迭代求解出適應(yīng)度全局極值以及對(duì)應(yīng)的最佳D 個(gè)頻率點(diǎn)。

2.2 適應(yīng)度函數(shù)

適應(yīng)度函數(shù)的好壞直接影響著智能優(yōu)化算法性能好壞［17］，時(shí)間乘以誤差絕對(duì)值積分是工程中常用的誤差積分準(zhǔn)則之一，由于其調(diào)節(jié)時(shí)間短，并且系統(tǒng)參數(shù)的細(xì)微改變可以使得指標(biāo)有明顯變化，能夠很好地反映控制系統(tǒng)性能，因此具有良好的選擇性和實(shí)用性［18］。ITAE 是時(shí)間乘以誤差絕對(duì)值積分的性能指標(biāo)，本文中的ITAE 值代表著曲線的擬合度，值越小代表擬合程度越高，以ITAE 為參考指標(biāo)對(duì)求得的參數(shù)進(jìn)行分析，ITAE 表達(dá)式如式（17）所示。

采用累加和形式，在采樣時(shí)間（0，T）內(nèi)將積分區(qū)間等分成n 個(gè)子區(qū)間，可近似表示成如式（18）所示。

式（18）中，y1 表示原始系統(tǒng)的階躍響應(yīng)，y2 表示辨識(shí)后模型的階躍響應(yīng)，δi為[ti-1,ti](t0=0,tn=T)中的某個(gè)數(shù)，Δti為一個(gè)極小值。

3 仿真實(shí)驗(yàn)

為了驗(yàn)證通過(guò)本文改進(jìn)方法得到的辨識(shí)參數(shù)具有更好的準(zhǔn)確性和穩(wěn)定性，考慮文獻(xiàn)［19］給出的串級(jí)系統(tǒng)外環(huán)高階時(shí)滯對(duì)象如式（19）所示。

假定輸入信號(hào)u（t）為階躍信號(hào)，通過(guò)上述提到的頻率辨識(shí)方法，可求出G（s）的穿越頻率為0.030 8，即蟻群算法中個(gè)體位置范圍為［0，0.030 8］。由于蟻群算法受其參數(shù)選取的嚴(yán)重影響，文獻(xiàn)［20-21］對(duì)蟻群算法的參數(shù)進(jìn)行了優(yōu)化，并對(duì)比了多組實(shí)驗(yàn)，得出最佳參數(shù)組合的取值范圍。依據(jù)此文獻(xiàn)，對(duì)各參數(shù)進(jìn)行選值：蟻群規(guī)模N 為50，D 為5，Rho 為0.7，Q 為100，p0 為0.2，Alpha 為2，Belta 為5，最大迭代次數(shù)為100。依次運(yùn)行5 次后，每次得到5 個(gè)最優(yōu)頻率點(diǎn)，計(jì)算出其平均值和方差如表1 所示。

Table 1 Algorithm data表1 算法數(shù)據(jù)

由表1 方差可得，該算法求出的最優(yōu)點(diǎn)比較穩(wěn)定，選取表1 中5 個(gè)點(diǎn)的平均值，計(jì)算出二階時(shí)滯模型各參數(shù)如表2 所示。

Table 2 Model parameter表2 模型參數(shù)

對(duì)于高階系統(tǒng)，工業(yè)中普遍采用一階時(shí)滯模型進(jìn)行建模，利用最小二乘法計(jì)算模型參數(shù)［22］，也常使用基于1/1Pade 逼近的二階模型擬合法［9］。表3 分別為采用傳統(tǒng)一階模型辨識(shí)法、1/1Pade 逼近法、傳統(tǒng)二階加時(shí)滯模型法及由本文方法所得出的傳遞函數(shù)模型。求出ITAE 指標(biāo)，本文的ITAE 值代表曲線擬合度，ITAE 值越小擬合程度越高。

Table 3 Results of four identification methods表3 4 種辨識(shí)方法結(jié)果

圖1 為上述4 種辨識(shí)方法所得出的Nyquist 圖，其中黑色實(shí)線是模型原型，藍(lán)色虛線為傳統(tǒng)一階方法，綠色實(shí)線為1/1Pade 逼近法，紅色實(shí)線為傳統(tǒng)二階加時(shí)滯方法，紅色虛線是本文方法，圖2 為其部分放大圖（彩圖掃OSID 碼可見(jiàn)）。

Fig.1 Nyquist diagram圖1 Nyquist 圖

Fig.2 Partial enlarged view of Nyquist圖2 部分Nyquist 放大圖

同樣，考慮文獻(xiàn)［19］中給出的串級(jí)系統(tǒng)內(nèi)環(huán)高階對(duì)象如式（20）所示。

采用上述同樣的方法，可求出G（s）的穿越頻率為0.025，即蟻群算法中個(gè)體位置范圍為［0，0.025］。分別為采用傳統(tǒng)一階模型辨識(shí)法、1/1Pade 逼近法、傳統(tǒng)二階加時(shí)滯模型法及本文方法，所得出的傳遞函數(shù)模型如表4 所示，并在此基礎(chǔ)上求出ITAE 性能指標(biāo)。

Table 4 Results of four identification methods表4 4 種辨識(shí)方法結(jié)果

圖3 為上述4 種辨識(shí)方法所得出的Nyquist 圖，其中黑色實(shí)線是模型原型，藍(lán)色虛線為傳統(tǒng)一階方法，綠色實(shí)線為1/1Pade 逼近法，紅色實(shí)線為傳統(tǒng)二階加時(shí)滯方法，紅色虛線是本文方法，圖4 為其部分放大圖（彩圖掃OSID 碼可見(jiàn)）。

Fig.3 Nyquist diagram圖3 Nyquist 圖

Fig.4 Partial enlarged view of Nyquist圖4 部分Nyquist 放大圖

由圖1—圖4 的Nyquist 圖及ITAE 指標(biāo)可知，本文基于二階加時(shí)滯模型的改進(jìn)方法與傳統(tǒng)一階方法、Pade 逼近法、二階時(shí)滯方法相比，與實(shí)際對(duì)象模型的擬合精度更高，在某些頻率上能夠達(dá)到完全擬合的效果，并且有效避免了傳統(tǒng)二階時(shí)滯方法中特征頻率選取點(diǎn)的隨機(jī)性對(duì)最終結(jié)果的影響，可以通過(guò)使用較少的頻率點(diǎn)獲得更優(yōu)的模型參數(shù)。

4 結(jié)語(yǔ)

本文將蟻群算法運(yùn)用在高階時(shí)滯系統(tǒng)的辨識(shí)中，是對(duì)普通二階時(shí)滯模型方法的一種改進(jìn)，根據(jù)系統(tǒng)正常運(yùn)行時(shí)產(chǎn)生的數(shù)據(jù)，求取對(duì)象的重要頻率段，再通過(guò)蟻群算法，獲取較少的特征頻率點(diǎn)，最終通過(guò)幅頻特性確定最終模型，并使得最終模型的參數(shù)更接近最優(yōu)值。Matlab 仿真結(jié)果表明，該方法對(duì)高階時(shí)滯對(duì)象具有更好的辨識(shí)效果，與傳統(tǒng)的一階方法和二階時(shí)滯方法相比，具有更高的模型精度，可應(yīng)用于實(shí)際工業(yè)生產(chǎn)模型辨識(shí)。