Brigitte Chebel-Morello

Simon Malinowski
email: simon.malinowski@femto-st.fr

H Senoussi
email: senoussih@yahoo.fr

Feature Selection for fault detection systems : application to the Tennessee Eastman Process

Keywords:

In fault detection systems, massive amount of data gathered from the life-cycle of equipment is often used to learn models or classifiers that aims at diagnosing different kind of errors or failures. Among this huge quantity of information, some features (or sets of features) are more correlated with the kind of failures than others. The presence of irrelevant features might affect the performance of the classifier. To improve the performance of a detection system, feature selection is hence a key step. We propose in this paper an algorithm named STRASS, that aims at detecting relevant features for classification purposes. In certain cases, when there exists a strong correlation between some features and the associated class, classical feature selection algorithms fail at selecting the most relevant features. In order to cope with this problem, STRASS algorithm makes use of k-way correlation between features and the class to select relevant features. To assess the performance of STRASS, we apply it on simulated data collected from the Tennessee Eastman chemical plant simulator. The Tennessee Eastman process (TEP) has been used in many fault detection studies and three specific faults are not well discriminated with classical algorithms. The results obtained by STRASS are compared to those obtained with reference feature selection algorithms. We show that the features selected by STRASS always improve the performance of a classifier compared to the whole set of original features and that the obtained classification is better than with most of the other feature selection algorithms.

Introduction

Fault detection has been extensively studied over the last few decades using various techniques. Tyan et al. [START_REF] Tyan | An application on intelligent control using neural network and fuzzy logic[END_REF] give a retrospective of the different methods such as parameter estimation, state observation schemes, pattern recognition techniques and artificial intelligence methods. Fault detection techniques can be divided into three categories: modelbased, knowledge-based [START_REF] Riverol | Integration of fault diagnosis based on case-based reasoning principles in brewing[END_REF][START_REF] Morello | A knowledge discovery process for a flexible manufacturing system[END_REF] and data-driven methods [START_REF] Verron | Fault detection and identification with a new feature selection based on mutual information[END_REF][START_REF] Tyan | An application on intelligent control using neural network and fuzzy logic[END_REF][START_REF] Sugumaran | Feature selection using decision tree and classification through proximal support vector machine for fault diagnostics of roller bearing[END_REF][START_REF] Torkkola | Sensor selection for maneuver classification[END_REF][START_REF] Widodo | Application of nonlinear feature extraction and support vector machines for fault diagnosis of induction motors[END_REF]. The approach proposed in this paper is data-driven and is based on machine learning tools.

A fault detection system, as shown in Fig. 1, can be identified from the knowledge data discovery process (KDD) where the output gives the state of the system (faulty or nonfaulty for instance). The input data correspond to recorded sensor measurements that are considered as features.

The data-mining component relies heavily on classical techniques from the fields of machine learning, pattern recognition, and statistics to find relevant patterns from the data and transform them into useful task-oriented knowledge.

[Fig. 1 about here.] Knowledge discovery and data mining have emerged as some of the most significant and fast expanding research areas and have found many successful real-world applications in a variety of disciplines like fault detection. For instance, Casimira at al. [7] develop a Knearest neighbors algorithm to identify stator and rotors faults in induction motors. From the thirty-one features extracted by time-frequency analysis of the stators currents and voltages, six relevant features are selected by a sequential backward algorithm from the initial subset. Experimental results demonstrated the effectiveness of this method in condition monitoring. Sugumara at al. [START_REF] Sugumaran | Feature selection using decision tree and classification through proximal support vector machine for fault diagnostics of roller bearing[END_REF] focused particularly on fault conditions in roller bearings of a rotary machine. They used vibration signals from different functional mode (good bearing, bearing with inner race fault, bearing with outer race fault, and inner and outer race fault). First, a set of eleven features were extracted by time-frequency analysis. Among these, four best features were selected from a given set of samples using C4.5 decision algorithm. Secondly, Proximal Support Vector Machine (PSVM), was used to efficiently classify the faults. Yang et al [START_REF] Yang | Support Vector Machine for Machine Fault Diagnosis and Prognosis[END_REF] presents a survey of fault diagnosis using Support Vector Machines (SVM) combined with other methods. In a similar study on fault detection of roller bearings, Jack at al. [START_REF] Jack | Genetic algorithms for feature selection in machine condition monitoring with vibration signals[END_REF] used Genetic Algorithms to select an optimal feature subset for SVM and artificial neural network based classifiers.

In chemical process industry large amounts of variables are monitored, which makes feature selection an important topic for that kind of applications. The Tennessee Eastman Process (TEP) benchmark has been the object of many studies in the literature [START_REF] Downs | A plant-wide industrial process control problem[END_REF][START_REF] Ricker | Decentralized control of the tennessee eastman challenge process[END_REF]. L. Wang and J. Yu have proposed in [START_REF] Wang | Fault feature selection based on modified binary pso with mutation and its application in chemical process fault diagnosis[END_REF] a binary Particle Swarm Optimization with mutation (MBPSOM) combined with Support Vector Machine (SVM) to select the most pertinent for fault diagnosis. Chiang et al. [9] applied Fisher Discriminant Analysis (FDA), Discriminant Partial Least Squares (DPLS) and Principal Component Analysis (PCA). Nashalji et al. [START_REF] Nashalji | Fault detection of the tennessee eastman process using improved pca and neural classifier[END_REF] used Genetic Algorithm and PCA to determine the main principle components and then used a neural networks-based classifier to detect faults during the operations of the industrial process. Verron et al. [START_REF] Verron | Fault detection and identification with a new feature selection based on mutual information[END_REF] proposed a fault diagnosis procedure based on discriminant analysis and mutual information. In order to obtain accurate classification performance, feature selection is performed with an algorithm based on the mutual information between variables. P. Cui et al. [START_REF] Cui | Improved kernel principal component analysis for fault detection[END_REF] applied Kernel principal component analysis (KPCA) for analysis of the TEP, and improved the fault detection of KPCA by a Fisher discriminant analysis scheme. In the TEP benchmark, three kinds of fault (4,9 and 11) are difficult to distinguish because they are strongly correlated. To improve the efficiency of fault detection systems with strongly correlated data, we develop in this paper a detection system based on a feature selection algorithm devoted to detect interactions between the features and the class. The presence of irrelevant and/or redundant features affects the speed and accuracy of learning or data mining algorithms. Therefore the selection of relevant information is very important for the development of a comprehensive and robust model and also for speeding up the training phase, hence reducing the cost as well as the data collection time for the classifier [START_REF] Almuallim | Learning with many irrelevant features[END_REF][START_REF] Blum | Selection of relevant features and examples in machine learning[END_REF][START_REF] Dash | Consistency based feature selection[END_REF][START_REF] Liu | Feature Selection for Knowledge Discovery and Data Mining[END_REF]. This paper takes up the feature interaction challenge. We are specifically interested in detecting partial correlations among variables. Two criteria are proposed (see section 3) respectively based on the weak and the strong relevance: the discriminating capacity measure (DC) and the discriminating capacity gain measure (DCG). These criteria are designed to focus on detecting k-way interactions in the data (i.e., interactions between sets of k features and the class) and therefore on detecting features that have the exclusiveness to discriminate concepts (also calles unavoidable features). These criteria are established in a greedy type algorithm named STRASS (STrong Relevant Algorithm of Subset Selection). This algorithm has the ability to treat partially correlated data. In order to highlight the effectiveness of the proposed algorithm, and assess its capability to detect partial correlations, STRASS is first tested on artificial data sets which are well known for their feature interactions. Indeed, those data sets have challenged many feature selection algorithms. STRASS is then applied to the TEP benchmark data, by feeding different classifiers with the set of features selected by STRASS. Experimental results are compared with reference algorithms as CFS, FCBF, mRMR, INTERACT, LASSO, ReliefF and SVM-RFE. They highlight the ability of the proposed algorithm to select relevant features for classification purposes.

The rest of the paper is organized as follows: Section 2 overviews some of the feature selection methods used for fault detection. Section 3 introduces the selection criteria uponn which STRASS is built. A comprehensive description of STRASS is presented in Section 4. In Section 5, STRASS is evaluated on synthetic data sets and on the TEP benchmark and compared with well-known feature selection algorithms. Conclusions are drawn in Section 6.

Overview of feature selection methods

Different feature selection strategies have been proposed and used in the fault detection context. Liu [START_REF] Liu | Feature Selection for Knowledge Discovery and Data Mining[END_REF][START_REF] Liu | Toward integrating feature selection algorithms for classification and clustering[END_REF], Blum and Langley [START_REF] Blum | Selection of relevant features and examples in machine learning[END_REF] compared different approaches dealing with data selection and emphasized four major points: the starting point in the feature space, the search strategy, the evaluation criterion and the stopping criterion.

-The starting point in the feature space could begin with no features with a forward search, all features with a backward search, or a random subset of features with a bidirectional search. Consequently, features can be successively added or removed by a certain procedure. -The search strategy for feature subsets can be done by random heuristics or complete procedure. -The evaluation criteria is an important component of any feature selection method, as it measures the goodness of a specific subset. An evaluation criteria can be categorized into three main groups based on their dependency on mining algorithm: filters, wrappers and embedded methods. Filters operate independently of any mining algorithm contrary to wrapper methods which use the performance of the mining algorithm. Embedded methods are built upon a data mining or a classification algorithm to perform the feature selection.

We propose to categorize feature selection algorithms depending of their evaluation criteria into three main groups based on how the interaction between features is treated [START_REF] Michaut | Filterign and variable selection in learning processes[END_REF]:

-The myopic measures which estimate the feature quality independently of the others features. Most of the existing measures belong to the first category that is why Kira and Rendel [START_REF] Kira | The feature selection problem: Traditional methods and a new algorithm[END_REF] and Kononenko [START_REF] Kononenko | Estimating attributes: Analysis and extensions of relief[END_REF] underlined the difficulties for classifiers to work with correlated data.

-Semi-contextual measures consider low order 2-way (one feature and the class) and 3way (two features and the class) interactions.

-Contextual measures consider k-way (k features and the class, k > 2) interactions.

The myopic measures

Algorithms based on myopic criteria Relief [START_REF] Kira | The feature selection problem: Traditional methods and a new algorithm[END_REF], B&B [START_REF] Narendra | A branch and bound algorithm for feature subset selection[END_REF], LVF [START_REF] Dash | Consistency based feature selection[END_REF], FOCUS [START_REF] Almuallim | Learning boolean concepts in the presence of many irrelevant features[END_REF] do not detect the correlations between the features and the class, unlike those using semi-contextual or contextual criteria. Indeed, most of the works in statistics make the features independence assumption. Relief, the most powerful individual feature selection algorithm, scores individual features rather than scoring feature subsets, those features with scores exceeding a user-specified threshold are selected for the final subset. The actual challenge in feature selection is to study feature interactions with relevant measures to select the optimal subset with maximum relevance and minimum redundancy.

Semi-contextual measures

CFS [START_REF] Hall | Correlation-based feature selection for discrete and numeric class machine learning[END_REF], mRMR [START_REF] Peng | Feature selection based on mutual information: criteria of maxdependency, max-relevance, and min-redundancy[END_REF], FCBF [START_REF] Yang | Support Vector Machine for Machine Fault Diagnosis and Prognosis[END_REF] and ReliefF [START_REF] Kononenko | Estimating attributes: Analysis and extensions of relief[END_REF]19] use semi contextual measures. CFS calculates a feature subset merit, thus detecting the best feature combination. The algorithm is powerful as long as the interaction between features is not too large. mRMR feature selection algorithm selects features that should be both minimally redundant among themselves and maximally relevant to the target classes. The optimal subset maximizes the distance between the two profits. FCBF uses a correlation measure based on the information gain to detect the redundancy between features. They chose the Symmetrical Uncertainty (SU).

The algorithm involves two steps: the first one calculates the SU value for each features, selects and orders relevant ones according to a predefined threshold, and the second one selects predominant features. However none of these algorithms treats high order interactions (k features and the class) like a contextual measure. Genetic Algorithm (GA) based feature selection algorithms [START_REF] Casillas | Genetic feature selection in a fuzzy rule-based classification system learning process for high dimensional problems[END_REF][START_REF] Lanzi | Fast feature selection with genetic algorithms: a filter approach[END_REF] have been successfully used as a variable subset selection and attempt to address the variable associations and various effects. Feature interaction is indirectly taken into account via the selection of a set of variables which is determined by a fitness function, generally based on the classifier accuracy. However, it should be noted that this method is a wrapper method unlike other methods. The LASSO method [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] is an embedded algorithm. It constructs a model and penalizes coefficients, shrinking many of them to zero. The principle of LASSO is linked with L1-norm regularization techniques which aim at penalizing complex models, upon which is also built the Elastic Net approach [START_REF] Zou | Table 3 Accuracy performance of different classifiers associated with feature selection algorithms[END_REF]. SVM-RFE [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] is also one of the most famous embedded method. Weights are assigned to features while building the model, and the ones with smaller weights are recursively eliminated.

Contextual measures

Two more recent algorithms treat and capture k-way interacting features : INTERACT [START_REF] Yu | Efficient feature selection via analysis of relevance and redundancy[END_REF] and STRASS [START_REF] Morello | A knowledge discovery process for a flexible manufacturing system[END_REF][START_REF] Michaut | Filterign and variable selection in learning processes[END_REF][START_REF] Senoussi | A new contextual based feature selection[END_REF]. To do so, INTERACT combines an information measure and a consistency measure. In the first part of the algorithm, the features are ranked in descending order based on their symmetrical uncertainty (SU) values. In the second part, features are evaluated according to their C-contribution which relies on the calculation of inconsistency rate. The features are evaluated one by one starting from the end of the ranked feature list.

The strong points of these algorithms are their effectiveness to deal with diverse problems like, modal, continuous, noisy and correlated data.

Criteria of relevance and redundancy

In this section, we define two new criteria for feature selection, that are elaborated from the class discriminatory power, in a pair-wise data representation approach. The STRASS feature selection algorithm (cf Section 4) is built upon these criteria. These criteria are designed to take into account k-way interactions in the data and hence lead to a contextual measure.

Data representation

Let the input data Ω consist of n samples ω 1 ; : : : ω n . Every sample in Ω is composed of r features. The set of features is denoted x = fx 1 ; : : : ; x r g. Every sample in Ω is labeled with a class c P C = fc 1 ; : : : ; c M g. In the following, the notation x k (ω i ) represents the value of feature x k in ω i and C(ω i ) represents the class of sample ω i . Let us associate to a feature x k the Boolean function φ k i j ; 1 i; j n; 1 k r :

φ k i j : Ω ¢ Ω 3 f0;1g (ω i ; ω j ) U 3 1 if x k (ω i ) = x k (ω j ) 0 otherwise: (1)
Let us also define the function φ C i j :

φ C i j : Ω ¢ Ω 3 f0;1g (ω i ; ω j ) U 3 1 if C(ω i ) = C(ω j )
0 otherwise:

(2)

Weak Relevance measure

The weak relevance of a set of features is defined by the number of all pairs of samples who have at least one discriminating variable and different labels or different distributions of labels. According to that definition, let us define the discriminating capacity (DC) measure of a feature set L = (x 1 ; : : : ; x m ) with the following formula :

DC(L; Ω ) = n ∑ i=1 n ∑ j=1 m ∏ k=1 φ k i; j φ C i; j (3)

Strong Relevance measure

To measure the exclusiveness of a feature to describe a concept, the equivalent of a "relevance gain" is defined as the measure related to a feature compared to a subset of features. The strong relevance (SR) of a feature x k on pairs of instances is defined as the relevance of a feature x k compared to a relevant preselected features subset L = (x 1 ; : : : x m ). This measure is given by:

SR(x k ; L; ω i ; ω j ) = φ C i; j φ k i; j m ∏ l=1 φ l i; j (4)
The aggregation of the Strong Relevance (SR) expression on the whole pairs obtained by the sample Ω of n patterns will define the Discriminating Capacity Gain (DCG) as:

DCG(x k ; L; Ω ) = n ∑ i=1 n ∑ j=1 SR(x k ; L; ω i ; ω j ) (5)
The DCG of a feature x k for a set of n objects compared to a set L of features is equal to the number of object couples discriminated by only x k and no other features.

Redundancy of a feature

A feature x k is said to be redundant in a feature subset L if the discriminating capacity measure of the set L fx k g is the same as the one of L. In other words, x k is said to be redundant if

DC(L; Ω ) = DC(L fx k g;Ω) (6)

Interest of the two criteria

The two suggested contextual criteria upon which STRASS is built are complementary. They can detect not only strongly relevant features but also weakly relevant ones. The first one calculates the discriminating capacity of a set of variables, and aims at extracting a subset of variables with the same DC as the entire set. The second one evaluates the discriminating capacity gain of one feature relatively to a set of features. This contextual criterion aims at detecting either the relevance or the redundancy of a feature compared to a subset of features.

Features with the largest gains are integrated into the selected set, while the redundant ones (with a null gain) are discarded. These criteria have the particularity to be calculated on a restricted set of object pairs and variables. It can hence detect partial correlations between features, and study the k-way interactions between them and the class. The criteria which we have developed under pair data set allow us to explore three aspects of the correlation:

1. the feature correlation on a pair-wise data set, i.e. the feature capacity to discriminate a part of the studied population. 2. the partial feature correlation relatively to a set of features. Redundant features (that play the same role of another feature) are searched to be excluded. On the other hand, a feature can be considered weakly relevant to the class when evaluated alone, but becomes very relevant when combined with other features. This correlation is also searched with these criteria. 3. the feature capacity to be the only one to discriminate a population subset. Such features are called strongly relevant (or predominant).

STRASS (STrong Relevant Algorithm of Subset Selection

Description of the algorithm

STRASS is based on the contextual criteria established in Michaut thesis [START_REF] Michaut | Filterign and variable selection in learning processes[END_REF][START_REF] Morello | A knowledge discovery process for a flexible manufacturing system[END_REF], and described in the previous section. The first criterion measure is a discriminating capacity (DC) of a set of variables and the second criterion is DCG (Discriminating Capacity Gain) measure. These criteria when associated with a greedy algorithm allow:

-To capture k-way interacting features -To detect the partially redundant features. Partial correlation exists when only a feature combination can discriminate the class -To rank the variables selected with the complementary criteria (DC: discriminating capacity).

The STRASS algorithm proposed here belongs to the greedy type category of algorithms. The research is a sequential bidirectional generation, i.e. a core of features is composed from an empty set S f which is built gradually until a subset having the same degree of relevance as the starting subset noted is obtained. The feature subset is progressively computed and re-evaluated at every feature addition. The algorithm breaks up into three steps depending on its initialization:

Step 1: The features are ranked in descending order based on their discriminating capacity gain and a subset of strongly relevant features or predominant features is selected.

Step 2: The remaining features are evaluated one by one starting from the top of the remaining features list and the weakly relevant features which have the largest discriminating capacity are combined with the previously selected features S f if the resulting overall discriminating power is increased. In fact a feature may have a little correlation with the class, but when it is combined with some others features, the resulting subset can be strongly correlated with the class.

Step 3: Suppression of redundant features. At this stage, backward elimination is employed to detect the features that become redundant compared to the preselected features subset S f when adding a new feature in the second stage of the algorithm. For a predefined threshold ρ; 0 < ρ < 1 features having a discriminating capacity DC < ρ ¢ DCtot are removed. Therefore, we obtain a subset of low cardinal and with no redundancy in the selected set of features.

The complete algorithm is given in Figure 2.

Strengths and weaknesses of STRASS

Amongst feature selection algorithms, STRASS is based on contextual criteria. These criteria are established in a greedy type algorithm and select from a learning set a minimal set of relevant features. The learning set must be consistent, not having missing data and consists of symbolic data and/or numeric. Noisy data have a negative impact on the associated performance. To our knowledge, the criteria used by STRASS are the first to take into account the different aspects of the correlation detailed above. STRASS identifies the k-way interaction between features and the partial correlations and partial redundancy in a set of pairwise data. Another interesting remark is the fact that STRASS is a filter algorithm which is less computationally intensive than wrapper techniques, and embedded methods.

STRASS is computationally efficient on databases with reasonable sizes (in the order of thousands entries and hundreds of features). However for large databases, the pairwise data representation is inherently a combinatorial problem, and is not adapted. In order to reduce the complexity of STRASS, we plan to simplify the criteria and express them under a contingency form. The transformation of the pairwise criteria to contingency criteria with the help of Marchotorchino [START_REF] Marcotorchino | Utilisation des comparaisons par paires en statistique des contingences[END_REF] might be of interest and will be the object of a future work.

[Fig. 2 about here.] 5 Experiments and results

Implementation

STRASS algorithm was implemented in MATLAB 7.5 environment. For the filtering algorithms and classifiers, existing tools in WEKA machine learning platform [START_REF] Hall | The weka data mining software: An update[END_REF] have been used. The experiments were run using WEKA with its default values. For evaluation purposes we will compare STRASS with other algorithms that consider feature interaction and correlation using contextual and semi-contextual measures, such as CFS, mrMR, CFS-FCBF, INTERACT, ReliefF, SVM-RFE, LASSO, Elastic Nets. The mRMR feature selection algorithm used in this work has been adapted from [START_REF] Peng | Feature selection based on mutual information: criteria of maxdependency, max-relevance, and min-redundancy[END_REF] and has been downloaded online 1 . Same conditions are used in Matlab and Weka.

Evaluation of the algorithm

STRASS performance was assessed in two complementary ways: (i) Direct evaluation through artificial data sets: Led, Monk, Bool, Parity, Corral and Agrawal's functions (Appendix B). These data sets highlight the behavior of our algorithm when the descriptive characteristics interact. Langley and Sage [START_REF] Langley | Computational learning theory and natural learning systems[END_REF] stressed the importance of artificial fields. (ii) Indirect evaluation on TEP benchmark allows us to study the classification performance with and without the filtering phase, the classifiers accuracy and the number of features removed by the filtering algorithm. The results obtained with STRASS are compared with filter methods, as mRMR, CFS (with best first search), FCBF (with threshold SU set to 0), INTERACT, ReliefF algorithms, and with embedded methods, as the LASSO principle [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], the Elastic Net principle [START_REF] Zou | Table 3 Accuracy performance of different classifiers associated with feature selection algorithms[END_REF] the SVM-RFE algorithm [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF]. To compare STRASS with mRMR and ReliefF, we take the same number of features for mRMR as the one selected by STRASS.

Synthetic data with known feature interaction

[Table 1

about here.]

We have applied STRASS on some synthetic data sets widely used for feature selection. Results are given in Table 1. Only STRASS is able to determine relevant features in all data sets. The INTERACT algorithm selects relevant features nine times out of thirteen, which is more than for the other algorithms. It can be seen from Table 1 that for LED Display Domain (Led, Led 24) all the algorithms fail to detect the relevance and the sufficiency of the five segments except STRASS. For the Parity data set, STRASS and INTERACT establish that the features x 1 ; x 2 ; x 3 are relevant whereas the others are useless. For the three MONKs data sets, STRASS was able to find true relevant features using a loss of 1% of the DCtot for MONK-3. STRASS, INTERACT and mRMR detected the redundancy of feature x 6 in Corral data set. This feature being correlated to 75% with the feature class is considered to be relevant alone by CFS and FCBF because the algorithms cannot evaluate the k-way interactions. In the case of Agrawal's functions (F1 to F4), STRASS gives the relevant features for the four functions. Thus, STRASS is the most powerful algorithm on these data sets. This is explained by the presence of large interactions in the data. Indeed CFS, mRMR and FCBF detect the pair-wise feature-feature correlation (inter-correlation) but they cannot identify neither interactions between subsets of features and the class nor unavoidable features.

Application on the TEP Benchmark

[Table 2 about here.] TEP benchmark [START_REF] Downs | A plant-wide industrial process control problem[END_REF] is described in Appendix B. Instances in this data set are composed of 52 variables and labeled with a fault type between 1 and 15. Three types of faults denoted as fault 4, fault 9 and fault 11 are considered here, because they are the most difficult to classify. The problem of their identification is due to a great interaction between the features for these faults, not to the classifier used. To solve this problem and improve the classifier performance, we need define a feature filtering taking account feature interactions. Table 2 presents the features that are selected by the different algorithms compared here. For STRASS, the most discriminating features are f51;41;38;40;37;50;19;18;9;20g.

STRASS uses DCG (discriminatory capacity gain) to rank each feature are sorted according to its contribution relatively with other features. Fig. 3 shows the DCG for each predominant feature.

[Fig. 3

about here.]

[Table 3

about here.]

For the classification task, five different classifiers have been used : a decision tree (C4.5), 1-nearest-neighbor (IB1), Naive Bayes (NB), a multilayer perceptron (MLP) and Support Vector Machines (SVM). We have applied these classifiers to the features selected by the feature selection algorithms described above. Results are obtained with 10-fold cross validation. We have compared the results with semi-contextual and contextual methods : mRMR, CFS, FCBF, INTERACT, ReliefF, LASSO, Elastic net and the SVM-RFE feature selection algorithms. To analyze the results obtained in this study, we have employed two performance measures: accuracy and Cohen's Kappa. The accuracy is defined as the number of successful hits relative to the total number of classifications. The Kappa is a statistical measure of inter-rater agreement. Cohen's Kappa can be adapted to classification tasks and it is also used in some well-known software packages, such as WEKA. Tables 3 and4 show the performance in terms of accuracy and Kappa respectively. The best results for each classifier are highlighted in bold. The symbols + and respectively indicates an improvement or a degradation in terms of performance compared with the ones obtained for the whole set of features. Let us examine the effect of feature selection on classification performance. Classification accuracy and kappa score is calculated before and after filtering. These results are reported in Tables 3 and4. An interesting point is that the features selected by STRASS always improve the performance of a classifier compared to the complete set of features. This is due to the fact that the feature subset selected by STRASS has the same discriminating capacity as the full set of features and hence the classification performance are equal or better than the one of the full set. The performance obtained with STRASS is always one of the two best feature selection algorithms apart from when using with the Multi-Layer perceptron. Its average position is 2:6 which is the best among the other algorithms. The best performance in terms of accuracy is obtained with IB1 and ReliefF algorithms (98:91%), but STRASS is just behind that performance (98:56%). Moreover, the feature selection made by the ReliefF algorithm does not always improve the performance compared to the whole features on the contrary to STRASS.

In terms of average classification performance, the best algorithms are STRASS, mRMR and INTERACT. These results hence highlight the benefit of detecting k-way correlations (particularly for k < 2, as STRASS and INTERACT) compared to only detecting featurefeature intercorrelation, as most of the semi-contextual and embedded consider. The characteristic of detecting k-way correlations lead to improve classification performance each time on these experiments.

[Table 4 about here.] [Fig. 4 about here.] Fig 4 gives the accuracy results obtained with STRASS with ordered selected features in decision tree (C4.5) and nearest neighbors classifier (IB1). For IB1 the best results are obtained with all the selected features whereas for C4.5 the seven first selected features give the best result.

The confusion matrices obtained with STRASS + IB1 and ReliefF + IB1 for their associate best feature set are given respectively in Table 5 and6. In all these cases, the fault 11 is less discriminated because this fault overlaps with the two others.

To conclude this experiment section, we have seen that STRASS gives better or equivalent performance than most of the compared feature selection algorithms when combined with different classifier. An interesting point is that the features selected by STRASS always improve the performance of a classifier compared to the complete set of features, which is not always the case for the other feature selection algorithms. These results hence highlight the benefit of the k-way feature selection process of STRASS.

[Table 5

Conclusion

This paper describes STRASS, a contextual-based feature selection algorithm for classication purposes able to detect the interaction between features and the class and select a minimum relevant feature subset. The efficiency and effectiveness of STRASS to handle large interactions are demonstrated through a comparative study with other representative feature selection algorithms on synthetic data known for their correlation. The proposed feature selection algorithm was then applied to a well-known fault detection benchmark: the Tennessee Eastman Process (TEP). STRASS has demonstrated its ability to reduce the dimensionality of data sets while maintaining or improving the performances of learning algorithms. In fact for TEP process, the features selected by STRASS decreased the data correlation and the overall misclassification for the testing set using 1-nearest neighbor decreased further to 1.4%. STRASS was also compared to other reference feature selection algorithms. The results of STRASS outperformed or lead to equivalent performance to those obtained with those reference algorithms. In addition, the features selected by STRASS always improve the performance of a classifier compared to the whole set of features.

A The Tennessee Eastman Process

[Table 7 about here.] The Tennessee Eastman Process (TEP) is a chemical process, created by the Eastman Chemical Company to provide a realistic industrial process in order to evaluate process control and monitoring methods [START_REF] Downs | A plant-wide industrial process control problem[END_REF]. This process was simulated on Matlab by Ricker [START_REF] Ricker | Decentralized control of the tennessee eastman challenge process[END_REF]. The simulator was used to generate overlapping data sets to evaluate the classification performance. Figure 5 shows a flow sheet of TEP. There are four unit operations: an exothermic two phase reactor, a flash separator, a re-boiler striper, and a recycle compressor. The TEP process produces two products (G and H) and one (undesired) by-product F from four reactants (A, C, D and E). This process has 12 input variables and 41 output variables. Only 52 variables are taken into account in this problem because one of the input variables (the reactor agitator speed) is constant. The system has fifteen types of identified faults. In this paper, we considered only three types of fault : fault 4, 9 and 11. These faults are described in Table 7.

[Fig. 5 about here.]

B Synthetic data

We describe in this appendix the synthetic data used in this paper for simulation purposes. The LED display domain data set is available on the UCI data set repository [START_REF] Bache | UCI machine learning repository[END_REF]. The MONK's problems [START_REF] Thrun | The MONK's problems: A performance comparison of different learning algorithms[END_REF] are composed of three target concepts : MONK-1 : (x 1 = x 2 ) (x 3 = 1) MONK-2 : exactly two of :

fx 1 = 1; x 2 = 1; x 3 = 1; x 4 = 1; x 5 = 1; x 6 = 1g MONK-3 : (x 5 = 3 x 4 = 1) (x 5 T = 4 x 2 T = 3)
The BOOL data set is composed of a function of 6 Boolean features giving a Boolean class, for instance : y class = (x 1 ¨x2 ) (x 3 x 4 ) (x 5 x 6 ). Six other randomly generated Boolean features are added to these features.

The Parity data set is composed of a of a function of 3 Boolean features y class = x 1 ¨x2 ¨x3 . Seven randomly generated Boolean features are added. This data set is particularly interesting because no relevant features taken separately can be distinguished from irrelevant ones. The Parity2 data set is the same as the Parity data set to which 2 redundant features are added : x 11 = x 1 and x 12 = x 2 . This data set allows testing the algorithms ability to work with redundant features. The Coral data set is composed of six binary features x 1 to x 6 among which x 5 is irrelevant and x 6 is correlated to 75% with the feature class y class = (x 1 x 2 ) (x 3 x 4 ).

Agrawal's functions are a series of classification functions of increasing complexity that uses 9 features to classify people into different groups. More details can be found in [START_REF] Agrawal | An interval classifier for database mining applications[END_REF]. 3 Accuracy performance for ordered selected features with STRASS and two classifiers : IB1 and C4.5 q q q q q q q q q q 2 4 6

Fig. 1 1 .Fig. 2

 112 Fig. 1 Knowledge data discovery process

Fig.

  Fig. 3 Accuracy performance for ordered selected features with STRASS and two classifiers : IB1 and C4.5

Fig. 4

 4 Fig.4Accuracy performance for ordered selected features with STRASS and two classifiers : IB1 and C4.5

Fig. 5

 5 Fig. 5 Process flow sheet of TEP

  about here.] [Table 6 about here.]

Table 1

 1 Features selected by different algorithms on synthetic data sets ; x 4:5 x 2 ; x 4:5x 2 ; x 5:6 x 2 ; x 5:6 x 2 ; x 5:6 x 1:2 ; x 4:5x 1:6 Parityx 1:3 x 1:3 x 2 ; x 6 ; x 10 x 5 ; x 8 ; x 10 ; x 9 x 1:2 ; x 9x 1:2 ; x 8 x 9x 1 ; x 9 x 1:2 ; x 8 x 1:2 ; x 9

	data sets	Relevant	STRASS	mRMR	CFS	FCBF	INTERACT	ReliefF
		features
	LED7	x 1:5	x 1:5	x 2:4 ; x 6:7	x 1:7	x 1:7	x 1:7	x 1:7
	LED 24	x 1:5	x 1:5	x 2:5 ; x 7	x 1:8 ; x 14	x 2:7 ; x 11	x 1:9 ; x 11	x 1:7
					x 16 ; x 19 ;	x 12 ; x 14 ;	x 14:20 ;
					x 21	x 15 ; x 19:21	x 22:23
	MONK1	x 1:2 ; x 5	x 1:2 ; x 5	x 1 ; x 4:5	x 1 ; x 3:5	x 1 ; x 3:5	x 1:2 ; x 5	x 1:5
	MONK2	x 1:6	x 1:6	x 1:6	x 4:6	x 4:6	x 1:6	x 1:6
	MONK3	x 2 x 10	x 1:3	x 1:3
	Parity 2	x 1:3	x 1:3	x 2 ; x 6 ; x 10	x 5 ; x 8 ; x 10	x 10	x 1:3	x 1:3 ; x 11:12
	Corral	x 1:4	x 1:4	x 1:4	x 1:4 ; x 6	x 1:4 ; x 6	x 1:4	x 1:4 ; x 6
	Bool	x 1:6	x 1:6	x 3:6 ; x 7 ;	x 1:6	x 3:6 ; x 10	x 1:6	x 1:6 ; x 7 ;
				x 10		x 12		x 12
	F1	x 3	x 3	x 3	x 3	x 3	x 3	x 3
	F2	x 1 ; x 3	x 1 ; x 3	x 1 ; x 8	x 1	x 1	x 1	x 1:3
	F3	x 1 ; x 3:4	x 1 ; x 3:4	x 2 ; x 4 ; x 8	x 2:4	x 2:4	x 2:4	x 1:9
	F4	x 1:2

Table 5

 5 Confusion matrix obtained with STRASS+IB1``````````T

	rue Class	Estimated class	Fault 4	Fault 9	Fault 11
		Fault 4	1430	6	4
		Fault 9	1	1435	4
		Fault 11	27	20	1393

Table 6

 6 Confusion matrix obtained with ReliefF+IB1``````````T

	rue Class	Estimated class	Fault 4	Fault 9	Fault 11
		Fault 4	1432	1	7
		Fault 9	1	1435	4
		Fault 11	20	14	1406

Table 7

 7 Description of the faults used in this paper

	Fault number	Description
	4	Step change in the reactor cooling water
		inlet temperature
	9	Random variation in D feed temperature
	11	Random variation in the reactor cooling
		water inlet temperature

http://penglab.janelia.org/software/

List of Figures


Simon Malinowski
email: simon.malinowski@femto-st.fr

Brigitte Chebel-Morello

Noureddine Zerhouni
email: noureddine.zerhouni@femto-st.fr

Remaining Useful Life estimation based on discriminating shapelet extraction

In the Prognostics and Health Management (PHM) domain, estimating the remaining useful life (RUL) of critical machinery is a challenging task. Various research topics as data acquisition, fusion, diagnostics, prognostics and decision are involved in this domain. This paper presents an approach to estimate the Remaining Useful Life of equipment based on shapelet extraction. This approach makes use, in an offline step, of a history of run-to-failure data to extract discriminative rul-shapelets, i.e. patterns that are correlated with the RUL of the considered equipment. A library of rul-shapelets is hence extracted at this step. Then, in an online step, these rul-shapelets are compared to testing units and the ones that match these units are used to estimate their RULs. Therefore, the RUL estimation of a testing unit is based on patterns that have been selected for their high correlation with the RUL. This approach is hence different from classical similaritybased approaches that attempt to match complete testing units (or only late instants of testing units) with training ones to estimate the RUL. The performance of our approach is assessed with a case study on the remaining useful life estimation of turbofan engines and performance is compared with other similarity-based approaches.

Introduction

Remaining Useful Life estimation is one of the main tasks in the Prognostics and Health Management domain. The aim of any RUL estimation technique is to provide accurate prediction of the time after which an equipment will not be able to meet its operating requirements. RUL estimation is hence very important for industrial purposes as it can help in planning maintenance strategies, maximizing the useful operational life of equipment, reducing maintenance costs and avoiding breakdowns that might have critical impacts.

RUL estimation techniques in the literature are separated into three families : model-based, data-driven and hybrid approaches. Model-based approaches rely on building a physical model describing the degradation behavior of the equipment. For instance, stochastic filtering [START_REF] Mytyri | Application of stochastic filtering for lifetime prediction[END_REF], particle filtering [START_REF] Cadini | Model-based monte carlo state estimation for condition-based component replacement[END_REF] have been used to model degradation in the case of fatigue crack growth. This kind of approaches is very accurate but requires the knowledge of the physical degradation of the system to be available, which is not always the case. In addition, in most of the cases, the equipment is composed of many components for which a different model needs to be defined. Furthermore, model-based approaches are closely related to the considered application. They hence lack genericity. Data-driven approaches make use of available run-to-failure data to extract relevant information based on a learning process. These approaches do not attempt to derive an analytical model from the data, but attempt to capture the relationship between sensory data and the degradation level of an equipment, in order to estimate its RUL. These approaches are usually easier to obtain and implement, but are often less accurate than model-based ones. They hence offer a trade-off between accuracy, complexity and applicability. Most of the data-driven approaches in the literature are based on machine learning, probabilistic or statistical tools. A good survey of machine learning and statistical techniques for prognostics can be found in [START_REF] Sikorska | Prognostic modelling options for remaining useful life estimation by industry[END_REF]. Neural Networks have been widely considered to estimate the RUL [START_REF] Gebraeel | Residual life predictions from vibration-based degradation signals: a neural network approach[END_REF][START_REF] Huang | Residual life predictions for ball bearings based on selforganizing map and back propagation neural network methods[END_REF][START_REF] Vassilopoulos | Artificial neural networks in spectrum fatigue life prediction of composite materials[END_REF][START_REF] Heimes | Recurrent neural networks for remaining useful life estimation[END_REF]. Many other approaches based on neural networks techniques have been proposed. The authors of [START_REF] Peng | A modified echo state network based remaining useful life estimation approach[END_REF] have used echo state networks for RUL estimation, Javed et al. [START_REF] Javed | Novel failure prognostics approach with dynamic thresholds for machine degradation[END_REF] developed a method combining wavelet based neural network and fuzzy clustering tools, deep belief networks are proposed in [START_REF] Tamilselvan | Failure diagnosis using deep belief learning based health state classification[END_REF] for health state classification, and complex valued neural networks are used in [START_REF] Fink | Predicting component reliability and level of degradation with complexvalued neural networks[END_REF]. Other machine learning tools have also been applied for RUL estimation, for instance support vector regression in [START_REF] Loutas | Remaining useful life estimation in rolling bearings utilizing datadriven probabilistic e-support vectors regression[END_REF]. Probabilistic and statistical tools have also been used for prognostics applications. Bayesian learning techniques [START_REF] Wang | A generic probabilistic framework for structural health prognostics and uncertainty management[END_REF], Wiener processes [START_REF] Son | Remaining useful life estimation based on stochastic deterioration models: A comparative study[END_REF], copulas [START_REF] Xi | A copula-based sampling method for data-driven prognostics and health management[END_REF], Kalman filters [START_REF] Baraldi | A kalman filter-based ensemble approach with application to turbine creep prognostics[END_REF], Hidden Markov models [START_REF] Medjaher | Remaining useful life estimation of critical components with application to bearings[END_REF][START_REF] Moghaddass | An integrated framework for online diagnostic and prognostic health monitoring using a multistate deterioration process[END_REF] and autoregressive-moving-average models [START_REF] Wu | Prognostics of machine health condition using an improved arima-based prediction method[END_REF][START_REF] Caesarendra | Combined probability approach and indirect data-driven method for bearing degradation prognostics[END_REF] have also been considered in the literature. Hybrid approaches, for instance in [START_REF] Niu | Intelligent condition monitoring and prognostics system based on data-fusion strategy[END_REF][START_REF] Maio | Combining relevance vector machines and exponential regression for bearing residual life estimation[END_REF][START_REF] Liu | A data-model-fusion prognostic framework for dynamic system state forecasting[END_REF][START_REF] Zhao | Uncertainty quantification in gear remaining useful life prediction through an integrated prognostics method[END_REF], combine model-based and data-driven approaches. Similarly to model-based approaches, knowledge about the physical degradation of components is required for hybrid-ones, which lowers the applicability of such approaches.

[Figure 1

about here.]

In this paper, we focus on similarity-based approaches which are particular cases of data-driven ones. Fig. 1 depicts the general framework of such approaches. Sensory data (often multidimensional) is first processed : noise filtering, feature extraction, data fusion... Note that the kind of processing depends on the type of data. After these steps, data can sometimes be transformed into an 1-dimensional indicator, denoted health indicator (HI), that models the degradation of an equipment with a time-series (or trajectory). If the original sensory data already captures the health status evolution of the equipment, data can also be kept multidimensional. In Fig. 1, this step is depicted by the data processing box. After this step, data is formalized into instances. The kind of instances depend on the method : instances can represent complete HI trajectories, blocks of HI trajectories... Hence, in an offline phase, a library of instances is constructed from the training sensory data. In an online phase, the same processing and formalization operations are applied to testing sensory data (whose RUL is to be predicted) in order to obtain a new instance. This instance is then compared to the library of instances in the retrieval step and the most similar ones are selected in order to predict the RUL of the testing instance (RUL estimation box in Fig. 1).

The authors of [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] have won the 2008 PHM challenge with a similarity-based approach where instances were the HI's obtained after fusion of different sensor measurements and that relied on Euclidean distance between training and testing HI's to select most similar instances. In this method, the whole test trajectory is used to be matched to the library of training units. Another approach, in [START_REF] Ramasso | Joint prediction of continuous and discrete states in time-series based on belief functions[END_REF], only considers last parts of testing units HIs as instances for the matching, as last parts are more likely to be correlated with the degradation level. In this paper, we propose an alternative method for the instance formalization, retrieval and RUL estimation steps (gray boxes in Fig. 1), based on the extraction of relevant patterns from the training data.

Shapelets have been introduced in [START_REF] Ye | Time series shapelets: a new primitive for data mining[END_REF] and [START_REF] Xing | Extracting interpretable features for early classification on time series[END_REF] for classification and early classification of time series. We extend here this notion and define a new kind of shapelets called rul-shapelets, that correspond to patterns correlated to the remaining useful life of an equipment. In this paper, in the formalization step, discriminative rul-shapelets are first extracted from a training set of trajectories representing run-to-failure data of an equipment. This step produces a library of rulshapelets that will be used for RUL estimation of testing units. In the online phase, the testing instance is created and rul-shapelets from the library are searched in this testing instance. The rul-shapelets that are found in the testing instance provides an estimation about the RUL. Hence, the RUL estimation relies here on finding similar behaviors between patterns of the testing time series and the shapelets that have been extracted for their correlation with the remaining useful life. Hence, a major difference with other traditional approaches is that the RUL estimation is based on some parts of the test unit (and not the whole testing unit or only last instants), these parts being chosen because of their high correlation with the RUL.

This approach is compatible with any applications satisfying the following assumptions:

• Run-to-failure data is available

• Trajectories capturing the health status evolution can be obtained from sensory data

• Testing components are assumed to go through the same degradation process as training ones.

The rest of this paper is organized as follows. Section 2 gives an overview of research work on similarity-based RUL estimation technique and positions our proposed approach within this context. Section 3 describes how discriminative shapelets are extracted from a set of time series representing run-to-failure data. Section 4 explains how these shapelets are used to perform RUL estimation on testing units, and Section 5 evaluates the performance of the proposed approach on a case study.

Related work

In this section, related work about similarity-based approaches for RUL prediction is presented. We focus on explaining in details how the formalization, retrieval and RUL prediction steps are performed in the literature. Similarity-based approaches are relatively new for prognostics applications. Wang et al. [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] proposed a similarity measure between testing and training trajectories based on Euclidean distance. Instances are represented as 1-dimensional health indicator trajectories obtained by linear regression of the sensory data. In the online step, for a given training trajectory, the Euclidean distance between a testing trajectory and every sub-trajectories (shifted by different number of cycles) of the training one is computed and the minimum distance is retained together with the optimal shift. This procedure is repeated for every training trajectories. The best trajectories are selected according to minimum distances, up to a certain threshold. RUL is estimated from the RUL of the selected trajectories, using a weighted average. Xue et al. [START_REF] Xue | An instance-based method for remaining useful life estimation for aircraft engines[END_REF] proposed a fuzzy instance based approach for RUL prediction. The approach starts by building local fuzzy models for testing engines. The fuzzy model defines a cluster of peers in which each of these peers is a similar instance to the given engine with comparable operational characteristics. The final RUL estimate is obtained by aggregating the RULs of similar training instance via a similarity weighted sum. Zio et al. [START_REF] Zio | A data-driven fuzzy approach for predicting the remaining useful life in dynamic failure scenarios of a nuclear system[END_REF] proposed a similarity-based approach for prognostics using a fuzzy point-wise similarity defined for degradation trajectories. The distance score between two trajectories is based on a fuzzy membership function that maps the difference between trajectories elements into membership. Weights are inversely proportional to the distance score. RUL is then obtained as a weighted sum of RULs of similar instances. Ramasso et al. [START_REF] Ramasso | Joint prediction of continuous and discrete states in time-series based on belief functions[END_REF] proposed a method that jointly predicts observations (continuous states) and health states (discrete states) in order to estimate the remaining useful life of components. In this paper, the authors worked with multidimensional sensory data and the instances were chosen as multidimensional blocks of data. The instances hence are subsequences of the original sensory data. Instead of aggregating RULs of the most similar instances, the observation trajectories are aggregated to predict the future observations. Those observations are classified as health states and RUL is predicted as the time transition from the degrading to fault state. The retrieval phase of the algorithm is based on a Euclidean distance measure where only the last block of the testing trajectory is considered and compared to blocks of the library. In [START_REF] Khelif | Rul prediction based on a new similarityinstance based approach[END_REF], Khelif et al. proposed a new similarity measure that aims at giving more weight to last observations (as they are more likely to be correlated with the degradation level), while using the total testing HI trajectory as instances.

In these works, the retrieval step is hence based on trying to match testing trajectories (or testing sensory data), or blocks of trajectories with training ones in order to estimate the RUL. This implies that every piece of trajectories in the training set carry relevant information about their correlation with the degradation level, which might not be the case. In this paper, we first rely on extracting, from the training set, patterns that are correlated with the remaining useful life of an equipment. These patterns are then matched to testing trajectories to estimate the RUL.

Shapelet extraction and selection

Let T be a training set composed of |T | time series, T 1 , . . . , T |T | . In this set, the time series may have different lengths. The length of the time series T i is denoted by l(T i ). With this notation, a time series T i in the set T can be written T i = t i 1 ,t i 2 , . . . ,t i l(T i ) . For the sake of clarity, we assume here that the time series T i are univariate, i.e.

t i j ∈ R, ∀ 1 ≤ j ≤ l(T i ).
The feature extraction process described in this section can be extended to multivariate time series by applying it to every dimension of the time series.

As we focus in this paper on prognostics (estimation of the remaining useful life of equipment before failure), time series (or trajectories) in the set T represent a run-to-failure health indicator trajectory of an equipment. In this section, we aim at extracting, from the set T of time series, features that are correlated with the remaining length of a time series (time period between the instant when the feature is met and the failure of the equipment).

In the following, we consider features under the form of time series of small length (relatively to the average length of the time series in T ), that will be denoted rul-shapelets.

Definition 1: A rul-shapelet is defined by a tuple f = (S, δ , µ), where S = s 1 , . . . , s l(S) is a time series and δ is a distance threshold. µ represents the average value of the remaining length of a time series that is matched by f (cf. Definition 2).

Definition 2: A rul-shapelet f = (S, δ , µ) is said to match a time series T if there exists a subsequence T of T (of length l(S)) such that the Euclidean distance between S and T is less than or equal to δ . In other words, f matches T if

∃ k ∈ [1, l(T ) -l(S) + 1], s. t. l(S) ∑ j=1 (s j -t k+ j-1 ) 2 ≤ δ . ( 1
)
The match is hence defined at a time instant k. If more than one k satisfies (1), the instant of the match is defined as the one that leads to the minimal distance.

According to Definitions 1 and 2, when a time series T is matched by a rul-shapelet f = (S, δ , µ) at a time instant k, we can estimate that the RUL of T (from time k) is µ, i.e. we estimate the length of the time series T to be equal to k + µ.

In this section, we describe how to extract rul-shapelets from a set of time-series T and select the ones that convey sufficiently relevant information about the RUL.

Shapelet extraction

[Figure 2 about here.]

To extract a rul-shapelet f , we first need to extract a time series of small length (which represents the feature S of the rul-shapelet) from the set T , and we then need to estimate the other three features associated with f : δ and µ. For that purpose, all the subsequences of lengths l 1 , . . . , l N from the set T are first extracted and stored as it can be seen in Figure 2. The lengths l i are parameters chosen by the user. Depending on the number of time series in T and their lengths, the number of subsequences extracted here can be very high. In order to keep a reasonable amount of shapelets, subsequences of length l i , ∀ 1 ≤ i ≤ N can be quantized using the K-means algorithm into a smaller number of subsequences (only the centroids of the obtained clusters are kept). After this step, a set F = { f 1 , . . . , f |F | } of shapelets is obtained. For all these shapelets, only the first feature S is known for the moment. In this case, a rul-shapelet f can also be written f = (S, ?, ?). We explain in the following how to obtain the other two features that completely determine a rul-shapelet.

Shapelet selection

[Figure 3 about here.] Definition 3: Let f = (S, ?, ?) be a rul-shapelet and T a time series from T . The best-match features (BMF) between f and T is defined as the pair (d, ρ), where :

               d = min 1≤ j≤l(T )-l(S)+1 l(S) ∑ i=1 (s i -t i+ j-1 ) 2 ρ = l(T ) -arg min 1≤ j≤l(T )-l(S)+1 l(S) ∑ i=1 (s i -t i+ j-1 ) 2 .
(

) 2
In other words, d is the minimum Euclidean distance between S and a subsequence of T of the same length and ρ is the remaining time before the end of T when this minimum distance d is met.

In the following, d will also be denoted by best-match distance and ρ by best-match RUL. We can, according to Definition 3, compute the BMF between a rul-shapelet f and every time series of the set T . Definition 4: For a rul-shapelet f = (S, ?, ?) and a set T of time series, we define the best-match features list between f and T as the list :

L f = bm f 1 = (d 1 , ρ 1 ), . . . , bm f |T | = (d |T | , ρ |T | ) , (3)
where bm f i , 1 ≤ i ≤ T is the BMF between f and the i th time series of T . This list of pairs is then sorted so that the d i 's are in an increasing order (

d 1 ≤ d 2 ≤ • • • ≤ d T ).
On Figure 3 are shown the two first elements of a BMF list (bm f 1 and bm f 2 ) for a given example shapelet and a given training set. The left picture corresponds to the time series for which the best match distance d 1 is the smallest amongst the other time series. The shapelet is depicted by circles. The instant of the match is 68 here. This leads to a best-match rul ρ 1 of 120 (as the total length of T 1 is 188 here). On the right picture, we can see the second time series (in terms of best match distance). The best-match distance d 2 for this second time series is higher than d 1 but lower than the ones obtained for the other training time series of this example. The instant of the match is 83 for T 2 , which leads to a best-match rul ρ 2 = 112.

We are now interested, for a rul-shapelet f = (S, ?, ?), in finding the parameters δ and µ such that, when f matches a time series T (i.e. the best-match distance between f and T is lower than δ ), we have a high confidence in estimating that the remaining time (from the instant of the match) before the end of the series T is mu in average. For that purpose, we use the best-match features list L f between f and the set of time series T . From this list, we first extract the list of best-match RULs: R = ρ 1 , . . . , ρ |T | . This list R is normalized so that its average value equals 0 and its variance equals 1. We compute the index i (2 ≤ i ≤ |T |) defined by:

i = arg min 2≤ j≤|T | var[ρ 1 , . . . , ρ j ], (4)
where var denotes the statistical variance. This equation also means that we are searching for the i first elements of R of minimum variance. In order to select only the most discriminative rulshapelets, the shapelets yielding to a minimal partial variance (computed as in ( 4)) that is above a threshold τ can be discarded. When the values of R are normalized, the variance of the whole list is 1 (whatever the considered rul-shapelet). As τ gets close to 0, the shapelets become more discriminative.

Then, the parameters of the selected rul-shapelet f are computed as :

δ = d i µ = (ρ 1 + • • • + ρ i )/i (5)
where the values ρ j correspond to the ones before normalization of R.

[Figure 4 about here.]

Example 1:

Figure 4 shows the steps described above to estimate the parameters of a rul-shapelet f . In this example, a training set of 100 training time series is used. Figure 4-(a) shows the values of the best-match RULs between f and the 100 trainingtime series. Note that these values are sorted according to the best match distances as explained above. The right part of Figure 4 shows the partial variances var[ρ 1 , . . . , ρ j ], 2 ≤ j ≤ 100. From these values, the index i = 11 is chosen according to Equation ( 4). The 11 st best-match distance (d 11 in the best-match features list) is selected as the parameter δ and the 11 first values of the left image are selected (shown by triangles) to estimate µ according to Equation [START_REF] Huang | Residual life predictions for ball bearings based on selforganizing map and back propagation neural network methods[END_REF].

After these steps, the set F = { f 1 , . . . , f |T | } of rul-shapelets is filled, i.e. every f i is defined with its three parameters.

Shapelet-based RUL estimation

In this section, we explain how we perform RUL estimation using a set of rul-shapelets F obtained as described in the previous section. The context is the following: we are monitoring the behavior of a component and our aim is to predict when this component is likely to face a failure, based on previous experiences. The monitoring of the component is modeled by a time series U = u 1 , . . . , u l(U) . This time series is obtained as explained in the introduction after data processing and instance formalization (cf. Figure 1). No information is given about the last monitoring instant l(U): it can be at an early stage of the life of the component, a late stage or any stage in between. The RUL estimation process described here aims at predicting the time difference between l(U) and the failure of the testing equipment.

Extracting the rul-shapelets that match U

The first step of the RUL estimation described in this section is to find the rul-shapelets in F that match the testing time series U, and the time instant of the match when applicable. Let f = (S, δ , µ) be a rul-shapelet of F . The best-match distance between f and U is computed (i.e. the minimum Euclidean distance between S and a subsequence of U of length l(S)). If this bestmatch distance is lower than δ , then according to Definition 2, f matches U. The time instant idx f of the match (i.e. the time index of the beginning of the subsequence of U that leads to the minimal distance) is stored together with f . If the best-match distance between f and U is greater than δ , then f is discarded. The same operation is repeated for all the rul-shapelets of F , leading to a set Match(U) = {( f 1 , idx 1 ), . . . , ( f k , idx k )}, where f i , 1 ≤ i ≤ k is a rul-shapelet and idx i the time instant when f i matches U. This set contains all the shapelets that match the testing time series U, and is used to estimate the RUL of U.

For some testing units, it might happen that no rul-shapelets match U. In that case, the set Match(U) is empty and no RUL estimation can be made at that point. To prevent this from happening, the condition of the match (Definition 2) can be relaxed when Match(U) is found to be empty, by adding a coefficient β > 1 to Equation (1) :

∃ k ∈ [1, l(T ) -l(S) + 1], s. t. l(S) ∑ j=1 (s j -t k+ j-1 ) 2 ≤ β × δ . (6)
Relaxing this condition enables to perform a RUL estimation even when no match is found for a test instance U.

RUL estimation

Every rul-shapelet in Match(U) conveys an information about the RUL of U given by the parameter µ of the rul-shapelet. Let f i = (S i , δ i , µ i ) be a rul-shapelet in Match(U) and let idx i be its matching time. According to this rul-shapelet, the estimated total length of U is µ i + idx i .

Taking into account all the information brought by the rul-shapelets of Match(U), the estimated total length lU of U is given by lU

= 1 k k ∑ i=i (idx i + µ i ) (7)
Note that weights can also be inserted in Equation [START_REF] Heimes | Recurrent neural networks for remaining useful life estimation[END_REF] to favor for instance rul-shapelets that match U at late time instants (as late time instants are closer to failures than early ones). We use in this paper a weight that corresponds to the ratio between the instant when the shapelet is matched and the length of the time series where it is matched. Ratios around 1 mean that the shapelet is matched close to the late instants.

Case study

To assess the performance of the approach proposed in this paper on a prognostics case study, we used turbofan engines data available on the NASA Prognostic Data Repository1 [START_REF] Saxena | C-MAPSS Data Set[END_REF]. We first describe the data, then give the performance metrics that we use, and finally present some results.

Turbofan data set

[Table 1 about here.] The data available from [START_REF] Saxena | C-MAPSS Data Set[END_REF] consists of multivariate sensory time series that have been gathered from turbofan engine dynamic simulation process. Run-to-failure data has been obtained using C-MAPSS, a simulation program for turbofan engines. Four different experimental setups have been conducted and have led to four different data sets. The settings of these experiments are given in Table 1. In this paper, we use two of these datasets : the #1 and the #4. The #1 is the easiest one as it considers only one type of fault and one operating mode, whereas dataset #4 is the most difficult one with two types of faults and six operating modes.

Run-to-failure data for the training units are composed of 21 sensor measurements at each cycle. The measurements start at a similar level of degradation, which is considered as healthy and stop when the equipment has reached a level that is considered not sufficient to meet its operating requirements. From these 21 measures, only 7 are kept here as they have been shown to be the most pertinent ones in [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF]. These sensors used here are the number 2,3,4,7,11,12 and 15. This sensory data is corrupted by noise. A third-order polynomial curve is used to smooth the sensor values and keep only the trends of the time series. For illustration purposes, the 100 training time series of dataset #1 given by sensor 2 after smoothing are shown in Fig. 5-(a).

We have considered in this paper two different approaches for shapelet-based RUL estimation : the first one consists in working directly with the sensory data (the process explained in Sections 3 and 4 is done for each sensor) and the second one consists, as explained in the introduction, in computing a 1-dimensional trajectory called health indicator from the sensory data. For dataset #1, these two approaches can be considered and compared, whereas for dataset #4 only the second approach can be considered as the six operating conditions prevents from working directly with the sensory data. To obtain the HI trajectories, we have used linear regression. The linear regression maps the seven sensor values at each cycle into a real number between 0 and 1 (approximately). The coefficients of the linear regression are inferred using training samples. The training samples are chosen as follows : sensory data from the last 5 cycles of every training units are selected and assigned the value 0 (which means bad health state), and the sensory data from the first 5 cycles of the training units whose total life is above 240 cycles are selected and assigned the value 1 (healthy state). The coefficients of the regression are inferred from the training samples with a least square optimization. These coefficients are stored and will be used to construct HI trajectories of the testing units accordingly. For the dataset #4, six models are learned, one for each operating mode. The 100 training HI trajectories obtained for the dataset #1 are shown in Fig. 5-(b).

Note that, as explained in the introduction, our contribution in this paper concerns the retrieval and RUL estimation steps in a similarity-based RUL estimation scheme. Our approach makes use of HI trajectories in order to estimate the RUL of an equipment. In this paper, the HI trajectories are obtained as explained above from sensory data using linear regression. This method was used by Wang et al. [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] in their winning approach in the 2008 PHM challenge. They have then optimized techniques to obtain HI trajectories with more complex tools, as principal component analysis, kernel smoothing and the use of neural networks (cf. [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF] for a complete study). As these approaches rely on heavy parameter settings and given that our contribution concerns the retrieval and RUL estimation steps, we used the approach presented in [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] and focus our results on showing the efficiency of our method given HI trajectories. Better overall performance might be obtained with optimized HI trajectories.

In the online phase, the RUL of the testing units are estimated using the testing sensory data. This data is incomplete: the measurements are stopped at a certain instant and the RUL from this instant needs to be predicted. Sensory data is smoothed and converted into a HI trajectory (if needed, depending on the considered approach) using the learned regression model(s), and the RUL estimation will be made as detailed in Section 4.

[Figure 5 about here.]

Performance metrics

We present and explain the different performance metrics that we use in this paper to assess the performance of the proposed RUL estimation method. These metrics are either taken from [START_REF] Saxena | Metrics for offline evaluation of prognostic performance[END_REF] or were defined by the PHM Challenge community. Let t be a time instant, rt the RUL prediction made at t and r t the real RUL at time t. Let also t E represent the time at which the considered equipment reaches end of life (t E = t + r t ).

Prediction horizon

Let α h ∈ R be a coefficient between 0 and 1 that represent the tolerance to errors. Let t p be the first time instant when the absolute value of the difference between the predicted RUL and the real RUL is less than α h times t E . This time instant t p is hence defined as follows:

t p = min t, s.t. r t -α h × t E ≤ rt ≤ r t + α h × t E . ( 8
)
As it can be seen in Fig. 6-(a), the prediction horizon PH is then defined as PH = t Et p . The higher PH is, the better the prediction approach is. In addition, the prediction horizon should increase with α h .

[Figure 6 about here.]

Rate of acceptable predictions

This metric evaluates the rate of predictions that fall in a cone-shaped region of acceptable prediction errors. This metric has two parameters : a time instant t H from which the predictions are evaluated and α r a coefficient in [0, 1] that represents the tolerance to errors (similarly to α h above). The rate of acceptable prediction RAP is defined as:

RAP = Card({t ∈ {t H , . . . ,t E }, s.t. (1 -α r ) × r t ≤ rt ≤ (1 + α r ) × r t }) t E -t H + 1 , (9)
where Card() represent the cardinal of a set. RAP is a real number ∈ [0, 1]. High values of RAP are associated with good estimation performance. Fig. 6-(b) illustrates how the RAP metric is computed. In this figure, t H is set to 50 and α r to 0.2. Estimations that fall in the cone-shaped region are depicted by circle points. The RAP metric is equal here to 19/43. Note that estimations were done every 5 cycles in this example.

Relative accuracy

The relative accuracy (RA) evaluates the mean absolute percentage errors of RUL predictions for every time instant between t H (set by the user) and t E . It is formally defined as :

RA = 1 -Mean({ |r t -rt | rt , t H ≤ t ≤ t E }). ( 10
)
Values of RA close to 1 are associated with good estimation performance.

Prediction score

This metric was defined by the PHM community for the 2008 PHM challenge. It assigns a score PS to a prediction based on the difference between the prediction and the real RUL as follows.

PS = e -(r t -r t )/10 -1 if r t -rt ≤ 0 e (r t -r t )/13 -1 if r t -rt > 0. ( 11
)
The lower the score PS is, the better the prediction is. Note that late predictions are more penalized than early ones, as late prediction can have harmful impact on the equipment.

In order to compute the first three metrics presented here (PH, RAP and RA), a complete run-to-failure testing HI trajectory is needed. Different subtrajectories are used from this complete trajectory in order to compute these three metrics. Hence, the testing data available in the Turbofan datasets are not compatible with these metrics. Indeed, in the testing data, trajectories are stopped at a given time and the rest of the trajectory is unknown (only the RUL is given for validation purposes). In the following, we present some results in terms of these measures (for the data set #4). To obtain these results, we had to split the training dataset into two parts : a real training part and a validation part in which the complete HI trajectories are used as testing ones. On the contrary, the fourth metric (prediction score) can be used with any kind of testing data. For this metric, we use the testing data available in the Turbofan dataset and we made RUL predictions from the available testing trajectory (after the last measurement of this trajectory).

Experimental settings and results

We present in this section the performance obtained with our proposed RUL estimation technique on two Turbofan datasets : the #1 and the #4. Three parameters influence the performance of the proposed technique : the lengths l 1 , . . . , l N of the extracted shapelets (Section 3.1), the number of centroids n 1 , . . . , n N used to quantize the shapelets (Section 3.1), and the discriminative threshold τ used to select the discriminative shapelets (Section 3.2). The lengths of the extracted shapelets needs to be chosen according to the distribution of lengths of the training units. Short lengths need to be taken in order to be able to cope with short testing units. In the considered application, the length of the training units vary from about 120 of 400. In our experiments, the lengths 10, 20, 30, 40 and 50 have been chosen to extract the shapelets. The total number of extracted shapelets before quantization is here more than 16, 000 for each l i . In the first experiment on dataset #1 (Section 5.3.1), we analyse the performance of our RUL estimation technique with varying values for the parameters n 1 , . . . , n 5 and τ. Conclusions about this are drawn in Section 5.3.1.

Dataset #1

In this first experiment, different values for the parameters n 1 , . . . , n 5 and τ are used. The values for n i , i ∈ 1, . . . , 5 are taken in the set {50, 100, 150, 200, 250, 300}. Note that we took n 1 = • • • = n 5 in this experiment. The parameter τ is taken in the set {0.15, 0.25, 0.35, 0.45, 0.55}. In this experiment, the health indicator trajectories obtained from the sensory data are used. Table 2 sums up the performance obtained in terms of the average prediction score for the 100 testing units. We can observe that the performance generally improves when the value of the parameters n i increases. Too few shapelets lead to bad estimation scores. When n i ≥ 200, the performance is somehow stabilized. We can also observe that when τ is very low (τ < 0.25), the obtained scores are bad for every values of n i . This is mainly because too many shapelets are discarded when τ is low. When τ ≥ 0.25, the performance also somewhat stabilizes. It can be seen from these experiments that when the parameters l i and τ do not take very low values, the performance of the proposed RUL estimation technique is stable, and that the predictions are quite accurate. The best performance is obtained for l i = 300 and τ = 0.45 and is equal to 6.52. These parameters are used in the following.

Figure 7 represents the histograms of the prediction errors for three different approaches : (a) an estimation approach based on [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF], (b) the proposed approach using the 7 dimensional time series and extracting rul-shapelets on each dimension and (c) the proposed approach using the health indicator time series obtained by linear regression. We can observe that the prediction errors are more concentrated around zero for the most right histogram. The prediction error range is [-39, 62] for Fig. 7-(c). Authors of [START_REF] Ramasso | Joint prediction of continuous and discrete states in time-series based on belief functions[END_REF] have also applied their technique to this dataset. A similar error histogram is given, in which the prediction errors are more spread than the ones obtained with our approach. The range of prediction errors is about [-80, 120] in [START_REF] Ramasso | Joint prediction of continuous and discrete states in time-series based on belief functions[END_REF].

[Table 2 about here.]

[Figure 7 about here.]

[Table 3 about here.] Table 3 gives the performance of the proposed approach and compares it to a RUL estimation technique based on [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] in terms of the average score given in Equation [START_REF] Fink | Predicting component reliability and level of degradation with complexvalued neural networks[END_REF]. In the method proposed by Wang et al. [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF], the RUL estimation depends on a parameter that determines the number of nearest neighbors kept. A training trajectory is kept if the Euclidean distance between this training trajectory and a testing one is less than γ times the minimum distance between the testing trajectory and all the training ones. We fixed this γ here to 3 as it gives the best results.

From this table, we can see that using the health indicator approach yields better results than working directly with the sensors. In addition, a better score is obtained using the shapelet-based approach than the one based on [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF].

Dataset #4

We have also applied the shapelet-based RUL estimation approach to the dataset #4, which is the most difficult one. In order to compare our approach with the optimized approach of [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF], in a first experiment we have split the training dataset into two parts, as explained above : 150 instances are kept for training and 99 for validation purposes. This splitting enables us to compute the three metrics PH, RAP and RA and compare the performance with the ones given in [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF] for the same dataset and same configuration. The performance in terms of the metrics PH, RAP and RA of our approach is given in Table 4. In this table, we also report the best results obtained by Wang on a testing set of 99 trajectories (cf. results in the thesis [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF]). The parameters for metric evaluation have been chosen according to [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF], i.e. t H = 80, α h = 0.2 and α r = 0.2. The values reported in this table are the median values obtained for the 99 testing instances. It can be seen that the performance obtained by our approach is very close to the ones obtained in [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF]. Note that, as explained above, the modeling of the HI trajectory in [START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF] has been optimized with many different tools (PCA, Kernel smoothing and neural networks), whereas we have considered here a simple modeling with linear regression.

[Table 4

about here.]

In a last experiment, we have also considered the real testing set of dataset #4. In this testing set, 248 RULs have to be predicted from incomplete run-to-failure data. Table 3 reports the performance of the proposed approach and compares it to a RUL estimation technique based on [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] in terms of the average score given in Equation [START_REF] Fink | Predicting component reliability and level of degradation with complexvalued neural networks[END_REF]. The average score obtained by our approach is significantly better than the one obtained from a RUL estimation based on [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF].

[Table 5 about here.]

Conclusion

We have proposed in this paper a RUL estimation technique based on shapelet extraction. The shapelet extraction process aims at selecting, from a training set of run-to-failure data, patterns (under the form of small time series) that are correlated with the remaining useful life of an equipment. These extracted patterns convey each an information about the RUL of the equipment from the instant when they are found. They are then used in an online step to estimate the RUL of testing units (units for which the remaining useful life is not known). Therefore, the proposed RUL estimation of a testing unit is based on patterns that have been selected for their high correlation with the RUL. This approach is a similarity-based approach : it makes use of available run-tofailure data (in which the health state evolution is captured) to estimate the RUL of equipment. It is hence compatible with techniques aiming at constructing health indicators from sensory data. We assessed the performance of the proposed RUL estimation technique on Turbofan datasets. The RUL estimation performance is shown to be effective and competitive compared to other similarity-based approaches. The proposed method is more adapted to applications where many training instances are available, in order to select relevant shapelets and to be able to provide an estimation of the RUL even at early life stages. Extending this method to limited training sets will focus our attention in a near future. This paper also opens up other possibilities of using shapelets in the field of Prognostics and Health Management. We are planning to extend the use of shapelets for diagnostic purposes by linking the apparition of shapelets with different kind of faults, or for predicting the future evolution of the health status of components. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q qq q q q q q q q q q q q q q q q q q q qq q qq q qqq q qq q q qqq q qqq qq qqqqqqqqq qq qqq qqqqqqqq q q q qq qq q qqq qqqq q q qqq q q q qqqq qqq q 0.2 Time RUL q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Prediction Horizon = 135 Time RUL q q q q q q q q q q q q q q q q q q q

List of Figures

Method

Average Prediction score of Eqn. (11) Proposed approach (health indicator) 37.7097 Estimation based on Wang et al. [START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] 69.2928

Figure 3 :

 3 Figure 3: Shapelet selection process.

Figure 4 :Figure 5 :

 45 Figure 4: (a) Best-match RULs obtained between a rul-shapelet and a set of 100 training time series, ordered according to the best-match distances (Equation (2)) -(b) Partial variances of the best-match RULs of (a). The index of the minimum partial variance is equal to 11 here.

Figure 6 :Figure 7 :

 67 Figure 6: (a) Prediction horizon and (b) Rate of acceptable predictions.

1

  General framework of similarity-based RUL estimation schemes. . . . . . . . . . 17 2 Shapelet extraction process : all subsequences of length l i are extracted and stored. The number of such subsequences is then reduced by the K-means algorithm. This process is repeated for different lengths l i if needed. . . . . . . . . . . . . . . . . 18 3 Shapelet selection process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4 (a) Best-match RULs obtained between a rul-shapelet and a set of 100 training time series, ordered according to the best-match distances (Equation (2)) -(b) Partial variances of the best-match RULs of (a). The index of the minimum partial variance is equal to 11 here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 5 (a) Evolution of the measures of sensor 2 and (b) health indicators obtained by linear regression for the 100 training time series of dataset #1. These curves are obtained after smoothing with a third-order polynomial curve. . . . . . . . . . . 21 6 (a) Prediction horizon and (b) Rate of acceptable predictions. . . . . . . . . . . . 22 7Histograms of the prediction errors (actual RULpredicted RUL) obtained by applying (a) an estimation approach based on[START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF], the proposed approach on (b) the 7-dimensional time series and on (c) the health indicator obtained after linear regression on the original time series. . . . . . . . . . . . . . . . . . . . . . . . 23Figure2: Shapelet extraction process : all subsequences of length l i are extracted and stored. The number of such subsequences is then reduced by the K-means algorithm. This process is repeated for different lengths l i if needed.

					T 1		Extracted shapelets
				0.8
				0.6				. .
								.
				0.4	T1			T2
	0.6 0.8	Training sensory data qqqq qqqq qqqq qqq qqq qqq qqq qqq qqq qq qq qq multi-dim. 0.0 0.2 qq qq qq qq qq qq qq 0 d1= 0.0088	Data processing (optional) 50 ρ1=120	single or dim. Time 100 multi-	Instance formalization 150 200 0.6 0.8	library of instances . qqqq qqqq qqqq qqq qqq qqq qqq qqq qqq qq qq qq qq qq qq qq qq qq qq d2= 0.014	ρ2=112
	0.2 0.4	Testing sensory data Offline stage	(optional) processing Data	T k . . .	formalization 0.2 Instance 0.4	Retrieval step new instance K-means	RUL estimation
	0.0			0.6 Online stage		0.0	. . .
		68	80	Figure 1: General framework of similarity-based RUL estimation schemes. 0.4 100 120 140 160 188 83 100 120 140 160	180	195
					Time			Time
								Library of shapelets
				0.2
				0.0
				0	50	100	150
						Time

Table 1 :

 1 Experimental settings for turbofan data sets.

	Experiment number	#1	#2	#3	#4
	Number of fault modes	1	1	2	2
	Number of operating conditions	1	6	1	6
	Number of training units	100 260 100 249
	Number of testing units	100 259 100 248

Table 2 :

 2 Prediction score of the proposed RUL estimation technique with different values of the parameters l i and τ.

	P P τ	P P # centroids P P P P P P	50	100	150	200	250	300
	0.15	15.35 17.19 19.39 16.97 17.82 12.34
	0.25	12.36 10.30 7.12	7.38	7.36	7.64
	0.35	13.30 10.94 7.19	6.88	6.87	6.75
	0.45	17.49 10.34 7.00	6.80	6.71	6.52
	0.55	17.49 10.90 7.19	6.91	6.79	6.57

Table 3 :

 3 Performance evaluation for the dataset #1 in terms of the average prediction score of the proposed approach and comparison with[START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF].

	Method	Average Prediction score of Eqn (11)
	Proposed approach (7 sensors)	8.07
	Proposed approach (health indicator)	6.52
	Estimation based on Wang et al. [25]	7.91

Table 4 :

 4 Performance evaluation in terms of the metrics PH, RAP and RA of the proposed approach and comparison with[START_REF] Wang | Trajectory similarity based prediction for remaining useful life estimation[END_REF] for dataset #4.

	Method	Prediction horizon Rate of acceptable prediction Relative accuracy
	Proposed approach (health indicator)	176	0.3750	0.6819
	Results from Wang et al. [33]	177	0.40	0.7124

Table 5 :

 5 Performance evaluation in terms of the average prediction score of the proposed approach and comparison with[START_REF] Wang | A similarity-based prognostics approach for remaining useful life estimation of engineered systems[END_REF] for dataset #4.

http://ti.arc.nasa.gov/tech/dash/pcoe/prognostic-data-repository/

Acknowledgment

The authors would like to acknowledge the European Regional Development Fund (FEDER) who funded this work as a part of a project named ALTIDE: Aid to Lifecycle Traceability for Intelligently Developed Equipment.

List of Tables