{'fraction_non_alphanumeric': 0.03005434978224094, 'fraction_numerical': 0.013173523377605011, 'mean_word_length': 4.523658051689861, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A recurrent neural network approach for remaining useful life prediction utilizing a novel trend features construction method

Data-driven methods for remaining useful life (RUL) prediction normally learn features from a fixed window size of a priori of degradation, which may lead to less accurate prediction results on different datasets because of the variance of local features. This paper proposes a method for RUL prediction which depends on a trend feature representing the overall time sequence of degradation. Complete ensemble empirical mode decomposition, followed by a reconstruction procedure, is created to build the trend features. The probability distribution of sensors' measurement learned by conditional neural processes is used to evaluate the trend features. With the best trend feature, a data-driven model using long short-term memory is developed to predict the RUL. To prove the effectiveness of the proposed method, experiments on a benchmark C-MAPSS dataset are carried out and compared with other state-of-the-art methods. Comparison results show that the proposed method achieves the smallest root mean square values in prediction of all RUL.


Introduction
Contrary to the time-based maintenance strategy which checks mechanical systems on a regularly scheduled basis, conditionbased maintenance (CBM) reduces the uncertainty of maintenance which making decisions according to the real-time information through condition monitoring. Remaining useful life (RUL) prediction is a key aspect of CBM, which estimates the duration of a component or a system from the current time to the end of its lifetime. An accurate RUL prediction improves the productivity, reliability and safety of a mechanical system while minimizing its maintenance cost. With the technology maturing and costs dropping, large amounts of sensor information have been captured by industry manufacturers and research laboratories, for the purposes of fault and degradation analysis.
Given its importance, RUL prediction has been the subject of many investigations in the field of reliability for modern manufacturing [1,2]. Generally, in these papers, methods for RUL prediction can be classified into model-based and data-driven approaches. Model-based approaches develop mathematical models based on the knowledge of the failure mechanism of the models. However, for some complex systems, such as the engine turbine fan, accurate failure mechanism models are hard to build precisely due to system uncertainties and environment noises.
Data-driven approaches can be divided into statistical modelbased approaches and artificial intelligence approaches. Statistical model-based approaches [3] include methods based on the Wiener and Gamma process [4], the Kalman filter method [5][6][7], and particle filters [8,9]. Statistical models are built according to physical laws or probability distributions. Given an example, in [7], Cui et al. established models for each kind of bearing operation and incorporated unscented Kalman filter into the Bayesian estimation method to calculate the probability of each state. Instead, the artificial intelligence approaches [10][11][12] build models without knowing the state of the system. They train models based on the historical degradation information by intelligent learning algorithms.
In terms of the intelligent learning algorithms, feature extraction is first performed to capture the key information and also to reduce the dimension of the original signal. Traditional features, such as mean square value in time domain or fast Fourier transform in the frequency domain [13], are often highly noised. In addition, features captured for a mechanical system are commonly not the best representations for another mechanical system. Therefore, a lot of human effort is expended on finding proper features in a case-by-case manner [14][15][16]. Few of these features consider the uncertainty of the measurement. In fact, industrial time series https://doi.org/10.1016/j.measurement.2019.06.004 0263-2241/Ó 2019 Elsevier Ltd. All rights reserved. signals acquired by sensors are intensely disturbed by the unpredictable environment and the uncertainties of the physical systems. Therefore, each measurement point of the time series signal is uncertain itself. However, the tendency of the physical systems during system degradation is relatively certain. In this study, we regard the measurement of each point during degradation as a trend-centered probability distribution. In the Bayesian machine literature, Gaussian Processes (GPs) are popular choices for modeling a distribution over regression functions. Yet the design of the appropriate prior for GPs is hard, especially when modeling the distribution over different sensors' measurements. Conditional Neural Processes (CNPs) [17] are recently developed network systems which solve this problem. Structured as neural networks, CNPs are able to learn appropriate prior from data themselves. Thus, it is natural to use the CNPs to map the sensors' measurements in probability distribution. But as CNPs are still time consuming, the more effective trend features construction is still meaningful. In [18], M. Torres et al. proposed a complete ensemble empirical mode decomposition (CEEMD) method for analyzing transient events, which is an enhancement of the extended empirical mode decomposition (EEMD) [19] and empirical mode decomposition (EMD) [20]. CEEMD decomposes the original signal into a residue and a number of intrinsic mode functions (IMFs), representing low to fast oscillation components of the signal. In [21], Wu, J. et al. used three IMFs to track the degradation evolution of rolling bearings. In [22], M. Niu et al. use three IMFs and the residue of the CEEMD as individual features for four separate SVR models to predict the RUL. In this work, we found that the residual of the CEEMD and the extracted IMFs can also be reconstructed in an optimal way to build general a trend feature for mechanical systems.
When the trend feature of a testing machine is extracted, advanced learning algorithms are commonly followed to train prediction models. The most well-known advanced learning algorithms for RUL prediction modeling include neural networks (NN) [11] and support vector machines (SVM) [23]. SVM is trained to solve a quadratic optimization problem by introducing bound constraints and a linear equality constraint. SVM has been well investigated for RUL prediction, see successful implementations in [24][25][26]. Recently, deep learning methods (i.g., convolutional neural networks [CNN] and recurrent neural networks [RNN]) have merged into research and industry fields, and yield very accurate results in speech recognition [27] and image recognition [28] that surpass other traditional NN. The methods based on CNN and SVM learn features from a fixed time window size of a priori of degradation, which may lead to less accurate prediction results on a different dataset because of the variance of local features. To improve the prediction performance, the information of the entire degradation process should contribute to constructing a more accurate model. RNN is a method which has been found well-suited for this task due to its ability to remember past input information over an period of time. However, RNN itself suffers a drawback called the 'vanishing error problem' [29] that results in vanishing or exploding gradients when learning long-term dependencies. This drawback raises the problem of which information is valuable to remember and which to forget. [29] proposed the long shortterm memory unit (LSTM) RNN which adds memory controllers to the conventional RNN architecture to decide when to select related information for prediction and to discard redundant information over a long period without error vanishing.
In this paper, we apply a signal decomposition technique CEEMD to obtain a number of IMFs and a residue. Several trend features are reconstructed based on the IMFs and the residual. The trend features are evaluated using the probability distribution of sensors' measurement learned by CNPs. The best trend feature is chosen and used as the input feature for a model developed using LSTM (denoted as DLSTM) to predict the RUL. The learning capabilities and prediction accuracy of the proposed DLSTM RUL prediction method, with respect to other state-of-the-art approaches, are evaluated through extensive experiments using a benchmark turbine engine dataset. The proposed trend features with the DLSTM prediction method are characterized by the following main features: 1. Proposing a novel trend feature construction method by using the CEEMD method. 2. Proposing a new perspective on the measurement of the sensors, as well as conducting experiments that evaluate several trend features using the distribution of the measurement learned by the neural processes. 3. LSTM RNN, using the novel trend features as inputs, performing more precise RUL predictions when compared to the Kalman filter, SVR, and CNN methods.
The rest of the paper is organized as follows: The methodology of developing the overall RUL prediction framework is explained in Section 2. Section 3 presents the experimental studies of RUL prediction on a turbofan engine dataset and provides the evaluation and comparison of our proposed RUL method, with respect to other existing research results on the same dataset. Finally, Section 4 concludes this paper.

Methodology
The central idea of this project is illustrated in Fig. 1, which consists of a trend feature extraction process described in Section 2.1 and a LSTM modeling process described in Section 2.2. . . . ; x n ð Þ ½ 2 R n , which can be decomposed into several intrinsic mode functions (IMFs) (or modes) by empirical mode decomposition (EMD). The IMFs have to possess the following two properties: 1) Over the entire time range, the function must have the same number of local extreme points and zero crossing, or only have difference by one; and 2) The envelope of the local maximum (upper envelope) and local minimum (lower envelope) must be zero on average.
As is defined by EEMD [19], the k-th 'true' IMF (denoted IMF k ) means the average of IMFs obtained by executing EMD on an ensemble of traits generated by adding different uniformly distributed Gaussian white noises to the original dataset D. The IMF k can be calculated through the following expressions: where where w i t ð Þ 2 R n denotes the i-th realization of the white noises for i ¼ 1; Á Á Á ; I. EMD (Á) represents the original EMD function proposed in [20], which could generate k number of IMFs for i-th realization.

Complete ensemble empirical mode decomposition
CEEMD [18] works better than EEMD with lower computing cost and smaller reconstruction error. It realizes these advantages by adding Guassian noises ð Þ denotes the k-th mode of CEEMD. Now compute the g IMF k t ð Þ and the first residual r 1 : Eqs. (4)-(7) are defined as the first step of CEEMD. It can be viewed that g denotes the noise standard deviation (NSD) for k-th residual and e 0 denotes the NSD for the original signal x t ð Þ. In Eq. (5), EMD 1 Á ð Þ represents the first step of the EMD in which the first IMF mode is created.
Then, define the second step of CEEMD by the following equations to calculate k-th residual, r k , and k-th IMF, g IMF k t ð Þ: g Eqs. (8)-(10) are repeated until the residual r k t ð Þ can not be decomposed any more (i.e., r k t ð Þ does not have at least two extremes). The third step is to obtain the final residual by This means that the original signal can be reconstructed by adding up the final residual and g

Trend feature construction
As has been mentioned that the g IMF k t ð Þ is the fluctuation signal of x t ð Þ; g IMF k t ð Þ with smaller k contains more information with higher frequency oscillation but also more noise. In order to balance the containment of information and noise, the trend feature signal Q t ð Þ is finally reconstructed, where v ¼ K þ 1 À m, and m 2 1; 2; . . . ; K þ 1 f g controls the usage of the g IMF k t ð Þ.

Model training and predicting
Using the constructed trend features as the inputs, in this section, a DLSTM framework for RUL prediction is proposed. Section 2.2.1 describes how the LSTM network makes predictions over a long time without the impact of local noises and the problem of error vanishing. Section 2.2.1 and 2.2.3 explain two techniques that help to train the network. Section 2.2.4 summarizes the DLSTM algorithm.

Concept of Vanilla LSTM
Vanilla LSTM is the most widely used LSTM architecture variation. The architecture of the vanilla LSTM is shown in Fig. 2.
To summarize, the core idea of LSTM lies in two information flows and includes three gates, which function is to add and remove information. More specifically, two information flows are detailed as illustrated in Fig. 2: The lower green flow controls the combination of the short-term state of the last cell (h tÀ1 ) and the input of the current cell (Q t ). During this process, three gates named forget gate, input gate, and output gate are used to generate three signals named forget control signal, input control signal, and the output control signal (denoted f t ; i t , and o t respectively); the upper brown flow shows the update of the current cell's long-term state (C t ). Two branches are introduced into the upper flow, including the long-term state of the last cell (C tÀ1 ) controlled by f t and the current input value controlled by i t . More details of LSTM are explained below: In terms of the forget gate, the state information (f t ) is captured from the short-term state of last cell (h tÀ1 ) and the current input (Q t ), and goes through a sigmoid function r Á ð Þ: where W f ; R f are the weights, and b f is the bias parameter.
For the input gate, this gate generates the input information (i t ), which is captured from the short-term state of last cell (h tÀ1 ) and the current input (Q t ), and goes through a sigmoid function to obtain i t : where W i ; R i , and b f are the weights and bias parameters. Then, the long-term state of current cell (C t ) is updated by the forget control signal (f t ) and input control signal (i t ) obtained in Eq. (14) and Eq. (15): HereC t is the candidate memory cell gate, W c and R c are the weights for Q t and h tÀ1 , and b c is the bias value.
The bottom data flow goes through the output gate where the output control signal o t is generated by the following equation: where W o and R o are the weights for h tÀ1 and Q t , and b o is the bias value.
The output of the current cell (Z t ) and the short-term cell state of the current cell (h t ) is generated by multiplying the current long-term cell state (C t ) by the output control signal (o t ) and it is also transferred to the next cell through the bottom data flow,  [30] which aims to solve the problem of overfitting. Each neuron in the dropout layer is randomly omitted from the network with a probability P c within (0,1) while training. Networks vary through training cases, but all networks share the same weights. All neurons are activated in testing cases, where the weights are calculated by multiplying the training weight and P c .

Masking layer
Generally, the proposed model should fit and make predictions on different samples with different lengths, however, the length of the RNN is a pre-defined value. A masking layer is introduced to tackle this problem. Samples are first extended to the longest sample length by padding in front with zeros. The masking layer will work by skipping the points where all the input features are zeros.

DLSTM
According to the above descriptions, the proposed DLSTM framework (see Fig. 3) for RUL prediction can be summarized into Algorithm 1.

Data description
In this paper, NASA's turbofan engine RUL prediction problem [32] is employed to evaluate the proposed method, which was produced on the commercial modular aero-propulsion system simulation (C-MAPSS) test-bed. The detailed information of the C-MAPSS dataset are presented in the Table 1. The dataset contains four subsets, named as FD001, FD002, FD003, and FD004. Engine units for each subset for training and testing are listed in the first and second row of Table 1 respectively. FD001 and FD002 are with fault mode 1, while FD003 and FD004 are with fault mode 2.
The time-series measurement signals used for training and prediction are acquired from 21 sensors under three operational conditions associated with different sea levels. The operational status of each engine unit is healthy in the early stage and begins to degrade over time until a failure occurs. The predicted RUL indicates the time periods over which those engines are likely to operate before they require repair or replacement.
The training datasets contain the whole evolution of the engines' degradation until the end of their lifetimes, which means the RUL value at the last time point is 0. In testing datasets, only the first partial degradation is given of each engine, hence, the RUL at the last point is not 0. The objective of this study is to precisely predict the RUL during the system degradation till the last point on the testing datasets. In this project, models are trained individually based on different subsets. For each modeling training process, the training data is used as the input and the corresponding RUL is used as the label. Then the trained models are applied to the test datasets for prediction and evaluation.
Programming languages are Python 3.6 with deep learning library Keras, and Matlab 2018a. Both programs are accelerated by GPU.

Data pre-processing
C-MAPSS dataset contains time-series measurements from 21 sensors. However, measurements of some sensors hardly change over time, indicating there are no useful information in these measurements. We abandon these measurements and keep the remaining ones. Finally, the available sensor datasets are x norm n ð Þ 2 R n . The normalized measurement data x j norm t ð Þ are computed by: where x j t ð Þ is the t-th measurement of the j-th sensor. x j max and x j min denote the maximum and minimum values of measurements of the j-th sensor, respectively. Training set augmentation: Sequences in the training sets are running to the end of engines' lifetime (i.e., RULs equivalent to 0), however, sequences in the testing sets always end with nonzero RULs. To solve this problem, one step of augmentation is proposed. In this step, the last 1 to 100 points are cropped out from the end of different sequence in the training set, thus, we get some training sequences end with the RUL from 1 to 100.
Set the RUL for the training set: Considering the fact that the engine units work in its good health condition in a very long time period, the degradation in this period is negligible to estimate. In this work, we follow the suggestion in [33], in which indicates that the faults usually occur at the time step around 120$130. We set the RUL bigger than 130 to be 130.

Performance metrics
In this study two metrics are used to evaluate the performance of the proposed prognostic method, i.e. score function and root mean square error (RMSE).  The score function is defined by 2008 PHM Data Challenge Competition [34]. The score function is given as follows: s ¼ X P p¼1 s p ; s p ¼ :e À dp 13 À 1; for d p < 0 e dp 10 À 1; for d p ! 0 where d p ¼ d RUL p À RUL p ; d RUL p and RUL p denote the predicted RUL and the label of p-th engine unit, respectively. RMSE for each sub-set is obtained by the following equation: 3.4. Trend feature construction using CEEMD CEEMD is used to obtain clear trend features for measurement data of each engine unit. For CEEMD implementation, the noise standard deviation is set to 0.02, the number of realization is set to be 100, and the maximum iterations is set to 5000.
To choose the optimal selection of m, three trend features Q j 0 ; Q j 1 , and Q j 2 for j ¼ 2; 8; 12; 15 of engine unit 1 in FD001 are reconstructed by substituting m ¼ K þ 1; K; K À 1 into Eq. (13), respectively. Fig. 4 shows the evolutions of Q j 0 ; Q j 1 ; Q j 2 and D j norm of engine unit 1 in FD001 sub-set. It is observed that the Q j 0 trend feature in Fig. 4(a) with m ¼ K þ 1 offers the most clear trend (with less fluctuation) compared with x j norm in Fig. 4(d), Q j 1 in Fig. 4(b), and Q j 2 in Fig. 4(c). The trend features will be further evaluated in subSection 3.5.

Trend features evaluation
CNPs are employed to learn the probability density for the measurement of the sensors. By setting the probability distribution as the baseline for tendency features, we numerically compare Q 0 ; Q 1 ; Q 2 ; D norm with other three tendency features constructed by EEMD (i.e., E 0 ; E 1 , and E 2 ), one popular time-domain statistical features (mean). E 0 ; E 1 , and E 2 are constructed through the following expression: where R 0 t ð Þ and IMF k t ð Þ denote the residual and the intrinsic mode functions extracted based on x norm t ð Þ using EEMD method (more technique details of EEMD can refer to [19]).
mean value is computed through: where Two evaluation metrics including the mean probability density (MPD) and the mean absolute relative error (MARE) are employ to assess these features: where x c;t is the trend feature for engine unit c and time period t, l c;t is the median for the probability distribution. r c;t is the standard derivative for the probability distribution. C is the number of engine units, T c is the number of time periods for engine unit c.   demonstrating that Q 0 is the best choice of representing the degradation tendency of this mechanical system. The temporal evolution of tendency feature Q 0 , normalized original signal D norm , 95% probability interval (PI) distribution, and the median curve of the distribution from Sensor 2, 8, 12, and 15 are also shown in Fig. 5. It is observed that Q 0 is almost overlapped with the median curve of the distribution. We can conclude that the measurements along the Q 0 are very close to the highest probability of measurement points for system degradation.

DLSTM modeling and optimization
To optimize the network architecture, we train and test several models with different architectures on subset FD001. The most suitable architecture for this task is comprised of one masking layer, two LSTM layers, one dropout layer, and one fullyconnected layer. Moreover, RUL prediction performance is also varied by changing the output number of the LSTM layer. Table 4 presents the score and RMSE of several configurations. It can be seen that the configuration with the first LSTM layer of 128 outputs and the second LSTM layer with 100 outputs show better prediction performance than other configurations. Therefore, experiments in this project use this configuration.   Fig. 8 shows RUL prediction results of engine unit numbers 5, 33, 43, and 77 which have been chosen at random from 100 engine units in the FD001 subset. In Fig. 8, during the healthy condition of the engine unit (flatten stage of the label), the proposed method shows precise RUL prediction where the red lines are almost overlapping with the blue lines. When degradation begins, the prediction translates smoothly and then decreases almost linearly with time until the end of the period. All of the four engine units show precise and accurate RUL predictions when close to the end of the period. Near the point of degradation, it appears that the evolution of the predictions does not match the given labels. This is due to the fact that not all of the engine units begin to degrade when RUL equals to 130 and not all of them degrade suddenly without transition. Meanwhile, for the implementation of the proposed method in industry, the accuracy of the prediction is of more importance when the engines work well and when they go to fault. The accurate predictions can reduce the maintenance cost as well as guarantee the safety of the system.

Prognostic results
To further investigate the prediction of all the engine units in one subset, we summarize the last RUL prediction values of all the 100 engine units in FD001 in Fig. 9, where the last point values are sorted by ascending number. It can be seen that the prediction values is closer to the label values when the engine units are close to the end of their lives (i.e., RUL < 50). This result indicates that the less the RUL there is, the higher of the prediction accuracy    achieved. This is because when the engine unit is about to break down, the trend features of the degradation are enhanced and are captured more easily by the DLSTM.

Comparison
The performances of the different methods evaluated on the FD001 test set are presented in Table 5. The proposed method (DLSTM) is compared with the following methods: 1. RNN with Gate Recurrent Unit (GRU): Besides the difference that the GRU does not use a memory unit to control the flow of the information like the Vanilla LSTM, the RNN with GRU has the same architecture as the DLSTM in Fig. 3 In Fig. 10 we show the comparison of the RUL predictions of DLSTM and CNN on the 58th engine unit in the FD001 test set. Since the CNN method is based on a time window with the length of 25 [35], the prediction of CNN starts at the time of 24. The DLSTM prediction has a more continuous prediction (a smoother curve) than CNN. Meanwhile, the DLSTM has a better prediction when the engine unit is still in healthy condition. These performances of DLSTM benefit from considering of the overall trend as a learning feature. With the overall trend, the prediction of the RUL will not be influenced by the noises on local features. Table 6 compares the method proposed in this paper with several existing methods including switching Kalman filter ensemble [36], SVR [15], deep CNN [38], and vanilla LSTM [37]. From Table 6, last row, we see that the proposed method shows better performance with its RMSE and Score metrics are both having smaller values for all the four subsets in C-MAPSS dataset. More specifically, compared with SVR and CNN where predictions are based on the local features, the proposed method considers the overall trend feature of the sequence and results in a smoother degradation curve; Compared with the model-based method such as the Kalman filter, the proposed method offers more flexibility in adjusting the architecture and has better performance in learning the hidden relationship between features and the prediction; Furthermore, the trend feature extracted by CEEMD in this project helps to prevent the problem of the overfitting of vanilla LSTM which uses noised feature to train the model.

Conclusion
In this paper, we proposed a novel RUL prediction method by using complete ensemble empirical mode decomposition (CE-EMD) for trend feature construction and long short-term memory recurrent neural network (LSTM RNN) for prediction, which are able to prevent the impact of signal fluctuation as well as minimize the prediction error. The trend feature construction associated with mechanical system degradation, which is the major challenge in this study, is constructed according to the combination of intrinsic mode functions (IMFs) and the signal residue. The proposed RUL prediction method is examined on a benchmark dataset of aircraft turbofan engines, and compared to the Kalman filter, support vector regression, deep CNN, and vanilla LSTM. Results verified that the proposed method achieved smaller prediction errors.
As a future project, the adoption of the proposed method in more complex mechanical systems will be considered. Furthermore, we aim to find what hidden information can be extracted from different oscillations to further improve RUL prediction accuracy.