多元退化系統維修與備件訂購策略優化模型

楊志遠， 趙建民， 程中華， 郭馳名， 李俐瑩

(1. 陸軍工程大學石家莊校區 裝備指揮與管理系，石家莊 050003；2. 河北科技大學 信息科學與工程學院， 石家莊 050000)

現代復雜系統通常裝備有狀態監測單元如傳感器、儀表、計算設備[1]，這為基于狀態的維修(CBM)實施提供了系統狀態數據.相對于基于時間的維修方式，CBM能夠有效控制“維修不足”和“維修過度”的問題[2].通常來說，獲取系統狀態信息的方式可分為定期檢測和連續監測兩類，相比于定期檢測，連續監測能夠實時獲取系統狀態信息，從而更為及時有效地進行CBM決策，降低系統故障風險和維修費用.隨著傳感器技術的發展，連續監測的應用場景越來越多，特別是對于一些關鍵部件或系統，如發電機組，其本身價值高、停機損失大，對系統狀態進行連續監測更為適用.同時，復雜系統通常存在多個退化過程，且由于共同的運行環境、應力等因素，退化過程間往往是相關的，在CBM決策中需要充分考慮這些因素[3]，這也增加了多元退化系統CBM決策的難度.

對于有多個相關退化過程的系統而言，CBM決策的基礎是對不同退化過程間的相關性進行建模描述.目前，常用的方法包括多維退化模型[4-5]、退化率相關模型[6-7]和基于Copula函數的退化相關模型[8-9].文獻[4-5]分別采用多元正態分布和二元Birnbaum-Saunders(B-S)分布來描述不同退化過程間的相關關系，而分布形式也限制了此類方法的應用.文獻[6-7]將不同退化過程間的相關關系表示為一個部件退化水平對其余部件退化速率的影響.在工程實際中，這種影響關系及程度很難采用定量化方法描述和確定.相對而言，Copula函數在定量描述隨機變量相關性方面具有很高的靈活性，并且提供了分析相關結構對系統可靠性影響的參數化方法.基于不同種類的Copula相關結構，文獻[10]利用仿真方法對定期檢測策略下的系統維修費用進行分析，其中預防性維修閾值為固定值.文獻[11]建立了多部件系統檢測維修策略優化模型，其中部件故障間的相關性采用Gaussian Copula函數來描述.文獻[12]采用Levy Copula函數描述系統二元退化間的相關關系，并在此基礎上提出系統定期檢測與機會維修策略，利用仿真方法獲得了最優預防性維修閾值和機會維修概率閾值.在現有的基于Copula函數建立多元退化模型的文獻中，很少研究在連續監測條件下的CBM決策優化問題，且沒有考慮系統備件約束，而對一些價值較高的重要系統而言，很難保證時刻有可用備件，通常需要提前預定.在連續監測條件下，文獻[13]研究了單一退化過程的預防性更換閾值優化問題.在此基礎上，文獻[14]考慮二元相關Levy退化過程，以系統長期運行條件下的期望維修費用率和可用度為準則，計算獲得在不同系統結構下的最優預防性更換閾值.上述研究均是在維修延遲(維修資源的到達需要一定時間，此處的維修資源也可指備件)條件下，對預防性更換閾值進行決策，沒有綜合考慮備件訂購閾值以及預防性維修閾值的優化問題.

針對上述問題，本文考慮系統不同退化過程間具有的相關性，利用Copula函數建立系統多元退化模型.在此基礎上，提出系統維修及備件訂購策略，分析直接維修費用、備件庫存費用和故障停機損失，建立系統維修費用模型.以系統期望維修費用率準則，引入人工蜂群(ABC)算法搜索最優維修及備件訂購策略，從而提高維修工作的經濟性，為科學制定系統維修及備件訂購計劃提供方法支撐.

1 系統描述

系統性能會隨著工作時間的增加不斷退化，當退化超過閾值時會導致系統出現故障.假設系統存在n個相關性能退化過程，任意退化過程超過特定閾值時系統即會出現故障.

現有的退化過程模型主要包括退化軌跡模型、退化量分布模型和基于隨機過程的退化模型等.相對于前兩種模型，基于隨機過程的退化模型能夠描述性能退化在時間軸上的不確定性，更符合工程實際.在基于隨機過程的退化模型中，Gamma過程由于其良好的性質，廣泛應用于連續非減的退化過程建模，包括材料的磨損、侵蝕、裂紋、老化等.本文采用Gamma過程來建立系統退化過程模型.

令Xi(t),i=1,2,…,n表示系統第i個性能退化過程，根據Gamma過程平穩獨立增量性質可知，對于任意t>u>0，Xi(t)-Xi(u)=Xi(t-u)～Γ(αi(t-u),βi)，其中：Γ(·,·)為Gamma分布；αi>0和βi>0分別為形狀參數和尺度參數.假設系統初始退化量均為0，那么Xi(t)的概率密度函數和概率密度分布函數可表示為

(1)

如上所述，本文采用Copula函數描述退化過程的相關性.為便于實際應用，這里假設：在不同時間間隔內，不同退化過程的退化過程增量間的相關性可忽略.在文獻[8]中也采用了類似的假設分析多元Wiener退化過程的相關性.令t>u>0，根據Sklar定理，系統n個退化過程在時間區間[u,t]上退化增量Xi(t-u)的聯合概率分布函數可表示為

Ht-u(x)=C(FX1(t-u)(x1),FX2(t-u)(x2),…,

FXn(t-u)(xn);θ)

(2)

其對應的概率密度函數可表示為

c(FX1(t-u)(x1),FX2(t-u)(x2),…,

(3)

c(FX1(t-u)(x1),FX2(t-u)(x2),…,FXn(t-u)(xn);θ)=

(4)

式中：C為Copula函數；c為Copula密度函數；x=[x1x2…xn]，xi≥0；FXi(t-u)(xi)為第i個退化過程在[u,t]內增量的分布函數；θ=[θ1θ2…θm]為Copula函數中的參數向量，其取值直接影響著隨機變量間相關關系強弱.對于Copula函數參數估計問題，可基于退化數據通過最大似然估計或邊際推斷(IFM)法分兩步獲得.在此基礎上，可采用赤池信息準則(AIC)選擇合適的Copula函數描述系統退化相關性.關于Copula函數的參數估計方法參考文獻[15].

在此基礎上，假設系統退化故障閾值向量為QL=[QL1QL2…QLn]，可獲得系統可靠度函數為

Rsys(t)=P(X1(t)≤QL1,X2(t)≤

QL2,…,Xn(t)≤QLn)=

Ht(QL)=C(FX1(t)(QL1),

FX2(t)(QL2),…,FXn(t)(QLn);θ)

(5)

式中：P為概率.

由于Gamma函數和不完全Gamma函數的影響，FXi(t)(QLi)的計算較為復雜，可采用B-S分布對其進行近似計算[16].FXi(t)(QLi)可近似表示為

FXi(t)(QLi)≈1-FBS(t;QLi)=

(6)

Rsys(t)=

C(FX1(t)(QL1),FX2(t)(QL2),…,FXn(t)(QLn);θ)≈

C(Φ(U1(t)),Φ(U2(t)),…,Φ(Un(t));θ)

(7)

2 系統維修與備件訂購策略

本文的系統維修和備件訂購采用控制限策略(CLP)，即根據系統退化狀態確定相應的備件訂購和維修活動.備件訂購閾值為QA=[QA1QA2…QAn]，預防性更換閾值為QM=[QM1QM2…QMn].為建立系統維修與備件訂購決策優化模型，對相關策略做出如下假設.

(1) 對系統在運行過程中的退化狀態進行連續監測，系統各性能退化指標水平可以通過連續監測得到；單位時間狀態監測費用記為CI.

(2) 初始備件庫存為0，當系統任意退化過程達到備件訂購閾值時，系統發出訂貨信號，備件經過時間τ后到達.此時，若系統不進行更換，會產生相應的備件庫存費用；備件訂購費用記為CO，單位時間庫存費用記為CS.

(3) 當系統任意退化過程達到預防性更換閾值后，系統發出維修信號，如果此時庫存中有備件，直接對系統進行更換；否則，在備件到達后對系統進行更換.更換后，系統修復如新，且之后系統退化過程與更換前的系統無關.如果備件是在系統故障后到達，即會產生停機損失，系統單位時間內的停機損失記為CLO.一般來說，系統維修更換費用與其退化量相關，退化量高的系統通常所需維修費用也較多.系統退化量是由多個退化過程共同決定的，由于各退化過程的故障閾值和量綱可能會不同，進而將影響系統退化量的計算.因此，首先定義各退化過程的相對退化量為

(8)

由式(8)可以看出，當任意χi(t)≥1時，系統處于故障狀態.多退化過程引發的系統故障為競爭失效過程，系統狀態由多退化過程中的最大相對退化量決定.基于此，系統在t時刻的維修更換費用可表示為

CMAN=CR+ρEsup(χt)

(9)

式中：Esup(χt)=max{E(χ1(t)),E(χ2(t)),…,E(χn(t))}，E(χi(t))為χi(t)的期望；CR為系統更換固定費用；ρ為更換費用隨系統退化量變化系數.

(4) 維修時間相對于系統壽命周期來說非常短，在建模過程中忽略不計.

(5) 為使備件訂購閾值和系統預防性更換閾值有意義，令QAi≤QMi≤QLi,i=1,2,…,n.

以單退化過程為例，系統退化過程和更新過程如圖1所示.其中：tA、tM和tL分別為退化過程到達備件訂購閾值、預防性更換閾值及故障閾值的時間，且tL≥tM≥tA>0.

圖1 系統退化和更新過程

對于Gamma過程而言，其退化路徑為跳躍過程，因此有P(tM-tA=0)>0，P(tL-tA=0)>0，如圖1中系統第2個更新周期內的退化過程所示.事實上，對于很多左極右連的Levy過程(包含Gamma過程)來說，均存在上述性質.

系統不同維修時機如圖2所示.根據備件到達的時間不同，系統維修情況可劃分為以下3種：

圖2 系統的不同維修時機

(1) 當tA+τ

(2) 當tM≤tA+τ

(3) 當tA+τ≥tL時，在備件到達時對系統進行故障更換.

3 費用模型

本文以系統長期運行情況下單位時間平均費用(即期望費用率)為指標，對系統預防性維修和備件訂購策略進行評價.由前文分析可以看出，備件訂購閾值和預防性更換閾值的高低會影響不同維修方式發生概率、維修費用、庫存費用、更新周期長度，從而影響系統期望維修費用率.

3.1 期望維修費用率計算

考慮上述3種不同的系統維修情況，在不同情況下系統在一個更新周期內的維修費用C(T)及其相應的更新周期長度T可表示為

(1) 當tA+τ

(10)

(2) 當tM≤tA+τ

(11)

(3) 當tA+τ≥tL時，

(12)

在此基礎上，根據更新報酬原理可得系統的期望維修費用率為

γ∞=[CSE(max{tM-tA,τ})-

CLOE(min{tL-tA,τ})+(CLO-CS)τ+

ρ(Esup(χtA+τ)-Esup(χtM))P(tM-tA≤τ)+

ρEsup(χtM)+CR+CO]/

[E(max{tM-tA,τ})+E(tA)]+CI

(13)

3.2 具體表達式推導

τGM(τ)+τ-τGM(τ)+

(14)

(15)

P(tM-tA≤τ)=1-P(tM-tA>τ)=

(16)

由所建立的系統退化過程模型可得tA的生存函數為

P(tA>t)=P(X1(t)

QA2,…,Xn(t)

(17)

式中:Ht(QA)可由式(2)獲得.在此基礎上可得期望E(tA)的表達式為

(18)

在上述推導結果基礎上，當τ為常數時，可得系統期望維修費用率表達式為

(19)

(20)

(21)

對于Esup(χtM)和Esup(χtA+τ)，由定義可知：

E(Xi(tM))=αiβiE(tM)

(22)

式中：E(tM)與式(18)有類似表達式，此處不再贅述.在此基礎上，可得：

Esup(χtM)=

(23)

同理可得：

Esup(χtA+τ)=

(24)

將式(21)、(23)和(24)代入式(19)即可得到系統期望維修費用率的解析表達式.

3.3 期望維修費用率估算

通過前文獲得的系統期望維修費用率模型，可以得到相關的數值計算結果.但由于模型中涉及到x=[x1x2…xn]的高維積分，雖然在MATLAB等科學計算軟件中集成了高維積分的數值算法，可以直接使用，但算法的計算復雜度較高.為簡化計算，下面給出模型的近似表達式.

假設系統各退化過程在訂貨時間點tA的退化量已知，記為

XtA=[X1(tA)X2(tA) …Xn(tA)]，那么容易得到:

P(Xs

(25)

由于退化過程均為隨機過程，所以在實際中XtA為隨機變量.因此，這里用XtA的期望值E(XtA)=[E(X1(tA))E(X2(tA)) …E(Xn(tA))]代替XtA.由此可得：

Hs(QM-E(XtA))

(26)

同理可得：

(27)

(28)

4 優化模型及算法

在上述費用模型基礎上，為制定最優系統維修與備件訂購策略，以系統長期運行期望維修費用率最小為優化目標，對預防性維修閾值和備件訂購閾值進行優化.基于此，建立如下優化模型

(29)

s.t.0

在上述決策模型中，可以看出目標函數TC∞具有非線性、不可微的特點，且模型決策變量較多，當系統存在n個退化過程時，模型有2n個決策變量，因此，難以得到模型的解析解.

為解決這一問題，采用ABC算法搜索系統最優維修與備件訂購策略.ABC算法是一種基于蜂群搜索蜜源行為的群體智能優化算法，在智能優化算法中屬于比較新的算法，該算法控制參數少、易于實現、全局收斂性能好，并且其優化性能優于其他一些傳統智能優化算法[17]，如遺傳算法(GA)、粒子群優化(PSO)算法、差分進化(DE)算法等，特別是在非線性優化函數求解方面具有良好的性能.ABC算法已廣泛應用于解決各類優化問題，關于ABC算法的具體機理和優化過程可參考文獻[18]，這里不再贅述.在案例分析中，ABC算法采用MATLAB實現.

5 案例分析

為驗證上述決策優化模型的有效性，以兩部件系統為例，即n=2，對系統維修及備件訂購策略進行優化，并分析相關參數對模型優化結果的影響.考慮兩種不同的情形，在案例1中假設系統兩個退化過程參數相同，單位時間退化量均值為2，方差為4，故障閾值也相同；在案例2中假設系統兩個退化過程參數存在差異，單位時間內退化量均值分別為2和1.5，方差分別為4和1，相應的故障閾值也不同.不同案例的退化過程及故障閾值參數如表1所示.

表1 系統退化過程及故障閾值參數

對于模型中的費用參數，為方便說明，假定單位時間系統狀態監測費用為1，其余各項費用值均為與狀態監測費用的比值.系統維修費用參數值如表2所示.

表2 系統維修費用參數

在Copula函數選擇方面，以Gaussian Copula函數為例，描述系統兩個退化過程間的相關關系，二元Gaussian Copula函數的分布函數為

(30)

式中：Φ-1(·)為標準正態分布逆函數.當表示兩個變量間為正相關關系時，θ∈[0,1]，在本案例中令θ=0.7.當然，這里也可以采用其他類型的Copula函數描述退化間相關關系，其分析過程相同.

5.1 系統維修及備件訂購策略優化

在以上參數設置基礎上，假設備件訂購交付時間為τ=1.為驗證期望維修費用率式(13)的準確性，在不同的維修及備件訂購策略下采用Monte-Carlo(MC)方法對維修費用進行仿真，具體過程為：首先對系統退化過程進行仿真，得到備件訂購時間、系統更換、系統故障時間，進而計算系統庫存費用、更換費用及停機損失，最后得到費用具體值.在這里重復上述過程104次，得到系統維修費用的104個實現值.在此基礎上，求得平均值并與解析結果進行比較，如表3所示.

表3 解析結果與仿真結果對比

由表3可以看出，解析方法得到的期望費用率與仿真結果很接近，并且在MC仿真結果95%置信區間內，由此說明了上述系統期望維修費用率表達式的準確性.

下面分析當備件訂購閾值(QA1,QA2)為常數時，系統期望維修費用率γ∞隨預防性更換閾值(QM1,QM2)的變化情況，結果如圖3所示.

圖3 γ∞隨(QM1,QM2)變化情況

由圖3可以看出，無論對于案例1還是案例2，當(QA1,QA2)為常數時，隨著(QM1,QM2)的增大，γ∞總體上來說呈現先減小后增大的趨勢，這是由庫存費用、維修費用增加以及系統更新周期增長共同作用的結果.在給定的備件訂購閾值下，存在最優的(QM1,QM2)使得γ∞達到最小.同理，下面固定預防性更換閾值(QM1,QM2)，分析系統期望維修費用率γ∞隨備件訂購閾值(QA1,QA2)的變化情況，結果如圖4所示.

由圖4可以看出，無論對于案例1還是案例2，當(QM1,QM2)為常數時，隨著(QA1,QA2)的增大，γ∞總體上來說呈現先減小后增大的趨勢，這主要是受預防性更換概率減小、故障更換概率增大以及庫存費用減少的影響.同樣，在給定的預防性更換閾值下，存在最優的(QA1,QA2)使得γ∞達到最小.

圖4 γ∞隨(QA1,QA2)的變化情況

下面基于ABC算法對模型進行求解，得到最優維修及備件訂購策略.在ABC算法中，設置種群數量為10，適應度函數為1/γ∞(QA,QM).以案例1為例，考慮到算法中隨機因素的影響，運行ABC算法5次，獲得的系統維修及備件訂購策略優化過程如圖5所示，其中：K為迭代次數；上標*表示最優.

圖5 ABC算法優化過程

由圖5可以看出，當迭代次數達到100次后，最優期望費用率值不再變化，且各優化過程均收斂于相同的最優值，對案例2也有類似的結果.當退化過程增多時，可適當增加算法中種群數量的設置，以提高算法收斂速度及精度.在此基礎上，得到的系統最優維修及備件訂購策略如表4所示.此外，假如在這里采取文獻[14]中的策略，即不單獨考慮系統備件訂購閾值，而是在預防性更換閾值(QM1,QM2)處訂購備件，即QAi=QMi.在此策略下，所獲得的系統最優維修策略見表4.

表4 系統最優維修及備件訂購策略

圖6 不同(QA1,QA2)下和γ∞對比

5.2 靈敏度分析

圖7 備件交付時間對優化結果的影響

下面分析θ的變化對系統最優維修及備件訂購策略的影響，在案例1和案例2規定的情形下，計算結果如表5所示.

表5 不同θ下的系統最優維修及備件訂購策略

6 結語

多元退化系統CBM決策通常需要考慮不同退化過程間的相關性.在維修策略中同時考慮預防性更換閾值和備件訂購閾值，在推導系統退化狀態到達各閾值時間聯合分布基礎上，得到了連續監測條件下的系統期望維修費用率精確和近似表達式.在案例分析中，通過ABC算法得到了最優維修和備件訂購策略，驗證了模型和算法的有效性.通過和已有策略的對比表明了所提系統維修和備件訂購策略的優勢.同時，在案例中也檢驗了系統期望維修費用近似表達式的精度，分析了備件交付時間和退化相關程度參數對決策結果的影響.值得注意的是，本文采用Gamma過程建立退化過程模型，事實上，對于包含Gamma過程的Levy過程來說，上述模型也同樣適用.此外，對于一些非Copula相關的多元退化系統在連續監測下的維修閾值優化問題，可參考上述方法建立優化模型.