基于差分進化算法的二階時滯系統參數辨識

李敏花,柏 猛

(山東科技大學電氣信息系,山東 濟南 250031)

0 引言

時滯現象廣泛存在于各類生產過程中。由于時滯環節對系統控制性能有重要影響,時滯現象在過程控制領域受到廣泛關注[1]。在實際應用中,通常將系統模型簡化為一階或二階時滯系統模型,以便進行控制器設計[2]。獲取時滯系統模型參數,特別是滯后時間參數,對控制系統的分析和設計具有重要意義[3]。目前,二階時滯系統參數辨識方法主要包括時域法[4]和頻域法[5]這2類。時域法通過對系統輸入輸出數據進行擬合,以實現參數估計。這類方法簡單、有效,適用范圍廣,但辨識結果易受觀測噪聲影響。頻域法對系統輸入特定的、滿足辨識要求的信號,根據系統頻域特性獲取系統參數。這類方法抗噪性強、算法復雜,對系統輸入/輸出數據要求高。系統階躍響應辨識系統模型是1種簡單、有效的時域辨識方法,廣泛應用于控制領域。目前,由階躍響應辨識系統傳遞函數的方法主要包括兩點法、相良節夫法、半對數法和面積法等。這些方法可以在一定程度上解決二階時滯系統的參數辨識問題。但總體而言,這些方法的辨識結果易受觀測噪聲影響,其參數估計精度有待提高。

對此,本文提出1種基于差分進化(differential evolution,DE)算法的二階時滯系統參數辨識方法。該方法將系統參數估計問題轉換為階躍響應數據的擬合問題,通過采用DE算法進行優化,可有效提高參數估計的準確性和魯棒性。

1 二階時滯系統參數估計方法

1.1 問題提出

二階時滯系統模型可表示為:

(1)

式中:a1和a2為模型參數;k為系統增益;τ為滯后時間。

根據系統阻尼比是否大于1,可將系統分為二階過阻尼時滯系統和二階欠阻尼時滯系統。

①二階過阻尼時滯系統。

典型的二階過阻尼時滯系統傳遞函數可表示為:

(2)

式中:T1和T2為時間常數。

式(2)與式(1)有參數關系,即a1=T1+T2、a2=T1T2。

零初始狀態下,二階過阻尼時滯系統單位階躍響應可表示為:

H(t-τ)

(3)

當T1=0或T2=0時,式(2)可轉化為一階時滯系統,即:

(4)

一階時滯系統的單位階躍響應可表示為:

(5)

②二階欠阻尼時滯系統。

二階欠阻尼時滯系統傳遞函數可表示為:

(6)

式中:ξ為阻尼比,0<ξ<1;wn為自然頻率。

在零初始狀態下,二階欠阻尼時滯系統單位階躍響應可表示為:

(7)

假設在ti時刻,二階時滯系統階躍響應y(ti)的觀測值為z(ti),即:

z(ti)=y(ti)+v(ti)

(8)

式中:v(ti)為觀測噪聲,i=1,2,…,m,m為數據采集個數。

(9)

式中:J(θ)為目標函數。

本文取J(θ)為:

(10)

顯然,目標函數J(θ)的極小化問題是典型的非線性最優化問題。經典求解方法,如Guass-Newton方法、Levenberg-Marquarat方法等,大多屬于局部優化方法[6]。這類方法需給定合適的參數初值,否則無法收斂到全局最優解。為獲得式(10)的全局最優解,本文采用DE算法對J(θ)進行優化。

1.2 DE算法

DE算法是由Storn和Price于1995年提出的1種智能進化算法。該算法采用實數編碼,通過引入變異、交叉和選擇這3種算子實現種群進化。研究結果表明,DE算法在求解非線性最優化問題時具有良好的全局收斂性[7]。同時,DE算法還具有算法簡單、易于實現、可靠性和魯棒性強的特點,非常適用于求解數值優化問題,因此被廣泛應用于各類優化問題的求解[8]。

針對式(10)所示的J(θ)優化問題,為便于描述,本文令X=0。對于式(4)所示的一階時滯系統,θ=[k,T,τ],即X=[x1,x2,x3]。對于式(2)所示的二階過阻尼和欠阻尼系統,θ分別取θ1=[k,T1,T2,τ]和θ2=[k′,ξ,wn,τ],即X=[x1,x2,x3,x4]。采用DE算法估計未知參數θ的主要步驟如下。

①種群初始化。

本文假設X的維數為D(對于二階時滯系統D=4),并且X取值的上限和下限分別為Xmax={x1,max,x2,max,…,xD,max}和Xmin={x1,min,x2,min,…,xD,min},則在Xmin和Xmax范圍內,隨機生成N個D維向量{X1,X2,…,XN}作為初始種群。對于任一向量Xi中的第j個元素xi,j,其生成方法如下。

xi,j=xj,min+rand[a,b]×(xj,max-xj,min)

(11)

式中:rand[a,b]為在區間[a,b]服從均勻分布的隨機數。

式(11)取a=0、b=1。

在上述過程中,N為種群規模(數量),一般取N∈[5D,10D]。綜合考慮算法效率和收斂性,本文取N=10D。

②變異操作。

對當前種群X1,X2,…,XN中的任一向量Xi產生相應的變異個體Vi:

Vi=Xr1+Fi(Xr2-Xr3)

(12)

式中:Xr1、Xr2和Xr3為在當前種群中隨機選擇的互不相同且與Xi不同的3個向量;Fi∈(0,1)為縮放因子。

③交叉操作。

將當前種群中的任一向量Xi與其變異個體Vi進行交叉,生成新的測試個體Ui。Ui中的第j個元素ui,j為:

(13)

式中:Ri,j=rand[0,1]為ui,j對應的隨機數;Ci∈[0,1]為Ui對應的交叉概率。

④選擇操作。

(14)

通過上述過程,生成新一代種群{X′1,X′2,…,X′N}繼續進行迭代,直到滿足迭代終止條件。

在DE算法中,縮放因子Fi和交叉概率Ci的設置對算法性能有重要影響,Fi和Ci的選擇與目標函數特性緊密相關。由于不同二階時滯系統所對應目標函數J(θ)的特性不同,為減小Fi和Ci對參數估計收斂性的影響,本文采用“either-or”策略[9]取代變異和交叉算子。采用“either-or”策略產生Ui的方法如下。

(15)

式中:pF∈[0,1]為變異概率;Ki為設定的參數,一般取Ki=0.5(Fi+1)。

研究結果表明[9],當Fi∈[0.5,1]、pF在大范圍內取值時,DE算法均可取得較好的優化結果。基于此,為簡化二階時滯系統參數辨識方法的參數選擇過程,本文取pF=0.5、Fi=rand[0.5,1]。

1.3 模型參數辨識方法

綜上所述,本文提出的采用DE算法估計二階時滯系統未知參數θ的主要步驟如下。

①判斷系統類型,確定未知參數。根據系統階躍響應是否有超調判斷系統類型,將二階過阻尼和欠阻尼系統未知參數分別設定為θ1=[k,T1,T2,τ]和θ2=[k′,ξ,wn,τ]。

②設定DE算法迭代次數最大值Gmax和較小的迭代終止參數ε>0,確定未知參數θ的取值范圍θmin和θmax,將式(10)作為目標函數,以產生初始種群。

③通過式(15)所示“either-or”策略產生試驗個體,并采用式(14)所示的選擇操作生成新一代種群。

④取當前種群中使目標函數最小的向量θbest作為解向量,判斷迭代終止條件:當J[θbest]小于ε或迭代次數G大于Gmax時,迭代終止。若不滿足迭代終止條件,則返回步驟③繼續迭代。

2 仿真結果與分析

為檢驗本文所提參數估計方法的性能,本文分別對二階過阻尼時滯系統、二階欠阻尼時滯系統和高階時滯系統進行參數估計試驗。試驗采用式(10)所示的J(θ)作為目標函數。為便于比較,試驗在采用目標函數J(θ)作為評價函數的同時采用歸一化均方根誤差(normalized root mean square error,NRMSE)進行比較。NRMSE的值f可表示為[10]:

(16)

顯然,參數估計結果越接近真值,則f值越大。通過比較準則函數值的大小,即可檢驗參數估計算法的性能。

2.1 無觀測噪聲的情況

本節主要檢驗本文方法在無觀測噪聲情況下的參數估計性能。試驗中,DE算法最大迭代次數取Gmax=200、種群規模為40、采樣周期Ts=0.01 s。試驗分別取基于繼電器反饋(relay feedback,RF)的辨識方法[11]、頻域辨識(frequency identification,FI)方法[12]和本文方法的參數估計結果進行比較。

①二階過阻尼時滯系統。

試驗采用的二階過阻尼時滯系統傳遞函數為:

(17)

系統待估計參數分別為a2=10、a1=11、k=1、τ=2。試驗采用本文方法進行參數估計時,以式(3)作為擬合模型,0≤t≤70s,參數向量取θ=[k,T1,T2,τ],且有a1=T1+T2、a2=T1T2,參數取值范圍為θmin=[0,0,0,0]和θmax=[5,15,15,5]。采用不同方法得到的G1參數估計結果如表1所示。

表1 G1參數估計結果

由表1可知:本文方法在估計二階過阻尼時滯系統參數時具有很高的估計精度;FI方法估計效果較好;RF方法估計效果較差。G1參數估計收斂特性如圖1所示。

圖1 G1參數估計收斂特性

由圖1可知,采用本文方法進行參數估計時迭代次數少于50次即可達到最優值,收斂速度較快。由此可見,本文方法可有效解決二階過阻尼時滯系統的參數估計問題,具有較高的估計精度。

②二階欠阻尼時滯系統。

試驗采用的二階欠阻尼時滯系統的傳遞函數為:

(18)

為便于與其他算法估計結果進行比較,本文將式(18)寫成以下形式:

(19)

系統未知參數分別為a2=1、a1=1.5、k=1、τ=0.5。試驗采用本文算法進行參數估計時,系統待估計參數向量取θ=[k′,ξ,wn,τ]。試驗采用式(7)作為擬合模型,0≤t≤10s,參數取值范圍為θmin=[0,0,0,0]和θmax=[2,1,5,2]。采用不同方法得到的G2參數估計結果如表2所示。

表2 G2參數估計結果

由表2可知,3種參數估計方法都可有效估計二階欠阻尼時滯系統的未知參數。相比較而言,本文方法的參數估計精度優于其他2種方法。G2參數估計收斂特性如圖2所示。

圖2 G2參數估計收斂特性

由圖2可知,本文方法具有較快的收斂速度,迭代次數少于50次即可達到最優值。由此可見,本文方法可有效解決二階欠阻尼時滯系統的參數估計問題。

③高階系統近似模型參數估計。

在實際應用中,常將高階系統近似為二階系統以便設計控制器。高階系統近似模型的參數估計是設計控制系統時需解決的首要問題。試驗對以下高階系統進行近似:

(20)

由試驗系統階躍響應可知,該系統為過阻尼系統。因此,試驗采用式(2)所示二階過阻尼時滯系統模型對上述系統進行近似。在本文方法進行參數估計時,本文采用式(3)作為擬合模型,0≤t≤10s,參數向量取θ=[k,T1,T2,τ],參數范圍取θmin=[0,0,0,0]和θmax=[2,5,5,5]。采用不同方法得到的G3參數估計結果如表3所示。

表3 G3參數估計結果

由表3可知,針對試驗高階系統的近似問題:RF方法辨識得到的近似系統效果較差;FI方法可獲得較好的近似效果;本文方法辨識得到的近似系統與實際系統輸出差距最小。G3參數估計收斂特性如圖3所示。

圖3 G3參數估計收斂特性

由圖3可知,本文方法具有較快的收斂性,迭代次數在50次內即可達到最優值。由此可見,本文方法可有效解決高階系統的低階近似問題。

2.2 含有觀測噪聲的情況

為檢驗本文方法在不同噪聲條件下的參數估計性能,本節采用式(17)所示二階過阻尼時滯系統進行參數估計試驗。試驗中,系統模型在每種噪聲強度下均進行30次參數估計試驗,取估計結果的均值和均方差作為參數估計結果。系統觀測噪聲取高斯白噪聲,輸出信號噪信比(noise-to-signal ratio,NSR)定義為:

(21)

式中:{y(ti)}和{v(ti)}分別為系統階躍響應輸出信號和觀測噪聲的均方差。

試驗中,DE算法的種群規模為40、最大迭代次數為tmax=400、采樣周期Ts=0.01 s,其他參數關系及設置均與無噪聲試驗時相同。不同觀測噪聲情況下的參數估計結果如表4所示。

表4 不同觀測噪聲情況下的參數估計結果

由表4可知,當系統階躍響應含有觀察噪聲時,本文方法依然可以獲得較好的參數估計結果。其中,參數估計評價指標f值與NSR近似成線性關系,即f≈100-10×R。同時,由試驗結果可以看出,當NSR較小(δ<0.2)時,參數辨識精度較高,方差較小;而當NSR較大(δ>0.2)時,參數估計結果誤差較大。由此可知,當系統階躍響應數據無觀測噪聲或噪聲較小時,采用本文方法可獲得滿意的辨識結果;而當觀察噪聲較大時,參數辨識結果雖有一定偏差,但依然可為控制系統設計提供模型參考。

3 結論

本文提出的基于DE算法的二階時滯系統參數辨識方法是1種根據系統階躍響應估計系統未知參數的時域辨識方法。試驗結果表明,該方法將參數估計問題轉換為系統階躍響應數據擬合問題,采用DE算法可有效解決二階時滯系統的參數估計問題。與其他方法相比,本文方法只需根據系統階躍響應即可辨識二階系統未知參數,對系統輸入輸出數據要求較低,具有算法簡單、可靠性高、適用性強和參數估計精度高的特點,非常適合工業應用。后續工作將主要研究多變量時滯系統的參數辨識問題和強觀測噪聲條件下低階時滯系統的參數辨識問題。