王鐵寧， 吳龍濤， 楊 帆

(裝甲兵工程學(xué)院技術(shù)保障工程系， 北京 100072)

基于貝葉斯法的定時截尾小樣本指數(shù)型裝備器材需求預(yù)測

王鐵寧， 吳龍濤， 楊 帆

(裝甲兵工程學(xué)院技術(shù)保障工程系， 北京100072)

針對高新裝備器材故障數(shù)據(jù)少、需求規(guī)律不明確的問題，提出了一種定時截尾小樣本條件下指數(shù)型裝備器材的需求預(yù)測方法?；谪惾~斯法進(jìn)行了裝備器材壽命分布參數(shù)估計，討論了先驗分布和損失函數(shù)的選擇問題。引入了K-S檢驗法對壽命分布模型進(jìn)行擬合優(yōu)度檢驗，并設(shè)計了壽命分布參數(shù)估計的卡方檢驗方法。綜合考慮故障更換和定時更換，提出了部隊裝備器材年度需求預(yù)測方法，并采用蒙特卡洛法進(jìn)行了仿真驗證。結(jié)果表明:基于貝葉斯法的參數(shù)估計結(jié)果能夠順利通過檢驗，需求預(yù)測結(jié)果與仿真結(jié)果一致。

定時截尾小樣本; 裝備器材； 指數(shù)分布； 需求預(yù)測； 貝葉斯法

隨著部隊高新裝備列裝步伐的加快和換件修理方式的廣泛開展，裝備器材需求結(jié)構(gòu)逐漸發(fā)生變化，如對光電類及總成類裝備器材的需求量顯著增加。由于高新裝備列裝時間短，當(dāng)前高新裝備器材保障普遍存在消耗規(guī)律掌握不清、保障不及時的難題。

目前，我軍主要采用定額計算[1]方法管理庫存，各級倉庫一般依據(jù)庫存標(biāo)準(zhǔn)，并結(jié)合現(xiàn)有庫存來制定需求計劃。但由于庫存標(biāo)準(zhǔn)本身就是粗略、不精確的估計[2]，因此定額計算方法難以精確測算裝備器材需求。同時，由于高新裝備服役時間短，可利用的歷史數(shù)據(jù)較少，故難以應(yīng)用傳統(tǒng)的基于大樣本的方法，如時間序列分析、回歸分析和神經(jīng)網(wǎng)絡(luò)方法等進(jìn)行器材需求預(yù)測。灰色預(yù)測模型[3]是小樣本預(yù)測經(jīng)常采用的方法，但其預(yù)測精度較低，且無法計算裝備器材滿足率。

基于可靠性的需求預(yù)測方法根據(jù)失效機理選擇裝備器材壽命分布類型[4]，并以此建立裝備器材需求計算模型。然而，由于裝備承制單位在交付裝備時往往無法提供其壽命分布的具體形式，且不同部隊的裝備使用和維護(hù)水平以及使用環(huán)境差異明顯，裝備承制單位提供的數(shù)據(jù)只能作為參考[5]。因此，筆者針對定時截尾小樣本，提出了基于貝葉斯法[6]的指數(shù)型裝備器材的需求預(yù)測方法，利用貝葉斯法估計裝備器材的壽命分布參數(shù)，然后綜合考慮故障更換和定時更換，給出了裝備器材的需求預(yù)測方法，并進(jìn)行了仿真驗證。

1 指數(shù)型壽命分布參數(shù)估計

1.1指數(shù)分布的一般形式

指數(shù)分布的重要特性之一是故障率λ(t)=λ為常數(shù)，即產(chǎn)品在下一時刻故障失效的概率與已使用時間無關(guān)，這也被稱作指數(shù)分布的“無后效性”[7]。設(shè)某裝備器材壽命T服從指數(shù)分布，其概率密度函數(shù)、累積分布函數(shù)和產(chǎn)品的可靠度函數(shù)分別為

f(t)=λexp(-λt),(t,λ>0)，

(1)

F(t)=1-exp(-λt),

(2)

R(t)=exp(-λt)。

(3)

E(T)=1/λ，為指數(shù)分布的期望，即產(chǎn)品的期望壽命(Mean Time To Failure，MTTF)。由指數(shù)分布的一般形式可知：對指數(shù)型裝備器材進(jìn)行可靠性估計就是估計故障率λ。極大似然估計法是一種在大樣本條件下簡單有效的參數(shù)估計方法，但在截尾小樣本情形下偏差較大，且不易收斂[8]。貝葉斯法可充分利用專家經(jīng)驗、裝備承制單位的參考數(shù)據(jù)等先驗信息來提高參數(shù)估計的準(zhǔn)確性，同時隨著新信息的出現(xiàn)可不斷更新結(jié)果，是解決小樣本數(shù)據(jù)預(yù)測的有效方法[9]。

1.2故障率λ先驗分布的選擇

應(yīng)用貝葉斯法進(jìn)行壽命分布參數(shù)估計時，首先要為參數(shù)選擇合適的先驗分布。通?？蛇x擇無信息先驗分布和有信息的共軛先驗分布。

1)當(dāng)關(guān)于待估參數(shù)的先驗信息很少時，一般選擇無信息先驗分布。根據(jù)最常用的確定無信息先驗分布的Jeffreys法則[9]，λ的先驗分布為

π(λ)=λ-1。

(4)

2)共軛先驗分布是指與后驗分布屬于同一類分布的先驗分布。伽馬分布為指數(shù)分布故障率λ常用的共軛先驗分布，即λ～Γ(α,β)，其概率密度函數(shù)為

(5)

式中：α,β>0，為先驗分布的超參數(shù)，其值通常根據(jù)先驗信息使用矩估計法來確定；Γ(α)為伽馬函數(shù)。

1.3故障率λ的貝葉斯估計

在裝備維修過程中，裝備器材的更換原因大致分為2類：因故障失效被更換和到了指定的維修間隔期進(jìn)行等級維修時被更換。因此，可將裝備器材的更換時間樣本視為定時截尾樣本，截尾時間為定時維修間隔期。

設(shè)裝備的定時維修間隔期為t*，裝備器材更換時間樣本t的容量為n，故障失效時間數(shù)量為z，第i(1≤i≤z)個故障失效時間為ti，則樣本的聯(lián)合分布密度函數(shù)為

(6)

1)若選擇無信息先驗，則樣本t和λ的聯(lián)合概率分布為

(7)

(8)

根據(jù)貝葉斯定理，在給定樣本t下，λ的后驗分布密度函數(shù)為

(9)

即

(10)

2)若選擇共軛伽馬先驗，λ的后驗分布密度函數(shù)為

(11)

(12)

(13)

2 壽命分布參數(shù)估計檢驗

2.1壽命分布模型的擬合優(yōu)度檢驗

通過貝葉斯估計確定了裝備器材壽命分布的形式后，就需要通過壽命分布的擬合優(yōu)度檢驗來驗證分布模型是否適用于樣本。χ2檢驗是一種常用的擬合優(yōu)度檢驗方法，但只適用于大樣本(n≥50)完全數(shù)據(jù)的情形，針對小樣本定時截尾數(shù)據(jù)，筆者引入了K-S檢驗法[5]。

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

H0:Ψ(t)=F(t)；H1:Ψ(t)≠F(t)。

K-S檢驗的統(tǒng)計量為

(14)

式中：Ψz(t)為壽命T關(guān)于給定失效樣本的經(jīng)驗分布函數(shù)。由于Ψz(t)是小樣本，其計算公式[11]為

(15)

2.2參數(shù)估計的假設(shè)檢驗

若默認(rèn)裝備器材壽命T服從指數(shù)分布，則只需要對λ的估計結(jié)果進(jìn)行假設(shè)檢驗。且檢驗假設(shè)為

3 指數(shù)型裝備器材需求預(yù)測

圖1為部隊裝備器材需求發(fā)生的一般過程。

圖1 部隊裝備器材需求發(fā)生一般過程

假設(shè)某部隊有A型現(xiàn)役裝備m臺，每臺裝備有1件器材a。根據(jù)滾動式循環(huán)動用原則，一般每臺裝備的年度動用時間不同，設(shè)Ti(1≤i≤m)為第i臺裝備的年度動用時間。為了便于分析計算，進(jìn)行如下假設(shè)：

1)裝備器材故障失效時采取換件維修，換件時間忽略不計；

2)裝備的動用時間即為裝備器材的動用時間；

3)各裝備中器材a的壽命相互獨立且服從同一指數(shù)分布。

在計算裝備器材年度需求時，需要分別考慮故障更換需求和定時更換需求。已知壽命服從指數(shù)分布的器材a在年度動用時間TS內(nèi)的故障次數(shù)N(Ti)服從泊松分布P(λTi)，即

(16)

式中：k∈N。根據(jù)泊松分布的可加性可得m件器材a的年度故障次數(shù)為N(TS)=N(T1)+N(T2)+…+N(Tm)，服從泊松分布N(TS)～P(λT1+λT2+…+λTm)，即

(17)

若規(guī)定器材a的年度裝備器材滿足率為ρ，則因故障更換產(chǎn)生的裝備器材需求量S1滿足

(18)

根據(jù)裝備維修計劃，一般每年都會安排摩托小時已到達(dá)定時維修間隔期的裝備進(jìn)行等級維修。在等級維修時，器材a可能采取不更換、必須更換或視情更換3種維修策略[2]，對于視情更換，通常指定一個更換概率p。假設(shè)在m臺裝備中有m0臺需要進(jìn)行等級維修，如器材a不更換，其定時更換需求量S2=0；如必須更換，S2=m0；如視情更換，則S2服從二項分布S2～B(m0,p)，同式(18)，滿足率為ρ時的S2為

(19)

裝備器材a的年度需求量為S=S1+S2。

4 算例分析

以某型裝甲車輛的“火控計算機通信板”(簡稱“通信板”)為例，利用本文提出的方法進(jìn)行可靠性統(tǒng)計分析和需求預(yù)測。已知單臺裝甲車輛安裝一塊通信板，通過裝備承制單位獲知通信板的工作壽命“大約為1 200 h”。某部隊在裝備動用過程中得到一組更換時間樣本為81 h，230 h，274 h，387 h，600 h，688 h，798 h，974 h，1 000*h，1 000*h(帶*的為截尾數(shù)據(jù))。

4.1參數(shù)估計

首先，計算通信板故障率λ的后驗分布。若選擇無信息Jeffreys先驗，由式(10)可得λ的后驗分布為Γ(8,6 032)。若選擇伽馬先驗，以裝備承制單位提供的數(shù)據(jù)作為先驗信息，設(shè)通信板壽命均值為1 200 h，標(biāo)準(zhǔn)差也設(shè)為1 200 h，當(dāng)λ～Γ(α,β)時，MTTF服從IΓ(α,β)分布。應(yīng)用矩估計法可得

(20)

由式(20)可得先驗分布的超參數(shù)α=3，β=2 400。再由式(11)可得λ的后驗分布為Γ(11,8 432)。圖2 為λ先驗分布和后驗分布的概率密度函數(shù)對比分析?？梢钥闯觯?種后驗分布得到的故障率都高于裝備承制單位給出的故障率；同時，由Jeffreys先驗得到的后驗故障率高于由伽馬先驗得到的后驗故障率，這是因為Jeffreys先驗為無信息先驗，其后驗分布取決于樣本。一方面，由于裝備承制單位無法預(yù)知部隊的裝備使用環(huán)境及裝備使用、維修人員的水平，因而高估了通信板的實際使用壽命；另一方面，該樣本僅為1次小樣本抽樣，存在一定的偶然性和不確定性。由伽馬先驗得到的后驗分布綜合考慮了裝備承制單位先驗信息和樣本的新信息，因此更加合理。

圖2 λ先驗分布和后驗分布概率密度對比

表1 K-S模型擬合優(yōu)度檢驗結(jié)果

4.2需求預(yù)測

假定該單位有該型裝甲車輛100臺，年度動用計劃為(40,15)、(30,25)、(20,60)、(10,150)，并計劃安排10臺車輛進(jìn)行等級維修，等級維修時通信板采取視情維修策略，更換概率為0.2，需預(yù)測滿足率為0.95時通信板的需求量。

為了驗證預(yù)測結(jié)果的準(zhǔn)確性，采用蒙特卡洛法對裝甲車輛通信板需求預(yù)測過程進(jìn)行仿真。圖3為基于MATLAB平臺的裝備器材需求預(yù)測仿真流程。

圖3 裝備器材需求預(yù)測仿真流程

圖4 10 000次仿真的年度裝甲車輛通信板需求預(yù)測樣本

圖5 仿真樣本的經(jīng)驗分布函數(shù)

5 結(jié)論

筆者以指數(shù)壽命型高新裝備器材為對象，結(jié)合實際情況，對小樣本條件下裝備器材的需求預(yù)測問題進(jìn)行了研究。使用貝葉斯法綜合了廠家先驗信息和使用過程中產(chǎn)生的樣本信息，對裝備器材的壽命分布參數(shù)進(jìn)行了估計，并給出了分布模型和分布參數(shù)的檢驗方法。在此基礎(chǔ)上，綜合考慮故障更換和定時更換，提出了裝備器材年度需求預(yù)測方法，并采用蒙特卡洛法進(jìn)行了仿真驗證，為解決當(dāng)前高新裝備器材消耗規(guī)律掌握不清、保障不及時的問題提供了一種可行的方法。

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

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

[3] 郭曉君，劉思峰，吳利豐. 污染物減排預(yù)測的灰色Markov組合模型與算法[J]. 計算機應(yīng)用研究，2013，20(12)：3670-3673.

[4] 李源，楊建軍.導(dǎo)彈武器裝備備件需求量分析與計算方法[J].兵工自動化，2007，26(8)：22-23.

[5] 蔡景，左洪福，王華偉.基于成本的民用航空發(fā)動機維修方案優(yōu)化研究[J].機械科學(xué)與技術(shù)研究，2007，26(2)：167-171.

[6] 曾志國.貝葉斯可靠性[M].北京：國防工業(yè)出版社，2014：67-70.

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

[8] 馮自立，陳曉陽，顧家銘.Weibull分布下基于矩法的小樣本量定時截尾數(shù)據(jù)的估計[J].四川兵工學(xué)報，2010，31(3)：93-97.

[9] SARHAN A.M. Empirical bayes estimates in exponential reliability model [J]. Applied mathematics and computation 2003，135(1)：319-332.

[10] 尹慧琳，楊筱菡，陸恒.Wiener過程性能退化產(chǎn)品可靠性評估新Bayes方法[J].同濟大學(xué)學(xué)報(自然科學(xué)版)，2015,43(8)：1234-1237.

[11] EBELING C.E.An introduction to reliability and maintainability engineering [M].New York：Waveland Pr Inc,2009:52-53.

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

(責(zé)任編輯: 王生鳳)

DemandForecastingforEquipmentMaterialswithExponentialLifeDistributionBasedonBayesianEstimationunderTypeICensoredSmallSample

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

(Department of Technical Support Engineering, Academy of Armored Force Engineering, Beijing100072, China)

In order to address the problem that the demand discipline of new and high-tech equipment materials cannot be mastered well because of lack of failure data, a method of demand forecasting is proposed for equipment materials with exponential distribution under small failure samples. In a view of small samples, the parameter of equipment material life distribution is estimated by Bayesian estimation, the option of prior distribution and loss function is discussed. Then, statistical tests for goodness-of-fit of the life distribution and parameter estimation are performed by K-S test and Chi-Square test respectively. In consideration of fault replacement and timely replacement, the annual demand forecasting method of army equipment is put forward and devised with a Monte-Carlo simulation test. As the example shows, the estimation result performs well in the test, and the forecasting value is identical to that from the simulation.

type I censored small sample; equipment materials; exponential distribution; demand forecasting; Bayesian estimation

1672-1497(2017)04-0029-06

2017-03-28

軍隊科研計劃項目

王鐵寧(1962-)，男，教授，博士。

E92;N945.24

：ADOI：10.3969/j.issn.1672-1497.2017.04.006