劉 鑫

系統(tǒng)辨識作為一種有效的系統(tǒng)建模手段,在過去的一段時間里受到國內(nèi)外學(xué)者的廣泛關(guān)注[1-5].與傳統(tǒng)的第一原理建模方法相比,系統(tǒng)辨識方法直接通過可測的系統(tǒng)數(shù)據(jù)提煉出與實(shí)際系統(tǒng)等價的數(shù)學(xué)模型,無需深入窺探系統(tǒng)的內(nèi)部復(fù)雜機(jī)理,因而使得建模過程變得簡單、高效[6-11].在系統(tǒng)辨識領(lǐng)域,由于長時間的數(shù)據(jù)傳輸、數(shù)據(jù)獲取方式的不同等因素,常常導(dǎo)致不同的過程變量之間存在未知的時間延遲(時滯).時滯作為一個常見的現(xiàn)實(shí)問題,很容易降低辨識數(shù)據(jù)的數(shù)據(jù)質(zhì)量并且對辨識算法的設(shè)計提出更高的要求[11-13].在辨識算法設(shè)計的過程中,忽略時滯因素的影響會導(dǎo)致有偏的參數(shù)估計,甚至錯誤的辨識結(jié)果[14-16].

在實(shí)際的工業(yè)過程中,過程變量之間的時滯大體上可以分為兩類:固定時滯和時變時滯[13-18].對于包含固定時滯的系統(tǒng)辨識問題,傳統(tǒng)的方法一般是分為兩個步驟進(jìn)行:1)對輸入-輸出數(shù)據(jù)利用相關(guān)分析法,先求解未知的時滯參數(shù);2)基于時滯參數(shù)求解未知的模型參數(shù)[17-21].但是這樣的解決方案實(shí)質(zhì)上是將時滯估計和參數(shù)估計獨(dú)立分開執(zhí)行,很容易造成誤差累計,進(jìn)而影響辨識結(jié)果的精度[17].為了解決誤差累計的問題,文獻(xiàn)[22]在統(tǒng)計學(xué)的框架下并基于極大似然判據(jù)提出了一類有效的迭代辨識算法.利用期望最大化(Expectation maximization,EM)算法首先將辨識問題公式化,將未知的時滯變量當(dāng)作隱含變量來處理.在辨識過程中,不斷計算時滯的后驗(yàn)概率密度函數(shù)(Probability density function,PDF),通過比較后驗(yàn)概率密度函數(shù)的值來判斷未知時滯的取值.該方法的主要優(yōu)點(diǎn)包括:1)給出了顯性的時滯估計公式;2)能夠同時給出時滯參數(shù)和模型參數(shù)的極大似然估計,避免了誤差累計降低辨識結(jié)果的精度.隨后該辨識思想被重新利用在線性變參數(shù)時滯系統(tǒng)過程中,文獻(xiàn)[23]采用多模型方法對線性變參數(shù)時滯系統(tǒng)設(shè)計有效的局部辨識算法.首先在每一個工作點(diǎn)處,選取線性自回歸時滯模型為局部模型,每一個局部模型的時滯參數(shù)各不相同;然后利用期望最大化算法同時估計各子模型的時滯參數(shù)、模型參數(shù)以及有效寬度參數(shù);最后利用輸出插值策略估計系統(tǒng)的全局輸出.在每一個采樣時刻,將系統(tǒng)全局輸出寫成每個子模型輸出數(shù)據(jù)的加權(quán)組合的形式.此外,文獻(xiàn)[24]結(jié)合了冗余法則以及遞歸最小二乘算法針對一類線性時滯系統(tǒng)進(jìn)行了辨識研究.核心思想是擴(kuò)展原系統(tǒng)的維數(shù),然后再辨識擴(kuò)展系統(tǒng)的模型參數(shù).如果某項(xiàng)系數(shù)趨近于零,則認(rèn)為該項(xiàng)在原系統(tǒng)中不存在,并依據(jù)于此來判斷時滯參數(shù).這種辨識思路無形中增加了算法的復(fù)雜度.

針對時變時滯系統(tǒng)辨識而言,首先要確定未知時滯的變化機(jī)制.在現(xiàn)有的文獻(xiàn)中,最常用的變化機(jī)制是首先假定時滯可能的取值范圍,然后假定時滯在取值范圍內(nèi)服從均勻分布,即在每個采樣時刻,時滯取每個值的概率相同[3,14,16,18-19].文獻(xiàn)[25]針對多率采樣時滯系統(tǒng)提出了一類有效的辨識方法.輸入、輸出變量分別采用快率、慢率采樣,且輸入和輸出變量之間存在服從均勻分布的時變時滯.利用期望最大化算法實(shí)現(xiàn)上述系統(tǒng)辨識問題,將未知的時滯變量當(dāng)作隱含變量來處理,上述辨識問題可歸納為:在慢率輸出基礎(chǔ)上同時估計時滯參數(shù)和模型參數(shù).另外,基于辨識得到的有限脈沖響應(yīng)模型還可以進(jìn)一步促進(jìn)輸出丟失情況下的輸出誤差模型和傳遞函數(shù)模型的辨識問題.文獻(xiàn)[14]采用多模型局部辨識的思想針對線性變參數(shù)輸出誤差模型辨識進(jìn)行深入研究.能夠同時給出各子模型的時滯參數(shù)和模型參數(shù)的估計公式,由于目標(biāo)代價函數(shù)與模型參數(shù)呈非線性關(guān)系,采用阻尼牛頓法來迭代搜索模型參數(shù)的最優(yōu)解.此外,文獻(xiàn)[18]在輸入輸出數(shù)據(jù)雙率采樣的條件下考慮了線性狀態(tài)空間時滯模型的辨識問題.首先利用離散化技術(shù)得到待辨識系統(tǒng)的慢率采樣模型;利用數(shù)學(xué)工具將慢率采樣模型轉(zhuǎn)化成能觀標(biāo)準(zhǔn)型;利用卡爾曼濾波算法估計未知的狀態(tài)變量;最后利用隨機(jī)梯度算法和遞歸最小二乘算法辨識模型參數(shù).在辨識過程中考慮了狀態(tài)變量與輸出變量之間存在未知的時變時滯,且假定時滯在可能的取值范圍內(nèi)服從均勻分布.

然而,對于時變時滯系統(tǒng)而言,認(rèn)為時滯在可能的取值范圍內(nèi)服從均勻分布這個假設(shè)通常過于嚴(yán)格,在很多情況下很難得到滿足.為了解決這個問題,本文在時滯取值概率未知的條件下研究線性時滯系統(tǒng)辨識問題.仍舊假定時滯可能的取值范圍先驗(yàn)已知,但是時滯取每個值的概率各不相同,即在時滯的取值范圍內(nèi),時滯取每一個值的比重各不相同.不難看出,時滯服從均勻分布可以看作本文研究內(nèi)容的一個特例,因此本文的研究更加貼近實(shí)際應(yīng)用.采用極大似然算法搭建辨識框架,將時滯當(dāng)作隱含變量來處理.在每個采樣時刻不斷計算時滯的后驗(yàn)概率密度函數(shù),且該函數(shù)作為權(quán)重自適應(yīng)地分配給每一個數(shù)據(jù)點(diǎn),因此提高時滯的估計精度也能促進(jìn)模型參數(shù)估計精度的提升.同時給出模型參數(shù)、噪聲參數(shù)、控制概率向量中的每一個元素和未知的時滯的估計公式,最后利用一個數(shù)值例子驗(yàn)證算法的有效性.

1 辨識問題闡述

考慮以下線性雙率采樣時滯系統(tǒng)

其中,u1:N={uk}k=1,···,N表示快速率采樣的輸入數(shù)據(jù)且假設(shè)采樣周期為Δt;xk表示理想的無噪聲快率采樣輸出;表示慢速率采樣的真實(shí)輸出數(shù)據(jù),它的采樣周期假設(shè)為快率采樣的整數(shù)倍,即fΔt,其中f表示一個任意的正整數(shù);τ1:M={τi}i=1,···,M表示未知的時變時滯;eTi表示輸出測量噪聲且服從零均值的高斯分布,即G(z-1)表示系統(tǒng)的多項(xiàng)式傳遞函數(shù)并且有

其中,m表示系統(tǒng)的階次.結(jié)合式(1)和式(2),系統(tǒng)的模型可改寫為

其中

并且待辨識的模型參數(shù)θ定義為

此外,系統(tǒng)輸出量測過程則可改寫為

由式(6)不難看出yTi服從以下高斯分布

在本文中,將未知的時變時滯τi當(dāng)作隱含變量來處理.首先假設(shè)時滯的取值范圍已知,即[1,J].在每個采樣時刻,假設(shè)時滯的取值為[1,J]中的任意一個整數(shù),但取值概率不同且未知,用數(shù)學(xué)語言描述如下:

并且

不難看出,服從均勻分布的時變時滯系統(tǒng)可看作是本文的一個特例,即

在本文中,慢率采樣的輸出數(shù)據(jù)和快率采樣的輸入數(shù)據(jù)構(gòu)成了可用辨識數(shù)據(jù)集Cobs={yT1:TM,u1:N};利用未知的時變時滯構(gòu)造隱含數(shù)據(jù)集Cmis={τ1:M};待辨識的參數(shù)集構(gòu)造為Θ={θ,β1:J,r}.因此本文的主要任務(wù)可以進(jìn)一步總結(jié)為:基于可用數(shù)據(jù)集Cobs設(shè)計有效的辨識算法,同時估計參數(shù)集Θ和未知的時變時滯τ1:M.

2 基于期望最大化(EM)算法的辨識算法推導(dǎo)

2.1 EM 算法簡介

在系統(tǒng)辨識領(lǐng)域,EM 算法對于解決數(shù)據(jù)丟失問題或者隱含變量問題非常有效.EM 算法是一個迭代優(yōu)化算法,通過不斷重復(fù)執(zhí)行期望步驟(E-step)和最大化步驟(M-step)以實(shí)現(xiàn)待辨識參數(shù)的估計.在E-step 中首先構(gòu)造系統(tǒng)對數(shù)似然函數(shù)為

由于隱含數(shù)據(jù)集不可用,使得計算上述似然函數(shù)變得尤為困難.這里通過計算其相對于隱含變量的期望來構(gòu)造代價函數(shù)如下:

其中,Θs表示在第s次迭代中得到的參數(shù)估計值.在M-step 中,通過優(yōu)化目標(biāo)代價函數(shù)以得到新的模型參數(shù)估計值

重復(fù)執(zhí)行E-step 和M-step,直至參數(shù)收斂[26].期望最大化算法實(shí)質(zhì)上是通過不斷迭代優(yōu)化L1的條件期望函數(shù)來代替直接優(yōu)化L1函數(shù),最終得到所需的模型參數(shù)估計值.算法1 總結(jié)了EM 算法的具體實(shí)施步驟.

算法 1.EM 算法

2.2 基于EM 算法的辨識過程推導(dǎo)

1)E-step

首先構(gòu)造系統(tǒng)的對數(shù)似然函數(shù)如下:

其中,C1=logp(u1:N|Θ)表示一個與參數(shù)無關(guān)的常數(shù)項(xiàng).基于式 (4)和式(6)不難看出,在每個采樣時刻,慢采樣輸出yTi與之前一段時間的輸入變量、未知時滯τi以及模型參數(shù)Θ有關(guān),因此等式右邊的概率密度函數(shù)可進(jìn)一步化簡為

再結(jié)合式 (7),可以進(jìn)一步化簡為

此外,由于時滯τi與模型參數(shù)以及輸入變量無關(guān),因此logp(τ1:M|u1:N,Θ)可以進(jìn)一步化簡為

結(jié)合式 (16)和式(17),系統(tǒng)的對數(shù)似然函數(shù)最終化簡為

其中,C2=C1-(M/2)log 2π是參數(shù)無關(guān)項(xiàng).此時再基于式 (12),可求得

由于時滯變量τi屬于離散變量,因此

觀察式 (20)可得,想要優(yōu)化代價函數(shù)Q(Θ|Θs),首先需要計算時滯取值的后驗(yàn)概率密度函數(shù)p(τi=j|Cobs,Θs).根據(jù)文獻(xiàn)[17]中給出的計算公式可得

為了方便優(yōu)化,可將代價函數(shù)Q(Θ|Θs)拆分成以下形式

其中,

2)M-step

在此步驟中通過優(yōu)化上述代價函數(shù)來得到參數(shù)估計.首先Q1(θ,r)屬于一個復(fù)合函數(shù),它同時是模型參數(shù)θ以及噪聲參數(shù)r的函數(shù),只能采取循環(huán)優(yōu)化的方式得到優(yōu)化的參數(shù)估計.當(dāng)優(yōu)化θ時,令r=rs為定值;反之,當(dāng)優(yōu)化r時,令θ=θs為定值.將函數(shù)Q1(θ,r)對r求偏導(dǎo)并令導(dǎo)數(shù)等于零,得

由于模型參數(shù)θ與代價函數(shù)Q1(θ,r)并非呈線性關(guān)系,因此無法直接利用導(dǎo)數(shù)的方式求得θ的估計公式.結(jié)合式 (3)、(4)和(23),可將代價函數(shù)Q1(θ,r)改寫為

其中,向量?Ti-j定義為

此時,函數(shù)Q1(θ,r)與模型參數(shù)θ呈線性關(guān)系.再進(jìn)一步對Q1(θ,r)求偏導(dǎo)可得

此外,根據(jù)函數(shù)Q2(βj)求解βj的迭代估計公式,需要構(gòu)造拉格朗日函數(shù).首先假設(shè)Lβ為拉格朗日因子,結(jié)合式 (9)可構(gòu)造所需的拉格朗日函數(shù)如下:

通過求解上述拉格朗日函數(shù)可得

并且在每個采樣時刻,未知時變時滯的迭代估計公式如下:

至此,基于期望最大化算法的未知參數(shù)和未知時變時滯估計公式全部推導(dǎo)完成.值得注意的是,在參數(shù)循環(huán)迭代的過程中,只用到了前一次的模型參數(shù)θs和噪聲參數(shù)rs.因此在執(zhí)行算法初始化時,只需設(shè)置θ和r的初值.本文提出的辨識算法可總結(jié)如算法2 所示.

算法 2.時滯取值概率未知下的線性時滯系統(tǒng)辨識算法

3 辨識算法驗(yàn)證

3.1 算法性能驗(yàn)證

考慮以下雙率采樣的線性時滯系統(tǒng)

其中,輸入信號選取為幅值在 [-1,1]之間的高斯白噪聲,即uk=-1+2×N(0,1).輸出數(shù)據(jù)采用慢率采樣且采樣周期是輸入的5 倍,即f=5.在每個慢率采樣時刻,時滯取值范圍分布在[1,4]中,并且按照概率[0.35,0.2,0.25,0.2]取值,即

輸出測量噪聲選取為方差r=0.01 的高斯噪聲,且輸出噪聲與輸入信號相互獨(dú)立.對于輸入數(shù)據(jù)生成N=2 000 個數(shù)據(jù)點(diǎn),則收集到的辨識數(shù)據(jù)如圖1 所示.

圖1 輸入輸出辨識數(shù)據(jù)Fig.1 The input-output identification data

在執(zhí)行辨識算法之前,首先需要將待辨識的參數(shù)初始化.模型參數(shù)θ的初始值選取為

噪聲方差r的初始值選取為

然后執(zhí)行算法2,經(jīng)過50 次迭代之后,得到的時滯分布向量估計值如下:

同時,得到模型參數(shù)的收斂曲線如圖2 所示.

圖2 經(jīng)過50 次迭代得到的模型參數(shù)估計值Fig.2 The estimated model parameters after 50 iterations

通過圖2 可以看出,本文提出的辨識算法具有非常良好的收斂性,經(jīng)過少量的迭代之后,每個模型參數(shù)都能精確地收斂到真實(shí)值附近.此時,未知的噪聲參數(shù)迭代曲線如圖3 所示.

圖3 噪聲方差的收斂曲線Fig.3 The convergence curve of the identified noise variance

此外,結(jié)合式 (31),得到未知時滯的估計結(jié)果如圖4 所示.

圖4 時滯的估計值與真實(shí)值對比Fig.4 The comparison between the estimated and true delays

在圖4 中,藍(lán)色的實(shí)線表示真實(shí)的時滯;紅色的虛線表示估計的時滯.經(jīng)過統(tǒng)計,時滯估計的正確率在91% 以上.由此可見,本文提出的辨識算法能夠同時估計模型參數(shù)以及未知的時變時滯.

為了更加系統(tǒng)地驗(yàn)證本算法的有效性,本文還設(shè)計了蒙特卡洛模擬仿真實(shí)驗(yàn)對算法進(jìn)行更加系統(tǒng)的測試.在接下來的測試中,運(yùn)行100 次獨(dú)立的蒙特卡洛仿真,迭代次數(shù)設(shè)置為50 次.與文獻(xiàn)[1]中類似,每次獨(dú)立的蒙特卡洛仿真的初始值都從以真實(shí)參數(shù)為中心的對稱區(qū)間中選取未知的初值θ0.得到的蒙特卡洛辨識結(jié)果如圖5～7 所示.

圖5 參數(shù)g1和g2的蒙特卡洛辨識結(jié)果Fig.5 The identifiedg1andg2in Monte Carlo simulations

圖6 參數(shù)g3和g4的蒙特卡洛辨識結(jié)果Fig.6 The identifiedg3andg4in Monte Carlo simulations

圖7 噪聲參數(shù)的蒙特卡洛辨識結(jié)果Fig.7 The identified noise variance in Monte Carlo simulations

此外,本文還在不同的噪聲條件下測試了算法的性能.首先定義信噪比 (Signal-to-noise rate,SNR)的計算式為

其中,var(·)表示方差操作符.分別在信噪比等于10 dB,15 dB,20 dB,25 dB,30 dB時執(zhí)行算法2,并且得到辨識結(jié)果,如圖8、圖9 和表1 所示.

表1 不同信噪比下的未知時滯估計精度Table 1 The estimation accuracy of the unknown delays under different SNRs

圖8 不同噪聲下的g1和g2估計結(jié)果Fig.8 The identifiedg1andg2under different noise levels

圖9 不同噪聲下的g3和g4估計結(jié)果Fig.9 The identifiedg3andg4under different noise levels

從上述辨識結(jié)果可以得出以下結(jié)論:

1)從圖2～4 可以看出,本文提出的辨識算法能夠同時給出模型參數(shù)、噪聲方差和時變時滯的估計值.經(jīng)過少數(shù)迭代之后,模型參數(shù)和噪聲參數(shù)可以精確地收斂到真實(shí)值;同時時滯的估計值和真實(shí)值之間的匹配精度達(dá)到91%以上,證明了所提算法的有效性.

2)從圖5～7 可以看出,在蒙特卡洛仿真中,本文提出的辨識算法依然穩(wěn)定,在不同的初始值條件下得到的模型參數(shù)估計以及時滯參數(shù)估計都能精確地收斂到真實(shí)值,這進(jìn)一步證明了所提算法的穩(wěn)定性.對于迭代過程中出現(xiàn)的數(shù)值差異,是由于參數(shù)初始值不同而造成的.

3)從圖8、圖9 和表1 可以看出,當(dāng)信噪比不斷增大時,意味著噪聲不斷減少且數(shù)據(jù)質(zhì)量不斷變好,模型參數(shù)和時滯參數(shù)的估計精度不斷提高,并且隨著數(shù)據(jù)質(zhì)量不斷提升,模型參數(shù)收斂到真實(shí)值的速度也不斷加快.上述辨識結(jié)果也再次驗(yàn)證了所提算法的有效性.

3.2 對比實(shí)驗(yàn)驗(yàn)證

為了進(jìn)一步驗(yàn)證算法2 的有效性,本文還設(shè)計了比較實(shí)驗(yàn).將算法2 分別與兩個基于遞推最小二乘 (Recursive least squares,RLS)辨識方法進(jìn)行比較.在第1 種方法中,假設(shè)時滯參數(shù)精確已知,將這種算法命名為 RLS-accurate;在第2 種方法中,忽略未知的時滯,將這種算法命名為 RLS-ignore.輸出量測噪聲采用信噪比為20 dB的高斯白噪聲,通過執(zhí)行上述三種辨識算法得到的結(jié)果如表2 所示.

表2 對比實(shí)驗(yàn)結(jié)果Table 2 The identification results of the comparison tests

由上述辨識結(jié)果可以看出:

1)本文提出的辨識算法在參數(shù)估計精度上與RLS-accurate 方法相當(dāng).再結(jié)合表1 的辨識結(jié)果可以進(jìn)一步看出,本文提出的辨識算法在估計時滯參數(shù)的同時,也能夠精確地給出模型參數(shù)的估計結(jié)果.

2)由RLS-ignore 方法的辨識結(jié)果可以看出,如果忽略時滯因素的影響會造成有偏的參數(shù)估計,這從側(cè)面證明了本文提出的算法2 的有效性.

4 結(jié)束語

針對時滯取值概率未知的情況,本文基于期望最大化算法提出了一類有效的線性時滯系統(tǒng)辨識算法.將時滯的取值概率當(dāng)成未知的參數(shù)來處理,在辨識過程中所提出的算法能夠有效地同時估計模型參數(shù)、噪聲參數(shù)和時滯參數(shù).仿真結(jié)果表明,本文提出的辨識算法具有快速的收斂速度以及較強(qiáng)的穩(wěn)定性.基于本文的研究內(nèi)容,未來的研究工作可注重以下兩方面:

1)跳過高斯假設(shè),在更加復(fù)雜的數(shù)據(jù)條件下,例如在非高斯噪聲、自相關(guān)噪聲以及量化的輸出條件下研究時滯系統(tǒng)辨識方法;

2)進(jìn)一步考慮相鄰時刻時滯間的相關(guān)性,對于包含相關(guān)時滯 (例如:馬爾科夫時滯)的系統(tǒng)辨識問題進(jìn)行深入研究.此方向的研究甚至可以擴(kuò)展到馬爾科夫切換系統(tǒng)辨識方法的研究.