基于巴氏系數(shù)的非高斯隨機(jī)系統(tǒng)性能評估

黃國豆,張金芳,黃俊雄,翟宇卓

(華北電力大學(xué)控制與計算機(jī)工程學(xué)院,北京 102206)

1 引言

控制回路在自動化系統(tǒng)中起著最重要的作用。性能評估的意義在于實現(xiàn)、恢復(fù)和保持控制回路的最佳性能。因此,在工業(yè)運行中快速、準(zhǔn)確地評估系統(tǒng)的性能具有重要的應(yīng)用價值[1]。

實際工業(yè)過程中的噪聲擾動大都不符合高斯分布,低階統(tǒng)計量(均方誤差)所包含的信息顯然無法完全反映變量的高階統(tǒng)計特性。對于非高斯隨機(jī)系統(tǒng),應(yīng)該在均值和方差之外選擇一個更合適的指標(biāo)來反映變量的高階統(tǒng)計特性[2]。信息熵作為一種不確定性度量,相較于均值或方差,對于任何隨機(jī)變量都有更為一般的含義。所有高階矩(不僅僅是二階)的最小化都可以通過最小化熵來實現(xiàn)[3,4]。因此,近年來一種基于熵的控制系統(tǒng)性能評估方法發(fā)展迅速。

鐘振芳[5]等通過計算最小熵控制下的概率密度函數(shù)(Probability Density Function,PDF),以此為基準(zhǔn)構(gòu)建了系統(tǒng)的性能評估指標(biāo)。Zhang[6]使用最小Renyi熵作為評估指標(biāo),將閉環(huán)控制下的輸出分解為只與擾動相關(guān)的反饋不變量和反饋變量。Zhou[4]提出了基于有理熵的系統(tǒng)性能評價指標(biāo),并給出了該指標(biāo)理論基準(zhǔn)值和估計基準(zhǔn)值的計算方法。在文獻(xiàn)[7,8]中,熵指標(biāo)被用于評估典型串級控制系統(tǒng)的性能。然而,熵指標(biāo)具有均值漂移的缺點,且熵基準(zhǔn)計算過程十分復(fù)雜,所耗時間較長,因此缺乏實際工程意義。

針對以上問題,本文提出了一種基于巴氏系數(shù)的非高斯隨機(jī)系統(tǒng)控制回路性能評估指標(biāo)。在統(tǒng)計學(xué)中,巴氏距離(Bhattacharyya distance)多被用來測量兩個離散或連續(xù)概率分布之間的相似性,常用于性能監(jiān)測[9],圖像處理[10],神經(jīng)網(wǎng)絡(luò)[11]等領(lǐng)域。由于該指標(biāo)具有無界性,因此并不適合直接用作系統(tǒng)性能評價指標(biāo),但與其密切相關(guān)的巴氏系數(shù)卻是對稱且有界的,在反映兩個分布之間的相似性方面并無太大差異,且計算更為簡便。通過計算系統(tǒng)基準(zhǔn)輸出與實際輸出之間的巴氏系數(shù)就可以得到系統(tǒng)的性能指標(biāo)。

為了更快速準(zhǔn)確地估計性能基準(zhǔn),結(jié)合摸石頭過河算法(Wading Across Stream Algorithm,WSA)的優(yōu)點,提出了一種改進(jìn)的混合分布估計算法(Hybrid Estimation of Distribution Algorithm,HEDA)進(jìn)行系統(tǒng)參數(shù)尋優(yōu)和噪聲PDF估計。通過對高斯和非高斯噪聲擾動下的單回路反饋控制系統(tǒng)的仿真,驗證了改進(jìn)算法和新指標(biāo)的有效性。

2 性能評估指標(biāo)

2.1 巴氏系數(shù)

巴氏系數(shù)是衡量兩個統(tǒng)計樣本之間重疊程度的度量,與巴氏距離密切相關(guān)。同時,它還可以用來確定兩個樣本是否相近及可分離性[11]。

對于定義在同一域X上的概率分布P和Q,巴氏距離的定義為

DB(P,Q)=-ln(BC(P,Q))

(1)

其中,BC(P,Q)即為概率分布P和Q的巴氏系數(shù)。對于離散的概率分布,P和Q的巴氏系數(shù)定義如下

(2)

且滿足條件0≤BC≤1,0≤D≤∞。

對于連續(xù)概率分布,P和Q的巴氏系數(shù)定義為

(3)

兩個概率分布之間的巴氏系數(shù)值會隨著兩個樣本相同部分的增多而變大。對于任意兩個概率分布,巴氏系數(shù)越接近1,意味著兩個概率分布越相似;反之,表示兩個概率分布越不相似。

2.2 巴氏系數(shù)性能指標(biāo)

為便于研究,考慮下圖1所示的單輸入單輸出反饋控制回路。其中y(k),u(k)和v(k)分別代表過程的輸出、輸入和未知隨機(jī)噪聲擾動。

圖1 SISO反饋控制回路結(jié)構(gòu)框圖

圖2 點分布離散近似

假定通過尋優(yōu)算法獲取的基準(zhǔn)數(shù)據(jù)yb為系統(tǒng)性能理想時的數(shù)據(jù),實際輸出的yt是待評估數(shù)據(jù),那么基于巴氏系數(shù)的性能指標(biāo)IBC

(4)

式中,p(x)和q(x)分別為yb與yt的概率分布。基于巴氏系數(shù)的性能指標(biāo)具有如下性質(zhì):

1) 有界性:0≤IBC≤1。當(dāng)基準(zhǔn)與實際輸出概率分布越相近時,性能指標(biāo)越接近于1,系統(tǒng)越接近于理想性能;反之,性能指標(biāo)越接近于0,系統(tǒng)性能越差,需要改進(jìn)。

2) 對稱性:IBC(p,q)=IBC(q,p)。

2.3 IBC性能指標(biāo)計算

如果p,q同時滿足高斯分布,p～N(μp,σp),q～N(μq,σq),那么其巴氏距離可以通過提取它們的均值和方差進(jìn)行簡便計算

根據(jù)式(1)可推出巴氏距離與巴氏系數(shù)的關(guān)系如下

BC(p,q)=e-DB(p,q)

(6)

此時的性能指標(biāo)IBC即可由式(5)和(6)求得

IBC=e-D(p,q)

(7)

對于實際系統(tǒng)中更為常見的非高斯輸出分布,僅依靠均值和方差顯然不能完全描述該隨機(jī)變量的統(tǒng)計特性。此時,將基準(zhǔn)輸出yb和實際輸出yt進(jìn)行離散化處理。

首先對隨機(jī)變量進(jìn)行固定頻率采樣,記錄收集到的有效數(shù)據(jù)并確定其落入的區(qū)間,將區(qū)間均勻且不重疊地劃分為多個子區(qū)間,統(tǒng)計落入每個子區(qū)間的數(shù)據(jù)個數(shù),將頻率記為每個子區(qū)間中心值的概率,即用不同子區(qū)間的中心值概率來近似逼近連續(xù)隨機(jī)變量的概率分布。

對于概率p的計算,直接利用樣本的直方圖對其進(jìn)行離散化處理,具體的實現(xiàn)步驟如下:

1) 以固定頻率采樣x,收集N個有效數(shù)據(jù)點,將所有數(shù)據(jù)都放在區(qū)間S=[E1,EM]內(nèi);

2) 將區(qū)間S平均分為M個子空間(S1∪S2∪…∪SM),每個子空間互不交叉,并將每個子空間的中心點記ei(i=1,2,3,…,M);

3) 統(tǒng)計每個子區(qū)間內(nèi)采樣數(shù)據(jù)的個數(shù)l作為頻數(shù),記為(l1,l2,…,lM);

4) 將頻數(shù)與數(shù)據(jù)總量N(N?M)相比得到頻率,記為P：

(8)

BC=∑Mi=1i(x)i(x)

(9)

3 基準(zhǔn)輸出數(shù)據(jù)的獲取

3.1 基準(zhǔn)輸出的定義

假定圖1系統(tǒng)中的設(shè)定值為0,則系統(tǒng)輸出為

(10)

Gv(z-1)=F(z-1)+z-τR(z-1)=(1+n1z-1+n2z-2+…+nτ-1z-(τ-1))+z-τR(z-1)

(11)

式中,F(z-1)是擾動傳遞函數(shù)Gv的脈沖響應(yīng)系數(shù),R(z-1)是滿足恒等式(11)的剩余函數(shù)。

y(t)=Fv(t)+Lv(t-τ)=

(12)

對于線性非高斯隨機(jī)系統(tǒng),最小熵控制器的目標(biāo)就是最小化輸出變量的熵[12,13]。

H(yt)=H(Fvt+Lvt-τ)≥H(Fvt)

(13)

由式(13)可知,當(dāng)且僅當(dāng)L=0時,系統(tǒng)輸出數(shù)據(jù)的熵達(dá)到最小值,此時控制系統(tǒng)的性能達(dá)到最優(yōu),相應(yīng)的反饋不變量部分即為基準(zhǔn)輸出數(shù)據(jù),記:

yb(t)=Fv(t)=(1+n1z-1+n2z-2+…+nvz-(τ-1))v(t)

(14)

實際輸出數(shù)據(jù)yt可由Simulink仿真平臺采集得到。

3.2 系統(tǒng)參數(shù)辨識和噪聲PDF估計

實際工業(yè)過程中的控制系統(tǒng)相當(dāng)于一個“黑箱”,內(nèi)部的結(jié)構(gòu)參數(shù)均是未知的,而根據(jù)式(14),計算基準(zhǔn)輸出數(shù)據(jù)必須知道模型參數(shù)[n1,n2,…,nv],因此就需要對給定系統(tǒng)進(jìn)行參數(shù)辨識。

圖1所示的隨機(jī)線性系統(tǒng)離散時間模型可描述為

A(z-1)y(k)=z-τB(z-1)u(k)+C(z-1)v(k)

(15)

其中A,B,C可表示為

(16)

其中ana,bnb,cnc和τ是模型的結(jié)構(gòu)參數(shù);na,nb,nc分別是A(z-1),B(z-1)和C(z-1)的階次,τ為系統(tǒng)延遲。

非高斯隨機(jī)系統(tǒng)參數(shù)辨識問題本質(zhì)上是高維參數(shù)空間優(yōu)化的問題,因此系統(tǒng)參數(shù)可以通過優(yōu)化算法獲得。為了獲得系統(tǒng)參數(shù),首先要知道給定系統(tǒng)的階次和時延。

對于遲延τ大小的估計,最簡單常用方法是通過分析輸入ut和輸出yt信號之間的相關(guān)性,即

(17)

應(yīng)用Akaike信息準(zhǔn)則[14]得到系統(tǒng)的階次。對于給定的CARMA模型,AIC準(zhǔn)則為

(18)

3.2.1 改進(jìn)的混合EDA算法

EDA算法是一種隨機(jī)搜索啟發(fā)式算法,用于創(chuàng)建解空間的概率模型,該模型基于模型采樣的解的質(zhì)量進(jìn)行迭代更新[15]。在之前的研究中,EDA算法已經(jīng)暴露出在優(yōu)化過程中存在計算量大、局部搜索能力較差等缺點。本文針對以上缺點進(jìn)行了改進(jìn)。

1) 獲取參數(shù)辨識空間

在參數(shù)未知的情況下,參數(shù)空間過大會影響算法收斂速度,特別是對于高維空間無法保證均勻采樣的情況。因此,在初始化參數(shù)空間之前,需要對模型參數(shù)進(jìn)行初步估計,以確保參數(shù)空間盡可能覆蓋真實參數(shù)。對于非高斯系統(tǒng),沒有很好的參數(shù)初始化方法,因此采用高斯系統(tǒng)參數(shù)初始化方法進(jìn)行粗略估計。

首先將CARMA模型寫成最小二乘形式

y(k)=φT(k)θ+v(k)

(19)

采用遞推增廣最小二乘算法(RELS)對系統(tǒng)參數(shù)進(jìn)行初步估計。目標(biāo)是求目標(biāo)函數(shù)J()達(dá)到極小值時的參數(shù),J()定義如下

J()＾U3]2

(20)

2) WSA算法思想

為了提高算法的局部搜索能力,在傳統(tǒng)EDA算法的中引入了WSA算法。WSA算法[16]是受“摸石頭過河”思想的啟發(fā),首先在“岸邊”慎重考察一下選擇初始起點,向該起點附近鄰域隨機(jī)搜索若干解,找出其中最優(yōu)解作為迭代結(jié)果;然后在以這個點為起點再向附近鄰域隨機(jī)搜索若干個解,選出最優(yōu)解作為第三次迭代結(jié)果,以此類推直到滿足終止條件[12]。

a)起點

在參數(shù)辨識空間已知的情況下,初始選擇的解將會對整個算法產(chǎn)生較大的影響。找初始解的思路就是先在參數(shù)空間內(nèi)均勻隨機(jī)取值產(chǎn)生R個解,找出其中的最優(yōu)解作為起點解。對于(20)中的待辨識參數(shù)向量要在區(qū)間[-3v,+3v]內(nèi)均勻隨機(jī)取值,產(chǎn)生R個解A=[1;2; …;R],并計算這R個解的目標(biāo)函數(shù)值(誤差熵值),根據(jù)目標(biāo)函數(shù)值高低找出最優(yōu)解*記為初始起點解,為方便描述,將起點解記*=[x*1,x*2,…x*n],n=na+nb+nc+1。

b)鄰域搜索策略

根據(jù)“摸石頭過河”的思想,當(dāng)摸到一個“石頭”后必然要以該“石頭”為起點向周圍摸索其它“石頭”,以此類推進(jìn)行搜索。因此當(dāng)?shù)玫狡瘘c解后就在起點解的鄰域半徑Lk內(nèi)隨機(jī)搜索m個領(lǐng)域解組成B。

B=[′1;′2;…;′m]

(21)

3)交叉操作

傳統(tǒng)尋優(yōu)算法得到的最優(yōu)解只包含一個解的信息,這在一定程度上忽略了其它優(yōu)秀解中所包含的信息,改進(jìn)方法就是交叉操作。根據(jù)適應(yīng)度值對以上m個個體進(jìn)行排序,選出前N*(N*?m)個優(yōu)秀解(適應(yīng)度值最優(yōu)的個體),D=[′1;′2;…;′N*]。計算D的中心點式(22)所示

(22)

*=a*+(1-a)θ

(23)

式中:a為[0,1]的隨機(jī)數(shù)。

4)適應(yīng)度值的選擇

對于本文的輸入輸出模型,由式(19)可知,k時刻的誤差可以定義為:

(24)

對于估計的誤差序列e=[e1,e1,…,eL],使用二階Renyi熵來估計過程參數(shù),表示為:

(25)

通過最小化熵的誤差序列,能夠得到最優(yōu)模型參數(shù)和輸出與擾動之間關(guān)系。當(dāng)誤差熵達(dá)到最小值時對應(yīng)的參數(shù)也達(dá)到最優(yōu)。即:

opt=arg minθ∈ΩWH2(e)

(26)

式中,ΩW即為參數(shù)尋優(yōu)空間。

5)算法總結(jié)

HEDA優(yōu)化算法的程序步驟如下:

算法1尋優(yōu)算法的程序步驟:

①用輸入輸出模型描述系統(tǒng),通過分析u(t)和y(t)之間的相關(guān)性來估計系統(tǒng)延遲τ;確定模型階次na,nb,nc;

②初步估計。由RELS算法得到初步的模型參數(shù),并以此確定尋優(yōu)空間ΩW為±3v;

③慎重選擇初始起點。在參數(shù)空間ΩW內(nèi)隨機(jī)生成R個個體A=[1;2; …;R],計算每個個體的適應(yīng)度值(誤差熵),根據(jù)適應(yīng)度值(對應(yīng)的誤差熵最小)選擇最優(yōu)解*0;

WhileI≤nmax

a)計算上述m個解的適應(yīng)度值并按大小進(jìn)行排序CI=[(I)*1;(I)*2;…;(I)*m];

b)從CI中篩選出最優(yōu)個體*I和前N*個適應(yīng)度值最優(yōu)的個體DI=[(I)*1;(I)*2;…;N*(I)*];

c)修改最優(yōu)個體*I。計算DI的中心點(均值)交叉操作并修改最優(yōu)個體的位置為*I=a*I+(1-a)θI,a=rand;

d)If 滿足終止條件

兩次相鄰迭代誤差熵之差小于0.0001,結(jié)束循環(huán);

Else

LI=αLI-1;

I=I+1;

End

End

3.2.2 基于參數(shù)靈敏度分析的算法驗證

考慮如下ARMA系統(tǒng),通過分析改進(jìn)算法初始參數(shù)的敏感性,進(jìn)一步比較HEDA和傳統(tǒng)EDA。

(27)

上述仿真實例的待尋優(yōu)參數(shù)向量可表示為θ=[-1.7,0.7,1.5,0.9]。大多數(shù)的算法初始參數(shù)是基于前人的工作[2-4,12-14]。一般情況下,EDA中初始種群的個體數(shù)為N=1000,每次迭代的最佳個體數(shù)m=200,最大迭代次數(shù)是nmax=120。在WSA中,交叉操作的優(yōu)秀個體數(shù)為N*=30,R=800。本文重點討論收縮系數(shù)α和初始鄰域半徑L0的初始值的敏感性。

表1 基于初始參數(shù)敏感性分析的模型參數(shù)辨識結(jié)果

綜上可得,獲取系統(tǒng)基準(zhǔn)輸出數(shù)據(jù)主要分為如下兩步:

仿真驗證

在Matlab/Simulink中搭建圖1所示的典型SISO系統(tǒng),考慮如下CARMA模型:

(28)

控制器的傳遞函數(shù)為:

(29)

表2 不同噪聲擾動下的參數(shù)辨識結(jié)果

圖3 實際和估計的高斯噪聲PDF

圖4 實際和估計的非高斯噪聲PDF

圖3,4可以看出,無論對于高斯噪聲還是非高斯噪聲,HEDA算法估計出的擾動分布與真實的系統(tǒng)擾動都更加接近;表2中的測試結(jié)果也進(jìn)一步驗證了改進(jìn)混合算法在尋優(yōu)速度和精度上的優(yōu)越性。

當(dāng)?shù)玫较到y(tǒng)參數(shù)和噪聲PDF的估計值后,可以很容易計算出性能指標(biāo),如表3所示。BC是通過優(yōu)化算法求得的基于巴氏系數(shù)的性能指標(biāo);IBC為其理論值;IME是基于最小有理熵的性能指標(biāo);IMV為最小方差指標(biāo)。根據(jù)表3,無論對于高斯還是非高斯噪聲,新指標(biāo)都能準(zhǔn)確地評估系統(tǒng)性能,且與理論值的差異很小。對比最小熵指標(biāo),新指標(biāo)在計算過程中更為簡便,實時性更強(qiáng),更重要的是有效避免了熵值平移不變性的缺點,具有更好的實際意義。

表3 不同噪聲擾動下的性能指標(biāo)

4 總結(jié)

針對最小熵性能指標(biāo)的缺陷,克服巴氏距離無界性的問題,提出了一種基于巴氏系數(shù)的隨機(jī)性能指標(biāo)IBC。高斯擾動下,該指標(biāo)可通過輸出數(shù)據(jù)的均值和方差直接求得;非高斯情況下,利用離散化處理獲得概率分布后求出指標(biāo)近似值。為了更加快速準(zhǔn)確地獲得基準(zhǔn)輸出,對傳統(tǒng)EDA算法進(jìn)行了改進(jìn),通過引入初步估計和WSA算法思想,極大地提高了算法的效率。仿真結(jié)果也驗證了HEDA算法在尋優(yōu)速度和精度上的優(yōu)越性。下一步研究可將新指標(biāo)和尋優(yōu)算法用于串級、非線性等系統(tǒng)的性能評估中。