Prediction Method of Hot Metal Silicon Content in Blast Furnace Based on Optimal Smelting Condition Migration
-
摘要: 高炉铁水硅含量是铁水品质与炉况的重要表征, 冶炼过程关键参数频繁波动及大时滞特性给高炉铁水硅含量预测带来了巨大挑战. 提出一种基于最优工况迁移的高炉铁水硅含量预测方法. 首先, 针对过程变量频繁波动问题, 提出基于邦费罗尼指数的自适应密度峰值聚类算法, 实现对高炉冶炼过程变量的工况划分, 并建立不同工况硅含量预测子模型. 其次, 针对冶炼过程的大时滞特性, 定义相邻时间节点间的硅含量工况迁移代价函数, 并提出多源路径寻优算法, 实现冶炼过程中硅含量最优工况迁移路径及当前时刻硅含量最优预测值的求解. 最后, 基于工业现场数据验证了所提方法的有效性与准确性.Abstract: The hot metal silicon content in blast furnace can characterize the hot metal quality and the condition of blast furnace. It poses a great challenge to the online prediction of silicon content because of the frequent fluctuation of smelting parameters and the existence of large time delay during the ironmaking process. This paper proposes an algorithm for predicting the hot metal silicon content in blast furnace based on optimal smelting condition migration. Firstly, arming at the frequent fluctuation of smelting process variables, an adaptive density peak clustering algorithm based on the Bonferroni index to dynamically cluster the process variables of blast furnace ironmaking process is proposed, which can obtain clusters of different smelting conditions, and establish sub-models for different smelting conditions. Secondly, to mitigate the large time delay of blast furnace ironmaking process, this paper defines the silicon content migration cost function between adjacent time nodes, and proposes a multi-source path optimization algorithm to solve the optimal migration path of silicon content during the smelting process and the optimal prediction value of silicon content at the current time. Finally, the effectiveness and accuracy of the proposed method are verified based on the industrial field data.
-
随着高新技术的迅猛发展, 现代工业设备正朝着大型化、复杂化和智能化趋势快速发展. 这类设备在运行过程中由于受到内部和外部因素的随机影响, 性能和健康状态不可避免地呈现下降趋势乃至退化失效, 导致无法完成正常任务和功能, 进而引发严重事故, 造成环境破坏和人员伤亡[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} $$ -
图 2 基于最优工况迁移的建模策略
Fig. 2 Modeling strategy based on optimal smelting condition migration
图 5 过程变量工况隶属度与模型预测
Fig. 5 Process variable membership degree matching and model prediction
图 6 铁水硅含量工况迁移图
Fig. 6 Smelting condition migration diagram of hot metal silicon content
图 13 最优工况迁移模型硅含量预测结果
Fig. 13 Prediction of silicon content in hot metal based on optimal smelting condition migration model
表 1 过程变量MIC相关性系数
Table 1 MIC correlation coefficient of process variables
过程变量 MIC 系数 过程变量 MIC 系数 富氧率 0.291 总压差 0.204 透气性指数 0.270 炉腹煤气指数 0.278 标准风速 0.275 热风压力 0.268 富氧流量 0.218 实际风速 0.173 冷风流量 0.264 冷风温度 0.209 鼓风动能 0.204 热风温度 0.213 设定喷煤量 0.241 顶温下降管 0.209 理论燃烧温度 0.248 铁水红外温度 0.291 顶压 0.195 顶温 0.292 富氧压力 0.229 鼓风湿度 0.179 冷风压力 0.197 阻力系数 0.204 
下载: 导出CSV
表 2 聚类中心截断标志
Table 2 Cluster center truncation flag
序号 1 2 3 4 5 6 截断系数 3.00 4.02 42.30 52.50 28.02 24.34
下载: 导出CSV
表 3 寻优算法耗时对比
Table 3 Comparison of the time consumption of optimization algorithms
寻优算法 节点数 40 80 120 160 200 Floyd 算法
耗时 (ms)3.20 × 104 2.72 × 105 8.96 × 105 2.09 × 106 4.05 × 106 本文算法
耗时(ms)3 8 11 13 18
下载: 导出CSV
表 4 模型性能对比
Table 4 Model performance comparison
模型类别 性能指标 数值预测
命中率 (%)趋势预测
准确率 (%)预测均方误差 工况迁移预测模型 88 82 0.0043 Elman 网络 79 69 0.0069 Elman-Adaboost 85 71 0.0054 FEEMD-Adaboost-Elman 86 74 0.0049
下载: 导出CSV
-
[1] Biswas A K. Principles of Blast Furnace Ironmaking —— Theory and Practice. Brisbane: Cootha Publishing House, 1981. 1−12 [2] Zhou Ping, Zhang Shuai, Dai Peng. Recursive Learning Based Bilinear Subspace Identification for Online Modeling and Predictive Control of a Complicated Industrial Process. IEEE Access, 2020, 8: 6253-62541. [3] Saxén H, Gao Chuan-Hou, Gao Zhi-Wei. Data-driven time discrete models for dynamic prediction of the hot metal silicon content in the blast furnace—A review. IEEE Transactions on Industrial Informatics, 2012, 9(4): 2213-2225. [4] Onorin O P, Spirin N A. Features of blast furnace transient processes. Metallurgist, 2017, 61(1): 121-126. [5] Spirin N A, Onorin O P, Istomin A S, Lavrov V V, Gurin I A. Study of transition processes of blast-furnace smelting by the mathematical model method. In: Proceedings of the 2018 IOP Conference Series, Materials Science and Engineering. Novokuznetsk, Russia: Institute of Physics Publishing, 2018. 012−073 [6] Spirin N, Onorin O, Alexander I. Prediction of blast furnace thermal state in real-time operation. Solid State Phenomena, 2020, 299: 518-523. doi: 10.4028/www.scientific.net/SSP.299.518 [7] Spirin N. A, Polinov A A, Gurin I A, Pishnograev SN. Information System for Real-Time Prediction of the Silicon Content of Iron in a Blast Furnace. Metallurgist, 2020, 63(9): 898-905. [8] Östermark R, Saxen H. VARMAX-modelling of blast furnace process variables. European Journal of Operational Research, 1996, 90(1): 85-101. doi: 10.1016/0377-2217(94)00304-1 [9] Saxen H, Östermark R. State realization with exogenous variables-A test on blast furnace data. European journal of operational research, 1996, 89(1): 34-52. doi: 10.1016/S0377-2217(96)90050-8 [10] Bhattacharya T. Prediction of silicon content in blast furnace hot metal using partial least squares. ISIJ international, 2005, 45(12): 1943-1945. doi: 10.2355/isijinternational.45.1943 [11] Li Jun-Peng, Hua Chang-Chun, Yang Yan-Na, Guan Xin-Ping. Bayesian block structure sparse based T–S fuzzy modeling for dynamic prediction of hot metal silicon content in the blast furnace. IEEE Transactions on Industrial Electronics, 2017, 65(6): 4933-4942. [12] Xu Xia, Hua Chang-Chun, Tang Ying-Gan, Guan Xing-Ping. Modeling of the hot metal silicon content in blast furnace using support vector machine optimized by an improved particle swarm optimizer. Neural Computing and Applications, 2016, 27(6): 1451-1461. doi: 10.1007/s00521-015-1951-7 [13] Han Y, Li J, Yang X L, Liu W X, Zhang Y Z. Dynamic prediction research of silicon content in hot metal driven by big data in blast furnace smelting process under hadoop cloud platform. Complexity, DOI: 10.1155/2018/8079697 [14] Xu X, Hua C C, Tang Y G, Guan X P. Wiener model identification of blast furnace ironmaking process based on laguerre filter and linear programming support vector regression. In: Proceedings of the 2014 International Joint Conference on Neural Networks. Beijing, China: IEEE Press, 2014. 2198−2204 [15] Zhou Ping, Li Wen-Peng, Wang Hong. Robust Online Sequential RVFLNs for Data Modeling of Dynamic Time-Varying Systems With Application of an Ironmaking Blast Furnace. IEEE Transactions on Cybernetics, 2020, 50(11): 4783-4795. doi: 10.1109/TCYB.2019.2920483 [16] 郜传厚, 渐令, 陈积明, 孙优贤. 复杂高炉炼铁过程的数据驱动建模及预测算法. 自动化学报, 2009, 35(06): 725-730. doi: 10.3724/SP.J.1004.2009.00725Gao Chuan-Hou, Jian Ling, Chen Ji-Ming, Sun You-Xian. Data-driven modeling and prediction algorithm for complex blast furnace ironmaking process. Acta Automatica Sinica, 2009, 35(06): 725-730. doi: 10.3724/SP.J.1004.2009.00725 [17] David S F, David F F, Machado M L P. Artificial neural network model for predict of silicon content in hot metal blast furnace. Materials Science Forum, 2016, 869: 572-577. doi: 10.4028/www.scientific.net/MSF.869.572 [18] 宋菁华, 杨春节, 周哲. 改进型 EMD-Elman 神经网络在铁水硅含量预测中的应用. 化工学报, 2016, 67(3): 729-735.Song Jing-Hua, Yang Chun-Jie, Zhou Zhe. Application of improved EMD-Elman neural network in prediction of silicon content in molten iron. CIESC Journal, 2016, 67(3): 729-735. [19] Jiang Ke, Jiang Zhao-Hui, Xie Yon-Fang, Chen Zhi-Peng. Classification of silicon content variation trend based on fusion of multilevel features in blast furnace ironmaking. Information Sciences, 2020, 521: 32-45. doi: 10.1016/j.ins.2020.02.039 [20] 周平, 张丽, 李温鹏, 戴鹏, 柴天佑. 集成自编码与PCA的高炉多元铁水质量随机权神经网络建模. 自动化学报, 2018, 44(10): 1799-1811.Zhou Ping, Zhang Li, Li Wen-Peng, Dai Peng, Chai Tian-You. Modeling of blast furnace multi-element molten iron quality with random weight neural network based on self-encoding and PCA. Acta Automatica Sinica, 2018, 44(10): 1799-1811. [21] Jian Ling, Song Yun-Quan, Shen Shu-Qian, Wang Yan, Yin Hai-Qing. Adaptive least squares support vector machine predictor for blast furnace ironmaking process. ISIJ International, 2015, 55(4): 845-850. doi: 10.2355/isijinternational.55.845 [22] Zeng Jiu-Sun, Liu Xiang-Guan, Gao Chuan-Hou, Luo Shi-Hua, Jian Ling. Wiener model identification of blast furnace ironmaking process. ISIJ International, 2008, 48(12): 1734-1738. doi: 10.2355/isijinternational.48.1734 [23] 蒋朝辉, 董梦林, 桂卫华, 阳春华, 谢永芳. 基于Bootstrap的高炉铁水硅含量二维预报. 自动化学报, 2016, 42(05): 715-723.Jiang Zhao-Hui, Dong Meng-Lin, Gui Wei-Hua, Yang Chun-Hua, Xie Yon-Fang. Two-dimensional prediction for silicon content of hot metal of blast furnace based on bootstrap. Acta Automatica Sinica, 2016, 42 (5): 715-723. [24] 李温鹏, 周平. 高炉铁水质量鲁棒正则化随机权神经网络建模. 自动化学报, 2020, 46(04): 721-733.Li Wen-Peng, Zhou Ping. Wen Liang. Blast furnace hot metal quality robust regularization random weight neural network modeling. Acta Automatica Sinica, 2020, 46(04): 721-733. [25] 温亮, 周平. 基于多参数灵敏度分析与遗传优化的铁水质量无模型自适应控制. 自动化学报.Wen Liang, Zhou Ping. Model-free adaptive control of molten iron quality based on multi-parameter sensitivity analysis and genetic optimization. Acta Automatica Sinica, to be published. [26] Rodriguez, A, Laio A. Clustering by fast search and find of density peaks. Science, 2014, 344(6191): 1492-1496. doi: 10.1126/science.1242072 [27] Martin B, Elena, Silber J. The Bonferroni index and the measurement of distributional change. Metron, 2017, 75(1): 1-16. doi: 10.1007/s40300-016-0105-8 [28] 孙甜, 凌卫新. 基于模拟退火的 Levenberg-Marquardt 算法在神经网络中的应用. 科学技术与工程, 2008(18): 5189−5192Sun Tian, Ling Wei-Xin. Application of Levenberg-Marquardt algorithm based on simulated annealing in neural network. Science Technology and Engineering, 2008(18): 5189−5192 [29] Reshef, D N, Reshef Y A. Detecting novel associations in large data sets. Science, 2011, 334(6062): 1518-1524. doi: 10.1126/science.1205438 [30] Pan Dong, Jiang Zhao-Hui, Chen Zhi-Peng. Temperature measurement and compensation method of blast furnace molten iron based on infrared computer vision. IEEE Transactions on Instrumentation and Measurement, 2018, 68(10): 3576-3588. [31] 庄田, 杨春节. 基于Elman-Adaboost强预测器的铁水硅含量预测方法. 冶金自动化, 2017, 41(04): 1-6+17.Zhuang Tian, Yang Chun-Jie. Prediction method of silicon content in molten iron based on Elman-Adaboost strong predictor. Metallurgical Automation, 2017, 41(04): 1-6+17. [32] 王凯, 毕贵红, 高晗, 蒲娴怡, 陈仕龙. 基于改进快速集合经验模态分解和Elman-Adaboost的短期风速预测方法. 电力科学与工程, 2020, 36(05): 32-39. doi: 10.3969/j.ISSN.1672-0792.2020.05.005Wang Kai, Bi Gui-Gong, Gao Han, Pu Xian-Yi, Chen Shi-Long. Short-term Wind Speed Prediction Method Based on Improved Fast Ensemble Empirical Mode Decomposition and Elman-Adaboost. Electric Power Science and Engineering, 2020, 36(05): 32-39. doi: 10.3969/j.ISSN.1672-0792.2020.05.005 [33] 蒋珂, 蒋朝辉, 谢永芳, 潘冬, 桂卫华. 大型高炉铁水硅含量变化趋势的智能预报. 控制工程, 2020, 27(03): 540-546.Jiang Ke, Jiang Zhao-Hui, Xie Yon-Fang, Pan Dong, Gui Wei-Hua. Intelligent prediction of silicon content change trend in molten iron of large blast furnace. Control Engineering, 2020, 27 (03): 540-546. 期刊类型引用(14)
1. 潘冬，许川，龚芃旭，蒋朝辉，桂卫华. 基于红外与可见光视觉的高炉铁口铁水温度场在线检测. 自动化学报. 2025(02): 343-355 . 
本站查看2. 方怡静，蒋朝辉，黄良，桂卫华，潘冬. 基于知识与AW-ESN融合的烧结过程FeO含量预测. 自动化学报. 2024(02): 282-294 . 本站查看3. 王帅，张朝晖，邢相栋，惠佳豪，折媛. 基于大数据的IWOA-KELM铁水硅含量预测模型. 钢铁研究学报. 2024(02): 188-198 . 百度学术4. 刘然，高子扬，李宏扬，刘小杰，吕庆，陈树军. 高炉智能优化系统综述与展望. 中国冶金. 2024(03): 13-24+36 . 百度学术5. 贺锦峰，张屹龙，曹益春，谭涛，陈斌. 基于Attention-LSTM的低硅冶炼硅含量预测研究. 机电工程技术. 2024(09): 56-59 . 百度学术6. 邱国兴，蔡明冲，张毅，苏炳瑞，杨永坤，李小明. 基于灰色关联分析和机器学习的高炉铁水硅含量预测. 材料导报. 2024(20): 256-261 . 百度学术7. 蒋珂，蒋朝辉，谢永芳，潘冬，桂卫华. 基于时序关联矩阵的高炉冶炼过程多重关联时延估计方法. 自动化学报. 2023(02): 329-342 . 本站查看8. 蒋珂，蒋朝辉，谢永芳，潘冬，桂卫华. 基于动态注意力深度迁移网络的高炉铁水硅含量在线预测方法. 自动化学报. 2023(05): 949-963 . 本站查看9. 刘然，赵伟光，刘颂，刘小杰，李宏扬，吕庆. 高炉冶炼智能化的发展与探讨. 钢铁. 2023(05): 1-10 . 百度学术10. 黄建才，蒋朝辉，桂卫华，潘冬，许川，周科. 基于状态识别的高炉料面视频关键帧提取方法. 自动化学报. 2023(11): 2257-2271 . 本站查看11. 季玉璋，张卫军，刘馨，池中源. 大数据技术推动高炉炼铁智能化发展的研究现状及展望. 冶金经济与管理. 2023(06): 17-20 . 百度学术12. 黄山文. 铁水硅含量对铁钢系统生产平衡的影响分析. 冶金与材料. 2023(12): 163-165+168 . 百度学术13. 蒋珂，蒋朝辉，谢永芳，潘冬，桂卫华. 高炉铁水质量信息在线检测方法综述. 冶金自动化. 2022(02): 19-33+45 . 百度学术14. 廖亚楠，王业林，李萌，肖清泰，王华. 基于CEEMDAN-RVM-EC的还原冶炼温度预报. 控制理论与应用. 2022(11): 2177-2184 . 百度学术其他类型引用(14)
-
下载:
计量
- 文章访问数: 1513
- HTML全文浏览量: 337
- PDF下载量: 298
- 被引次数: 28
