考虑不同维修效果的多状态可修系统可靠性模型

狄鹏，黎放，陈童

(海军工程大学 管理工程系，湖北 武汉430033)

0 引言

随着军事装备功能与结构的日益复杂，系统功能与故障模式表现出多样性的特点，将系统状态简单分为“运行”和“故障”的传统可靠性分析方法已经显现出很大的局限性［1-2］。因此，多状态系统于20世纪70年代被提出［3］，到20世纪80年代初步建立了多状态系统的可靠性理论。

多状态系统可靠性解析模型研究作为一个重要方向，受到国内外学者的广泛关注［4］。如Amico 等［5］利用非齐次半马尔可夫过程研究了具有M+1 个工作状态的系统可靠性问题。Lisnianski 等［6］利用多状态马尔可夫模型研究了具有多性能状态的大型发电机组的可靠性规律，其中假设系统工作状态的转移服从指数分布，采用参数估计法获取各转移速率。Levitin 等［7］利用通用函数法和可靠性框图研究了多状态串并联系统的可靠性规律。Moghaddass等［8］采用参数估计法对具有连续状态劣化情况的环境监测设备的可靠性规律进行了研究。王丽花等［9］利用马尔可夫过程研究了离散状态退化系统的可靠性问题，维修策略选择小修和一般型更换，部件更换时间服从一般分布。Yuan 等［10］利用几何过程研究了存在性能退化的可修系统可靠性问题，模型假设系统有一个多重休假模式的维修工，系统故障时间和维修工休假时间服从指数分布，维修时间为一般分布。

上述研究具有如下特点:解析建模过程复杂，模型假设条件严格，模型只适用于特定对象，重用性不佳。为了改善这种情况，Phase-type(PH)分布被引入到多状态可修系统可靠性解析建模工作中。Frostig 等［11］研究了具有连续劣化和随机冲击的可修系统，冲击烈度服从PH 分布。Ruiz-Castro 等［12］研究了具有多状态部件的n 中取k 系统，部件寿命和维修时间均服从离散PH 分布。Segovia 等［13］研究了具有外部冲击和内部磨损的多状态系统可靠性，冲击到达间隔时间和冲击损坏效果均服从连续PH 分布。

PH 分布作为指数分布的矩阵形式推广，保持了指数分布的易处理性，其密度函数、拉普拉斯变换和各阶矩等都具有简洁的矩阵表示和概率意义。而由PH 分布在非负实数轴上全体概率分布函数中的稠密性所决定，对于任何非负随机变量，总可以找到一个PH 分布把该随机变量逼近到任意需要的精度。即对于未知分布规律，假设为PH 分布是合理的。因此采用PH 分布作为解析模型假设条件能够有效提升模型的描述能力。此外，PH 分布在很多情况下使模型分析中复杂的数值积分转化为矩阵运算，有利于提升模型的可计算性［14］。

本文对多状态可修系统的一类常见现象展开研究，即这类系统通常可以根据实际需要选择不同方式的维修，维修效果也随之不同。本文考虑修理后恢复到修前状态、维修效果较差和维修效果好3 种效果，假设系统状态可以划分为完好、一般和故障3 种状态集，其中系统在完好和一般状态的停留时间为系统工作时间。模型假设系统工作时间和维修时间均服从连续PH 分布，通过建立系统状态马尔可夫转换的无穷小生成元矩阵，求得各状态的稳态概率向量，给出系统稳态可用度、稳态故障率的解析表达式，并通过算例验证模型的有效性和适用性。

下面首先对连续PH 分布作简要介绍。

1 连续PH 分布

定义［15］［0，+∞)上的概率分布函数F(·)称为PH 分布，并具有m 阶(α，T)表示，当且仅当它是一个有限状态马尔可夫过程的吸收时间分布，有分布函数F(x)=1 -αexp (Tx)eT.

该马尔可夫过程具有状态集{1，…，m，m +1}，状态1，…，m 都是非常返的(瞬态)，状态m +1 是吸收态。初始概率为(α，αm+1)，其中α=(α1，…，αm)且满足αeT+αm+1=1，e 为元素值均为1 的行向量，eT为e 的转置。该马尔可夫过程无穷小生成元Q 可写成下列分块形式:式中:T =(Tij)为m 阶方阵，Tij≥0 (i≠j;i、j =1，…，m)表示从瞬态i 至瞬态j 的转移率，且满足Tii＜0;T0=(T01，…，T0m)T是非负列向量，其中T0i 表示从瞬态i至吸收态的转移率，满足TeT+T0=0. PH 分布中的每一个瞬态称为相位。

2 问题描述

某多状态系统的状态集为S={1，…，n，n+1}.根据系统工况，可以将S 分为3 类:完好状态，即系统性能保持在设计指标附近，有状态集G ={1，…，k};一般状态，即系统可以正常运行，但性能较差，有状态集B={k +1，…，n};故障状态，指系统发生停机故障，有状态集F ={n +1}，此时必须经过修理，系统才能恢复到状态集G 或B.

由于使用和环境等因素影响，系统会从完好和一般状态直接进入故障状态。而本文考虑3 种维修效果:

1)修后恢复修前状态，即系统经过维修后，返回到发生故障前一刻的状态，如图1所示。

图1 修后恢复修前状态示意图Fig.1 System state being the same as good/old after repair

2)维修效果较差，即系统经过维修，只能恢复到一般状态集B，如图2所示。

图2 维修效果较差的系统状态变化示意图Fig.2 Transition of system state for imperfect repair

3)维修效果好，即系统经过维修，恢复到完好状态集G，如图3所示。

图3 维修效果好的系统状态变化示意图Fig.3 Transition of system state for perfect repair

对该系统做进一步假设:

1)系统工作时间，即系统在状态集G 和B 的停留时间服从PH 分布，(α，T)为该PH 分布的n 阶不可约表示，其中α 和T 可以表示为分块形式:α=

矩阵TGG、TGB、TBG和TBB分别表示系统在状态集G 和B 内部和相互之间的转移速率。通常不经过维修时，系统不能从一般状态转换到完好状态，故TBG=0.

工作时间吸收速率矩阵T0亦可写为分块形式:

2)系统维修时间，即系统在状态集F 的停留时间服从PH 分布，(β，S)为该PH 分布的m 阶不可约表示。

因此，有PGGe+PGBe=e，PBGe+PBBe=e.

3 模型分析

3.1 无穷小生成元矩阵

根据系统状态变化情况，需定义4 个宏状态:OG、OB、RG、RB. 系统状态空间为Ω=OG∪OB∪RG∪RB，其中:OG={(i)，i∈G}表示系统处于完好状态的相位i;OB={(i)，i∈B}表示系统处于一般状态的相位i;RG={(i，j)，i∈G，j=1，…，m}表示系统故障前处于完好状态的相位i，维修状态处于相位j;RB={(i，j)，i∈B，j=1，…，m}表示系统故障前处于一般状态的相位i，维修状态处于相位j.

因此，该马尔可夫链的无穷小生成元矩阵可以表示为

式中运算符⊗为Kronecker 积［14］，下面分别从4 个方面说明Q 中各元素的构成。

3.1.1 宏状态OX向OY(X，Y∈{G，B})的转移

以宏状态OG向OG的转移为例，系统在OG内部进行转移，因此相位i(i∈G)之间的转移率矩阵为TGG. 同理，可知宏状态OB内部，以及宏状态OG和OB之间的转移率矩阵。

3.1.2 宏状态OX向RY(X，Y∈{G，B})的转移

以宏状态OG向RG的转移为例，说明系统进入故障状态的前一刻处于G，并以初始概率β 进入维修相位，因此有diag(T0G)⊗β，其中diag(T0G)为将单位矩阵的对角线元素换为T0G各分量。同理，可得OB向RB的转移率矩阵。

以宏状态OB向RG的转移为例，因为系统进入故障状态的前一刻处于B，不可能进入宏状态RG，因此为0 矩阵。同理，可得OG向RB的转移率矩阵。

3.1.3 宏状态RX向OY(X，Y∈{G，B})的转移

以宏状态RG向OG的转移为例，说明维修活动进入了吸收态，同时以概率矩阵PGG返回G，因此有PGG⊗S0. 同理，可得RG向OB、RB向OG、RB向OB的转移率矩阵。

3.1.4 宏状态RX向RY(X，Y∈{G，B})的转移

以宏状态RG向RG的转移为例，系统在维修状态内部发生相位转移，即I⊗S. 同理，可得RB向RB的转移率矩阵。

因为不可能发生RG向RB、RB向RG的转移，故相应的转移矩阵均为0 矩阵。

3.2 稳态概率向量

当系统进入稳态时，由连续时间马尔可夫过程稳态概率向量定义［15］，无穷小生成元Q 矩阵中各个宏状态所对应的概率组成了稳态概率向量π=(πOG，πOB，πR，πRB)，并且π 满足如下方程组:

将上述方程组展开，得(2)式～(6)式:

由(4)式和(5)式可知:

将(7)式、(8)式代入(2)式、(3)式和(6)式，即可求得π.

3.3 系统特性

3.3.1 系统稳态可用度

在获得π 后，可以确定系统稳态可用度，即系统处于宏状态OG和OB的概率:

3.3.2 系统稳态故障率

故障率是系统分别从完好和一般状态进入维修状态的速率，当系统进入稳态后，稳态故障率为

3.4 维修效果对系统的影响

下面考虑模型的特例:系统只存在一种维修效果的情况。

3.4.1 修后恢复修前状态

令向量δG为从维修状态进入完好状态各相位的概率，向量δB为从维修状态进入一般状态各相位的概率。则PGG=eδG，PGB=0，PBG=0，PBB=eδB，并有

3.4.2 维修效果较差

此时PGG=0，PGB=eδB，PBG=0，PBB=eδB，并有

3.4.3 维修效果好

此时PGG=eδG，PGB=0，PBG=eδG，PBB=0，并有

4 算例

4.1 模型有效性验证

某型舰用电站系统具有5 种工作状态，如表1所示。

表1 电站系统工作状态表Tab.1 Operating states of power station

将该型电站系统功率在50%以上的工况划分为完好状态，其余工况归为一般状态，故障停机时为故障状态。该系统初始状态概率向量为α=(1，0，0，0，0)。对该型电站系统长期的运行和维修记录等数据进行分析，可得各工况之间的状态转移率矩阵为

维修时间分布为

系统维修后状态转移概率矩阵为

由3.2 节得出系统稳态概率向量后，根据(9)式和(10)式，可得系统稳态可用度A =0.961 6，系统稳态故障率r=0.255 7 次/1 000 h.

该算例说明利用本文模型能够有效获得该类多状态可修系统的可靠性参数。下面通过改变模型输入的随机分布类型，验证模型的适用性。

4.2 模型适用性验证

当系统维修时间服从指数分布、Weibull 分布等典型分布时，可以将这些分布表示为PH 形式［16］，再利用本文模型即可求得系统可靠性。首先假设系统维修时间服从修复率为μ 的指数分布，则β =(1)，S=(-μ). 令μ 在区间(0，20］中变化，可以计算得到修复率分别与系统稳态可用度、系统稳态故障率之间的关系，如图4和图5所示。

图4 修复率与系统稳态可用度关系图Fig.4 Relationship between repair rate and stationary availability

图5 修复率与系统稳态故障率关系图Fig.5 Relationship between repair rate and stationary failure rate

图4说明修复率取值的增大能够有效增加系统稳态可用度。对于本文研究的多状态可修系统，在系统进入稳态之前，其故障率并非常数，只有当系统进入稳态后，在给定修复率下，系统的稳态故障率才会固定。而随着修复率取值的增大，系统稳态故障率也逐步上升，并最终趋于稳定，如图5所示。这是因为修复率直接影响到系统在完好状态和一般状态的停留时间，即稳态概率πOG和πOB随着修复率取值的增大而增大，结合(10)式可知，r 也将增大。图5正是反映了r 与μ 之间的这种正相关关系。而r 最终趋向于某一稳定值，这正是在忽略维修时间情况下系统长期运行时，系统固有的故障率。

当维修时间为Weibull 分布时，令a 为形状参数，b 为尺度参数。采用文献［16］方法将各Weibull分布拟合为PH 分布形式，然后带入本文模型可得对应系统可靠性参数如表2所示。

从表2可以看到，当给定形状参数a 后，尺度参数b 与系统稳态可用度A 和系统稳态故障率r 均表现出正相关关系，分别如图6和图7所示，可以清楚地看到:b 取不同的值时，a 对A、r 的影响机理并不会因为b 值不同而发生改变，显现出本文算法的稳定性。

算例说明本文模型能够方便地用于不同随机分布类型的输入，适用范围有所拓展。

表2 维修时间服从Weibull 分布时的系统可靠性Tab.2 System reliability measures with Weibull repair time

图6 形状参数a 对系统稳态可用度的影响Fig.6 Influence of shape parameter a on system stationary availability

图7 形状参数a 对系统稳态故障率的影响Fig.7 Influence of shape parameter a on system stationary failure rate

5 结束语

本文采用PH 分布研究了具有不同维修效果的多状态可修系统可靠性规律，在保证良好解析特性的同时，有效提升了模型的描述能力。本文模型在计算时主要涉及矩阵运算，而目前高性能计算机和矩阵解析理论的应用能对大型矩阵的运算提供良好的支持，算例充分演示了模型的良好计算性和适用性。

References)

［1］ Zio E. Reliability engineering:old problems and new challenges［J］. Reliability Engineering and System Safety，2009，94:125 -141.

［2］ Ying K G，Jing L. Multi-state system reliability:a new and systematic review［J］. Procedia Engineering，2011，29:531 -536.

［3］ Barlow R E，Wu A S. Coherent systems with multi-state components［J］，Mathematics of Operations Research，1978，3:275 -281.

［4］ Lisnianski A，Frenkel I，Ding Y. Multi-state system reliability analysis and optimization for engineers and industrial managers［M］，London:Springer，2010:8 -19.

［5］ Amico G D，Biase G D，Manca R. A customer’s utility measure based on the reliabilityof multi-state systems［J］. Decisions in Economics Finance，2011，34(1):1 -20.

［6］ Lisnianski A，Elmakias D，Laredo D，et al. A multi-state Markov model for a short-term reliability analysis of a power generatingunit［J］. Reliability Engineering and System Safety，2012，98(1):1 -6.

［7］ Levitin G，Xing L. Reliability and performance of multi-state systems with propagated failure shaving selective effect［J］. Reliability Engineering and System Safety，2010，95:655 -661.

［8］ Moghaddass R，Zuo M. A parameter estimation method for a condition-monitored device under multi-state deterioration［J］. Reliability Engineering and System Safety，2012，106:94 -103.

［9］ 王丽花，岳德权，张静，等. 具有小修和一般型更换策略的多状态退化系统［J］. 辽宁工程技术大学学报:自然科学版，2011，30(6):935 -938.WANG Li-hua，YUE De-quan，ZHANG Jing，et al. A multi-state degraded system with minimal repairs and the replacement policy in general distribution［J］. Journal of Liaoning Technical University:Natural Science，2011，30(6):935 -938.(in Chinese)

［10］ Yuan L，Xu J. A deteriorating system with its repairman having multiple vacations［J］. Applied Mathematics and Computation，2011，217(10):4980 -4989.

［11］ Frostig E，Kenzin M. Availability of inspected systems subject to shocks-a matrix algorithmic approach［J］. European Journal of Operational Research，2009，193(1):168 -183.

［12］ Ruiz-Castro J E，Li Q L. Algorithm for a general discrete k-outof-n:G system subject to several types of failure with an indefinite number of repairpersons［J］. European Journal of Operational Research，2011，211(1):97 -111.

［13］ Segovia M C，Labeau P E. Reliability of a multi-state system subject to shocks using phase-type distributions［J］. Applied Mathematical Modelling，2013，37(7):4883 -4904.

［14］ 陈童. 基于PH 分布和马尔可夫到达过程的装备备件需求与库存模型研究［D］. 长沙:国防科学技术大学，2010.CHEN Tong.Research on demand and inventory models of equipment spare parts based on phase-type distribution and Markovian arrival process［D］. Changsha:National University of Defense Technology，2010.(in Chinese)

［15］ 田乃硕. 休假随机服务系统［M］. 北京:北京大学出版社，2001:5 -7.TIAN Nai-shuo. Queuing systems with sever vacations［M］，Beijing:Peking University Press，2001:5 -7.(in Chinese)

［16］ 黄卓.Phase-type 分布数据拟合方法及其应用研究［D］. 长沙:国防科学技术大学，2007.HUANG Zhuo. Research on phase-type distributions data fitting method and its applications［D］. Changsha:National University of Defense Technology，2007.(in Chinese)