Jaouher Ben Ali

Lotfi Saidi

Aymen Mouelhi

Brigitte Chebel-Morello

Farhat Fnaiech

Linear Feature Selection and Classification Using PNN and SFAM Neural Networks for a Nearly Online Diagnosis of Bearing Naturally Progressing Degradations

Keywords: Empirical Mode Decomposition (EMD), Linear Discriminant Analysis (LDA), Principal Component Analysis (PCA), Probabilistic Neural Network (PNN), Scatter matrix, Simplified Fuzzy Adaptive resonance theory Map (SFAM)

In this work, an effort is made to characterize seven bearing states depending on the energy entropy of Intrinsic Mode Functions (IMFs) resulted from the Empirical Modes Decomposition (EMD). Three run-to-failure bearing vibration signals representing different defects either degraded or different failing components (roller, inner race and outer race) with healthy state lead to seven bearing states under study. Principal Component Analysis (PCA) and Linear Discriminant Analysis (LDA) are used for feature reduction. Then, six classification scenarios are processed via a Probabilistic Neural Network (PNN) and a Simplified Fuzzy Adaptive resonance theory Map (SFAM) neural network. In other words, the three extracted feature data bases (EMD, PCA and LDA features) are processed firstly with SFAM and secondly with a combination of PNN-SFAM. The computation of classification accuracy and scattering criterion for each scenario shows that the EMD-LDA-PNN-SFAM combination is the suitable strategy for online bearing fault diagnosis. The proposed methodology reveals better generalization capability compared to previous works and it's validated by an online bearing fault diagnosis. The proposed strategy can be applied for the decision making of several assets.

Introduction

Anomalous operating condition and fault detection of industrial assets is the core of process engineers in industry. Effective detection ensures the product quality improvement and reduces the cost and the time of repair. Consequently, industrial assets safety and reliability are very important to ensure the continuity of production, to improve process operation, to increase plant throughput, to reduce process downtime and to comply with increasingly stringent environmental rules and safety regulations [START_REF] Lau | Fault diagnosis of Tennessee Eastman process with multi-scale PCA and ANFIS[END_REF].

Rolling Element Bearing (REB) is the most used solution for industrial machinery guiding in rotation.

This critical component requires running with high reliability to decrease fault occurrences and fatal breakdowns of machineries [START_REF] Lei | Fault diagnosis of rotating machinery based on multiple ANFIS combination with Gas[END_REF][START_REF] Konar | Bearing fault detection of induction motor using wavelet and Support Vector Machines (SVMs)[END_REF]. It is well known that bearing defects are one of the most common fault sources in induction machineries (about 40-50%) (Widodo and Yang, 2007;[START_REF] Frosini | Stator current and motor efficiency as indicators for different types of bearing faults in induction motors[END_REF]. Failure surveys by the electric power research institute indicate that bearing-related faults are about 40% among the most frequent faults in induction motors [START_REF] Bellini | Advances in diagnostic techniques for induction machines[END_REF]. Therefore, accurate strategies for REB defect detection are a trivial task nowadays in the industry.

Over the past few years, various signal processing methods have been proposed to detect and diagnose bearing defects. These methods may be roughly classified into vibration and acoustic measurements, temperature measurements and wear debris analysis [START_REF] Nandi | Condition monitoring and fault diagnosis of electrical motors-A review[END_REF]. The mechanical bearing vibration signals remain the most immediate, simplest, and richest source of information for understanding phenomena related to bearing defects [START_REF] Seungdeog | Implementation of a fault-diagnosis algorithm for induction machines based on advanced digital-signal-processing techniques[END_REF][START_REF] Žvokelj | Non-linear multivariate and multiscale monitoring and signal denoising strategy using Kernel Principal Component Analysis combined with Ensemble Empirical Mode Decomposition method[END_REF].

Certainly, vibration analysis is a powerful source to extract interesting information and to apply appropriate techniques. However, REB accelerations are considered as non-stationary and nonlinear [START_REF] Randall | Rolling element bearing diagnostics-A review[END_REF]. Besides, noises present a serious trouble in the study of this type of signals [START_REF] Zhang | Rolling element bearing fault diagnosis based on the combination of genetic algorithms and fast kurtogram[END_REF]. Moreover, the relatively weak bearing signals are always affected by quite stronger signals (gear meshing, unbalance, bars, …) [START_REF] Randall | Rolling element bearing diagnostics-A review[END_REF].

Several efforts based on vibration signals have been recently

proposed in the open literature to deal with the diagnosis of bearing faults. When a fault occurs on one bearing surface, the acceleration is characterized by the presence of periodic repetitive sharp peaks and further modulated by a number of harmonic frequencies [START_REF] Zhang | Rolling element bearing fault diagnosis based on the combination of genetic algorithms and fast kurtogram[END_REF]. However, interesting information is often contained in the periodicity of impacts, rather than in the rest of the signal frequency content. Thereby, vibration signals are non-stationary and always masked by machine noise [START_REF] Randall | Rolling element bearing diagnostics-A review[END_REF]. For this reason, conventional vibration analysis techniques, which were successfully employed for stationary signals such as fast Fourier transform, haven't provided good and significant results based on REB vibration signals [START_REF] Braun | Decomposition of non-stationary signals into varying times scales: Some aspects of the EMD and HVD methods[END_REF].

In the literature, non-stationary or transient signals are often decomposed into sub-bands where stationary and linear characteristics could be obtained and fault patterns are easily extracted [START_REF] Priestley | Nonlinear and Nonstationary Time Series Analysis[END_REF]. Consequently, Wavelet Transforms (WT) [START_REF] Wang | Detection of weak transient signals based on wavelet packet transform and manifold learning for rolling element bearing fault diagnosis[END_REF], Empirical Mode Decomposition (EMD) [START_REF] Liu | Bearing faults diagnostics based on hybrid LS-SVM and EMD method[END_REF], Local Characteristic-scale Decomposition (LCD) (Zheng et al., 2013a), Ensemble Empirical Mode Decomposition [START_REF] Wang | Feature extraction of rolling bearing's early weak fault basedn on EEMD and tunable Q-factor wavelet transform[END_REF], and Generalized Empirical Mode Decomposition (Zheng et al., 2013b) were widely used for the processing of bearing vibration signals in the presence of nonlinear and non-stationary data. By analyzing some series of decomposed REB signals based on WT and EMD, it was demonstrated that the energy entropy changes in different frequency bands when a bearing fault occurs. Besides, EMD is more highlighted than WT in terms of mean fault characterization [START_REF] Tavakkoli | A ball bearing fault diagnosis method based on wavelet and EMD energy entropy mean[END_REF].

After decomposing bearing accelerations into some linear and stationary signals, it is recommended to extract the most useful features. Bearing fault diagnosis is not an easy task, it is essentially a problem of pattern recognition. The most effective features and accurate classifiers are needed to obtain higher diagnostic accuracy [START_REF] Li | Supervised locally linear embedding projection for machinery fault diagnosis[END_REF].

Several previous works based on feature extraction and classification techniques have been proposed for bearing fault diagnosis. Yang et al have used the discrete WT to decompose vibration and acoustic signals into different frequency levels. Then, some statistical measures were extracted from the first-four levels. Also, an Artificial Neural Networks (ANNs) combined with a Support Vector Machine (SVM) were used to classify the faults of small reciprocating compressor used in refrigerators [START_REF] Yang | Intelligent fault diagnosis of rolling element bearing based on SVMs and fractal dimension[END_REF]. In [START_REF] Lin | Fault diagnosis of natural gas compressor based on EEMD and Hilbert marginal spectrum[END_REF], the EEMD was applied to decompose vibration signals into different Intrinsic Mode Functions (IMFs) and then the Hilbert spectrum was computed to extract bearing fault features. In [START_REF] Bin | Early fault diagnosis of rotating machinery based on wavelet packets-empirical mode decomposition feature extraction and neural network[END_REF], authors have demonstrated that the WT generally has some shortcomings and difficulties in selecting suitable basic function because there are no standards or general selection rules for different tasks. However, the EMD is an automated method that extracts easily failure patterns thanks to the use of a wide scale of frequency. In [START_REF] Saidi | Bi-spectrum based-EMD applied to the non-stationary vibration signals for bearing faults diagnosis[END_REF], a combination of EMD and bi-spectrum was presented. The computation of the bi-spectrum was based on the first IMF which contains more information. The proposed combination was promising but it needs always an expert intervention for the analyses of the bi-spectrum plot. Unfortunately, all these methods presented an offline bearing fault detection. Generally, some holes were created in the tested REB and subsequently the fault detection efficiency of the proposed method was given. In this case, the bearing vibration signals increase and the fault detection task become effortless. In an industrial environment, REBs are in continuous rotation where fault characteristics are always submerged and hidden by noises. So, it is recommended to investigate bearing defects online by using run-to-failure histories.

The online bearing monitoring is very important to avoid unexpected breakdowns where the breakdown of a single machine can halt the production of the entire process chain. However, investigations related to naturally induced and naturally progressed defects of REBs are relatively scarce. In fact, online detection of REB degradation for naturally progressing defect in damage stages is not available in the literature [START_REF] Shakya | A novel methodology for online detection of bearing health status for naturally progressing defect[END_REF]. Recently, some approaches have been proposed for the online bearing health monitoring. In [START_REF] Pan | Machine health condition prediction via online dynamic fuzzy neural networks[END_REF], the authors have proposed a new approach based on dynamic fuzzy neural networks where only consequent parameters were updated online. In [START_REF] Shakya | A novel methodology for online detection of bearing health status for naturally progressing defect[END_REF], various damage stages for naturally progressing REB defect were given online. For this, time domain, frequency domain and time-frequency domain were together under study. The Mahalanobis distance was used as a trend parameter for bearing fault diagnosis. In [START_REF] Loutas | The combined use of vibration, acoustic emission and oil debris on-line monitoring towards a more effective condition monitoring of rotating machine[END_REF], the combination of vibration, acoustic emission and oil debris measures is used in an online condition monitoring of bearings and gearboxes. Several techniques were combined such as fusion technique, Principal Component Analysis (PCA) and Independent Component Analysis (ICA) to identify different damage modes. The computational time was not discussed in this work and with respect the definition of failure thresholds was arbitrary. In (Ben Ali et al., 2015b), the combination of traditional statistical features and EMD energy entropy were jointed with the back propagation neural network for an online bearing fault diagnosis. The classification results were relatively low (93%) that is why the authors have proposed a health index for a more reliable online monitoring.

In this paper, an attempt is made to diagnose the state of bearings using a new nearly online method. The proposed method is based on the EMD feature extraction, linear feature reduction and feature classification using neural networks. The validation of the proposed method is based on bearing run-tofailure histories. This method is able to detect online perfectly the state of the bearing. The remainder of this paper is organized as follows: Section 2 is devoted to describing the used techniques in this work. This section describes the different steps of the EMD method. Also, a brief description of linear discriminant analysis and principal component analysis algorithms is given. Besides, this section details the different steps to implement the used neural networks in this work. Section 3 is dedicated to the experimental results by analyzing the proposed feature extraction, feature reduction and feature classification techniques. Section 4 gives a numerical analysis of the most effective features based on a mathematical criterion. Also, a good discussion and analysis of the experimental results by comparing the performances of the proposed method with some previous works is provided in this section. Finally, conclusions and prospects of this work are given in Section 5.

Methods

Empirical Mode Decomposition (EMD)

EMD algorithm is developed from the simple assumption that any signal consists of different simple intrinsic modes of oscillations called Intrinsic Mode Functions (IMFs). This adaptive decomposition method is especially applicable to the analysis of nonlinear and non-stationary signals [START_REF] Yu | A roller bearing fault diagnosis method based on EMD energy entropy and ANN[END_REF].

By the definition, any signal x(t) can be decomposed as detailed in Table 1 [START_REF] Huang | An intelligent fault diagnosis method of high voltage circuit breaker based on improved EMD energy entropy and multi-class support vector machine[END_REF].

The step (G) in Table 1 presents a stoppage criterion for the EMD which is one of the weaknesses of this algorithm [START_REF] Lei | A review on empirical mode decomposition in fault diagnosis of rotating machinery[END_REF]. Huang et al proposed a robust way to access the IMF component of a signal as "Sifting process". The sifting process is stopped by limiting the size of the standard deviation (S D ), computed from two consecutive results as [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis[END_REF]:

  0.2 0.3 2 T k -1 k D 2 t=0 k -1 h (t)-h (t) S = h (t)   (1)
Numerous publications of the EMD for bearing fault diagnosis have shown several advantages and very encouraging results. Nowadays, EMD is one of the most powerful advanced signal processing techniques and it has been extensively studied and widely applied in fault diagnosis of rotating machinery. To achieve better results, it is recommended to combine the EMD with other signal processing techniques or artificial intelligent techniques for an effective strategy for bearing diagnosis [START_REF] Lei | A review on empirical mode decomposition in fault diagnosis of rotating machinery[END_REF].

Brief description of principal component analysis

In literature, Principal Component Analysis (PCA) is the most used technique for feature reduction which means transforming the original features into a lower dimensional space. PCA, as a linear technique, is a quantitatively rigorous method for achieving data dimensionality reduction of the extracted features. The method generates a new set of variables, called Principal Components (PCs), which maximize the variance of the projected vectors. Each PC is a linear combination of the original variables. All PCs are orthogonal to each other, so there is no redundant information. The PCs are a whole of an orthogonal basis for the data space. Thus, the first PC consists of the high variability, the second PC consists of the next highest variability and so on for other directions. The first few components are kept and others with less variability are discarded (Widodo et al., 2007). In this work, the PCA algorithm will be used to eliminate the non significant features in the recorded vibration signals. The step by step procedure of PCA is provided below. Let's say we have N observations of Mdimensional data.

 Step 1: To perform PCA we should first compute the covariance matrix ΣR M  M ; Σ=E{(xx )(x-

x ) T }, where x is the given signal, x is its mean vector, and E{.} is an average over an ensemble of N samples.



Step 2: Find matrix Λ of eigenvectors and diagonal elements of matrix D as eigenvalues of covariance matrix Σ as given by, Λ -1 ΣΛ= D.

 Step 3: Sort the eigenvectors of PCs in the decreasing order of importance of eigenvalues in D.



Step 4: Project the data into the directions of sorted eigenvectors by taking the dot product between the given data and eigenvectors.



Step 5: Select the first few principal components depending on the containment of a given percentage of variability (like 95% or 98% depending on the problem).

Brief look into linear discriminant analysis

Linear Discriminant Analysis (LDA) is a mathematical transformation based on a linear discriminant combination. As the PCA method, LDA is used to reduce features by transforming the original features into a lower dimensional space. Typically, LDA is used to analyze the original dataset and to produce K-1 independent components where K is the number of classes. Generally, LDA dataset preprocessing is more suitable for classification tasks than PCA [START_REF] Mouelhi | Automatic image segmentation of nuclear stained breast tissue sections using color active contour model and an improved watershed method[END_REF]: PCA tends to minimize feature variances of the same class although LDA tends to minimize feature variances of the same class and to maximize feature variances of different classes. The aim of LDA is to create new variables by maximizing discrimination between all classes. New variables maximize the differences between the predefined classes and minimize differences in the same class [START_REF] Subasi | EEG signal classification using PCA, ICA, LDA and support vector machines[END_REF].

Probabilistic neural network

The first formulation of the Probabilistic Neural Network (PNN) was by [START_REF] Specht | Probabilistic Neural Networks[END_REF]. The different steps of PNN algorithm are well described in [START_REF] Specht | Probabilistic Neural Networks[END_REF]. Brief, PNN is based on the computing of Bayes statistics for two categories (A and B). Consequently, for each p dimensional vector x, the state θ could either be θ A or θ B . The state decision is based on the Bayes rule as follow:

 If     A A A B B B h l f x h l f x  then   A dx   : State A.  Contrarily, if     A A A B B B h l f x h l f x  then   B dx   : State B.
where f A (x) and f B (x) are respectively the probability density functions of the category A and B, l A and l B are the loss functions associated respectively with the decision d(x)=θ B and d(x)=θ A , h A is the prior probability of category A pattern occurrences and vice-versa for h B . We should note that the boundary of density probabilities A and B is given by:

    BB AB AA hl f x f x hl  (2)
As shown in Fig. 1, PNN consists of four layers: the input layer, the pattern layer, the summation layer and the decision layer. The p nodes in the input layer consist of the p input vector components. The pattern layer is designed to receive all components of the input vector x. For each patter node i, the product Z i =xW i is computed. To deal with the nonlinear feature distributions for any classification task, a nonlinear function, where x and i W should be normalized to the unit length, is introduced as follow:

    T ii i 2 W X W X Z1 exp exp 22            (3)
After computing the nonlinear function for all nodes, this function will be transferred to the summation layer. Afterwards, summation units sum pattern unit inputs which correspond to the desired category.

The decision layer consists of two neurons that produce binary outputs taking into account the variable weight C k with the following mathematical expression:

k k k k k k B B A k A A B h l n C h l n  (4)
where n A and n B is the number of training patterns from the category A and B respectively. More details are given by [START_REF] Specht | Probabilistic Neural Networks[END_REF].

Simplified Fuzzy Adaptive Resonance Theory Map neural network (SFAM)

The Simplified Fuzzy Adaptive Resonance Theory map (SFAM) neural network is a fast, incremental, supervised learning system for analog inputs [START_REF] Ali | Identification of early stage Alzheimer's disease using SFAM neural network[END_REF]. In 1993, Kasuba proposed the SFAM network which is a simplification of Fuzzy ARTMAP [START_REF] Venkatesan | Classification of Renal Failure Using Simplified Fuzzy Adaptive Resonance Theory Map[END_REF]. This neural network is a step ahead of Fuzzy ARTMAP in reducing the computation and architectural redundancy.

As illustrates Fig. 2, the SFAM is formed of three layers; the input layer F 1 , the competitive layer F 2 and the output layer F ab .

The activity patterns of layers F 1 and F 2 are fully interconnected. Each neuron is initially connected to every neuron on the other layer. All connections are weighted by synaptic weights lying between 0 and 1.

Each neuron of the F 2 layer represents one category formed by the network and it is characterized by its weight vector W j (j is the index of the neuron). The weight vector's size is equal to the dimension M of layer F 1 . In the training process, if the weights of a neuron are not modified, this neuron is called uncommitted. Inversely, once a neuron's weights have been modified, this neuron is named committed.

SFAM is based on a form of normalization called complement coding. The operation consists on taking the input vector and concatenating it with its complement. The resulting vector should be presented to the layer F 1 . Therefore, the dimension M of layer F 1 is equal to the double of the input vector's dimension (M=2L').

The F ab layer is connected with the competitive layer F 2 by the W ab weights to locate the appropriate class t of the input vector in the test mode.

The SFAM has a self-determining architecture: the training vectors are presented to the network, the output category nodes are grown and assigned to a particular class. If there is a new category of input training example, a new output category node is grown and its weights are set equal to the elements of the complement-coded input vector I.

For further applied training vectors, the subset hood of the weight vectors to the complement coded input is computed using the activation function:

  j j j IW TI W     (5)
where :

 ^ is the fuzzy AND operator;

 α is a very small value (typically 5 10  );

 j is the number of the current output category node.

The node with the highest T j (I), called the "winner", is initially chosen (G) to represent the input vector.

To know if associating W G with I is a good match, the matching function MF(I,W G ) between the two is calculated as :

  G G IW MF I, W I   (6)
The matching function should be compared to the vigilance parameter ρ.

 If MF(I,W G ) ≤ ρ, the winner G is disabled and the next most activated node is selected.

 If MF(I,W G )>ρ, a category agreement test between the winner and the input example is performed.

If the categories disagree, the vigilance ρ is increased by a small amount α and the winner is disabled. If they agree the SFAM is considered in a state of resonance and the weight vector W G is updated using the following formula:

    new old old G G G W 1 W I W      ( 7
)
where β is the learning rate introduced to select the mode of learning as follow:

 If β=1, the fast learning mode is selected with just one iteration of the training data.

 Whereas if 0<β<1, slow learning is selected and it requires a number of iterations for the training.

Finally, when classifying, each test vector is applied in turn and its predicted class is that of the node with the highest activation. More details are given in [START_REF] Ali | Identification of early stage Alzheimer's disease using SFAM neural network[END_REF]Ben Ali et al., 2015a).

3 Experimental results

Experimental setup and data recording

The o'clock is used. We should also note that the last file in the first test is uninformative because the record was made in the stoppage phase of the machine.

Table 2 focuses on some important points. The REB life is a highly nonlinear phenomenon. The 12 used bearings in the experimental setup were tested in the same conditions of speed and load. Despite that, only four bearings were broken and one bearing was slightly degraded (bearing 3 of testing 3).

Moreover, the failure time and the failure nature of the three broken bearings were not the same. All these confirm that online bearing fault assessment is a very hard task.

Feature extraction

In this section, the three broken bearings in the IMS database are explored. These bearings are bearing 4 and 3 of testing 1 ending respectively with a failure roller and a failure inner race and bearing 1 of testing 2 ending with a failure outer race. The run-to-failure vibration signals of these bearings are shown in the ii. Decompose just the first part of the extracted signal (with 0.1 second as length) into some IMFs. If the EMD algorithm still running after 10s, stop it and consider the second part of the signal (from 0.1 second to 0.2 second) and so on.

iii. Compute the total energy E i of the first 12 IMFs as (where i is the IMF number) :

2048 2 i ij j=1 E = c  (8)
iv. Compute the total energy of all IMFs:

12 i i=1 E = E  (9)
v. Construct the feature vector as:

  1 2 12 en en en 1 en 2 en 12 E E E H , , ,…, = H ,H IMF , H IMF , ,H IMF E E E       (10)
where H en is EMD energy entropy of the whole original signal and H en IMF G is EMD energy entropy of the IMF number G. The EMD energy entropy is designated as [START_REF] Hong | An adaptive method for health trend prediction of rotating bearings[END_REF]:

  n en i i i=1 H = -p log p  (11)
Where p i =E i /E is the percent of the energy of c i (t) in the whole signal energy.

By considering just the first part of the vibration signal, the number of the treated points by the EMD algorithm is reduced from 20480 to 2048 and the computational time is reduced impressively. Besides, if the computational time of the EMD algorithm exceeds 10s, we consider the next part of the vibration signals. This is a very smart step to ensure a quick convergence of the proposed feature extraction strategy. We note that for all the considered vibration signals, just the first three parts of the signal were considered and thereby the maximum computational was time less than 30s. Note that the reported time is on a typical PC configured with CPU i3 2.3 GHz, 4GB RAM, and running Windows 7.

The energy entropy of the original signal and the first twelve IMFs of the bearing 1 of testing with an outer race defect (shown in Fig. 4.c) are shown in Fig. 5. In this plot, feature number zero is the energy entropy of the original signal and the other features (from 1 to 12) are the energy entropies of IMFs (from 1 to 12 respectively). Through Fig. 5 it is obvious that the energy entropy is near zero for the last IMFs.

In [START_REF] Yu | A roller bearing fault diagnosis method based on EMD energy entropy and ANN[END_REF], Yu et al. demonstrated that the EMD energy entropy changes with the energy variation of bearing fault vibration signals when the REB operates with different faults. Thereby, we propose to use the energy entropy of the first 12 IMFs as features to characterize REB defects. We select just the first 12 IMFs because it is the minimum IMF number found using the EMD. Thus, we are sure that any REB vibration signals produce more than 12 IMFs.

Feature reduction and classification

To show the repartition of the extracted features, 50 vectors were selected arbitrarily from three classes (H, DR, FR). Fig. 6 shows the energy reparation of the first eight IMF energy entropies and the original signal energy entropy. This plot clearly demonstrates that the first IMFs consist of the highest discrimination. For more clarification, Fig. 7 and Fig. 8 show two bi-dimensional plots of the highest features' discrimination (shown in Fig. 6.a) and lowest features' discrimination (shown Fig. 6.g) respectively. Fig. 7 and Fig. 8 reveal that the overlapping of features increases for IMFs higher order than 4.

As shown previously, plotting EMD energy entropy features indicates the necessity of the preprocessing step them to generate separable features more suitable for the classification task. Disorder distribution of original features tends to decrease the performance of the used classifier. So, the combination of these features directly with a classifier is not recommended.

To avoid this drawback, Principal Component Analysis (PCA) and Linear Discriminant Analysis (LDA) are proposed to extract and to reduce the feature dimensionality based on the eigenvalues of the covariance matrix. Firstly, because the overlap between normal-degradation zones and degradationfaulting zones, the number of records of the broken three bearings is reduced from 5394 to 3911. This reduction is made with a very great caution to ensure the reliability of measures. Then, the new data set is divided into two sub-sets: one for training (containing 60% of the samples), and the other for test (containing 40% of the samples).

By using the PCA for the training dataset, the first eight principal components have been selected to replace the original feature vector (which contains 13 components). As shown in Fig. 9, the number of features is reduced from 13 to 8 where principal components 9, 10, 11, 12 and 13 are discarded.

By using LDA technique, original features of the training dataset are reduced to automatically to K-1

where K is the number of classes. In this paper, we have used seven REB states. Thereby, the original features are reduced from 13 to 6 by applying the LDA. Fig. 10.b and Fig. 10.c shows the distribution of the first three features extracted using PCA and LDA respectively. To evaluate the importance of PCA and LDA techniques, in Fig. 10 we have used the same features shown in Fig. 6. To have a good interpretation, Fig. 10.a shows the optimal and the desired feature distribution where classes are perfectly separated.

Compared to previous works, this paper is based on a large data set which contains 3911 records extracted from run-to-failure REB vibration signals. For each record, the EMD energy entropy is used to extract 13 features. Then, the PCA and LDA are employed to reduce features to 8 and to 6 respectively.

To obtain a good analysis of the proposed method, six scenarios are proposed: In the first three scenarios, the SFAM neural network is used to classify seven REB states (H, DIR, FIR, DOR, FOR, DR, FR). However, in the last three scenarios, the PNN is introduced to classify REB states in two classes: healthy and not healthy. If the PNN decision shows that the processed state is unhealthy, then the SFAM is used to detect the type of failure (DIR, FIR, DOR, FOR, DR, FR).

For each scenario, the training-bearing dataset is divided into ten equal sized sub-datasets (ten folds).

Then, ten iterations of training and validation are performed. One iteration presents nine folds for the training and one fold for the validation. Also, ten folds cross-validation experiments are done to validate the classification accuracy of each scenario.

After the validation phase, we have projected the testing dataset using PCA and LDA techniques. Then, the classification accuracy (CA) is computed for each scenario. The CA is the ratio between the total number of correctly classified test samples to the total number of test samples and it can be expressed as follow:

(number of correctly classified samples ) CA[%] = 100 (total number of samples in testing dataset ) (12)

Table 3 shows the CA for each scenario where the introduction of PNN is highlighted compared to the direct use of the SFAM neural network. The PNN neural network is able to classify the whole dataset in two classes (healthy and not healthy) with a CA about 100%. Compared to the use of SFAM only, the proposed combination of PNN-SFAM is very promising for classification tasks. Besides, LDA is better than PCA in terms of feature discrimination. This result is validated when comparing the CA of scenarios (2,3) and (5,6). Better experimental results are insured by the combination of LDA, PNN and SFAM as illustrated by the 6 th scenario summarized in Fig. 11.

Nearly Online bearing fault diagnosis

As mentioned above, the online diagnosis is a very hard task for scientists. With the recent advances in modern technology, industry and scientist researchers are progressing toward providing enhanced strategies to detect online REB anomalies.

In this section, we will present the application of the proposed EMD-LDA-PNN-SFAM combination for an online diagnosis task. We recall that the three used bearings (see Fig. To improve the quality of the proposed diagnosis approach, we propose a smart step to minimize the false alarm rate: It is a sliding window that filters external tests which were not used in the training mode and it can be summarized as follow:

-Select a sliding window (sliding step (ss=1), window length (wl=16)) to browse the nearly online classification results.

-Set k DIR =k FIR =k DOR =k FOR =k DR =k FR =0, where k i is the occurrence of the defect i along the sliding window.

-For each point of the window, if the classification result is i then k i =k i +1;

-

Discussion

Feature performance measures

The ultimate goal in pattern classification is to find the best feature dataset. It is a common practice to get a feature dataset with high performance measures. To achieve this goal, a mathematical treatment is presented in this section. From various measures known in literature, we discuss the interclass/intraclass distance well known as scatter matrix. This measure is applied to the multi-class case and it is well known as the J criterion. In this paper, the J criterion is used to evaluate original extracted features, PCA features and LDA features. As clarified previously, seven REB states are defined in this work. For each class, 50 extracts are used to be evaluated. Also, these features are projected with PCA and LDA to be evaluated. In this way seven classes have been used (K = 7) and 50 extracts are evaluated (N = 50). The procedure of the J criterion is summarized as follows [START_REF] Ali | Importance of the fourth and fifth intrinsic mode functions for bearing fault diagnosis[END_REF]:

i. Compute the average feature vectors of the k th class:

N k,n k n=1 1 µ = x N  (13)
where x k,n the n th extract of the k th class.

ii. Compute the total average feature vectors of all classes:

1 1 K c k k μ K     (14)
iii. Compute the intra-class variance matrix of average dispersion coefficients:

   , , intra 11 1 KN t k n k n k k kn S x µ x µ KN      (15)
iv. Compute the inter-class variance of average dispersion between different classes:

   inter 1 1 K t k c k c k S µ µ µ µ K      (16) v.
Finally, the J degree of the average scattering within classes is given by:

  1 intra inter J trace S .S   (17)
The computation of the J criterion based on original features, PCA features and LDA features is given respectively in Table 4, Table 5 andTable 6. These three tables improve the importance of the LDA feature reduction technique. First, the use of LDA reduces the 13 original features to 6 features which decrease the training computational time of PNN and SFAM. In addition, the first component LD 1 shows the greatest scattering criterion compared to all features. Besides, the scattering criterion of all LDA features is more discriminant compared to original and PCA features. In conclusion, the comparison of the six proposed scenarios confirms the important of EMD-LDA-PNN-SFAM combination.

Comparison with some previous works

Classification accuracy is much more improved for all bearing diagnosis conditions by the use of methodologies based on feature extraction, feature selection (or reduction) and classification tools. The classification accuracy results are all greater than 92%. To get a good classification, the majority of scientific researches in the previous works have used four bearing states (H, FIR, FOR and FR). The superiority of this work is highlighted in the introduction of the default severity which was not held before. Severity consists of three types: normal, degraded and faulty. Thereby, seven bearing states have been used in our work. These added modifications allow the real time monitoring degradation. Table 7 summarizes previous works on automated identification of bearing faults.

The importance of our work over others is the use of real data which results from real time vibrations describing the dynamic response of defective REBs. However, previous works are based on simulation defects by creating holes with different diameters and depths. The proposed approach produces good classification accuracy despite the important number of classes and records. Moreover, our proposed methodology is tested online however the majority of the previous works do not satisfy this aim and scope.

Compared to the work proposed in (Ben Ali et al., 2015b), this paper presents a strategic improvement.

As showed in Table 7, the number of features used in (Ben Ali et al., 2015b) was 18 and the used classifier was the back propagation neural network (which is perhaps, with respect, even a bit outdated).

The structure of the used neural network was important (four layers and 60 nodes) which increase remarkably the computational time. Besides, the use of statistical features certainly degrades the quality classification: it is well known that they are redundant and consequently we don't need all of them as they provide the same information [START_REF] Tandon | A review of vibration and acoustic measurement methods for the detection of defects in rolling element bearings[END_REF]. Also, statistical features present strange behavior: their values come down to the level of an undamaged bearing when the damage is well advanced [START_REF] Tandon | A review of vibration and acoustic measurement methods for the detection of defects in rolling element bearings[END_REF]. In this paper, we have applied the LDA technique to select best features derived from the EMD method and guarantee an optimal separation between features of different classes. Moreover, we have proposed a new combination of PNN and SFAM neural networks. This combination is very promising and it is through her that we obtained 97.82% as classification accuracy. The combination of PNN and SFAM don't require structure initialization, however, the back propagation neural network needs to be compiled many times to confirm a structure which can be not optimal. In the literature there are some empirical recommendations to define the structure of the back propagation neural network, however, no mathematical justification was given [START_REF] Capocchi | Wound-Rotor Induction Generator Short-Circuit Fault Classification Using a New Neural Network Based on Digital Data[END_REF].

Although that classification accuracy, of the previous works given in Table 7, was very satisfactory, these works were based on synthetic bearing defects which were always introduced by the user (holes with different diameters, depths and numbers). In reality, this is not the case of a bearing defect in industrial environment. Many factors can cause deterioration of bearings such as:

-Contamination and corrosion, -Lack of lubrication causing heating and abrasion, -Defect of bearing's mounting, by improperly forcing the bearing onto the shaft, -Misalignment defects.

REB degradation is a highly nonlinear phenomenon; using the same type of bearing in the same experimental conditions never produces the same kind of failure in terms of time, type, and severity. This property is confirmed by Intelligent Maintenance Systems (IMS) bearing run-to-failure data set.

Relatively, only few papers discussed the bearing diagnosis based on run-to-failure vibration signals. This is due to the difficulties in processing these signals. The major problem is the noise; it is very difficult to detect degradation when it is smaller than the noise measurements. This form of degradation is often known as "naturally progressing degradations". Even though many researchers have performed fault bearing detection, they do not perform the diagnosis and identification of naturally progressing degradations.

The twelve used bearings in this work confirm that bearing degradation is a highly nonlinear phenomenon. The failure time and the failure nature of the three broken bearings were not the same. All these confirm that bearing fault assessment is a very hard task. In conclusion, the proposed method is very adequate in terms of fault accuracy detection.

Statistical validation

The confusion matrix is a useful tool for analyzing how well the classifier can recognize tuples of different groups, which contains information about real and predicted classifications done by a classification system (Fawett, 2001).

From the confusion matrix of the proposed method given in Table 8 and presented in Fig. 16 by a compact box plot, it can be seen that the EMD-LDA-PNN-SFAM combination recognize the type and the severity of defects effectively.

For more statistical significance and better generalization, some series of experiments should be conducted to support the author's claims. In fault detection problems, the performance of a detection algorithm usually depends on the trade-off between robustness and sensitivity. The sensitivity and robustness of the proposed method need to be explored by running a series of experiments. A receiver operating characteristic (ROC) curve will make the results more convincing. The ROC curves are often the only valid method of evaluation (Fawett, 2006). An ROC curve is a detection performance evaluation methodology that shows how successfully a certain detector can separate two groups in a quantitative manner (Fawett, 2001). An ROC curve gives the trade-off between the probability of detection or true positives rate, also called sensitivity and recall versus the probability of false alarm or false positives rate (Fawett, 2006).

In order to compute the performance of the developed diagnosis methodology, the ten folds crossvalidation experiments mentioned in section 3.3 are re-used in this validation phase. Based on the correct classification of patterns in the test set (in each fold), the performance measures are computed.

The average of these performance measures is computed to report the overall performance of the proposed method. To handle the problem of multiple classes, we used the one-against-all strategy: Each class was treated as appositive class in turn, while all remaining six classes were combined to form a negative class.

In each fold, based on correct classification of test patterns, the following parameters are defined:

- Based on TP, TN, FP and FN, the sensitivity (Sn), and specificity (Sp), were defined as follows.

21 TP Sn TP FN   (18) TN Sp TN FP   (19)
Fig. 17 shows one of the found ROC performances of our proposed strategy. It can be seen that PNN-SFAM classifier provides good accuracies of sensitivity (91.7%) and specificity (93.4%) and an area under curves equal to 0.847. We should note that this curve has the least important features among all ROC curves found which validates the strength of our strategy.

Conclusion and prospects

In this work, a new online automated bearing fault diagnosis methodology is proposed. This effective method is introduced by feature extraction, characterization, reduction and classification based on nonstationary vibration signals. Compared to previous works, experimental results demonstrate that the proposed method is extremely promising and it can be summarized as follows:

-EMD is a self-adaptive signal processing method applied to decompose nonlinear and non-stationary vibration signals into some stationary IMFs.

-The diagnosis method is based on eigenvalues of LDA method to investigate a higher classification accuracy compared to the use of EMD energy entropy and PCA features.

-The combination of PNN and SFAM neural networks presents a better categorization than the use of SFAM only.

Finally, the proposed method accurately identifies bearing faults online and a better generalization capability is shown compared to previous works. As future works, we will focus on the variable operating conditions of speed and charge which can generalise the proposed methodology. (1) Initialize: r 0 =x(t), and i=1 (2) Extract the i th IMF c i A) Initialize: h i (k-1)=r i -1, k=1 B) Extract the local maxima and minima of h i (k-1) C) Interpolate the local maxima and the minima by cubic spline lines to form upper and lower envelops of h i (k-1) D) Compute the mean m i (k-1)of the upper and lower envelops of h i (k-1) E) Let h ik =h i (k-1)-m i (k-1) F) If h ik is an imf then set c i =h ik else go to step (B) with k=k+1define the residue r i+1 =r ic G) If r i+1 still has least 2 extrema then go to step (2) with i=i+1 else the decomposition process is finished and r i+1 is the monotonic residue of the signal.

File Recording Interval

Every 10 minutes (except the first 43 files were taken every 5 minutes).

Every 10 minutes Every 10 minutes.

Description

At the end of the testto-failure experiment, inner race defect occurred in bearing 3 and roller element defect in bearing 4.

At the end of the testto-failure experiment, outer race failure occurred in bearing 1.

At the end of the testto-failure experiment, outer race failure occurred in bearing 3.

  vibration signals used in this paper were provided by the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati, USA, in collaboration with the National Aeronautics and Space Administration (NASA). A schematic of experimental test rig is shown in Fig. 3. Four Rexnord ZA-2115 double row bearings are installed on the shaft. Each bearing contains 16 rollers (for each row), a pitch diameter of 2.815in, a roller diameter of 0.331in, and a tapered contact angle of 15.17° (Qiu et al., 2006). Three tests were made. Each test is an experience of 4 bearings. In this way 12 bearings are used but only 4 bearings have reached failure with known defects. Each data set describes a run-to-failure experiment. It consists of individual files that are 1-second vibration signal snapshots recorded at specific intervals (every 10min). Each file consists of 20480 points with the sampling rate set at 20 kHz. Records (row) in the data ASCII files are data points. Data collection is provided by NI DAQ Card 6062E (Ben Ali et al., 2015b). The three tests have been carried out for 35 days until a significant amount of metal debris is found on the magnetic plug of the test bearing. Larger intervals of time stamps (showed in file names) indicate resumption of the experiment in the next working day.

Fig. 4 .

 4 Fig. 4. Each run-to-failure history is divided into three zones: the healthy zone (number 1), the degradation zone (number 2) and the faulting zone (number 3). This partition was done based on the magnitude of vibration signals by proposing two thresholds for each history. Consequently, seven REB states are considered: healthy (H), degraded roller (DR), failure roller (FR), degraded inner race (DIR), failure inner race (FIR), degraded outer race (DOR), and failure outer race (FOR).

i. Scenario 1 :

 1 Combination of original features with SFAM. ii. Scenario 2: Combination of PCA and SFAM. iii. Scenario 3: Combination of LDA and SFAM. iv. Scenario 4: Combination of original features, PNN and SFAM. v. Scenario 5: Combination of PCA, PNN and SFAM. vi. Scenario 6: Combination of LDA, PNN and SFAM.

  4) contain 5394 records were just 3911 are used in the training and testing modes. Consequently, the other 1483 records are unknown for the validated SFAM and PNN neural networks. Experimental results, given in Fig. 12, show that the application of the proposed strategy is very promising for a nearly online bearing diagnosis. The online evaluation of the proposed strategy for the diagnosis of naturally bearing progressing defects is very encouraging. The proposed strategy is able to indicate qualitatively a clear distinction between the healthy and damaged states. In addition, this approach proposes automatic fault detection without any human intervention. The application of the proposed methodology perfectly separates the healthy stage, the degraded stage and the faulting stage. Besides, the type and the severity of defects are given for all tests which provide an idea about the remaining useful life of the tested bearing. However, the three tested bearings are already modeled by the PNN and SFAM neural networks and another test is needed. Another bearing is needed where records were not yet treated by the neural networks. Hence, we thought to use the bearing 3 of testing 3 where any record was used in the training phase. Experimental results shown in Fig. 13 confirm that the proposed methodology is really interesting. The probability density of failure increases at the end of the experiment although that the magnitude of vibration signals remains substantially constant. For more clarity and better resolution, experimental results are enlarged in Fig. 14.

Fig. 15

 15 Fig.15shows the adverse effect of the proposed smart filtering phase on the bearing 3 of testing 3.

  True Positive (TP): It is the number of actually abnormal test patterns which were classified as abnormal by our automated classifiers. -True Negative (TN): It is the number of normal patterns in the test set which were classified into normal class by the automated classifiers. -False Positive (FP): It is the number of normal test patterns classified as abnormal by the automated classifiers. -False Negative (FN): It is the number of abnormal test patterns which are classified as normal by the automated classifiers.

Fig. 1 :

 1 Fig. 1: Probabilistic Neural Network architecture.

Fig. 3 :Fig. 4 :Fig. 5 :

 345 Fig. 3: Experimental test rig, (a) photo from the IMS center (Qiu et al., 2006), (b) simplified schematic of the test rig.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Original feature distributions.

Fig. 9 :Fig. 10 :Fig. 11 :Fig. 12 :

 9101112 Fig. 9: PCA eigenvalue: (a) Entire distribution of covariance matrix, (b) Zoom in the area of the last features

Fig. 13 :

 13 Fig. 13: Application of the proposed method for an online diagnosis of bearing 3 of testing 3.

Fig 16 :Fig 17 :

 1617 Fig 16: Compact box plot of the confusion matrix.

  After 16 sliding steps, determine the highest occurrence as occ=max(k DIR , k FIR , k DOR , k FOR , k DR , k FR ); -Determine the fault with the highest occurrence as occtype= argmax(k DIR , k FIR , k DOR , k FOR , k DR ,

	k FR );

-If occ≥2 then the fault type is occtype else the state of bearing is healthy.

-Run the next 16 classification results.

Table(s) Click here to view linked References Fig

. 2. SFAM neural network architecture.

Text Figure(

s) Table(s) Click here to view linked References

Table(s) Click here to view linked References

	a)												b)
													0.4	H
		0.35										DOR
	IMF 2	0.25 0.3				DOR FOR H					H en IMF 3	0.1 0.2 0.3	0.06 0.08	FOR
	H en		0.2		0.08							0.04
		0.15		0.04 0.06							0.4 0	H en IMF 8	0 0.02
			0.8 0.1	-0.02 0.7 0 0.02 H en IMF 8	0.6						1.4	1.6	1 0.5 H e n F 0.6 I M	0.7	-0.04 -0.02	0.25	0.3	0.35
					H en IMF 1 -0.04	0.5		0.4	0.8	1	H en 1.2	0.8	0.1	0.15	0.11 H en IMF 2 0.12 0.2 -0.06
					-0.06
					0.12							0.1
	c)						0.11						-0.4	d) -0.2	0.09	-0.8
									0.1				0.08	-1
													H	IMF
													-0.6	en	7	H
													en
											0.09
													-0.8
												0.08	-1
								H en IMF	7			H en IMF	6	0.05	H
			0.12				H						0.04	DOR
	H en IMF 4		0.06 0.08 0.1		0.06 0.08	FOR DOR					H en IMF 5	0.01 0.02 0.03	0.04 0.06 0.08	FOR
			0.3 0.02 0.04	H en IMF 8	0 0.02 0.04							0.3 H 0.4 0	0.2	H en IMF 8	-0.02 0 0.02
			0.2 -0.02 0.25						0.3	0.35	e n I F M	0.1	-0.04	0.04	0.02
													3	0.06
													0.25	0.08
					-0.04	0.15						0.1	-0.06
				-0.06 H en IMF 3	0.1			0.15	0.2	H en IMF 2	0	0.12	H en IMF 4 0.12
									0.05	0.1
													0.11
					0.12
													0.1
	e)						0.11		0.1				-0.4	f)	-0.2	0.09	-0.8
													-0.6
													0.08	-1
											0.09		H	IMF
													-0.8	en	7
												0.08	-1
								H	IMF
	H en IMF 6	0.01 0.02 0.03					en	7		0.08	FOR DOR H	en H IMF	6	H en IMF 7	0.02 0.04 0.06	0.06 0.08	FOR DOR H
												0.06
													0.04
													0
			0									0.04
			0.15 0.2							H en IMF 8	0 0.02	0	0.01	0.02	H en IMF 8	0 0.02
													H	-0.02	0.025
			H		0.1						-0.02	0.03 F M I n e	-0.04	0.015	0.02
			e n I M F	0.05					-0.04	0.01	0	5	0.04	0.05 -0.06	0	0.005	0.01	H en IMF 6
												-0.06	0.12
					4								0.02
												0.12
												0.03
							0	0.05		0.04		H en IMF 5 0.11	0.11 -0.2	0.1	-0.4
													-0.4
													0.1	-0.6
													0.09
													0.09	-0.8 g)	-0.6	0.08	-1	-0.8
													H	IMF
													0.08	-1	en	7	H	IMF
													H	IMF	en	6
													en	7	H	IMF
													en	6

Text Figure(s) Table(s) Click here to view linked References

Text Figure(s) Table(s) Click here to view linked References

	a)												b)
		0.8 0.045											0.3	0.015
		0.75 0.04 0.035						H H DOR DOR FOR FOR				0.26 0.28					H DOR FOR		H DOR FOR
		0.7 0.03				0.08 0.08								0.01			0.08		0.08
	H en IMF 1 IMF 7 H en	0.55 0.6 0.65 0.01 0.015 0.02 0.025			H en IMF 8 H en IMF 8	-0.04 -0.04 -0.02 0 0.02 0.04 0.06 0.06 0.04 -0.02 0 0.02						H en IMF 2	0.16 0.18 0.2 0.22 0.24 IMF 8 H en	0.005		H en IMF 8	-0.04 -0.02 0 0.02 0.04 0.06	H en IMF 8	-0.04 -0.02 0.06 0.04 0 0.02
		0.8 0.4 0.45 0.5 0.005 0 0	0.9 0.005	1	1.1 H en IMF 6 H en -0.06 -0.06 1.2 0.12 0.12 0.01 0.015	0.11 0.11	1.3 H en IMF 7 0.1 0.02 0.1 H en IMF 7	1.4 0.09 0.025 0.09	1.5 0.08 0.08	-1 -1	0.8 0.1 0.12 0.14 -0.8 -0.8	0.9 -0.6 H en IMF 6 0 0 -0.6 H en IMF 6	1 -0.4 0.005 -0.4	1.1 0.12 H en -0.06 -0.2 0.01 -0.2 H en IMF 6 1.2 0.11 0.015 H en IMF 7 1.3 0.1 -0.06 0.12	0.09 0.02 0.11	1.4 0.08 0.025 0.1 H en IMF 7	1.5 -1 0.09	-0.8 0.08	-1	-0.6 H en IMF 6	-0.8	-0.4 H en IM -0
												c)
							0.3 0.015
							0.28
																H
							0.26								DOR
																H FOR
							0.24 0.01							0.08	DOR FOR
															0.06
							0.22
						H en IMF 2 IMF 8 H en	0.18 0.2						H en IMF 8	0.04 0 0.02	0.04 0.08 0.06
							0.16 0.005							-0.04 H en IMF 8 -0.02	0 0.02
							0.14							-0.06	-0.02
															0.12	-0.04
							0.4 0.1 0.12 0 0	0.45 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.5 0.55 0.6 0.65 0.7 0.75 0.8 H en IMF 1 0.11 0.11 -0.06 0.1 0.12 0.09 0.1 H en IMF 7 H en IMF 7	0.09 0.08	-1		-0.8	-0.8 H en IMF 6 -0.6	-0.6	-0.4	-0.4	-0.2	-0.2
																	H en IMF 7		0.08	-1	H en IMF 6
		Fig. 8: Bi-dimensional distributions of the lowest features' discriminations (shown in
											Fig.6.g).

Text Figure(s) Table(s) Click here to view linked References

Table(s) Click here to view linked References Fig 15:

  Enhanced nearly online diagnosis of bearing 3 of testing 3 using a smart average technique.

	Other decision
	Online decision	DOR FOR				False arlarms: 1.17%
		H	0	500	1000	1500	2000	2500	3000	3500	4000	4500
							Time (10min)

Text Figure(s) Table(s) Click here to view linked References

Table 1 :

 1 EMD algorithm.

	Text Figure(s) Table(s)
	Click here to view linked References

Table 2 :

 2 IMS bearing dataset description.

	Text Figure(s) Table(s)
	Click here to view linked References

Table 3 :

 3 Testing accuracies of the proposed methodologies.

	Text Figure(s) Table(s)
	Click here to view linked References

Table 4 :

 4 J criterion of original features.

	Text Figure(s) Table(s)
	Click here to view linked References

Table 5 :

 5 J criterion of PCA features.

	Text

Figure(s) Table(s) Click here to view linked References

	LDs	J
	LD 1	3.8574
	LD 2	0.3629
	LD 3	0.228
	LD 4	0.3539
	LD 5	0.5555
	LD 6	0.0513
	All features 6.8044

Table 6 :

 6 J criterion of LDA features.

	Text

Figure(s) Table(s) Click here to view linked References Literature Features/techniques used

			Classifier	REB	Accuracy	Online
				classes	[%]	diagnosis
	(Gryllias	Frequency domain based	Support	3: H, FIR,	100	No
	and	vibration energy features.	Vector	FOR.
	Antoniadis, 2012)		Machine (SVM)
	(Yu et al.,	EMD energy entropy of	Artificial	3: H, FIR,	93	No
	2006)	the first eight IMFs	neural	FOR
			networks
			(ANN)
	(Liu et al.,	EMD features extraction	Particle	4: H, FIR,	97.5	No
	2013)	in time and frequency	swarm	FOR,
		domains	optimization	FR.
			(PSO) and
			SVM
	(Xu et al.,	EMD energy entropy of	Improved	3: H, FIR,	97	No
	2013)	the first six IMFs	PSO and	FOR.
			least squares
			SVM
	(Zarei et	Removing non-bearing	ANN	4: H, FIR,	97.66	No
	al., 2014)	fault component (RNFC)		FOR (one
		filter		hole),
				FOR
				(double
				holes).
	(Pandya et	Hilbert Huang transform	Asymmetric	5: H, FIR,	96.67	No
	al., 2013)	(HHT) based on EMD	proximity	FOR,
			function and	FR,
			K-nearest	combined
			neighbor	defect.
	(Ben Ali et	EMD energy entropy of	ANN	7: H,	93	Yes
	al., 2015b)	the first seven IMFs and		DIR,
		classical statistical		DOR,
		features		DR, FIR,
				FOR, FR.
	This work EMD energy entropy and	PNN and	7: H,	97.82	Yes
		LDA	SFAM	DIR,
				DOR,
				DR, FIR,
				FOR, FR.

Table 7 :

 7 Summary of some studies reporting bearing fault diagnosis.

	Text Figure(s

) Table(s) Click here to view linked References

	Estimated
		H	DR	FR	DIR	FIR	DOR	FOR
	Real
	H	2900	0	0	0	0	0	0
	DR	3	139	0	4	5	0	0
	FR	6	8	423	1	1	7	5
	DIR	3	0	1	140	2	4	1
	FIR	2	0	1	2	46	2	3
	DOR	2	0	0	12	1	130	6
	FOR	1	0	0	0	1	1	48

Table 8 :

 8 Confusion matrix of the proposed method.

	Text Figure(s

) Table(s) Click here to view linked References

Acknowledgment

Authors wish to thank the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati, USA, for his permission to use his bearing data. Also, authors wish to thank Dr Barrie William Jervis, former professor in Sheffield Hallam University (U.K) for his relevant advices during the implementation of PNN and SFAM neural network.


Julien Alexandre
email: alexandre@ensta.fr

Douglas Sandretto

Alexandre Picciani De Suza

Chapoutot
email: chapoutot@ensta.fr

Appropriate Design Guided by Simulation: An Hovercraft Application

A design methodology based on simulation of dynamical behavior is presented in this paper. The particularity of our method is that it exploits the set-membership simulation. Indeed, this method allows one to consider an interval of values for each parameter of the dynamical model. Finding the parameters validating the requirements is then a filter, based on a Branch and Prune algorithm. We prefer an approach of appropriate design that is a design which validates the physical constraints coming from requirements, to an optimal design which does not lead to satisfy imperative requirements. Our method is described and applied on the complex problem of design of an hovercraft under dynamic requirements.

Introduction

During the design process of a complex system such as a vehicle (aircraft, ship, etc.), a chemistry device or a production chain, it is usual to define constraints on the behavior of the futur system (see [START_REF] Diez | Multidisciplinary robust optimization for ship design[END_REF] for a list of design applications), i.e., on its temporal evolution or its dynamics. These constraints come from the requirements defined during the engineering step. The most common approach to model the behavior of a complex system uses differential equations such as Ordinary Differential Equations (ODEs). On contrary, there is no standard method to validate the requirements over the dynamics of the system.

One well known method, consisting in Monte-Carlo simulations for estimating the behavior of a system, is often used to solve design problems. This method cannot try all the possible configurations with continuous valued parameters and quickly leads to a combinatorial problem. Moreover, it is not easy to manipulate uncertainties with this method, which are important in design of complex systems [START_REF] Delaurentis | Uncertainty modeling and management in multidisciplinary analysis and synthesis[END_REF].

In a Multidisciplinary Design Optimization (MDO) process, one may consider the influence of the physical process variability, information uncertainty and the choice of modeling on the design. If it is different than a behavior constrained approach, it focuses on the uncertainties [START_REF] Batill | Modeling and simulation uncertainty in multidisciplinary design optimization[END_REF]. A common solution to the issue of uncertainties is to use a probabilistic approach [START_REF] Mavris | A probabilistic approach to multivariate constrained robust design simulation[END_REF]. Unfortunately, this approach does not provide guarantee on the requirement validation.

Finally, two main approaches can be used to solve a design problem

• Optimal design, which tries to find the best design w.r.t. a given cost, even if not feasible in the sense of requirements;

• Appropriate design, which tries to find one feasible design, but does not consider the cost of the solution.

We can notice that some methods mix both approaches such as [START_REF] Diez | Multidisciplinary robust optimization for ship design[END_REF]. If optimization approach is not directly linked to the possible values of design parameters [START_REF] Kroo | Multidisciplinary optimization methods for aircraft preliminary design[END_REF], the appropriate design based on interval analysis searches design parameters into intervals [START_REF] Merlet | Appropriate Design of Parallel Manipulators[END_REF].

Another similar method uses intervals and Branch and Bound algorithm to find a feasible and optimal solution [START_REF] Stuber | Robust simulation and design using parametric interval methods[END_REF]. Both of these interval-based methods do not consider the dynamics of the studied system. We propose in this paper an interval-based method using set-membership simulation and constraint satisfaction problem (CSP) formalism. Our method is tested on the design of an hovercraft, which is an application very close to [START_REF] Nguyen | A multidisciplinary robust optimization framework for UAV conceptual design[END_REF].

After some preliminary definitions, we present our design methodology in three sections. We describe two tools in sections 2 and 3, before presenting the complete algorithm in Section 4. Then we apply our method on the design of an hovercraft and present the results in Section 5. Section 6 concludes.

Notation and Definitions

The main underlying tool, used in our design methodology, is interval analysis [START_REF] Moore | Interval Analysis[END_REF]. In the following, we will often use the notation [x] ∈ IR (the set of intervals with real bounds) with [x] = [x, x] = {x ∈ R | x x x} denoting an interval. By abuse of notation, [x] will also denote a vector of intervals, i.e., a Cartesian product of intervals, a.k.a. a box. We now introduce two fundamental definitions.

Definition 1 Consider an Ordinary Differential Equation (ODE) with a given initial condition • X = {x 1 , . . . , x n } is a set of variables;

ẋ(t) = f (t, x(t), p) with x(0) ∈ [x 0 ], p ∈ [p], (1)
with f : R + × R n × R m → R n assumed
• D = {[x 1 ], . . . , [x n ]} is a set of domains ([x i ] contains all possible values of x i ); • C = {c 1 , . . . , c m } is a set of constraints. A valuation (v 1 , . . . , v n ), v i ∈ [x i ] with i = 1, . . . , n
, of a CSP can be classified as i) accepted if it satisfies all the constraints; ii) rejected if it does not satisfies at least one constraint; or iii) uncertain if it partially satisfies, in a set-theoretic sense, at least one constraint and satisfies all the others constraints. Many methods have been developed to solve CSPs, such as [START_REF] Lhomme | Consistency techniques for numeric CSPs[END_REF] and [START_REF] Rueher | Solving continuous constraint systems[END_REF].

In order to be able to consider ODEs in CSP formalism, two main approaches can be distinguished. First, [START_REF] Goldsztejn | Including ordinary differential equations based constraints in the standard CP framework[END_REF] has shown that it is possible to consider ODE directly in CSP in an elegant manner by considering the solution x(t) of (1) as a function involved into a constraint c. Second, [START_REF] Cruz | Constraint reasoning over differential equations[END_REF] has defined a new extended formalism named Constraint Satisfaction Differential Problem (CSDP). In this formalism, the differential constraint is clearly identified. We place our approach in the first and classical formalism of CSP, and we do not give more details, the implementation being out of the scope of this paper.

Set-based Simulation

Several methods exist to solve [START_REF] Dit Sandretto | Validated Explicit and Implicit Runge-Kutta Methods[END_REF] involving interval values such as [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF][START_REF] Dit Sandretto | Validated Explicit and Implicit Runge-Kutta Methods[END_REF]. In this section, we describe our approach for validated simulation based on Runge-Kutta methods [START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF][START_REF] Dit Sandretto | Validated Explicit and Implicit Runge-Kutta Methods[END_REF].

A numerical integration method computes a sequence of approximations (t n , x n ) of the solution x(t; x 0 ) of the IVP defined in Equation (1) such that x n ≈ x(t n ; x n-1 ). The simplest method is Euler's method in which t n+1 = t n +h for some step-size h and x n+1 = x n +h×f (t n , x n , p); so the derivative of x at time t n , f (t n , x n , p), is used as an approximation of the derivative on the whole time interval to perform a linear interpolation. This method is very simple and fast, but requires small step-sizes. More advanced methods coming from the Runge-Kutta family use a few intermediate computations to improve the approximation of the derivative. The general form of an s-stage Runge-Kutta formula, that is using s evaluations of f , is

x n+1 = x n + h s i=1 b i k i , (2a)
k i = f t n + c i h, x n + h s j=1 a ij k j , p , i = 1, 2, . . . , s . (2b)
The coefficients c i , a ij and b i fully characterize the method. To make Runge-Kutta validated, the challenging question is how to compute a bound on the distance between the true solution and the numerical solution, defined by x(t n ; x n-1 ) -x n . This distance is related to the local truncation error (LTE) of the numerical method.

To bound the LTE, we rely on order condition [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF] respected by all Runge-Kutta methods. This condition states that a method of this family is of order O iff the O + 1 first coefficients of the Taylor expansion of the solution and the Taylor expansion of the numerical methods are equal. In consequence, LTE is proportional to the Lagrange remainders of Taylor expansions. Formally, LTE is defined by (see [START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF]):

x(t n ; x n-1 ) -x n = h O+1 (O + 1)! f (O) (ξ, x(ξ; x n-1 ), p) - d O+1 φ dt O+1 (η) ξ ∈]t n , t n+1 [ and η ∈]t n , t n+1 [ . (3)
The function f (n) stands for the n-th derivative of function f w.r.t. time t that is

d n f dt n and h = t n+1 -t n is the step-size. The function φ : R → R n is defined by φ(t) = x n + h s i=1 b i k i (t)
where k i (t) are defined as Equation (2b).

The challenge to make Runge-Kutta integration schemes safe w.r.t. the true solution of IVP is then to compute a bound of the result of Equation ( 3). In other words we have to bound the value of f (O) (ξ, x(ξ; x n-1 ), p) and the value of d O+1 φ dt O+1 (η). The latter expression is straightforward to bound because the function φ only depends on the value of the step-size h, and so does its (O + 1)-th derivative. The bound is then obtain using interval analysis tools.

However, the expression f (O) (ξ, x(ξ; x n-1 ), p) is not so easy to bound as it requires to evaluate f for a particular value of the IVP solution x(ξ; x n-1 ) at an unknown time ξ ∈]t n , t n+1 [. The solution used is the same as the one found in [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF] and it requires to bound the solution of IVP on the interval [t n , t n+1 ]. This bound is usually computed using the Banach's fixpoint theorem applied with the Picard-Lindelöf operator, see [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF]. This operator is used to compute an enclosure of the solution [x] of IVP over a time interval [t n , t n+1 ], that is for all t ∈ [t n , t n+1 ], x(t; x n-1 ) ∈ [x]. We can hence bound f (O) substituting x(ξ; x n-1 ) by [x].

From requirements to CSP

During the process of design a system, it is common to define some requirements on its behavior, in particular for functional requirements. One requirement on the behavior of a system can be decomposed as one scenario to follow and one or more criteria to fulfill (see theory of System Engineering [START_REF] Blanchard | Systems engineering and analysis[END_REF]). Our approach exploits the fact that a scenario can be simulated while each criterion can be translated into a constraint. Then, a design problem is equivalent to find the design parameters which validate the requirement by simulating the related scenario and by verifying a CSP. Figure 1 shows an example of the use of constraints on a solution of a validated simulation. The square in black represents the hovercraft position, computed by validated simulation, at some instants. We set constraints which are represented by green boxes (for valid sets) and red boxes (for invalid sets). The constraints can then be of different types, for example: (OK-1) Hovercraft stays inside a valid set ; (OK-2) Hovercraft finishes in a valid set; (NOK-3) Hovercraft avoids an obstacle; and (NOK-4) Hovercraft is not in an invalid state at a given time. Constraints OK-1, NOK-3 and NOK-4 are associated to safety constraints while Constraint OK-2 is associated to a functional requirement.

Appropriate design algorithms

The modeling parameters can be split into two parts: the design parameters and the other ones, such that p = {p design , p model }. The model parameters p model are fixed and can not be changed by the design. Nevertheless, they can be uncertain and then bounded in intervals (due to uncertainties on construction or modeling). The design parameters p design , taken inside intervals of possible values (depending on construction constraints, part available on the market, etc.) will be filtered till they validate all the requirements. For that, we use the classical Branch and Prune algorithm which is efficient for CSP solving [START_REF] Lebbah | Accelerating filtering techniques for numeric CSPs[END_REF][START_REF] Benhamou | Clp (intervals) revisited[END_REF].

OK-1 NOK-4 NOK-3 OK-2
Figure 1: Example of solution by simulation with interval analysis and some allowed constraints in green and prohibited constraints in red.

Branch and Prune

A Branch and Prune algorithm consists on alternatively branching and pruning to produce two subpavings L and S, with L the boxes too small to be bisected and S the solution boxes. We are then sure that all the solutions are included in L ∪ S and that every point in S is solution.

Literally, this algorithm browses a list of boxes and for each one i) Prune: the CSP is evaluated (or contracted) on the current box, if it is a solution it is added to S, otherwise ii) Branch: if the box is large enough, it is bisected and the two boxes resulting are added to the main list, otherwise the box is added to L.

Propagation

If a variable domain has been reduced, the reduction is propagated to all constraints on that variable, which allows one to narrow the other variables domains. This process is repeated till a fixed point is reached (or gain too small w.r.t. computation). In our case, it is important to note that if a parameter design p design is reduced, all the simulation need to be computed again. A diagram of the proposed design methodology is presented in Figure 2 for only one requirement. A requirement comes with a scenario and some constraints to verify over this scenario. This part represents the engineering process of design. Another part, the mathematical one, provides a model of the considered system. This model comes with some parameters that we want to compute in order to obtain a design which validates the requirement. This model is used in the simulation process to obtain the behavior of the system over the scenario. Then, CSPs are applied to this behavior to verify the constraints defined by the requirement. If the result of CSPs is positive then the set of design parameters is accepted, if the solution does not respect one of the CSPs then the set of design parameters is rejected. In other case, we classify it as uncertain and bisect the set of design parameters to perform a new analysis on each part. It is a classical Branch and Prune algorithm supplemented with set-membership simulation to deal with ODEs.

Design Algorithm

Algorithm 1 The design methodology. A general algorithm to handle more than one requirement is presented in Algorithm 1. In this algorithm, we use both CSP for acceptance and CSP for rejection to speed up the algorithm. Note that a CSP for rejection can be seen as the negation b of a Boolean formula b, i.e., the values that satisfy b are not satisfy b and reciprocally. Hence, if we found values (v 1 , . . . , v n ) satisfying a CSP rejection we know that (v 1 , . . . , v n ) does not satisfy the CSP and so we can look for an other valuation.

Require: Stack = ∅, Stackacc = ∅, Stackrej = ∅, Stackunc = ∅, [par]0 ⊂ R n Push [

An hovercraft design

We are currently developing many autonomous small vehicles. One of them is an hovercraft, which has many advantages but also a complex behavior. In order to build an efficient platform, we performed a design step with our method previously presented.

Problem Definition

The starting point to model this particular vehicle is the study of surface ships [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF]. It helps to define a generic model of the dynamics of a rigid-body with 6 degrees of freedom. Applying some physical constraints, such as lying on a plane surface (altitude fixed to zero), and assuming that the body-fixed reference frame is placed at the center of gravity, the mathematical model is given by Equation ( 4) that is

u = vr + 1 m X (4a) v = -ur + 1 m Y (4b) ṙ = 1 I z N (4c) ẋ = cos(ψ)u -sin(ψ)v (4d) ẏ = sin(ψ)u + cos(ψ)v (4e) ψ = r (4f)
The function u stands for rectilinear acceleration in x-axis while v stands for rectilinear acceleration in the y-axis. r stands for the angular acceleration while ψ is the angular velocity. The positions of the hovercraft are given by x and y. The forces and moments considered are related to the motor propulsion, the air resistance, the friction on the ground, damping and rudder. They act on the system by X, Y, N as exposed in equations ( 5) to [START_REF] Delaurentis | Uncertainty modeling and management in multidisciplinary analysis and synthesis[END_REF].

X = F u -µ sec N f (u) u √ u 2 + v 2 - 1 2 ρC d S u u 2 + v 2 u -D 11 u (5)
Y = 1 4A C L S v F u δ -µ sec N f (v) v √ u 2 + v 2 - 1 2 ρC d S v u 2 + v 2 v -D 22 v (6)
N = -L 1 4A C L S g F u δ -D 33 r ( 7
)
where f is a regularization function avoiding discontinuity between the static and the dynamic dry frictions. Table 1 gives the physical meanings of all the parameters used in ( 4) to [START_REF] Delaurentis | Uncertainty modeling and management in multidisciplinary analysis and synthesis[END_REF] with their SI units. Design parameters, i.e., those which will be used in CSP, are in bold font.

Parameters configuration

Table 2 shows the chosen values, points or intervals, of each of the mathematical model parameters. The design parameters m, I z , D, S g , L imply that our research space is in R 5 . The domains adopted by these parameters determine a space region where one will perform the task to find a point in this region with the values that lead to a validated design configuration, w.r.t. requirements. In this application, some parameters are taken in continuous domains (i.e., the mass m) and some other are discrete modeling.

Discrete values represent a finite number of physical components that can be used to implement the hovercraft such as the propulsion propeller diameter D. Remark that our algorithm straightforwardly handles both kinds of parameter values. For example, if we pick two different points in the design parameter space, we are able to simulate a scenario which provides two very different behaviors, see Figure 3.

Requirements on Behavior

In this section, we show how requirements on the behavior of the hovercraft can be translated into a simulation scenario and a CSP. For the hovercraft case study, we will consider four requirements each of them is described in a separated section.

In the rest of this section, a simulation scenario is defined by i) the initial condition associated to (4); ii) the control law used to obtain the desired behavior of the hovercraft; and iii) a simulation duration. Each CSP will be described by the list of variables, their domains and a set of constraints.

Global constraints

Before describing the main requirements we considered for the hovercraft case study, we note that a dynamical systems is constrained by global constraints depending on its intrinsic characteristics, e.g., its maximal possible speed. In other words, we take into account constraints which represents some limitations on the behavior of the hovercraft for all its possible trajectories. These constraints will naturally be added into all the CSP associated to the requirements we considered and can be also seen as an independent CSP.

We will assume that the maximum speed reached by the hovercraft will be V max = 10m/s and its maximum angular speed will be R max = πrad/s. Moreover, we define some constants useful to describe simulation scenario such as the cruising speed V c = 2m/s and the maximum travel distance D max = V max × t f where t f is the simulation duration depending on each requirement.

In Table 3 we give the additional constraints which have to be added in each CSP associated to requirements on the behavior of the hovercraft. More precisely, Constraint C 0 t,u means that the rectilinear acceleration in the x-axis is positive and bounded, i.e., we cannot go backward. Constraint C 0 t,v means that the acceleration in the y-axis of the hovercraft is bounded. Constraint C 0 t,r means that the angular acceleration is bounded. Constraints C 0 t,x and C 0 t,y give a bound on the position w.r.t. the duration of the displacement.

Requirement 1: Cruising Speed

The first requirement on the behavior of the hovercraft is give in Table 4. It represents that the system may have an acceleration allowing to reach the cruising speed V c before 1 second. Note that constraints C 1 t,y and C 1 t,ψ mean that the hovercraft follow a line along the x-axis, i.e., their values are closed to zero within a distance of ε > 0.

Requirement 2: Braking Distance

The second requirement is to have the hovercraft stopping within a distance between 3m and 7m, once the engine is turned off after a run in a straight line at cruising speed. Table 5 gives the mathematical formulation in terms of simulation scenario and CSP. Constraint C 2 t,u means that the rectilinear speed along x-axis is positive and bounded by 0.1m/s at the end of the simulation. Constraints C 2 t,x and C 2 t,y mean that the position of the hovercraft are bounded during the simulation.

Requirement 3: Curve

The third requirement specifies that the hovercraft can follow a curve without deriving too much from the desired trajectory. The scenario, given in Table 6, consists to run in a straight line at the cruising speed and turn at a given time. The requirement imposes that the hovercraft keeps a good behavior in a sense of the CSP. Constraints C 3 t,x and C 3 t,y mean that at the end of the simulation the position of the hovercraft is strictly positive in x and strictly negative in y. In other terms, theses constraints mean that the hovercraft has turn to the right. Note that this requirement is hard to express in terms of radius of curvature that is why we simplified it.

Requirement 4: Perturbation of Symmetry

The fourth requirement consists in verifying the robustness of the design to a possible uncertainty on the mass distribution (expressed by an uncertainty on the gravity center position). For that, the gravity center in y-coordinate is set to y G = [-0.2, 0.2] which represents 10% of the width of the hovercraft. The Simulation scenario Initial condition (u0, v0, r0, x0, y0, ψ0) = (0.1, 0, 0, 0, 0, 0) Control input

: R+ × R 6 Constraints: {C 0 t,u , C 0 t,v , C 0 t,r , C 0 t,x , C 0 t,y } C 0 t,u : (t, u(t)) ∈ [0, t f ] × [0, Vmax] C 0 t,v : (t, v(t)) ∈ [0, t f ] × [-Vmax, Vmax] C 0 t,r : (t, r(t)) ∈ [0, t f ] × [-Rmax, Rmax] C 0 t,x : (t, x(t)) ∈ [0, t f ] × [x0 -D, x0 + D] C 0 t,y : (t, y(t)) ∈ [0, t f ] × [y0 -D, y0 + D]
W (t) = WV max , ∀t ∈ [0, t f ] δ(t) = 0 Simulation duration t f = 1 s CSP Variables {t, u, v, r, x, y, ψ} Domains R+ × R 6 Constraints C 1 t,y : (t, y(t)) ∈ [0, t f ] × [-ε, ε] C 1 t,ψ : (t, ψ(t)) ∈ [0, t f ] × [-ε, ε] C 1 t,u : (t = t f ∧ u(t f ) = Vc) Table 5: Braking distance requirement Requirement 1 -Braking distance Simulation scenario Initial condition (u0, v0, r0, x0, y0, ψ0) = (Vc, 0, 0, 0, 0, 0) Control input W (t) = 0, ∀t ∈ [0, t f ] δ(t) = 0 Simulation duration t f = 10 s CSP Variables {t, u, v, r, x, y, ψ} Domains R+ × R 6 Constraints C 2 t,u : (t = t f ) ∧ (u(t) ∈ [0, 0.1]) C 2 t,x : (t, x(t)) ∈ [0, t f ] × [3, 7] C 2 t,y : (t, y(t)) ∈ [0, t f ] × [-0.5, 0.5]
scenario, given in Table 7, imposes to run in a straight line along x-axis at the cruising speed. After 3 seconds, the requirement is that the hovercraft did not drift more than 15cm from x-axis. Constraints C 4 t,y and C 4 t,ψ mean that for the duration of the simulation the y coordinate and the headings of the hovercraft are bounded.

Results

Applying this methodology to the hovercraft design, a set of appropriate solutions verifying the exigences was found. A particular design solution was chosen to start the hovercraft construction. This latter is chosen in the middle of the validated region, in order to be as robust as possible to the building/assembly uncertainties. Simulation scenario Initial condition (u0, v0, r0, x0, y0, ψ0) = (V, 0, 0, 0, 0, 0) Control input W (t) = To display the set of all solutions found, a projection from R 5 to R 2 is performed by fixing two discrete parameters D and S g . We then obtain three figures gathered in Figure 4.

Conclusion

We presented in this paper a tool for the design of complex systems, described by differential equations with bounded uncertainties. In accordance with the requirements, expressed under the constraint satisfaction problem formalism, our method is able to find design parameters. These parameters validate, even in presence of uncertainties, the correspondence between behavior simulation and constraints defined by requirements. This approach has been applied to the design of an hovercraft. The hovercraft is now in construction and a step of validation will then be done. Obviously, the model being a simplification of the real behavior of the hovercraft, we are waiting for a difference between the expected behavior and the reality. Anyway, this difference can be surpassed by a parameter identification process, which may allow to match the model and the real behavior of the hovercraft. This will be the futur work to validate our methodology.

Figure 2 :

 2 Figure 2: Schematic of the design methodology for one requirement.

Figure 3 :

 3 Figure 3: Two simulation solutions with different punctual values of the design parameters.

Figure 4 :

 4 Figure 4: Solution set for m × L, after projection with I z = [0.0575, 0.065], D = 0.1778 and S g = 0.0176 (a); S g = 0.0235 (b); and S g = 0.0352 (c)

  par]0 in Stack while Stack = ∅ do Pop a [par] from Stack ∀Si ∈ {S1, ...Sm} compute Soli solution using Simulation Si with [par] parameters if ∀j ∈ {1, ..., mi}, CSP i j accept(Soli) then Push [par] in Stack accepted else if ∃j ∈ {1, ..., mi} : CSP i j reject(Soli) then Push [par] in Stack rejected else if width([par]) > lim then ([par lef t ], [par

right ]) = Bisect([par]) Push [par lef t ] in Stack Push [par right ] in Stack else Push [par] in Stackuncertain end if end while

Table 1 :

 1 Model and design (in bold) parameters

	Param. Description	IS Units
	m	Mass	kg
	Iz	Inertia moment in relation	kg.m 2
		to the axis z
	C d	Air resistance coefficient	none
	Su	Frontal surface area	m 2
	Sv	Side surface area	m 2
	ρ	Air density	kg.m -3
	µsec	Ground friction coefficient	none
	D11	Damping in relation to the axis	kg.s -1
		x
	D22	Damping in relation to the axis	kg.s -1
		y
	D33	Damping in relation to the axis	kg.s -1
		z
	xG	x-Coordinate of the gravity cen-	m
		ter
	yG	y-Coordinate of the gravity cen-	m
		ter
	CF u	Propulsion coefficient	none
	D	Propulsion propeller diam-	m
		eter
	CL	Rudder lift coefficient	rad -1
	Sg	Rudder area	m 2
	L	Distance between CG and	m
		the rudder aerodynamic
		center

Table 2 :

 2 Values of the parameters

		Model parameters
	C d	0.6
	Su	0.1276
	Sv	0.2146
	ρ	1.225
	D11 0
	D22 0.2
	D33 0.1
	xG	0
	yG	0
	CF u 1.339
	CL	1.891
		Design parameters
	m	[1.0, 1.3]
	Iz	[0.05, 0.08]

D { 0.127, 0.1524, 0.1778, 0.2032, 0.2286 } Sg { 0.0042, 0.005, 0.0059, 0.0067, 0.0075, 0.0084, 0.0101, 0.0117, 0.0126, 0.0134, 0.014, 0.0151, 0.0168, 0.0176, 0.0196, 0.0201, 0.0224, 0.0226, 0.0235, 0.0251, 0.0268, 0.0279, 0.0302, 0.0335, 0.0352, 0.0391, 0.0402, 0.0419, 0.0447, 0.0453, 0.0503, 0.0587, 0.0671, 0.0754 } L [0.25 , 0.35]

Table 3 :

 3 CSP of global constraints

	Variables: {t, u, v, r, x, y, ψ}
	Domains

Table 4 :

 4 Cruising speed requirement

	Requirement 1 -Cruising Speed

Table 6 :

 6 Curve requirement

	Requirement 3

  1 3 WV max , ∀t ∈ [0, t f ] δ(t) = π 12 1 + sin π 2 t

	Simulation duration
	t f = 3s
	CSP
	Variables
	{t, u, v, r, x, y, ψ}
	Domains
	R+ × R 6
	Constraints
	C 3 t,x : (t = t f ) ∧ (x(t) ∈ [1, ∞])
	C 3 t,

y : (t = t f ) ∧ (y(t) ∈ [-∞, -1])

Table 7 :

 7 Symmetry robustness requirement

	Requirement 4
	Simulation scenario
	Initial condition
	(u0, v0, r0, x0, y0, ψ0) = (Vc, 0, 0, 0, 0, 0)
	Control input
	W (t) = 1 4 WV max , ∀t ∈ [0, t f ]
	δ(t) = 0
	Simulation duration
	t f = 3s
	CSP
	Variables
	{t, u, v, r, x, y, ψ}
	Domains
	R+ × R 6
	Constraints
	C 4 t,

y : (t, y(t)) ∈ [0, t f ] × [y0 -0.15, y0 + 0.15] C 4 t,ψ : (t, y(t)) ∈ [0, t f ] × [ψ0 -π 16 , ψ0 + π 16 ]