基于MCMC的導(dǎo)彈多階段可靠性評估*

朱會傳 徐廷學(xué) 董 琪

(1.海軍航空工程學(xué)院研究生管理大隊 煙臺 264001)(2.海軍航空工程學(xué)院兵器科學(xué)與技術(shù)系 煙臺 264001)

?

基于MCMC的導(dǎo)彈多階段可靠性評估*

朱會傳1徐廷學(xué)2董 琪1

(1.海軍航空工程學(xué)院研究生管理大隊 煙臺 264001)(2.海軍航空工程學(xué)院兵器科學(xué)與技術(shù)系 煙臺 264001)

針對導(dǎo)彈在研制試驗中存在多階段、樣本量小而造成傳統(tǒng)評估誤差大的問題,建立了基于模糊理論、Bayes方法和馬爾可夫鏈蒙特卡羅法(MCMC)的可靠性評估模型。根據(jù)專家信息和模糊理論確定可靠性先驗分布;并利用Bayes公式得到可靠度后驗分布函數(shù);結(jié)合試驗多階段的特點,通過MCMC模型對不同階段導(dǎo)彈的可靠性進行評估,其中利用Gibbs抽樣算法對模型進行求解。通過算例分析,驗證了模型的有效性,為可靠性增長評估提供了新的思路。

可靠性評估; 仿真; 馬爾可夫; 蒙特卡羅

Class Number TP182

1 引言

導(dǎo)彈系統(tǒng)是復(fù)雜的武器系統(tǒng),對可靠性要求較高,對其進行可靠性評估是研制工作的客觀需求和重要內(nèi)容[1]。

由于導(dǎo)彈中主要部件像發(fā)動機和引信這類結(jié)構(gòu)復(fù)雜的系統(tǒng),研制成本較高,在試驗中會產(chǎn)生樣本量不大的專家信息和成敗型試驗數(shù)據(jù)。在對這樣的系統(tǒng)進行可靠性評估時,傳統(tǒng)的可靠性理論適用性不強。

Bayes方法在可靠性多源信息融合中有著廣泛的應(yīng)用。現(xiàn)有的Bayes可靠性增長模型主要為Smith模型[2]、Bar-Scheuer模型[3]、折合因子法[4]等,而且同樣運用Bayes方法,文獻[5]主要融合了產(chǎn)品驗前信息和現(xiàn)場信息進行進精度鑒定;文獻[6]則注重多維動態(tài)參數(shù)的融合來分析飛行器的試驗中的可靠性;文獻[7]綜合利用專家信息和各實驗階段的信息,并基于狄氏分布和可靠性增長來分析固體火箭發(fā)動機可靠性;文獻[8]主要利用性能退化數(shù)據(jù)進行可靠性評定;文獻[9]利用發(fā)動機運行數(shù)據(jù)和地面試驗數(shù)據(jù)進行可靠性分析。

不同的裝備性能的變化有著不同的規(guī)律,需要用不同的模型來進行分析。本文針對導(dǎo)彈在研制過程出現(xiàn)的試驗數(shù)據(jù)的特點,給出一種綜合利用模糊理論、Bayes方法和MCMC(Markov Chain Monte Carlo)方法,融合可靠性先驗信息和小樣本現(xiàn)場試驗信息的可靠性建模與分析方法,充分利用研制階段的多重信息,在很大程度上提高了可靠性評估的精度。本文研究思路如下:以導(dǎo)彈的關(guān)鍵部件發(fā)動機為例,首先根據(jù)利用模糊理論對專家信息進行賦值,確定先驗分布參數(shù),進而擬合出成敗型系統(tǒng)先驗分布函數(shù);然后結(jié)合Bayes、可靠度聯(lián)合先驗密度函數(shù)及可靠度似然函數(shù)得到可靠度后驗分布函數(shù);最后,利用MCMC算法中最有代表性的Gibbs抽樣算法對模型進行求解計算。

2 確定先驗參數(shù)

先驗參數(shù)的求解是確定先驗分布的關(guān)鍵,也是利用Bayes進行可靠性評估的難點和重點[10]。通常情況下,先驗參數(shù)的確定依據(jù)來自專家經(jīng)驗或者相似產(chǎn)品信息等,而且應(yīng)優(yōu)先參考相似產(chǎn)品信息,當(dāng)無相似產(chǎn)品信息時,再由專家經(jīng)驗給出先驗信息。但是由于導(dǎo)彈在研制階段,缺少相似產(chǎn)品的參考信息,所以采用專家經(jīng)驗來計算先驗信息。由于不同的專家給出的先驗信息的可信度不同,本文對不同專家給出的先驗信息附加不同的權(quán)重值,使獲得的先驗信息更科學(xué)、合理。

根據(jù)專家給出的可靠度區(qū)間,宜采用均勻分布進行描述,求出與該均勻分布最為接近的Beta分布,作為先驗分布。步驟如下:

然后,利用模糊理論對專家信息進行權(quán)重分配,可先構(gòu)造模糊判斷矩陣?？梢詮膶＜覐氖孪嚓P(guān)研究的時間、職稱、職務(wù)等方面綜合考慮,并采用0～1七標(biāo)度構(gòu)造互補模糊矩陣,0～1七標(biāo)度表示不同專家兩兩相對重要程度,含義如表1所示,其中X、Y表示兩個不同的專家。

表1 模糊判斷矩陣標(biāo)度及其含義

接著,計算相應(yīng)權(quán)重向量。方法如下:

1) 將模糊互補矩陣M(mik)l×l按式(6)求和

(1)

2) 通過式(7)得到模糊一致矩陣MT(tij)l×l

(2)

3) 由式(3)對矩陣MT(tij)l×l行和歸一化即可得到不同專家先驗信息的影響權(quán)重排序向量w=(w1,w2,…,wl)T。其中第i個專家的權(quán)重值為

(3)

最后,根據(jù)權(quán)重排序和l個專家的先驗信息,利用式(3)經(jīng)過加權(quán)綜合得到第k個試驗階段的可靠性先驗信息。

(4)

該試驗階段可靠性先驗信息的均值和方差為

(5)

令先驗信息均勻分布的均值和方差與先驗Beta分布的均值和方差等價,通過解方程可確定先驗參數(shù)值。因此需要求出在已知前k-1個試驗階段可靠度的條件下,第k個試驗階段可靠度條件均值和條件方差為

(6)

由于在實際計算得到的先驗參數(shù)有時會出現(xiàn)負值,不能滿足先驗參數(shù)都為正值的要求,針對這種情況,可建立最優(yōu)化模型,如式(7)～式(9),將均值作為約束,方差作為目標(biāo),確定第k個試驗階段最優(yōu)先驗參數(shù)。

目標(biāo)函數(shù)為

min(vk-VRk)2

(7)

約束條件為

μk=ERk

(8)

ak>0bk>0

(9)

其中式(6)中Rk-1可由第k-1個試驗階段的條件均值μk-1代替。

3 基于MCMC預(yù)測模型的建立

馬爾可夫鏈預(yù)測模型對未來狀態(tài)的預(yù)測,不需要尋求系統(tǒng)復(fù)雜因素之間的相互規(guī)律,只需考慮系統(tǒng)本身歷史狀態(tài)的演變特點,是一種動態(tài)隨機數(shù)學(xué)模型。

3.1 模型建立

由調(diào)研和分析知道,導(dǎo)彈在研制過程中,多是在出現(xiàn)故障后進行故障的維修,而在兩次故障期間,可靠性水平依次降低,該趨勢稱為序化約束模型。如式(10)所示:

0≤R1≤R2≤…≤Rk≤…≤Rm≤1

(10)

其中,Rk為第k個時間點的可靠性水平。

構(gòu)造第k個時間階段的可靠度先驗分布函數(shù),即fk(Rk|Rk-1) =fk(Rk|Rk-1;ak,bk)

×(1-Rk-1)1-ak-bkI(Rk-1,1)(Rk)

(11)

式中B(ak,bk)是先驗參數(shù)為ak和bk的Beta函數(shù),且ak>0,bk>0,令R0=0,Rm+1=1。

綜合各試驗階段的可靠度分布函數(shù)可得聯(lián)合先驗密度函數(shù)

(12)

式中R=(R1,R2,…,Rm)為各試驗階段可靠度估計值。

根據(jù)不同試驗階段的試驗數(shù)據(jù)可以確定可靠度似然函數(shù),若在第k個試驗階段現(xiàn)場數(shù)據(jù)為(nk,sk),其中nk表示第k個試驗階段的總試驗數(shù),sk表示第k個試驗階段的總成功數(shù)。則可靠度Rk的似然函數(shù)為

(13)

整個試驗階段可靠度R的似然函數(shù)為

(14)

3.2 可靠度后驗推斷

確定先驗參數(shù)后即可確定先驗分布函數(shù),由式(12)和式(14),結(jié)合Bayes公式可得R的后驗分布為

×(1-Rk-1)1-ak-bk

(15)

在計算后驗分布和后驗估計時,會面臨難以求解的高維積分問題。而MCMC方法已成為一種處理復(fù)雜統(tǒng)計問題的有效工具,更經(jīng)常應(yīng)用于需要復(fù)雜高維積分的Bayes分析領(lǐng)域。Gibbs抽樣是一種MCMC算法[11],它將隨機過程中的馬爾可夫鏈應(yīng)用到蒙特卡羅模擬中。本文采用Gibbs抽樣算法進行后驗積分和后驗估計的計算。利用Gibbs抽樣方法產(chǎn)生隨機樣本的步驟如下:

1) 抽樣時,可由式(15)經(jīng)過推導(dǎo)得R的后驗密度的核

×(1-Rk)bk+nk-sk-ak+1-bk+1

(16)

2) 確定完全條件分布。對于后驗密度的核,在給定Rk-1的情況下,fk(Rk|Rk-1;ak,bk)僅為Rk的分布函數(shù),因此,fk(Rk|Rk-1;ak,bk)顯然是參數(shù)Rk的完全條件分布。

…

將需要計算的后驗估計寫成函數(shù)φ(Ri)關(guān)于后驗分布fk(Rk|Rk-1;ak,bk)的期望:

E[φ(Ri)]=∫φ(Ri)fk(Rk|Rk-1;ak,bk)dRi

(17)

4) 判斷收斂性。模型是否收斂常根據(jù)迭代圖形來判斷,即同時輸入多組初始值經(jīng)過一段時間的Gibbs抽樣后,通過迭代圖形結(jié)果判斷,若迭代圖形趨于重合,則認為Gibbs抽樣收斂。

根據(jù)后驗分布式(17)得到完全條件分布為

×(1-Rk)bk+nk-sk-ak+1-bk+1

(18)

此時,對完全條件分布式(17)進行標(biāo)準(zhǔn)取舍抽樣,令

×(1-Rk)bk+nk-sk-ak+1-bk+1

(19)

g(Rk)∝(Rk-Rk-1)ak-1(1-Rk)ak+1-1

(20)

可見g(Rk)為區(qū)間(Rk-1,1)上的截尾Beta分布,易于從g(Rk)抽樣。

利用Gibbs算法得到R的抽樣值后,就可以利用這些抽樣值對Rk的后驗估計和后驗區(qū)間進行推斷。

4 實例分析

在導(dǎo)彈發(fā)動機研制過程中,使用人員制定了四個可靠性增長試驗階段。首先,由于未找到相似產(chǎn)品信息,現(xiàn)只提取專家信息。本次專家信息提取的調(diào)查對象,是對該裝備使用過程和可靠性較熟悉的專家,參加調(diào)查的共有五位專家,分別來自研制部門、基層部隊和相關(guān)院校。在職稱方面,高級工程師三名,工程師二名;學(xué)歷方面,博士二名,碩士三名。工作年限,高級工程師一般在21～30年,工程師一般在11年～20年。

本次專家信息的提取方式采用調(diào)查問卷方式。其中,專家信息主要是統(tǒng)計學(xué)描述方式,即置信區(qū)間形式,如表2所示。

表2 發(fā)動機可靠度專家信息調(diào)查表

本次調(diào)查表一共發(fā)放20份,收回有效問卷20份,下面對專家信息進行處理。

首先,根據(jù)專家信息調(diào)查表構(gòu)造模糊互補矩陣M:

根據(jù)式(2)得到模糊一致矩陣MT:

由式(3)對矩陣MT行和歸一化得到專家先驗信息影響權(quán)重排序向量w=(0.2188,0.2188,0.1875,0.1875,0.1875)T。

其次,對專家給出的可靠性信息進行加權(quán)得到不同階段可靠性先驗信息,并結(jié)合試驗數(shù)據(jù)如表3所示。

表3 發(fā)動機先驗信息和試驗數(shù)據(jù)

再次,利用式(7)～式(9)可以得到各試驗階段與均勻分布等價的Beta分布的先驗參數(shù)ak和bk,如表4所示。

表4 Beta分布先驗參數(shù)

確定先驗參數(shù)后,即可將試驗數(shù)據(jù)與先驗參數(shù)代入式(16),利用Gibbs抽樣計算導(dǎo)彈發(fā)動機可靠度,采用三個不同的初始值并產(chǎn)生三條Markov鏈,迭代次數(shù)取10000次,然后觀測三條Markov鏈的抽樣值,并判斷收斂性。

圖1 R[4]的三條Markov鏈的抽樣值

如圖1所示,第四試驗階段發(fā)動機可靠度R[4]的三條Markov鏈迭代圖形結(jié)果趨于重合,說明抽樣值隨著迭代次數(shù)的增加趨于穩(wěn)定。圖2～圖4顯示了R[4]的迭代軌跡、分位數(shù)、核密度等抽樣信息,而如圖5所示,R[4]的自相關(guān)函數(shù)很快接近于0,也證明迭代過程已收斂。

圖2 R[4]的迭代軌跡圖

圖3 R[4]的分位數(shù)

圖4 R[4]的核密度函數(shù)圖

圖5 R[4]的自相關(guān)函數(shù)

通過Gibbs抽樣,得到各試驗階段可靠性的后驗估計結(jié)果,如表5所示。

表5 發(fā)動機各試驗階段可靠度的后驗估計

若僅利用第四階段的試驗數(shù)據(jù),得到的可靠度估計值為0.9286,略低于0.93的要求值,且受樣本量影響較大。由表5可見,導(dǎo)彈發(fā)動機在試驗結(jié)束后的可靠度估計值為0.9433,基本滿足了0.93的可靠性要求,這是因為本文采用的方法利用先驗信息以及前階段的可靠性信息,利用的可靠性信息更多,得出的結(jié)果更真實可靠。

5 結(jié)語

針對導(dǎo)彈研制過程成敗型數(shù)據(jù)樣本量小的特點,給出一種綜合利用模糊理論、Bayes方法和MCMC方法的可靠性建模與分析方法,充分利用了可靠性先驗信息和小樣本現(xiàn)場試驗信息,在一定程度上提高了可靠性的計算精度。實例分析表明,該方法對研制階段導(dǎo)彈的可靠性評估是有效和可行的。

[1] 周忠寶,厲海濤,劉學(xué)敏,等.航天長壽命產(chǎn)品可靠性建模與評估的Bayes信息融合方法[J].系統(tǒng)工程理論與實踐,2012,32(7):2517-2522.

[2] Smith A F M. A Bayesian notes on reliability growth during a development testing program[J]. IEEE Transactions on reliability,1977,(26):346-347.

[3] Barlow R E, Scheuer E M. reliability growth during a development testing program[J]. Technometrics,1966,(8):53-60.

[4] 田國梁.二項分布的可靠性增長模型[J].宇航學(xué)報,1992,(1):55-61.

[5] 張金槐.多源信息的Bayes融合精度鑒定方法[J].國防科技大學(xué)學(xué)報,2001,23(3):93-97.

[6] 張金槐.多維動態(tài)參數(shù)的多層Bayes融合設(shè)計[J].國防科技大學(xué)學(xué)報,2003,25(4):97-101.

[7] 劉飛,竇毅芳,張為華.基于狄氏先驗分布的固體火箭發(fā)動機可靠性增長Bayes分析[J].固體火箭技術(shù),2006,29(4):239-242.

[8] 馮靜,劉琦,周經(jīng)倫,等.基于性能退化數(shù)據(jù)的液體火箭發(fā)動機可靠性Bayes分析[J].航空計算技術(shù),2003,33(3):6-10.

[9] 劉琦,馮靜,周經(jīng)倫.火箭發(fā)動機性能可靠性的Bayes分析[J].中國空間科學(xué)技術(shù),2003,23(6):46-51.

[10] 白璐.可靠性評估中的一種Bayes方法[J].質(zhì)量與可靠性,2007,129(3):20-23.

[11] 張云安.小樣本復(fù)雜系統(tǒng)的多階段可靠性增長評估與可靠性鑒定[D].長沙:國防科技大學(xué),2008:100-102.

Multi-stage Reliability Assessment of Missile Based on MCMC

ZHU Huichuan1XU Tingxue2DONG Qi1

(1. Graduate Students’ Brigade, Naval Aeronartical and Astronautical University, Yantai 264001)(2. Department of Weapon Science and Technology, Naval Aeronartical and Astronautical University, Yantai 264001)

The reliability evaluation model based on fuzzy theory, Bayesian analysis method and MCMC is established for solving the problem of the small sample of missile reliability and mulit-stage during the development of the missile. According to the expert information and fuzzy theory, the prior distribution of reliability is determined. The posterior distribution is given by making use of the Bayesian formula. The reliability assessment of the every stage of the missile equipment is verified by sampling simulation of Markov Chain Monte Carlo and Gibbs Calculate method. Finally, an example is used to illustrate the effective of the model, and provide a new method for the reliability-growth assessment.

reliability assessment, simulation, Markov, Monte Carlo

2014年9月1日,

2014年10月19日

朱會傳,男,碩士研究生,研究方向:武器裝備綜合保障理論與技術(shù)。徐廷學(xué),男,博士,教授,研究方向:武器裝備綜合保障理論與技術(shù)。董琪,男,博士研究生,研究方向:武器裝備綜合保障理論與技術(shù)。

TP182

10.3969/j.issn1672-9730.2015.03.027