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

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.


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. e 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]. e 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. e 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. e artificial intelligence methods achieve RUL prognostic results via data fitting mainly including machine learning. e 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]. e 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. e 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. ey 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. erefore, 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]. is 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. e 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. e equipment parallel simulation diagram is shown in Figure 1. In the equipment parallel simulation, the actual equipment and simulation system exchange information through sensors and actuators. e 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. e actual equipment information provided by the sensors can be divided into two categories, namely, observable state information S t and behavior information B t at time t. Equipment can be controlled by the control information C t sent by the parallel simulation system. According to the implementation of control information, it can be divided into automatic control information AC t and manual control information MC t .
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.
is 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. e 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]. e latter refers to dividing the actual complex problem into several subproblems handled simultaneously [29]. eoretical 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]. e 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. e failure of these components will cause the entire equipment to be shut down. erefore, the performance degradation process of these components is often used to measure the RUL of mechanical equipment. e 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. e structure of the paper is as follows. Based on the modeling analysis, the proposed simulation model is constructed in Section 2. e 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. en, the RUL real-time prediction method of hybrid degradation equipment based on parallel simulation is presented in Section 4. e 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.

Modeling Analysis.
e 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. erefore, 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. e equipment degradation SSM includes a degradation state equation and an observation equation. e 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. e stochastic process is suitable for describing the randomness and uncertainty of the degradation  process.
e Wiener process and Gamma process are the most commonly used stochastic processes, but the application conditions of the latter are too harsh. e 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. erefore, 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.

Construction of HWSSM.
In order to construct HWSSM, two equations need to be developed, i.e., hybrid degradation state equation and observation equation. e 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. e continuous degradation state equation can be formulated as where x(t), t ≥ 0 { } is the continuous degradation process driven by the standard Brownian motion B(t), t ≥ 0 { } and B(t)∼N(0, t); x(0) is an initial degradation state; η and σ are the drift coefficient and diffusion coefficient of the standard Brownian motion [37]. en, equation (1) is transformed by Euler discretization and the degradation state equation at discrete time points t k (k � 1, 2, . . .) without considering the unknown discrete shocks is yielded as where τ k � t k − t k−1 is the sensor sampling interval and x k � x(t k ) is the simulated degradation state at time t k . Time t k is short for time, k and ϖ k is the noise sequence which obeys the standard normal distribution. e 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 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., where M(t k ) 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 data y(t) are obtained by sensor measurement. e random relationship between the simulated degradation state and the observation data can be described by the observation equation, i.e., where π(t)∼N(0, ϕ 2 ) and ϕ 2 denotes the variance of the measurement noise. It is also assumed that π(t) is independent of the Brownian motion. en, equation (5) is transformed by Euler discretization, and the observation equation at discrete time points t k (k � 1, 2, . . .) is yielded as where y k � y(t k ) 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

HWSSM Evolution
HWSSM belongs to the state space model with hidden degradation state. e parallel simulation system realizes the evolution of HWSSM by the following two ways. e first is the model soft switch based on the interactive multiple model (IMM) filtering [40,41]. rough 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. e model parameters are updated by utilizing the latest observation data via the parameter evolution. e model soft switch and parameter evolution are not executed in isolation, but they are iterative to each other. rough 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.

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. e 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. Y k � y 1 , y 2 , . . . , y k and X k � x 0 , x 1 , x 2 , . . . , x k represent the observation data vector and simulated degradation state vector until time k, respectively. m represents the number of simulation models. e prior probability G u of the simulation model u is defined and calculated by where e hybrid state estimation x 0u k−1|k−1 of simulation model u is determined by e hybrid covariance estimation P 0u k−1|k−1 of simulation model u is determined by 3.1.2. Kalman Filtering (Model u). Kalman filtering can be divided into two steps, i.e., the prediction step and the updating step. x u k|k−1 and P u k|k−1 express the simulated degradation state estimation and covariance estimation, respectively, conditioned on C u k and the previous k − 1 measured values. en, the prediction step refers to achieve the predicted results x u k|k−1 and P u k|k−1 according to the degradation equations and the results of input interaction. e prediction step can be expressed as e updating step refers to obtain the posterior state estimation x u k|k and covariance estimation P u k|k of model u based on the prognostic results in the prediction step and the observed data y k . e updating step can be expressed as where y u k denotes the new information, S u k is the variance of y u k , K u k is the Kalman gain, and I represents one-dimensional unit matrix.

Model Probability Calculation.
At this stage, the likelihood function Λ u k of the model u is used to calculate the model probability and is defined as Furthermore, according to equation (16), Λ u k can be given by en, the probability of the simulation model u is Computational Intelligence and Neuroscience 5 where G u is a normalization constant that satisfies G u � m u�1 Λ u k G u .

Output
Interaction. According to the model probability μ u k and the estimation results of each model in the Kalman filtering stage, the degradation state estimation x k|k and covariance estimation P k|k are achieved through the weighted sum and can be given by In addition, in the initial time of degraded state estimation, there exists Once a new measurement value y k+1 is obtained, the model probability can be updated to realize model soft switch according to equations (8)-(20).

Model Parameter Evolution.
e 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 where E(·) denotes the mathematical expectation operator, L(X k , Y k , 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 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 loglikelihood function. e model parameters can be estimated online via the iteration of E-step and M-step.

E-
Step. Based on the Markov property and the multiplicative equation of conditional probability, the joint log-likelihood function of the simulated state vector X k , the observed data vector Y k , and the model indicating matrix R 1:k at time k is given as According to HWSSM, there exists where en, the derivation of equation (23) is as follows: Based on the above derivation, the final expression of the joint log-likelihood function is formulated as Next, Q(θ, θ (j) k ) is denoted as the mathematical expectation of the joint log-likelihood function at time k and the jth iteration of EM algorithm, i.e., rough neglecting the irrelevant items that are independent of the parameter θ, Q(θ, θ (j) k ) can be expressed as where · represents Y k , θ (j−1) k . According to the property of mathematical expectation, the following intermediate variables are defined as Based on the definition of the above variables, the derivation of the four conditional mathematical expectations in equation (28) is as follows: After deriving the above four parts of the conditional mathematical expectation, Q(θ, θ (j) k ) can be given by Computational Intelligence and Neuroscience where . According to equation (31), in order to calculate Q(θ, θ (j) k ), the values of the variables x i|k , x i−1|k , P i|k , P i−1|k , P i,i−1|k , and ω v i should be acquired, which belong to the smoothed variables. is 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. e operation flow of backward filtering is similar to the forward IMM filtering mentioned in Section 3.1, but there are also obvious differences between them. e 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 onestep prediction equation for backward filtering is given by where x B,u i|i+1 denotes the backward one-step prediction of the model u and P B,u i|i+1 represents the covariance of backward one-step prediction error for the model u. e 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 where ω u i represents the smoothed model probability, ω v|u i+1|k denotes the smoothed hybrid probability, and Λ uv i denotes the likelihood function. e definition and calculation of the three variables are given by where G u and G are all the normalized constants. ey can be acquired by In equation (33), with respect to p(x i |C u i , C v i+1 , Y k ), the mixed smoothing state estimation x uv i|k , covariance estimation P uv i|k , and the mixed interaction covariance estimation P uv i,i−1|k can be obtained by In equation (33), with respect to p(x i | C u i , Y k ), the mixed smoothing state estimation x u i|k , covariance estimation P u i|k , and the mixed interaction covariance estimation P u i,i−1|k of model u can be obtained by

Computational Intelligence and Neuroscience
In particular, if i � k, there exists Finally, after the smoothed state fusion of each model, the smooth estimation x i|k , covariance estimation P i|k , and interactive covariance estimation P i,i−1|k of the simulated degradation state can be acquired by In particular, the initial smoothed estimations are given by

M-
Step. e estimation of the model parameters θ at the jth step of the EM algorithm can be obtained by differentiating equation (31), i.e., e updated parameters are given by where ω 2 i 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]. e function is an effective method to search for the minimum value of the multidimensional function.
e specific solving process is as follows. Equation (42) is taken into equation (31) to get an expression that only contains the parameters D and η. en, 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(θ, θ

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 is w, the equipment RUL T is defined as the time when the degradation process first passes the failure threshold [47], i.e., In order to obtain the analytical expression of the PDF of RUL, the derivation of the RUL calculation at time k 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 x(t) needs to be transformed. At the current monitoring time k, the simulated degradation state is x k . According to the Markov characteristic of the Brownian motion, the Wiener process can be rewritten as en the two state equations of HWSSM are merged into one equation. In other words, the damages caused by n Poisson shocks are integrated into the Wiener process. A merged state equation is expressed as 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, x k , θ, and Y k obeys the inverse Gaussian distribution, i.e., where T k represents the RUL at time k. Its mean is given by Equation (47) does not take into account the simulated degradation state distribution obtained by the parallel simulation system. e 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(c − ξ 2 ) and A, B, and C are all constant,
e 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 process x(t), t ≥ 0
{ }, the RUL distribution at time k can be given by According to the form of the specific HWSSM and IMM filtering algorithm, x k follows the normal distribution, i.e., x k ∼N(x k|k , P k|k ). Let p(x k | Y k ) express the conditional PDF about Y k of x k . Based on the total probability theorem, the simulated degradation state distribution is integrated into the inverse Gaussian distribution and gives c � x k|k , and ξ 2 � P k|k and calculate the mathematical expectation about en, the weighted sum with the occurrence probability of the n Poisson shocks is calculated. As a result, equation (49) is obtained. is completes the proof of eorem 1.
According to the property of the mathematical expectation and equation (49), the expected value of the RUL can be obtained by Obviously, it is difficult to directly obtain the analytical expression of RUL's mathematical expectation by using the above equation.
erefore, 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., According to the equations (49) and (52), the parallel simulation system can support maintenance decisionmaking by calculating the probability density function of the RUL, and its expected value in an online and real-time manner.

Data Introduction.
e performance degradation of bearings of mechanical equipment is a typical hybrid degradation process with shock characteristic. is 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. ese data are provided by the FEMTO-ST Institute in France. e 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. e 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. e 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. e root mean square (RMS) value of the vibration signal is a commonly used degradation feature and it is calculated by where N is the sampling point number and satisfies N � 2560 and e i denotes the vibration acceleration signal at the ith sampling point. e 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. e RMS begins to change significantly after the 1500th monitoring time, which is regarded as the starting time of RUL prediction oriented parallel simulation. e 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.

Model Evolution and RUL Real-Time Prediction.
e initial simulation configurations include x 0 � 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.6 0.4] T , i.e., μ 1 0 � 0.6 and μ 2 0 � 0.4. e model transition probability matrix is chosen as P � [0.5 0.5; 0.6 0.4]. e 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 where r(r � 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. e model probability is shown in Figure 4. e Poisson shocks characteristic is not significant in the monitoring period from t 1500 to t 1765 , and the linear degradation characteristic is more obvious. e 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. e 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 t 1766 , t 1827 , t 1877 , t 2130 , t 2234 , etc. e 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.
With the dynamic injection of the observed degradation data, the parallel simulation system uses the IMM-EM algorithm to perform model evolution. e estimated results of model parameters are shown in Figure 5. 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 t 1766 , t 1827 , etc., reflecting the accelerated degradation rate of bearings 1-3. e 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, Computational Intelligence and Neuroscience 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. e 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.
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 t 1600 to t 2300 with the prediction interval of 100 monitoring points.
According to Figure 6, 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. e 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. is 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.

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. e equations of model parameters evolution and RUL prediction for the single-model method are given in Appendix.
e RUL prognostic results of the single-model method and the proposed method are shown in Figure 7. e 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 "rightbiased" 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 Figure 8, 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 singlemodel 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.
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 time k, the definition of the relative prediction accuracy of the RUL is given by where T k denotes the true RUL at time k. On the basis of RA k , the average relative accuracy MRA is defined as 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 Based on the definition of the TMSE, the smaller TMSE implies the more accurate prediction result. e calculated comparison results are shown in Table 1. 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. e TMSE of the proposed method is much smaller than that of the single-model method, implying the smaller RUL prediction error.

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. erefore, the forward IMM filtering is utilized to dynamically calculate the probabilities of different models, achieving the model soft switch. en, 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. rough 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. en, 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. e 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. is 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.
ere 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. e 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.