# Hybrid Degradation Equipment Remaining Useful Life Prediction Oriented Parallel Simulation considering Model Soft Switch

Authors: Chenglong Ge; Yuanchang Zhu; Yanqiang Di
Journal: Computational Intelligence and Neuroscience
 (2019)

---

## Abstract

Equipment parallel simulation is an emerging simulation technology in recent years, and equipment remaining useful life (RUL) prediction oriented parallel simulation is an important branch of parallel simulation. An important concept in equipment parallel simulation is the model evolution driven by real-time data, including model selection and model parameter evolution. The current research on equipment RUL prediction oriented parallel simulation mainly focuses on a single continuous degradation mode, such as linear degradation and nonlinear degradation. Under this degradation condition, the model parameter evolution methods in parallel simulation can effectively predict equipment RUL. However, in practice, most of the equipment degradation processes exhibit a mixture of continuous degradation and discrete shock. So this requires adaptive selection of simulation models based on real-time degradation data. In this paper, the hybrid degradation equipment RUL prediction oriented parallel simulation considering model soft switch is studied. Firstly, under the modeling framework of the state space model (SSM), two kinds of degradation simulation models are established using the Wiener process and Poisson effect. Driven by the real-time degradation data, the model probability is calculated by using the forward interactive multiple model filtering algorithm to realize the model soft switch and data assimilation. On the basis of model soft switch, the expectation maximization algorithm is utilized to achieve model parameter evolution. Through the iteration between model soft switch and model parameter evolution, the simulation fidelity can be effectively improved and the actual equipment degradation state is continuously approached. According to the full probability theorem and the concept of first hitting time, the simulated degradation state distribution is integrated into the inverse Gaussian distribution. Then the analytical expression of the RUL probability density function is obtained to achieve RUL real-time prediction. Finally, a case study was conducted by using a bearing degradation data. The results show that the parallel simulation can effectively model the hybrid degradation process of the bearing. Compared with the single-model method that only considers the model parameter evolution, the RUL obtained by the method proposed in this paper has higher prediction accuracy and smaller uncertainty.


---

## Body

## 1. Introduction

As the equipment complexity increases and the time-varying property of the equipment-operating environment enhances, equipment maintenance has become increasingly complicated. Due to external shock, wear, fatigue, corrosion, and other reasons, the equipment’s performance will be inevitably degraded, eventually causing equipment failure or even causing serious accidents. If we can determine the best opportunity for equipment maintenance and formulate the corresponding spare parts management and ordering plan based on the predicted RUL in the initial degradation stage, the reliability of equipment will be improved and the operational risks and operating costs of the equipment will be reduced effectively. In recent years, the prognostic and health management (PHM) technology has received more and more attention and has become an active investigation area in the reliability field [1]. PHM aims to accurately predict the equipment RUL. Accordingly, the reasonable maintenance and management of equipment are executed to guarantee the safety, reliability of equipment operation. The main tasks of PHM include RUL prediction and health management. Pecht and Chinnam [2, 3] both believe that RUL prediction is the core content of PHM and the RUL prognostic results provides a scientific basis for maintenance activities such as maintenance replacement and spare parts ordering. RUL prediction includes the probability density function (PDF) of RUL and its mathematical expectation. Since the RUL probability density function characterizes the uncertainty of RUL prediction, the probability density function is the primary predictor [4]. The current RUL prediction approaches can be broadly classified as failure physical model methods, statistics-based methods, and artificial intelligence methods [5]. For complex equipment, its failure mechanism is difficult to obtain, so the latter two methods get more attention. Statistics-based methods achieve RUL prognostic results via data fitting on the basis of statistics models. The commonly used statistics-based methods include the Markov chain [6], Bayesian method [7], inverse Gaussian process [8], Gamma process [9], Wiener process [10, 11], etc. The artificial intelligence methods achieve RUL prognostic results via data fitting mainly including machine learning. The commonly used artificial intelligence methods include neural network [12], support vector machine [13], etc. However, the artificial intelligence methods are limited by the two deficiencies in the practical application. One is that it usually cannot achieve the analytical expression of the probability density function of RUL, which restricts the real-time application of the methods [14]. The other is that it also requires a large amount of training samples, which usually cannot be satisfied in practice [15]. Contrastively, statistics-based methods do not depend on a large amount of training samples and are more flexible because of its independency regarding specific application objects, making it widely utilized in degradation modeling and RUL prediction. Among these statistics-based methods, the Wiener process is the most popular method because of the following advantages. On the one hand, compared with other methods (e.g., inverse Gaussian process and Gamma process), the Wiener process can model not only monotonous degradation paths but also nonmonotonous degradation paths. On the other hand, the Wiener process has excellent physical interpretation and mathematical property (i.e., its first hitting time distribution is inverse Gaussian distribution), which is helpful to obtain the analytical expression of the PDF of RUL.Recent advances about Wiener process-based RUL prediction methods were discussed in the most recent review paper [16]. In this review, Zhang et al. systematically reviewed the conventional Wiener process-based RUL prediction methods and many generalizations and variants from the conventional Wiener process-based RUL prediction methods by considering the factors of nonlinearity [17], multisource variability [18], covariates [19], and multivariates [20]. In the above numerous applications, Wiener process-based methods showed a flexible modeling capability and has undergone extensive development. But there still exists a typical issue deserving further research. The issue is that it is usually assumed that the equipment degradation mode is a fixed continuous degradation in the entire degradation period among the Wiener process-based prediction methods. As a result, a single and fixed model is used for RUL prediction. However, besides the continuous degradation, the degradation mode may present a hybrid of multiple degradation modes in practice, which requires multiple models. Specially, the damage caused by the random shock is also an important reason accounting for equipment failure. In this paper, the degradation mode incorporating the continuous degradation and the random shock is called the hybrid degradation mode. Until now, in order to solve the RUL prediction issue of the hybrid degradation equipment, several researchers have attempted to establish the prognostic models on the basis of the Wiener process.Si et al. [21] presented a new prognostic model to characterize the hybrid degradation equipment. In the proposed model, the linear Wiener model is used to describe the continuous degradation process and a compound Poisson process is utilized to characterize the randomly arriving shock. And each randomly arriving shock is described by a random variable which obeys the normal distribution with two parameters. Inspired by the research of Si et al. [21], Zhang et al. [22] regarded the hybrid degradation process as the state switching process and utilized a nonlinear Wiener process with state switching depicted by a continuous time Markov chain (CTMC) to describe the continuous degradation process. And also a random variable which obeys the normal distribution is used to characterize each shock introduced in the state switching. Regrettably, the model parameters estimation method was not investigated. On the basis of the previous studies [21, 22], Zhang et al. [23] proposed a hybrid degradation model with a continuous model described by a nonlinear Wiener process and a randomly arriving shock model characterized by a nonhomogeneous compound Poisson process. They not only obtained the approximated analytical lifetime under the concept of FHT, but also obtained parameters updating formulas by combining the expectation conditional maximization (ECM) algorithm and maximum likelihood estimation (MLE). Additionally, a numerical example and a case study of furnace wall were studied. Furthermore, similar to the previous study [23], Zhang et al. [24] utilized a specific nonlinear Wiener process with the power low model to describe the continuous degradation and the same method to describe the randomly arriving shock. On the basis of the above modeling mechanism, Du et al. [25] proposed a more generalized hybrid degradation model consisting of trend term and stochastic fluctuating term. It is worth noting that all the above latest research studies use one model to model the hybrid degradation process. Considering the unknown characteristic of arriving time and the amplitudes for the shock, it is a reasonable choice to construct multiple models to describe the hybrid degradation process and replace the model dynamically through model soft switch. Therefore, in order to accurately predict the RUL of hybrid degradation equipment, a novel real-time prediction method is needed urgently considering both the model soft switch and online evolution of model parameters. Equipment parallel simulation has become a possible solution [26, 27]. This method is called hybrid degradation equipment RUL prediction oriented parallel simulation.Equipment parallel simulation is an emerging simulation paradigm which aims to combine the simulation system with the actual equipment. It was originally defined at the AsiaSim/SCS AutumnSim conference in 2016. The simulation system benefits from the online acquisition of equipment information to update the simulation model. Conversely, the equipment benefits from the simulation results of the simulation system to improve the equipment performance. Particularly, the simulation system running in this mode is called the parallel simulation system. The equipment parallel simulation diagram is shown in Figure1. In the equipment parallel simulation, the actual equipment and simulation system exchange information through sensors and actuators. The sensor provides the actual equipment information to the parallel simulation system, and the actuator lets the parallel simulation system to perform control and other operations on the actual equipment. The actual equipment information provided by the sensors can be divided into two categories, namely, observable state information St and behavior information Bt at time t. Equipment can be controlled by the control information Ct sent by the parallel simulation system. According to the implementation of control information, it can be divided into automatic control information ACt and manual control information MCt.Figure 1
Overview of the equipment parallel simulation.An important concept in equipment parallel simulation is model evolution driven by the real-time data, including two aspects of model adaptive selection and model parameter evolution. This is considered to be a typical feature that distinguishes it from previous simulation techniques. In the past simulation technologies, the simulation model focused on one-time construction. After the simulation operation, the simulation model and model parameters no longer change, i.e., there is no model evolution process. Although the word “parallel” in parallel simulation translates as “parallel” in English, it is essentially different from “parallel” in parallel simulation proposed for many years. The former refers to expand the actual complex problem into the virtual space, and the actual complex problem is handled by the interaction between virtual space and actual space. It has the same connotation as “parallelism” in parallel system theory [28]. The latter refers to dividing the actual complex problem into several subproblems handled simultaneously [29]. Theoretical origins of parallel simulation are closely related to parallel system theory [28, 30], dynamic data drive application system (DDDAS) [31, 32], symbiotic simulation [33, 34], and online simulation [35]. The literature [29] has been reviewed in detail, but there are also significant differences. Parallel system theory emphasizes agent-based modeling. DDDAS introduces the idea of control theory into the simulation field, which stresses utilizing the simulation results to control the measurement process. Symbiotic simulation underlines using what-if analysis (WIA) to execute multiscenario simulation. Online simulation emphasizes the connecting relationship between the actual system and the simulation system, which is a contrary concept to offline simulation.According to its technical principle and typical characteristics, equipment parallel simulation provides an effective way to solve the RUL prediction issue of hybrid degradation equipment. In the field of mechanical equipment, components such as bearings and gearboxes are widely used and they are all critical components. The failure of these components will cause the entire equipment to be shut down. Therefore, the performance degradation process of these components is often used to measure the RUL of mechanical equipment. The degradation process of these critical components is a typical hybrid degradation process. In this paper, with the background of RUL prediction issue of hybrid degradation equipment, the simulation model construction, model soft switch, model parameter evolution, and RUL prediction of hybrid degradation equipment RUL prediction oriented parallel simulation are investigated.The structure of the paper is as follows. Based on the modeling analysis, the proposed simulation model is constructed in Section2. The evolution method of the parallel simulation model is put forward in Section 3, which includes model probability based model soft switch and expectation maximization algorithm based model parameter evolution. Then, the RUL real-time prediction method of hybrid degradation equipment based on parallel simulation is presented in Section 4. The parallel simulation method given in this paper is validated by using a bearing degradation data, and comparative study is also conducted in Section 5. Finally, some conclusions and future perspectives are discussed in Section 6.


## 2. Modeling of Parallel Simulation

### 2.1. Modeling Analysis

The simulation model and its evolution belong to the model theory category of equipment parallel simulation, and it is also the basic issue in the research of equipment parallel simulation. Due to the strong field correlation of the model coupled with evolution characteristic, it becomes a difficult research issue in equipment parallel simulation. Therefore, it is the primary basic problem to determine the modeling method and model form in hybrid degradation equipment RUL prediction oriented parallel simulation.Literatures [27, 36] pointed out that building a state space model for equipment performance degradation is an reasonable modeling direction for equipment RUL prediction oriented parallel simulation. In the SSM modeling approach, the simulation output is a hidden degradation state. The equipment degradation SSM includes a degradation state equation and an observation equation. The former describes the relationship between degradation states at adjacent moments, and the latter describes the relationship between the observation and degradation state. In the equipment degradation SSM, the dynamic and time-varying characteristics for the degradation process are both taken into consideration, which is helpful for the simulated degradation state estimation and RUL prediction. Considering that statistics-based methods are easier to obtain the analytical expression of the RUL probability density function, this paper combines the stochastic process approach with SSM to develop the parallel simulation model. The stochastic process is suitable for describing the randomness and uncertainty of the degradation process. The Wiener process and Gamma process are the most commonly used stochastic processes, but the application conditions of the latter are too harsh. The Gamma process is only suitable for describing the degradation processes with strictly monotonic characteristic. Contrarily, the Wiener process is suitable for describing the nonmonotonic degradation process, which has a more relaxed application conditions. Therefore, it is an advisable choice to construct the Wiener state space model (WSSM) by combining the SSM modeling method and the Wiener process in the parallel simulation modeling [36]. In particular, in order to describe the hybrid degradation process with unknown discrete shocks, a hybrid Wiener state space model (HWSSM) should be established as the parallel simulation model, which includes the continuous degradation model and the degradation model with unknown discrete shocks.


### 2.2. Construction of HWSSM

In order to construct HWSSM, two equations need to be developed, i.e., hybrid degradation state equation and observation equation. The Wiener process and shock effect are used to construct the hybrid degradation state equation with two forms, including the continuous degradation state equation and the degradation state equation with unknown discrete shock. The continuous degradation state equation can be formulated as(1)xt=x0+ηt+σBt,where xt,t≥0 is the continuous degradation process driven by the standard Brownian motion Bt,t≥0 and Bt∼N0,t; x0 is an initial degradation state; η and σ are the drift coefficient and diffusion coefficient of the standard Brownian motion [37]. Then, equation (1) is transformed by Euler discretization and the degradation state equation at discrete time points tkk=1,2,… without considering the unknown discrete shocks is yielded as(2)xk=xk−1+ητk+στkϖk,where τk=tk−tk−1 is the sensor sampling interval and xk=xtk is the simulated degradation state at time tk. Time tk is short for time, k and ϖk is the noise sequence which obeys the standard normal distribution.The effect of unknown discrete shocks on equipment performance can be represented by a shock variable. According to the discrete characteristic of the shocks, the damage caused by the shocks is integrated into equation (2), and a degradation state equation with discrete shocks is obtained by(3)xk=xk−1+ητk+στkϖk+D,where D is the damage caused by the shocks. We assume that the arrival of the shocks obeys the Poisson process [38, 39], i.e.,(4)PMtk+Δt−Mtk=n=ρΔtnn!e−ρΔt,where Mtk denotes the total number of shocks that appear from the initial time until the moment k and ρ denotes the shocks arrival rate. Specifically, considering that the sampling interval of the sensor is small, the probability of more than one shock occurring within one sampling interval is very small. In other words, the probability of more than one shock occurring within one sampling interval can be neglected, and n is zero or one. In addition, it is assumed that the Poisson process is independent of the Brownian motion.Equipment degradation observation datayt are obtained by sensor measurement. The random relationship between the simulated degradation state and the observation data can be described by the observation equation, i.e.,(5)yt=xt+πt,where πt∼N0,ϕ2 and ϕ2 denotes the variance of the measurement noise. It is also assumed that πt is independent of the Brownian motion. Then, equation (5) is transformed by Euler discretization, and the observation equation at discrete time points tkk=1,2,… is yielded as(6)yk=xk+ϕζk,where yk=ytk denotes the degraded observation data at time k and ζk is the noise sequence that obeys the standard normal distribution. Furthermore, ζk is independent of ϖk.According to equations (1)–(6), the HWSSM at the discrete time point can be formulated as(7)xk=xk−1+ητk+στkϖk,ifMtk−Mtk−1=0,xk−1+ητk+στkϖk+D,ifMtk−Mtk−1=1,yk=xk+ϕζk.


## 2.1. Modeling Analysis

The simulation model and its evolution belong to the model theory category of equipment parallel simulation, and it is also the basic issue in the research of equipment parallel simulation. Due to the strong field correlation of the model coupled with evolution characteristic, it becomes a difficult research issue in equipment parallel simulation. Therefore, it is the primary basic problem to determine the modeling method and model form in hybrid degradation equipment RUL prediction oriented parallel simulation.Literatures [27, 36] pointed out that building a state space model for equipment performance degradation is an reasonable modeling direction for equipment RUL prediction oriented parallel simulation. In the SSM modeling approach, the simulation output is a hidden degradation state. The equipment degradation SSM includes a degradation state equation and an observation equation. The former describes the relationship between degradation states at adjacent moments, and the latter describes the relationship between the observation and degradation state. In the equipment degradation SSM, the dynamic and time-varying characteristics for the degradation process are both taken into consideration, which is helpful for the simulated degradation state estimation and RUL prediction. Considering that statistics-based methods are easier to obtain the analytical expression of the RUL probability density function, this paper combines the stochastic process approach with SSM to develop the parallel simulation model. The stochastic process is suitable for describing the randomness and uncertainty of the degradation process. The Wiener process and Gamma process are the most commonly used stochastic processes, but the application conditions of the latter are too harsh. The Gamma process is only suitable for describing the degradation processes with strictly monotonic characteristic. Contrarily, the Wiener process is suitable for describing the nonmonotonic degradation process, which has a more relaxed application conditions. Therefore, it is an advisable choice to construct the Wiener state space model (WSSM) by combining the SSM modeling method and the Wiener process in the parallel simulation modeling [36]. In particular, in order to describe the hybrid degradation process with unknown discrete shocks, a hybrid Wiener state space model (HWSSM) should be established as the parallel simulation model, which includes the continuous degradation model and the degradation model with unknown discrete shocks.


## 2.2. Construction of HWSSM

In order to construct HWSSM, two equations need to be developed, i.e., hybrid degradation state equation and observation equation. The Wiener process and shock effect are used to construct the hybrid degradation state equation with two forms, including the continuous degradation state equation and the degradation state equation with unknown discrete shock. The continuous degradation state equation can be formulated as(1)xt=x0+ηt+σBt,where xt,t≥0 is the continuous degradation process driven by the standard Brownian motion Bt,t≥0 and Bt∼N0,t; x0 is an initial degradation state; η and σ are the drift coefficient and diffusion coefficient of the standard Brownian motion [37]. Then, equation (1) is transformed by Euler discretization and the degradation state equation at discrete time points tkk=1,2,… without considering the unknown discrete shocks is yielded as(2)xk=xk−1+ητk+στkϖk,where τk=tk−tk−1 is the sensor sampling interval and xk=xtk is the simulated degradation state at time tk. Time tk is short for time, k and ϖk is the noise sequence which obeys the standard normal distribution.The effect of unknown discrete shocks on equipment performance can be represented by a shock variable. According to the discrete characteristic of the shocks, the damage caused by the shocks is integrated into equation (2), and a degradation state equation with discrete shocks is obtained by(3)xk=xk−1+ητk+στkϖk+D,where D is the damage caused by the shocks. We assume that the arrival of the shocks obeys the Poisson process [38, 39], i.e.,(4)PMtk+Δt−Mtk=n=ρΔtnn!e−ρΔt,where Mtk denotes the total number of shocks that appear from the initial time until the moment k and ρ denotes the shocks arrival rate. Specifically, considering that the sampling interval of the sensor is small, the probability of more than one shock occurring within one sampling interval is very small. In other words, the probability of more than one shock occurring within one sampling interval can be neglected, and n is zero or one. In addition, it is assumed that the Poisson process is independent of the Brownian motion.Equipment degradation observation datayt are obtained by sensor measurement. The random relationship between the simulated degradation state and the observation data can be described by the observation equation, i.e.,(5)yt=xt+πt,where πt∼N0,ϕ2 and ϕ2 denotes the variance of the measurement noise. It is also assumed that πt is independent of the Brownian motion. Then, equation (5) is transformed by Euler discretization, and the observation equation at discrete time points tkk=1,2,… is yielded as(6)yk=xk+ϕζk,where yk=ytk denotes the degraded observation data at time k and ζk is the noise sequence that obeys the standard normal distribution. Furthermore, ζk is independent of ϖk.According to equations (1)–(6), the HWSSM at the discrete time point can be formulated as(7)xk=xk−1+ητk+στkϖk,ifMtk−Mtk−1=0,xk−1+ητk+στkϖk+D,ifMtk−Mtk−1=1,yk=xk+ϕζk.


## 3. HWSSM Evolution

HWSSM belongs to the state space model with hidden degradation state. The parallel simulation system realizes the evolution of HWSSM by the following two ways. The first is the model soft switch based on the interactive multiple model (IMM) filtering [40, 41]. Through the IMM filtering, the probabilities of different simulation models are calculated dynamically, and the data assimilation between the observation data and the simulation output is achieved. As a result, the simulation output is updated and the estimation of the simulated degradation state is obtained. Another way is the model parameter evolution based on parameter online estimation. The model parameters are updated by utilizing the latest observation data via the parameter evolution. The model soft switch and parameter evolution are not executed in isolation, but they are iterative to each other. Through the iteration of the two ways, the simulation output is continuously approaching the equipment’s actual degradation state. Finally, a high-fidelity simulation model is provided for accurately predicting the equipment RUL.


### 3.1. IMM Filtering-Based Model Soft Switch

In the HWSSM, the unknown characteristic of the shocks make it impossible to know whether there is a shock arrival from the observed data, which leads to no way to inform the conditions of switching the simulation model. Consequently, it is necessary to calculate the probabilities of different simulation models. In the SSM modeling method, the IMM filtering provides an effective way for calculating the probability of different models, and the Kalman filter is utilized in IMM filtering. The simulation model soft switch includes four stages, i.e., model input interaction, Kalman filtering, model probability calculation, and model output interaction.For convenience, the following definitions are given.Yk=y1,y2,…,yk and Xk=x0,x1,x2,…,xk represent the observation data vector and simulated degradation state vector until time k, respectively. m represents the number of simulation models. μku represents the probability of simulation model u at time k. pvu represents the transition probability from the model v to the model u. P=pvum×m represents the Markov probability transition matrix. Cku denotes that the valid model is model u in time interval k−1,k. x^kku denotes the estimated mean of the simulated degradation state at time k conditioned on Cku and the previous k measured values. Pkku denotes the estimated covariance of the simulated degradation state at time k conditioned on Cku and the previous k measured values. x^kk0u denotes the hybrid estimated mean of the simulated degradation state at time k conditioned on Ck+1u and the previous k measured values. Pkk0u denotes the hybrid estimated covariance of the simulated degradation state at time k conditioned on Ck+1u and the previous k measured values. Nx;0,κ denotes the Gaussian distribution with value x, mean 0, and variance κ.


#### 3.1.1. Input Interaction (Modelu)

At this stage, the hybrid state estimationx^k−1k−10u and the covariance estimation Pk−1k−10u are obtained by the state estimation x^k−1k−1v and model probability μk−1v of model v at time k−1. x^k−1k−10u and Pk−1k−10u are both used as the initial state of the Kalman filtering at the time k. The specific steps are as follows.The prior probabilityG¯u of the simulation model u is defined and calculated by(8)G¯u=ΔpCkuYk−1=∑v=1mpCkuCk−1v,Yk−1pCk−1vYk−1=∑v=1mpvuμk−1v,where pvu=ΔPCkuCk−1v,Yk−1.The hybrid probabilityμk−1k−1vu transferred from the model v to the model u is defined and calculated by(9)μk−1k−1vu=ΔpCk−1vCku,Yk−1=pCkuCk−1v,Yk−1pCk−1vYk−1pCkuYk−1=pvuμk−1vG¯u.The hybrid state estimationx^k−1k−10u of simulation model u is determined by(10)x^k−1k−10u=∑v=1mx^k−1k−1vμk−1k−1vu.The hybrid covariance estimationPk−1k−10u of simulation model u is determined by(11)Pk−1k−10u=∑v=1mμk−1k−1vuPk−1k−1v+x^k−1k−1v−x^k−1k−10u⋅x^k−1k−1v−x^k−1k−10uT.


#### 3.1.2. Kalman Filtering (Modelu)

Kalman filtering can be divided into two steps, i.e., the prediction step and the updating step.x^kk−1u and Pkk−1u express the simulated degradation state estimation and covariance estimation, respectively, conditioned on Cku and the previous k−1 measured values. Then, the prediction step refers to achieve the predicted results x^kk−1u and Pkk−1u according to the degradation equations and the results of input interaction. The prediction step can be expressed as(12)x^kk−11=x^k−1k−101+ητk,(13)x^kk−12=x^k−1k−102+ητk+D,(14)Pkk−1u=Pk−1k−10u+σ2τk.The updating step refers to obtain the posterior state estimationx^kku and covariance estimation Pkku of model u based on the prognostic results in the prediction step and the observed data yk. The updating step can be expressed as(15)y˜ku=yk−x^kk−1u,Sku=Pkk−1u+ϕ2,Kku=Pkk−1uSku−1,x^kku=x^kk−1u+Kkuy˜ku,Pkku=I−KkuPkk−1u,where y˜ku denotes the new information, Sku is the variance of y˜ku, Kku is the Kalman gain, and I represents one-dimensional unit matrix.


#### 3.1.3. Model Probability Calculation

At this stage, the likelihood functionΛku of the model u is used to calculate the model probability and is defined as(16)Λku=ΔpykCku,Yk−1≈Nyk;x^kk−1u,Pkk−1u+ϕ2.Furthermore, according to equation (16), Λku can be given by(17)Λku=12π1/2Sku1/2exp−12y˜kuTSku−1y˜ku.Then, the probability of the simulation modelu is(18)μku=ΔpCkuYk=ΛkuG¯uGu,where Gu is a normalization constant that satisfies Gu=∑u=1mΛkuG¯u.


#### 3.1.4. Output Interaction

According to the model probabilityμku and the estimation results of each model in the Kalman filtering stage, the degradation state estimation x^kk and covariance estimation Pkk are achieved through the weighted sum and can be given by(19)x^kk=∑u=1mx^kkuμku,(20)Pkk=∑u=1mμkuPkku+x^kku−x^kkx^kku−x^kkT.In addition, in the initial time of degraded state estimation, there existsx00u=x00v and pC0u=pC0v conditioned on ∀u,v∈m,u≠v. Once a new measurement value yk+1 is obtained, the model probability can be updated to realize model soft switch according to equations (8)–(20).


### 3.2. Model Parameter Evolution

The model parametersθ of HWSSM include μ0, Σ0, η, σ, D, and ϕ, where μ0 and Σ0 denote the mean and covariance of the initial simulated degradation state, respectively. Since HWSSM is a state space model with hidden state, it is considered to use the expectation maximization (EM) algorithm to achieve model parameter evolution. EM algorithm is suitable for estimating model parameters of SSM with hidden state [42, 43]. Based on the model soft switch, this paper uses EM algorithm to realize the evolution of model parameters.According to the EM algorithm, the estimations of the model parametersθ at time k and the jth iteration of EM algorithm can be obtained by(21)θ^kj=argmaxθEXk,R˜1:kYk,θ^kj−1LXk,Yk,R˜1:kθYk,θ^kj−1,where E⋅ denotes the mathematical expectation operator, LXk,Yk,R˜1:kθ is the joint log-likelihood function, and R˜1:k represents the model indicating matrix. R˜k represents the model indicating vector at time k and satisfies R˜k=R˜k1,R˜k2. Its element R˜kv is defined as(22)R˜kv=1,conditioned onCkv,0,others.EM algorithm consists of two steps, i.e., E-step and M-step. E-step refers to calculate the mathematical expectation of the joint log-likelihood function, and M-step refers to maximize the mathematical expectation of the joint log-likelihood function. The model parameters can be estimated online via the iteration of E-step and M-step.


#### 3.2.1. E-Step

Based on the Markov property and the multiplicative equation of conditional probability, the joint log-likelihood function of the simulated state vectorXk, the observed data vector Yk, and the model indicating matrix R˜1:k at time k is given as(23)LXk,Yk,R˜1:kθ=lnpXk,Yk,R˜1:kθ=lnpYkXk,θpXk,R˜1:kθ=lnpYkXk,θ︸①+lnpXk,R˜1:kθ︸②.According to HWSSM, there exists(24)pyixi,θ=Nyi;xi,ϕ2;pxixi−1,R˜iv,R˜1:i−1,θ=Nxi;xi−1+χiv,σ2τi;px0θ=Nx0;μ0,Σ0,where χi1=ητi and χi2=ητi+D. Then, the derivation of equation (23) is as follows:(25)①=ln∏i=1kpyixi,θ=∑i=1klnpyixi,θ=∑i=1kln12πϕ2exp−yi−xi22ϕ2=−k2ln2π−k2lnϕ2−12ϕ2∑i=1kyi−xi2,②=ln∏i=2kpxi,R˜ixi−1,R˜1:i−1,θpx1,R˜1x0,θpx0θ=lnpx0θ∏i=1k∏v=12μivNxi;xi−1+χiv,σ2τiR˜iv=lnpx0θ︸③+ln∏i=1k∏v=12μivNxi;xi−1+χiv,σ2τiR˜iv︸④,③=lnpx0θ=−12ln2π−12lnΣ0−x0−μ022Σ0,④=∑i=1k∑v=12lnμivNxi;xi−1+χiv,σ2τiR˜iv=∑i=1k∑v=12R˜ivlnμivNxi;xi−1+χiv,σ2τi=∑v=12∑i=1kR˜ivlnμiv+lnNxi;xi−1+χiv,σ2τi=∑v=12∑i=1kR˜ivlnμiv+ln12πτi−lnσ−12σ2τixi−xi−1+χiv2.Based on the above derivation, the final expression of the joint log-likelihood function is formulated as(26)LXk,Yk,R˜1:kθ=−k2ln2π−k2lnϕ2−12ϕ2∑i=1kyi−xi2+∑v=12∑i=1kR˜ivlnμiv+ln12πτi−lnσ−12σ2τixi−xi−1+χiv2−12ln2π−12lnΣ0−x0−μ022Σ0.Next,Qθ,θkj is denoted as the mathematical expectation of the joint log-likelihood function at time k and the jth iteration of EM algorithm, i.e.,(27)Qθ,θkj=EXk,R˜1:kYk,θ^kj−1LXk,Yk,R˜1:kθYk,θ^kj−1.Through neglecting the irrelevant items that are independent of the parameterθ, Qθ,θkj can be expressed as(28)Qθ,θkj∝−∑i=1k12lnϕ2+12ϕ2EXk,R˜1:k⋅yi−xi2⋅︸⑤+∑v=12∑i=1kEXk,R˜1:k⋅R˜iv⋅︸⑥lnμiv+ln12πτi−lnσ−EXk,R˜1:k⋅xi−xi−1−χiv2⋅︸⑦2σ2τi−12lnΣ0−12Σ0EXk,R˜1:k⋅x0−μ02Yk,θ^kj−1︸⑧,where ⋅ represents Yk,θ^kj−1. According to the property of mathematical expectation, the following intermediate variables are defined as(29)x^ik=EXk,R˜1:k⋅xi⋅,Pik=varXk,R˜1:k⋅xi2⋅=EXk,R˜1:k⋅xi2⋅−EXk,R˜1:k⋅xi⋅2=EXk,R˜1:k⋅xi2⋅−x^ik2,Pi,i−1k=covXk,R˜1:k⋅xixi−1⋅=EXk,R˜1:k⋅xixi−1⋅−EXk,R˜1:k⋅xi⋅EXk,R˜1:k⋅xi−1⋅=EXk,R˜1:k⋅xixi−1⋅−x^ikx^i−1k.Based on the definition of the above variables, the derivation of the four conditional mathematical expectations in equation (28) is as follows:(30)⑤=EXk,R˜1:k⋅yi2+xi2−2yixi⋅=yi2+EXk,R˜1:k⋅xi2⋅−2yiEXk,R˜1:k⋅xi⋅=yi2+Pik+x^ik2−2yix^ik,⑥=EXk,R˜1:k⋅R˜iv⋅=pR˜iv=1⋅=ωiv,⑦=EXk,R˜1:k⋅xi−xi−1−χiv2⋅=EXk,R˜1:k⋅χiv2+xi−xi−12−2χivxi−xi−1⋅=χiv2+x^ik2+Pik+x^i−1k2+Pi−1k−2x^ikx^i−1k+Pi,i−1k−2χivx^ik−x^i−1k,⑧=EXk,R˜1:k⋅x0−μ02⋅=EXk,R˜1:k⋅x02+μ02−2μ0x0⋅=x^0k2+P0k+μ02−2μ0x^0k.After deriving the above four parts of the conditional mathematical expectation,Qθ,θkj can be given by(31)Qθ,θkj∝−∑i=1k12lnϕ2+12ϕ2yi2+Pik+x^ik2−2yix^ik+∑v=12∑i=1kωivlnμiv+ln12πτi−lnσ−χiv2+Γi−2χivx^ik−x^i−1k2σ2τi−12lnΣ0−12Σ0x^0k2+P0k+μ02−2μ0x^0k,where Γi=x^ik2+Pik+x^i−1k2+Pi−1k−2x^ikx^i−1k+Pi,i−1k. According to equation (31), in order to calculate Qθ,θkj, the values of the variables x^ik, x^i−1k, Pik, Pi−1k, Pi,i−1k, and ωiv should be acquired, which belong to the smoothed variables. This paper uses the IMM backward smoothing (IMMBS) algorithm to obtain the values of the above six smoothed variables [44, 45].IMMBS algorithm consists of backward filtering and model state fusion. Backward filtering refers to the backward filtering from the latest measured value. The operation flow of backward filtering is similar to the forward IMM filtering mentioned in Section3.1, but there are also obvious differences between them. The forward IMM filtering performs the input interaction before performing one-step prediction, while the backward filtering performs one-step prediction before performing the input interaction. Backward filtering can be divided into five steps, including backward one-step prediction, backward input interaction, backward filtering update, backward model probability calculation, and backward output fusion. In particular, the backward one-step prediction equation for backward filtering is given by(32)x^ii+1B,1=x^i+1i+1B,1−ητi+1,x^ii+1B,2=x^i+1i+1B,2−ητi+1−D,Pii+1B,u=Pi+1i+1B,u+σ2τi+1,where x^ii+1B,u denotes the backward one-step prediction of the model u and Pii+1B,u represents the covariance of backward one-step prediction error for the model u. The other steps of backward filtering can be found in [44]. In the stage of model state fusion for the IMMBS algorithm, according to the full probability theorem, the smoothing estimation of the simulated state can be expressed as(33)pxiYk=∑u=1mωiupxiCiu,Yk=∑u=1mωiu∑v=1mω˜i+1kvupxiCiu,Ci+1v,Yk=∑u=1mωiu∑v=1mω˜i+1kvuΛiuv,where ωiu represents the smoothed model probability, ω˜i+1kvu denotes the smoothed hybrid probability, and Λiuv denotes the likelihood function. The definition and calculation of the three variables are given by(34)ω˜i+1kvu=ΔpCi+1vCiu,Y1k=1G˜upCi+1vCiu,YipYi+1kCi+1v,Ciu,Yi=1G˜upuvpYi+1kCi+1v,Ciu,Yi,Λiuv=ΔpYi+1kCi+1v,Ciu,Yi≈pxii+1B,vCi+1v,Ciu,xiiu=Nx^ii+1B,v−x^iiu;0,Pii+1B,v+Piiu,ωiu=ΔpCiuY1k=G˜uGpCiuYi=G˜uGμiu,where G˜u and G are all the normalized constants. They can be acquired by(35)G˜u=∑v=1mpuvΛiuv,G=∑u=1mG˜uμiu.In equation (33), with respect to pxiCiu,Ci+1v,Yk, the mixed smoothing state estimation x^ikuv, covariance estimation Pikuv, and the mixed interaction covariance estimation Pi,i−1kuv can be obtained by(36)x^ikuv=PikuvPiiu−1x^iiu+Pii+1B,v−1x^ii+1B,v,Pikuv=Piiu−1+Pii+1B,v−1−1,Pi,i−1kuv=Piiu−1+Pii+1B,v−1−1+Pi−1i−1u−1+Pi−1iB,v−1−1−1.In equation (33), with respect to pxiCiu,Yk, the mixed smoothing state estimation x^iku, covariance estimation Piku, and the mixed interaction covariance estimation Pi,i−1ku of model u can be obtained by(37)x^iku=∑u=1mω˜i+1kvux^ikuv,Piku=∑u=1mω˜i+1kvuPikuv+x^ikuv−x^ikux^ikuv−x^ikuT,Pi,i−1ku=∑u=1mω˜i+1kvuPi,i−1kuv+x^ikuv−x^ikux^i−1kuv−x^i−1kuT.In particular, ifi=k, there exists(38)Pk,k−1ku=I−KkuPk−1k−1u.Finally, after the smoothed state fusion of each model, the smooth estimationx^ik, covariance estimation Pik, and interactive covariance estimation Pi,i−1k of the simulated degradation state can be acquired by(39)x^ik=∑u=1mωiux^iku,Pik=∑u=1mωiuPiku+x^iku−x^ikx^iku−x^ikT,Pi,i−1k=∑u=1mωiuPi,i−1ku+x^iku−x^ikx^i−1ku−x^i−1kT.In particular, the initial smoothed estimations are given by(40)x^0k=∑u=1mω1ux^01B,u,P0k=∑u=1mω1uP01B,u+x^01B,u−x^0kx^01B,u−x^0kT.


#### 3.2.2. M-Step

The estimation of the model parametersθ at the jth step of the EM algorithm can be obtained by differentiating equation (31), i.e.,(41)∂Qθ,θkj∂θ=0.The updated parameters are given by(42)μ^0=x^0k,Σ^0=P0k,D^=∑i=1kωi2x^ik−x^i−1k−ητi/τi∑i=1kωi2/τi,η^==∑i=1kx^ik−x^i−1k−Dωi2∑i=1kτi,σ^2=1k∑i=1kητi2+ωi22ητiD+D2+Γi−2ητi+ωi2Dx^ik−x^i−1kτi,ϕ^2=1k∑i=1kyi2+Pik+x^ik2−2yix^ik,where ωi2 denotes the smoothed probability of the model formulated by equation (3). Note that the updated equation of parameter D contains parameter η, the updated equation of parameter η contains parameter D, and the updated equation of parameter σ contains parameters η and D, which makes it impossible to directly use equation (42) to obtain the estimated values of parameters D, η, and σ. For this reason, the simplex method of multidimensional search in [36] is utilized to achieve the estimated values of three parameters, which is integrated as the “fminsearch” function in MATLAB [46]. The function is an effective method to search for the minimum value of the multidimensional function. The specific solving process is as follows. Equation (42) is taken into equation (31) to get an expression that only contains the parameters D and η. Then, the “fminsearch” function is used to perform a two-dimensional search starting from the initial value of the parameters D and η. When the expression −Qθ,θkj gets the minimum value, the expression Qθ,θkj obtains the maximum value. In this case, the corresponding parameter values are the estimated values of D and η. Finally, the estimated values of D and η are brought into the updated equation of the parameter σ in equation (42), and the estimated values of σ is achieved.The abovementioned solving method can obtain the estimated values of model parameters at thejth iteration of EM algorithm. In other words, the one iteration from θ^kj−1 to θ^kj is completed. Then, the estimated values are taken into IMM filtering-based model soft switch to update the model probabilities and the simulated degradation state and EM algorithm is performed again until a criterion of convergence is satisfied.


## 3.1. IMM Filtering-Based Model Soft Switch

In the HWSSM, the unknown characteristic of the shocks make it impossible to know whether there is a shock arrival from the observed data, which leads to no way to inform the conditions of switching the simulation model. Consequently, it is necessary to calculate the probabilities of different simulation models. In the SSM modeling method, the IMM filtering provides an effective way for calculating the probability of different models, and the Kalman filter is utilized in IMM filtering. The simulation model soft switch includes four stages, i.e., model input interaction, Kalman filtering, model probability calculation, and model output interaction.For convenience, the following definitions are given.Yk=y1,y2,…,yk and Xk=x0,x1,x2,…,xk represent the observation data vector and simulated degradation state vector until time k, respectively. m represents the number of simulation models. μku represents the probability of simulation model u at time k. pvu represents the transition probability from the model v to the model u. P=pvum×m represents the Markov probability transition matrix. Cku denotes that the valid model is model u in time interval k−1,k. x^kku denotes the estimated mean of the simulated degradation state at time k conditioned on Cku and the previous k measured values. Pkku denotes the estimated covariance of the simulated degradation state at time k conditioned on Cku and the previous k measured values. x^kk0u denotes the hybrid estimated mean of the simulated degradation state at time k conditioned on Ck+1u and the previous k measured values. Pkk0u denotes the hybrid estimated covariance of the simulated degradation state at time k conditioned on Ck+1u and the previous k measured values. Nx;0,κ denotes the Gaussian distribution with value x, mean 0, and variance κ.


### 3.1.1. Input Interaction (Modelu)

At this stage, the hybrid state estimationx^k−1k−10u and the covariance estimation Pk−1k−10u are obtained by the state estimation x^k−1k−1v and model probability μk−1v of model v at time k−1. x^k−1k−10u and Pk−1k−10u are both used as the initial state of the Kalman filtering at the time k. The specific steps are as follows.The prior probabilityG¯u of the simulation model u is defined and calculated by(8)G¯u=ΔpCkuYk−1=∑v=1mpCkuCk−1v,Yk−1pCk−1vYk−1=∑v=1mpvuμk−1v,where pvu=ΔPCkuCk−1v,Yk−1.The hybrid probabilityμk−1k−1vu transferred from the model v to the model u is defined and calculated by(9)μk−1k−1vu=ΔpCk−1vCku,Yk−1=pCkuCk−1v,Yk−1pCk−1vYk−1pCkuYk−1=pvuμk−1vG¯u.The hybrid state estimationx^k−1k−10u of simulation model u is determined by(10)x^k−1k−10u=∑v=1mx^k−1k−1vμk−1k−1vu.The hybrid covariance estimationPk−1k−10u of simulation model u is determined by(11)Pk−1k−10u=∑v=1mμk−1k−1vuPk−1k−1v+x^k−1k−1v−x^k−1k−10u⋅x^k−1k−1v−x^k−1k−10uT.


### 3.1.2. Kalman Filtering (Modelu)

Kalman filtering can be divided into two steps, i.e., the prediction step and the updating step.x^kk−1u and Pkk−1u express the simulated degradation state estimation and covariance estimation, respectively, conditioned on Cku and the previous k−1 measured values. Then, the prediction step refers to achieve the predicted results x^kk−1u and Pkk−1u according to the degradation equations and the results of input interaction. The prediction step can be expressed as(12)x^kk−11=x^k−1k−101+ητk,(13)x^kk−12=x^k−1k−102+ητk+D,(14)Pkk−1u=Pk−1k−10u+σ2τk.The updating step refers to obtain the posterior state estimationx^kku and covariance estimation Pkku of model u based on the prognostic results in the prediction step and the observed data yk. The updating step can be expressed as(15)y˜ku=yk−x^kk−1u,Sku=Pkk−1u+ϕ2,Kku=Pkk−1uSku−1,x^kku=x^kk−1u+Kkuy˜ku,Pkku=I−KkuPkk−1u,where y˜ku denotes the new information, Sku is the variance of y˜ku, Kku is the Kalman gain, and I represents one-dimensional unit matrix.


### 3.1.3. Model Probability Calculation

At this stage, the likelihood functionΛku of the model u is used to calculate the model probability and is defined as(16)Λku=ΔpykCku,Yk−1≈Nyk;x^kk−1u,Pkk−1u+ϕ2.Furthermore, according to equation (16), Λku can be given by(17)Λku=12π1/2Sku1/2exp−12y˜kuTSku−1y˜ku.Then, the probability of the simulation modelu is(18)μku=ΔpCkuYk=ΛkuG¯uGu,where Gu is a normalization constant that satisfies Gu=∑u=1mΛkuG¯u.


### 3.1.4. Output Interaction

According to the model probabilityμku and the estimation results of each model in the Kalman filtering stage, the degradation state estimation x^kk and covariance estimation Pkk are achieved through the weighted sum and can be given by(19)x^kk=∑u=1mx^kkuμku,(20)Pkk=∑u=1mμkuPkku+x^kku−x^kkx^kku−x^kkT.In addition, in the initial time of degraded state estimation, there existsx00u=x00v and pC0u=pC0v conditioned on ∀u,v∈m,u≠v. Once a new measurement value yk+1 is obtained, the model probability can be updated to realize model soft switch according to equations (8)–(20).


## 3.1.1. Input Interaction (Modelu)

At this stage, the hybrid state estimationx^k−1k−10u and the covariance estimation Pk−1k−10u are obtained by the state estimation x^k−1k−1v and model probability μk−1v of model v at time k−1. x^k−1k−10u and Pk−1k−10u are both used as the initial state of the Kalman filtering at the time k. The specific steps are as follows.The prior probabilityG¯u of the simulation model u is defined and calculated by(8)G¯u=ΔpCkuYk−1=∑v=1mpCkuCk−1v,Yk−1pCk−1vYk−1=∑v=1mpvuμk−1v,where pvu=ΔPCkuCk−1v,Yk−1.The hybrid probabilityμk−1k−1vu transferred from the model v to the model u is defined and calculated by(9)μk−1k−1vu=ΔpCk−1vCku,Yk−1=pCkuCk−1v,Yk−1pCk−1vYk−1pCkuYk−1=pvuμk−1vG¯u.The hybrid state estimationx^k−1k−10u of simulation model u is determined by(10)x^k−1k−10u=∑v=1mx^k−1k−1vμk−1k−1vu.The hybrid covariance estimationPk−1k−10u of simulation model u is determined by(11)Pk−1k−10u=∑v=1mμk−1k−1vuPk−1k−1v+x^k−1k−1v−x^k−1k−10u⋅x^k−1k−1v−x^k−1k−10uT.


## 3.1.2. Kalman Filtering (Modelu)

Kalman filtering can be divided into two steps, i.e., the prediction step and the updating step.x^kk−1u and Pkk−1u express the simulated degradation state estimation and covariance estimation, respectively, conditioned on Cku and the previous k−1 measured values. Then, the prediction step refers to achieve the predicted results x^kk−1u and Pkk−1u according to the degradation equations and the results of input interaction. The prediction step can be expressed as(12)x^kk−11=x^k−1k−101+ητk,(13)x^kk−12=x^k−1k−102+ητk+D,(14)Pkk−1u=Pk−1k−10u+σ2τk.The updating step refers to obtain the posterior state estimationx^kku and covariance estimation Pkku of model u based on the prognostic results in the prediction step and the observed data yk. The updating step can be expressed as(15)y˜ku=yk−x^kk−1u,Sku=Pkk−1u+ϕ2,Kku=Pkk−1uSku−1,x^kku=x^kk−1u+Kkuy˜ku,Pkku=I−KkuPkk−1u,where y˜ku denotes the new information, Sku is the variance of y˜ku, Kku is the Kalman gain, and I represents one-dimensional unit matrix.


## 3.1.3. Model Probability Calculation

At this stage, the likelihood functionΛku of the model u is used to calculate the model probability and is defined as(16)Λku=ΔpykCku,Yk−1≈Nyk;x^kk−1u,Pkk−1u+ϕ2.Furthermore, according to equation (16), Λku can be given by(17)Λku=12π1/2Sku1/2exp−12y˜kuTSku−1y˜ku.Then, the probability of the simulation modelu is(18)μku=ΔpCkuYk=ΛkuG¯uGu,where Gu is a normalization constant that satisfies Gu=∑u=1mΛkuG¯u.


## 3.1.4. Output Interaction

According to the model probabilityμku and the estimation results of each model in the Kalman filtering stage, the degradation state estimation x^kk and covariance estimation Pkk are achieved through the weighted sum and can be given by(19)x^kk=∑u=1mx^kkuμku,(20)Pkk=∑u=1mμkuPkku+x^kku−x^kkx^kku−x^kkT.In addition, in the initial time of degraded state estimation, there existsx00u=x00v and pC0u=pC0v conditioned on ∀u,v∈m,u≠v. Once a new measurement value yk+1 is obtained, the model probability can be updated to realize model soft switch according to equations (8)–(20).


## 3.2. Model Parameter Evolution

The model parametersθ of HWSSM include μ0, Σ0, η, σ, D, and ϕ, where μ0 and Σ0 denote the mean and covariance of the initial simulated degradation state, respectively. Since HWSSM is a state space model with hidden state, it is considered to use the expectation maximization (EM) algorithm to achieve model parameter evolution. EM algorithm is suitable for estimating model parameters of SSM with hidden state [42, 43]. Based on the model soft switch, this paper uses EM algorithm to realize the evolution of model parameters.According to the EM algorithm, the estimations of the model parametersθ at time k and the jth iteration of EM algorithm can be obtained by(21)θ^kj=argmaxθEXk,R˜1:kYk,θ^kj−1LXk,Yk,R˜1:kθYk,θ^kj−1,where E⋅ denotes the mathematical expectation operator, LXk,Yk,R˜1:kθ is the joint log-likelihood function, and R˜1:k represents the model indicating matrix. R˜k represents the model indicating vector at time k and satisfies R˜k=R˜k1,R˜k2. Its element R˜kv is defined as(22)R˜kv=1,conditioned onCkv,0,others.EM algorithm consists of two steps, i.e., E-step and M-step. E-step refers to calculate the mathematical expectation of the joint log-likelihood function, and M-step refers to maximize the mathematical expectation of the joint log-likelihood function. The model parameters can be estimated online via the iteration of E-step and M-step.


### 3.2.1. E-Step

Based on the Markov property and the multiplicative equation of conditional probability, the joint log-likelihood function of the simulated state vectorXk, the observed data vector Yk, and the model indicating matrix R˜1:k at time k is given as(23)LXk,Yk,R˜1:kθ=lnpXk,Yk,R˜1:kθ=lnpYkXk,θpXk,R˜1:kθ=lnpYkXk,θ︸①+lnpXk,R˜1:kθ︸②.According to HWSSM, there exists(24)pyixi,θ=Nyi;xi,ϕ2;pxixi−1,R˜iv,R˜1:i−1,θ=Nxi;xi−1+χiv,σ2τi;px0θ=Nx0;μ0,Σ0,where χi1=ητi and χi2=ητi+D. Then, the derivation of equation (23) is as follows:(25)①=ln∏i=1kpyixi,θ=∑i=1klnpyixi,θ=∑i=1kln12πϕ2exp−yi−xi22ϕ2=−k2ln2π−k2lnϕ2−12ϕ2∑i=1kyi−xi2,②=ln∏i=2kpxi,R˜ixi−1,R˜1:i−1,θpx1,R˜1x0,θpx0θ=lnpx0θ∏i=1k∏v=12μivNxi;xi−1+χiv,σ2τiR˜iv=lnpx0θ︸③+ln∏i=1k∏v=12μivNxi;xi−1+χiv,σ2τiR˜iv︸④,③=lnpx0θ=−12ln2π−12lnΣ0−x0−μ022Σ0,④=∑i=1k∑v=12lnμivNxi;xi−1+χiv,σ2τiR˜iv=∑i=1k∑v=12R˜ivlnμivNxi;xi−1+χiv,σ2τi=∑v=12∑i=1kR˜ivlnμiv+lnNxi;xi−1+χiv,σ2τi=∑v=12∑i=1kR˜ivlnμiv+ln12πτi−lnσ−12σ2τixi−xi−1+χiv2.Based on the above derivation, the final expression of the joint log-likelihood function is formulated as(26)LXk,Yk,R˜1:kθ=−k2ln2π−k2lnϕ2−12ϕ2∑i=1kyi−xi2+∑v=12∑i=1kR˜ivlnμiv+ln12πτi−lnσ−12σ2τixi−xi−1+χiv2−12ln2π−12lnΣ0−x0−μ022Σ0.Next,Qθ,θkj is denoted as the mathematical expectation of the joint log-likelihood function at time k and the jth iteration of EM algorithm, i.e.,(27)Qθ,θkj=EXk,R˜1:kYk,θ^kj−1LXk,Yk,R˜1:kθYk,θ^kj−1.Through neglecting the irrelevant items that are independent of the parameterθ, Qθ,θkj can be expressed as(28)Qθ,θkj∝−∑i=1k12lnϕ2+12ϕ2EXk,R˜1:k⋅yi−xi2⋅︸⑤+∑v=12∑i=1kEXk,R˜1:k⋅R˜iv⋅︸⑥lnμiv+ln12πτi−lnσ−EXk,R˜1:k⋅xi−xi−1−χiv2⋅︸⑦2σ2τi−12lnΣ0−12Σ0EXk,R˜1:k⋅x0−μ02Yk,θ^kj−1︸⑧,where ⋅ represents Yk,θ^kj−1. According to the property of mathematical expectation, the following intermediate variables are defined as(29)x^ik=EXk,R˜1:k⋅xi⋅,Pik=varXk,R˜1:k⋅xi2⋅=EXk,R˜1:k⋅xi2⋅−EXk,R˜1:k⋅xi⋅2=EXk,R˜1:k⋅xi2⋅−x^ik2,Pi,i−1k=covXk,R˜1:k⋅xixi−1⋅=EXk,R˜1:k⋅xixi−1⋅−EXk,R˜1:k⋅xi⋅EXk,R˜1:k⋅xi−1⋅=EXk,R˜1:k⋅xixi−1⋅−x^ikx^i−1k.Based on the definition of the above variables, the derivation of the four conditional mathematical expectations in equation (28) is as follows:(30)⑤=EXk,R˜1:k⋅yi2+xi2−2yixi⋅=yi2+EXk,R˜1:k⋅xi2⋅−2yiEXk,R˜1:k⋅xi⋅=yi2+Pik+x^ik2−2yix^ik,⑥=EXk,R˜1:k⋅R˜iv⋅=pR˜iv=1⋅=ωiv,⑦=EXk,R˜1:k⋅xi−xi−1−χiv2⋅=EXk,R˜1:k⋅χiv2+xi−xi−12−2χivxi−xi−1⋅=χiv2+x^ik2+Pik+x^i−1k2+Pi−1k−2x^ikx^i−1k+Pi,i−1k−2χivx^ik−x^i−1k,⑧=EXk,R˜1:k⋅x0−μ02⋅=EXk,R˜1:k⋅x02+μ02−2μ0x0⋅=x^0k2+P0k+μ02−2μ0x^0k.After deriving the above four parts of the conditional mathematical expectation,Qθ,θkj can be given by(31)Qθ,θkj∝−∑i=1k12lnϕ2+12ϕ2yi2+Pik+x^ik2−2yix^ik+∑v=12∑i=1kωivlnμiv+ln12πτi−lnσ−χiv2+Γi−2χivx^ik−x^i−1k2σ2τi−12lnΣ0−12Σ0x^0k2+P0k+μ02−2μ0x^0k,where Γi=x^ik2+Pik+x^i−1k2+Pi−1k−2x^ikx^i−1k+Pi,i−1k. According to equation (31), in order to calculate Qθ,θkj, the values of the variables x^ik, x^i−1k, Pik, Pi−1k, Pi,i−1k, and ωiv should be acquired, which belong to the smoothed variables. This paper uses the IMM backward smoothing (IMMBS) algorithm to obtain the values of the above six smoothed variables [44, 45].IMMBS algorithm consists of backward filtering and model state fusion. Backward filtering refers to the backward filtering from the latest measured value. The operation flow of backward filtering is similar to the forward IMM filtering mentioned in Section3.1, but there are also obvious differences between them. The forward IMM filtering performs the input interaction before performing one-step prediction, while the backward filtering performs one-step prediction before performing the input interaction. Backward filtering can be divided into five steps, including backward one-step prediction, backward input interaction, backward filtering update, backward model probability calculation, and backward output fusion. In particular, the backward one-step prediction equation for backward filtering is given by(32)x^ii+1B,1=x^i+1i+1B,1−ητi+1,x^ii+1B,2=x^i+1i+1B,2−ητi+1−D,Pii+1B,u=Pi+1i+1B,u+σ2τi+1,where x^ii+1B,u denotes the backward one-step prediction of the model u and Pii+1B,u represents the covariance of backward one-step prediction error for the model u. The other steps of backward filtering can be found in [44]. In the stage of model state fusion for the IMMBS algorithm, according to the full probability theorem, the smoothing estimation of the simulated state can be expressed as(33)pxiYk=∑u=1mωiupxiCiu,Yk=∑u=1mωiu∑v=1mω˜i+1kvupxiCiu,Ci+1v,Yk=∑u=1mωiu∑v=1mω˜i+1kvuΛiuv,where ωiu represents the smoothed model probability, ω˜i+1kvu denotes the smoothed hybrid probability, and Λiuv denotes the likelihood function. The definition and calculation of the three variables are given by(34)ω˜i+1kvu=ΔpCi+1vCiu,Y1k=1G˜upCi+1vCiu,YipYi+1kCi+1v,Ciu,Yi=1G˜upuvpYi+1kCi+1v,Ciu,Yi,Λiuv=ΔpYi+1kCi+1v,Ciu,Yi≈pxii+1B,vCi+1v,Ciu,xiiu=Nx^ii+1B,v−x^iiu;0,Pii+1B,v+Piiu,ωiu=ΔpCiuY1k=G˜uGpCiuYi=G˜uGμiu,where G˜u and G are all the normalized constants. They can be acquired by(35)G˜u=∑v=1mpuvΛiuv,G=∑u=1mG˜uμiu.In equation (33), with respect to pxiCiu,Ci+1v,Yk, the mixed smoothing state estimation x^ikuv, covariance estimation Pikuv, and the mixed interaction covariance estimation Pi,i−1kuv can be obtained by(36)x^ikuv=PikuvPiiu−1x^iiu+Pii+1B,v−1x^ii+1B,v,Pikuv=Piiu−1+Pii+1B,v−1−1,Pi,i−1kuv=Piiu−1+Pii+1B,v−1−1+Pi−1i−1u−1+Pi−1iB,v−1−1−1.In equation (33), with respect to pxiCiu,Yk, the mixed smoothing state estimation x^iku, covariance estimation Piku, and the mixed interaction covariance estimation Pi,i−1ku of model u can be obtained by(37)x^iku=∑u=1mω˜i+1kvux^ikuv,Piku=∑u=1mω˜i+1kvuPikuv+x^ikuv−x^ikux^ikuv−x^ikuT,Pi,i−1ku=∑u=1mω˜i+1kvuPi,i−1kuv+x^ikuv−x^ikux^i−1kuv−x^i−1kuT.In particular, ifi=k, there exists(38)Pk,k−1ku=I−KkuPk−1k−1u.Finally, after the smoothed state fusion of each model, the smooth estimationx^ik, covariance estimation Pik, and interactive covariance estimation Pi,i−1k of the simulated degradation state can be acquired by(39)x^ik=∑u=1mωiux^iku,Pik=∑u=1mωiuPiku+x^iku−x^ikx^iku−x^ikT,Pi,i−1k=∑u=1mωiuPi,i−1ku+x^iku−x^ikx^i−1ku−x^i−1kT.In particular, the initial smoothed estimations are given by(40)x^0k=∑u=1mω1ux^01B,u,P0k=∑u=1mω1uP01B,u+x^01B,u−x^0kx^01B,u−x^0kT.


### 3.2.2. M-Step

The estimation of the model parametersθ at the jth step of the EM algorithm can be obtained by differentiating equation (31), i.e.,(41)∂Qθ,θkj∂θ=0.The updated parameters are given by(42)μ^0=x^0k,Σ^0=P0k,D^=∑i=1kωi2x^ik−x^i−1k−ητi/τi∑i=1kωi2/τi,η^==∑i=1kx^ik−x^i−1k−Dωi2∑i=1kτi,σ^2=1k∑i=1kητi2+ωi22ητiD+D2+Γi−2ητi+ωi2Dx^ik−x^i−1kτi,ϕ^2=1k∑i=1kyi2+Pik+x^ik2−2yix^ik,where ωi2 denotes the smoothed probability of the model formulated by equation (3). Note that the updated equation of parameter D contains parameter η, the updated equation of parameter η contains parameter D, and the updated equation of parameter σ contains parameters η and D, which makes it impossible to directly use equation (42) to obtain the estimated values of parameters D, η, and σ. For this reason, the simplex method of multidimensional search in [36] is utilized to achieve the estimated values of three parameters, which is integrated as the “fminsearch” function in MATLAB [46]. The function is an effective method to search for the minimum value of the multidimensional function. The specific solving process is as follows. Equation (42) is taken into equation (31) to get an expression that only contains the parameters D and η. Then, the “fminsearch” function is used to perform a two-dimensional search starting from the initial value of the parameters D and η. When the expression −Qθ,θkj gets the minimum value, the expression Qθ,θkj obtains the maximum value. In this case, the corresponding parameter values are the estimated values of D and η. Finally, the estimated values of D and η are brought into the updated equation of the parameter σ in equation (42), and the estimated values of σ is achieved.The abovementioned solving method can obtain the estimated values of model parameters at thejth iteration of EM algorithm. In other words, the one iteration from θ^kj−1 to θ^kj is completed. Then, the estimated values are taken into IMM filtering-based model soft switch to update the model probabilities and the simulated degradation state and EM algorithm is performed again until a criterion of convergence is satisfied.


## 3.2.1. E-Step

Based on the Markov property and the multiplicative equation of conditional probability, the joint log-likelihood function of the simulated state vectorXk, the observed data vector Yk, and the model indicating matrix R˜1:k at time k is given as(23)LXk,Yk,R˜1:kθ=lnpXk,Yk,R˜1:kθ=lnpYkXk,θpXk,R˜1:kθ=lnpYkXk,θ︸①+lnpXk,R˜1:kθ︸②.According to HWSSM, there exists(24)pyixi,θ=Nyi;xi,ϕ2;pxixi−1,R˜iv,R˜1:i−1,θ=Nxi;xi−1+χiv,σ2τi;px0θ=Nx0;μ0,Σ0,where χi1=ητi and χi2=ητi+D. Then, the derivation of equation (23) is as follows:(25)①=ln∏i=1kpyixi,θ=∑i=1klnpyixi,θ=∑i=1kln12πϕ2exp−yi−xi22ϕ2=−k2ln2π−k2lnϕ2−12ϕ2∑i=1kyi−xi2,②=ln∏i=2kpxi,R˜ixi−1,R˜1:i−1,θpx1,R˜1x0,θpx0θ=lnpx0θ∏i=1k∏v=12μivNxi;xi−1+χiv,σ2τiR˜iv=lnpx0θ︸③+ln∏i=1k∏v=12μivNxi;xi−1+χiv,σ2τiR˜iv︸④,③=lnpx0θ=−12ln2π−12lnΣ0−x0−μ022Σ0,④=∑i=1k∑v=12lnμivNxi;xi−1+χiv,σ2τiR˜iv=∑i=1k∑v=12R˜ivlnμivNxi;xi−1+χiv,σ2τi=∑v=12∑i=1kR˜ivlnμiv+lnNxi;xi−1+χiv,σ2τi=∑v=12∑i=1kR˜ivlnμiv+ln12πτi−lnσ−12σ2τixi−xi−1+χiv2.Based on the above derivation, the final expression of the joint log-likelihood function is formulated as(26)LXk,Yk,R˜1:kθ=−k2ln2π−k2lnϕ2−12ϕ2∑i=1kyi−xi2+∑v=12∑i=1kR˜ivlnμiv+ln12πτi−lnσ−12σ2τixi−xi−1+χiv2−12ln2π−12lnΣ0−x0−μ022Σ0.Next,Qθ,θkj is denoted as the mathematical expectation of the joint log-likelihood function at time k and the jth iteration of EM algorithm, i.e.,(27)Qθ,θkj=EXk,R˜1:kYk,θ^kj−1LXk,Yk,R˜1:kθYk,θ^kj−1.Through neglecting the irrelevant items that are independent of the parameterθ, Qθ,θkj can be expressed as(28)Qθ,θkj∝−∑i=1k12lnϕ2+12ϕ2EXk,R˜1:k⋅yi−xi2⋅︸⑤+∑v=12∑i=1kEXk,R˜1:k⋅R˜iv⋅︸⑥lnμiv+ln12πτi−lnσ−EXk,R˜1:k⋅xi−xi−1−χiv2⋅︸⑦2σ2τi−12lnΣ0−12Σ0EXk,R˜1:k⋅x0−μ02Yk,θ^kj−1︸⑧,where ⋅ represents Yk,θ^kj−1. According to the property of mathematical expectation, the following intermediate variables are defined as(29)x^ik=EXk,R˜1:k⋅xi⋅,Pik=varXk,R˜1:k⋅xi2⋅=EXk,R˜1:k⋅xi2⋅−EXk,R˜1:k⋅xi⋅2=EXk,R˜1:k⋅xi2⋅−x^ik2,Pi,i−1k=covXk,R˜1:k⋅xixi−1⋅=EXk,R˜1:k⋅xixi−1⋅−EXk,R˜1:k⋅xi⋅EXk,R˜1:k⋅xi−1⋅=EXk,R˜1:k⋅xixi−1⋅−x^ikx^i−1k.Based on the definition of the above variables, the derivation of the four conditional mathematical expectations in equation (28) is as follows:(30)⑤=EXk,R˜1:k⋅yi2+xi2−2yixi⋅=yi2+EXk,R˜1:k⋅xi2⋅−2yiEXk,R˜1:k⋅xi⋅=yi2+Pik+x^ik2−2yix^ik,⑥=EXk,R˜1:k⋅R˜iv⋅=pR˜iv=1⋅=ωiv,⑦=EXk,R˜1:k⋅xi−xi−1−χiv2⋅=EXk,R˜1:k⋅χiv2+xi−xi−12−2χivxi−xi−1⋅=χiv2+x^ik2+Pik+x^i−1k2+Pi−1k−2x^ikx^i−1k+Pi,i−1k−2χivx^ik−x^i−1k,⑧=EXk,R˜1:k⋅x0−μ02⋅=EXk,R˜1:k⋅x02+μ02−2μ0x0⋅=x^0k2+P0k+μ02−2μ0x^0k.After deriving the above four parts of the conditional mathematical expectation,Qθ,θkj can be given by(31)Qθ,θkj∝−∑i=1k12lnϕ2+12ϕ2yi2+Pik+x^ik2−2yix^ik+∑v=12∑i=1kωivlnμiv+ln12πτi−lnσ−χiv2+Γi−2χivx^ik−x^i−1k2σ2τi−12lnΣ0−12Σ0x^0k2+P0k+μ02−2μ0x^0k,where Γi=x^ik2+Pik+x^i−1k2+Pi−1k−2x^ikx^i−1k+Pi,i−1k. According to equation (31), in order to calculate Qθ,θkj, the values of the variables x^ik, x^i−1k, Pik, Pi−1k, Pi,i−1k, and ωiv should be acquired, which belong to the smoothed variables. This paper uses the IMM backward smoothing (IMMBS) algorithm to obtain the values of the above six smoothed variables [44, 45].IMMBS algorithm consists of backward filtering and model state fusion. Backward filtering refers to the backward filtering from the latest measured value. The operation flow of backward filtering is similar to the forward IMM filtering mentioned in Section3.1, but there are also obvious differences between them. The forward IMM filtering performs the input interaction before performing one-step prediction, while the backward filtering performs one-step prediction before performing the input interaction. Backward filtering can be divided into five steps, including backward one-step prediction, backward input interaction, backward filtering update, backward model probability calculation, and backward output fusion. In particular, the backward one-step prediction equation for backward filtering is given by(32)x^ii+1B,1=x^i+1i+1B,1−ητi+1,x^ii+1B,2=x^i+1i+1B,2−ητi+1−D,Pii+1B,u=Pi+1i+1B,u+σ2τi+1,where x^ii+1B,u denotes the backward one-step prediction of the model u and Pii+1B,u represents the covariance of backward one-step prediction error for the model u. The other steps of backward filtering can be found in [44]. In the stage of model state fusion for the IMMBS algorithm, according to the full probability theorem, the smoothing estimation of the simulated state can be expressed as(33)pxiYk=∑u=1mωiupxiCiu,Yk=∑u=1mωiu∑v=1mω˜i+1kvupxiCiu,Ci+1v,Yk=∑u=1mωiu∑v=1mω˜i+1kvuΛiuv,where ωiu represents the smoothed model probability, ω˜i+1kvu denotes the smoothed hybrid probability, and Λiuv denotes the likelihood function. The definition and calculation of the three variables are given by(34)ω˜i+1kvu=ΔpCi+1vCiu,Y1k=1G˜upCi+1vCiu,YipYi+1kCi+1v,Ciu,Yi=1G˜upuvpYi+1kCi+1v,Ciu,Yi,Λiuv=ΔpYi+1kCi+1v,Ciu,Yi≈pxii+1B,vCi+1v,Ciu,xiiu=Nx^ii+1B,v−x^iiu;0,Pii+1B,v+Piiu,ωiu=ΔpCiuY1k=G˜uGpCiuYi=G˜uGμiu,where G˜u and G are all the normalized constants. They can be acquired by(35)G˜u=∑v=1mpuvΛiuv,G=∑u=1mG˜uμiu.In equation (33), with respect to pxiCiu,Ci+1v,Yk, the mixed smoothing state estimation x^ikuv, covariance estimation Pikuv, and the mixed interaction covariance estimation Pi,i−1kuv can be obtained by(36)x^ikuv=PikuvPiiu−1x^iiu+Pii+1B,v−1x^ii+1B,v,Pikuv=Piiu−1+Pii+1B,v−1−1,Pi,i−1kuv=Piiu−1+Pii+1B,v−1−1+Pi−1i−1u−1+Pi−1iB,v−1−1−1.In equation (33), with respect to pxiCiu,Yk, the mixed smoothing state estimation x^iku, covariance estimation Piku, and the mixed interaction covariance estimation Pi,i−1ku of model u can be obtained by(37)x^iku=∑u=1mω˜i+1kvux^ikuv,Piku=∑u=1mω˜i+1kvuPikuv+x^ikuv−x^ikux^ikuv−x^ikuT,Pi,i−1ku=∑u=1mω˜i+1kvuPi,i−1kuv+x^ikuv−x^ikux^i−1kuv−x^i−1kuT.In particular, ifi=k, there exists(38)Pk,k−1ku=I−KkuPk−1k−1u.Finally, after the smoothed state fusion of each model, the smooth estimationx^ik, covariance estimation Pik, and interactive covariance estimation Pi,i−1k of the simulated degradation state can be acquired by(39)x^ik=∑u=1mωiux^iku,Pik=∑u=1mωiuPiku+x^iku−x^ikx^iku−x^ikT,Pi,i−1k=∑u=1mωiuPi,i−1ku+x^iku−x^ikx^i−1ku−x^i−1kT.In particular, the initial smoothed estimations are given by(40)x^0k=∑u=1mω1ux^01B,u,P0k=∑u=1mω1uP01B,u+x^01B,u−x^0kx^01B,u−x^0kT.


## 3.2.2. M-Step

The estimation of the model parametersθ at the jth step of the EM algorithm can be obtained by differentiating equation (31), i.e.,(41)∂Qθ,θkj∂θ=0.The updated parameters are given by(42)μ^0=x^0k,Σ^0=P0k,D^=∑i=1kωi2x^ik−x^i−1k−ητi/τi∑i=1kωi2/τi,η^==∑i=1kx^ik−x^i−1k−Dωi2∑i=1kτi,σ^2=1k∑i=1kητi2+ωi22ητiD+D2+Γi−2ητi+ωi2Dx^ik−x^i−1kτi,ϕ^2=1k∑i=1kyi2+Pik+x^ik2−2yix^ik,where ωi2 denotes the smoothed probability of the model formulated by equation (3). Note that the updated equation of parameter D contains parameter η, the updated equation of parameter η contains parameter D, and the updated equation of parameter σ contains parameters η and D, which makes it impossible to directly use equation (42) to obtain the estimated values of parameters D, η, and σ. For this reason, the simplex method of multidimensional search in [36] is utilized to achieve the estimated values of three parameters, which is integrated as the “fminsearch” function in MATLAB [46]. The function is an effective method to search for the minimum value of the multidimensional function. The specific solving process is as follows. Equation (42) is taken into equation (31) to get an expression that only contains the parameters D and η. Then, the “fminsearch” function is used to perform a two-dimensional search starting from the initial value of the parameters D and η. When the expression −Qθ,θkj gets the minimum value, the expression Qθ,θkj obtains the maximum value. In this case, the corresponding parameter values are the estimated values of D and η. Finally, the estimated values of D and η are brought into the updated equation of the parameter σ in equation (42), and the estimated values of σ is achieved.The abovementioned solving method can obtain the estimated values of model parameters at thejth iteration of EM algorithm. In other words, the one iteration from θ^kj−1 to θ^kj is completed. Then, the estimated values are taken into IMM filtering-based model soft switch to update the model probabilities and the simulated degradation state and EM algorithm is performed again until a criterion of convergence is satisfied.


## 4. Parallel Simulation-Based Equipment RUL Real-Time Prediction

According to the widely used concept of first hitting time (FHT), the definition of the equipment RUL is given. Assuming that the equipment failure threshold isw, the equipment RUL T is defined as the time when the degradation process first passes the failure threshold [47], i.e.,(43)Tw=inft:xt≥wx0<w.In order to obtain the analytical expression of the PDF of RUL, the derivation of the RUL calculation at timek is performed below. Considering that there may be two degraded state equations at a particular moment in the HWSSM, the RUL probability density distribution cannot be obtained directly, so the degradation process xt needs to be transformed. At the current monitoring time k, the simulated degradation state is xk. According to the Markov characteristic of the Brownian motion, the Wiener process can be rewritten as(44)xt=xk+ηt−tk+σBt−Btk=xk+ηt−tk+σBt−tk.Then the two state equations of HWSSM are merged into one equation. In other words, the damages caused byn Poisson shocks are integrated into the Wiener process. A merged state equation is expressed as(45)xt=xk+nD+ηt−tk+σBt−tk.Equation (45) is a variant of the Wiener process. According to the Wiener process property and the definition determined by equation (43), the RUL probability density function conditioned on n Poisson shocks, xk, θ, and Yk obeys the inverse Gaussian distribution, i.e.,(46)Tkn,xk,θ,Yk∼IGw−xk−nDη,w−xk−nD2σ2,where Tk represents the RUL at time k. Its mean is given by ETkn,xk,θ,Yk=w−xk−nD/η and its PDF is formulated by(47)fTkn,xk,θ,Yk=w−xk−nD2πTk3σ2exp−w−ηTk−xk−nD22σ2Tk.Equation (47) does not take into account the simulated degradation state distribution obtained by the parallel simulation system. The distribution reflects the uncertainty of the simulated degradation state. Integrating it into the RUL distribution determined by equation (47) can improve the prediction accuracy and enhance the prediction rationality. However, it will involve complex integral operation. So a lemma is utilized to obtain the probability density function of the RUL.Lemma 1.
Supposing thatΩ∼Nγ−ξ2 and A, B, and C are all constant,(48)EΩA−Ωexp−B−Ω22C=Cξ2+CA−ξ2B+γCξ2+Cexp−B−γ22ξ2+C.The proof of Lemma 1 can refer to literature [48]. On the basis of Lemma 1, a theorem is proposed to achieve the RUL distribution.Theorem 1.
For the hybrid degradation processxt,t≥0, the RUL distribution at time k can be given by(49)fTTkθ,Yk=∑n=0+∞ρTknσ2w−x^kk−nD+ηPkk2πPkk+σ2Tk3⋅n!⋅exp−ρTk−w−ηTk−x^kk−nD22Pkk+σ2Tk.According to the form of the specific HWSSM and IMM filtering algorithm,xk follows the normal distribution, i.e., xk∼Nx^kk,Pkk. Let pxkYk express the conditional PDF about Yk of xk. Based on the total probability theorem, the simulated degradation state distribution is integrated into the inverse Gaussian distribution and gives(50)fTTkθ,Yk=∫−∞+∞fTTkxk,θ,YkpxkYkdxk=∑n=0+∞∫−∞+∞fTTkn,xk,θ,YkpxkYkdxkpn=∑n=0+∞ExkfTTkn,xk,θ,Ykpn.According to Lemma 1, letA=w−nD, B=w−ηTk−nD, C=σ2Tk, Ω=xk, γ=x^kk, and ξ2=Pkk and calculate the mathematical expectation about xk of fTkxk,θ,Yk. Then, the weighted sum with the occurrence probability of the n Poisson shocks is calculated. As a result, equation (49) is obtained. This completes the proof of Theorem 1.According to the property of the mathematical expectation and equation (49), the expected value of the RUL can be obtained by(51)ETkθ,Yk=Exkθ,YkETTkxk,θ,Yk=Exkθ,YkET∑n=0+∞fTTkn,xk,θ,Ykpn.Obviously, it is difficult to directly obtain the analytical expression of RUL’s mathematical expectation by using the above equation. Therefore, according to the definition of mathematical expectation, the parallel simulation system uses numerical integration to calculate the mathematical expectation value of the RUL, i.e.,(52)ETTk=∫0+∞Tk⋅fTTkθ,YkdTk.According to the equations (49) and (52), the parallel simulation system can support maintenance decision-making by calculating the probability density function of the RUL, and its expected value in an online and real-time manner.


## 5. A Case Study

### 5.1. Data Introduction

The performance degradation of bearings of mechanical equipment is a typical hybrid degradation process with shock characteristic. This paper uses the life test data of a bearing of the IEEE PHM 2012 Prediction Competition [49] to verify the parallel simulation method considering model soft switch. These data are provided by the FEMTO-ST Institute in France. The life test is conducted on the PRONOSTIA platform shown in Figure 2(a), and the test data have been widely used in method validation in the reliability field in recent years. The life test is divided into 3 different working conditions, while the first working condition is rotating speed of 1800 rpm with a load of 4000 N. The record time of the test data is from 2010/11/17 08:33:01 to 2010/11/17 15:08:41, and the sampling frequency of the vibration signal is 25.6 kHz with the sampling interval of 10 s. A total of 2375 data samples are collected in the first working condition. In the case study, the life test data of the third bearing under the first working condition are used to validate the method. Particularly, the bearing is called bearing 1–3.Figure 2
Test platform and degradation data. (a) PRONOSTIA platform. (b) RMS of bearing 1–3.The root mean square (RMS) value of the vibration signal is a commonly used degradation feature and it is calculated by(53)RMS=1N∑i=1Nei2,where N is the sampling point number and satisfies N=2560 and ei denotes the vibration acceleration signal at the ith sampling point. The RMS of bearing 1–3 is shown in Figure 2(b). It can be seen that the shock characteristic of the performance degradation process for bearing 1–3 is obvious, indicating that the data are suitable for verifying the method. The RMS begins to change significantly after the 1500th monitoring time, which is regarded as the starting time of RUL prediction oriented parallel simulation. The failure criterion of the bearing is that the vibration intensity of the original signal reached 20 g at the 2341th monitoring time, and the corresponding RMS value is 4.7145. As a result, the failure threshold is set to the RMS value at the 2341th monitoring time, i.e., w=4.7145.


### 5.2. Model Evolution and RUL Real-Time Prediction

The initial simulation configurations includex0=0.2, η=0.02, σ=0.5, D=0.02, ρ=0.5, τ=1, and ϕ=0.1. Furthermore, the vector of initial model probability is μ0=0.60.4T, i.e., μ01=0.6 and μ02=0.4. The model transition probability matrix is chosen as P=0.50.5;0.60.4. The comparison of the degraded trajectories is shown in Figure 3. It shows that the difference between the simulated degradation trajectory and the actual degradation trajectory is minimal, indicating that the simulation output can effectively approach the actual degradation process driven by the real-time degraded data. In order to quantify the comparison results, the root mean square error (RMSE) is given by(54)RMSE=1r∑k=1ryk−x^kkyk2×100%,where rr=842 denotes the number of monitoring points. After calculation, the RMSE of the simulated degradation trajectory and the actual observed degradation trajectory is only 3.497%, which fully demonstrates that the parallel simulation considering model soft switch can effectively model and simulate the performance degradation process with discrete shocks of bearing 1–3.Figure 3
Comparison of the degraded traces.The model probability is shown in Figure4. The Poisson shocks characteristic is not significant in the monitoring period from t1500 to t1765, and the linear degradation characteristic is more obvious. The probability of the linear degradation model which is noted as Model 1 is obviously higher than that of the degradation model with discrete Poisson shocks which is noted as Model 2. The model probabilities of the two models are about 0.74 and 0.26, respectively. In this monitoring period, the dominant model is Model 1. However, as time passes, the Poisson shock characteristic becomes more and more prominent, especially at the moments t1766, t1827, t1877, t2130, t2234, etc. The probability of Model 2 generally shows a dynamic upward trend. Conversely, the probability of Model 1 shows a dynamic decline trend. In the late degradation stage, the probability of Model 2 surpasses the probability of Model 1, which shows that the former model is more suitable for describing the current degradation process. Above all, the parallel simulation method considering model soft switch can effectively meet the needs of model suitability for RUL prediction. It is worth noting that the probability curves of the two models are symmetric about the probability μ=0.5.Figure 4
Model probability.With the dynamic injection of the observed degradation data, the parallel simulation system uses the IMM-EM algorithm to perform model evolution. The estimated results of model parameters are shown in Figure5. It shows that the drift coefficient η fluctuates around 0.004 with the fluctuation range [0, 0.012]. In addition, the fluctuations are larger at the monitoring times t1766, t1827, etc., reflecting the accelerated degradation rate of bearings 1–3. The diffusion coefficient σ can converge quickly. When the shocks characteristic is obvious, σ fluctuates greatly and reaches a new convergence state, which is conducive to obtaining a stable remaining life probability density function. Furthermore, the convergence value of σ is rather small, which is helpful to obtain a narrower PDF of the RUL and improve the accuracy of RUL prediction. The damage D caused by Poisson shocks fluctuates dynamically within the interval [0.03, 0.07]. Considering the damage into the RUL prediction, it can effectively reduce the occurrence of “lack maintenance.” Parameter ϕ converges faster and the convergence value is about 0.01, reflecting that the fluctuation of measurement error is gradually stable.Figure 5
Parameters evolution of HWSSM. (a) Parameterη evolution. (b) Parameter σ evolution. (c) Parameter D evolution. (d) Parameter ϕ evolution.
(a)
(b)
(c)
(d)With the execution of model soft switch and parameters evolution, the parallel simulation system uses equations (49) and (52) to predict the RUL of bearings 1–3, including the PDF of the RUL and its mathematical expectation. Figure 6 shows the PDF curves predicted at eight different monitoring times from t1600 to t2300 with the prediction interval of 100 monitoring points.Figure 6
PDF of the RUL for bearing 1–3 at different monitoring time.According to Figure6, at each monitoring time of predicting the RUL, the PDF curve of the RUL can effectively cover the actual RUL. As the degradation data of bearing 1–3 continuously accumulates, the PDF curve of the RUL gradually narrows with the weaker “right-biased” characteristic and the stronger “normal” characteristic. The prognostic results indicate that the model matching degree is gradually improved, and the model parameters are more accurate. As a result, the uncertainty of the RUL prediction is getting smaller. This is due to the model soft switch and parameters online estimation. In addition, the error between the mathematical expectation of the RUL and the actual RUL is small. And also, the mathematical expectation of the RUL is close to the peak of the PDF curve, indicating that the uncertainty of the PDF is small, and the prognostic results can provide an important basis for maintenance decision.


### 5.3. Comparative Study

To further verify the validity of the parallel simulation method considering the model soft switch, the comparative study is performed by comparing with the method without considering the model soft switch which can be found in [50]. In that case, the single model is used to only execute model parameters evolution, i.e., the state equation equation (2) and the observation equation (6). So the method is called the single-model method. The equations of model parameters evolution and RUL prediction for the single-model method are given in Appendix. The RUL prognostic results of the single-model method and the proposed method are shown in Figure 7.Figure 7
Comparative results of the RUL prediction.The PDF curves of the single-model method are flatter than that of the proposed method, which indicates that it has stronger uncertainty. Moreover, the “right-biased” characteristic of the PDF curves of the single-model method is more obvious and there is a long “tail.” Although the PDF curves of the single-model method can also cover the true RUL, the flat distribution of the RUL and the obvious “right-biased” characteristic cause the prognostic results to be unfavorable to the maintenance decision. In addition, at each monitoring time of predicting the RUL, the mathematical expectation of the RUL obtained by the single-model method is all greater than that of the proposed method, implying the larger prediction error. On the contrary, the proposed method has better performance, and the prognostic results at the 1900th monitoring time are taken as examples for analysis. As shown in Figure8, the PDF curves obtained by the proposed method is more compact with less uncertainty, and the corresponding peak value is 2.24 × 10−3, which is bigger than the peak value 1.31 × 10−3 achieved by the single-model method. Besides, the RUL corresponding to the peak value of the proposed method and the single-model method are 292 cycles and 170 cycles, respectively. Considering that the true RUL at the 1900th monitoring time is 441 cycles, it shows that the RUL corresponding to the peak value of the proposed method is closer to the true RUL than that of the single-model method. Additionally, the mathematical expectation of the RUL obtained by the proposed method and the single-model method are 465.92 cycles and 553.67 cycles, respectively, also illustrating that the prediction error obtained by the proposed method is smaller than that of the single-model method.Figure 8
Comparative results of the RUL prediction at the 1900th monitoring time.In order to further quantify the comparison results of the RUL prediction, two quantitative indicators, the average relative accuracy and the total mean square error (TMSE), are introduced. At the monitoring timek, the definition of the relative prediction accuracy of the RUL is given by(55)RAk=1−T˜k−TkT˜k,where T˜k denotes the true RUL at time k. On the basis of RAk, the average relative accuracy MRA is defined as(56)MRA=1p∑k=1pRAk,where p represents the number of monitoring time points used to predict the RUL. According to the definition of MRA, it satisfies MRA∈0,1 and the bigger MRA indicates the higher prediction accuracy with respect to the RUL. In addition, the TMSE between the actual RUL and the mathematical expectation of the predicted RUL is defined as(57)TMSE=∑k=1pETk−T˜k2=∑k=1p∫0∞Tk−T˜k2fTTkθ,YkdTk.Based on the definition of the TMSE, the smaller TMSE implies the more accurate prediction result. The calculated comparison results are shown in Table1. According to Table 1, the MRA obtained by the proposed method is obviously larger than that of the single-model method, indicating the higher relative prediction accuracy. The TMSE of the proposed method is much smaller than that of the single-model method, implying the smaller RUL prediction error.Table 1
Comparative results of RUL prediction.
Methods
Indicators
MRA
TMSE (×104)
Single-model method
0.6427
13.5037
Proposed method
0.8607
5.4611


## 5.1. Data Introduction

The performance degradation of bearings of mechanical equipment is a typical hybrid degradation process with shock characteristic. This paper uses the life test data of a bearing of the IEEE PHM 2012 Prediction Competition [49] to verify the parallel simulation method considering model soft switch. These data are provided by the FEMTO-ST Institute in France. The life test is conducted on the PRONOSTIA platform shown in Figure 2(a), and the test data have been widely used in method validation in the reliability field in recent years. The life test is divided into 3 different working conditions, while the first working condition is rotating speed of 1800 rpm with a load of 4000 N. The record time of the test data is from 2010/11/17 08:33:01 to 2010/11/17 15:08:41, and the sampling frequency of the vibration signal is 25.6 kHz with the sampling interval of 10 s. A total of 2375 data samples are collected in the first working condition. In the case study, the life test data of the third bearing under the first working condition are used to validate the method. Particularly, the bearing is called bearing 1–3.Figure 2
Test platform and degradation data. (a) PRONOSTIA platform. (b) RMS of bearing 1–3.The root mean square (RMS) value of the vibration signal is a commonly used degradation feature and it is calculated by(53)RMS=1N∑i=1Nei2,where N is the sampling point number and satisfies N=2560 and ei denotes the vibration acceleration signal at the ith sampling point. The RMS of bearing 1–3 is shown in Figure 2(b). It can be seen that the shock characteristic of the performance degradation process for bearing 1–3 is obvious, indicating that the data are suitable for verifying the method. The RMS begins to change significantly after the 1500th monitoring time, which is regarded as the starting time of RUL prediction oriented parallel simulation. The failure criterion of the bearing is that the vibration intensity of the original signal reached 20 g at the 2341th monitoring time, and the corresponding RMS value is 4.7145. As a result, the failure threshold is set to the RMS value at the 2341th monitoring time, i.e., w=4.7145.


## 5.2. Model Evolution and RUL Real-Time Prediction

The initial simulation configurations includex0=0.2, η=0.02, σ=0.5, D=0.02, ρ=0.5, τ=1, and ϕ=0.1. Furthermore, the vector of initial model probability is μ0=0.60.4T, i.e., μ01=0.6 and μ02=0.4. The model transition probability matrix is chosen as P=0.50.5;0.60.4. The comparison of the degraded trajectories is shown in Figure 3. It shows that the difference between the simulated degradation trajectory and the actual degradation trajectory is minimal, indicating that the simulation output can effectively approach the actual degradation process driven by the real-time degraded data. In order to quantify the comparison results, the root mean square error (RMSE) is given by(54)RMSE=1r∑k=1ryk−x^kkyk2×100%,where rr=842 denotes the number of monitoring points. After calculation, the RMSE of the simulated degradation trajectory and the actual observed degradation trajectory is only 3.497%, which fully demonstrates that the parallel simulation considering model soft switch can effectively model and simulate the performance degradation process with discrete shocks of bearing 1–3.Figure 3
Comparison of the degraded traces.The model probability is shown in Figure4. The Poisson shocks characteristic is not significant in the monitoring period from t1500 to t1765, and the linear degradation characteristic is more obvious. The probability of the linear degradation model which is noted as Model 1 is obviously higher than that of the degradation model with discrete Poisson shocks which is noted as Model 2. The model probabilities of the two models are about 0.74 and 0.26, respectively. In this monitoring period, the dominant model is Model 1. However, as time passes, the Poisson shock characteristic becomes more and more prominent, especially at the moments t1766, t1827, t1877, t2130, t2234, etc. The probability of Model 2 generally shows a dynamic upward trend. Conversely, the probability of Model 1 shows a dynamic decline trend. In the late degradation stage, the probability of Model 2 surpasses the probability of Model 1, which shows that the former model is more suitable for describing the current degradation process. Above all, the parallel simulation method considering model soft switch can effectively meet the needs of model suitability for RUL prediction. It is worth noting that the probability curves of the two models are symmetric about the probability μ=0.5.Figure 4
Model probability.With the dynamic injection of the observed degradation data, the parallel simulation system uses the IMM-EM algorithm to perform model evolution. The estimated results of model parameters are shown in Figure5. It shows that the drift coefficient η fluctuates around 0.004 with the fluctuation range [0, 0.012]. In addition, the fluctuations are larger at the monitoring times t1766, t1827, etc., reflecting the accelerated degradation rate of bearings 1–3. The diffusion coefficient σ can converge quickly. When the shocks characteristic is obvious, σ fluctuates greatly and reaches a new convergence state, which is conducive to obtaining a stable remaining life probability density function. Furthermore, the convergence value of σ is rather small, which is helpful to obtain a narrower PDF of the RUL and improve the accuracy of RUL prediction. The damage D caused by Poisson shocks fluctuates dynamically within the interval [0.03, 0.07]. Considering the damage into the RUL prediction, it can effectively reduce the occurrence of “lack maintenance.” Parameter ϕ converges faster and the convergence value is about 0.01, reflecting that the fluctuation of measurement error is gradually stable.Figure 5
Parameters evolution of HWSSM. (a) Parameterη evolution. (b) Parameter σ evolution. (c) Parameter D evolution. (d) Parameter ϕ evolution.
(a)
(b)
(c)
(d)With the execution of model soft switch and parameters evolution, the parallel simulation system uses equations (49) and (52) to predict the RUL of bearings 1–3, including the PDF of the RUL and its mathematical expectation. Figure 6 shows the PDF curves predicted at eight different monitoring times from t1600 to t2300 with the prediction interval of 100 monitoring points.Figure 6
PDF of the RUL for bearing 1–3 at different monitoring time.According to Figure6, at each monitoring time of predicting the RUL, the PDF curve of the RUL can effectively cover the actual RUL. As the degradation data of bearing 1–3 continuously accumulates, the PDF curve of the RUL gradually narrows with the weaker “right-biased” characteristic and the stronger “normal” characteristic. The prognostic results indicate that the model matching degree is gradually improved, and the model parameters are more accurate. As a result, the uncertainty of the RUL prediction is getting smaller. This is due to the model soft switch and parameters online estimation. In addition, the error between the mathematical expectation of the RUL and the actual RUL is small. And also, the mathematical expectation of the RUL is close to the peak of the PDF curve, indicating that the uncertainty of the PDF is small, and the prognostic results can provide an important basis for maintenance decision.


## 5.3. Comparative Study

To further verify the validity of the parallel simulation method considering the model soft switch, the comparative study is performed by comparing with the method without considering the model soft switch which can be found in [50]. In that case, the single model is used to only execute model parameters evolution, i.e., the state equation equation (2) and the observation equation (6). So the method is called the single-model method. The equations of model parameters evolution and RUL prediction for the single-model method are given in Appendix. The RUL prognostic results of the single-model method and the proposed method are shown in Figure 7.Figure 7
Comparative results of the RUL prediction.The PDF curves of the single-model method are flatter than that of the proposed method, which indicates that it has stronger uncertainty. Moreover, the “right-biased” characteristic of the PDF curves of the single-model method is more obvious and there is a long “tail.” Although the PDF curves of the single-model method can also cover the true RUL, the flat distribution of the RUL and the obvious “right-biased” characteristic cause the prognostic results to be unfavorable to the maintenance decision. In addition, at each monitoring time of predicting the RUL, the mathematical expectation of the RUL obtained by the single-model method is all greater than that of the proposed method, implying the larger prediction error. On the contrary, the proposed method has better performance, and the prognostic results at the 1900th monitoring time are taken as examples for analysis. As shown in Figure8, the PDF curves obtained by the proposed method is more compact with less uncertainty, and the corresponding peak value is 2.24 × 10−3, which is bigger than the peak value 1.31 × 10−3 achieved by the single-model method. Besides, the RUL corresponding to the peak value of the proposed method and the single-model method are 292 cycles and 170 cycles, respectively. Considering that the true RUL at the 1900th monitoring time is 441 cycles, it shows that the RUL corresponding to the peak value of the proposed method is closer to the true RUL than that of the single-model method. Additionally, the mathematical expectation of the RUL obtained by the proposed method and the single-model method are 465.92 cycles and 553.67 cycles, respectively, also illustrating that the prediction error obtained by the proposed method is smaller than that of the single-model method.Figure 8
Comparative results of the RUL prediction at the 1900th monitoring time.In order to further quantify the comparison results of the RUL prediction, two quantitative indicators, the average relative accuracy and the total mean square error (TMSE), are introduced. At the monitoring timek, the definition of the relative prediction accuracy of the RUL is given by(55)RAk=1−T˜k−TkT˜k,where T˜k denotes the true RUL at time k. On the basis of RAk, the average relative accuracy MRA is defined as(56)MRA=1p∑k=1pRAk,where p represents the number of monitoring time points used to predict the RUL. According to the definition of MRA, it satisfies MRA∈0,1 and the bigger MRA indicates the higher prediction accuracy with respect to the RUL. In addition, the TMSE between the actual RUL and the mathematical expectation of the predicted RUL is defined as(57)TMSE=∑k=1pETk−T˜k2=∑k=1p∫0∞Tk−T˜k2fTTkθ,YkdTk.Based on the definition of the TMSE, the smaller TMSE implies the more accurate prediction result. The calculated comparison results are shown in Table1. According to Table 1, the MRA obtained by the proposed method is obviously larger than that of the single-model method, indicating the higher relative prediction accuracy. The TMSE of the proposed method is much smaller than that of the single-model method, implying the smaller RUL prediction error.Table 1
Comparative results of RUL prediction.
Methods
Indicators
MRA
TMSE (×104)
Single-model method
0.6427
13.5037
Proposed method
0.8607
5.4611


## 6. Conclusions and Future Perspectives

Equipment parallel simulation is an emerging simulation paradigm, which has an important concept of model evolution. Based on the background of RUL prediction for the hybrid degradation equipment with continuous degradation and discrete shocks, the hybrid degradation equipment RUL prediction oriented parallel simulation considering model soft switch is studied, including parallel simulation modeling, model evolution, RUL prediction, and case study. Under the modeling framework of SSM, two different models are constructed by using the Wiener process and the effect of Poisson shocks. One model describes a continuous degradation process, and the other model expresses a degradation process with discrete shocks. Owing to the discrete, unknown characteristics of the shocks, it is not possible to directly determine the model morphology at a specific time. Therefore, the forward IMM filtering is utilized to dynamically calculate the probabilities of different models, achieving the model soft switch. Then, the model probability-based weighted summation is performed to obtain the simulated degradation state estimation. On the basis of model soft switch, the evolution of model parameters based on the EM algorithm is studied. Through the iteration between model soft switch and model parameters evolution, the output of the simulation model can dynamically approach the actual degradation process. In order to realize RUL real-time prediction, the unknown Poisson shocks is firstly integrated into the Wiener process. According to the concept of the first hitting time and the mathematical property of the Wiener process, the RUL distribution which neglects the simulated degradation state distribution is achieved and subjects to the inverse Gaussian distribution. Then, based on the total probability theorem, the simulated degradation state distribution is integrated into the inverse Gaussian distribution to obtain the analytical expression of the PDF of the RUL. A bearing degradation data with typical hybrid degradation characteristic is regarded as the data-driven source to verify the proposed method considering model soft switch. The results show that the proposed method can effectively model the bearing performance degradation process. Comparative study shows that the PDF of the RUL acquired by the proposed method has a less uncertainty and higher prediction accuracy than that of the single-model method. This paper researches and perfects the modeling method, model evolution mechanism, and RUL prediction method of parallel simulation in the field of equipment RUL prediction, which is helpful to promote the application of parallel simulation in actual equipment maintenance support.There are several underlying directions deserving further implementation. First, this paper only considers the hybrid degradation with linear continuous degradation and discrete shocks. It is a challenging work further considering the hybrid degradation with nonlinear continuous degradation and discrete shocks. The mechanisms of model soft switch and parameters evolution will be more complicated and difficult. Additionally, the multistage case can be considered, and the corresponding model dynamically evolution is deserved to research.



---

*Source: 9179870-2019-03-12.xml*