退化相關多部件系統預防性維修決策模型

楊志遠, 趙建民, 程中華

(陸軍工程大學石家莊校區裝備指揮與管理系, 河北 石家莊 050003)

0　引　言

隨著技術的發展,現代工業設備結構越來越復雜,相應的價值和停機造成的損失也越來越高,維修對于降低系統故障率和減少使用成本具有重要的意義。在復雜系統中,各部件間的退化失效往往具有相互影響,這就使得這類系統的維修決策變得更為復雜,如果將此類系統的維修決策問題按照獨立退化過程來處理,將會影響維修決策的科學性乃至正確性。

部件間的相關關系可劃分為經濟相關、結構相關和故障相關,近年來,在復雜系統維修決策研究中,最受關注的是對故障相關性的研究。故障相關也稱隨機相關,是指一個部件的狀態會影響其他部件的狀態。這里的“狀態”包括壽命、故障率、故障狀態等各種狀態度量標準。目前,在故障相關研究中最常見的是考慮故障率之間的相關關系,這種方法是建立在多部件系統故障失效隨機數據基礎上的。在關于故障相關的早期研究中,文獻[1-2]提出了多部件系統中存在的3種類型的故障相關;文獻[3]考慮兩部件并聯系統在交叉檢測下的可靠度問題,其中一個部件故障后將以概率p影響另一個部件的故障概率。文獻[4]建立了一個解析模型用于定量分析部件間的相關故障率。文獻[5]假設兩部件可修復系統中故障相關是單向的,在此基礎上,建立了周期性檢測優化模型。另外共因失效也是一種特殊的故障相關表現,目前關于共因失效的研究也可以認為是基于故障率相關的,共因失效增大了系統各失效模式間的聯合失效概率[6],關于其研究可參考文獻[7-8]。故障率交互的方法是基于故障失效數據的,在實際中,有些設備的可靠度較高因而能獲得的故障數據較少[9],會限制該方法的應用,而且由于復雜系統各個部件的獨立故障率難以通過失效數據統計得到,因此也很難推導故障率之間的相互影響關系。

故障相關性另一種研究方法是基于退化過程相關開展的,通過各部件退化過程間的相互影響來對部件間的相關性進行分析研究。目前關于退化過程相關的研究文獻還較少,比較常用的是用兩變量退化分布模型[10-12]和Copula函數[13-14]來表現復雜系統不同部件退化量間的相關關系。使用兩變量退化分布模型和Copula函數難點在于確定相關函數的形式,就Copula法來說,常見的Copula函數有Gumbel Copula、Clayton Copula、Frank Copula、Gaussian Copula、t-Copula等[15],確定系統適合的Copula函數及相關參數比較困難。退化過程相關也可用部件退化速率之間的關系表示,目前采用此種方法的研究還較少且相關維修策略較為單一,在維修決策模型中也通常只考慮兩部件相關,文獻[16]運用退化速率相關的方法研究了多部件系統部件退化間的相關關系模型并提出了相應的基于狀態維修策略,但模型采用數據驅動的方法,沒有建立關于部件退化過程相互影響的解析模型,也未得到相應的退化過程顯式模型,這限制了對該方法的深入研究和使用。

基于此,本文在非線性Wiener退化過程模型基礎上,考慮退化速率影響建立了退化相關多部件系統退化過程解析模型和可靠度模型,并提出了相應的預防性維修策略,以單位時間維修費用最小為目標建立了預防性維修決策優化模型,進而得到了系統最優預防性維修策略,以提高這類系統維修活動的經濟效益。

1　退化過程模型

在現有退化過程模型中,常用的包括退化軌跡模型、退化量分布模型、累積損傷模型、Gamma過程模型以及Wiener過程模型,其中Gamma過程模型與Wiener過程模型是常用的隨機過程模型。由于Wiener過程模型能夠描述多種典型產品的性能退化過程,并且具有良好的計算性質,因而是目前性能可靠性建模領域較常用的退化模型之一。本文中采用非線性Wiener過程模型來建立部件退化模型,令X(t)表示部件在時刻t的退化量,那么X(t)可表示為

(1)

式中,λ(t;?)是關于時間t的連續非減函數,表示退化過程的平均趨勢,?為影響退化過程的相關參數;B(t)為標準布朗運動,表示退化過程中的隨機性。可以看出,如果令λ(t;?)=λ,那么模型就變為了線性Wiener過程模型;如果令σB=0,那么模型將退化為真實退化路徑模型。為便于分析,令X(0)=0,即在開始時刻0部件的退化量為0,用X(t)表示部件的累計退化量,則有

(2)

在此基礎上,進一步設定

(3)

式中,μ(t;?)表示的是期望退化量隨時間的變化情況。由此,部件退化非線性Wiener過程模型可表示為

X(t)=μ(t;?)+σBB(t)

(4)

假設系統由n個部件組成,部件退化過程中存在相關關系。根據建立的退化模型,n個部件退化過程模型可表示為

X(t)=μ(t;?)+σBB(t)

(5)

式中,X(t),μ(t;?)和σB均為n維列向量,表示不同部件的退化模型。可以看出,退化量X(t)由兩部分組成,第二部分表示的是部件在退化過程中的隨機性因素,假設對于不同的部件隨機性因素是相互獨立的,所以在本文中采用μ(t;?)描述不同部件退化過程間的相互關系。

在系統運行過程中,部件的退化會對與其具有相關關系部件的退化速率產生影響,這是由于某些部件的退化造成了整個系統在非正常狀況下運行,這會影響相關部件的退化速率。以齒輪箱為例,軸承的磨損將會造成過量的振動,因此相關的齒輪退化速度也會加快,反之,齒輪的磨損也會導致軸承的加速退化。為方便實際應用,本文中采用退化量期望值的變化定義部件在t時刻的退化速率,表示為

(6)

基于此,借鑒文獻[16]中的退化率相關法和文獻[18]中相關故障率的概念,建立部件退化相關解析模型為

?)=Θn×nμ(t;?)+Q(t)

(7)

式中,Θn×n是部件間的相關系數矩陣,表示部件退化間相關關系的強弱;Q(t)表示部件的固有退化率。以上參數可由系統或類似系統歷史退化數據或實驗數據估計得到,以較小時間間隔Δt對各部件的退化狀態進行取樣,令ΔXt=Xt+Δt-Xt,根據式(5),可得到ΔXt=μ(t+Δt;?)-μ(t;?)+σBB(t+Δt)-σBB(t),由Wiener過程可得

(8)

基于此,可采用最大似然估計法對Θn×n和Σ的值進行估計;對于Q(t)可參照文獻[16]中所采用的回歸方法得到,這里不再贅述。

式(7)表示部件退化速率受到其自身和具有相關關系部件退化量的影響。為便于理解,將式(7)展開為

(9)

式中,θij為部件j累積退化量對部件i退化速率的影響,θij=0表示部件j對部件i無影響。如果θij為常數,式(7)可以看作一個常系數線性非齊次微分方程組,根據線性非齊次微分方程組解的存在與唯一定理,給定相應的初始條件,那么式(7)存在唯一解。通過求解式(7),就可以求得相應的μ(t;?),進而得到各部件退化過程模型X(t)。為求解式(7),引入線性微分算子L,L表示為

(10)

引入線性微分算子L后,式(7)可簡化表示為

L[μ(t;?)]=Q(t)

(11)

其對應的齊次常微分方程組形式為

L[μ(t;?)]=0

(12)

根據微分方程組解得結構性質,式(11)的解可表示為

(13)

根據Cayley-Hamilton定理[19]可以得到φ(t)的表達式為

(14)

式中,P0=E;Pj=(Θ-λjE)(Θ-λj-1E)…(Θ-λ1E);j=1,2,…,n。其中,λ1,λ2,…,λn是矩陣Θ的特征根(不一定相異),而函數γj(t)(j=1,2,…,n)是微分方程組式(15)初值問題的解。

(15)

(16)

由此可得部件退化相關模型式(7)的通解為

(17)

給定系統各部件的初始退化量μ(t0;?)=μ0,即可得到具有退化相關關系各部件的期望退化量函數為

(18)

在此基礎上,將式(18)代入式(5)中,可得部件退化過程模型為

(19)

2　可靠度模型

由于本文在退化建模中采用的是非線性Wiener過程模型,然而由于退化模型的復雜性,難以得到部件精確的故障分布函數。為此,考慮到在實際中部件退化過程一般是非減的,做出假設:對于任意t>s≥0,部件退化量X(t)

2.1　部件可靠度

根據式(19)可得到退化量在t時刻的分布Gt(x),定義

Gt(x)=Pr{X(t)≤x}

(20)

用T(x)表示退化量達到一定水平x時的隨機時間,定義

Fx(t)=Pr{T(x)≤t}

(21)

基于以上假設,容易得到當且僅當T(x)>t時,X(t)≤x,因此Gt(x)和Fx(t)之間存在關系為

Gt(x)=1-Fx(t)

(22)

假設退化過程的故障閾值為D,在本文中,部件故障僅由退化引起,不考慮偶然故障影響,也就是說,當退化水平X(t)超過閾值D時，部件出現故障。那么根據假設和可靠度的定義,部件的可靠度函數可表示為

RD(t)=Gt(D)=1-FD(t)

(23)

(24)

式中,Φ(·)表示標準正態分布。

Ψt(x;u)=Pr{ΔXt(u)≤x}=

(25)

2.2　系統可靠度

根據部件退化相關模型可以看出,部件退化間的相互關系是通過μ(t;?)來表示的。而退化模型中σBB(t)只與時間相關,不同部件之間σBB(t)是相互獨立的,當時間一定時,μ(t;?)是固定值,所以不同部件間是否發生故障(退化量是否超過故障閾值)可認為相互獨立。因此,系統可靠度函數可表示為

(26)

3　系統預防性維修模型

3.1　基本假設與預防性維修策略

在建立預防性維修模型之前,做出以下基本假設:

假設1任意部件故障都會造成系統故障,系統故障是由不同部件故障引起的競爭失效過程,系統一旦發生故障則會立即停止運行并被發現。

假設2部件的退化狀態只能通過檢測得到,且檢測是完善的,即通過檢測能夠精確地得到部件退化狀態。

假設3系統由1個核心部件和n-1個輔助部件組成,任何一個部件的故障都將導致系統出現功能故障,核心部件是指價值較高、功能重要且發生故障造成損失較大的部件。需要注意，本文中的“部件”既可以是單個部件,也可以指一個子系統,本文中所指的核心部件是符合上述特點,有必要進行定期檢測的部件組成的子系統。

圖1　系統預防性維修策略

在以上預防性維修策略下,由于部件退化間的相互關系,部件在每個間隔期內的退化過程有所差異。假設部件在第i個預防性維修間隔內t時刻期望退化量記為μi(t;?),為了便于以下模型推導,令每個預防性維修間隔期均從0時刻開始,即0≤t≤Ti。在任意間隔期內,部件退化間的相互關系模型是相同的,均由式(7)給出,所以式(16)給出的通解在每個預防性維修間隔期內均適用,只需改變相應的初始條件即可。對部件進行更換使得部件修復如新,即退化量變為0,因此,在第i個預防性維修間隔期內模型初始條件可表示為

向量中第一項是對于核心部件來說的,由于對其進行檢測不影響退化量,在第i個間隔期初始退化量等于第i-1個間隔期結束時的退化量,其中對核心部件退化模型存在μ1(0;?)=μ0(T0;?)=0。由此,可得到在第i個預防性維修間隔期內部件的退化過程模型為

(27)

在式(27)基礎上,根據式(23)可得到各部件在各預防性維修間隔期內的可靠度和故障分布函數。

已經給出了部件可靠度和故障分布函數,但對于核心部件而言,由于其采用功能檢測策略,所以其發生故障更新或預防性更換的概率不能由可靠度函數簡單得到,需要分情況進行討論。

當i=1時,即在第一個預防性維修間隔期內,發生故障更新或預防性更換的概率較為簡單,根據部件可靠度函數及功能檢測策略,其分別為FD,1(T1)和GT1,1(D)-GT1,1(xp)。符號Gt,i(x)=Pr{Xi(t)≤x},表示在第i個預防性維修間隔期內的Gt(x)函數,其他相關符號表示相同,下面不再贅述。

當i≥2時,核心部件首次更新發生在第i個預防性維修間隔期內t時刻前故障更新的概率為

HD,i(t)=Pr{Xi(0)

Pr{Xi(0)

(28)

式中,g0,i(x)=dG0,i(x)/dx。

核心部件在第i個預防性維修間隔期內t時刻前未出現部件更新的概率為

Pr{Xi(0)

Pr{Xi(0)

(29)

核心部件在第i次檢測時出現預防性更換更新的概率為

Pr{Xi(0)

Pr{Xi(0)

(30)

(31)

式中,EC表示每個更新周期內系統維修費用期望值；EL表示系統期望更新周期長度。

在預防性維修策略中,系統有兩種可能的更新方式,一種是由于任意部件故障引起的系統故障更換,稱為故障更新;另一種是由于核心部件在檢測時，退化量在預防性更換閾值xp和故障閾值D之間引起的系統預防性更換,稱為檢測更新。

3.2　系統維修總費用

(32)

當i≥2時,根據系統故障更新是否由核心部件(假設核心部件編號為1)發生故障引起,考慮兩種情況,和首次預防性維修間隔期內故障更新概率計算原理類似,可得在第i個預防性維修間隔期內,系統故障更新概率為

(33)

(34)

基于此,可得到系統更新周期內故障更新期望費用為

(35)

系統檢測更新是由核心部件預防性更換引起的,而核心部件第i次檢測時發生預防性更換的概率已在之前推導給出,如式(29)所示,綜合所有可能的檢測點(i=1,2,…)即可得到系統檢測更新發生的概率,在此基礎上,可得到系統更新周期內檢測更新期望費用為

(36)

綜合系統更新周期內故障更新和檢測更新期望費用,可得到每個更新周期內系統維修費用期望值

EC=ECfs+ECps

(37)

3.3　系統更新周期長度

在計算期望更新周期長度時,和維修費用模型類似,也可劃分為故障更新和檢測更新,更新概率在維修模型中已得到。在此基礎上,可得到故障更新期望周期長度為

(38)

檢測更新期望周期長度為

(39)

綜合系統更新周期內故障更新和檢測更新期望周期長度,可得到系統期望更新周期長度為

EL=ELfs+ELps

(40)

將式(37)和式(40)代入式(31)中即可得到系統單位時間維修費用Cavg(T1,T,xp)的表達式。

(41)

4　算例分析

假設系統由兩個部件(子系統)組成,不失一般性,令部件1為核心部件,部件2為輔助部件,兩部件退化間存在相關關系,任意部件故障都會引起系統失效,系統故障是競爭失效過程。為減少系統維修費用,采用所提預防性維修策略對系統進行維修,退化和維修模型相關參數設定如表1所示,為方便分析,費用單位統一設置為元,部件退化量單位為毫米(mm),時間單位為月,在以下敘述中不再重復標明。

基于退化過程相關模型和預防性維修策略,可利用非線性Wiener過程的獨立增量特性生成部件退化過程的仿真軌跡。當T1=25,T=10時,對兩部件的退化過程獨立仿真5次,得到的仿真軌跡如圖2所示。

表1　模型參數設置

圖2　預防性維修策略下部件仿真退化軌跡

可以看出,兩個部件的退化過程均呈現明顯的遞增趨勢,且兩部件退化間有明顯的相關關系,如隨著部件1的退化量增加,部件2的退化速率也明顯加快,在相同時間內的退化量增加;對部件2進行預防性更換時,部件1的退化速率也會相應變慢,退化仿真軌跡出現轉折。如果不考慮部件間退化的相互影響,即任意部件維修活動不會影響其他部件的退化過程,將會對維修決策優化結果產生影響,在下面的維修費用分析中將詳細說明。

在得到退化過程模型基礎上,為得到部件可靠度函數,將部件退化近似做非減退化過程處理,為驗證該方法的合理性和精確性,在以上參數設置下對部件1和部件2首次到達故障閾值時間進行仿真,分別取樣5 000次,做出首次故障時間的直方圖,得到其近似分布,與解析故障密度函數進行對比,如圖3所示。

圖3　首次故障時間仿真數據和解析分布對比

由結果圖仿真數據和解析分布對比可以看出,所提方法可較為精確地得到部件的故障分布函數,也驗證了當σB較小時,部件非減退化過程假設的合理性。在此基礎上,對系統首次故障時間進行仿真,同樣取樣5 000次,得到系統故障時間分布仿真數據與解析分布對比圖,結果顯示系統可靠度函數式(25)可較為精確地描述系統真實故障分布情況。在用迭代法搜尋模型最優值的過程中,為方便預防性維修策略的實施,預防性維修間隔期的搜索步長取1,即只在整數范圍內求預防性維修間隔期的最優值。當預防性維修間隔期(T1,T)確定時,根據優化模型可得到系統單位時間維修費用Cavg(T1,T,xp)隨預防性更換閾值xp的變化情況,如圖4所示。

圖4　xp與Cavg(T1,T,xp)的關系

圖4中列舉了4種不同預防性維修間隔期下單位時間維修費用隨xp的變化情況,顯示了對于固定的(T1,T),存在最優的預防性維修閾值xp使得單位時間維修費用最小。從整體上看,隨著預防性維修間隔期的減小,xp的最優值是逐漸增大的,這是由于較小的預防性維修間隔期可以更為及時地檢測系統退化量,從而可以在部件達到故障閾值前設置更大的預防性維修閾值以充分利用部件使用壽命。較小的xp值會引起系統過早更換,而較大的xp值則可相對減少過早更換，但同時會增加系統故障風險,這都會導致維修費用增加。同時通過曲線的對比可以說明區分首次預防性維修間隔期與重復預防性維修間隔期的必要性。

當預防性更換閾值xp固定時,系統單位時間維修費用Cavg(T1,T,xp)與預防性維修間隔期(T1,T)的關系如圖5所示(以xp=5為例)。由系統維修費用率變化三維圖和等高線圖可以看出,對于固定的xp,從整體上看,無論是隨著T1或T的增大,Cavg(T1,T,xp)呈現先降低后升高的趨勢,這是預防性維修費用減少和故障更換費用增加共同作用結果。存在最優的(T1,T)使得系統單位時間維修費用最小,例如當xp=5時,最優預防性維修間隔期為T1=16,T=8,此時系統單位時間維修費用為528.27元。

在不考慮部件間相互影響的條件下,兩部件的退化過程不會隨著另一個部件的維修發生變化,所以此時兩部件退化過程可由第一次預防性維修間隔期內的退化模型表示,圖5(b)中的T1=25,T=4點表示在這種情況下對應的優化結果,可以看出如果部件間相互影響,所得到的實際費用和預防性維修間隔期并不是最優值。

圖5　(T1,T)與Cavg(T1,T,xp)的關系(xp=5)

在以上分析基礎上,對不同xp取值下的模型進行優化,得到相應的最優預防性維修間隔期和系統單位時間維修費用值,如圖6所示。

圖6　不同xp下系統最小單位時間維修費用

表2　不同xp對應優化結果

在預防性維修策略實施過程中,由于工作任務或者其他一些因素影響,可能不能嚴格按照最優預防性維修間隔期對部件實施檢測和預防性更換工作,這時就需要分析預防性維修間隔期變化對系統單位時間維修費用產生的影響。

表3　不同(T1,T)下系統單位時間維修費用

5　結　論

本文在非線性Wiener退化過程模型基礎上提出了基于退化速率的退化相關多部件系統退化過程模型和可靠度模型,在此基礎上,建立了系統預防性維修決策優化模型。在模型算例分析部分,通過對不同部件的退化過程仿真說明了部件間的相關關系;通過仿真數據和解析分布的對比,驗證了可靠度模型的精確性;通過優化得到了系統最優預防性維修策略,對比說明了最優預防性維修策略的經濟效益。最后分析了決策變量變化對維修費用的影響,為預防性維修策略的實施提供參考。

參考文獻：

[1] MURTHY D, NGUYEN D. Study of two-component system with failure interaction[J].Naval Research Logistics,1985,32(2): 239-247.

[2] MURTHY D, NGUYEN D. Study of a multi-component system with failure interaction[J]. European Journal of Operational Research, 1985, 21(3): 330-338.

[3] ZEQUEIRA R I, BERENGUER C. On the inspection policy of a two-component parallel system with failure interaction[J]. Reliability Engineering and System Safety, 2005, 88(1): 99-107.

[4] SUN Y, MA L, MATHEW J. An analytical model for interactive failures[J]. Reliability Engineering and System Safety, 2006, 91(5): 495-504.

[5] GOLMAKANI H R, MOAKEDI H. Periodic inspection optimization model for a two-component repairable system with failure interaction[J].Computers and Industrial Engineering,2012,63(3): 540-545.

[6] 周金宇, 謝里陽, 王學敏. 多狀態系統共因失效分析及可靠性模型[J]. 機械工程學報, 2005, 41(6): 66-70.

ZHOU J Y, XIE L Y, WANG X M. Analysis for common cause failure and reliability model in multi-state systems[J]. Journal of Mechanical Engineering, 2005, 41(6): 66-70.

[7] LEVITIN G. Incorporating common-cause failures into non-repairable multistate series-parallel system analysis[J]. IEEE Trans.on Reliability, 2001, 50(4): 380-388.

[8] 李春洋,陳循,易曉山,等.基于馬爾可夫過程的k/n(G)系統共因失效分析[J].系統工程與電子技術,2009,31(11):2793-2796.

LI C Y, CHEN X, YI X S, et al. Analysis of k-out-of-n: G systems subject to common cause failures based on Markov process[J]. Systems Engineering and Electronics, 2009, 31(11): 2793-2796.

[9] HSIEH M H, JENG S L, SHEN P S. Assessing device reliability based on scheduled discrete degradation measurements[J]. Probabilistic Engineering Mechanics, 2009, 24(2): 151-158.

[10] SARI J K, NEWBY M J, BROMBACHER A C. Bivariate constant stress degradation model: LED lighting system reliability estimation with two-stage modeling[J]. Quality and Reliability Engineering International, 2009, 25(8): 1067-1084.

[11] PAN Z, BALAKRISHNAN N. Reliability modeling of degradation of products with multiple performance characteristics based on gamma processes[J]. Reliability Engineering and System Safety, 2011, 96(8): 949-957.

[12] PAN Z, BALAKRISHNAN N, SUN Q, et al. Bivariate degradation analysis of products based on wiener processes and copulas[J]. Journal of Statistical Computation and Simulation,2013,83(7): 1316-1329.

[13] KAISHEV V K, DIMITROVA D S, HAVERMAN S. Modeling the joint distribution of competing risks survival times using copula functions[J]. Insurance, Mathematics and Economics, 2007, 41(3): 339-361.

[14] WANG Y P, PHAM H. Modeling the dependent competing risks with multiple degradation processes and random shock using time-varying copulas[J]. IEEE Trans.on Reliability, 2012, 61(1): 13-22.

[15] 郭馳明.基于獨立增量過程的視情維修優化方法研究[D].長沙: 國防科技大學, 2013.

GUO C M. Condition-based maintenance optimization with independent increments processes[D]. Changsha: National University of Defense Technology, 2013.

[16] RASMEKOMEN N, PARLIKAD A K. Condition based maintenance of multi-component systems with degradation state-rate interactions[J]. Reliability Engineering and System Safety, 2016, 148(1): 1-10.

[17] 劉次華.隨機過程[M].5版.武漢:華中科技大學出版社,2014: 21-22.

LIU C H. Stochastic process[M]. 5th ed. Wuhan: Huazhong University of Science and Technology Press, 2014: 21-22.

[18] SUN Y, MA L, MATHEW J. Failure analysis of engineering systems with preventive maintenance and failure interactions[J]. Computer and Industrial Engineering, 2009, 57(1): 539-549.

[19] 焦寶聰.常微分方程[M].北京:清華大學出版社,2008:177-178.

JIAO B C. Ordinary differential equations[M]. Beijing: Tsinghua University Press, 2008: 177-178.

[20] BAE S J, KUO W, KVAM P H. Degradation models and implied lifetime distributions[J]. Reliability Engineering and System Safety, 2007, 92(5): 601-608.

[21] LIU B, XU Z G, XIE M, et al. A value-based preventive maintenance policy for multi-component system with continuously degrading components[J]. Reliability Engineering and System Safety, 2014, 132(132): 83-89.

[22] WANG X L, BALAKRISHNAN N, GUO B. Residual life estimation based on a generalized wiener degradation process[J].Reliability Engineering and System Safety,2013,124(4):13-23.