二元性能相關多階段退化系統可靠性建模與分析

秦紅志, 孫明娟, 王偉偉, 董慶來

(延安大學 數學與計算機科學學院,陜西 延安 716000)

0 引言

隨著現代技術的不斷進步，產品(或設備)越來越復雜，而產品的質量與可靠性問題日益突顯出來，如何分析產品的可靠性已經成為實際生產過程中必須面對的一項重要工作[1,2]。傳統的基于失效數據的可靠性建模與評估方法面臨著失效數據不足等嚴重問題[3]。基于性能退化數據的可靠性建模與評估方法為解決這類問題提供了一條有效途徑，已經成為可靠性理論與工程領域的研究熱點之一。

基于性能退化的可靠性建模與評估的前提是假設系統存在一個性能指標能充分反映系統的健康狀態，但是許多系統通常具有多個性能指標，即系統的健康狀態需由多個相關性能指標共同反映[4]。例如，可用容量和可用能量是表征鋰離子電池退化的兩個重要指標[5]。因此，如何準確描述具有多個相關性能指標系統的退化過程并計算相關可靠性指標，給學者們帶來了巨大的挑戰[6,7]。近年來，學者們將Copula函數引入到退化建模中，描述多個性能指標退化過程的相關性[8,9]。Li等[10]基于Copula函數建立了二元退化模型，分析了產品的可靠性并給出了剩余壽命預測的方法。周義蛟[11]和劉小平[12]研究了基于Wiener退化過程和Copula函數的二元相關性能可靠性模型。張志鵬[13]在其碩士論文中給出了系統多元相關退化過程的建模和可靠性評估方法，詳細介紹歸納了基于Copula函數的建模方法。賈旭杰[14]介紹了隨機Copula模型及其模型參數估計方法，基于該模型提出了串并聯動態相依系統的可靠度計算方法。鮑兆偉[15]提出了確定邊緣退化量分布函數和Copula函數的方法，進一步提高了可靠性模型的精確性。

上述文獻大多假設系統遵循單階段隨機退化過程。然而在實踐中，系統的退化過程會受系統內部和環境等因素的綜合影響，系統的退化過程將呈現出兩階段甚至多階段的特性[16～18]。例如，鋰離子電池在開始經歷一個平穩退化期，隨著充放電次數的增加，由于內部損耗以及環境溫度變化等，導致鋰電池容量在第二階段迅速退化[19]。對于此類存在變點的多階段隨機退化系統的研究，已取得一定進展。張鵬[20]研究了考慮隨機效應下的兩階段退化系統剩余壽命預測。黃金波[3]研究了多階段可校正系統退化模型，并給出了系統可靠性的評估方法。Gao等[21]考慮到在整個生命周期中，系統可能會經歷各種環境，提出了一種基于Wiener過程的多階段退化模型來描述失效閾值和漂移系數的突變現象。高洪達[22]在其博士論文中詳細闡述了具有變動閾值、固定變點、隨機變點及競爭失效的多階段一維隨機退化系統的可靠性建模與計算問題。

不難發現，現有成果以單階段多元退化系統和多階段一元退化系統為研究對象，而具有多個相關性能指標且呈現多階段退化特性的復雜系統的可靠性建模與評估研究卻鮮有報道。在工程實際中，鋰離子電池的退化過程會呈現緩慢退化、加速退化、瀕臨失效退化的多階段退化過程[23]，鋰離子電池容量隨著工作時間增加而不斷減少，同時鋰離子電池阻抗也將隨著鋰離子電池容量的減少而發生變化[13]。本文將針對這一問題，提出一類多階段二元隨機退化模型，采用Copula函數描述性能指標間的相關關系，并推導系統的可靠度函數。

1 Copula函數

Copula函數由Sklar[24]在1959年提出，它將一個n(n≥2)維隨機變量的聯合分布函數分解成n個邊緣分布函數和一個Copula函數，該函數描述了變量間的相關關系。由于Copula函數將聯合分布函數與它們各自的邊緣分布函數連接在一起，因此Copula函數也被稱為連接函數[25]。

本文使用的是二元Copula函數。令c(u,v)表示一個二元Copula函數，則相關隨機變量X和Y的聯合分布函數F(X,Y)可用一個二元Copula函數來描述：

F(X,Y)=C(u;vθ)

(1)

式中，u=Fs(X),v=Fy(F)分別為變量X和Y的邊緣分布函數，且u,v均服從均勻分布，θ為Copula函數的參數。Copula密度函數為:

(2)

2 二元多階段退化模型

本節將構建具有固定變點和隨機變點兩種情形下的二元多階段退化模型。

為建立固定變點下的二元多階段退化模型，下面首先給出幾點假設。

假設1系統存在兩個相關性能指標，每個性能指標的退化過程包含n=(n=1,2,…)個階段。令τ1,τ2,…,τn-1分別表示退化過程的n-1個固定的變點，且滿足0<τ1<τ2<…<τn。在變點τi(i=1,2,…,n)之前，系統在第i-1階段運行，如果在第i-1階段系統未失效，則系統進入第i階段。

假設2每個性能指標在每個階段的退化過程服從Wiener過程，即對于k=1,2和t∈[τi-1,τi],i=1,2,…,n, 有

(3)

其中，Xk(t)為性能指標k在時刻t的退化量，Xi,k(t)為性能指標k在第i階段時刻t的退化量，Xi-1,k(τi-1)為性能指標k在第i階段的初始退化量，即性能指標k在第i-1階段處于變點τi-1時的退化量，且規定τ0≡0,X0,k≡0, 即假設兩個性能指標的初始退化量為0;μi,k和σi,k則分別表示性能指標k在第i階段的漂移系數和擴散系數，{B(t),t≥0}為標準布朗運動。

假設3兩個性能指標的退化過程是相關的，其相關性用Copula函數來描述。

假設4當系統的任意一個性能指標的累積退化量超過其失效閾值時，系統失效，記兩個性能指標的失效閾值分別為D1,D2。因此，系統的壽命為Ti=min(Ti,1,Ti,2)。其中Ti,k=inf{t≥0|Xi,k(t)≥Dk}表示性能指標k首次達到其失效閾值的時間。

基于以上假設，建立多階段二元相關Wiener過程退化模型:

ΔXi,k(t)=Xi,k(t+Δt)-Xi,k(t)

Fi(ΔXi,1(t),ΔXi,2(t))

=C(Fi,1(ΔXi,1(t)),Fi,2(ΔXi,2(t));θ)

其中i=1,2,…,n;k=1,2;t∈[τi-1,τi],這里所描述的二元相關退化過程建立在相同增量區間上，在不同增量區間[t,t+Δt]上假設是獨立的。

注：假設Copula函數為Frank Copula函數，當θ→0時，Fi,1(ΔXi,1(t)),Fi,2(ΔXi,1(t))趨向于獨立[12]，最終模型變為兩個性能指標獨立的退化模型，即

F(ΔXi,1(t),(ΔXi,2(t))

=C(Fi,1(ΔXi,1(t)),Fi,2(ΔXi,1(t));θ)

=Fi,1(ΔXi,1(t))·Fi,2(ΔXi,2(t))

在具有固定變點的二元多階段退化模型的基礎上，提出了隨機變點二元多階段退化可靠性模型，該模型進一步貼近了工程實際中的應用情況。保持其他假設條件不變的情形下，將具有固定變點的二元多階段退化模型的假設1修改為:

假設1系統存在兩個相關性能指標，每個性能指標的退化過程包含n(n=1,2,…)個階段。令Γ1,Γ2,…,Γn-1分別表示退化過程的n-1個隨機變點的發生時刻，令Xi=Γi-Γi-1(i=1,2,…,n)表示兩相鄰變點間的時間間隔，其中Xi相互獨立且服從均值為λ的指數分布。

4 系統可靠度

根據假設，系統可靠度函數R(t)為:

Ri(t)=P(Ti>t)=p(min(Ti,1,Ti,2)>t)

=P(Ti,1>t,Ti,2>t)

=1-P(Ti,1≤t)-P(Ti,2≤t)+P(Ti,1≤t,Ti,2≤t)

=Ri,1(t)+Ri,2(t)-1+C(Fi,1(t),Fi,2(t);θ)

(4)

其中Ri,1(t)和Ri,2(t)分別表示兩個性能指標在第i階段時刻t的可靠度,Fi,1(t)和Fi,2(t)則分別表示兩個性能指標在第i階段首次達到其失效閾值時間的分布函數。

由于系統退化過程包含多個階段, 下面根據時間t與變點的關系進行分類討論, 推導系統的可靠度函數。

如果0

R1(t|x0,μ1,σ1,D)

(5)

其中μ1=(μ1,1,μ1,2),σ1=(σ1,1,σ1,2),D=(D1,D2),則該系統壽命的概率密度函數為

(6)

如果τ1

R2(t|X(τ1),μ2,σ2,D)

Dk,X1,1(τ1)=x1,1X1,2(τ1)=x1,2}×

(7)

其中X(τ1)=(X1,1(τ1),X1,2(τ1))=(x1,1,x1,2),μ2=(μ2,1,μ2,2),σ2=(σ2,1,σ2,2)由文獻[26]可知:

(8)

由于Wiener過程的強馬爾科夫性，有

g1,k(x1,k|x0,μ1,k,σ1,k,τ1,Dk)dx1,k,dx2,k

(9)

如果τn-1

Rn(t|X(τn-1),μn,σn,D)

(10)

其中式(18)中x(τn-1)=(xn-1,1,xn-1,2),μn=(μn,1,μn,2),σn=(σn,1,σn,2),且有gn-1,k(xn-1,k)=gn-1,k(xn-1,k|xn-2,k,μn-1,k,σn-1,k,τn-1,Dk),而

gn-1,k(xn-1,k|x0,μn-1,k,σn-1,k,τn-1,Dk)dxn-1,k

(11)

在固定變點情形下的模型可靠度計算公式的基礎上，給出具有隨機變點的二元多階段隨機退化系統的可靠度計算公式。由于系統在不同階段退化速率是不同的，下面對時刻進行分類討論。

當t≤Γ1，時刻之前未產生變點，此時系統可靠度為:

R(t)=R1(t|x0,μ1,σ1,D)×P{N(t)=0}

(12)

當t>Γ1，即時刻t之前至少有一個變點，如果此時Γ2≥t≥Γ1，系統可靠度為:

P{N(t)=2}×f(t1,t2|N(t)=2)dt1dt2

(13)

綜上所述，可知在隨機變點情形下的系統可靠度為:

f(t1,t2,…,tn|N(t)=m)dt1dt2…dtm

(14)

其中Rj(t|·)(j=1,2,…,n)表示系統在第j階段的系統可靠度，其解析表達式為(10)。P{N(t)=m}表示系統在t時刻之前有m個變點的概率，{N(t),t≥0}為強度是λ的Poisson過程，而f(·|N(t)=m)表示t時刻前有m個變點的條件下，隨機變點分布的聯合條件概率密度函數。

5 數值模擬

通過數值模擬驗證本文提出的二元相關多階段退化模型的有效性。使用蒙特卡洛方法(MC)產生相關的仿真樣本數據，將仿真模擬結果與解析結果進行比對，最終分析判斷模型的準確性。由模型假設，首先要仿真得到相關的退化增量數據，然后對增量數據依次累加，最終得到相關的退化量數據。定義Copula函數的偏導數Cu(v)為:

(15)

(16)

(17)

進一步，設定退化量初始值x0,退化變點時刻τ=(τ0,τ1,τ2,…,τn),性能指標1在各階段的漂移系數μ1=(μ1,1,μ2,1,…,μn,1)和擴散系數σ1=(σ1,1,σ2,1,…,σn,1),步長Δt。代入Xi,k(t)便可得到性能指標1的退化量為:

X1=(X1,1(0),X1,1(Δt),…,X2,1(τ1),

X2,1(τ1+Δt),…,Xn,1(τn-1),Xn,1(τn-1+Δt),…)

此時,性能指標1的退化增量為:

ΔX1=(ΔX1,1(τ0),ΔX1,1(τ0+Δt),…)

而退化增量的分布函數表達式為:

(18)

將性能指標1的退化增量分布函數u與服從(0,1)上均勻分布的s代入上式，便可得到性能指標2的退化增量分布函數v。同理可知性能指標2的分布函數為:

(19)

則性能指標2的退化增量表達式為:

(20)

其中Φ-1(v)、Δt、μi,2和σi,2已知，通過累加退化增量最終可得到性能指標2的退化量為:

X2=(X1,2(0),X1,2(Δt),X1,2(2Δt),…

X2,2(τ1),X2,2(τ1+Δt),X2,2(τ1+2Δt),…

Xn,2(τn-1),Xn,2(τn-1+Δt),Xn,2(τn-1+2Δt),…)

利用上述方法模擬出系統設備的退化數據后，需要根據設備的失效閾值對設備進行可靠性模擬求值。為此，提出以下可靠度模擬算法：

步驟1設定模擬次數N,構造兩個空矩陣X1和X2。獲取第n次模擬的兩組不同性能的退化數據，分別記錄在矩陣X1和X2中的第n列中；

步驟2判斷N與n的大小，若n

步驟3令n=n+1,執行步驟1；

步驟4得到矩陣X1和X2分別為:

其中aTN的表示性能1第N次模擬對應時刻的退化量,bTN表示性能2第N次模擬對應時刻的退化量；

步驟5遍歷矩陣元素，判斷是否超出閾值D,若超出則記錄對應位置元素為1, 否則記錄對應位置元素為0；

步驟6矩陣X1和X2的相同位置取邏輯且運算，將結果保存在矩陣A對應位置中；

步驟7對矩陣A每一行求和得到失效總次數并與總模擬次數做比值，該比值即為系統的模擬失效概率。

取定退化量初始值為0, 退化變點時刻τ=(1,2,3),性能指標1和2在各階段上的漂移系數μ1=(4,5,6),μ2=(5,6,7)和擴散系數σ1=(1,1,1),σ2=(1,1,1),步長Δt=0.1,Copula函數為Frank Copula函數,同時假設各階段參數θ均為10, 模擬次數N=10000,閾值D=(10,10)。表1中所示為部分退化模擬數據。圖1為其對應的退化路徑。

圖1 二元相關性能指標的退化路徑圖

表1 退化樣本數據

由于式(15)、式(18)和式(19)涉及多元函數的多重積分，且被積分的函數表達式結果難以用顯示表達式表達，只能采用數值方法進行計算，下面將采用蒙特卡洛模擬方法給出近似解[28]。

最終得到可靠度解析解和模擬解的曲線如圖2所示。

圖2 模擬解與解析解的可靠度曲線

利用解析表達式和模擬算法可以得到可靠度模擬解和解析解的曲線如圖4所示。從圖4可以發現基于解析解和模擬解的系統可靠度曲線基本一致。此外，還可以看出在運行初始階段，系統可靠度保持在非常高的狀態，隨著系統的運行，系統可靠度逐漸下降。

為得到隨機變點情形下系統的可靠度曲線，首先需要生成一系列隨機變點Γi，然后代入到固定變點情形下的算例當中。設置初始時刻Γ0=0，生成服從指數分布的X1，得到一階段變點Γ1=Γ0+X1,隨后生成服從指數分布的X2，得到二階段變點Γ2=Γ1+X2。以此類推，最終得到一系列服從指數分布的隨機變點Γj(j=1,2,…,n-1)。

取定退化量初始值為0，產生3階段的2個隨機變點，性能指標1和2在各階段上的漂移系數μ1=(4,5,6),μ2=(5,6,7)和擴散系數σ1=(1,1,1),σ2=(1,1,1),步長Δ=0.01,Copula函數為Frank Copula函數,參數θ為10和30,模擬次數N=10000,閾值D=(10,10)。此時系統可靠度解析解的計算涉及無窮積分，采用文獻[29]中的逐項積分定理及相關算法，系統的可靠度的結果如圖5所示。

根據結果可知:參數θ值越小，系統的可靠度下降的越早，這是因為參數θ決定了相關程度的度量，當θ趨向于0時，兩個性能的可靠度愈趨向于獨立，由文獻[14]可知兩個退化量如果相互獨立會低估系統可靠度。同時隨著值越大，可靠度曲線下降的越晚，這是因為變點產生的時間間隔隨著λ值增加而增加，進入退化速率較快的階段較晚，因此失效時間靠后。

圖3 隨機變點情況下的可靠度曲線

根據結果，為保證系統設備運行可靠性，在設備性能隨時間、季度、暖熱等情況發生變化更替時要及時對系統進行維修替換保養，從而避免損失。

6 結論

本文針對工程設備退化過程存在兩個相關性能指標且表現出多階段特征的問題，利用基于Copula連接函數的方法，建立了多階段二元Wiener過程退化模型，得到了系統的可靠度函數，通過蒙特卡洛模擬數據驗證了其可行性和有效性。

本文研究的多階段二元性能相關退化模型是多階段多元相關的基礎，基于多元相關多階段模型的退化建模與可靠性分析問題可在本文研究結論的基礎上拓展推廣得到。同時，在可靠度解析解的計算過程中，涉及大量繁瑣的復雜多重積分計算，算法時間復雜度的降低是進一步的研究工作核心。但相較于可靠度模擬計算方法來說，本文給出的可靠度計算方法精度較高，運算次數較少。針對模型的參數估計問題和應用性評價是可靠性重要的研究環節，未來將以此模型為基礎，進行退化失效數據的分析，利用對應的統計方法對模型進行參數估計。