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

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

?

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

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

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

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

可靠性評(píng)估; 仿真; 馬爾可夫; 蒙特卡羅

Class Number TP182

1 引言

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

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

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

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

2 確定先驗(yàn)參數(shù)

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

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

然后,利用模糊理論對(duì)專家信息進(jìn)行權(quán)重分配,可先構(gòu)造模糊判斷矩陣。可以從專家從事相關(guān)研究的時(shí)間、職稱、職務(wù)等方面綜合考慮,并采用0～1七標(biāo)度構(gòu)造互補(bǔ)模糊矩陣,0～1七標(biāo)度表示不同專家兩兩相對(duì)重要程度,含義如表1所示,其中X、Y表示兩個(gè)不同的專家。

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

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

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

(1)

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

(2)

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

(3)

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

(4)

該試驗(yàn)階段可靠性先驗(yàn)信息的均值和方差為

(5)

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

(6)

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

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

min(vk-VRk)2

(7)

約束條件為

μk=ERk

(8)

ak>0bk>0

(9)

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

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

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

3.1 模型建立

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

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

(10)

其中,Rk為第k個(gè)時(shí)間點(diǎn)的可靠性水平。

構(gòu)造第k個(gè)時(shí)間階段的可靠度先驗(yàn)分布函數(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)是先驗(yàn)參數(shù)為ak和bk的Beta函數(shù),且ak>0,bk>0,令R0=0,Rm+1=1。

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

(12)

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

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

(13)

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

(14)

3.2 可靠度后驗(yàn)推斷

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

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

(15)

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

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

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

(16)

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

…

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

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

(17)

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

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

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

(18)

此時(shí),對(duì)完全條件分布式(17)進(jìn)行標(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的抽樣值后,就可以利用這些抽樣值對(duì)Rk的后驗(yàn)估計(jì)和后驗(yàn)區(qū)間進(jìn)行推斷。

4 實(shí)例分析

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

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

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

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

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

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

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

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

表3 發(fā)動(dòng)機(jī)先驗(yàn)信息和試驗(yàn)數(shù)據(jù)

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

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

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

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

如圖1所示,第四試驗(yàn)階段發(fā)動(dòng)機(jī)可靠度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抽樣,得到各試驗(yàn)階段可靠性的后驗(yàn)估計(jì)結(jié)果,如表5所示。

表5 發(fā)動(dòng)機(jī)各試驗(yàn)階段可靠度的后驗(yàn)估計(jì)

若僅利用第四階段的試驗(yàn)數(shù)據(jù),得到的可靠度估計(jì)值為0.9286,略低于0.93的要求值,且受樣本量影響較大。由表5可見,導(dǎo)彈發(fā)動(dòng)機(jī)在試驗(yàn)結(jié)束后的可靠度估計(jì)值為0.9433,基本滿足了0.93的可靠性要求,這是因?yàn)楸疚牟捎玫姆椒ɡ孟闰?yàn)信息以及前階段的可靠性信息,利用的可靠性信息更多,得出的結(jié)果更真實(shí)可靠。

5 結(jié)語(yǔ)

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

[1] 周忠寶,厲海濤,劉學(xué)敏,等.航天長(zhǎng)壽命產(chǎn)品可靠性建模與評(píng)估的Bayes信息融合方法[J].系統(tǒng)工程理論與實(shí)踐,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] 田國(guó)梁.二項(xiàng)分布的可靠性增長(zhǎng)模型[J].宇航學(xué)報(bào),1992,(1):55-61.

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

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

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

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

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

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

[11] 張?jiān)瓢?小樣本復(fù)雜系統(tǒng)的多階段可靠性增長(zhǎng)評(píng)估與可靠性鑒定[D].長(zhǎng)沙:國(guó)防科技大學(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日

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

TP182

10.3969/j.issn1672-9730.2015.03.027