考虑非完美维修的多部件系统劣化状态空间划分建模

石冠男，张晓红，曾建潮,2+

(1.太原科技大学 工业与系统工程研究所，山西 太原 030024； 2.中北大学 计算机与控制工程学院，山西 太原 030051)

0 引言

维修是使设备保持、恢复或改善其运行状态的全部技术和管理活动[1]。优化维修策略的目的是在保证系统可靠性的前提下，尽可能降低系统运行成本，以获得更好的经济效益。如今工业生产设备的结构越来越复杂，各组成部件之间存在复杂的关联关系，为了使维修决策模型更加符合实际需要，考虑部件自身特点和部件间依赖关系的多部件系统维修决策模型研究具有重大的理论价值和现实意义[2]。

在多部件系统的最优维修决策建模中，系统的维修需求表示及其概率计算是建立维修优化目标与决策变量联系的基础，而系统劣化状态的概率密度函数又是计算维修需求分组概率的先决条件。2015年，Zhang等[3-4]针对由相同和不相同部件组成的多部件系统，先后提出统一的多部件系统维修建模方法，给出了基于视情机会维修的维修需求划分模型及基于该模型的维修需求概率计算通式。

目前，已有的多部件系统联合劣化状态空间划分模型在建立过程中假设部件经预防性维修后“修复如新”，这种维修方式通常被称作完美维修，而现实中，大多数设备在维修后是不能修复如新的，这种维修方式通常被称作非完美维修[5]。受维修器材、维修能力、维修人员等诸多因素影响，非完美维修更贴切实际并已被实际应用，如煤炭运输系统保养[6]、自卸卡车运输车队维护[7]、发电厂维修养护[8]、汽车产品保证服务[9]等。非完美维修的研究工作最早可追溯到Chan等[10]的最坏预防维修(worst PM)； Chaudhuri等[11]在1977年首次提出非完美预防维修(imperfect PM)概念； Pham等[12]针对非完美维修做了大量研究和文献综述，将之前的非完美维修模型总结为(p,q)法则、(p(t),q(t))法则、改善因子、虚拟年龄法则、冲击振动模型、(α,β)法则、复合(p,q)法则及其他8类，其中，(p,q),(p(t),q(t))和复合(p,q)法则利用劣化状态回退概率来体现维修效果，而改善因子、虚拟年龄法则、冲击振动模型则是假设部件经维修后表征部件劣化状态的特征参数恢复到某一特定值，其受多种因素共同影响，如维修方式、维修次数或者劣化程度等。近年来，针对非完美维修的研究也有很多，例如Li等[13]建立了基于可靠度的非完美维修模型；Zhang等[14]将基于可靠性的非完美预防维修和机会维修相结合，研究了风力发电机组的维修与维护；Nodem等[15]用固定系数体现非完美维修效果，即设备维修后恢复至现在劣化水平的固定百分比； Liu等[16]采用同时考虑退化状态和退化速度的改善因子建立非完美维修模型。在以往的非完美维修研究中，大部分文献只针对两部件系统[17]进行讨论，而对于部件数量较多或者需要同时考虑其他影响因素的复杂问题，多借助仿真方法[18]解决，因此如何用比较合理的解析方法分析考虑非完美维修的多部件系统中不同视情维修策略的组合优化，是一个值得深入研究的问题。

考虑非完美维修的多部件系统劣化状态空间划分模型，能够更准确地描述实际维修情况而更具有参考应用价值，本文在前人研究基础上，假设机会维修和预防维修区域的维修效果均为非完美的，系统部件的劣化状态在维修后介于“全新”和“如旧”之间，针对连续劣化多部件系统，建立了考虑非完美维修的多部件系统视情机会维修策略的劣化状态空间划分模型，通过分析劣化状态空间划分与维修需求分组之间的对应关系，给出了维修需求分组的概率计算通式，最后重新推导了考虑非完美维修时相同多部件系统的联合稳态概率密度函数并给出数值解法，利用数值实验验证了上述模型的正确性及对后续维修优化研究的有效性。

1 系统的特征假设及定义

1.1 系统描述

假设某一多部件系统是由M个相同且相互独立的单部件组成，每个部件均独立遭受连续的累积劣化，其表征部件劣化状态的特征参数(机械部件磨损量、齿轮裂纹长度等)在系统工作过程中逐渐增加，此类连续劣化系统在生活中十分常见，如大型风力发电场、重型商业运输车队[19]等。基于连续劣化系统的特性，相关假设如下：

(5)设备检测和维修时间均忽略不计。

1.2 维修策略

对部件采取视情机会维修，部件的整个寿命周期内可能的维修需求有不需要维修(正常运行)、非完美维修(包括机会维修和预防性维修)和故障后更换3类。为系统中所有相同的部件均定义相同的机会、预防和故障后维修阈值Do，Dp，Df，通常认为0

(4)对不符合上述3种情况的部件，不采取任何维修措施。

如此，部件i从初始状态开始正常运行，可能经过若干次机会或预防性维修，最后由故障后更换重新回到初始状态。

2 考虑非完美维修的联合劣化状态空间划分

2.1 劣化状态空间划分

图1所示为仅考虑完美维修时部件i的劣化过程示意图及维修区域划分。其中，设定的3条维修阈值将其劣化状态空间划分为4个维修活动区域：正常运行区、机会维修区、预防维修区和故障后更换区。在每次检测后根据当前劣化状态信息所落区域，采取相应的维修活动，即不需要维修、机会维修、预防维修、故障后更换。

在考虑非完美维修的单部件系统维修建模中，同样为部件i设定3个维修阈值将其劣化状态空间分为4个维修活动区域：正常运行区域、机会维修区域、预防维修区域、故障后更换区域。不同的是，为了使研究更符合实际维修需求，当设备的劣化状态在检测点落在机会维修和预防维修区域时，对部件采取非完美维修。图2所示为考虑非完美维修时部件i可能的劣化示意图，从全新初始状态开始，其可能接受数次机会维修或预防维修，维修效果为非完美的，最终经由故障后更换重新恢复至初始状态。

图3所示为考虑非完美维修的单部件系统劣化状态空间划分简化图，其中：U0表示部件没有维修需求，U1表示部件在等待机会维修，P表示部件需要预防维修，C表示部件需要故障后更换。图4所示为相同两部件系统联合劣化状态空间划分示意图,在M个部件组成的多部件系统中，各部件的维修阈值相互交叉，将其联合劣化状态空间划分为更多区域，不同区域对应不同的维修需求组合。

类似地，可以得出M=3,4,…时相同部件系统联合劣化状态空间划分。对于不相同部件系统而言，区别仅在于各部件的维修阈值不同，其劣化状态空间划分思想不变。

在具体的维修优化建模中，研究者通常关注的是求取多部件系统的维修需求组合和相应的概率，而系统劣化状态的稳态概率密度函数是计算维修需求分组概率的先决条件，下面在假设系统联合劣化状态的概率密度函数存在且为π(x1,x2,…,xM)的情况下，给出M个相同部件系统中所有维修需求组合的概率计算通式。

2.2 考虑非完美维修的多部件系统维修需求分组及其概率计算

(1)当i=0,j=0时，

(1)

(2)当i=0,1≤j≤M时，

(2)

(3)当1≤i≤M,j=0时,

(3)

(4)当1≤i

(4)

3 考虑非完美维修的劣化状态回退

(5)

至此，得到相同多部件系统中所有维修需求分组及其概率计算通式，并对非完美维修进行了合理建模。在相同多部件系统维修中引入非完美维修虽然不会影响其维修需求分组，但是会改变部件劣化状态的概率密度函数，进而影响系统联合劣化状态的概率密度函数，因此需要对考虑非完美维修的多部件系统联合劣化状态的概率密度函数进行推导求解。

4 考虑非完美维修的多部件系统稳态联合概率密度计算

在单部件维修决策过程中常利用更新理论进行建模，但在基于视情机会维修策略的多部件系统维修中，特别是引入非完美维修后，所有部件在某一个检测时刻全部恢复至初始全新状态的可能几乎不存在，寻找多部件系统的完全更新周期无从谈起。然而，考虑到系统中部件的后续劣化仅与当前检测时刻的状态有关，与系统的历史状态无关，符合Markov特性，可利用半更新特性将系统检测周期作为一个半更新周期，分析一个半更新周期内联合劣化状态的稳态概率密度函数，可近似得到系统联合劣化状态在整个劣化空间的稳态概率密度函数[23-24]。

已知π(x1,x2,…,xM)为M相同多部件系统的联合劣化状态的概率密度函数，设任意部件i的劣化状态稳态概率密度函数为Ω(xi)，因为各部件的劣化过程相互独立，所以只需求出任意部件i的劣化状态稳态概率密度函数Ω(xi)，即可得出系统联合劣化状态的稳态概率密度函数π(x1,x2,…,xM)=Ω(x1)Ω(x2)…Ω(xM)，而部件i的劣化状态稳态概率密度函数Ω(xi)可通过分析半更新周期内部件i所有可能的维修及其状态转移情况得到。

在多部件系统维修中分析任意部件i的维修需求，需要在决策时刻考虑是否有其他部件为部件i提供维修机会。因此，在M相同的多部件系统中，分析部件i需要考虑以下两种情形：

情形1当其他M-1个部件均没有维修需求，即其他M-1个部件部件的劣化状态均落在区域[0,Dp)时，由式(1)可得其发生概率为

(6)

在上述基础上，假设半更新周期开始前的部件劣化状态为yi，继续分析部件i所有可能的维修情况及其转移状态，在半更新周期开始时部件i又有维修和没有被维修两种情况(如图5)：

(1)部件i在半更新周期开始时没有被维修过

(2)部件i在半更新周期开始时进行了维修

可得，在其他M-1个部件不会为部件i提供机会维修机会时，部件i的劣化状态为xi的稳态概率密度函数可表示为

(7)

情形2当其他M-1个部件至少有一个有维修需求时，其会为部件i提供机会维修机会，其发生概率为

dx1dx2…dxM-1。

(8)

类似情形1，继续分析此时部件i所有可能的维修及其转移状态情况，如下：

(1)部件i在半更新周期开始时没有被维修过

(2)部件i在半更新周期开始时被维修过

在其他M-1个部件会为部件i提供维修机会时，部件i的劣化状态为xi的稳态概率密度函数可表示为

(9)

综合情形1和情形2的发生概率及上述所有状态的转移可能，合并式(7)和式(9)，可得部件i的劣化状态为xi的稳态概率密度

(10)

式(10)关于Ω(xi)的表达式是一个隐式积分方程，很难精确求出其精确解析解，下文尝试采用离散方式求其数值的近似解。

5 数值求解

5.1 数值解法

首先对式(10)进行合并项简化：

(11)

f(x-y*))。

(12)

根据定积分的意义，每一个积分和都可以看作是定积分的一个近似值，即

(13)

这在几何意义上，是用一系列小矩形面积近似小曲边梯形的结果，将所示的正交近似数值求解规则应用于式(12)的求解，同时采用Dmax对式(12)中的∞进行截尾，每一个正交点的近似方程可表示为

wherex/h=i,y/h=j,z/h=k,

Do/h=a,Dp/h=b,Df/h=c,Dmax/h=d。

(14)

将要求解的稳态概率密度函数用矩阵形式表示为Ω=[Ω(h)Ω(2h) …Ω(ih)]T，然后根据矩阵相乘的原则将式(14)改写为式(15)的矩阵形式：

Ω=hK1+h2K21K22+hK3+

(hK4)M-1(hK5+h2K61K62)。

(15)

式中K1,K21,…的具体定义见附录。如此可通过求解Ω(x)的线性方程组得到其数值近似解。

5.2 数值实验

5.2.1 数值求解

目前,大多数研究多采用Markov过程(如复合泊松分布、Gamma分布、布朗漂移运动等)描述系统的劣化过程，Gamma分布最适合用来分析非负增量、严格单调的劣化过程[25]，并有利于进行数学上的解析推导，文献[26-27]以Gamma过程描述系统的劣化过程，研究了考虑非完美维修劣化系统的视情机会维修优化模型；李玲等[28]将Gamma过程与几何过程相结合，研究了加速劣化系统的最优视情维修策略。因此，本文假设部件单位时间内的劣化状态变化量Δxi服从Gamma分布。

设部件单位时间内的劣化状态变化量Δxi服从概率密度函数为式(16)的Γ(α,β)分布，则t个单位时间内的累积增量服从Γ(tα,β)分布。

(16)

式中：α为形状参数，β为尺度参数。

令M=10，Do=8，Dp=12，Df=20，用5Df对积分上限∞作截尾处理。图6和图7所示为β=4对应不同α和α=0.04对应不同β时，求解非线性方程组得到的Gamma函数f(x)和稳态概率密度函数Ω(x)的近似数值解。

5.2.2 正确性分析

表1 不同劣化参数对应的R

≈1，

(17)

即π(x1,x2,…,xM)在整个联合劣化状态空间的积分同样近似为1，满足概率密度函数的一般特性。

5.2.3 有效性分析

在实际维修优化建模时，通常选取系统的半更新周期和维修阈值作为优化的决策变量，不同的半更新周期和维修阈值会影响稳态概率密度函数的求解，下面分析半更新周期t和各维修阈值对稳态概率密度函数的影响程度。选取Do=8，Dp=12，Df=20，α=0.04，β=4，t=100为标准对照参照组，每次仅改变参照组中一个参数，以分析各参数对稳态概率密度函数的影响。

图8所示为当半更新周期t发生改变，其余参数不变时，劣化增量分布Gamma函数f(x)和稳态概率密度函数Ω(x)的示意图。从图中可以看出，随着半更新期周期t的增加，t时间内累积劣化增量分布Γ(tα,β)的形状参数明显增加，故Gamma函数f(x)和稳态概率密度函数Ω(x)的均值明显右移。造成该现象的原因是Gamma分布的特性：在固定β值时，随着α值的增大，其超过某定值的概率会随之增加，对于本文所建立的视情机会维修策略，体现为提高检测周期，会提高设备劣化状态到达失效阈值的可能性。

图9～图11所示分别为相同劣化增量分布Gamma函数f(x)，以及不同的机会维修阈值、预防维修阈值和故障后维修阈值对应的稳态概率密度函数Ω(x)示意图。

从图9～图11的分析来看，选择不同的检测周期、机会维修阈值、预防维修阈值和故障更换阈值将对应不同的维修概率，从而对后续维修优化模型的维修成本产生影响，因此该理论方法可以为后续的维修优化决策提供理论依据。实际维修中，检测周期过长使检测间隔内发生停机故障的概率增大，增加了故障维修成本，检测周期过短又会对系统干预频繁，从而大大增加检测成本。同样，在故障更换阈值确定的情况下，减小机会和预防维修区域，会缩短部件的生命周期，使预防维修成本过高，反之，增大机会和预防维修区域则会增大周期内故障维修的概率和成本。

6 结束语

本文主要针对设备维修后不能修复如新这一实际情况，将已有的劣化状态空间划分方法与非完美维修模型相结合，建立了考虑非完美维修的多部件系统劣化状态空间划分模型，在假设系统联合概率密度函数存在的情况下，给出了各维修需求分组的概率计算式，然后利用半更新理论求得系统的联合稳态概率密度函数表达式，并给出了相应的数值解法，最后通过数值实验分析验证了上述稳态概率密度函数定义及解法的正确性。

本文还分析了该方法对后续维修优化建模的有效性，后续将以此为基础对相同多部件系统的最优维修优化决策展开进一步研究;同时,本文仅考虑定期检测下的维修策略劣化状态空间划分，事实上，由于劣化程度的加深，系统劣化速度会加快，采用序贯检测更加合理，后续将以本文为基础进行深入研究。