吴龙涛， 王铁宁， 杨帆

(陆军装甲兵学院 技术保障工程系， 北京 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