一种两阶段衰退与冲击竞争风险的综合预防维修方法

摘要

本发明涉及一种两阶段衰退与冲击竞争风险的综合预防维修方法，包括：同时考虑稳态与非稳态两阶段衰退与冲击竞争风险，建立产品失效过程模型；采用步骤S1建立的失效过程模型进行仿真，得到样本在稳态阶段的失效时间或稳态阶段的区间长度，以及样本在非稳态阶段的失效时间：根据步骤S2获得的数据，得到产品失效时间累积分布函数；利用步骤S3得到的产品失效时间累积分布函数建立预防维修策略并优化。本发明实施例方法全面考虑了两阶段衰退过程，对有效开展产品预防维修具有很强的指导意义。

说明书

技术领域

本发明涉及产品可靠性工程技术领域，尤其涉及一种两阶段衰退与冲击竞争风险的综合预防维修方法。

背景技术

各类产品在运行或使用过程中可能经历多种相互竞争的失效过程，任何一种失效过程满足一定条件时均可能导致产品失效。两类主要的产品失效过程是衰退过程与随机冲击。各类产品经常因为磨损、疲劳、腐蚀、老化等因素逐渐衰退或失去既定的功能。另一方面，冲击、过载或其他外部应力也可能导致产品停止工作。这些失效过程可能是独立的，也可能是相互关联的。目前关于多失效过程的研究成果大多假设产品自投入运行后衰退的程度就不断加大，很少考虑多阶段的衰退过程。例如，据观察很多产品的衰退过程就存在稳态与非稳态两个阶段。在稳态阶段，产品的衰退程度随时间没有明显的增加，而在非稳态阶段产品的衰退程度不断增加直至失效。最近，Yang等(Li Yang,Xiaobing Ma,Rui Peng,Qingqing Zhai,Yu Zhao.A preventive maintenance policy based on dependent two-stage deterioration and external shocks.Reliability Engineering and SystemSafety,2017,160(4),201-211)提出了一个两阶段衰退与外部冲击竞争风险的预防维修策略，考虑产品存在稳态与非稳态两阶段衰退过程，且两阶段中产品所有外部冲击的概率有所不同。然而，(Li Yang,Xiaobing Ma,Rui Peng,Qingqing Zhai,Yu Zhao.A preventivemaintenance policy based on dependent two-stage deterioration and externalshocks.Reliability Engineering and System Safety,2017,160(4),201-211)中假设产品所受的外部冲击均将导致产品的失效，并制定了基于时间与状态观测的预防维修策略，一旦观察到产品进入非稳态阶段则开展预防维修。实际中，很多产品虽然在恶劣工况中遭受外部冲击，但未必会立即失效，而是导致衰退程度的恶化或衰退速率的加快；此外，很多产品进入非稳态阶段，很可能仍能开展有效的工作，可以考虑让其继续工作一段时间再开展预防维修。可见现有技术中关于两阶段衰退中预防维修的方法具有局限性。

发明内容

鉴于上述的分析，本发明旨在提供一种两阶段衰退与冲击竞争风险的综合预防维修方法，在产品衰退进入非稳态阶段后仍允许其工作一段时间，以提高产品的利用率。

本发明的目的主要是通过以下技术方案实现的：

在基于本发明方法的一个实施例中，提供了一种两阶段衰退与冲击竞争风险的综合预防维修方法，包括：

S1、同时考虑稳态与非稳态两阶段衰退与冲击竞争风险，建立产品失效过程模型；

S2、采用步骤S1建立的失效过程模型进行仿真，得到样本在稳态阶段的失效时间或稳态阶段的区间长度，以及样本在非稳态阶段的失效时间：

S3、根据步骤S2获得的数据，得到产品失效时间累积分布函数；

S4、利用步骤S3得到的产品失效时间累积分布函数建立预防维修策略并优化。

在基于本发明方法的另一个实施例中，步骤S1具体包括：

产品包括稳态阶段和非稳态阶段两个阶段，

当冲击力大小大于冲击失效阈值D_hf时为致命冲击，产品立即失效；

当冲击力大小小于冲击失效阈值D_hf时为非致命冲击，此时如果产品处于稳态阶段，则产品衰退变量值X_s(t)保持稳定，并在一常值附近上下波动；此时如果产品处于非稳态阶段，则产品衰退变量值X_s(t)随非致命冲击次数增加而增加，直到衰退值大于衰退损伤阈值时，产品失效。

在基于本发明方法的另一个实施例中，产品在非稳态阶段，当冲击力大小小于产品冲击失效阈值D_hf时，判断冲击力大小是否大于状态转移阈值D_e：

否，则冲击导致衰退程度增加；

是，则冲击导致衰退程度和衰退速率增加；

其中衰退程度和衰退速率增加导致衰退变量值增加。

在基于本发明方法的另一个实施例中，当产品在稳态阶段，

其中，为正态分布的截距参数，ε(t)为正态分布的过程噪声；

当产品在非稳态阶段，

其中，t_sw为产品衰退过程由稳态阶段向非稳态阶段的切换时刻，β_k为第k次状态转移之前，第k-1次状态转移之后的衰退速率参数；t_ck为产品进入非稳态阶段后发生第k次状态转移的时间；Y_i对应第i次冲击导致的衰退增加程度；ε(t)为正态分布的过程噪声；n(t)为产品衰退进入非稳态阶段后至时刻t发生的冲击数；nc(t)为产品衰退进入非稳态阶段后至时刻t发生的状态转移数，β_nc(t)+1是第nc(t)+1次状态转移之前，第nc(t)次状态转移之后的衰退速率；t_nc(t)为第nc(t)次状态转移发生的时间；

当产品发生第k次状态转移时，产品衰退速率由β_k变为β_k+1，定义：

β_k+1＝β_k+η_k(3)

其中η_k为正态分布的正随机变量；

在产品衰退稳态阶段各冲击的到达时间符合到达率为λ₁(t)的齐次非齐次Poisson分布，而在产品衰退非稳态阶段各冲击的到达时间符合到达率为λ₂(t)的齐次非齐次Poisson分布，且有λ₂(t)≥λ₁(t)；

W_i为一系列独立同分布的非负随机变量，且有相同的累积分布函数G(w)，G(w)＝P(W_i＜w)，其中G(w)与D_e、D_hf无关，w为一个代表冲击大小的随机变量；

根据冲击大小的累积分布函数，则对于第i次冲击产品不发生硬失效的概率为：

P(W_i＜D_hf)(4)

另一方面，Y_i与W_i相关，记为：

Y_i＝Q(W_i)(5)

其中，Q(W_i)为一递增的映射函数。

在基于本发明方法的另一个实施例中，步骤S2中得到样本在稳态阶段的失效时间或稳态阶段的区间长度具体包括：

S211、确定仿真样本数量N，时间轴长度T_s1，仿真步长Δt；

其中，时间轴长度T_s1使得全部样本在T_s1内能进入非稳态阶段；仿真步长△t满足在△t内不会发生两次冲击；

S212、选取样本n，初始值n＝1；仿真产生样本n的稳态阶段的区间长度t_swn；

S213、基于强度为λ₁(t)的Poisson过程产生样本n的冲击序列，对应时间为t₁、t₂、…t_M，大小为W₁、W₂、…W_M，其中t_M≤T_s1；

S214、定义冲击时间t_Z，其中t_z＜t_M，使得t_z≤t_swn≤t_z+1；

其中t_swn为样本n衰退过程由稳态阶段向非稳态阶段的切换时刻；

S215、判断Max(W₁、W₂、…W_z)是否大于等于D_hf，否，则IND_n＝0，样本序号n＝n+1，并执行步骤S218；是，则执行步骤S216；

其中IND_n用于指示样本n是否在稳态区间失效，如果是，则IND_n＝1；如果否，则IND_n＝0；

S216、确定第一个导致失效的冲击时刻t_j，1≤j≤z，该冲击的大小为W_j；

S217、记录样本n的失效时间t_f1＝t_j，置IND_n为1；样本序号n＝n+1；

S218、判断n是否大于N，否，则执行步骤S212；是，则结束仿真。

在基于本发明方法的另一个实施例中，步骤S2中得到样本在非稳态阶段的失效时间具体包括：

S221、确定仿真样本数量N，时间轴长度T_s2，仿真步长△t；

其中，时间轴长度T_s2使得全部样本在T_s2以内失效；仿真步长△t满足在△t内不会发生两次冲击；

S222、选取样本n，初始值n＝1；

S223、基于强度为λ₂(t)的Poisson过程产生样本n的冲击序列，对应时间为t₁、t₂、…t_M，大小为W₁、W₂、…W_M，其中t_M≤T_s2，冲击次数标识i，初始值i＝1；状态转移标识k，初始值k＝1；

S224、进入仿真，仿真产生ε(t)、W_i、β_k、

S225、判断t是否大于等于t_i，否，X_s(t)＝X_s(t-Δt)+β_kΔt+ε(t)，并执行S226；是，判断W_i是否大于等于D_hf；是，则执行步骤S228，否，执行S227，其中t_i代表第i次冲击的时间；

S226、判断Xs(t)是否大于等于D_sf；否，t＝t+Δt，返回步骤S224，是，则执行S228；

S227、判断W_i是否大于等于D_e；

是，则仿真产生η_k，i＝i+1、k＝k+1，β_k＝β_(k-1)+η_(k-1)；

X_s(t)＝X_s(t-Δt)+β_k(t-t_k)Δt+β_(k-1)(t_k-t+Δt)+Y_i+ε(t)并执行步骤S226；

否，则i＝i+1，Xs(t)＝Xs(t-Δt)+β_kΔt+Y_i+ε(t)，并执行步骤S226；

S228、记录样本n的失效时间t_f2＝t，样本序号n＝n+1；并判断n是否大于N；否，则执行步骤S222，是，则仿真结束。

在基于本发明方法的另一个实施例中，步骤S3具体包括：

根据步骤S2中获取的产品在稳态阶段的区间长度t_swn或失效时间t_f1，以及产品在非稳态阶段的失效时间t_f2，样本n的最终失效时间定义为：

根据TF_n以及数据经验累积分布函数估计方法估计产品在全周期的失效时间累积分布函数F₁(t)，或根据t_f2以及数据经验累积分布函数估计方法估计产品非稳态阶段的失效时间累积分布函数F₂(t)。

在基于本发明方法的另一个实施例中，步骤S4具体包括：

S41、根据产品失效时间的累积分布函数计算和

其中，是针对产品全寿命周期的最优预防维修时间，是产品在非稳态阶段的最优预防维修时间；

S42、产品开始工作，t＝0；

S43、每隔Δt_in观测一次产品衰退状态；

Δt_in为产品状态观测间隔，当一次状态观测的费用较高时，可先根据历史样本数据，采用交叉校验方法优化选取Δt_in；

S44、判断产品是否进入非稳态阶段，是，则在后开展预防维修；否，则执行步骤S45；

S45、判断当前时刻是否大于等于时刻，否，则执行步骤S43；是，则立即开展预防维修。

在基于本发明方法的另一个实施例中，步骤S41中计算和的具体方法为：

根据步骤S3中得到的产品在全周期的失效时间累积分布函数F₁(t)，或产品非稳态阶段的失效时间累积分布函数F₂(t)，分别对应计算产品在全周期或者在非稳态阶段的预期维修费用率，分别得到使产品预期维修费用率最低的预防维修时间，即为产品在全周期或者在非稳态阶段的最优预防维修时间。

本发明有益效果如下：与现有技术不同，本发明实施例全面考虑了两阶段衰退过程中，冲击可能对衰退过程造成的各种影响，同时该策略中在观测到产品进入非稳态阶段后仍允许产品工作一段时间，对于经历两阶段衰退与竞争风险且在衰退非稳态阶段仍可正常工作的产品，本发明实施例的方法对有效开展产品预防维修具有很强的指导意义。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分的从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。

图1为本发明实施例方法的流程图；

图2为两类相关竞争失效过程：(a)软失效：衰退过程；(b)硬失效：冲击失效；

图3为样本稳态阶段的区间长度或失效时间获取流程；

图4为样本非稳态阶段失效时间获取流程；

图5为综合预防维修策略执行流程；

图6为产品全寿命周期内失效时间累积分布函数；

图7为产品非稳态阶段内失效时间累积分布函数；

图8为各种预防维修时刻情况下产品全寿命周期的预期维修费用率；

图9为仅考虑产品非稳态阶段的各种预防维修时刻的预期维修费用率；

图10为不同状态观测间隔情况下两种预防维修策略表现。

具体实施方式

下面结合附图来具体描述本发明的优选实施例，其中，附图构成本发明一部分，并与本发明的实施例一起用于阐释本发明的原理。

本发明提供了一种两阶段衰退与冲击竞争风险的综合预防维修方法，如图1所示，具体包括：

S1、同时考虑产品稳态阶段与非稳态阶段两阶段衰退与冲击竞争风险并以此建立产品失效过程模型；

建立两阶段衰退与冲击竞争风险失效过程模型，即建立考虑产品稳态阶段与非稳态阶段两阶段衰退、冲击两种失效因素，同时考虑冲击对衰退过程的可能影响的失效过程模型。如图2所示，演示了两类相关失效过程：由不断衰退导致的软失效和由随机冲击导致的硬失效。

如图2所示，产品受到第i次冲击时，当W_i大于D_hf则产品立即失效，即发生了硬失效，其中，t_i为产品自投入使用始受到第i次冲击的时间，W_i为冲击力大小，D_hf为冲击失效阈值。

t_sw为产品衰退过程由稳态阶段向非稳态阶段的切换时刻，假设产品稳态阶段与非稳态阶段的区间长度无关，t_sw为均值为μ_sw，方差为的正态分布的正数。X_s(t)为产品在时刻t的衰退变量值。通常而言，产品在稳态阶段其衰退变量值保持稳定，并在一常值附近上下波动，而产品在非稳态阶段的衰退变量呈现趋势明显的变化。考虑最常见的情况，定义产品在稳态阶段的衰退变量值：

其中，为均值为μ₀、方差为的正态分布的截距参数，ε(t)为均值为零，方差为σ²的正态分布的过程噪声。

在产品进入非稳态阶段后，非致命冲击则影响产品衰退变量值X_s(t)，本实施例采用随机线性衰退模型刻画X_s(t)，即：

其中，β_k为第k次状态转移之前，第(k-1)次状态转移之后的衰退速率参数。特别地，β₁为第一次状态转移之前的均值为μ₁、方差为σ₁²的正态分布衰退速率参数衰退速率，t_ck为产品进入非稳态阶段后发生第k次状态转移的时间；Y_i对应第i次冲击导致的衰退增加程度；ε(t)为均值为零，方差为σ²的正态分布的过程噪声；n(t)为产品衰退进入非稳态阶段后至时刻t发生的冲击数；nc(t)为产品衰退进入非稳态阶段后至时刻t发生的状态转移数，β_nc(t)+1是第nc(t)+1次状态转移之前，第nc(t)次状态转移之后的衰退速率；t_nc(t)为第nc(t)次状态转移发生的时间。

当产品进入非稳态阶段后，非致命失效以两种方式影响着产品的衰退过程：当W_i小于状态转移阈值D_e时，导致产品衰退过程的衰退程度增加Y_i；当W_i大于D_e且小于D_hf时，记为一次非致命严重冲击(状态转移)，导致Y_i、β_k增加，其中，Y_i对应第i次冲击导致的衰退增加程度；β_k对应第k次状态转移之前的均值为μ₁、方差为的正态分布衰退速率参数，β₁为产品进入非稳态阶段后的初始衰退速率；软失效方面，当衰退值X_s(t)大于衰退损伤阈值D_sf则产品失效。

当产品发生第k次状态转移时，产品衰退速率由β_k变为β_k+1，定义：

β_k+1＝β_k+η_k(3)

其中η_k为均值为μ_η、方差为正态分布的正随机变量；

在产品衰退稳态阶段各冲击的到达时间符合到达率为λ₁(t)的齐次非齐次Poisson分布，而在产品衰退非稳态阶段各冲击的到达时间符合到达率为λ₂(t)的齐次非齐次Poisson分布，且有λ₂(t)≥λ₁(t)。

W_i为一系列独立同分布的非负随机变量，且有相同的累积分布函数G(w)，G(w)＝P(W_i＜w)，其中G(w)与D_e、D_hf无关，w为一个代表冲击大小的随机变量。

根据冲击大小的累积分布函数，则对于第i次冲击产品不发生硬失效的概率为：

P(W_i＜D_hf)(4)

另一方面，Y_i与W_i相关，记为：

Y_i＝Q(W_i)(5)

其中，Q(W_i)为一递增的映射函数。

这样，给定G(w)、Q(W_i)的具体形式，及式(1)～(5)中各参数的具体取值，即可开展产品在衰退冲击竞争衰退过程的可靠性分析。

S2、采用步骤S1建立的失效过程模型进行仿真，得到样本在稳态阶段的失效时间或稳态阶段的区间长度，以及样本在非稳态阶段的失效时间：

本实施例采用大量样本取样方法分析评估产品的可靠性，即基于大量样本的实际失效过程数据分析产品的可靠性，其中每一个样本均经历从开始运行直至失效的过程。对于本发明实施例考察的两阶段衰退过程，由于产品稳态阶段与非稳态阶段的冲击与衰退均相对独立，提出采用并行仿真的样本失效时间获取方法。

S21、采用步骤S1建立的失效过程模型进行仿真，得到样本在稳态阶段的失效时间或稳态阶段的区间长度：

S211、确定仿真样本数量N，时间轴长度T_s1，仿真步长Δt；

仿真样本数量N一般应较大以使得整个仿真结果具有较高的准确性，时间轴长度T_s1一般应较长以使得全部样本在T_s1内能进入非稳态阶段，即，T_s1＞max(t_swn)；仿真步长△t一般应远小于T_s1同时满足在△t内不会发生两次冲击为宜。

S212、选取样本n，初始值n＝1；仿真产生样本n的稳态阶段的区间长度t_swn；

S213、基于强度为λ₁(t)的Poisson过程产生样本n的冲击序列，对应时间为t₁、t₂、…t_M，大小为W₁、W₂、…W_M，其中t_M≤T_s1；

S214、定义冲击时间t_Z，其中t_z＜t_M，使得t_z≤t_swn≤t_z+1；

其中t_swn为样本n衰退过程由稳态阶段向非稳态阶段的切换时刻；

S215、判断Max(W₁、W₂、…W_z)是否大于等于D_hf，否，则IND_n＝0，样本序号n＝n+1，并执行步骤S218；是，则执行步骤S216；

其中IND用于指示产品是否在稳态区间失效，如果是，则IND_n＝1；如果否，则IND_n＝0。

S216、确定第一个导致失效的冲击时刻t_j，1≤j≤z，该冲击的大小为W_j；

S217、记录样本n的失效时间t_f1＝t_j，置IND_n为1；样本序号n＝n+1；

S218、判断n是否大于N，否，则执行步骤S212；是，则结束仿真，得到样本在稳态阶段的失效时间t_f1或稳态阶段的区间长度t_swn。

更优地，本发明是基于Monte Carlo仿真的可靠性分析评估，如图3给出了基于Monte Carlo仿真的流程。通常而言，需选取较大的样本量(例如：N＝1000)，此外仿真时间轴长度T_s1应足够长也确保所有仿真样本均在仿真结束前失效。

S22、采用步骤S1建立的失效过程模型进行仿真，得到样本在非稳态阶段的失效时间：

S221、确定仿真样本数量N，时间轴长度T_s2，仿真步长△t；

仿真样本数量N一般应较大以使得整个仿真结果具有较高的准确性，时间轴长度T_s2一般应较长以使得全部样本在T_s2以内失效；仿真步长△t一般应远小于T_s2同时满足在△t内不会发生两次冲击为宜。

S222、选取样本n，初始值n＝1；

S223、基于强度为λ₂(t)的Poisson过程产生样本n的冲击序列，对应时间为t₁、t₂、…t_M，大小为W₁、W₂、…W_M，其中t_M≤T_s2，冲击次数标识i，初始值i＝1；状态转移标识k，初始值k＝1；

S224、进入仿真，仿真产生ε(t)、W_i、β_k、

S225、判断t是否大于等于t_i，否，X_s(t)＝X_s(t-Δt)+β_kΔt+ε(t)，并执行S226；是，判断W_i是否大于等于D_hf；是，则执行步骤S228，否，执行S227；

t_i代表第i次冲击的时间，判断当前仿真时刻是否达到产生第i次冲击时刻。

S226、判断Xs(t)是否大于等于D_sf；否，t＝t+Δt，返回步骤S224，是，则执行S228；

S227、判断W_i是否大于等于D_e；

是，则仿真产生η_k，i＝i+1、k＝k+1，β_k＝β_(k-1)+η_(k-1)；

X_s(t)＝X_s(t-Δt)+β_k(t-t_k)Δt+β_(k-1)(t_k-t+Δt)+Y_i+ε(t)并执行步骤S226；

否，则Xs(t)＝Xs(t-Δt)+β_kΔt+Y_i+ε(t)，i＝i+1并执行步骤S226；

判断冲击是否为大于D_e的非致命严重冲击，并按是与不是非致命严重冲击的情况分别考虑产品衰退的发展情况。

S228、记录样本n的失效时间t_f2＝t，样本序号n＝n+1；并判断n是否大于N；否，则执行步骤S222，是，则仿真结束，得到样本在非稳态阶段的失效时间。

更优地，本发明是基于Monte Carlo仿真的可靠性分析评估，如图4所示，给出了基于Monte Carlo仿真的流程。通常而言，需选取较大的样本量(例如：N＝1000)，此外仿真时间轴长度T_s2应足够长也确保所有仿真样本均在仿真结束前失效。

S3、根据步骤S2获得的数据，得到产品失效时间累积分布函数；

基于图3中获取的稳态阶段的区间长度t_swn或失效时间t_f1，以及图4中获取的非稳态阶段失效时间t_f2，样本n的最终失效时间定义为：

在获得了这些样本的失效时间数据后，即可根据常用的经验累积分布函数估计方法(即Kaplan-Meier估计方法)估计产品失效时间的累积分布函数F₁(t)，则产品的可靠性函数即为1-F₁(t)。当然，仅基于图4中获取的产品非稳态阶段的失效时间数据，同样可根据常用的经验累积分布函数估计方法，估计得到产品非稳态阶段的失效时间累积分布函数F₂(t)。

S4、利用步骤S3得到的产品失效时间累积分布函数建立预防维修策略并优化；

本实施例中，预防维修要素共包括四个：

维修目标：最小化产品的维修费用率；

维修效果：预防维修与事后维修均“修复如新”；

维修策略：顺序型维修，即产品的相邻两次维修时间没有限制，完全根据优化结果确定；

维修限制：无，即不存在备件等候，产品可用度指标等约束。

具体包括以下步骤：

S41、根据产品失效时间的累积分布函数计算和

基于上述假设，产品在一个维修周期内的预期维修费用率为：

t_pm1为规划预防维修时间，C_pm、C_cm分别为开展一次预防维修与事后维修的费用，E[t_o]为一个维修周期内产品预期工作时间。

这样，根据式(7)，使得产品预期维修费用率最低的预防维修时间为：

其中，是针对产品全寿命周期的最优预防维修时间。随着状态监测技术的不断发展，有望对产品的衰退过程予以观测并指导维修策略的制定。由于产品非稳态阶段往往更接近产品失效时刻，在观测到产品进入非稳态阶段后，专门针对非稳态阶段制定预防维修策略有望获得更好的维修效果，因此本发明实施例在开展产品全寿命周期预防维修规划的基础上，针对非稳态阶段仍开展预防维修规划。

产品在非稳态阶段的预期维修费用率为：

其中，t_pm2为产品衰退进入非稳态阶段后的规划预防维修时间。同样地，使得产品进入非稳态阶段后，维修费用率最低的预防维修时间为：

其中，是产品在非稳态阶段的最优预防维修时间。

S42、产品开始工作，t＝0；

S43、每隔Δt_in观测一次产品衰退状态；

Δt_in为产品状态观测间隔，当一次状态观测的费用较高时，可先根据历史样本数据，采用交叉校验方法优化选取Δt_in。

S44、判断产品是否进入非稳态阶段，是，则在后开展预防维修；否，则执行步骤S45；

对于在时刻之前就进入非稳态阶段的产品，则自观测到产品进入非稳态时刻起，在时间后开展预防维修。由于当Δt_in较大时无法确知产品进入非稳态阶段的具体时刻，扣除Δt_in项是因此而采取的一种保守方案。

具体地，可根据工程经验设定衰退变量阈值来判断产品是否进入非稳态阶段。此外，执行预防维修策略的过程中，如果产品突然失效则立即开展事后维修。

S45、判断当前时刻是否大于等于时刻，否，则执行步骤S43；是，则立即开展预防维修。

对于某些长时间处于稳态阶段而未进入非稳态阶段的产品，为避免冲击导致的硬失效，仍按全寿命周期预防维修策划优化的结果在时刻开展预防维修。

在一个具体的实施例中，定义μ₀＝4，σ₀＝1，μ₁＝4，σ₁＝1，σ＝1，μ_η＝0，σ_η＝1，λ₁(t)＝8，λ₂(t)＝5，D_e＝0.7，D_hf＝0.95，D_sf＝200，X_i为均值为0，方差为1的非负正态分布变量，Y_i＝X_i。此外，μ_SW＝60，σ_SW＝20。选取仿真样本N_s＝1000，仿真步长Δt＝1。按照图1中步骤2产品可靠性分析评估方法，产品全寿命周期内失效时间累积分布函数F₁(t)如图6所示，产品非稳态阶段内失效时间累积分布函数F₂(t)如图7所示。

为制定产品预防维修策略，定义参数C_pm＝1000，C_cm＝4000，式(7)与式(9)分别给出了各种预防维修时刻情况下产品全寿命周期的预期维修费用率与仅考虑产品非稳态阶段的各种预防维修时刻的预期维修费用率，由图8与图9中结果可见，

基于图8中的结果，若采用传统的基于产品全寿命周期的预防维修策略(即所有产品均在开展预防维修)，则所有样本的实际平均维修费用率为：

其中，N_pm为预防维修策略中接收预防维修的产品数，N_cm为预防维修策略中接收事后维修的产品数，T_oj为产品j的实际工作时间，δ_j为产品j的预防维修指针，若产品j接收预防维修则δ＝1，反之则δ＝0。T_pmj与T_fj分别为产品j的预防维修时间与失效时间。

基于式(11)，计算得到传统的基于产品全寿命周期的预防维修策略下，产品实际平均维修费用率为：r_ac1＝36.56。

假设一次状态观测的费用C_in＝10，则不同状态观测间隔情况下，基于图5中的综合维修策略的实际平均维修费用率为：

其中，N_in为全部产品的状态观测总次数。基于式(12)，图10给出了不同状态观测间隔情况下，本发明实施例提出的综合预防维修策略的表现。作为比较，图10还给出了预防维修策略2(即一旦观测到产品进入非稳态阶段或在进入非稳态阶段前即达到时开展预防维修)的表现。由图10中的结果可见，选择恰当的状态观测间隔，本发明实施例提出的综合预防维修策略表现优于传统的基于产品全寿命周期的预防维修策略，同样优于预防维修策略2的表现，证明了本专利提出的综合预防维修策略的有效性。这也说明，对于那些进入非稳态阶段但仍能开展一段时间有效工作的产品，在观测到产品进入非稳态阶段后仍允许产品适当工作一段时间是有益的。

本发明有益效果如下：

与现有技术不同，本发明全面考虑了两阶段衰退过程中，冲击可能对衰退过程造成的各种影响，同时该策略中在观测到产品进入非稳态阶段后仍允许产品工作一段时间，对于经历两阶段衰退与竞争风险且在衰退非稳态阶段仍可正常工作的产品，本发明实施例的方法对有效开展产品预防维修具有很强的指导意义。

本领域技术人员可以理解，实现上述实施例方法的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于计算机可读存储介质中。其中，所述计算机可读存储介质为磁盘、光盘、只读存储记忆体或随机存储记忆体等。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。