董慶來(lái), 王偉偉, 司書賓

(1.延安大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院, 陜西 延安 716000; 2.西北工業(yè)大學(xué) 機(jī)電學(xué)院, 陜西 西安 710072)

關(guān)于隨機(jī)退化系統(tǒng)的研究可以追溯到20世紀(jì)60年代，但是直到20世紀(jì)90年代才引起國(guó)內(nèi)外學(xué)者的重視，近年來(lái)隨著傳感器技術(shù)和計(jì)算技術(shù)的發(fā)展，隨機(jī)退化系統(tǒng)的可靠性建模與分析已經(jīng)發(fā)展為可靠性理論與工程領(lǐng)域研究的熱點(diǎn)之一[1-4]。目前，在多元隨機(jī)退化系統(tǒng)的研究中，通常假設(shè)系統(tǒng)的退化程度由每個(gè)性能指標(biāo)的退化量直接進(jìn)行描述[5];除此之外，在實(shí)踐中，還經(jīng)常將性能指標(biāo)融合成一個(gè)變量間接地描述系統(tǒng)的健康狀態(tài)[6]。例如，船舶艉軸承是船舶動(dòng)力裝置的主要組成部分，艉軸承間隙是描述船舶艉軸承健康狀態(tài)的重要指標(biāo)，而艉軸承間隙可通過(guò)艉軸直徑磨損量和艉軸承襯磨損量來(lái)確定，當(dāng)艉軸直徑磨損量和艉軸承襯磨損量之和(即將二者融合為一個(gè)新的變量)達(dá)到失效閾值(即允許的間隙最大值)時(shí)，船舶艉軸承失效[7]。楊圓鑒[7]采用Wiener過(guò)程描述艉軸直徑磨損量和艉軸承襯磨損量的退化過(guò)程，對(duì)船舶艉軸承的可靠性進(jìn)行了評(píng)估，但是為了計(jì)算方便，假設(shè)二者是相互獨(dú)立的。實(shí)際上，由于受共同的工作環(huán)境及應(yīng)力影響，艉軸直徑磨損量和艉軸承襯磨損量的退化過(guò)程是相關(guān)的，并且系統(tǒng)性能指標(biāo)退化過(guò)程的相關(guān)性在許多系統(tǒng)中廣泛存在，因此在多元隨機(jī)退化系統(tǒng)的可靠性建模與分析中考慮性能指標(biāo)的相關(guān)性更加符合實(shí)際，但是多元性能指標(biāo)之間的相關(guān)性卻給退化建模帶來(lái)了巨大的挑戰(zhàn)[4,8]，其主要原因在于多元隨機(jī)過(guò)程(例如多元Wiener過(guò)程)的首達(dá)時(shí)分布是目前學(xué)術(shù)界無(wú)法解決的難題[9]。作為多元隨機(jī)退化系統(tǒng)的特殊情形，二元隨機(jī)退化系統(tǒng)是最基本、最重要的情形，是研究多元隨機(jī)退化系統(tǒng)的基礎(chǔ)，不僅能夠反映多元情形下的許多本質(zhì)問(wèn)題，而且相比多元情形在數(shù)學(xué)上更易于處理[10]，因此本文以二元隨機(jī)退化系統(tǒng)為主要研究對(duì)象。

在隨機(jī)退化系統(tǒng)的可靠性建模中，除了建立描述系統(tǒng)退化過(guò)程的退化模型，還要確定系統(tǒng)失效的準(zhǔn)則，兩者缺一不可。學(xué)者們通常將系統(tǒng)的壽命定義為退化過(guò)程關(guān)于失效閾值的首達(dá)時(shí)間，一般稱這種失效為退化失效[11-13]。事實(shí)上，許多系統(tǒng)的退化過(guò)程存在警戒閾值[14]，在退化量達(dá)到警戒閾值前，系統(tǒng)處于正常狀態(tài)，但是超過(guò)警戒閾值后，系統(tǒng)將處于缺陷狀態(tài)，例如警戒水位、腐蝕裕量等[15-16]，系統(tǒng)在2種狀態(tài)下的退化速率通常是不同的。此時(shí)，如果采用相同參數(shù)的隨機(jī)過(guò)程對(duì)系統(tǒng)兩階段退化過(guò)程進(jìn)行描述，可能會(huì)導(dǎo)致可靠性指標(biāo)計(jì)算的誤差[17]。此外，上述系統(tǒng)除了發(fā)生因退化導(dǎo)致的退化失效外，系統(tǒng)還會(huì)發(fā)生突發(fā)失效[18-19]，例如系統(tǒng)在缺陷狀態(tài)的持續(xù)時(shí)間過(guò)長(zhǎng)，系統(tǒng)也會(huì)失效。Dong和Cui[20]考慮退化失效和持續(xù)時(shí)間失效，構(gòu)建了基于非平穩(wěn)隨機(jī)過(guò)程的兩階段退化模型，給出了系統(tǒng)可靠度等可靠性指標(biāo)。Qiu和Cui[21]考慮不同失效準(zhǔn)則，給出了最優(yōu)維修策略。上述文獻(xiàn)均以一元隨機(jī)退化系統(tǒng)為研究對(duì)象，而在二元或多元隨機(jī)退化系統(tǒng)的可靠性建模與分析中均只考慮退化失效[22-26]，因此探討二元或多元隨機(jī)退化系統(tǒng)的競(jìng)爭(zhēng)失效分析問(wèn)題具有重要的現(xiàn)實(shí)意義。

基于上述分析，本文的主要研究?jī)?nèi)容包括：①針對(duì)二元隨機(jī)退化系統(tǒng)，考慮性能指標(biāo)的相關(guān)性、性能指標(biāo)融合以及退化失效和持續(xù)時(shí)間失效2類競(jìng)爭(zhēng)失效形式，構(gòu)建基于二元Wiener過(guò)程的可靠性模型。②在性能指標(biāo)不進(jìn)行融合以及性能指標(biāo)按照退化量的線性組合進(jìn)行融合情形下，推導(dǎo)持續(xù)時(shí)間失效閾值為常數(shù)和隨機(jī)變量時(shí)的系統(tǒng)可靠度解析表達(dá)式。③提出基于蒙特卡洛模擬的數(shù)值模擬算法,模擬性能指標(biāo)按照任意形式進(jìn)行融合情形下二元以及多元隨機(jī)退化系統(tǒng)的可靠度，并用于驗(yàn)證解析結(jié)果的正確性。

1 問(wèn)題描述與可靠性建模

考慮二元隨機(jī)退化系統(tǒng)，根據(jù)性能指標(biāo)是否進(jìn)行融合，構(gòu)建2個(gè)可靠性模型。性能指標(biāo)融合是指利用新變量描述性能指標(biāo)之間的關(guān)系，即系統(tǒng)的退化程度由性能指標(biāo)退化量的函數(shù)進(jìn)行描述。假設(shè)X(t)=(X1(t),X2(t))′,Xi(t)為系統(tǒng)的第i個(gè)性能指標(biāo)在時(shí)刻t的退化量(i=1,2),則在性能指標(biāo)融合情形下,系統(tǒng)在時(shí)刻t的退化程度由X(t)的函數(shù)H(x1,x2)進(jìn)行描述,其中x1和x2為X1(t)和X2(t)的樣本實(shí)現(xiàn)值。

模型1性能指標(biāo)不進(jìn)行融合情形下的可靠性模型。

1) 兩階段退化過(guò)程

假設(shè)系統(tǒng)的退化過(guò)程包括正常階段(或第一階段,即從新系統(tǒng)投入運(yùn)行到缺陷發(fā)生時(shí)刻)和缺陷階段(或第二階段,即從缺陷發(fā)生時(shí)刻到失效發(fā)生時(shí)刻),從而系統(tǒng)的狀態(tài)包括正常、缺陷和失效。令

(1)

式中:Xik(t)是系統(tǒng)的第i個(gè)性能指標(biāo)在第k個(gè)階段時(shí)刻t的退化量,i=1,2,k=1,2;Γ1是變點(diǎn),定義為2個(gè)性能指標(biāo)首達(dá)警戒閾值的時(shí)刻,若定義隨機(jī)變量γi=inf{t>0:Xi1(t)≥Di1},i=1,2,則Γ1=min{γ1,γ2}。

2) 二元Wiener過(guò)程退化模型

由于受共同的環(huán)境應(yīng)力,2個(gè)性能指標(biāo)的退化過(guò)程是相關(guān)的,在每個(gè)階段采用二元相關(guān)Wiener過(guò)程描述系統(tǒng)的退化程度,即

X(t)=

(2)

式中:x0=(x110,x210)′是系統(tǒng)的初始退化量,并假設(shè)x0=(0,0)′,即系統(tǒng)在時(shí)刻t=0是全新的;X(Γ1)=(X11(Γ1),X21(Γ1))′為系統(tǒng)在變點(diǎn)處的退化量向量。對(duì)于i=1,2,k=1,2,μk=(μ1k,μ2k)′為漂移系數(shù)向量;μik為系統(tǒng)第i個(gè)性能指標(biāo)在第k個(gè)階段的漂移系數(shù);Σk是第k個(gè)階段的擴(kuò)散矩陣,使得

為協(xié)方差矩陣,ρk是相關(guān)系數(shù),描述2個(gè)性能指標(biāo)在第k個(gè)階段的相關(guān)性,σ1k和σ2k為標(biāo)準(zhǔn)差;W(t)=(W1(t),W2(t))′,{W1(t),t≥0}和{W2(t),t≥0}是相互獨(dú)立的標(biāo)準(zhǔn)Wiener過(guò)程。具體的,如果記

3) 競(jìng)爭(zhēng)失效

當(dāng)一個(gè)性能指標(biāo)的退化量超過(guò)其退化失效閾值時(shí),或系統(tǒng)在缺陷階段的持續(xù)時(shí)間超過(guò)持續(xù)時(shí)間失效閾值τ時(shí),系統(tǒng)失效。系統(tǒng)的退化路徑如圖1所示,其中γ和y分別為Γ1和Y的實(shí)現(xiàn)值。

圖1 模型1下系統(tǒng)的退化路徑

因此,模型1下系統(tǒng)的壽命為

Y1=min{Λ1,τ}+Γ1

(3)

模型2性能指標(biāo)融合情形下的可靠性模型。

本文考慮H(x1,x2)的3種表達(dá)形式:

①H(x1,x2)=k1x1(t)+k2x2(t),k1,k2∈R;②H(x1,x2)=|x1(t)-x2(t)|;③一般形式。

與模型1的不同之處在于變點(diǎn)以及系統(tǒng)壽命的表達(dá)形式。在模型2下,變點(diǎn)定義為2個(gè)性能指標(biāo)退化量的函數(shù)關(guān)于警戒閾值L1的首達(dá)時(shí)刻,即

Γ2=inf{t≥0:H(X11(t),X21(t))≥L1}

(4)

則系統(tǒng)的壽命為

Y2=min{Λ2,τ}+Γ2,

(5)

2 系統(tǒng)可靠度

本節(jié)將給出系統(tǒng)可靠度的解析表達(dá)式,包括模型1和模型2下H(x1,x2)=k時(shí)2種情形,以及模型2下H(x1,x2)為一般函數(shù)時(shí)可靠度的模擬算法。

1) 模型1下系統(tǒng)可靠度

根據(jù)失效準(zhǔn)則,系統(tǒng)的可靠度為

R(t)=P{Y1>t}=P{min{Λ1,τ}+Γ1>t}

(6)

首先,假設(shè)持續(xù)時(shí)間失效閾值為常數(shù)。

注意到Γ1=min{γ1,γ2}以及γi=inf{t≥0:Xi(t)≥Di1},i=1,2,通過(guò)對(duì)一個(gè)性能指標(biāo)達(dá)到失效閾值的時(shí)刻及另一個(gè)未達(dá)到失效閾值的性能指標(biāo)在該時(shí)刻的數(shù)值取條件,得到

(7)

當(dāng)其中的一個(gè)性能指標(biāo)達(dá)到失效閾值時(shí),系統(tǒng)將進(jìn)入缺陷階段。根據(jù)Xi2(t)(i=1,2)的定義,當(dāng)

(8)

式中,f(λ1,λ2)(u1,u2|x0i)的表達(dá)式由文獻(xiàn)[27]的公式(2)給出。將(8)式代入(7)式,得到

(9)

在工程實(shí)踐中,持續(xù)時(shí)間失效閾值通常是隨機(jī)變量。下面假設(shè)持續(xù)時(shí)間失效閾值是具有概率密度函數(shù)fτ(t)和分布函數(shù)Fτ(t)的隨機(jī)變量,并考慮以下2種情形：一種情形是分布在[0,+∞)上的隨機(jī)變量;一種情形是分布在[τ1,τ2]上的隨機(jī)變量,其中τ1和τ2是正常數(shù)。

如果持續(xù)時(shí)間失效閾值τ是分布在[0,+∞)上的隨機(jī)變量,例如τ是服從指數(shù)分布或韋布爾分布的隨機(jī)變量,則(9)式變?yōu)?/p>

類似的,可以得到持續(xù)時(shí)間失效閾值τ是分布在[τ1,τ2]上的隨機(jī)變量時(shí)的可靠度。

2) 模型2下H(x1,x2)=k時(shí)系統(tǒng)可靠度

根據(jù)Γ2的定義,如果τ是常數(shù),系統(tǒng)可靠度為

(10)

式中:f(X1,Γ2)(x1,u)是X1(Γ2)和Γ2的聯(lián)合概率密度函數(shù),滿足

f(X1,Γ2)(x1,u)=fX11(u)(x1)fΓ2(u)

(11)

(12)

從而,要得到(11)式的解析表達(dá)式,只需得到變點(diǎn)Γ2的概率密度函數(shù)fΓ2(u)即可。然而,除了H(x1,x2)的一些特殊形式外,fΓ2(u)的解析表達(dá)式是很難獲得的。

下面,給出H(x1,x2)=k1x1+k2x2情形下的系統(tǒng)可靠度。由(2)式知

H(X11(t),X21(t))=(k1μ11+k2μ21)t+δ1B1(t)

(13)

式中，{B1(t)=[(k1a111+k2a211)W1(t)+(k1a121+k2a221)W2(t)]/δ1,t≥0}是標(biāo)準(zhǔn)Wiener過(guò)程，

由Γ2的定義,即Γ2=inf{t>0:H(X11(t),X21(t))≥L1},fΓ2(u)為{Z1(t),t≥0}關(guān)于L1的首達(dá)時(shí)的概率密度函數(shù),其中Z1(t)(k1μ11+k2μ21)t+δ1B1(t)。因此,由(13)式可以得到

(14)

將(12)式和(14)式代入(11)式,得到

(15)

接下來(lái),給出FΛ2(t-u|x1,x2)的解析表達(dá)式。注意到k1x1+k2x2=L1,與(13)式類似

L1+(k1μ12+k2μ22)t+δ2B2(t),t≥0

(16)

將(15)式和(16)式代入(10)式,得到H(x1,x2)=k1x1+k2x2情形下的系統(tǒng)可靠度

(17)

注：①如果令τ=0,(17)式將變?yōu)閱坞A段一元隨機(jī)退化系統(tǒng)的可靠度[11];②如果令τ→+∞,(17)式將變?yōu)閮呻A段一元隨機(jī)退化系統(tǒng)的可靠度[11];③如果τ是隨機(jī)變量,通過(guò)積分可以得到相應(yīng)的系統(tǒng)可靠度的解析表達(dá)式。

3) 模型2下H(x1,x2)為一般函數(shù)時(shí)系統(tǒng)可靠度的數(shù)值模擬

由(11)式可以發(fā)現(xiàn),性能指標(biāo)按照任意形式進(jìn)行融合情形下的系統(tǒng)可靠度表達(dá)式在形式上是簡(jiǎn)單的,但是對(duì)一般H(x1,x2)而言,Λ2和Γ2的分布函數(shù)的表達(dá)式并不容易求解。基于蒙特卡洛模擬的算法來(lái)模擬系統(tǒng)可靠度,具體的模擬算法流程如算法1所示。

算法1模型2下H(x1,x2)為一般函數(shù)時(shí)系統(tǒng)可靠度的模擬算法

步驟1設(shè)定參數(shù)值,包括初始時(shí)刻t=0,初始退化量x0,漂移參數(shù)向量μ1和μ2,擴(kuò)散矩陣Σ1和Σ2,警戒閾值L1,退化失效閾值L2,持續(xù)時(shí)間失效閾值τ和其他相關(guān)參數(shù)。令模擬次數(shù)和步長(zhǎng)分別為N和Δt。

步驟2對(duì)于i=1,2,…,N,執(zhí)行步驟3至4。

步驟5完成模擬,根據(jù)失效數(shù)據(jù)集A計(jì)算系統(tǒng)可靠度。

注：1) 如果持續(xù)時(shí)間失效閾值是隨機(jī)變量,只需在步驟4中增加在指定的分布中生成獨(dú)立同分布的隨機(jī)數(shù)即可。

2) 算法1還適用于基于性能指標(biāo)融合的多元隨機(jī)退化系統(tǒng)的可靠度計(jì)算,即隨機(jī)退化系統(tǒng)具有多個(gè)性能指標(biāo)(例如n個(gè)性能指標(biāo))時(shí),只需做如下修改:①在步驟3和步驟4中生成多元正態(tài)分布隨機(jī)數(shù),即擴(kuò)散矩陣Σ1和Σ2分別為n×n階矩陣;②采用H(x1,x2,…,xn)描述性能指標(biāo)融合,即將各性能指標(biāo)退化量的數(shù)值代入H(x1,x2,…,xn)描述系統(tǒng)的退化程度。

3 數(shù)值算例

本節(jié)給出一些數(shù)值算例來(lái)驗(yàn)證所構(gòu)建模型的有效性,并對(duì)模型中的參數(shù)進(jìn)行靈敏度分析,探討模型參數(shù)對(duì)隨機(jī)退化系統(tǒng)可靠性指標(biāo)的影響。

假設(shè)2個(gè)性能指標(biāo)退化量的警戒閾值和退化失效閾值分別為0.2和0.4,持續(xù)時(shí)間失效閾值τ=4。在第一階段,μ11=μ21=0.016 7,σ11=σ21=0.045 0,ρ1=-0.786 2;在第二階段,μ12=μ22=0.024 0,σ12=σ22=0.060 4,ρ2=-0.694 3,從而兩階段的協(xié)方差矩陣分別為

1) 模型1下系統(tǒng)可靠度

圖2 模型1下基于模擬解和解析解的系統(tǒng)度曲線

圖3 持續(xù)時(shí)間失效閾值服從β分布情形下系統(tǒng)可靠度曲線

2) 模型2下H(x1,x2)=k1x1+k2x2時(shí)系統(tǒng)可靠度

假設(shè)系統(tǒng)的警戒閾值、退化失效閾值和持續(xù)時(shí)間失效閾值分別為0.3，0.5和4.0,其余參數(shù)取值與圖3相同。為了描述2個(gè)性能指標(biāo)退化量的線性組合對(duì)系統(tǒng)可靠度的影響,考慮k1和k2取不同參數(shù)值情形下的系統(tǒng)可靠度,即:k1=1.5且k2=1.2、k1=1且k2=1、k1=0.5且k2=0.5、k1=0.4且k2=0.3、k1=1且k2=0。將可靠度曲線畫在同一個(gè)圖中,如圖4所示,其中k1=1且k2=0時(shí)的系統(tǒng)可靠度曲線可以看作一元隨機(jī)退化系統(tǒng)的可靠度曲線。

圖4 H(x1,x2)=k1x1+k2x2情形下系統(tǒng)可靠度曲線

從圖4可以發(fā)現(xiàn),k1和k2的值越大,系統(tǒng)的可靠度越低。這是因?yàn)閗1和k2反映了2個(gè)性能指標(biāo)退化量對(duì)系統(tǒng)整體退化量的貢獻(xiàn),它們的值越大,說(shuō)明系統(tǒng)的退化速度越快,從而系統(tǒng)的可靠度越低,這與我們的直覺是一致的。另一個(gè)現(xiàn)象是,當(dāng)k1和k2的取值較小時(shí),其可靠度在退化過(guò)程的早期階段高于單變量Wiener過(guò)程退化模型的結(jié)果,而在退化過(guò)程的后期結(jié)果卻相反,這可以用2個(gè)性能指標(biāo)的相關(guān)性(2個(gè)性能指標(biāo)是負(fù)相關(guān)的)來(lái)解釋。

3) 模型2下H(x1,x2)=|x1-x2|時(shí)系統(tǒng)可靠度分析

不等式|x1-x2|>L1曾被Cui等[8]用作平衡系統(tǒng)的失效準(zhǔn)則,即當(dāng)2個(gè)性能指標(biāo)退化量的距離大于失效閾值L1時(shí),系統(tǒng)失效;采用Markov過(guò)程描述2個(gè)性能指標(biāo)的退化過(guò)程,并假設(shè)它們是相互獨(dú)立的。下面,利用算法1給出2個(gè)性能指標(biāo)相關(guān)時(shí)采用函數(shù)H(x1,x2)=|x1-x2|描述系統(tǒng)退化程度情形下的系統(tǒng)可靠度曲線。假設(shè)τ=0,L1=0.1,即只考慮單階段退化過(guò)程。為了得到2個(gè)性能指標(biāo)之間的相關(guān)性對(duì)系統(tǒng)可靠度的影響,分別令ρ1=-0.8，-0.5，0，0.5和0.8,將可靠度曲線畫在同一個(gè)圖中,如圖5所示。其余參數(shù)取值與圖2相同。從圖5可以發(fā)現(xiàn),在其余參數(shù)值相同的情形下,相關(guān)系數(shù)ρ1的值越大,系統(tǒng)的可靠度越高。因此,2個(gè)性能指標(biāo)的退化過(guò)程是正相關(guān)(或負(fù)相關(guān))時(shí),如果假設(shè)它們彼此獨(dú)立,則會(huì)低估(或高估)系統(tǒng)可靠度。

圖5 模型2下H(x1,x2)=|x1-x2|時(shí)單階段隨機(jī)退化系統(tǒng)可靠度曲線

圖6 模型2下H(x1,x2)=|x1-x2|時(shí)兩階段隨機(jī)退化系統(tǒng)可靠度曲線

進(jìn)一步,為了探究持續(xù)時(shí)間失效閾值對(duì)系統(tǒng)可靠度的影響,假設(shè)L1=0.1,L2=0.3,其余參數(shù)取值與圖3相同,考慮τ取不同參數(shù)值情形下的系統(tǒng)可靠度,包括τ=0,4,7,10和50,并將可靠度曲線畫在同一個(gè)圖中,如圖6所示。從圖6可以發(fā)現(xiàn),持續(xù)時(shí)間失效閾值越小,系統(tǒng)的可靠度越低。這意味著系統(tǒng)處于缺陷階段的可承受時(shí)間越短,系統(tǒng)發(fā)生失效的可能性越大,這與文獻(xiàn)[20]的結(jié)論是一致的。

4 結(jié) 論

本文針對(duì)二元隨機(jī)退化系統(tǒng)，考慮性能指標(biāo)融合以及競(jìng)爭(zhēng)失效，提出了基于二元相關(guān)Wiener過(guò)程的退化模型，給出了性能指標(biāo)不進(jìn)行融合以及性能指標(biāo)按照退化量的線性組合進(jìn)行融合情形下，持續(xù)時(shí)間失效閾值為常數(shù)和隨機(jī)變量時(shí)的系統(tǒng)可靠度的解析表達(dá)式。當(dāng)性能指標(biāo)按照任意形式進(jìn)行融合時(shí)，提出了一個(gè)算法程序來(lái)模擬系統(tǒng)的可靠度。研究結(jié)果表明持續(xù)時(shí)間失效閾值越大，系統(tǒng)的可靠性越高。此外，當(dāng)性能指標(biāo)按照2個(gè)性能指標(biāo)退化量的距離進(jìn)行融合時(shí)，若2個(gè)性能指標(biāo)正相關(guān)(或負(fù)相關(guān))而假設(shè)它們彼此獨(dú)立，則會(huì)低估(或高估)系統(tǒng)可靠度。因此，在工程實(shí)踐中，如果存在2個(gè)(或多個(gè))性能指標(biāo)，應(yīng)充分考慮它們的相關(guān)性，從而提高可靠性指標(biāo)計(jì)算的精度。此外，本文給出的算法盡管能夠用于性能指標(biāo)按照任意形式進(jìn)行融合情形下二元以及多元隨機(jī)退化系統(tǒng)可靠度的模擬，但是由于該算法是基于蒙特卡洛模擬的，因此算法的精度強(qiáng)烈依賴于模擬次數(shù)，因此在今后研究中尋求系統(tǒng)可靠度的解析表達(dá)式以及其他算法將更加有意義。