基于差分進(jìn)化算法的二階時滯系統(tǒng)參數(shù)辨識

李敏花,柏 猛

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

0 引言

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

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

1 二階時滯系統(tǒng)參數(shù)估計方法

1.1 問題提出

二階時滯系統(tǒng)模型可表示為:

(1)

式中:a1和a2為模型參數(shù);k為系統(tǒng)增益;τ為滯后時間。

根據(jù)系統(tǒng)阻尼比是否大于1,可將系統(tǒng)分為二階過阻尼時滯系統(tǒng)和二階欠阻尼時滯系統(tǒng)。

①二階過阻尼時滯系統(tǒng)。

典型的二階過阻尼時滯系統(tǒng)傳遞函數(shù)可表示為:

(2)

式中:T1和T2為時間常數(shù)。

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

零初始狀態(tài)下,二階過阻尼時滯系統(tǒng)單位階躍響應(yīng)可表示為:

H(t-τ)

(3)

當(dāng)T1=0或T2=0時,式(2)可轉(zhuǎn)化為一階時滯系統(tǒng),即:

(4)

一階時滯系統(tǒng)的單位階躍響應(yīng)可表示為:

(5)

②二階欠阻尼時滯系統(tǒng)。

二階欠阻尼時滯系統(tǒng)傳遞函數(shù)可表示為:

(6)

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

在零初始狀態(tài)下,二階欠阻尼時滯系統(tǒng)單位階躍響應(yīng)可表示為:

(7)

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

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

(8)

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

(9)

式中:J(θ)為目標(biāo)函數(shù)。

本文取J(θ)為:

(10)

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

1.2 DE算法

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

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

①種群初始化。

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

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

(11)

式中:rand[a,b]為在區(qū)間[a,b]服從均勻分布的隨機(jī)數(shù)。

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

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

②變異操作。

對當(dāng)前種群X1,X2,…,XN中的任一向量Xi產(chǎn)生相應(yīng)的變異個體Vi:

Vi=Xr1+Fi(Xr2-Xr3)

(12)

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

③交叉操作。

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

(13)

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

④選擇操作。

(14)

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

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

(15)

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

研究結(jié)果表明[9],當(dāng)Fi∈[0.5,1]、pF在大范圍內(nèi)取值時,DE算法均可取得較好的優(yōu)化結(jié)果?；诖?為簡化二階時滯系統(tǒng)參數(shù)辨識方法的參數(shù)選擇過程,本文取pF=0.5、Fi=rand[0.5,1]。

1.3 模型參數(shù)辨識方法

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

①判斷系統(tǒng)類型,確定未知參數(shù)。根據(jù)系統(tǒng)階躍響應(yīng)是否有超調(diào)判斷系統(tǒng)類型,將二階過阻尼和欠阻尼系統(tǒng)未知參數(shù)分別設(shè)定為θ1=[k,T1,T2,τ]和θ2=[k′,ξ,wn,τ]。

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

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

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

2 仿真結(jié)果與分析

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

(16)

顯然,參數(shù)估計結(jié)果越接近真值,則f值越大。通過比較準(zhǔn)則函數(shù)值的大小,即可檢驗(yàn)參數(shù)估計算法的性能。

2.1 無觀測噪聲的情況

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

①二階過阻尼時滯系統(tǒng)。

試驗(yàn)采用的二階過阻尼時滯系統(tǒng)傳遞函數(shù)為:

(17)

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

表1 G1參數(shù)估計結(jié)果

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

圖1 G1參數(shù)估計收斂特性

由圖1可知,采用本文方法進(jìn)行參數(shù)估計時迭代次數(shù)少于50次即可達(dá)到最優(yōu)值,收斂速度較快。由此可見,本文方法可有效解決二階過阻尼時滯系統(tǒng)的參數(shù)估計問題,具有較高的估計精度。

②二階欠阻尼時滯系統(tǒng)。

試驗(yàn)采用的二階欠阻尼時滯系統(tǒng)的傳遞函數(shù)為:

(18)

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

(19)

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

表2 G2參數(shù)估計結(jié)果

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

圖2 G2參數(shù)估計收斂特性

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

③高階系統(tǒng)近似模型參數(shù)估計。

在實(shí)際應(yīng)用中,常將高階系統(tǒng)近似為二階系統(tǒng)以便設(shè)計控制器。高階系統(tǒng)近似模型的參數(shù)估計是設(shè)計控制系統(tǒng)時需解決的首要問題。試驗(yàn)對以下高階系統(tǒng)進(jìn)行近似:

(20)

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

表3 G3參數(shù)估計結(jié)果

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

圖3 G3參數(shù)估計收斂特性

由圖3可知,本文方法具有較快的收斂性,迭代次數(shù)在50次內(nèi)即可達(dá)到最優(yōu)值。由此可見,本文方法可有效解決高階系統(tǒng)的低階近似問題。

2.2 含有觀測噪聲的情況

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

(21)

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

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

表4 不同觀測噪聲情況下的參數(shù)估計結(jié)果

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

3 結(jié)論

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