蒙特卡罗方法在维修过程仿真中的应用研究

赵思桥，刘井泉

（清华大学工程物理系，北京100084）

鉴于核电站对安全的特殊要求，对核电站设备进行预防性维修已经成为必不可少的一个环节。目前在有关核电站维修策略制定过程中，主要应用解析计算的方法，然而鉴于电站设备寿命过程的随机性，解析方法已经难以满足电站的经济性、安全性对维修计划的要求。因此采用仿真方法求解维修策略已经逐渐成为研究的热点。目前并没有针对于存在维修条件下设备寿命过程的仿真方法相关的研究。本文对各类维修条件下设备的寿命过程模型进行了分析，对完全维修与最小维修的仿真方法进行了总结并阐述了仿真方法的特点，本文中同样对不完全维修条件下设备的寿命过程进行了建模，并对其仿真方法进行了构造。

根据传统分类可将维修过程分为三种类型：完全维修、最小维修和不完全维修。完全维修是指设备维修后恢复最初时刻的状态，完全维修条件下设备逐次的故障间隔时间分布相互独立同分布；最小维修是指设备在维修后仅恢复到维修前一时刻的状态，设备的寿命过程服从非齐次泊松过程；不完全维修是介于完全维修与最小维修之间的一种维修，也就是设备维修后恢复到“修复如新”与“修复如旧”之间的某一个状态。

对于维修过程的仿真，其主要任务在于对各种维修条件下设备的逐次失效间隔时间进行抽样。对于完全维修，鉴于每次的故障间隔时间分布相互独立同分布，只需重复对相同的寿命分布进行抽样即可。对于最小维修，其抽样就是对确定强度函数的非齐次泊松过程进行抽样，目前业界对非齐次泊松过程的抽样方法有一定的研究，其主要思路是根据逐次失效间隔时间的联合密度分布函数或非齐次泊松过程的特性构造相应的抽样方法，关于最小维修的仿真在第2节会有详细的说明。不完全维修不属于泊松过程，目前有多种描述不完全维修的模型，在核电行业中应用最多的为等效役龄回退因子模型所描述的不完全维修，本文中所提出的有关不完全维修仿真方法是在最小维修仿真方法的基础上进行改进而构造的，有关这方面的内容在第3节会有详细的阐述。

1 最小维修过程仿真

1．1 最小维修

最小维修的定义为：维修使得系统的故障率恢复到失效前的状态。在维修决策优化中，为简化模型通常考虑监测维修为最小维修。在核电站的维修中最小维修的维修时间相对一般较短，在不考虑维修时间的情况下，最小维修条件下的寿命过程是一种非齐次泊松过程，其存在每次故障间隔时间相互不独立的特点。

1．2 非齐次泊松过程

非齐次泊松过程是用来描述可修复系统的相邻故障间隔呈某种趋势的寿命过程，是齐次泊松过程的一种推广。

非齐次泊松过程的数学的定义为：

｛N（t）；t≥0｝若满足以下条件称作是一个非齐次泊松过程。

（1）N（0）＝0；

（2）N（t）具有独立增量；

（3）对任意的0≤t＜s，在（t，s］中的故障数N（s）－N（t）有参数的泊松分布：

非齐次泊松过程与齐次泊松过程的区别在于其相邻故障间隔时间既不独立亦不同分布，在某一区间上其故障参数还依赖于区间的起点。

泊松过程的强度函数决定了泊松过程，电厂的设备寿命分布多为威布尔分布，因此当累计强度函数呈威布尔分布形式时将作为主要的研究对象，这类非齐次泊松过程称为威布尔过程。

1．3 非齐次泊松过程仿真方法

对于最小维修过程的仿真主要即为对服从某一非齐次泊松过程的设备失效间隔时间进行抽样。目前业界存在几种关于非齐次泊松过程的仿真方法，但没有有关各种方法的适用范围和使用特点的说明，事实上其中有很多种是难以在实际应用中得以实施的。根据对目前主要应用的仿真方法进行比较，本人筛选出适合于核电设备的两种：间隔时间产生方法、尺度变换方法。

方法一：间隔时间产生方法

间隔时间产生方法是根据设备的失效密度函数直接抽样从而得到抽样样本，根据非齐次泊松过程的特点，设备第n次失效间隔时间的分布函数是与第n－1次失效的发生时间Sn－1相关，其关系式如定理1所述。

定理1［1］：设S0＝s0＝0，在Sn－1＝sn－1的条件下，tn＝Sn－Sn－1的条件密度函数为：

式中：Sn代表发生第n次失效的时间点，tn代表第n次失效发生的失效间隔时间，λ为非齐次泊松过程的强度函数。

根据定理1的叙述，在Sn－1＝sn－1条件下第n次的失效间隔时间的条件分布函数为：

采用间隔时时间产生方法其最大的限制在于需要对条件分布函数进行抽样，无论采用何种蒙卡的抽样方法都需要分布形式不能过于复杂，因此间隔时间产生方法更多适用于失效率形式相对简单的情况，例如威布尔分布（Λ（t）＝（t／α）β）等。

在设备寿命服从威布尔过程的条件下，其条件寿命分布函数的形式为：

式中：α是威布尔分布中的尺度参数，β为形状参数。

该方法是通过依次产生设备失效间隔时间来对设备的寿命过程进行仿真，在仿真过程中需要给出一个截止条件，可以选择为观测总时间。

方法二：尺度变换方法

尺度变换方法主要是根据非齐次泊松过程与齐次泊松过程之间的关系而构造出的一种抽样方法，其主要依据如定理2所述：

定理2［1］：设T1，…，Tn为强度函数为λ（t）的非齐次泊松过程头n个事件发生时刻，记，则Z1，…，Zn独立同分布，服从参数为1的指数分布。

根据定理2所构造的尺度变换方法的仿真步骤为：

1．产生λ＝1的齐次泊松过程前n个事件的发生时刻Z1，…，Zn；

2．计算Si＝Λ－1（Zi），则S1，…，Sn满足强度为λ（t）的非齐次泊松过程。

根据上面的步骤可以看出尺度变换方法是利用已经产生的齐次泊松过程的失效时间通过一定的时间尺度变换映射，得出非齐次泊松过程的失效时间。在这一过程之中最主要的计算步骤在于计算累计强度函数的逆函数，这一点就可以降低该方法相对于时间间隔产生方法对失效率形式复杂度的要求。对于形式较为复杂的情况，例如常用于描述电器类设备失效的对数正态分布，很难构造出适用于它的时间间隔产生方法的抽样方法，而采用尺度变换方法，可以通过数值积分的方法计算出变量t与累计强度函数Λ（t）之间的映射关系，在计算Λ－1的时候，只需从之前得到的映射中寻找对应值即可。

综上，尺度变换方法主要适用于失效率分布形势复杂的情况。在仿真的过程中需要给出一定的截止条件，由于在做尺度变换的前提是已经仿真得出所有的Z1，…，Zn，因此截止条件可设置为总的失效次数n。

1．4 实例计算

为了验证上面叙述的两种仿真方法是否正确，对一假设参数的算例进行抽样，并对结果进行拟合优度检验。为方便起见选择失效率为威布尔分布：λ（t）＝2×10－8t。适用于威布尔分布的检验方法可以采用Cramer Von－Mises检验。

1．4．1 Cramer Von－Mises拟合优度检验

Cramer Von－Mises检验方法［2］为：记n次故障的故障时间分别为T1，T2，…，Tn；

对定数截尾情况，统计量C2M计算表达式：

若统计量计算值大于某一置信度下的统计量临界值，则说明在该置信度下接受样本服从该过程，反之则舍弃该样本。

1．4．2 实例计算结果

分别采用时间间隔产生方法以及尺度变换方法对设备最小维修条件下的5次寿命过程进行了仿真计算，得到每一次寿命过程中的设备逐次失效时间点Ti的结果如表1所示。

表1 最小维修抽样结果Table 1 Sampled data of least maintenance

寿命过程1寿命过程2寿命过程3寿命过程4寿命过程5 9 947．92 7 671．81 16 130．96 12 770．2 9 597．52 T218 193．34 9 815．34 19 353．89 13 047．81 9 714．52 T322 337．03 13 172．31 19 371．01 14 229．92 19 895．92 T418 722．35 21 422．28 14 990．78 20 079．97 T525 777．95 23 853．56 19 294．04 20 080．48 T628 050．23 29 158．68 20 633．39 T721 010．94 T821 070．96 T924 942．03 T1025 388．5 T1 T11 25 409．91

采用Cramer Von－Mises检验：

通过拟合优度检验可以在较大的置信度下接受抽样样本服从威布尔分布，因此两种方法得到的结果是可信的，由此证明了两种抽样方法的可行性。

2 不完全维修过程仿真

不完全维修是最重要也是最难研究的一类维修过程，在电站的管理中，考虑设备的预防性维修为不完全维修。不完全维修的难点在于设备在该条件下维修后状态的无法确定，设备的寿命过程无法继续采用非齐次泊松过程进行分析。因此在解决不完全维修过程仿真方法的过程中，建立设备寿命过程模型是其中一个关键的问题。

在目前研究中，Pham和Wang将描述不完全维修的各种方法归纳为七种类型［3］：（p，q）法则、（p（t），q（t））法则、复合（p，q）法则、（α，β）法则、冲击模型、等效役龄、回退因子。在诸多描述方法中等效役龄与回退因子具有较强的灵活性和实用性，在核电站的数据分析中最多考虑。本章中将主要根据等效役龄与回退因子建立寿命过程模型并提出仿真方法。

2．1 等效役龄与回退因子

在Chan和Shaw［4］的研究中认为预防性维修过程改善了设备的可靠性，设备可靠性的改善体现在失效率与可靠度的减少，其减少的程度用回退因子η来描述，而回退因子通常是由设备的役龄与维修的次数等参数决定的。根据之前的研究结果，η可以写成预防性维修成本与维修次数的函数：

式中：a，b是回退因子的计算参数；Cpmi是第i次预防性维修的总成本，它的具体数值是与预防性维修固定成本Cpmf、预防性维修变动成本Cpmv以及设备役龄Ti均存在关系，计算关系式可以写为：

因此回退因子的表达式可以写为：

2．2 不完全维修过程仿真方法

目前还没有任何有关不完全维修仿真方法的研究，通过分析引入回退因子后设备寿命过程与非齐次泊松过程之间的关系，在以非齐次泊松过程的仿真方法为基础之上提出了不完全维修过程的仿真方法，填补了这方面的空白，为应用仿真方法求解维修策略提供了支持。

本文中提出的不完全维修过程仿真方法是在最小维修的两种仿真方法：间隔时间产生方法与尺度变换方法基础之上进行针对于回退因子的修改构造得出的。

方法一：间隔时间产生方法

若设备第n－1次维修后认为达到等效役龄Yn－1，则设备的第n次失效间隔时间分布函数可以写为：

因此在考虑设备的等效役龄后，其失效间隔时间分布密度函数只需要将非齐次泊松过程的寿命分布函数中的第n－1次失效发生时间sn－1替换为第n次失效发生时设备的等效役龄Yn－1。

同样，这种方法只适用于失效率分布形式相对简单的情况，例如对于失效率服从威布尔分布情况，条件密度函数服从：

方法二：时间尺度变换方法

根据定理3的描述，对于服从非齐次泊松过程的前n个失效时刻S1，…，Sn，对其进行尺度变换后的Z1，…，Zn服从强度为1的齐次泊松过程。在存在回退因子的情况下，显然不能简单地将Zn进行尺度变换的逆变换，这是由于第n－1次维修后设备的等效役龄不同于维修前设备的等效役龄。因此需要对n－1次维修前的役龄进行一定的变换使其与维修后的役龄相协调才可以对Z1，…，Zn进行时间尺度的逆变换。

这种方法的本质在于构造出符合于存在回退因子情况下经过时间尺度变换后服从于齐次泊松过程的逐次失效发生时间点Z＊1，…，Z＊n，再对其进行尺度逆变换从而得到需要的结果。

其具体的步骤为：

1）对Z1进行尺度变换产生第一次维修前的设备寿命；

2）设备寿命与回退因子η相乘产生第一次维修后的等效役龄；

5）以S2作为步骤二的输入，重复上述步骤，直到产生所有的样本S1，…，Sn。

其流程如图1所示。

图1 尺度变换方法流程图Fig．1 Flowchart of time scale conversion method

实例计算：

为了验证本文提出的两种不完全维修仿真方法是否正确，需要对实例进行计算，并对结果进行验证。需要说明的是目前并没有针对役龄回退因子模型适合的检验方法，为此在实例计算中采用与最小维修一节相同的算例，并对两者之间的结果进行比较，从而给出分析与结论。采用的抽样分布为威布尔分布：λ（t）＝2×10－8t；役龄回退因子与设备役龄之间的关系设定为η＝（0．1＋1×10－5T）2，对该条件下设备寿命过程进行仿真，逐次的失效时点抽样样本如表2所示。

表2 不完全维修抽样结果Table 2 Sampled data of incomplete maintenance

将上表中的结果与表1中的结果进行比较，由于采用的随机数完全一致，因此两种维修模型设备出现故障的时间应处在相同的分布密度函数区间内。通过比较可以看出对于每一次寿命过程除首次故障时间外设备的寿命均有所增加，出现这种情况恰是由于引入了役龄回退因子，改善了维修后的设备可靠性，增加了设备的维修间隔时间。

3 结论

在核电站的维修策略优化过程中，预防性维修是计算的主要对象。由于设备的失效存在概率性的特点，因此选用蒙特卡罗方法对存在维修的设备寿命过程进行模拟所得到的结果相对数值方法而言存在应用简便、易于扩展的特点。

在本文中分别针对最小维修与不完全维修提出了设备寿命过程的抽样方法。最小维修模型即非齐次泊松过程，在该模型下本文总结了两种设备寿命过程的抽样方法：间隔时间产生法、尺度变换方法。通过对两种方法的使用特点进行分析提出间隔时间产生方法局限于失效率形式相对简单的情况，尺度变换方法步骤相对复杂但可以处理失效率复杂的情况。通过实例计算，对设备的寿命过程进行仿真并对设备逐次的失效时间点进行拟合优度检验，可以证明两种抽样方法的正确性。不完全维修是一种复杂的过程，它不属于非齐次泊松过程但与其存在一定的联系。在通过对最小维修模型的两种抽样方法进行修改后本文构造出适用于等效役龄回退因子描述的不完全维修模型下设备寿命过程的仿真方法。通过分析实例计算得到的数据可以证明两种方法的正确性不完全维修过程仿真方法的提出完善了核电站设备寿命过程的仿真方法，为应用蒙卡仿真方法计算维修策略提供了支持。

［1］ 曹晋华，程侃．可靠性数学引论［M］．北京：高等教育出版社，2006．

［2］ Darling D A．The Kolmogorov－Smirnov，Cramer－von Mises Tests，Ann［J］．Math．Statist，1957，284，28（4）：823．

［3］ Chan J K，Shaw L．Modeling repairable systems with failure rates that depend on age and maintenance［J］．IEEE Trans．on Reliability，1993，42：566－570．

［4］ H Wang，H Pham．Some maintenance models and availability with imperfect maintenance in production systems［J］．Annals of Operations Research，1999，91，（0）．

［5］ 徐文静．不完全维修条件下的可用度与维修策略分析［D］．湖南：国防科技大学，2008．

［6］ 程志君，高大化．不完全维修条件下的视情维修优化模型［A］．系统工程与电子技术，2006，28（7）：1106－1108．

［7］ 徐文静，郭波．基于回复因子的不完全维修模型分析［A］．科技创新导报，2008（24）：74－75．

［8］ 张宁平，陈跃良．基于可靠度经济性约束的不完全维修周期研究［A］．新技术新工艺，2010，3：17－19．

［9］ Ebastian Martorell，Ana Sanchez，Vicente Serradell．Age－dependent reliability model considering effects of maintenance and working conditions［J］．Reliability Engineering and System Safety，1999，64，64（1）：19－31．

［10］ Pham H，Wang H．Optimal maintenance policies for several imperfect maintenance models［J］．International Journal of Systems Science，1996，12，12（2）：87－101．