Moncef Soualhi
email: moncef.soualhi@enit.fr

Thi Phuong

Khanh Nguyen

Kamal Medjaher

Denis Lebel
email: denis.lebel@metallicadour.fr

David Cazaban-Loustaunau

D Cazaban

Health indicator construction for system health assessment in smart manufacturing

Keywords: Smart manufacturing, smart monitoring, industry 4.0, health indicator construction, signal processing, tool condition monitoring

I. INTRODUCTION

Industry 4.0 is the evolution from automatic systems to cyber-physical systems (CPS) that are fully integrated, automated and optimized [START_REF] Lasi | Industry 4.0[END_REF]. This evolution makes industrial systems more complex and consequently leads to numerous challenges for Prognostics and Health Management (PHM). Indeed, the traditional degradation modeling methods, which are based on the analysis of physical characteristics and system dynamic behaviors, become unfeasible for complex systems. To cope with this complexity, data-driven health monitoring approaches have been developed. They use sensor data to track the system degradation without prior knowledge of the physical nature of the system. The performance of these methods strictly depends on signal processing and health indicator construction techniques.

Considering signal processing methods to construct health indicators, the available studies can be classified into three groups: time, frequency, and time-frequency domain.

In time domain, statistical features (root mean square (RMS), standard deviation (STD), skewness (SKEW), variance (VAR), kurtosis (KUR), crest factor (CF), shape factor (SF), etc.) are extracted from the raw signals to monitor the asset degradation process [START_REF] Gouriveau | From Prognostics and Health Systems Management to Predictive Maintenance 1: Monitoring and Prognostics[END_REF]. These features are widely used due to their easy implementation and their relation with the system's anomalies [START_REF] Nguyen | Analysis and comparison of multiple features for fault detection and prognostic in ball bearings[END_REF]. However, these techniques are not viable for noisy signals. They could require preprocessing tasks to enhance the manipulated raw signals before the construction of the health indicators.

The frequency domain analysis allows visualizing the frequency amplitudes by decomposing the signal into a number of discrete frequencies/spectrum of continuous frequencies. For example, the fast Fourier transform (FFT) is used to observe the harmonic frequency amplitudes characterizing faults in planetary gearboxes [START_REF] Zhang | Stator Current Analysis From Electrical Machines Using Resonance Residual Technique to Detect Faults in Planetary Gearboxes[END_REF]. This method is easy to implement and detect faults that can be represented by characteristic frequencies (e.g. bearings and gears). However, these frequencies cannot be generalized to all systems.

Finally, the time-frequency domain methods use both time and frequency analysis to construct the health indicator. They are good candidates for non-stationary signals. The most effective techniques used in the literature are the short Fourier transform (SFT), Wigner-Ville distribution (WVD), wavelet transform (WT), and Hilbert Huang transform (HHT) [START_REF] Soualhi | Early detection of bearing faults by the hilbert-huang transform[END_REF] to extract the useful information of the system state. These approaches combine the advantages of the time domain analysis and the spectrum content in the frequency analysis to reduce the uncertainties and false alarms.

Although the traditional features in time, frequency, and time-frequency domain are widely proposed in numerous studies of the literature, the robustness of these indicators was not verified under different operating conditions. In addition, to our humble knowledge, the existing works in the literature do not investigate the health indicator performance with different signal types. This paper aims to fill this literature gap. It proposes a new health indicator that is constructed by a combination of different features in order to efficiently characterize different system states. The robustness and performance of the proposed health indicator are tested on different sensor measurements, and also on various operating conditions.

The remainder of this paper is structured as follows: Section II presents the global methodology used for the system's health assessment and details the main steps to construct the health indicator. In section III, the performance of the proposed health indicator is highlighted through numerous experimental tests carried out on a robot cutting tool. Finally, the conclusion and perspectives of this work will be presented in section IV.

II. METHODOLOGY FOR SYSTEM HEALTH ASSESSMENT

This section aims to present the main steps to monitor and asses the health state of a given system. The methodology goes from system analysis to fault detection and diagnostics [START_REF] Atamuradov | Prognostics and health management for maintenance practitioners-review, implementation and tools evaluation[END_REF][START_REF] Gouriveau | From Prognostics and Health Systems Management to Predictive Maintenance 1: Monitoring and Prognostics[END_REF] as shown in Fig. 1.

Ω1 Ω2 Ω3 Ω4 Indi 2 Indi 3 1 Indi 1

System Sensor instrumentation Data acquisition

Feature extraction Health indicator construction Fault detection and diagnostics

Ratio! = Peak_to_peak (y(!")) max (fft) (y(!")) indi= Rms!(ratio) StD!(ia) Y(ijh) Y(ijn)
Fig. 1: General overview of the proposed methodology.

First, it is necessary to analyze the system to identify the critical components leading to the system failure [START_REF] Gouriveau | From Prognostics and Health Systems Management to Predictive Maintenance 1: Monitoring and Prognostics[END_REF]. Once the critical components are identified, it is important to define the appropriate physical parameters to be observed and the instrumentation of the sensors to track the component degradation process. The data recorded will then be injected into signal processing algorithms to extract relevant features and build health indicators. Finally, the health indicators will be used to create a model for fault detection and diagnostics (FDD) of the studied system.

Among the system health assessment process, this paper focuses on the health indicator construction step. The developed health indicator is a combination of different features, that are extracted from both time and frequency domains. Its performance is verified on different operating conditions and also on various signal types. In details, the current (i a , i b , i c ), the vibration (acc a , acc b , acc c ), the force (F a , F b , F c ) and the torque signals (T a , T b , T c ) are used to illustrate the health indicator effectiveness. The main steps of the proposed methodology are presented hereafter. 2) Signal processing: Considering Fig. 3, each signal segment is divided by the max value of the segment spectrum, max(FFT(y ijh )). This step is inspired by the results obtained in [START_REF] Ondel | Fault detection and diagnosis in a set inverter-induction machine through multidimensional membership function and pattern recognition[END_REF]. It allows limiting features dispersion caused by load variations and then separating the classes that represent different load levels of the operating mode.

Fig. 3: Time and frequency domain analysis.

3) Feature extraction from every signal window y ijh : This step aims to extract the statistical features from the splitting signals to construct the health indicator (HI). Its expression is denoted by the following equation:

HI = KUR (y ijh ) × VAR (Y ij ) 2 (1)
The KUR and VAR are respectively the kurtosis and the variance values. The kurtosis measures the segmented signal flattening, which indicates the impulsive property of the signal by the centered moment of order 4 of the segmented signal y ijh [START_REF] Soualhi | Data fusion for fault severity estimation of ball bearings[END_REF]. It reduces signal's noises due to the sensor measurement. The variance measures the dispersion of the signal Y ij [START_REF] Soualhi | Data fusion for fault severity estimation of ball bearings[END_REF].

KUR = 1 ne ne n=1 (y ijh -y ijh ) 4 1 ne ne n=1 (y ijh -y ijh ) 2 2
(2)

VAR = 1 N e N e n=1 Y ij -Y ij 2 (3)
where ȳ ijh and Ȳij are respectively the mean values of the segmented signal and the total signal.

4) Construction of health indicators:

The obtained results by ( 1) are used to construct the health indicators to detect the system's anomalies and also to classify the different system's states. Note that the use of one sensor measurement, for example only the phase a of the current, is not sufficient to detect different fault types. The combination of several health indicators obtained from sensor measurements, for example, all the three phases of the current, is necessary to clearly separate different system's states. The health indicator presentation is therefore changed from one dimension to three dimension space. In detail, the obtained vector from each signal is used to build the health indicator matrix as illustrated hereafter.

indi ijh =           indi ij1 indi ij1 . . . indi ij1 indi ij2 indi j2 indi ij2 indi ij3 indi j3 indi ij3 . . . . . . . . . indi ijN e indi ijN e . . . indi ijN e          

III. CASE STUDY

This section presents the case study used to highlight the performance of the proposed health indicator. The test bench is a six-axis robot used for machining an aluminum blank, see Fig. 4. This robot is composed of six servo-drive motors to control the robot-arms positions. At the sixth axis of the robot, a tool is placed for machining the workpiece. The machining tool contains the unit (flat-end mill) that is the most critical component in the machining processes. To assess the component's health state, current, vibration, force, and torque measurements are the most suitable parameters to monitor this kind of process. First, the subsection (III-A) describes the experimental process. Then, the subsection (III-B) verifies the performance and robustness of the proposed health indicator on different operating conditions. Finally, the fault diagnostics results are presented in the subsection (III-C).

A. Test bench description

The overall scheme of the test bench is presented in Fig. 4. It is installed at METALLICADOUR, a resource and technology transfer center in the south-west of France.

In detail, a three-phase synchronous motor drives a cutting tool. This latter is used for machining aluminum blanks by 4 level steps of 5 mm as marked in green color in Fig. 4. Note that the cutting tool contains three cutting edges. In the case study, four components are used to achieve the experimental tests. The first one represents a new tool. The second, third and fourth tools represent different degraded states. These are respectively the first, the second cutting tool defect, and the broken tooth state, see Fig. 5.

Healthy

Surface damage level 1

Surface damage level 2 Teeth broken

Fig. 5: Illustration of different tool experiences.

Concerning the machining process, the manufacturing of an aluminum blank is investigated. This object is a small part used in aircraft. The machining parameters, which correspond to two different operating conditions (i.e two cutting depth level) of the machining process, are summarized in table I.

Regarding the data acquisition part, firstly, three-axis force and torque sensors are placed at the sixth axis of the robot. Secondly, the three axes vibration sensors are placed as near as possible to the cutting tool. Finally, the three-phase current sensors are placed at the output of the inverter of the machining tool. The sensor data are then recorded by a National Instrument software (Labview) with a sampling frequency of 25.6 kHz and a duration of 40 seconds. These data are saved in csv files of 5 seconds of data in each file.

B. Investigation of the proposed methodology performance

In this subsection, the performance of the proposed health indicator is highlighted when compared with the classical ones proposed in the literature, such as (STD and KUR). The results are shown in Fig. 6. Fig. 8 shows also a negligible dispersion of health indicator observations within a group and a large distance separation between different groups representing the tools health states.This proposed methodology offers a promising way for assessing the system's health, e.g. when using the current sensors, which are not expensive and are easy to mount on the machine without perturbing the machining process.

C. Fault classification

In this section, to verify the fault detection and diagnostics performance, the constructed health indicators matrix is divided into two sets: training (50%) and testing sets (50%). In this application, the number of considered system's states is equal to 4 denoted by s, where s = 1 represents the healthy state and the rest of the values correspond to three fault types. The training set is used to train an ANFIS model [START_REF] Jang | Anfis: adaptive-network-based fuzzy inference system[END_REF] to classify tool's states.

The classification results obtained by ANFIS are verified by considering two operating conditions, see table I. It can be seen that, in all cases, the accuracy score given by the ANFIS classifier is 100%.

TABLE II: Accuracy score (%) using ANFIS at different operating modes. 2 operating modes: C1 (5 mm, 1890 mm/mn, 14000 tr/mn), C2 (10 mm, 1890 mm/mn, 14000 tr/mn).

Case

C1 C2 ANFIS 100 100

IV. CONCLUSION

In this paper, a new health indicator has been presented for a system's health assessment. This HI is a combination of statistical features extracted from the time and the frequency domains. Then, the proposed HI is fed into the ANFIS model to classify the different health states of the monitored system. The methodology performance was highlighted when investigating a machining robot test bench, where the studied system was a cutting tool (end flat mill), that is used for milling an aluminum workpiece. The experiment tests were performed in different operating conditions and monitored by different parameters (current, vibration, force, and torque signals) in three-dimensional space. The obtained results highlighted the robustness of the proposed health indicator regardless of the system operating modes and of the signal types.

As a perspective of this work, the proposed health indicator will be applied in different complex systems with various operating conditions to investigate its robustness and effectiveness.

1 )

 1 Load and split every measured signal into N e observations of length L: This step, illustrated in Fig.2, allows reducing the data size for signal processing and takes only the relevant features. It splits the signals into N e segments denoted y ijh of length L, where i and j characterize the signal types (current, vibration, force and torque) and the axes (i a , i b , i c , acc a , acc b , acc c , F a , F b , F c , T a , T b , T c ) respectively with h representing the segmented signal samples (h ∈ [1...N e]).

Fig. 2 :

 2 Fig. 2: Sampling signals of different measurements.

Fig. 4 :

 4 ABB 6 axis robot

Fig. 6 :

 6 Fig. 6: Health indicators construction with traditional features.

Fig. 6

 6 Fig.6shows clearly that the traditional features (STD and KUR) proposed in the literature presents a large dispersion of the observations in the same group. Thus, one can conclude that when using these traditional features, different system health states cannot be detected. Next, the performance of the proposed health indicator is then proved when considering all signal types: the current, the vibration, the force, and the torque signals. First, the proposed HI is used to diagnose different states of the cutting tool that works in the first operating mode (5 mm, 1890 mm/mn, 14000 tr/mn).

Fig. 7 :Fig. 8 :

 78 Fig. 7: Health indicators construction for TCM in first operating condition: 5 mm, 1890 mm/mn, 14000 tr/mn.From Fig.7, it can be seen that the proposed methodology allows clearly separating different health states of the cutting tool, regardless of the signal type. In addition, another operating mode (10 mm, 1890 mm/mn, 14000 tr/mn) is investigated to highlight the HI performance.

TABLE I :

 I Machining parameters for experimental tests.

	Cutting depth Feed rate Speed Health state
	(mm)	(mm/mn) (rpm)	HS
	5	1890	14000	Healthy
	5	1890	14000	Faulty 1
	5	1890	14000	Faulty 2
	5	1890	14000	Faulty 3
	10	1890	14000	Healthy
	10	1890	14000	Faulty 1
	10	1890	14000	Faulty 2
	10	1890	14000	Faulty 3

ACKNOWLEDGMENT

The project has been 65% cofinanced by the European Regional Development Fund (ERDF) through the Interreg V-A Spain France Andorra programme (POCTEFA 2014-2020). POCTEFA aims to reinforce the economic and social integration of the French-Spanish-Andorran border. Its support is focused on developing economic, social and environmental cross-border activities through joint strategies favouring sustainable territorial development.


Ferhat Tamssaouet

Tp Khanh Nguyen

Kamal Medjaher

Marcos E Orchard

Uncertainty Quantification in System-level Prognostics: Application to Tennessee Eastman Process

This paper addresses the problem of uncertainty quantification in system-level prognostics. To this purpose, a three-step methodology, based on the inoperability input-output model, is presented. The first step concerns the estimation of the system inoperability, using a new adapted particle filtering method, while considering the interactions between its components. The second step focuses on the long-term prediction of the system inoperability in order to determine its evolution. Finally, in the third step, a method for calculating the remaining useful life of the system, based on the system configuration, is formulated. The proposed methodology is applied on data obtained from the Tennessee Eastman Process simulations to predict the shutdown due to violation of process constraints.

I. INTRODUCTION

Prognostics and Health Management (PHM) is essential to ensure safe, reliable and correct operation of complex technical systems. Among the key elements of PHM, prognostics allows predicting the remaining useful life (RUL) of components, subsystems or systems before they become inoperable. Based on these predictions, effective actions can be taken to minimize losses, optimize maintenance, and extend components life.

According to practical requirements, prognostics received a great attention in literature. However, it has often been approached from a component view without considering interactions with other system components and the environment [START_REF] Atamuradov | Prognostics and health management for maintenance practitioners-review, implementation and tools evaluation[END_REF], [START_REF] Rodrigues | Remaining useful life prediction for multiplecomponent systems based on a system-level performance indicator[END_REF]. Hence, for complex engineering systems, it is necessary to study the concept of failure prognostics at system-level by considering the mutual interactions between its elements.

Moreover, notwithstanding the increasing accuracy and precision of prognostics algorithms, their objects of study, i.e. degradation and failure mechanisms, remain stochastic phenomena and, therefore, the uncertainty cannot be eliminated totally [START_REF] Saha | Uncertainty management for diagnostics and prognostics of batteries using bayesian techniques[END_REF]. Indeed, various sources contribute to make the estimation and prediction of one system state uncertain. The number of the uncertainty sources will rapidly increase when considering prognostics at system-level.

To provide a solution to the problems stated above, a new method that allows quantifying the uncertainty in the system remaining useful life (SRUL) prediction is proposed in this paper. It introduces, as a first contribution, a novel systemlevel prognostics framework based on the inoperability inputoutput model (IIM). This model allows tackling the issue 1 F. Tamssaouet, TPK Nguyen and K. Medjaher are with Laboratoire Génie de Production, LGP, Université de Toulouse, INPT-ENIT, Tarbes, France ferhat.tamssaouet@enit.fr related to the interactions between the system components. As a second contribution, a methodology to quantify uncertainty in the SRUL predictions based on the particle filtering is proposed. This methodology will be applied to the wellknown Tennessee Eastman Process (TEP) to demonstrate its effectiveness.

The remainder of the paper is organized as follows. Section II presents the new system degradation model which is based on the inoperability input-output model. Section III describes the proposed methodology for uncertainty quantification and the SRUL determination. Section IV deals with uncertainty modeling and quantification for the Tennessee Eastman Process. A comparison of the obtained results with process real data, were made in order to show the effectiveness of the proposed method. Finally, Section V concludes the paper and gives some future works.

II. INOPERABILITY INPUT-OUTPUT MODEL

One of the main challenges for the SRUL prediction is to develop a model that allows taking into account the mutual components interactions and effects of the mission profile on the degradation evolution. For this purpose, a new model for the degradation of multi-component systems is proposed in this section, that is based on the inoperability input-output model (IIM). The proposed model is then used in Section III to quantify the uncertainty when predicting the SRUL.

The proposed model is inspired by the IIM which is an extension of the input-output model developed by Leontief Wassily in 1936 [START_REF] Leontief | Quantitative input and output relations in the economic systems of the united states[END_REF]. The IIM model and its variants are usually used to investigate the global effects of negative events on highly interdependent infrastructures or multisector economies [START_REF] Haimes | Leontief-based model of risk in complex interconnected infrastructures[END_REF], [START_REF] Santos | Modeling the demand reduction input-output (i-o) inoperability due to terrorism of interconnected infrastructures[END_REF]. This is achieved by using the concept of inoperability, which is defined as the inability of a system to perform its intended functions. The aptitude of the IIM to consider mutual interactions between numerous elements offers a promising perspectives when applying it in PHM domain.

The IIM adapted to prognostics is proposed firt in [START_REF] Tamssaouet | System-level prognostics based on inoperability input-output model[END_REF] and is represented by the following formula:

q(t) = K(t).[A.q(t -1) + c(t)] (1)
where:

• q(t) is a vector representing the overall inoperability of the system components at time t. Each component of this vector is a value between 0 and 1, where q i (t) = 0 corresponds to a healthy component (with an ideal performance) and q i (t) = 1 to a faulty component (no longer able to perform its tasks).

• A is a matrix representing the interdependencies between the system components. Each element a i j of the matrix corresponds to the influence of the inoperability of component j on the inoperability of component i.

The bigger a i j is, the greater is the influence of j on i. • c(t) is a vector representing the internal inoperabilities of the system components at time t, i.e. the degradation of the component due to wear, corrosion or any other failure mechanism. The parameter c i (t) can be obtained by normalizing the health indicator of component i to its failure threshold. More details about the normalization of the component's health indicator can be found in [START_REF] Tamssaouet | System-level prognostics based on inoperability input-output model[END_REF]. • A.q(t) represents the inoperability of a component due to its interdependencies. This quantity informs about the degradations caused by the interactions between components.

• K(t) is a diagonal matrix representing the factors influencing the inoperabilities of components at time t with respect to the system inputs (mission profiles and environment conditions). Each element k i is specific to only one component i. As one can notice in (1), the degradation of component i, characterized by an inoperability q i (t), depends on its inherent natural degradation mechanisms expressed by c i (t) and the degradation induced by the interactions with other components through the term A.q(t). By integrating these two types of degradation, IIM can estimate the health state of systems more accurately.

The advantages of using the IIM to model systems are multiple, among which: 1) the normalization of health indicators to obtain inoperability allows modeling systems with heterogeneous components (different health indicators, range values, degradation patterns and failure thresholds); 2) IIM describes a direct relationship between the mission profile effects and the degradation evolution, which eases the adaptation of the mission profile to extend the system life; 3) multiplying the inoperability by 100 gives a percentage of the component degradation relative to its failure threshold, which facilitates communication with the decision-makers.

III. UNCERTAINTY QUANTIFICATION IN SYSTEM-LEVEL PROGNOSTICS

The methodology proposed in this paper, and illustrated by Fig. 1, combines the estimated and the predicted system inoperabilities to compute the SRUL. The computational process requires the component-level degradation models, the interactions between the components, the thresholds related to each failure mode and the distributions associated to the uncertainties. The details of the three main steps of the proposed methodology is explained in the next subsections.

A. Inoperability uncertainty estimation

The objective of the first step of the methodology is to estimate the inoperability posterior density of the M system components at each time instant k given the observations y k . To do that, the particle filtering, which is a popular Fig. 1: Uncertainty quantification methodology for systemlevel prognostics.

technique explored by several works in prognostics domain [START_REF] Acuña | Particle-filtering-based failure prognosis via sigma-points: Application to lithium-ion battery stateof-charge monitoring[END_REF], [START_REF] Orchard | A particle filtering-based framework for on-line fault diagnosis and failure prognosis[END_REF], is used. This tool can be applied to systems with nonlinear dynamics and non-Gaussian noise. However, contrary to traditional utilization, in this paper a particle is considered as a vector representing the state of health (inoperability) of the system components. Thus, the weight associated to a particle represents the approximation of the inoperability probabilities of all the M components at the same time, as shown in Fig. 2. The process of estimating the inoperability state of a system at time k is explained below.

Firstly, using the IIM presented in Section II, the prior probability density distributions PDFs of the system components inoperabilities p(q k |q k-1 ) at time k are predicted based on the ones at the previous time k -1:

p(q k |q k-1 ) ∼ IIM(q k-1 ) (2)
Next, given new observations yk i at time k for a component i, i ∈ {0, 1, ..., M}, the system posterior PDFs inoperabilities are updated by the particle filtering. In detail, considering a set of N particles {q (l) } l=1,...,N , their associated normalized weights {w (l) } l=1,...,N are evaluated by the likelihood functions p(yk

i |qk i ) using the importance distribution functions π(qk i |qk-1 i , y1:k i ): w (l) k ∝ w (l) k-1 M ∏ i p(yk i |q (l) k i )p(q (l) k i |q (l) k-1 i ) π(qk i |qk-1 i , y1:k i ) (3)
Finally, to overcome the degeneracy problem, a resampling Fig. 2: Inoperability PDFs of a system with 3 components. process is applied in each time step to replace particles having low importance weights with particles that have higher importance weights. The posterior PDFs of the system inoperability at time k (Fig. 1) can be approximated before the resampling step by:

p(q k |y 0:k ) ≈ N ∑ l=1 w (l) k δ (l) q k (q k ) (4)
where δ (•) denotes the Dirac delta function. The estimation procedure is repeated at every instant k, k ∈ {1, 2, ..., k p }, where k p is the starting time of the prediction step presented in the next subsection.

B. Inoperability uncertainty prediction

Prognostics, and thus generation of long-term predictions, is a problem that goes beyond the scope of filtering problem, since it involves future time horizons in which no measurements are available for the Bayesian updating through the equation (3). Thus, the particle filtering, which is more suitable for estimation problems, needs to be adapted to use it for predictions.

In this work, to reduce the computation requirement, we suggest to follow the procedure proposed in [START_REF] Doucet | On sequential monte carlo sampling methods for bayesian filtering[END_REF] and which is based on the assumption that the particle weights are constant from time k p to time k. According to this procedure, the predicted PDF of the inoperability of the system's components at time k (i.e., p(q k |y 1:k p )) can be obtained by applying recursively (2) to q (l)

k p .
Once the prediction of the future system inoperability is done, it will be used to determine the system remaining useful life (SRUL), as explained in the next subsection.

C. SRUL determination

The SRUL provides information related to the time when the whole system fails (i.e., when the combined failures of individual components lead to system failure) [START_REF] Rodrigues | Remaining useful life prediction for multiplecomponent systems based on a system-level performance indicator[END_REF] or when a system reaches performance level that is considered unacceptable. However, the consequence of the degradation of one or more components depends on the considered architecture (e.g. parallel or series). Therefore, the SRUL must be calculated according to system configuration.

Assuming that the system is healthy at time k pth, moment when the prediction algorithm is launched, the SRUL can be computed as follows:

SRUL = τ F -k p (5)
with τ F is the system time-of-failure ToF (or the system end-of-life (EOL)). ToF is chosen in this work because it is a more general concept which can be used in multiple applications [START_REF] Acuña | A theoretically rigorous approach to failure prognosis[END_REF].

τ F = in f (k ∈ N : system f ailure at k) (6)
In practice, and given the complexity of industrial systems, it is important to consider the uncertainty associated with the ToF. To do this, the notations and the new paradigms proposed in [START_REF] Acuña | Particle-filtering-based failure prognosis via sigma-points: Application to lithium-ion battery stateof-charge monitoring[END_REF], [START_REF] Acuña | A theoretically rigorous approach to failure prognosis[END_REF] are used in the remainder of this paper.

Let's denote a healthy system (with no occurrence of catastrophic failure) and a faulty system (with occurrence of catastrophic failure) at kth by H k and F k , respectively. Let's also consider

H k p :k = (H k p , H k p + 1, • • • , H k )
as the sample space that determines all possible sequences where a system has not catastrophically failed until the time k. Then, according to the definition of the conditional probability, the failure probability at kth is given by:

P(F k ) = P(F k , H k p :k-1 ) P(H k p :k-1 |F k ) (7)
As the system can only fail once (without maintenance), given that the failure has occurred at time k, the probability of staying healthy until time k -

1 is P(H k p :k-1 |F k ) = 1. P(F k ) = P(F k , H k p :k-1 ) = P(F k |H k p :k-1 )p(H k p :k-1 ); ∀k > k p (8) where P(F k |H k p :k-1
) is given by:

P(F k |H k p :k-1 ) = R nq p( f ailure|q k )p(q k |y 1:k p )d q k (9)
The second term of (8), p(H k p :k-1 ), stands for the probability that one component is healthy from k p -th until time (k -1)th, which can be expressed as:

p(H k p :k-1 ) = k-1 ∏ h=k p +1 p(H h |h k p :h-1 ) (10)
As F k and H k are exclusive events, the failure event can be modeled through a Bernoulli stochastic process: p(H j |H k p : j-1 ) = 1p(F j |H k p : j-1 ). It follows that:

p(H k p : j-1 ) = k-1 ∏ h=k p +1 (1 -p(F h |H k p :h-1 )) (11)
The expressions presented in ( 8) and ( 11) are valid whether for prognostics of a single component or complex systems. However, when considering a multi-components system, the way of characterizing p(F k |H k p :k-1 ) will change according to the system configuration.

For exemple, in series configuration of M components, the probability that a system will fail at time k, conditional that it is healthy at k -1, is a finite union of the components failure events. As only one component failure can appear at an instantaneous moment, the components failure events can be considered as incompatible. Then, the system failure probability can be written as:

p(F k |H k p :k-1 ) = M ∑ i=1 p(F i k |H k p :k-1 ) (12)
where p(F i k |H k p :k-1 ) is the probability that component i will fail at time k, conditional that the system is healthy at k -1. Then:

p(F k |H k p :k-1 ) = M ∑ i=1 q k ∈R nq p( f ailure i |q i k )p(q i k |y i 1:kp )d q k ( 13
)
Fig. 3: P&ID of the Tennessee Eastman Process [START_REF] Downs | A plant-wide industrial process control problem[END_REF].

IV. TENNESSEE EASTMAN PROCESS

In this section, the proposed methodology is applied to solve the failure prognostics issue of the Tennessee Eastman Process (TEP).

A. Process description

The Tennessee Eastman Process built by Eastman Chemical Company has been widely used as a realistic benchmark for process control optimization, fault diagnostics and, to a lesser extent, for component-level prognostics. Downs and Vogel [START_REF] Downs | A plant-wide industrial process control problem[END_REF] described it in detail and provided its simulation where the components, kinetics, and operating conditions have been modified for proprietary reasons. The TEP involves five major units (working in open-loop) including a two-phase reactor, a partial condenser, a separator, a stripper, and a compressor. The schematic piping and instrumentation diagram (P&ID) of the TEP is shown in Fig. 3.

In the TEP, the gaseous reactants {A,C, D, E} are fed to the reactor where the liquid products {G, H} and the byproduct {F} are formed through the following exothermic reactions:

         A(g) +C(g) + D(g) → G(l), A(g) +C(g) + E(g) → H(l), A(g) + E(g) → F(l), 3D(g) → 2F(g).
The process has a total of 53 measured variables, of which 22 variables are continuous process measurements (such as: temperatures, pressures, flow rates, levels), 19 variables are composition measurements and remaining 12 are manipulated variables. 28 faults can be injected in the process [START_REF] Bathelt | Revision of the tennessee eastman process model[END_REF], which can be related to step changes, drifts, and random variation of variables, etc.

B. Problem Statement

In this case study, we consider a failure as interruption of the operational continuity resulting from violation of the variables shutdown limits. Therefore, only components with shutdown constraints are considered, i.e. reactor, stripper and separator. Each of these components is monitored by a single parameter: pressure for the reactor, and level for the stripper and the separator. Table I lists the specific operational constraints related to the system parameters that the control system should respect. TABLE I: Process operating constraints [START_REF] Downs | A plant-wide industrial process control problem[END_REF].

Normal operating limits

Shut down limits Process variables Low limit High limit In Matlab R simulations of the TEP [START_REF] Ricker | Tennessee eastman challenge archive[END_REF], three disturbances predefined in [START_REF] Bathelt | Revision of the tennessee eastman process model[END_REF] were injected. Two disturbances are deviations in the reactor and the stripper parameters that are, respectively: deviation in the reactor cooling water flow (c 1 (t)) and deviation in the heat transfer of the heat exchanger of the stripper (c 2 (t)). These deviations c 1 (t) and c 2 (t) are assumed to follow equations ( 14) and (15), respectively.

c

1 (t) = α.c 1 (t -1) + β (14) c 2 (t) = γ.c 2 (t -1) (15)
with α, β and γ are the parameter of the two models. As shown in Fig. 4, when these two faults are injected, the system shuts down at 0.8 hour.

The third injected fault is a random variation of the reactor cooling water inlet temperature and can be considered as process noise.

The purpose, here, is to estimate the SRUL when taking into account the process uncertainty, which is characterized by a random variation of the reactor cooling water inlet temperature. To do this, the proposed IIM, Eq.( 1) where c(t) is given by Eq.( 14) and (15), is used to model the system degradation. The IIM parameters are estimated and adjusted thanks to data collected from the TEP simulation. In the next subsection, construction of the IIM is detailed.

C. Inoperability input-output model of TEP

The raw data acquired from the TEP simulation are normalized to the initial state and the failure threshold, as Fig. 4: Components inoperabilities evolution when considering two faults.

shown in Fig. 4. They represent the inoperability of the three components: reactor q 1 (t), stripper q 2 (t) and separator q 3 (t). Their initial states correspond to the parameters base value [START_REF] Downs | A plant-wide industrial process control problem[END_REF] and the failure thresholds to the shutdown limits values (Table I).

The IIM of the TEP is built from data obtained after injection of these faults. It should be noted here that the IIM model does not model the physical phenomena operating in the process, but is more used to fit the actual data of the process. For this case study, the interdependence matrix (A) is estimated and adjusted to make the proposed model better fit the process data. Then, we obtain the following result:

A =   0 0.1 0 10 -5 0 0 0.4 0.3 0   (16)
Regarding the matrix K, it has been assumed that the environmental conditions have no effects on the evolution of the components inoperabilities. Therefore, its diagonal is equal to 1.

D. Inoperability estimation and prediction

After application of the random variation of the reactor cooling water inlet temperature, the data obtained from the process and the built IIM are used in the particle filtering to estimate the components inoperabilities. To evaluate the inoperabilities densities, 200 particles were used with the initial distributions of the components inoperabilities considered as Gaussian. The selection of the particles to be retained after each filtering step was done by using residual resampling.

When the inoperability of one component exceeds the normal operating limits (as indicated in Table I), the prediction step will be launched (at time k p ).

The results of the estimation and prediction of the inoperabilities uncertainty of the system components are shown in Fig. 5. The reactor is the first component to go out of its normal operating limits after 0.34 hour. This time corresponds to the time where the long-term inoperability prediction is launched. Also, it is the reactor pressure, that triggers the system shutdown (system failure) at 0.44 hour.

From these results, we can highlight the power of estimation and prediction of the proposed methodology. Indeed, before the adding of the random variation (the process noise), the system shuts down at 0.8 hour and after the fault injection, it shuts down at 0.44 hour. Despite this high process variability, the particle filtering was able to estimate the actual inoperability of the system (as shown in Fig. 5). Also, one can notice that the uncertainty related to the predicted inoperability increases for k > k p . This is due to the fact that no measurements were received and, therefore, there is neither updating of the particle weights nor resampling.

E. SRUL determination

For this case study, the operability of the studied system depends on the operability of its components since they all contribute to realization of the system function (G and H production). Therefore, one can conclude that the system has a series architecture. Consequently, the system ToF probability will be determined using the equation ( 13).

Fig. 6 shows the PMF of the SRUL. The mean value of the SRUL is equal to 0.1 hour. The true SRUL, which is equal to 0.102 hour, is within the 95% confidence interval of the predicted SRUL distribution. We can conclude that Fig. 5: Components real and predicted inoperabilities evolution when considering three faults. the predicted SRUL is close to reality and is slightly pessimistic, which does not put the system, its operators and the environment in danger.

In order to discuss the results of the proposed methodology, a study of the impact of the prognostics horizon on uncertainty intervals is performed by considering the α-accuracy metric, which determines whether a prediction falls within an α% interval. In fact, α-accuracy is a useful metric to judge if one prognostics algorithm converges to the true value as more information are accumulated over time. Indeed, a faster convergence is desired to achieve a high confidence, keeping the prediction horizon as large as possible. To that end, in this study, the prognostics time has been varied within the interval shown in Fig. 7 and the accuracy is defined with α = 10%. The figure shows the mean values and uncertainties of the predicted SRUL distributions compared with the true SRUL. As it can be seen, the prediction of the SRUL becomes more accurate each time the measurements are obtained. One can also notice that for t ≤ 0.6, the SRUL remains almost constant because the values of the monitored parameters do not change a lot, as it can be seen in Fig. 4.

V. CONCLUSION

A methodology for the uncertainty quantification at system-level prognostics is proposed in this paper. This methodology results in three main contributions. The first concerns the modeling of the system degradation by using the Fig. 7: SRUL prediction performance with α = 0.1. inoperability input-output model. The second deals with the health state estimation based on an adapted particle filtering. Finally, the third contribution is related to the calculation of the system remaining useful life, based on the recent developments and achievements proposed at componentlevel prognostics and generalized in this paper to systemlevel prognostics.

The proposed methodology was applied on the Tennessee Eastman process in order to predict the shutdown time caused by violation of the process constraints. The obtained results show the effectiveness of the methodology in estimating and predicting the system remaining useful life and characterization of uncertainty.

For future work, a systematic method to determine the IIM parameters from a real data will be proposed. Also, it would be worthwhile to apply the other predefined disturbances in the TEP in order to test the robustness of the methodology.

Fig. 6 :

 6 Fig. 6: SRUL probability distribution of the system.

M. Orchard is with the department of Electrical Engineering, University of Chile, Santiago, 8370451, Chile.


Houda Sarih
email: houda.sarih@enit.fr

Ayeley P Tchangani
email: ayeley.tchangani@enit.fr

Kamal Medjaher
email: kamal.medjaher@enit.fr

Eric Pere
email: e.pere@worldcastsystems.com

Data preparation and preprocessing for broadcast systems monitoring in PHM framework

Keywords: Prognostics and Health Management, Data collection, Data cleaning, Data preprocessing, Useful information

Nowadays, companies producing goods use production systems that are equipped by different sensors in order to monitor efficiently their behavior. Most of the time, the information collected by these sensors is mainly used for production monitoring rather than to analyzing the state of health of the production system. By so doing, these companies have a large and growing amount of data at their disposal. These data make it possible to extract information and knowledge for a better control of the system in order to improve its efficiency and reliability. With the emergence of Prognostics and Health Management (PHM) paradigm few years ago, it has become possible to study the state of health of an equipment and predict its future evolution. Globally, the principle of PHM is to transform a set of raw data gathered on the monitored equipment into one or more health indicators. In this framework, the present paper addresses issues related to raw data. A generic approach is proposed for obtaining monitoring data that are reliable and exploitable in a PHM application. The proposed approach is based on 2 steps: collecting data and preprocessing data. This approach will be applied to a real world case in broadcast industry to show its feasibility.

I. INTRODUCTION

In recent decades, manufacturers have begun to focus on activities related to remote monitoring, fault diagnostic and, more broadly, predictive maintenance solutions. This interest is motivated by the operators desire to reduce operational expenses, such as operating and maintenance costs. In this context, efforts are being made to develop monitoring and diagnostic systems and decision support solutions based on Prognostics and Health Management paradigm. The state of art of these approaches show that various works aiming at the integration of predictive maintenance are published and some solutions already exist on the market. In 2006, ALSTOM launched a preventive maintenance tool called TrainTracer, which collects useful data on the train and captures them from several sources of sensors to continuously monitor the health status of various components. A few years later, the same company launched a predictive tool called "HealthHub" that integrates TrainTracer and automates certain tasks such as alerts and anticipating repairs [START_REF] Lu | Survey on the Railway Telematic System for Rolling Stocks[END_REF]. Another solution for the implementation of predictive maintenance, known as "Orbita", has been developed by Bombardier Transportation, which combines data from rolling stock with Bombardier's engineering know-how and fleet knowledge. This capability improves active vehicle management and infrastructure maintenance, performance and operation [START_REF] Dub | predictive asset management[END_REF]. In 2015, SNCF also introduced predictive maintenance for its trains and pays particular attention to the relevance of the data collected [START_REF] Versailles | Digitalization in Maintenance at SNCF[END_REF]. In 2017, as part of the Connect project, Air Liquide created a remote control center for all its French plants, whose challenges are: the development of predictive maintenance, data recording and analysis and the measurement of the slightest discrepancies that could be a sign of the failure of a particular installation [START_REF][END_REF]. At the same year, Air France presented its Prognos range of two dedicated predictive maintenance solutions at the Paris Air Show 2017: "Prognos for Engine" for motors, and "Prognos for Aircraft" for planes. This analysis tool makes it possible to identify and understand the wear of various components, plan their replacement before they fail and avoid costly maintenance operations [START_REF]La maintenance prédictive avec Prognos[END_REF]. Table I summarizes some companies that have established predictive maintenance tools, their industrial sectors and the names of the tools. The targeted application of the approach being developed in this paper is a broadcast system from WorldCast Systems company. At WorldCast Systems, which is a truly global actor offering a multitude of product lines for the broadcast market, experts and engineers created a proactive service for their products that sends to the user regular reports on the status of key parameters such as temperature, current and voltage, as well as information on components performance and life. Condition monitoring (CM) is one of WorldCast Systems historic businesses and such functions are integrated into all the equipment of various ranges. However, even if the latter have knowledge of the physical behavior of the systems, they encounter difficulties when it comes to exploit this knowledge and implement PHM applications. Because this process requires to transform the raw data into useful information, which makes it possible to estimate and predict the systems health state at every instant. These few analyzed systems show the interest in predictive maintenance along with customers requirements that are more and more demanding (Higher performances, low prices, etc.). Moreover, with the complexity of nowadays systems and the constraints the users may face, it is necessary to establish a scientifically proven approach to deal with their maintenance. The main purpose of this paper is to position the industrial problem within a more general scientific framework. The research work that will be presented in this contribution aims to preprare and clean data in order to have exploitable and reliable data to be used for monitoring, diagnostic and prognostic of failures. The objective is to enable anomalies to be detected and faults to be anticipated by taking into account the context and conditions of use (environment, localization, etc.). This research work can thus be positioned in the PHM framework where the physical system is a Frequency Modulation (FM) transmitter monitored by different sensors. In the field of PHM, three approaches can be used for the monitoring purpose [START_REF] Gouriveau | From Prognostics and Health Systems Management to Predictive Maintenance 1: Monitoring and Prognostics[END_REF]: physics-based (also called model-based) approach, data-driven approach and hybrid approach that combines the two first ones. The pluridisciplinarity and complexity of the physical phenomena characterizing industrial systems, added to the relatively high cost of the experimental considerations to generate degradation models, make it difficult to implement the approach based on physical model. This has been reinforced by various collaborations with industrial partners. In this paper, and given the knowledge and the data available, the data-driven approach would be the most appropriate to implement. Indeed, given the complexity of the studied system, which is located in different places in the world and subjected to different conditions and contexts of use, it would be difficult to use the approach based on a physical model. For example, in practice, it is very difficult to obtain an exploitable model of nominal behavior including damage to an FM transmitter. Indeed, the FM transmitter considered here is a complex system involving different physical phenomena (electrical, mechanical, etc.) and consequently, its modeling requires multidisciplinary skills. Furthermore, for the generated model to be usable, it must be simple, if possible linear, with stationary parameters and take into account variable operating and environmental conditions. Thus, the difficulty of constructing easily exploitable models, the relatively excessive cost and the low applicability of the model-based approach led us to choose the data-driven approach. Also, this choice was motivated, among others, by the simplicity of implementation of the approach on industrial systems and its low cost. Each of the aforementioned approaches has its advantages and disadvantages (see table II), which limit their applicability.

In this context, this paper deals with the description of a general approach for collecting and preprocessing data to extract useful information for monitoring, which are representative of the critical degradation mechanisms and which allow extracting relevant features for fault detection, diagnostic and prognostic. For this purpose, two main steps are identified:

• Collecting data. This step consists in specifying technical characteristics of the acquisition system and defining the storage format in order to be used by PHM algorithms.

• Preprocessing data. At this step, it is a question of removing noise from raw data and extracting useful information for PHM algorithms. One can notice that these tasks are basically at the system comprehension and the engineering levels, but they are nevertheless crucial for the success of the PHM process. Indeed, a poor choice of components, of parameters to monitor, or a poor specification of the acquisition system may lead to biased results in detection, diagnostics, prognostics and decision. Furthermore, implementing the approach on a complex system is potentially a challenging task. All this requires multidisciplinary skills and close collaboration with the manufacturer and/or the operator of the system who have precious information. Moreover, in practice, and most of the time, the raw data collected from the system are not directly exploitable and need some preprocessing actions (data cleaning, denoising, resampling, missing data estimation, wrong data substitution, etc.), which is the purpose of this paper. Figure 1 shows the main idea of the contribution.

The remainder of this paper is organized as following. The second section deals with the proposed methodology by presenting its purpose as well as its different steps; the third section is devoted to the application of the methodology to a real world case in the domain of broadcast products in order to show its effectiveness; and finally, the fourth section concludes the paper and presents some future works.

II. METHODOLOGY

The main aim of this paper is to obtain exploitable and reliable data. The purpose of data acquisition is to collect monitoring data. The data collected must make it possible to monitor the health assessment of the system over time and to understand its degradation phenomena. Data can come from different sensors installed on/in the critical components or from inspections performed by operators working on the system. The choice of the sensors depends on the components to be monitored and the type of data to be recovered. This choice should take into account the use of a minimum number of dedicated sensors to collect representative and relevant data and deploy appropriate acquisition facilities. These data are recovered via dedicated acquisition cards or recording files and stored in a secure database. They are then preprocessed to eliminate errors in data reception that might affect the interpretation of degradation phenomena to ensure that they are complete, reliable and ready for further processing. Furthermore, it is necessary to save the data in formats easily accessible by processing algorithms. In a report published by Forbes [START_REF] Gill | Cleaning Big Data: Most Time-Consuming, Least Enjoyable Data Science Task, Survey Says[END_REF], data preparation accounts for about 80% of a data scientists role, with 60% of that time spent with the cleaning and organizing of that data for analysis; see for instance Figure 2 for most occupation of a data scientist.

The overall PHM framework that we are concerned with can be schematically summarized by the Figure 3 where the goal is to obtain monitoring, reliable and exploitable data which represent the degradation mechanisms and which are necessary in a PHM application.

In the following subsections each of the aforementioned steps will be described more precisely.

A. Data collection

In a company there are many products. The first step consists in choosing the system to analyze. One possible criterion to do this is to proceed by the amount of the collected data. The choice will then be made according to the product which has a maximum representative data over a long period [START_REF] Sarih | Critical components identification based on experience feedback data in the framework of PHM[END_REF]. The problem is that the company may have several data sources. Burn-in data that provides information on youthful defects, customer online data and data in the SD card that provides information on the system in normal operation and the wear period (Figure 4). All these information are important for system monitoring. The first solution to this problem is to create a single database to group the information disseminated in the databases, organize them and store them in the same format to help decision-making in the company and facilitate the processing of these data for use in PHM algorithms.

B. Data preprocessing

Data Preprocessing is a technique that is used to convert raw data into clean data set that are reliable and exploitable in a PHM or any other application. In other words, whenever the data are gathered from different sources, they are collected in raw format, which is not suitable for the analysis. Real-world data are often incomplete, inconsistent and lacking in certain behaviors or trends and are likely to contain many errors. Data preprocessing is a proven technique of resolving such issues. This technique prepares the raw data for further processing [START_REF] Sagar | Data Mining Technique: Fundamental Concept and Statistical Analysis[END_REF]. ISO definition for data preprocessing can be found in [START_REF]Condition monitoring and diagnostics of machines Data processing, communication and presentation Part 1: General guidelines[END_REF].

Different steps are involved for Data Preprocessing. These steps are described below [START_REF] Sagar | Data Mining Technique: Fundamental Concept and Statistical Analysis[END_REF]:

Data cleaning: This is the first step which is implemented in Data Preprocessing. In this step, the primary focus is on handling missing and noisy data, detection, and removal of outliers, minimizing duplication and computed biases within the data.

Data integration: In this step, the data, which can be gathered from various sources, are combined. The combined data are then cleaned and used for processing and analysis.

Data transformation: This step is used to convert the raw data into a specified format according to the need of the model. For this purpose, different options exist, such as the following ones: Data Reduction: This step is used to minimize the amount of data that needs to be stored in a data storage environment. After the transformation and scaling, duplicated (or redundant) data are removed and efficiently organized. The creation of a single database was the solution to both problems: the problem of several data sources and also duplicate data. Each line in the table in our database corresponds to one record. One of the main objectives of a database is to avoid data redundancy or duplication. A duplicate is two (or more) identical records. To be careful about this, a primary key is used as a unique identifier for each line.

The overall data preprocessing procedure performed in this work is depicted in Figure 5.

III. CASE STUDY

The methodology described in the previous section is applied to a case study in the field of radio and television broadcasting. This case study consists of analyzing on a FM transmitter from the WorldCast Systems company.

A. Description of the considered system

The FM transmitter is a complex system composed of different subsystems (please refer to subsection 3.1 System selection in [START_REF] Sarih | Critical components identification based on experience feedback data in the framework of PHM[END_REF] for more details)

In the following subsection, we will present the first task of the proposed methodology that is identifying parameters to measure and necessary sensors to this end.

B. Parametres to measure and sensors to use

The parameters which are to be monitored in a PHM implementation can be selected based on their relationships to functions that are crucial for safety, likely to be implicated in catastrophic failures, essential for mission completeness or that can result in long downtimes. The selection is also based on the knowledge of the critical parameters established by past experiences [START_REF] Pecht | Prognostics and Health Management of Electronics: Fundamentals, Machine Learning, and the Internet of Things[END_REF], on field failure data taken from similar products and on qualification testing. For the considered FM transmitter, 14 parameters have been considered because sensors for this purpose are already integrated during the design process of the FM transmitter.

The monitoring data collected by different sensors from several sources, where initially saved as data files in different formats (usually .txt, .csv or .tdms). In our case, the initial data where in XML format. A typical data set are shown on the left side of Figure 6 from which we extract the csv file of the right side of that Figure that is more convenient for analysis by data analysis software such as Matlab.

In order to be able to manage all the company's data from sensors and have only one data source instead of searching for information in several services, the solution was to create a single database by mySQL Workbench and send all the data from different sources into this database. The preprocessing of the latter takes 2 months Data cleaning processes are then used to obtain exploitable data. For instance, in the case of missing data, techniques such as replacement by the value of previous not missed data, average value, or ignoring missed data where used. For missing data, 3 methods are tested with an example of 100 fan speed measurements and 15% missing values. Figure 7 shows an example of these procedures for the fan speed parameter. In our case, the preprocessing of the FM transmitter data took 2 months.

C. Data import and data mining

The studied example of the company WorldCast systems contain a CSV file (comma-separated values) with measurements taken every minute. The typical approach to working with data distributed across multiple files and sources is to create a database, explore it to identify the data values to be analyzed, and then import these values for the entire data set. The MATLAB import tool allows us to select columns in a CSV file and import the selected data into various MATLAB data structures, including vectors, matrices, cell tables and tables. With the import tool, we select CSV file columns and a target format. We can import the data directly from the example file or generate a MATLAB function that imports all files corresponding to the format of the example file. Once the data are imported, preliminary plots are generated to identify patterns and make conversions.

Most real-world data contain missing or incorrect values and, before the data can be explored, they must be identified and processed. After the graphical representation of the example of the measured supply temperature data, we found peaks outside the normal cyclical fluctuations (Figure 8).

There are several ways to automate peak identification. For example, one can apply a smoothing function on Matlab (Figure 9) and locate the peaks by calculating the difference between the smoothed and the original curves (Figure 10).

After removing abnormal peaks from the data, one can notice the absence of missing values (Figure 11). Indeed, a moving average filter smooths the data by averaging a number of points from the input signal to produce each point in the [START_REF] Smith | The scientist and engineer's guide to digital signal processing[END_REF]. This process is equivalent to low-pass filtering with the the response of the smoothing. In equation form [START_REF] Lu | Survey on the Railway Telematic System for Rolling Stocks[END_REF], this is written:

y[i] = 1 M M -1 j=0 x[i + j] (1)
where:

x[ ] is the input signal.

y[ ] is the output signal.

M is the number of points used in the moving average.

By applying data filtering and smoothing methods, it was possible to remove the outliers and the missing data and get reliable and exploitable data for further processing (feature extraction, fault detection, fault diagnosis and fault prognostics).

IV. CONCLUSION AND PERSPECTIVES

In this paper, a methodology has been proposed to obtain reliable and exploitable monitoring data for a PHM application in order to assure predictive analysis and reliability objectives demanded by industrial companies. The developed approach consists in analyzing the system, identifying the parameters to measure along with necessary sensors for this purpose, and the pre-processing tasks to do to obtain PHM exploitable data. The application of this methodology to a real world problem with data stored in different sources and different formats has shownits promising results. The next step will consist in developing algorithms to extract knowledge from these data; this knowledge will be the basics for PHM algorithms to be developed. On this basis, the PHM algorithms development (diagnostics, prognostics and decision) can be considered. For this purpose, it is necessary to give meanings to the preprocessed data by processing them; this will reveal their content and allow estimating the health status of the components or of the system at every instant.

Fig. 1 .

 1 Fig. 1. Idea around the contribution of the paper.

Fig. 2 .

 2 Fig.2. What data scientists spend the most time doing[START_REF] Gill | Cleaning Big Data: Most Time-Consuming, Least Enjoyable Data Science Task, Survey Says[END_REF].

Fig. 3 .

 3 Fig. 3. Collecting and preprocessing data methodology flowchart.

Fig. 4 .

 4 Fig. 4. Collecting data and creating of a single database process.

Fig. 5 .

 5 Fig. 5. Preprocessing data

Fig. 6 .

 6 Fig. 6. Examples of data file formats: on the left a formal XML, on the right a CSV format.

Fig. 7 .

 7 Fig. 7. a) Example fan speed measurements with missing values, b) Method 1: Replacing by the previous value not missing, c) Method 2: Replacing missing values with the average value and d) Method 3: Ignore missing values (Listwise deletion).

Fig. 8 .

 8 Fig. 8. Graphical representation of the "Supply temperature" measurement.

Fig. 9 .

 9 Fig. 9. Graph of the actual and smoothed "Supply temperature" measurement.

Fig. 10 .

 10 Fig. 10. Plotting the difference between actual and smoothed values

Fig. 12 .

 12 Fig. 12. Data without missing values and outliers

TABLE I SOME

 I PREDICTIVE MAINTENANCE SOLUTIONS USED IN INDUSTRIES.

	Company	Sector	maintenance solutions
	Alstom	Railway transport	HealthHub
	Bombardier	Railway transport	Orbita
	Air liquid	Industrial Gas	A remote operation and optimization center
	Air France	Aerospace	Prognos
	General Electric	Energy	Equipment insight

TABLE II SUMMARY

 II OF ADVANTAGES AND DISADVANTAGES OF THE THREE PROGNOSTIC APPROACHES.

• Normalization -In this method, the numerical data are converted into a specified range, for instance to be ranged between 0 and 1 so that scaling can be performed.• Aggregation -The concept can be derived from the word itself, this method is used to combine some groups of data into one group.• Generalization -In this case, lower level attributes are converted to a higher standard.


Saba Amirdehi
email: sabaamirdehi@alstomgroup.com

Baptiste Trajin

Paul-Etienne Vidal

Johana Vally

Didier Colin

Frédéric Rotella

Open Archive Toulouse Archive Ouverte

In order to study the effect of pantograph bounce on railway traction chain a global model including the power supply, catenary, transformer and effects of pantograph bounce has to be established. In this paper, pantograph bounce model considering arc model including all significant parameters such as relative airflow, variable disconnection length is presented. Pantograph arc phenomenon is analyzed in the shape of global model of power input of traction. Model is established using bond graph in order to express differential system of equations that can be easily implemented on numerical simulators. Simulation results are compared with experimental results from real train.

Introduction

Electrical arcing during pantograph bounce in electrified railway is a common phenomenon. Due to mechanical oscillations of train, pantograph and catenary, irregularity of contact wires and influence of high speed airflow, pantograph and catenary can be disconnected for a period of time [START_REF] Habedank | On the mathematical-description of arc behavior in the vicinity of current zero[END_REF], [START_REF] Hongyi | Extended black-box model of pantograph arcing considering varying pantograph detachment distance[END_REF]. Electrical arc can be modeled through a variable and nonlinear resistance. It influences the available voltage for train traction chain and current of catenary as well [START_REF] Hongyi | Extended black-box model of pantograph arcing considering varying pantograph detachment distance[END_REF]. Arcing phenomenon was modeled in detail in previous papers, [START_REF] Hiroshi | Development of a new pantograph contact strip for ultrahigh-speed operations[END_REF], [START_REF] Wang | The pantograph-catenary arc test system for high-speed railways[END_REF], [START_REF] Cassie | An Improved arc model before current zero based on the combined Mayr and Cassie arc models[END_REF]. Most arc models are presented based on the results from laboratory experiments which are not exactly representative for this phenomenon on complete traction chain [START_REF] Wei | Study on Pantograph Arcing in a Laboratory Simulation System by High-Speed Photography[END_REF], [START_REF] Li | Pantograph Arcing's Impact on Locomotive Equipments[END_REF], [START_REF] Midya | DC Component rom Pantograph Arcing in AC Traction System -Influencing Parameters, Impact, and Mitigation Techniques[END_REF]. In order to analyze the global influence of pantograph bounce and electrical arcing, it is needed to have a complete model of power supply, catenary and train traction chain. In this paper, alternating current power chain including pantograph bounce effects is modeled using bond graph. Bond graph allows considering and representing physical phenomenon. Finally state equations of the system are deduced for simulation purposes. In order to validate the proposed model, simulation results are given and compared with real train measurements.

Model of the traction chain

Power part

In order to analyze the effect of electrical arcing during pantograph bounce on the traction chain, it is needed to have an accurate model of whole power supply and traction chain of train. In this approach, power chain from substation to secondary part of transformer is particularly under concern. Electromechanical part and associated power converters of the traction chain are simplified through equivalent voltage sources based on experimental steady state measurements (! "# , ! "$ ). Figure 1 shows the corresponding schematic of the traction chain. At the secondary of the transformer, there are two windings for electromechanical parts and associated converters of traction chain. While train is moving, there are some variations of P model of catenary parameters: resistance, capacitance and inductance of catenary depend of distance between the train and the substation [START_REF] Grainger | Power systems analysis[END_REF]. Note that these parameters may also depend on environmental conditions such as temperature of humidity. However, these effects won't be studied in this paper. Only general values of parameters will be indicated and should be adapted. Electrical arc can be modelled as a nonlinear resistance [START_REF] Hongyi | Extended black-box model of pantograph arcing considering varying pantograph detachment distance[END_REF]. Therefore, a variable resistance (% &'( ) represents connection or disconnection of pantograph and catenary including arcing phenomenon, in the circuit of figure 1. Table 1 shows all the parameters of circuit in figure 1. In order to have a model of the whole power chain (Figure 1) bond graph method is chosen. This method illustrates energy transfers in systems by considering power links with two energy variables: generalized effort (voltage) and generalized flow (current) [START_REF] Gandanegara | Modeling and Multi-time Scale Analysis of Railway Traction Systems Using Bond Graphs[END_REF]. All physical phenomena, including integral causality, are modeled by some elements linked by junctions [START_REF] Karnopp | System dynamics: a unified approach[END_REF]. Figure 2 depicts the bond graph of whole system shown in Figure 1. It can be deduced that the inputs of the system are ! ) , ! "# , ! "$ .

Parameter

Based on bond graph definitions, all the components of the power chain are modelled using elements and junctions as described in table 2:

Parameter Definitions 2 ) 3
Effort source representing voltage source 3 % Resistance Bond graph directly allows identifying inputs of the system (generalized effort of effort source) but also state variables (generalized flow of I elements and effort of C elements). Consequently, following junction equations, state equations of the whole system can be directly achieved considering energy variables and causality:

; = # F GHIJ 4 (&. @AB C # F GHIJ 4 / @AB!!! 4 / ; = K GHIJ @.B * L C M> HNG @.BO> L P * L 4 / @AB C K ? @.B * L !!!!!!!(1)! 4 0 ; @AB = K ? @.B * Q !!
4 "# ; = @<> R S T RS @UBO0 S K ? @UBV<K RS @UBB * RS !! 4 "$ ; = @<> R J T RJ @UBO0 J K ? @UBV<K RJ @UBB * RJ !! ! * @AB = % 1 @4 W @AB C X # 4 "# @AB C X $ 4 "$ @AB C 4 0 B!!!

(1)
where m 1 and m 2 are primary to secondary ratios of transformer. State equations are useful for simulation purposes due to the possibility to solve the system along time using numerical integration algorithms such as Euler, Heun or Runge-Kutta methods [START_REF] Süli | An Introduction to Numerical Analysis[END_REF][START_REF] Press | Numerical Recipes in C, The art of scientific computing[END_REF].

Arc Model

Generally, there are four different conditions during pantograph bounce: no arc, arc drawn, arc extinguished and arc approach. Arc drawn happens while distance between pantograph and catenary is sufficiently small and increases. In this case, due to the existence of arc, there is still electrical connection between pantograph and catenary. While this disconnection length increases until there is no enough energy to maintain the arc, it will be extinguished and the electrical connection of pantograph and catenary will be loosed. Following these steps, pantograph goes back to the catenary. While distance between pantograph and catenary decreases, electrical field increases. So that, arc phenomenon will occur gain. This step is called arc approaching. Note that the maximum arc length in this step is less than the one in the step of arc drawn [START_REF] Wang | The pantograph-catenary arc test system for high-speed railways[END_REF]. Figure 3 shows all the steps of arc drawn, arc extinguished and arc approaching. During this process, variation of disconnection length along time is considered as a constant [START_REF] Hongyi | Extended black-box model of pantograph arcing considering varying pantograph detachment distance[END_REF]. Note that some of these three steps could be neglected due to different case of disconnection on reality. Indeed, arc extinguishing may not appear if the arc length does not exceed its maximum value.

As assumed before, electrical arc is modelled through a nonlinear resistance R arc that depends on the relative position of catenary and pantograph. Obviously, % &'( = 8 if there is a contact between catenary and pantograph and no arc is created. Moreover, % &'( = Y if there is disconnection between catenary and pantograph without arc phenomenon (arc extinguished). In this paper, Habedank's equation is used to describe arc phenomenon. Indeed, in case of train study, it allows to take into account effects of airflow speed that can be assumed to be similar to train speed [START_REF] Habedank | On the mathematical-description of arc behavior in the vicinity of current zero[END_REF]. Moreover, the model describes variations of length between pantograph and catenary denoted in the following + &'( @AB [START_REF] Cassie | An Improved arc model before current zero based on the combined Mayr and Cassie arc models[END_REF]. Equation [START_REF] Hongyi | Extended black-box model of pantograph arcing considering varying pantograph detachment distance[END_REF] gives arc resistance along time when there is disconnection between catenary and pantograph with arc phenomenon depending on g c (t) and g m (t), respectively Cassie and Mayr conductance. Thus, Habedank's model is a series association of Cassie's and Mayr's models. Note that Cassie's model is more adapted for high current values and Mayr's model for arc with low current values [START_REF] Yuan | Simulation of Fault Arc Using Conventional Arc Models[END_REF][START_REF] Tseng | An Experimentally Verified Hybrid Cassie-Mayr Electric Arc Model for Power Electronics Simulations[END_REF].

% &'( @AB = # Z HNG @.B = # Z G @.B D # Z Q @.B (2)
Conductance's g c (t) and g m (t) are governed by nonlinear differential equations ( 3) and ( 4)

[Z G [. = # \ J ] T L V@.B J 9 G @.B J Z G @.B C ^(@AB_VVVVVVVVVVV (3)
where `(@AB = -+ &'( @AB,

[Z Q [. = # \ S ] T L @.B J W a @.B C ^0@AB_ (4) with:
-Vb c @AB = d^& '( @AB e + &'( @AB f if g h i c

-b c @AB = j+ &'( @ABk # @g D lBm4 /,'0" @AB if g n i c 4 /,'0" @AB is the RMS value over a period of 4 / @AB with 4 /,'0" @AB = o # p q r 4 / $ @`Bs.

.<p t, u, Vv # , Vv $ , g, Vk # , d,are respectively constant arc dissipation power, parameter depending on arc phenomenon, time constant, airflow speed, correlation coefficient of dissipated power for arc, pyroelectric coefficient and ratio between + &'( @AB and `(. V 0 is a given airflow speed.

For the parameters k and β, the ranges of their values are k>0 and 0<t<3, according to measured data in the pantograph arcing test system developed in [START_REF] Wang | An Extended Habedank's Equation-Based EMTP Model of Pantograph Arcing Considering Pantograph-Catenary Interactions and Train Speeds[END_REF], [START_REF] Gao | Study on the mechanism and suppression method of transient overvoltage in the state of fast train[END_REF].

It can be assumed that the minimum values of conductance, i.e. the maximum values of arc resistance occur for I p (t) close to 0. Thus, voltage spikes appear on primary of transformer for zero crossing current. Considering the elements in the power model, such as inductance and capacitances that change the phase shift of current have to be precisely known at each time in order to have accurate values of voltage. Finally, voltage spikes time location depends on train position regarding substation, catenary parameters and also electromechanical chain that modify phase shift of current on the secondary part of the transformer.

Simulation results

In this section, some simulation results will be detailed. First of all, numerical integration algorithm will be detailed. Then, for a given set of parameters, state variables will be studied in case of pantograph approach with arc phenomenon.

Numerical integration

The solving method chosen for differential equations is Heun method [START_REF] Süli | An Introduction to Numerical Analysis[END_REF]. Let's consider differential equation in [START_REF] Cassie | An Improved arc model before current zero based on the combined Mayr and Cassie arc models[END_REF]. From a time instant j, an approximation of state X at time instant j+1 is given by ( 6) with dt the time step between the two time instants.

w ; = x@w, !B (5) where X is the vector of state variables and U the vector of inputs.

d # = sAq xMw y , ! y P, ! d $ = sAq xMw y D d # , ! yO# P (6)
w yO# = w y D 8qz { @d # D d $ B

State variables

Based on the explanations in section 3.1, the whole circuit of figure 1 is simulated using Heun method. Table 3 shows the value of parameters of circuit that are used in simulation. Figure 4 shows some of the state variable signals that were obtained by simulation. Arc parameters values were extracted from [START_REF] Hongyi | Extended black-box model of pantograph arcing considering varying pantograph detachment distance[END_REF]. It can be observed that there is an amplitude decrease of the current signals (Catenary and primary of transformer) due to the arc phenomenon. There are no changes in the catenary voltage. Indeed, catenary capacitor may be seen as a voltage source, thus its voltage amplitude may not change during arc. However, as global resistance of the system increases due to arc phenomenon, at constant voltage, catenary and transformer currents may decrease. Some effects also appear on the derivative of currents.

Parameter

+ (&. 1 mH ! "# V 920 V % / 0.14 W ! "$ 920 V

Experimental validation

In order to be the most accurate as possible, simulations are performed with similar conditions to tests on train. The detachment of pantograph and catenary was done very fast and no arc was created. Then pantograph is reattached to catenary slowly and the arc phenomenon can be observed. Therefore, all the steps of disconnection between pantograph and catenary are neglected except arc approach. In the step of simulation, disconnection length is reducing linearly as a function of time. During the test, substation voltage is 15 kV (RMS) and frequency of 7~, • € • Hz. Train speed is 100 km/h. Therefore, we have g ' i c .

Based on the investigation of [START_REF] Hongyi | Extended black-box model of pantograph arcing considering varying pantograph detachment distance[END_REF], equation of L arc (t) for simulation part will be as below:

+ &'( @AB = C + V0&ƒ "A @A C A $ B D + V0&ƒ (7)
where:

-L max =5cm the maximum distance between pantograph and catenary, -Dt=0.5s the approach duration, -A $ = €7…q8zs is the start time of pantograph approach. For the parameters of equation ( 3), based on the investigation [START_REF] Press | Numerical Recipes in C, The art of scientific computing[END_REF], the arc waveform is close to experimental results with a=0.81, b=10, K 1 =2, v # = v $ =0.1 ms and C=15.

All the arc parameters are identified by internal software based on experimental results. Due to the confidential aspects, it is not possible to go in detail about method of identification. Figure 5 shows the comparison between the train result and simulation result based on the equations ( 1) and (2). It can be seen that effects of arc phenomenon on primary voltage of transformer in simulation is representative of real measurements on train. As explained in section 2, there are two voltage spikes per period due to the zero crossing point of current. Zero crossing point of current are perfectly predicted thanks to modeling of the whole power chain. Due to measurements parameters such as sample time, the comparison of spikes high frequency oscillations are difficult. However, this is balanced by the fact that the worst effect of arcing lies in the values of spikes that can induce overvoltage and then damages in the traction chain.

In this paper, pantograph bounce with arc phenomenon was modeled inside a complex power traction chain using bond graph. This model is achieved to lead to state equations that can be easily simulated. The major effects of pantograph bounce are on the input of traction chain: primary voltage and current of transformer. Comparison between real measurements and simulation results shows that presented model is sufficiently accurate to predict the behavior of the whole system during pantograph bounce phenomenon. In future work, arc behavior regarding parameters of the power chain such as inductance and resistance of catenary and primary of transformer will be studied. A sensitivity analysis will be performed to clearly understand the effect of train distance and speed regarding substation.

Acknowledge

This research was supported partially by national association research technology (ANRT).

Fig. 1 :

 1 Fig. 1: Power supply, catenary and traction chain including arc resistance.

  Generalized efforts are identical on each power link attached to the junction

Fig. 2 :

 2 Fig. 2: Bond graph of power supply, catenary and traction chain including R arc

Fig. 3 :

 3 Fig. 3: Pantograph bounces process

Fig. 4 :

 4 Fig. 4: Simulation result of state variable: a) Catenary voltage, b) Current of catenary, c) Primary current of transformer.

Fig. 5 :

 5 Fig. 5: Primary voltage of transformer on train and simulation in presence of arc phenomenon

Table 1 :

 1 Parameters of the circuit

	! ) % * + * -(&.# , -(&.$ % (&. + (&. % &'( % / + / + 0 % 1 % "# , % "$ + "# , + "$	Definitions Voltage of substation Resistance of substation Inductance of substation Capacitors of catenary Resistance of catenary Inductance of catenary Nonlinear resistance of arc Leakage resistance of primary of input train transformer Leakage inductance of primary of input train transformer Magnetizing inductance of input train transformer Parallel magnetizing resistance of input train transformer Resistances of secondary of input train transformer Leakage inductances of secondary of input train
	! "# , ! "$	transformer Equivalent voltage sources for electromechanical traction chain

Table 2 :

 2 Bond graph elements definitions used in figure3

Table 3 :

 3 Values of circuit parameters used in simulation