屏蔽數據并聯系統的可靠性分析

師小琳,師義民,段俊杰

(1.西安郵電大學 電子工程學院, 陜西 西安　710121；2.西北工業大學 應用數學系, 陜西 西安　710072;3.西北大學 研究生院,陜西 西安　710069)

對于多個部件組成的系統，為了解系統的性能和可靠性指標,通常要對系統進行壽命試驗。在理想狀態下,試驗不僅可以觀測到系統的失效時間,還可以確定引起其失效的原因(即系統的失效是由哪一個部件失效而引起的)。然而對某些系統,由于其結構的復雜性以及試驗經費、設備、技術等各種條件的限制,人們往往無法獲知引起系統失效的確切部件,而只能將其失效原因歸結為一些部件組成的集合,從而系統的失效原因就被屏蔽了。這種試驗數據稱為屏蔽數據,相應的系統稱為屏蔽數據系統。近年來,關于屏蔽數據系統的可靠性分析已有較多研究成果。文獻[1]首次在屏蔽數據下給出觀測樣本的似然函數,并利用極大似然估計法對指數部件的可靠性進行分析;Sarhan[2-3]分別討論了串聯系統中部件服從Weibull分布及幾何分布時,部件參數和可靠性指標的極大似然估計和Bayes估計。Kuo等[4]基于屏蔽數據研究了威布爾部件串聯系統中未知參數的極大似然估計(MLE)和Bayes估計(BE)。相關研究可參見文獻[5-9]。

上述研究主要針對屏蔽數據串聯系統。關于屏蔽數據并聯系統的研究相對較少。文獻[10-11]分別研究了廣義指數分布、廣義瑞利分布并聯系統的可靠性分析。文獻[12]研究了兩部件并聯系統中補充指數型分布參數的MLE與BE。相關的研究可見文獻[13-14]。但是，在已有的文獻中,關于屏蔽數據并聯系統的研究是在完全樣本或Ⅰ型、Ⅱ型截尾試驗下完成的。與上述采樣方式不同,逐步II型截尾試驗允許在試驗的不同階段從試驗中移走一些未失效的樣品,以供研究產品的某些性能(例如退化性能)。這種截尾方式可用于縮短試驗時間和保證一定量失效數據二者間的權衡。本文在逐步II型截尾試驗下,研究補充指數分布三部件屏蔽數據并聯系統的可靠性分析問題。討論部件參數及可靠度函數的MLE和BE，并推導未知參數的最大后驗密度置信區間。最后給出數值例子,對估計結果進行驗證比較。

1　試驗模型與似然函數

(T1,S1),(T2,S2),…,(Tm,Sm),m

對于上述試驗,我們假設在同一時刻有且僅有一個系統失效,且屏蔽的發生與失效原因無關,即

P(Si=si|Ti=ti,Ki=j)=

P(Si=si|Ti=ti,Ki=j′),

?j,j′∈si,j≠j′

若以隨機變量Tij表示系統i中部件j的壽命,則Ti=max(Ti1,Ti2,…,TiJ),i=1,2,…,n,j=1,…,J。

根據上述假設,觀測樣本的似然函數可由下述定理給出。

定理1設并聯系統中第j(j=1,…,J)個部件的密度函數和分布函數分別為fj(t),Fj(t),則在逐步Ⅱ型截尾試驗下,觀測樣本的似然函數為

(1)

證明由文獻[12]知,系統i的壽命Ti的密度函數為

在m個受試系統全部失效且試驗中途無移走的情形下,樣本的似然函數為

若有m個并聯系統參與逐步II型截尾壽命試驗,在第i個失效時刻Ti,從未失效系統中任意移走Ri個,則由文獻[15]可得如下形式的似然函數

(2)

本文假設并聯系統中各部件相互獨立,第j(j=1,2,…,J)個部件的壽命服從參數為λj的補充指數分布,其概率密度函數與分布函數分別為:

fj(t)=λjt-2exp{-λj/t},t>0,λj>0，

(3)

Fj(t)=exp{-λj/t},t>0,λj>0。

(4)

2　參數及可靠度函數的估計

下面考慮3個部件組成的串聯系統(J=3),我們以mj表示由于第j個部件失效而引發的系統失效個數, 即失效原因為部件集合Sj={j},j=1,2,3時系統的失效數;以m12,m13,m23，m123分別表示當失效原因為部件集合{1,2},{1,3},{2,3},{1,2,3}時系統的失效數。顯然有

2.1　極大似然估計

下面考慮并聯系統中補充指數部件參數及可靠度函數的MLE。令λ=(λ1,λ2,λ3),將式(3),(4)代入式(1)中,即可得似然函數:

(λ1+λ2)m12·(λ1+λ3)m13·

(λ2+λ3)m23·(λ1+λ2+λ3)m123

(5)

m12ln(λ1+λ2)+m13ln(λ1+λ3)+

由此可得下列方程組:

(6)

2.2　貝葉斯估計

當系統的失效原因被完全屏蔽時無法得到極大似然估計。這時考慮部件參數及可靠度函數的貝葉斯估計。記|Ε|為集合E的元素個數。對于觀測數據

D={(T1,S1),(T2,S2),…,(Tm,Sm)},記Aj={i||si|=1,Ki=j},j=1,2,3

以Sl分別表示屏蔽集{1,2}{1,3}{2,3}{1,2,3},記Ml={i||si|>1,si=Sl},l=1,2,3,4。此時似然函數可表示為

P(Ki=j|Ti=ti,Si=si)=

(7)

I(0,∞)(λj),j=1,2,3

I(0,∞)(λj)}

(8)

利用貝葉斯公式及式(7) 及(8),可以得到λ的聯合后驗密度函數為

(9)

由于λ的聯合后驗密度函數計算復雜。故可應用Metropolis-Hastings(M-H)抽樣算法[16]模擬參數(λ1,λ2,λ3)的后驗樣本。M-H算法計算步驟如下:

5)q=q+1;

并聯系統部件j可靠度函數在時刻t的貝葉斯估計為

(j=1,2,3)。

2.3　區間估計

對于給定的置信度1-α,從后驗分布獲得的置信區間不止一個,常用的方法是把α平分,用α/2和1-α/2分位數來獲得λ1,λ2,λ3的置信區間,此即為實際中常用的等尾置信區間。但最理想的置信區間應是區間長度最短的,需要把具有最大后驗密度的點都包含在區間內,而在區間外的點上的后驗密度函數值不超過區間內的后驗密度函數值。此時的區間稱為最大后驗密度(Highest posterior density,即HPD)置信區間[17]。

構造參數λ1,λ2,λ3的HPD置信區間的具體步驟如下:

λj(1)≤λj(2)≤…≤λj(L),j=1,2,3;

3) 計算λ1,λ2,λ3的100(1-α)%置信區間,即[λj(q),λj(q+(1-α)*L)],j=1,2,3,;q=1,2,…,L-(1-α)L；

4) 在[λj(q),λj(q+(1-α)*L)]，q=1,2,…,L-(1-α)L中尋找區間長度最短的區間,將其作為參數λj，的HPD置信區間(j=1,2,3)。

3　數值模擬

假定有60個完全相同的系統參加逐步II截尾壽命試驗 (n=60), 每個系統由3個相互獨立的部件并聯而成。部件壽命分別服從參數為λ1=1,λ2=1.2,λ3=0.7的補充指數分布,并設定在第40個失效發生時停止試驗,即m=40。試驗中采用下列3種不同的逐步移走方案。

方案1:R1=R2=…=Rm-3=0,Rm-2=Rm-1=1,Rm=n-m-2;

方案2:R1=R2=…=0,Rm-3=n-m-2,Rm-2=0,Rm-1=Rm=1;

方案3:R1=R2=…=RK+1=1,RK+2=Rm-1=0,Rm=n-m-1-K。

其中K=[m/4]+1,[m/4]表示不超過m/4的最大整數。

重復上述模擬步驟2 000次,計算出部件參數及可靠度函數MLE的平均值(仍記為MLE),BE的平均值(仍記為BE)及相應的均方誤差(MSE),并將上述計算結果列于表1中。均方誤差的計算公式為

按照上一節給出的方法,在不同移走方案(CS)及屏蔽水平p下計算參數置信度為95%的HPD置信區間(HCI)、區間長度(CIL)及覆蓋率(CP),并將計算結果列于表2中。

表1不同CS及屏蔽水平p下參數和可靠度函數的估計及均方誤差(括號中為MSE, “---”表示無法計算)

Tab.1Estimation and MSE of parameters and reliability functions for different CS and masking levelsp(The number in brackets is MSE,"---" means that it cannot be calculated)

參數CS　　　　　　　p=02　　　　　　　　　　　　　p=06　　　　　　　　　　　　p=1　　　　　　MLE(MSE)BE(MSE)MLE(MSE)BE(MSE)MLE(MSE)BE(MSE)λ1[1]10380(00298)08214(00179)11034(00328)10143(00528)---11027(00763)[2]12408(00276)08957(00248)08815(00284)09061(00262)---10879(00535)[3]08973(00252)08373(00239)09241(00295)09426(00256)---09513(00425)λ2[1]14283(00701)10646(00642)15482(00836)09124(00746)---13568(00827)[2]13418(00739)13092(00685)09122(00898)12645(00860)---13466(00946)[3]12771(00631)14105(00502)13885(00768)12642(00802)---12348(00925)λ3[1]08204(00798)04858(00699)05293(00889)05727(00773)---07312(00879)[2]07244(00711)07906(00694)09967(00856)08504(00768)---07312(00895)[3]09107(00706)08136(00659)09374(00890)07548(00743)---07257(00859)R1[1]03914(00082)04342(00074)03140(00096)02780(00087)---02758(00098)[2]03605(00082)02721(00071)02546(00094)02727(00078)---02940(00089)[3]02350(00075)02565(00067)03262(00083)02871(00074)---02860(00089)

續表1

參數CS　　　　　　　p=02　　　　　　　　　　　　　p=06　　　　　　　　　　　　p=1　　　　　　MLE(MSE)BE(MSE)MLE(MSE)BE(MSE)MLE(MSE)BE(MSE)R2[1]04111(00090)03380(00075)04034(00102)02516(00087)---03315(00095)[2]03505(00087)03364(00066)02642(00109)03249(00078)---03478(00085)[3]03547(00084)03631(00064)03705(00097)03173(00083)---03209(00097)R3[1]02491(00051)03126(00045)02619(00063)02378(00053)---02989(00076)[2]02345(00053)02207(00042)02827(00075)02179(00051)---02195(00075)[3]02718(00051)02408(00035)02684(00063)02161(00049)---02082(00068)

表2　不同CS及屏蔽水平p下參數95% 的HCI、區間長度(CIL)及覆蓋率(CP)Tab.2　95% HCI,CIL and CP for different CS and masking levels p

由表1及表2中的數據可以看出:

1) 在相同屏蔽水平和移走方案下,Bayes估計的均方誤差(MSE)較MLE的MSE小,故Bayes估計比MLE的精度更高。

2) 給定移走方案CS,屏蔽水平p越大,參數λ1,λ2,λ3及可度函數的MLE和BE的均方誤差越大;給定屏蔽水平p,采用移走方案3時,參數λ1,λ2,λ3及可度函數的MLE和BE的MSE相對較小。采用移走方案1與2, 對參數及可度函數的MLE和BE的MSE無顯著影響。

3) 當失效原因完全被屏蔽(p=1)時,無法計算參數λ1,λ2,λ3及可靠性指標的MLE,但可以獲得其BE,且精度較高。

4) 從各置信區間的長度及覆蓋率數據可以看出,各置信區間的覆蓋率接近于0.95,覆蓋率較大。隨著屏蔽水平的提高,置信區間長度增加覆蓋率減小。

4　結　語

本文研究了補充指數部件屏蔽數據并聯系統的可靠性估計問題。基于逐步Ⅱ型截尾試驗下屏蔽系統壽命數據,推導出試驗樣本的似然函數,進而獲得了部件參數及可靠度函數的MLE與BE。在此基礎上,利用M-H抽樣算法給出了部件參數的最大后驗密度置信區間。最后給出隨機模擬例子,并計算出各種估計的均方誤差及參數置信區間的長度及覆蓋率,結果表明MLE和BE的均方誤差都隨屏蔽水平的增大而增大,且Bayes估計效果更優。

參考文獻：

[1]USHER J S, HODGSON T J.Maximum likelihood analysis of component reliability using masked system life data [J].IEEE Transactins on Reliability, 1988, 37(5):550-555.

[2]SARHAN A M.Estimation of system components reliabilities using masked data [J].Applied Mathematics and Computation, 2003, 136:79-92.

[3]SARHAN A M, KUNDU D.Bayes estimators for reliability measures in geometric distribution model using masked system life test data [J].Computational Statistics & Data Analysis, 2008, 52(4):1821-1836.

[4]KUO L, YANG T Y. Bayesian reliability modeling for masked system lifetime data[J].Statistics & Probability Letters, 2000,47: 229-241.

[5]SARHAN A M, KUNDU D. Bayes estimators for reliability measures in geometric distribution model using masked system life test data[J].Computational Statistics and Data Analysis, 2008,52(4):1821-1836.

[6]SARHAN A M. Estimations of parameters in Pareto reliability model in the presence of masked data[J].Reliability Engineering and System Safety, 2003,82(1):75-83.

[7]CAI J, SHI Y, LIU B. Inference for a series system with dependent masked data under progressively interval censoring[J].Journal of Applied Statistics, 2017, 44(1): 3-15.

[8]FAN T H, HSU T M. Accelerated life tests of a series system with masked interval data under exponential lifetime distributions[J].IEEE Trans. on Reliability, 2012,61(3):798-808.

[9]CAI J, SHI Y, YUE H. Accelerated life tests for Log-normal series system with dependent masked data under Type-I progressive hybrid censoring[J].Communications in Statistics-Simulation and Computation, 2017, 46(2):1628-1646.

[10] ZHANG M, SHI Y, YANG Y. Reliability estimation for the generalized exponential components in parallel system using masked data[J].Information & Control, 2011, 40(4): 483-496.

[11] GU X, SHI Y. Estimation of generalized rayleigh components reliability in a parallel system using dependent masked data [J].Journal of Physical Science, 2011, 15: 33-42.

[12] AMMAR M, SARHAN A, AHMED H.Estimation of components reliability in a parallel system using masked system life data [J].Applied Mathematics and Computation, 2003, 138: 61-75.

[13] XU X, WANG Y, WANG R, et al.Statistical analysis of parallel system for masked data under successive censored life test [J].Journal of Shanghai Normal University (Natural Sciences),2016,45(3):378-385.

[14] XU X, WANG R, GU B.Statistical analysis of new model for masked data under parallel system [J].Journal of Sichuan Ordnance,2017,38(4)164-169.

[15] 李鳳,師義民,田亞愛.逐步增加II型截尾下Weibull分布的Bayes估計[J].工程數學學報, 2008, 25(4):641-650.

[16] FAN T H,HSU T M. Accelerated life tests of a series system with masked interval data under exponentiallifetime distributions[J].IEEE Transactions on Reliability, 2012,61(3):798-803.

[17] 茆詩松.貝葉斯統計[M].北京：中國統計出版社, 1999.

[18] MUKHOPADHYAY C, BASU A P.Bayesian analysis of incomplete time and cause of failure data [J].Journal of Statistical Planning and Inference, 1997, 59: 79-100.