Remaining Useful Life Prediction for Complex Systems With Multiple Indicators Based on Particle Filter and Parameter Correlation

In practical applications, the failure of large-scale complex equipment is often caused by the simultaneous degradation of multiple components. It is necessary to predict the remaining useful life (RUL) of the equipment with multiple degradation indicators. This article proposes a new joint-RUL-prediction method in the presence of multiple degradation indicators based on parameter correlation. The stochastic process model is established for each degradation indicator, and the model parameters are estimated by kernel smoothing particle filter (KS-PF) and maximum likelihood estimation (MLE). Meanwhile, to facilitate the dependencies between multiple degradation indicators, correlations of the degradation model parameter between multiple degradation indicators are established in KS-PF. In addition, optimal tuning (OT) is introduced to choose the best kernel parameter. A case study on the Commercial Modular Aero-Propulsion System Simulation (C-MAPSS) dataset is applied to verify the proposed method, the experiment shows that the proposed joint-RUL-prediction method based on parameter correlation possesses a superior prediction performance compared with that by using a single degradation indicator.


I. INTRODUCTION
As an essential part of prognostics and health management (PHM), remaining useful life (RUL) prediction can effectively help to improve the reliability of equipment, avoid the huge loss caused by equipment failure, and reduce the cost of equipment maintenance. The failure of large-scale complex equipment is often caused by the simultaneous degradation of its multiple components [1], [2]. The multiple degradations of different components are identified as multiple degradation indicators, and the RUL prediction based on multiple degradation indicators is expected to capture the health status of equipment effectively.
The RUL prediction method can be classified into physics model-based methods, empirical model-based methods, artificial intelligent (AI) methods, and hybrid methods [3]. The empirical model-based methods are the most popular one The associate editor coordinating the review of this manuscript and approving it for publication was Yu Liu . among the four categories and have been widely used in recent years. When the empirical model-based methods are used for RUL prediction, an empirical model describing the degradation process of equipment is established first, and then the model parameters and degradation states are updated using the degradation measurements, consequently deriving the RUL. Stochastic process-based models such as Wiener process, Gamma process, and Markov chain-based method are commonly used [4], [5] as the degradation models and they are more interpretable compared with other models [6]- [9]. The Bayesian filter methods such as Kalman filter (KF) and particle filter (PF) are widely used for model updating in recent years [10]- [17]. However, most of the existing empirical model-based methods only involve a single degradation indicator, which is insufficient for characterizing the complex degradation process of equipment in most cases. Consequently, the empirical model-based RUL prediction involving multiple degradation indicators should be developed necessarily for practical applications. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ For equipment with multiple performance characteristics, or the systems which have multiple components, the multiple performance characteristics or components usually degrade simultaneously when the equipment or system degenerate, and the multiple indicators can be identified to indicate the degradation degree of different performance characteristics or components. And there may be dependencies between these indicators since they contain the information of the same degenerate equipment or system. Here the failure definition of the series systems is considered in this article, i.e., for equipment with multiple degradation indicators, the first hitting time when any degradation indicator exceeds the corresponding failure threshold is regarded as the Endof-Life (EoL) [1], [18], [19]. Based on this idea, this article intends to use the stochastic process model and PF to realize the RUL prediction of equipment with multiple degradation indicators. A multivariate degradation model needs to be established first, and the challenging part is how to characterize the dependencies between the multiple degradation indicators reasonably.
The dependencies are established on the stochastic term of the degradation process in some literature [1], [4], [18]- [20]. Peng et al. [1] established multivariate degradation models with copula functions and multiple stochastic processes, where the copula functions were used to establish the dependencies on the increments of degradation states. Rodriguez-Picon et al. [20] established different multivariate degradation models considering copula functions and different combinations of stochastic processes as marginal distributions, the dependencies were established on the degradation states, and then the product reliability under heterogeneous models was studied on the crack propagation data. Fang et al. [18] proposed a bivariate stochastic process model, with the dependencies establishing on the degradation states, and the Bayesian method was used for parameter estimation. Xi et al. [19] characterize the dependencies among different degradation indicators using a diffusion coefficient matrix, and a multivariate degradation model for multiple indicators was established. Then the Monte Carlo simulation method was used to predict the RUL. In addition to the above methods in which the dependencies are established on the multivariate degradation models for multiple indicators, Sun et al. [21] directly established the dependencies on multiple RUL probability distributions estimated in terms of different degradation indicators.
Since the model parameters can affect the long-term degradation trend of indicators, we consider establishing the dependencies on the model parameters in this article, and the long-term dependencies among different indicators are expected to be characterized accordingly. With the establishment of correlations of model parameters, the degradation states of indicators which are estimated using model parameters can capture the long term dependencies among different indicators effectively. Therefore, a joint-RUL-prediction method based on the stochastic process model and PF is proposed in this article, characterizing the dependencies between multiple degradation indicators as the correlations between parameters of the models built for the corresponding indicators. The long-term dependencies among different indicators are taken into consideration using the correlation between parameters, and the multivariate degradation modeling in the proposed method can consequently provide accurate information of dependencies for RUL prediction.
In the proposed method, parameters in the degradation model can be divided into 1) the parameters of degradation trajectory and 2) the fluctuation parameters, according to their physical implications. Maximum likelihood estimation (MLE) is used on the training units to estimate the fluctuation parameters offline. Then kernel smoothing particle filter (KS-PF) [22], [23] is used to estimate the degradation trajectory parameters as well as the degradation states simultaneously online. In addition, optimal tuning (OT) [24] is adopted in KS-PF, i.e., OTKS-PF to obtain the optimal kernel parameters at each monitoring time when the degradation state and the parameters are estimated. To characterize the dependencies between multiple degradation indicators, the correlations between the model parameters are introduced to generate the degradation trajectory parameters in OTKS-PF. Then the joint probability distribution of the multidimensional degradation state of indicators is estimated, and the RUL can be consequently predicted with the prespecified failure thresholds. The respective RUL predicted based on a single degradation indicator and multiple degradation indicators are comprehensively compared by using the Commercial Modular Aero-Propulsion System Simulation (C-MAPSS) dataset for the illustration of the effectiveness and superiority of the proposed method.
The remainder of this article is organized as follows. Section II describes the establishment of the marginal degradation models for each degradation indicator and the estimation method for fluctuation parameters. OTKS-PF is introduced in section III, to simultaneously estimate the degradation states and model parameters on a single degradation indicator. The joint-RUL-prediction method based on parameter correlation with multiple degradation indicators is proposed in section IV. Section V performs a case study on the C-MAPSS dataset, and the RUL predicted based on a single degradation indicator and multiple degradation indicators respectively are compared. Section VI gives a conclusion.
The measurement variability represents the uncertainty caused by measurement error. Therefore, the following stochastic process models are established for each marginal degradation process: where k = 1, 2, . . . , K represents the kth degradation indicator. X k (t) is the degradation state of the kth degradation indicator at time t. µ k (θ k , t) is the mean function of the stochastic process X k (t), and σ k B (t) is the stochastic term of X k (t). µ k (θ k , t) is a differentiable and monotonical function which is defined as the degradation trajectory, and it may be nonlinear. θ k = [a k , b k , · · · ] T is a vector of the degradation trajectory parameters, which determines the degradation rate of different units. Therefore, µ k (θ k , t) describes the nonlinear variability of the degradation process, and θ k characterizes the unit-to-unit variability by specifying the degradation trajectories of different equipment. B (t) is a standard Brownian motion process and σ k is the volatility parameter.
, and it describes the temporal variability by representing the inherent uncertainty of degradation process (1). Y k (t) is the measurement of the kth degradation indicator. ε k is a zero-mean Gaussian random variable with the variance ν 2 k , which can be seen as the measurement noise representing the measurement variability. Then the marginal degradation models are established for each degradation indicator and the models can characterize all the categories of uncertainty that may exist in equipment degradation processes.
It is considered that σ 2 k and ν 2 k are only related to the operating conditions and measuring conditions of equipment. In other words, σ 2 k and ν 2 k are the same and θ k is different for different equipment which runs under a single operating condition. Therefore, σ 2 k and ν 2 k are defined as fluctuation parameters which are different from θ k , and they are estimated using different methods in this article.
Let State transition equation f x k,j | x k,j−1 : Measurement equation h y k,j | x k,j : When (3) and (4) are used on multiple degradation indicators, let x j = x 1,j , · · · , x K ,j T denote the states of all degradation indicators at time t j , and y j = y 1,j , · · · , y K ,j T denote the measurements. Then the state transition equation and the measurement equation are converted to f x j | x j−1 and h y j | x j respectively. For equipment with K degradation indicators, if the kth degradation indicator X k (t) exceeds its failure threshold H k first, and that time is regarded as the failure time. The estimated RUL at time t j is defined as

B. ESTIMATION OF FLUCTUATION PARAMETERS
The parameters θ k , σ 2 k , and ν 2 k in (3) and (4) for each degradation indicator need to be estimated. For parameter estimation, the fluctuation parameters σ 2 k and ν 2 k are estimated offline using MLE based on the training units. Then PF is used to estimate θ k and the degradation state X k (t) online simultaneously for each testing unit. The RUL prediction result based on the kth degradation indicator can be obtained after the twostep parameter estimation. The estimation of the fluctuation parameters is described in detail as follows.
Here the superscript ' * ' is used to represent the variables for the training units. Let represent the states, measurements, and degradation trajectories of the kth degradation indicator respectively for training units. Then the distributions of X * k and Y * k can be derived by using the properties of standard Brownian motion: where Q = min t i , t j 1≤i,j≤M and I M is an identity matrix of order M . The log-likelihood functions of the fluctuation parameters are expressed as follows: and The measurements in Y * k of each training unit are collected, while the degradation states in X * k are usually implicit. In order to estimate the fluctuation parameters, the degradation state X * k can be distinguished through some appropriate methods first, such as the simple moving average (SMA), the VOLUME 8, 2020 weighted moving average (WMA), exponential moving average (EMA), etc, and then ν 2 k can be estimated by maximizing the log-likelihood function (9). After a proper functional expression µ k (·) is used to characterize the degradation trajectory, the least-squares method is used to fit the degradation state X * k of training units to obtain the estimated θ * k as well as the degradation trajectory. σ 2 k can be finally estimated by maximizing (8).

III. SIMULTANEOUS ESTIMATION OF STATE AND PARAMETER WITH OTKS
To present the proposed method clearly, the RUL prediction with a single degradation indicator is discussed in this section first, and OTKS-PF is introduced. When the kth degradation indicator is used for RUL prediction, after the fluctuation parameters σ 2 k and ν 2 k are estimated, θ k and X k (t) are jointly updated according to the new measurements using KS-PF. Moreover, the kernel parameter that is used in KS-PF is optimized online by introducing OT. Then the RUL can be predicted based on the estimated parameters and states.

A. PARTICLE FILTER
PF is a Bayesian method for state estimation. Let x j be the state of equipment at time t j that needs to be estimated in PF, and y j be the measurement at time t j . After establishing the state transition equation and the measurement equation for equipment, the posterior probability p x j | y 1:j can be obtained through two steps: prediction and update. For prediction, and for update, where the normalizing constant is Since the integral needs to be calculated in (10) and (12), it is difficult to implement the above optimal Bayesian filtering algorithms for nonlinear and non-gaussian systems. To solve the problem, the Monte Carlo sampling is introduced into PF to represent the probability distribution by particles that sampled from it, so the posterior probability of x j can be estimated as a set of weighted particles (samples) and N is the number of particles. x i j is the ith particle of x j and w i j is the weight associated with it. The posterior probability distribution of x j is estimated as is the Dirac delta function.

B. KERNEL SMOOTHING
The degradation state x k,j and the degradation trajectory parameter θ k,j can be jointly updated according to the measurement y k,j , where θ k,j = a k,j , b k,j , · · · T are the estimated parameters of θ k at time t j . Here the most commonly used approach is the augment PF, in which the parameter θ k,j and the system state x k,j are all considered as the hidden states which are estimated in PF, and the artificial evolution is introduced to realize the evolution of parameters in the particles. However, this approach will increase the variance of particles, and the convergence of PF will be delayed [22]. In this article, kernel smoothing (KS) is introduced into PF to solve the problem. The evolution of parameters in the particles θ i k,j | θ i k,j−1 is achieved by two steps: shrinkage and perturbation. For the shrinkage, and for the perturbation, where is the particle set of θ k with the particle weights at time t j−1 .
is a diagonal matrix in which the entries outside the main diagonal are all zero. s ∈ [0, 1] is the kernel parameter.
After the shrinkage, each particle of parameters shrinks towards the mean of the parameters' particles at the previous moment. The larger the kernel parameter s is, the greater the volume of the shrinkage is. When s = 1, each particle of parameters will shrink to the mean and all the particles will be the same, so the variance of the particles will decrease due to the shrinkage. Then the artificial random walk is added to ensure the diversity of particles in the perturbation step. Finally, the variance of particles keeps unchanged and the diversity of particles is also kept after the shrinkage and perturbation. Let q = √ 1 − s 2 , and the state transition equation (3) then can be deduced as Using the sequential importance resampling (SIR) framework with the state transition equation (15) and the measurement equation (4), X k (t) and θ k are jointly estimated online by the following steps: 1. Initialization: Let j = 0, and the initial distributions are set as p(x k,0 ) and p(θ k,0 ).

State estimation:
Let j = j + 1. If the measurement y k,j exists, perform the following steps, otherwise terminate this procedure.

1) Prediction: the particles {x
are sampled from the proposal distributions p(x k,j |x k,j−1 ) and p(θ k,j |θ k,j−1 ), and the prior probability of x k,j and θ k,j can be represented by the weighted particles.
2) Update: the weights of the particles are updated according to the measurement y k,j at time t j : and the weights are normalized aŝ 3) Resampling: the number of effective particles is approxi- If N neff is less than a given threshold, the resample method is implemented, otherwise resampling is not performed. After the particle resampling, the updated particles are obtained as , and the posterior probability distribution of x k,j and θ k,j can be consequently represented by the weighted particles. 4) Iteratively perform Steps 2-3.

C. OPTIMAL TUNING
The value of the kernel parameter s affects the scale of shrinkage and perturbation of the particles, so it is important to select an appropriate value for s. s is usually set as a constant in the present literature [1], [22]. Instead, we adopt OT here, to select the optimal kernel parameter s j at each monitoring time t j . The estimation accuracy of PF is related to the proposal distribution. The optimal proposal distribution should minimize the variance of the updated weights of particles, and it can also be interpreted that a good proposal distribution should make the prior distribution p x j |y 1:j−1 as similar as possible to the posterior distribution p x j |y 1:j [24], [25]. KL divergence is an index to measure the matching degree of two probability distributions, so the optimal s j is chosen in OT to make the KL divergence between the prior distribution and the posterior distribution minimal [24].
According to PF, the prior distribution of the state of the kth degradation indicator at time t j is estimated as and the posterior distribution is The KL divergence between the prior distribution and the posterior distribution is We can directly derive to (21) by substituting (11), (12), and (18) into (20) using the sifting property of the Dirac delta function.
Substitute (16) into (17), and the following equation can be derived: therefore, substitute (22) into as (21) as Finally, the optimal kernel parameter for the kth degradation indicator at monitoring time t j is derived as

IV. RUL PREDICTION WITH MULTIPLE DEGRADATION INDICATORS BASED ON PARAMETER CORRELATION
To characterize the dependencies between multiple degradation indicators, we establish the correlations between parameters of the state transition equation, and the dependencies are also taken into consideration in OTKS-PF for state estimation of multiple indicators. The commonly used RUL prediction methods with multiple degradation indicators establish the dependencies on the stochastic term of the degradation process, i.e., it is considered that there are correlations between σ 1 B 1 (t) , · · · , σ K B K (t) in (1) which cause the dependencies between X 1 (t), · · · , X K (t). Here we consider that the hidden dependencies between different degradation indicators are caused by the correlations between θ 1 , · · · , θ K , and the RUL prediction method with multiple degradation indicators based on parameter correlation is proposed.
The state estimation based on any single degradation indicator using OTKS-PF has been introduced in section III, and then the RUL prediction can be conducted based on the estimated states and parameters. Consider the RUL prediction method with multiple degradation indicators, the prediction framework is still the same as that with a single indicator. After the establishment of marginal degradation models and the estimation of fluctuation parameters, the proper multivariate state transition equation and the measurement equation can be established according to the marginal degradation models. Then OTKS-PF can be expanded considering the correlations between θ 1 , · · · , θ K for the case of multiple indicators, specifically, the correlations are introduced in the generating process of parameters' particles in the kernel smoothing step. The parameters θ 1 , · · · , θ K and the system states X 1 (t), · · · , X K (t) are jointly estimated online using OTKS-PF, and the joint probability distribution p x 1j , · · · , x Kj of the degradation state x j is obtained. At the prediction stage, the last obtained particles (the joint probability distribution) are propagated using the state transition equation f x j | x j−1 until equipment fails, i.e., any propagated indicator exceeds its corresponding failure threshold. Finally, the probability distribution of RUL can be obtained and the joint RUL prediction with multiple degradation indicators based on OTKS-PF is consequently derived.

A. JOINT STATE ESTIMATION OF MULTIPLE DEGRADATION INDICATORS BASED ON PARAMETER CORRELATION
It is considered that the θ 1 , θ 2 , · · · , θ K in (1) are linearly correlated with each other, and these correlations cause the dependencies between different degradation indicators when equipment degrades. And it should be noted that the linear or nonlinear correlations between multiple degradation indicators X 1 (t), X 2 (t), · · · , X K (t) can be characterized using a combination of the linearly correlated parameters and the mapping relationship between θ k and X k (t) defined in (1). Therefore, compared with the RUL prediction method with a single indicator, the proposed joint estimation method needs to estimate ont only the parameters θ k , σ 2 k , and ν 2 k , but also needs to estimate the Pearson correlation coefficients between θ 1 , θ 2 , · · · , θ K . Assume that the models for all degradation indicators are set as µ k (θ k , t) = a k exp (b k · t), where θ k = [a k , b k ] T . The Pearson correlation coefficients between a k 1 and a k 2 , i.e., ρ a k 1 ,k 2 and that between b k 1 and b k 2 , i.e., ρ b k 1 ,k 2 can be estimated on the training units respectively, To characterize the correlations between parameters, this approach generates the parameters' particles of different indicators from a multivariate Gaussian distribution directly to realize the evolution of parameters. This multivariate Gaussian distribution is the proposal distribution of the parameters with the correlations between each component. Let 1,j , · · · , b i K ,j ] T , and x i j = [x i 1,j , · · · , x i K ,j ] T be the ith particle of parameter a, parameter b and the state of multiple degradation indicators x at time t j , respectively, where a = [a 1 , · · · , a K ] T , b = [b 1 , · · · , b K ] T , and x = [X 1 (t), · · · , X K (t)] T . The state transition equation f x j | x j−1 of the proposed joint estimation method with multiple indicators is where It should be noted that there is no correlation between a i j and b i j in (25). It can be seen that when particles a i j and b i j are generated from the proposal distributions in the KS step, the proposal distributions are two multivariate Gaussian distributions with non-zero covariance, i.e., ρ a k 1 ,k 2 and ρ b k 1 ,k 2 which represent the correlations between indicators are contained in the covariance matrices of the multivariate Gaussian distributions. In this way, there are correlations of the parameters' particles between different indicators that are sampled from the proposal distribution, and then the weights of particles can be updated based on the newly collected measurements.
When the state transition equation (25) and the measurement equation (4) are used to estimate the degradation states and the parameters simultaneously with OTKS-PF, the steps are still the same as that in the section III-B, while the difference is that the update equation (16) is derived as and (24) is derived as so as to determine the optimal kernel parameter s j in the OT step.

The weighted particles {x
can be obtained at time t j using this joint estimation method with multiple degradation indicators based on parameter correlation, and the obtained particle x i j = [x i 1,j , · · · , x i K ,j ] T representing the estimated states of multiple indicators is multidimensional with each dimension represents the state of the corresponding degradation indicator. There are correlations between the obtained x i k 1 ,j and x i k 2 ,j , a i k 1 ,j and a i k 2 ,j , b i k 1 ,j and b i k 2 ,j for any k 1 , k 2 ∈ [1, K ]. Therefore, the dependencies between different degradation indicators are established on the corresponding model parameters using the proposed method, and the joint probability distribution of degradation indicators p x 1j , · · · , x Kj is obtained.

B. RUL PREDICTION WITH MULTIPLE DEGRADATION INDICATORS
According to the joint estimation method with multiple indicators using OTKS-PF in Section IV-A, the particle set can be obtained at time t j . Then the joint probability distribution of the state of multiple degradation indicators is The RUL corresponding to the ith particle which is predicted at monitoring time t j is and the probability density function (PDF) of RUL that is predicted at time t j is

V. EXPERIMENTAL VERIFICATION
The turbine engine data (C-MAPSS) that is published by NASA is applied to demonstrate the effectiveness of the proposed method for series systems, which contains four subsets (named FD001 ∼ FD004). Details of these subsets can be accessed in [26], and the subset involving a single failure mode and a single operating condition (FD001) is adopted here. The proposed joint-RUL-prediction method with multiple degradation indicators (OTKS-PF-Joint) is used to predict the RUL of the engines, and the prediction result is compared with the RUL result that obtained with a single degradation indicator (OTKS-PF-Single).

A. DATASET DESCRIPTION AND PREPROCESSING
The dataset contains degradation data from 100 simulated engines, each of them runs from a specific state to system failure under a single operating condition and single fault mode. The degradation status of each engine is measured by 21 sensors. The degradation data are preprocessed first, which includes standardization and dimension reduction. Z-score standardization is carried out on the original data for all the 21 sensors, and PCA is used for dimension reduction to remove features that provide little degradation information. The obtained two features using PCA, i.e., the first two principal components, retain 87.6% of the variance of the original data, so the first two principal components are chosen as the two degradation indicators which characterize the degradation of engines. Then the proposed joint-RUL-prediction method with two degradation indicators is conducted to predict the RUL of engines, and the RUL prediction result based on a single indicator is compared with that based on two indicators.
The stochastic process models are established as (1) and (2) for each degradation indicator, where k = 1, 2 represents the kth degradation indicator. Then SMA is used to obtain the smoother curves regarded as the states of the degradation indicators X * k (t) , k = 1, 2. Figure. 1 shows the measurements and the degradation states of the engine #1.
To demonstrate the effectiveness of the proposed method for series systems, the failure thresholds of indicator #1 and indicator #2 are set as H 1 = 1.22 and H 2 = 0.85, respectively. And the failure time of the engine is defined as the first hitting time that any indicator exceeds the corresponding VOLUME 8, 2020 failure threshold, which is consistent with the failure definition of the series systems. For all the 100 engines, 5-fold cross-validation is adopted to divide the testing units and the training units and five experiments are conducted using different training units and testing units. The run-to-failure data of training units are already collected, while partial data of each testing unit are removed to simulate the practical cases, where there are only partial data for the practicer to derive RUL. In each experiment, the RUL of the testing units is predicted at different monitoring time. The following analysis takes the first experiment as an example due to the limit of the length of this article, and the results of all the five experiments are also illustrated.

B. DEGRADATION MODELING AND PARAMETER ESTIMATION
The marginal degradation models are established according to (1) and (2), and the degradation trajectory µ k (·) should be determined first. Liu et al. [27] predicted the RUL on this dataset with the exponential model, here we also establish the exponential model for degradation indicators as Two strategies are considered when the model (31) is used for testing units, and the commonly used strategy (strategy I) is that a k , b k and c k are all considered as unknown parameters for testing units which should be estimated by PF. However, the estimation error of parameters will be large because of the interaction between b k and c k in PF [22], and considering the consistency of different units given the same operating condition and failure mode [3], the degradation processes of different units are expected to be similar. Therefore, strategy II can be chosen by defining b k or c k as constants which are calculated using training units for the testing units.
Degradation data of training units are used to verify strategy I first. Consistent with section V-A, the measurements Y * k (t), k = 1, 2 of the training units are processed by SMA, and the degradation states are obtained as X * k (t) , k = 1, 2. Matlab curve fitting toolbox is used to fit the degradation states of the two indicators with the model (31) to get the fitting parameters a * k , b * k , and c * k using the least-squares method, and Fig. 2 is the boxplot that illustrates the distribution of the fitting parameters.
It can be seen from Fig. 2 that, different from the parameters a * k and c * k , there are a number of outliers in the distributions of b * 1 and b * 2 , which indicates b * k may not reflect the consistency of different units under the same failure mode and operating condition accurately. Moreover, according to [28], model parameters can be categorized as fixedeffect and random effect parameters, and the former describes the homogeneity among the units. Therefore, strategy II is considered by defining the parameters b 1 and b 2 as two mixed-effect parameters which are estimated as the median of the fitting parameters: b 1 = 0.085, b 2 = 0.056, and a * k and c * k are random for testing units to describe the heterogeneity among the units. For this experiment (the first experiment),   model (31) can be redefined by using strategy II: Figure. 3 shows the fitting result of two engines in the training set using strategy II, and it can be seen that the model fits the degradation data well using strategy II.
As a result, strategy II is used for the testing units in this experiment. Then (8) and (9) are used to estimate the fluctuation parameters σ 2 k and ν 2 k on the training units according to the method in section II-B. Table 1 shows the parameters estimated based on the training units.

C. JOINT RUL PREDICTION BASED ON PARAMETER CORRELATION WITH TWO INDICATORS
After the establishment of the degradation models and the estimation of the fluctuation parameters, OTKS-PF is used for joint state estimation of the two degradation indicators involving parameter correlation. According to (25), the following state transition equation f x j , c j | x j−1 , c j−1 is established for the two-dimensional degradation indicators of the engine: where k = 1, 2 represents the kth degradation indicator.
wherec * k = 80 u=1 c * k,u /80, k = 1, 2, and it is calculated as ρ c 1,2 = 0.9873 according to (34). Figure. 4 shows the scatter plot of the parameters c * 1 and c * 2 of all the training units. It can be seen that there is a significant correlation between the parameters c 1 and c 2 from Fig. 4 and the calculated value of ρ c 1,2 , which indicates that the dependence between degradation indicators is evident and can be characterized in the correlation of parameters.
The OTKS-PF method in Section IV is used to conduct the joint estimation of X 1 (t), X 2 (t), c 1 , and c 2 online using the state transition equation (33) and the measurement equation (4). The initial distribution of states and parameters p(x 0 ) and p(c 0 ) are set as Gaussian distributions, i.e.,   optimal fitted values (POFV). After all the measurements are collected, the POFVs are obtained by fitting the degradation states estimated by using SMA. It can be seen that the estimation result of OTKS-PF-Joint is similar to that of OTKS-PF-Single. The estimated parameters can gradually converge to the POFV with the increase of the number of collected measurements. Figure. 6 and Fig. 7 respectively show the RUL prediction results of engine #19 and engine #68 in the testing set. Figure. 6-a and Fig. 7-a show the predicted PDFs and RULs estimated using indicator #1, indicator #2, and OTKS-PF-Joint respectively. It can be seen from Fig. 6-a that the predicted PDF of RUL with OTKS-PF-Joint keeps closer tracking to the predicted PDF based on indicator #2 than that based on indicator #1. Such a situation changes and the predicted PDF of RUL with OTKS-PF-Joint keeps closer to the indicator #1 in Fig. 7-a. These are because the degradation indicator #2 of engine #19 exceeds its corresponding failure threshold first, while the degradation indicator #1 of engine #68 exceeds the failure threshold first. The prediction results in Fig. 6-a and Fig. 7-a demonstrate that no matter which indicator exceeds its corresponding failure threshold first, the PDF of RUL predicted by the proposed OTKS-PF-Joint can track the PDF that predicted based on that indicator in an adaptive way. As a result, the proposed method will provide accurate RUL prediction for equipment according to the definition of equipment failure, i.e., the first hitting time that any indicator exceeds the corresponding failure threshold is regarded as the EoL.
Subsequently, the above prediction results of the proposed method which are shown in Fig. 6-a and Fig. 7-a are analyzed in detail. At the prediction stage, each particle is propagated using the state transition equation, and the predicted failure time of the equipment based on each particle is defined as the first hitting time that any dimension of the particle reaches the corresponding failure threshold. Therefore, the predicted PDF of RUL with multiple degradation indicators is obtained and the predicted RUL value will be accurate. Figure. 6-b and Fig. 7-b respectively show the RUL prediction result of engine #19 and engine #68 using OTKS-PF-Joint. The widths of the 95% confidence intervals become narrower as the monitored equipment running, which means that the uncertainty of prediction decrease gradually with the number of collected measurements increasing. It also can be seen that the predicted RUL is getting closer to the real one as the monitored equipment running, which indicates that the accuracy of RUL prediction increases when the monitored equipment runs closer to its EoL.
To compare the RUL prediction results of the proposed OTKS-PF-Joint and the OTKS-PF-Single intuitively, the RUL of testing units is predicted at each cycle after the engine runs beyond two-thirds of its whole life so that several illustrative prediction results can be obtained for each testing unit. Then the mean of the several prediction results is calculated. Figure. 8 shows the mean prediction error and the mean width of the 95% confidence interval for the testing units.
It can be seen that the prediction errors of the proposed method are generally lower than that of the RUL prediction via OTKS-PF-Single. The 95% confidence interval widths  obtained by the proposed method are also narrower than that with a single indicator. This result confirms the validity of the proposed joint-RUL-prediction method with multiple indicators based on parameter correlation. This result confirms the validity of the proposed joint-RUL-prediction method with multiple indicators based on parameter correlation.
For each experiment, the mean prediction errors and the mean widths of 95% confidence interval for all testing units are averaged, and Table 2 gives the average prediction result using the proposed OTKS-PF-Joint and the OTKS-PF-Single respectively for all 5 experiments of the 5-fold cross-validation. The results demonstrate that the proposed method is superior to OTKS-PF-Single in the accuracy and the uncertainty of the prediction.

VI. CONCLUSION AND DISCUSSION
In this article, a joint-RUL-prediction method with multiple degradation indicators based on parameter correlation is proposed by combining the stochastic process model and PF. The dependencies between different degradation indicators are characterized as the correlations between their parameters, and the covariance matrices of model parameters between different degradation indicators are introduced into KS-PF to quantify such a correlation. Besides, the OT approach is used to choose the optimal kernel parameter in KS. Then the joint probability distribution of multidimensional degradation state can be estimated via OTKS-PF considering parameter correlation, and the RUL of equipment is consequently predicted with multiple degradation indicators. A case study regarding engine degradation is performed, and the result shows that the proposed method is superior to the prediction method with a single indicator in prediction accuracy and prediction uncertainty.
The proposed method is designed for the product with multiple performance characteristics, or the system with multiple components, in which the dependencies between multiple degradation indicators can be characterized as the correlations of model parameters. In addition to the series systems, the proposed method can also be applied to other systems with differnt failure structures, such as the parallel systems or the combination systems, since the degradation states of multiple degradation indicators can be estimated simultaneously by the proposed method.