-
摘要: 针对退化过程呈现两阶段特征的一类随机退化设备, 现有剩余寿命预测方法不适用于测量间隔分布不均匀、监测数据的测量频率与历史数据频率不一致的情况, 并且忽略了自适应漂移的可变性. 鉴于此, 提出了一种新的考虑个体差异性的两阶段自适应Wiener过程剩余寿命预测模型与方法. 首先, 基于自适应Wiener过程分阶段构建随机退化模型, 在首达时间意义下推导出寿命和剩余寿命解析式. 然后, 结合Kalman滤波技术和期望最大化算法进行参数自适应更新, 同时利用赤池信息准则实现退化模型变点的辨识. 最后, 通过蒙特卡洛仿真和锂电池实例, 验证了本文所提方法的有效性和实用价值.
-
关键词:
- 剩余寿命 /
- 两阶段自适应Wiener过程 /
- 期望最大化算法 /
- 赤池信息准则
Abstract: For a class of stochastic degradation equipments with two-stage degradation process, the existing remaining useful life (RUL) prediction method is not suitable for the situation that non-uniform measurement interval, inconsistent measurement frequency of monitoring data and historical data, and ignored it the variability of adaptive drift. In view of this, the paper proposes a new RUL prediction model and method which considers individual differences two-stage adaptive Wiener process. Firstly, the stochastic degradation model is constructed in stages based on the adaptive Wiener process, and the analytical expressions of life and RUL are derived in the sense of first arrival time. Then, the Kalman filtering technology and expectation maximization (EM) algorithm are combined to update the parameters adaptively, and the Schwarz information criterion is used to identify the change points of the degraded model. Finally, the effectiveness and practical value of the proposed method are verified by an example of lithium-ion battery data and the Monte Carlo simulation. -
随着高新技术的迅猛发展, 现代工业设备正朝着大型化、复杂化和智能化趋势快速发展. 这类设备在运行过程中由于受到内部和外部因素的随机影响, 性能和健康状态不可避免地呈现下降趋势乃至退化失效, 导致无法完成正常任务和功能, 进而引发严重事故, 造成环境破坏和人员伤亡[1-3]. 如果能在设备性能退化初期对其进行剩余寿命预测(Remaining useful life, RUL), 并基于预测结果确定维修决策的最佳时机, 制定相应的备件订购或替换策略, 将有效提高设备运行可靠性、降低运行成本. 近年来, 预测与健康管理技术(Prognostics and health management, PHM)得到广泛关注和应用, 已经成为可靠性领域的热点研究方向, 而PHM技术的关键在于预测运行设备的剩余寿命. 因此, 如何得到精确且符合实际情况的剩余寿命, 对切实保障系统的运行安全性、可靠性与经济性具有重要的意义[4-7].
经过几十年的发展, RUL预测取得了丰硕的理论成果并得到广泛应用. 袁烨等[8]将寿命预测研究主要分为模型方法和数据驱动方法. 数据驱动方法主要包括统计方法和机器学习方法, 基于设备大量退化数据推导退化模型, 进而判断超过失效阈值的时间预测剩余寿命, 如: 最小二乘方法、支持向量机方法和深度学习方法等. 模型方法分为物理退化模型(机理建模)和经验退化模型. 相比于物理退化模型, 经验退化模型能通过随机模型对监测数据建模, 进而得到寿命或剩余寿命的概率分布, 便于量化寿命或剩余寿命的不确定性, 从而为健康管理奠定基础, 更加适用于现代复杂工业设备. 而在经验退化模型中, Wiener过程和Gamma过程是两种最常用随机过程退化建模的方法. Gamma过程是一种增量非负的单调过程, 主要适用于单调退化过程, 如磨损过程、疲劳扩展过程. 但在工程实际中, 设备的退化过程大多为非单调, 如锂电池容量退化、惯性平台陀螺仪的退化等. Wiener过程作为非单调退化过程, 凭借良好的数学特性, 在RUL预测和健康管理领域得到广泛应用.
近年来, 大多数基于Wiener过程退化建模方法普遍假设系统在退化过程中是一种遵循单一阶段的随机退化模型. 但在工程实际中, 由于受到内部因素(如: 退化机理突变)或外部因素(如: 动态环境、状态切换等)的影响, 许多设备的退化特性呈现出两阶段乃至多阶段退化特征[9]. 例如锂电池[10]开始时经历一个平稳退化期, 随着充放电的进行, 固体电解质层在电极上的生长以及副反应导致的活性材料的损失, 导致锂电池容量在后一阶段迅速衰落; 液力耦合器[11]开始时经历一个快速退化期, 当到达某一时刻(变点)退化速度明显下降, 与之类似的还有半导体激光器[12]、等离子显示板[13-14]等.
对于这种存在变点、呈现两阶段退化特性的设备进行退化建模和RUL预测, 已有不少学者进行了研究和拓展. Ng[12]根据退化数据的两阶段特性, 提出一种基于单个变点独立增量的两阶段随机退化模型, 并采用期望最大化(Expectation maximization, EM)算法对模型参数进行估计. Yan等[15]基于两阶段Wiener过程模型对液力耦合器进行可靠性校验, 并根据赤池信息准则(Schwarz information criterion, SIC)对变点进行辨识. Chen等[16]改进两阶段线性对数模型来描述滚球轴承的分阶段退化过程, 并用贝叶斯方法更新模型参数进行寿命估计. Wang等[17]提出了一种两阶段退化模型用于轴承退化数据的建模, 在第一段假设处于健康状态, 在第二段结合卡尔曼滤波和EM算法进行RUL估计. Peng等[18]为了提高RUL预测的鲁棒性和效率, 开发了一种半解析预测模型, 该模型可以避免RUL预测的大幅度波动, 自动跟踪不同的退化阶段, 并自适应地更新超参数. Zhang等[19]在两阶段Wiener过程退化模型的框架下, 推导出基于首达时间意义的寿命分布, 该模型优势在于充分考虑并量化变点处退化量的不确定性同时能够推广至更具有一般性的多阶段退化模型中.
尽管两阶段以及多阶段退化模型已经取得了一些理论与实际应用成果, 但仍存在一些问题有待解决. 目前大多数两阶段退化模型(如: Zhang等[19])都是基于Wang等[20]所提出的一阶自回归模型进行建模, 但该模型存在三点不足: 1) 假设噪声项是独立且均匀分布, 并且仅适用于均匀测量间隔. 由于不是自动测量或根据某些设计方案进行测量等原因, 在工程实际中设备退化过程的测量间隔往往是不均匀的; 2) 当使用多组同类型退化设备的历史数据或先验信息估计模型未知参数时, 必须要求监测数据的测量频率与历史数据中的测量频率相同. 否则, 历史数据将不再适用; 3) 该模型退化建模存在一个潜在假设, 即在后一时刻估计的随机参数与前一时刻的随机参数的后验估计完全相等, 并且当该模型用于RUL预测时, 使用最新的监测值来更新漂移系数, 该漂移系数从最后监测点开始保持不变, 直到系统发生故障. 这意味着该模型假设可以根据实时监测数据自适应更新漂移系数, 但在未来的RUL预测中忽略这种自适应漂移可变性.
针对上述问题, 本文提出了一种基于自适应Wiener过程的两阶段退化模型, 突破测量间隔固定和采样频率一致的要求限制, 同时考虑对表征退化个体差异性的漂移系数实现自适应更新. 在首达时间意义下, 推导出两阶段自适应Wiener过程模型的RUL分布解析式, 结合EM算法和Kalman滤波技术对模型参数进行估计和更新, 并基于SIC实现退化变点辨识, 最后通过锂电池的实例研究验证了本文所提方法可有效实现两阶段退化设备的RUL预测.
1. 两阶段随机退化过程建模
建立两阶段退化模型主要针对退化过程表现为两阶段特征的设备, 且整个退化时间内存在一个变点, 变点前后的退化率呈现显著差异性.
1.1 两阶段线性Wiener过程退化建模
令
$X(t)$ 表示设备在运行中的退化过程, 常用的Wiener过程模型具有以下形式[21]$$X(t) = {x_0} + \lambda t + \sigma B(t)$$ (1) 其中,
$\lambda $ 和$\sigma $ 分别表示退化过程的漂移系数和扩散系数,$B(t)$ 是一个标准Brownian运动,${x_0}$ 为退化初值, 通常假设${x_0} = 0$ .基于上述模型和假设条件, 对于存在变点的随机退化系统可以建立两阶段Wiener过程退化模型:
$$X(t) = \left\{ \begin{aligned} &{x_0} + {\lambda _1}t + {\sigma _1}B(t),&{{ 0 < t}} \leq \tau \\ &{x_\tau } + {\lambda _2}(t - \tau ) + {\sigma _2}B(t - \tau ),&{{ t > }}\tau \qquad\end{aligned} \right.$$ (2) 其中,
${x_0}$ 表示退化过程初值,${x_\tau }$ 表示第二阶段退化初值即变点处的退化量,$\tau $ 表示变点发生时间;${\lambda _1}$ 和${\sigma _1}$ 分别表示第一阶段退化过程的漂移系数和扩散系数,${\lambda _2}$ 和${\sigma _2}$ 分别表示第二阶段退化过程的漂移系数和扩散系数.1.2 两阶段自适应Wiener过程退化建模
由于上述的Wiener过程模型存在测量间隔不均匀、测量频率不一致以及在RUL预测中没有利用实时监测数据自适应更新漂移系数三点不足. 因此, 本文考虑基于Zhai等[22]所提如下自适应Wiener过程模型建模:
$$\begin{split} & \lambda (t) = {\lambda _0} + kW({{t}}) \\ & X(t) = \int_0^t {\lambda (\tau )} {\rm{d}}\tau + \sigma B(t) \end{split} $$ (3) 其中,
$\lambda (t)$ 是一个遵循Wiener过程且随时间变化的漂移系数.${\lambda _0}$ 是初始漂移率,$k$ 是自适应漂移率的扩散系数,$W(t)$ 是一个独立于$B(t)$ 的标准Brownian运动.根据式(1)、(3), 可以推导出两者各自的离散状态空间模型:
$$\left[ {\begin{array}{*{20}{c}} {{\lambda _i}} \\ {{X_i}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&0 \\ {\Delta {t_i}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\lambda _{i - 1}}} \\ {{X_{i - 1}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\eta _i}} \\ {\sigma \Delta {B_i}} \end{array}} \right]$$ (4) $$\begin{split}\left[ {\begin{array}{*{20}{c}} {{\lambda _i}} \\ {{X_i}} \end{array}} \right] =\;& \left[ {\begin{array}{*{20}{c}} 1&0 \\ {\Delta {t_i}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\lambda _{i - 1}}} \\ {{X_{i - 1}}} \end{array}} \right] +\\ &\left[ {\begin{array}{*{20}{c}} {k\Delta {W_i}} \\ {k\displaystyle\int_{{t_{i - 1}}}^{{t_i}} {W(\tau ){\rm{d}}\tau } + \sigma \Delta {B_i}} \end{array}} \right] \end{split}$$ (5) 其中,
${\eta _i}\sim {\rm{N}}(0,Q)$ ,$\Delta {B_i} = B({t_i}) - B({t_{i - 1}})$ ,$\Delta {W_i} = $ $ W({t_i}) - W({t_{i - 1}})$ .从上式中, 可以看到自适应Wiener过程离散模型的噪声项比式(4)的Wiener过程离散模型, 增加了一个自适应漂移项
$k\displaystyle\int_{{t_{i - 1}}}^{{t_i}} {W(\tau ){\rm{d}}\tau }$ , 当使用最新的监测值${X_i}$ 更新漂移系数时, 自适应Wiener过程的漂移项从最后监测点开始仍可以动态变化, 直到系统发生故障. 此外, 如式(4)所示的Wiener过程模型假定两个监测点之间的漂移系数${\lambda _i}$ 为一个随机游走模型, 并且依赖于前一个时刻的漂移系数, 当测量间隔发生变化时, 可能导致模型参数估计不准确. 鉴于此, 自适应Wiener过程增加了一个随时间变化的漂移项. 当测量间隔发生变化时, 模型漂移部分能动态相应变化, 适用于不等间隔测量下的RUL预测.为了进一步说明, 本文给出式(1)和式(3)的两步预测模型:
$$\begin{split} \left[ {\begin{array}{*{20}{c}} {{\lambda _i}} \\ {{X_i}} \end{array}} \right] =\;& \left[ {\begin{array}{*{20}{c}} 1&0 \\ {\Delta {t_i}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\lambda _{i - 1}}} \\ {{X_{i - 1}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\eta _i}} \\ {\sigma \Delta {B_i}} \end{array}} \right] =\\ &\left[ {\begin{array}{*{20}{c}} 1&0 \\ {\Delta {t_i} + \Delta {t_{i - 1}}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\lambda _{i - 2}}} \\ {{X_{i - 2}}} \end{array}} \right] +\\ &\left[ {\begin{array}{*{20}{c}} {{\eta _i} + {\eta _{i - 1}}} \\ {{\eta _{i - 1}}\Delta {t_i} + \sigma (\Delta {B_i} + \Delta {B_{i - 1}})} \end{array}} \right] \\[-15pt]\end{split} $$ (6) $$\begin{split} \left[ {\begin{array}{*{20}{c}} {{\lambda _i}} \\ {{X_i}} \end{array}} \right] =\;& \left[ {\begin{array}{*{20}{c}} 1&0 \\ {\Delta {t_i}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\lambda _{i - 1}}} \\ {{X_{i - 1}}} \end{array}} \right] + \\ &\left[ {\begin{array}{*{20}{c}} {k\Delta {W_i}} \\ {k\displaystyle\int_{{t_{i - 1}}}^{{t_i}} {W(\tau ){\rm{d}}\tau } + \sigma \Delta {B_i}} \end{array}} \right] = \\ &\left[ {\begin{array}{*{20}{c}} 1&0 \\ {\Delta {t_i} + \Delta {t_{i - 1}}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\lambda _{i - 2}}} \\ {{X_{i - 2}}} \end{array}} \right] + \\ &\left[ {\begin{array}{*{20}{c}} {k(\Delta {W_i} + \Delta {W_{i - 1}})} \\ \begin{array}{l} k\Delta {W_{i - 1}}\Delta {t_i} + k\displaystyle\int_{{t_{i - 2}}}^{{t_i}} {W(\tau ){\rm{d}}\tau }+ \\ \sigma (\Delta {B_i} + \Delta {B_{i - 1}}) \end{array} \end{array}} \right] \approx \\ &\left[ {\begin{array}{*{20}{c}} 1&0 \\ {\Delta {t_i} + \Delta {t_{i - 1}}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\lambda _{i - 2}}} \\ {{X_{i - 2}}} \end{array}} \right] + \\ &\left[ {\begin{array}{*{20}{c}} {k(\Delta {W_i} + \Delta {W_{i - 1}})} \\ {k\displaystyle\int_{{t_{i - 2}}}^{{t_i}} {W(\tau ){\rm{d}}\tau } + \sigma (\Delta {B_i} + \Delta {B_{i - 1}})} \end{array}} \right] \end{split} $$ (7) 从式(6)中观察到, Wiener过程两步预测模型的噪声项中第一个元素
${\eta _i} + {\eta _{i - 1}}$ 的方差是一步预测模型的2倍, 然而两步预测模型的噪声项中第二个元素除了$\sigma (\Delta {B_i} + \Delta {B_{i - 1}})$ 还有一个附加项${\eta _{i - 1}}\Delta {t_i}$ , 与一步预测模型的噪声项第二个元素相比, 方差不能构成2倍关系. 因此, 当式(1)所示的Wiener过程应用于测量间隔不均匀的数据时, 可能造成估计值不准确, RUL预测准确度随之下降. 相比之下, 自适应Wiener过程的两步预测模型增加了一个漂移项$k\displaystyle\int_{{t_{i - 2}}}^{{t_i}} {W(\tau ){\rm{d}}\tau }$ , 导致附加项$k\Delta {W_{i - 1}}\Delta {t_i}$ 的影响可忽略不计内容. 因此, 自适应Wiener过程两步预测与一步预测可以相互兼容, 当测量间隔变化时能解决此类问题带来的影响.基于上述结论, 本文对存在变点的随机退化系统建立两阶段自适应Wiener过程退化模型:
$$X(t) = \left\{ \begin{aligned} &{x_0} + \int_0^t {{\lambda _1}\left( s \right)} {\rm{d}}s + {\sigma _1}B(t),\;\;\qquad\quad{{0 < t}} \le \tau \\ &{x_\tau } + \int_0^{t - \tau } {{\lambda _2}\left( s \right)} {\rm{d}}s + {\sigma _2}B(t - \tau ),\;\;{{t > }}\tau \end{aligned} \right.$$ (8) 其中,
${x_0}$ 表示退化初值,${x_\tau }$ 表示第二阶段退化初值即变点处的退化量,$\tau $ 表示变点发生时间;${\lambda _1}(s)$ 和${\sigma _1}$ 分别表示第一阶段退化过程的漂移系数和扩散系数,${\lambda _2}(s)$ 和${\sigma _2}$ 分别表示第二阶段退化过程的漂移系数和扩散系数.2. 两阶段自适应Wiener过程剩余寿命预测
为描述同批次设备中某一个体的退化过程, 体现个体差异性, 将退化模型的漂移系数随机化[23], 即
${\lambda _1} \sim {\rm{N}}({u_a},\;\sigma _a^2)$ 和${\lambda _{\rm{2}}} \sim {\rm{N}}({u_b},\;\sigma _b^2)$ . 若发生变点时的退化量已知, 根据文献[22]中自适应Wiener过程的寿命分布, 可以推导出两阶段自适应Wiener过程寿命分布的概率密度函数(Probability density function, PDF), 如式(9)所示.其中,
${X_0}$ 表示退化初值,$D$ 表示设备退化的失效阈值,$\tau $ 表示变点发生时间.实际中, 在变点出现前, 变点处退化量准确值
${X_\tau }$ 是未知的, 为了得到寿命估计值, 首先要得到首达时间意义下${X_\tau }$ 的分布形式, 即在${X_\tau } < D$ 条件下经过时间$\tau $ , 退化量从0到${X_\tau }$ 的转移概率${g_\tau }({X_\tau })$ . 因此, 要计算退化过程在$ ({X}_{\tau },\;\infty )$ 失效概率, 需保证退化过程在$ (0,\;{X}_{\tau })$ 上未超过失效阈值. 如果${g_\tau }({X_\tau })$ 的解析表达式可以得到, 则基于全概率公式可推导出寿命分布的PDF.引理 1[19, 22]. 若退化过程为两阶段自适应Wiener过程模型, 且漂移系数
${\lambda _1}$ 、${\lambda _2}$ 服从正态分布${\lambda _1} \sim {\rm{N}}({u_a},\;\sigma _a^2)$ 、${\lambda _{\rm{2}}} \sim {\rm{N}}({u_b},\;\sigma _b^2)$ . 如果变点时间$\tau $ 给定, 那么在首达时间意义下的寿命分布的PDF, 如式(10)、(11)所示.$$\begin{split} & {f_T}(t) = \left\{ \begin{aligned} & \frac{1}{{\sqrt {2\pi ({\varphi _1}(t) + \sigma _a^2{t^2})} }}\left[ {\frac{{{\varphi _1}'(t)}}{{{\varphi _1}(t)}}(D - {X_0}) + \left(1 - \frac{{{\varphi _1}'(t)}}{{{\varphi _1}(t)}}t\right)\frac{{{u_a}{\varphi _1}(t) + (D - {X_0})\sigma _a^2t}}{{{\varphi _1}(t) + \sigma _a^2{t^2}}}} \right] \cdot\\ &\qquad\exp \left\{ { - \frac{{{{(D - {X_0} - {u_a}t)}^2}}}{{2({\varphi _1}(t) + \sigma _a^2{t^2})}}} \right\},\qquad\qquad\qquad{{ 0 < t}} \leq \tau \\ &\frac{1}{{\sqrt {2\pi ({\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2})} }}\left[ \begin{aligned} &\frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(D - {X_\tau }) + \left(1 - \frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(t - \tau )\right) \\ &\frac{{{u_b}{\varphi _2}(t - \tau ) + (D - {X_\tau })\sigma _b^2(t - \tau )}}{{{\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2}}} \end{aligned} \right] \cdot\\ &\qquad\exp \left\{ { - \frac{{{{(D - {X_\tau } - {u_b}(t - \tau ))}^2}}}{{2({\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2})}}} \right\},\;\;\qquad{{ t > }}\tau \end{aligned} \right. \\ &{\varphi _1}(t)=\frac{{\rm{1}}}{{\rm{3}}}k_1^2{t^3} + \sigma _1^2t,{\rm{ }}{\varphi _1}'(t)=k_1^2{t^2} + \sigma _1^2,\;{\varphi _2}(t - \tau )=\frac{{\rm{1}}}{{\rm{3}}}k_2^2{(t - \tau )^3} + \sigma _2^2(t - \tau ),{\rm{ }}{\varphi _2}'(t - \tau )=k_2^2{(t - \tau )^2} + \sigma _2^2 \end{split} $$ (9) $${f_T}(t) = \left\{ \begin{aligned} & \frac{1}{{\sqrt {2\pi ({\varphi _1}(t) + \sigma _a^2{t^2})} }}\left[ \begin{aligned} &\frac{{{\varphi _1}'(t)}}{{{\varphi _1}(t)}}(D - {X_0}) + \left(1 - \frac{{{\varphi _1}'(t)}}{{{\varphi _1}(t)}}t\right) \\ & \frac{{{u_a}{\varphi _1}(t) + (D - {X_0})\sigma _a^2t}}{{{\varphi _1}(t) + \sigma _a^2{t^2}}} \end{aligned} \right] \cdot \exp \left\{ { - \frac{{{{(D - {X_0} - {u_a}t)}^2}}}{{2({\varphi _1}(t) + \sigma _a^2{t^2})}}} \right\}, \;\;\;\;\;\; {{0 < t}} \leq \tau \\ &\int_{ - \infty }^D \frac{1}{{\sqrt {2\pi ({\varphi _2}(t - \tau ) + \sigma _2^2{{(t - \tau )}^2})} }}\left[ \begin{aligned} &\frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(D - {X_\tau }) + \left(1 - \frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(t - \tau )\right) \\ & \frac{{{u_b}{\varphi _2}(t - \tau ) + (D - {X_\tau })\sigma _b^2(t - \tau )}}{{{\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2}}} \end{aligned} \right]\cdot\\ &\qquad\exp \left\{ { - \frac{{{{(D - {X_\tau } - {u_b}(t - \tau ))}^2}}}{{2({\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2})}}} \right\} \cdot {g_\tau }({X_\tau }|{{{u}}_a},{\sigma _a}){\rm{d}}{X_\tau } \approx {{{A}}_1}{\rm{ - }}{{{B}}_1},{\rm{ }} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\qquad\qquad {{ t > }}\tau \end{aligned} \right.$$ (10) 推论 1. 若已知当前时刻
${t_k}$ 的退化状态${x_k}$ , 用${l_k}$ 表示设备剩余寿命,${f_L}({l_k})$ 表示设备RUL分布的PDF, 在随机退化速率${\lambda _1}$ 和${\lambda _2}$ 的影响下, 可获得首达意义下两阶段自适应Wiener过程模型RUL的PDF, 其形式与首达意义下得到的寿命分布PDF, 即与式(10)、(11)类似, 具体可分为以下两种情况:情况 1. 当前时刻
${t_k}$ 位于变点前, 即${t_k} < \tau $ , 此时随机设备退化失效又存在两种情况: 1) 失效阈值位于变点前, 即${l_k} + {t_k} \leq \tau$ ; 2) 失效阈值位于变点后, 即${l_k} + {t_k} > \tau $ , 此时RUL的PDF, 如式(12)所示.情况 2. 当前时刻
${t_k}$ 位于变点后, 即${t_k} > \tau $ , 此时退化设备RUL的PDF为$$\begin{split} &{A_1} = \frac{{\varphi' _2(t - \tau )}}{{{\varphi _2}(t - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a1}^2 + \sigma _{b1}^2)}}} \exp \left[ { - \frac{{{{({u_{a1}} - {u_{b1}})}^2}}}{{2(\sigma _{a1}^2 + \sigma _{b1}^2)}}} \right] \cdot \left\{\frac{{{u_{b1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}\Phi \left( {\frac{{{u_{b1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right)+ \right.\\ &\qquad\left. \sqrt {\frac{{\sigma _{a1}^2\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}} \phi \left( {\frac{{{u_{b1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right) \right\} \\ &{B_1} = \exp \left\{ {\frac{{2{u_a}D}}{{\varphi' _1(\tau )}} + \frac{{2({D^2}\sigma _a^4\tau + {D^2}\sigma _a^2\varphi' _1(\tau ))}}{{(\varphi' _1(\tau ) + \sigma _a^2\tau )\varphi' _1{{(\tau )}^2}}}} \right\}\frac{{\varphi' _2(t - \tau )}}{{{\varphi _2}(t - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a1}^2 + \sigma _{b1}^2)}}} \exp \left[ { - \frac{{{{({u_{a1}} - {u_{c1}})}^2}}}{{2(\sigma _{a1}^2 + \sigma _{b1}^2)}}} \right] \cdot\\ &\qquad\left\{ {\frac{{{u_{c1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}\Phi \left( {\frac{{{u_{c1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right) + \sqrt {\frac{{\sigma _{a1}^2\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}} \phi \left( {\frac{{{u_{c1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right)} \right\} \\ &{u_{a1}} = {u_b}(t - \tau ),\;\;{\rm{ }}{u_{b1}} = D - {u_a}\tau ,\;\;{\rm{ }}{u_{c1}} = - D - {u_a}\tau - \frac{{\sigma _a^2\tau }}{{\varphi _1^{'}(\tau )}} \\ & \sigma _{a1}^2 = {\varphi _2}(t - \tau ) + \sigma _b^2{(t - \tau )^2},\;\;{\rm{ }}\sigma _{b1}^2 = {\varphi _1}(\tau ) + \sigma _a^2{\tau ^2} \end{split} $$ (11) $$\begin{split} &{f_L}({l_k}) = \left\{ \begin{aligned} &\frac{1}{{\sqrt {2\pi ({\varphi _{\rm{1}}}({l_k}) + \sigma _a^2{l_k}^2)} }}\left[ \frac{{{\varphi' _{\rm{1}}}({l_k})}}{{{\varphi _{\rm{1}}}({l_k})}}(D - {X_k}) + \left(1 - \frac{{{\varphi' _{\rm{1}}}({l_k})}}{{{\varphi _{\rm{1}}}({l_k})}}{l_k}\right) \frac{{{u_a}{\varphi _{\rm{1}}}({l_k}) + (D - {X_k})\sigma _a^2{l_k}}}{{{\varphi _{\rm{1}}}({l_k}) + \sigma _a^2{l_k}^2}} \right] \cdot\\ &\qquad\exp \left\{ { - \frac{{{{(D - {X_k} - {u_a}{l_k})}^2}}}{{2({\varphi _{\rm{1}}}({l_k}) + \sigma _a^2{l_k}^2)}}} \right\},\;\;\;\;\;\, {\rm{ }}{l_k} + {t_k} \leq \tau \; \\ &{{{A}}_2}{\rm{ - }}{{{B}}_2}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\quad{\rm{ }}{l_k} + {t_k} > \tau \end{aligned} \right. \\ &{A_2} = \frac{{\varphi' _2({l_k} + {t_k} - \tau )}}{{{\varphi _2}({l_k} + {t_k} - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a2}^2 + \sigma _{b2}^2)}}} \cdot \exp \left[ { - \frac{{{{({u_{a2}} - {u_{b2}})}^2}}}{{2(\sigma _{a2}^2 + \sigma _{b2}^2)}}} \right] \cdot \left\{ \frac{{{u_{b2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}\Phi \left( {\frac{{{u_{b2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right)+\right. \\ & \qquad\left.\sqrt {\frac{{\sigma _{a2}^2\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}} \phi \left( {\frac{{{u_{b2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right) \right\} \\ & {B_2} = \exp \left\{ {\frac{{2{u_a}(D - {X_k})}}{{\varphi _1^{'}(\tau - {t_k})}} + \frac{{2({{(D - {X_k})}^2}\sigma _a^4\tau + {{(D - {X_k})}^2}\sigma _a^2\varphi' _1(\tau - {t_k}))}}{{(\varphi' _1(\tau - {t_k}) + \sigma _a^2(\tau - {t_k}))\varphi' _1{{(\tau - {t_k})}^2}}}} \right\} \cdot \frac{{\varphi' _2({l_k} + {t_k} - \tau )}}{{{\varphi _2}({l_k} + {t_k} - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a2}^2 + \sigma _{b2}^2)}}}. \\ &\;\;\qquad\exp \left[ { - \frac{{{{({u_{a2}} - {u_{c2}})}^2}}}{{2(\sigma _{a2}^2 + \sigma _{b2}^2)}}} \right] \cdot \left\{ \frac{{{u_{c2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}\Phi \left( {\frac{{{u_{c2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right) +\right.\\ &\;\;\qquad\left.\sqrt {\frac{{\sigma _{a2}^2\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}} \phi \left( {\frac{{{u_{c1}}\sigma _{a2}^2 + {u_{a1}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right) \right\} \\ &{u_{a2}} = {u_b}({l_k} + {t_k} - \tau ),\;\;\;\;{\rm{ }}{u_{b2}} = D - {X_k} - {u_a}(\tau - {t_k}),\;\;\;\;{\rm{ }}{u_{c2}} = - D + {X_k} - {u_a}(\tau - {t_k}) - \frac{{\sigma _a^2(\tau - {t_k})}}{{\varphi' _1(\tau - {t_k})}} \\ &\sigma _{a2}^2 = {\varphi _2}({l_k} + {t_k} - \tau ) + \sigma _b^2{({l_k} + {t_k} - \tau )^2},\;\;\;\;{\rm{ }}\sigma _{b2}^2 = {\varphi _1}(\tau - {t_k}) + \sigma _a^2{(\tau - {t_k})^2} \end{split} $$ (12) $$\begin{split} {f_L}({l_k}) =\;& \frac{1}{{\sqrt {2\pi ({\varphi _{\rm{2}}}({l_k}) + \sigma _b^2{l_k^2})} }} .\\ &\left[ \begin{aligned} &\frac{{{\varphi' _{\rm{2}}}({l_k})}}{{{\varphi _{\rm{2}}}({l_k})}}(D - {X_k}) + \left(1 - \frac{{{\varphi _{\rm{2}}}'({l_k})}}{{{\varphi _{\rm{2}}}({l_k})}}{l_k}\right) \\ & \frac{{{u_b}{\varphi _{\rm{2}}}({l_k}) + (D - {X_k})\sigma _b^2{l_k}}}{{{\varphi _{\rm{2}}}({l_k}) + \sigma _b^2{l_k}^2}} \end{aligned} \right] \cdot\\ &\exp \left\{ { - \frac{{{{(D - {X_k} - {u_b}{l_k})}^2}}}{{2({\varphi _{\rm{2}}}({l_k}) + \sigma _b^2{l_k^2})}}} \right\}\\[-15pt] \end{split} $$ (13) 在推论1中, 变点发生时间为某一固定值, 只适用于事先预设情况. 在实际中, 利用监测数据进行预测时, 通常两阶段变点位置在不同情况或不同个体下存在差异性. 因此, 假设变点时间
$\tau $ 为随机变量来描述这种差异性. 在这种情况下, 随机退化设备的寿命和剩余寿命PDF为$$\left\{ \begin{aligned} & {f_T}(t) = \int_0^{ + \infty } {{f_T}(t|\tau )p(\tau ){\rm{d}}\tau } \\ & {f_L}({l_k}) = \int_{{t_k}}^{ + \infty } {{f_L}({l_k}|\tau )p(\tau ){\rm{d}}\tau } \\[-10pt]\end{aligned} \right. $$ (14) 其中,
$p(\tau )$ 为变点发生时间的PDF. 由于寿命与剩余寿命分布形式比较复杂, 上述的积分难以得到具体解析表达式, 故本文考虑采用数值积分方法求解.3. 模型参数估计与更新
3.1 潜在退化状态估计
当前时间可定义为
${t_k}$ , 而当前运行设备从时间${t_0}\sim {t_k}$ 获取的退化数据为${{\boldsymbol{x}}_{0:k}} = \{ {x_0},{x_1}, \cdots ,{x_k}\} $ , 如果此时变点未出现, 即${t_k} \leq \tau$ , 也就是说退化设备仍处于第一阶段且尚无第二阶段退化数据, 那么仅需根据收集的退化数据来更新第一阶段模型参数; 反之, 若变点已经出现, 即${t_k} > \tau $ , 那么仅需要更新第二阶段模型参数.根据建立的两阶段模型, 可将漂移系数
${\lambda _{\rm{1}}}$ 和${\lambda _{\rm{2}}}$ 视作“隐含状态”, 因此基于实时监测数据${{\boldsymbol{x}}_{0:k}}$ , 可以利用Kalman滤波进行状态估计. 在此, 定义${\lambda _{\rm{1}}}$ 和${\lambda _{\rm{2}}}$ 的均值和方差分别为${\widehat \lambda _1} = {\rm{E}}({\lambda _1}|{{\boldsymbol{x}}_{0:k}})$ 、${\widehat \lambda _2} = $ $ {\rm{E}}({\lambda _2}|{{\boldsymbol{x}}_{\tau + 1:k}})$ 和${P_1}_{k|k} = {\rm{var}} ({\lambda _1}|{{\boldsymbol{x}}_{0:k}})$ 、${P_2}_{k|k} = {\rm{var}} ({\lambda _2}| $ $ {{\boldsymbol{x}}_{\tau + 1:k}})$ .算法1. Kalman滤波算法
第一阶段(
${t_k} \leq \tau$ ):初始化
$${\widehat \lambda _{_{{\rm{10}}}}} = {a_{10}},{P_{10}}$$ 状态估计
$$\begin{split} &{P_1}_{k|k - 1} = {P_1}_{k - 1|k - 1} + k_1^2\Delta {t_k} \\ &{K_1}(k) = {P_1}_{k|k - 1}\Delta {t_k}{({P_1}_{k|k - 1}\Delta t_k^2 + {\varphi _1}(\Delta {t_k}))^{ - 1}} \\ &{\widehat \lambda _{1k}} = {\widehat \lambda _{1k - 1}} + {K_1}(k)({x_k} - {x_{k - 1}} - {\widehat \lambda _{1k - 1}}\Delta {t_k}) \end{split} $$ 方差更新
$${P_1}_{k|k} = {P_1}_{k|k - 1} + {P_1}_{k|k - 1}{K_1}(k)\Delta {t_k}$$ 类似地, 若
${t_k} > \tau $ , 可利用当前运行设备退化数据更新参数${\lambda _{\rm{2}}}$ , 由于第一阶段数据与第二阶段模型无关, 因此仅需要数据${{\boldsymbol{x}}_{{\boldsymbol{\tau}} :{\boldsymbol{k}}}} = \{ {x_{\tau + 1}},{x_{\tau + 2}}, \cdots ,{x_k}\} $ 用于更新.第二阶段(
${t_k} > \tau $ ):初始化
$${\widehat \lambda _{_{{\rm{20}}}}} = {a_{20}},{P_{20}}$$ 状态估计
$$\begin{split} &{P_2}_{k|k - 1} = {P_2}_{k - 1|k - 1} + k_2^2\Delta {t_k} \\ &{K_2}(k) = {P_2}_{k|k - 1}\Delta t{({P_2}_{k|k - 1}\Delta t_k^2 + {\varphi _2}(\Delta {t_k}))^{ - 1}} \\ & {\widehat \lambda _{2k}} = {\widehat \lambda _{2k - 1}} + {K_2}(k)({x_k} - {x_{k - 1}} - {\widehat \lambda _{2k - 1}}\Delta {t_k}) \end{split} $$ 方差更新
$${P_2}_{k|k} = {P_2}_{k|k - 1} + {P_2}_{k|k - 1}{K_2}(k)\Delta {t_k}$$ 当两阶段自适应Wiener模型用于RUL预测时, 模型参数
${a_{10}},{a_{20}},{P_{10}},{P_{20}}$ ,$k_1^2,k_2^2,\sigma _1^2,\sigma _2^2$ 均是未知的. 对此, 本文采用EM算法对参数自适应估计, 使得估计的寿命更好地反映设备当前健康状态.3.2 基于EM算法的自适应估计
假设对某个退化设备进行状态监测, 监测点为
$m$ 个, 即${\boldsymbol{x}} = \{ {x_1},{x_2}, \cdots ,{x_m}\} $ , 其各自对应的监测时间为$\{ {t_1},{t_2}, \cdots ,{t_m}\} $ . 同时, 本文假设变点发生时间已知, 即$\tau \in \{ {t_1},{t_2}, \cdots ,{t_m}\} $ , 那么$\{ {x_1},{x_2}, \cdots , {x_\tau }\}$ 表示设备第一阶段的退化数据,$\{ {x_{\tau + 1}},{x_{\tau + 2}}, \cdots , $ $ {x_m}\}$ 表示设备第二阶段的退化数据.首先, 令
${{\boldsymbol{\theta }}_1} = {({a_{10}},{P_{10}},k_1^2,\sigma _1^2)^{\rm{T}}}$ 表示第一阶段未知参数向量, 因此在未知参数${{\boldsymbol{\theta }}_1}$ 条件下${t_k}$ 时刻的监测数据${{\boldsymbol{x}}_{0:k}}$ 的对数似然函数为$${L_{1k}}({{\boldsymbol{\theta }}_1}) = \ln p({{\boldsymbol{x}}_{0:k}}|{{\boldsymbol{\theta }}_1})$$ (15) 式中,
$p({{\boldsymbol{x}}_{0:k}}|{{\boldsymbol{\theta }}_1})$ 为监测数据${{\boldsymbol{x}}_{0:k}}$ 的联合PDF.然后, 基于监测数据
${{\boldsymbol{x}}_{0:k}}$ 的似然函数,${{\boldsymbol{\theta }}_1}$ 的极大似然估计值${\widehat {\boldsymbol{\theta }}_{_{1k}}}$ 可通过最大化似然函数得到, 表示为$${\widehat {\boldsymbol{\theta }}_{1k}} = \arg \;\mathop {\max }\limits_{{{\boldsymbol{\theta }}_{\boldsymbol{1}}}} {L_{1k}}({{\boldsymbol{\theta }}_1})$$ (16) 在本文中, 由于漂移系数
${\lambda _1}$ 被视作一个“隐含状态”. 无法使未知参数${{\boldsymbol{\theta }}_1}$ 最大化. 而EM算法可通过最大化联合似然函数$p({\lambda _{1k}},{{\boldsymbol{x}}_{0:k}}|{{\boldsymbol{\theta }}_1})$ 来估计逼近参数的极大似然估计, 估计值可以通过迭代以下两步实现.1) E-步骤
$$\ell ({{\boldsymbol{\theta }}_1}|\widehat {\boldsymbol{\theta }}_{1k}^{(i)}) = {{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}},\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}\{ p({\lambda _{1k}},{{\boldsymbol{x}}_{0:k}}|{{\boldsymbol{\theta }}_1})\} $$ (17) 其中,
$\widehat {\boldsymbol{\theta }}_{1k}^{(i)} = [a_{10k}^{(i)},P_{10k}^{(i)},k_{{1}k}^{{2}(i)},\sigma _{{1}k}^{{2}(i)}]$ 表示基于监测数据${{\boldsymbol{x}}_{0:k}}$ 在第$i$ 步估计的参数值.2) M-步骤
$$\widehat {\boldsymbol{\theta }}_{1k}^{(i)} = \arg \;\mathop {\max }\limits_{{{\boldsymbol{\theta }}_{{1}}}} \{ \ell ({{\boldsymbol{\theta }}_{{1}}}|\widehat {\boldsymbol{\theta }}_{1k}^{(i)})\} $$ (18) 通过不断迭代E-步骤和M-步骤直到满足某一收敛条件截止, 由此得到对应的参数估计值, 一般来说随着迭代次数增加, 得到的参数估计值会越来越好.
对于第一阶段的随机退化模型, EM算法中的联合对数似然函数可以表示为
$$\begin{split} {\ell _{1k}}({{\boldsymbol{\theta }}_1}) =\;& \ln p({\lambda _{10}}|{{\boldsymbol{\theta }}_1}) + \ln \prod\limits_{j = 1}^k {p(} {\lambda _{1j}}|{\lambda _{1j - 1}},{{\boldsymbol{\theta }}_1}) +\\ &\ln \prod\limits_{j = 1}^k {p({x_j}|{\lambda _{1j - 1}},{{\boldsymbol{\theta }}_1})}\\[-15pt] \end{split} $$ (19) 根据线性状态空间模型, 有
$$\begin{split} & {\lambda _{1k}}|{\lambda _{1k - 1}} \sim {\rm{N}}({\lambda _{1k - 1}},{k^2}_1\Delta {t_k}) \\ & {x_k}|{\lambda _{1k - 1}} \sim {\rm{N}}({x_{k - 1}} + {\lambda _{1k - 1}}\Delta {t_k},{\varphi _1}(\Delta {t_k}) \\ & {\lambda _{10}} \sim {\rm{N}}({a_{10}},{P_{10}}) \end{split} $$ (20) 因此, 由E-步骤计算
${\ell _{1k}}({\boldsymbol{\theta }}_1)$ 的条件期望, 即$\ell ({\boldsymbol{\theta }}_1|\widehat {\boldsymbol{\theta }}_{1k}^{(i)}) .$ $$\begin{split} &2\ell ({\boldsymbol{\theta }}_1|\widehat {\boldsymbol{\theta }}_{1k}^{(i)}) = {{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}},\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}[2{\ell _{1k}}({\boldsymbol{\theta }}_1)] = {{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}},\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}\\ &\left[ \begin{array}{l} - \ln {P_{10}} - {({\lambda _{10}} - {a_{10}})^2}/{P_{10}} \\ - \sum\limits_{j = 1}^k {(\ln k_1^2\Delta {t_j} + (} {\lambda _{1j}} - {\lambda _{1j - 1}})k_1^2\Delta {t_j}) \\ - \sum\limits_{j = 1}^k {(\ln \sigma _1^2} + {({x_j} - {x_{j - 1}} - {\lambda _{1j - 1}}\Delta {t_j})^2}/\sigma _1^2\Delta {t_j}) \end{array} \right] \end{split} $$ (21) 显然, 要想计算
${{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}},\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}({\lambda _{1j}})$ 、${{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}}{{,}}\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}} $ $ ({\lambda ^2}_{1j})$ 和${{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}}{{,}}\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}({\lambda _{1j}}{\lambda _{1j - 1}})$ , 可由Rauch-Tung-Striebel (RTS)最优平滑算法来实现内容算法2. RTS最优平滑算法
$1)$ 后向迭代$$\begin{split} &{S_j} = {P_{j|j}}P_{j{\rm{ + 1|}}j}^{ - 1} \\ & {\widehat \lambda _{j|k}} = {\widehat \lambda _j} + {S_j}({\widehat \lambda _{j + 1|k}} - {\widehat \lambda _j}) \\ &{P_{j|k}} = {P_{j|j}} + {S_j}^2({P_{j + 1|k}} - {P_{j + 1|j}}) \end{split} $$ (22) 2)初始化
$${M_{k|k}} = (1 - K(k)\Delta {t_k}){P_{k - 1|k - 1}}$$ (23) 3)协方差后向迭代
$${M_{j|k}} = {P_{j|j}}{S_{j - 1}} + {S_j}({M_{j + 1|k}} - {P_{j|j}}){S_{j - 1}}$$ (24) 引理2 [24]. 基于当前时刻估计的未知参数
$\widehat {\boldsymbol{\theta }}_{1k}^{(i)}$ 和监测数据${{\boldsymbol{x}}_{0:k}}$ , 有$$\begin{split} &{{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}}{{,}}\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}({\lambda _{1j}}) = {\widehat \lambda _{j|k}} \\ & {{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}}{{,}}\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}({\lambda ^2}_{1j}) = {\widehat \lambda ^2}_{j|k} + {P_{j|k}} \\ &{{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}}{{,}}\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}({\lambda _{1j}}{\lambda _{1j - 1}}) = {\widehat \lambda _{j|k}}{\widehat \lambda _{j - 1|k}} + {M_{j|k}} \end{split} $$ (25) 结合式(21)和式(25), 对未知变量
${{\boldsymbol{\theta }}_{1}}$ 求偏导, 且偏导等于0, 可得到第$i + 1$ 步的参数估计值$\widehat {\boldsymbol{\theta }}_{1k}^{(i + 1)}$ 如下:$$ \begin{split} & {({a_1})_{0k}^{(i + 1)} = {{({a_1})}_{0|k}}}\\ & {({P_1})_{0k}^{(i + 1)} = {{({P_1})}_{0|k}}}\\ & (k_1^2)_k^{(i + 1)} = \frac{1}{k}\sum\limits_{j = 1}^k \frac{{\ln ({t_j} - {t_{j - 1}})}}{{{t_j} - {t_{j - 1}}}} \cdot\\ &\qquad\qquad\quad({C_{j|k}} - 2{C_{j,j - 1|k}} + {C_{j - 1|k}}) \\ &(\sigma _1^2)_k^{(i + 1)} = \frac{1}{k}\sum\limits_{j = 1}^k \frac{1}{{{t_j}-{t_{j - 1}}}}\cdot\\ &\qquad\Big[ ({x_j} - {x_{j - 1}})^2 - 2\hat\lambda_{1{j - 1|k}}{({x_j} - {x_{j - 1}})}\cdot \\ &\qquad{({t_j} - {t_{j - 1}}) + ({t_j} - {t_{j - 1}})^2C_{j - 1|k}} \Big] \end{split}$$ (26) 式中:
${C_{j|k}} = {{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}}{{,}}\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}({\lambda ^2}_{1j}),\; {C_{j,j - 1|k}} = $ $ {{\rm{E}}_{{\lambda _{1k}}|{{\boldsymbol{x}}_{0:k}}{{,}}\widehat {\boldsymbol{\theta }}_{1k}^{(i)}}}({\lambda _{1j}}{\lambda _{1j - 1}}) .$ 第二阶段参数估计方法同上, 不再赘述. 在退化实验中, 监测到的设备性能退化数据一般为离散值, 变点
$\tau $ 的值通常未知. SIC是Akaike信息准则的改进, 对变点检测效果良好, 下面通过SIC, 确定变点$\tau $ 的值[15].3.3 变点检测
SIC是由Schwarz. 于1978年提出, 可以判断模型是否存在变点问题. 其原理是如果待检测序列存在变点, 其样本的熵要大于不存在变点的样本的熵[25]. 利用SIC来估计变点的个数和位置是较为简单的, 对变点的检测效果良好. 其定义为
$${\rm{SIC}} = - 2{\rm{ln}}L (\widehat {\boldsymbol{\theta }}) + p\ln m$$ (27) 式中,
$L(\widehat {\boldsymbol{\theta }})$ 是模型的极大似然函数,$\widehat {\boldsymbol{\theta }}$ 是${\boldsymbol{\theta }}$ 的极大似然估计;$p$ 是模型中的自由参数个数;$m$ 是样本大小. 依据SIC原则, 为了确定变点所在位置, 对本文做出以下假设:原假设
${H_0}$ : 各参数值相等, 表示模型中不存在变点.备择假设
${H_1}$ : 存在一个变点$\tau $ , 在$\tau $ 之前一个阶段按${X_1}(t;{\lambda _1},\sigma _1^2)$ 退化, 在$\tau $ 之后一个阶段按${X_2}(t;{\lambda _2},\sigma _2^2)$ 退化.根据式(27), 基于原假设
${H_0}$ 下的${\rm{SIC}}(m)$ 值为$$\begin{split} &{\rm{SIC}}(m) = m\ln 2\pi + m {\rm{ln}} \sum\limits_1^m {(\Delta {x_i}} - \Delta \overline x {)^{\rm{2}}} +\\ &\qquad\;\;\qquad m + (2 - m)\ln m \\ &\Delta \overline x = \frac{1}{m}\sum\limits_1^m {\Delta {x_i}} \end{split} $$ (28) 基于备择假设
${H_1}$ 下的${\rm{SIC}}(k)$ 值为$$\begin{split} &{\rm{SIC}}(k) = m\ln 2\pi + k\ln\frac{1}{k}\sum\limits_1^k {(\Delta {x_i}} - \Delta \overline {{x_1}} {)^2}+ \\ & \quad\qquad 4\ln m + (m - k)\ln\frac{1}{m}\sum\limits_{k + 1}^m {(\Delta {x_i}} - \Delta \overline {{x_2}} {)^2} - m \\ &\Delta \overline {{x_1}} = \frac{1}{k}\sum\limits_1^k {\Delta {x_i}} ,{\rm{ }}\Delta \overline {{x_2}} = \frac{1}{{m - k}}\sum\limits_{k + 1}^m {\Delta {x_i}} \end{split} $$ (29) 根据SIC, 如果
${\rm{SIC}}(m) > \mathop {\min }\limits_{2 < k \leq m - 2} {\rm{SIC}}(k)$ , 则拒绝原假设${H_0}$ , 认为存在变点. 同时估计的变点值$\widehat \tau = {t_{\widehat k}}$ 为$${\rm{SIC}}(\widehat k) = \mathop {\min }\limits_{2 < k \leq m - 2} {\rm{SIC}}(k)$$ (30) 为了描述不同设备变点时间的个体差异性, 本文假设变点时间
$\tau $ 服从随机分布. 相比于其他分布, Gamma分布能包含其他常见分布, 如指数分布等, 并且形状参数$\alpha $ 越大, Gamma分布越逼近正态分布, 计算也较为容易. 因此, 本文假设变点时间$\tau $ 服从Gamma分布, 且形状参数为$\alpha $ 、尺度参数为$\beta $ .通过SIC方法, 利用退化先验信息可离线估计设备的变点发生时间, 通过Gamma分布的统计形式, 进而得到变点时间
$\tau $ 的分布参数及概率分布函数$p(\tau )$ .4. 实例验证
在本节, 通过引入蒙特卡洛仿真, 验证本文基于自适应Wiener过程所提出的模型优越性. 然后, 将所提出的模型应用于锂电池容量退化数据中.
4.1 模拟仿真
为了验证本文所提模型能够解决现有两阶段模型不能刻画测量间隔不均匀、测量频率不一致的问题, 在此考虑将所提出模型与Zhang等[19]所提的两阶段模型作比较.
首先, 在这里主要考虑漂移参数的随机效应, 利用两种退化模型分别对仿真得到的数据进行退化建模, 并利用第四节模型参数估计方法, 可求得两个阶段模型参数估计值, 最后得到两种方法各自预测的RUL结果. 在此, 本文增加一个仿真验证的例子, 设定仿真参数值为
$({\lambda _1},{\lambda _2},{\sigma _1},{\sigma _2},{k_1},{k_2}) = (2,1, $ $ 1,1,1,0.5)$ , 并且生成一组间隔为1、次数为200的监测数据. 为了证明本文所提模型比现有模型适用于测量间隔不均匀、测量频率不一致情况, 考虑将数据变为以下情况: 存在测量间隔为1、2、4的测量间隔不均匀的混合数据. 其中, 变点发生时间设为120, 失效阈值为155.图1为本文所提模型和Zhang等所提的两阶段模型在不均匀退化数据下的RUL结果. 从图中, 可以看到RUL预测的PDF随着时间逐渐变窄, 表明预测的不确定性越来越小. 对于测量间隔不均匀的数据, 所提模型的PDF曲线更加窄而尖锐, 说明在寿命预测方面可以取得更好的预测结果.
为了进一步量化测量间隔不均匀条件下两种模型的预测结果, 本文采用可靠性领域常用的性能指标: 绝对误差(Absolute error, AE)和相对误差(Relative error, RE)指标.
从图2和表1中, 可以看出随着退化数据的逐渐累积, 两种模型各自预测的误差值也在随之下降. 总体上看, 本文所提模型的预测准确度要优于Zhang等所提的模型. 当退化数据呈不均匀分布时, 相比较Zhang等所提的两阶段模型, 本文所提模型能较好解决这种情况带来的影响, 更为准确地估计参数, 进而提高预测准确度.
表 1 不同监测点相对误差的比较结果Table 1 Comparison results of relative error at different monitoring points方法 不同监测点相对误差 20 40 60 80 提出模型 20.8 % 19.6 % 15.2 % 12.3 % Zhang等[19]模型 27.2 % 25.3 % 20.1 % 17.6 % 综上, 本文基于自适应Wiener过程所建立的模型相比于Zhang等[19]的两阶段退化模型, 可以克服测量间隔不均匀、测量频率不一致问题, 较为准确预测设备RUL.
4.2 锂电池实例验证
本节中, 基于马里兰大学Pecht[26]课题组的锂电池容量退化数据进行实例验证. 该数据是在室温条件下通过充放电实验得到的, 记录了电池状态信息(包括容量)随充放电次数的变化. 由于复杂的化学反应, 锂电池容量开始时经历一个平稳退化期, 随着充放电的进行, 由于固体电解质层在电极上的生长以及副反应导致的活性材料的损失, 导致锂电池容量在后一阶段迅速降低[9]. 编号为CS2-35、CS2-36、CS2-37、CS2-38四组电池容量退化数据如图3所示, 图中退化过程呈现出明显的两阶段特性, 这里采用CS2-37锂电池数据进行RUL预测验证, 其余三组(即CS2-35、CS2-36和CS2-38)用作变点时间离线参数估计.
首先, 对CS2-37锂电池利用SIC进行变点检测, 确定变点所在时刻. 根据SIC, 分别计算相应的SIC值, 对应的SIC值变化趋势如图4所示.
从图中可以观测到, 第96监测点数值最小, 由于本文采用两阶段随机退化建模, 从起始退化到第96监测点, 监测时间较短不适宜作为变点, 因而, 对于此情况忽略不计. 此时,
${\rm{SIC}}(1\,006) = - 16\,735 > $ $ {\rm{SIC}}(736)$ , 且${\rm{SIC}}(736) = - 26\,891$ 最小, 所以CS2-37锂电池退化中存在变点, 且变点发生在第736监测点. 同理, 其余三组电池数据辨识到的变点分别为623、681和753. 由于变点时间$\tau $ 服从Gamma分布, 可利用辨识得到的数据进行拟合, 可求得形状参数和尺度参数分别为$\alpha = 106$ 、$\beta = 7.98$ . 将变点引入到参数估计中, 再结合退化数据, 可得两个阶段漂移系数的均值和方差估计值分别为${u_a} = $ $ {\rm{8}}.{\rm{0}}56 \times {10^{ - 4}}$ 、${\sigma _a} = 2.94 \times {10^{ - {\rm{5}}}}$ 、${u_b} = 0.0022{\rm{1}}$ 、${\sigma _b} = $ $ {\rm{7}}{\rm{.41}} \times {10^{ - {\rm{5}}}}$ . 利用上述估计值, 结合CS2-37锂电池数据和卡尔曼滤波技术实现漂移系数在线更新. 图5、图6展示了隐含状态即漂移系数的在线更新过程. 结果表明, 两个阶段容量的退化率相差较大, 如果只用单一阶段的退化过程进行退化建模, 误差将会比较大. 由于变点发生时刻在第736监测点, 因此第一阶段模型参数在变点出现后不再更新, 第二阶段模型参数在变点之前不进行更新.为了验证所提模型的有效性, 用建立模型对CS2-37锂电池退化数据进行拟合. 图7为符合Wiener过程的锂电池容量退化预测情况, 从图中可以看出所建立的模型能较好地反映锂电池退化过程.
由于变点在第736监测点, 而一般电池容量的失效阈值定义为电容量损失到初始容量的百分比. 在本文中设定失效阈值为初始电容量的45 %. 为了说明本文方法的有效性和合理性, 将估计的模型参数代入到式(8)中, 可得到RUL的PDF和预测值, 如图8所示.
从图8中可以看出, 对于锂电池退化数据的RUL预测, 所提模型与单一阶段线性模型相比较, 前者更为准确. 与Zhang等提出的两阶段模型相比, 本文的方法能取得更好的结果, 且随着监测数据的增加, RUL预测结果不确定性越来越小, 精度越来越高. 为了更加直观说明本文方法有效性, 给出几种方法的RUL预测绝对误差和
$\alpha {\rm{ - }}\beta $ 性能指标对预测结果进行验证.从图9、图10中可以看出, 与单一阶段线性模型相比, 本文方法考虑变点前后呈现两阶段特征, 即变点前后的退化速率存在明显差异进行建模, 且考虑同批产品个体差异性的影响, 其模型更符合实际退化情况. 与Zhang等所提的两阶段模型相比, 本文方法考虑Wiener过程模型存在测量间隔不均匀、测量频率不一致以及在RUL预测中忽略了自适应漂移的可变性等三点不足, 结果表明在监测前期退化数据较少时, 所提模型能取得较好的预测结果. 其原因是本文所提模型与Zhang等所提模型相比, 噪声项增加了一个自适应漂移项
$k\int_0^{{l_k}} {W(\tau ){\rm{d}}\tau }$ , 它是一个随时间变化的过程. 在设备监测初期, 剩余寿命${l_k}$ 值比较大, 无法忽略自适应漂移可变性的影响, 因此所提模型预测结果优于现有模型. 随着电池充放电循环在寿命将尽时, 自适应漂移项的影响降到最低, 此时两种退化模型结构相似, 进而提供近似的预测结果.此外, 通过引入相对误差指标进一步量化预测准确度, 在寿命的35 %、55 %、75 %和95 %分位点给出两种方法相对误差的比较结果.
从表2可以看出, 本文方法可以有效减小RUL预测的RE, 进而提高预测精度, 尤其在95 %分位点处, RUL预测的RE仅为0.66 %. 综上, 本文所提出的两阶段自适应Wiener过程模型和方法, 预测精度高, 并且具有一定的适用性, 可为后续设备的备件订购、最优替换等维护决策提供依据.
表 2 不同寿命分位点相对误差的比较结果Table 2 Comparison results of relative error at different life quantiles方法 寿命分位点相对误差 35 % 55 % 75 % 95 % 提出模型 6.04 7.65 7.91 0.66 Zhang 等[19]模型 11.4 9.37 9.68 1.46 5. 结论
本文针对工程实际中存在阶段性退化特征的一类设备, 建立了两阶段自适应Wiener过程退化模型进行RUL预测, 重点阐述了模型参数估计和隐含状态更新方法, 最后通过锂电池退化数据验证所提模型的可靠性和有效性. 具体结论如下:
1)建立自适应Wiener过程模型, 克服了一般Wiener过程模型存在的测量间隔不均匀、测量频率不一致以及在RUL预测中没有利用实时监测数据自适应更新漂移系数三点不足, 提出一种新的考虑个体差异性两阶段预测模型.
2)在首达时间意义下, 推导出两阶段自适应Wiener过程模型寿命和剩余寿命PDF的解析表达式, 实现RUL的实时估计.
3)基于Kalman滤波算法和EM算法进行参数估计和自适应更新, 实现设备的实时可靠性评估, 为维护决策提供依据. 最后通过锂电池实例验证了本文所提方法的有效性.
本文主要适用于退化数据呈两阶段随机退化设备研究, 但对于大型复杂设备, 由于环境以及工作任务变换的影响, 可能存在多个工况或状态切换的现象, 进而导致多阶段情况发生. 下一步可深入拓展多状态多阶段复杂随机系统的退化建模、RUL预测与维护决策的问题研究.
附录A
A.1 两阶段自适应Wiener过程寿命PDF
A.1.1 变点处退化量已知
根据文献[21]
$, $ 线性自适应Wiener过程寿命分布的PDF为$$\begin{split} {f_{T/\lambda }}(t/\lambda ) =\;& \frac{1}{{\sqrt {2\pi \varphi (t)} }} \cdot \exp \left\{ { - \frac{{{{(D - {X_0} - {\lambda _n}t)}^2}}}{{2\varphi (t)}}} \right\} \cdot\\ &\left[ {\frac{{D - {X_0} - {\lambda _n}t}}{{\varphi (t)}}\varphi '(t) + {\lambda _n}} \right]\\[-10pt] \end{split}\tag{A1} $$ 引理A1. 若
$p\sim {\rm{N}}(u,{\sigma ^2})$ , A、B都为常数, C为正实数, 则下式成立$$\begin{split} &{{\rm{E}}_p}\left[ {({A} - p)\exp \left( { - \frac{{{{(B - p)}^2}}}{{2C}}} \right)} \right] = \\ &\sqrt {\frac{C}{{{\sigma ^2} + C}}} \left( {A - \frac{{{\sigma ^2}B + uC}}{{{\sigma ^2} + C}}} \right) \cdot \exp \left( { - \frac{{{{(B - u)}^2}}}{{2({\sigma ^2} + C)}}} \right)\\[-10pt] \end{split}\tag{A2}$$ 当考虑设备个体差异性时, 漂移参数
$\lambda $ 服从正态分布, 即$\lambda \sim {\rm{N}}({u_\lambda },\sigma _\lambda ^2)$ . 此时, 基于引理A1和全概率公式可求得设备寿命分布的PDF为$$\begin{split} &{f_T}(t) = {{\rm{E}}_{\lambda |{{\boldsymbol{X}}_{1:n}}}}({f_{T/\lambda }}(t/\lambda )= \\ & \frac{1}{{\sqrt {2\pi (\varphi (t) + \sigma _\lambda ^2{t^2})} }}\left[ \begin{aligned} &\frac{{\varphi '(t)}}{{\varphi (t)}}(D - {X_0}) + (1 - \frac{{\varphi '(t)}}{{\varphi (t)}}t) \\ &\frac{{{u_\lambda }\varphi (t) + (D - {X_0})\sigma _\lambda ^2t}}{{\varphi (t) + \sigma _\lambda ^2{t^2}}} \end{aligned} \right]\cdot \\ & \exp \left\{ { - \frac{{{{(D - {X_0} - {u_\lambda }t)}^2}}}{{2(\varphi (t) + \sigma _\lambda ^2{t^2})}}} \right\}\\[-10pt] \end{split}\tag{A3} $$ 将式(A3)扩展到两阶段模型. 在此, 本文考虑最简单的一种情况: 假定变点发生时间
$\tau $ 和变点处的退化量${X_\tau }$ 已知.当寿命满足
$t < \tau $ 时, 此时寿命分布的PDF与单一阶段自适应Wiener过程的寿命分布一致,$$\begin{split} {f_T}(t) =\;& \frac{1}{{\sqrt {2\pi ({\varphi _1}(t) + \sigma _a^2{t^2})} }}.\\ &\left[ \begin{aligned} &\frac{{{\varphi _1}'(t)}}{{{\varphi _1}(t)}}(D - {X_0}) + (1 - \frac{{{\varphi _1}'(t)}}{{{\varphi _1}(t)}}t) \\ &\frac{{{u_a}{\varphi _1}(t) + (D - {X_0})\sigma _a^2t}}{{{\varphi _1}(t) + \sigma _a^2{t^2}}} \end{aligned} \right] \cdot \\ &\exp \left\{ { - \frac{{{{(D - {X_0} - {u_a}t)}^2}}}{{2({\varphi _1}(t) + \sigma _a^2{t^2})}}} \right\} \end{split}\tag{A4} $$ 当寿命满足
$t > \tau $ 时, 不用考虑变点退化量${X_\tau }$ 的影响, 此时寿命分布的PDF为,$$\begin{split} & {f_T}(t) = \frac{1}{{\sqrt {2\pi ({\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2})} }}\cdot\\ &\left[ \begin{aligned} &\frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(D - {X_\tau }) + (1 - \frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(t - \tau )) \\ &\frac{{{u_b}{\varphi _2}(t - \tau ) + (D - {X_\tau })\sigma _b^2(t - \tau )}}{{{\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2}}} \end{aligned} \right] \cdot\\ &\exp \left\{ { - \frac{{{{(D - {X_\tau } - {u_b}(t - \tau ))}^2}}}{{2({\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2})}}} \right\}\\[-15pt] \end{split}\tag{A5} $$ A.1.2 变点处退化量未知
实际上, 在变点时间发生前, 变点处退化量的值
${X_\tau }$ 通常未知. 因此, 为了得到寿命估计值, 首先要得到首达意义下${X_\tau }$ 的分布形式, 即在${X_\tau } < D$ 条件下经过时间$\tau $ , 退化量从0到${X_\tau }$ 的转移概率${g_\tau }({X_\tau })$ .将式(A3)扩展到两阶段模型. 当寿命满足
$t < \tau $ 时, 仅受随机参数${\lambda _1}$ 影响, 此时寿命分布的PDF为$$\begin{split} &{f_T}(t) = \int_{ - \infty }^{ + \infty } {\left( \begin{aligned} & \frac{1}{{\sqrt {2\pi {\varphi _1}(t)} }} \cdot \frac{1}{{\sqrt {2\pi \sigma _a^2} }}\\ &\exp \left( { - \frac{{{{({\lambda _1} - {u_a})}^2}}}{{\sigma _a^2}}} \right) \\ &\left[ {\frac{{D - {X_0} - {\lambda _1}t}}{{{\varphi _1}(t)}}{\varphi' _1}(t) + {\lambda _1}} \right] \cdot\\ &\exp \left\{ { - \frac{{{{(D - {X_0} - {\lambda _1}t)}^2}}}{{2{\varphi _1}(t)}}} \right\} \end{aligned} \right)} {\rm{d}}{\lambda _1}= \\ &\frac{1}{{\sqrt {2\pi ({\varphi _1}(t) + \sigma _a^2{t^2})} }}\left[ \begin{aligned} &\frac{{{\varphi' _1}(t)}}{{{\varphi _1}(t)}}(D - {X_0}) + (1 - \frac{{{\varphi _1}'(t)}}{{{\varphi _1}(t)}}t) \\ &\frac{{{u_a}{\varphi _1}(t) + (D - {X_0})\sigma _a^2t}}{{{\varphi _1}(t) + \sigma _a^2{t^2}}} \end{aligned} \right] \cdot \\ &\exp \left\{ { - \frac{{{{(D - {X_0} - {u_a}t)}^2}}}{{2({\varphi _1}(t) + \sigma _a^2{t^2})}}} \right\} \\[-15pt]\end{split}\tag{A6} $$ 引理A2. 对于线性两阶段自适应Wiener过程, 如果漂移系数
$\lambda $ 服从正态分布, 即$\lambda \sim {\rm{N}}({u_a},\sigma _a^2)$ , 那么首达意义下的转移概率为$$\begin{split} &{g_\tau }({X_\tau }|{{{u}}_a},{\sigma _a})=\int_{ - \infty }^{ + \infty } {{g_\tau }({X_\tau }|{\lambda _1})p(} {\lambda _1}){\rm{d}}{\lambda _1}= \\ & \int_{ - \infty }^{ + \infty } {{g_\tau }({X_\tau }|{\lambda _1})\frac{1}{{\sqrt {2\pi \sigma _a^2} }}\exp \left( { - \frac{{{{({\lambda _1} - {u_a})}^2}}}{{\sigma _a^2}}} \right)} {\rm{d}}{\lambda _1} =\\ & \left[ {1{\rm{ - }}\exp \left( {\frac{{4{w^2} - 4{X_\tau }w}}{{2{\varphi _1}(\tau )}}} \right)} \right] \cdot \frac{1}{{\sqrt {2\pi ({\varphi _1}(\tau ) + \sigma _a^2{\tau ^2})} }} \cdot\\ &\exp \left[ { - \frac{{{{({X_\tau } - {u_a}\tau )}^2}}}{{2({\varphi _1}(\tau ) + \sigma _a^2{\tau ^2})}}} \right]\\[-15pt] \end{split} \tag{A7}$$ 当寿命满足
$t > \tau $ 时, 在这种情况下受随机参数${X_\tau }$ 和${\lambda _2}$ 的影响, 此时寿命分布的PDF为$$\begin{split} &{f_T}(t) =\\ & \int_{ - \infty }^D {\left( \begin{aligned} &\frac{1}{{\sqrt {2\pi ({\varphi _2}(t - \tau ) + \sigma _2^2{{(t - \tau )}^2})} }} \\ & \left[ \begin{aligned} &\frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(D - {X_\tau }) + \\ &\left(1 - \frac{{{\varphi _2}'(t - \tau )}}{{{\varphi _2}(t - \tau )}}(t - \tau )\right) \\ &\frac{{{u_b}{\varphi _2}(t - \tau ) + (D - {X_\tau })\sigma _b^2(t - \tau )}}{{{\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2}}} \end{aligned} \right] \cdot\\ &\exp \left\{ { - \frac{{{{(D - {X_\tau } - {u_b}(t - \tau ))}^2}}}{{2({\varphi _2}(t - \tau ) + \sigma _b^2{{(t - \tau )}^2})}}} \right\} \cdot\\ &\qquad{g_\tau }({X_\tau }|{{u}_a},{\sigma _a}) \end{aligned} \right){\rm{d}}{X_\tau }{\kern 1pt} } \approx\\ &{A_1} - {B_1}\\[-15pt] \end{split} \tag{A8}$$ 其中:
$$\begin{split}& {A_1} = \frac{{\varphi _2^{'}(t - \tau )}}{{{\varphi _2}(t - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a1}^2 + \sigma _{b1}^2)}}} \exp \left[ { - \frac{{{{({u_{a1}} - {u_{b1}})}^2}}}{{2(\sigma _{a1}^2 + \sigma _{b1}^2)}}} \right] \cdot \\ &\qquad\left\{ \frac{{{u_{b1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}\Phi \left( {\frac{{{u_{b1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right)+\right. \\ &\qquad\left.\sqrt {\frac{{\sigma _{a1}^2\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}} \phi \left( {\frac{{{u_{b1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right) \right\} \\ &{B_1} = \exp \left\{ {\frac{{2{u_a}D}}{{\varphi _1^{'}(\tau )}} + \frac{{2({D^2}\sigma _a^4\tau + {D^2}\sigma _a^2\varphi _1^{'}(\tau ))}}{{(\varphi _1^{'}(\tau ) + \sigma _a^2\tau )\varphi _1^{'}{{(\tau )}^2}}}} \right\} \\ &\qquad\frac{{\varphi _2^{'}(t - \tau )}}{{{\varphi _2}(t - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a1}^2 + \sigma _{b1}^2)}}} \exp \left[ { - \frac{{{{({u_{a1}} - {u_{c1}})}^2}}}{{2(\sigma _{a1}^2 + \sigma _{b1}^2)}}} \right] \cdot \\ &\qquad\left\{\frac{{{u_{c1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}\Phi \left( {\frac{{{u_{c1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right) +\right.\\ & \qquad\left.\sqrt {\frac{{\sigma _{a1}^2\sigma _{b1}^2}}{{\sigma _{a1}^2 + \sigma _{b1}^2}}} \phi \left( {\frac{{{u_{c1}}\sigma _{a1}^2 + {u_{a1}}\sigma _{b1}^2}}{{\sqrt {\sigma _{a1}^2\sigma _{b1}^2(\sigma _{a1}^2 + \sigma _{b1}^2)} }}} \right) \right\} \\ &{u_{a1}} = {u_b}(t - \tau ),\;\;{u_{b1}} = D - {u_a}\tau \\ &{u_{c1}} = - D - {u_a}\tau - \frac{{\sigma _a^2\tau }}{{\varphi _1^{'}(\tau )}} \\& \sigma _{a1}^2 = {\varphi _2}(t - \tau ) + \sigma _b^2{(t - \tau )^2}\;\;\;\\&\sigma _{b1}^2 = {\varphi _1}(\tau ) + \sigma _a^2{\tau ^2} \end{split}\tag{A9} $$ A.2 两阶段自适应Wiener过程剩余寿命PDF
令当前时刻
${t_k}$ 的退化状态${{{x}}_k}$ , 用${l_k}$ 表示设备剩余寿命,${f_L}({l_k})$ 表示设备RUL分布的PDF, 在随机退化速率${\lambda _1}$ 和${\lambda _2}$ 的影响下, 可获得首达意义下两阶段自适应Wiener过程模型RUL的PDF.情况1. 当前时刻
${t_k}$ 位于变点前, 即${t_k} < \tau $ , 此时随机设备退化失效又存在两种情况:1) 失效阈值位于变点前, 即
${l_k} + {t_k} \leq \tau$ , 在此种情况下仅受随机参数${\lambda _1}$ 影响, RUL分布的PDF为$$\begin{split} &{f_L}({l_k}) = \\ &\int_{ - \infty }^{ + \infty } \left\{ \frac{1}{{\sqrt {2\pi {\varphi _1}({l_k})} }} \cdot \exp \left[ { - \frac{{{{(D - {X_k} - {\lambda _1})}^2}}}{{2{\varphi _1}({l_k})}}} \right]\cdot\right. \\ &\left( {\frac{{D - {X_k} - {\lambda _1}{l_k}}}{{{\varphi _1}({l_k})}}{\varphi _1}'({l_k})} \right)\cdot \\ &\left.\frac{1}{{\sqrt {2\pi \sigma _a^2} }}\exp \left( { - \frac{{{{({\lambda _1} - {u_a})}^2}}}{{\sigma _a^2}}} \right) \right\} {\rm{d}}{\lambda _1}= \\ &\frac{1}{{\sqrt {2\pi ({\varphi _{\rm{1}}}({l_k}) + \sigma _a^2{l_k^2})} }} \cdot \exp \left\{ { - \frac{{{{(D - {X_k} - {u_a}{l_k})}^2}}}{{2({\varphi _{\rm{1}}}({l_k}) + \sigma _a^2{l_k^2)}}}} \right\} \cdot\\ &\left[ \begin{aligned} &\frac{{{\varphi '_{\rm{1}}}({l_k})}}{{{\varphi _{\rm{1}}}({l_k})}}(D - {X_k}) + \left(1 - \frac{{{\varphi' _{\rm{1}}}({l_k})}}{{{\varphi _{\rm{1}}}({l_k})}}{l_k}\right) \\ &\frac{{{u_a}{\varphi _{\rm{1}}}({l_k}) + (D - {X_k})\sigma _a^2{l_k}}}{{{\varphi _{\rm{1}}}({l_k}) + \sigma _a^2{l_k^2}}} \end{aligned} \right]\\[-20pt] \end{split}\tag{A10} $$ 2) 失效阈值位于变点后, 即
${l_k} + {t_k} > \tau $ , 此时可以看作受两个随机参数的影响, 即${\lambda _2}$ 和${X_\tau }$ , RUL的PDF可以用两重积分表示$$\begin{split} {f_L}({l_k}) =\;& \int_{ - \infty }^D \int_{ - \infty }^{ + \infty } \Bigg\{ \frac{1}{{\sqrt {2\pi {\varphi _2}({l_k} + {t_k} - \tau )} }}{g_\tau }({X_\tau }|{{\rm{u}}_a},{\sigma _a})p({\lambda _2}) \cdot \\ &\exp \left[ { - \frac{{{{(D - {X_k} - {\lambda _2}({l_k} + {t_k} - \tau ))}^2}}}{{2{\varphi _2}(({l_k} + {t_k} - \tau ))}}} \right]\cdot \\ &\Bigg( \frac{{D - {X_k} - {\lambda _2}({l_k} + {t_k} - \tau )}}{{{\varphi _2}(({l_k} + {t_k} - \tau ))}} \\ &{\varphi' _2}(({l_k} + {t_k} - \tau )) \Bigg) \Bigg\} {\rm{ }}{\rm{d}}{\lambda _2}{\rm{d}}{X_\tau }=\\ & \int_{ - \infty }^D \int_{ - \infty }^{ + \infty }\Bigg\{ \frac{1}{{\sqrt {2\pi {\varphi _2}({l_k} + {t_k} - \tau )} }} \cdot\\ &\Bigg( \frac{{D - {X_k} - {\lambda _2}({l_k} + {t_k} - \tau )}}{{{\varphi _2}(({l_k} + {t_k} - \tau ))}} {\varphi' _2}(({l_k} + {t_k} - \tau )) \Bigg) \cdot \\ &\exp \left[ { - \frac{{{{(D - {X_k} - {\lambda _2}({l_k} + {t_k} - \tau ))}^2}}}{{2{\varphi _2}(({l_k} + {t_k} - \tau ))}}} \right]. \\ &\left[ {1{\rm{ - }}\exp \left( {\frac{{4{w^2} - 4{X_\tau }w}}{{2{\varphi _1}(\tau )}}} \right)} \right] \cdot\frac{1}{{\sqrt {2\pi ({\varphi _1}(\tau ) + \sigma _a^2{\tau ^2})} }} \cdot\\ &\exp \left[ { - \frac{{{{({X_\tau } - {u_a}\tau )}^2}}}{{2({\varphi _1}(\tau ) + \sigma _a^2{\tau ^2})}}} \right] \cdot\\ &\frac{1}{{\sqrt {2\pi \sigma _b^2} }}\exp \left( { - \frac{{{{({\lambda _2} - {u_b})}^2}}}{{\sigma _b^2}}} \right) \Bigg\} {\rm{d}}{\lambda _2}{\rm{d}}{X_\tau } \approx\\ & {A_2} - {B_2}\\[-10pt] \end{split} \tag{A11}$$ 其中:
$$\begin{split} & {A_2} = \frac{{\varphi '_2({l_k} + {t_k} - \tau )}}{{{\varphi _2}({l_k} + {t_k} - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a2}^2 + \sigma _{b2}^2)}}} \cdot \\ &\qquad\exp \left[ { - \frac{{{{({u_{a2}} - {u_{b2}})}^2}}}{{2(\sigma _{a2}^2 + \sigma _{b2}^2)}}} \right] \cdot \\ &\qquad\left\{ \frac{{{u_{b2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}\Phi \left( {\frac{{{u_{b2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right)+ \right.\\ &\qquad\left. \sqrt {\frac{{\sigma _{a2}^2\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}} \phi \left( {\frac{{{u_{b2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right) \right\}\\ &{B_2} = \frac{{\varphi '_2({l_k} + {t_k} - \tau )}}{{{\varphi _2}({l_k} + {t_k} - \tau )}}\sqrt {\frac{1}{{2\pi (\sigma _{a2}^2 + \sigma _{b2}^2)}}}. \\ &\qquad\exp \left[ { - \frac{{{{({u_{a2}} - {u_{c2}})}^2}}}{{2(\sigma _{a2}^2 + \sigma _{b2}^2)}}} \right] .\\ &\qquad\exp \left\{ \frac{{2({{(D - {X_k})}^2}\sigma _a^4\tau + {{(D - {X_k})}^2}\sigma _a^2\varphi _1^{'}(\tau - {t_k}))}}{{(\varphi' _1(\tau - {t_k}) + \sigma _a^2(\tau - {t_k}))\varphi _1^{'}{{(\tau - {t_k})}^2}}}+ \right. \end{split} $$ $$\begin{split} &\qquad\left. \frac{{2{u_a}(D - {X_k})}}{{\varphi' _1(\tau - {t_k})}} \right\} \cdot\\ &\qquad\left\{ \frac{{{u_{c2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}\Phi \left( {\frac{{{u_{c2}}\sigma _{a2}^2 + {u_{a2}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right)+\right. \\ &\left.\qquad \sqrt {\frac{{\sigma _{a2}^2\sigma _{b2}^2}}{{\sigma _{a2}^2 + \sigma _{b2}^2}}} \phi \left( {\frac{{{u_{c1}}\sigma _{a2}^2 + {u_{a1}}\sigma _{b2}^2}}{{\sqrt {\sigma _{a2}^2\sigma _{b2}^2(\sigma _{a2}^2 + \sigma _{b2}^2)} }}} \right)\right\} \\ &{u_{a2}} = {u_b}({l_k} + {t_k} - \tau ),\;\;{u_{b2}} = D - {X_k} - {u_a}(\tau - {t_k})\;\;\;\;\;\; \\ &{u_{c2}} = - D + {X_k} - {u_a}(\tau - {t_k}) - \frac{{\sigma _a^2(\tau - {t_k})}}{{\varphi' _1(\tau - {t_k})}} \\ &\sigma _{a2}^2 = {\varphi _2}({l_k} + {t_k} - \tau ) + \sigma _b^2{({l_k} + {t_k} - \tau )^2}{\rm{ }} \\ & \sigma _{b2}^2 = {\varphi _1}(\tau - {t_k}) + \sigma _a^2{(\tau - {t_k})^2}\\[-10pt] \end{split}\tag{A12} $$ 情况2. 当前时刻
${t_k}$ 位于变点后, 即${t_k} > \tau $ , 此时仅受随机参数${\lambda _2}$ 影响, RUL的PDF为$$ \begin{split} {f_L}({l_k}) = \;&\frac{1}{{\sqrt {2\pi ({\varphi _2}({l_k}) + \sigma _b^2{l_k^2})} }} \cdot\\ & \left[ \begin{aligned} &\frac{{{\varphi _2}'({l_k})}}{{{\varphi _2}({l_k})}}(D - {X_k}) + (1 - \frac{{{\varphi _2}'({l_k})}}{{{\varphi _2}({l_k})}}{l_k}) \\ &\frac{{{u_b}{\varphi _2}({l_k}) + (D - {X_k})\sigma _b^2{l_k}}}{{{\varphi _2}({l_k}) + \sigma _b^2{l_k^2}}} \end{aligned} \right] \cdot\\ &\exp \left\{ { - \frac{{{{(D - {X_k} - {u_b}{l_k})}^2}}}{{2({\varphi _2}({l_k}) + \sigma _b^2{l_k^2})}}} \right\} \end{split}\tag{A13} $$ -
[1] 陆宁云, 陈闯, 姜斌, 邢尹. 复杂系统维护策略最新研究进展: 从视情维护到预测性维护. 自动化学报, 2021, 47(1): 1-17.Lu Ning-Yun, Chen Chuang, Jiang Bin, Xing Yin. Latest progress on maintenance strategy of complex system: from condition-based main-tenance to predictive maintenance. Acta Automatica Sinica, 2021, 47(1): 1-17. [2] 喻勇, 司小胜, 胡昌华, 崔忠马, 李洪鹏. 数据驱动的可靠性评估与寿命预测研究进展: 基于协变量的方法. 自动化学报, 2018, 44(2): 216-227.Yu Yong, Si Xiao-Sheng, Hu Chang-Hua, Cui Zhong-Ma, Li Hong-Peng. Data driven reliability assessment and life-time prognostics: a review on covariate models. Acta Automatica Sinica, 2018, 44(2): 216-227. [3] Omshi E M, Grall A, Shemehsavar S. A dynamic auto-adaptive predictive maintenance policy for degradation with unknown parameters. European Journal of Operational Research, 2020, 282(1): 81-92. doi: 10.1016/j.ejor.2019.08.050 [4] Hu J W, Sun Q Z, Ye Z S, Zhou Q. Joint modeling of degradation and lifetime data for RUL prediction of deteriorating products. IEEE Transactions on Industrial Informatics. 2020, 17(7): 4521-4531. [5] 周东华, 魏慕恒, 司小胜. 工业过程异常检测、寿命预测与维修决策的研究进展. 自动化学报, 2013, 39(6): 711−722.Zhou Dong-Hua, Wei Mu-Heng, Si Xiao-Sheng. A survey on anomaly detection, life prediction and maintenance decision for industrial processes. Acta Automatica Sinica, 2013, 39(6): 711−722. [6] Wu S, Castro I T. Maintenance policy for a system with a weighted linear combination of degradation processes. European Journal of Operational Research, 2020, 280(1): 124-133. doi: 10.1016/j.ejor.2019.06.048 [7] 王泽洲, 陈云翔, 蔡忠义, 项华春, 王莉莉. 基于比例关系加速退化建模的设备剩余寿命在线预测. 系统工程与电子技术, 2021, 43(2): 584-592. doi: 10.12305/j.issn.1001-506X.2021.02.34Wang Ze-Zhou, Chen Yun-Xiang, Cai Zhong-Yi, Xiang Hua-Chun, Wang Li-Li. Equipment remaining useful lifetime online prediction based on accelerated degradation modeling with the proportion relationship. Systems Engineering and Electronics, 2021, 43(2): 584-592. doi: 10.12305/j.issn.1001-506X.2021.02.34 [8] 袁烨, 张永, 丁汉. 工业人工智能的关键技术及其在预测性维护中的应用现状. 自动化学报, 2020, 46(10): 2013-2030.Yuan Ye, Zhang Yong, Ding Han. Research on key technology of industrial artificial intelligence and its application in predictive maintenance. Acta Automatica Sinica, 2020, 46(10): 2013-2030. [9] Zhou R S, Serban N, Gebraeel N. Degradation based residual life prediction under different environments. The Annals of Applied Statistics, 2014, 8(3): 1671-1689. [10] Burgess W L. Valve regulated lead acid battery float service life estimation using a Kalman filter. Journal of Power Sources, 2009, 191(1): 16-21. doi: 10.1016/j.jpowsour.2008.12.123 [11] Wang X L, Jiang P, Guo B, Cheng Z J. Real-time reliability evaluation for an individual product based on change-point Gamma and Wiener process. Quality and Reliability Engineering International, 2014, 30(4): 513-525. doi: 10.1002/qre.1504 [12] Ng T S. An application of the EM algorithm to degradation modeling. IEEE Transactions on Reliability, 2008, 57(1): 2-13. doi: 10.1109/TR.2008.916867 [13] Yuan T, Bae S J, Zhu X. A Bayesian approach to degradation-based burn-in optimization for display products exhibiting two-phase degradation patterns. Reliability Engineering & System Safety, 2016, 155(11): 55-63. [14] Bae S J, Kvam P H. A change-point analysis for modeling incomplete burn-in for light displays. IIE Transactions, 2006, 38(3): 489-498. [15] Yan W A, Song B W, Duan G L, Shi Y M Real-time reliability evaluation of two-phase Wiener degradation process. Communications in Statistics-Theory and Methods, 2017, 46(1): 176-188. doi: 10.1080/03610926.2014.988262 [16] Chen N, Tsui K L. Condition monitoring and remaining useful life prediction using degradation signals: Revisited. IIE Transactions, 2013, 45(9): 939-952. doi: 10.1080/0740817X.2012.706376 [17] Wang Y, Peng Y, Zi Y, Jin X H, Tsui K L. A two-stage data-driven-based prognostic approach for bearing degradation problem. IEEE Transactions on Industrial Informatics, 2016, 12(3): 924–932. doi: 10.1109/TII.2016.2535368 [18] Peng Y, Wang Y, Zi Y. Switching state-space degradation model with recursive filter/smoother for prognostics of remaining useful life. IEEE Transactions on Industrial Informatics, 2018, 15(2): 822-832. [19] Zhang J X, Hu C H, He X, Si X S, Liu Y, Zhou D H. A novel lifetime estimation method for two-phase degrading systems. IEEE Transactions on Reliability, 2018, 68(2): 689-709. [20] Wang. W, Carr M, Xu W, Kobbacy K. A model for residual life prediction based on Brownian motion with an adaptive drift. Microelectronics Reliability, 2011, 51(2): 285-293. [21] Si X S, Wang W, Hu C H, Zhou D H. Remaining useful life estimation-a review on the statistical data driven approaches. European Journal of Operational Research, 2011, 213(1): 1-14. doi: 10.1016/j.ejor.2010.11.018 [22] Zhai Q Q, Ye Z S. RUL prediction of deteriorating products using an adaptive Wiener process model. IEEE Transactions on Industrial Informatics, 2017, 13(6): 2911-2921. doi: 10.1109/TII.2017.2684821 [23] Si X S, Wang W, Chen M Y, Hu C H, Zhou D H. A degradation path-dependent approach for remaining useful life estimation with an exact and closed-form solution. European Journal of Operational Research, 2013, 226(1): 53-66. doi: 10.1016/j.ejor.2012.10.030 [24] Shumway R H, Stoffer D S. Times Series Analysis and Its Applications. New York: Springer, 2011. 326−344 [25] 夏正敏. 基于分形的网络流量分析及异常检测技术研究 [博士学位论文], 上海交通大学, 2012Xia Zheng-Min. Research on network traffic analysis and abnormal detection based on fractal theory [Ph. D. Dissertation], Shanghai jiao-tong University, 2012 [26] He W, Williard N, Osterman M, Pecht M. Prognostics of lithium-ion batteries based on Dempster–Shafer theory and the Bayesian Monte Carlo method. Journal of Power Sources, 2011, 196(23): 10314–10321. doi: 10.1016/j.jpowsour.2011.08.040 -
下载:
下载:
计量
- 文章访问数: 2364
- HTML全文浏览量: 611
- PDF下载量: 280
- 被引次数: 0
