基于马氏链蒙特卡洛方法

1基于马氏链蒙特卡洛方法的数控系统可靠性评估李斌全，戴怡（天津职业技术师范大学，机械工程学院，天津300222）摘要：数控系统属于高可靠性产品，应用贝叶斯理论，已经成为其可靠性评估的重要方法。由于数控系统寿命试验数据服从Weibull分布，导致后验分布出现数值大、高维复杂的情况，贝叶斯计算依靠的数值积分方法难以实施。根据马氏链蒙特卡洛（MCMC）方法思想，建立后验平稳分布的马尔科夫链，对分布参数的贝叶斯估计进行求解，解决了数值积分问题，保证可靠性评估的有效实施。通过与BUGS软件结果进行比较，表明提高了模型计算的稳健性、有效性及精度。关键词：数控系统；马氏链蒙特卡洛（MCMC）方法；可靠性评估；Weibull分布；贝叶斯理论中图分类号：TH17ThereliabilityevaluationofNCsystembasedonMarkovchainMonteCarlomethodLiBin-quan，DaiYi(SchoolofMechanicalEngineering,TianjinUniversityofTechnologyandEducation,Tianjin300222，China)Abstract:NCsystemwashigh-reliabilityproduct,sotheapplicationofBayesiantheoryhadbecomeanimportantmeansofitsreliabilityevaluation.ThefailuredataofNCsystemobeyedWeibulldistributionwhichhadcomplexforms.SotheposteriordistributionbecamemorecomplicatedandthenumericalintegrationwhichBayesiancomputingdependedonwasnotavailable.TheproblemofnumericalintegrationwassolvedbyMarkovchainMonteCarlo(MCMC)method.Itensurestheeffectiveimplementationofreliabilityevaluation.ComparedwithBUGSpackage,theresultofBayesianestimationshowsthatitincreasestherobustnessandeffectivenessofthecalculation.Keywords:NCsystem；MarkovchainMonteCarlo(MCMC)method；reliabilityevaluation；Weibulldistribution；Bayesiantheory0引言可靠性工程中，对产品进行可靠性研究或可靠性特征量预测，最基本的工作是要确定其寿命分布参数。据多位学者的工作，数控系统寿命服从两参数Weibull分布[1,2]。其分布参数的贝叶斯估计需要对后验分布进行二重积分。根据贝叶斯公式，后验分布为基金项目：国家自然科学基金（50875186）；国家数控机床科技重大专项（2009ZX04014-013）00),(),(),(),()|,(dmdmmLmmLTm（1）其中，m和为Weibull分布两参数，),(mL为样本似然函数，),(m为先验分布。则分布参数的贝叶斯估计为00)|,()(dmdTmmmE（2）00)|,()(dmdTmE（3）由于数控系统寿命较长且Weibull分2布形式复杂，使得被积函数数值巨大。在实际评估中，计算（2）、（3）式的内层积分时就已经出现浮点数溢出情况，导致依靠被积函数值计算的数值积分方法失效随着马氏链蒙特卡洛（MCMC）方法的不断发展，贝叶斯后验分布计算有望得到彻底解决。通过MCMC方法，复杂的后验分布)|,(Tm被直观的模拟出来，这就使得参数求解运用的数值积分方法转化成从简单的分布中抽样并推断。基于抽样所得分布参数样本，后验分布的统计计算简单易行，各种可靠性特征量求解的可操作性显著提高。目前贝叶斯计算多用国外的BUGS软件进行[3,4]，但是该软件未公布其核心算法，这给掌握算法实质以及进行算法改进带来困难。另一方面BUGS软件采用的Gibbs抽样是MCMC方法的特例，所以自主掌握MCMC方法，并应用于可靠性评估有重要意义。本文利用MCMC方法对参数的贝叶斯估计进行求解，并与BUGS软件结果进行比较，表明其适用于数控系统的可靠性评估。1基于MCMC方法的贝叶斯可靠性1.1截尾样本的Weibull似然函数构建对数控系统进行有替换定时截尾寿命试验，得到),,,,,(2121srLLLtttT，其中，rttt,,,21为完全失效数据样本，sLLL,,,21为截尾数据样本。两参数Weibull分布概率密度函数)(tf为])/(exp[)/(1mmmtmt，生存函数)(St为])/(exp[mt，其中m为形状参数，为尺度参数。则截尾样本似然函数为sjmjrimirimirmrsjjriiLttmLtfmL111111)(exp)(exp)(S)(),(（4）1.2后验分布形式的确立取参数m的先验分布为伽马分布),(~)(baGam，其中a为形状参数，b为尺度参数。参数的先验分布为逆伽马分布),(~)(vuIGa，其中u为形状参数，v为尺度参数。则两参数的联合先验分布为/exp/1)()exp()(),(),(),(11vuvbmmabvuIGabaGamuuaaΓΓ（5）由贝叶斯公式得到后验分布：00111111(,)(,)(,|)(,)(,)exp()()/rsjmmiijrarmirmuiLmmmTLmmdmdLtbmvmt（6）由于后验分布)|,(Tm数值巨大，高维复杂，贝叶斯计算依靠的数值积分方法难以实施。但是可以利用MCMC方法思想解决上述问题。2后验平稳分布的马尔科夫链建立贝叶斯后验量的计算可归结为关于后验分布的积分[5]。MCMC方法思想是通过建立后验平稳分布)|,(Tm的马尔科夫链[6]，并从中抽样，来获取后验分布样本进行统计推断，这就不需要对被积函数进行数值计算，从而保证求解过程的有效实施和稳定性。令(,)(,)(,)pxyqxyxy，表示从易于抽样的建议分布(|)qx产生一个潜在的转移xy，然后根据概率(,)xy决定是否转移。如果链在时刻t处于状态x，即()tXx，从U(0,1)nif中抽一个随机数u，3(1),(,),(,)tyuxyXxuxy（7）则此时，(,)pxy为马氏链一个转移核。在选取),(q后，根据（6）式，概率(,)xy取为(,|)(,)min1,(,|)(,)yyxxmTqyxmTqxy（8）这时，转移核(,)pxy以)|,(Tm为马氏链平稳分布。选取均匀分布分布为建议分布，此时，11111111(,)min{1,exp[()()/()()/](/)}/yyxxyxyxrmiiysrmjmiyyjiyxsjmxxjxrarmmyxirmurmuiyxtxyLtbmvLbmvmmt（9）其中，(,)xy为转移概率，ym和y是在一步迭代中，从建议分布中抽取的建议值，xm和x是一步迭代的初始值。过程如下:1,...,:tN（1）(1)(1)01(,)(,)ttm（2）从建议分布U(0.5,0.5)nif中产生候选点（3）令，即让建议值随机游动一个距离（4）根据式（9）计算接受概率(,)（5）以概率(,)接受()t，否则令()t通过上述过程，后验分布)|,(Tm的马尔科夫链被建立起来，并得到Weibull分布参数m和的样本，基于所得样本就可以对分布参数进行各种统计推断。3仿真分析通过有替换定时截尾寿命试验，获得数控系统寿命的现场样本如表1所示，表1数控系统寿命现场样本（单位：h）完全失效数据176256339455495117612601896236221342133322431853158151120134102438013604251截尾数据11369315232231245625883322102242131021250246614503210310250116011145141522332根据表1数据，用上述方法编程进行计算，得到结果如表2表2MCMC方法计算结果建议分布接受率E(m)E(η)sd(m)sd(η)均匀分布0.511.3540990.171.8BUGS是进行贝叶斯推断的软件包，用该软件处理表1数据，得到结果如表3表3应用BUGS软件计算结果nodemeansdstartsamplem1.320.18100157000η41114.1100157000通过表2计算结果，表明MCMC方法的有效性与稳定性，与表3进行比较，可以看出其计算精度也有所提高。图1、2分别是m和η的全部迭代轨迹和后验分布的密度直方图。0200005000030005000Iterations0200005000030005000IterationsValuesofMTBF020000500000.81.42.0Iterationsm（a）参数m40200005000030005000Iterations0200005000030005000IterationsValuesofMTBF020000500000.81.42.0Iterationsm（b）参数η图1参数m和η的迭代轨迹图mDensity1.01.52.00.01.02.0Density3000500070000e+006e-04（a）参数mmDensity1.01.52.00.01.02.0Density3000500070000e+006e-04（b）参数η图2后验分布的密度直方图4数控系统可靠性特征量的计算用建议分布为均匀分布时建立的马氏链对数控系统可靠性进行估计，Weibull分布概率密度函数0.351.351.35()1.35/4099exp[(/4099)]fttt（）1(1)3758.7MTBFm采用传统方法计算MTBF值的区间估计是复杂而且难以实施的。利用MCMC方法可以有效地予以估计。图3是MTBF值的迭代轨迹。mDensity1.01.52.00.01.02.0Density3000500070000e+006e-040200005000030005000IterationsMTBF图3MTBF值的迭代轨迹图可以得到，MTBF值的90%置信区间为[3666.76，3989.30]。5结论（1）应用贝叶斯理论进行可靠性评估是一种有前途的方法。由于以往不能很好地解决高维数值积分问题，从而制约了该理论在实际中的应用。MCMC方法解决了高维数值积分问题带来的不便，有利于贝叶斯理论在可靠性评估中的推广。（2）由于后验分布)|,(Tm被模拟构造出来，后验分布的各种统计推断都可从简单分布中抽样进行。基于抽样所得分布参数样本，可靠度、失效率等各种可靠性特征量的求解都可以有效实施，不再受数值积分的限制，提高了模型计算的稳定性、可操作性与适应性[参考文献][1]戴怡，周云飞，陈学东等.加工中心故障分布规律及其研究方法[J].系统工程与电子技术，2004，26(3)：413-415.[2]张海波，贾亚洲，周广文.数控系统故障间隔时间分布模型的研究[J].哈尔滨工业大学学报