朱會傳 徐廷學 董 琪
(1.海軍航空工程學院研究生管理大隊 煙臺 264001)(2.海軍航空工程學院兵器科學與技術系 煙臺 264001)
?
基于MCMC的導彈多階段可靠性評估*
朱會傳1徐廷學2董 琪1
(1.海軍航空工程學院研究生管理大隊 煙臺 264001)(2.海軍航空工程學院兵器科學與技術系 煙臺 264001)
針對導彈在研制試驗中存在多階段、樣本量小而造成傳統(tǒng)評估誤差大的問題,建立了基于模糊理論、Bayes方法和馬爾可夫鏈蒙特卡羅法(MCMC)的可靠性評估模型。根據專家信息和模糊理論確定可靠性先驗分布;并利用Bayes公式得到可靠度后驗分布函數;結合試驗多階段的特點,通過MCMC模型對不同階段導彈的可靠性進行評估,其中利用Gibbs抽樣算法對模型進行求解。通過算例分析,驗證了模型的有效性,為可靠性增長評估提供了新的思路。
可靠性評估; 仿真; 馬爾可夫; 蒙特卡羅
Class Number TP182
導彈系統(tǒng)是復雜的武器系統(tǒng),對可靠性要求較高,對其進行可靠性評估是研制工作的客觀需求和重要內容[1]。
由于導彈中主要部件像發(fā)動機和引信這類結構復雜的系統(tǒng),研制成本較高,在試驗中會產生樣本量不大的專家信息和成敗型試驗數據。在對這樣的系統(tǒng)進行可靠性評估時,傳統(tǒng)的可靠性理論適用性不強。
Bayes方法在可靠性多源信息融合中有著廣泛的應用。現有的Bayes可靠性增長模型主要為Smith模型[2]、Bar-Scheuer模型[3]、折合因子法[4]等,而且同樣運用Bayes方法,文獻[5]主要融合了產品驗前信息和現場信息進行進精度鑒定;文獻[6]則注重多維動態(tài)參數的融合來分析飛行器的試驗中的可靠性;文獻[7]綜合利用專家信息和各實驗階段的信息,并基于狄氏分布和可靠性增長來分析固體火箭發(fā)動機可靠性;文獻[8]主要利用性能退化數據進行可靠性評定;文獻[9]利用發(fā)動機運行數據和地面試驗數據進行可靠性分析。
不同的裝備性能的變化有著不同的規(guī)律,需要用不同的模型來進行分析。本文針對導彈在研制過程出現的試驗數據的特點,給出一種綜合利用模糊理論、Bayes方法和MCMC(Markov Chain Monte Carlo)方法,融合可靠性先驗信息和小樣本現場試驗信息的可靠性建模與分析方法,充分利用研制階段的多重信息,在很大程度上提高了可靠性評估的精度。本文研究思路如下:以導彈的關鍵部件發(fā)動機為例,首先根據利用模糊理論對專家信息進行賦值,確定先驗分布參數,進而擬合出成敗型系統(tǒng)先驗分布函數;然后結合Bayes、可靠度聯合先驗密度函數及可靠度似然函數得到可靠度后驗分布函數;最后,利用MCMC算法中最有代表性的Gibbs抽樣算法對模型進行求解計算。
先驗參數的求解是確定先驗分布的關鍵,也是利用Bayes進行可靠性評估的難點和重點[10]。通常情況下,先驗參數的確定依據來自專家經驗或者相似產品信息等,而且應優(yōu)先參考相似產品信息,當無相似產品信息時,再由專家經驗給出先驗信息。但是由于導彈在研制階段,缺少相似產品的參考信息,所以采用專家經驗來計算先驗信息。由于不同的專家給出的先驗信息的可信度不同,本文對不同專家給出的先驗信息附加不同的權重值,使獲得的先驗信息更科學、合理。
根據專家給出的可靠度區(qū)間,宜采用均勻分布進行描述,求出與該均勻分布最為接近的Beta分布,作為先驗分布。步驟如下:
然后,利用模糊理論對專家信息進行權重分配,可先構造模糊判斷矩陣??梢詮膶<覐氖孪嚓P研究的時間、職稱、職務等方面綜合考慮,并采用0~1七標度構造互補模糊矩陣,0~1七標度表示不同專家兩兩相對重要程度,含義如表1所示,其中X、Y表示兩個不同的專家。
表1 模糊判斷矩陣標度及其含義
接著,計算相應權重向量。方法如下:
1) 將模糊互補矩陣M(mik)l×l按式(6)求和
(1)
2) 通過式(7)得到模糊一致矩陣MT(tij)l×l
(2)
3) 由式(3)對矩陣MT(tij)l×l行和歸一化即可得到不同專家先驗信息的影響權重排序向量w=(w1,w2,…,wl)T。其中第i個專家的權重值為
(3)
最后,根據權重排序和l個專家的先驗信息,利用式(3)經過加權綜合得到第k個試驗階段的可靠性先驗信息。
(4)
該試驗階段可靠性先驗信息的均值和方差為
(5)
令先驗信息均勻分布的均值和方差與先驗Beta分布的均值和方差等價,通過解方程可確定先驗參數值。因此需要求出在已知前k-1個試驗階段可靠度的條件下,第k個試驗階段可靠度條件均值和條件方差為
(6)
由于在實際計算得到的先驗參數有時會出現負值,不能滿足先驗參數都為正值的要求,針對這種情況,可建立最優(yōu)化模型,如式(7)~式(9),將均值作為約束,方差作為目標,確定第k個試驗階段最優(yōu)先驗參數。
目標函數為
min(vk-VRk)2
(7)
約束條件為
μk=ERk
(8)
ak>0bk>0
(9)
其中式(6)中Rk-1可由第k-1個試驗階段的條件均值μk-1代替。
馬爾可夫鏈預測模型對未來狀態(tài)的預測,不需要尋求系統(tǒng)復雜因素之間的相互規(guī)律,只需考慮系統(tǒng)本身歷史狀態(tài)的演變特點,是一種動態(tài)隨機數學模型。
3.1 模型建立
由調研和分析知道,導彈在研制過程中,多是在出現故障后進行故障的維修,而在兩次故障期間,可靠性水平依次降低,該趨勢稱為序化約束模型。如式(10)所示:
0≤R1≤R2≤…≤Rk≤…≤Rm≤1
(10)
其中,Rk為第k個時間點的可靠性水平。
構造第k個時間階段的可靠度先驗分布函數,即fk(Rk|Rk-1) =fk(Rk|Rk-1;ak,bk)
×(1-Rk-1)1-ak-bkI(Rk-1,1)(Rk)
(11)
式中B(ak,bk)是先驗參數為ak和bk的Beta函數,且ak>0,bk>0,令R0=0,Rm+1=1。
綜合各試驗階段的可靠度分布函數可得聯合先驗密度函數
(12)
式中R=(R1,R2,…,Rm)為各試驗階段可靠度估計值。
根據不同試驗階段的試驗數據可以確定可靠度似然函數,若在第k個試驗階段現場數據為(nk,sk),其中nk表示第k個試驗階段的總試驗數,sk表示第k個試驗階段的總成功數。則可靠度Rk的似然函數為
(13)
整個試驗階段可靠度R的似然函數為
(14)
3.2 可靠度后驗推斷
確定先驗參數后即可確定先驗分布函數,由式(12)和式(14),結合Bayes公式可得R的后驗分布為
×(1-Rk-1)1-ak-bk
(15)
在計算后驗分布和后驗估計時,會面臨難以求解的高維積分問題。而MCMC方法已成為一種處理復雜統(tǒng)計問題的有效工具,更經常應用于需要復雜高維積分的Bayes分析領域。Gibbs抽樣是一種MCMC算法[11],它將隨機過程中的馬爾可夫鏈應用到蒙特卡羅模擬中。本文采用Gibbs抽樣算法進行后驗積分和后驗估計的計算。利用Gibbs抽樣方法產生隨機樣本的步驟如下:
1) 抽樣時,可由式(15)經過推導得R的后驗密度的核
×(1-Rk)bk+nk-sk-ak+1-bk+1
(16)
2) 確定完全條件分布。對于后驗密度的核,在給定Rk-1的情況下,fk(Rk|Rk-1;ak,bk)僅為Rk的分布函數,因此,fk(Rk|Rk-1;ak,bk)顯然是參數Rk的完全條件分布。
…
將需要計算的后驗估計寫成函數φ(Ri)關于后驗分布fk(Rk|Rk-1;ak,bk)的期望:
E[φ(Ri)]=∫φ(Ri)fk(Rk|Rk-1;ak,bk)dRi
(17)
4) 判斷收斂性。模型是否收斂常根據迭代圖形來判斷,即同時輸入多組初始值經過一段時間的Gibbs抽樣后,通過迭代圖形結果判斷,若迭代圖形趨于重合,則認為Gibbs抽樣收斂。
根據后驗分布式(17)得到完全條件分布為
×(1-Rk)bk+nk-sk-ak+1-bk+1
(18)
此時,對完全條件分布式(17)進行標準取舍抽樣,令
×(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ū)間進行推斷。
在導彈發(fā)動機研制過程中,使用人員制定了四個可靠性增長試驗階段。首先,由于未找到相似產品信息,現只提取專家信息。本次專家信息提取的調查對象,是對該裝備使用過程和可靠性較熟悉的專家,參加調查的共有五位專家,分別來自研制部門、基層部隊和相關院校。在職稱方面,高級工程師三名,工程師二名;學歷方面,博士二名,碩士三名。工作年限,高級工程師一般在21~30年,工程師一般在11年~20年。
本次專家信息的提取方式采用調查問卷方式。其中,專家信息主要是統(tǒng)計學描述方式,即置信區(qū)間形式,如表2所示。
表2 發(fā)動機可靠度專家信息調查表
本次調查表一共發(fā)放20份,收回有效問卷20份,下面對專家信息進行處理。
首先,根據專家信息調查表構造模糊互補矩陣M:
根據式(2)得到模糊一致矩陣MT:
由式(3)對矩陣MT行和歸一化得到專家先驗信息影響權重排序向量w=(0.2188,0.2188,0.1875,0.1875,0.1875)T。
其次,對專家給出的可靠性信息進行加權得到不同階段可靠性先驗信息,并結合試驗數據如表3所示。
表3 發(fā)動機先驗信息和試驗數據
再次,利用式(7)~式(9)可以得到各試驗階段與均勻分布等價的Beta分布的先驗參數ak和bk,如表4所示。
表4 Beta分布先驗參數
確定先驗參數后,即可將試驗數據與先驗參數代入式(16),利用Gibbs抽樣計算導彈發(fā)動機可靠度,采用三個不同的初始值并產生三條Markov鏈,迭代次數取10000次,然后觀測三條Markov鏈的抽樣值,并判斷收斂性。
圖1 R[4]的三條Markov鏈的抽樣值
如圖1所示,第四試驗階段發(fā)動機可靠度R[4]的三條Markov鏈迭代圖形結果趨于重合,說明抽樣值隨著迭代次數的增加趨于穩(wěn)定。圖2~圖4顯示了R[4]的迭代軌跡、分位數、核密度等抽樣信息,而如圖5所示,R[4]的自相關函數很快接近于0,也證明迭代過程已收斂。
圖2 R[4]的迭代軌跡圖
圖3 R[4]的分位數
圖4 R[4]的核密度函數圖
圖5 R[4]的自相關函數
通過Gibbs抽樣,得到各試驗階段可靠性的后驗估計結果,如表5所示。
表5 發(fā)動機各試驗階段可靠度的后驗估計
若僅利用第四階段的試驗數據,得到的可靠度估計值為0.9286,略低于0.93的要求值,且受樣本量影響較大。由表5可見,導彈發(fā)動機在試驗結束后的可靠度估計值為0.9433,基本滿足了0.93的可靠性要求,這是因為本文采用的方法利用先驗信息以及前階段的可靠性信息,利用的可靠性信息更多,得出的結果更真實可靠。
針對導彈研制過程成敗型數據樣本量小的特點,給出一種綜合利用模糊理論、Bayes方法和MCMC方法的可靠性建模與分析方法,充分利用了可靠性先驗信息和小樣本現場試驗信息,在一定程度上提高了可靠性的計算精度。實例分析表明,該方法對研制階段導彈的可靠性評估是有效和可行的。
[1] 周忠寶,厲海濤,劉學敏,等.航天長壽命產品可靠性建模與評估的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].宇航學報,1992,(1):55-61.
[5] 張金槐.多源信息的Bayes融合精度鑒定方法[J].國防科技大學學報,2001,23(3):93-97.
[6] 張金槐.多維動態(tài)參數的多層Bayes融合設計[J].國防科技大學學報,2003,25(4):97-101.
[7] 劉飛,竇毅芳,張為華.基于狄氏先驗分布的固體火箭發(fā)動機可靠性增長Bayes分析[J].固體火箭技術,2006,29(4):239-242.
[8] 馮靜,劉琦,周經倫,等.基于性能退化數據的液體火箭發(fā)動機可靠性Bayes分析[J].航空計算技術,2003,33(3):6-10.
[9] 劉琦,馮靜,周經倫.火箭發(fā)動機性能可靠性的Bayes分析[J].中國空間科學技術,2003,23(6):46-51.
[10] 白璐.可靠性評估中的一種Bayes方法[J].質量與可靠性,2007,129(3):20-23.
[11] 張云安.小樣本復雜系統(tǒng)的多階段可靠性增長評估與可靠性鑒定[D].長沙:國防科技大學,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日
朱會傳,男,碩士研究生,研究方向:武器裝備綜合保障理論與技術。徐廷學,男,博士,教授,研究方向:武器裝備綜合保障理論與技術。董琪,男,博士研究生,研究方向:武器裝備綜合保障理論與技術。
TP182
10.3969/j.issn1672-9730.2015.03.027