李敏花,柏 猛

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

0 引言

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

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

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

1.1 問(wèn)題提出

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

(1)

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

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

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

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

(2)

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

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

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

H(t-τ)

(3)

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

(4)

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

(5)

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

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

(6)

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

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

(7)

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

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

(8)

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

(9)

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

本文取J(θ)為:

(10)

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

1.2 DE算法

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

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

①種群初始化。

本文假設(shè)X的維數(shù)為D(對(duì)于二階時(shí)滯系統(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個(gè)D維向量{X1,X2,…,XN}作為初始種群。對(duì)于任一向量Xi中的第j個(gè)元素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。

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

②變異操作。

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

Vi=Xr1+Fi(Xr2-Xr3)

(12)

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

③交叉操作。

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

(13)

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

④選擇操作。

(14)

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

在DE算法中,縮放因子Fi和交叉概率Ci的設(shè)置對(duì)算法性能有重要影響,Fi和Ci的選擇與目標(biāo)函數(shù)特性緊密相關(guān)。由于不同二階時(shí)滯系統(tǒng)所對(duì)應(yīng)目標(biāo)函數(shù)J(θ)的特性不同,為減小Fi和Ci對(duì)參數(shù)估計(jì)收斂性的影響,本文采用“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)取值時(shí),DE算法均可取得較好的優(yōu)化結(jié)果。基于此,為簡(jiǎn)化二階時(shí)滯系統(tǒng)參數(shù)辨識(shí)方法的參數(shù)選擇過(guò)程,本文取pF=0.5、Fi=rand[0.5,1]。

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

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

①判斷系統(tǒng)類(lèi)型,確定未知參數(shù)。根據(jù)系統(tǒng)階躍響應(yīng)是否有超調(diào)判斷系統(tǒng)類(lèi)型,將二階過(guò)阻尼和欠阻尼系統(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)生初始種群。

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

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

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

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

(16)

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

2.1 無(wú)觀測(cè)噪聲的情況

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

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

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

(17)

系統(tǒng)待估計(jì)參數(shù)分別為a2=10、a1=11、k=1、τ=2。試驗(yàn)采用本文方法進(jìn)行參數(shù)估計(jì)時(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ù)估計(jì)結(jié)果如表1所示。

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

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

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

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

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

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

(18)

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

(19)

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

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

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

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

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

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

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

(20)

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

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

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

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

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

2.2 含有觀測(cè)噪聲的情況

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

(21)

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

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

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

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

3 結(jié)論

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