時滯取值概率未知下的線性時滯系統辨識方法

劉 鑫

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

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

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

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

1 辨識問題闡述

考慮以下線性雙率采樣時滯系統

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

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

其中

并且待辨識的模型參數θ定義為

此外,系統輸出量測過程則可改寫為

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

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

并且

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

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

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

2.1 EM 算法簡介

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

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

其中,Θs表示在第s次迭代中得到的參數估計值.在M-step 中,通過優化目標代價函數以得到新的模型參數估計值

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

算法 1.EM 算法

2.2 基于EM 算法的辨識過程推導

1)E-step

首先構造系統的對數似然函數如下:

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

再結合式 (7),可以進一步化簡為

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

結合式 (16)和式(17),系統的對數似然函數最終化簡為

其中,C2=C1-(M/2)log 2π是參數無關項.此時再基于式 (12),可求得

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

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

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

其中,

2)M-step

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

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

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

此時,函數Q1(θ,r)與模型參數θ呈線性關系.再進一步對Q1(θ,r)求偏導可得

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

通過求解上述拉格朗日函數可得

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

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

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

3 辨識算法驗證

3.1 算法性能驗證

考慮以下雙率采樣的線性時滯系統

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

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

圖1 輸入輸出辨識數據Fig.1 The input-output identification data

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

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

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

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

圖2 經過50 次迭代得到的模型參數估計值Fig.2 The estimated model parameters after 50 iterations

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

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

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

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

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

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

圖5 參數g1和g2的蒙特卡洛辨識結果Fig.5 The identifiedg1andg2in Monte Carlo simulations

圖6 參數g3和g4的蒙特卡洛辨識結果Fig.6 The identifiedg3andg4in Monte Carlo simulations

圖7 噪聲參數的蒙特卡洛辨識結果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時執行算法2,并且得到辨識結果,如圖8、圖9 和表1 所示.

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

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

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

從上述辨識結果可以得出以下結論:

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

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

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

3.2 對比實驗驗證

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

表2 對比實驗結果Table 2 The identification results of the comparison tests

由上述辨識結果可以看出:

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

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

4 結束語

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

1)跳過高斯假設,在更加復雜的數據條件下,例如在非高斯噪聲、自相關噪聲以及量化的輸出條件下研究時滯系統辨識方法;

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