考虑竞争失效的多阶段二元随机退化系统可靠性建模与分析

Reliability Modeling and Analysis of Considering Competition Failure Multi-Stage Bivariate Stochastic Degradation Systems

摘要: 针对具有两个相关性能指标的随机退化系统，考虑竞争失效和退化过程的多阶段特征，采用Copula函数描述性能指标间的相关关系，建立考虑竞争失效的多阶段二元随机退化系统模型。推导了系统可靠度的解析表达式，提出基于蒙特卡洛方法模拟系统可靠度的流程步骤，得到系统可靠度的模拟解，并验证了解析解的正确性。最后，通过数值算例验证模型的有效性。

Abstract: In the paper, for stochastic degenerate systems with two correlated performance indexes, considering the multi-stage characteristics of competition failure and degradation process, and by using Copula function to describe the correlation between performance indicators, a multi-stage binary stochastic degenerate system model considering competition failure is established. Explicit expression of the system reliability is given, and a numerical simulation algorithm based on Monte Carlo simulation is proposed, which verifies the correctness of the analytical results. Finally, some numerical examples are given to show the effectiveness of the proposed model.

1. 引言

现代产品的制造工艺不断进步，其内部结构越来越复杂，受到应用环境的动态性和工作模式变化性的影响，传统的基于失效数据的可靠性建模与评估方法用于分析此类产品的可靠度存在一定困难 [1]。目前，基于产品退化数据的可靠性建模与评估方法为解决这类问题提供了一条有效途径，已经成为当前可靠性理论与工程领域的研究热点之一。

随着研究的深入，学者们发现单一性能指标不足以反映复杂系统整体上的运行状态。目前，学界主要使用Copula函数来描述多个性能指标退化过程的相关性。杨志远等 [2] 构建了多性能非线性Wiener退化过程模型，利用Copula函数描述相关性，并分析了系统可靠度。韩玉成等 [3] 基于Copula函数给出了混合随机退化过程的可靠性建模方法。胡启国等 [4] 采用Vine Copula函数研究了多失效模式下的随机退化可靠性模型。周义蛟 [5] 和刘小平 [6] 研究了基于Wiener退化过程和Copula函数的二元相关性能可靠性模型。鲍兆伟 [7] 提出了确定边缘退化量分布函数和Copula函数的方法，进一步提高了可靠性模型的精确性。

另外，产品还会遭受各种随机冲击的影响，当这些冲击对系统的损伤超出一定范围时，就会导致系统失效。根据冲击类型的不同目前主要分为累积冲击过程、极值冲击过程和 $\delta$ -冲击过程等 [8]。刘汉葱 [9] [10] 建立了不同随机冲击影响下的二元相关随机退化系统可靠性模型。杨志远等 [11] 基于Gamma过程和Copula函数，考虑冲击过程对退化过程的影响，退化过程间以及冲击与退化过程间的相关关系，建立了基于退化相关性分析的竞争失效系统可靠性模型。黄文平等 [12] 在阈值变动情况下，假定退化过程为线性退化过程，分别建立了考虑极值冲击、累积冲击和 $\delta$ -冲击的可靠性模型。

同时，很多产品的退化过程还会表现出多阶段性，表现为在不同退化阶段产品的退化速率不同。例如，锂电池容量随着汽车运行环境的变化，充放电次数的增加，电容和阻抗等性能指标的退化曲线前期退化趋势较缓和，后期的退化趋势则较激烈迅速 [13]。Dong等 [14] [15] 利用Wiener过程描述退化，建立了动态环境下的多阶段随机退化系统的可靠性模型，并给出了可靠度解析表达式。Gao等 [16] 提出了一种基于Wiener过程来描述失效阈值和漂移系数存在突变现象的多阶段随机退化模型。Shi等 [17] 利用一般退化模型描述退化，建立了多阶段随机退化系统的可靠性模型。Farid等 [18] 针对多阶段线性退化过程中的变点分布进行了细致研究。

不难发现，现有考虑竞争失效的退化系统可靠性模型以单阶段多元随机退化系统和多阶段一元随机退化系统的可靠性建模与评估为研究对象，而考虑竞争失效的多阶段多元随机退化系统的可靠性建模与评估研究却鲜有报道。在工程实际中，例如，锂离子电池的退化过程会呈现多阶段退化过程 [13] ，同时，电池容量会随着工作时间的增加而不断减少，电池阻抗也将随着锂离子电池容量的减少而发生变化 [19]。本文将针对上述问题，提出一类考虑竞争失效的多阶段二元多阶段随机退化模型，并推导系统的可靠度函数。

2. 模型描述

假设系统的退化过程由两个相关性能指标表征，在不同阶段退化速率不同。同时外界的冲击分为致命冲击和非致命冲击两种模式，致命冲击会引起系统直接失效，非致命冲击会引起性能退化量的增加，这时系统实际的退化量由退化过程的退化量和冲击引起的累积增量共同组成。因此，系统具有极值冲击引发的硬失效以及累积冲击引发的软失效两种失效模式，下面给出此类系统的具体模型。

2.1. 二元多阶段退化模型

假设各性能指标的退化轨迹为多阶段退化轨道，令 $X\left(t\right)$ 表示 $t$ 时刻的系统退化量。由模型描述可知系统存在两个性能相关的退化过程，记 $X_{\text{1}}\left(t\right)$ 和 $X_{\text{2}}\left(t\right)$ 分别表示性能指标1和性能指标2的退化过程。则有：

$\{\begin{array}{l}{X}_i\left(t\right)=X_{i,\text{1}}\left(t\right)I_{[{\tau }_{\text{0}},{\tau }_{\text{1}})}\left(t\right)+X_{i,\text{2}}\left(t\right)I_{[{\tau }_{\text{1}},{\tau }_{\text{2}})}\left(t\right)+\cdots +X_{i,j}\left(t\right)I_{[{\tau }_{j-\text{1}},{\tau }_j)}\left(t\right),\hfill \\ X_{i,j}\left(t\right)=X_i\left({\tau }_{j-\text{1}}\right)+{\beta }_{i,j}\left(t-{\tau }_{j-\text{1}}\right),\hfill \end{array}$ (1)

其中 $X_i\left(t\right)$ 为性能指标 $i$ 在时刻 $t$ 的退化量， $X_{i,j}\left(t\right)$ 为性能指标 $i\left(i=1,2\right)$ 在第 $j\left(j\in N^{*}\right)$ 阶段时刻 $t$ 的退化量。 ${\tau }_j$ 为退化过程中退化性能从第 $j$ 阶段到第 $j+\text{1}$ 阶段的变点时刻，其中 ${\tau }_{\text{0}}=\text{0}$。 $X_i\left({\tau }_{j-\text{1}}\right)$ 为性能指标 $i$ 在时刻 ${\tau }_{j-\text{1}}$ 的退化量，即性能指标 $i$ 在第 $j$ 阶段的初始退化量，且性能指标 $i$ 在整个退化过程中的初始退化量均为0，即 $X_i\left({\tau }_{\text{0}}\right)=\text{0}$。 $I_{[{\tau }_{j-\text{1}},{\tau }_j)}\left(t\right)$ 表示示性函数，当 $t\in \left[{\tau }_{j-\text{1}},{\tau }_j\right)$ 时， $I_{[{\tau }_{j-1},{\tau }_j)}\left(t\right)=1$，否则 $I_{[{\tau }_{j-1},{\tau }_j)}\left(t\right)=0$。 ${\beta }_{i,j}$ 表示性能指标 $i$ 在第 $j$ 阶段的退化速率，假设 ${\beta }_{i,j}~N\left({\mu }_{i,j},{\sigma }_{i,j}^{\text{2}}\right)$。

2.2. 相关性描述

如上所述，本文采用Copula函数来描述退化过程的相关性。根据二维Sklar定理 [20] 可知， $X_{\text{1}}\left(t\right)$ 和 $X_{\text{2}}\left(t\right)$ 在时间区间 $\left[t,t+\Delta t\right]\subset \left[{\tau }_{j-\text{1}},{\tau }_j\right)$ 上的退化增量 $\Delta X_{\text{1}}\left(t\right)$ 和 $\Delta X_{2,j}\left(t\right)$ 的联合概率分布函数可表示为：

$H\left(x_{1,j},x_{2,j}\right)=C\left(F_{\Delta X_{1,j}(t)}\left(x_{1,j}\right),F_{\Delta X_{2,j}(t)}\left(x_{2,j}\right);\theta \right),$ (2)

其中 $H(\cdot )$ 为两个性能的联合分布函数， $C(\cdot )$ 表示Copula函数， $\theta$ 为描述相关性强弱关系的Copula函数参数。其对应的概率密度函数为：

$h\left(x_{1,j},x_{2,j}\right)=\frac{{\partial }^{2}H\left(x_{1,j},x_{2,j}\right)}{\partial x_{1,j}\partial x_{2,j}}=c\left(F_{\Delta X_{1,j}(t)}\left(x_{1,j}\right),F_{\Delta X_{2,j}(t)}\left(x_{2,j}\right);\theta \right)\times {\displaystyle \underset{i=1}{\overset{2}{\prod }}{f}_{\Delta X_{i,j}(t)}\left(x_{i,j}\right)}.$ (3)

式中

$c\left(F_{\Delta X_{1,j}(t)}\left(x_{1,j}\right),F_{\Delta X_{2,j}(t)}\left(x_{2,j}\right);\theta \right)=\frac{{\partial }^{2}C\left(F_{\Delta X_{1,j}(t)}\left(x_{1,j}\right),F_{\Delta X_{2,j}(t)}\left(x_{2,j}\right);\theta \right)}{\partial F_{\Delta X_{1,j}(t)}\left(x_{1,j}\right)\partial F_{\Delta X_{2,j}(t)}\left(x_{2,j}\right)}.$ (4)

常见的二维Copula函数有Frank Copula、Gumbel Copula、Clayton Copula和Gaussian Copula函数。以Frank Copula函数为例，其分布函数和密度函数为：

$C\left(u,v;\theta \right)=-\frac{1}{\theta }\ln \left(1+\frac{\left(\exp \left(-\theta u\right)-1\right)\left(\exp \left(-\theta v\right)-1\right)}{\exp \left(-\theta \right)-1}\right),$ (5)

$c\left(u,v;\theta \right)=\frac{-\theta \left(\exp \left(-\theta \right)-1\right)\left(\exp \left(-\theta \left(u+v\right)\right)\right)}{{\left[\left(\exp \left(-\theta u\right)-1\right)\left(\exp \left(-\theta v\right)-1\right)+\exp \left(-\theta \right)-1\right]}^2}.$ (6)

2.3. 冲击过程模型

假设外界冲击为只有1种冲击类型的极值冲击模型，冲击到达的强度为速率 $\lambda$ 的齐次泊松过程，同时冲击对处于同一时刻的性能指标1和性能指标2的影响相同且独立。记 $W_k\left(k\in N\right)$ 表示第 $k$ 次冲击量的大小，该冲击量为独立同分布的随机变量，假设 $W_k~N\left({\mu }_w,{\sigma }_w^2\right)$， $D_w$ 表示该模型冲击对应的失效阈值。则系统在第 $k$ 次冲击中不失效的概率为：

$P\left\{W_k<D_w\right\}=\Phi \left(\frac{D_w-{\mu }_w}{{\sigma }_w}\right),$ (7)

其中 $\Phi \left(x\right)={\displaystyle {\int }_{-\infty }^x\exp \left(-z^2/2\right)\text{d}z}/\sqrt{2\pi }$ 为标准正态分布的分布函数。则系统此时的可靠度为：

$R_w\left(t\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}P\left\{W_k<D_w|N\left(t\right)=k\right\}P\left\{N\left(t\right)=k\right\}}={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\Phi \left(\frac{D_w-{\mu }_w}{{\sigma }_w}\right)\frac{{\left(\lambda t\right)}^k{\text{e}}^{-\lambda t}}{k!}}.$ (8)

假设外界冲击为只有1种冲击类型的累积冲击模型，冲击到达的强度为速率 $\lambda$ 的齐次泊松过程，同时冲击对处于同一时刻的性能指标1和性能指标2的影响相同。记 $Y_n\left(n\in N\right)$ 表示第 $k$ 次冲击导致的性能退化量，该退化量为独立同分布的随机变量，假设 $Y_n~N\left({\mu }_y,{\sigma }_y^2\right)$。则 $t$ 时刻外界冲击对系统的累积退化量为：

$S_i\left(t\right)=X_i\left(t\right)+Y\left(t\right)={\displaystyle \underset{m=1}{\overset{j}{\sum }}{X}_{i,m}\left(t\right)I_{[{\tau }_{m-\text{1}},{\tau }_m)}\left(t\right)}+{\displaystyle \underset{k=0}{\overset{N(t)}{\sum }}{Y}_n}.$ (9)

注：

1) 如果 $t\in \left[{\tau }_0,{\tau }_1\right)$，则上述模型变为考虑竞争失效的单阶段二元随机退化系统模型，文献 [10] [11] 对此进行了研究。

2) 如果 $t\in \left[{\tau }_0,{\tau }_2\right)$ 且 $v=1$ 或 $u=1$，则上述模型变为考虑竞争失效的两阶段一元随机退化系统模型，文献 [21] 对此进行了研究。

3. 可靠性分析

3.1. 固定变点情形下考虑竞争失效的二元多阶段随机退化系统的可靠度

本节将建立固定变点情形下的可靠度模型，仅考虑累积冲击过程，倘若系统的任意一个性能指标的总的性能退化量大于等于失效阈值时，系统失效。记两个性能指标的失效阈值分别 $D_1$ 和 $D_2$，故性能指标i的可靠度表达式为：

$\begin{array}{c}{R}_i\left(t\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{S_i\left(t\right)<D_i|N\left(t\right)=n\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{{\displaystyle \underset{m=1}{\overset{j}{\sum }}{X}_{i,m}\left(t\right)I_{[{\tau }_{m-\text{1}},{\tau }_m)}\left(t\right)}+{\displaystyle \underset{k=0}{\overset{N(t)}{\sum }}{Y}_n}<D_i\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\Phi \left(\frac{D_i-\left({\displaystyle {\sum }_{m=0}^j\left(X_i\left({\tau }_{m-\text{1}}\right)+\left(t-{\tau }_{j-\text{1}}\right){\mu }_{i,m}\right)I_{[{\tau }_{m-\text{1}},{\tau }_m)}\left(t\right)}+\left(n+1\right){\mu }_y\right)}{\sqrt{{\displaystyle {\sum }_{m=0}^j{\left(\left(t-{\tau }_{j-\text{1}}\right){\sigma }_{i,m}\right)}^2{I}_{[{\tau }_{m-\text{1}},{\tau }_m)}\left(t\right)}+\left(n+1\right){\sigma }_y^2}}\right)\frac{{\left(\lambda t\right)}^n{\text{e}}^{-\lambda t}}{n!}},\end{array}$ (10)

令 $T_i=\inf \left\{t\ge 0|Y_i\left(t\right)\ge D_i\right\}$ 表示性能指标 $i$ 首次到达失效阈值的时间，则系统整个的寿命 $T=$ $\min \left\{T_1,T_2\right\}$，故系统的可靠度函数 $R\left(t\right)$ 为：

$\begin{array}{c}R\left(t\right)=P\left\{T>t\right\}=P\left\{T_1>t,T_2>t\right\}=1-P\left\{T_1\le t\right\}-P\left\{T_2\le t\right\}+P\left\{T_1\le t,T_2\le t\right\}\\ =R_{1,j}\left(t\right)+R_{2,j}\left(t\right)-1+P\left\{S_1\left(t\right)>D_1,S_2\left(t\right)>D_2\right\}\\ =R_{1,j}\left(t\right)+R_{2,j}\left(t\right)-1+C\left(F_{1,j}\left(t\right),F_{2,j}\left(t\right);\theta \right).\end{array}$ (11)

本节假设系统产生变点的时间固定已知，则当系统运行时刻 $t$ 处在不同的阶段时，其初始退化量为随机变量。

如果 $t\in \left[0,{\tau }_1\right)$，则该阶段的系统可靠度为：

$R_1\left(t\right)={\displaystyle \underset{i=1}{\overset{2}{\sum }}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\Phi \left(\frac{D_i-\left(X_i\left({\tau }_0\right)+t{\mu }_{i,1}+\left(n+1\right){\mu }_y\right)}{\sqrt{t^2{\sigma }_{i,1}^2+\left(n+1\right){\sigma }_y^2}}\right)}}\frac{{\left(\lambda t\right)}^n{\text{e}}^{-\lambda t}}{n!}-1+C\left(F_{1,1}\left(t\right),F_{2,1}\left(t\right);\theta \right),$ (12)

其中 $F_{i,1}\left(t\right),i=1,2$ 的表达式为：

$F_{i,1}\left(t\right)=1-{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\Phi \left(\frac{D_i-\left(X_i\left({\tau }_0\right)+t{\mu }_{i,1}+\left(n+1\right){\mu }_y\right)}{\sqrt{t^2{\sigma }_{i,1}^2+\left(n+1\right){\sigma }_y^2}}\right)}.$ (13)

如果 $t\in \left[{\tau }_1,{\tau }_2\right)$，则该阶段性能指标 $i$ 的可靠度：

$\begin{array}{c}{R}_{i,2}\left(t\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{S_i\left(t\right)<D_i,\underset{0\le u<{\tau }_1}{\sup }{S}_i\left(u\right)<D_i|N\left(t\right)=n\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{S_i\left(t\right)<D_i,\underset{0\le u<{\tau }_1}{\sup }{S}_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{S_i\left(t\right)<D_i|\underset{0\le u<{\tau }_1}{\sup }{S}_i\left(u\right)<D_i\right\}P\left\{\underset{0\le u<{\tau }_1}{\sup }{S}_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_1}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}}P\left\{S_i\left(t\right)<D_i|\underset{0\le u<{\tau }_1}{\sup }{S}_i\left(u\right)<D_i\right\}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_1}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}{\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t\right)<D_i|\underset{0\le u<{\tau }_1}{\sup }{S}_i\left(u\right)<D_i,S_i\left({\tau }_1\right)=s_{i,1}\right\}{f}_{S_i({\tau }_1)}\left(s_{i,1}\right)\text{d}{s}_{i,1}}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_1}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}{\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t-{\tau }_1\right)<D_i-s_{i,1}\right\}{f}_{S_i({\tau }_1)}\left(s_{i,1}\right)\text{d}{s}_{i,1}}},\end{array}$ (14)

其中 $P\left\{S_i\left(u\right)<D_i\right\}$， $P\left\{N\left(t\right)=n\right\}$ 和 $P\left\{S_i\left(t-{\tau }_1\right)<D_i-s_{i,1}\right\}$ 的表达式分别为：

$P\left\{S_i\left(u\right)<D_i\right\}=\Phi \left(\frac{D_i-\left(u{\mu }_{i,1}+\left(n+1\right){\mu }_y\right)}{\sqrt{u^2{\sigma }_{i,1}^2+\left(n+1\right){\sigma }_y^2}}\right),$ (15)

$P\left\{N\left(t\right)=n\right\}\text{=}\frac{{\left(\lambda t\right)}^n{\text{e}}^{-\lambda t}}{n!},$ (16)

$P\left\{S_i\left(t-{\tau }_1\right)<D_i-s_{i,1}\right\}=\Phi \left(\frac{D_i-s_{i,1}-\left(\left(t-{\tau }_1\right){\mu }_{i,2}+\left(n+1\right){\mu }_y\right)}{\sqrt{{\left(t-{\tau }_1\right)}^2{\sigma }_{i,1}^2+\left(n+1\right){\sigma }_y^2}}\right).$ (17)

由于 $E\left[S_i\left({\tau }_1\right)\right]=E\left[X_i\left({\tau }_1\right)+Y\left({\tau }_1\right)\right]={\tau }_1{\mu }_{i,1}+\left(n+1\right){\mu }_y$， $Var\left[S_i\left({\tau }_1\right)\right]={\tau }_1{}^2{\sigma }_{i,1}^2+\left(n+1\right){\sigma }_y^2$，故：

$f_{S_i({\tau }_1)}\left(s_{i,1}\right)=\frac{1}{\sqrt{2\pi \left({\tau }_1{}^2{\sigma }_{i,1}^2+\left(n+1\right){\sigma }_y^2\right)}}\exp \left(-\frac{{\left(s_{i,1}-{\tau }_1{\mu }_{i,1}-\left(n+1\right){\mu }_y\right)}^2}{2\left({\tau }_1{}^2{\sigma }_{i,1}^2+\left(n+1\right){\sigma }_y^2\right)}\right).$ (18)

则该阶段的可靠度为：

$\begin{array}{c}{R}_2\left(t\right)={\displaystyle \underset{i=1}{\overset{2}{\sum }}{R}_{i,2}\left(t\right)}-1+C\left(F_{1,2}\left(t\right),F_{2,2}\left(t\right);\theta \right)\\ {\displaystyle \underset{i=1}{\overset{2}{\sum }}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_1}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}}}\\ {\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t-{\tau }_1\right)<D_i-s_{i,1}\right\}{f}_{S_i({\tau }_1)}\left(s_{i,1}\right)\text{d}{s}_{i,1}}-1+C\left(F_1\left(t\right),F_2\left(t\right);\theta \right),\end{array}$ (19)

其中 $F_{i,2}\left(t\right),i=1,2$ 的表达式为：

$F_{i,2}\left(t\right)=1-{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_1}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}{\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t-{\tau }_1\right)<D_i-s_{i,1}\right\}{f}_{S_i({\tau }_1)}\left(s_{i,1}\right)\text{d}{s}_{i,1}}}.$ (20)

如果 $t\in \left[{\tau }_{j-1},{\tau }_j\right),j>2$，则该阶段性能指标 $i$ 的可靠度为：

$\begin{array}{c}{R}_{i,j}\left(t\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{S_i\left(t\right)<D_i,{\displaystyle \underset{m=2}{\overset{1}{\cap }}\left\{\underset{{\tau }_{m-2}\le u<{\tau }_{m-1}}{\sup }{S}_i\left(u\right)<D_i\right\}}|N\left(t\right)=n\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{S_i\left(t\right)<D_i,\underset{0\le u<{\tau }_{j-1}}{\sup }{S}_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}P\left\{S_i\left(t\right)<D_i|\underset{0\le u<{\tau }_{j-1}}{\sup }{S}_i\left(u\right)<D_i\right\}P\left\{\underset{0\le u<{\tau }_{j-1}}{\sup }{S}_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_{j-1}}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}}P\left\{S_i\left(t\right)<D_i|\underset{0\le u<{\tau }_{j-1}}{\sup }{S}_i\left(u\right)<D_i\right\}\end{array}$

$\begin{array}{c}={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_{j-1}}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}{\displaystyle {\int }_{-\infty }^{D_i}\cdots {\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t\right)<D_i|\underset{0\le u<{\tau }_{j-1}}{\sup }{S}_i\left(u\right)<D_i,{\displaystyle \underset{m=2}{\overset{j}{\cap }}{S}_i\left({\tau }_{m-1}\right)=s_{i,m-1}}\right\}}}}\\ f_{S_i({\tau }_1),\cdots S_i({\tau }_{j-1})}\left(s_{i,1},\cdots s_{i,j-1}\right)\text{d}{s}_{i,1}\cdots \text{d}{s}_{i,j-1}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_{j-1}}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}{\displaystyle {\int }_{-\infty }^{D_i}\cdots {\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t-{\displaystyle \underset{m=2}{\overset{j}{\sum }}{\tau }_{m-1}}\right)<D_i-{\displaystyle \underset{m=2}{\overset{j}{\sum }}{s}_{i,m-1}}\right\}}}}\\ {\displaystyle \underset{k=1}{\overset{j-1}{\prod }}{f}_{S_i({\tau }_k)}\left(s_{i,k}\right)}\text{d}{s}_{i,1}\cdots \text{d}{s}_{i,j-1},\end{array}$ (21)

其中 $P\left\{S_i\left(u\right)<D_i\right\}$， $P\left\{N\left(t\right)=n\right\}$ 的表达式参见公式(15)到(17)，由于各阶段退化过程独立，故在变点处的退化量的联合分布独立，因此

$f_{S_i({\tau }_1),\cdots S_i({\tau }_{j-1})}\left(s_{i,1},\cdots s_{i,j-1}\right)={\displaystyle \underset{k=1}{\overset{j-1}{\prod }}{f}_{S_i({\tau }_k)}\left(s_{i,k}\right)},$ (22)

$f_{S_i({\tau }_k)}\left(s_{i,k}\right)=\frac{1}{\sqrt{2\text{π}\left({\tau }_k{}^2{\sigma }_{i,k}^2+\left(n+1\right){\sigma }_y^2\right)}}\exp \left(-\frac{{\left(s_{i,k}-{\tau }_k{\mu }_{i,k}-\left(n+1\right){\mu }_y\right)}^2}{2\left({\tau }_1{}^2{\sigma }_{i,k}^2+\left(n+1\right){\sigma }_y^2\right)}\right).$ (23)

则该阶段的可靠度为：

$\begin{array}{c}{R}_j\left(t\right)={\displaystyle \underset{i=1}{\overset{2}{\sum }}{R}_{i,j}\left(t\right)-1+C\left(F_{1,j}\left(t\right),F_{2,j}\left(t\right);\theta \right)}\\ {\displaystyle \underset{i=1}{\overset{2}{\sum }}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_{j-1}}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}{\displaystyle {\int }_{-\infty }^{D_i}\cdots {\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t-{\displaystyle \underset{m=2}{\overset{j}{\sum }}{\tau }_{m-1}}\right)<D_i-{\displaystyle \underset{m=2}{\overset{j}{\sum }}{s}_{i,m-1}}\right\}}}}}\\ {\displaystyle \underset{k=1}{\overset{j-1}{\prod }}{f}_{S_i({\tau }_k)}\left(s_{i,k}\right)}\text{d}{s}_{i,1}\cdots \text{d}{s}_{i,j-1}-1+C\left(F_{1,j}\left(t\right),F_{2,j}\left(t\right);\theta \right),\end{array}$ (24)

其中 $F_{i,j}\left(t\right),i=1,2$ 的表达式为：

$\begin{array}{c}{F}_{i,j}\left(t\right)=1-{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\underset{u\to {\tau }_{j-1}}{\lim }P\left\{S_i\left(u\right)<D_i\right\}P\left\{N\left(t\right)=n\right\}{\displaystyle {\int }_{-\infty }^{D_i}\cdots {\displaystyle {\int }_{-\infty }^{D_i}P\left\{S_i\left(t-{\displaystyle \underset{m=2}{\overset{j}{\sum }}{\tau }_{m-1}}\right)<D_i-{\displaystyle \underset{m=2}{\overset{j}{\sum }}{s}_{i,m-1}}\right\}}}}\\ {\displaystyle \underset{k=1}{\overset{j-1}{\prod }}{f}_{S_i({\tau }_k)}\left(s_{i,k}\right)}\text{d}{s}_{i,1}\cdots \text{d}{s}_{i,j-1}.\end{array}$ (25)

综上可得仅考虑累积冲击过程的二元性能相关多阶段退化可靠度表示为：

$R_s\left(t\right)=R_1\left(t\right)I_{[{\tau }_{\text{0}},{\tau }_{\text{1}})}\left(t\right)+R_2\left(t\right)I_{[{\tau }_{\text{1}},{\tau }_{\text{2}})}\left(t\right)+\cdots +R_j\left(t\right)I_{[{\tau }_{j-\text{1}},{\tau }_j)}\left(t\right).$ (26)

当累积冲击过程引起的软失效和极值冲击引起的硬失效任意一种失效模式发生时，系统都会发生失效。因此，同时考虑系统的软失效和硬失效，此时的可靠度为：

$R\left(t\right)={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}P\left\{{\displaystyle \underset{i=1}{\overset{2}{\cap }}{S}_i\left(u\right)<D_i},{\displaystyle \underset{k=1}{\overset{N(t)}{\cap }}{W}_k<D_w}|N\left(t\right)=m\right\}P\left\{N\left(t\right)=m\right\}.}$ (27)

假设累积冲击过程的系统退化量 $Y_n$ 和极值冲击过程的性能退化量 $W_k$ 是独立的。在已知冲击次数的条件下，系统发生的两种失效模式是独立的。故此模型的可靠度表达式为：

$R\left(t\right)={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}P\left\{{\displaystyle \underset{i=1}{\overset{2}{\cap }}{S}_i\left(u\right)<D_i}|N\left(t\right)=m\right\}{\left(P\left\{{\displaystyle \underset{k=1}{\overset{N(t)}{\cap }}{W}_k<D_w}\right\}\right)}^{m}P\left\{N\left(t\right)=m\right\}},$ (28)

其中 $P\left\{{\displaystyle \underset{i=1}{\overset{2}{\cap }}{S}_i\left(u\right)<D_i}|N\left(t\right)=m\right\}$ 的表达式为 $R_s\left(t\right)$。结合公式(7)和公式(26)，公式(28)的表达式为：

$R\left(t\right)={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}{R}_s\left(t\right){\left(\Phi \left(\frac{D_w-{\mu }_w}{{\sigma }_w}\right)\right)}^m\frac{{\left(\lambda t\right)}^m{\text{e}}^{-\lambda t}}{m!}}.$ (29)

3.2. 随机变点情形下考虑竞争失效的二元多阶段随机系统的可靠度

本节将在固定变点情形下的可靠度模型基础上，建立变点随机情形下的可靠度模型。令 $T_j$ 表示退化过程中退化性能从第 $j$ 阶段到第 $j+1$ 阶段的随机变点的发生时刻，点过程 $\left\{N_r\left(t\right),t\ge 0\right\}$ 表示 $t$ 时刻之前的变点个数，是强度为 ${\lambda }_r$ 的Poisson过程。

如果 $t\in \left[0,T_1\right)$，表示 $t$ 时刻之前未产生变点，此时的系统可靠度为：

$R_{r,0}\left(t\right)=R\left(t\right)P\left\{N_r\left(t\right)=0\right\}.$ (30)

如果 $t\in \left[T_1,T_2\right)$，表示 $t$ 时刻之前存在一个变点，此时的系统可靠度为：

$R_{r,1}\left(t\right)={\displaystyle {\int }_0^{t}R\left(t\right)P\left\{N_r\left(t\right)=1\right\}f\left(t_1|N_r\left(t\right)=1\right)\text{d}{t}_1}.$ (31)

同理，如果 $t\in \left[T_m,T_{m+1}\right)$，表示 $t$ 时刻之前存在 $m$ 个变点，此时的系统可靠度为：

$R_{r,m}\left(t\right)={\displaystyle {\int }_{t_m}^t\cdots {\displaystyle {\int }_{t_2}^{t_3}{\displaystyle {\int }_0^{t_2}R\left(t\right)P\left\{N_r\left(t\right)=m\right\}f\left(t_1,t_2,\cdots ,t_m|N_r\left(t\right)=m\right)\text{d}{t}_1}\text{d}{t}_2}\cdots \text{d}{t}_m}.$ (32)

综上所述，可知在随机变点情形下的系统可靠度为：

$\begin{array}{c}{R}_r\left(t\right)=R_{r,0}\left(t\right)+R_{r,1}\left(t\right)+\cdots +R_{r,j}\left(t\right)=R_{r,0}\left(t\right)+{\displaystyle \underset{m=1}{\overset{j}{\sum }}{R}_{r,m}\left(t\right)}\\ =R\left(t\right)P\left\{N_r\left(t\right)=0\right\}+{\displaystyle \underset{m=1}{\overset{j}{\sum }}{\displaystyle {\int }_{t_{m\text{-}1}}^t\cdots {\displaystyle {\int }_{t_2}^{t_3}{\displaystyle {\int }_0^{t_2}R\left(t\right)P\left\{N_r\left(t\right)=m\right\}f\left(t_1,t_2,\cdots ,t_m|N_r\left(t\right)=m\right)\text{d}{t}_1}\text{d}{t}_2}\cdots \text{d}{t}_m}}\\ =R\left(t\right){\text{e}}^{-\lambda t}+{\displaystyle \underset{m=1}{\overset{j}{\sum }}{\displaystyle {\int }_{t_{m\text{-}1}}^t\cdots {\displaystyle {\int }_{t_2}^{t_3}{\displaystyle {\int }_0^{t_2}R\left(t\right){\lambda }^m{\text{e}}^{-\lambda t}\text{d}{t}_1}\text{d}{t}_2}\cdots \text{d}{t}_m}},\end{array}$ (33)

其中 $R\left(t\right)$ 的表达式参见公式(29)。 $f\left(\cdot |N\left(t\right)=m\right)$ 表示 $t$ 时刻之前存在 $m$ 个变点的条件下，各随机变点的联合条件概率密度函数。由文献 [22] 可知 $f\left(t_1,t_2,\cdots ,t_m|N\left(t\right)=m\right)=\frac{m!}{t^m}$。

4. 数值模拟

本节将通过数值算例来验证本文提出的考虑竞争失效的多阶段二元随机退化模型的有效性。使用蒙特卡洛方法产生相关的仿真样本数据，将仿真模拟结果与解析结果进行比对，最终分析判断模型的准确性。

首先需要仿真得到两组相关的性能退化数据，由二维Copula的一般形式

$H\left(u,v\right)=C\left(u,v\right),$ (34)

可知Copula函数对 $u$ 的偏导数 $C_u\left(v\right)$ 为

$C_u\left(v\right)=\frac{\partial C\left(u,v\right)}{\partial u}.$ (35)

令 $s=C_u\left(v\right)$，则 $v=C_u^{-1}\left(s\right)$。此时有 $C^{-1}:\left(u,s\right)\to v$，输入 $u$， $s$ 即可得到 $v$，这样得到的 $u$， $v$ 就是服从联合分布的 $C\left(u,v\right)$ 的随机变量。以Frank Copula函数为例 [23] ，可得

$C_u\left(v\right)=\frac{\exp \left(-\theta u\right)\left(\exp \left(-\theta v\right)-1\right)}{\left(\exp \left(-\theta u\right)-1\right)\left(\exp \left(-\theta v\right)-1\right)+\exp \left(-\theta \right)-1},$ (36)

$C_u^{-1}\left(s\right)=-\frac{1}{\theta }\ln \left(\frac{\left(1-s\right)\exp \left(-\theta u\right)+s\exp \left(-\theta \right)}{\left(1-s\right)\exp \left(-\theta u\right)+s}\right).$ (37)

给定退化量初始值，退化变点时刻，性能指标1的各阶段退化速率，迭代步长 $\Delta t$ 等参数，可得到性能指标1的退化量为：

$X_1=\left(X_{1,1}\left(0\right),X_{1,1}\left(\Delta t\right),\cdots ,X_{1,2}\left({\tau }_1\right),X_{1,2}\left({\tau }_1+\Delta t\right),\cdots X_{1,j}\left({\tau }_{j-1}\right),X_{1,j}\left({\tau }_{j-1}+\Delta t\right)\right).$ (38)

此时，性能指标1的退化增量为：

$\Delta X_1=\left(\Delta X_{1,1}\left(\Delta t\right),\cdots ,\Delta X_{1,2}\left({\tau }_1\right),\Delta X_{1,2}\left(\Delta t\right),\cdots \Delta X_{1,j}\left({\tau }_{j-1}\right),\Delta X_{1,j}\left(\Delta t\right)\right).$ (39)

而退化增量的分布函数表达式为：

$F_{\Delta X_{1,j}(t)}\left(x_{1,j}\right)=\Phi \left(\frac{\Delta X_{1,j}\left(t\right)-{\mu }_{1,j}\Delta t}{{\sigma }_{1,j}\sqrt{\Delta t}}\right)=u.$ (40)

将性能指标1的退化增量分布函数 $u$ 与服从(0, 1)上均匀分布的 $s$ 代入Frank Copula函数的偏导函数当中，便可得到性能指标2的退化增量分布函数 $v$。同理，可知性能指标2的分布函数为：

$F_{\Delta X_{2,j}(t)}\left(x_{2,j}\right)=\Phi \left(\frac{\Delta X_{2,j}\left(t\right)-{\mu }_{2,j}\Delta t}{{\sigma }_{2,j}\sqrt{\Delta t}}\right)=v.$ (41)

则性能指标2的退化增量表达式为：

${\mu }_{2,j}\Delta t+{\Phi }^{-1}\left(v\right){\sigma }_{2,j}\sqrt{\Delta t}=\Delta X_{2,j}\left(t\right),$ (42)

其中 ${\Phi }^{-1}\left(v\right)$ 、 $\Delta t$ 、 ${\mu }_{2,j}$ 和 ${\sigma }_{2,j}$ 已知，累加退化增量最终可得到性能指标2的退化量为：

$X_2=\left(X_{2,1}\left(0\right),X_{2,1}\left(0\right)+\Delta X_{2,1}\left(0\right),\cdots ,X_{2,j}\left({\tau }_{j-1}\right),X_{2,j}\left({\tau }_{j-1}\right)+\Delta X_{2,j}\left({\tau }_{j-1}\right)\right).$ (43)

模拟出系统设备的退化数据后，需要进一步模拟冲击过程。由于冲击过程为泊松过程，因此可以根据泊松过程的到达率得到在 $[t,t+\Delta t]$ 范围内冲击到达的时间分布，求和之后便可以得到 $[t,t+\Delta t]$ 发生的总冲击次数 $N\left(t+\Delta t\right)-N\left(t\right)$。定义初始时刻 $t=0,N\left(t\right)=0$，迭代相加后便可得到对应时刻发生的冲击次数 $m$。

随后模拟生成 $m$ 个服从一定分布的该时刻冲击退化量，若冲击过程为累积冲击过程，则将 $m$ 个冲击退化量相加求和。若冲击过程为极值冲击过程，则取 $m$ 个冲击退化量中的最大值。综上，系统可靠度的模拟仿真算法如下：

步骤1 (输入基本参数)：给定迭代步长 $\Delta t$，模拟次数 $N$，参数 ${\mu }_{i,j}$， ${\sigma }_{i,j}$， $\theta$， $\lambda$， ${\mu }_w$， ${\sigma }_w$， $D_w$， ${\mu }_y$， ${\sigma }_y$， $D_i$ 的数值，并初始化模拟次数 $n\left(=1,2,\cdots ,N\right)=1$，四个空矩阵 $X1$ 、 $X2$ 、 $Dw$ 和 $Dy$，分别用于记录第 $n$ 次模拟的两组性能指标的退化数据、极值冲击过程数据和累积冲击过程数据，其中矩阵的行数表示模拟次数，列数表示依序排列的时刻。

步骤2 (开始循环)：判断当前模拟次数 $n$ 与 $N$ 的大小，若 $n<N$，则初始化当前时刻 $t=0$，然后执行步骤3~4，否则进入步骤5。

步骤3 (模拟退化过程)：利用输入的参数 ${\mu }_{i,j}$， ${\sigma }_{i,j}$， $\theta$ 得到退化速率 ${\beta }_{i,j}$，得到第 $n$ 次模拟的两组性能指标的退化数据，并分别记录在矩阵 $X1$ 、 $X2$ 的第 $n$ 行当中。

步骤4 (模拟冲击过程)：初始化冲击次数 $m=0$，利用输入的参数 $\lambda$， ${\mu }_w$， ${\sigma }_w$， ${\mu }_y$， ${\sigma }_y$ 分别得到累积冲击过程的累积退化量和极值冲击过程的最大退化量，并分别记录在矩阵 $Dy$ 和 $Dw$ 的第 $n$ 行当中。然后令 $n=n+1$，返回执行步骤2。

步骤5 (判断是否失效)：将基于退化过程数据得到的矩阵 $X1$ 、 $X2$ 分别与基于累积冲击过程数据得到的矩阵 $Dy$ 相加并遍历相加后的两个矩阵 $X1+Dy$ 和 $X2+Dy$，判断对应的性能指标i相加后的矩阵中的元素是否超过阈值 $D_i$，若超出则记录对应位置元素为1，否则记录对应位置元素为0。同理，遍历基于极值冲击过程数据得到的矩阵 $Dw$，判断矩阵中的元素是否超过阈值 $D_w$，若超出则记录对应位置元素为1，否则记录对应位置元素为0。

步骤6 (计算失效概率)：随后，矩阵 $X1+Dy$ 、 $X2+Dy$ 和 $Dw$ 的相同位置元素做逻辑或运算，将得到的结果保存在矩阵当中。对矩阵进行列求和得到次模拟中不同时刻失效的次数，并分别与总模拟次数做比值，得到的比值向量即为系统的模拟失效概率。

本算例中，所有的参数值为：，，，，，，，，，，，。模拟次数，步长，变点时间，Copula函数为Frank Copula函数。表1给出了次模拟后的部分样本退化均值，图1给出了基于表1数据得到的不同情况下的退化路径图像及其累积退化增量分布直方图。

Table 1. Degraded sample data

表1. 退化样本数据

Figure 1. General degradation path and cumulative impact degradation path. (a) General degradation paths; (b) Considering the degradation path of cumulative impact surge; (c) The degenerate path function considering the cumulative impact surge; (d) Histogram of cumulative impact increment

图1. 一般退化路径与累积冲击退化路径。(a) 一般退化路径；(b) 考虑累积冲击突增的退化路径；(c) 考虑累积冲击突增的退化路径函数；(d) 累积冲击增量直方图

另外，利用解析表达式和模拟算法可以得到变点固定情况下的极值冲击和累积冲击的退化过程可靠度模拟解和解析解的曲线，如图2所示。从图2中可以看到基于解析解与模拟解的系统可靠度曲线基本一致。图3则给出了基于模拟解得到的不同冲击类型下的系统可靠度。

Figure 2. System reliability with fixed change point

图2. 变点固定情况下的系统可靠度

Figure 3. The reliability curves of different impact types with fixed change points

图3. 变点固定情况下不同冲击类型的可靠度曲线

图3中可以看到，不考虑冲击影响的系统可靠度曲线与仅考虑累积过程和同时考虑两种冲击过程的系统可靠度曲线的相对位置不同，同时考虑两种冲击的系统，其可靠度曲线呈下降趋势的时间更早，仅考虑累积的次之，没有冲击的可靠度曲线的呈下降趋势的时间则最晚。符合日常系统或产品的一般失效规律。

接下来，考虑变点随机情况下的系统可靠度曲线。假设系统存在两个变点，呈现三阶段退化过程。而变点产生的时间间隔服从参数为的指数分布。本算例中，，，，，和30，，，，，，，，步长，和1。得到的可靠度曲线如图4所示。

根据结果可知，随着值的增大，可靠度曲线呈现下降的趋势更晚，由于变点产生的时间间隔服从指数分布，值越大，间隔越大，故进入退化速率较快的第二、三阶段的时间较晚，因此系统的失效时间延后。另外，由图5可知，值越大，可靠度曲线呈现下降的趋势更晚，由于Copula函数的性质可知，当趋于0时，两个性能的可靠度趋向于独立，而独立的两个退化量会降低整体的可靠度 [24]。

Figure 4. The reliability curve of random change points

图4. 变点随机情况下的可靠度曲线

Figure 5. The reliability curve of random change points

图5. 变点随机情况下的可靠度曲线

5. 结论

本文针对考虑竞争失效的多阶段二元随机退化系统的可靠性建模问题，利用Copula连接函数的方法，建立了可靠性模型，得到了系统的可靠度解析表达式，通过数值算例验证模型的有效性。

在可靠度解析解的计算过程中，涉及大量繁琐的复杂多重积分计算，算法时间复杂度的降低是进一步的研究工作核心。但相较于可靠度模拟计算方法来说，本文给出的可靠度计算方法精度较高，运算次数较少。针对模型的参数估计问题和应用性评价是可靠性重要的研究环节，未来将以此模型为基础，进行退化失效数据的分析，利用对应的统计方法对模型进行参数估计。