基于Gamma退化過(guò)程的剩余壽命預(yù)測(cè)及維修決策優(yōu)化模型研究

谷玉波，賈云獻(xiàn)，張英波

(軍械工程學(xué)院 裝備指揮與管理系，石家莊 050003)

機(jī)械設(shè)備在使用過(guò)程中，某些部件的運(yùn)行狀態(tài)會(huì)逐漸退化。狀態(tài)的退化可以表現(xiàn)為機(jī)械部件的磨損、裂紋的增長(zhǎng)、腐蝕程度的加深等，這些都是一系列理化效應(yīng)的結(jié)果。產(chǎn)品狀態(tài)退化的最終結(jié)果會(huì)導(dǎo)致功能故障的發(fā)生，而且在產(chǎn)品的退化過(guò)程中能夠直接表征退化狀態(tài)的狀態(tài)指標(biāo)(磨損、裂紋、腐蝕等)經(jīng)常是不可測(cè)量或難于直接進(jìn)行測(cè)量的。為獲得設(shè)備狀態(tài)退化程度的數(shù)據(jù)，可以運(yùn)用CBM狀態(tài)監(jiān)控技術(shù)對(duì)設(shè)備的狀態(tài)數(shù)據(jù)進(jìn)行提取分析，間接地反映出設(shè)備的健康狀態(tài)，實(shí)現(xiàn)對(duì)設(shè)備剩余壽命的預(yù)測(cè)，最終實(shí)現(xiàn)維修決策的優(yōu)化，使維修效果達(dá)到最佳。在此，運(yùn)用Gamma退化過(guò)程狀態(tài)空間模型來(lái)描述設(shè)備狀態(tài)性能的退化，通過(guò)該模型實(shí)現(xiàn)了剩余壽命的預(yù)測(cè)，并運(yùn)用剩余壽命預(yù)測(cè)結(jié)果對(duì)維修決策模型進(jìn)行優(yōu)化。

1 Gamma退化過(guò)程的狀態(tài)空間模型

在常用的分布中，Gamma分布可以用來(lái)描述連續(xù)累積磨損過(guò)程，并且具有退化過(guò)程建模所需的所有屬性，即非負(fù)、平穩(wěn)增長(zhǎng)和獨(dú)立增量過(guò)程[1]。因此，運(yùn)用Gamma分布建立狀態(tài)空間模型，可以描述設(shè)備的退化過(guò)程，并實(shí)現(xiàn)剩余壽命的預(yù)測(cè)。

1.1 狀態(tài)空間模型

狀態(tài)空間模型(State Space Model，SSM)是一種用來(lái)描述系統(tǒng)演化過(guò)程的模型，可以用來(lái)建立產(chǎn)品的狀態(tài)退化模型。狀態(tài)空間模型假設(shè)系統(tǒng)狀態(tài)隨時(shí)間的演化可以由一個(gè)不可觀測(cè)或難以觀測(cè)的狀態(tài)序列確定，與狀態(tài)序列伴隨的是一個(gè)可以觀測(cè)的測(cè)量序列，兩者之間的關(guān)系可以由狀態(tài)空間模型來(lái)描述，其中應(yīng)用較為普遍的狀態(tài)空間模型是典型相關(guān)(Canonical Correlation)方法[2]。通過(guò)建立狀態(tài)方程和觀測(cè)方程，狀態(tài)空間模型為描述系統(tǒng)的動(dòng)態(tài)演化特性建立了一個(gè)有效的模型框架，一些復(fù)雜的問(wèn)題就可以用比較簡(jiǎn)單的形式表示。狀態(tài)空間模型的基本方程為

xi+1=Fi(xi,ui,ζ,θ1)，

(1)

yi=Hi(xi,ui,ε,θ2)，

(2)

式中：Fi為狀態(tài)轉(zhuǎn)移函數(shù)，Hi為測(cè)量函數(shù)，一般而言Fi和Hi假設(shè)是已知的；ui為系統(tǒng)的輸入變量或控制變量；θ1,θ2分別為狀態(tài)方程和觀測(cè)方程的靜態(tài)參數(shù)；ζ和ε分別為系統(tǒng)噪聲和觀測(cè)噪聲，一般假設(shè)ζ和ε相互獨(dú)立，且服從某一特定分布。

1.2 Gamma過(guò)程

Gamma過(guò)程是一個(gè)具有獨(dú)立、非負(fù)增量的隨機(jī)過(guò)程。運(yùn)用Gamma過(guò)程對(duì)一些退化過(guò)程建立模型，能夠使建模過(guò)程中的數(shù)學(xué)計(jì)算更加明確易懂[3]，適合描述在時(shí)域上[4]單調(diào)累積微小增量的退化過(guò)程，比如磨損、腐蝕、裂紋增長(zhǎng)等。

根據(jù)Gamma過(guò)程，定義非負(fù)隨機(jī)過(guò)程{X(t);t≥0}，其概率密度函數(shù)為

(3)

2 Gamma-SSM退化模型

2.1 模型的建立

假設(shè)設(shè)備狀態(tài)x的退化符合Gamma過(guò)程，y表示與x相對(duì)應(yīng)的觀測(cè)變量，xi表示設(shè)備ti時(shí)刻所對(duì)應(yīng)狀態(tài)，則建立的狀態(tài)空間模型的狀態(tài)方程和觀測(cè)方程為

xi-xi-1～Gamma(α(ti)-α(ti-1),λ)，

(4)

yi=Hi(xi)+ε。

(5)

假設(shè)當(dāng)系統(tǒng)狀態(tài)x達(dá)到故障閾值Xf時(shí)系統(tǒng)發(fā)生故障，如圖1所示。則系統(tǒng)狀態(tài)從0首次達(dá)到故障閾值所需的時(shí)間為

圖1 系統(tǒng)退化過(guò)程

Tf=inf{t:x=Xf,t>0}。

(6)

為簡(jiǎn)化問(wèn)題，假設(shè)觀測(cè)噪聲是可以加的。假設(shè)設(shè)備狀態(tài)退化過(guò)程為平穩(wěn)Gamma過(guò)程，形狀參數(shù)α(t)為線性函數(shù)，即α(t)=a·t，并且假設(shè)觀測(cè)量y與狀態(tài)量x之間的關(guān)系為y=c·x+ε，則在獲得觀測(cè)序列y0:n的條件下，退化模型可以表示為

xi+1-xi～Gamma(α·(ti+1-ti),λ)，

(7)

yi=c·xi+ε。

(8)

由假設(shè)可知，觀測(cè)噪聲ε服從均值為0的正態(tài)分布。所以，y～N(c·x,σ)。因此只要獲得參數(shù)a,λ,c,σ的值，就可以確定該退化模型的形式。

2.2 模型參數(shù)估計(jì)

在對(duì)上述模型中參數(shù)進(jìn)行求解的過(guò)程中，會(huì)出現(xiàn)新的需要估計(jì)的參數(shù)，導(dǎo)致運(yùn)算量增大。為解決此問(wèn)題，運(yùn)用文獻(xiàn)[5]提到的經(jīng)驗(yàn)最大化(Experience Maximization，EM)算法與粒子濾波(Particle Filtering，PF)算法相結(jié)合的方法求解模型中出現(xiàn)的參數(shù)。

EM算法是進(jìn)行極大似然估計(jì)的一種有效方法，簡(jiǎn)單來(lái)說(shuō)，這一算法的每一次迭代主要包括2個(gè)步驟：(1)求期望(Expectation Step)，稱為E步；(2)求極大值(Maximization)，稱為M步。EM算法通過(guò)假設(shè)潛在變量的存在，來(lái)簡(jiǎn)化似然函數(shù)極大估計(jì)[3]。

PF算法是一種基于Monte Carlo方法和遞推Bayes估計(jì)的統(tǒng)計(jì)濾波方法[6]，建立在序貫重要性采樣(Sequence Importance Sampling，SIS)和Bayes理論的基礎(chǔ)上，對(duì)于解決非線性、非Gauss問(wèn)題有很好的效果[7]。PF算法的基本思想是：首先根據(jù)系統(tǒng)狀態(tài)向量的經(jīng)驗(yàn)分布在狀態(tài)空間產(chǎn)生一組隨機(jī)樣本集合(即粒子)，然后根據(jù)觀測(cè)數(shù)據(jù)不斷調(diào)整粒子的權(quán)重和位置，通過(guò)調(diào)整后的粒子信息修正最初的經(jīng)驗(yàn)分布。其實(shí)質(zhì)是用由粒子及權(quán)重組成的離散隨機(jī)測(cè)度近似相關(guān)的概率分布，并根據(jù)算法遞推更新離散隨機(jī)測(cè)度。當(dāng)粒子容量足夠大時(shí)，就近似于狀態(tài)變量真實(shí)的后驗(yàn)概率分布。

由于EM-PF算法兼有兩種算法各自的優(yōu)點(diǎn)，可以更加有效地對(duì)模型中的參數(shù)進(jìn)行求解。模型參數(shù)的求解過(guò)程如圖2所示。

圖2 模型參數(shù)求解過(guò)程

3 剩余壽命預(yù)測(cè)

剩余壽命是指產(chǎn)品從被檢測(cè)的某一時(shí)刻起到該產(chǎn)品發(fā)生故障的時(shí)間長(zhǎng)度[8]，設(shè)備剩余壽命的預(yù)測(cè)是維修管理中的一個(gè)重要環(huán)節(jié)，由于剩余壽命預(yù)測(cè)對(duì)產(chǎn)品使用過(guò)程中的安全性、經(jīng)濟(jì)性和任務(wù)性具有重要的影響，因此它是制定維修決策的重要依據(jù)。

通過(guò)確定狀態(tài)空間模型參數(shù)，就能得到Gamma狀態(tài)空間模型的具體形式，進(jìn)而得到產(chǎn)品的剩余壽命分布函數(shù)。令Ns為粒子個(gè)數(shù)；Tf為設(shè)備故障時(shí)間；wi為第i個(gè)粒子的權(quán)重；Xf為故障閾值；y0:c為直至當(dāng)前時(shí)刻所獲得的觀測(cè)序列；y0:i為當(dāng)前獲得的觀測(cè)值序列；ti為第i次狀態(tài)檢測(cè)時(shí)間；xi為ti時(shí)刻的狀態(tài)值；x1:n為n次檢測(cè)所獲得的狀態(tài)序列，即x1:n={x1,x1,x2,…,xn} ；yi為ti時(shí)刻的退化量觀測(cè)值；y0:n為n次檢測(cè)所獲得的觀測(cè)序列，即y0:n={y0,y1,y2,…,yn},則ti時(shí)刻的剩余壽命τi分布函數(shù)為

(9)

p(τi+ti≥Tf|xi)=p(x(τi+ti)≥

(10)

(11)

另外，預(yù)測(cè)tk時(shí)刻的狀態(tài)概率密度函數(shù)可表示為

(12)

(13)

(14)

(15)

4 維修決策優(yōu)化模型

構(gòu)建剩余壽命概率密度函數(shù)的最終目的是根據(jù)剩余壽命與實(shí)際需求做出維修決策，確定合理的維修或更換時(shí)間。平時(shí)對(duì)設(shè)備進(jìn)行維修管理的過(guò)程中，人們總是希望設(shè)備的維修管理費(fèi)用盡可能達(dá)到最低。以費(fèi)用最小為目標(biāo)的維修決策模型可以根據(jù)設(shè)備當(dāng)前的運(yùn)行狀態(tài)確定在下一間隔期內(nèi)是否進(jìn)行維修以及何時(shí)進(jìn)行維修，可最大程度地避免故障的發(fā)生。

在維修策略中，產(chǎn)品有2種維修方式：(1)修復(fù)性維修，即產(chǎn)品在達(dá)到預(yù)定維修時(shí)刻之前發(fā)生故障而進(jìn)行的維修或更換活動(dòng)，也稱故障維修或故障更換；(2)預(yù)防性維修，即產(chǎn)品在達(dá)到預(yù)定的維修時(shí)刻進(jìn)行的維修或更換活動(dòng)。因此按照維修決策的目標(biāo)，并根據(jù)剩余壽命概率密度函數(shù)，建立以最小費(fèi)用為目標(biāo)的維修決策模型，其決策過(guò)程如圖3所示。

圖3 以費(fèi)用最小為目標(biāo)的維修決策示意圖

4.1 假設(shè)和參數(shù)

(1)每隔Δt進(jìn)行一次狀態(tài)檢測(cè)并獲取當(dāng)前時(shí)刻的狀態(tài)數(shù)據(jù)。

(2)如果在檢測(cè)間隔期內(nèi)發(fā)生故障，則立刻進(jìn)行修復(fù)性維修。

(3)設(shè)備的預(yù)防性維修費(fèi)用小于故障后的修復(fù)性維修費(fèi)用，即Tp

(4)預(yù)防性維修和修復(fù)性維修可以使設(shè)備恢復(fù)初始狀態(tài)，即修復(fù)如新。

(5)參數(shù)：E(C)為更新周期內(nèi)的期望總費(fèi)用；E(T)為期望更新周期長(zhǎng)度；Δt為狀態(tài)檢測(cè)間隔期，Δt=ti-ti-1；f(τi|y0:i)為ti時(shí)刻狀態(tài)信息為y0:i時(shí)產(chǎn)品的剩余壽命概率密度分布函數(shù)；Pp為部件進(jìn)行預(yù)防性維修的概率；Pf為部件進(jìn)行修復(fù)性維修的概率；cp為預(yù)防性維修費(fèi)用；cf為修復(fù)性維修費(fèi)用；TR為最佳維修更換時(shí)間；C(TR)為部件維修更換時(shí)間為TR時(shí)，長(zhǎng)期使用下的單位時(shí)間費(fèi)用；ti為當(dāng)前時(shí)刻的狀態(tài)監(jiān)控點(diǎn)，t0=0；Tp為預(yù)防性維修所需的平均時(shí)間；Tc為修復(fù)性維修所需的平均時(shí)間；ci為每次實(shí)施狀態(tài)檢測(cè)的費(fèi)用。

4.2 以費(fèi)用為目標(biāo)的決策模型

設(shè)備部件在長(zhǎng)期使用下的單位時(shí)間費(fèi)用可表示為

(16)

E(C)=cfPf+cpPp+nci=cff(τi

ti|y0:i)+cp[1-f(τi

cp+(cf-cp)f(τi

(17)

(18)

(19)

將(17)～(18)式代入(16)式中，則費(fèi)用模型為

(20)

其中，nci為更新周期內(nèi)監(jiān)控所產(chǎn)生的費(fèi)用，由于狀態(tài)檢測(cè)中檢測(cè)費(fèi)用大多為設(shè)備投入費(fèi)用，nci常?？梢院雎圆挥?jì)，因此在實(shí)際決策中，近似取nci=0。如果Tp和Tc相對(duì)于TR來(lái)說(shuō)很小，在費(fèi)用模型中也可以將其忽略。

利用數(shù)值方法最小化(20)式，得出設(shè)備進(jìn)行基于狀態(tài)的維修時(shí)可達(dá)到的最小單位費(fèi)用及所對(duì)應(yīng)的維修時(shí)刻TR。當(dāng)獲取新的狀態(tài)信息值y0:i+1時(shí)，需要將數(shù)據(jù)代入(20)式重新計(jì)算更新結(jié)果。

4.3 最佳維修時(shí)間確定方法

通過(guò)費(fèi)用決策模型可以得到各檢測(cè)時(shí)刻對(duì)應(yīng)的維修時(shí)間TR，使得單位時(shí)間維修費(fèi)用最小。當(dāng)TR-ti>Δt時(shí)，不對(duì)部件進(jìn)行維修，直到下次正常檢測(cè)時(shí)刻；當(dāng)首次出現(xiàn)TR-ti≤Δt時(shí)，則到達(dá)TR時(shí)對(duì)部件進(jìn)行維修或更換，此時(shí)的TR即為基于狀態(tài)維修的最佳維修時(shí)間。

5 實(shí)例分析

試驗(yàn)采用3套R(shí)exnord ZA-2155雙列滾子軸承，對(duì)2套軸承進(jìn)行全壽命試驗(yàn)，用于模型參數(shù)估計(jì)和模型驗(yàn)證，對(duì)第3套軸承進(jìn)行截尾壽命試驗(yàn)，用于維修決策，分別記為軸承1、軸承2和軸承3。試驗(yàn)得到2組全壽命試驗(yàn)數(shù)據(jù)和1組截尾壽命數(shù)據(jù)。

全壽命試驗(yàn)過(guò)程中，軸承的轉(zhuǎn)速和載荷保持在2 000 r/min和10 000 LB(約44 500 N)，并采用PCB加速度傳感器對(duì)軸承的振動(dòng)信號(hào)進(jìn)行采集，每1 h采集1次信號(hào)，采樣頻率為20 kHz，長(zhǎng)度為1 s。

試驗(yàn)過(guò)程中，軸承1和軸承2分別在966 h,982 h出現(xiàn)強(qiáng)烈的噪聲，此時(shí)認(rèn)為軸承故障。用振動(dòng)信號(hào)的能量均值作為間接狀態(tài)觀測(cè)數(shù)據(jù)y，磨損量為狀態(tài)數(shù)據(jù)x，2組全壽命數(shù)據(jù)的振動(dòng)信號(hào)能量均值在整個(gè)壽命周期上的變化如圖4所示。

圖4 試驗(yàn)軸承平均能量的變化過(guò)程

應(yīng)用軸承1的間接狀態(tài)觀測(cè)數(shù)據(jù)對(duì)模型參數(shù)進(jìn)行估計(jì)，并采用EM-PF算法對(duì)剩余壽命模型中的參數(shù)進(jìn)行計(jì)算，結(jié)果見(jiàn)表1。

表1 模型參數(shù)計(jì)算結(jié)果

確定模型參數(shù)之后就能得到基于Gamma退化過(guò)程的狀態(tài)空間模型，通過(guò)計(jì)算求得軸承剩余壽命概率密度函數(shù)f(τi|y0:i)。對(duì)于軸承2的觀測(cè)數(shù)據(jù)，運(yùn)用Matlab 7.0對(duì)剩余壽命預(yù)測(cè)方法進(jìn)行實(shí)現(xiàn)，不同狀態(tài)檢測(cè)時(shí)刻(500～980 h)的剩余壽命概率密度函數(shù)如圖5所示，其中*表示剩余壽命的估計(jì)值，?表示剩余壽命真實(shí)值。剩余壽命預(yù)測(cè)值與實(shí)際值見(jiàn)表2。

圖5 不同檢測(cè)時(shí)刻剩余壽命概率密度

表2 剩余壽命預(yù)測(cè)值與實(shí)際值比較(軸承2) h

由上述結(jié)果可以看出，該模型對(duì)于此類試驗(yàn)背景下的軸承剩余壽命預(yù)測(cè)具有一定的實(shí)用性，下面利用該模型對(duì)軸承3進(jìn)行維修決策，確定使得單位時(shí)間費(fèi)用最小的最佳維修時(shí)間。假設(shè)cp=600元，cf=1 200元，在得到軸承的剩余壽命概率密度函數(shù)之后，將f(τi|y0:i)代入到(2)式中，計(jì)算各維修時(shí)刻的單位時(shí)間費(fèi)用，得到每次檢測(cè)時(shí)費(fèi)用最低的更換時(shí)間TR，如圖6所示。圖中可以看出，在某一檢測(cè)時(shí)刻，隨著更換時(shí)間的增加，單位時(shí)間費(fèi)用先減少后增加，在極值點(diǎn)處可以取得該次檢測(cè)的最佳維修時(shí)間。圖中?表示各檢測(cè)時(shí)刻使單位時(shí)間費(fèi)用最低的維修時(shí)間。

圖6 不同檢測(cè)時(shí)刻所確定的維修更換時(shí)間

結(jié)合最佳維修時(shí)間確定的方法，通過(guò)Matlab計(jì)算，可以得出首次滿足TR-ti≤Δt條件時(shí)，TR=701 h。所以在701 h對(duì)軸承進(jìn)行更換，可以使單位時(shí)間費(fèi)用C(TR)達(dá)到最低,約為289.6元。

6 結(jié)束語(yǔ)

為解決維修決策的實(shí)際問(wèn)題，重點(diǎn)研究了剩余壽命預(yù)測(cè)方法和維修決策優(yōu)化模型。首先根據(jù)產(chǎn)品狀態(tài)的退化過(guò)程的特征建立了基于Gamma 退化過(guò)程的狀態(tài)空間模型，并對(duì)模型參數(shù)的估計(jì)方法進(jìn)行了說(shuō)明，進(jìn)而通過(guò)該模型實(shí)現(xiàn)剩余壽命的預(yù)測(cè)，得到產(chǎn)品的剩余壽命概率密度函數(shù)。最后以滾子軸承壽命試驗(yàn)得到的數(shù)據(jù)為例，驗(yàn)證了剩余壽命預(yù)測(cè)模型，并以費(fèi)用最小為目標(biāo)建立了維修決策模型，對(duì)另一軸承進(jìn)行了維修決策的優(yōu)化，實(shí)現(xiàn)了剩余壽命預(yù)測(cè)與維修決策的有機(jī)結(jié)合，實(shí)例證明了所建立模型的有效性和可行性。