基于貝葉斯法和蒙特卡洛仿真的威布爾型裝備器材需求預測

吳龍濤， 王鐵寧， 楊帆

(陸軍裝甲兵學院 技術保障工程系， 北京 100072)

基于貝葉斯法和蒙特卡洛仿真的威布爾型裝備器材需求預測

吳龍濤， 王鐵寧， 楊帆

(陸軍裝甲兵學院 技術保障工程系， 北京 100072)

為了解決裝備器材歷史需求數據少、需求規律不明確的問題，提出一種基于貝葉斯法和蒙特卡洛仿真的威布爾型裝備器材需求預測方法。針對威布爾分布尺度參數未知以及形狀和尺度參數均未知兩種情況，分別基于貝葉斯方法通過解析求解和數值模擬的方式進行了參數估計，并引入柯爾莫哥洛夫- 斯米爾諾夫檢驗法對壽命分布模型進行擬合優度檢驗；綜合考慮修復性維修、預防性維修和裝備器材的已使用時間，提出了基于蒙特卡洛仿真的部隊裝備器材年度需求預測方法。算例分析表明：小樣本下通過貝葉斯估計得到的壽命分布模型擬合度高，基于仿真的需求預測方法簡單、有效。

兵器科學與技術； 裝備器材需求預測； 威布爾分布； 貝葉斯法； 蒙特卡洛仿真

0 引言

裝備器材需求數量的確定是裝備器材保障工作的一個重要環節。隨著部隊高新裝備列裝步伐的加快和換件修理方式的開展，裝備器材需求結構正逐漸發生變化。然而，由于裝備列裝時間短、裝備器材歷史需求數據少，需求規律一時難以掌握，給裝備器材保障工作造成了一定的困難。

我軍目前主要采用定額計算[1]的方法管理庫存，各級裝備器材保障部門通常依據庫存標準并結合現有庫存來制定需求計劃。但由于庫存標準本身是粗略、不精確的估計[2]，定額計算法難以精確計算裝備器材的需求量。同時，由于可以利用的需求數據較少，傳統的基于大樣本方法如時間序列分析、回歸分析和神經網絡等也難以發揮作用。灰色預測模型[3]是小樣本預測經常采用的方法，但其預測精度有限、適用范圍較小。

基于可靠性的需求預測方法可以精確計算裝備器材的需求量，此類方法根據失效機理選擇裝備器材壽命所服從的分布類型[4]，并以此建立裝備器材需求計算模型。但承制單位在交付裝備時往往難以提供裝備器材的壽命分布，且不同部隊的裝備使用和維護水平以及使用環境不同，裝備器材實際的使用可靠性也各不相同。因此，在需求預測前，首先要確定裝備器材的壽命分布形式。威布爾分布是可靠性分析中廣泛應用的壽命分布方法，具有明顯耗損期的裝備器材壽命大都服從威布爾分布。針對威布爾型裝備器材的需求預測，目前常見的模型大多是通過復雜的卷積計算得到一個多重無窮級數[5]，或以指數分布等效代替威布爾分布以簡化計算[6-7]。前一種方法解析計算困難，不便于工程應用；后一種方法存在較大的誤差；而且在計算裝備器材年度需求時，大都未考慮裝備器材上一年度的已使用時間，從而會產生一定的誤差。本文圍繞上述問題，首先給出了小樣本條件下基于貝葉斯法的威布爾分布單參數估計和雙參數估計方法；然后結合實際，綜合考慮修復性維修和預防性維修以及裝備器材的已使用時間，提出了基于蒙特卡洛仿真的威布爾型裝備器材年度需求預測方法，并進行了算例驗證和分析。

1 威布爾分布參數估計

設裝備器材壽命T服從威布爾分布，其概率密度函數、累積分布函數和裝備器材的可靠度函數分別為

f(t)=mλtm-1exp(-λtm)，

(1)

F(t)=1-exp(-λtm)，

(2)

R(t)=exp(-λtm)，

(3)

式中：m為形狀參數；λ為尺度參數；t為失效時間；m、λ和t均大于0. 當m<1時，失效率隨時間遞減，適用于裝備列裝的早期階段；當m>1時，失效率隨時間遞增，適用于裝備器材的耗損階段；當m=1時，威布爾分布退化為失效率恒定的指數分布。由此可見，威布爾分布具有很好的通用性。

常見的威布爾分布參數估計方法包括圖估計法、最小二乘估計法[5]、矩估計法和極大似然估計法[8]等。其中：前3種方法都是粗略的估計，精度較低；極大似然估計法通過迭代法求解超越方程，在大樣本下估計精度較高，但有時不易收斂，且面對截尾小樣本時偏差較大。貝葉斯法[9]可以充分利用專家經驗、廠家數據等先驗信息提高參數估計的準確性；同時，隨著使用過程中新信息的出現不斷更新結果，是解決小樣本參數估計的有效方法。

裝備器材在使用過程中被更換的原因大致有兩個：一是因故障失效被更換，二是在預防性維修時被更換。據此，可將裝備器材的更換時間樣本視為裝備器材壽命的定時截尾樣本，截尾時間即為預防性維修間隔期。設Ω為來自同一威布爾分布的定時截尾樣本，容量為n，截尾時間為t0，t1≤t2≤…≤tr為r個失效時間。下面對定時截尾小樣本下的威布爾分布參數m和λ進行估計。

1.1 參數m已知而參數λ未知

通常，形狀參數m可由專家經驗或相似產品法[10]直接給出，故只對尺度參數λ進行估計。運用貝葉斯法進行分布參數估計時，首先要為參數選擇合適的先驗分布。通常可以選擇以下兩種先驗分布：無信息先驗分布和有信息的共軛先驗分布。

1.1.1 無信息先驗分布

當待估參數的先驗信息很少時，一般選擇在參數空間內分布較為廣泛的無信息先驗分布。根據Jeffreys法則[11]，尺度參數λ的無信息先驗分布為

π(λ)=λ-1.

(4)

由(1)式和(3)式，在尺度參數λ給定的條件下樣本Ω的聯合分布密度函數為

(5)

(6)

式中：Γ(r)為伽馬函數。由(6)式可知，λ服從伽馬分布，即λ|Ω～Ga(r,τ).

用損失函數L(λ,)表示尺度參數真值為λ時，選擇作為估計值所造成的損失。平方損失L(λ,)=(λ-)2是一種常用的損失函數，平方損失最小時參數的貝葉斯估計值為后驗分布期望[9]。為了便于計算，通常也令η=λ-1/m表示尺度參數，由 (6) 式可得，參數η的貝葉斯估計值為

=1/m=1/m.

(7)

1.1.2 共軛先驗分布

共軛先驗分布是指與后驗分布屬于同一類分布的先驗分布。威布爾分布參數λ常用的共軛先驗分布是伽馬分布λ～Ga(α，β)，其概率密度函數為

(8)

式中：α>0，β>0為先驗分布的超參數，其值通常根據先驗信息使用矩估計法確定。將(8)式代入(6)式可知，λ的后驗分布仍為伽馬分布，且λ～Ga(α+r，β+τ). 則平方損失最小時參數η的貝葉斯估計值為

=1/m=1/m.

(9)

1.2 參數m和參數λ未知

當形狀參數和尺度參數均未知時，威布爾分布不存在適用的共軛先驗分布，因此暫且只討論無信息先驗分布。根據Jeffreys法則，參數m和參數λ的無信息先驗分布為

π(m,λ)=(mλ)-1.

(10)

實際中，通常可以根據經驗給定參數m的取值范圍[m1,m2]. 將(10)式代入(6)式，通過計算可得到參數m和參數λ的聯合后驗分布為

(11)

在此基礎上分別對λ和m積分，即可得到參數的后驗密度函數。基于平方損失最小原則，可得參數m和參數λ的貝葉斯估計值分別為

=,

(12)

=.

(13)

由于 (12) 式和(13)式難以直接解析求解，本文采用馬爾可夫鏈- 蒙特卡洛(MCMC)數值模擬方法中的隨機游走Metropolis-Hastings算法[12]求解。參數m候選點m*的建議密度函數[9]為區間[m1,m2]上的均勻分布，參數η候選點η*的建議密度函數為

η*=η(j-1)exp(v),v～N(0,σ2)，

(14)

式中：j為當前模擬次數；η(j)為η的第j次模擬值；σ為常數。候選點m*和η*的接收概率[8]分別為

(15)

式中：m(j)為m的第j次模擬值，j=0時的初始值m(0)和n(0)可根據實際情況自行設定。

1.3 參數估計檢驗

確定裝備器材壽命分布的具體形式后，需要通過壽命分布的擬合優度檢驗來驗證分布模型是否適用于樣本。χ2檢驗是一種常用的擬合優度檢驗方法，但只適用于大樣本完全數據。針對定時截尾小樣本的情況，本文引入了柯爾莫哥洛夫- 斯米爾諾夫(K-S)檢驗法[13]。

設裝備器材壽命T的真實分布函數為Ψ(t)，則檢驗假設為

H0∶Ψ(t)=F(t)?H1∶Ψ(t)≠F(t).

K-S檢驗統計量為

(16)

式中：φ(t)為壽命T關于給定失效樣本的經驗分布函數，其計算公式為

(17)

2 裝備器材需求預測

實際中部隊裝備器材需求發生的一般過程如圖1所示。假設某部隊有某型裝備z臺，每臺裝備安裝有a器材1件。根據裝備滾動式循環動用原則，不同裝備a器材的年初已使用時間T′i(i=1,2,…,z)和年度計劃使用時間Ti一般不同。為了便于分析計算，作如下假設：裝備故障失效時采取換件修理，換件時間忽略不計；所有a器材的壽命相互獨立，且服從同一威布爾分布。

由文獻[5，7]可知，多裝備同時動用時的威布爾型裝備器材需求模型是一個多重卷積計算過程，且未考慮裝備器材的已使用時間和預防性維修需求。下面給出基于蒙特卡洛仿真的裝備器材需求預測方法。

2.1 修復性維修需求仿真

實行換件修理的裝備器材故障間隔時間即為其使用壽命，可根據其壽命分布模型隨機確定。采用逆變法[14]生成威布爾分布的隨機故障間隔時間為

(18)

式中：μ為(0, 1)區間內服從均勻分布的隨機數。然而，與指數分布的“無后效性”不同，威布爾型裝備器材的故障間隔時間與已使用時間是相關的。通常，已使用時間越長，發生故障的可能性越大。因此，特別是對于可靠性較高的裝備器材，在計算需求時需要考慮其已使用時間。對于任一臺裝備i的a器材，給定已使用時間T′i的條件下，其剩余壽命分布函數為

(19)

由此可得其剩余壽命隨機數生成公式為

(20)

采用事件驅動機制模擬裝備器材故障更換需求過程，具體步驟如下：對于裝備i，首先利用 (20) 式生成隨機時間ΔTi1，然后利用(18)式連續生成隨機時間ΔTi2，ΔTi3,…，按此步驟模擬裝備i上a器材“使用- 故障- 更換”的循環過程，直到前k次故障累積時間ΔTi1+ΔTi2+…+ΔTik首次大于裝備i的年度使用時間Ti，則裝備i對a器材的年度需求N(Ti)=k-1. 對每臺裝備重復上述過程，可得z臺裝備的a器材年度故障更換需求量Xf=N(T1)+N(T2)+…+N(Tz).

2.2 預防性維修需求仿真

在需求預測時，除修復性維修需求外，還需要考慮預防性維修需求。根據裝備維修計劃，一般每年都會安排已到達定時維修間隔期的裝備進行預防性維修。在預防性維修時，可能對a器材采取不更換、必須更換或視情更換3種維修策略[2]，且對于視情更換，通常指定一個更換概率p. 設在z臺裝備中有z0臺需要進行等級維修：若a器材不更換，則其定時更換需求量Xp=0；若a器材必須更換，則Xp=z0；若a器材視情更換，則Xp服從二項分布B(z0,p)，可通過生成二項分布隨機數的方式確定Xp的值。

2.3 基于蒙特卡洛仿真的需求預測

部隊多裝備動用情形下的裝備器材需求仿真流程如圖2所示。通過仿真，可以得到一個容量為c的需求仿真樣本S={s1,s2,…,sc}. 對樣本值進行統計，記需求量不超過l的頻率為l，則根據大數定律，當仿真次數足夠多時，可將l視為需求量不超過l的概率Pl的估計。進而可得滿足率為ω時的裝備器材需求量d滿足d-1<ω≤d.

顯然，為了增加估計精度應當增加仿真次數c. 根據中心極限定理，對于給定值ε>0，誤差|l-Pl|<ε的概率為

Pr{|l-Pl|<ε}≈2Φ-1.

(21)

(22)

例如，當ε=0.01、γ=0.95時，通過計算可得c≥9 604.

3 算例分析

下面以某型裝甲裝備的行星變速箱為例，利用本文方法基于MATLAB平臺進行可靠性的統計分析和需求預測。已知廠家給出的變速箱更換周期為500 h. 現得到某旅級單位該裝備動用過程中一組容量為10的變速箱更換時間樣本：229 h、388 h、444 h、468 h、500 h*、500 h*、500 h*、500 h*、500 h*、500 h*(帶*的為截尾數據)。

3.1 參數估計

3.1.1 參數m已知而參數λ未知

首先利用定時截尾的極大似然估計法[8]擬合威布爾分布，得到參數=3.70. 若選擇無信息先驗分布，則由 (7) 式可得參數=600.44. 若選擇共軛先驗分布，則以廠家提供的更換周期作為先驗信息，設變速箱壽命均值和標準差均為500 h. 當λ～Ga(α,β)時，1/λ服從逆伽馬分布IGa(α,β)，通過矩估計法可得

(23)

經求解，α=3，β=2×5003.70，則=548.72.

3.1.2 參數m和參數λ均未知

首先利用區間估計法得到參數m的置信度為95%的取值區間[2.46,6.38]。設置模擬次數為10 000次，m(0)=(2.46+6.38)/2=4.42，η(0)=500，σ=5. 為了消除數據之間的相關性，只取后5 000次模擬結果進行統計分析，如圖3所示。根據平方損失最小原則，以樣本均值作為估計結果，得到=6.19，3=600.63.

圖4所示分別為3種方法得到的變速箱壽命概率密度函數。若以可靠度為0.5時的中位壽命表示變速箱的可靠性，則3種方法得到的變速箱壽命分布模型的中位壽命分別為544 h、496 h和566 h. 由圖4可以看出：形狀參數和尺度參數均未知且采用無信息先驗分布得到的變速箱壽命模型W3的可靠性最高；形狀參數已知且采用無信息先驗分布得到的壽命模型W1的可靠性次之；形狀參數已知且采用共軛先驗得到的壽命模型W2的可靠性最低。

由于無信息先驗分布只根據樣本信息進行參數估計，得到的變速箱的可靠性更高；共軛先驗分布綜合了廠家數據和樣本信息，因此得到的變速箱可靠性較低。由此可見，從樣本數據來看，廠家低估了變速箱的使用可靠性，部隊可適當延長更換周期；另一方面，從廠家數據來看，此樣本僅為一次小樣本抽樣，存在一定的偶然性，需要繼續觀察。

使用K-S法對上述3種方法得到的壽命分布模型進行擬合優度檢驗，檢驗結果如表1所示。3個分布模型的K-S檢驗統計量的值分別為D(1)=0.04，D(2)=0.15，D(3)=0.16. 令檢驗水平α=0.05，通過查表得到臨界值D10,0.05=0.41. 由表1可以看出，各壽命分布檢驗統計量的觀測值均小于臨界值，故均不能拒絕原假設；其中第1個壽命分布模型的檢驗統計量最小。因此，選擇形狀參數m=3.70、尺度參數λ=600.44的威布爾分布作為變速箱的壽命分布模型，并進行需求預測。

3.2 需求預測

得到變速箱的壽命分布后，即可根據裝備動用信息對其需求量進行預測。假定該單位有100輛此型裝備，變速箱的已使用時間和下一年度動用計劃如表2所示。此外，該單位計劃下一年度安排10臺裝備進行預防性維修，預防性維修時對變速箱采取視情維修策略，更換概率為0.2. 下面運用3.1節得到的壽命分布模型計算滿足率為0.95時該單位的變速箱年度需求量。

運行仿真程序執行10 000次仿真，得到的仿真樣本及其經驗分布函數分別如圖5和圖6所示。由圖6可知，8<0.95<9，故滿足率為0.95時變速箱年度需求量為9臺。圖6同時顯示了不考慮變速箱已使用時間時的需求仿真結果，同樣滿足率下變速箱的年度需求量僅為4臺。由此不難發現，已使用時間對變速箱需求預測的影響顯著，不可忽略。

為了進一步研究已使用時間對不同裝備器材需求預測的影響，保持其他數據不變，僅等比例改變裝備器材壽命分布模型的尺度參數η和已使用時間，分析忽略已使用時間造成的需求量絕對百分比誤差的變化規律，結果如圖7所示。從圖7中可以看出：對于低可靠性裝備器材，在計算年度需求時已使用時間可忽略不計；高可靠性裝備器材的年度需求受已使用時間的影響更加顯著，若忽略則會產生較大誤差。

4 結論

本文針對裝備器材需求數據少、需求規律難以掌握的問題，以威布爾型裝備器材為對象，綜合利用貝葉斯法和蒙特卡洛仿真方法對小樣本條件下裝備器材需求預測問題進行了研究。得出以下結論：

1)貝葉斯法能夠綜合廠家的先驗信息和部隊使用過程中產生的樣本信息對裝備器材可靠性進行估計，尤其在預測非穩態需求時非常有效。

2)使用蒙特卡洛仿真方法可以綜合考慮各種因素對威布爾型器材需求進行簡單有效的預測。

3)在對可靠性較高的裝備器材進行需求預測時，忽略已使用時間會產生較大誤差。

4)本文的研究成果可為裝備器材保障部門更好地掌握高新裝備器材消耗規律、提高裝備器材精確化保障水平提供一定的參考。

)

[1] 劉旭陽，吳龍濤，周萬里.基于ARIMA模型的裝備器材需求預測方法[J].裝甲兵工程學院學報，2016，30(6)：21-25.

LIU Xu-yang, WU Long-tao, ZHOU Wan-li. Equipment material demand forecasting method based on ARIMA model [J]. Journal of Academy of Armored Force Engineering, 2016, 30(6): 21-25. (in Chinese)

[2] 劉慎洋，高崎，李志偉，等.基于使用壽命的裝備備件消耗預測模型[J].裝甲兵工程學院學報，2015，29(4)：27-30.

LIU Shen-yang, GAO Qi, LI Zhi-wei, et al. Forecasting model of equipment spare parts consumption based on service lives [J]. Journal of Academy of Armored Force Engineering, 2015, 29(4): 27-30. (in Chinese)

[3] 潘顯俊，張煒，趙田，等.分數階離散灰色模型及其在備件需求預測中的應用[J]. 兵工學報，2017，38(4)：785-792.

PAN Xian-jun, ZHANG Wei, ZHAO Tian, et al. Fractional order discrete grey model and its application in spare parts demand forecasting [J]. Acta Armamentarii, 2017, 38(4): 785-792.(in Chinese)

[4] 徐宗昌，張永強，呼凱凱，等.備件攜行量研究方法綜述[J].航空學報，2016，37(9)：2623-2633.

XU Zong-chang, ZHANG Yong-qiang, HU Kai-kai, et al. Survey on amount configuration methods of carrying spare parts [J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(9): 2623-2633. (in Chinese)

[5] 張瑋，楊善林，劉婷婷.基于威布爾分布的MRO企業估算零部件維修更換率的方法[J].計算機集成制造系統，2013,19(9)：2237-2243.

ZHANG Wei, YANG Shan-lin, LIU Ting-ting. Estimation method of maintenance replacement rate in MRO enterprise based on Weibull distribution [J]. Computer Integrated Manufacturing Systems, 2013, 19(9): 2237-2243. (in Chinese)

[6] 劉天華，張志華，梁勝杰，等.一種Weibull型備件需求量的改進算法[J].系統工程理論與實踐，2012，32(5)：1124-1128.

LIU Tian-hua, ZHANG Zhi-hua, LIANG Sheng-jie, et al. An improved method for the spare demand of the Weibull-distribution [J]. Systems Engineering-Theory & Practice, 2012, 32(5): 1124-1128. (in Chinese)

[7] 劉天華，張志華，李慶民，等.威布爾型多不可修部件備件需求確定方法[J].系統工程理論與實踐，2012，32(9)：2010-2015.

LIU Tian-hua, ZHANG Zhi-hua, LI Qing-min, et al. Determination method of the spare demand for multiple components with Weibull distribution [J]. Systems Engineering-Theory & Practice, 2012, 32(9):2010-2015. (in Chinese)

[8] 劉曉舟，李再興.截尾情形下Weibull分布的最大似然估計[J].數學理論與應用，2014，34(1)：125-128.

LIU Xiao-zhou, LI Zai-xing. Maximum likelihood estimates of the Weibull distribution under censoring samples [J]. Mathematical Theory and Applications, 2014, 34(1):125-128. (in Chinese)

[9] Michael S H, Alyson G W, Reese C S, et al. Bayesian reliability [M]. NY，US: Springer, 2014:67-70.

[10] 董驍雄，陳云翔，項華春，等.基于SST和Bayes的初始備件需求確定方法[J].北京航空航天大學學報，2017，43(5): 1186-1192.

DONG Xiao-xiong, CHEN Yun-xiang, XIANG Hua-chun, et al. Determination method of initial spares requirement based on similarity system theory and Bayesian theory [J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(5): 1186-1192. (in Chinese)

[11] 張春曉，周美茵，王志平.基于貝葉斯更新的飛機結構腐蝕可靠度模型[J].航空學報，2014，35(7)：1931-1940.

ZHANG Chun-xiao, ZHOU Mei-yin, WANG Zhi-ping. Reliability model of aircraft structure corrosion based on Bayesian updated [J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(7):1931-1940. (in Chinese)

[12] 李蕓，易志強.大規模MIMO系統的改進MCMC檢測算法[J].電子技術應用, 2016,42(5)：101-104.

LI Yun, YI Zhi-qiang. An improved MCMC detection algorithm for massive MIMO systems [J]. Application of Electronic Technique, 2016, 42(5):101-104. (in Chinese)

[13] 任偉建，任新元，王磊，等. 基于威布爾分布的油田機采井故障率研究[J].信息與控制，2015，44(6)：722-728.

REN Wei-jian, REN Xin-yuan, WANG Lei, et al. Research on failure rate of pump recovery well in oil fields based on the Weibull distribution [J]. Information and Control, 2015, 44(6): 722-728. (in Chinese)

[14] 張永強，徐宗昌，孫寒冰，等.基于蒙特卡洛仿真和并行粒子群優化算法的攜行備件優化[J].兵工學報，2016，37(1)：122-130.

ZHANG Yong-qiang, XU Zong-chang, SUN Han-bing, et al. Optimization of carried spare parts based on Monte Carlo simulation and parallel particle swarm optimization algorithm [J]. Acta Armamentarii, 2016, 37(1):122-130.(in Chinese)

DemandForecastingofEquipmentandMaterialsbyWeibullDistributionBasedonBayesianEstimationandMonteCarloSimulation

WU Long-tao, WANG Tie-ning, YANG Fan

(Department of Technical Support Engineering，Academy of Army Armored Force，Beijing 100072，China)

The demand of new equipment and materials cannot be mastered well because of less historical demand data and undefined demand. To address this problem, a demand forecasting method based on Weibull distribution is proposed for equipment and materials in the case of small failure samples. The parameters of equipment and material life distribution are estimated by Bayes estimation and MCMC simulation for K-S goodness-of-fit test, including scale and shape parameters. A Monte Carlo simulation-based forecasting method for the annual demand of equipment and materials is presented, in which repairing maintenance, preventive maintenance and service time of equipment and materials are considered. The analysis of examples shows that the life distribution model derived from Bayes estimation has higher degree of fitting in the case of few samples, and the Monte Carlo simulation-based forecasting method is simple and effective.

ordnance science and technology； equipment and material demand forecasting； Weibull distribution； Bayes estimation； Monte Carlo simulation

E92

A

1000-1093(2017)12-2447-08

10.3969/j.issn.1000-1093.2017.12.019

2017-05-18

軍隊科研計劃項目 ( 2016ZB31 )

吳龍濤(1989—)， 男， 博士研究生。 E-mail： tony9076@163.com

王鐵寧(1962—)， 男， 教授， 博士生導師。 E-mail： wtn0728@163.com