王澤洲，陳云翔，蔡忠義，王莉莉，項(xiàng)華春

(1.空軍工程大學(xué) 裝備管理與無(wú)人機(jī)工程學(xué)院，陜西 西安 710051;2.中國(guó)人民解放軍93920部隊(duì)，陜西 漢中 723200)

隨著科技的進(jìn)步，航空、航天、核能等諸多高新技術(shù)領(lǐng)域得到快速發(fā)展，同時(shí)也對(duì)相關(guān)裝備提出了更高的安全性要求。為了切實(shí)提升裝備的可靠性，減少系統(tǒng)運(yùn)營(yíng)維護(hù)費(fèi)用，預(yù)測(cè)與健康管理(Prognostics and Health Management,PHM)技術(shù)應(yīng)運(yùn)而生，并取得了良好效益[1-3]。PHM技術(shù)的核心是通過監(jiān)測(cè)系統(tǒng)關(guān)鍵組成設(shè)備的性能退化數(shù)據(jù)來(lái)預(yù)測(cè)其剩余使用壽命，并提出相應(yīng)的維修保障策略，從而提升裝備運(yùn)行的安全性與可靠性，降低故障或失效帶來(lái)的風(fēng)險(xiǎn)。

傳統(tǒng)的維修決策研究多基于隨機(jī)過程描述設(shè)備的退化規(guī)律，進(jìn)而依據(jù)設(shè)備的可靠性指標(biāo)進(jìn)行維修策略優(yōu)化。Kallen等[4]采用Gamma過程描述設(shè)備的退化過程，并分析了檢測(cè)周期與預(yù)防性維修閾值對(duì)失效風(fēng)險(xiǎn)和費(fèi)用的影響。van der Weide等[5]通過幾何過程預(yù)測(cè)設(shè)備的退化趨勢(shì)，并優(yōu)化了設(shè)備的維修時(shí)間。Ho等[6]則基于Markov過程構(gòu)建設(shè)備性能退化狀態(tài)轉(zhuǎn)移方程，并以費(fèi)用為目標(biāo)函數(shù)建立優(yōu)化模型，從而確定最優(yōu)檢測(cè)時(shí)間和維修時(shí)間。然而，上述研究對(duì)設(shè)備實(shí)際退化過程的描述并不完善，均是基于經(jīng)驗(yàn)公式構(gòu)建退化模型，而未能充分利用設(shè)備運(yùn)行過程中的狀態(tài)監(jiān)測(cè)數(shù)據(jù)來(lái)進(jìn)行退化建模，這將導(dǎo)致所建模型難以準(zhǔn)確反映設(shè)備的真實(shí)退化規(guī)律，進(jìn)而影響維修決策結(jié)果的科學(xué)性與合理性。

針對(duì)傳統(tǒng)維修決策方法存在的不足，利用狀態(tài)監(jiān)測(cè)數(shù)據(jù)準(zhǔn)確構(gòu)建設(shè)備退化模型，并基于設(shè)備剩余壽命預(yù)測(cè)信息進(jìn)行維修決策的方法開始逐步得到研究者的關(guān)注[7-9]。Guo等[7]通過引入殘余退化量描述維修活動(dòng)對(duì)設(shè)備退化量的影響，并推導(dǎo)出對(duì)應(yīng)的剩余壽命函數(shù)，進(jìn)而通過建立優(yōu)化模型確定最優(yōu)預(yù)防性維護(hù)閾值。Zhang等[8]通過引入改善因子描述維修活動(dòng)對(duì)設(shè)備退化速率的影響，并以退化速率為對(duì)象分析了設(shè)備的最優(yōu)維修策略。然而，文獻(xiàn)[7]和文獻(xiàn)[8]僅單獨(dú)研究了退化速率或退化量對(duì)維修決策的影響，未能分析二者對(duì)維修決策的綜合作用。為了進(jìn)一步提升維修決策的科學(xué)性，裴洪等[9]在綜合考慮維修活動(dòng)對(duì)退化速率和退化量影響的基礎(chǔ)上，基于設(shè)備的剩余壽命預(yù)測(cè)信息，構(gòu)建了維修決策模型，從而確定了最優(yōu)檢測(cè)時(shí)間和預(yù)防性維護(hù)閾值，降低了設(shè)備的壽命周期費(fèi)用。然而，文獻(xiàn)[9]所提維修決策方法僅適用于可修復(fù)設(shè)備，難以滿足不可修復(fù)設(shè)備制定最優(yōu)維修策略的需求；且該方法在剩余壽命預(yù)測(cè)過程中認(rèn)為設(shè)備性能退化失效閾值滿足固定值假設(shè)，忽略了不確定失效閾值對(duì)剩余壽命預(yù)測(cè)的影響，這可能降低剩余壽命預(yù)測(cè)的準(zhǔn)確性[10-12]，進(jìn)而對(duì)制定科學(xué)合理的維修策略產(chǎn)生消極影響。此外，目前針對(duì)不確定失效閾值分布系數(shù)確定方法的研究尚不充分，文獻(xiàn)[13]提出采用極大似然估計(jì)(Maximum Likelihood Estimation,MLE)對(duì)不確定失效閾值的分布系數(shù)進(jìn)行估算，然而該方法在似然函數(shù)不存在解析形式時(shí)，僅能通過泰勒級(jí)數(shù)展開得到近似解，降低了參數(shù)估計(jì)的準(zhǔn)確性，也對(duì)設(shè)備最優(yōu)維修決策造成了不利影響。

針對(duì)上述問題，本文以不可修復(fù)設(shè)備為研究對(duì)象，開展了基于剩余壽命預(yù)測(cè)數(shù)據(jù)與不確定失效閾值的設(shè)備最優(yōu)替換策略研究，主要?jiǎng)?chuàng)新點(diǎn)有：

1)基于期望最大(Expectation Maximization,EM)算法提出一種新型不確定失效閾值分布系數(shù)估計(jì)法，相較于傳統(tǒng)MLE方法，能有效提升參數(shù)估計(jì)的準(zhǔn)確性；

2)基于設(shè)備的剩余壽命預(yù)測(cè)數(shù)據(jù)，依據(jù)更新報(bào)酬理論建立維修決策模型，分析不確定失效閾值對(duì)最優(yōu)替換時(shí)機(jī)的影響。

1 設(shè)備剩余壽命預(yù)測(cè)

1.1 退化建模

Wiener過程可以準(zhǔn)確描述具有非單調(diào)退化特征設(shè)備的退化規(guī)律，且具備良好的數(shù)學(xué)特性，現(xiàn)已被廣泛應(yīng)用于退化建模研究[14]。

理想狀態(tài)下，Wiener過程可表示為：

X(t)=X(0)+at+bB(t)

(1)

其中：X(t)表示設(shè)備在t時(shí)刻的性能退化量；a為漂移系數(shù)；b為擴(kuò)散系數(shù)；B(t)表示標(biāo)準(zhǔn)布朗運(yùn)動(dòng)；X(0)表示初始時(shí)刻的性能退化量，通常認(rèn)為X(0)=0。

在實(shí)際運(yùn)行過程中，受設(shè)備內(nèi)外部應(yīng)力影響以及非理想測(cè)量手段制約，通常難以準(zhǔn)確描述設(shè)備的退化過程。為了解決該問題，現(xiàn)有研究多通過將非線性、個(gè)體差異以及測(cè)量誤差融入退化建模過程，形成改進(jìn)的Wiener退化模型，其具體表達(dá)式如式(2)所示。

Y(t)=aγ(t,θ)+bB(t)+ε

(2)

1.2 參數(shù)估計(jì)

1.2.1 退化模型參數(shù)估計(jì)

(3)

令Y=[ΔY1,ΔY2,…,ΔYN]，則Y表示全部性能退化監(jiān)測(cè)數(shù)據(jù)，其對(duì)數(shù)似然函數(shù)可表示為：

(4)

lnL(Y)=

(5)

(6)

(7)

(8)

(9)

1.2.2 不確定失效閾值分布系數(shù)估計(jì)

在真實(shí)使用環(huán)境下，同類設(shè)備不同個(gè)體間往往存在差異性，這種個(gè)體差異在退化過程中多體現(xiàn)為設(shè)備退化速率(漂移系數(shù)a)的隨機(jī)性，而在失效過程中則體現(xiàn)為失效閾值ω的不確定性。例如，彈簧的伸縮極限、銑刀的磨損上限、風(fēng)機(jī)的振動(dòng)閾值等，不同個(gè)體間失效閾值相似但不完全相等，難以采用一個(gè)固定值進(jìn)行明確，因而多采用具有不確定性的區(qū)間值進(jìn)行描述。為了科學(xué)分析失效閾值的不確定性，現(xiàn)有研究多采用隨機(jī)變量來(lái)描述設(shè)備的不確定失效閾值，其中正態(tài)隨機(jī)變量成了當(dāng)前的研究熱點(diǎn)[11-12]。

在設(shè)備的實(shí)際退化過程中，易知設(shè)備性能退化的失效閾值應(yīng)大于設(shè)備失效前一時(shí)刻的性能退化量X(t*),0

f(ω)=

(11)

其中：Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)。

由于Φ(·)不存在解析表達(dá)式，難以采用傳統(tǒng)MLE法來(lái)對(duì)其進(jìn)行參數(shù)估計(jì)，為此，本文提出一種基于EM算法的不確定失效閾值分布系數(shù)估計(jì)方法。EM算法針對(duì)缺失/隱含數(shù)據(jù)情形下的參數(shù)估計(jì)具有良好效果，因此適用于估算截?cái)嗾龖B(tài)分布的參數(shù)估計(jì)值。

E步：對(duì)式(12)求虛擬失效閾值ω′的期望，可得

(13)

由截?cái)嗾龖B(tài)分布的性質(zhì)可知，對(duì)于任意ω′i均滿足：

(14)

(15)

(16)

(17)

(18)

令式(17)、式(18)等于零，可得

(19)

由于R未知，因此還需計(jì)算E(R)。文獻(xiàn)[15]給出了E(R)的計(jì)算方法，如式(21)所示。

(21)

將式(21)代入式(19)、式(20)即可得到M步的迭代公式：

(22)

其中：

(24)

1.3 剩余壽命分布推導(dǎo)

設(shè)備的壽命通常被定義為性能退化量首次達(dá)到失效閾值的時(shí)間，也被稱為首達(dá)時(shí)間(First Hitting Time，F(xiàn)HT)?；谏鲜龆x，設(shè)備的壽命可表示為：

T=inf{t:X(t)≥ω|X(0)<ω}

(25)

對(duì)于式(2)所描述的退化模型，在不考慮測(cè)量誤差ε的前提下，可證明其壽命T近似服從逆高斯分布，對(duì)應(yīng)的概率分布如式(26)所示[14]。

(26)

進(jìn)一步可推導(dǎo)出tk時(shí)刻設(shè)備剩余壽命lk概率分布如式(27)所示。具體證明過程詳見文獻(xiàn)[16]。

fLk|ω,X1:k(lk|ω,X1:k)≈

(27)

其中：xk表示tk時(shí)刻目標(biāo)設(shè)備所對(duì)應(yīng)的真實(shí)性能退化數(shù)據(jù)；X1:k則表示直至tk時(shí)刻所獲取的全部真實(shí)性能退化數(shù)據(jù)；

ψ(lk)=γ(tk+lk,θ)-γ(tk,θ)

(28)

(29)

(30)

引理1的證明過程詳見文獻(xiàn)[16]。

則基于全概率公式，可得：

fLk|ω,Y1:k(lk|ω,Y1:k)=

P(X1:k|Y1:k)P(a|Y1:k)dxkda=

Ea{Exk[fLk|ω,a,X1:k(lk|ω,a,X1:k)]}

(31)

其中：Y1:k表示直至tk時(shí)刻所獲取的全部性能退化數(shù)據(jù)監(jiān)測(cè)值；P(·)為求概率。

令xk=D1,a=D2,E=β(lk),F=ψ(lk),G=b2lk，利用引理1，則可求出固定失效閾值條件下設(shè)備剩余壽命的PDF為：

(32)

為了推導(dǎo)考慮不確定失效閾值條件下設(shè)備剩余壽命的PDF，本文給出如下定理1。

定理1若D～TN(μ,σ2),E,F∈R,G∈R+,則

(33)

定理1的證明過程可由文獻(xiàn)[16]中引理1的證明經(jīng)變形得到，這里不再進(jìn)行單獨(dú)推導(dǎo)。

基于全概率公式，若Y1:k已知，則考慮不確定失效閾值條件下設(shè)備的剩余壽命可表示為：

fLk|Y1:k(lk|Y1:k)

=Eω[fLk|ω,Y1:k(lk|ω,Y1:k)]

(34)

若令ω=D,H4=E,H2=F,H1=G，將其代入式(33)，可得設(shè)備剩余壽命的PDF為：

fLk|Y1:k(lk|Y1:k)≈

(35)

其中：

(36)

H2=yk+μaψ(lk)

(37)

(38)

(39)

基于上述分析，可得考慮不確定失效閾值條件下設(shè)備剩余壽命的累積分布函數(shù)為：

(40)

2 最優(yōu)替換策略決策模型

基于更新報(bào)酬理論[17]建立決策模型，進(jìn)而確定最優(yōu)替換時(shí)間。具體決策模型可表示為：

(41)

其中，C(τ)表示τ時(shí)刻進(jìn)行替換操作對(duì)應(yīng)的期望壽命周期費(fèi)用率，EC表示設(shè)備運(yùn)行壽命周期總費(fèi)用的期望，ET表示設(shè)備運(yùn)行總時(shí)間的期望。進(jìn)一步分析可得：

EC=c1P(t>τ-tk|Y1:k)+c2P(t<τ-tk|Y1:k)

(42)

(43)

其中：c1表示預(yù)防性替換的費(fèi)用，c2表示失效性替換的費(fèi)用，tk為設(shè)備當(dāng)前運(yùn)行時(shí)間，τ為替換時(shí)刻，f(l)為變量l的函數(shù)且l=τ-tk。

基于Wiener過程的性質(zhì)，易知lk=τ-tk，則f(l)等價(jià)于第1.3小節(jié)中設(shè)備剩余壽命的概率密度函數(shù)fLk|Y1:k(lk|Y1:k)，進(jìn)一步分析可知P(t<τ-tk|Y1:k)等價(jià)于FLk|Y1:k(lk|Y1:k)。進(jìn)而可將替換策略決策模型改寫為：

(44)

具體證明過程如下：

EC=c1P(t>τ-tk|Y1:k)+c2P(t<τ-tk|Y1:k)

=c1[1-P(t<τ-tk|Y1:k)]+c2P(t<τ-tk|Y1:k)

=c1+(c2-c1)P(t<τ-tk|Y1:k)

=c1+(c2-c1)FLk|Y1:k(l|Y1:k)

(45)

=tk+lk[1-FLk|Y1:k(lk|Y1:k)]+lkFLk|Y1:k(lk|Y1:k)-

(46)

由式(45)除以式(46)即可證明式(44)。

通過求解式(44)即可得到設(shè)備的最優(yōu)替換時(shí)機(jī)。

3 算例分析

3.1 仿真退化數(shù)據(jù)分析

圖1 仿真退化軌跡Fig.1 Simulation of degradation path

(a) μω

為了驗(yàn)證前文所提基于EM算法的不確定失效閾值分布系數(shù)估計(jì)法較傳統(tǒng)MLE方法更具優(yōu)勢(shì)，本文引入均方誤差(Mean Squared Error,MSE)作為判別標(biāo)準(zhǔn)進(jìn)行分析。MSE的定義式為：

(47)

為了進(jìn)一步消除仿真結(jié)果的隨機(jī)性，在原有仿真參數(shù)的基礎(chǔ)上，再分別仿真出30組、300組退化數(shù)據(jù)，對(duì)應(yīng)得到30個(gè)、300個(gè)仿真退化失效閾值數(shù)據(jù)，并采用EM與MLE方法分別進(jìn)行參數(shù)估計(jì)，得到參數(shù)估計(jì)結(jié)果如表1所示。

表1 不確定失效閾值分布參數(shù)估計(jì)Tab.1 Estimation of uncertain failure threshold distribution parameters

表1中MLE對(duì)應(yīng)的失效閾值分布系數(shù)估計(jì)值由MATLAB軟件中normfit命令求出。由表1可知，同等仿真數(shù)據(jù)量條件下，EM算法對(duì)應(yīng)的參數(shù)估計(jì)值較MLE方法更貼近于實(shí)際仿真參數(shù)，且MSE值更小，表明EM算法具有更高的估計(jì)準(zhǔn)確性。進(jìn)一步分析可以發(fā)現(xiàn)，MLE方法對(duì)仿真數(shù)據(jù)量較為敏感，隨著仿真數(shù)據(jù)量的增多，MLE方法的參數(shù)估計(jì)值逐步接近于真實(shí)值，且對(duì)應(yīng)MSE值逐步減??；而EM算法對(duì)仿真數(shù)據(jù)量變化的穩(wěn)定性更好，隨著仿真數(shù)據(jù)的增多，EM算法估計(jì)結(jié)果波動(dòng)較小，參數(shù)估計(jì)誤差變化不明顯?；谏鲜龇治隹梢宰C明，在中、小樣本條件下，EM算法的準(zhǔn)確性要明顯優(yōu)于傳統(tǒng)的MLE方法。而在工程實(shí)際中，退化數(shù)據(jù)往往具有小樣本特性，進(jìn)一步說(shuō)明了EM算法具有較高的工程應(yīng)用價(jià)值。

為了進(jìn)一步分析失效閾值的不確定性對(duì)設(shè)備剩余壽命預(yù)測(cè)與維修決策的影響，選取4號(hào)設(shè)備為目標(biāo)設(shè)備進(jìn)行研究(目標(biāo)設(shè)備壽命為1 h，失效閾值為2.45)。為便于描述，記考慮不確定失效閾值的最優(yōu)替換策略模型為M0，考慮固定失效閾值的最優(yōu)替換策略模型為M1。則針對(duì)M0與M1模型，不同狀態(tài)監(jiān)測(cè)時(shí)刻(0.2 h、0.4 h、0.6 h、0.8 h)對(duì)應(yīng)剩余壽命預(yù)測(cè)情況如圖3所示。

(a) tk=0.2 h

由圖3可知，在不同狀態(tài)監(jiān)測(cè)時(shí)刻，M0對(duì)應(yīng)的剩余壽命PDF均可以包含設(shè)備的真實(shí)剩余壽命，而M1對(duì)應(yīng)的PDF均無(wú)法包含設(shè)備的真實(shí)剩余壽命，表明M0較M1的剩余壽命預(yù)測(cè)準(zhǔn)確性更高，體現(xiàn)了在剩余壽命預(yù)測(cè)過程中考慮不確定失效閾值的必要性。若假設(shè)c1=55元、c2=95元，結(jié)合上述分析得到的剩余壽命預(yù)測(cè)數(shù)據(jù)，將其代入式(44)，即可確定最優(yōu)替換時(shí)機(jī)，具體結(jié)果如圖4所示。

圖4 M0、M1對(duì)應(yīng)的最優(yōu)替換策略Fig.4 Optimal replacement strategy under M0 and M1

由圖4可知，M1模型的最優(yōu)替換時(shí)機(jī)為τ=0.7 h，對(duì)應(yīng)的最小期望壽命周期費(fèi)用率為8.538元/h；M0模型的最優(yōu)替換時(shí)機(jī)為τ=0.8 h，對(duì)應(yīng)的最小期望壽命周期費(fèi)用率為7.655元/h。由此可以說(shuō)明，在確定設(shè)備最優(yōu)替換策略的過程中，考慮不確定失效閾值將有助于延長(zhǎng)設(shè)備的運(yùn)行時(shí)間，減少維修保障的費(fèi)用消耗。

3.2 真實(shí)退化數(shù)據(jù)分析

基于NASA公開數(shù)據(jù)源中的銑刀退化數(shù)據(jù)進(jìn)行分析(如圖5所示)[18]。銑刀是典型的不可修復(fù)組件，當(dāng)銑刀的磨損量超過一定閾值時(shí)，銑刀發(fā)生失效。為避免銑刀失效造成的產(chǎn)品不合格問題與生產(chǎn)安全問題，需要實(shí)時(shí)監(jiān)測(cè)其退化狀態(tài)，并適時(shí)采取更換策略，確保生產(chǎn)的安全性與經(jīng)濟(jì)性。

圖5 銑刀退化數(shù)據(jù)Fig.5 Milling degradation data

假設(shè)銑刀預(yù)防性替換費(fèi)用為c1=10元，失效性替換費(fèi)用為c2=100元，且當(dāng)前目標(biāo)銑刀已完成了5次銑削操作。基于本文第2節(jié)建立的最優(yōu)維修決策模型，可得目標(biāo)銑刀的最優(yōu)替換時(shí)機(jī)，詳見圖6。

圖6 銑刀最優(yōu)替換策略Fig.6 Optimal replacement strategy of milling

由圖6可知：M0求解得到的銑刀最優(yōu)替換時(shí)機(jī)為第13次，對(duì)應(yīng)壽命周期費(fèi)用率為0.769 2元/次；而M1得到的銑刀最優(yōu)替換時(shí)機(jī)為第12次，對(duì)應(yīng)壽命周期費(fèi)用率為1.106 3元/次。由此可知，本文所提最優(yōu)維修決策方法能夠延長(zhǎng)銑刀的可靠使用壽命并降低其壽命周期費(fèi)用消耗。該結(jié)論與仿真結(jié)果相一致，進(jìn)一步說(shuō)明了本文所提方法的優(yōu)越性。

4 結(jié)論

通過研究，建立了基于剩余壽命預(yù)測(cè)數(shù)據(jù)與不確定失效閾值的維修決策優(yōu)化模型，確定了不可維修設(shè)備的最優(yōu)替換時(shí)機(jī)。主要結(jié)論有：

1)設(shè)備失效閾值的不確定性真實(shí)存在，在設(shè)備剩余壽命預(yù)測(cè)研究中考慮不確定失效閾值具有合理性，能夠有效提升設(shè)備剩余壽命預(yù)測(cè)的準(zhǔn)確性。

2)針對(duì)不確定失效閾值分布系數(shù)的估計(jì)問題，EM算法較MLE算法在中、小樣本條件下適用性更強(qiáng)且準(zhǔn)確性更高，更能滿足工程應(yīng)用要求。

3)基于剩余壽命預(yù)測(cè)數(shù)據(jù)與不確定失效閾值的維修決策模型可以實(shí)現(xiàn)設(shè)備的最優(yōu)替換策略。同不考慮不確定失效閾值的決策結(jié)果相比，考慮不確定失效閾值的可靠使用壽命顯著延長(zhǎng)，設(shè)備壽命周期費(fèi)用率明顯降低。

在研究中，不確定失效閾值分布下限值κ多基于工程經(jīng)驗(yàn)人為給定，這可能引入主觀誤差，不利于實(shí)現(xiàn)科學(xué)決策。因此，未來(lái)應(yīng)著重針對(duì)κ的估計(jì)方法展開研究，以進(jìn)一步提升維修決策的準(zhǔn)確性。