CN112487694B - 一种基于多退化指标的复杂设备剩余寿命预测方法 - Google Patents

一种基于多退化指标的复杂设备剩余寿命预测方法 Download PDF

Info

Publication number
CN112487694B
CN112487694B CN202011366448.7A CN202011366448A CN112487694B CN 112487694 B CN112487694 B CN 112487694B CN 202011366448 A CN202011366448 A CN 202011366448A CN 112487694 B CN112487694 B CN 112487694B
Authority
CN
China
Prior art keywords
degradation
equation
parameters
indexes
index
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011366448.7A
Other languages
English (en)
Other versions
CN112487694A (zh
Inventor
陈绍炜
王美男
赵帅
温鹏飞
汪盛悦
窦智
李毅
薛峰
李勇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202011366448.7A priority Critical patent/CN112487694B/zh
Publication of CN112487694A publication Critical patent/CN112487694A/zh
Application granted granted Critical
Publication of CN112487694B publication Critical patent/CN112487694B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Economics (AREA)
  • Physics & Mathematics (AREA)
  • Strategic Management (AREA)
  • Theoretical Computer Science (AREA)
  • Entrepreneurship & Innovation (AREA)
  • General Physics & Mathematics (AREA)
  • Development Economics (AREA)
  • Marketing (AREA)
  • Game Theory and Decision Science (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Educational Administration (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Testing And Monitoring For Control Systems (AREA)

Abstract

本发明提供了一种基于多退化指标的复杂设备剩余寿命预测方法，对多个退化指标建立多元退化模型，将该多元退化模型中的模型参数划分外部参数和退化轨迹参数，并且在该多元退化模型中使用协方差矩阵建立多个退化指标的退化轨迹参数之间的相关关系，使用最大似然估计对外部参数进行离线估计，且使用基于核平滑的粒子滤波方法对模型中的退化轨迹参数与多个指标的退化状态进行在线估计，从而实现基于多退化指标的设备剩余寿命预测。本发明实现了对多退化指标的合理建模，从有噪声的量测数据中对多退化指标的真实退化状态以及退化轨迹参数进行在线估计，基于最后一个时刻的多退化指标状态估计的概率分布，预测得到设备的剩余使用寿命。

Description

一种基于多退化指标的复杂设备剩余寿命预测方法
技术领域
本发明涉及系统可靠性领域，更具体的涉及复杂设备的剩余寿命预测方法。
背景技术
由于设备不同部件的多个退化过程包含了相同的退化信息，这多个退化过程往往具有一定的依赖关系。因此对于具有多个部件的大型复杂系统，如何同时对其多个部件的退化情况进行监测和状态估计，并且合理的描述多个退化过程之间的依赖关系，从而实现基于多退化指标的设备剩余寿命(remaining useful life,RUL)预测具有重要的工程意义。
通过对目前的文献检索发现，现有技术大多仅对单一的退化指标进行分析，仅有少部分文献同时分析多个退化指标。大部分基于多退化指标的预测方法将多个退化指标之间的依赖关系建立在多个退化过程的内部随机性当中，如Xi等人在《Remaining usefullife prediction for multicomponent systems with hidden dependencies》中为每个退化指标建立维纳过程的退化模型，并且使用一个扩散系数矩阵(diffusion coefficientmatrix)将多个退化指标之间的依赖关系建立在维纳过程的随机项中；Peng等人所著的《Joint Online RUL Prediction for Multivariate Deteriorating Systems》一文，以及Fang等人所著的《Copula-based reliability analysis of degrading systems withdependent failures》一文中，均为每个退化指标建立随机过程的退化模型，并且使用copula函数将多个指标之间的依赖关系建立在随机过程模型的随机项中。
由于不同设备的多个退化过程之间可能具有不同的依赖结构，例如，某些设备的不同部件的退化速度之间可能存在一定相关性，这种相关性可能是导致不同退化指标依赖关系的原因。当为一个退化指标建立对应的随机过程模型时，模型中的退化轨迹参数在一定程度上反映了该指标的退化速率，因此针对这种依赖结构，将依赖关系体现在随机过程模型的参数相关性上更为合理。然而当前的基于多指标的预测方法均将多指标之间的依赖关系建立在随机过程模型的随机项中，对于这种依赖关系体现为不同指标退化速率相关性的多部件复杂设备，尚未出现一种能够合理解释依赖关系的多退化指标设备剩余寿命预测方法，使得对于这类设备的健康状态估计与剩余寿命预测不能满足工程应用的要求。
发明内容
为了克服现有技术的不足，本发明提供一种基于多退化指标的复杂设备剩余寿命预测方法。对多个退化指标建立多元退化模型，将该多元退化模型中的模型参数根据其物理意义分划分外部参数和退化轨迹参数，并且在该多元退化模型中使用协方差矩阵建立多个退化指标的退化轨迹参数之间的相关关系。使用最大似然估计对外部参数进行离线估计，且使用基于核平滑的粒子滤波(kernel smoothing–particle filter,KS-PF)方法对模型中的退化轨迹参数与多个指标的退化状态进行在线估计，从而实现基于多退化指标的设备剩余寿命预测。此外，本发明还采用最优调整(optimal tuning,OT)的方法对KS-PF中的核参数进行在线的自适应调整。
本发明解决其技术问题所采用的技术方案的详细步骤如下：
步骤1：边缘退化模型的建立；
假设设备有K个退化指标，为每个退化指标建立维纳过程的边缘退化模型：
Xk(t)＝μk(θk,t)+σkB(t) (1)
Figure BDA0002805559250000023
公式(1)为状态转移方程，公式(2)为量测方程，其中，k＝1,2…K代表第k个退化指标，Xk(t)是第k个退化指标在t时刻的退化状态，μk(θk,t)是一个单调可微函数，被称为退化轨迹，为线性或非线性函数；θk＝[ak,bk,…]T为退化轨迹参数，决定指标的退化速率，个体不确定性指的是同一批次的不同设备之间的退化多样性，通过对退化指标的状态值进行回归得到，B(t)是一个标准布朗运动过程，σk是波动系数；σkB(t)具有正态增量，即有
Figure BDA0002805559250000021
其中0≤t2≤t1，描述了第k个指标在退化过程中的内部不确定性，即指标在退化过程中本身的退化随机性；Yk(t)代表第k个退化指标的量测值，其中
Figure BDA0002805559250000022
是服从高斯分布的随机变量εk的方差，εk代表量测噪声随机变量，体现了量测不确定性，服从一个零均值的正态分布，对退化指标建立随机过程模型，则退化模型中的每个部分均代表了一定的物理含义；
令xk,j＝Xk(tj)，yk,j＝Yk(tj)，Δt＝tj-tj-1，其中0≤j≤M，M为待预测设备中退化指标量测点的个数，仅考虑第k个退化指标，则边缘退化模型的公式(1)和公式(2)被转化为粒子滤波方法需要的状态转移方程和量测方程，其中，状态转移方程为f(xk,j|xk,j-1)：
Figure BDA0002805559250000031
量测方程h(yk,j|xk,j)为：
Figure BDA0002805559250000032
当公式(3)和(4)被用于多个退化指标时，令xj＝[x1,j,…,xK,j]T表示所有退化指标在tj时刻的状态值，令yj＝[y1,j,…,yK,j]T代表所有退化指标在tj时刻的量测值，则状态转移方程和量测方程转化为f(xj|xj-1)和h(yj|xj)；此时的多元退化模型f(xj|xj-1)和h(yj|xj)未建立多个退化指标之间的依赖关系；
在多个退化指标同时退化的情况中，假设设备的失效时刻是先到达对应失效阈值的某退化指标到达该退化指标的失效阈值的时刻，因此对一个具有K个退化指标的复杂设备，每个退化指标的失效阈值定义为Hk,k＝1,2,…,K，该设备在tj时刻开始预测时的剩余使用寿命为：
Figure BDA0002805559250000033
步骤2：外部参数离线估计；
外部参数
Figure BDA0002805559250000034
和
Figure BDA0002805559250000035
使用最大似然估计方法利用训练集数据进行离线估计；为了区分训练集变量以及测试集变量，训练集的变量用上标'*'表示，令
Figure BDA0002805559250000036
Figure BDA0002805559250000037
及
Figure BDA0002805559250000038
分别表示训练集中的设备的第k个退化指标的退化状态、量测值以及退化轨迹，则
Figure BDA0002805559250000039
和
Figure BDA00028055592500000310
的分布由标准布朗运动过程的性质推得：
Figure BDA00028055592500000311
Figure BDA00028055592500000312
其中，Q＝[min{ti,tj}]1≤i,j≤M，IM是一个M阶单位矩阵；
外部参数的似然函数如下：
Figure BDA00028055592500000313
Figure BDA00028055592500000314
为了估计外部参数，首先将退化状态
Figure BDA00028055592500000315
识别出来，然后通过公式(9)使用最大似然估计得到参数
Figure BDA0002805559250000041
选择一个退化趋势与实际量测值的趋势一致的函数μk(·)描述第k个退化指标的退化轨迹，然后使用μk(·)对训练集的退化状态
Figure BDA0002805559250000042
进行最小二乘拟合，估计得到训练集每个设备的退化轨迹参数
Figure BDA0002805559250000043
参数
Figure BDA0002805559250000044
通过使公式(8)的似然函数最大来估计得到；
步骤3：测试集设备退化轨迹参数与退化指标状态的同时在线估计；
步骤3.1：基于参数相关性的多元退化模型建立；
假设模型(1)中的退化轨迹参数θ1,θ2,…,θK之间线性相关，则不同退化指标X1(t),X2(t),…,XK(t)之间的复杂依赖关系由线性相关的退化轨迹参数以及θk和Xk(t)之间的线性或非线性的映射共同构建和描述；首先估计θ1,θ2,…,θK之间的皮尔逊相关系数，假设退化轨迹被设置为μk(θk,t)＝akexp(bk·t)，θk＝[ak,bk]T，则
Figure BDA0002805559250000045
和
Figure BDA0002805559250000046
之间的皮尔逊相关系数
Figure BDA0002805559250000047
以及
Figure BDA0002805559250000048
和
Figure BDA0002805559250000049
之间的皮尔逊相关系数
Figure BDA00028055592500000410
由训练集数据
Figure BDA00028055592500000411
估计得到，而
Figure BDA00028055592500000412
由最小二乘估计得到，其中1≤k1,k2≤K；
使用基于核平滑的粒子滤波方法对测试集设备的多个退化指标及其对应的退化轨迹参数进行在线估计，首先建立基于参数间相关性的状态转移方程，通过在状态转移方程中构建一个协方差矩阵来实现建立参数之间相关性的目的，从一个多元高斯分布中进行采样，得到多个退化指标的退化轨迹参数的粒子群，而该多元高斯分布的协方差矩阵体现了多个退化指标参数之间的相关关系；
步骤3.2：KS-PF实现退化轨迹参数与多退化指标的同时在线估计
使用粒子滤波中的贯序重要性重采样(sequential importance resampling,SIR)框架，按照状态转移方程(10)和量测方程(4)，对测试集设备的多指标退化状态x＝[X1(t),…,XK(t)]T和退化轨迹参数θk,j＝[ak,j,bk,j,…]T,k＝1,…,K进行在线估计；
步骤3.3：最优调整OT方法实现核参数s的自适应在线调整；
采用最优调整(Optimized Tuning,OT)方法，在每个状态估计时刻tj，均选择一个最优的核参数sj；
在tj时刻，设备多个退化指标的多元退化状态的先验分布Qj＝p(xj|y1:j-1)表示为公式(11)，后验分布Pj＝p(xj|y1:j)表示为公式(15)，先验分布与后验分布的KL散度计算为：
Figure BDA00028055592500000413
贝叶斯滤波中的更新公式为：
Figure BDA0002805559250000051
将公式(18)和公式(11)带入公式(17)，并且使用狄拉克函数δ(·)的挑选性质，得到：
Figure BDA0002805559250000052
将公式(13)带入公式(14)，得到：
Figure BDA0002805559250000053
最后，将公式(20)带入公式(19)，得到：
Figure BDA0002805559250000054
因此，在每个tj时刻，在SIR框架中用到的状态转移方程(10)里的核参数s均使用最优的核参数：
Figure BDA0002805559250000055
用状态转移方程(10)和量测方程(4)，按照粒子滤波的SIR框架，即在每个tj时刻均使用公式(13)、(14)得到更新后的粒子集
Figure BDA0002805559250000056
得到的粒子
Figure BDA0002805559250000057
代表了K个退化指标的状态估计粒子值，得到的
Figure BDA0002805559250000058
和
Figure BDA0002805559250000059
对于任意k1,k2∈[1,K]都不是独立的，具有一定的依赖关系，多个退化指标的多元联合分布p(x1,j,…,xK,j|y1:j)即可得到；
步骤4：多退化指标的复杂设备剩余寿命预测；
在每个tj时刻均得到粒子集
Figure BDA00028055592500000510
则多个退化指标的多元联合分布可被估计为：
Figure BDA00028055592500000511
在tj时刻开始预测时，将估计得到的每个粒子代表的参数和系统状态均代入模型(3)中，则通过设备失效定义得到基于每个粒子预测的剩余使用寿命；基于第i个粒子预测得到的设备剩余使用寿命为：
Figure BDA0002805559250000061
在tj时刻开始预测时，预测得到的设备的剩余使用寿命的概率密度函数(probability density function,PDF)为：
Figure BDA0002805559250000062
即求得剩余使用寿命的概率密度函数，即式(25)中的p(lj|y1:j)。
所述步骤3.1的基于参数相关性的多元退化模型建立的具体步骤如下：
令Δuj＝[Δμ1(tj),…,ΔμK(tj)]T，ωj＝[ω1,j,…,ωK,j]T，并且令
Figure BDA0002805559250000063
Figure BDA0002805559250000064
分别为参数a＝[a1,…,aK]T，参数b＝[b1,…,bK]T和多个退化指标x＝[X1(t),…,XK(t)]T在tj时刻的第i个粒子。令θj＝[θ1,j,…,θK,j]T，θk,j＝[ak,j,bk,j]T，则基于参数相关性的多个退化指标的状态转移方程f(xj,θj|xj-1,θj-1)如下：
Figure BDA0002805559250000065
其中，
Figure BDA0002805559250000066
Figure BDA0002805559250000067
是tj时刻x，a和b的粒子与其对应的粒子权值，
Figure BDA0002805559250000068
和
Figure BDA0002805559250000069
分别是
Figure BDA00028055592500000610
的加权均值和加权方差，
Figure BDA00028055592500000611
和
Figure BDA00028055592500000612
分别是
Figure BDA00028055592500000613
的加权均值和加权方差；
当粒子
Figure BDA00028055592500000614
从粒子滤波中的建议分布(proposal distribution)中采样时，建议分布为两个多元高斯分布，并且表征参数之间相关关系的
Figure BDA00028055592500000615
和
Figure BDA00028055592500000616
均包含在该多元高斯分布的协方差矩阵Va,j-1和Vb,j-1中，从建议分布中采样得到的参数粒子之间是具有相关关系的。
所述步骤3.2中KS-PF实现退化轨迹参数与多退化指标的同时在线估计的具体步骤如下：
1)初始化：
令j＝0，设置t0时刻的设备的多指标退化状态和退化轨迹参数的初始分布p(x0)和p(θ0)；
2)状态估计：
令j加1，若存在量测向量yj，顺序执行以下步骤，否则终止该程序；
a)预测：从建议分布p(xj,θj|xj-1,θj-1)中按照状态转移方程(10)采样得到粒子集
Figure BDA0002805559250000071
xj和θj的先验分布，用该粒子集表示为：
Figure BDA0002805559250000072
Figure BDA0002805559250000073
b)更新：使用tj时刻新到来的退化指标量测向量yj按照量测方程(4)对粒子权值进行更新：
Figure BDA0002805559250000074
并且该权值被归一化为：
Figure BDA0002805559250000075
c)重采样：有效粒子数被计算为
Figure BDA0002805559250000076
若Nneff小于一个给定的阈值，则执行重采样步骤，否则不执行，得到的粒子集为
Figure BDA0002805559250000077
xj和θj的后验分布用该粒子集表示为：
Figure BDA0002805559250000078
Figure BDA0002805559250000079
d)迭代执行步骤a)-步骤c)，即得到状态估计。
所述步骤2中，退化状态
Figure BDA00028055592500000710
识别采用简单滑动平均(simple moving average,SMA)或加权滑动平均(weighted moving average,WMA)。
所述
Figure BDA00028055592500000711
s是处于0到1之间的核参数。
本发明的有益效果在于通过为多个退化指标建立随机过程的多元退化模型，并将多个退化指标之间的依赖关系建立在模型退化轨迹参数的相关性上，实现了对多退化指标的合理建模；使用最大似然方法实现了模型中外部参数的离线估计；构建了一种最优调整核参数的基于核平滑的粒子滤波(OTKS-PF)方法，从有噪声的量测数据中对多退化指标的真实退化状态以及退化轨迹参数进行在线估计；最后基于最后一个时刻的多退化指标状态估计的概率分布，预测得到设备的剩余使用寿命。
附图说明
图1是本发明的实现框架图。
图2是第一组实验中训练集第1号发动机的两个退化指标的退化状态及退化状态量测值。其中图2(a)为训练集第1号发动机的量测值及退化指标1，图2(b)为训练集第1号发动机的量测值及退化指标2。
图3是第一组实验中所有训练集发动机的退化轨迹参数箱线图。
图4是第一组实验中训练集第1号和8号发动机的退化状态及模型拟合曲线。
图5是第一组实验中所有训练集发动机的退化轨迹参数散点图。
图6是第一组实验中测试集第5号发动机的退化轨迹参数的在线估计结果。其中图6(a)为采用退化指标1的估计结果，图6(b)为采用退化指标2的估计结果。
图7是第一组实验中测试集第19号发动机分别使用单指标和多指标联合预测得到的RUL预测结果。其中图7(a)显示了RUL点预测结果及概率密度函数预测结果；图7(b)显示了RUL点预测结果及置信区间预测结果。
图8是第一组实验中测试集第68号发动机分别使用单指标和多指标联合预测得到的RUL预测结果。其中图8(a)显示了RUL点预测结果及概率密度函数预测结果；图8(b)显示了RUL点预测结果及置信区间预测结果。
图9是第一组实验所有20个测试集发动机的RUL平均预测结果。其中图9(a)显示了预测的平均误差结果，图9(b)显示了预测的平均95％置信区间宽度结果。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
本发明的实现框架如图1所示，由离线模型训练和模型的在线估计两部分构成。在离线模型训练部分，使用一组训练集设备对退化模型进行参数估计，以得到能够反映这组设备的共同特征的模型参数(模型的外部参数)。而在线估计部分，则针对特定的测试集设备，对其进行个性化。具体的，利用该设备的量测信息去在线估计其特有的退化信息(退化轨迹参数和退化指标的状态值)，从而实现对该测试集设备的RUL预测。
本发明的具体步骤如下：
步骤1：边缘退化模型的建立
假设设备有K个退化指标，为每个退化指标建立维纳过程的边缘退化模型：
Xk(t)＝μk(θk,t)+σkB(t) (1)
Figure BDA0002805559250000091
公式(1)为状态转移方程，公式(2)为量测方程，其中，k＝1,2…K代表第k个退化指标，Xk(t)是第k个退化指标在t时刻的退化状态，μk(θk,t)是一个单调可微函数，被称为退化轨迹，为线性或非线性函数；θk＝[ak,bk,…]T为退化轨迹参数，决定指标的退化速率，μk(θk,t)体现了退化过程的个体不确定性以及非线性不确定性；个体不确定性指的是同一批次的不同设备之间的退化多样性，通过对退化指标的状态值进行回归得到，B(t)是一个标准布朗运动过程，σk是波动系数；σkB(t)具有正态增量，即有
Figure BDA0002805559250000092
其中0≤t2≤t1，描述了第k个指标在退化过程中的内部不确定性，即指标在退化过程中本身的退化随机性；Yk(t)代表第k个退化指标的量测值，其中
Figure BDA0002805559250000093
是服从高斯分布的随机变量εk的方差，εk代表量测噪声随机变量，体现了量测不确定性，它服从一个零均值的正态分布，对退化指标建立随机过程模型，则退化模型中的每个部分均代表了一定的物理含义；
假设
Figure BDA0002805559250000094
和
Figure BDA0002805559250000095
仅与工况和量测条件有关，对于处于单工况相同量测条件的一批不断发生退化的设备而言，参数
Figure BDA0002805559250000096
是相同的，而参数θk对于不同的退化设备是不同的。因此，区别于退化轨迹参数θk，
Figure BDA0002805559250000097
和
Figure BDA0002805559250000098
被定义为外部参数，并且两种参数用在线和离线的方法分别估计得到，
Figure BDA0002805559250000099
和
Figure BDA00028055592500000910
的估计方法是最大似然估计，似然函数为公式(8)和公式(9)，具体方法参见步骤2；θk的在线估计方法参见步骤3)。
令xk,j＝Xk(tj)，yk,j＝Yk(tj)，Δt＝tj-tj-1，其中0≤j≤M，M为待预测设备中退化指标量测点的个数，仅考虑第k个退化指标，则边缘退化模型的公式(1)和公式(2)被转化为粒子滤波方法需要的状态转移方程和量测方程，其中，状态转移方程为f(xk,j|xk,j-1)：
Figure BDA0002805559250000101
量测方程h(yk,j|xk,j)为：
Figure BDA0002805559250000102
当公式(3)和(4)被用于多个退化指标时，令xj＝[x1,j,…,xK,j]T表示所有退化指标在tj时刻的状态值，令yj＝[y1,j,…,yK,j]T代表所有退化指标在tj时刻的量测值，则状态转移方程和量测方程转化为f(xj|xj-1)和h(yj|xj)；但此时的多元退化模型f(xj|xj-1)和h(yj|xj)未建立多个退化指标之间的依赖关系；
在多个退化指标同时退化的情况中，假设设备的失效时刻是：先到达对应失效阈值的某退化指标到达该退化指标的失效阈值的时刻，因此对一个具有K个退化指标的复杂设备，每个退化指标的失效阈值定义为Hk,k＝1,2,…,K，该设备在tj时刻开始预测时的剩余使用寿命为：
Figure BDA0002805559250000103
步骤2：外部参数离线估计；
外部参数
Figure BDA0002805559250000104
和
Figure BDA0002805559250000105
使用最大似然估计方法利用训练集数据进行离线估计；为了区分训练集变量以及测试集变量，训练集的变量用上标'*'表示，令
Figure BDA0002805559250000106
Figure BDA0002805559250000107
及
Figure BDA0002805559250000108
分别表示训练集中的设备的第k个退化指标的退化状态、量测值以及退化轨迹，则
Figure BDA0002805559250000109
和
Figure BDA00028055592500001010
的分布由标准布朗运动过程的性质推得：
Figure BDA00028055592500001011
Figure BDA00028055592500001012
其中，Q＝[min{ti,tj}]1≤i,j≤M，IM是一个M阶单位矩阵；
外部参数的似然函数如下：
Figure BDA00028055592500001013
Figure BDA00028055592500001014
训练集设备中
Figure BDA00028055592500001015
里的量测值是已知的，然而
Figure BDA00028055592500001016
中的退化状态是隐含的。为了估计外部参数，首先需要用合适的方法将退化状态
Figure BDA00028055592500001115
识别出来，例如简单滑动平均(simplemoving average,SMA)，加权滑动平均(weighted moving average,WMA)等等，然后即可通过公式(9)使用最大似然估计得到参数
Figure BDA0002805559250000111
选择一个退化趋势与实际量测值的趋势一致的函数μk(·)描述第k个退化指标的退化轨迹，然后即可使用μk(·)对训练集的退化状态
Figure BDA0002805559250000112
进行最小二乘拟合，估计得到训练集每个设备的退化轨迹参数
Figure BDA0002805559250000113
参数
Figure BDA0002805559250000114
即可通过使公式(8)的似然函数最大来估计得到。
步骤3：测试集设备退化轨迹参数与退化指标状态的同时在线估计；
步骤3.1：基于参数相关性的多元退化模型建立；
对于不同退化指标之间的依赖关系是由其退化速率之间的相关性导致的复杂设备，假设模型(1)中的退化轨迹参数θ1,θ2,…,θK之间线性相关，则不同退化指标X1(t),X2(t),…,XK(t)之间的复杂依赖关系由线性相关的退化轨迹参数以及θk和Xk(t)之间的线性或非线性的映射(公式(1)中的线性或非线性函数μk(θk,t))共同构建和描述。因此，首先需要估计θ1,θ2,…,θK之间的皮尔逊相关系数，假设退化轨迹被设置为μk(θk,t)＝akexp(bk·t)，θk＝[ak,bk]T，则
Figure BDA0002805559250000115
和
Figure BDA0002805559250000116
之间的皮尔逊相关系数
Figure BDA0002805559250000117
以及
Figure BDA0002805559250000118
和
Figure BDA0002805559250000119
之间的皮尔逊相关系数
Figure BDA00028055592500001110
需要由训练集数据
Figure BDA00028055592500001111
估计得到，而
Figure BDA00028055592500001112
由最小二乘估计得到，其中1≤k1,k2≤K；
使用基于核平滑的粒子滤波方法对测试集设备的多个退化指标及其对应的退化轨迹参数进行在线估计，首先需要建立合理的基于参数间相关性的状态转移方程，通过在状态转移方程中构建一个协方差矩阵来实现建立参数之间相关性的目的。直接从一个多元高斯分布中进行采样，得到多个退化指标的退化轨迹参数的粒子群，而该多元高斯分布的协方差矩阵体现了多个退化指标参数之间的相关关系。
令Δuj＝[Δμ1(tj),…,ΔμK(tj)]T，ωj＝[ω1,j,…,ωK,j]T，并且令
Figure BDA00028055592500001113
Figure BDA00028055592500001114
分别为参数a＝[a1,…,aK]T，参数b＝[b1,…,bK]T和多个退化指标x＝[X1(t),…,XK(t)]T在tj时刻的第i个粒子。令θj＝[θ1,j,…,θK,j]T，θk,j＝[ak,j,bk,j]T，则基于参数相关性的多个退化指标的状态转移方程f(xj,θj|xj-1,θj-1)如下：
Figure BDA0002805559250000121
其中：
Figure BDA0002805559250000122
s是处于0到1之间的核参数，
Figure BDA0002805559250000123
Figure BDA0002805559250000124
Figure BDA0002805559250000125
是tj时刻x，a和b的粒子与其对应的粒子权值。
Figure BDA0002805559250000126
和
Figure BDA0002805559250000127
分别是
Figure BDA0002805559250000128
的加权均值和加权方差，
Figure BDA0002805559250000129
和
Figure BDA00028055592500001210
分别是
Figure BDA00028055592500001211
的加权均值和加权方差。
从公式(10)中看到，当粒子
Figure BDA00028055592500001212
从粒子滤波中的建议分布(proposaldistribution)中采样时，建议分布为两个多元高斯分布，并且表征参数之间相关关系的
Figure BDA00028055592500001213
和
Figure BDA00028055592500001214
均包含在该多元高斯分布的协方差矩阵Va,j-1和Vb,j-1中。因此，从建议分布中采样得到的参数粒子之间是具有相关关系的。
步骤3.2：KS-PF实现退化轨迹参数与多退化指标的同时在线估计
使用粒子滤波中的贯序重要性重采样(sequential importance resampling,SIR)框架，按照状态转移方程(10)和量测方程(4)，对测试集设备的多指标退化状态x＝[X1(t),…,XK(t)]T和退化轨迹参数θk,j＝[ak,j,bk,j,…]T,k＝1,…,K进行在线估计，具体步骤如下：
1)初始化：
令j＝0，设置t0时刻的设备的多指标退化状态和退化轨迹参数的初始分布p(x0)和p(θ0)；
2)状态估计：
令j加1，若存在量测向量yj，顺序执行以下步骤，否则终止该程序；
2.1)预测：从建议分布p(xj,θj|xj-1,θj-1)中按照状态转移方程(10)采样得到粒子集
Figure BDA00028055592500001215
xj和θj的先验分布，可用该粒子集表示为：
Figure BDA0002805559250000131
Figure BDA0002805559250000132
2.2)更新：使用tj时刻新到来的退化指标量测向量yj按照量测方程(4)对粒子权值进行更新：
Figure BDA0002805559250000133
并且该权值被归一化为：
Figure BDA0002805559250000134
2.3)重采样：有效粒子数被计算为
Figure BDA0002805559250000135
若Nnef小于一个给定的阈值，则执行重采样步骤，否则不执行。此时得到的粒子集为
Figure BDA0002805559250000136
xj和θj的后验分布即可用该粒子集表示为：
Figure BDA0002805559250000137
Figure BDA0002805559250000138
2.4)迭代执行步骤2.1)，即得到状态估计；
步骤3.3：最优调整OT方法实现核参数s的自适应在线调整；
公式(10)中的参数
Figure BDA0002805559250000139
s即为KS-PF中的核参数，取值为0到1之间的常数。s的取值影响了核平滑步骤中参数粒子收缩和扰动的程度，即s的不同取值会导致从建议分布中采样得到不同的参数粒子集，从而影响KS-PF的在线状态估计准确度。采用最优调整(Optimized Tuning,OT)方法，在每个状态估计时刻tj，均选择一个最优的核参数sj。
粒子滤波方法的在线估计准确度与建议分布的选择有关，一个最优的建议分布应该使得更新后的粒子权值方差最小，这也可以理解为：一个好的建议分布应该使得先验分布与后验分布尽可能的相似。KL散度(Kullback-Leibler divergence)是一种度量两个概率分布相似程度的指标，因此OT方法在每个tj时刻，均选择使得先验分布与后验分布的KL散度最小的核参数sj。
在tj时刻，设备多个退化指标的多元退化状态的先验分布Qj＝p(xj|y1:j-1)表示为公式(11)，后验分布Pj＝p(xj|y1:j)表示为公式(15)，先验分布与后验分布的KL散度计算为：
Figure BDA0002805559250000141
贝叶斯滤波中的更新公式为：
Figure BDA0002805559250000142
将公式(18)和公式(11)带入公式(17)，并且使用狄拉克函数δ(·)的挑选性质，得到：
Figure BDA0002805559250000143
将公式(13)带入公式(14)，得到：
Figure BDA0002805559250000144
最后，将公式(20)带入公式(19)，得到：
Figure BDA0002805559250000145
因此，在每个tj时刻，在SIR框架中用到的状态转移方程(10)里的核参数s均使用最优的核参数：
Figure BDA0002805559250000146
这样，使用该部分介绍的OTKS-PF方法，用状态转移方程(10)和量测方程(4)，按照粒子滤波的SIR框架，即可在每个tj时刻均使用公式(13)、(14)得到更新后的粒子集
Figure BDA0002805559250000147
得到的粒子
Figure BDA0002805559250000148
代表了K个退化指标的状态估计粒子值，并且这样得到的
Figure BDA0002805559250000149
和
Figure BDA00028055592500001410
对于任意k1,k2∈[1,K]都不是独立的，具有一定的依赖关系。多个退化指标的多元联合分布p(x1,j,…,xK,j|y1:j)即可得到；
步骤4：多退化指标的复杂设备剩余寿命预测
根据步骤3.2和步骤3.3描述的“最优调整核参数的核平滑(OTKS)的基于核平滑的粒子滤波(OTKS-PF)实现多个退化指标的状态在线估计方法，在每个tj时刻均可得到粒子集
Figure BDA00028055592500001411
则多个退化指标的多元联合分布可被估计为：
Figure BDA0002805559250000151
在tj时刻开始预测时，将估计得到的每个粒子代表的参数和系统状态均代入模型(3)中，则可以通过设备失效定义得到基于每个粒子预测的剩余使用寿命；基于第i个粒子预测得到的设备剩余使用寿命为：
Figure BDA0002805559250000152
在tj时刻开始预测时，预测得到的设备的剩余使用寿命的概率密度函数(probability density function,PDF)为：
Figure BDA0002805559250000153
即求得剩余使用寿命的概率密度函数，即上式中的p(lj|y1:j)。
1.发动机退化数据介绍与数据预处理
使用NASA发布的C-MAPSS软件仿真的飞机涡轮发动机的退化数据进行多退化指标方法的验证和分析。该发动机退化数据集包含100个同系列单工况单故障模式的发动机的退化数据，每个发动机均是从某个特定状态一直运行到系统故障。每个发动机的退化状态均通过21个传感器得到的时间序列量测数据来反映。
(1)数据预处理及多退化指标的构造
首先，由于发动机原始的测量信号量纲不一样，因此对原始测量数据进行z-score标准化。之后需要构造发动机的退化指标，原始数据中有21个传感器的量测数据，对原始数据使用主成分分析(principal components analysis,PCA)进行降维处理以去除噪声数据保留有用信息。
经过PCA处理及分析，发现仅需二个维度即可保留原始数据87.6％的信息量，因此将原始数据降为二维，并且将降维得到的二维数据作为退化指标1和退化指标2来共同反映该发动机的退化情况。
对每个退化指标建立如公式(1)和(2)的边缘退化模型，由于在该发动机数据中使用两个退化指标进行设备的RUL预测，因此公式(1)和(2)中的k＝1,2代表第k个退化指标。
这样，经过PCA降维后得到数据即作为该发动机的退化指标1和退化指标2的量测值Yk(t),k＝1,2，再将降维后得到的退化指标1和退化指标2的退化状态量测值均进行简单滑动平均滤波，得到相对更平滑的曲线，滑动窗口大小设置为25；将滤波后的曲线作为发动机的两个退化指标的退化状态值Xk(t),k＝1,2。图2展示了1号发动机通过PCA降维后得到的两个退化指标的量测值Yk(t)及经过滑动平均后得到的退化状态值Xk(t)。
(2)失效阈值及设备失效的定义
定义Hk,k＝1,2分别为两个退化指标的失效阈值。这里，对退化指标1设定失效阈值H1＝1.22，对退化指标2设定失效阈值H2＝0.85。定义一个发动机的真实失效时刻tEOL为：发动机发生退化时，先退化至其对应的失效阈值的指标的退化状态值Xk(t)到达其阈值的时刻。
在实验验证部分，取一部分发动机作为训练集，其他发动机作为测试集，训练集发动机的全寿命的退化数据均视为已知，而测试集发动机则视为仅已知部分退化数据。对于所有100个发动机，采用k-fold来验证我们的方法，k值取5，即将100个发动机随机分为5组，每组20个发动机。将实验重复五次，称为5组实验，每次取其中一组的20个发动机数据做测试集，其余80个发动机为训练集。在每组实验中，均使用训练集数据进行参数估计，并对测试集发动机在不同时刻进行多次的剩余寿命预测。由于篇幅原因，以下分析均以第一组实验为例，我们也给出了所有5组实验的预测结果，并进行了结果对比分析。
2.边缘退化模型的建立及外部参数估计；
对于建立的公式(1)和(2)的边缘退化模型，首先要确定退化轨迹μk(θk,t)的形式。通过对两个退化指标状态值的退化轨迹进行观察，看到其呈现指数形式的退化，因此本文为两个退化指标建立如下的指数模型：
μk(θk,t)＝ak+bk·exp(ck·t) (47)
对训练集发动机已知的两个退化指标的量测值Yk(t),k＝1,2进行滑动窗口为25的滑动平均处理，得到其退化状态值Xk(t),k＝1,2。图3为使用公式(26)的模型对所有训练集发动机的两个退化指标的状态值进行最小二乘拟合后，对拟合的退化轨迹参数进行统计得到的箱线图。因此对参数b1,b2取其各自的中位数作为固定值：b1＝0.085,b2＝0.056。这样，模型(26)变为：
Figure BDA0002805559250000161
图4展示了第一组实验中的2个训练集发动机(第1号和第8号发动机)使用公式(27)进行拟合时的模型拟合效果。从图4中可以看到，使用公式(27)的退化轨迹模型对两个退化指标的退化状态进行拟合，得到的拟合结果较好。之后，使用训练集数据以及公式(8)和(9)，对退化模型的外部参数
Figure BDA0002805559250000171
进行离线估计。表1展示了第一组实验中，使用80个训练集数据得到的外部参数的估计结果。
表1第一组实验的外部参数估计结果
Figure BDA0002805559250000172
3.测试集发动机退化轨迹参数在线估计及剩余寿命预测
按照公式(10)的形式，建立两个退化指标的状态转移方程f(xj,cj|xj-1,cj-1)如下：
Figure BDA0002805559250000173
其中，k＝1,2代表第k个退化指标。
Figure BDA0002805559250000174
b1＝0.085，b2＝0.056(对第一组实验)，
Figure BDA00028055592500001715
是tj-1时刻得到的参数向量c＝[c1,c2]T的粒子集及每个粒子的对应权值。,
Figure BDA0002805559250000177
和
Figure BDA0002805559250000178
分别是ck,j-1的基于粒子集
Figure BDA0002805559250000179
的加权均值向量和加权协方差矩阵。
从公式(28)可以看到在状态转移的过程中仅包含了状态的增量，因此虽然此时θk＝[ak,ck]T，但在OTKS-PF的状态参数同估计时，仅需估计参数ck。参数c1和c2之间的皮尔逊相关系数由训练集计算为
Figure BDA00028055592500001710
图5展示了训练集的参数
Figure BDA00028055592500001711
和
Figure BDA00028055592500001712
的散点图。从计算得到的参数间皮尔逊相关系数以及图5中可以看出，不同退化指标之间的依赖关系是明显的，并且这种依赖关系可以通过参数之间的相关关系来描述。
使用OTKS-PF方法对多个退化指标的状态值X1(t),X2(t)及退化轨迹参数c1,c2进行在线估计，状态和参数的初始分布被设置为高斯分布，即：x0～N(μx,Σx)，c0～N(μc,Σc)，其中
Figure BDA00028055592500001713
Σx＝diag(σx,1,σx,2)，Σc＝diag(σc,1,σc,2)。
Figure BDA00028055592500001714
和σx,k分别是训练集设备退化状态
Figure BDA0002805559250000181
的均值和方差，这里的
Figure BDA0002805559250000182
是t1时刻第k个退化指标的退化状态。
Figure BDA0002805559250000183
和σc,k分别是训练集设备退化轨迹参数
Figure BDA0002805559250000184
的均值和方差。
图6展示了在第一组实验中测试集第5号发动机使用该方法进行双指标联合预测及单指标预测的参数估计效果。可以看到，该基于参数相关性的双指标联合预测方法与单指标的预测方法效果相近；在OTKS-PF的状态参数同估计中，随着已知观测点的增多，估计得到的参数均能逐渐收敛到参数的真实值。
图7和图8分别展示了第一组实验中测试集第19号发动机和测试集第68号发动机的寿命预测结果。其中，位于“RUL”轴与“开始预测时刻”轴组成的平面内的黄线代表设备的真实RUL值，(a)图中的蓝线，黑线，红线分别代表使用退化指标1，退化指标2以及使用基于参数的联合预测方法预测得到的RUL的概率密度曲线，而蓝点，黑点和红点则是对应的预测得到的RUL期望值。
图7和图8中的(a)图展示了分别使用退化指标1、退化指标2以及使用基于参数的联合预测方法预测得到的RUL的概率密度曲线以及预测值；(b)图展示了基于参数相关性的多指标联合估计的RUL预测结果及95％置信区间。观察图7和图8可以得到以下分析及结论：
(1)基于参数相关性的粒子滤波多退化指标设备RUL预测方法可以同时估计两个退化指标的状态值，并且跟踪到先到达其失效阈值的退化指标，从而预测得到准确的设备RUL的后验概率。从图7和图8的(a)图中可以看到，对于19号发动机，联合估计方法得到的预测PDF随着已知量测点的增加，不断地跟踪到退化指标2预测得到的RUL的PDF；而对于第68号发动机，联合估计的预测PDF不断的跟踪到指标1预测得到的RUL的PDF。
这是因为对于第19号发动机的两个退化指标，退化指标2首先退化至其对应的失效阈值，该时刻是第19号发动机的真实失效时刻。而对于第第68号发动机，退化指标1首先退化至其对应的失效阈值，该时刻是第68号发动机的真实失效时刻。因此，对于同一批设备，无论其哪个退化指标先退化至其对应的失效阈值，本发明提出的联合预测的方法均能跟踪到该退化指标的预测结果，这证明了本发明提出的基于参数相关性的多指标联合预测方法的有效性。
(2)基于参数相关性的粒子滤波多退化指标设备RUL预测方法随着开始预测时刻的增长，其RUL预测的精确度不断提升。从图7和图8的(b)图可以看到，随着已知观测点的增多，基于参数相关性的联合估计方法预测得到的RUL的95％置信区间的宽度不断减小，预测的精度不断提高，并且预测的准确度也较高。这符合使用粒子滤波实现RUL预测的逻辑，因为已知的量测值越多，则得到的先验知识越多，因此预测也会更准确，更精确。这说明了本发明所提出的方法具有实际的应用价值。
为了更直观的对比使用单个退化指标实现RUL预测与双指标联合估计实现RUL预测的预测效果，在每组实验中，对所有的测试集发动机从其寿命终止时刻tEOL的三分之二时刻开始预测，之后每隔1个时刻点预测一次。这样，每个测试集发动机均能得到若干个寿命预测的结果，再对该结果取平均。图9展示了第一组实验的预测误差与95％置信区间的平均值。由图9可以看出，在第一组实验的20个测试发动机中，基于参数相关性的联合估计的预测误差基本低于退化指标1和退化指标2的单个估计的预测误差，这符合分析的结果。并且联合估计的95％置信区间宽度也均低于单指标估计得到的95％置信区间宽度。
表2五组实验的单指标预测与多指标联合预测结果
Figure BDA0002805559250000191
对于每一组实验，对所有20个发动机的平均预测误差和平均95％置信区间宽度算平均值，表2展示了所有5组实验的预测结果对比情况。从表2中可以看到，当考虑同一批次具有个体差异性的不同测试集设备时，本专利提出的基于参数相关性的多退化指标设备RUL预测方法在预测准确度(由预测误差体现)和预测精度(由95％置信区间体现)上都比单指标的预测方法有较大的优势，这说明本发明提出的方法能够对具有多个退化指标的大型复杂设备进行有效且准确的RUL预测。

Claims (5)

1.一种基于多退化指标的复杂设备剩余寿命预测方法，其特征在于包括下述步骤：
步骤1：边缘退化模型的建立；
假设设备有K个退化指标，为每个退化指标建立维纳过程的边缘退化模型：
Xk(t)＝μk(θk,t)+σkB(t) (1)
Figure FDA0002805559240000011
公式(1)为状态转移方程，公式(2)为量测方程，其中，k＝1,2…K代表第k个退化指标，Xk(t)是第k个退化指标在t时刻的退化状态，μk(θk,t)是一个单调可微函数，被称为退化轨迹，为线性或非线性函数；θk＝[ak,bk,…]T为退化轨迹参数，决定指标的退化速率，个体不确定性指的是同一批次的不同设备之间的退化多样性，通过对退化指标的状态值进行回归得到，B(t)是一个标准布朗运动过程，σk是波动系数；σkB(t)具有正态增量，即有
Figure FDA0002805559240000012
其中0≤t2≤t1，描述了第k个指标在退化过程中的内部不确定性，即指标在退化过程中本身的退化随机性；Yk(t)代表第k个退化指标的量测值，其中
Figure FDA0002805559240000013
是服从高斯分布的随机变量εk的方差，εk代表量测噪声随机变量，体现了量测不确定性，服从一个零均值的正态分布，对退化指标建立随机过程模型，则退化模型中的每个部分均代表了一定的物理含义；
令xk,j＝Xk(tj)，yk,j＝Yk(tj)，Δt＝tj-tj-1，其中0≤j≤M，M为待预测设备中退化指标量测点的个数，仅考虑第k个退化指标，则边缘退化模型的公式(1)和公式(2)被转化为粒子滤波方法需要的状态转移方程和量测方程，其中，状态转移方程为f(xk,j|xk,j-1)：
Figure FDA0002805559240000014
量测方程h(yk,j|xk,j)为：
Figure FDA0002805559240000015
当公式(3)和(4)被用于多个退化指标时，令xj＝[x1,j,…,xK,j]T表示所有退化指标在tj时刻的状态值，令yj＝[y1,j,…,yK,j]T代表所有退化指标在tj时刻的量测值，则状态转移方程和量测方程转化为f(xj|xj-1)和h(yj|xj)；此时的多元退化模型f(xj|xj-1)和h(yj|xj)未建立多个退化指标之间的依赖关系；
在多个退化指标同时退化的情况中，假设设备的失效时刻是先到达对应失效阈值的某退化指标到达该退化指标的失效阈值的时刻，因此对一个具有K个退化指标的复杂设备，每个退化指标的失效阈值定义为Hk,k＝1,2,…,K，该设备在tj时刻开始预测时的剩余使用寿命为：
Figure FDA0002805559240000021
步骤2：外部参数离线估计；
外部参数
Figure FDA0002805559240000022
和
Figure FDA0002805559240000023
使用最大似然估计方法利用训练集数据进行离线估计；为了区分训练集变量以及测试集变量，训练集的变量用上标'*'表示，令
Figure FDA0002805559240000024
Figure FDA0002805559240000025
及
Figure FDA0002805559240000026
分别表示训练集中的设备的第k个退化指标的退化状态、量测值以及退化轨迹，则
Figure FDA0002805559240000027
和
Figure FDA0002805559240000028
的分布由标准布朗运动过程的性质推得：
Figure FDA0002805559240000029
Figure FDA00028055592400000210
其中，Q＝[min{ti,tj}]1≤i,j≤M，IM是一个M阶单位矩阵；
外部参数的似然函数如下：
Figure FDA00028055592400000211
Figure FDA00028055592400000212
为了估计外部参数，首先将退化状态
Figure FDA00028055592400000213
识别出来，然后通过公式(9)使用最大似然估计得到参数
Figure FDA00028055592400000214
选择一个退化趋势与实际量测值的趋势一致的函数μk(·)描述第k个退化指标的退化轨迹，然后使用μk(·)对训练集的退化状态
Figure FDA00028055592400000215
进行最小二乘拟合，估计得到训练集每个设备的退化轨迹参数
Figure FDA00028055592400000216
参数
Figure FDA00028055592400000217
通过使公式(8)的似然函数最大来估计得到；
步骤3：测试集设备退化轨迹参数与退化指标状态的同时在线估计；
步骤3.1：基于参数相关性的多元退化模型建立；
假设模型(1)中的退化轨迹参数θ1,θ2,…,θK之间线性相关，则不同退化指标X1(t),X2(t),…,XK(t)之间的复杂依赖关系由线性相关的退化轨迹参数以及θk和Xk(t)之间的线性或非线性的映射共同构建和描述；首先估计θ1,θ2,…,θK之间的皮尔逊相关系数，假设退化轨迹被设置为μk(θk,t)＝akexp(bk·t)，θk＝[ak,bk]T，则
Figure FDA0002805559240000031
和
Figure FDA0002805559240000032
之间的皮尔逊相关系数
Figure FDA0002805559240000033
以及
Figure FDA0002805559240000034
和
Figure FDA0002805559240000035
之间的皮尔逊相关系数
Figure FDA0002805559240000036
由训练集数据
Figure FDA0002805559240000037
估计得到，而
Figure FDA0002805559240000038
由最小二乘估计得到，其中1≤k1,k2≤K；
使用基于核平滑的粒子滤波方法对测试集设备的多个退化指标及其对应的退化轨迹参数进行在线估计，首先建立基于参数间相关性的状态转移方程，通过在状态转移方程中构建一个协方差矩阵来实现建立参数之间相关性的目的，从一个多元高斯分布中进行采样，得到多个退化指标的退化轨迹参数的粒子群，而该多元高斯分布的协方差矩阵体现了多个退化指标参数之间的相关关系；
步骤3.2：KS-PF实现退化轨迹参数与多退化指标的同时在线估计
使用粒子滤波中的贯序重要性重采样框架，按照状态转移方程(10)和量测方程(4)，对测试集设备的多指标退化状态x＝[X1(t),…,XK(t)]T和退化轨迹参数θk,j＝[ak,j,bk,j,…]T,k＝1,…,K进行在线估计；
步骤3.3：最优调整OT方法实现核参数s的自适应在线调整；
采用最优调整(Optimized Tuning,OT)方法，在每个状态估计时刻tj，均选择一个最优的核参数sj；
在tj时刻，设备多个退化指标的多元退化状态的先验分布Qj＝p(xj|y1:j-1)表示为公式(11)，后验分布Pj＝p(xj|y1:j)表示为公式(15)，先验分布与后验分布的KL散度计算为：
Figure FDA0002805559240000039
贝叶斯滤波中的更新公式为：
Figure FDA00028055592400000310
将公式(18)和公式(11)带入公式(17)，并且使用狄拉克函数δ(·)的挑选性质，得到：
Figure FDA00028055592400000311
将公式(13)带入公式(14)，得到：
Figure FDA00028055592400000312
最后，将公式(20)带入公式(19)，得到：
Figure FDA0002805559240000041
因此，在每个tj时刻，在SIR框架中用到的状态转移方程(10)里的核参数s均使用最优的核参数：
Figure FDA0002805559240000042
用状态转移方程(10)和量测方程(4)，按照粒子滤波的SIR框架，即在每个tj时刻均使用公式(13)、(14)得到更新后的粒子集
Figure FDA0002805559240000043
得到的粒子
Figure FDA0002805559240000044
代表了K个退化指标的状态估计粒子值，得到的
Figure FDA0002805559240000045
Figure FDA0002805559240000046
和
Figure FDA0002805559240000047
对于任意k1,k2∈[1,K]都不是独立的，具有一定的依赖关系，多个退化指标的多元联合分布p(x1,j,…,xK,j|y1:j)即可得到；
步骤4：多退化指标的复杂设备剩余寿命预测；
在每个tj时刻均得到粒子集
Figure FDA0002805559240000048
则多个退化指标的多元联合分布可被估计为：
Figure FDA0002805559240000049
在tj时刻开始预测时，将估计得到的每个粒子代表的参数和系统状态均代入模型(3)中，则通过设备失效定义得到基于每个粒子预测的剩余使用寿命；基于第i个粒子预测得到的设备剩余使用寿命为：
Figure FDA00028055592400000410
在tj时刻开始预测时，预测得到的设备的剩余使用寿命的概率密度函数为：
Figure FDA00028055592400000411
即求得剩余使用寿命的概率密度函数，即式(25)中的p(lj|y1:j)。
2.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述步骤3.1的基于参数相关性的多元退化模型建立的具体步骤如下：
令Δuj＝[Δμ1(tj),…,ΔμK(tj)]T，ωj＝[ω1,j,…,ωK,j]T，并且令
Figure FDA0002805559240000051
Figure FDA0002805559240000052
分别为参数a＝[a1,…,aK]T，参数b＝[b1,…,bK]T和多个退化指标x＝[X1(t),…,XK(t)]T在tj时刻的第i个粒子。令θj＝[θ1,j,…,θK,j]T，θk,j＝[ak,j,bk,j]T，则基于参数相关性的多个退化指标的状态转移方程f(xj,θj|xj-1,θj-1)如下：
Figure FDA0002805559240000053
其中，
Figure FDA0002805559240000054
Figure FDA0002805559240000055
是tj时刻x，a和b的粒子与其对应的粒子权值，
Figure FDA0002805559240000056
和
Figure FDA0002805559240000057
分别是
Figure FDA0002805559240000058
的加权均值和加权方差，
Figure FDA0002805559240000059
和
Figure FDA00028055592400000510
分别是
Figure FDA00028055592400000511
的加权均值和加权方差；
当粒子
Figure FDA00028055592400000512
从粒子滤波中的建议分布中采样时，建议分布为两个多元高斯分布，并且表征参数之间相关关系的
Figure FDA00028055592400000513
和
Figure FDA00028055592400000514
均包含在该多元高斯分布的协方差矩阵Va,j-1和Vb,j-1中，从建议分布中采样得到的参数粒子之间是具有相关关系的。
3.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述步骤3.2中KS-PF实现退化轨迹参数与多退化指标的同时在线估计的具体步骤如下：
1)初始化：
令j＝0，设置t0时刻的设备的多指标退化状态和退化轨迹参数的初始分布p(x0)和p(θ0)；
2)状态估计：
令j加1，若存在量测向量yj，顺序执行以下步骤，否则终止该程序；
a)预测：从建议分布p(xj,θj|xj-1,θj-1)中按照状态转移方程(10)采样得到粒子集
Figure FDA0002805559240000061
xj和θj的先验分布，用该粒子集表示为：
Figure FDA0002805559240000062
Figure FDA0002805559240000063
b)更新：使用tj时刻新到来的退化指标量测向量yj按照量测方程(4)对粒子权值进行更新：
Figure FDA0002805559240000064
并且该权值被归一化为：
Figure FDA0002805559240000065
c)重采样：有效粒子数被计算为
Figure FDA0002805559240000066
若Nneff小于一个给定的阈值，则执行重采样步骤，否则不执行，得到的粒子集为
Figure FDA0002805559240000067
xj和θj的后验分布用该粒子集表示为：
Figure FDA0002805559240000068
Figure FDA0002805559240000069
d)迭代执行步骤a)-步骤c)，即得到状态估计。
4.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述步骤2中，退化状态
Figure FDA00028055592400000611
识别采用简单滑动平均或加权滑动平均。
5.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述
Figure FDA00028055592400000610
s是处于0到1之间的核参数。
CN202011366448.7A 2020-11-29 2020-11-29 一种基于多退化指标的复杂设备剩余寿命预测方法 Active CN112487694B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011366448.7A CN112487694B (zh) 2020-11-29 2020-11-29 一种基于多退化指标的复杂设备剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011366448.7A CN112487694B (zh) 2020-11-29 2020-11-29 一种基于多退化指标的复杂设备剩余寿命预测方法

Publications (2)

Publication Number Publication Date
CN112487694A CN112487694A (zh) 2021-03-12
CN112487694B true CN112487694B (zh) 2022-09-23

Family

ID=74936802

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011366448.7A Active CN112487694B (zh) 2020-11-29 2020-11-29 一种基于多退化指标的复杂设备剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN112487694B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114626449A (zh) * 2022-03-03 2022-06-14 南京航空航天大学 一种基于最小二乘确权法的对抗环境下目标识别网络退化分析方法
CN116738862B (zh) * 2023-07-12 2024-05-31 南方电网电力科技股份有限公司 一种锅炉结垢故障预测方法、装置、设备和介质
CN117291445B (zh) * 2023-11-27 2024-02-13 国网安徽省电力有限公司电力科学研究院 一种综合能源系统下基于状态转移的多目标预测方法
CN117829002B (zh) * 2024-03-05 2024-05-14 深圳市明谋科技有限公司 电力线缆的老化诊断监测方法及系统

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111160666B (zh) * 2020-01-02 2023-06-23 西北工业大学 强噪声与非周期状态监测的健康状态与可靠性评估方法
CN111368403B (zh) * 2020-02-24 2022-03-08 西安交通大学 一种自适应非线性退化剩余寿命预测方法

Also Published As

Publication number Publication date
CN112487694A (zh) 2021-03-12

Similar Documents

Publication Publication Date Title
CN112487694B (zh) 一种基于多退化指标的复杂设备剩余寿命预测方法
US8725456B1 (en) Decomposition technique for remaining useful life prediction
Ahwiadi et al. An adaptive particle filter technique for system state estimation and prognosis
US7580812B2 (en) Trending system and method using window filtering
Zimmer et al. Safe active learning for time-series modeling with gaussian processes
CN111144458A (zh) 子空间嵌入特征分布对齐的不同工况下机械故障识别方法
CN109917777B (zh) 基于混合多采样率概率主成分分析模型的故障检测方法
CN112800616B (zh) 基于比例加速退化建模的设备剩余寿命自适应预测方法
CN110757510B (zh) 一种机器人剩余寿命预测方法及系统
CN114819054A (zh) 一种基于物理信息神经网络的电力电子系统状态监测方法
Lim et al. A novel time series-histogram of features (TS-HoF) method for prognostic applications
CN112613617A (zh) 基于回归模型的不确定性估计方法和装置
CN116384224A (zh) 一种基于条件化参数动态卷积神经网络的航空发动机寿命预测方法
Kotenko et al. An approach for intelligent evaluation of the state of complex autonomous objects based on the wavelet analysis
Lin et al. A novel product remaining useful life prediction approach considering fault effects
Liao et al. Nonparametric and semi-parametric sensor recovery in multichannel condition monitoring systems
CN115392056A (zh) 一种高压架空输电线路运行状态监测及预警方法及装置
CN111309973A (zh) 基于改进马尔可夫模型和改进k最近邻的缺失值填补方法
Kakati et al. Remaining useful life predictions for turbofan engine degradation using online long short-term memory network
CN113468720B (zh) 一种数模联动的随机退化设备寿命预测方法
CN115375038A (zh) 一种飞机发动机失效模式识别和寿命预测方法
CN104200269A (zh) 一种基于在线学习最小嵌入维网络的实时故障诊断方法
CN111160464B (zh) 基于多隐层加权动态模型的工业高阶动态过程软测量方法
CN110060374B (zh) 一种飞机燃油系统异常检测方法及装置
Wang et al. Health indicator forecasting for improving remaining useful life estimation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant