基于比例故障率模型的柴油發(fā)動機(jī)視情維修決策

王海龍， 張 營， 左付山

(南京林業(yè)大學(xué)汽車與交通工程學(xué)院， 南京 210037)

柴油發(fā)動機(jī)作為汽車的動力部件，其結(jié)構(gòu)復(fù)雜，運(yùn)動部件多的特點(diǎn)，使得柴油機(jī)維修診斷較為困難。目前，對于車輛柴油機(jī)的維修主要采用定時(shí)維修和事后維修的策略，定時(shí)維修具有容易執(zhí)行等優(yōu)點(diǎn)，然而并沒有考慮發(fā)動機(jī)當(dāng)前運(yùn)行狀態(tài)；事后維修則是在故障發(fā)生之后再進(jìn)行相應(yīng)的維修，維修成本高，經(jīng)濟(jì)損失大。為了提高經(jīng)濟(jì)效益和社會效益，避免維修不足或者維修過度，相關(guān)學(xué)者提出了視情維修。視情維修是通過監(jiān)測設(shè)備的狀態(tài)參數(shù)來反映設(shè)備的性能狀態(tài)，然后根據(jù)相應(yīng)的狀態(tài)在故障發(fā)生之前采取一定的維修措施[1]，其最主要的就是建立設(shè)備運(yùn)行狀態(tài)參數(shù)與設(shè)備完好程度之間的對應(yīng)關(guān)系[2]，建立模型的方法包括基于設(shè)備可靠性、概率性以及卡爾曼濾波等建模方法。

董思辰等[3]選擇振動信號中的冗余參數(shù)作為協(xié)變量，利用威布爾比例故障率模型(Weibull proportional hazards model，WPHM)進(jìn)行齒輪箱的運(yùn)行可靠性評估，得到可靠度衰減時(shí)間與振動信號故障時(shí)間相近。劉璐[4]選擇航空機(jī)載設(shè)備電機(jī)軸承的振動信號，通過特征提取后，基于WPHM進(jìn)行了其可靠性評估和壽命預(yù)測。蔣文博等[5]提出了一種基于比例風(fēng)險(xiǎn)模型與機(jī)器學(xué)習(xí)的混合方法，有效的預(yù)測了電梯的剩余壽命。通過上述的研究發(fā)現(xiàn)，對于機(jī)械設(shè)備來說，WPHM能綜合考慮設(shè)備運(yùn)行狀態(tài)及工作時(shí)間，有效反映設(shè)備的狀態(tài)好壞，從而實(shí)現(xiàn)壽命預(yù)測與視情維護(hù)，因此柴油機(jī)作為機(jī)械設(shè)備，狀態(tài)參數(shù)復(fù)雜并且維修決策較為困難，通過建立WPHM進(jìn)行視情維修將減少大量維修成本并且提高工作效率。

綜上所述，構(gòu)建一種基于威布爾比例故障率模型的發(fā)動機(jī)視情維修策略，首先對監(jiān)測到的柴油發(fā)動機(jī)在各個(gè)階段的潤滑油金屬含量數(shù)據(jù)進(jìn)行了主成分分析(principal component analysis，PCA)數(shù)據(jù)降維并通過極大似然估計(jì)求出了所建立的模型參數(shù)；為綜合考慮成本和可用度，采用最小維修成本和最大可用度的綜合維修策略，求得相應(yīng)的最優(yōu)預(yù)防時(shí)間間隔以及失效率閾值，從而根據(jù)當(dāng)前監(jiān)測到柴油機(jī)金屬含量數(shù)據(jù)判斷是否需要對其進(jìn)行預(yù)防性維修。

1 威布爾比例故障率模型

1.1 比例故障率模型

比例故障率模型前期主要應(yīng)用于醫(yī)學(xué)領(lǐng)域[6]，其綜合考慮了設(shè)備運(yùn)行時(shí)間和當(dāng)前狀態(tài)信息之間的關(guān)系[7]，從而實(shí)現(xiàn)了對于設(shè)備健康狀況的評估。該模型的性質(zhì)是不同個(gè)體的故障率函數(shù)成比例[8]，并且故障率函數(shù)不隨時(shí)間t的變化而變化。其基本形式為

h(t;Z)=h0(t)exp(rZ)

(1)

式(1)中：h(t;Z)為故障率函數(shù)；h0(t)為基本故障率函數(shù)，只與時(shí)間有關(guān)[9]；rZ=r1Z1+r2Z2+…+rpZp，Z為協(xié)變量，反映裝備運(yùn)行狀態(tài)，按作用可分為外部協(xié)變量與內(nèi)部協(xié)變量[10]，如振動幅度、金屬含量等，r是協(xié)變量系數(shù)，反映了與之對應(yīng)的協(xié)變量對裝備故障率影響的嚴(yán)重程度[11]，ri的值越大，表明對應(yīng)的協(xié)變量對故障率函數(shù)影響越大。

1.2 威布爾比例故障率模型的建立

在汽車可靠性研究中，指數(shù)分布、威布爾分布、正態(tài)分布、對數(shù)正態(tài)分布等都是較為常用的理論分布類型。由于柴油發(fā)動機(jī)為機(jī)械部件，一般機(jī)械部件的故障數(shù)據(jù)分布規(guī)律均服從威布爾分布，同時(shí)其對各種類型的數(shù)據(jù)具有很強(qiáng)的擬合能力，所以應(yīng)用廣泛，因此上述基本故障率函數(shù)選擇較為經(jīng)典和簡單的兩參數(shù)威布爾分布故障率函數(shù)作為比例故障率模型的基礎(chǔ)函數(shù)[12]，兩參數(shù)威布爾故障率函數(shù)具體形式為

(2)

式(2)中：β為形狀參數(shù)，可改變曲線的陡峭程度；η為尺度參數(shù)[13]，可改變曲線中心線位置。

由此，得到了威布爾比例故障率函數(shù)[14]為

(3)

故障概率密度函數(shù)f(t)、可靠度函數(shù)R(t)以及故障率函數(shù)h(t)定義如下：f(t)是指所有部件在時(shí)間間隔(t，t+Δt)內(nèi)發(fā)生故障的個(gè)數(shù)占總數(shù)的比例；R(t)是指規(guī)定時(shí)間內(nèi)完成規(guī)定功能的概率，即部件壽命T大于某一時(shí)刻t的概率；h(t)是指定義為時(shí)間(0，t)內(nèi)不發(fā)生故障的條件下，下一個(gè)單位時(shí)間(t+Δt)內(nèi)發(fā)生故障的概率。由此，可以得到[15]

(4)

(5)

(6)

式中：N為部件總個(gè)數(shù)；n(t)為到t時(shí)刻發(fā)生故障的個(gè)數(shù)。

對式(6)左右兩邊同時(shí)在0～t上積分，得可靠度函數(shù)R(t)，即

(7)

(8)

同理可得，相應(yīng)的威布爾比例故障率可靠度函數(shù)R(t;Z)為

(9)

1.3 模型參數(shù)求解

為了通過威布爾比例故障率模型進(jìn)行柴油發(fā)動機(jī)的可靠性評估，最關(guān)鍵的是對于模型中的三個(gè)參數(shù)進(jìn)行求解β、η、r。對于三個(gè)參數(shù)估計(jì)的方法主要有最大似然估計(jì)法、最小二乘法、矩估計(jì)法、區(qū)間估計(jì)法等。由于極大似然估計(jì)法對于設(shè)備的不完全壽命數(shù)據(jù)是處理具有獨(dú)特的優(yōu)勢，使得模型更加精確，因此本文中選擇極大似然估計(jì)法對于模型的三個(gè)參數(shù)進(jìn)行估計(jì)。

極大似然估計(jì)的基本原理是通過估計(jì)的參數(shù)使得樣本出現(xiàn)的概率更大，似然函數(shù)的一般形式為

(10)

針對威布爾比例故障率模型，需要通過極大似然法求解，即求參數(shù)θ=(β，η，r1，r2，…，rp)何值時(shí)，使得似然函數(shù)的值最大。根據(jù)柴油發(fā)動機(jī)的失效密度函數(shù)f(t;Z)=h(t;Z)R(t;Z)，威布爾比例故障率模型似然函數(shù)為

exp(rZ)dt]}

(11)

式(11)中：F為所測數(shù)據(jù)中故障集；q為故障數(shù)據(jù)個(gè)數(shù)；N為數(shù)據(jù)總集。

為了求得似然函數(shù)的最大值，需要先對其進(jìn)行取對數(shù)，得到對數(shù)似然函數(shù)，然后分別對于β、η、r進(jìn)行一階求偏導(dǎo)并令其等于0，可以得到對數(shù)似然函數(shù)方程組。似然函數(shù)及其一階方程組[16]為

(12)

(13)

(14)

對于上述對數(shù)似然函數(shù)方程組求解方法主要有牛頓迭代法、遺傳算法[17]、粒子群法。

牛頓迭代法主要是通過泰勒一階展開，從初始值開始，分別對下一個(gè)值進(jìn)行泰勒展開，不斷循環(huán)，直到達(dá)到理想精度或最大迭代次數(shù)停止計(jì)算[18]，從而求得模型參數(shù)β、η、r。

原理為

f(x)=f(x0)+(x-x0)f′(x0)

(15)

步驟為

(16)

具體求解：由于上述得到了似然函數(shù)一階方程組，所以要求得似然函數(shù)方程組的解，還需要對于其進(jìn)行二次求導(dǎo)[19]，即

(17)

根據(jù)牛頓迭代法可得

(18)

式(18)中：β0、η0、r0為初始值。

通過牛頓迭代法求解得到3個(gè)待估參數(shù)之后，便得到了威布爾比例故障率模型以及相應(yīng)的可靠度模型，繼而就可以選擇最小成本維修策略或者最大可用度維修策略建立決策模型，求取最優(yōu)預(yù)防時(shí)間間隔。

2 綜合維修優(yōu)化模型

2.1 最小維修成本

針對汽車的柴油發(fā)動機(jī)，由于經(jīng)濟(jì)性在維修中的是首先被考慮的，其故障后維修費(fèi)用遠(yuǎn)大于預(yù)防性維修費(fèi)用，所以選擇最小維修成本法作為決策目標(biāo)，其具體模型求解如下所示。

當(dāng)以單位時(shí)間內(nèi)柴油發(fā)動機(jī)維修成本最小為目標(biāo)，制定相應(yīng)地最優(yōu)維修策略。單位時(shí)間維修成本為E，其計(jì)算公式[20-22]為

(19)

式(19)中：CPM為部件的預(yù)防性維修成本；CCM為部件故障的時(shí)候維修成本；N(t)=(t/η)β為0～t內(nèi)部件的平均失效數(shù)。

為了求得E的最小值，對t進(jìn)行一階求導(dǎo)，并等于0，得到最優(yōu)預(yù)防時(shí)間間隔T為

(20)

2.2 最大可用度法

可用度A(t)表示t時(shí)刻系統(tǒng)處于正常狀態(tài)的概率，其計(jì)算公式[23]為

(21)

(22)

求A(t)的最大值，即轉(zhuǎn)化為求t取何值時(shí)，α可以取到最小值。

2.3 綜合維修優(yōu)化模型

由于單方面只考慮減少柴油發(fā)動機(jī)維修成本，有可能會造成可用度無法保障，為了綜合考慮柴油發(fā)動機(jī)的維修成本最小以及可用度最大，為此提出了以可用度為約束，以單位維修最少為目標(biāo)，建立柴油發(fā)動機(jī)綜合維修優(yōu)化模型[24]，即

(23)

2.4 維修決策模型

根據(jù)綜合維修優(yōu)化模型得到最優(yōu)預(yù)防維修時(shí)間間隔T，在故障數(shù)據(jù)中找出一組與之最接近的失效時(shí)間與伴隨時(shí)間變量數(shù)據(jù)，帶入威布爾比例故障率模型，得到失效率閾值下限h。在部件工作正常時(shí)，故障率低于h，即

(24)

然后根據(jù)監(jiān)測的周期Δt,繼續(xù)代入模型計(jì)算，得到失效率模型閾值上限h*，當(dāng)設(shè)備當(dāng)前狀態(tài)及壽命數(shù)據(jù)代入模型后大于h*,表明設(shè)備已經(jīng)處于維修區(qū)，即

(25)

將式(25)經(jīng)過取對數(shù)等變形后，得到了對數(shù)風(fēng)險(xiǎn)度rZ，并將維修區(qū)域分為正常工作區(qū)、加強(qiáng)監(jiān)測區(qū)、維修區(qū)[25-26]，即

(26)

基于綜合優(yōu)化維修決策流程如圖1所示。

3 實(shí)例分析

3.1 模型選擇

對于柴油發(fā)動機(jī)的壽命數(shù)據(jù)，它有可能服從幾種分布類型，如指數(shù)分布、對數(shù)正態(tài)分布、威布爾分布等[27]。為了判斷壽命數(shù)據(jù)最符合哪種分布類型，首先需要通過最大似然法估計(jì)出各種分布的參數(shù)值，然后通過擬合優(yōu)度檢驗(yàn)方法進(jìn)行擬合優(yōu)度檢驗(yàn)[28]。然而，有時(shí)經(jīng)過分布的擬合優(yōu)度檢驗(yàn)判斷出符合條件的分布類型不止一種，所以引入了赤池信息準(zhǔn)則(akaike information criterion，AIC)進(jìn)行分布模型的最終選擇[29]。AIC信息準(zhǔn)則是基于最大熵原理，能夠檢驗(yàn)出不同模型間差異的顯著性，可通過計(jì)算錯(cuò)判概率進(jìn)行模型選擇分析，其計(jì)算公式為

(27)

為了驗(yàn)證上文中選擇威布爾模型的有效性，基于柴油發(fā)動機(jī)壽命數(shù)據(jù)對于對數(shù)分布、對數(shù)正態(tài)分布以及威布爾分布的AIC值進(jìn)行了計(jì)算與判定。柴油發(fā)動機(jī)的故障壽命數(shù)據(jù)如表1所示。

表1 故障壽命數(shù)據(jù)

首先采用極大似然估計(jì)的方法得到了3種分布的參數(shù)值，繼而選擇擬合優(yōu)度檢驗(yàn)方法中的K-S檢驗(yàn)[30]對于3種分布的p值和H值進(jìn)行了求解，求解結(jié)果如表2所示。

從表2可以看出，3種分布的K-S檢驗(yàn)的H值均為0，表示接受原假設(shè)，檢驗(yàn)結(jié)果均服從原假設(shè)分布。為了進(jìn)一步確定最佳分布，下面根據(jù)AIC公式計(jì)算得到了3種分布的AIC值如意表3所示。

表2 參數(shù)估計(jì)值和K-S檢驗(yàn)結(jié)果

表3 3種分布類型AIC值

根據(jù)AIC準(zhǔn)則，AIC值最小的分布類型作為最佳的分布類型[31],所以推斷出威布爾分布為故障數(shù)據(jù)的最佳分布類型。

3.2 模型的建立與求解

在不同階段，柴油發(fā)動機(jī)的運(yùn)行狀態(tài)與柴油發(fā)動機(jī)的潤滑油中的各種金屬含量有著直接關(guān)系，因此本文選擇其中鐵、鈣、鎂、鉛在潤滑油中的含量作為模型的協(xié)變量。在監(jiān)測的數(shù)據(jù)中，包括柴油發(fā)動機(jī)工作時(shí)間、換油時(shí)間、4種金屬含量等數(shù)據(jù)。其部分?jǐn)?shù)據(jù)如表4所示。

通過表4可以發(fā)現(xiàn)，在同一時(shí)刻，潤滑油中4種金屬含量數(shù)據(jù)量級差距較大，為了避免較大數(shù)據(jù)淹沒其他數(shù)，提高模型精度，需要先將原始數(shù)據(jù)歸一化。同時(shí)，為了減少協(xié)變量的個(gè)數(shù)，并得到一組更加全面的有效的特征協(xié)變量，可以對于數(shù)據(jù)進(jìn)一步的處理，即通過主成分分析得到新的協(xié)變量[32-33]，在MATLAB種將主成分貢獻(xiàn)度設(shè)置為90%，便得到了三個(gè)相互獨(dú)立的新的特征協(xié)變量Z1、Z2、Z3，回歸變量的相關(guān)系數(shù)矩陣和相關(guān)系數(shù)矩陣的特征值與貢獻(xiàn)率分別如表5、表6所示。

表5 回歸變量的相關(guān)系數(shù)矩陣

表6 相關(guān)系數(shù)矩陣的特征值與貢獻(xiàn)率

將處理后的監(jiān)測數(shù)據(jù)帶入對數(shù)似然函數(shù)中，并使用牛頓迭代法借助MATLAB對于模型的3個(gè)參數(shù)進(jìn)行求解，得到參數(shù)β=2.25，η=14 846，r1=-2.883，r2=0.741，r3=-1.633 ,因此威布爾比例故障率模型為

(28)

3.3 維修閾值確定

以綜合維修策略作為決策方法，通過調(diào)查及歷史經(jīng)驗(yàn)，柴油發(fā)動機(jī)的預(yù)防性維修費(fèi)用CPM平均為100元，故障后維修費(fèi)用CCM平均為2 000元，因此CPM/CCM=0.05,將上述結(jié)果代入式(19)中，畫出時(shí)間t與單位時(shí)間維修成本E的曲線如圖2所示，并求出相應(yīng)的最優(yōu)預(yù)防時(shí)間間隔。

圖2 平均時(shí)間維修成本曲線

通過圖2以及相應(yīng)的計(jì)算可知，當(dāng)T=3 486時(shí)，單位時(shí)間維修成本可以達(dá)到最小值，并滿足可用度要求，因此T=3 486為最優(yōu)預(yù)防時(shí)間間隔。然后，在歷史監(jiān)測的壽命與油液金屬含量數(shù)據(jù)中，找到一組與壽命時(shí)間最近的一組數(shù)據(jù)及其對應(yīng)的3個(gè)新的協(xié)變量，將其帶入WPHM模型后，可以得到相應(yīng)的故障率閾值h=3.32×10-5，再將h代入式(26)中，即可得到?jīng)Q策閾值下限曲線為

0.741×(-0.075 2)-1.636×

(-0.083 3)]=3.32×10-5

(29)

為了進(jìn)一步確定維修決策閾值上限，設(shè)置監(jiān)測時(shí)間間隔時(shí)間為1 260 h，即得到時(shí)間T=4 746，同樣找到一組與壽命時(shí)間最近的一組數(shù)據(jù)及其對應(yīng)的3個(gè)新的協(xié)變量，代入模型后便可得到對應(yīng)的故障率閾值h*=6.75×10-5，以及決策閾值上限曲線。決策閾值上下限曲線如圖3所示。

圖3 決策閾值上下限曲線

在后續(xù)的維修決策中，將監(jiān)測到的柴油發(fā)動機(jī)壽命時(shí)間和潤滑油金屬含量數(shù)據(jù)經(jīng)過數(shù)據(jù)處理后，帶入決策模型中，如果得到的數(shù)值處于決策曲線下方，表明工作正常，繼續(xù)進(jìn)行監(jiān)測；如果得到的數(shù)值處于上限與下限之間，則表明需要加強(qiáng)監(jiān)測，縮短監(jiān)測的時(shí)間間隔；如果求得的數(shù)值處于決策曲線的上方，則需要對于柴油發(fā)動機(jī)采取相應(yīng)的維修措施。

3.4 模型驗(yàn)證

為了驗(yàn)證上述建立的視情決策維修模型的有效性與可用性，選擇監(jiān)測到的一組同類型的柴油發(fā)動機(jī)的潤滑油中金屬含量數(shù)據(jù)以及對應(yīng)的運(yùn)行狀況進(jìn)行了分析驗(yàn)證。該發(fā)動機(jī)從1994年1月開始運(yùn)行，在運(yùn)行到1994年12月出現(xiàn)了故障，總共運(yùn)行了7 523 h。將柴油發(fā)動機(jī)的金屬含量數(shù)據(jù)經(jīng)過PCA后，帶入維修決策模型，得到相應(yīng)的對數(shù)風(fēng)險(xiǎn)值并在圖中畫出，如圖4所示。

圖4 柴油發(fā)動機(jī)狀態(tài)維修曲線圖

從圖4可以看出，該柴油發(fā)動機(jī)在運(yùn)行到2 516 h的時(shí)候，其最接近決策曲線下限，這時(shí)需要加強(qiáng)監(jiān)測；在運(yùn)行6 196 h后，其狀態(tài)接近決策曲線上限，此時(shí)應(yīng)當(dāng)采取一定的預(yù)防性維修措施。根據(jù)該發(fā)動機(jī)的實(shí)際狀態(tài)，由于沒有采取相應(yīng)的維修措施，其在運(yùn)行到7 523 h發(fā)生了故障，該實(shí)驗(yàn)結(jié)果表明本文提出的決策方法有效合理，可用于實(shí)際的該型發(fā)動機(jī)視情維修決策。

4 結(jié)論

基于收集的柴油發(fā)動機(jī)壽命數(shù)據(jù)及其對應(yīng)的潤滑油金屬含量數(shù)據(jù)，在對其進(jìn)行PCA后，選取了累計(jì)貢獻(xiàn)率達(dá)到90%的3個(gè)特征協(xié)變量，使用最大似然估計(jì)的方法估計(jì)得到了威布爾比例危險(xiǎn)率模型參數(shù)，建立了威布爾比例故障率模型；進(jìn)而以最小維修成本為目標(biāo)，以最大可用度作為約束，計(jì)算求得了最優(yōu)預(yù)防時(shí)間間隔，建立了同類型柴油發(fā)動機(jī)維修決策模型，為后續(xù)的狀態(tài)維修提供了有效的決策方法，不僅最大可能地降低了維修成本，又保證了柴油發(fā)動機(jī)的安全運(yùn)行和最大可用時(shí)間。最后以同類型發(fā)動機(jī)的狀態(tài)及壽命數(shù)據(jù)為例，驗(yàn)證了該模型可以有效地通過壽命和金屬含量數(shù)據(jù)判斷柴油發(fā)動機(jī)是否需要進(jìn)行維修，具有較大的應(yīng)用于視情維修決策的實(shí)用價(jià)值。