CN112949026B - Age and state dependence considered degradation equipment residual life prediction method - Google Patents
Age and state dependence considered degradation equipment residual life prediction method Download PDFInfo
- Publication number
- CN112949026B CN112949026B CN202110068209.1A CN202110068209A CN112949026B CN 112949026 B CN112949026 B CN 112949026B CN 202110068209 A CN202110068209 A CN 202110068209A CN 112949026 B CN112949026 B CN 112949026B
- Authority
- CN
- China
- Prior art keywords
- state
- degradation
- parameter
- time
- algorithm
- 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
Links
- 238000006731 degradation reaction Methods 0.000 title claims abstract description 108
- 230000015556 catabolic process Effects 0.000 title claims abstract description 90
- 238000000034 method Methods 0.000 title claims abstract description 52
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 41
- 238000001914 filtration Methods 0.000 claims description 21
- 238000005259 measurement Methods 0.000 claims description 21
- 239000013598 vector Substances 0.000 claims description 15
- 230000000694 effects Effects 0.000 claims description 13
- 238000009499 grossing Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 9
- 238000012544 monitoring process Methods 0.000 claims description 8
- 238000009792 diffusion process Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000007476 Maximum Likelihood Methods 0.000 claims description 4
- 230000005653 Brownian motion process Effects 0.000 claims description 3
- 230000003044 adaptive effect Effects 0.000 claims description 3
- 238000005537 brownian motion Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 238000010845 search algorithm Methods 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims description 3
- 230000036541 health Effects 0.000 description 3
- 238000012360 testing method Methods 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000001747 exhibiting effect Effects 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 208000011580 syndromic disease Diseases 0.000 description 1
- 238000012549 training Methods 0.000 description 1
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
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/04—Ageing analysis or optimisation against ageing
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
- Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
- Y02P90/30—Computing systems specially adapted for manufacturing
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)
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a method for predicting the residual life of degradation equipment considering age and state dependence. Then, the degradation state and unknown parameters are adaptively estimated and updated using a Kalman Filter (KF) and Expectation Maximization (EM) algorithm. Deducing a residual life distribution analysis form which is based on the existing observation value and can be updated in real time based on the estimated degradation state and the model parameters; the residual life prediction accuracy is effectively improved.
Description
Technical Field
The invention belongs to the technical field of reliability engineering, and particularly relates to a degradation equipment residual life prediction method considering age and state dependence.
Background
As one of the key issues of Predictive and Health Management (PHM), estimating the remaining life (RUL) of a degraded device may provide a useful information reference for state-based maintenance decisions, optimal detection intervals, and spare part ordering. To obtain an accurate estimate of the RUL distribution, a suitable degradation model is needed that adequately characterizes the actual degradation process.
Nonlinearity and randomness are important factors that must be considered when estimating RUL in a random degradation model framework, and random process models are typically used to describe the randomness of degradation. At the same time, individual variability is also an important factor to consider in the estimation of RUL. Individual variability may be defined as different individuals in the same batch of equipment exhibiting different rates of degradation due to different operating environments and health conditions. Furthermore, due to the complexity of the system or the high cost of directly observing the degradation state, implicit or partially observable degradation processes are often encountered in engineering practice. The above-mentioned problems have been studied to some extent in the prior art, and it is noted that in the prior art and the invention, the degradation process is aimed at a device which depends only on age, but there is still a large number of degradation processes which are not only affected by the age of the device but also on the health state of the device itself. Taking fatigue crack propagation as an example, the degradation rate at the initial stage of the fatigue crack is low. As fatigue cracks grow, the stress strength will gradually increase, resulting in an increased degradation rate. Currently, there is no related study of a degradation apparatus remaining life prediction method that considers both age and state in an individual difference and implicit degradation state.
Disclosure of Invention
The invention aims to provide a degradation device with age and state dependence in a degradation process, and provides a degradation device residual life prediction method with age and state dependence in consideration of the individual variability of the device and the influence of an implicit degradation state.
The technical scheme adopted by the invention is as follows:
a method for predicting remaining life of a degraded device considering age and state dependence, comprising the steps of:
step 4, deducing a residual life distribution analysis form which is based on the existing observation value and can be updated in real time based on the degradation state and the model parameters estimated in the step 3;
and 5, when a new observed value is generated, updating parameters in the state space model in real time, and substituting the parameters into an expression of the residual life distribution of the equipment, so that life prediction of the degraded equipment is realized.
Preferably, in the step 1, the degradation model building process is as follows:
let X (t) denote the degradation amount of the sample at time t, then the degradation process { X (t), t.gtoreq.0 } based on the diffusion process can be expressed as:
dX(t)=μ(X(t),t;θ)dt+σ B dB(t) (1)
where μ (X (t), t; θ) represents a nonlinear drift coefficient function while depending on age and state, μ (X, t; θ) =ax+bh (t; ζ), where θ= (a, b, ζ), and h (t) is a function with respect to time t, ζ is an unknown parameter vector in the function h (t; ζ); wherein, the parameter a is a fixed constant, and the parameter b is a random effect coefficient representing individual variability, namely b:μ b is the mean value of parameter b,/->Is the variance of parameter b; sigma (sigma) B Is a diffusion coefficient, θ is a parameter vector, B (t) is a standard Brownian motion, and B (t): N (0, t); without loss of generality, when t=0, X (0) =0;
the relation between the uncertain observed value and the potential degradation state at the moment t is described by { Y (t), t is more than or equal to 0}, and the expression is as follows:
Y(t)=X(t)+ε (2)
for the random degradation process given by equation (1), in the first time sense, the lifetime T of the device can be defined as:
T=inf{t:X(t)≥ω|X(0)<ω} (3)
wherein εFor measuring errors, the service life T is a random variable, and omega is a failure threshold value of equipment; the probability density function of the hit T can be f T (t) represents;
to utilize the state monitoring information, it is assumed that the device is at discrete time point 0=t 0 <t 1 <t 2 <L<t k Monitoring to t k By the time a set of degradation measurements is recorded as Y 1:k ={y 1 ,y 2 ,L,y k}, wherein yk =Y(t k ) Representing t k The observed value at the moment, the real degradation state corresponding to the observed value is expressed as X 1:k ={x 1 ,x 2 ,L,x k}, wherein xk =X(t k ) The method comprises the steps of carrying out a first treatment on the surface of the Defining t according to the concept of the time of arrival k Time of day, remaining life L k Is that
L k =inf{l k >0:X(l k +t k )≥ω|x k <ω} (4)
Preferably, in the step 2, the state space model building process is as follows:
the residual life estimation takes into account the influence of individual variability and measurement errors during degradation, and therefore only the observed value Y 1:k Available, actual degradation state x k Cannot be directly utilized;
random effect parameter presence update procedure b k =b k-1 +α, wherein α:b 0 =b:/>the initial distribution is obtained by using the observation value to obtain posterior distribution of the parameter b; at the same time, the state equation and the measurement equation are converted into discrete time equations for state estimation when new observation data is available, and then the state equationThe equation (1) and the observation equation (2) may be at discrete points in time t k The transition at k=1, 2, l is a state space model as follows:
wherein ,v k =σ B [B(t k )-B(t k-1 )],ε k for epsilon at t k A specific value of time; { v k } k≥1 and {εk } k≥1 Is an independent co-distributed noise sequence, and v k :/>ε k :/>
Let z k =x k -ax k Then x k =z k /(1-a)=βz k The method comprises the steps of carrying out a first treatment on the surface of the Thus, a new potential degradation state z k And a random parameter b k Can be regarded as a new implicit state, requiring dependency on the observed value Y 1:k Estimating; based on the state space model (5), a new implicit state can be estimated by using Kalman filtering, and the state space model (5) can be rewritten as:
wherein ,sk ∈R 2×1 ,η k ∈R 2×1 ,A k ∈R 2×2 ,C∈R 1×2 ,η k :N(0,Q k ) In particular, to
Preferably, in step 3, implicit degradation state estimation:
based on a state space model (6), implicit states are estimated using a Kalman filtering algorithm, first s is defined k Based on t k The conditional expectation and variance of the time-of-day available observations are as follows:
accordingly, the expected sum variance obtained by one-step prediction is respectively
According to the definition above, the Kalman filtering algorithm may be summarized as follows:
(1) State estimation
(2) Variance update
P k|k =P k|k-1 -K(k)Cr′ k|k-1 P k|k-1
Wherein, the initial state of the Kalman filtering is set as follows:
from formula (6), s k The conditional probability density function of (2) is a bivariate gaussian distribution, i.e. s k :
Then, the potential degradation state z k And random effect parameter b k The posterior distribution of (a) is the posterior distribution of (b) with time t k Relatedly, based on the property of the double-variable Gaussian distribution, it is possible to obtain
wherein ,
preferably, in step 3, adaptive parameter estimation:
representing unknown parameters in the state space model (6) as vectorsBy iterative calculation and maximization of conditional expectation of log likelihood function, expectation maximization algorithm can generate parameter estimation sequence converging to maximum likelihood estimation of parameter, let ∈ ->Estimated value representing unknown parameters in jth iteration of the expectation maximization algorithm, +.>Representing a conditional expectation operator;
the expectation maximization algorithm comprises the steps of:
(1) Calculate the j+1st iteration value
wherein ,l(s1:k ,Y 1:k |Θ)=lnp(s 1:k ,Y 1:k |Θ);
for the state space model (4), t is cut off k Time of day, implicit degradation state sequence s 1:k And observing sequence Y 1:k The joint log likelihood function of (2) can be expressed as:
next, the conditional expectation of the log-likelihood function (14) is calculated, and it is possible to obtain:
to calculateThe conditional expectation of each term in equation (17) needs to be deduced; first, the definition is as follows:
the method can be obtained through complex algebraic operation:
wherein ,
obviously, for the purpose of calculationNeed to get +.>P i|k ,P i-1|k and Mi|k The method comprises the steps of carrying out a first treatment on the surface of the These conditions are expected to be calculated by a kalman filter algorithm and an RTS smoothing algorithm; the specific steps of the RTS smoothing algorithm are as follows: />
(2) Backward smoothing iterations, for i=k, k-1, l,1, backward recursion
(3) Initialization of
M k|k =[I-K(k)C]A k-1 P k-1|k-1 (20)
(4) Updating covariance matrix
In order to reduce the complexity of the parameter estimation algorithm,can be divided into three parts, the first part only contains the initial value of the Kalman filtering algorithm +.>Can be expressed as:
respectively regarding the formula (22) as s 0 and P0|0 The partial derivative is calculated and is set to 0, so that the estimated value of the j+1st iteration can be obtained:
unknown parametersThe estimation can be performed by a section likelihood function method, and in order to facilitate the solution of the partial derivative, the formula (24) is rewritten as:
assuming that ζ is fixed, the formula (25) is respectively related to and />Taking the partial derivative and letting it be 0, one can get:
substituting equations (26) and (27) into (25) to obtain a section likelihood function about xi, and obtaining an estimated value of j+1th iteration of xi by using a search algorithm;
the third part contains the parameter vector Θ in the observation equation 3 ={β,σ ε -can be expressed as:
unknown parameter theta 3 ={β,σ ε Estimation method and apparatusSimilarly, assuming β is fixed, formula (28) is given for +.>Taking the partial derivative and letting it be 0, one can get:
substituting formula (29) into formula (28) to maximize the cross-sectional likelihood function with respect to beta, thereby obtaining an estimated valueSince β=1/(1-a), it is possible to obtain +.>
Preferably, in step 4, for the degradation model (1) equation, while taking into account the personalized differences between the devices and the measurement errors, t is defined according to the remaining lifetime definition (4) equation k Probability density function of time remaining lifeCan be expressed as:
wherein ,
the invention has the beneficial effects that: the invention utilizes a nonlinear diffusion process to describe a degradation process of degradation equipment with age and state dependent simultaneously, and provides a degradation modeling and residual life prediction method which simultaneously considers individual variability and implicit degradation state of the equipment. First, by constructing a state space model, the relationship between the implicit degradation state and the observed value is described. Then, the degradation state and unknown parameters are adaptively estimated and updated using a Kalman Filter (KF) and Expectation Maximization (EM) algorithm. Based on the estimated degradation state and model parameters, a residual life distribution analysis form which is based on the existing observation value and can be updated in real time is deduced; the prediction accuracy of the residual life is effectively improved.
Drawings
In order to more clearly illustrate the embodiments of the invention or the technical solutions in the prior art, the drawings that are required in the embodiments or the description of the prior art will be briefly described, it being obvious that the drawings in the following description are only some embodiments of the invention, and that other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a graph of gyroscope degradation under step acceleration stress;
FIG. 2 is a graph of probability density of remaining life of a # 2 bearing;
FIG. 3 is a graph of probability density of remaining life of a 3# bearing;
fig. 4 is a 4# bearing remaining life probability density function.
Detailed Description
For the purpose of making the objects, technical solutions and advantages of the embodiments of the present invention more apparent, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention, and it is apparent that the described embodiments are some embodiments of the present invention, but not all embodiments of the present invention.
Thus, the following detailed description of the embodiments of the invention, as presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
The invention specifically provides a degradation equipment residual life prediction method considering age and state dependence, which comprises the following steps:
the degradation model building process is as follows:
let X (t) denote the performance degradation of the sample at time t, then the degradation model { X (t), t.gtoreq.0 } based on the diffusion process can be expressed as:
dX(t)=μ(X(t),t;θ)dt+σ B dB(t) (1)
where μ (X (t), t; θ) represents a nonlinear drift coefficient function while depending on age and state, μ (X, t; θ) =ax+bh (t; ζ), where θ= (a, b, ζ), and h (t) is a function with respect to time t, ζ is an unknown parameter vector in the function h (t; ζ); wherein, the parameter a is a fixed constant, and the parameter b is a random effect coefficient representing individual variability, namely b:μ b is the mean value of parameter b,/->Is the variance of parameter b; sigma (sigma) B Is a diffusion coefficient, θ is a parameter vector, B (t) is a standard Brownian motion, and B (t): N (0, t); without loss of generality, when t=0, X (0) =0;
in addition, perfect measurement of the potential degradation state is often difficult to achieve, and measurement errors caused by external interference, noise, non-ideal instruments and other factors are inevitably generated in the measurement process. In this case, the observation only partially reflects the potential degradation state. To characterize the effect of measurement uncertainty, the relationship between the uncertainty observations and the potential degradation state at time t is described by { Y (t), t.gtoreq.0 } as follows:
Y(t)=X(t)+ε (2)
for the random degradation process given by equation (1), in the first time sense, the lifetime T of the device can be defined as:
T=inf{t:X(t)≥ω|X(0)<ω} (3)
where ε is the measurement error, life T is a random variable, ω is the failure threshold of the device, and is typically determined by industry standards and expert experience. Probability Density Function (PDF) of lifetime T is available f T (t) represents.
The main object of the present invention is to estimate and update the RUL distribution of a single device in operation based on real-time observations of the degradation process. To utilize the state monitoring information, it is assumed that the device is at discrete time point 0=t 0 <t 1 <t 2 <L<t k Monitoring to t k By the time a set of degradation measurements is recorded as Y 1:k ={y 1 ,y 2 ,L,y k}, wherein yk =Y(t k ) Representing t k Observations of time of day. The true degradation state corresponding to this is denoted as X 1:k ={x 1 ,x 2 ,L,x k}, wherein xk =X(t k ). Defining t according to the concept of the time of arrival k Time of day, remaining life L k Is that
L k =inf{l k >0:X(l k +t k )≥ω|x k <ω} (4)
for a practical system, our goal is to at a specific point in time t k Estimating remaining lifetime (RUL), which is the kth monitored data from the start time, and the degradation state is x k =X(t k ). To obtain the monitoring t of the system k The remaining lifetime (RUL) at the moment is introduced as followsAnd (5) managing.
dU(l k )=μ(X(l k +t k ),l k +t k ;θ)dl k +σ B dB(l k )
Given random effect parameters b and t k Implicit degradation observations of time x k Residual life ofThe PDF of (C) can be approximated as +.>
For simplicity, let eta (u, l k ;θ)=μ(X(l k +t k ),l k +t k ;θ)=μ(U(l k )+x k ,l k +t k The method comprises the steps of carrying out a first treatment on the surface of the θ). Then and />Can be expressed as
In the above derivation, it is assumed that no random effect exists, and an implicit degradation state can be directly observed. However, there are individual differences in the different devices and it is often not possible to measure the potential degradation state perfectly in practical engineering.
For the effect of individual variability and measurement errors in the degradation process in the estimation of the remaining lifetime (RUL), in this case only the observed value Y 1:k Available, actual degradation state x k Cannot be directly utilized. The main purpose of the invention is based on the observed value Y 1:k Deriving Probability Density Function (PDF) for remaining life (RUL)
Random effect parameter b presence update procedure b k =b k-1 +α, wherein α:b 0 =b:/>the initial distribution is used to obtain a posterior distribution of the parameter b. Meanwhile, in order to recognize an implicit degradation state, it is necessary to convert a state equation and a measurement equation into a discrete time equation in order to perform state estimation when new observation data is available. Then, the state equation (1) and the observation equation (2) may be at discrete time points t k ,k=1,2,The transition at L is to a state space model as follows:
wherein ,v k =σ B [B(t k )-B(t k-1 )],ε k for epsilon at t k A specific value of time; { v k } k≥1 and {εk } k≥1 Is an independent co-distributed noise sequence, and v k :/>ε k :/>
Let z k =x k -ax k Then x k =z k /(1-a)=βz k The method comprises the steps of carrying out a first treatment on the surface of the Thus, a new potential degradation state z k And a random parameter b k Can be regarded as a new implicit state, requiring dependency on the observed value Y 1:k Estimating; let z k =x k -ax k Then x k =z k /(1-a)=βz k The method comprises the steps of carrying out a first treatment on the surface of the Thus, a new potential degradation state z k And random parameter b k Can be regarded as a new implicit state, requiring dependency on the observed value Y 1:k Estimating; based on the state space model (5), a new implicit state can be estimated by using Kalman filtering, and the state space model (5) can be rewritten as:
wherein ,sk ∈R 2×1 ,η k ∈R 2×1 ,A k ∈R 2×2 ,C∈R 1×2 ,η k :N(0,Q k ) In particular, to
implicit degradation state estimation:
based on a state space model (6), implicit states are estimated using a Kalman filtering algorithm, first s is defined k Based on t k The conditional expectation and variance of the time-of-day available observations are as follows:
accordingly, the expected sum variance obtained by one-step prediction is respectively
According to the definition above, the Kalman filtering algorithm may be summarized as follows:
(1) State estimation
(2) Variance update
P k|k =P k|k-1 -K(k)Cr′ k|k-1 P k|k-1
Wherein, the initial state of the Kalman filtering is set as follows:
from formula (6), s k The conditional probability density function of (2) is a bivariate gaussian distribution, i.e. s k :
Then, the potential degradation state z k And random effect parameter b k The posterior distribution of (a) is the posterior distribution of (b) with time t k Relatedly, based on the property of the double-variable Gaussian distribution, it is possible to obtain
wherein ,
and (3) adaptive parameter estimation:
for convenience, the unknown parameters in the model (6) are represented as vectorsEstimating Θ using the Expectation Maximization (EM) algorithm provides a basic framework for solving the Maximum Likelihood Estimation (MLE) problem of the existence of implicit states.
By iterative calculation and maximization of conditional expectation of log likelihood function, expectation maximization algorithm can generate parameter estimation sequence converging to maximum likelihood estimation of parameter, letEstimated value representing unknown parameters in the jth iteration of the expectation maximization algorithm, ++>Representing a conditional expectation operator.
In summary, the EM algorithm includes the following two steps:
(1) E, step E: calculate the j+1st iteration value
wherein ,l(s1:k ,Y 1:k |Θ)=ln p(s 1:k ,Y 1:k |Θ);
To be used forThe above two steps are iterated repeatedly until the convergence condition is satisfied. />
For the state space model (6), t is cut off k Time of day, implicit degradation state sequence s 1:k And observing sequence Y 1:k The joint log likelihood function of (2) can be expressed as:
next, the conditional expectation of the log-likelihood function (16) is calculated, and it is possible to obtain:
to calculateThe conditional expectation of each term in equation (17) needs to be deduced; first, the definition is as follows:
the method can be obtained through complex algebraic operation:
wherein ,
obviously, for the purpose of calculationNeed to get +.>P i|k ,P i-1|k and Mi|k The method comprises the steps of carrying out a first treatment on the surface of the These conditions are expected to be calculated by a kalman filter algorithm and an RTS smoothing algorithm; the specific steps of the RTS smoothing algorithm are as follows:
(2) Backward smoothing iterations, for i=k, k-1, l,1, backward recursion
(3) Initialization of
M k|k =[I-K(k)C]A k-1 P k-1|k-1 (20)
(4) Updating covariance matrix
In order to reduce the complexity of the parameter estimation algorithm,can be divided into three parts, the first part only contains the initial value of the Kalman filtering algorithm +.>Can be expressed as:
respectively regarding the formula (22) as s 0 and P0|0 The partial derivative is calculated and is set to 0, so that the estimated value of the j+1st iteration can be obtained:
unknown parametersThe estimation can be performed by a section likelihood function method, and in order to facilitate the solution of the partial derivative, the formula (24) is rewritten as: />
Assuming that ζ is fixed, the formula (25) is respectively related to and />Taking the partial derivative and letting it be 0, one can get:
substituting equations (26) and (27) into (25) to obtain a section likelihood function about xi, and obtaining an estimated value of j+1th iteration of xi by using a search algorithm;
the third part contains the parameter vector Θ in the observation equation 3 ={β,σ ε -can be expressed as:
unknown parameter theta 3 ={β,σ ε Estimation method and apparatusSimilarly, assuming β is fixed, formula (28) is given for +.>Taking the partial derivative and letting it be 0, one can get:
substituting formula (29) into formula (28) to maximize the cross-sectional likelihood function with respect to beta, thereby obtaining an estimated valueSince β=1/(1-a), it is possible to obtain +.>
Step 4, deducing a residual life distribution analysis form which is based on the existing observation value and can be updated in real time based on the degradation state and the model parameters estimated in the step 3
wherein ,p(sk |Y 1:k ) Is s k |Y 1:k Is the same as the condition expectation of p (z k |b k ,Y 1:k) and p(bk |Y 1:k ) Z respectively k |b k ,Y 1:k and bk |Y 1:k Is desirable for the conditions of (2).
To derive the distribution of remaining life, the following quotients are introduced.
The quotients 2 there is Z: N (μ, σ) 2 ) And ω, β, L, M εi, Q εi + Can obtain
From this, the Probability Density Function (PDF) of the remaining life in the time-of-arrival sense can be derived using a full probability formula, the specific remaining life expression being as follows:
for the degradation process (1), the individual variability and measurement errors between devices are taken into account; according to the definition of remaining life (RUL), L k =inf{l k >0:X(l k +t k )≥ω|x k < ω }; then t k Probability Density Function (PDF) of time remaining life (RUL)Can be expressed as:
wherein ,
and (3) proving: PDF with residual life can be obtained according to the quotation 1Can be approximated as
Will x k =βz k Substituted into the above formula to obtainAnd simplifying the method, and the result is as follows: />
wherein ,
from formula (12), z k |b k ,Y 1:k :Based on the sum of the full probability formulas, using the lemma 2, we can get:
Further, from the formula (12) and the formula (13), b is known k |Y 1:k :Based on the sum of the full probability formulas, using the lemma 2, we can get:
wherein ,
the syndrome is known.
Example application analysis
Taking degradation data of a certain type of bearing as an example, the effectiveness of the method provided by the invention is verified. The test data for the four samples are shown in figure 1.
The last 25 observations of the 2 nd-4 th bearing are used as a test data set, and the rest of degradation data are used as training data sets. According to the parameter estimation method proposed in step 3, an estimated value of an unknown parameter in the model can be obtained, as shown in table 1.
TABLE 1 degradation model parameter estimation
To illustrate the effectiveness of this method in residual life (RUL) prediction, samples may be subjected to residual life (RUL) prediction based on initial estimates of model parameters in table 1. In the RUL prediction process, once new degradation observations are available, the model parameters are updated, resulting in a more accurate residual life (RUL) distribution. In fig. 2-4, a Probability Density Function (PDF) of the remaining life (RUL), a predicted average remaining life (RUL), and an actual remaining life (RUL) at that time are shown for comparison.
For three bearings, the predicted residual life of the method is very close to the actual residual life, and the probability density function of the predicted residual life becomes higher and sharper along with the updating of model parameters, which shows that the uncertainty of the predicted result is gradually reduced.
The effectiveness and the superiority of the equipment degradation modeling and residual life prediction method considering age and state dependence are further verified.
The foregoing is merely illustrative of the present invention and not restrictive, and other modifications and equivalents thereof may occur to those skilled in the art without departing from the spirit and scope of the present invention.
Claims (4)
1. A method for predicting remaining life of a degraded device in consideration of age and state dependence, comprising the steps of:
step 1, a degradation model is established, dependency of equipment degradation, age and state is described by using a drift function, and individuation difference and measurement error between equipment in the degradation process are considered;
step 2, on the basis of the step 1, a state space model is established, and the relation between the implicit degradation state and the observed value is described;
step 3, estimating an implicit degradation state by using a Kalman filtering algorithm, and estimating parameters in a state space model by using an expected maximization algorithm;
step 4, deducing a residual life distribution analysis form which is based on the existing observation value and can be updated in real time based on the degradation state and the model parameters estimated in the step 3;
step 5, when a new observation value is generated, updating parameters in the state space model in real time, and substituting the parameters into an expression of the residual life distribution of the equipment so as to realize life prediction of the degraded equipment;
in the step 1, the degradation model building process is as follows:
let X (t) denote the degradation amount of the sample at time t, then the degradation process { X (t), t.gtoreq.0 } based on the diffusion process can be expressed as:
dX(t)=μ(X(t),t;θ)dt+σ B dB(t) (1)
where μ (X (t), t; θ) represents a nonlinear drift coefficient function while depending on age and state, μ (X, t; θ) =ax+bh (t; ζ), where θ= (a, b, ζ), and h (t) is a function with respect to time t, ζ is an unknown parameter vector in the function h (t; ζ); wherein, the parameter a is a fixed constant, and the parameter b is a random effect coefficient representing individual variability, namely b:μ b is a ginsengMean value of number b>Is the variance of parameter b; sigma (sigma) B For diffusion coefficients, θ is a parameter vector, B (t) is a standard Brownian motion, and B (t): n (0, t); without loss of generality, when t=0, X (0) =0;
the relation between the uncertain observed value and the potential degradation state at the moment t is described by { Y (t), t is more than or equal to 0}, and the expression is as follows:
Y(t)=X(t)+ε (2)
for the random degradation process given by equation (1), in the first time sense, the lifetime T of the device can be defined as:
T=inf{t:X(t)≥ω|X(0)<ω} (3)
wherein epsilon is a measurement error, the service life T is a random variable, and omega is a failure threshold of the equipment; the probability density function of the hit T can be f T (t) represents;
to utilize the state monitoring information, it is assumed that the device is at discrete time point 0=t 0 <t 1 <t 2 <…<t k Monitoring to t k By the time a set of degradation measurements is recorded as Y 1:k ={y 1 ,y 2 ,…,y k}, wherein yk =Y(t k ) Representing t k The observed value of the moment, the real degradation state corresponding to the observed value is expressed as X 1:k ={x 1 ,x 2 ,…,x k}, wherein xk =X(t k ) The method comprises the steps of carrying out a first treatment on the surface of the Defining t according to the concept of the time of arrival k Time of day, remaining life L k Is that
L k =inf{l k >0:X(l k +t k )≥ω|x k <ω} (4)
In the step 2, the state space model building process is as follows:
residual life estimationThe influence of individual variability and measurement errors during degradation is taken into account, so that only the observed value Y 1:k Available, actual degradation state x k Cannot be directly utilized;
random effect parameter presence update procedure b k =b k-1+α, wherein ,the initial distribution is obtained by using the observation value to obtain posterior distribution of the parameter b; at the same time, the state equation and the measurement equation are converted into discrete-time equations for state estimation when new observation data is available, and then the state equation (1) and the observation equation (2) can be obtained at discrete time points t k K=1, 2, …), where the transition is to a state space model, as follows:
wherein ,v k =σ B [B(t k )-B(t k-1 )],ε k for epsilon at t k A specific value of time; { v k } k≥1 and {εk } k≥1 Is an independent and equidistributed noise sequence, and +.>
Let z k =x k -ax k Then x k =z k /(1-a)=βz k The method comprises the steps of carrying out a first treatment on the surface of the Thus, a new potential degradation state z k And random parameter b k Can be regarded as a new implicit state, requiring dependency on the observed value Y 1:k Estimating; based on the state space model (5), a new implicit state can be estimated by using Kalman filtering, and the state space model (5) can be rewritten as:
wherein ,sk ∈R 2×1 ,η k ∈R 2×1 ,A k ∈R 2×2 ,C∈R 1×2 ,η k :N(0,Q k ) In particular, to
2. A method for predicting remaining life of a degraded device taking into account age and state dependence as claimed in claim 1, wherein in step 3, the degraded state estimation is implicit:
based on a state space model (6), implicit states are estimated using a Kalman filtering algorithm, first s is defined k Based on t k The conditional expectation and variance of the time-of-day available observations are as follows:
accordingly, the expected sum variance obtained by one-step prediction is respectively
According to the definition above, the Kalman filtering algorithm may be summarized as follows:
(1) State estimation
(2) Variance update
P kk =P kk-1 -K(k)Cr k ′ k-1 P kk-1
Wherein, the initial state of the Kalman filtering is set as follows:
from formula (6), s k The conditional probability density function of (2) is a bivariate gaussian distribution, i.eThen, the potential degradation state z k And random effect parameter b k The posterior distribution of (a) is the posterior distribution of (b) with time t k Relatedly, based on the nature of the bivariate Gaussian distribution, one can obtain
wherein ,
3. a method for predicting remaining life of a degraded equipment taking into account age and state dependency as claimed in claim 2, wherein in step 3, the adaptive parameter estimates:
representing unknown parameters in the state space model (6) as vectorsBy iterative calculation and maximization of conditional expectation of log likelihood function, expectation maximization algorithm can generate parameter estimation sequence converging to maximum likelihood estimation of parameter, let ∈ ->Representing an estimate of the unknown parameter in the jth iteration of the expectation maximization algorithm,representing a conditional expectation operator; />
The expectation maximization algorithm comprises the steps of:
(1) Calculate the j+1st iteration value
wherein ,l(s1:k ,Y 1:k |Θ)=lnp(s 1:k ,Y 1:k |Θ);
for the state space model (4), t is cut off k Time of day, implicit degradation state sequence s 1:k And observing sequence Y 1:k The joint log likelihood function of (2) can be expressed as:
next, the conditional expectation of the log-likelihood function (14) is calculated, and it is possible to obtain:
to calculateThe conditional expectation of each term in equation (17) needs to be deduced; head partFirst, the definition is as follows:
the method can be obtained through complex algebraic operation:
wherein ,
obviously, for the purpose of calculationNeed to get +.>P i|k ,P i-1|k and Mi|k The method comprises the steps of carrying out a first treatment on the surface of the These conditions are expected to be calculated by Kalman filteringCalculating by a method and an RTS smoothing algorithm; the specific steps of the RTS smoothing algorithm are as follows:
(2) Backward smoothing iterations, for i=k, k-1, …,1, backward recursion
(3) Initialization of
M k|k =[I-K(k)C]A k-1 P k-1|k-1 (20)
(4) Updating covariance matrix
In order to reduce the complexity of the parameter estimation algorithm,can be divided into three parts, the first part only contains the initial value of the Kalman filtering algorithm +.>Can be expressed as:
respectively regarding the formula (22) as s 0 and P0|0 The partial derivative is calculated and is set to 0, so that the estimated value of the j+1st iteration can be obtained:
Unknown parametersThe estimation can be performed by a section likelihood function method, and in order to facilitate the solution of the partial derivative, the formula (24) is rewritten as:
assuming that ζ is fixed, the formula (25) is respectively related to and />Taking the partial derivative and letting it be 0, one can get:
substituting equations (26) and (27) into (25) to obtain a section likelihood function about ζ, and obtaining an estimated value of j+1st iteration of ζ by using a search algorithm;
the third part contains the parameter vector Θ in the observation equation 3 ={β,σ ε -can be expressed as:
unknown parameter theta 3 ={β,σ ε Estimation method and apparatusSimilarly, assuming β is fixed, formula (28) is given for +.>Taking the partial derivative and letting it be 0, one can get:
substituting formula (29) into formula (28) to maximize the cross-sectional likelihood function with respect to beta, thereby obtaining an estimated valueSince β=1/(1-a), it is possible to obtain +.>
4. A method for predicting the remaining life of a degradation device in consideration of age and state dependency as claimed in claim 3, wherein in step 4, for the degradation model (1) equation, while taking into consideration the personalized difference and measurement error between devices, t is defined according to the remaining life definition (4) equation k Probability density function of time remaining lifeCan be expressed as:
wherein ,
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110068209.1A CN112949026B (en) | 2021-01-19 | 2021-01-19 | Age and state dependence considered degradation equipment residual life prediction method |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110068209.1A CN112949026B (en) | 2021-01-19 | 2021-01-19 | Age and state dependence considered degradation equipment residual life prediction method |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112949026A CN112949026A (en) | 2021-06-11 |
CN112949026B true CN112949026B (en) | 2023-05-23 |
Family
ID=76235676
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110068209.1A Active CN112949026B (en) | 2021-01-19 | 2021-01-19 | Age and state dependence considered degradation equipment residual life prediction method |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112949026B (en) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113656746B (en) * | 2021-07-21 | 2022-06-17 | 东南大学 | Travel mode chain selection method considering group heterogeneity under dynamic structure |
CN114063456B (en) * | 2021-11-15 | 2023-06-02 | 哈尔滨工业大学 | Fault prediction and early warning method using autoregressive model and Kalman filtering algorithm |
CN114417686B (en) * | 2022-01-20 | 2023-02-03 | 哈尔滨工业大学 | Self-adaptive online residual service life prediction method for single lithium ion battery |
CN114779088B (en) * | 2022-04-20 | 2023-05-12 | 哈尔滨工业大学 | Battery remaining service life prediction method based on maximum expected-unscented particle filtering |
CN114896801B (en) * | 2022-05-24 | 2024-07-09 | 河南科技大学 | Method for predicting residual life of knuckle bearing by considering state increment |
CN116029164B (en) * | 2023-03-30 | 2023-07-18 | 中国人民解放军火箭军工程大学 | Method and system for determining degradation degree of device performance, electronic device and storage medium |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107145720A (en) * | 2017-04-19 | 2017-09-08 | 浙江大学 | It is continuous to degenerate and the unknown equipment method for predicting residual useful life impacted under collective effect |
CN107145645A (en) * | 2017-04-19 | 2017-09-08 | 浙江大学 | The non-stationary degenerative process method for predicting residual useful life of the uncertain impact of band |
CN109977552A (en) * | 2019-03-28 | 2019-07-05 | 中国人民解放军火箭军工程大学 | A kind of equipment method for predicting residual useful life and system considering that state-detection influences |
WO2019174142A1 (en) * | 2018-03-14 | 2019-09-19 | 山东科技大学 | Multi-mode degradation process modelling and remaining service life prediction method |
CN111368403A (en) * | 2020-02-24 | 2020-07-03 | 西安交通大学 | Self-adaptive non-linear degradation residual life prediction method |
CN111460692A (en) * | 2020-04-26 | 2020-07-28 | 中国人民解放军火箭军工程大学 | Equipment residual life prediction method and system considering degradation rate mutual influence |
CN112231925A (en) * | 2020-10-27 | 2021-01-15 | 山东科技大学 | Residual life prediction method considering state dependence time lag |
-
2021
- 2021-01-19 CN CN202110068209.1A patent/CN112949026B/en active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107145720A (en) * | 2017-04-19 | 2017-09-08 | 浙江大学 | It is continuous to degenerate and the unknown equipment method for predicting residual useful life impacted under collective effect |
CN107145645A (en) * | 2017-04-19 | 2017-09-08 | 浙江大学 | The non-stationary degenerative process method for predicting residual useful life of the uncertain impact of band |
WO2019174142A1 (en) * | 2018-03-14 | 2019-09-19 | 山东科技大学 | Multi-mode degradation process modelling and remaining service life prediction method |
CN109977552A (en) * | 2019-03-28 | 2019-07-05 | 中国人民解放军火箭军工程大学 | A kind of equipment method for predicting residual useful life and system considering that state-detection influences |
CN111368403A (en) * | 2020-02-24 | 2020-07-03 | 西安交通大学 | Self-adaptive non-linear degradation residual life prediction method |
CN111460692A (en) * | 2020-04-26 | 2020-07-28 | 中国人民解放军火箭军工程大学 | Equipment residual life prediction method and system considering degradation rate mutual influence |
CN112231925A (en) * | 2020-10-27 | 2021-01-15 | 山东科技大学 | Residual life prediction method considering state dependence time lag |
Non-Patent Citations (3)
Title |
---|
Jing Xu ; Yang Lei ; Bin Liu ; Chao Ji ; Lijun Nan.Application of PHM Technology in the Design of Tank Fire Control System.《2018 Prognostics and System Health Management Conference (PHM-Chongqing)》.2019,第249-253页. * |
不完美维护活动干预下的设备剩余寿命估计;胡昌华等;《中国惯性技术学报》;20161015(第05期);全文 * |
基于机器学习的设备剩余寿命预测方法综述;裴洪胡昌华司小胜张建勋庞哲楠张鹏;《机械工程学报》;20190401;第55卷(第8期);第1-13页 * |
Also Published As
Publication number | Publication date |
---|---|
CN112949026A (en) | 2021-06-11 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112949026B (en) | Age and state dependence considered degradation equipment residual life prediction method | |
CN107145645B (en) | Method for predicting residual life of non-stationary degradation process with uncertain impact | |
CN107765347B (en) | Short-term wind speed prediction method based on Gaussian process regression and particle filtering | |
CN110851980B (en) | Method and system for predicting residual life of equipment | |
CN103488881B (en) | Equipment residual service life prediction method under the condition of uncertain degradation measured data | |
CN109241609B (en) | Bayesian dynamic prediction method based on Markov chain Monte Carlo | |
CN108304685A (en) | A kind of non-linear degradation equipment method for predicting residual useful life and system | |
CN105205288B (en) | Satellite based on the Evolution Modes prediction technique of state in orbit for a long time | |
CN112800616B (en) | Equipment residual life self-adaptive prediction method based on proportional acceleration degradation modeling | |
CN104615866B (en) | A kind of life-span prediction method based on physical-statistical model | |
Si et al. | A general stochastic degradation modeling approach for prognostics of degrading systems with surviving and uncertain measurements | |
CN104462015B (en) | Process the fractional order linear discrete system state updating method of non-gaussian L é vy noises | |
Badescu et al. | A marked Cox model for the number of IBNR claims: estimation and application | |
CN109918776B (en) | Fatigue crack prediction method of mechanical product based on two-step least square method | |
Schirru et al. | Particle filtering of hidden gamma processes for robust predictive maintenance in semiconductor manufacturing | |
CN112883550A (en) | Degradation equipment residual life prediction method considering multiple uncertainties | |
CN111460692A (en) | Equipment residual life prediction method and system considering degradation rate mutual influence | |
CN106874634A (en) | Residual life Bayesian forecasting method based on inverse Gauss degradation model | |
CN111523727B (en) | Method for predicting remaining life of battery by considering recovery effect based on uncertain process | |
Wu et al. | Remaining useful life estimation based on a nonlinear Wiener process model with CSN random effects | |
CN105674943A (en) | General multipoint non-linear overall deformation prediction method | |
CN112100574A (en) | Resampling-based AAKR model uncertainty calculation method and system | |
CN114564487B (en) | Meteorological raster data updating method combining forecast prediction | |
Wang et al. | Noise-dependent ranking of prognostics algorithms based on discrepancy without true damage information | |
WO2022161069A1 (en) | Anomaly detection method and apparatus for dynamic control system, and computer-readable medium |
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 |