尹东亮,吴 晶,黄晓颖

(海军工程大学 作战运筹与规划系, 湖北 武汉 430033)

军事装备或大型舰船、飞机关键系统结构复杂、可靠性要求高,通常采用冗余储备结构提高系统效能,如舰船动力系统、飞机电气系统等,往往采用n中取k这一冗余结构[1-4]。

但新装备组成单元繁杂、集成程度高、维修经验缺乏,导致故障后维修难度大。同时,大量新装备的列装使用,对装备维修保障带来了较大的压力,所需维修人员专业性较强,培训要求较高,维修力量较为紧缺。装备部门的修理所及军工厂维修力量通常不会“一对一”专职保障,往往是一组维修力量同时负责维护若干台装备[5],有时装备出现故障后不能第一时间得到有效检修与维护,影响装备使用效能的提升。因此,在装备数量与维修力量难以平衡的维修环境下,应立足于实际情况制定严谨科学且具有良好可行性的维修力量调度方案,这样既能够有效改善装备使用可靠性,还能够使装备的运行周期进一步延长,从而提高装备的战斗力。

关于冗余储备系统的维修力量调度问题,Takagi[6]于1991年提出修理工可变休假策略将其应用于系统可靠性建模中。该策略是指一旦维修台空出,修理工要求进行连续H次休假,H为确定的任意正整数。若在这H次休假内系统所有部件始终处于完好状态,第H次休假结束后修理工进入空闲期,直到系统内有故障件到达转入维修期。若系统在修理工H次休假内出现故障件,则修理工将在该次休假结束后终止而转入维修期,“休假”是指系统部件故障后不能立即得到维修的现象。该策略体现了维修调度管理者在修理工维修工作与休假期从事的辅助工作之间权衡得失来确定休假次数思想,更具一般性。

在考虑该种策略的系统可靠性建模研究中,Ke[7]针对M/G/1排队系统,假设服务台计划休假次数为J次,以可变休假为休假策略开展了一系列研究。Wang等[8]假设服务到达服从伯努利到达过程,服务台在基本休假后最多可休假J-1次,利用通用生成函数得出了一些装备性能指标,演示了以最小成本确定J最佳值的成本模型。Guha等[9]研究了服务台执行多个单一工作休假策略的情况,在服务时间、休假时间这两大参数皆满足指数分布要求的前提下,科学合理地确定了计划休假次数的均衡阈值。Jain等[10]针对含M个工作部件、k个备用部件的多部件制造系统,假设执行N次休假策略的修理工休假流服从泊松过程,部件寿命服从指数分布,得出了该系统的稳态概率向量。

在上述研究过程中,为简化分析,往往会假设部件寿命等相关指标满足指数分布要求等,而这容易造成模型计算复杂,同时也会降低模型的适用性。为克服此不足,Neuts等[11]率先将PH(Phase-type)分布以指数分布的一般形式进行了研究和论述。即便是在进行一系列复杂的运算后,PH分布族依旧表现出较强的封闭性,可以比较准确合理地拟合非负半轴上的各种概率分布情况,同时还可将各种不同的典型分布统一描述为PH分布形式,既能降低计算复杂度,还能增强适用性,令其在各类随机模型中得到灵活合理的应用。

Goswami等[12]研究了单服务台工作休假排队系统,服务员服务速率随顾客数变化而变化,假设服务时间、休假时间等均服从PH分布,分析了该系统的一系列性能指标。Banik[13]以考虑服务员可变休假的单服务台排队系统为研究对象,利用批量马尔可夫到达过程(batch Markovian arrival process, BMAP)表示顾客到达流,假设服务员服务时间、休假时间满足PH分布要求,科学合理地确定了系统稳态队长,并对服务台单位时间平均费用进行了优化。Liu等[14]分析了修理工工作休假且休假可中断的冷贮备可修系统,连续休假期表示为马尔可夫到达过程(Markovian arrival process, MAP),假设部件寿命、工作休假时间、维修时间服从PH分布,得出了系统稳态可用度等可靠性指标。Sheng等[15]基于休假策略研究了具有不确定性的多状态n中取k系统,引入重要性测度评估了系统可靠性参数。

因此,本文考虑单一维修台的修理工可变休假,选取较为典型的n中取k系统对大型复杂装备可靠性模型进行研究。假设系统工作部件寿命、备用部件寿命、修理工休假时间、维修时间以及空闲时间均服从不同的连续PH分布,构建了具有广泛适用性的系统可靠性解析模型,提供了包括系统可靠度、稳态故障频度等在内的一系列可靠性指标;最后通过算例验证了模型的适用性,模拟分析了修理工有无休假、修理工休假速率等相关可变因子对系统运行可靠性指标产生的各种影响。

1 连续PH分布

定义1[16]如果非负随机变量X满足连续PH分布要求,其分布函数可描述如下:

(1)

其中:t≥0,e是一种列向量,与公式相匹配的各类元素均为1;α为次随机矩阵,存在m个非负元素的行向量,m为正整数,αe≤1;T为m阶矩阵,其对角线元素均呈负,其他元素为正,矩阵可逆。

定义2[16]连续时间马尔可夫链{I(t),t≥0},其状态空间包含m+1个状态{1,2,3,…,m,m+1},无穷小生成元如下所示:

(2)

其中,T是PH分布的生成元。因所有行的元素之和皆是零,可推导出T0=-Te。状态m+1的转移可能性为零,状态m+1为吸收态。PH分布的一般描述形式为m阶的(α,T)表示。

定义3[17]m×n阶矩阵A和p×q阶矩阵B的Kronecker积被描述如下:

(3)

通过式(3)能够科学合理地推导出Kronecker积存在下述特征:

C(A⊗B)=(CA)⊗B=A⊗(CB)

(4)

(A1⊗B1)(A2⊗B2)=(A1A2)⊗(B1B2)

(5)

定义4[17]m阶矩阵A和n阶矩阵B的Kronecker和被描述如下:

A⊕B=A⊗In+Im⊗B

(6)

式中,Im、In分别为m、n阶的单位矩阵。

2 模型假设

结合装备保障中的实际问题背景构建模型,如舰船完成任务靠港后提交一系列装备维修清单,由岸上修理所进行维修保障,而修理所保障多艘舰船必定会出现维修调度问题,导致待修装备使用可靠性的降低和维修延迟。

对此类问题进行抽象假设,假设某n中取k系统包含n个相同部件,若k个或更多部件正常运行时,系统处于工作状态,否则系统停机;当系统工作时某工作部件故障,则由备用件替换,系统继续运行,故障件进入维修序列,系统内有单一维修台,进行可变休假。下面将此假设展开深入细致的探讨和分析:

1)在系统中,n个相同部件只存在两种状态:一种是完好无损状态,另一种则是异常状态,即故障或失效状态,工作部件(编号1,…,k)寿命独立同分布,满足PH分布,通过m1阶的(α,T)进行描述;备用部件(编号k+1,…,n)寿命独立同分布,满足PH分布,通过m2阶的(β,S)进行描述。

2)系统内有单一维修台,始终保持着完好的运行状态,修理工进行可变休假,即一旦维修台空出,修理工要求进行连续H次休假。若在这H次休假内系统n个部件始终处于完好状态,第H次休假结束后修理工进入空闲期,直到系统内有故障件到达,则修理工转入维修期。若系统在修理工H次休假内出现故障件,则修理工将在该次休假结束后终止而转入维修期,具体修理工调度流程如图1所示。

图1 修理工调度流程——可变休假Fig.1 Repairman scheduling process—variant vacations

3)维修台在运行过程中,遵循先到先服务的维修规则,故障件在经过一系列维修后恢复正常。修理工维修时间满足PH分布情况,通过l阶的(δ,U)进行描述,每次休假时间独立同分布,满足PH分布,通过g阶的(γ,R)进行描述,空闲时间服从PH分布,通过w阶的(θ,L)进行描述。

4)在整个系统中,如果运行正常的部件未达到k个,系统会自动停机,其中,运行正常的部件保持当前状态,不会出现异常情况,系统恢复正常后进入运行状态。

5)不考虑工作部件的替换时间。

6)工作部件寿命、维修时间等相关变量互相独立。

3 模型构建

3.1 系统状态空间

假设i指的是系统内存在的故障件数目,运行无恙的部件相位为O(t)={o1,o2,…,ok,ok+1,…,on},修理工的休假所处相位为V(t),工作所处相位为J(t),空闲所处相位为W(t),那么系统是一个多维连续时间轴上的马尔可夫过程。在对系统状态空间进行深入研究后,划分六种运行状态,状态空间为Ω={H1,H2,H3,H4,H5,H6}。

H1={i=0,(o1,o2,…,ok,ok+1,…,on),v(t)}表示系统处于正常运行状态,n个部件均完好,修理工处于H次以内的某次休假期。其中:完好工作部件所处相位为o1,o2,…,ok,1≤o1,o2,…,ok≤m1;完好备用部件所处相位为ok+1,…,on,1≤ok+1,…,on≤m2;修理工休假所处相位是v(t),1≤v(t)≤g。

H2={i=0,(o1,o2,…,ok,ok+1,…,on),w(t)}表示修理工经过了H次休假后系统仍处于n个部件均完好的正常运行状态,此时修理工处于H次休假后的空闲期。其中:完好工作部件所处相位为o1,o2,…,ok,1≤o1,o2,…,ok≤m1;完好备用部件所处相位为ok+1,…,on,1≤ok+1,…,on≤m2;修理工空闲期所处相位是w(t),1≤w(t)≤w。

H3={i≤n-k,(o1,o2,…,ok,ok+1,…,on-i),v(t)}指的是系统保持稳定可靠的工作状态,i个部件出现故障,i≤n-k,修理工处于H次内的某次休假期。完好无损的工作部件相位为o1,o2,…,ok,1≤o1,o2,…,ok≤m1;完好备用部件所处相位为ok+1,…,on-i,1≤ok+1,…,on-i≤m2;修理工休假所处相位是v(t),1≤v(t)≤g。

H4={i≤n-k,(o1,o2,…,ok,ok+1,…,on-i),j(t)}表示系统处于正常运行状态,有i个部件故障,i≤n-k,修理工H次以内的某次休假结束对故障件进行维修。其中,完好工作部件所处相位为o1,o2,…,ok,1≤o1,o2,…,ok≤m1;完好备用部件所处相位为ok+1,…,on-i,1≤ok+1,…,on-i≤m2;修理工维修工作所处相位是j(t),1≤j(t)≤l。

H5={i≥n-k+1,v(t)}指的是系统停机状态,有i个部件运行异常故障,i≥n-k+1。修理工处于H次以内的某次休假期,休假所处相位是v(t),1≤v(t)≤g。

H6={i≥n-k+1,j(t)}指的是系统停机状态,有i个部件运行异常故障,i≥n-k+1。修理工H次以内的某次休假结束,结合实际情况,按照相关规定和要求对故障件进行检修维护,维修工作所处相位是j(t),1≤j(t)≤l。

3.2 系统状态转移矩阵

基于前文论述,接下来对系统空间Ω={H1,H2,H3,H4,H5,H6}内部发生的各种转移情况展开详细分析:

H2内部转移情况分析:此次转移是修理工在进行了H次休假后,系统依旧保持稳定可靠的工作状态,各个部件皆完好无损,也未出现任何异常,n个部件自身转移,修理工处于空闲期。内部件和修理工空闲无法同步实现相位转移,故其转移矩阵可描述为Tk⊕Sn-k⊕L。

H3内部转移情况分析:在此次转移过程中,系统保持稳定可靠的工作状态,不过有i个部件出现异常,无法稳定安全地运行,i≤n-k,修理工正在休假。n-i个部件自身转移,也就是说,k个工作部件自身转移,n-k-i个备用部件自身转移,内部件无法和修理工休假同步实现相位转移,故其转移矩阵可描述为Tk⊕Sn-k-i⊕R。

H4内部转移情况分析:在此次转移过程中,系统保持稳定可靠的工作状态,不过有i个部件出现异常,无法稳定安全地运行,i≤n-k,修理工负责对所有异常部件的检修,此次维修在H次休假计划之内。n-i个部件自身转移,也就是说,k个工作部件自身转移,n-k-i个备用部件自身转移,内部件无法和修理工休假同步实现相位转移,故其转移矩阵可描述为Tk⊕Sn-k-i⊕U。

H5内部转移情况分析:在此次转移过程中,系统保持停机状态,i个部件出现异常,无法稳定安全地运行,i≥n-k+1,修理工在H次休假计划之内的休假期。只有修理工休假相位转移,其转移矩阵才可描述为R。

H6内部转移情况分析:在此次转移过程中,系统保持停机状态,i个部件出现异常,无法稳定安全地运行,i≥n-k+1,修理工负责对所有异常部件的检修,此次维修在H次休假计划之内。修理工维修工作相位转移,其转移矩阵可描述为U。

对{H1,H2,H3,H4,H5,H6}中所有状态对应的转移情况展开深入研究,基于此得出系统状态转移图,如图2所示。

图2 系统状态转移图Fig.2 State transition of system

根据图2各状态之间的转移关系,对系统状态空间Ω={H1,H2,H3,H4,H5,H6}各状态之间的转移情况进行说明:

从H1到H2的状态转移分析:该次转移是修理工经过了H次计划休假后系统n个部件始终处于完好状态,此时修理工进入空闲期,转移矩阵为(R0γ)H-1R0θ⊗[Tk⊕Sn-k]。

同理,可得其余各状态之间的转移矩阵,令0v表示修理工休假时系统完好无故障件,0id表示系统运行安全可靠,未出现任何异常部件而修理工H次休假后当前处于空闲期,iv表示修理工休假阶段异常部件的数目,ij表示修理工开展维修工作阶段异常部件的数目,此时,系统状态可用马尔可夫过程描述Ω={0v,0id,1v,1j,2v,…,(n-k)v,(n-k)j,(n-k+1)v,(n-k+1)j}。按照PH分布的概念,基于上文对系统状态转移的分析,可确定系统状态转移矩阵为:

其中,系统转移矩阵Q中各个分块矩阵的含义分别如下:

Aid=Tk⊕Sn-k⊕L

B=(R0γ)H-1R0θ⊗[Tk⊕Sn-k]

Ci=Tk⊕Sn-k-i⊕R,i=1,2,…,n-k

Cn-k+1=R

Dn-k+1=R0δ

Ei=[Tk-1⊕Sn-k-1-i]⊗[T0α⊗S0]⊗Ig,i=0,1,2,…,n-k-1

Fi=Tk⊕Sn-k-i⊕U,i=1,2,…,n-k

Fn-k+1=U

G0=[Tk-1⊕Sn-k-1]⊗[T0α⊗S0]⊗L0δ

Gi=[Tk-1⊕Sn-k-1-i]⊗[T0α⊗S0]⊗Il,i=1,2,…,n-k-1

3.3 系统稳态概率向量

系统在长期运行之后逐步保持稳定可靠的运行状态,此时,系统各状态在不同转移中停留在各状态的可能性构成系统稳态概率向量π,与Ω相匹配,能够进一步分解为π=(π0v,π0id,π1v,π1j,…,π(n-k)v,π(n-k)j,π(n-k+1)v,π(n-k+1)j),基于稳态概率向量的有关特性推导出其符合以下方程组要求:

(7)

将式(7)展开:

对上述方程进行求解,可得到系统稳态概率向量π。

4 系统可靠性指标

4.1 系统可靠度函数

系统可靠度主要指的是系统在某时间段t内稳定可靠运行的可能性。在n中取k系统中,可靠度指的是系统在某时间段t内停留状态集H1∪H2∪H3∪H4的可能性,结合PH分布的概念,将系统停机状态H5∪H6看作是吸收态,由此推导出下述状态转移矩阵:

该矩阵主要反映了系统在稳定可靠地工作时各种不同状态之间所发生的转移情况,根据系统可靠度能够推导出可靠度函数为。

R(t)=(αk⊗βn-k⊗γ,0,0,0)exp(Q*t)e

(8)

系统故障密度函数为:

f(t)=(αk⊗βn-k⊗γ,0,0,0)exp(Q*t)(Q*)0

(9)

4.2 系统稳态故障频度

故障频度主要指的是在某时间段里,系统或者其内部部件发生异常或者无法运行的频度。结合n中取k系统运行规律,推导出单一部件与整个系统的稳态故障频度。

1)系统工作部件稳态故障频度。工作部件在运行过程中如果出现异常或者无法运行,备用部件会立即安全规范地进入系统并保持稳定可靠的运行状态,若不对发生异常的部件进行维修,第n-k个备用部件迅速进入系统并安全规范地运行后又出现了异常,此时,已无可用备件,系统只能停机处理。故关于工作部件的稳态故障频度,其稳态概率可描述为(π1v,π1j,…,π(n-k)v,π(n-k)j,π(n-k+1)v,π(n-k+1)j)。由此推导出其稳态故障频度为:

monline=π0v(T0⊗e)+π0id(T0⊗e)+

(10)

2)系统稳态故障频度。在系统中,系统停机的稳态概率为(π(n-k+1)v,π(n-k+1)j),其状态集为H5∪H6。由此可推导出系统稳态故障频度为:

msystem=π(n-k)v(T0⊗e)+π(n-k)j(T0⊗e)

(11)

4.3 系统其他可靠性指标

1)系统稳态可用度。系统稳态可用度主要指的是系统在保持正常运行模式下的稳态可能性,换言之,是指系统处于状态集H1∪H2∪H3∪H4的稳态可能性,关于系统稳态可用度A,可将其描述如下:

A=1-(π(n-k+1)veg+π(n-k+1)jel)

(12)

2)系统平均故障间隔时间。系统故障间隔时间指的是系统在保持运行的过程中,两次连续停机时的工作时长[18]。基于此含义,可将系统停机状态H5∪H6看作是吸收态,由于PH分布具有封闭性,系统故障间隔时间亦满足PH分布,通过(φ,Q*)进行描述,φ=(αk⊗βn-k⊗γ,0,0,0)。基于以上论述可推导出MTBF=-φQ*-1e。

3)维修台忙期稳态概率。在系统停留状态集H4∪H6的情况下,维修台持续检修异常部件。故可将维修台忙期稳态概率Pj视为系统状态集H4∪H6内各子状态稳态概率之和。

(13)

按照上述分析,分别推导出维修台在休假期、空闲期相对应的稳态概率:

(14)

Pid=π0ide

(15)

5 算例

算例主要由三部分构成:①结合模型问题,提出假设条件即n=6、k=4,修理工休假次数H=4,工作部件、备用部件寿命均满足相应的PH分布要求,基于此构建系统可靠性模型,以此完成对模型适用性广的检验;②将n值作为调整参数,探讨n值增大对系统可靠性产生的影响;③假设修理工维修时间、休假时间的PH分别为[δ=1U=-λ]、[γ=1R=-φ],研究维修速率λ、休假速率γ这两大参数和系统可靠性之间的关系。

5.1 模型适用性验证

假设某n中取k系统,n=6、k=4,修理工计划休假次数H=4,工作部件、备用部件寿命均满足相应的PH分布要求,详情如下:

工作部件寿命分布:

备用部件寿命分布:

维修时间分布:

休假时间分布:

空闲时间分布:

H=4

θ=[0.43,0.23,0.34]

根据已掌握的各种信息及提出的模型假设,计算以上模型,基于此推导出系统稳态概率向量π,进一步确定下述各个重要可靠性指标:A=0.875 9,monline=0.461 1,msystem=0.261 5,Pj=0.493 8,Pv=0.369 6,Pid=0.136 6,MTBF=0.818 8单位时间。

结合可靠性相关定义,模型所得数据均符合可靠性要求,表明该模型可以较好地适用于复杂分布输入的系统可靠性建模问题。

将上述已知条件和模型假设中的工作部件、备用部件寿命以及维修时间的PH分布表示代入文献[19]中,对该文献中不考虑修理工休假的n中取k系统可靠性模型进行求解得出相关可靠性指标数据为:A=0.927 5,monline=0.351 2,msystem=0.213 7,Pj=0.232 6,MTBF=0.957 2单位时间。

通过将两者模型所得可靠性指标数据进行对比,可得出雷达图如图3所示。由对比雷达图可以看出,本文模型的一系列可靠性指标劣于不考虑修理工休假的文献[19],这表明修理工的调度与分配策略可能会造成系统部件故障后无法及时得到维修,进而使系统稳态可用度小和平均故障间隔时间较短,同时无法及时维修会造成系统稳态故障频度较高,导致维修台忙期稳态概率有一定的提高。两者可靠性指标的对比分析验证了该模型的实用价值,对于维修调度管理者而言,为弥补维修调度分配对系统可靠性的影响,须合理、科学安排维修计划,加大人员培训,提高维修效率,减少装备故障后待修时间,从而有效提高装备使用时限。

图3 与文献[19]模型所得可靠性指标对比雷达图Fig.3 Radar chart of system reliability features contrast to the model of reference[19]

5.2 n值变化对系统可靠性的影响

在该模型中以n值为调整参数,基于5.1节模型假设条件,分析n值单一变量增加对系统可靠性的影响,演示n取何值时系统可靠性最佳。具体分析如下:

令k=4,4≤n≤50,H=4,分析n值逐渐增大时系统可靠性指标变化趋势,如图4～7所示。

图4 系统稳态可用度与n的关系Fig.4 Relationship of n and system stationary availability

图5 系统稳态故障频度与n的关系Fig.5 Relationship of n and system steady-state failure frequency

图6 维修台忙期稳态概率与n的关系Fig.6 Relationship of n and busy steady-state probability of repair facility

图7 MTBF与n的关系Fig.7 Relationship of n and MTBF

由图4～7可以看出,随着系统部件的增加,系统可靠性能够有所增强,但是维修台忙期稳态概率则出现了不同的变化情况,主要表现为一开始下降,而后升高,最后逐步保持平缓状态,而这也反映了部件数目的增多,维修台忙期稳态概率不会持续下滑,如果n保持在30～35范围内,可获得最高值,其他可靠性指标也保持良好状态。所以,在修理工可变休假的n中取k系统中,应从可靠性、经济性等多方面入手进行全面细致的分析,以此确定最优部件数目。

5.3 修理工休假速率和维修速率对系统可靠性的影响

采用5.1节中系统工作部件、备用部件寿命所服从PH分布的相关数据,假设n=6、k=4,修理工计划休假次数H=4,修理工维修时间和休假时间的PH表示分别为[δ=1U=-λ]、[γ=1R=-φ]。结合系统可靠性相关定义,修理工休假速率和维修速率分别是修理工维修时间和休假时间的倒数,由此能够明确维修速率λ、休假速率γ这两大重要参数和系统运行性能间的关系,具体变化趋势如图8～10所示。

图8 系统稳态可用度与维修速率、休假速率的关系Fig.8 Relationship of maintenance rate, vacation rate and system stationary availability

图9 系统稳态故障频度与维修速率、休假速率的关系Fig.9 Relationship of maintenance rate, vacation rate and system steady-state failure frequency

图10 MTBF与维修速率、休假速率的关系Fig.10 Relationship of maintenance rate, vacation rate and MTBF

根据图8～10可知,若维修速率在2～3范围内,休假速率在4～5范围内,无论是系统稳态可用度,还是系统稳态故障频度等,均能够保持较高水平。由此可知,在系统运行过程中,结合实际情况,设定维修速率等相关参数,不仅有助于提高维修资源的利用水平,还能够减少维修费用,最重要的是能够增强系统可靠性,具有一定的现实意义。

6 结论

关于大型复杂装备维修力量管理调度这一重要现实问题,本文在借鉴前人研究成果的基础上,尝试以修理工可变休假策略为切入点对上述问题进行准确合理的描述。通过连续PH分布的相关特性,构建了n中取k系统可靠性模型,由此获得了包括系统稳态可用度等在内的一系列重要可靠性指标,借助算例分析对模型的适用性、通用性进行了检验,同时,深入分析了修理工休假、部件数目等相关参数的变动与系统运行性能的关系。针对大型复杂装备n中取k系统,本文认为维修调度管理人员在开展检修工作时,应分别从经济、安全角度等入手深入细致地分析修理工计划休假次数的调整对系统运行可靠性产生的具体影响,基于此,科学有效地部署维修力量,以期在增强系统运行可靠性的同时,促进维修效率大幅提升。

致谢

海军工程大学管理工程与装备经济系胡涛教授、陈童副教授在模型构建思路方面提供了帮助和指导,谨致谢意!