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.

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 Figure 1. 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 and behavior information at time . Equipment can be controlled by the control information sent by the parallel simulation system. According to the implementation of control information, it can be divided into automatic control information and manual control information .

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 Section 2. 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 aswhere is the continuous degradation process driven by the standard Brownian motion and ; 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 without considering the unknown discrete shocks is yielded aswhere is the sensor sampling interval and is the simulated degradation state at time . Time is short for time, and 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 bywhere is the damage caused by the shocks. We assume that the arrival of the shocks obeys the Poisson process [38, 39], i.e.,where denotes the total number of shocks that appear from the initial time until the moment 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 is zero or one. In addition, it is assumed that the Poisson process is independent of the Brownian motion.

Equipment degradation observation data 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.,where and denotes the variance of the measurement noise. It is also assumed that is independent of the Brownian motion. Then, equation (5) is transformed by Euler discretization, and the observation equation at discrete time points is yielded aswhere denotes the degraded observation data at time and is the noise sequence that obeys the standard normal distribution. Furthermore, is independent of .

According to equations (1)–(6), the HWSSM at the discrete time point can be formulated as

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. and represent the observation data vector and simulated degradation state vector until time , respectively. represents the number of simulation models. represents the probability of simulation model at time . represents the transition probability from the model to the model . represents the Markov probability transition matrix. denotes that the valid model is model in time interval . denotes the estimated mean of the simulated degradation state at time conditioned on and the previous measured values. denotes the estimated covariance of the simulated degradation state at time conditioned on and the previous measured values. denotes the hybrid estimated mean of the simulated degradation state at time conditioned on and the previous measured values. denotes the hybrid estimated covariance of the simulated degradation state at time conditioned on and the previous measured values. denotes the Gaussian distribution with value , mean 0, and variance .

3.1.1. Input Interaction (Model )

At this stage, the hybrid state estimation and the covariance estimation are obtained by the state estimation and model probability of model at time . and are both used as the initial state of the Kalman filtering at the time . The specific steps are as follows.

The prior probability of the simulation model is defined and calculated bywhere .

The hybrid probability transferred from the model to the model is defined and calculated by

The hybrid state estimation of simulation model is determined by

The hybrid covariance estimation of simulation model is determined by

3.1.2. Kalman Filtering (Model )

Kalman filtering can be divided into two steps, i.e., the prediction step and the updating step. and express the simulated degradation state estimation and covariance estimation, respectively, conditioned on and the previous measured values. Then, the prediction step refers to achieve the predicted results and according to the degradation equations and the results of input interaction. The prediction step can be expressed as

The updating step refers to obtain the posterior state estimation and covariance estimation of model based on the prognostic results in the prediction step and the observed data . The updating step can be expressed aswhere denotes the new information, is the variance of , is the Kalman gain, and represents one-dimensional unit matrix.

3.1.3. Model Probability Calculation

At this stage, the likelihood function of the model is used to calculate the model probability and is defined as

Furthermore, according to equation (16), can be given by

Then, the probability of the simulation model iswhere is a normalization constant that satisfies .

3.1.4. Output Interaction

According to the model probability and the estimation results of each model in the Kalman filtering stage, the degradation state estimation and covariance estimation are achieved through the weighted sum and can be given by

In addition, in the initial time of degraded state estimation, there exists and conditioned on . Once a new measurement value 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 , , , , , and , where and 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 and the th iteration of EM algorithm can be obtained bywhere denotes the mathematical expectation operator, is the joint log-likelihood function, and represents the model indicating matrix. represents the model indicating vector at time and satisfies . Its element is defined as

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 vector , the observed data vector , and the model indicating matrix at time is given as

According to HWSSM, there existswhere and . Then, 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, is denoted as the mathematical expectation of the joint log-likelihood function at time and the th iteration of EM algorithm, i.e.,

Through neglecting the irrelevant items that are independent of the parameter , can be expressed aswhere represents . 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, can be given bywhere . According to equation (31), in order to calculate , the values of the variables , , , , , and 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 Section 3.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 bywhere denotes the backward one-step prediction of the model and represents the covariance of backward one-step prediction error for the model . 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 aswhere represents the smoothed model probability, denotes the smoothed hybrid probability, and denotes the likelihood function. The definition and calculation of the three variables are given bywhere and are all the normalized constants. They can be acquired by

In equation (33), with respect to , the mixed smoothing state estimation , covariance estimation , and the mixed interaction covariance estimation can be obtained by

In equation (33), with respect to , the mixed smoothing state estimation , covariance estimation , and the mixed interaction covariance estimation of model can be obtained by

In particular, if , there exists

Finally, after the smoothed state fusion of each model, the smooth estimation , covariance estimation , and interactive covariance estimation of the simulated degradation state can be acquired by

In particular, the initial smoothed estimations are given by

3.2.2. M-Step

The estimation of the model parameters at the th step of the EM algorithm can be obtained by differentiating equation (31), i.e.,

The updated parameters are given bywhere denotes the smoothed probability of the model formulated by equation (3). Note that the updated equation of parameter contains parameter , the updated equation of parameter contains parameter , and the updated equation of parameter contains parameters and , which makes it impossible to directly use equation (42) to obtain the estimated values of parameters , , 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 and . Then, the “fminsearch” function is used to perform a two-dimensional search starting from the initial value of the parameters and . When the expression gets the minimum value, the expression obtains the maximum value. In this case, the corresponding parameter values are the estimated values of and . Finally, the estimated values of 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 the th iteration of EM algorithm. In other words, the one iteration from to 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 is , the equipment RUL 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 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 needs to be transformed. At the current monitoring time , the simulated degradation state is . According to the Markov characteristic of the Brownian motion, the Wiener process can be rewritten as

Then the two state equations of HWSSM are merged into one equation. In other words, the damages caused by 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 Poisson shocks, , , and obeys the inverse Gaussian distribution, i.e.,where represents the RUL at time . Its mean is given by and its PDF is formulated by

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 and A, B, and C are all constant,

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 process , the RUL distribution at time can be given by

According to the form of the specific HWSSM and IMM filtering algorithm, follows the normal distribution, i.e., . Let express the conditional PDF about of . Based on the total probability theorem, the simulated degradation state distribution is integrated into the inverse Gaussian distribution and gives

According to Lemma 1, let , , , , , and and calculate the mathematical expectation about of . Then, the weighted sum with the occurrence probability of the 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

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.,

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.

The root mean square (RMS) value of the vibration signal is a commonly used degradation feature and it is calculated bywhere is the sampling point number and satisfies and denotes the vibration acceleration signal at the th 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., .

5.2. Model Evolution and RUL Real-Time Prediction

The initial simulation configurations include , , , , , , and . Furthermore, the vector of initial model probability is , i.e., and . The model transition probability matrix is chosen as . 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 bywhere 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.

The model probability is shown in Figure 4. The Poisson shocks characteristic is not significant in the monitoring period from to , 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 , , , , , 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 .

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 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 , , 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 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 to 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. 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.

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 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 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.

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 , the definition of the relative prediction accuracy of the RUL is given bywhere denotes the true RUL at time . On the basis of , the average relative accuracy is defined aswhere represents the number of monitoring time points used to predict the RUL. According to the definition of , it satisfies and the bigger 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. The 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. The TMSE of the proposed method is much smaller than that of the single-model method, implying the smaller RUL prediction error.

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.

Appendix

According to the single-model method, it can be known that and . The joint log-likelihood function can be formulated as

Based on the Bayesian theorem and the Markov property and omitting independent parts, can be expanded aswhere .

Through the derivation of the M-step, the parameters evolution of the single-model method at the th iterative step is determined by

Via the iteration of Kalman filtering and EM algorithm, the model parameters can be obtained by equation (A.3). Then, the RUL distribution at time of the single-model method can be calculated by

The proof of equation (A.4) is similar to Theorem 1. The RUL expected value can be formulated as follows:

Data Availability

The bearing degradation data supporting the findings of this study are from previously reported studies and datasets, which have been cited. The processed data are available at “http://www.femto-st.fr/en/Researchdepartments/AS2M/Research-groups/PHM/IEEE-PHM-2012-Data-challenge” or in “P. Nectoux, R. Gouriveau, K. Medjaher, E. Ramasso, B. Morello, N. Zerhouni, C. Varnier, “PRONOSTIA: an experimental platform for bearings accelerated life test,” IEEE International Conference on Prognostics and Health Management, Denver, CO, USA, 2012.”

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.