基于高斯混合聚類的切換系統(tǒng)的辨識

柴秀俊 ,王宏偉,,王 林 ,嵇薪儒

(1.新疆大學電氣工程學院,新疆烏魯木齊 830047;2.大連理工大學控制科學與工程學院,遼寧大連 110024)

1 引言

混雜系統(tǒng)是一類復雜系統(tǒng),包含相互影響并相互作用的連續(xù)動力學和離散動力學[1].切換系統(tǒng)是一類典型的混雜系統(tǒng),它在通訊、電力系統(tǒng)和機器人控制等領域中有廣泛的應用[2–3].目前,在研究切換系統(tǒng)的結構特性上已經(jīng)獲得了許多重要的研究成果,比如在切換系統(tǒng)的控制器設計[4–6]與穩(wěn)定性分析[7–9]等方面.而在切換系統(tǒng)的辨識領域,還存在很多問題需要解決.

切換系統(tǒng)一般由若干個子系統(tǒng)和決定它們之間切換的規(guī)則組成,這就導致了從切換系統(tǒng)收集的輸入輸出數(shù)據(jù)來自于不同的子系統(tǒng).因此對于切換系統(tǒng)的辨識,通常需要通過模式檢測計算每個采樣數(shù)據(jù)的運行模式號以確定切換系統(tǒng)的切換時間與子系統(tǒng)的駐留時間,然后再對采樣數(shù)據(jù)進行參數(shù)辨識.

對于切換系統(tǒng)的辨識方法,一般分為子空間方法和基于輸入輸出(input and output,IO)模型的方法.文獻[10]首先引入子空間方法用于分段線性系統(tǒng)的辨識,文中假設切換規(guī)則已知,主要研究局部模式的辨識.對于未知切換規(guī)則的切換系統(tǒng)辨識,文獻[11]在子空間框架中使用變化檢測技術對切換系統(tǒng)進行辨識,但是這種方法需要假設系統(tǒng)只在穩(wěn)態(tài)下進行切換.其它的子空間方法還包括投影子空間方法[12],結構化子空間方法[13]等.對于子空間方法,最大的問題就是要求兩個連續(xù)的系統(tǒng)切換之間存在一個最小停頓時間[10].基于IO模型的方法是另一種簡單而有效的切換系統(tǒng)辨識方法,比如自組織模型法[14]、變分貝葉斯法[15]、聚類方法[16–17]、隨機迭代法[18]和橢球法[19]等.基于IO模型的方法適用于具有任意切換規(guī)則的切換系統(tǒng),很好的解決了子空間方法中要求最小停頓時間問題.

上述文獻提出的辨識方法中大都假設子系統(tǒng)的數(shù)量是先驗已知的,而一般情況下具有未知規(guī)則的切換系統(tǒng),只有輸入數(shù)據(jù)與輸出數(shù)據(jù)是可以采樣的,子系統(tǒng)的數(shù)量、切換時間、子系統(tǒng)的駐留時間往往是未知的.在此基礎上,本文提出一種將高斯混合聚類與遞推增廣最小二乘法結合的兩階段辨識方法對切換系統(tǒng)進行辨識.本文將具有未知切換規(guī)則與未知子系統(tǒng)數(shù)量的切換系統(tǒng)作為研究對象,在模式檢測階段,首先建立高斯混合模型表示采樣數(shù)據(jù)的分布,并使用輪盤法選擇較為合適的初始模型參數(shù).其次,分別計算采樣數(shù)據(jù)屬于各個子系統(tǒng)的后驗概率,同時通過極大似然估計算法迭代更新高斯混合模型中的模型參數(shù),使高斯混合模型最大化地擬合所有采樣數(shù)據(jù)的分布.在此基礎上,通過貝葉斯信息準則(Bayesian information criterion,BIC)確定切換系統(tǒng)中子系統(tǒng)的數(shù)量,同時根據(jù)最大后驗概率準則計算采樣數(shù)據(jù)的運行模式號,從而估計切換系統(tǒng)的切換規(guī)則.在參數(shù)辨識階段,通過遞推增廣最小二乘法估計每個子系統(tǒng)的參數(shù)向量.最后,利用含有有色噪聲的切換系統(tǒng)來驗證所提辨識方法的有效性.

2 問題描述

在存在有色噪聲的情況下,切換系統(tǒng)可以表示為

其中:u(k)和y(k)分別代表第k次均勻采樣的切換系統(tǒng)的輸入和輸出;v(k)表示均值為零、方差為δ2的高斯白噪聲;w(k)為與時間相關的分段切換函數(shù);N代表子系統(tǒng)的數(shù)量;z?1代表單位延遲算子,z?1y(k)=y(k?1);Ai(z?1),Bi(z?1)和Ci(z?1)分別表示對應的系統(tǒng)多項式,具體形式為

在k次采樣時,切換系統(tǒng)被切換至第i個子系統(tǒng),式(1)可以轉換成以下線性回歸模型:

其中:

其中θi和φi分別表示第i個子系統(tǒng)的參數(shù)向量和信息向量.

由于切換系統(tǒng)不提供任何的切換信息,所以必須在辨識子系統(tǒng)參數(shù)之前對切換系統(tǒng)進行模式檢測,以確認切換信息[20].因此切換系統(tǒng)的辨識可以分為兩個階段:1)在模式檢測階段,通過檢測技術估計切換規(guī)則,并通過目標函數(shù)法確定切換模式總數(shù);2)在參數(shù)辨識階段,通過辨識算法獲得每個子系統(tǒng)的參數(shù)估計值.本文將在第3節(jié)中討論模式檢測問題,在第4節(jié)中討論參數(shù)辨識問題.

3 基于高斯混合聚類的模式檢測

聚類分析的算法可以分為劃分法、層次法、基于密度的方法、基于網(wǎng)格的方法以及基于模型的方法[21].其中,基于模型的方法就是給每一個聚類假定一個模型,然后去尋找能夠很好的滿足這個模型的數(shù)據(jù)集.高斯混合聚類就是借助高斯混合模型來表示數(shù)據(jù)分布的一種基于統(tǒng)計模型的聚類方法.本文將通過高斯混合聚類對切換系統(tǒng)進行模式檢測.

由于只有輸入數(shù)據(jù)與輸出數(shù)據(jù)是可以采樣的,所以將輸入數(shù)據(jù)與輸出數(shù)據(jù)定義為IO向量ψ(k):

假設對切換系統(tǒng)進行采樣共得到m組IO向量,構成IO向量 集ψ={ψ(1),ψ(2),···,ψ(m)},則ψ的 分布可由高斯混合模型表示:

其中:N代表子系統(tǒng)的數(shù)量;αi代表第i個子系統(tǒng)對應簇類的混合系數(shù),滿足條件代表切換系統(tǒng)在第k次采樣得到的第i個子系統(tǒng)的IO向量;pζ(ψ(k))代表根據(jù)ψ的分布建立的高斯混合模型.由于系統(tǒng)的輸入信號為具有零均值,單位方差的不相關隨機序列信號,輸出信號與輸入信號呈線性關系,所以ψi(k)的分布由概率密度函數(shù)p(ψi(k)|μi,Σi)表示:

3.1 高斯混合模型初始模型參數(shù)的選擇

由式(3)–(4)可以看出高斯混合模型由混合系數(shù),均值向量與協(xié)方差矩陣所決定.對于高斯混合聚類,其聚類的效果十分依靠高斯混合模型中初始均值向量的選擇,本文將通過輪盤法尋找一組合適的初始均值向量.

從IO向量集ψ中隨機選取一組IO向量作為第一個初始聚類中心,分別計算其它的IO向量與當前聚類中心的IO向量之間的距離:

其中:ψe為當前聚類中心的IO向量,ψf為其它的IO向量,λ為距離系數(shù).

在聚類算法中,為了保證聚類的準確性,通常選擇相互距離較大的一組初始聚類中心,但在ψ中存在異常IO向量的情況下,如果直接通過式(5)選擇相互距離最大的一組初始聚類中心,將導致收斂速度慢或者局部收斂.為了解決這個問題,本文通過計算IO向量的選擇概率來選擇初始聚類中心,IO向量之間的距離越大,越有可能被選為下一個聚類中心.在式(5)的基礎上,IO向量的選擇概率為

在每一輪聚類中心的選擇中隨機生成一個大于0且逼近于0的數(shù)τ,按順序分別令選擇概率p(ψf)減去τ,直到滿足條件p(ψf)?τ>0時選擇該組IO向量作為下一個初始聚類中心.當存在兩個或兩個以上的聚類中心時,式(5)中選取其它的IO向量與當前所有聚類中心的IO向量之間最大的距離.之后重復計算式(5)–(6),直到找到N個聚類中心,并將其構建為高斯混合模型中的初始均值向量.

3.2 運行模式號的計算

通過上式構建的高斯混合模型,用于計算IO向量集ψ中m組IO向量的運行模式號.假設ψj(k)是ψ中的一組IO向量,其中j=1,2,···,N代表需要計算的運行模式號.由貝葉斯定理可知,ψj(k)屬于第i個子系統(tǒng)的后驗概率為

為了使高斯混合模型最大化的擬合樣本集ψ,需要迭代更新模型參數(shù)αi,μi和Σi.由式(3)可得樣本集ψ的極大似然函數(shù):

由于模型參數(shù)在每一輪迭代更新過程中,需要重復使用后驗概率pζ(i|ψj(k)),為了減小算法復雜度,對pζ(i|ψj(k))進行0–1離散化處理,若

則令pζ(i|ψj(k))=1,表示將第k個IO向量ψj(k)歸為第i個子系統(tǒng);否則令pζ(i|ψj(k))=0.至此,將樣本集ψ化分為N個數(shù)據(jù)集C1,C2,···,CN.

對于極大似然函數(shù)L(ψ),分別令

可以得到μi與Σi的更新公式:

其中:ψij(k)表示屬于第i類數(shù)據(jù)集Ci的IO向量;mi表示數(shù)據(jù)集Ci中IO向量的總數(shù),滿足條件m1+m2+···+mN=m.

對于αi的更新公式,由上述推論可知

在迭代更新模型參數(shù)的過程中,當L(ψ)的增長小于給定閾值γ(γ>0)時,模型參數(shù)αi,μi與Σi停止更新.此時根據(jù)式(7),由最大后驗概率原則計算出每個IO向量ψj(k)的運行模式號:

由式(13)計算可得j=i,表示IO向量ψj(k)屬于第i個子系統(tǒng).

3.3 切換模式總數(shù)的確定

在上述的計算過程中,切換模式總數(shù)或者子系統(tǒng)的數(shù)量N是未知的,本文通過比較模型擬合采樣數(shù)據(jù)的優(yōu)良性估計切換模式的總數(shù).

赤池信息準則(Akaike information criterion,AIC)是衡量統(tǒng)計模型擬合優(yōu)良性的一種標準,它提供了權衡估計模型復雜度和擬合數(shù)據(jù)優(yōu)良性的標準[22].通常情況下,AIC 準則定義為

其中:LN(ψ)表示子系統(tǒng)數(shù)量分別為N(N=1,2,···,Q,Q為設置的最大子系統(tǒng)數(shù)量)時IO向量集ψ的極大似然函數(shù),m代表IO向量的總數(shù).

BIC準則是在AIC準則的基礎上,增大了懲罰項,在樣本數(shù)量過多時,可有效的防止模型精度過高造成的模型復雜度過高[22].假設子系統(tǒng)的數(shù)量分別為N=1,2,···,Q,根據(jù)貝葉斯信息準則,可以得到IO向量集ψ的BIC值:

一般而言,式(15)中模型復雜度mln(N)與似然函數(shù)項LN(ψ)會隨著子系統(tǒng)數(shù)量N的增加而增加,當N較小時,似然函數(shù)項的增速大于模型復雜度的增速,從而導致BIC變小.當N過大時,模型過于復雜會出現(xiàn)過擬合的現(xiàn)象,似然函數(shù)項LN(ψ)增速減緩,導致BIC增大.所以從N=1,2,···,Q中選擇使BIC值最小的N,故有

通過式(16)可以計算出子系統(tǒng)的數(shù)量N.

為了便于理解,上述的建模及計算過程,可以總結為圖1.

圖1 模式檢測流程圖Fig.1 Flow chart of mode detection

4 基于遞推增廣最小二乘法的參數(shù)辨識

通過第3節(jié)中的高斯混合聚類方法,可以計算出IO向量的運行模式號,在此之后,本文使用遞推增廣最小二乘法[23–24]獲得式(2)中參數(shù)向量的估計值.由于信息向量φi(k)中的v(k)不可測,所以通過其估計值?v(k)代替,即

則遞推增廣最小二乘法的公式如下:

步驟1初始化參數(shù)PPPi(0),,k=1,其中i=1,2,···,N;

步驟2采集數(shù)據(jù)y(k)和u(k),通過式(21)構造信息向量

步驟3依次計算式(18)–(19)和式(20),獲得參數(shù)向量θi的估計值,并通過式(22)計算?v(k);

步驟4判斷k>m? 若是則完成參數(shù)辨識過程,若否則進入下一步驟,其中m為信息向量的總數(shù);

步驟5k=k+1,返回步驟2進行新一輪的迭代.

5 仿真結果

為了驗證所提方法的有效性,本文以如下包含3個子系統(tǒng)的切換系統(tǒng)作為研究對象:

其中:

在仿真中,將具有零均值,單位方差的不相關隨機序列信號作為輸入u(k);將均值為零,方差為δ2=0.1的高斯白噪聲作為噪聲信號v(k).設置距離系數(shù)λ=2,最大子系統(tǒng)數(shù)量Q=10,似然函數(shù)增長閾值γ=10?5,初始協(xié)方差矩陣Σi=III(其中i=1,···,N,···,Q,III為n階單位矩陣),初始混合系數(shù)

本文采用周期性隨機切換信號,每50個樣本為一個周期,對1000組樣本進行了仿真.圖2展示了在一次獨立實驗中得到的BIC曲線,由圖可以看出當N=3時,BIC為最小值,表明切換系統(tǒng)由3個子系統(tǒng)構成.圖3為模式檢測過程中IO向量對應于每個子系統(tǒng)的后驗概率與實際切換過程的比較圖,為了更加直觀的表示,圖中令所有小于10?9的后驗概率都以10?9表示,對于IO向量,使其后驗概率最大的子系統(tǒng)即為此IO向量所屬的子系統(tǒng).圖4為實際的切換信號與通過模式檢測重建的估計切換信號的對比圖,圖中兩者的匹配率為93.7%.

圖2 不同子系統(tǒng)數(shù)量對應的BIC值Fig.2 The BIC value corresponds to the number of different subsystems

圖3 IO向量對應于每個子系統(tǒng)的后驗概率與實際切換過程Fig.3 The IO vector corresponds to the posterior probabilities of each subsystems and the actual switched process

圖4 實際切換信號及其估計Fig.4 Actual switched signal and its estimation

圖5為通過50次重復實驗統(tǒng)計,將本文所提的方法與投影子空間法[12]和聚類方法[16]進行對比得到的分類錯誤率箱式圖,從圖中可以看出本文所提方法的分類錯誤率主要集中在7%到11%之間,低于其它兩種算法.表1為通過遞推增廣最小二乘法的參數(shù)辨識的統(tǒng)計結果,其中第1列和第2列的項分別是平均值和方差.可以看出,遞推增廣最小二乘法可以較為準確地估計切換系統(tǒng)的參數(shù)向量.

圖5 不同方法分類錯誤率的箱式圖Fig.5 Box plot of classification error rate of different algorithms

圖6展示了在噪聲信號的方差δ2分別為{0.1,0.2,···,1}時,切換系統(tǒng)的噪信比和模式檢測的平均分類錯誤率.從圖中可以看出分類錯誤率會隨著噪信比的增大而增大,當噪聲信號的方差為0.1時,分類錯誤率大致為8.75%,當噪聲信號的方差增大到為1時,分類錯誤率增大到26.7%.結果表明該方法在噪信比較低的情況下表現(xiàn)良好,隨著噪信比的增加,分類錯誤率也會略微增加.

圖6 噪聲信號方差不同時的分類錯誤率與噪信比的統(tǒng)計圖Fig.6 Statistical chart of classification error rate and noise signal ratio when the variance of noise signal is different

從上面的仿真中可以得出以下結論:1)對于具有未知切換規(guī)則的切換系統(tǒng)的辨識,本文所提的方法不僅能夠確定切換系統(tǒng)中子系統(tǒng)的數(shù)量,還能較為準確地計算出采樣數(shù)據(jù)的模式運行號,獲得整個切換規(guī)則運行情況;2) 在對切換系統(tǒng)進行模式檢測之后,利用遞推增廣最小二乘法能夠準確地得到各個子系統(tǒng)與噪聲項的估計參數(shù),也反應出模式檢測過程的有效性.

6 結論

針對具有未知切換規(guī)則與未知子系統(tǒng)數(shù)量的切換系統(tǒng),提出一種基于高斯混合聚類與遞推增廣最小二乘法的二階段辨識方法.本文借助高斯混合模型表示采樣數(shù)據(jù)的分布,并根據(jù)最大后驗概率原則計算采樣數(shù)據(jù)的運行模式號,以確定切換系統(tǒng)的切換時間與子系統(tǒng)的駐留時間.由于在高斯混合聚類的過程中初始均值向量對聚類效果影響較大,所以利用輪盤法選擇一組較為合適的聚類中心作為高斯混合模型的初始均值向量.同時借助采樣數(shù)據(jù)的極大似然函數(shù),根據(jù)BIC準則確定切換子系統(tǒng)的數(shù)量.最后通過遞推增廣最小二乘法獲得了較為準確的參數(shù)估計值.

切換系統(tǒng)的辨識在工業(yè)煉鐵中硅含量的預測[14],配水網(wǎng)絡的故障檢測[25]等場景具有廣泛應用.然而,對于切換系統(tǒng)的辨識依然還存在很多問題,比如在切換系統(tǒng)中含有非穩(wěn)定子系統(tǒng),但可以通過平均駐留時間理論來保證全局穩(wěn)定的情況下的參數(shù)辨識,以及非線性切換系統(tǒng)的辨識問題與收斂性分析.這些問題依然值得思考,接下來我們將進行更深入的研究.