受环境、工作载荷与自身材料、结构性能老化等因素的影响, 系统的性能会随着使用发生不可避免的退化.系统的性能退化最终会导致失效, 而失效的发生将带来经济、环境等多方面的损失, 甚至造成重大的安全事故.如果能够在系统完全失效之前进行有效预测, 并及时地进行维护和管理, 就可以有效地降低运行风险与成本, 提高系统的安全性与可靠性^[1-2].因此, 预测与健康管理(Prognostics and health management, PHM)成为当前可靠性领域的研究热点问题, 基于退化现象的分析方法在解决系统剩余寿命(Remaining useful lifetime, RUL)预测问题时取得了一系列的成果^[3-5].然而, 在剩余寿命预测时考虑控制作用的影响、对控制系统进行剩余寿命预测的研究几乎是空白.

隐含退化过程难以检测是闭环控制系统剩余寿命预测研究的一个难点.得益于特殊的闭环反馈结构, 控制系统可以补偿实际工作中产生的故障和扰动^[6], 具有一定的容错能力.当采用容错控制(Fault tolerant control, FTC)^[7-8]策略时, 闭环控制系统对故障的容错能力将进一步提升, 即使内部已经发生比较严重的性能退化, 系统输出可能仍然在正常范围之内.这使得内部部件的性能退化呈现隐含退化的特点, 给传统的寿命预测方法带来一定困难^[9].在现有文献中, 针对隐含退化过程的剩余寿命预测研究通常采用基于随机滤波的方法^[10]、基于协变量风险模型的方法^[11]和基于隐Markov或隐半Markov模型的方法^[12-13], Si等^[5]对此进行了详细的综述.无论哪种方法, 都需要获得观测值与退化状态、失效率或寿命之间的关系模型, 由于控制系统的解析模型已知, 基于随机滤波的方法不需要利用一定的历史寿命数据进行模型辨识, 相比于其他方法更具优势.文献[14]研究了考虑执行器性能退化的控制系统故障检测问题, 提出一种自适应卡尔曼滤波算法在线估计出执行器故障等级.然而, 这类采用控制效能损失^[14-18]程度对执行器退化水平建模的方法没有区分控制电路故障和执行器硬件性能退化这种早期微小故障^[19]的区别, 在此基础上设计的容错控制方法可能会增大执行器负载、加重退化程度. Langeron等^[20]将这种执行器故障模型^[14-18]的不足总结为:与控制算法相独立、与任何形式的物理退化过程相独立、故障发生于确定的时间点. Xu等^[21-22]针对动态系统提出一种实时可靠性预测方法, 选择系统模型参数作为表征系统性能的隐含退化量并采用粒子滤波算法在线辨识隐含的随机退化过程模型.然而, 由于没有考虑反馈控制的作用, 因而也就没有将退化过程模型与系统负载、控制算法联系在一起.

除了上述难点之外, 阈值大小的选择也是一个难点.由于退化过程隐含且不易观测, 因此缺乏历史失效时刻退化状态或观测值大小的经验数据, 难以根据工程经验直接给出阈值大小.文献[23]对阈值不确定这一问题进行了详尽的分析, 文献[24-26]也指出使用固定阈值是存在问题的.这一问题在闭环控制系统应用中尤其突出.除了系统内部一些部件发生性能退化会降低系统的控制品质并可能导致系统故障^[27], 控制算法对当前运行环境下的随机干扰等不良影响的鲁棒性也是影响系统可靠性的重要因素.这两个因素互相耦合, 导致在不同的工作环境和工程需求下, 控制系统的失效时刻难以通过一个固定的失效阈值来判定, 而如果失效时刻本身不够准确, 那么在此基础上获得的历史寿命数据或剩余寿命预测结果也是不准确的.基于此, 本文考虑了系统隐含退化过程的失效阈值并非固定值的情况, 提出一种在线计算失效阈值的方法.

综上分析, 当控制系统的设计需要同时满足控制性能要求和可靠性要求, 尤其是当系统不具有执行器冗余度时, 就需要在控制性能和可靠性之间进行折中^[28], 然而现有模型无法准确描述控制系统的可靠性.鉴于此, 本文提出一种考虑执行器性能退化的控制系统剩余寿命预测方法, 采用随机过程对执行器的退化水平进行建模, 并假设退化速率是一个与控制负载相关的变量.针对闭环反馈作用下部件呈现隐含退化的特点, 采用权值优选粒子滤波算法对系统输入输出的采样值进行分析, 进而估计系统状态和执行器隐含退化量.同时, 在预测时刻采用蒙特卡洛(Monte Carlo, MC)仿真以选择合理的执行器失效阈值, 通过将执行器的退化水平带入退化模型中外推出系统剩余寿命分布, 为后期健康管理活动提供可信依据.最后的仿真结果证明了本文所提方法的有效性.

1 问题描述

考虑只有一个执行器发生性能退化的单输入单输出闭环控制系统

$ \begin{equation} \label{eq1} \begin{cases} {\dot{{x}}(t) }=\mathit{\boldsymbol{f}}\!\left( {t, \phi, \mathit{\boldsymbol{x}}(t), u(t) } \right)\mathit{\boldsymbol{+w}}\!\left( t \right) \\ y(t) =h\!\left( {t, \phi, \mathit{\boldsymbol{x}}(t) } \right)+v(t) \\ \end{cases} \end{equation} $

(1)

式中, $ \mathit{\boldsymbol{x}}(t) \in {\bf R}^{n} $为$ n $维状态向量, $ u(t) $为控制输入, $ y(t) $为系统输出, $ \phi $为执行器的控制效能, 是隐含的性能参数, 难以直接测得. $ \mathit{\boldsymbol{w}}(t) $和$ v(t) $为合适维度、互相独立的过程噪声向量和观测噪声向量, 服从高斯或非高斯分布. $ \mathit{\boldsymbol{f}} $为状态向量$ \mathit{\boldsymbol{x}}\!\left( t \right) $和输入$ u(t) $的非线性函数, $ h $为状态向量$ \mathit{\boldsymbol{x}}(t) $的线性或非线性函数.输出残差定义为$ e(t) =y^{{\rm ref}}(t) -y(t) $, 表示系统期望输出与实际输出的偏差, $ y^{{\rm ref}}(t) $为参考输出.

图 1展示了执行器性能退化过程与闭环控制系统失效过程的关系, 从图 1中可以看出, 执行器的控制效能连续缓慢下降, 导致执行器的实际输出受到影响. 图 1左上角输出残差$ e(t) $的变化曲线, 起初波动的幅度较小, 表示在闭环反馈的作用下, 系统尚能补偿部分外部干扰的影响, 输出残差仍然保持在在性能指标允许的范围内, 但此时系统的控制品质已经有所下降.随着隐含退化量$ \phi (t) $变化曲线不断下降, $ e(t) $曲线突然持续上升直至超过残差的允许上限, 表示当执行器性能退化累积到一定程度时, 系统会在外部随机干扰的影响下发生失效.本文选择合适的模型参数$ \phi (t) $作为隐含退化量, 相比于传统性能变量有以下几点优势: 1)可以通过监测系统输入输出数据辨识出来, 不需要额外增加传感设备去监测性能变量; 2)隐含退化量通过系统模型作用于系统输出, 其大小直接反映系统的实际工作状态和性能水平, 不需要再通过机理分析和大量实验去建立性能参数和系统可靠性之间的关系; 3)其退化阈值可以通过系统解析模型和输出性能指标进行确定, 不同工作环境下只需要相应地修改系统模型, 即可求得当前环境下的阈值, 而传统性能参数的阈值需要在固定的环境中对同类系统进行大量实验, 或者通过大量工程经验的积累给出, 且在新的工作环境下基于固定阈值无法得到准确的系统剩余寿命预测结果; 4)模型参数通过闭环反馈作用与控制器耦合, 选择模型参数作为退化量, 便于建立与控制算法、系统负载的关系.

由于本文仅考虑执行器的退化过程, 控制器的形式一经确定便不再改变.选择工业和过程控制领域常用的比例-积分-微分(Proportional integral differential, PID)控制算法作为闭环控制系统控制器的设计依据.定义$ u(t) $为$ t $时刻的控制器输出, 则PID控制算法具有如下标准形式

$ \begin{equation} \label{eq2} u(t) \!=\!K_{{\rm P}} \left[{e(t) \!+\!\frac{1}{T_{{\rm I}} }\int_0^t {e( \tau ){\rm d}\tau } +T_{{\rm D}} \frac{{\rm d}e(t) }{{\rm d}t}} \right] \end{equation} $

(2)

式中, $ K_{{\rm P}} $为比例系数, $ T_{{\rm I}} $为积分时间常数, $ T_{{\rm D}} $为微分时间常数.关于这三个参数的选择与整定问题在许多参考书籍中都有所介绍, 这里不再赘述.

控制信号$ u(t) $被用来驱动执行器产生响应的输出来维持系统平衡, 执行器的输出是一个状态量, 设为$ x_{{\rm a}} $, 可以是电机输出的力矩、阀门开闭的大小等. $ t $时刻执行器的输出定义如下

$ \begin{equation} \label{eq3} x_{{\rm a}} (t) =g\!\left( {u(t), \phi } \right) \end{equation} $

(3)

式中, $ g $是一个由系统状态方程决定的函数.由式(3)可见执行器的实际输出由控制律$ u\!\left( t \right) $和执行器自身的控制效能$ \phi $共同决定.

文献[29-30]将控制效能损失视作一类执行器故障, 由于在连续负载下工作, 执行器的机械参数或电气参数发生缓慢变化, 进而导致其输出效能是一个渐变的、连续的退化过程.因此, 采用连续的随机过程对执行器的性能退化过程建模, 是一个合理的选择.常见的退化过程模型包括扩散过程模型^[31-32]、Gamma过程模型^[33]以及隐Markov模型^[13]等.不论采用哪种退化模型都是对系统真实性能退化过程的近似描述, 其中的扩散过程作为一种连续的随机过程, 近年来被广泛应用于机电系统的退化建模, 适用于非单调、非线性的退化过程描述.因此, 本文采用扩散过程对执行器的控制效能退化建模, 并假设退化过程的参数已经通过辨识得到.令$ \phi (t) $表示$ t $时刻执行器的退化量, 则基于扩散过程的执行器退化过程$ \left\{ {\phi (t), t\geq 0} \right\} $可以表示为

$ \begin{equation} \label{eq4} \phi (t) =\phi ( 0 )+\int_0^t {\mu( {\tau, \mathit{\boldsymbol{x}}} )} {\rm d}\tau +\sigma_{{\rm B}} B(t) \end{equation} $

(4)

式中, $ \phi ( 0 ) $表示隐含退化量的初始值; $ \sigma_{{\rm B}} $为扩散系数; $ B(t) $表示标准Brownian运动, 且$ B(t)\sim {\rm N}( {0, t} ) $; $ \mu( {\tau, \mathit{\boldsymbol{x}}} ) $为漂移系数, 用来表征漂移速度, 其中的$ \mathit{\boldsymbol{x}} $表示漂移速度受到系统负载和工作状态的影响.

在进行剩余寿命预测之前首先要对系统失效进行定义.文献[34-35]将系统失效定义为系统丧失一定功能以至于无法继续完成预定任务, 认为这种任务需求可以通过某种阈值来表示, 而这个阈值可以是系统状态、参数或是输出的函数.同时, 文献[36]认为通过将系统真实输出与某种指标比较得到残差并对残差进行评价是得到故障检测结果的常用手段.对于图 1所示闭环跟踪系统, 预定任务就是使得系统输出有效跟踪预设的参考输出并将输出残差控制在一定范围之内, 而当执行器性能退化到一定程度, 系统从早期微小故障状态向显著故障^[19]转变的时刻即为失效时刻.基于上述分析, 本文方法采用残差上界$ e_{\max } $作为判断系统是否失效的性能指标, 并将系统的使用寿命定义为输出残差$ \left\{ {\left| {e( t )} \right|, t\geq 0} \right\} $达到$ e_{\max } $的首达时间

$ \begin{equation} \label{eq5} T_{{\rm e}} =\inf \left\{ {t:\left| {e(t) } \right|\geq e_{\max } \left| {\left| {e(0)} \right|<e_{\max } } \right.} \right\} \end{equation} $

(5)

式中, $ {\rm inf} $为下确界算子, $ e_{\max } $为残差阈值, 是一个固定值, 根据工程中的实际要求给出, 这里阈值的设置需要满足一定的大小要求, 避免输出残差在正常工作下产生波动达到阈值而造成误报和虚警的情况. $ \left| {e(t) } \right| $实际上是残差评价函数^[36]针对单输入单输出跟踪系统(1)的简化形式, 当系统为多输入多输出系统且具有其他性能指标时, 只需将$ \left| {e(t) } \right| $替换成其他合理的残差评价函数即可.

由图 1可知, 从开始时刻起很长一段时间内系统的输出残差变化较小, 不便于描述系统实时的健康水平, 因此在预测时对隐含退化量建模并提出一种基于变失效阈值的剩余寿命预测方法.基于该方法的系统寿命预测结果依赖于运行过程中对未来退化量的估计数据, 通常定义为未来退化量估计值$ \left\{ {\hat{{\phi }}(t), t\geq 0} \right\} $到达失效阈值$ \phi _{\min } $的首达时间

$ \begin{equation} \label{eq6} T_{\phi } =\inf \left\{ {t:\hat{{\phi }}(t) \leq \phi _{\min } \left| {\hat{{\phi }}(0)>\phi_{\min } } \right.} \right\} \end{equation} $

(6)

式中, $ \hat{{\phi }}(t) $表示系统在预测时刻对未来退化量的估计值.

这里出现了式(5)和式(6)两种寿命定义, 对此说明如下:式(5)用作真实退化过程中系统是否失效的判断依据, 而式(6)用作预测计算过程中系统是否失效的判断依据.从图 1中可以看出, 由于性能参数的缓慢退化会导致控制器在某个时刻无法补偿随机干扰, 系统发散.因此在失效时刻前后, 系统输出残差变化非常剧烈.反过来说, 一定程度的残差上界不准确并不会引起失效时刻产生太大偏差, 这也是本文选取式(5)作为判断系统是否失效或发生显著故障依据的原因.这样一来, 失效时刻相对比较固定, 在不同大小的失效阈值下得到的预测结果才有了一个可信的参照, 即真实失效时间(对应于真实剩余寿命), 本文基于变阈值预测方法与传统基于固定阈值预测方法下的预测结果才有了可比性.

基于以上分析, 本文主要研究以下2个问题:

1) 如何在已知控制系统执行器退化模型的情况下确定合理的失效阈值, 建立基于隐含退化量失效阈值的系统失效判断准则;

2) 如何利用在线观测输入输出数据对隐含状态进行估计并实时预测系统的剩余寿命分布.

2 剩余寿命预测方法

提出一种考虑退化速率受负载影响的基于扩散过程的非线性退化过程建模方法, 在已知系统状态模型和噪声模型的情况下, 在每一个预测时刻$ t_{1}, t_{2}, \cdots, t_{n} $实时计算系统的剩余寿命.该方法包含: 1)在每一个预测时刻利用MC仿真合理设置隐含退化量的失效阈值, 建立基于隐含退化量变失效阈值的失效判断准则; 2)基于权值优选粒子滤波算法, 根据在线监测的系统输入输出数据实时估计出系统的隐含状态; 3)在不同的预测时刻计算系统剩余寿命分布.

2.1 基于隐含退化量变失效阈值的失效判断准则

首先定义系统基于残差准则的失效区域为$ {\cal F} $, $ {\cal F} $包含系统所有状态变量.系统输出偏差$ e $首次进入区域$ {\cal F} $即意味着系统在满足一定使用性能的意义下发生失效.基于上述定义和式(5), 假设系统首次进入$ {\cal F} $的采样时刻为$ t_{{\rm est}} $, 则系统在时刻$ t $的剩余寿命可以描述如下:

$ \begin{equation} \label{eq7} \begin{aligned} T_{{\rm RUL}} &=\inf \left\{ {t_{{\rm est}} \geq t, e(t_{{\rm est}} )\in {\cal F}} \right\}-t =\\ &\quad~ T_{e} -t\left| {T_{e} } \right.>t \\ \end{aligned} \end{equation} $

(7)

基于隐含退化量失效阈值的系统失效区域为$ \Phi $, 系统隐含退化量$ \phi $首次进入区域$ \Phi $即意味着系统的可靠性大幅下降, 进入故障易发区.基于上述定义和式(6), 假设系统在MC仿真中首次进入$ \Phi $的采样时刻为$ t_{{\rm est}}^{\ast}$, 则系统在时刻 $ t $的剩余寿命可以描述如下:

$ \begin{equation} \label{eq8} \begin{aligned} T_{{\rm RUL}}^{\ast } &=\inf \left\{ {t_{{\rm est}}^{{\rm \ast }} \geq t, \phi {\rm (}t_{{\rm est}}^{\ast } {\rm )}\in \Phi } \right\}-t =\\ &\quad~ T_{\phi } -t\left| {T_{\phi } } \right.>t \\ \end{aligned} \end{equation} $

(8)

式中, $ t_{{\rm est}}^{\ast } $为所有MC仿真轨迹中首次达到残差上界$ e_{\max } $的时刻, 是确定变失效阈值的参考时刻.根据式(6)系统寿命的定义, 系统在采样时刻$ t_{{\rm est}}^{{\rm \ast }} $还未失效的累积概率可以描述如下:

$ %\begin{equation} % \label{eq9} \begin{align} &\Pr \!\left( {T_{\phi } >t_{{\rm est}}^{{\rm \ast }} \left| {Y_{1} =y_{1}, \cdots } \right., Y_{n} =y_{l} } \right) =\\\nonumber &\int {R_{\phi }\!\left( {t_{{\rm est}}^{{\rm \ast }} } \right)\mu_{y_{1}, \cdots, y_{l} }\!\left( {\rm d}\mathit{\boldsymbol{z}} \right)} \approx \sum\limits_{i=1}^{N_{{\rm p}} } {R_{\phi }\!\left( {t_{{\rm est}}^{\ast } } \right)} \mu_{y_{1}, \cdots, y_{l} } \end{align} %\end{equation} $

(9)

式中, $ N_{p} $表示优选粒子数, $ \mu_{y_{1}, \cdots, y_{l} }({d\mathit{\boldsymbol{z}}}) $表示在$ t_{{\rm est}}^{{\rm \ast }} $时刻获得观测数据$ y_{1}, \cdots, y_{l} $轨迹的重要性函数; $ R_{\phi } ( {t_{{\rm est}}^{{\rm \ast }} }) $表示系统在时间区间$ \left( {0, t_{{\rm est}}^{\ast } } \right] $内基于隐含退化量失效阈值准则未失效的概率.式(9)是本文基于隐含退化量失效阈值求解系统可靠度和剩余寿命概率密度函数(Probability density function, PDF)的理论依据.

在不同预测时刻进行MC仿真, 选取式(8)中首达定义下的参考时刻$ t_{{\rm est}}^{{\rm \ast }} $作为失效阈值$ \phi_{\min } $的对应时刻, 则$ \phi _{\min } $通过下式计算

$ \begin{equation} \label{eq10} \phi_{\min } =\sum\limits_{i=1}^{N_{{\rm traj}} }\frac{{\phi_{i}\!\left( {t_{{\rm est}}^{{\rm \ast }} } \right)} }{{N_{\rm traj} }} \end{equation} $

(10)

式中, $ N_{{\rm traj}} $表示仿真次数.由于不同预测时刻的初值不同, MC仿真轨迹也不同, 因此失效阈值不是一个固定值.下面通过图 3进一步具体说明在预测时采用变失效阈值的必要性.

真实的控制系统由于随机干扰带来的时变不确定性, 可能会出现不同的退化轨迹, 如图 2中退化轨迹1和退化轨迹2所示.由于预测时刻不同, MC仿真的初值不同, 系统未来退化轨迹的预期也就不同. 图 2中故障易发区就是在第2节问题描述中阐述的系统已经达到控制约束上界或鲁棒性较差的区域.如果采用固定的失效阈值$ \phi _{\min } $, 则在$ t_{n} $时刻对系统2的预测结果为$ RUL_{2} $. 可以看到, 系统2 已经处在一个加速退化的阶段, 有较大的概率在运行到 $ t_{{\rm est2} }$之前就已经失效, 因此这样预测得到的剩余寿命结果是不符合实际的.而采用变失效阈值可以根据在线监测数据调整系统未来退化轨迹的预期, 得到更加接近真实的预测结果.

2.2 隐含退化量和系统状态的联合估计

为了方便对系统状态进行滤波估计, 首先对控制系统(1)、控制律(3)和退化量(5)进行欧拉离散化处理

$ \begin{equation} \label{eq11} \begin{cases} \mathit{\boldsymbol{x}}_{k+1} =f( {\mathit{\boldsymbol{x}}_{k}, \phi_{k}, u_{k} } )+ \mathit{\boldsymbol{w}}_{k} \\ \phi_{k+1} =\phi_{k} +\mu( {\mathit{\boldsymbol{x}}_{k} } )+\sigma _{B} B(t_{k} ) \\ y_{k} =h( {\mathit{\boldsymbol{x}}_{k} } )+v_{k} \\ \end{cases} \end{equation} $

(11)

$ \begin{equation} \label{eq12} \begin{aligned} u_{k} &=u_{k-1} +K_{{\rm P}} \!\left[{\left( {1+\frac{\Delta t}{T_{{\rm I}} }+\frac{T_{{\rm D}} }{\Delta t}} \right)e_{k} } \right. +\\ &\quad~ \left.\left( {-1-\frac{2T_{{\rm D}} }{\Delta t}} \right)e_{k-1} +\frac{T_{{\rm D}} }{\Delta t}e_{k-2} \right] \\ \end{aligned} \end{equation} $

(12)

从式(11)可以看出, 隐含退化量$ \phi_{k} $原来是系统的一个结构参数, 现在系统经过扩维$ \phi_{k} $成为系统的一个状态量, 对隐含退化量的估计问题变成状态变量与参数的联合估计问题.

本文参照文献[37]将系统原状态量$ \mathit{\boldsymbol{x}} $与隐含性能变量$ \phi $扩展成为扩维状态向量$ \mathit{\boldsymbol{z}} $.

$ \begin{equation} \label{eq13} \mathit{\boldsymbol{z}}_{k} =\begin{pmatrix} \mathit{\boldsymbol{x}}_{k} \\ \phi_{k} \\ \end{pmatrix} \end{equation} $

(13)

考虑到退化过程是一个非线性扩散过程, 且退化速率与系统状态量之间有耦合, 因此本文采用基于递归贝叶斯估计理论的粒子滤波(Particle filter, PF)算法解决扩维状态变量的联合估计问题.粒子滤波对系统模型的要求较为宽松, 适用于强非线性非高斯系统, 不过粒子权值退化和多样性匮乏问题是难以完全消除的一对矛盾^[38].粒子滤波算法有许多种改进和变形, 本文采用在著名的序贯重要性重采样(Sequential importance resampling, SIR)算法^[39]基础上改进的基于权值优选的粒子滤波器^[40]以克服粒子贫化的问题.

基于权值优选的粒子滤波算法^[40]的基本思想是:从众多备选粒子中选出权值较大的粒子用来参与将来时刻系统状态估计的迭代运算, 这样剩下的大量粒子在不影响滤波精度的情况下可以保持粒子的多样性.算法流程如下:

1) 初始化

对$ \forall i=1, \cdots, N_{{\rm s}} $, 从初始分布(一般采用均匀分布作为初始分布) $ \pi\! \left( {x_{0} } \right) $中抽取$ N_{{\rm s}} $个粒子构成初始粒子群$ \left\{ {\mathit{\boldsymbol{z}}_{0}^{(i)}, w_{0}^{(i)} } \right\}_{i=1}^{N_{{\rm s}} } $;

2) 状态更新

已知$ \left\{ {\mathit{\boldsymbol{z}}_{k-1}^{(i)}, w_{k-1}^{(i)} } \right\}_{i=1}^{N_{{\rm s}} } $, 根据$ k-1 $时刻的粒子$ \mathit{\boldsymbol{z}}_{k-1}^{(i)} $和系统的状态方程(11)计算$ k $时刻$ {\tilde{{z}}}_{k}^{(i)} {\rm \sim }\, p\!\left( {\mathit{\boldsymbol{z}}_{k} \left| {\mathit{\boldsymbol{z}}_{k-1}^{(i)} } \right.} \right) $;

3) 权值更新

根据$ k $时刻的观测值$ \mathit{\boldsymbol{y}}_{k} $更新对应的$ N_{{\rm s} } $个粒子的权值

$ \begin{equation} \label{eq14} w_{k}^{(i)} =w_{k-1}^{(i)} \frac{p\!\left( {\mathit{\boldsymbol{y}}_{k} \left| {{\tilde{{z}}}_{k}^{(i)} } \right.} \right)p\!\left( {{\tilde{{z}}}_{k}^{(i)} \left| {\mathit{\boldsymbol{z}}_{k-1}^{(i)} } \right.} \right)}{q\!\left( {{\tilde{{z}}}_{k}^{(i)} \left| {\mathit{\boldsymbol{z}}_{0:k-1}^{(i)}, \mathit{\boldsymbol{y}}_{1:k} } \right.} \right)} \end{equation} $

(14)

式(14)展示了粒子滤波又被称为序贯MC方法中序贯的含义.

4) 粒子选优

对这$ N_{{\rm s}} $个粒子$ {\tilde{{z}}}_{k}^{(i)} $按照权值$ w_{k}^{(i)} $的大小进行排序, 选出前面的$ N_{{\rm p}} $个粒子;

5) 权值归一化

对选出的$ N_{{\rm p}} $个粒子进行权值归一化处理:

$ w_{k}^{(i)} =\frac{w_{k}^{(i)} }{\sum\limits_{i=1}^{N_{{ p}} } {w_{k}^{(i)} } } $

6) $ k $时刻状态估计

$ \begin{equation} \label{eq15} \mathit{\boldsymbol{z}}_{k}^{\ast } \approx \sum\limits_{i=1}^{N_{{\rm p}} } {w_{k}^{(i)} } \times {\tilde{{z}}}_{k}^{(i)} \end{equation} $

(15)

7) 权值恢复

先将选出来的$ N_{{\rm p}} $个粒子的权值恢复到归一化之前: $ w_{k}^{(i)} =w_{k}^{(i)} \times \sum\nolimits_{i=1}^{N_{{\rm p}} } {w_{k}^{(i)} } $, 再对所有$ N_{{\rm s}} $个粒子进行权值归一化:

$w_{k}^{(i)} =\frac{w_{k}^{(i)} }{\sum\limits_{i=1}^{N_{{\rm s}} } {w_{k}^{(i)} } } $

8) 重采样

如果$ \hat{{N}}_{{\rm eff}} =1/{\sum\nolimits_{i=1}^{N_{{\rm s}} }\! {\left( {w_{k}^{(i)} } \right)} }^{2}<N_{{\rm thresh}} $, 采样自举法对所有$ N_{{\rm s}} $个粒子进行重采样, 并将所有粒子的权值恢复到初始分布: $ \left\{ {{\tilde{{z}}}_{k}^{(i)}, w_{k}^{(i)} } \right\}_{i=1}^{N_{{\rm s}} } \Rightarrow \left\{ {\mathit{\boldsymbol{z}}_{k}^{(i)}, 1/{N_{{\rm s}} }} \right\}_{i=1}^{N_{{\rm s}}} $;

9) 求得状态量的概率分布

实时更新状态量 $ \mathit{\boldsymbol{z}}_{k} $的后验概率分布:

$ \begin{equation} \label{eq16} \begin{aligned} p\!\left( {\mathit{\boldsymbol{z}}_{k} \left| {\mathit{\boldsymbol{y}}_{1:k} } \right.} \right)&=\hat{{p}}\!\left( {\mathit{\boldsymbol{z}}_{k} \left| {\mathit{\boldsymbol{y}}_{1:k} } \right.} \right)\approx \\ & \sum\limits_{i=1}^{N_{{\rm p}} } {w_{k}^{(i)} } \delta _{\mathit{\boldsymbol{z}}_{k}^{(i)} } \!\left( {\mathit{\boldsymbol{z}}_{k} -\mathit{\boldsymbol{z}}_{k-1}^{(i)} } \right) \end{aligned} \end{equation} $

(16)

返回步骤2)进行循环迭代.

图 3展示了基于粒子滤波的扩维状态估计和系统剩余寿命在线预测的过程. 图 3中剩余寿命PDF是在$ t_{n} $时刻$ N_{{\rm p}} $个优选粒子的基础上进行$ N_{{\rm traj}} $次MC仿真的结果, 具体方法见第3.3节.

2.3 剩余寿命分布的求解

由于预测时刻进行MC仿真的计算量较大, 假设系统仿真步长和粒子滤波采样间隔为$ \Delta t_{{\rm s}} =0.001\, {\rm s} $, 那么可以每隔$ \Delta t_{{\rm est}} =2\, 000\Delta t_{{\rm s}} =2\, {\rm s} $进行剩余寿命在线预测.

假设初始时刻为$ t=0\, {\rm s} $, 在系统进入失效区域$ \Phi $之前, 预测时刻是间隔为$ \Delta t_{{\rm est}} $的等间隔采样点$ t_{1}, t_{2}, \cdots, t_{n} $. $ t_{n} $时刻系统剩余寿命的计算步骤如下:

1) 确定初始状态及仿真次数

通过式(15)计算得到$ t_{n} $时刻优选粒子的扩维状态估计值$ \mathit{\boldsymbol{z}}_{t_{n} }^{\ast } $, 假设所有优选粒子总共产生$ N_{{\rm traj}} $条MC仿真轨迹, 粒子权重$ w_{t_{n} }^{(i)} $越大, 每个粒子分配到的仿真次数$ N_{i} $也越大;

2) 利用自举法产生仿真轨迹并分配初始状态

For $ j=1, \cdots, N_{{\rm traj}} $

a) 生成随机数: $ u_{j} \, {\rm \sim }\, U\!\left( {0, 1} \right) $;

b) 初始化累计概率$ w_{j} =0 $;

c) For $ i=1, \cdots, N_{{\rm p}} $, $ w_{j} =w_{j} +w_{t_{n} }^{(i)} $, 当$ w_{j} \geq u_{j} $时, $ \mathit{\boldsymbol{z}}_{j} =\mathit{\boldsymbol{z}}_{t_{n} }^{(i)} $, 退出循环;

3) MC仿真

以$ \mathit{\boldsymbol{z}}_{j} $为状态初值, 对第$ j $条轨迹按照式(11)和式(12)进行仿真, 当$ \phi <\phi _{\min } $时, 保存此时的采样时刻, 最终所有轨迹的仿真结果即为系统剩余寿命的数值PDF.

3 仿真实验与结果分析 3.1 仿真实验

仿真采用惯性平台稳定回路控制系统^[41]对本文提出的剩余寿命预测方法进行验证.如图 4所示, 该控制系统的传感器由平台上的敏感元件组成, 如陀螺仪、加速度计等; 控制器由框架角传感器和位置补偿控制器共同组成; 执行器直接作用于产生系统输出的被控对象, 此处执行器为力矩电机.

系统模型如下^[41]:

$ \begin{equation} \label{eq17} \begin{cases} L_{{\rm m}} \dfrac{{\rm d}i_{{\rm m}} }{{\rm d}t}+R_{{\rm m}} i_{{\rm m}} +K_{{\rm e}} \omega_{{\rm p}} =U_{{\rm m}} \\ \left( {J_{{\rm m}} +J_{{\rm L}} } \right)\dfrac{{\rm d}\omega_{{\rm p}} }{{\rm d}t}=K_{{\rm m}} i_{{\rm m}} \\ \dfrac{{\rm d}\theta }{{\rm d}t}=\omega_{{\rm p}} \\ \end{cases} \end{equation} $

(17)

式中, $ L_{{\rm m}} $为电机电枢电感; $ R_{{\rm m}} $为电枢绕组; $ i_{{\rm m}} $为电机电枢电流; $ K_{{\rm e}} $为电机反电势系数; $ \omega_{{\rm p}} $为框架角速度; $ U_{{\rm m}} $为电机控制电压; $ J_{{\rm m}} $和$ J_{{\rm L}} $分别为电机转子转动惯量和负载转动惯量; $ K_{{\rm m}} $为电机的力矩系数; $ \theta $为框架角.

由于执行器在工作过程中存在控制效能损失的现象, 系统模型中的模型参数也会随之变化, 而力矩系数中包含了电机的内部结构特征和性能水平, 因此假设隐含退化量为力矩系数$ K_{{\rm m}} $, 则退化过程可以描述如下:

$ \begin{equation} \label{eq18} K_{{\rm m}} (t) =K_{{\rm m}} (0)+\int_0^t {\mu \!\left( {\tau ;i_{{\rm m}} } \right)} {\rm d}\tau +\sigma_{{\rm B}} B(t) \end{equation} $

(18)

式中, $ \mu \!\left( {\tau ;i_{{\rm m}} } \right) $是受电机的负载即电枢电流影响的退化漂移系数, 是一个变量.这种将执行器的工作应力与其所受负载联系起来的假设是符合实际的, 而负载越大退化速率越快也符合客观规律.文献[20]采用控制量平方来量化执行器的工作应力, 鉴于此, 假设退化漂移系数具有如下形式:

$ \begin{equation} \label{eq19} \mu \!\left( {\tau ;i_{{\rm m}} } \right)=\lambda i_{{\rm m}}^{2}\!\left( \tau \right) \end{equation} $

(19)

式中, $ \lambda $为退化速率系数, 由式(19)可知, 单位时间的退化量与电枢电流的平方成正比, 代入式(18)可得:

$ \begin{equation} \label{eq20} K_{{\rm m}} (t) =K_{{\rm m}} (0)+\lambda \int_0^t {i_{{\rm m}}^{2}\!\left( \tau \right)} {\rm d}\tau +\sigma_{{\rm B}} B(t) \end{equation} $

(20)

数字式稳定回路控制系统建立在采样数据的基础上, 利用采样数据可以方便地对隐含退化量进行粒子滤波估计.因此, 对稳定回路采用数字PID控制, 且对控制系统仿真和粒子滤波选择相同的采样周期.为了描述和编程方便, 令$ x_{1} =i_{{\rm a}} $, $ x_{2} =\omega_{{\rm p}} $, $ x_{3} =\theta $, $ x_{4} =K_{{\rm m}} $, $ u=U_{{\rm m}} $, $ y=x_{3} $, $ J=J_{{\rm m}} +J_{{\rm L}} $, 对系统(17)进行离散化, 则整个闭环控制系统的完整描述如下:

$ \begin{align} t_{k} &=t_{k-1} +\Delta t\\ y_{k} &=x_{3_{k} } \\\nonumber y_{k}^{{\rm ref}} &=\frac{\pi }{2}\sin \!\left( {0.001\cdot t_{k} } \right)\\ \end{align} $

(21)

$ \begin{align}e_{k} &=y_{k}^{{\rm ref}} -x_{3_{k} }\\ u_{k} &=u_{k-1} +K_{{\rm P}}\! \left[{\left( {1+\frac{\Delta t}{T_{{\rm I}} }+\frac{T_{{\rm D}} }{\Delta t}} \right)e_{k} } +\right. \\\nonumber &\quad~ \left. {\left( {-1-\frac{2T_{{\rm D}} }{\Delta t}} \right)e_{k-1} +\frac{T_{{\rm D}} }{\Delta t}e_{k-2} } \right] \\ \end{align} $

(22)

$ \begin{align} \label{eq22} u_{k} &=\begin{cases} u_{\max } :u_{k} >u_{\max } \\ u_{k} :u_{\min } \leq u_{k} \leq u_{\max } \\ u_{\min } :u_{k} <u_{\min } \\ \end{cases}\\ \end{align} $

(23)

$ \begin{align} \label{eq23} x_{1_{k+1} } &=\left( {-\frac{R_{{\rm m}} \Delta t}{L_{{\rm m}} }+\Delta t+1} \right)x_{1_{k} } + \\\nonumber &\quad~ \left( {-\frac{K_{{\rm e}} \Delta t}{L_{{\rm m}} }} \right)x_{2_{k} } +\frac{\Delta t}{L_{{\rm m}} }u_{k} +Q_{1_{k} } \\ \end{align} $

(24)

$ \begin{align} \label{eq24} x_{2_{k+1} } =\frac{\Delta t}{J}x_{1_{k} } x_{4_{k} } +x_{2_{k} } +Q_{2_{k} } \\ \end{align} $

(25)

$ \begin{align} \label{eq25} x_{3_{k+1} } =\Delta t\cdot x_{2_{k-1} } +x_{3_{k-1} } +Q_{3_{k} } \\ \end{align} $

(26)

$ \begin{align} \label{eq26} x_{4_{k+1} } =\lambda \Delta t\cdot \left( {x_{1_{k} } } \right)^{2}+x_{4_{k} } +\sigma_{{\rm B}} B_{k} \end{align} $

(27)

其中, 仿真步长取$ \Delta t=\Delta t_{{\rm s}} =0.001\, {\rm s} $; $ Q_{3} $既是状态量$ x_{3} $的过程噪声, 也是输出量$ y $的观测噪声, 因此式(21)中不再考虑额外的测量误差.过程噪声$ Q_{1} $、$ Q_{2} $、$ Q_{3} $以及扩散运动$ B $均为独立高斯随机变量且均值为0, 方差分别为$ q_{1} $、$ q_{2} $、$ q_{3} $: $ Q_{1} \, \sim \, {\rm N}( {0, q_{1} } )$, $Q_{2}\, \sim\, {\rm N}( {0, q_{2} } )$, $Q_{3}\, \sim \, {\rm N}( {0, q_{3} }) $.

在初始时刻, 首先按照表 1对各个参数以及状态变量进行初始化设置.在到达预测时刻之前的每一个采样时刻, 同时对系统进行仿真和粒子滤波, $ x_{1} $、$ x_{2} $、$ x_{3} $、$ x_{4} $既是系统仿真轨迹中的状态变量, 也是粒子滤波过程中的估计量, 对其按照式(23) $ \sim $ (26)分别进行迭代运算, 同时按照第2.2节步骤2)至步骤9)进行粒子滤波.此时, 每个状态量的粒子分布即是该状态量概率分布的数值逼近结果.当到达预测时刻时, 按照第2.3节步骤1) ~步骤3)对所有优选粒子进行MC仿真并得到每条仿真轨迹的剩余寿命, 即为系统剩余寿命的数值PDF.

实验中粒子规模取为$ N_{{\rm s}} =500 $, 优选粒子数取为$ N_{{\rm p}} =40 $, MC仿真次数取为$ N_{{\rm traj}} =1\, 000 $, 当输出残差$ \left| e \right|\geq 0.28\;{\rm rad} $时认为系统基于输出残差准则失效.

3.2 结果分析

图 5 (a)~(f)展示的是系统的真实失效过程以及粒子滤波采样结果.从图 5中可以看出, 在闭环反馈的作用下, 系统在40 s之前一直保持稳定工作, 系统输出能够有效跟踪参考信号, 然而力矩系数已经发生退化.此时, 通过对隐含退化量设置合理的阈值来预测失效时刻, 相比基于输出信号分析的寿命预测方法更加有效.从图 5 (c)中可以明显地看出电枢电流(即状态$ x_{1} $)受到三种过程噪声的共同影响, 粒子滤波值对含噪声的真实电流值跟踪效果良好.

最后一步是在间隔$ \Delta t_{{\rm est}} $的各个预测时刻$ t_{1}, t_{2}, \cdots, t_{n} $, 分别对系统的剩余寿命分布进行预测.为进一步比较不同时刻预测结果的准确性, 本文选取$ T_{{\rm prog}} =t_{10} =10\, {\rm s} $和$ T_{{\rm prog}} =t_{30} =30\, {\rm s} $两个时刻的剩余寿命预测结果进行对比.首先得到首达失效时刻$ t_{\phi } =44.2\, {\rm s} $和$ t_{\phi } =45.5\, {\rm s} $时退化量的分布, 并由式(10)计算得到的失效阈值, 分别为$ \phi _{\min } ={\rm 0.1614} $和$ \phi_{\min } =0.1379 $.

图 6表示在不同时刻系统剩余寿命分布的预测结果.由图 6可见, 距离系统失效时刻越近, 剩余寿命分布就越集中, 预测均值相对于系统真实剩余寿命的差距也越小, 意味着预测的不确定性降低, 精度提高.同时, 基于变阈值失效准则的系统剩余阈值寿命与基于输出残差失效准则的真实剩余寿命也越接近, 这是因为采用的变阈值方法充分利用了在线数据, 对失效阈值起到了在线调整以逼近真实失效过程的作用.

图 7是系统在不同预测时刻、不同失效阈值下的可靠度直方图, 可见预测时刻越接近失效时刻, 基于变阈值的可靠度曲线越陡峭, 这说明系统由于控制效能损失导致控制冗余度下降, 在$ t=30\, {\rm s } $比在 $ t=10\, {\rm s} $更易受到随机干扰的影响而发生失效.由图可知, 越往后可靠度斜率越大、退化过程不确定性越小, 退化速率越快, 这也与非线性退化过程的特点吻合.

为进一步定量比较不同时刻剩余寿命预测结果的精度, 引入剩余寿命的均方误差:

$ \begin{align} \label{eq27} MSE_{e} =\sum\limits_{i=1}^{N_{{\rm traj}} } \frac{{\left( {T_{_{{\rm RUL}_{K_{{\rm m}} } } }^{(i)} -T_{e} } \right)^{2}} }{{N_{{\rm traj}} }} \end{align} $

(28)

$ \begin{align} \label{eq28} MSE_{K_{{\rm m}} } =\sum\limits_{i=1}^{N_{{\rm traj}} } \frac{ {\left( {T_{_{{\rm RUL}_{K_{{\rm m}} } } }^{(i)} -T_{K_{{\rm m}} } } \right)^{2}} } {{N_{{\rm traj}} }} \end{align} $

(29)

式中, $ T_{e} $为系统基于输出残差失效准则的真实剩余寿命, $ T_{K_{{\rm m}} } $为系统基于预测时刻$ t_{k} $失效阈值准则的剩余阈值寿命; $ T_{_{{\rm RUL}_{k} } }^{(i)} $为第$ i $条仿真轨迹所得的剩余寿命预测值. $ MSE_{e} $用来表示基于隐含退化量失效阈值准则的预测结果在描述系统真实健康水平时的准确性, $ MSE_{K_{{\rm m}}} $ 用来表示在基于隐含退化量失效阈值准则的前提下, 预测结果与真实剩余阈值寿命的偏差, 表征预测精度. 为方便对比基于变阈值方法与基于固定阈值方法在控制系统剩余寿命预测应用中 的区别, 选择0.12205 作为固定失效阈值, 将4 个预测时刻的计算结果列于表 2 和表 3. 表中, 真实剩余寿命是指在输出残差准则下系统的真实剩余寿命, 剩余阈值寿命是指在失效阈值准则下系统 的剩余寿命, 是对系统真实剩余寿命的近似.

从表 2中可以看出, 基于固定阈值方法在不同时刻的预测精度依赖于所选固定阈值的大小和准确性, 其中表征预测精度的指标$ MSE_{K_{{\rm m}}} $ 不具有随着时间推移而减小的特性, 说明在固定阈值的框架下, 基于隐含退化量失效阈值的预测方法不能通 过在线积累数据而对未来退化轨迹进行有效修正, 或者说当失效阈值这个判断标准本身不准确时, 预期退化轨 迹精度提高的利好被不准确的判断标准抵消了. 而在表 3 中, $ MSE_{K_{{\rm m}} } $和$ MSE_{e } $都会随着时间的推移而逐渐减小, 说明变阈值方法通过不断获取在线数据对系统未来的退化轨迹预期做出 调整, 有利于提高预测精度, 这是基于固定阈值方法所不具备的优势. 同时, 表 3中 $ MSE_{e} $的大小总体上也小于表 2, 证明本文所提方法能够更加准确地描述系统的健康水平.至此, 本文提出的剩余寿命预测方法的有效性和优越性得到了验证.

4 总结与展望

本文针对考虑执行器性能退化的闭环控制系统, 提出一种基于解析模型的控制系统剩余寿命预测方法.首先, 采用权值优选粒子滤波算法, 利用在线输入输出数据估计出执行器的隐含退化量; 然后, 利用MC仿真设置合理的失效阈值, 建立基于变失效阈值的系统失效判断准则; 最后, 将隐含退化量的估计值代入退化模型中外推出剩余寿命分布.仿真实验表明, 本文提出的方法能够解决闭环控制系统的可靠性评估问题, 为系统失效前的健康管理活动提供可信依据, 提高了系统的安全性和可靠性.相比于传统的单部件系统预测方法, 本文方法考虑了控制作用的影响, 具有潜在的工程应用价值.

下一步可能的研究方向包括: 1)实际工程中的控制系统往往是由多执行器、多传感器、多控制器共同组成的, 这些部件都有可能发生性能退化而导致系统失效, 可以针对多部件隐含退化过程展开研究; 2)可以根据剩余寿命预测结果进行相应的维修、管理, 或是改进控制器的设计以延长系统寿命, 使该方法进一步应用到实际工程中.