A method for predicting the remaining life of degraded equipment considering multiple uncertainties

Abstract

To address the lack of sufficient degradation data for long-life and high-reliability equipment, this paper proposes a method for predicting the remaining life of degraded equipment that accounts for multiple uncertainties. This method includes a step-accelerated degradation model based on a nonlinear diffusion process. The advantage of this model is that it requires only a small sample size and a short testing time. The model accounts for multiple uncertainties caused by the inherent characteristics of the degradation model, individual variability, measurement equipment performance, and human bias in the degradation process. The model considers time-varying uncertainty, individual variability, and measurement uncertainty of performance degradation and covariates. To estimate the remaining life of degraded equipment, the present invention derives an analytical approximate solution for the nonlinear diffusion process crossing a predetermined threshold in terms of first arrival time. Combining the maximum likelihood estimation (MLE) method with the simulation-based extrapolation (SIMEX) method, an MLE-SIMEX method is derived for estimating unknown parameters in the model. The effectiveness of the proposed model is demonstrated through simulations and real-world cases. Results demonstrate that this method achieves higher remaining life estimation accuracy and has practical engineering value.

Description

Technical Field

[0001] The present invention belongs to the technical field of reliability engineering, and in particular relates to a method for predicting the remaining service life of degraded equipment taking multiple uncertainties into consideration. Background Art

[0002] With the continuous improvement of current design levels and manufacturing processes, the demand for high-quality products is also increasing. This has led to the emergence of an increasing number of devices with long service lives and high reliability, particularly in the aerospace and military sectors. Traditional remaining life prediction methods are difficult to effectively apply to such devices due to the long testing time and high testing costs. To obtain more degradation data, engineers accelerate device degradation by increasing the severity of the test environment, such as through high temperature, high pressure, and random vibration. Therefore, accelerated degradation testing has become an effective method for collecting product degradation data.

[0003] After obtaining degradation data, random process models can be used to model and analyze the degradation data to estimate the remaining life of the equipment. For complex systems or equipment, nonlinearity and uncertainty are two important factors that affect the accuracy of degradation process modeling. Multiple uncertainties exist in nonlinear accelerated degradation process models, including time-varying uncertainty, individual variability, and measurement uncertainty in performance degradation and covariates. The time-varying uncertainty of equipment degradation is addressed by establishing a degradation model based on random processes; the random effect of the drift coefficient in the degradation model is considered to describe the individual variability of degraded equipment; in addition, due to the limitations of measurement equipment performance and human bias in test data readings, measurement errors will exist in the measurement process of most covariates and degradation performance. Currently, existing methods only consider some of the multiple uncertainties, and there is a lack of relevant research on nonlinear accelerated degradation modeling and remaining life prediction methods that consider multiple uncertainties. Summary of the Invention

[0004] This paper addresses the lack of sufficient degradation data for long-life, high-reliability devices by proposing a step-acceleration degradation model based on a nonlinear diffusion process to estimate the remaining life of the device. This model accounts for the inherent characteristics of degradation models, individual variability, measured device performance, and multiple uncertainties caused by human bias in the degradation process. The model accounts for time-varying uncertainty, individual variability, and measurement uncertainty of performance degradation and covariates.

[0005] The technical solution adopted in the present invention is:

[0006] A method for predicting the remaining life of degraded equipment considering multiple uncertainties includes the following steps:

[0007] Step 1: For degraded equipment, a general nonlinear accelerated degradation process model is established based on the diffusion process, and the standard Brownian motion is used to describe the inherent uncertainty of the random degradation process over time;

[0008] Step 2: Use the Arrhenius model to describe the relationship between the drift coefficient of the degradation model and the accelerated stress. Considering the individual differences of the degraded equipment and the measurement uncertainty of the accelerated stress as a covariate, an accelerated degradation model considering multiple uncertainties is established.

[0009] Step 3: Considering the measurement uncertainty during the performance degradation process, based on the concept of first arrival time, derive the expression for the remaining life distribution of the degraded equipment;

[0010] Step 4: Use the MLE-SIMEX method to estimate the unknown parameters in the accelerated degradation model considering multiple uncertainties, and substitute them into the expression of the remaining life distribution of the degraded equipment to achieve the life prediction of the degraded equipment.

[0011] Preferably, in step 1, the nonlinear accelerated degradation process model is established as follows:

[0012] Assume that X(t) represents the performance degradation of the sample at time t, then the degradation process based on the diffusion process {X(t), t≥0} can be expressed as:

[0013]

[0014] Where μ(t;θ) is the drift coefficient, a nonlinear function of time t, which is used to represent the nonlinear characteristics of the model. The parameter vector θ = (a, b) is an unknown parameter, where a represents the differences between different devices and b represents the common characteristics of similar devices. σ is called the diffusion coefficient, B(t) is the standard Brownian motion, and B(t) ~ N(0, t).

[0015] For formula (1), under the concept of first arrival time, the life of the equipment T can be defined as:

[0016] T=inf{t:X(t)≥ω|X(0)<ω} (2)

[0017] When the randomness of parameters a and σ is not considered, the probability density function of equipment life can be approximated as:

[0018]

[0019] in, ω is the failure threshold.

[0020] Preferably, in step 2, the accelerated degradation model considering multiple uncertainties is established as follows:

[0021] Assume that there are K steps of step-accelerated degradation test data. The test starts with the lowest stress level S1 and the initial time t = 0. Each device will undergo a predetermined stress increase process at a predetermined time. When t = T k (k=1,2,…,K), the corresponding stress is S k , which is expressed as follows:

[0022]

[0023] The drift coefficient a is related to the acceleration stress. The Arrhenius model is used to describe the relationship between a and the acceleration stress S as follows:

[0024] a k =c a exp(-d a / S k ) (5)

[0025] Among them, a k Indicates the drift coefficient a at different stress levels S k The value below, c a is a constant related to the failure mode, d a is a constant related to the activation energy and the Boltzmann constant;

[0026] For the nonlinear degradation model driven by the standard Brownian motion {B(t), t≥0}, the inherent properties of the random process degradation model that change with time can describe the time-varying uncertainty; due to different working conditions, the degradation trajectory will show different degradation rates, and the drift coefficient can be used as a random parameter to describe individual differences. The drift coefficient determines the degradation trajectory of the equipment under different stresses. Considering a k The random effect of , in formula (5), c a can be used as a random parameter and assume μ c and Respectively represent c a The mean and variance of ; therefore, we can get μ ak and Respectively represent the mean and variance of ak, that is:

[0027]

[0028] In actual engineering, to describe the impact of measurement uncertainty, let {Y(t), t≥0} represent the measurement process of the equipment, which is used to describe the relationship between the hidden degradation state and the measured value, as follows:

[0029] Y(t)=X(t)+ε (7)

[0030] Where ε is the random measurement error, assuming represents the variance of ε;

[0031] Usually, the accelerated degradation test is carried out under a predetermined stress level. Due to the existence of measurement error, the equipment actually degrades under an unknown accelerated stress level. The measurement error of the accelerated stress exists in the measurement process of the covariate. The measurement error of the covariate can be expressed as a random variable. Indicates that, and express The variance of ; therefore, the actual accelerated stress level can be expressed as Right now:

[0032]

[0033] In summary, the accelerated degradation model considering multiple uncertainties can be expressed as:

[0034]

[0035] Preferably, in step 3, the expression for the distribution of the remaining life of the degraded equipment is derived as follows:

[0036] Under the first arrival time, considering multiple uncertainties, the equipment life T e It can be defined as:

[0037] T e =inf{t:Y(t)≥ω|Y(0)<ω} (10)

[0038] Correspondingly, the probability density function of equipment life is expressed as

[0039] First, we derive the

[0040] For the degradation process {Y(t), t≥0}, its lifespan is defined as shown in formula (10), then the lifespan T e of It can be expressed as

[0041]

[0042] Among them, B=ω-a·h(t,θ), C=σ 2 t;

[0043] In order to estimate the remaining life of a single device, it is necessary to obtain the degradation data of the device under normal working stress, Y(t k ) indicates that the device is at t k The measured value at the moment, based on the concept of first arrival time, defines t k The remaining life at the moment L k for

[0044] L k =inf{l k >0:Y(l k +t k )≥ω} (12)

[0045] Known t k The degradation state measurement value at the moment is Y(t k ), then the probability density function of the remaining life can be expressed as

[0046]

[0047] Among them, ω tk =ω-Y(t k ),Ω(l k )=h(l k +t k ,θ)-h(t k ,θ);

[0048] The randomness of the drift coefficient a is used to represent the individual differences of the degraded equipment, and a obeys the distribution form of formula (6);

[0049] For the degradation process {Y(t),t≥0}, it is known that And t k The degradation state measurement value at the moment is Y(t k ), then t k The probability density function of the remaining life at time t can be expressed as

[0050]

[0051] in, L=Ω(l k )-l k μ(l k +t k ;θ), M=Ω(l k )=h(l k +t k ,θ)-h(t k ,θ),

[0052] Preferably, in step 4, the model parameters are estimated using the MLE-SIMEX method as follows:

[0053] (1) Simulation steps

[0054] Given a positive integer G and a sequence V={v1,v2,…,v M}, where v1=0, v M is a given positive number, and v2,…,v M-1Obey the uniform distribution U(0,v M ); for the test stress S k , k=1,2,…,K, let u kg , g=1,2,…,G, represents K×G samples that obey the standard normal distribution; therefore, the stress level of the kth covariate can be simulated as follows

[0055]

[0056] Where m = 1, 2, ..., M;

[0057] (2) Maximum likelihood estimation steps

[0058] Assumed test stress There are N test samples, and the degradation state measurement value y(t ijk ) represents the i-th measurement data of the j-th sample under the k-th stress, and the corresponding measurement time is expressed as t ijk , where i = 1, 2, ..., n k , j=1,2,…,N, k=1,2,…,K,n k represents the number of measured values ​​under the kth stress, and the actual degradation state is expressed as x(t ijk );

[0059] Therefore, y(t ijk ) can be expressed as

[0060] y(t ijk )=c aj Λ(t ijk )+σB(t ijk )+ε (16)

[0061] in, c aj represents the drift parameter a j The coefficient of Λ(t ijk ) can be expressed as

[0062]

[0063] For ease of description, let represents the degradation measurement vector of the jth sample, Where (·)′ represents the transpose of the vector, and Y represents all accelerated degradation data, including Y j , j=1,2,…,N;

[0064] According to formula (16), Y j It obeys the multidimensional normal distribution, and its mean and covariance can be expressed as

[0065]

[0066] in, I l is the l-dimensional identity matrix,

[0067] Therefore, for the unknown parameters Θ(g,m)=(μ c ,σ c ,σ,σ ε ,θ,d a ) can be expressed as

[0068]

[0069] in,

[0070] Substitute the log-likelihood function (19) for μ c Find the first-order partial derivative, and we can get

[0071]

[0072] Let formula (20) be 0, we can get

[0073]

[0074] Then according to the estimated Parameter σ c ,σ,σ ε , θ and d a The profile likelihood function can be expressed as

[0075]

[0076] Parameter σ c ,σ,σ ε , θ and d a The maximum likelihood estimate of can be obtained by maximizing the profile likelihood function (22) using the multidimensional search method. and Substituting into formula (21), we can get The maximum likelihood estimate You can get it;

[0077] (3) Extrapolation steps

[0078] If the degradation data of individual equipment under working stress is known, the parameters of the degradation model can be updated according to the Bayesian method to obtain the real-time remaining life prediction value of the equipment;

[0079] Assume that there are κ degradation data under normal working stress, expressed as Y κ =[y1,y2,…,y κ ], where yκ Indicates that the device is at t κ The degradation state measurement value at the moment; assuming that the prior estimation distribution of the drift coefficient a obeys the normal distribution, that is, Then, the joint posterior distribution p(a|Y κ ) can be expressed according to the Bayesian method as

[0080]

[0081] Among them, p(Y κ |a) is the likelihood function, p(a) is the PDF of the prior distribution, Δy g =y g -y g-1 , Δt g =t g -t g-1 ,

[0082] p(Y κ |a) and p(a) are both normally distributed, then the posterior distribution p(a|Y κ ) and the prior distribution should belong to the same distribution family, that is,

[0083]

[0084] in,

[0085]

[0086] Substituting the updated parameters (25) into (16), we can obtain the probability density function of the remaining life of individual equipment and realize the real-time prediction of the remaining life.

[0087] The present invention utilizes a step-acceleration degradation test model based on a nonlinear diffusion process to describe the degradation process of a device. This model simultaneously accounts for time-varying uncertainty, individual variability, and measurement uncertainty in performance degradation and covariates. This model has the advantage of requiring a smaller sample size and shorter testing time. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0089] Figure 1 is the degradation curve of the gyroscope under step acceleration stress;

[0090] Figure 2 is the degradation curve of the gyroscope under working stress;

[0091] Figure 3 Update curves for model parameters;

[0092] Figure 4 is the probability density function of the remaining life of the gyroscope. DETAILED DESCRIPTION

[0093] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments.

[0094] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention as claimed, but rather merely represents selected embodiments of the present invention. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without creative effort shall fall within the scope of protection of the present invention.

[0095] The present invention specifically provides a method for predicting the remaining life of degraded equipment considering multiple uncertainties, comprising the following steps:

[0096] Step 1: For degraded equipment, a general nonlinear accelerated degradation process model is established based on the diffusion process, and the standard Brownian motion is used to describe the inherent uncertainty of the random degradation process over time;

[0097] The nonlinear accelerated degradation process model is established as follows:

[0098] Assume that X(t) represents the performance degradation of the sample at time t, then the degradation process based on the diffusion process {X(t), t≥0} can be expressed as:

[0099]

[0100] Where μ(t; θ) is the drift coefficient, a nonlinear function of time t, which is used to represent the nonlinear characteristics of the model. The parameter vector θ = (a, b) is an unknown parameter, a represents the difference between different devices, and b represents the common characteristics of similar devices. σ is called the diffusion coefficient, B(t) is the standard Brownian motion, and B(t) ~ N(0, t).

[0101] Different functional forms of μ(t; θ) can describe different forms of nonlinear stochastic degradation processes. Clearly, if μ(t; θ) = μ, then Equation (1) transforms into a linear stochastic degradation model, namely, a Wiener process. Therefore, the linear stochastic degradation model is a special case of the nonlinear stochastic degradation model discussed in this paper. Without loss of generality, this paper primarily considers the case where X(0) = 0.

[0102] For formula (1), under the concept of first arrival time, the life of the equipment T can be defined as:

[0103] T=inf{t:X(t)≥ω|X(0)<ω} (2)

[0104] When the randomness of parameters a and σ is not considered, the probability density function of equipment life can be approximated as:

[0105]

[0106] in, ω is the failure threshold.

[0107] Step 2: Use the Arrhenius model to describe the relationship between the drift coefficient of the degradation model and the accelerated stress, while considering the individual differences of the degraded equipment and the measurement uncertainty of the accelerated stress as a covariate, and establish an accelerated degradation model that considers multiple uncertainties.

[0108] The accelerated degradation model considering multiple uncertainties is established as follows:

[0109] In the present invention, it is assumed that there are K-step step-accelerated degradation test data. The test starts with the lowest stress level S1 and the initial time t=0. Each device will undergo a predetermined stress increase process at a predetermined time. When t=T k (k=1,2,…,K), the corresponding stress is S k , which is expressed as follows:

[0110]

[0111] The drift coefficient a is related to the acceleration stress. The Arrhenius model is used to describe the relationship between a and the acceleration stress S as follows:

[0112] a k =c a exp(-d a / S k ) (5)

[0113] Among them, a k Indicates the drift coefficient a at different stress levels S k The value below, c a is a constant related to the failure mode, d ais a constant related to the activation energy and the Boltzmann constant;

[0114] For the nonlinear degradation model driven by the standard Brownian motion {B(t), t≥0}, the inherent properties of the random process degradation model that change over time can describe the time-varying uncertainty; due to different working conditions, the degradation trajectory will show different degradation rates, and the drift coefficient can be used as a random parameter to describe individual differences. The drift coefficient determines the degradation trajectory of the device under different stresses. Considering a k The random effect of , in formula (5), c a can be used as a random parameter and assume μ c and Respectively represent c a The mean and variance of ; therefore, we can get μ ak and Respectively represent a k The mean and variance of , that is:

[0115]

[0116] In practical engineering, it is difficult to accurately measure the degradation state of a device. It is inevitably affected by noise interference and measurement uncertainty caused by non-ideal instruments. To describe the impact of measurement uncertainty, let {Y(t), t≥0} represent the measurement process of the device, which is used to describe the relationship between the hidden degradation state and the measured value, as follows:

[0117] Y(t)=X(t)+ε (7)

[0118] Where ε is the random measurement error, assuming represents the variance of ε;

[0119] Usually, the accelerated degradation test is carried out under a predetermined stress level. The stress is measured and controlled within the specified range during the test. Due to the existence of measurement error, the equipment actually degrades under an unknown accelerated stress level. The measurement error of the accelerated stress exists in the measurement process of the covariate. The measurement error of the covariate can be expressed as a random variable. Indicates that, and express The variance of ; therefore, the actual accelerated stress level can be expressed as Right now:

[0120]

[0121] In summary, the accelerated degradation model considering multiple uncertainties can be expressed as:

[0122]

[0123] Step 3: Considering the measurement uncertainty during the performance degradation process, based on the concept of first arrival time, derive the expression for the remaining life distribution of the degraded equipment;

[0124] Under the first arrival time, considering multiple uncertainties, the equipment life T e It can be defined as:

[0125]

[0126] Correspondingly, the probability density function of equipment life is expressed as

[0127] First, we derive the Introduce the following lemma.

[0128] Lemma 1. There exists Z~N(μ,σ 2 ),and Can get

[0129]

[0130] According to Lemma 1, the total probability formula can be used to obtain the lifespan T in the sense of first arrival time: e Probability density function (PDF) of The specific results are as follows:

[0131] For the degradation process {Y(t), t≥0}, its lifespan is defined as shown in formula (10), then the lifespan T e of It can be expressed as

[0132]

[0133] Among them, B=ω-a·h(t,θ), C=σ 2 t.

[0134] Proof: From formula (7), we can see that Y(t) = X(t) + ε, where Therefore, the lifespan T e That is, the time when {Y(t), t≥0} first reaches the threshold ω-ε. Since the measurement error is random, T can be obtained by the total probability formula e The probability density function (PDF) of .

[0135]

[0136] Let B = ω - a·h(t,θ), C = σ 2 t, Z = ε, according to Lemma 1 we can get

[0137]

[0138] The proof is complete.

[0139] In order to estimate the remaining life of a single device, it is also necessary to obtain the degradation data of the device under normal working stress, Y(t k ) indicates that the device is at t k The measured value at the moment, based on the concept of first arrival time, defines t k The remaining life at the moment L k for

[0140] L k =inf{l k >0:Y(l k +t k )≥ω} (12)

[0141] Known t k The degradation state measurement value at the moment is Y(t k ), then the probability density function of the remaining life can be expressed as

[0142]

[0143] in, Ω(l k )=h(l k +t k ,θ)-h(t k ,θ).

[0144] Proof: Given t k The degradation state measurement value at the moment is Y(t k ), for t≥t k , the degradation process can be expressed as

[0145]

[0146] Let l k =tt k , degradation process {Y(t),t≥t k} can be expressed as

[0147] Y(l k +t k )-Y(t k )=a·(h(l k +t k ,θ)-h(t k ,θ))+σB(l k )

[0148] in,

[0149] Therefore, t kThe remaining life of the equipment at the moment is equivalent to the first arrival time, and the degradation process {Z(l k ),l k ≥0} crosses the threshold time, where Z(l k )=Y(l k +t k )-Y(t k ), and Z(0)=0. Equation (52) can be rewritten as

[0150] Z(l k )=a·(h(l k +t k ,θ)-h(t k ,θ))+σB(l k )

[0151] According to the conclusion of Theorem 1, the expression of the remaining life PDF (50) can be obtained.

[0152] The proof is complete.

[0153] The randomness of the drift coefficient a is used to represent the individual differences of the degraded equipment, and a obeys the distribution form of formula (6);

[0154] For the degradation process {Y(t),t≥0}, it is known that And t k The degradation state measurement value at the moment is Y(t k ), then t k The probability density function of the remaining life at time t can be expressed as

[0155]

[0156] in, L=Ω(l k )-l k μ(l k +t k ;θ), M=Ω(l k )=h(l k +t k ,θ)-h(t k ,θ),

[0157] Proof: From the total probability formula we can get

[0158]

[0159] Let L = Ω(l k )-l k μ(l k +t k ;θ), M=Ω(lk )=h(l k +t k ,θ)-h(t k ,θ), From Lemma 2, we can get

[0160]

[0161] The proof is complete.

[0162] Step 4: Use the MLE-SIMEX method to estimate the unknown parameters in the accelerated degradation model considering multiple uncertainties, and substitute them into the expression of the remaining life distribution of the degraded equipment to achieve the life prediction of the degraded equipment.

[0163] Considering that test stresses such as temperature and humidity are usually measured by common instruments, it is necessary to estimate parameters from independent samples. is feasible and is not related to other unknown parameters. Therefore, the present invention assumes The present invention proposes an improved MLE-SIMEX method to handle the parameter estimation problem of multivariate nonlinear SSADT.

[0164] The model parameters are estimated using the MLE-SIMEX method:

[0165] (1) Simulation steps

[0166] Given a positive integer G and a sequence V={v1,v2,…,v M}, where v1=0, v M is a given positive number, and v2,…,v M-1 Obey the uniform distribution U(0,v M ); for the test stress S k , k=1,2,…,K, let u kg , g=1,2,…,G, represents K×G samples that obey the standard normal distribution; therefore, the stress level of the kth covariate can be simulated as follows

[0167]

[0168] Where m = 1, 2, ..., M;

[0169] (2) Maximum likelihood estimation steps

[0170] Assumed test stress There are N test samples, and the degradation state measurement value y(t ijk ) represents the i-th measurement data of the j-th sample under the k-th stress, and the corresponding measurement time is expressed as t ijk , where i = 1, 2, ..., n k, j=1,2,…,N, k=1,2,…,K,n k represents the number of measured values ​​under the kth stress, and the actual degradation state is expressed as x(t ijk );

[0171] Therefore, y(t ijk ) can be expressed as

[0172] y(t ijk )=c aj Λ(t ijk )+σB(t ijk )+ε (16)

[0173] in, c aj represents the drift parameter a j The coefficient of Λ(t ijk ) can be expressed as

[0174]

[0175] For ease of description, let represents the degradation measurement vector of the jth sample, Where (·)′ represents the transpose of the vector, and Y represents all accelerated degradation data, including Y j , j=1,2,…,N;

[0176] According to formula (16), Y j It obeys the multidimensional normal distribution, and its mean and covariance can be expressed as

[0177]

[0178] in, I l is the l-dimensional identity matrix,

[0179] Therefore, for the unknown parameters Θ(g,m)=(μ c ,σ c ,σ,σ ε ,θ,d a ) can be expressed as

[0180]

[0181] in,

[0182] Substitute the log-likelihood function (19) for μ c Find the first-order partial derivative, and we can get

[0183]

[0184] Let formula (20) be 0, we can get

[0185]

[0186] Then according to the estimated Parameter σ c ,σ,σ ε , θ and d a The profile likelihood function can be expressed as

[0187]

[0188] Parameter σ c ,σ,σ ε , θ and d a The maximum likelihood estimate of can be obtained by maximizing the profile likelihood function (22) using the multidimensional search method. and Substituting into formula (21), we can get The maximum likelihood estimate You can get it;

[0189] (3) Extrapolation steps

[0190] If the degradation data of individual equipment under working stress is known, the parameters of the degradation model can be updated according to the Bayesian method to obtain the real-time remaining life prediction value of the equipment;

[0191] Assume that there are κ degradation data under normal working stress, expressed as Y κ =[y1,y2,…,y κ ], where y κ Indicates that the device is at t κ The degradation state measurement value at the moment; assuming that the prior estimation distribution of the drift coefficient a obeys the normal distribution, that is, Then, the joint posterior distribution p(a|Y κ ) can be expressed according to the Bayesian method as

[0192]

[0193] Among them, p(Y κ |a) is the likelihood function, p(a) is the PDF of the prior distribution, Δy g =y g -y g-1 , Δt g =t g -t g-1 ,

[0194] p(Y κ|a) and p(a) are both normally distributed, then the posterior distribution p(a|Y κ ) and the prior distribution should belong to the same distribution family, that is,

[0195]

[0196] in,

[0197]

[0198] Substituting the updated parameters (25) into (16), we can obtain the probability density function of the remaining life of individual equipment and realize the real-time prediction of the remaining life.

[0199] Example application analysis

[0200] The effectiveness of the method proposed in the present invention is verified by taking the gyroscope in a certain type of inertial navigation system as an example. Considering that temperature is the main factor causing gyroscope drift, the relationship between drift degradation data and temperature stress conforms to the Arrhenius model. In the step stress accelerated degradation test, the threshold is set to 0.6° / h, and the step stress is 50°C, 60°C and 70°C. The gyroscope works continuously in the actual working environment, and the drift degradation data is automatically recorded once an hour, and each stress is recorded 6 times. The test results of the six samples are as follows: Figure 1 shown.

[0201] According to the parameter estimation method proposed in step 4, the estimated values ​​of the unknown parameters in the model can be obtained using the improved MLE-SIMEX method, as shown in Table 1.

[0202] Table 1 Estimated values ​​of degradation model parameters

[0203]

[0204] Gyroscopes of the same type and batch were tested under normal operating stress, and drift degradation data was recorded every 1.5 hours. A total of 14 data were obtained, such as Figure 2 As shown in the figure, it can be seen that after the 13th data, the degradation trajectory reaches the failure threshold for the first time. Therefore, it can be considered that the true life of the sample is L = 20.10 hours.

[0205] In getting the i-th(1≤i≤13) After the measurement values ​​are obtained, the remaining useful life (RUL) of the sample can be predicted. In the remaining useful life (RUL) prediction process, once the new degradation measurement values ​​are available, the model parameters can be updated using the Bayesian method. The specific parameter update process is as follows Figure 3 shown.

[0206] To demonstrate the effectiveness of this method for remaining useful life (RUL) prediction, we present RUL prediction results at different measurement points based on actual observed degradation data, using the initial estimates of the model parameters and their updated values ​​as the degradation data is applied, as shown in Table 1. The RUL prediction results and their corresponding relative errors are shown in Table 2.

[0207] Table 2 Remaining useful life (RUL) prediction values

[0208]

[0209] As can be seen from Table 2, as the model parameters are updated, the predicted RUL is closer to the actual RUL. After the third parameter update, the relative error gradually decreases and stabilizes at around 2%.

[0210] In addition, we Figure 4 The probability density function (PDF) of the remaining useful life (RUL), the predicted average remaining useful life (RUL) and the actual remaining useful life (RUL) at that moment are shown in for comparison. Figure 4 It can be seen that the remaining life predicted by the method of the present invention is very close to the actual remaining life, and the probability density function of the predicted remaining life becomes higher and sharper as the model parameters are updated, indicating that the uncertainty of the prediction results gradually decreases.

[0211] The effectiveness and superiority of the nonlinear accelerated degradation modeling and remaining life prediction method proposed in this invention considering multiple uncertainties are further verified.

[0212] The above description is only used to illustrate the technical solution of the present invention and is not intended to limit it. Other modifications or equivalent substitutions made to the technical solution of the present invention by ordinary technicians in this field should be included in the scope of the claims of the present invention as long as they do not depart from the spirit and scope of the technical solution of the present invention.

Claims

1. A method for predicting the remaining life of degraded equipment considering multiple uncertainties, characterized in that: The following steps are involved: Step 1: For degraded equipment, a general nonlinear accelerated degradation process model is established based on the diffusion process, and the standard Brownian motion is used to describe the inherent uncertainty of the random degradation process over time; Step 2: Use the Arrhenius model to describe the relationship between the drift coefficient of the degradation model and the accelerated stress. Considering the individual differences of the degraded equipment and the measurement uncertainty of the accelerated stress as a covariate, an accelerated degradation model considering multiple uncertainties is established. Step 3: Considering the measurement uncertainty during the performance degradation process, based on the concept of first arrival time, derive the expression for the remaining life distribution of the degraded equipment; Step 4: Use the MLE-SIMEX method to estimate the unknown parameters in the accelerated degradation model considering multiple uncertainties, and substitute them into the expression of the remaining life distribution of the degraded equipment to achieve the life prediction of the degraded equipment; In step 2, the accelerated degradation model considering multiple uncertainties is established as follows: Assume that there are K steps of step-accelerated degradation test data. The test starts with the lowest stress level S1 and the initial time t = 0. Each device will undergo a predetermined stress increase process at a predetermined time. When t = T k (k=1,2,…,K), the corresponding stress is S k , which is expressed as follows: The drift coefficient a is related to the acceleration stress. The Arrhenius model is used to describe the relationship between a and the acceleration stress S as follows: a k =c a exp(-d a / S k ) (5) Among them, a k Indicates the drift coefficient a at different stress levels S k The value below, c a is a constant related to the failure mode, d a is a constant related to the activation energy and the Boltzmann constant; For the nonlinear degradation model driven by the standard Brownian motion {B(t), t≥0}, the inherent properties of the random process degradation model that change with time can describe the time-varying uncertainty; due to different working conditions, the degradation trajectory will show different degradation rates, and the drift coefficient can be used as a random parameter to describe individual differences. The drift coefficient determines the degradation trajectory of the equipment under different stresses. Considering a k The random effect of , in formula (5), c a can be used as a random parameter and assume μ c and Respectively represent c a The mean and variance of ; therefore, we can get μ ak and Respectively represent a k The mean and variance of , that is: In actual engineering, to describe the impact of measurement uncertainty, let {Y(t), t≥0} represent the measurement process of the equipment, which is used to describe the relationship between the hidden degradation state and the measured value, as follows: Y(t)=X(t)+ε (7) Where ε is the random measurement error, assuming represents the variance of ε; Usually, the accelerated degradation test is carried out under a predetermined stress level. Due to the existence of measurement error, the equipment actually degrades under an unknown accelerated stress level. The measurement error of the accelerated stress exists in the measurement process of the covariate. The measurement error of the covariate can be expressed as a random variable. Indicates that, and express The variance of , therefore, the actual accelerated stress level can be expressed as Right now: In summary, the accelerated degradation model considering multiple uncertainties can be expressed as: In step 3, the remaining life distribution expression of degraded equipment is derived as follows: Under the first arrival time, considering multiple uncertainties, the equipment life T e It can be defined as: T e = inf{t : Y(t)≥ω|Y(0)<ω} (10) Correspondingly, the probability density function of equipment life is expressed as f Te (t); First, f is derived without considering the random effect of the drift coefficient Te (t) For the degradation process {Y(t), t≥0}, its lifespan is defined as shown in formula (10), then the lifespan T e f Te|a (t) can be expressed as Where, B=ω-a·h(t,θ), C=σ 2 t; In order to estimate the remaining life of a single device, it is necessary to obtain the degradation data of the device under normal working stress, Y(t k ) indicates that the device is at t k The measured value at the moment, based on the concept of first arrival time, defines t k The remaining life at the moment L k for L k = inf{l k >0 : Y(l k +t k )≥ω} (12) Known t k The degradation state measurement value at the moment is Y(t k ), then the probability density function of the remaining life can be expressed as Among them, oh tk =ω-Y(t k ),Ω(l k )=h(l k +t k ,θ)-h(t k ,i); The randomness of the drift coefficient a is used to represent the individual differences of the degraded equipment, and a obeys the distribution form of formula (6); For the degradation process {Y(t),t≥0}, it is known that And t k The degradation state measurement value at the moment is Y(t k ), then t k The probability density function of the remaining life at time t can be expressed as Among them, oh tk =ω-Y(t k ), L=Ω(l k )-l k μ(l k +t k ;θ), M=Ω(l k )=h(l k +t k ,θ)-h(t k ,i), 2. The method for predicting the remaining life of degraded equipment considering multiple uncertainties according to claim 1, characterized in that: In step 1, the nonlinear accelerated degradation process model is established as follows: Assume that X(t) represents the performance degradation of the sample at time t, then the degradation process based on the diffusion process {X(t), t≥0} can be expressed as: Where μ(t; θ) is the drift coefficient, which is a nonlinear function of time t and is used to represent the nonlinear characteristics of the model. The parameter vector θ = (a, b) is an unknown parameter, where a represents the differences between different devices and b represents the common characteristics of similar devices. σ B is called the diffusion coefficient, B(t) is the standard Brownian motion, and B(t) ~ N(0,t); For formula (1), under the concept of first arrival time, the life of the equipment T can be defined as: T= inf{t : X(t)≥ω|X(0)<ω} (2) When the parameters a and σ are not considered B When the randomness is , the probability density function of equipment life can be approximated as: in, ω is the failure threshold.

3. The method for predicting the remaining life of degraded equipment considering multiple uncertainties according to claim 1, characterized in that: In step 4, the model parameters are estimated using the MLE-SIMEX method as follows: (1) Simulation steps Given a positive integer G and a sequence V={v1,v2,…,v M }, where v1=0, v M is a given positive number, and v2,…,v M-1 Obey the uniform distribution U(0,v M ); for the test stress S k , k=1,2,…,K, let u kg , g=1,2,…,G, represents K×G samples that obey the standard normal distribution; therefore, the stress level of the kth covariate can be simulated as follows Where m = 1, 2, ..., M; (2) Maximum likelihood estimation steps Assumed test stress There are N test samples, and the degradation state measurement value y(t ijk ) represents the i-th measurement data of the j-th sample under the k-th stress, and the corresponding measurement time is expressed as t ijk , where i = 1, 2, ..., n k , j=1,2,…,N, k=1,2,…,K,n k represents the number of measured values ​​under the kth stress, and the actual degradation state is expressed as x(t ijk ); Therefore, y(t ijk ) can be expressed as y(t ijk )=c aj Λ(t ijk )+σB(t ijk )+ε (16) Among them, t nkjk =T k , c aj represents the drift parameter a j The coefficient of Λ(t ijk ) can be expressed as For ease of description, let represents the degradation measurement vector of the jth sample, Where (·)′ represents the transpose of the vector, and Y represents all accelerated degradation data, including Y j , j=1,2,…,N; According to formula (16), Y j It obeys the multidimensional normal distribution, and its mean and covariance can be expressed as in, I l is the l-dimensional identity matrix, Therefore, for the unknown parameters Θ(g,m)=(μ c ,σ c ,σ,σ ε ,θ,d a ) can be expressed as in, Substitute the log-likelihood function (19) for μ c Find the first-order partial derivative, and we can get Let formula (20) be 0, we can get Then according to the estimated Parameter σ c ,σ,σ ε , θ and d a The profile likelihood function can be expressed as Parameter σ c ,σ,σ ε , θ and d a The maximum likelihood estimate of can be obtained by maximizing the profile likelihood function (22) using the multidimensional search method. and Substituting into formula (21), we can get The maximum likelihood estimate You can get it; (3) Extrapolation steps If the degradation data of individual equipment under working stress is known, the parameters of the degradation model can be updated according to the Bayesian method to obtain the real-time remaining life prediction value of the equipment; Assume that there are κ degradation data under normal working stress, expressed as Y κ =[y1,y2,…,y κ ], where y κ Indicates that the device is at t κ The degradation state measurement value at the moment; assuming that the prior estimate distribution of the drift coefficient a obeys the normal distribution, that is, Then, the joint posterior distribution p(a|Y κ ) can be expressed according to the Bayesian method as Among them, p(Y κ |a) is the likelihood function, p(a) is the PDF of the prior distribution, Δy g =y g -y g-1 , Δt g =t g -t g-1 , p(Y κ |a) and p(a) are both normally distributed, then the posterior distribution p(a|Y κ ) and the prior distribution should belong to the same distribution family, that is, in, Substituting the updated parameter (25) into (14), the probability density function of the remaining life of the individual equipment can be obtained, and the real-time prediction of the remaining life can be achieved.

Images