趙 騫, 譚 堯, 吳 波, 蔣 平,*, 郭 波

(1. 國防科技大學(xué)系統(tǒng)工程學(xué)院, 湖南 長沙 410005; 2. 國防科技大學(xué)信息通信學(xué)院, 陜西 西安 710106)

0 引 言

冗余是提高系統(tǒng)可靠性的重要技術(shù),在工業(yè)和軍事裝備上被廣泛應(yīng)用[1-3],其可被進一步劃分為熱備、溫備與冷備3種。熱備指的是工作部件與備份部件處于相同的工作狀態(tài),冷備則認為備份部件在貯備期間不工作,失效率為零,而溫備的部件失效率介于熱備和冷備之間。n中取k冷備系統(tǒng)共由n個部件組成,至少需要k個部件正常工作才可以保證系統(tǒng)的正常運行。當(dāng)工作位置部件發(fā)生故障,備份部件依次替換失效的工作部件繼續(xù)工作,直至系統(tǒng)中能夠正常工作的部件總數(shù)小于k時,冷備系統(tǒng)失效。

從裝備質(zhì)量特性與管理的角度,常常會關(guān)心當(dāng)前正常工作的裝備產(chǎn)品還能繼續(xù)運行的時間,這就是產(chǎn)品的剩余壽命。準確的剩余壽命預(yù)測對于工業(yè)系統(tǒng)、裝備系統(tǒng)任務(wù)完成、維修以及更換而言均意義重大，因而對復(fù)雜裝備系統(tǒng)的剩余壽命預(yù)測具有極大的經(jīng)濟價值和戰(zhàn)略意義。然而，因冷備系統(tǒng)結(jié)構(gòu)復(fù)雜,對其剩余壽命估計的研究卻相對較少。

許雙偉等[4]通過分析部件狀態(tài)之間的關(guān)系,并基于最小割集理論提出了一種計算可靠度和平均故障間隔時間的仿真算法。同時,基于不確定性理論,推導(dǎo)了3種不可修冷備系統(tǒng)的可靠度和平均失效時間(mean time to failure, MTTF)的解析表達式[5]。利用極大似然原理、Bayes估計和經(jīng)驗Bayes估計,部件服從Gamma分布的冷備系統(tǒng)失效率、可靠度以及平均壽命的估計方法可參見文獻[6-7]。廣義近似置信限以及E-Bayes方法在文獻[8]中進行了研究。但是以上研究只針對n中取1冷備的特殊情況,并且缺少對剩余壽命估計的討論。更進一步,Eryilmaz[9]通過計算得到了一個冷備部件的表決系統(tǒng)平均剩余壽命,但是當(dāng)部件壽命分布復(fù)雜時,該方法同樣面臨局限性。應(yīng)用Farlie-Gumbel-Morgenstern Copula函數(shù)完成對可靠度和MTTF的估計可參見文獻[10-11]。事實上,這些研究也僅考慮了部件服從指數(shù)分布的情況。

雖然對于冷備系統(tǒng)剩余壽命的研究較少,但是對其可靠度的估計方法的研究相對較為成熟。周仲夏等[12]通過推導(dǎo)得到了n中取1的冷備系統(tǒng)可靠度的閉合表達式,類似研究也可參見文獻[13]。通過運用馬爾可夫過程,2n+1中取n的特殊冷備系統(tǒng)在文獻[14]中進行了研究,類似研究還可參見文獻[15-17]。這些方法的優(yōu)勢是可以獲得閉合表達式,但是以上結(jié)論必須基于部件均服從相同指數(shù)分布這一假設(shè)。

由于計算較為復(fù)雜,數(shù)值計算和Bayes理論等方法得到了進一步應(yīng)用。Levitin等[18]提出了基于通用生成函數(shù)獲取n中取k的冷備系統(tǒng)可靠度函數(shù)的方法。這也為本文進行剩余壽命估計奠定了重要基礎(chǔ)。事實上,在部件壽命已知的情況下,冷備系統(tǒng)可靠度可以通過卷積進行計算,但是往往由于計算復(fù)雜,無法或者難以通過推導(dǎo)得到冷備系統(tǒng)可靠度的解析表達式。Amari等[19]通過梯形法實現(xiàn)了卷積的近似計算。一些學(xué)者嘗試運用Bayes理論進行可靠度估計[20-22]。Jia等[23]基于馬爾可夫鏈-蒙特卡羅方法獲得可靠度Bayes估計。

綜上所述,對于n中取k的冷備系統(tǒng)剩余壽命的估計比較困難,特別是當(dāng)部件服從復(fù)雜分布時,相關(guān)研究更為缺乏?；诖?本文首先研究了剩余壽命估計與可靠度函數(shù)之間的關(guān)系,證明了剩余壽命計算的核心是獲得系統(tǒng)可靠度,然后通過研究給出解析方法、數(shù)值方法以及仿真方法,獲得了冷備系統(tǒng)任意時刻的可靠度,在此基礎(chǔ)上基于組合辛普森公式,提出一種基于任意時刻可靠度近似計算剩余壽命的方法,高效、準確地獲得了冷備系統(tǒng)剩余壽命的點估計與區(qū)間估計。

1 基于可靠度函數(shù)的剩余壽命預(yù)測方法

1.1 剩余壽命的概率密度函數(shù)

假設(shè)T為產(chǎn)品壽命,F(t)、f(t)和R(t) 分別表示產(chǎn)品失效時間的累積分布函數(shù)(cumulative distribution function, CDF)、概率密度函數(shù)以及可靠度函數(shù)。則該產(chǎn)品在工作到當(dāng)前時刻τ時,剩余壽命分布函數(shù)FL(t)[9]為

(1)

式(1)對t求導(dǎo),可得剩余壽命概率密度函數(shù)為

(2)

1.2 可靠度函數(shù)與剩余壽命

τ時刻,剩余壽命點估計Eτ(TL)為

(3)

將式(2)代入(3)可得

其中,

將其代入式(3)可得剩余壽命點估計為

(4)

假設(shè)剩余壽命100(1-α)%置信水平下的雙側(cè)區(qū)間為[μL,μH],則μL與μH需滿足:

(5)

也可記為

(6)

從式(4)和式(6)均不難看出,要得到冷備系統(tǒng)剩余壽命的點估計與區(qū)間估計,核心是獲得該系統(tǒng)可靠度。本文將在第2節(jié),針對冷備系統(tǒng)組成部件壽命分布的不同情況,依次介紹解析方法、數(shù)值方法以及仿真方法3種可靠度計算方法。

2 可靠度的計算

2.1 解析方法計算可靠度

對于部件服從相同指數(shù)分布的情況,根據(jù)已有文獻研究,可以通過解析推導(dǎo)得到冷備系統(tǒng)可靠度的閉合表達式。

例如n中取1的冷備系統(tǒng),假設(shè)n個部件的壽命為X1,X2,…,Xn,相互獨立且均服從失效率為λ0的指數(shù)分布,則冷備系統(tǒng)的可靠度可以表示為

R(t)=1-P(X1+X2+…+Xn≤t)

由拉普拉斯變換及逆變換[12]可得

(7)

更一般地,對于n中取k冷備系統(tǒng),部件仍服從失效率為λ0的指數(shù)分布,可靠度[24]為

(8)

文獻[25]推導(dǎo)了該種情況下τ時刻冷備系統(tǒng)剩余壽命的閉合表達式為

(9)

(10)

式(9)與式(10)的推導(dǎo)結(jié)果將在算例中進行運用,作為該種情況下的理論值與本文第2.2節(jié)所提方法進行對比。

2.2 數(shù)值方法計算可靠度

由于指數(shù)分布相對簡單,可以推導(dǎo)得到對應(yīng)的冷備系統(tǒng)剩余壽命閉合表達式,對于更為復(fù)雜的情況,數(shù)值方法提供了新的解決思路。設(shè)k個工作部件中第i個(1≤i≤k)工作部件和對應(yīng)該位置的備份部件壽命分布CDF分別為Fi(t)和Hi(t),令Gi, j(t)表示第i個部件所在的位置至少有j個部件已經(jīng)失效。不難得到

Gi, j(t)=Gi, j-1(t)*Hi(t)

(11)

式中：Gi,1(t)=Fi(t)。并且，

(12)

關(guān)于式(11)中的卷積,可以將積分區(qū)間劃分為m個子區(qū)間,采用梯形法則進行數(shù)值計算,計算精度隨著m的增大而提高[19]。令Pi, j(t)表示t時刻在i個部件所在位置上正好發(fā)生j個失效的概率,因而

Pi, j(t)=Gi, j(t)-Gi, j+1(t)

(13)

且Pi,0(t)=1-Gi,1(t)。

令Pj(t)表示t時刻系統(tǒng)中共有j個部件失效,容易看出Pj(t)可由所有位置的Pi, j(t)計算得到。則當(dāng)系統(tǒng)失效個數(shù)小于n-k+1時,系統(tǒng)正常工作,所以冷備系統(tǒng)的可靠度為

(14)

值得一提的是,上述數(shù)值計算方法要求同一位置上工作的部件需要服從相同分布,因而也面臨一定局限性。

2.3 仿真方法計算可靠度

對于更為復(fù)雜的情況,可采用仿真抽樣從而得到任意時刻的可靠度,抽樣算法如算法1所示。

算法1 冷備系統(tǒng)可靠度估計抽樣方法已知:給定系統(tǒng)n個部件壽命分布概率密度函數(shù)fi(t)(1≤i≤n) 和抽樣次數(shù)S,任意工作時刻τ,令count=0。步驟1 由fi(t)(1≤i≤n)分別抽樣得到部件壽命樣本L1,L2,…,Ln,令工作位置的累計壽命CLi=Li(1≤i≤k),p=1;步驟2 尋找到所有CLi (1≤i≤k)中的最小值,在該CLi上增加Lk+p,并令p=p+1;步驟3 重復(fù)步驟2,(n-k)次后,可以得到最終k個CLi值,則系統(tǒng)壽命Ls=miniCLi(1≤i≤k);步驟4 If Ls≥τ, then count+1;步驟5 重復(fù)步驟1～步驟4 S次,τ時刻的可靠度近似值為R(τ)=countS。

需要指出的是,仿真方法適用于分布更為復(fù)雜的情況,在實際應(yīng)用中,僅需要將算法進行對應(yīng)調(diào)整即可,本文在算例中也將進行詳細分析。

3 剩余壽命的近似計算

3.1 組合辛普森公式

第2節(jié)討論了不同情況下獲得冷備系統(tǒng)任意時刻可靠度的方法。由式(4)可知,要對分子中的積分進行計算，相對復(fù)雜,很多情況下難以推導(dǎo)得到閉合表達式,即使獲得了冷備系統(tǒng)可靠度的解析表達,由于其形式復(fù)雜,積分計算工作量也很大,因此本節(jié)基于組合辛普森公式引入近似計算方法。近似計算方法簡單高效,適用范圍廣,表現(xiàn)出了更大的優(yōu)勢。利用可靠度函數(shù)與剩余壽命的關(guān)系,在通過解析、數(shù)值、仿真方法獲得任意時刻系統(tǒng)可靠度的基礎(chǔ)上,可以通過求解得到剩余壽命的點估計與區(qū)間估計。

對于式(4),本文采用組合辛普森公式進行計算。

如果f(x)∈C4[a,b],則存在ξ∈[a,b],使得[a,b]上w個子區(qū)間的組合辛普森公式和截斷誤差項[26-27]滿足:

(15)

式中：h=(b-a)/w;f(4)(ξ)為f(x)在ξ的四階導(dǎo)數(shù)。

(16)

式中：M為足夠大的正數(shù);w為劃分的區(qū)間個數(shù),且h=(M-τ)/ω。顯然,M的取值對結(jié)果會產(chǎn)生影響,相關(guān)分析見第4.1.1節(jié)。

3.2 誤差分析

上述方法計算剩余壽命的誤差主要存在于采用梯形法則近似計算卷積產(chǎn)生的誤差,采用辛普森公式求解積分產(chǎn)生的截斷誤差以及將無限區(qū)間的積分轉(zhuǎn)化為有限區(qū)間積分所產(chǎn)生的截斷誤差。一般地,區(qū)間劃分個數(shù)m,ω越大,以及所選有限區(qū)間[τ,M]中的區(qū)間上限M越大,計算精度越高,誤差也越小。剩余壽命方差推導(dǎo)如下:

(17)

其中,

(18)

所以,隨機變量TL的方差為

(19)

4 算例分析

4.1 部件服從指數(shù)分布的冷備系統(tǒng)

以5中取2的冷備系統(tǒng)為研究對象,假設(shè)部件完全相同并且相互獨立,均服從失效率λ0=0.01h-1的指數(shù)分布。

4.1.1 區(qū)間的確定

由式(10)和式(19),可得壽命分布的期望為E(TL)=200 h,方差為var(TL)=1×104h2。表1展示了在積分區(qū)間上限M的不同設(shè)定下,壽命的點估計計算結(jié)果。

表1 不同積分區(qū)間下的相對誤差

由表1可知,一般情況下,取值區(qū)間上限M越大,計算所得到的壽命以及剩余壽命點估計值越準確。分析表1可知,積分區(qū)間上限為E(TL)+4σ時,誤差已經(jīng)明顯降到1%以下,因而可確定剩余壽命積分區(qū)間范圍為[τ,E(TL)+4σ]時,能基本滿足誤差要求,當(dāng)然在計算時間允許范圍內(nèi),M的取值越大,估計結(jié)果越精確。

4.1.2m的確定

在確定積分區(qū)間上限M的基礎(chǔ)上,需要研究使用梯形法則時,區(qū)間劃分個數(shù)m的確定方法和原則。為了更清晰地進行比較,首先利用式(1)和式(8)獲得當(dāng)前工作時刻為τ=20 h時剩余壽命的理論CDF,然后通過本文提出的數(shù)值方法,分別畫出m分別取10,100,500時,剩余壽命的CDF圖像，如圖1所示。

圖1 不同m設(shè)定下的CDF比較Fig.1 Comparison of CDFs with different m

需要說明的是,為了便于對比,圖1(a)和圖1(b)均插入了局部放大圖,且放大比例依次增大。通過分析圖1(a)可以看到,在3種m設(shè)定下,誤差均比較小。但在繼續(xù)局部放大得到圖1(b)后可以發(fā)現(xiàn),m=10所對應(yīng)的CDF圖像與理論CDF圖像存在較明顯的偏差,但是m=100和m=500所得到的兩條曲線仍然沒有明顯區(qū)別。因而一般情況下,m值越大，結(jié)果越準確,考慮計算成本,令m=100更好。

4.1.3w的確定

表2 不同w下的相對誤差表

4.1.4 剩余壽命驗證

第4.1.1至第4.1.3節(jié)先后研究了影響近似計算方法精度的因素,為了進一步驗證在部件服從相同指數(shù)分布情況下,n中取k的冷備系統(tǒng)剩余壽命估計方法的準確性,改變參數(shù)值,利用式(7)求得冷備系統(tǒng)可靠度函數(shù),進而利用式(16)計算出剩余壽命近似值,并將該近似計算結(jié)果與由文獻[25]推導(dǎo)的剩余壽命點估計閉合表達式(式(9))進行比較，比較結(jié)果如表3所示。

表3 不同參數(shù)下的剩余壽命估計結(jié)果

由表3可以看出,在部件壽命服從相同指數(shù)分布的情況下,本文提出的數(shù)值方法求解冷備系統(tǒng)剩余壽命的方法計算結(jié)果準確,誤差均在可接受范圍內(nèi)。

4.2 部件服從相同威布爾分布的冷備系統(tǒng)

為進一步驗證本文所提方法,假設(shè)部件壽命服從相同的威布爾分布,威布爾分布CDF如下:

(20)

式中:η,β為威布爾分布參數(shù)[28-30]。

針對部件服從威布爾分布的情況,運用第2.2節(jié)數(shù)值方法獲得任意時刻該冷備系統(tǒng)的可靠度,然后運用辛普森公式得到剩余壽命的近似估計。為了進行結(jié)果驗證,與仿真方法進行對比,設(shè)定抽樣次數(shù)S=10 000，驗證結(jié)果如表4所示。

表4 部件服從威布爾分布時不同參數(shù)下的剩余壽命估計結(jié)果

由于CDF圖像能夠反應(yīng)更詳細的比較信息,因而根據(jù)表4第1組實驗參數(shù)設(shè)置(n=5,k=2,η=100,β=2,τ=20),利用仿真方法和數(shù)值方法繪制的剩余壽命CDF圖像如圖2所示。

圖2 部件服從威布爾分布情況下仿真方法和數(shù)值方法 得到的CDF比較Fig.2 Comparison of CDFs between simulation method and numerical method with components following Weibull distribution

仍以表4第1組實驗為例,在得到任意時刻可靠度的基礎(chǔ)上,由式(9)可以分別得到兩種方法的剩余壽命區(qū)間估計,假設(shè)置信水平為90%,比較結(jié)果如表5所示。

表5 實驗1剩余壽命點估計與區(qū)間估計

通過表4、表5和圖2可以發(fā)現(xiàn),本文提出的數(shù)值方法無論從剩余壽命點估計、區(qū)間估計或是CDF圖像考量,結(jié)果都比較準確,誤差在可接受范圍內(nèi)。

4.3 更為復(fù)雜情況下的冷備系統(tǒng)

對于更為復(fù)雜的情況,第2.1節(jié)與第2.2節(jié)提出的解析方法與數(shù)值方法均無法獲得冷備系統(tǒng)可靠度,所以采用第2.3節(jié)所提出的仿真抽樣方法。以8中取4的冷備系統(tǒng)為例,假設(shè)4個位置部件均服從威布爾分布,分布參數(shù)(η,β)分別設(shè)定為(110,2.25), (120,2.50), (130,2.75), (140,3.00)。其工作到τ=20 h 時系統(tǒng)未失效,通過仿真抽樣,近似計算得到其剩余壽命的估計值為150.95 h,區(qū)間估計為(90.80,218.11)h。其剩余的CDF圖像如圖3所示。

圖3 不同位置部件服從不同威布爾分布時冷備系統(tǒng) 剩余壽命CDF圖像Fig.3 CDF image of residual life of cold standby system for components at different positions following different Weibull distributions

仍以上述情況為例,假設(shè)初始時,工作部件與冷備部件分別服從不同威布爾分布,分布參數(shù)(η,β)分別為 (100,2), (140,3), 其工作到τ=20 h 時系統(tǒng)未失效,則通過計算得到剩余壽命估計值為140.55 h,區(qū)間估計為 (71.11,219.08)h。剩余的CDF圖像如圖4所示。

圖4 部件服從不同威布爾分布時冷備系統(tǒng)剩余壽命CDF圖像Fig.4 CDF image of residual life of cold standby system for components following different Weibull distributions

通過上述實驗可以看出,對于更為復(fù)雜的情況,解析方法、數(shù)值方法面臨局限性,但是仿真方法在犧牲了可接受的時間成本的基礎(chǔ)上,彌補了前者的不足。

5 實例計算

冷備作為重要冗余技術(shù),在航空航天產(chǎn)品中得到了非常廣泛的應(yīng)用。例如，文獻[31]研究了一種工程中存在的典型冷備系統(tǒng),即指數(shù)-威布爾混合型冷備系統(tǒng)。該系統(tǒng)工作部件為電子設(shè)備,服從指數(shù)分布,而冷貯備部件為機電設(shè)備,服從威布爾分布。由于推導(dǎo)復(fù)雜,文獻[31]也只是研究了一備一系統(tǒng),即2中取1的冷備系統(tǒng)。本文針對文獻[31]所提到的開關(guān)完全可靠和開關(guān)壽命為指數(shù)分布兩種不同情況進行了對比驗證。

5.1 開關(guān)完全可靠

當(dāng)轉(zhuǎn)換開關(guān)完全可靠時,在模擬分析中,文獻[31]假設(shè)工作部件服從失效率為λ=0.15的指數(shù)分布,冷貯備部件服從η=2,β=1.5的威布爾分布。

運用本文第2.2節(jié)所提出的基于梯形法則的計算方法,可以快速得到任務(wù)時刻τ=2時,該冷備系統(tǒng)的可靠度為0.916 7。這與文獻[31]中的計算結(jié)果一致,說明本文所提方法準確合理。

文獻[32]基于Fiducial方法給出了冷備系統(tǒng)的置信下限,該方法同時也可以得到冷備系統(tǒng)可靠度的點估計。運用該方法計算得到上述系統(tǒng)可靠度點估計為0.915 1,與本文方法及文獻[31]中方法的計算結(jié)果一致,再一次證明了本文方法的準確性。同時,利用文獻[32]中方法計算得到可靠度95%置信區(qū)間的下限為0.904 7。多次利用本文第2.3節(jié)中的仿真方法得到該系統(tǒng)的多個可靠度點估計樣本,同樣也可以計算得到可靠度置信下限,為0.903 0,與文獻[32]的計算結(jié)果一致。說明本文方法也可以很好地解決置信下限點估計問題。

進一步,在獲得任意時刻該冷備系統(tǒng)可靠度后,利用辛普森公式,計算得到該系統(tǒng)在工作時刻τ=2時，剩余壽命為7.105 5,剩余壽命90%置信區(qū)間估計為(0.518 5,20.825 7)。

根據(jù)剩余壽命計算結(jié)果,可以合理安排該冷備系統(tǒng)的更換、維修等策略,與文獻[31]僅僅給出冷備系統(tǒng)可靠度相比,文獻[32]給出的可靠度置信下限更具意義。

5.2 開關(guān)壽命為指數(shù)分布

文獻[31]還研究了轉(zhuǎn)換開關(guān)不可靠的情形,并且假設(shè)轉(zhuǎn)換開關(guān)失效時間服從指數(shù)分布,失效率為λb=0.15。同時，工作部件服從失效率λ為 0.15的指數(shù)分布,冷貯備部件服從η=3,β=1.5的威布爾分布。此時,需要采用第2.3節(jié)的仿真方法獲取冷備系統(tǒng)可靠度。

根據(jù)仿真抽樣樣本,該冷備系統(tǒng)壽命樣本及其擬合概率密度曲線如圖5所示。

圖5 開關(guān)壽命為指數(shù)分布時冷備系統(tǒng)壽命概率密度曲線Fig.5 Probability density function curve of cold standby system life with exponentially distributed switch life

通過計算可以得到,任務(wù)時刻τ=2時的可靠度為0.918 7,這與文獻[31]中的計算結(jié)果(0.918 8)基本一致。同時,根據(jù)本文方法可以得到任意時刻可靠度,所以運用辛普森公式便可以計算得到開關(guān)壽命為指數(shù)分布時該冷備系統(tǒng)的剩余壽命為6.711 7,90%剩余壽命置信區(qū)間為(0.573 9,18.960 2)。

由于結(jié)構(gòu)復(fù)雜,文獻[31]的研究具有如下局限性:① 僅考慮了2中取1冷備系統(tǒng);② 工作部件和備份部件分布形式固定,分別服從指數(shù)分布與威布爾分布;③ 僅研究了可靠度估計,并未給出剩余壽命估計。本文所提方法解決了上述難題,并且通過對比實驗再次說明,本文所提出的通過仿真獲取任意時刻可靠度、進而估計剩余壽命的方法能夠解決部件分布復(fù)雜、開關(guān)不可靠等更加多樣、復(fù)雜的實際情形,具有很好的適用性。

6 結(jié) 論

冷備系統(tǒng)日益受到關(guān)注,但是因其結(jié)構(gòu)復(fù)雜,針對復(fù)雜冷備系統(tǒng)剩余壽命的研究較少?；诖?本文提出了n中取k冷備系統(tǒng)剩余壽命點估計與區(qū)間估計的近似計算方法。通過推導(dǎo)發(fā)現(xiàn),對于冷備系統(tǒng)剩余壽命估計的核心是獲得其可靠度函數(shù)。本文提供了不同情況下,獲得冷備系統(tǒng)可靠度的解析方法、數(shù)值方法以及仿真方法。在此基礎(chǔ)上,基于組合辛普森公式,利用上述方法得到任意時刻可靠度,便可完成冷備系統(tǒng)剩余壽命的點估計和區(qū)間估計。

算例分析部分討論了與估計誤差密切相關(guān)的積分上限M、卷積區(qū)間劃分個數(shù)m、積分區(qū)間劃分個數(shù)w的一般確定原則,并且通過部件服從相同指數(shù)分布、威布爾分布兩種情況對本文提出的近似計算方法進行了檢驗。當(dāng)部件服從更為復(fù)雜的分布時,本文所提方法仍可以解決剩余壽命點估計和區(qū)間估計問題,適用范圍廣。

同時,本文在最后的實例計算部分,采用已有論文數(shù)據(jù)進行計算,并與原論文計算結(jié)果進行了對比,說明了本文計算方法準確可靠,同時也解決了原論文研究中存在的局限,具有較大的工程應(yīng)用價值。