WO2019174142A1 - 一种多模式的退化过程建模及剩余寿命预测方法 - Google Patents
一种多模式的退化过程建模及剩余寿命预测方法 Download PDFInfo
- Publication number
- WO2019174142A1 WO2019174142A1 PCT/CN2018/090928 CN2018090928W WO2019174142A1 WO 2019174142 A1 WO2019174142 A1 WO 2019174142A1 CN 2018090928 W CN2018090928 W CN 2018090928W WO 2019174142 A1 WO2019174142 A1 WO 2019174142A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- mode
- degradation
- time
- change point
- slope
- Prior art date
Links
Images
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B23/00—Testing or monitoring of control systems or parts thereof
- G05B23/02—Electric testing or monitoring
- G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
- G05B23/0259—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterized by the response to fault detection
- G05B23/0283—Predictive maintenance, e.g. involving the monitoring of a system and, based on the monitoring results, taking decisions on the maintenance schedule of the monitored system; Estimating remaining useful life [RUL]
-
- 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
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B23/00—Testing or monitoring of control systems or parts thereof
- G05B23/02—Electric testing or monitoring
- G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
- G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
- G05B23/0243—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults model based detection method, e.g. first-principles knowledge model
- G05B23/0254—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults model based detection method, e.g. first-principles knowledge model based on a quantitative model, e.g. mathematical relationships between inputs and outputs; functions: observer, Kalman filter, residual calculation, Neural Networks
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B23/00—Testing or monitoring of control systems or parts thereof
- G05B23/02—Electric testing or monitoring
- G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
- G05B23/0259—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterized by the response to fault detection
- G05B23/0275—Fault isolation and identification, e.g. classify fault; estimate cause or root of failure
- G05B23/0281—Quantitative, e.g. mathematical distance; Clustering; Neural networks; Statistical analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
- G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
-
- 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/04—Ageing analysis or optimisation against ageing
Definitions
- the invention belongs to the technical field of health management, and particularly relates to a multi-mode degradation process modeling and residual life prediction method.
- the remaining life prediction of industrial equipment can provide effective information for the equipment to develop maintenance strategies and production decisions, thereby reducing the losses caused by equipment failure and ensuring the safety and reliability of the system.
- Degradation process modeling is a key step in residual life prediction.
- the degradation model should be as close as possible to the actual degradation.
- Most existing methods currently assume that the equipment does not have large operating conditions and maintenance operations throughout its life cycle. However, in actual industrial production, there may be multiple operating conditions during the operation of the equipment, and the equipment will be regularly or inspected for maintenance or maintenance activities. These activities can affect the degradation process of the device and present multiple patterns during the degradation process. At present, there is no residual life prediction method that can automatically identify multiple retreat modes.
- the multi-mode degradation process modeling and residual life prediction have the following difficulties: First, since the degeneration mode switching has no tags, it is necessary to identify the number of degradation modes and the degradation models in each mode based on historical data. Second, since the degradation model contains fractal Brownian motion, it is not a Markov process or a half-turn, so it is difficult to find the first-time distribution of the parsing. Third, the future degradation mode switching situation is unknown, and it is necessary to consider the possible degeneration mode switching in the future when predicting the remaining life.
- the present invention proposes a multi-mode degradation process modeling and residual life prediction method, which is reasonable in design, overcomes the deficiencies of the prior art, and has good effects.
- a multi-mode degradation process modeling and residual life prediction method includes the following steps:
- Step 1 Collect device degradation data x 0 , x 1 , x 2 , . . . , x k at equally spaced sampling instants t 0 , t 1 , t 2 , . . . , t k , respectively , where the sampling interval is ⁇ , k is the number of samples;
- Step 2 According to the change point detection method, the slope change point of the historical degradation process is detected, and is recorded as ⁇ 1 , ⁇ 2 , ... ⁇ j , ⁇ j+1 ...;
- Step 3 Obtain the degenerate segment by using the points ⁇ j and ⁇ j+1 obtained in step 2 as the end point, and calculate the slope of the degenerate segment according to the following formula And this slope The characteristic value of the jth degenerate segment;
- d c is the truncation distance
- d ji
- ⁇ ( ⁇ ) is defined as follows:
- Step 4 If ⁇ j and ⁇ j are respectively greater than the corresponding thresholds, the line segment obtained by using the slope change points ⁇ j and ⁇ j+1 as the end point is a cluster center, and according to this method, the number of cluster centers to be obtained is obtained.
- N that is, the line segment obtained by dividing the slope change point can be clustered into N categories according to the slope thereof, and the degradation mode of the sampling time u is recorded as ⁇ (u), and And record the time point of the degradation mode change as c 1 , c 2 ,...;
- Step 5 Establish a degradation model:
- X(0) is the initial value of the degradation process and ⁇ [ ⁇ (u)] is the drift term coefficient, assuming ⁇ H is the diffusion term coefficient, B H (t) is the standard fractal Brownian motion, and the degradation mode ⁇ (u) is a continuous-time Markov chain whose transfer rate matrix is Q;
- Step 6 According to Estimating q j and q ij in the transfer rate matrix Q of continuous-time Markov chains:
- j is the m t k until the time reaches the degraded mode and the number of the j-th mode
- m ij is the transition from the degraded mode times t i-th mode to the j-th travel mode before the time k, The dwell time for the ith time to reach the jth mode
- Step 7 Estimate the Hurst exponent H of the degradation process:
- ⁇ 1 , ⁇ 2 , . . . , ⁇ p are wavelet decomposition high-pass filter coefficients based on the Symles wavelet function, p is the wavelet function vanishing moment order; E( ⁇ ) is a mathematical expectation;
- Step 8 Estimate the estimated value ⁇ i of the drift term coefficient ⁇ [ ⁇ (u)] in the i-th segment degradation mode, respectively:
- I i is a column vector of c i+1 -c i dimensions, each element is 1, Is a c i+1 -c i -dimensional covariance matrix whose elements in the i-th row and j-th column are
- Step 9 Estimate the expectation of the drift term coefficient ⁇ [ ⁇ (u)] for each degradation mode And variance
- ⁇ j,i is an estimate of the drift term coefficient obtained when the i thth reaches the jth degradation mode
- Step 10 Estimate the diffusion term coefficient ⁇ H :
- Step 11 Order Using the Monte Carlo method, the numerical distribution of ⁇ [ ⁇ (t k ), l k ] is obtained.
- Step 12 For a given failure threshold ⁇ , the approximate distribution of the first arrival time of the degradation process is:
- ⁇ min min ⁇ (1), ⁇ (2),..., ⁇ (N) ⁇
- ⁇ max max ⁇ (1), ⁇ (2),..., ⁇ (N) ⁇
- ⁇ ( ⁇ ) is a Gamma function.
- step 2 the method specifically includes the following steps:
- Step 2.3 Calculate the change point detection index ⁇ (i) of the i-th point:
- step 11 the method specifically includes the following steps:
- Step 11.2 Generate n uniformly distributed random numbers r j obeying [0, 1];
- Step 11.4 Count the total length of time each Monte Carlo sample sequence stays in each mode during the (t k , t k + l k ) time interval, recorded as
- Step 11.5 Calculate the numerical distribution of ⁇ [ ⁇ (t k ), l k ]
- the invention provides a method for identifying a degradation mode from historical degradation data, and the method is applicable to a system with state switching, environmental change or maintenance operation, which is closer to the actual situation than the prior art; and different degradation is obtained according to historical data identification.
- a degradation model with state-switching based on fractal Brownian motion is established.
- the degeneration mode is described by a continuous-time Markov chain.
- the maximum likelihood method is used to estimate the state transition rate of the Markov chain and the drift coefficient in the degradation model.
- the present invention in order to obtain the first time of the degradation process, the present invention first obtains the numerical distribution of the diffusion term in the future time period by the Monte Carlo method, and then uses a time-space transformation to first the degradation process based on the fractal Brownian motion. The time is converted to the time when the standard Brownian motion first reaches the time-varying threshold with uncertainty, and then the distribution of the remaining life of the degradation process is obtained.
- FIG. 1 is a flow chart of a degradation process modeling and remaining life prediction method of the present invention
- FIG. 2 is a flow chart of a method for detecting a change point in the present invention
- Figure 3 is a flow chart of the Monte Carlo method of the present invention.
- Figure 4 is a schematic view showing the temperature degradation curve of the blast furnace wall in the example of the present invention.
- FIG. 5 is a schematic view showing a change point detection result of 1200 days before the temperature degradation curve of the blast furnace wall in the example of the present invention
- FIG. 6 is a schematic diagram showing the degradation mode identification result of the first 1200 days before the temperature degradation curve of the blast furnace wall in the example of the present invention
- Fig. 7 is a view showing the results of predicting the remaining life of the blast furnace wall in the example of the present invention.
- a multi-mode degradation process modeling and residual life prediction method including the following steps:
- Step 1 Collect device degradation data x 0 , x 1 , x 2 , . . . , x k at equally spaced sampling instants t 0 , t 1 , t 2 , . . . , t k , respectively , where the sampling interval is ⁇ , k is the number of samples;
- Step 2 According to the change point detection method (the flow is shown in FIG. 2), the slope change point of the historical degradation process is detected, and is recorded as ⁇ 1 , ⁇ 2 , .
- Step 3 Obtain the degenerate segment by using the points ⁇ j and ⁇ j+1 obtained in step 2 as the end point, and calculate the slope of the degenerate segment according to the following formula And this slope The characteristic value of the jth degenerate segment;
- d c is the truncation distance
- d ji
- ⁇ ( ⁇ ) is defined as follows:
- Step 4 If ⁇ j and ⁇ j are respectively greater than the corresponding thresholds, the line segment obtained by using the slope change points ⁇ j and ⁇ j+1 as the end point is a cluster center, and according to this method, the number of cluster centers to be obtained is obtained.
- N that is, the line segment obtained by dividing the slope change point can be clustered into N categories according to the slope thereof, and the degradation mode of the sampling time u is recorded as ⁇ (u), and And record the time point of the degradation mode change as c 1 , c 2 ,...;
- Step 5 Establish a degradation model:
- X(0) is the initial value of the degradation process and ⁇ [ ⁇ (u)] is the drift term coefficient, assuming ⁇ H is the diffusion term coefficient, B H (t) is the standard fractal Brownian motion, and the degradation mode ⁇ (u) is a continuous-time Markov chain whose transfer rate matrix is Q;
- Step 6 According to Estimating q j and q ij in the transfer rate matrix Q of continuous-time Markov chains:
- j is the m t k until the time reaches the degraded mode and the number of the j-th mode
- m ij is the transition from the degraded mode times t i-th mode to the j-th travel mode before the time k, The dwell time for the ith time to reach the jth mode
- Step 7 Estimate the Hurst exponent H of the degradation process:
- ⁇ 1 , ⁇ 2 , . . . , ⁇ p are wavelet decomposition high-pass filter coefficients based on the Symles wavelet function, p is the wavelet function vanishing moment order; E( ⁇ ) is a mathematical expectation;
- Step 8 Estimate the estimated value ⁇ i of the drift term coefficient ⁇ [ ⁇ (u)] in the i-th segment degradation mode, respectively:
- I i is a column vector of c i+1 -c i dimensions, each element is 1, Is a c i+1 -c i -dimensional covariance matrix whose elements in the i-th row and j-th column are
- Step 9 Estimate the expectation of the drift term coefficient ⁇ [ ⁇ (u)] for each degradation mode And variance
- ⁇ j,i is an estimate of the drift term coefficient obtained when the i thth reaches the jth degradation mode
- Step 10 Estimate the diffusion term coefficient ⁇ H :
- Step 11 Order Using the Monte Carlo method, the numerical distribution of ⁇ [ ⁇ (t k ), l k ] is obtained.
- Step 12 For a given failure threshold ⁇ , the approximate distribution of the first arrival time of the degradation process is:
- ⁇ min min ⁇ (1), ⁇ (2),..., ⁇ (N) ⁇
- ⁇ max max ⁇ (1), ⁇ (2),..., ⁇ (N) ⁇
- ⁇ ( ⁇ ) is a Gamma function.
- step 2 the following steps are specifically included (the process is as shown in FIG. 2):
- Step 2.3 Calculate the change point detection index ⁇ (i) of the i-th point:
- step 11 the following steps are specifically included (the flow is shown in FIG. 3):
- Step 11.2 Generate n uniformly distributed random numbers r j obeying [0, 1];
- Step 11.4 Count the total length of time each Monte Carlo sample sequence stays in each mode during the (t k , t k + l k ) time interval, recorded as
- Step 11.5 Calculate the numerical distribution of ⁇ [ ⁇ (t k ), l k ]
- the time to reach the temperature threshold of 530 ° C for the first time is the 1456 day.
- the remaining life predictions were performed at 1300, 1310, 1320, 1330, 1340, 1350, 1360, 1370, 1380, 1390, and 1400 days, respectively, and the results obtained are shown in FIG. It can be seen that the true value of the remaining life is located at a higher probability density of the prediction result, which verifies the effectiveness of the proposed method.
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Automation & Control Theory (AREA)
- Artificial Intelligence (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Evolutionary Biology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Probability & Statistics with Applications (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Computational Mathematics (AREA)
- Algebra (AREA)
- Geometry (AREA)
- Computer Hardware Design (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
一种多模式的退化过程建模及剩余寿命预测方法,属于健康管理的技术领域,包括以下步骤:首先,收集等间隔采样的退化数据;对退化数据进行变点检测;对变点分割得到的退化段以退化速率为特征进行聚类;建立包含模式切换的退化模型,模式切换通过一个连续时间Markov链描述;采用基于二次变分的方法估计退化过程的Hurst指数;利用极大似然法分别估计Markov链的状态转移率矩阵和各模式下的漂移项系数,以及扩散项系数;利用蒙特卡洛算法得到未来一段时间在状态切换的影响下漂移项服从的分布;然后在给定的阈值下,得到剩余寿命的分布。
Description
本发明属于健康管理的技术领域,具体涉及一种多模式的退化过程建模及剩余寿命预测方法。
工业设备的剩余寿命预测能够为设备制定维护策略和生产决策提供有效的信息,从而减少由于设备故障带来的损失,确保系统的安全性和可靠性的。
退化过程建模是剩余寿命预测的关键步骤。为了得到准确的剩余寿命预测结果,退化模型应该尽可能地符合实际退化情况。目前大多数现有的方法均假设设备再其整个寿命周期内不存在较大的工况变化及维护操作。然而,在实际的工业生产中,设备的运行过程可能存在多种工况,同时设备会进行定期或视情的检修或维护活动。这些活动会影响设备的退化过程,在退化过程中呈现出多种模式。目前,尚未出现能够自动辨识多种退模式的剩余寿命预测方法。
多模式的退化过程建模及剩余寿命预测主要有以下难点:第一,由于退化模式切换没有标签,需要根据历史数据辨识退化模式的数量和各个模式下的退化模型。第二,由于退化模型中包含分形布朗运动,它不是Markov过程也不是半鞅,因此难以求得解析的首达时间分布。第三,未来的退化模式切换情况未知,在预测剩余寿命时需要考虑未来可能发生的退化模式切换。
发明内容
针对现有技术中存在的上述技术问题,本发明提出了一种多模式的退化过程建模及剩余寿命预测方法,设计合理,克服了现有技术的不足,具有良好的效果。
为了实现上述目的,本发明采用如下技术方案:
一种多模式的退化过程建模及剩余寿命预测方法,包括如下步骤:
步骤1:分别在等间隔采样时刻t
0,t
1,t
2,...,t
k收集设备退化数据x
0,x
1,x
2,...,x
k,其中,采样间隔为τ,k为采样个数;
步骤2:根据变点检测方法,检测历史退化过程的斜率变化点,记为γ
1,γ
2,...γ
j,γ
j+1...;
计算每个退化段的特征值的局部密度ρ
j,并计算比局部密度更大的特征值的最小距离δ
j,其中局部密度ρ
j根据式(1)计算:
其中,d
c为截断距离,d
ji=|η
i-η
j|,函数χ(·)定义如下:
根据下式(3)计算最小距离δ
j:
步骤4:若ρ
j和δ
j分别大于相应的阈值,则以斜率变化点γ
j和γ
j+1为端点得到的线段是一个聚类中心,按照这种方法,将得到的聚类中心数量记为N,即由斜率变化点分割得到的线段按照其斜率可以聚为N类,将采样时刻u的退化模式记为Φ(u),并令
并将退化模式变化的时间点记为c
1,c
2,...;
步骤5:建立退化模型:
步骤7:估计退化过程的Hurst指数H:
其中,θ
1,θ
2,...,θ
p为基于Symlets小波函数的小波分解高通滤波器系数,p为该小波函数消失矩阶数;E(·)为数学期望;
步骤8:分别估计第i段退化模式下的漂移项系数λ[Φ(u)]的估计值λ
i:
其中,λ
j,i是第i次到达第j个退化模式时得到的漂移项系数估计值;
步骤10:估计扩散项系数σ
H:
步骤12:对于给定的失效阈值ω,退化过程首达时间的近似分布为:
其中,Γ(·)为Gamma函数。
优选地,在步骤2中,具体包括如下步骤:
步骤2.1:初始化参数,令γ
1=0,i=1,i
γ=1,选择两个变点间的最小间隔mτ,以及变点检测的阈值ω
β;
步骤2.2:计算从上一个变点到第i个点的退化段的斜率,如果i-i
γ>m,根据式(14)计算当前段的斜率η
i,并令i=i+1:
步骤2.3:计算第i个点的变点检测指标β(i):
步骤2.5:如果i≤k,则令i=i+1,然后执行步骤2.2。
优选地,在步骤11中,具体包括如下步骤:
步骤11.2:生成n个服从[0,1]上的均匀分布的随机数r
j;
步骤11.3:对于第j个蒙特卡洛样本序列,令i=i+1,v
i+1,j=s,其中,s满足
为经过时间τ退化模式由第v
i,j个模式变换到第w个模式的概率,如果i<l
k/τ,则返回步骤11.2,否则执行步骤11.4;
本发明所带来的有益技术效果:
本发明提出的一种从历史退化数据中辨识退化模式的方法,该方法适用于存在状态切换、环境变化或维护操作的系统,比现有技术更加贴近实际情况;根据历史数据辨识得到不同的退化模式,建立了基于分形布朗运动的含有状态切换的退化模型,退化模式的切换通过一个连续时间Markov链描述;利用极大似然方法分别估计Markov链的状态转移率以及退化模型中的漂移系数及扩散系数;为了求得退化过程的首达时间,本发明首先通过蒙特卡洛方法得到未来时间段扩散项的的数值分布,然后通过一个时间-空间变换,将基于分形布朗运动的退化过程的首达时间转换为标准布朗运动首达含有不确定性的时变阈值的时间,进而得到了退 化过程剩余寿命的分布。
图1是本发明退化过程建模及剩余寿命预测方法的流程图;
图2是本发明中变点检测方法的流程图;
图3是本发明中蒙特卡洛方法的流程图;
图4是本发明示例中,高炉炉壁温度退化曲线示意图;
图5是本发明示例中,高炉炉壁温度退化曲线前1200天的变点检测结果示意图;
图6是本发明示例中,高炉炉壁温度退化曲线前1200天的退化模式辨识结果示意图;
图7是本发明示例中,高炉炉壁剩余寿命预测结果示意图。
下面结合附图以及具体实施方式对本发明作进一步详细说明:
一种多模式的退化过程建模及剩余寿命预测方法,其流程如图1所示,包括如下步骤:
步骤1:分别在等间隔采样时刻t
0,t
1,t
2,...,t
k收集设备退化数据x
0,x
1,x
2,...,x
k,其中,采样间隔为τ,k为采样个数;
步骤2:根据变点检测方法(其流程如图2所示),检测历史退化过程的斜率变化点,记为γ
1,γ
2,...;
计算每个退化段的特征值的局部密度ρ
j,并计算比局部密度更大的特征值的最小距离δ
j,其中局部密度ρ
j根据式(1)计算:
其中,d
c为截断距离,d
ji=|η
i-η
j|,函数χ(·)定义如下:
根据下式(3)计算最小距离δ
j:
步骤4:若ρ
j和δ
j分别大于相应的阈值,则以斜率变化点γ
j和γ
j+1为端点得到的线段是一个聚类中心,按照这种方法,将得到的聚类中心数量记为N,即由斜率变化点分割得到的线段按照其斜率可以聚为N类,将采样时刻u的退化模式记为Φ(u),并令
并将退化模式变化的时间点记为c
1,c
2,...;
步骤5:建立退化模型:
步骤7:估计退化过程的Hurst指数H:
其中,θ
1,θ
2,...,θ
p为基于Symlets小波函数的小波分解高通滤波器系数,p为该小波函数消失矩阶数;E(·)为数学期望;
步骤8:分别估计第i段退化模式下的漂移项系数λ[Φ(u)]的估计值λ
i:
其中,λ
j,i是第i次到达第j个退化模式时得到的漂移项系数估计值;
步骤10:估计扩散项系数σ
H:
步骤12:对于给定的失效阈值ω,退化过程首达时间的近似分布为:
其中,Γ(·)为Gamma函数。
在步骤2中,具体包括如下步骤(其流程如图2所示):
步骤2.1:初始化参数,令γ
1=0,i=1,i
γ=1,选择两个变点间的最小间隔mτ,以及变点检测的阈值ω
β;
步骤2.2:计算从上一个变点到第i个点的退化段的斜率,如果i-i
γ>m,根据式(14)计算当前段的斜率η
i,并令i=i+1:
步骤2.3:计算第i个点的变点检测指标β(i):
步骤2.5:如果i≤k,则令i=i+1,然后执行步骤2.2。
在步骤11中,具体包括如下步骤(其流程如图3所示):
步骤11.2:生成n个服从[0,1]上的均匀分布的随机数r
j;
步骤11.3:对于第j个蒙特卡洛样本序列,令i=i+1,v
i+1,j=s,其中,s满足
为经过时间τ退化模式由第v
i,j个模式变换到第w个模式的概率,如果i<l
k/τ,则返回步骤11.2,否则执行步骤11.4;
本示例基于MATLAB工具,利用柳州钢铁2号高炉的数据对本发明进行说明,结合附图展示本发明的效果。
(1)收集高炉炉壁的温度传感器数据,选择位于高8.2米的热电偶,连续采样1506天的温度数据(取每天的温度平均值),记为x
0,x
1,x
2,...,x
1506,退化数据如图4所示;
(2)利用变点检测算法,检测t
k时刻之前的变点,记为γ
1,γ
2,...,变点检测结果如图5所示(以t
k=1300为例);
(3)利用式(14),分别计算由变点分割得到的线段斜率。
(4)根据式(1)、式(2)和式(3),分别计算每个线段的局部密度ρ
j和与比其密度更大的线段的最小距离δ
j;
(5)选择ρ
j>2且δ
j>26的线段作为聚类中心,共得到3个聚类中心,即N=3;
(6)分别计算距离每个线段最近的聚类中心,完成线段的聚类,得到的聚类结果如图6所示;
(7)利用式(5)和式(6)求得Markov链状态转移率矩阵估计值;
(8)利用式(7)求得退化过程的Hurst指数估计值H=0.8959;
(9)利用式(8)求得每段退化路径的漂移项系数估计值;
(10)利用式(9)和式(10)分别求得每种退化模式下漂移项系数的期望和方差;
(11)利用式(11)求得扩散项系数估计值σ
H=1.0599;
(10)假设阈值为530℃,根据式(12)得到第k个采样时刻剩余寿命分布的预测结果。
根据高炉炉壁的温度退化曲线可知,首次到达温度阈值530℃的时间为第1456天。分别在第1300、1310、1320、1330、1340、1350、1360、1370、1380、1390、1400天进行剩余寿命预测,得到的结果如图7所示。可以看出,剩余寿命真值均位于预测结果概率密度较高处,验证了本发明提出方法的有效性。
当然,上述说明并非是对本发明的限制,本发明也并不仅限于上述举例,本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换,也应属于本发明的保护范围。
Claims (3)
- 一种多模式的退化过程建模及剩余寿命预测方法,其特征在于:包括如下步骤:步骤1:分别在等间隔采样时刻t 0,t 1,t 2,...,t k收集设备退化数据x 0,x 1,x 2,...,x k,其中,采样间隔为τ,k为采样个数;步骤2:根据变点检测方法,检测历史退化过程的斜率变化点,记为γ 1,γ 2,...γ j,γ j+1...;计算每个退化段的特征值的局部密度ρ j,并计算比局部密度更大的特征值的最小距离δ j,其中局部密度ρ j根据式(1)计算:其中,d c为截断距离,d ji=|η i-η j|,函数χ(·)定义如下:根据下式(3)计算最小距离δ j:步骤4:若ρ j和δ j分别大于相应的阈值,则以斜率变化点γ j和γ j+1为端点得到的线段是一个聚类中心,按照这种方法,将得到的聚类中心数量记为N,即由斜率变化点分割得到的线段按照其斜率可以聚为N类,将采样时刻u的退化模式记为Φ(u),并令 并将退化模式变化的时间点记为c 1,c 2,...;步骤5:建立退化模型:步骤7:估计退化过程的Hurst指数H:其中,θ 1,θ 2,...,θ p为基于Symlets小波函数的小波分解高通滤波器系数,p为该小波函数消失矩阶数;E(·)为数学期望;步骤8:分别估计第i段退化模式下的漂移项系数λ[Φ(u)]的估计值λ i:其中,λ j,i是第i次到达第j个退化模式时得到的漂移项系数估计值;步骤10:估计扩散项系数σ H:步骤12:对于给定的失效阈值ω,退化过程首达时间的近似分布为:其中,Γ(·)为Gamma函数。
- 根据权利要求1所述的多模式的退化过程建模及剩余寿命预测方法,其特征在于:在步骤11中,具体包括如下步骤:步骤11.2:生成n个服从[0,1]上的均匀分布的随机数r j;步骤11.3:对于第j个蒙特卡洛样本序列,令i=i+1,v i+1,j=s,其中,s满足 为经过时间τ退化模式由第v i,j个模式变换到第w个模式的概率,如果i<l k/τ,则返回步骤11.2,否则执行步骤11.4;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US16/475,583 US20210048807A1 (en) | 2018-03-14 | 2018-06-13 | Method of modeling multi-mode degradation process and predicting remaining useful life |
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810207540.5 | 2018-03-14 | ||
CN201810207540.5A CN108629073B (zh) | 2018-03-14 | 2018-03-14 | 一种多模式的退化过程建模及剩余寿命预测方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
WO2019174142A1 true WO2019174142A1 (zh) | 2019-09-19 |
Family
ID=63706224
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
PCT/CN2018/090928 WO2019174142A1 (zh) | 2018-03-14 | 2018-06-13 | 一种多模式的退化过程建模及剩余寿命预测方法 |
Country Status (3)
Country | Link |
---|---|
US (1) | US20210048807A1 (zh) |
CN (1) | CN108629073B (zh) |
WO (1) | WO2019174142A1 (zh) |
Cited By (13)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111159846A (zh) * | 2019-12-03 | 2020-05-15 | 上海工程技术大学 | 基于分数levy稳定运动的滚动轴承剩余有效寿命预测方法 |
CN112685933A (zh) * | 2020-12-24 | 2021-04-20 | 中国人民解放军海军工程大学 | 一种滚轮丝杠副剩余使用寿命预测方法 |
CN112800616A (zh) * | 2021-02-05 | 2021-05-14 | 中国人民解放军空军工程大学 | 基于比例加速退化建模的设备剩余寿命自适应预测方法 |
CN112949026A (zh) * | 2021-01-19 | 2021-06-11 | 中国人民解放军火箭军工程大学 | 一种考虑年龄和状态依赖的退化设备剩余寿命预测方法 |
CN113011036A (zh) * | 2021-03-26 | 2021-06-22 | 中国人民解放军火箭军工程大学 | 一种针对线性维纳退化过程模型参数的无偏估计方法 |
CN113139276A (zh) * | 2021-03-23 | 2021-07-20 | 温州大学 | 一种印刷机组的两阶段退化分析方法及装置 |
CN113468801A (zh) * | 2021-06-07 | 2021-10-01 | 太原科技大学 | 一种齿轮核密度估计剩余寿命预测方法 |
CN113836741A (zh) * | 2021-09-30 | 2021-12-24 | 中国工程物理研究院研究生院 | 基于多功能系统退化过程的重构和可靠性评估的方法 |
CN114818345A (zh) * | 2022-05-05 | 2022-07-29 | 兰州理工大学 | 一种光伏组件剩余寿命预测方法及预测系统 |
CN114896801A (zh) * | 2022-05-24 | 2022-08-12 | 河南科技大学 | 一种考虑状态增量的关节轴承剩余寿命预测方法 |
CN116227366A (zh) * | 2023-05-08 | 2023-06-06 | 浙江大学 | 两阶段电机绝缘寿命预测方法 |
WO2023216015A1 (zh) * | 2022-05-07 | 2023-11-16 | 潍柴动力股份有限公司 | 一种尿素泵剩余寿命的预测方法及相关装置 |
WO2023240907A1 (zh) * | 2022-06-14 | 2023-12-21 | 广东粤海水务投资有限公司 | 基于分数布朗运动的管道建模方法 |
Families Citing this family (18)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111628901B (zh) * | 2019-02-28 | 2022-11-18 | 华为技术有限公司 | 一种指标异常检测方法以及相关装置 |
CN109992875B (zh) * | 2019-03-28 | 2020-11-17 | 中国人民解放军火箭军工程大学 | 一种切换设备剩余寿命的确定方法及系统 |
CN109960884A (zh) * | 2019-03-28 | 2019-07-02 | 中国人民解放军火箭军工程大学 | 对存在运行状态切换的工程设备寿命预测的方法及系统 |
EP3959571A1 (de) * | 2019-04-23 | 2022-03-02 | Volkswagen Aktiengesellschaft | Verfahren zum bestimmen von restnutzungszyklen, restnutzungszyklusbestimmungsschaltung, restnutzungszyklusbestimmungsvorrichtung |
TWI708197B (zh) * | 2019-04-26 | 2020-10-21 | 國立成功大學 | 生產機台組件的預測保養方法與其電腦程式產品 |
CN112036084B (zh) * | 2020-08-28 | 2022-08-02 | 北京航空航天大学 | 一种相似产品寿命迁移筛选方法和系统 |
CN112257904B (zh) * | 2020-09-29 | 2022-09-30 | 上海工程技术大学 | 一种基于长相关分数阶退化模型预测锂离子电池剩余有效寿命的方法 |
CN112836381B (zh) * | 2021-02-19 | 2023-03-14 | 震兑工业智能科技有限公司 | 一种基于多源信息的船舶剩余寿命预测方法及系统 |
CN113516299B (zh) * | 2021-05-31 | 2024-08-13 | 震兑工业智能科技有限公司 | 一种船舶故障状态预测方法及系统 |
CN113221252B (zh) * | 2021-05-31 | 2023-01-20 | 震兑工业智能科技有限公司 | 一种船舶寿命预测方法及系统 |
CN113610249B (zh) * | 2021-08-13 | 2022-05-20 | 中国石油大学(华东) | 一种全电控井下安全阀的视情维修方法 |
CN113689044B (zh) * | 2021-08-26 | 2024-06-21 | 北京航空航天大学 | 一种开关电源剩余使用寿命预测方法及系统 |
CN114707098B (zh) * | 2022-01-22 | 2024-03-08 | 西北工业大学 | 一种基于多源传感器状态距离的航空发动机性能退化状态评价方法 |
CN118575145A (zh) * | 2022-03-30 | 2024-08-30 | 西门子股份公司 | 计算电子系统的剩余使用寿命的方法、装置及计算机介质 |
CN114896861B (zh) * | 2022-04-02 | 2024-11-01 | 西安交通大学 | 基于平方根容积卡尔曼滤波的滚动轴承剩余寿命预测方法 |
CN114879046B (zh) * | 2022-04-24 | 2023-04-14 | 上海玫克生储能科技有限公司 | 一种基于卡尔曼滤波的锂电池寿命预测方法和系统 |
CN115099164B (zh) * | 2022-08-26 | 2022-12-30 | 中国科学技术大学 | 基于机器学习的锡基钙钛矿薄膜晶体管优化方法及其装置 |
CN115630519B (zh) * | 2022-10-31 | 2023-05-26 | 哈尔滨工业大学 | 基于永磁一致性的极化磁系统式继电器性能退化建模方法 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20070029974A1 (en) * | 2005-08-08 | 2007-02-08 | Toyota Jidosha Kabushiki Kaisha | Powertrain battery life predicting and warning apparatuses |
CN101870076A (zh) * | 2010-07-02 | 2010-10-27 | 西南交通大学 | 一种基于性能退化模型的数控机床导轨副寿命预测方法 |
CN102789545A (zh) * | 2012-07-12 | 2012-11-21 | 哈尔滨工业大学 | 基于退化模型匹配的涡轮发动机剩余寿命的预测方法 |
CN103678858A (zh) * | 2012-09-26 | 2014-03-26 | 中国人民解放军第二炮兵工程大学 | 一种存在竞争失效条件下的设备剩余寿命预测方法 |
CN106021826A (zh) * | 2016-07-11 | 2016-10-12 | 北京航空航天大学 | 一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法 |
CN107358347A (zh) * | 2017-07-05 | 2017-11-17 | 西安电子科技大学 | 基于工业大数据的装备集群健康状态评估方法 |
CN107480442A (zh) * | 2017-08-08 | 2017-12-15 | 山东科技大学 | 依赖于时间和状态的长程相关退化过程剩余寿命预测方法 |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104573881B (zh) * | 2015-02-10 | 2018-01-09 | 广东石油化工学院 | 一种基于退化数据建模的服役设备剩余寿命自适应预测方法 |
CN106484949B (zh) * | 2016-09-12 | 2019-08-16 | 西安理工大学 | 基于退化数据的动量轮可靠性分析与剩余寿命预测方法 |
CN107145645B (zh) * | 2017-04-19 | 2020-11-24 | 浙江大学 | 带不确定冲击的非平稳退化过程剩余寿命预测方法 |
CN107480440B (zh) * | 2017-08-04 | 2020-01-21 | 山东科技大学 | 一种基于两阶段随机退化建模的剩余寿命预测方法 |
-
2018
- 2018-03-14 CN CN201810207540.5A patent/CN108629073B/zh active Active
- 2018-06-13 US US16/475,583 patent/US20210048807A1/en not_active Abandoned
- 2018-06-13 WO PCT/CN2018/090928 patent/WO2019174142A1/zh active Application Filing
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20070029974A1 (en) * | 2005-08-08 | 2007-02-08 | Toyota Jidosha Kabushiki Kaisha | Powertrain battery life predicting and warning apparatuses |
CN101870076A (zh) * | 2010-07-02 | 2010-10-27 | 西南交通大学 | 一种基于性能退化模型的数控机床导轨副寿命预测方法 |
CN102789545A (zh) * | 2012-07-12 | 2012-11-21 | 哈尔滨工业大学 | 基于退化模型匹配的涡轮发动机剩余寿命的预测方法 |
CN103678858A (zh) * | 2012-09-26 | 2014-03-26 | 中国人民解放军第二炮兵工程大学 | 一种存在竞争失效条件下的设备剩余寿命预测方法 |
CN106021826A (zh) * | 2016-07-11 | 2016-10-12 | 北京航空航天大学 | 一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法 |
CN107358347A (zh) * | 2017-07-05 | 2017-11-17 | 西安电子科技大学 | 基于工业大数据的装备集群健康状态评估方法 |
CN107480442A (zh) * | 2017-08-08 | 2017-12-15 | 山东科技大学 | 依赖于时间和状态的长程相关退化过程剩余寿命预测方法 |
Cited By (21)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111159846B (zh) * | 2019-12-03 | 2023-10-31 | 上海工程技术大学 | 基于分数levy稳定运动的滚动轴承剩余有效寿命预测方法 |
CN111159846A (zh) * | 2019-12-03 | 2020-05-15 | 上海工程技术大学 | 基于分数levy稳定运动的滚动轴承剩余有效寿命预测方法 |
CN112685933B (zh) * | 2020-12-24 | 2022-10-18 | 中国人民解放军海军工程大学 | 一种滚轮丝杠副剩余使用寿命预测方法 |
CN112685933A (zh) * | 2020-12-24 | 2021-04-20 | 中国人民解放军海军工程大学 | 一种滚轮丝杠副剩余使用寿命预测方法 |
CN112949026A (zh) * | 2021-01-19 | 2021-06-11 | 中国人民解放军火箭军工程大学 | 一种考虑年龄和状态依赖的退化设备剩余寿命预测方法 |
CN112949026B (zh) * | 2021-01-19 | 2023-05-23 | 中国人民解放军火箭军工程大学 | 一种考虑年龄和状态依赖的退化设备剩余寿命预测方法 |
CN112800616A (zh) * | 2021-02-05 | 2021-05-14 | 中国人民解放军空军工程大学 | 基于比例加速退化建模的设备剩余寿命自适应预测方法 |
CN112800616B (zh) * | 2021-02-05 | 2023-07-18 | 中国人民解放军空军工程大学 | 基于比例加速退化建模的设备剩余寿命自适应预测方法 |
CN113139276A (zh) * | 2021-03-23 | 2021-07-20 | 温州大学 | 一种印刷机组的两阶段退化分析方法及装置 |
CN113011036A (zh) * | 2021-03-26 | 2021-06-22 | 中国人民解放军火箭军工程大学 | 一种针对线性维纳退化过程模型参数的无偏估计方法 |
CN113468801A (zh) * | 2021-06-07 | 2021-10-01 | 太原科技大学 | 一种齿轮核密度估计剩余寿命预测方法 |
CN113468801B (zh) * | 2021-06-07 | 2024-03-26 | 太原科技大学 | 一种齿轮核密度估计剩余寿命预测方法 |
CN113836741A (zh) * | 2021-09-30 | 2021-12-24 | 中国工程物理研究院研究生院 | 基于多功能系统退化过程的重构和可靠性评估的方法 |
CN113836741B (zh) * | 2021-09-30 | 2024-01-26 | 中国工程物理研究院研究生院 | 基于多功能系统退化过程的重构和可靠性评估的方法 |
CN114818345A (zh) * | 2022-05-05 | 2022-07-29 | 兰州理工大学 | 一种光伏组件剩余寿命预测方法及预测系统 |
CN114818345B (zh) * | 2022-05-05 | 2023-09-12 | 兰州理工大学 | 一种光伏组件剩余寿命预测方法及预测系统 |
WO2023216015A1 (zh) * | 2022-05-07 | 2023-11-16 | 潍柴动力股份有限公司 | 一种尿素泵剩余寿命的预测方法及相关装置 |
CN114896801A (zh) * | 2022-05-24 | 2022-08-12 | 河南科技大学 | 一种考虑状态增量的关节轴承剩余寿命预测方法 |
WO2023240907A1 (zh) * | 2022-06-14 | 2023-12-21 | 广东粤海水务投资有限公司 | 基于分数布朗运动的管道建模方法 |
CN116227366A (zh) * | 2023-05-08 | 2023-06-06 | 浙江大学 | 两阶段电机绝缘寿命预测方法 |
CN116227366B (zh) * | 2023-05-08 | 2023-08-11 | 浙江大学 | 两阶段电机绝缘寿命预测方法 |
Also Published As
Publication number | Publication date |
---|---|
US20210048807A1 (en) | 2021-02-18 |
CN108629073A (zh) | 2018-10-09 |
CN108629073B (zh) | 2019-05-03 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
WO2019174142A1 (zh) | 一种多模式的退化过程建模及剩余寿命预测方法 | |
CN110263846B (zh) | 基于故障数据深度挖掘及学习的故障诊断方法 | |
WO2018076571A1 (zh) | Lte网络中的异常值检测方法及系统 | |
CN111046564B (zh) | 两阶段退化产品的剩余寿命预测方法 | |
CN116625438B (zh) | 燃气管网安全在线监测系统及其方法 | |
WO2022126526A1 (zh) | 一种电池温度预测方法及系统 | |
CN110365647B (zh) | 一种基于pca和bp神经网络的虚假数据注入攻击检测方法 | |
CN106888205B (zh) | 一种非侵入式基于功耗分析的plc异常检测方法 | |
WO2019080367A1 (zh) | 一种机械设备健康状态评估方法 | |
Duan et al. | Product technical life prediction based on multi-modes and fractional Lévy stable motion | |
CN105405151A (zh) | 基于粒子滤波和加权Surf的抗遮挡目标跟踪方法 | |
Branisavljević et al. | Improved real-time data anomaly detection using context classification | |
CN103840988A (zh) | 一种基于rbf神经网络的网络流量测量方法 | |
CN109396375A (zh) | 一种基于特征向量和层次聚类的结晶器漏钢预报方法 | |
CN104267610B (zh) | 高精度的高炉冶炼过程异常数据检测及修补方法 | |
CN108446714B (zh) | 一种多工况下的非马尔科夫退化系统剩余寿命预测方法 | |
CN114048546B (zh) | 一种基于图卷积网络和无监督域自适应的航空发动机剩余使用寿命预测方法 | |
CN108436050A (zh) | 一种采用空间密度聚类dbscan预报连铸结晶器漏钢的方法 | |
CN116520236B (zh) | 一种智能电表的异常检测方法和系统 | |
CN116930609A (zh) | 一种基于ResNet-LSTM模型的电能计量误差分析方法 | |
CN117985077B (zh) | 利用无线通信实现铁路沿线施工列车接近智能预警方法 | |
CN105894014A (zh) | 基于多因素不一致度量的异常行为序贯检测方法 | |
CN102508997B (zh) | 一种地下水模型输出不确定性分析方法 | |
CN109101778B (zh) | 基于性能退化数据和寿命数据融合的Wiener过程参数估计方法 | |
CN117056815A (zh) | 基于对比预测编码与支持向量数据的窃电用户检测方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 18909618 Country of ref document: EP Kind code of ref document: A1 |
|
NENP | Non-entry into the national phase |
Ref country code: DE |
|
122 | Ep: pct application non-entry in european phase |
Ref document number: 18909618 Country of ref document: EP Kind code of ref document: A1 |