二元性能相關(guān)多階段退化系統(tǒng)可靠性建模與分析

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

(延安大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,陜西 延安 716000)

0 引言

隨著現(xiàn)代技術(shù)的不斷進(jìn)步，產(chǎn)品(或設(shè)備)越來越復(fù)雜，而產(chǎn)品的質(zhì)量與可靠性問題日益突顯出來，如何分析產(chǎn)品的可靠性已經(jīng)成為實(shí)際生產(chǎn)過程中必須面對(duì)的一項(xiàng)重要工作[1,2]。傳統(tǒng)的基于失效數(shù)據(jù)的可靠性建模與評(píng)估方法面臨著失效數(shù)據(jù)不足等嚴(yán)重問題[3]?；谛阅芡嘶瘮?shù)據(jù)的可靠性建模與評(píng)估方法為解決這類問題提供了一條有效途徑，已經(jīng)成為可靠性理論與工程領(lǐng)域的研究熱點(diǎn)之一。

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

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

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

1 Copula函數(shù)

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

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

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

(1)

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

(2)

2 二元多階段退化模型

本節(jié)將構(gòu)建具有固定變點(diǎn)和隨機(jī)變點(diǎn)兩種情形下的二元多階段退化模型。

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

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

假設(shè)2每個(gè)性能指標(biāo)在每個(gè)階段的退化過程服從Wiener過程，即對(duì)于k=1,2和t∈[τi-1,τi],i=1,2,…,n, 有

(3)

其中，Xk(t)為性能指標(biāo)k在時(shí)刻t的退化量，Xi,k(t)為性能指標(biāo)k在第i階段時(shí)刻t的退化量，Xi-1,k(τi-1)為性能指標(biāo)k在第i階段的初始退化量，即性能指標(biāo)k在第i-1階段處于變點(diǎn)τi-1時(shí)的退化量，且規(guī)定τ0≡0,X0,k≡0, 即假設(shè)兩個(gè)性能指標(biāo)的初始退化量為0;μi,k和σi,k則分別表示性能指標(biāo)k在第i階段的漂移系數(shù)和擴(kuò)散系數(shù)，{B(t),t≥0}為標(biāo)準(zhǔn)布朗運(yùn)動(dòng)。

假設(shè)3兩個(gè)性能指標(biāo)的退化過程是相關(guān)的，其相關(guān)性用Copula函數(shù)來描述。

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

基于以上假設(shè)，建立多階段二元相關(guān)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],這里所描述的二元相關(guān)退化過程建立在相同增量區(qū)間上，在不同增量區(qū)間[t,t+Δt]上假設(shè)是獨(dú)立的。

注：假設(shè)Copula函數(shù)為Frank Copula函數(shù)，當(dāng)θ→0時(shí)，F(xiàn)i,1(ΔXi,1(t)),Fi,2(ΔXi,1(t))趨向于獨(dú)立[12]，最終模型變?yōu)閮蓚€(gè)性能指標(biāo)獨(dú)立的退化模型，即

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))

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

假設(shè)1系統(tǒng)存在兩個(gè)相關(guān)性能指標(biāo)，每個(gè)性能指標(biāo)的退化過程包含n(n=1,2,…)個(gè)階段。令Γ1,Γ2,…,Γn-1分別表示退化過程的n-1個(gè)隨機(jī)變點(diǎn)的發(fā)生時(shí)刻，令Xi=Γi-Γi-1(i=1,2,…,n)表示兩相鄰變點(diǎn)間的時(shí)間間隔，其中Xi相互獨(dú)立且服從均值為λ的指數(shù)分布。

4 系統(tǒng)可靠度

根據(jù)假設(shè)，系統(tǒng)可靠度函數(shù)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)分別表示兩個(gè)性能指標(biāo)在第i階段時(shí)刻t的可靠度,Fi,1(t)和Fi,2(t)則分別表示兩個(gè)性能指標(biāo)在第i階段首次達(dá)到其失效閾值時(shí)間的分布函數(shù)。

由于系統(tǒng)退化過程包含多個(gè)階段, 下面根據(jù)時(shí)間t與變點(diǎn)的關(guān)系進(jìn)行分類討論, 推導(dǎo)系統(tǒng)的可靠度函數(shù)。

如果0

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

(5)

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

(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)由文獻(xiàn)[26]可知:

(8)

由于Wiener過程的強(qiáng)馬爾科夫性，有

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),且有g(shù)n-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)

在固定變點(diǎn)情形下的模型可靠度計(jì)算公式的基礎(chǔ)上，給出具有隨機(jī)變點(diǎn)的二元多階段隨機(jī)退化系統(tǒng)的可靠度計(jì)算公式。由于系統(tǒng)在不同階段退化速率是不同的，下面對(duì)時(shí)刻進(jìn)行分類討論。

當(dāng)t≤Γ1，時(shí)刻之前未產(chǎn)生變點(diǎn)，此時(shí)系統(tǒng)可靠度為:

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

(12)

當(dāng)t>Γ1，即時(shí)刻t之前至少有一個(gè)變點(diǎn)，如果此時(shí)Γ2≥t≥Γ1，系統(tǒng)可靠度為:

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

(13)

綜上所述，可知在隨機(jī)變點(diǎn)情形下的系統(tǒng)可靠度為:

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

(14)

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

5 數(shù)值模擬

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

(15)

(16)

(17)

進(jìn)一步，設(shè)定退化量初始值x0,退化變點(diǎn)時(shí)刻τ=(τ0,τ1,τ2,…,τn),性能指標(biāo)1在各階段的漂移系數(shù)μ1=(μ1,1,μ2,1,…,μn,1)和擴(kuò)散系數(shù)σ1=(σ1,1,σ2,1,…,σn,1),步長(zhǎng)Δt。代入Xi,k(t)便可得到性能指標(biāo)1的退化量為:

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

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

此時(shí),性能指標(biāo)1的退化增量為:

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

而退化增量的分布函數(shù)表達(dá)式為:

(18)

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

(19)

則性能指標(biāo)2的退化增量表達(dá)式為:

(20)

其中Φ-1(v)、Δt、μi,2和σi,2已知，通過累加退化增量最終可得到性能指標(biāo)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),…)

利用上述方法模擬出系統(tǒng)設(shè)備的退化數(shù)據(jù)后，需要根據(jù)設(shè)備的失效閾值對(duì)設(shè)備進(jìn)行可靠性模擬求值。為此，提出以下可靠度模擬算法：

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

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

步驟3令n=n+1,執(zhí)行步驟1；

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

其中aTN的表示性能1第N次模擬對(duì)應(yīng)時(shí)刻的退化量,bTN表示性能2第N次模擬對(duì)應(yīng)時(shí)刻的退化量；

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

步驟6矩陣X1和X2的相同位置取邏輯且運(yùn)算，將結(jié)果保存在矩陣A對(duì)應(yīng)位置中；

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

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

圖1 二元相關(guān)性能指標(biāo)的退化路徑圖

表1 退化樣本數(shù)據(jù)

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

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

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

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

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

取定退化量初始值為0，產(chǎn)生3階段的2個(gè)隨機(jī)變點(diǎn)，性能指標(biāo)1和2在各階段上的漂移系數(shù)μ1=(4,5,6),μ2=(5,6,7)和擴(kuò)散系數(shù)σ1=(1,1,1),σ2=(1,1,1),步長(zhǎng)Δ=0.01,Copula函數(shù)為Frank Copula函數(shù),參數(shù)θ為10和30,模擬次數(shù)N=10000,閾值D=(10,10)。此時(shí)系統(tǒng)可靠度解析解的計(jì)算涉及無窮積分，采用文獻(xiàn)[29]中的逐項(xiàng)積分定理及相關(guān)算法，系統(tǒng)的可靠度的結(jié)果如圖5所示。

根據(jù)結(jié)果可知:參數(shù)θ值越小，系統(tǒng)的可靠度下降的越早，這是因?yàn)閰?shù)θ決定了相關(guān)程度的度量，當(dāng)θ趨向于0時(shí)，兩個(gè)性能的可靠度愈趨向于獨(dú)立，由文獻(xiàn)[14]可知兩個(gè)退化量如果相互獨(dú)立會(huì)低估系統(tǒng)可靠度。同時(shí)隨著值越大，可靠度曲線下降的越晚，這是因?yàn)樽凕c(diǎn)產(chǎn)生的時(shí)間間隔隨著λ值增加而增加，進(jìn)入退化速率較快的階段較晚，因此失效時(shí)間靠后。

圖3 隨機(jī)變點(diǎn)情況下的可靠度曲線

根據(jù)結(jié)果，為保證系統(tǒng)設(shè)備運(yùn)行可靠性，在設(shè)備性能隨時(shí)間、季度、暖熱等情況發(fā)生變化更替時(shí)要及時(shí)對(duì)系統(tǒng)進(jìn)行維修替換保養(yǎng)，從而避免損失。

6 結(jié)論

本文針對(duì)工程設(shè)備退化過程存在兩個(gè)相關(guān)性能指標(biāo)且表現(xiàn)出多階段特征的問題，利用基于Copula連接函數(shù)的方法，建立了多階段二元Wiener過程退化模型，得到了系統(tǒng)的可靠度函數(shù)，通過蒙特卡洛模擬數(shù)據(jù)驗(yàn)證了其可行性和有效性。

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