Method for predicting residual life of mechanical system by considering dynamic environment influence

Abstract

The invention belongs to the technical field of system life prediction, and discloses a mechanical system residual life prediction method considering dynamic environment influence, and the specific technical scheme is as follows: step 1, system degradation modeling: firstly, according to received monitoring data representing an operation state and a degradation state, carrying out correlation analysis and processing on two types of obtained states, constructing a system dynamics evolution model and a part degradation model; step 2, according to the established system degradation model, estimating an observable degradation state and an observable operation state by using a dual-time scale Kalman filtering algorithm and a Rauch Tung Striebel filtering algorithm, and estimating model parameters by using a two-step expectation condition maximization algorithm; 3, defining the residual life of the system according to the first arrival time on the basis of considering the dual-time scale influence, the measurement error and the bidirectional correlation of the two states; the system life prediction precision is high, the effectiveness is good, and the system life prediction method has better general adaptability.

Description

Technical Field

[0001] The present invention belongs to the technical field of system life prediction, and particularly relates to a method for predicting the remaining life of a mechanical system considering the influence of a dynamic environment. Background Technique

[0002] With the continuous evolution of industrial Internet and industrial cloud storage / cloud computing technologies, the industrial big data platform has been continuously promoted in the fields of energy, equipment, manufacturing, etc., providing strong data support for the prediction research of systems, setting off an unprecedented research boom in prediction and predictive maintenance decision-making, and driving the extensive application demand for predictive maintenance. Prognostics and Health Management (PHM) aims to apply a large amount of online monitoring data, through system degradation modeling and reliability analysis, predict its future state or remaining life, and based on this, make predictive maintenance decisions and management, so as to realize remote real-time intelligent operation and maintenance services and health management for the entire life cycle of the system. With the continuous development of industrial big data technology, PHM has become a new research hotspot in the field of intelligent operation and maintenance and health management of current complex systems.

[0003] As the foundation and core of PHM technology, the prediction of Remaining Useful Life (RUL) has always been the focus of PHM research. At present, most of the research on system degradation modeling and remaining life prediction implicitly assumes that the operating state of the system is stable, that is, the operating state of the system has no influence on its degradation process. In fact, when the system operates in a dynamic environment, its operating state will change with the changes of the environment, load or operating mode, which brings challenges to the system degradation modeling and remaining life prediction. For example, the rotational speed of a fan changes continuously under the influence of wind speed changes and wind direction fluctuations, resulting in the fan failing prematurely before reaching its expected life. The vibration level and degradation rate of rotating components such as gears, bearings and pumps often increase with the increase of rotational speed. Therefore, a degradation model that ignores the system state cannot accurately predict the RUL of the system.

[0004] Currently, research on RUL prediction considering the impact of system state on degradation has been conducted in complex systems. Some current literature has proposed that as the charging time of the battery increases, the ion concentration in the battery decreases, and thus the battery capacity gradually decays. For a high-temperature furnace, as the temperature increases, the furnace wall is corroded by molten iron, generating new substances, resulting in an increase in the degradation rate of the high-temperature furnace. For an inertial platform system, due to the increasing disturbing torque generated by the motor shaft and the gyro output shaft, and the system being subjected to continuous loads as well as cumulative damage from internal stress and the external environment, mechanical or electrical parameters gradually change over time. This makes the degradation process of the actuator output efficiency gradual and random. In these cases, the state of the system is also crucial for its degradation process, and conversely, the degradation process of the system will also affect the overall operation of the system.

[0005] These studies have focused on the impact of component operating states on system degradation behavior, further constructing system evolution models and degradation models, and predicting the remaining life of the system. However, in the actual application process, the system degradation process is much slower compared to its state changes. Although the above studies have considered the problem of life estimation related to the system operating state and degradation state, there are no reports in the literature on the stochastic degradation modeling and RUL prediction under the above-mentioned related double time scales. From the perspective of the system operating state, it is a rapidly changing process, and from the perspective of the system degradation state, it is a slowly changing process. At the same time, these two types of states affect each other. In the research process of RUL prediction for degraded components, this problem is considered a double time scale. For example, the change in the operating speed of a bearing is a rapidly changing process, but the degradation, corrosion, and wear of the bearing are slowly changing processes. The ill-conditioned problem caused by the interaction between slow and fast dynamic modes needs to be solved.

[0006] To solve the problem of different time scales between the system state and the operating state, a time scale separation parameter is introduced to unify the system state change and the degradation state to the same time scale for measurement, and a state space model is used to construct the system model. After introducing the time separation scale parameter in this model, although the problem of different time scales between the degradation state is solved, due to the inclusion of extremely small parameters in the model, ill-conditioned problems may occur during the state estimation process. To address this problem, the phase space reconstruction method is adopted in some literature. The fast system operating state and the slow degradation state are transformed to an intermediate time scale for processing, which can avoid ill-conditioned problems during the state estimation process. The selection of this intermediate time scale is crucial. For a fast-changing system, an inappropriate intermediate time scale cannot accurately reflect the dynamic characteristics of the system, thereby affecting the accuracy of the correlation modeling between the system state and the degradation state. Therefore, it is crucial to consider the remaining life prediction related to the system operating state and degradation state under the influence of double time scales. Summary of the Invention

[0007] To solve the technical problems existing in the prior art, the present invention provides a method for predicting the remaining useful life of a mechanical system considering the influence of a dynamic environment, with high prediction accuracy.

[0008] To achieve the above object, the technical solution adopted by the present invention is as follows: A method for predicting the remaining useful life of a mechanical system considering the influence of a dynamic environment, the specific steps are as follows:

[0009] Step 1. System degradation modeling: First, based on the received monitoring data representing the operating state and degradation state, perform correlation analysis and processing on the obtained two types of states, and construct a system dynamics evolution model and a component degradation model;

[0010] Step 2. According to the established system degradation model, use the dual-time-scale Kalman filter algorithm and the Rauch-Tung-Striebel (RTS) filter algorithm to estimate the observable degradation state and operating state, and use the two-step expectation conditional maximization algorithm to estimate the model parameters;

[0011] Step 3. Based on considering the influence of the dual-time scale, measurement error, and the bidirectional correlation between the two types of states, define the remaining useful life of the system according to the first passage time.

[0012] System degradation modeling:

[0013] First, based on the received monitoring data representing the operating state and degradation state, perform correlation analysis and processing on the obtained two types of states, and construct its system dynamics evolution model and component degradation model;

[0014] (1). System dynamics evolution model:

[0015] X(t) = A(Z(T))X(t - 1) + B(Z(T))U(t - 1) + ω(t) (1)

[0016] Y(t) = C(Z(T))X(t) + ξ(t) (2)

[0017] Wherein, X(t) represents the system operating state, A(Z(T)) is the state matrix related to system degradation, B(Z(T)) is the input matrix related to system degradation, C(Z(T)) is the output matrix related to system degradation, Y(t) is the measured value of the operating state data, t is the time scale of the system operating state, T is the time scale of the system degradation state, U(t) is the external input quantity, ω(t) is the system noise and ω(t) ~ (0, Q), and ξ(t) is the measurement noise of the operating data and ξ(t) ~ (0, R).

[0018] (2). System degradation model:

[0019] For complex systems, degradation is mostly non - linear and non - monotonic. Therefore, the Wiener process is used to model the system degradation, that is:

[0020]

[0021] where \(Z(T)\) is the degradation process of the system, \(\sigma\) B \(B(T)\) is the diffusion process, \(\sigma\) B is the diffusion coefficient, \(B(T)\) is the Brownian motion, is the system degradation rate, \(\Omega\) is the parameter of the degradation rate function, and the average operating state is expressed as:

[0022]

[0023] where \(D\) nx1 reflects the influence of different components of the system state deviating from the steady - state value \(x\) ε on the degradation, which can be determined in advance according to the specific system, \(()\) T represents the transpose of the matrix.

[0024] To observe the system degradation state, an observation equation for its degradation state needs to be constructed:

[0025] \(V(T)=GZ(T)+\Lambda(T)\ (5)\)

[0026] where \(V(T)\) represents the measured value of the degradation data, \(Z(T)\) represents the degradation state of the system, \(\Lambda(T)\) represents the measurement noise of the degradation data, that is, \(\Lambda(T)\sim MVN(0,\sigma\) R ) and \(\sigma\) R represents the variance of the measurement noise, and \(G\) represents the measurement coefficient of the degradation state.

[0027] In an actual system, the system operating state changes continuously multiple times within the unit degradation time. Through analysis, the influence of the degradation state on the system operating state is at a certain point, and the influence of the system operating state on the degradation is over a certain period of time, referring to the comprehensive effect of the system state within that time period. Thus, the influence of the system operating state on the degradation rate can be reflected by the average effect of the system state change within the unit degradation time.

[0028] Without loss of generality, the drift term is a degradation function related to the system operating state and the degradation rate, which can have different forms, such as linear function, polynomial function, exponential function, and power function, etc. Specifically, it needs to be determined according to the physical characteristics of the system. \(\Omega\) represents the parameter of the degradation rate function, and this function is also related to the operating environment, individual differences, and load of the system.

[0029] In practical applications, the system operates in a dynamic environment, and its operating state changes continuously over time. The deviation λ(t) between any operating state X(t) at time t and the steady-state value x ω is defined as λ(t) = X(t) - x ω and is called the offset. The change in λ(t) will affect the degradation rate η(T) of the system. When λ(t) = 0, the operating state of the system is the steady-state value X(t) ≡ x ω . At this time, the degradation rate of the system is a constant.

[0030] Degradation state and parameter estimation:

[0031] According to the established system degradation model, the dual-time-scale Kalman filtering algorithm and the RTS filtering algorithm are used to estimate the partially observable degradation state and operating state, and the two-step expectation conditional maximization algorithm is used to estimate the model parameters θ = [Ω, σ B , σ R , G]. First, consider that in the small time scale, the system state changes, but the degradation state remains unchanged; in the large time scale, the corresponding degradation state changes and then affects the change of the system operating state.

[0032] First, given the initial values Z(T0), X(0), U(0), A(Z(T0)), B(Z(T0)), and C(Z(T0)), the filtering algorithm and the EM algorithm are used to simultaneously estimate the degradation state Z(T k-1 ) and the parameter θ k-1 of the system per unit degradation time T k-1 . According to the evolution model, calculate A(Z(T k-1 ), B(Z(T k-1 ), C(Z(T k-1 ). When the system state X(t) changes over time t ∈ [T k-1 , T k , estimate the system operating state X(t i-1 ) at time t i-1 and determine whether the system operating time reaches t p . If not, enter the next sampling moment until the end. At this time, Z(T k-1 ) remains unchanged, and the corresponding A(Z(T k-1 ), B(Z(T k-1 ), C(Z(T k-1 ) remain unchanged.

[0033] Secondly, when the system degradation time enters the next sampling moment, estimate the degradation state Z(T k-1 ) and the parameter θ in the degradation model according to the average state of the system obtained in the previous degradation period T k ) and parameter θ k, the corresponding degradation affects the system, making A(Z(T k-1 ))、B(Z(T k-1 )) and C(Z(T k-1 )) become A(Z(T k ))、B(Z(T k )) and C(Z(T k ))。

[0034] Determine whether the system degradation sampling time has reached T e . If not, enter the next sampling time and continue to execute the above operations until the sampling ends.

[0035] Degradation state estimation:

[0036] Let the sampling period of the system state be Δt, and the period of the degradation state be ΔT with ΔT >> Δt. Denote Assume that at the end of time T k , the observed data of the system operating state at t ∈ [T k-1 , T k is The corresponding system degradation state set is Z 0:k = {Z0, Z1..., Z k}; the set of measured values of the degradation data is V 0:k = {V0, V1..., V k}. Assume that in t ∈ (T, T + ΔT), the system degradation state Z(T) remains unchanged, and design a filter to update the system state and the degradation state; in equations (1)-(5), assume that the structure and parameters of the system are determined, then the actual operating state X i at t ∈ [T k-1 , T k is estimated using the system operating data Y 0:ρ , and the actual degradation state Z k of the system is estimated through the system degradation data V 0:k . For the convenience of calculation, represents a nonlinear function, b is a nonlinear parameter, and the specific double-time-scale Kalman filtering algorithm is as follows:

[0037] 1), Small-time-scale one-step prediction, including operating state prediction and covariance prediction:

[0038]

[0039] One-step update, including operating state update and covariance update:

[0040]

[0041] Among them, represents the estimated value of the operating state at time i, represents the operating state gain at time i, represents the covariance matrix of the predicted operating state X at time i, represents the covariance matrix of the estimated operating state X at time i, Y i represents the measurement data of the operating state at time i.

[0042] Finally, it is obtained that and the unit degradation period T k the average expected operating state of and covariance are:

[0043]

[0044] Among them, tr() represents the trace of the matrix, represents the summation function.

[0045] 2) One-step prediction on a large time scale, including degradation state prediction and covariance prediction:

[0046]

[0047] The second step of updating, including degradation state update and covariance update:

[0048]

[0049] Among them, represents the estimated value of the degradation state at time k, represents the degradation state gain at time k, represents the covariance matrix of the predicted degradation state Z at time k, represents the covariance matrix of the estimated degradation state Z at time k, V k represents the measurement data of the degradation state at time k. Finally, it is obtained that

[0050] Parameter estimation:

[0051] If the observed data V 0:k of the system degradation state is known, based on Equation (3), the system performance degradation state Z k and the observed data V 0:k can be obtained, and the joint log-likelihood function is:

[0052] ln f(Z k , V 0:k ; θ) = ln f(Z k ; θ) + ln f(V 0:k |Z k; θ) (18)

[0053] where, ln f(Z k ; θ) represents the probability density function of the system degradation state, ln f(V 0:k |Z k ; θ) represents the probability density function of the measurement data; ln() represents the logarithmic function.

[0054] Based on the assumptions of the model in Equation (18), we obtain:[[]]

[0055]

[0056] where,[[]] represents the product of functions,[[]] represents the summation of functions from j = 1 to j = k; τ represents the time interval and τ = T k -T k-1 , Z j-1 and Z j represent the values of the state Z at times j-1 and j respectively. Similarly,[[]]

[0057]

[0058] Substituting Equations (19) and (20) into the equation, we can obtain:[[]]

[0059]

[0060] Based on the l-th estimation result θ l , by solving the expectation of the likelihood function l n(Z k , V 0:k ; θ) with respect to the hidden state Z, we can obtain E[ln f(Z k , V 0:k ; θ)]. When maximizing the logarithmic likelihood function of Equation (21), since there are many hidden variables in Equation (21), it cannot be directly maximized.[[]]

[0061] To calculate the conditional expectation of the logarithmic likelihood function, the expectation is obtained using the Kalman smoothing filter algorithm. First, we get:[[]]

[0062]

[0063] where, j = k..., 1, 0,[[]] and can be obtained through the RTS backward filtering algorithm.[[]]

[0064]

[0065] where,[[]] is the mean of the backward filtering,[[]] is the variance of the inverse filtering, is the cross-covariance matrix of the inverse filtering, M j is the filtering gain of the inverse filtering.

[0066] Therefore, the present invention substitutes equations (25)-(28) into equation (21) to obtain:

[0067]

[0068] Further arranging equation (29) gives:

[0069]

[0070] wherein, represents the mean value of the operating state estimated by the inverse filtering at time j.

[0071] Equation (31) can be used to implement the iterative estimation of all parameters:

[0072]

[0073] wherein, argmax() represents obtaining the maximum value, represents the value of the unknown parameter at the (l + 1)-th step of iteration at time k.

[0074] Let After estimating by setting the partial derivative of the function to 0, we have:

[0075]

[0076] For the parameter b, taking the partial derivative according to equation (30) gives:

[0077]

[0078] By using the fminsearch function in MATLAB, the parameter estimation value is obtained by solving the non-linear equation

[0079] Remaining useful life prediction:

[0080] On the basis of considering the influence of the double time scale, measurement error, and two-way correlation of the two types of states, the remaining useful life L of the system is defined according to the first passage time k :

[0081]

[0082] wherein, inf{} represents the infimum, l k is the RUL of the system at time T k Z(T k + l k ) represents Tk +l k The degradation state of the time system, where w is the system failure threshold, and the probability density function of the remaining life of the system is:

[0083]

[0084] where exp() represents the exponential function.

[0085] To predict the remaining life of the system in the implicit degradation state, define and Through the entire monitoring sequence V k at time T 1:k =(V0, V1,..., V k ) to obtain the conditional expectation k and variance of state Z That is Further, the PDF h(Z k |V k |) of the degradation state Z 0:k is:

[0086]

[0087] By the total probability formula, the k at the current time T, the of the RUL of the system in the implicit degradation state is:

[0088]

[0089] where represents the integral function.

[0090] Substitute equations (38) and (39) into equation (40) for calculation and simplification, and further organize to obtain:

[0091]

[0092] On this basis, considering the system operating state distribution That is Obtain the PDF f k (l k |V 0:k ) of the RUL of the system under the influence of the double-time scale as:

[0093]

[0094] Substitute equation (41) into (42) to further obtain:

[0095]

[0096] Among them,

[0097] Based on the existing technology, the present invention introduces factors such as the two-way correlation between the system operation state and the degradation state, the time-varying degradation rate, and the hidden state, and constructs a stochastic degradation model and remaining life prediction of the system under the influence of a two-time-scale, which will provide a decision-making basis for the research on RUL prediction and maintenance of mechanical systems considering the influence of dynamic environments. Description of the Drawings

[0098] Figure 1 It is a relationship diagram of the system operation state and the degradation state under the influence of a two-time-scale. Figure 1 (a) is the degradation state relationship diagram. Figure 1 (b) is the operation state relationship diagram.

[0099] Figure 2 It is a flow chart of the filtering algorithm based on the influence of a two-time-scale.

[0100] Figure 3 It is a curve diagram of the system operation state estimation.

[0101] Figure 4 It is a curve diagram of the system degradation state estimation.

[0102] Figure 5 It is a curve diagram of the system parameter estimation result.

[0103] Figure 6 It is a degradation trajectory diagram considering the system state and not considering the system state.

[0104] Figure 7 It is a remaining life prediction result diagram of different methods. Figure 7 (a) is the prediction result diagram of the system degradation model without considering the influence of the operation state. Figure 7 (b) is the prediction result diagram of the system degradation model without considering the influence of the time-scale. Figure 7 (c) is the prediction result diagram of the present invention.

[0105] Figure 8 It is a comparison result diagram of the RUL means of different prediction methods.

[0106] Figure 9 It is a comparison diagram of the MSE of different prediction methods.

[0107] Figure 10 It is a curve diagram of the armature current operation state.

[0108] Figure 11 It is a curve diagram of the torque coefficient degradation state.

[0109] Figure 12Prediction result graph of the method proposed in the present invention. Detailed implementation manners

[0110] In order to make the technical problems, technical solutions and beneficial effects to be solved by the present invention more clear and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0111] The method for predicting the remaining useful life of a mechanical system considering the influence of a dynamic environment is as follows:

[0112] System degradation modeling:

[0113] First, based on the monitoring data that can be received to characterize the operating state and degradation state, the correlation analysis and processing of the two types of states are carried out, and its system dynamics evolution model and component degradation model are constructed;

[0114] (1) System dynamics evolution model:

[0115] X(t) = A(Z(T))X(t - 1) + B(Z(T))U(t - 1) + ω(t) (1)

[0116] Y(t) = C(Z(T))X(t) + ξ(t) (2)

[0117] Among them, X(t) represents the system operating state, A(Z(T)) is the state matrix related to system degradation, B(Z(T)) is the input matrix related to system degradation, C(Z(T)) is the output matrix related to system degradation, Y(t) is the measured value of the operating state data, t is the time scale of the system operating state, T is the time scale of the system degradation state, U(t) is the external input quantity, ω(t) is the system noise and ω(t) ∼ (0, Q), and ξ(t) is the measurement noise of the operating data and ξ(t) ∼ (0, R).

[0118] (2) System degradation model:

[0119] For the degradation of complex systems, it is mostly non-linear and non-monotonic. Therefore, the Wiener process is used to model the degradation of the system, that is:

[0120]

[0121] Among them, Z(T) is the degradation process of the system, σ B B(T) is the diffusion process, σ B is the diffusion coefficient, B(T) is the Brownian motion, is the system degradation rate, Ω is the parameter of the degradation rate function, and the average operating state is expressed as:

[0122]

[0123] Among them, D nx1 reflects the deviation of different components of the system state from the steady-state value x ε The impact on degradation can be determined in advance according to the specific system, () T represents the transpose of the matrix.

[0124] In order to observe the degradation state of the system, it is necessary to construct an observation equation for its degradation state:

[0125] V(T) = GZ(T) + Λ(T) (5)

[0126] Among them, V(T) represents the measured value of the degradation data, Z(T) represents the degradation state of the system, and Λ(T) represents the measurement noise of the degradation data, that is, Λ(T) ~ MVN(0, σ R ) and σ R represents the variance of the measurement noise, and G is the measurement coefficient of the degradation state.

[0127] In an actual system, the operating state of the system will change continuously multiple times within a unit degradation time. Through analysis, the impact of the degradation state on the operating state of the system is at a certain point, and the impact of the operating state of the system on degradation is over a certain period of time, referring to the comprehensive effect of the system state within that time period. Therefore, the impact of the operating state of the system on the degradation rate can be reflected by the average effect of the change in the system state within a unit degradation time.

[0128] Without loss of generality, the drift term is a degradation function related to the operating state of the system and the degradation rate, which can have different forms, such as linear function, polynomial function, exponential function, and power function, etc. Specifically, it needs to be determined according to the physical characteristics of the system. Ω represents the parameter of the degradation rate function, and this function is also related to the operating environment, individual differences, and load of the system.

[0129] In practical applications, the system operates in a dynamic environment, and its operating state changes continuously with time. The deviation λ(t) between the arbitrary operating state X(t) at time t and the steady-state value x ω is defined as λ(t) = X(t) - x ω and is called the offset. The change of λ(t) will affect the degradation rate η(T) of the system. When λ(t) = 0, the operating state of the system is the steady-state value X(t) ≡ x ω . At this time, the degradation rate of the system is a constant.

[0130] According to the correlation analysis of the system dynamics evolution model and the degradation model, a state space model of a single-component system is established simultaneously. Based on the relationship between the operating state and the degradation state of the system under the influence of double time scales as Figure 1 shown.

[0131] Degradation state and parameter estimation:

[0132] According to the established system degradation model, the dual-time-scale Kalman filtering algorithm and the RTS filtering algorithm are used to estimate the partially observable degradation state and operating state, and the two-step expectation conditional maximization algorithm is used to estimate the model parameters θ = [Ω, σ B , σ R , G]. First, consider that in the small time scale, the system state changes, but the degradation state remains unchanged; in the large time scale, the corresponding degradation state changes and then affects the change of the system operating state. Specifically, the filtering algorithm process based on the dual-time-scale influence is as Figure 2 shown.

[0133] First, given the initial values Z(T0), X(0), U(0), A(Z(T0)), B(Z(T0)) and C(Z(T0)), the filtering algorithm and the EM algorithm are used to simultaneously estimate the degradation state Z(T k-1 ) and the parameter θ k-1 of the system unit degradation time T k-1 . According to the evolution model, calculate A(Z(T k-1 ), B(Z(T k-1 ), C(Z(T k-1 ). When the system state X(t) changes with time t ∈ [T k-1 , T k , estimate the system operating state X(t i-1 ) at time t i-1 , and judge whether the system operating time reaches t p . If not, enter the next sampling time until the end. At this time, Z(T k-1 ) remains unchanged, and the corresponding A(Z(T k-1 ), B(Z(T k-1 ), C(Z(T k-1 ) remain unchanged.

[0134] Secondly, when the system degradation time enters the next sampling time, according to the system average state X(T k-1 ) obtained in the previous degradation period T k-1 , estimate the degradation state Z(T k ) and the parameter θ k in the degradation model. The corresponding degradation affects the system so that A(Z(T k-1 ), B(Z(T k-1 ), and C(Z(T k-1 ) become A(Z(T k ), B(Z(T k ), and C(Z(T k ).

[0135] Determine whether the system degradation sampling time reaches T e , if not, enter the next sampling time and continue to execute the above operations until the sampling ends.

[0136] Degradation state estimation:

[0137] Let the sampling period of the system state be Δt, the period of the degradation state be ΔT and ΔT >> Δt, and denote Suppose at T k moment, the observed data of the system operating state at t ∈ [T k-1 , T k is the observed value corresponding to The corresponding system degradation state set is Z 0:k = {Z0, Z1..., Z k}; the set of measured values of the degradation data is V 0:k = {V0, V1..., V k}. Suppose that in t ∈ (T, T + ΔT), the system degradation state Z(T) remains unchanged, and a filter is designed to update the system state and the degradation state; in equations (1)-(5), it is assumed that the structure and parameters of the system are determined, then the actual operating state X i at t ∈ [T k-1 , T k is estimated using the system operating data Y 0:ρ , and the actual degradation state Z k of the system is estimated through the system degradation data V 0:k . For the convenience of calculation, represents a nonlinear function, b is a nonlinear parameter, and the specific two-time-scale Kalman filtering algorithm is as follows:

[0138] 1), One-step prediction on the small time scale, including operating state prediction and covariance prediction:

[0139]

[0140] One-step update, including operating state update and covariance update:

[0141]

[0142] Among them, represents the estimated value of the operating state at the i-th moment, represents the operating state gain at the i-th moment, represents the covariance matrix of the prediction of the operating state X at the i-th moment, represents the covariance matrix of the estimation of the operating state X at the i-th moment, Y iMeasurement data representing the operating state at time i.

[0143] Finally obtain And obtain the unit degradation period T k The expected value of the average operating state And covariance Is:

[0144]

[0145] Where tr() represents the trace of the matrix, Represents the summation function.

[0146] 2), One-step prediction on a large time scale, including degradation state prediction and covariance prediction:

[0147]

[0148] The second step of updating, including degradation state update and covariance update:

[0149]

[0150] Where Represents the estimated value of the degradation state at time k, Represents the degradation state gain at time k, Represents the covariance matrix of the predicted degradation state Z at time k, Represents the covariance matrix of the estimated degradation state Z at time k, V k Represents the measurement data of the degradation state at time k, and finally obtain

[0151] Parameter estimation:

[0152] If the observed data V of the system degradation state is known 0:k , Based on Equation (3), the system performance degradation state Z can be obtained k And the observed data V 0:k Of the joint log-likelihood function:

[0153] ln f(Z k ,V 0:k ; θ) = ln f(Z k ; θ) + ln f(V 0:k |Z k ; θ) (18)

[0154] Where, ln f(Z k ; θ) represents the probability density function of the system degradation state, ln f(V 0:k |Z k ; θ) represents the probability density function of the measurement data; ln() represents the logarithmic function.

[0155] It is obtained according to the assumptions of the model in Equation (18):

[0156]

[0157] where denotes the product of functions, denotes the summation of functions from j = 1 to j = k; τ represents the time interval and τ = T k -T k-1 Z j-1 and Z j represent the values of state Z at times j - 1 and j respectively. Similarly,

[0158]

[0159] Substituting Equation (19) and Equation (20) into the equation gives:

[0160]

[0161] Based on the l-th estimated result θ l , by solving the expectation of the likelihood function ln(Z k , V 0:k ; θ) with respect to the hidden state Z, we can obtain E[ln f(Z k , V 0:k ; θ)]. When maximizing the log-likelihood function of Equation (21), since there are many hidden variables in Equation (21), it cannot be directly maximized.

[0162] To calculate the conditional expectation of the log-likelihood function, the expectation is obtained using the Kalman smoothing filter algorithm. First, we get:

[0163]

[0164] where j = k..., 1, 0, and can be obtained through the RTS backward filtering algorithm.

[0165]

[0166]

[0167] where is the mean of the backward filtering, is the variance of the backward filtering, and is the cross-covariance matrix of the backward filtering, M j is the filtering gain of the backward filtering.

[0168] Therefore, by substituting equations (25)-(28) into equation (21), the present invention can obtain:

[0169]

[0170] By further arranging equation (29), we get:

[0171]

[0172] Among them, represents the mean value of the running state of the backward filtering estimation at time j.

[0173] Equation (31) can be used to achieve the iterative estimation of all parameters:

[0174]

[0175] Among them, argmax() represents finding the maximum value, represents the value of the unknown parameter at the (l + 1)-th step of iteration at time k.

[0176] Let After estimating by setting the partial derivative of the function to 0, we have:

[0177]

[0178]

[0179] For the parameter b, taking the partial derivative according to equation (30) gives:

[0180]

[0181] By using the fminsearch function in MATLAB, the parameter estimation value is obtained by solving the non-linear equation

[0182] Remaining useful life prediction:

[0183] On the basis of considering the influence of the double time scale, measurement error, and two-way correlation of the two types of states, the remaining useful life L of the system is defined according to the first passage time k :

[0184]

[0185] Among them, inf{} represents the infimum, l k is the RUL of the system at time T k , Z(T k + l k ) represents the degradation state of the system at time T k + l k , w is the system failure threshold, and the probability density function of the remaining useful life of the system is:

[0186]

[0187] wherein, exp() represents the exponential function.

[0188] In order to predict the remaining useful life of the system in the implicit degradation state, define and Through the entire monitoring sequence V k at time T 1:k =(V0, V1,..., V k ) to obtain the conditional expectation k and variance of state Z That is Further, the PDF h(Z k |V k | 0:k ) is:

[0189]

[0190] Through the total probability formula, at the current time T k , the of the RUL of the system in the implicit degradation state is:

[0191]

[0192] wherein, represents the integral function.

[0193] Substitute equations (38) and (39) into equation (40) for calculation and simplification, and further organize to obtain:

[0194]

[0195] On this basis, considering the system operating state distribution That is obtain the PDF f k (l k |V 0:k ) of the RUL of the system under the influence of the double-time scale is:

[0196]

[0197] Substitute equation (41) into (42), and further obtain:

[0198]

[0199] wherein,

[0200] To verify the correctness and effectiveness of the system remaining life prediction method under the influence of double time scales proposed in the present invention, numerical experiments and case experiments are used for verification and analysis.

[0201] To verify the effectiveness of the proposed method, a numerical simulation experiment is carried out. The specific evolution model parameters are the operating state C = [0 1], Q = [0 (0.25) 2 T , R = 0.0625; the degradation function of the system degradation model is and the corresponding parameters are a = 1.6, b = 0.8, The system operating state is set at a periodic interval of Δt = 0.001, and the degradation state sampling interval is ΔT = 10 3 *Δt. Here, D = [0 1] is set T , that is, the change of the system state X2(t) affects the degradation state, and the steady-state value of the system state X2(t) is set to X 2w = 5, and the failure threshold of the system degradation state is w = 63.1.

[0202] The present invention uses the double time scale Kalman filter algorithm and the EM algorithm for joint estimation of the operating state, degradation state and degradation parameters. The operating state of the system at t ∈ [T k-1 , T k The curve is as Figure 3 shown. The system degradation state estimation curve during the entire degradation period T is as Figure 4 shown.

[0203] From Figure 3 and Figure 4 it can be seen that the operating state and degradation state of the system change continuously with time, and the predicted two types of states match well with the actual states of the system, verifying the effectiveness of the double time scale Kalman filter algorithm. The results of parameter estimation are as Figure 5 shown. Figure 6 The degradation trajectories considering the system state and not considering the system state are given. Among them, the first arrival time of the degradation trajectory without considering the system state is approximately 510 sampling points, and the first arrival time of the degradation curve considering the system state is approximately the 480th sampling point.

[0204] The parameters of the model proposed in the present invention are close to the given values, as shown in Figure 5 . From Figure 6 it can be seen that the time-varying operating state of the system affects its degradation process. Specifically, compared with not considering the time-varying operating state, the system degradation considering the time-varying operating state is accelerated.

[0205] ​The system degradation model without considering the influence of the operating state is denoted as M1, the system degradation model without considering the influence of the time scale is M2, and the method proposed in the present invention is denoted as M3. And it is compared with the remaining life results without considering the correlation and the time scale. Taking the monitoring moments from 400 to 500 as an example, the remaining life prediction results of different methods are as Figure 7 shown, and the comparison results of the remaining life means are as Figure 8 shown.

[0206] In order to verify the correctness of the prediction method of the present invention, the mean square error (MSE) evaluation index is used to compare the three prediction methods, and the results are shown in Figure 9 .

[0207] It can be seen from Figure 9 that compared with other prediction methods, the MSE error of the prediction method proposed in the present invention is the smallest and the prediction result is more accurate. Therefore, it is of great significance to consider the influence of the double time scale for the remaining life prediction.

[0208] In this implementation case, a simulation inertial platform stabilization loop system is used to verify the proposed remaining life prediction method, and the applicability of the method proposed in the present invention is proved. For the inertial platform, due to the phenomenon of efficiency loss in the actuator during operation, the system model will change, and the torque coefficient includes the internal structure characteristics and performance level of the motor. Therefore, it is assumed that the implicit degradation quantity is the torque coefficient K m , then its degradation process Z(T) is modeled as:

[0209]

[0210] In the formula, is the degradation drift coefficient affected by the load of the motor, that is, the armature current i m . It can be understood that the working stress of the actuator is related to its load, and the greater the load, the faster the system degradation rate. During the operation of the system, the performance of the actuator will decrease over time. The system failure threshold is w = 0.14 Nm / A, ω p and are the frame angular velocity and angle respectively. Let X1 = i m , X2 = ω p , The corresponding dynamic evolution model of the operating state of the entire inertial platform system is:

[0211]

[0212] At the same time, the corresponding degradation state model of the entire inertial platform system is:

[0213]

[0214] V k = GZ k + Λ k (50)

[0215] Wherein, represents the m-th operating state at the i-th moment, R m is the armature winding, L m is the armature inductance of the motor, K e is the back electromotive force coefficient of the motor, J is the moment of inertia, Q1, Q2 and Q3 are the process noises of the operating state system, and the operating state period is Δt = 0.001, and the degradation state period is ΔT = 10 3 *Δt, ξ k is the system noise of the degradation state, V k represents the measurement data, Λ k represents the measurement noise, and the steady-state value of the system state X1(t) is set to X 1w = 18 / A.

[0216] Using the dual-time-scale Kalman filtering algorithm proposed by the present invention, the operating state of the armature current and the degradation state of the torque coefficient of the inertial platform are estimated simultaneously. For the operating state of the system at t ∈ [T k-1 , T k , the curve is as Figure 10 shown. The results of the system degradation state estimation curve within the entire degradation period T are shown in Figure 11 shown.

[0217] It can be seen from Figure 10 and Figure 11 that by adopting the dual-time-scale Kalman filtering algorithm proposed by the invention, the operating state of the armature current can better track its true state, and at the same time, the estimated system actuator degradation state can track its true degradation state. As Figure 11 shown, the proposed time-varying degradation process with a parameter update process is very well matched with the implicit degradation estimated by the Kalman filter and the simulated actual degradation. Therefore, using the learning process in Equation (44) for remaining life prediction will decouple the mutual influence between the implicit degradation and the system state. The expectation of doing so is that the prognosis accuracy will not be greatly affected, and the complexity of predicting future degradation progress will be reduced. Figure 12 shows the remaining life probability density function curve based on the model proposed by the invention.

[0218] In order to further quantitatively compare the accuracy of the remaining life prediction results at different times, the mean absolute error (Mean Absolute Error, MRE) and root mean square error (Root Mean Square Error, RESE) of the remaining life are introduced, and the results are shown in Table 1.

[0219] Table 1 Error Comparison of Different Prediction Methods

[0220]

[0221] In Table 1, by comparing with different prediction methods, the relative errors of the MRE and RMSE of the method proposed in the invention are smaller than those of the methods that do not consider the time scale and the correlation. Further, it can be concluded that without considering the influence of the double time scale, it leads to an ill-conditioned problem in state estimation and cannot fully reflect the health level of the complex system under the influence of the correlation between the system operating state and the degradation state, which proves the effectiveness and superiority of the remaining useful life prediction method proposed in the invention.

[0222] The present invention aims at a complex degradation system that can be monitored in real time. Under the consideration of the influence of the dynamic environment, a system degradation modeling and remaining useful life prediction method based on the influence of double time scales is studied. A system state space model is established on the basis of considering the two-way correlation between the system operating state and the degradation state. This model not only considers the inherent evolution model and degradation model of the system itself, but also considers the influence of the system state on its degradation process, so it has more general adaptability. On the basis of the proposed model, an approximate analytical expression of the RUL probability density function is given, and the effectiveness of the proposed method is verified by using a case experiment of an inertial platform.

[0223] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, and improvements made within the spirit and principle of the present invention shall be included within the scope of the present invention.

Claims

1. A method for predicting the remaining useful life of a mechanical system considering the influence of a dynamic environment, characterized in that, The specific steps are as follows: Step 1. System degradation modeling: First, based on the monitor data that can be received to represent the operating state and degradation state, perform correlation analysis and processing on the two types of obtained states, and construct a system dynamics evolution model and a component degradation model; Step 2. According to the established system degradation model, adopt the dual-time-scale Kalman filtering algorithm and the Rauch-Tung-Striebel filtering algorithm to estimate the observable degradation state and operating state, and use the two-step expectation conditional maximization algorithm to estimate the model parameters; Step 3. On the basis of considering the dual-time-scale influence, measurement error, and two-way correlation of the two types of states, define the remaining life of the system according to the first passage time.

2. The method for predicting the remaining useful life of a mechanical system considering the influence of a dynamic environment according to claim 1, characterized in that In Step 1, the specific steps of system degradation modeling are: The system dynamics evolution model is: X(t) = A(Z(T))X(t - 1) + B(Z(T))U(t - 1) + ω(t) (1) Y(t) = C(Z(T))X(t) + ξ(t) (2) Where, X(t) represents the system operating state, A(Z(T)) is the state matrix, B(Z(T)) is the input matrix, C(Z(T)) is the output matrix, Y(t) is the measured value of the operating state data, t is the time scale of the system operating state, T is the time scale of the system degradation state, U(t) is the external input quantity, ω(t) is the system noise and ω(t)~(0,Q), Q is the covariance of the system noise, ξ(t) is the measurement noise of the operating data and ξ(t)~(0,R), R is the covariance of ξ(t); The system degradation model is: Model the degradation of the system to obtain: Among them, Z(T) is the degradation process of the system, and σ B B(T) is the diffusion process, and σ B is the diffusion coefficient, B(T) is the Brownian motion, is the system degradation rate, Ω is the parameter of the degradation rate function, and the average operating state is expressed as: Among them, D nx1 reflects the deviation of different components reflecting the system state from the steady-state value x ε on the impact of degradation, () T represents the transpose of a matrix; Set the corresponding degradation state observation equation: V(T) = GZ(T) + Λ(T) (5) Among them, V(T) represents the measured value of the degraded data, Z(T) represents the degradation state of the system, and Λ(T) represents the measurement noise of the degraded data, that is, Λ(T) ∼ MVN(0, σ R ) and σ R represents the variance of the measurement noise, and G represents the measurement coefficient of the degradation state; The deviation λ(t) = X(t) - x between any operating state X(t) at time t and the steady-state value x ω is called the offset. The change in λ(t) will affect the degradation rate η(T) of the system. When λ(t) = 0, the operating state of the system is the steady-state value X(t) ≡ x ω ; ω ; According to the correlation analysis of the system dynamics evolution model and the degradation model, establish a single-component system state space model.

3. The method for predicting the remaining useful life of a mechanical system considering the influence of a dynamic environment according to claim 2, wherein In Step 2, the degradation state and parameter estimation are: Given the initial values Z(T0), X(0), U(0), A(Z(T0)), B(Z(T0)), and C(Z(T0)), simultaneously estimate the system unit degradation time T using a filtering algorithm and the Expectation-Maximization (EM) algorithm k-1 of the degradation state Z(T k-1 ) and the parameter θ k-1 , calculate A(Z(T k-1 ))、B(Z(T k-1 ))、C(Z(T k-1 )) according to the evolution model. The change of the system state X(t) over time is for t ∈ [T k-1 , T k . Estimate the system operating state X(t i-1 ) at time t i-1 and determine whether the system operating time has reached t p . If it has not reached t p , then enter the next sampling time until it ends. At this time, Z(T k-1 ) remains unchanged, and the corresponding A(Z(T k-1 ))、B(Z(T k-1 ))、C(Z(T k-1 )) remain unchanged; When the system degradation moment enters the next sampling moment, according to the previous degradation period T k-1 obtained average system state X(T k-1 ), then estimate the degradation state Z(T k ) and parameter θ k . The corresponding impact of degradation on the system makes A(Z(T k-1 ))、B(Z(T k-1 )) and C(Z(T k-1 )) become A(Z(T k ))、B(Z(T k )) and C(Z(T k )); Determine whether the system degradation sampling time has reached T e If it has not reached T e then enter the next sampling time and continue to execute the above operations until the sampling ends; Degradation state estimation: Let the adoption period of the system state be Δt, the period of the degradation state be ΔT and ΔT >> Δt, and denote Suppose at the end of time T k the observed data of the system operating state at time t ∈ [T k-1 , T k the corresponding observed value is the corresponding system degradation state set is Z 0:k = {Z0, Z1..., Z k}; the set of measured values of the degradation data is V 0:k = {V0, V1..., V k}, assuming that at t ∈ (T, T + ΔT), the system degradation state Z(T) remains unchanged, design a filter to update the system state and the degradation state; In equations (1)-(5), assuming that the structure and parameters of the system are determined, the actual operating state X of the system i At t ∈ [T k-1 , T k , the actual degradation state Z of the system is estimated using the system operation data Y 0:ρ and is obtained by estimating the system degradation data V k . Let 0:k represent a nonlinear function and b be a nonlinear parameter. The specific double-time-scale Kalman filter algorithm is as follows: ​ 1). One-step prediction on the small time scale, including operating state prediction and covariance prediction: One-step update, including operating state update and covariance update: Among them, represents the estimated value of the operating state at time i, represents the operating state gain at time i, represents the covariance matrix of the predicted operating state X at time i, represents the covariance matrix of the estimated operating state X at time i, Y i represents the measurement data of the operating state at time i; Finally obtain and obtain the unit degradation period T k the expected average operating state and covariance as follows: Among them, tr() represents the trace of a matrix; represents the summation function; 2). One-step prediction on the large time scale, including degradation state prediction and covariance prediction: One-step update, including degradation state update and covariance update: Among them, represents the estimated value of the degradation state at time k, represents the degradation state gain at time k, represents the covariance matrix of the predicted degradation state Z at time k, represents the covariance matrix of the estimated degradation state Z at time k, V k represents the measurement data of the degradation state at time k, and finally obtains Parameter estimation: If the observed data V of the system degradation state is known 0:k , the system performance degradation state Z can be obtained based on Equation (3) k and the joint log-likelihood function lnf(Z 0:k , V k ; θ) of the observed data V 0:k : lnf(Z k ,V 0:k ; θ) = lnf(Z k ; θ) + lnf(V 0:k |Z k ; θ) (18) where, lnf(Z k ; θ) represents the probability density function of the system degradation state, lnf(V 0:k |Z k ; θ) represents the probability density function of the measurement data; ln() represents the logarithmic function; Obtained according to the assumptions of the model in Equation (18): Among them, represents the product of functions, represents the summation of functions from time j = 1 to j = k; τ represents the time interval and τ = T k -T k-1 , Z j-1 and Z j represent the values of state Z at times j - 1 and j respectively. Similarly, Substitute Equations (19) and (20) into the equation to obtain: Based on the l-th estimated result θ l , by solving the expectation of the log-likelihood function l n(Z k , V 0:k ; θ) with respect to the hidden state Z, we obtain E[lnf(Z k , V 0:k ; θ)]; In order to calculate the conditional expectation of the log-likelihood function, use the Kalman smoothing filtering algorithm to obtain the expectation: where j = k..., 1, 0, and is obtained by the RTS inverse filtering algorithm; Among them, is the mean of the inverse filtering, is the variance of the inverse filtering, is the cross-covariance matrix of the inverse filtering, M j is the filtering gain of the inverse filtering; Substitute Equations (25)-(28) into Equation (21) to obtain: Further organize Equation (29) to obtain: Among them, represents the mean of the running state of the inverse filtering estimation at time j; Use Equation (31) to realize the iterative estimation of all parameters: where argmax() represents obtaining the maximum value, represents the value of the unknown parameter at the (l + 1)-th iteration at time k; Let Estimate by setting the partial derivative of the function to 0 to obtain: For parameter b, take its partial derivative according to Equation (30) to obtain: Solve the non-linear equation to obtain the parameter estimate value 4. The method for predicting the remaining useful life of a mechanical system considering the influence of a dynamic environment according to claim 3, wherein In Step 3, the specific steps of remaining life prediction are: Based on considering the influence of double time scales, measurement errors, and the two-way correlation of two types of states, the remaining life L of the system is defined according to the first passage time k : where inf{} represents the infimum, and l k is the RUL of the system at time T k , and Z(T k +l k ) represents the degradation state of the system at time T k +l k . w is the system failure threshold, and the probability density function of the remaining life of the system is as follows: Among them, exp() represents the exponential function; In order to predict the remaining life of the system under the implicit degradation state, we define and By T k The entire monitoring sequence at time V 1:k =(V0,V1,...,V k )Get state Z k The conditional expectation of and variance Right now Further given the degenerate state Z k PDFh(Z k |V 0:k )for: The current T is obtained through the total probability formula k At the moment, the RUL of the system in the implicit degradation state is as follows: Among them, represents an integral function; Substitute Equations (38) and (39) into Equation (40) for calculation and simplification, and further organize to obtain: On this basis, consider the system operation state distribution That is Obtain the PDF f of the system RUL under the influence of double time scales k (l k |V 0:k ) is as follows: Substitute Equation (41) into (42) to further obtain: Among them,

Images