杨皓洁, 吕建伟, 徐一帆

(海军工程大学管理工程与装备经济系, 湖北 武汉 430033)

0 引 言

具有复杂冗余结构、时间相关故障的多状态系统是当前可靠性领域的热点研究对象之一。复杂冗余结构中的元件常表现出动态故障行为,相应元件在启用前后具有不同的失效率;而时间相关故障意味着故障率具有时变的特点,状态转移时间服从非指数分布。在该条件下故障率是时间的函数,因此故障与时间相关。

考虑这些因素会使系统故障演化规律和任务过程更为复杂,增大了系统可靠性及任务成功性评估的难度。系统可靠性评估的研究方法大致可分为解析法和仿真法。其中,解析法有组合模型法[1-2]、状态空间法[3-4]、概率近似法[5]等,具有计算精度高、计算结果理想可靠等特点。组合模型法中,以故障树分析为基础的多值决策图[6-8](布尔扩展模型法)、序列二元决策图等[9-10]是解决多状态系统和动态系统可靠性的较好方法,但缺点是单独使用此类方法时,不能处理可修系统,对于系统中的可修部分必须结合状态空间法。另一类主流方法是通用生成函数[11-12],该方法通过引入z变换函数、通用生成算子等工具,再结合系统元件的状态概率分布,计算得到系统的可靠性指标。通用生成函数是解决多状态系统可靠性问题的高效方法,但当系统动态规律较为复杂时,仍存在建模和求解困难的问题。当面对时间相关故障问题时,需要利用半马尔可夫模型进行建模,但该模型在建模和计算求解上都较为复杂。此外,解析法在面对时间相关故障问题时同样存在建模和求解复杂的问题。仿真法作为一类评估复杂系统可靠性及任务成功性的统计类方法[13-15],主要包括基于离散事件仿真[16]、贝叶斯理论[17-18]、petri网[19]等建模方法。其优点是易于建模且模型通用性较强。但由于仿真抽样效率在针对具体问题时各不相同,计算精度和仿真次数之间的矛盾难以兼顾。目前主要通过特定的方差减小方法[20-23]来提高计算结果的精度,但相关研究中均以元件状态转移时间服从指数分布为前提，该条件下不需考虑状态转移时间的剩余分布问题,并且系统状态转移时间也服从指数分布,容易实现系统状态转移时间的直接抽样,在简化问题的同时限制了方法的使用范围。一旦考虑元件状态转移时间不服从指数分布,系统状态转移时间的概率分布将变得相当复杂,导致在对其进行抽样时,相关计算较为困难。

为了在考虑时间相关故障、多状态性能输出、系统冗余储备等因素的同时,改进抽样效率和估值精度,本文提出了一种基于蒙特卡罗离散事件仿真,冗余多状态系统可靠性及任务成功性的评估方法,用以针对故障率具有时变特点的时间相关故障。通过设计不服从指数分布状态转移随机抽样的数值算法,实现故障事件的随机抽样。在此基础上,结合受迫转移、失效偏倚等方差减小方法[24-27],改善小概率事件抽样效率和估值精度,为该问题提出一种较为完善的解决方案。

1 问题描述及分析

某多状态系统S由n台不同功能的设备组成,其中设备j有mj种性能状态。任意设备的状态变化都会引起系统状态转移,因此在任务过程中的某一时刻t,系统状态S(t)及设备性能状态Gj(t)的关系如下:

S(t)={G1(t),G2(t),…,Gn(t)}

(1)

(2)

(3)

式中:s为系统状态空间;gj为设备的性能状态空间。实际上,系统状态就是由各设备状态组成的向量,但不同的系统状态可能有相同的性能输出表现,因此定义系统性能状态G(t)为

G(t)=Φ(S(t))∈g={g1,g2,…,gm}

(4)

式中:Φ(·)为某具体系统的结构和运行方式确定的构型函数;g为系统的性能状态空间。设备的性能状态转移关系如图1所示。

图1 设备性能状态转移图

图1中各类状态及对应编号如下:完好为n、完好到故障间的中间状态为(n-1)到2、一般故障为1、致命故障为0。其中完好对应性能输出为100%。中间状态对应性能输出为0%～100%。一般故障、致命故障对应性能输出均为0%,其区别在于一般故障为可修复故障,致命故障为无法修复或任务时间内不具有维修条件的故障,也标为吸收态。λi,i′和μi,i′分别为该设备从状态i′→i的故障率和修复率。

为达到提高系统可靠性的目的,系统中存在一些冗余结构,包括冷、温、热储备。由冗余结构和系统运行产生设备的工作状态大致可分为3类:工作、储备待机和维修。设备工作状态根据系统结构、备用设备启动规则和设备当前状态来确定。假设x、y分别为一个冗余结构中的主设备和备用设备,图2展示了冷、温储备结构下的主设备与备用设备工作状态的转移关系。热储备结构较为简单(即均视为工作),在此不再赘述。

图2 冷、温储备设备工作状态转移规则图

图2中Gx(t)和Gy(t)表示设备x和y在t时刻的性能状态,数值为图1中性能状态的编号;Sysstate表示系统状态;down表示系统宕机,available表示系统可用。当主设备性能状态低于备用设备或发生一般故障时应当转入维修,同时备用设备转入工作状态;主设备维修完成后立即转入工作状态,同时备用设备转入待机状态;当系统中其他部分性能退化导致系统宕机时,主、备单元均转入待机状态。根据此规则并结合可靠性框图即可确定所有设备的工作状态。在动态系统中，某设备的状态由工作状态和性能状态共同定义。其意义在于确定系统中所有的状态转移速率,也是蒙特卡罗随机事件抽样的前提。

2 基于蒙特卡罗仿真的模型描述

2.1 仿真流程设计

本文所采用的蒙特卡罗离散事件仿真方法流程如图3所示。在单次任务仿真过程中,需进行多次离散事件抽样循环。一次完整的离散事件抽样步骤如下:

步骤 1在系统最底层(设备级)进行系统状态转移随机抽样,得到故障或维修事件发生的时间和转入状态,这是仿真的核心步骤,如图3中灰色部分所示;

图3 总体仿真流程图

步骤 2以底层设备状态变化驱动上层系统状态变化;

步骤 3根据系统状态变化更新任务相关变量;

步骤 4根据任务相关变量和系统状态判断任务是否失败。

在完成仿真后整理数据得到统计结果。设任务成功概率为D,系统瞬时可用度为A(t),则有

D=(N-Nd)/N

(5)

A(t)=[N-Nf(t)]/N

(6)

式中:Nd为任务失败频数;Nf(t)为t时刻系统处于不可用状态频数。Nd和Nf(t)的计数与系统状态转移的权重wt有关。wt是一个与时间相关的变量,其含义为t时刻前所有随机事件对之后事件发生的贡献程度。因此,wt代表t时刻前每一次状态转移累积的权重,同时代表t时刻某随机事件的仿真计数值。任务开始时令w0=1。具体wt的计算与系统状态转移随机抽样在第2.2节、第2.3节一并展开。当任务失败时更新任务失败频数:

Nd←Nd+wT0

(7)

当t时刻系统处于不可用状态时更新系统不可用状态频数:

Nf(t)←Nf(t)+wt

(8)

不同的抽样方法所产生的随机事件的权重可能是不同的,wt可能不为1,因此Nd和Nf(t)也可能不是整数。

2.2 系统状态转移随机抽样

蒙特卡罗方法通过生成随机事件,模拟系统的动态行为。本研究使用的基本抽样方法是蒙特卡罗间接仿真法,在此基础上才能够使用方差减小方法来提高仿真效率及减小统计误差。蒙特卡罗间接仿真法的基本步骤是:先随机生成系统当前状态的驻留时间(或状态转移的发生时间),再通过轮盘赌等方式随机确定触发事件对应的设备和该设备转移时进入的状态,如图4所示。

图4 多状态系统状态转移的蒙特卡罗间接仿真法抽样过程

对于多状态系统而言,分系统、设备存在多个性能状态,各状态之间的转移概率服从特定概率分布。系统在某一给定时刻发生状态转移并进入新的状态是由一个状态概率转移核(下文简称为“转移核”)决定的,转移核可以全面地描述系统在时域内的随机行为。转移核本质上是一个概率密度函数。多状态系统的转移核可以表示为

K((k,t)|(k′,t′))=f(t|k′,t′)q(k|k′,t)

(9)

这里用f(t|k′,t′)表示在时刻t′发生下一次状态转移时刻为t的概率密度函数,用q(k|k′,t)表示在时刻t时状态k′→k的转移概率,且有

(10)

(11)

(12)

式中：γk′(t)表示当前处于时刻t时系统转移出当前状态k′的转移速率;γkk′(t)表示当前处于时刻t时状态k′→k的转移速率;k′和k均属于系统状态空间s。系统经过多次状态转移形成了一条状态的随机游动序列ΓM:(k0,t0),(k1,t1),…,(kM,tM),以实现离散事件抽样的循环。

在转移核的基础上就可以对转移发生时刻t和转移进入状态k进行抽样,蒙特卡罗间接仿真在正常状态转移空间的概率分布下进行抽样,因此系统状态转移权重为

wt=wt′

(13)

(1)转移发生时刻t的抽样

发生状态转移时刻t计算公式为

(14)

式中：ξ为0～1之间均匀分布的随机数。为了描述的简洁,这里将故障、降级、退化等系统状态变化统称为故障类随机事件,将维护、修理等系统状态变化统称为修复性随机事件。根据常规标识习惯,当转移后系统发生故障或退化时,转移速率用λ表示;当转移后系统修复时,转移速率用μ表示;Ok′表示状态k′的可达状态集;用Λk′表示状态k′可达的故障或降级状态集,Ψk′表示可达的修复状态集:

(15)

Ok′=Λk′∪Ψk′

(16)

由于本文假设设备当前状态的驻留时间均服从非指数分布,因此系统状态转移的抽样带有一定复杂性,具体内容在下文中展开。

(2)转移进入状态k的抽样

用轮盘赌的方式可以确定转移到达的具体状态:

(17)

式中：ζ为0～1之间均匀分布的随机数。当ζ落在式(17)区间内,则系统从状态k′转移至状态k,也就是确定了发生状态转移的设备和该设备转入的状态。蒙特卡罗间接仿真法利用两个随机数就可以产生一次状态转移,而且不需要比较等步骤,一定程度上提高了计算效率,在此基础上可使用方差减小方法进行加速抽样。

2.3 方差减小方法

(18)

(19)

(1)受迫转移

(20)

(21)

可见受迫转移法只改变转移时间的抽样空间,不改变转移状态的抽样空间。然后按照蒙特卡罗间接仿真法进行转移发生时刻t的抽样,此时式(14)变为

(22)

受迫转移法使状态转移时刻发生在任务结束时间内,客观上起到了强制提升随机事件发生概率的作用,从而增加了随机事件的抽样实现次数。为使估计值保持无偏性,修正仿真抽样计数：

(23)

(2)失效偏倚

(24)

(25)

式中:x为失效偏倚系数,用于调节倾向故障类转移及修复类转移之间的相对程度。可见失效偏倚法只改变转移状态抽样空间,不改变转移时间抽样空间。x的具体含义是本次抽样中故障类转移占所有状态转移的比例;(1-x)是本次抽样中维修类转移占所有状态转移的比例,通常取值范围在0.5～0.7[21]。这时系统转移到故障状态的转移速率将远大于不偏倚抽样时的转移速率。然后根据式(17)确定转移进入的状态k。为了使估计值保持无偏性应当修正仿真抽样计数:

wt=

(26)

(3)抽样方法的使用条件

每1次随机事件抽样只能从蒙特卡罗间接仿真、受迫转移、失效偏倚3种方法中选用一种,根据每种方法的特点可在仿真中做以下规定:① 每次任务初始化后第1次状态转移(k0,t0)→(k1,t1)使用受迫转移法;② 当系统中存在维修类转移时使用失效偏倚法;③ 其余情况使用蒙特卡罗间接仿真法。

3 非指数分布随机抽样的数值算法

由于系统中存在服从非指数分布的状态转移,其转移率为时间的函数,因此系统当前状态的驻留时间无“记忆性”,是典型的非齐次马尔可夫过程。在离散事件仿真中,随着每次仿真时间的推进,一直正常工作或一直处于维修状态的设备的当前状态驻留时间应按其剩余分布进行抽样。以某设备的故障转移为例,假设其故障时间T的初始分布为F(t),则在累计工作tc后,其剩余分布函数[28]为

(27)

作为减小方差的基础方法,蒙特卡罗间接仿真法首先需要对系统当前状态驻留时间进行抽样,则考虑剩余分布后的系统状态驻留时间所服从的分布函数可表示为

FS-RD(t)=P{min(T1,T2,…,TK)

(28)

步骤 1选取一个足够大的数ts,使得FS-RD(ts)=1。尽管FS-RD(ts)在理论上只能无限接近于1,但在计算机上可以实现等于1。

步骤 2在[0,tu]上均匀选择若干点组成向量Ts,并根据式(28)计算相应时间的概率值Ps=FS-RD(Ts)。由此组成离散样本点组合(Ts,Ps)=[(t0,p0),(t1,p1), …,(ts,ps)]。

步骤 3对组合(Ts,Ps)进行补充和删选。确定向量Ts中使得FS(tu)>0.99的最小元素tu,并在[0,tu]上选择若干点以步骤2相同的方式补充离散样本点组合(Ts,Ps)。Ps中可能存在重复元素,只保留Ps中不重复元素和Ts中对应的元素,删除其余元素,为插值法做准备。

步骤 4利用插值法进行随机抽样,以概率Ps为自变量,时间Ts为因变量。生成0～1之间均匀分布的随机数ξ,若此次抽样方法为受迫转移法,则生成一个[0,FS-RD(T0-t′)]间的随机数ξ。在组合(Ts,Ps)中选取自变量与ξ相邻的两个点(ti,pi)、(ti+1,pi+1),再进行线性插值计算得到系统当前状态驻留时间Δt：

(29)

该算法不需大量样本点就可以完整地展示分布函数全貌,在有限的计算资源下提高了插值法精度,还可以降低计算量。并且状态转移时间分布函数FS-RD(t)是单调递增函数,能保证上述算法Δt有唯一值。

以威布尔分布为例,在仿真开始时令累计工作时间tc=0,正常工作及热储备设备的tc应累加仿真推进时间,发生故障需维修的设备的tc清零,同时更新该设备相应的状态转移率:

γ(tc)=bη-b(tc)b-1

(30)

式中：b为形状参数;η为尺度参数。根据式(27),得到该失效机理发生时间剩余分布为

(31)

系统中存在相互独立、相互竞争的多种服从威布尔分布的失效机理和维修转移,bkk′和ηkk′为相关转移的形状参数和尺度参数。根据式(28)可得系统的当前状态驻留时间分布函数。

(32)

图5展示了该算法的抽样结果。在随机抽样生成系统驻留时间Δt后,根据式(17)进行轮盘赌确定最终发生状态转移的设备和具体转移状态。

图5 非指数分布随机抽样数值算法结果

4 算例分析

4.1 算例背景

本算例以舰船及其装备系统为对象进行分析,由于篇幅所限,这里仅以舰船动力、电力系统在航渡任务阶段的任务可靠性和任务成功率评估为例,其他系统及任务的相应计算也同样适用。设某舰动力、电力系统的系统结构组成如图6所示。相关设备的性能状态转移关系如图7所示。属于第I类设备的有柴油机、燃气轮机,这类设备有多种性能输出表现,且可能发生致命故障;属于第II类设备的有柴油发电机组、轴系、螺旋桨,这类设备只有完好和失效两种性能输出表现,且可能发生致命故障;属于第Ⅲ类设备的有监控设备、辅助设备。在本文应用背景中,假设这类设备不具有影响任务成败的致命故障,可将这类设备简化为可修的二态设备。系统中的冗余结构包括冷储备结构(柴油机1,柴油机2和燃气轮机)和热储备结构(轴系和螺旋桨1,螺旋桨2),其中冷储备结构按照图2所示规则启用备用设备。

图6 航渡任务中舰船系统结构

图7 相关设备性能状态转移图

该舰在航渡阶段的主要任务目标[29]为:① 按时抵达目标海域;② 当航渡任务结束时系统仍处于可工作(非故障)状态。目标①是对完成任务时间和空间的要求,指在规定的任务时间T0里完成规定的任务航程S0。舰船既不可提前到达,也不可晚于规定时间到达。提前到达可能导致本舰脱离协同兵力、提前暴露行踪;晚于预定时间到达则难以实现兵力协同、完成既定作战任务。目标②是对系统下一阶段完成作战任务的要求。总之,当且仅当以上两项目标均满足时,方可判定任务成功(实际情况下的任务成功性问题判据可能更为复杂,此处略去)。据此可知,任务目标是否达成应结合任务成功性判定规则,根据特定的依据给予评判,如表1所示。在表1中,航速超限是指由于经历系统降级或故障维修,导致该舰在剩余航程中即使以最大航速也不能完成任务。维修超时是指任务结束时由于设备故障正在维修,导致该舰总体处于宕机、不可用的状态。当发生维修超时时,舰船通常还未到达目标区域。致命故障是指由于设备致命故障导致该舰总体宕机且无法修复,从而无法完成剩余任务。

表1 航渡任务成功判据

4.2 模型假设条件

根据上文分析对该案例相关分系统做进一步假设。

假设 1在航渡阶段,其他分系统不影响该舰总体状态。

假设 2分系统和舰船的可用状态包括完好和降级状态,系统降级时会影响该舰航速。

假设 3所有状态转移均为独立状态转移。

假设 4设备维修为完全维修(修旧如新),修理后状态均返回完好状态,且维修完成后故障率回归初始状态,即λ(t)=λ(0)。

假设 5不考虑维修资源限制。

4.3 任务结果判定

表1中3种任务失败类型可通过以下方法进行判定。

(1)航速超限

设S0为任务航程;St为当前累计航程;T0为任务总时间;t为当前时刻;ϑ为航程松弛系数[30](0<ϑ<1);vmax为最大可持续航速;vp为剩余航程计划航速;v为实际航速。根据松弛系数定义确定任务航程S0=ϑvmaxT0。通过上次随机事件确定的实际航速v与此次随机事件对应的时间间隔Δt(抽样得到的状态驻留时间)的乘积来更新累计航程St。

t=t′+Δt

(33)

St=St′+vΔt

(34)

(35)

结合实际情况,舰船通常具有特定的经济巡航航速。航速与动力分系统的功率输出情况和航速规则有关,当剩余航程计划航速vp低于veco时,则以veco航速航行,此次随机事件发生后的实际航速v为

(36)

式中：θ为功率输出系数(0<θ<1),表示舰船在某状态下的最大功率占完好状态下满功率的比值,具体取值根据动力系统状态确定。根据vp与vmax的关系判定是否航速超限。需要指出的是,在仿真中当动力分系统功率输出降级时,可能暂时出现vmax>vp>v的情况,此时暂不判定失败。当修理完成后,功率输出恢复完好状态,使得vpvmax,则认定为航速超限,判定任务失败。

(2)维修超时

若本次随机事件的发生时间超出任务时间,即t>T0,并且此时系统处于宕机维修状态,则认定为维修超时,判定任务失败。

(3)致命故障

若某设备发生致命故障,导致系统进入吸收态,即系统无当场修复的可能,则认定为致命故障,判定任务失败。

由于存在多种任务失败类型,因此可扩展Nd为向量，Nd=[NOS,NUM,NFF],其中元素分别对应航速超限、维修超时和致命故障的频数。按照下式计算不同失败类型概率,以便于观察其对任务成功性的影响。

Dd=Nd/N

(37)

4.4 算例数据设置

(1)系统参数

通常机械设备磨损、电子元件老化等退化失效模式是造成设备一般故障和降级退化的主要原因,由此可认为其失效、退化时间服从形状参数大于1的威布尔分布。而发生致命故障的随机性较大,故此算例将其失效时间设为指数分布(即形状参数等于1的威布尔分布)。设备的维修时间通常为正态分布,但正态分布下可能得到负数的抽样值,这与维修时间非负性的要求不符。因此,通常用威布尔分布或对数正态分布来近似模拟正态分布。本文假设维修时间服从威布尔分布。当威布尔分布的形状参数大于3时,形态逼近于正态分布。例如某随机变量X1～N(50,100)可用X2～Wbl(54, 5.8)代替,其概率密度函数如图8所示。

图8 正态分布与威布尔分布概率密度函数对比

系统中各设备状态转移时间所服从威布尔分布的尺度参数η和形状参数b如表2所示。当系统处于可用状态时,式(36)中功率输出系数θ主要由柴油机、燃气轮机状态及功能结构决定,设备性能状态与功率输出系数关系如表3所示。可见系统也会表现出介于完好与故障间的降级退化状态,这时最大航速只有完好情况下的60%。

表2 设备状态转移分布的参数设置

表3 动力系统设备性能状态与功率输出系数关系

(2)任务及仿真参数

设经济航速veco=18 kn,满功率最大航速vmax=25 kn,航程松弛系数ϑ=0.8,任务时间T0=50 h,仿真次数N=105,失效偏倚系数x=0.5。

4.5 仿真结果分析

(1)仿真模型有效性检验

本文通过解析方法对照验证仿真模型的有效性。由于本文研究的问题具有一定的复杂性,所以解析法建模和计算的成本会很高,因此对条件进行一定的简化,使得解析法过程不过于复杂。简化假设条件并不影响仿真模型的有效性。一是令状态转移时间服从指数分布,即令表2中所有形状参数b=1;二是令系统中储备方式均为热储备，其余参数均不改变。

解析法只能求得系统的可靠性指标(可用度),所以只能通过检验仿真法中系统可用度是否与解析法一致的方式验证其有效性。由于在检验仿真法中,系统可用度和任务成功概率为同时输出,且任务成功性实则是系统可靠性在任务方面的表现。因此也能够验证任务成功性指标的有效性。

本文采用马尔可夫方程与通用生成函数相结合的解析方法求解系统可用度。在图9展示的结果中,仿真法下可用度曲线与解析法下基本一致,从而验证了仿真方法的有效性。

图9 仿真法与解析法下的系统可用度对比

(2)抽样方法效率对比

表4和图10展示了仿真次数均为105次时,使用方差减小方法(受迫转移&失效偏倚)和仅使用蒙特卡罗间接仿真法的任务成功概率仿真结果。在系统存在冗余设计和任务时间较短的高可靠性系统背景下,方差减小方法的优势显著。该方法在一次仿真中能够抽得更多的故障样本，在总体仿真中能够抽得更多的任务失败的样本,因此任务成功性指标收敛速度明显更快,如图10所示。方差减小方法所得结果的精确小数位数不受仿真次数的限制,105次仿真结果的精确小数位数不止5位。

表4 不同抽样方法下任务成功性仿真结果(T0=50 h)

图10 任务成功性指标收敛速度对比

图11比较了相同仿真次数N=105情况下蒙特卡罗间接仿真法与方差减小方法的系统瞬时可用度,并以样本方差展现了瞬时可用度估值的波动程度。与蒙特卡罗间接仿真法相比,图11中方差减小方法得到的可用度曲线波动幅度和方差范围减小的效果显著。两种方法得到的仿真结果趋势基本吻合,而方差减小方法的结果精度更高。这说明方差减小方法能够有效结合本文所设计的随机抽样数值算法,达到了非齐次马尔可夫系统加速仿真进程和减小方差的目的。

图11 系统瞬时可用度(T0=50 h)

(3)任务过程的仿真结果

图12和图13选取了两次较为典型的仿真过程,能够体现系统多状态特性。并且能够明显地展示出系统状态、系统的任务表现与任务目标之间的实时关系。计划航速表示目标完成任务所需要达到的航速。系统功率输出系数表示当前舰船状态下的性能输出能力,其代表的舰船状态将直接影响实际航速,从而影响该舰能否完成任务目标。需要指出的是,功率输出系数表明系统输出的能力高低,该系数为1时,舰船可通过动力及传动装置在经济航速至最大航速之间的航速航行。实际航速将依据计划航速的需要而定,其规则如式(36)所示。在图12所示的任务过程中,该舰先后经历了2次宕机和1次降级退化,即便在25 h时已修复为完好状态,但是此时已无法满足vp>vmax条件,无法按时到达指定海域,因此认定为航速超限,任务失败。图13所示的过程与图12类似,但宕机、降级所延误的时间相对较短,且此后一直保持完好的性能状态,最终任务成功。

图12 任务失败仿真过程

图13 任务成功仿真过程

为了比较计划航速、实际航速和系统功率输出系数变化趋势之间的关系,图14展示了上述指标值的样本均值曲线。需要强调的是,图14中各指标值的每个数据点并非对应相应时刻单个样本的实际数据,而是统计意义下的平均水平。尽管数据的变化幅度较小,但仍能展现各指标变化趋势之间的相互关系。随着任务推进、系统性能不断退化,系统宕机、降级概率逐渐增加,平均功率输出系数呈下降趋势。因此,需要更高的计划航速才能完成航渡任务,表现为平均计划航速随着任务时间的推进而上升。实际航速在满足功率输出条件的情况下,应不低于计划航速才能完成航渡任务。因此,实际航速随着计划航速的增长而提升。但幅度低于计划航速的增幅,这是因为系统状态的劣化导致任务后期实际航速可能无法达到计划航速的需求。

图14 实际航速、计划航速、功率输出系数的均值对比

(4)状态转移分布的影响

按照第4.4节中的任务参数设置,对非指数分布与指数分布状态转移下的仿真结果进行对比。为保证仿真具有可比性,应使不同状态转移分布具有相同期望值。根据表2中威布尔分布的各项参数,计算各类转移时间的期望值作为指数分布条件下的平均故障、维修时间,可表示为

E(T)=ηΓ(1+1/b)

(39)

式中:Γ(·)为伽马函数;E(T)为威布尔分布的期望值。任务成功性与可靠性结果如表5和图15(a)所示,可见指数分布下的任务成功概率和可用度相对较低,这是由于在系统运行之初指数分布的失效率更高,这种现象在任务时间较短的情况下较为明显。说明若按照指数分布进行建模可能会低估系统可靠性及任务成功性。在本算例参数设置下,当任务时间较长时,系统中可修部分对可用度的影响已趋于稳定。可用度的变化主要受不可修部分(致命故障)影响,由于不同状态转移分布下致命故障均视为指数分布,因此可用度的变化趋势基本相同,但在相同时刻指数分布下的估值相对较低,如图15(b)所示。

表5 不同状态转移概率分布下任务成功性仿真结果(T0=50 h)

图15 不同状态转移概率分布下的系统瞬时可用度

5 结 论

本文研究了考虑时间相关状态转移的系统任务可靠性和任务成功性问题。针对非齐次马尔可夫过程在离散事件仿真中的剩余分布抽样,设计了服从非指数分布状态转移随机抽样的数值算法,并结合受迫转移、失效偏倚等方差减小方法,改善了抽样效率和估值精度。以舰船航渡任务为算例,评估了其任务可靠性和任务成功率,结果表明本文建立的仿真模型可以反映系统状态、任务目标和系统性能表现间的变化关系,验证了本文方法对评估复杂系统行为和故障机制条件下系统可靠性及任务成功性的合理性与有效性。本文方法不仅适用于马尔可夫系统和具有解析形式状态转移分布的非齐次马尔可夫系统,还可适用于其他复杂系统行为及任务剖面组合的情况,具有一定的通用性。对于竞争性故障、级联故障、共因失效等更为复杂的相关故障机制,以及等待维修状态、维修资源可用性等维修方面问题,将在进一步研究中加以考虑。