CN113642158B - 基于Box-Cox变换的非线性设备剩余寿命预测方法 - Google Patents

基于Box-Cox变换的非线性设备剩余寿命预测方法 Download PDF

Info

Publication number
CN113642158B
CN113642158B CN202110825429.4A CN202110825429A CN113642158B CN 113642158 B CN113642158 B CN 113642158B CN 202110825429 A CN202110825429 A CN 202110825429A CN 113642158 B CN113642158 B CN 113642158B
Authority
CN
China
Prior art keywords
degradation
box
cox
equipment
random
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
CN202110825429.4A
Other languages
English (en)
Other versions
CN113642158A (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.)
Rocket Force University of Engineering of PLA
Original Assignee
Rocket Force University of Engineering of PLA
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 Rocket Force University of Engineering of PLA filed Critical Rocket Force University of Engineering of PLA
Priority to CN202110825429.4A priority Critical patent/CN113642158B/zh
Publication of CN113642158A publication Critical patent/CN113642158A/zh
Application granted granted Critical
Publication of CN113642158B publication Critical patent/CN113642158B/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • 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)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明属于寿命预测技术领域，具体公开了一种基于Box‑Cox变换的非线性设备剩余寿命预测方法，包括历史随机退化设备的未知参数估计和在线服役设备的剩余使用寿命预测两部分；本发明利用Box‑Cox变换对非线性退化数据做变换处理，提高了随机退化设备退化数据的线性程度，利用已有的线性Wiener过程建模该变换后的退化数据可得到剩余寿命的精确解析解，避免了已有时间尺度变换函数选择受制于随机退化设备失效数据的约束。

Description

基于Box-Cox变换的非线性设备剩余寿命预测方法
技术领域
本发明属于剩余寿命预测技术领域，涉及一种关键设备剩余寿命预测方法，具体涉及一种基于Box-Cox变换的非线性设备剩余寿命预测方法。
背景技术
剩余寿命预测与健康管理技术是现代复杂工程系统、重大产品、重大设施提高运行可靠性、安全性、可维护性的关键技术，可为重大装备的长周期安全可靠运行提供重要保障。随着先进传感和状态监测技术的发展，获取能够反映设备健康状态的性能退化过程监测数据已成为可能。在此背景需求下，数据驱动的基于随机过程建模的随机退化设备剩余寿命预测技术在过去十余年得到了广泛关注和蓬勃发展，其中以标准布朗运动驱动的Wiener过程为代表的随机过程，由于其良好的数学特性和适于描述非单调退化特征而被广泛应用于剩余寿命预测。
对于复杂的随机退化设备，其退化过程表现出明显的非线性。针对非线性退化数据建模，一般采取两种思路：一是采用对数变换或时间尺度变换将非线性退化数据近似处理成线性退化数据，然后采用线性Wiener过程建模退化数据，进而可得到剩余寿命概率密度函数的精确解析解，常见时间尺度变换函数有幂函数、多项式函数等，具体函数形式的选择取决于具体对象的退化数据，只有当设备退化及寿命数据充足时是可行的，然而对于高可靠性装备，要获取大量退化及寿命数据是不现实的，同时基于对数变换的非线性退化数据的线性化技术，要求退化数据具有指数特征，因此其适用对象受限；二是基于一个一般性的非线性时变漂移系数函数，通过对非线性随机退化模型进行时间-空间变换，将非线性随机过程首达固定边界的问题转化为标准布朗运动首达时变边界的问题，区别于基于时间尺度变换、对数变换等非线性数据重构技术，然而基于该一般性的非线性Wiener过程模型只能得到剩余寿命的近似解。
发明内容
针对现有技术中存在的问题，本发明的目的在于提供一种基于Box-Cox变换的非线性设备剩余寿命预测方法，根据提出的基于Box-Cox变换，用线性Wiener过程实现非线性退化数据建模，同时得到剩余寿命概率密度函数精确解析解，提高预测精度。
为了达到上述目的，本发明采用以下技术方案予以实现。
基于Box-Cox变换的非线性设备剩余寿命预测方法，包括以下步骤：
步骤1，对于在线服役运行的非线性随机退化设备，收集在线服役设备从开始运行时刻t0到当前时刻tk的退化监测值x1，x2，…，xk，记为在线服役设备原始监测退化数据X1：k＝{x1，x2，...，xk}；使用Box-Cox变换对X1：k进行处理，得到经Box-Cox变换后的退化数据
Figure BDA0003173392380000021
其中，
Figure BDA0003173392380000022
为Box-Cox变换参数，其值通过历史随机退化设备的Box-Cox变换退化监测数据估计得到；
步骤2，采用线性Wiener随机过程对在线服役设备经Box-Cox变换后的退化数据进行建模，得到在线服役设备的随机退化模型，模型参数中退化量初值β和漂移系数θ服从正态分布，其超参数取值通过历史随机退化设备的退化监测数据估计得到；
步骤3，采用Bayesian方法对所述在线服役设备随机退化模型中的退化量初值β和漂移系数θ的后验分布进行更新，得到退化量初值β和漂移系数θ的联合后验分布，从而得到退化量初值β和漂移系数θ的联合后验概率密度函数；
步骤4，在所述在线服役设备经Box-Cox变换后的随机退化模型的基础上，基于首达时间，对在线服役设备在当前时刻tk的剩余寿命Lk进行描述；基于退化量初值β和漂移系数θ的联合后验概率密度函数，得到在线服役设备剩余使用寿命的概率密度函数和概率分布函数；随着所述在线服役设备经Box-Cox变换后的随机退化模型中参数集合
Figure BDA0003173392380000031
的更新，基于全概率公式得到在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解；将Lk的均值作为服役设备剩余使用寿命的预测值
Figure BDA0003173392380000032
进一步地，所述使用Box-Cox变换对X1：k进行处理，具体为：
使用Box-Cox变换对{x(t)，t≥0}进行处理，对于每一个具体的时间t，有
Figure BDA0003173392380000033
其中，x(t)＝xt；
经过Box-Cox变换后的在线服役设备在tk时刻的退化监测量记为
Figure BDA0003173392380000034
即
Figure BDA0003173392380000035
则在线服役设备经Box-Cox变换的监测序列记为
Figure BDA0003173392380000036
进一步地，所述在线服役设备随机退化模型为：
Figure BDA0003173392380000037
其中，σ为扩散系数，B(t)为是标准布朗运动，且有σB(t)～N(0，σ2t)。
进一步地，所述采用Bayesian方法对所述在线服役设备随机退化模型中的退化量初值β和漂移系数θ的后验分布进行更新，具体为：
3.1，由于β和θ是相互独立的，因此β和θ的联合先验分布为p(β，θ)＝p(β)p(θ)；对于给定的β和θ，基于标准布朗运动{B(t)，t≥0}的独立增量性和马尔科夫性，
Figure BDA0003173392380000041
服从多元正态分布，其联合概率密度函数为：
Figure BDA0003173392380000042
3.2，基于Bayesian定理，已知
Figure BDA0003173392380000043
前提下β和θ的联合后验分布为：
Figure BDA0003173392380000044
3.3，由于
Figure BDA0003173392380000045
和p(β，θ)是共轭的，从而使得后验估计
Figure BDA0003173392380000046
也是高斯的，即
Figure BDA0003173392380000047
则其联合后验概率密度函数为：
Figure BDA0003173392380000048
其中，
Figure BDA0003173392380000049
上式中的
Figure BDA00031733923800000410
参数值通过历史随机退化设备的Box-Cox变换退化监测数据估计得到，即
Figure BDA00031733923800000411
进一步地，所述基于首达时间，对在线服役设备在当前时刻tk的剩余寿命Lk进行描述为：
Figure BDA0003173392380000051
其中，
Figure BDA0003173392380000052
是在线服役设备经过Box-Cox变换后的失效阈值，定义为
Figure BDA0003173392380000053
ω为在线服役设备的失效阈值，其取值根据设计指南、工业标准或专家经验事先确定；lk为tk时刻的剩余寿命对应的随机变量；inf表示下确界。
进一步地，所述在线服役设备剩余使用寿命的概率密度函数和概率分布函数的获取过程为：
已知在线服役设备预测时刻为tk，同时给定θ、β和到tk时刻的经Box-Cox变换的退化监测数据集，利用标准布朗运动的性质，将剩余寿命Lk在首达时间下的描述改写为：
Figure BDA0003173392380000054
其中，W(lk)＝B(tk+lk)-B(tk)；
由于{W(lk)，lk≥0}仍然是一个标准布朗运动，因此{S(tk+lk，λ)，lk≥0}是一个带有线性漂移系数θlk的Wiener过程，其初始值为
Figure BDA0003173392380000055
首达失效阈值时间意义下带有线性漂移系数的Wiener过程描述的在线服役设备剩余寿命服从逆高斯分布，在已知β，θ和
Figure BDA0003173392380000056
的前提下，其剩余使用寿命的概率密度函数
Figure BDA0003173392380000057
和概率分布函数
Figure BDA0003173392380000058
分别为：
Figure BDA0003173392380000059
Figure BDA0003173392380000061
随着所述在线服役设备经Box-Cox变换后的随机退化模型中参数集合
Figure BDA0003173392380000062
的更新，基于全概率公式得到在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解为：
Figure BDA0003173392380000063
Figure BDA0003173392380000064
式中，
Figure BDA0003173392380000065
Figure BDA0003173392380000066
中各参数的取值为式(1*)经Bayesian参数在线更新得到的估计值；
最后，基于在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解公式得到在线服役设备剩余使用寿命的均值估计结果。
进一步地，还包括：基于在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解公式计算Lk的方差和置信区间，以量化剩余寿命预测结果的不确定性。
与现有技术相比，本发明的有益效果为：
本发明一是利用Box-Cox变换对非线性退化数据做变换处理，提高了随机退化设备退化数据的线性程度，利用已有的线性Wiener过程建模该变换后的退化数据可得到剩余寿命的精确解析解，避免了已有时间尺度变换函数选择受制于随机退化设备失效数据的约束；二是本发明基于Box-Cox变换对非线性随机退化设备剩余使用寿命预测进行研究，得到了非线性随机退化设备剩余寿命的精确解析解，避免直接采用非线性Wiener过程进行非线性退化数据建模，只能得到剩余使用寿命的近似解的问题。
附图说明
下面结合附图和具体实施例对本发明做进一步详细说明。
图1是本发明仿真实验的CS2-35、CS2-36、CS2-37和CS2-38原始监测退化轨迹；
图2是本发明仿真实验的基于Box-Cox变换后的CS2-35、CS2-36、CS2-37、CS2-38线性化退化轨迹；其中，(a)CS2-35，(b)CS2-36，(c)CS2-37，(d)CS2-38；
图3是本发明仿真实验的最后监测时刻为tk＝100个循环周期时Box-Cox变换退化路径拟合效果对比图；其中，(a)经过Box-Cox变换，(b)未经过Box-Cox变换；
图4是本发明仿真实验的最后监测时刻为tk＝400个循环周期时经过Box-Cox变换和未经过变换的退化路径拟合效果对比图；其中，(a)经过Box-Cox变换，(b)未经过Box-Cox变换；
图5是本发明仿真实验的基于本发明提出方法得到的CS2-36号锂电池剩余使用寿命概率密度函数及均值图。
具体实施方式
下面将结合实施例对本发明的实施方案进行详细描述，但是本领域的技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。
本发明提供的一种基于Box-Cox变换的非线性设备剩余寿命预测方法，包括：
第一部分：历史随机退化设备的参数估计
(1)假设有N个同类历史随机退化设备，每个退化设备有Mj(j＝1，2，...，N)个监测时间点个数，且将第j(1≤j≤N)个随机退化设备对应监测时刻表示为
Figure BDA0003173392380000081
则第j(1≤j≤N)个随机退化设备在
Figure BDA0003173392380000082
时刻的退化监测量记为
Figure BDA0003173392380000083
Figure BDA0003173392380000084
是第j(1≤j≤N)个随机退化设备的第qj个测量时刻，因此qj＝1，2，…，Mj。
Figure BDA0003173392380000085
记第j(1≤j≤N)个随机退化设备原始退化监测数据为
Figure BDA0003173392380000086
则N个同类随机退化设备原始退化监测数据记为
Figure BDA0003173392380000087
(2)历史退化数据Box-Cox变换
特殊地，{Xj(t)，t≥0}表示第j(1≤j≤N)个历史随机退化设备随着监测时间t推进的非线性退化过程，为了提高非线性退化过程的线性程度，本专利使用Box-Cox变换对{Xj(t)，t≥0}进行处理，对于每一个具体的时间t，有
Figure BDA0003173392380000088
其中，λ是Box-Cox变换参数。
经过Box-Cox变换后的j(1≤j≤N)个随机退化设备在
Figure BDA0003173392380000089
时刻的退化监测量记为
Figure BDA00031733923800000810
即
Figure BDA00031733923800000811
j＝1，2，...，N；qj＝1，2，...，Mj，则第j(1≤j≤N)个历史随机退化设备经Box-Cox变换的监测序列记为
Figure BDA00031733923800000812
同理，N个同类历史随机退化设备经Box-Cox变换的监测序列记为
Figure BDA00031733923800000813
(3)线性Wiener随机退化建模
用{S(t，λ)，t≥0}表示随机退化设备经过Box-Cox变换后的退化过程，基于线性Wiener建模该退化过程可得：
S(t，λ)＝β+θt+σB(t) (2)
其中，θ是漂移系数，用来表征设备退化率；σ是扩散系数；β是退化量初值；{B(t)，t≥0}是标准布朗运动，且有σB(t)～N(0，σ2t)。
考虑到历史随机退化设备的个体差异性，令β和θ为随机变量，且分别服从正态分布，即
Figure BDA0003173392380000091
基于以上描述，为实现基于线性Wiener随机过程的Box-Cox变换后退化监测数据建模，核心任务是对式(2)的未知参数集合进行估计，记未知参数集合为
Figure BDA0003173392380000092
(4)未知参数极大似然估计
由标准布朗运动的数学性质可得，不失一般性，第j(1≤j≤N)个历史随机退化设备经Box-Cox变换的监测序列
Figure BDA0003173392380000093
服从多元正态分布，即
Figure BDA0003173392380000094
其均值为：
Figure BDA0003173392380000095
其中，
Figure BDA0003173392380000096
相应的方差为：
Figure BDA0003173392380000097
其中，
Figure BDA0003173392380000098
假设N个历史随机退化设备的退化过程是相互独立的，可得该N个历史随机退化设备经Box-Cox变换的历史监测数据
Figure BDA0003173392380000099
的对数似然函数为：
Figure BDA00031733923800000910
其中，C是一个常数，其取值大小不影响随机退化设备剩余使用寿命估计结果。
为方便数学推导，对以下参数进行重构，有
Figure BDA0003173392380000101
基于式(6)参数重构后的未知参数集合记为
Figure BDA0003173392380000102
其联合对数似然函数为：
Figure BDA0003173392380000103
基于式(7)分别对μβ，μθ和
Figure BDA0003173392380000104
求偏导数有
Figure BDA0003173392380000105
Figure BDA0003173392380000106
Figure BDA0003173392380000107
为求得
Figure BDA0003173392380000108
的极大似然估计值，首先固定
Figure BDA0003173392380000109
值，让式(8)、(9)、(10)等于零，解偏微分方程组可得
Figure BDA00031733923800001010
的极大似然估计值，即：
Figure BDA00031733923800001011
其中，
Figure BDA0003173392380000111
接下来，将基于式(11)、式(12)得到的
Figure BDA0003173392380000112
的极大似然估计值
Figure BDA0003173392380000113
知
Figure BDA0003173392380000114
代入式(7)可得到
Figure BDA0003173392380000115
的边缘似然函数为：
Figure BDA0003173392380000116
同理，将式(13)所示边缘似然函数分别对
Figure BDA0003173392380000117
求偏导数并令其等于零，解方程即可得到到
Figure BDA0003173392380000118
的极大似然估计值。然后再基于
Figure BDA0003173392380000119
变换，可得到式(2)所描述的线性Wiener过程且考虑个体差异性的参数集合估计值，记为
Figure BDA00031733923800001110
第二部分：在线服役设备的剩余使用寿命的自适应预测
步骤1，对于在线服役运行的非线性随机退化设备，收集在线服役设备从开始运行时刻t0到当前时刻tk的退化监测值x1，x2，…，xk，记为在线服役设备原始监测退化数据X1：k＝{x1，x2，...，xk}；使用Box-Cox变换对X1：k进行处理，得到经Box-Cox变换后的退化数据
Figure BDA00031733923800001111
其中，
Figure BDA00031733923800001112
为Box-Cox变换参数，其值通过历史随机退化设备的Box-Cox变换退化监测数据估计得到；
对于在线服役设备，假设其退化过程在0＜t1＜...＜tk时刻被离散监测，得到的退化监测值记为x1，x2，…，xk，不失一般性，监测时间间隔为等间隔，即Δt＝tk-tk-1，xk＝X(tk)表示最后监测时刻tk得到的退化量监测值。基于步骤(1)(2)的退化数据描述和Box-Cox变换，记在线服役设备原始监测退化数据为X1：k＝{x1，x2，...，xk}，经Box-Cox变换后的数据为
Figure BDA0003173392380000121
其中
Figure BDA0003173392380000122
为步骤(4)通过极大似然估计得到的参数估计结果。
步骤2，采用线性Wiener随机过程对在线服役设备变换后的退化数据进行建模，得到在线服役设备经Box-Cox变换后的随机退化模型；
同理基线性Wiener退化过程可建模
Figure BDA0003173392380000123
Figure BDA0003173392380000124
其中各参数物理内涵同式(2)。
步骤3，采用Bayesian方法对所述在线服役设备经Box-Cox变换后的随机退化模型中的退化量初值β和漂移系数θ的后验分布进行更新，得到退化量初值β和漂移系数θ的联合后验分布，从而得到退化量初值β和漂移系数θ的联合后验概率密度函数；
为体现在线服役随机退化设备与历史随机退化设备之间的个体差异性，需要利用
Figure BDA0003173392380000125
对参数β和θ进行更新，即求取β和θ后验分布。本专利采用Bayesian方法来对β和θ的后验分布进行更新。
考虑到β和θ是相互独立的，因此β和θ的联合先验分布为p(β，θ)＝p(β)p(θ)。对于给定的β和θ，基于标准布朗运动{B(t)，t≥0}的独立增量性和马尔科夫性，
Figure BDA0003173392380000126
服从多元正态分布，其联合概率密度函数为：
Figure BDA0003173392380000127
然后，基于Bayesian定理，已知
Figure BDA0003173392380000128
前提下β和θ的联合后验分布为：
Figure BDA0003173392380000131
由于
Figure BDA0003173392380000132
和p(β，θ)是共轭的，从而使得后验估计
Figure BDA0003173392380000133
也是高斯的，即
Figure BDA0003173392380000134
其联合后验概率密度函数为：
Figure BDA0003173392380000135
其中，
Figure BDA0003173392380000136
式(18)中的参数值利用历史N个随机退化设备Box-Cox变换退化监测数据基于第一部分的第(4)步的未知参数极大似然估计结果作为
Figure BDA0003173392380000137
的值，即
Figure BDA0003173392380000138
步骤4，在所述在线服役设备经Box-Cox变换后的随机退化模型的基础上，基于首达时间，对在线服役设备在当前时刻tk的剩余寿命Lk进行描述；基于退化量初值β和漂移系数θ的联合后验概率密度函数，得到在线服役设备剩余使用寿命的概率密度函数和概率分布函数；随着所述在线服役设备经Box-Cox变换后的随机退化模型中参数集合
Figure BDA0003173392380000141
的更新，基于全概率公式得到在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解；将Lk的均值作为服役设备剩余使用寿命的预测值
Figure BDA00031733923800001411
在式(14)描述的线性Wiener随机退化过程基础上，基于首达时间的概念，在线服役随机退化设备在tk时刻的剩余使用寿命可由下式求得：
Figure BDA0003173392380000142
其中，
Figure BDA0003173392380000143
是经过Box-Cox变换后的失效阈值，定义为
Figure BDA0003173392380000144
ω为在线服役随机退化设备失效阈值，一般可根据设计指南、工业标准以及专家经验事先确定。
已知在线服役设备预测时刻为tk，同时给定θ，β和到tk时刻的经Box-Cox变换的退化监测数据集，利用标准布朗运动的性质，可将式(19)改写为：
Figure BDA0003173392380000145
其中，W(lk)＝B(tk+lk)-B(tk)。
易证明{W(lk)，lk≥0}仍然是一个标准布朗运动，因此{S(tk+lk，λ)，lk≥0}是一个带有线性漂移系数θlk的Wiener过程，其初始值为
Figure BDA0003173392380000146
首达失效阈值时间意义下带有线性漂移系数的Wiener过程描述的在线服役设备剩余寿命服从逆高斯分布，在已知β，θ和
Figure BDA0003173392380000147
的前提下，其剩余使用寿命的概率密度函数
Figure BDA0003173392380000148
和概率分布函数
Figure BDA0003173392380000149
分别为：
Figure BDA00031733923800001410
Figure BDA0003173392380000151
在得到
Figure BDA0003173392380000152
随着参数的更新
Figure BDA0003173392380000153
基于全概率公式可得在线服役设备tk剩余使用寿命的概率密度函数
Figure BDA0003173392380000154
和概率分布函数
Figure BDA0003173392380000155
封闭解析解为：
Figure BDA0003173392380000156
Figure BDA0003173392380000157
式(23)、(24)中的参数取值为式(18)经Bayesian参数在线更新得到的估计值。
进而，基于式(23)(24)可以得到在线服役设备剩余使用寿命的均值、方差、置信区间等估计结果，完成在线服役设备剩余使用寿命的自适应预测。
仿真实验
下面通过仿真数据处理结果进一步说明本发明的正确性和有效性。
仿真内容：本实施例以锂电池为研究对象，众所周知，锂电池是现代工业系统中一个重要的退化部件，且其退化数据表现出非常明显的非线性。具体的，本专利采用马里兰大学先进生命周期工程中心发布的锂电池退化数据，选取CS2-35、CS2-36、CS2-37和CS2-38共四组相同退化测试条件锂电池退化数据作为退化建模和剩余使用寿命预测验证数据，锂电池的退化通常用电池的容量来表示，图1给出了CS2-35、CS2-36、CS2-37和CS2-38退化趋势数据。
由图1可知，CS2-35、CS2-36、CS2-37和CS2-38四组退化轨迹均表现出明显的非线性。接下来，应用Box-Cox变换对原始具有明显非线性的退化数据进行处理，其中CS2-35、CS2-37和CS2-38三组数据被用来进行式(2)所示线性Wiener随机退化过程参数的极大似然估计，即，N＝3，求得参数估计值记为
Figure BDA0003173392380000161
其中
Figure BDA0003173392380000162
基于
Figure BDA0003173392380000163
经Box-Cox变换后的原始退化监测数据(如图1所示)线性化效果示于图2中。
在图2中，实线为经过Box-Cox变换后的锂电池退化轨迹，虚线为本发明提出的基于公式(2)所示线性Wiener随机过程且其参数由极大似然估计得到进而拟合出的线性退化轨迹，由图2可以看出，经Box-Cox变换后的退化轨迹能很好的跟踪线性拟合退化轨迹，说明本发明提出的Box-Cox变换可以有效的把非线性退化轨迹变换成近似线性退化轨迹。
进一步，以CS2-36号锂电池为在线服役随机退化设备，不失一般性，假设最后监测时刻分别为第100个和第400个循环周期，利用本发明提出的方法对实际非线性退化轨迹做Box-Cox变换，其变换后的退化轨迹和预测得到的退化轨迹如图3、图4所示，同时图3、图4还给出了不经Box-Cox变换的实际非线性退化轨迹以及针对该非线性退化轨迹采用线性Wiener建模得到的拟合退化轨迹图。
分析图3、图4可知，基于本发明提出的Box-Cox变换的实际退化路径能很好的跟踪预测退化趋势，而不经Box-Cox变换的非线性退化轨迹不能很好的匹配线性拟合退化趋势。进一步说明了本发明提出的Box-Cox变换在处理非线性退化数据进而提升线性预测趋势精度上的有效性。
最后，基于Box-Cox变换后的退化轨迹进行CS2-36锂电池的剩余使用寿命预测，锂电池的退化通常用电池的额定容量来表示，本专利定义当电池的额定容量低于45％时即认为电池失效，基于此，经Box-Cox变换后的退化轨迹对应的失效阈值为-0.215，对应于CS2-36的第772个循环周期，基于本专利提出的历史随机退化设备随机退化线性Wiener随机模型参数极大似然估计以及基于Bayesian算法的参数自适应更新，得到CS2-36的剩余使用寿命概率密度函数如图5所示。
分析图5可知，随着在线服役设备监测时刻的进行，其剩余使用寿命的概率密度函数被实时更新，并且剩余使用寿命预测均值逐渐接近于剩余使用寿命真值，即剩余使用寿命预测精度越来越大。
本发明基于Box-Cox变换提出了一种更一般性的非线性到线性的数据变换技术，以便用线性Wiener过程实现非线性退化数据建模，同时能得到剩余寿命概率密度函数精确解析解，提高预测精度。
虽然，本说明书中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (9)

1.一种基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，包括以下步骤：
步骤1，对于在线服役运行的非线性随机退化设备，收集在线服役设备从开始运行时刻t0到当前时刻tk的退化监测值x1，x2，…，xk，记为在线服役设备原始监测退化数据X1：k＝{x1，x2，…，xk}；使用Box-Cox变换对X1：k进行处理，得到经Box-Cox变换后的退化数据
Figure FDA0003173392370000011
其中，
Figure FDA0003173392370000012
为Box-Cox变换参数，其值通过历史随机退化设备的退化监测数据估计得到；
步骤2，采用线性Wiener随机过程对在线服役设备经Box-Cox变换后的退化数据进行建模，得到在线服役设备随机退化模型，模型参数中退化量初值β和漂移系数θ服从正态分布，其超参数取值通过历史随机退化设备的退化监测数据估计得到；
步骤3，采用Bayesian方法对所述在线服役设备随机退化模型中的退化量初值β和漂移系数θ的后验分布进行更新，得到退化量初值β和漂移系数θ的联合后验分布，从而得到退化量初值β和漂移系数θ的联合后验概率密度函数；
步骤4，在所述在线服役设备随机退化模型的基础上，基于首达时间，对在线服役设备在当前时刻tk的剩余寿命Lk进行描述；基于退化量初值β和漂移系数θ的联合后验概率密度函数，得到在线服役设备剩余使用寿命的概率密度函数和概率分布函数；基于全概率公式得到在线服役设备tk时刻剩余使用寿命的概率密度函数和概率分布函数的封闭解析解；将Lk的期望作为服役设备剩余使用寿命的预测值
Figure FDA0003173392370000013
2.根据权利要求1所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，所述使用Box-Cox变换对X1：k进行处理，具体为：
使用Box-Cox变换对{x(t)，t≥0}进行处理，对于每一个具体的时间t，有
Figure FDA0003173392370000021
其中，x(t)＝xt；
经过Box-Cox变换后的在线服役设备在tk时刻的退化监测量记为
Figure FDA0003173392370000022
即
Figure FDA0003173392370000023
则在线服役设备经Box-Cox变换的监测序列记为
Figure FDA0003173392370000024
3.根据权利要求2所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，所述在线服役设备随机退化模型为：
Figure FDA0003173392370000025
其中，σ为扩散系数，B(t)为是标准布朗运动，且有σB(t)～N(0，σ2t)。
4.根据权利要求3所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，所述采用Bayesian方法对所述在线服役设备经Box-Cox变换后的随机退化模型中的退化量初值β和漂移系数θ的后验分布进行更新，具体为：
3.1，由于β和θ是相互独立的，因此β和θ的联合先验分布为p(β，θ)＝p(β)p(θ)；对于给定的β和θ，基于标准布朗运动{B(t)，t≥0}的独立增量性和马尔科夫性，
Figure FDA0003173392370000026
服从多元正态分布，其联合概率密度函数为：
Figure FDA0003173392370000027
3.2，基于Bayesian定理，已知
Figure FDA0003173392370000028
前提下β和θ的联合后验分布为：
Figure FDA0003173392370000031
3.3，由于
Figure FDA0003173392370000032
和p(β，θ)是共轭的，从而使得后验估计
Figure FDA0003173392370000033
也是高斯的，即
Figure FDA0003173392370000034
则其联合后验概率密度函数为：
Figure FDA0003173392370000035
其中，
Figure FDA0003173392370000036
上式中的μβ，
Figure FDA0003173392370000037
μθ，
Figure FDA0003173392370000038
参数值通过历史随机退化设备的退化监测数据估计得到，即
Figure FDA0003173392370000039
5.根据权利要求4所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，所述基于首达时间，对在线服役设备在当前时刻tk的剩余寿命Lk进行描述为：
Figure FDA00031733923700000310
其中，
Figure FDA00031733923700000311
是在线服役设备经过Box-Cox变换后的失效阈值，定义为
Figure FDA00031733923700000312
ω为在线服役设备的失效阈值，其取值根据设计指南、工业标准或专家经验事先确定；lk为tk时刻的剩余寿命对应的随机变量；inf表示下确界。
6.根据权利要求5所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，所述在线服役设备剩余使用寿命的概率密度函数和概率分布函数的获取过程为：
已知在线服役设备预测时刻为tk，同时给定θ、β和到tk时刻的经Box-Cox变换的退化监测数据集，利用标准布朗运动的性质，将剩余寿命Lk在首达时间下的描述改写为：
Figure FDA0003173392370000041
其中，W(lk)＝B(tk+lk)-B(tk)；
由于{W(lk)，lk≥0}仍然是一个标准布朗运动，因此{S(tk+lk，λ)，lk≥0}是带有线性漂移系数θlk的Wiener过程，其初始值为
Figure FDA0003173392370000042
首达失效阈值时间意义下带有线性漂移系数的Wiener过程描述的在线服役设备剩余寿命服从逆高斯分布，在已知β，θ和
Figure FDA0003173392370000043
的前提下，其剩余使用寿命的概率密度函数
Figure FDA0003173392370000044
和概率分布函数
Figure FDA0003173392370000045
分别为：
Figure FDA0003173392370000046
Figure FDA0003173392370000047
随着所述在线服役设备经Box-Cox变换后的随机退化模型中参数集合
Figure FDA0003173392370000048
的更新，基于全概率公式得到在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解为：
Figure FDA0003173392370000051
式中，
Figure FDA0003173392370000052
Figure FDA0003173392370000053
中各参数的取值为式(1)经Bayesian参数在线更新得到的估计值；
最后，基于在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解公式得到在线服役设备剩余使用寿命的均值估计结果。
7.根据权利要求1所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，还包括：基于在线服役设备tk剩余使用寿命的概率密度函数和概率分布函数的封闭解析解公式计算Lk的方差和置信区间，以量化剩余寿命预测结果的不确定性。
8.根据权利要求1-7任一项所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，所述通过历史随机退化设备的Box-Cox变换退化监测数据估计得到参数
Figure FDA0003173392370000054
的具体过程为：
(a)设有N个同类历史随机退化设备，每个退化设备有Mj，j＝1，2，...，N个监测时间点个数，且将第j个随机退化设备对应监测时刻表示为
Figure FDA0003173392370000055
则第j个随机退化设备在
Figure FDA0003173392370000056
时刻的退化监测量记为
Figure FDA0003173392370000057
Figure FDA0003173392370000058
是第j个随机退化设备的第qj个测量时刻，因此qj＝1，2，…，Mj；
Figure FDA0003173392370000059
记第j个随机退化设备的原始退化监测数据为
Figure FDA00031733923700000510
则N个同类随机退化设备原始退化监测数据记为
Figure FDA00031733923700000511
(b)使用Box-Cox变换对{Xj(t)，t≥0}进行处理，对于每一个具体的时间t，有
Figure FDA0003173392370000061
其中，λ是Box-Cox变换参数；
经过Box-Cox变换后的第j个随机退化设备在
Figure FDA0003173392370000062
时刻的退化监测量记为
Figure FDA0003173392370000063
即
Figure FDA0003173392370000064
则第j个历史随机退化设备经Box-Cox变换的监测序列记为
Figure FDA0003173392370000065
同理，N个同类历史随机退化设备经Box-Cox变换的监测序列记为
Figure FDA0003173392370000066
(c)采用{S(t，λ)，t≥0}表示随机退化设备经过Box-Cox变换后的退化过程，基于线性Wiener建模该退化过程得：
S(t，λ)＝β+θt+σB(t)
其中，θ是漂移系数，σ是扩散系数；β是退化量初值；{B(t)，t≥0}是标准布朗运动，且有σB(t)～N(0，σ2t)；
令β和θ为随机变量，且分别服从正态分布，即
Figure FDA0003173392370000067
基于以上描述，记未知参数集合为
Figure FDA0003173392370000068
(d)设N个历史随机退化设备的退化过程是相互独立的，则该N个历史随机退化设备经Box-Cox变换的历史监测数据
Figure FDA0003173392370000069
的对数似然函数为：
Figure FDA00031733923700000610
其中，C是常数，第j个历史随机退化设备经Box-Cox变换的监测序列
Figure FDA00031733923700000611
服从多元正态分布，即
Figure FDA00031733923700000612
其均值为
Figure FDA00031733923700000613
方差为：
Figure FDA00031733923700000614
Figure FDA00031733923700000615
采用极大似然估计求解所述对数似然函数
Figure FDA00031733923700000616
即可得到参数
Figure FDA0003173392370000071
的估计值。
9.根据权利要求8所述的基于Box-Cox变换的非线性设备剩余寿命预测方法，其特征在于，所述采用极大似然估计求解所述对数似然函数
Figure FDA0003173392370000072
具体求解过程为：
对参数进行重构，有
Figure FDA0003173392370000073
Figure FDA0003173392370000074
Figure FDA0003173392370000075
将参数重构后的未知参数集合记为
Figure FDA0003173392370000076
其联合对数似然函数为：
Figure FDA0003173392370000077
基于式(3)分别对μβ，μθ和
Figure FDA0003173392370000078
求偏导数有：
Figure FDA0003173392370000079
Figure FDA00031733923700000710
Figure FDA00031733923700000711
为求得
Figure FDA00031733923700000712
的极大似然估计值，首先固定
Figure FDA00031733923700000713
λ值，令式(4)、(5)、(6)分别等于零，解偏微分方程组，得到μβ，μθ，
Figure FDA00031733923700000714
的极大似然估计值：
Figure FDA0003173392370000081
其中，
Figure FDA0003173392370000082
接下来，将基于式(7)、式(8)得到的μβ，μθ，
Figure FDA0003173392370000083
的极大似然估计值
Figure FDA0003173392370000084
和
Figure FDA0003173392370000085
代入式(2)，得到
Figure FDA00031733923700000812
λ的边缘似然函数为：
Figure FDA0003173392370000087
同理，将式(9)的边缘似然函数分别对
Figure FDA0003173392370000088
λ求偏导数并令其等于零，解方程即可得到
Figure FDA0003173392370000089
λ的极大似然估计值；
然后再基于
Figure FDA00031733923700000810
变换，得到未知参数集合估计值
Figure FDA00031733923700000811
CN202110825429.4A 2021-07-21 2021-07-21 基于Box-Cox变换的非线性设备剩余寿命预测方法 Active CN113642158B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110825429.4A CN113642158B (zh) 2021-07-21 2021-07-21 基于Box-Cox变换的非线性设备剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110825429.4A CN113642158B (zh) 2021-07-21 2021-07-21 基于Box-Cox变换的非线性设备剩余寿命预测方法

Publications (2)

Publication Number Publication Date
CN113642158A CN113642158A (zh) 2021-11-12
CN113642158B true CN113642158B (zh) 2022-09-30

Family

ID=78417910

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110825429.4A Active CN113642158B (zh) 2021-07-21 2021-07-21 基于Box-Cox变换的非线性设备剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN113642158B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112883550A (zh) * 2021-01-19 2021-06-01 中国人民解放军火箭军工程大学 一种考虑多重不确定性的退化设备剩余寿命预测方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112883550A (zh) * 2021-01-19 2021-06-01 中国人民解放军火箭军工程大学 一种考虑多重不确定性的退化设备剩余寿命预测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Bayesian更新与EM算法协作下退化数据驱动的剩余寿命估计方法;司小胜等;《模式识别与人工智能》;20130415(第04期);全文 *
考虑随机失效阈值的设备剩余寿命在线预测;王泽洲等;《系统工程与电子技术》;20190319(第05期);全文 *

Also Published As

Publication number Publication date
CN113642158A (zh) 2021-11-12

Similar Documents

Publication Publication Date Title
CN109472110B (zh) 一种基于lstm网络和arima模型的航空发动机剩余使用寿命预测方法
Le Son et al. Remaining useful lifetime estimation and noisy gamma deterioration process
CN110851980B (zh) 一种设备剩余寿命预测方法及系统
CN112926273B (zh) 一种多元退化设备剩余寿命预测方法
Huang et al. Remaining useful life prediction for an adaptive skew-Wiener process model
CN107765347B (zh) 一种高斯过程回归和粒子滤波的短期风速预测方法
Pang et al. A Bayesian inference for remaining useful life estimation by fusing accelerated degradation data and condition monitoring data
Hu et al. Online performance assessment method for a model-based prognostic approach
Djeziri et al. Data-driven approach augmented in simulation for robust fault prognosis
Peng et al. The transformed inverse Gaussian process as an age-and state-dependent degradation model
CN110955963A (zh) 一种航空电缆剩余寿命预测方法
CN111458661A (zh) 一种配电网线变关系诊断方法、装置及系统
CN116306798A (zh) 一种超短时风速预测方法及系统
CN115128427B (zh) Mos器件寿命预测方法、装置、电子设备、介质及程序产品
CN111523727B (zh) 基于不确定过程的考虑恢复效应的电池剩余寿命预测方法
Ge et al. An improved PF remaining useful life prediction method based on quantum genetics and LSTM
Gao et al. A novel remaining useful life prediction method for capacity diving lithium-ion batteries
CN113642158B (zh) 基于Box-Cox变换的非线性设备剩余寿命预测方法
Lin An integrated procedure for bayesian reliability inference using MCMC
Song et al. A sliding sequence importance resample filtering method for rolling bearings remaining useful life prediction based on two Wiener-process models
Zhao et al. Rolling bearing remaining useful life prediction based on wiener process
CN112541243B (zh) 一种数据驱动的系统剩余使用寿命预测方法
CN115577854A (zh) 一种基于eemd-rbf组合的分位数回归风速区间预测方法
CN110895624A (zh) 基于最大熵谱估计的加速贮存与自然贮存退化数据一致性检验法
CN109992875B (zh) 一种切换设备剩余寿命的确定方法及系统

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