溫東賓， 張正道， 潘 豐

(江南大學(xué)輕工過程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室，江蘇無錫 214122)

近年來，隨機(jī)分布系統(tǒng)的故障檢測(cè)與隔離問題成為研究的焦點(diǎn)，也提出了許多有效的解決方法。大體可分為如下兩類［1-2］:一是使用統(tǒng)計(jì)學(xué)理論，例如似然函數(shù)和貝葉斯理論，再結(jié)合諸如蒙特卡羅方法、粒子濾波器等的數(shù)學(xué)計(jì)算，對(duì)故障(或系統(tǒng)參數(shù)突變)進(jìn)行預(yù)估;二是基于特定方法的濾波器，例如最優(yōu)化熵濾波的最小化思想，能使估計(jì)誤差滿足各項(xiàng)性能指標(biāo)［3］。

在研究線性高斯系統(tǒng)的故障檢測(cè)與隔離時(shí)典型的方法是卡爾曼濾波器(KF)及改進(jìn)(EKF)。而對(duì)于非線性系統(tǒng)，濾波器的設(shè)計(jì)還需要一定的變換或逼近，以使得誤差方程滿足一定的線性結(jié)構(gòu)［4］。盡管如此，系統(tǒng)會(huì)受到高斯噪聲的影響，可能產(chǎn)生非高斯的輸出信號(hào)。對(duì)于一個(gè)非高斯的參數(shù)，均值和方差均不能有效表征其動(dòng)態(tài)特性。此時(shí)，改進(jìn)的卡爾曼濾波器也只能在一定范圍內(nèi)使用，因?yàn)樗鼉H僅給出了一個(gè)最小化參數(shù)估計(jì)［5］。熵能權(quán)衡隨機(jī)特性，廣泛應(yīng)用于信息、熱力學(xué)和控制理論中。在文獻(xiàn)［6-7］中，用B樣條原理逼近系統(tǒng)輸出的概率密度函數(shù)，隨機(jī)故障檢測(cè)與診斷問題轉(zhuǎn)變成由線性系統(tǒng)支配的動(dòng)態(tài)權(quán)重選擇問題。文獻(xiàn)［8］給出了一種使用可測(cè)的輸出概率密度函數(shù)的魯棒性故障檢測(cè)與診斷方法，解決了一類非線性系統(tǒng)的動(dòng)態(tài)權(quán)重模型。但是，文獻(xiàn)［9］的計(jì)算過程過于繁瑣，結(jié)論需要滿足假設(shè):(1)輸出信號(hào)可測(cè);(2)輸出的概率密度函數(shù)能用B樣條原理近似逼近;(3)控制輸入量與在線權(quán)重之間的動(dòng)態(tài)權(quán)重可建模。

從實(shí)際情況出發(fā)，有兩類故障隔離問題需要考慮:第一類是輸出信號(hào)可測(cè)情況;第二類是輸出的數(shù)據(jù)有部分缺失情況。針對(duì)第一種情況的方法，假設(shè)分析的采樣數(shù)據(jù)集是完整的，可以通過構(gòu)建一個(gè)殘差濾波器;然后，基于熵最優(yōu)化準(zhǔn)則(EOP)對(duì)該濾波器的殘差故障檢測(cè)性能優(yōu)化，其目的是將目標(biāo)故障與非目標(biāo)故障進(jìn)行分離;最后，為了簡(jiǎn)化過程的計(jì)算量，引入輔助映射和Renyi熵來計(jì)算誤差的概率密度函數(shù)。第二類情況反映實(shí)際工業(yè)過程當(dāng)中，采集到的數(shù)據(jù)不完整。引起缺失數(shù)據(jù)的原因很多，如觀測(cè)器故障、觀測(cè)器觀察范圍之外的誤差、數(shù)據(jù)采集系統(tǒng)故障、計(jì)算機(jī)系統(tǒng)崩潰、中斷的輸電線路、本地?cái)?shù)據(jù)的錯(cuò)誤格式化、管理軟件的數(shù)據(jù)錯(cuò)誤以及數(shù)據(jù)備份的故障等。有時(shí)，數(shù)據(jù)在關(guān)鍵位置的值無法獲得，其原因歸咎于觀測(cè)器或采樣的問題，例如分析儀器的采樣頻率與溫度、流速、壓力的測(cè)量頻率相比過低。對(duì)于這類情況下的故障隔離問題，目前的成果依然較少。

文中主要提出一種非線性非高斯帶有數(shù)據(jù)缺失的系統(tǒng)故障隔離方法。首先利用EM迭代算法修補(bǔ)缺失的數(shù)據(jù)，然后基于數(shù)據(jù)驅(qū)動(dòng)模型的思想設(shè)計(jì)最優(yōu)化熵濾波器，實(shí)現(xiàn)目標(biāo)故障的隔離作用。利用一個(gè)數(shù)據(jù)缺失示例，驗(yàn)證了該方法的有效性。

1 對(duì)象模型和故障隔離問題

1．1 對(duì)象模型

構(gòu)建非線性非高斯動(dòng)態(tài)離散系統(tǒng)的差分狀態(tài)方程和輸出動(dòng)態(tài)方程:

其中:xk為系統(tǒng)k時(shí)刻的狀態(tài)量;δk為需要隔離的目標(biāo)故障;ηk為非目標(biāo)故障，由系統(tǒng)方程建模誤差引起;ωk為干擾的噪聲信號(hào);vk為干擾輸出yk的噪聲信號(hào)。δk，ηk，ωk，vk為有界的、互為獨(dú)立的隨機(jī)變量。假設(shè)1 f(·)和h(·)是兩個(gè)已知的波爾可測(cè)和非線性平滑函數(shù)，初值 f(0，0，0)=0;

假設(shè) 2 外部輸入 δk，ηk，ωk，vk是有界、互為獨(dú)立的隨機(jī)變量，且各自的概率密度函數(shù)已知。初始時(shí)刻x0是一個(gè)概率密度函數(shù)已知的隨機(jī)變量。

1．2 動(dòng)態(tài)誤差和故障隔離濾波器

依據(jù)非線性非高斯系統(tǒng)參數(shù)模型，構(gòu)建如下的濾波器參數(shù)模型:

其中，Uk為濾波器控制增益，狀態(tài)誤差表達(dá)式:ek=xk－k，則

定義系統(tǒng)殘差:êk=h(xk，vk)－ h(k，0)。

由假設(shè)1的內(nèi)容可知，系統(tǒng)殘差可寫成:

其中，Lk為確定的常量;ρk(xk，k)為由狀態(tài)誤差引起的不確定量;êk與誤差ek成線性關(guān)系。ρk(xk，k)∈［M1，M2］，因此殘差的熵滿足不等式:

通常情況下，yk是可測(cè)的，誤差ek+1可看成是Uk的參數(shù)表達(dá)式。設(shè)計(jì)適當(dāng)?shù)腢k，使誤差ek+1受目標(biāo)故障 δk的影響最大化，而受變量 ηk，ωk，vk的影響最小，實(shí)現(xiàn)目標(biāo)故障δk隔離。本中與參考文獻(xiàn)不同之處在于:yk是部分?jǐn)?shù)據(jù)缺失的，且yk概率密度函數(shù)未知。文中主要考慮離線系統(tǒng)的故障隔離問題，先修補(bǔ)缺失數(shù)據(jù)，然后計(jì)算yk的概率密度函數(shù)，為最優(yōu)化熵濾波器設(shè)計(jì)做準(zhǔn)備。依據(jù)引言的兩類實(shí)際情況進(jìn)行分析討論:

Case 1:yk任意時(shí)刻可測(cè)。

方程式3的誤差ek+1可以看作兩部分組成:

由假設(shè) 2 已知，xk，δk，ηk，ωk是互為獨(dú)立的變量，其概率密度函數(shù)已知，所以xk+1的概率密 度函數(shù)表達(dá)式:

假設(shè)2給定x0、目標(biāo)故障、非目標(biāo)故障及噪聲的概率密度函數(shù)，經(jīng)過k次迭代積分，獲得xk+1的概率密度函數(shù):

定義積分區(qū)域:假設(shè)τ是給定的常數(shù):

因此，狀態(tài)誤差ek的概率密度函數(shù)寫為

Case 2:yk數(shù)據(jù)缺失。

首先，需要對(duì)缺失數(shù)據(jù)yk使用EM算法，然后計(jì)算概率密度函數(shù)γyk。EM算法需要以下幾個(gè)步驟:

1)初始化分布參數(shù);

2)重復(fù)E-Step和M-Step直至收斂。

E步驟:估計(jì)未知參數(shù)的期望值，給出當(dāng)前的參數(shù)估計(jì)。

M步驟:重新估計(jì)分布參數(shù)，以使得數(shù)據(jù)的似然性最大，給出未知變量的期望估計(jì)。

在迭代收斂過程，缺失數(shù)據(jù)使用估計(jì)值修補(bǔ)完整。令

由方程式5的關(guān)系可知:

在上述條件下，動(dòng)態(tài)誤差的概率密度函數(shù)的表達(dá)式為 γek(τ)= γdk(γ － θ2k)。

證 狀態(tài)誤差ek+1=dk+1+θ2k。

方程式(5)給出了 dk+1是 f(xk，δk，ηk，ωk)和Ukyk兩個(gè)隨機(jī)變量的求和。利用方程式(6)可得γdk+1(τ)，又因?yàn)?ek+1=dk+1+ θ2k，所以

γek(τ)= γdk(τ － θ2k)。

證明完畢。

顯然γek(τ)可以用5個(gè)互為獨(dú)立的概率密度函數(shù) γδk(τ)，γηk(τ)，γωk(τ)，γvk(τ)，γxk－1(τ)和 Uk來表示，是一個(gè)在其積分區(qū)域Ωk上的多重積分。

2 最優(yōu)化熵濾波器

2．1 最優(yōu)化熵原則

定義狀態(tài)誤差ek的熵為

因?yàn)?γek(τ)是可用 γδk(τ)，γηk(τ)，γωk(τ)，γvk(τ)，γxk－1(τ)和Uk表示，所以 H(ek)也可用它們的熵來表示。實(shí)際上，H(ek)是關(guān)于5個(gè)變量的條件概率。

因?yàn)闉V波器的目的是隔離目標(biāo)故障，所有狀態(tài)誤差分成兩部分:目標(biāo)故障引起的誤差e1k+1和非目標(biāo)故障引起的誤差e2k+1。先考慮Case 1，兩部分對(duì)應(yīng)的參數(shù)表達(dá)式:

定理1

積分區(qū)域:假設(shè)給定常數(shù)τ。

Π1:={(τ1，τ2)|f(τ1，τ2，0，0)≤ τ}，

Π2:={(τ1，τ2，τ3)|f(τ1，0，τ2，τ3)≤ τ}

對(duì)于Case 2，同理可得，由方程式(6)可知:

定理2 Case 2對(duì)應(yīng)的概率密度表達(dá)式為

其中

2．2 優(yōu)化算法

假 設(shè) 3 xk∈ ［a，b］ 取 任 意 值 時(shí):?f(xk，δk，ηk，ωk)/?xk≠0 均成立。

給定一個(gè)常數(shù) τ，對(duì)于 f(τ1，τ2，0，0)，f(τ1，0，τ2τ3)，f(τ1，τ2)，f(τ1，τ2，τ3)有可逆函數(shù)，其中:

由可逆函數(shù)的定義可知:

構(gòu)建如下的多維輔助映射:

定理3 對(duì)于Case 1，利用假設(shè)1～3，則在t=k時(shí)，的概率密度表達(dá)式:

其中

得出誤差的概率密度表達(dá)式為

證

由方程組(16)的第一個(gè)方程可知:

將上式代入式(20)即可證明。

定理4 對(duì)于Case 2，利用假設(shè)1～3，在t=k時(shí)，狀態(tài)誤差的表達(dá)式為

ψ1(τ，τ1)，ψ2(τ，τ1，τ2)在定理 3 中已經(jīng)給出。先代入式(11)，然后代入式(21)，求得誤差的概率密度函數(shù)。

2．3 最優(yōu)化熵濾波器設(shè)計(jì)

Renyi熵定義:

當(dāng)α=2時(shí)，

最大則須H2(ek)取最小化。因此，考慮如下的瞬時(shí)代價(jià)函數(shù):

令?Jk/?Uk=0獲得的控制增益Uk，進(jìn)行如下的公式變形:

在Uk－1處對(duì)Jk(Uk)進(jìn)行泰勒級(jí)數(shù)展開:

其中

定理5 對(duì)于Case 1的故障隔離，最優(yōu)濾波器控制增益:

設(shè)定初始值，在任意k時(shí)刻的控制增益，狀態(tài)誤差的概率密度函數(shù)通過以下實(shí)驗(yàn)仿真步驟獲得:

1)定義初始時(shí)刻的參數(shù)值x0，0，U0及初始化分布參數(shù);

2)重復(fù)3)和4)步驟直至收斂，跳至5);

3)E步驟:估計(jì)未知參數(shù)的期望值，給出當(dāng)前的參數(shù)估計(jì);

4)M步驟:重新估計(jì)分布參數(shù)，以使得數(shù)據(jù)的似然性最大，給出未知變量的期望估計(jì);

6)根據(jù)定理4，求得該時(shí)刻的誤差表達(dá)式，關(guān)于Uk的一個(gè)多項(xiàng)表達(dá)式;

7)式(27)，(29)，(30)及

求得t=k時(shí)刻的Uk，ΔUk;

8)k=k+1≤N，返回到5)，否則到9);

9)繪制狀態(tài)誤差ek。

3 實(shí)驗(yàn)仿真

為了故障隔離的仿真，非線性非高斯參數(shù)模型:

其中目標(biāo)故障δk的概率密度函數(shù)(見圖1)表達(dá)式:

非目標(biāo)故障ηk的概率密度函數(shù)(見圖2)表達(dá)式:

干擾信號(hào)wk的概率密度函數(shù)表達(dá)式:

vk為服從N(0，1)的高斯變量。

根據(jù)方程式2，故障隔離濾波器設(shè)計(jì)為

在仿真中，權(quán)重系數(shù)R1=1，R2=20，設(shè)定的初始值為 U0=0，0=0，x0是一個(gè)在區(qū)間［0，1］上的均勻分布，即概率密度函數(shù)為1。按照2．3節(jié)中濾波器設(shè)計(jì)的迭代步驟，得到圖3，4。圖3為只存在目標(biāo)故障δk的系統(tǒng)狀態(tài)誤差的動(dòng)態(tài)響應(yīng);圖4為只存在ηk影響下系統(tǒng)狀態(tài)誤差的動(dòng)態(tài)響應(yīng)。

圖3 存在目標(biāo)故障時(shí)的誤差響應(yīng)Fig．3 Response when the target fault occurs

圖4 只存在非目標(biāo)故障時(shí)的誤差響應(yīng)Fig．4 Response when the nuisance fault occurs

圖3，4仿真是針對(duì)Case 1的仿真結(jié)果。顯然在數(shù)據(jù)完整的情況下，故障隔離目標(biāo)故障的效果較好。對(duì)于 Case 2，預(yù)計(jì)對(duì)y值依次采用5%，10%，20%的缺失情況進(jìn)行EM修補(bǔ)并故障隔離仿真，圖5，6，7是目標(biāo)故障隔離的仿真。

故障的發(fā)生時(shí)間、振動(dòng)幅值可以清晰觀測(cè)得到，為了分析數(shù)據(jù)方便，構(gòu)建表1進(jìn)行對(duì)比描述。從表1可以看出，缺失5% 的數(shù)據(jù)時(shí)，系統(tǒng)故障隔離仍能發(fā)揮作用，但是有1．112 s的時(shí)間延遲，這是數(shù)據(jù)修補(bǔ)算法過程帶來的時(shí)間消耗。圖6反應(yīng)因?yàn)槿笔?shù)值過大，誤差相應(yīng)的幅值開始變大，時(shí)間延遲更加增大(延遲2．451 s)。圖7因?yàn)槿笔?shù)值達(dá)到了20%，故障隔離已經(jīng)無法準(zhǔn)確地判斷故障發(fā)生的時(shí)刻，且因?yàn)檎`差振幅過大，影響故障信號(hào)的檢測(cè)。圖8，9，10是非目標(biāo)故障仿真示意圖，通過詳細(xì)觀察對(duì)比可以發(fā)現(xiàn)幅值和信號(hào)跳變時(shí)刻的改變。

表1 圖5，6，7的仿真結(jié)果比較Tab．1 Comparative simulation result of Fig．5，6，7

為了便于說明，建立表2分別對(duì)故障的發(fā)生時(shí)刻與振幅變動(dòng)加以對(duì)比分析。從表2的數(shù)據(jù)可以看出，圖8為5% 的數(shù)據(jù)缺失，非目標(biāo)故障發(fā)生時(shí)刻有1．176 s的延遲，且平穩(wěn)狀態(tài)時(shí)的幅值波動(dòng)增大;圖9為10% 的數(shù)據(jù)缺失，故障檢測(cè)的時(shí)延增大為2．578 s，誤差波動(dòng)的幅值更加劇烈;圖10為20%的缺失數(shù)據(jù)，誤差的幅值波動(dòng)過大，與非目標(biāo)故障的幅值接近，無法進(jìn)行故障檢測(cè)。

非目標(biāo)故障的隔離在y缺失5%，10%，20% 時(shí)的仿真可以看出，與目標(biāo)故障隔離相同的實(shí)驗(yàn)結(jié)論，隨著缺失數(shù)據(jù)的增大，誤差幅值漸漸增大，故障檢測(cè)延遲越來越大，直到數(shù)據(jù)缺失20%，非目標(biāo)故障抑制作用失效。

表2 圖8，9，10的仿真結(jié)果比較Tab．2 Comparative simulation result of Fig．8，9，10

4 結(jié)語(yǔ)

在現(xiàn)有故障隔離方法的基礎(chǔ)上，提出了在輸出信號(hào)Y存在數(shù)據(jù)缺失的故障隔離問題。首先，使用EM算法，通過E-Step和M-Step迭代直至收斂，得到的估計(jì)值填補(bǔ)到缺失集合內(nèi)。其次，基于修補(bǔ)完整的數(shù)據(jù)驅(qū)動(dòng)模型，構(gòu)建濾波器，可實(shí)現(xiàn)對(duì)系統(tǒng)狀態(tài)的估計(jì)。然后，狀態(tài)誤差由非線性非高斯系統(tǒng)模型和濾波器模型來表達(dá)，其關(guān)系用概率密度函數(shù)表達(dá)，再對(duì)概率密度函數(shù)用最優(yōu)化準(zhǔn)則進(jìn)行優(yōu)化，通過迭代優(yōu)化，可使得目標(biāo)故障存在時(shí)狀態(tài)誤差概率密度函數(shù)的熵最大化，同時(shí)保證在只有非目標(biāo)故障時(shí)最小化，實(shí)現(xiàn)目標(biāo)故障的隔離作用。

文中仍然具有一定的局限性，EM算法是假設(shè)缺失以隨機(jī)缺失機(jī)制，修補(bǔ)方法的好壞，影響著基于數(shù)據(jù)驅(qū)動(dòng)的濾波器設(shè)計(jì)，進(jìn)而影響故障隔離的效果。

［1］周靖林，岳紅，王宏．基于有理平方根B樣條模型的概率密度函數(shù)形狀控制［J］．自動(dòng)化學(xué)報(bào)，2005，31(3):343-351．ZHOU Jing-ling，YUE Hong，WANG Hong．Shaping of output PDF based on the rational square-root B-spline model［J］．Acta Automatica，2005，31(3):343-351．(in Chinese)

［2］ZHOU J L，WANG H．Optimal tracking control of the output probability density functions:square root B-spline model［J］．Control Theory and Applications，2005，22(3):369-376．

［3］YUE H，WANG H．Recent developments in stochastic distribution control:a review［J］．J Meas Control，2003，36(7):209-216．

［4］GUO L，YIN L P，WANG H，et al．Entropy optimization filtering for fault isolation of nonlinear non-gaussian stochastic systems［J］．IEEE Transactions of Automatic Control，2009，54(4):804-810．

［5］GUO L，WANG H．Minimum entropy filtering for multivariate stochastic systems with non-gaussian noises［J］．IEEE Transactions on Automatic Control，2006，51(4):695-700．

［6］WANG H．Bounded stochastic distributions control for pseudo ARMAX stochastic systems［J］．IEEE Transactions on Automatic Control，2001，46(3):486-490．

［7］WANG H．Fault diagnosis and fault tolerant control for non-gaussian stochastic systems with random parameters［J］．Springer Tracts in Advanced Robotics，2003(1):59-84．

［8］XU Zheng-guo，JI Yin-dong，ZHOU Dong-hua．A new real-time reliability prediction method for dynamic systems based on on-line fault prediction［J］．IEEE Transactions on Reliability，2009，58(3):523-538．

［9］LI T，GUO L．Optimal fault-detection filtering for non-Gaussian systems via output PDFs［J］．IEEE Transactions on Systems，Man and Cybernetics A，2009，39(2):476-481．