Saba Amirdehi
email: saba.amirdehi@yahoo.com

Alstom Transport

Johana Vally
email: johana.vally@alstomgroup.com

Baptiste Trajin
email: baptiste.trajin@enit.fr

Didier Colin
email: didier.colin@alstomgroup.com

Paul-Etienne Vidal
email: paul-etienne.vidal@enit.fr

Sylvain Tastet
email: sylvain.tastet@alstomgroup.com

Identification of PEM Fuel cell system model for traction chain based on profile mission measurement

Keywords: PEMFC system model, polarization curve, train traction chain, parameter identification

Testing of proton exchange membrane fuel cell system is one of the most important issues due to the cost and security reasons in railway applications. To reduce the cost, it would be necessary to have a representative simulation model for this case. Having an accurate model of proton exchange membrane fuel cell system requires considering several auxiliary circuits. However, in industrial application, the accurate model of proton exchange membrane fuel cell system and corresponding information are not often provided by the designer. Moreover, the simplicity and representativity of model is required to be compatible with overall traction chain model. The model has to be suitable for system simulation and not design purpose. In this paper, a simplified model of PEMFC system is proposed using the experimental profile mission. The impact of all the environment system of proton exchange membrane fuel cell stack is estimated using the average of polarization curve on several operational points. The parameter identification process is applied, and the obtained parameters are then inserted on the overall system of fuel cell stack. The experimental data obtained through train is compared with the consecutive simulation results.

I. INTRODUCTION

Fuel cells (FCs) are one of the most promising renewable energy sources. They are environmentally clean, and they have low emission of oxides of nitrogen and sulfur. There are wide range of FC technologies classified into different categories based on fuel used, operating temperature and type of electrolyte. Among different types of FCs, Proton Exchange Membrane Fuel Cell (PEMFC) has lower temperature (under 100°C), fast connection time and easiest manufacturing process compare to the other types of FC. It has fast dynamic response in the condition of load variation, and it is compact, lightweight and noiseless [START_REF] Blunier | Pile Combustible (French)[END_REF].

The open circuit voltage of PEMFC is obtained by the differential potential between oxidant (Oxygen coming from air) and fuel (Hydrogen) in cathode and anode, respectively [START_REF] Dicks | [END_REF]. The maximum differential potential provided by PEMFC is about 1.23 V. In railway applications, to have sufficient amount of power (around 250 KW) to supply the whole traction chain, PEM fuel cell system is used. It consists of several PEMFC stacks in different configuration in series, parallel or combination of series and parallel. To manage an adequate gas, it is required to have several auxiliaries such as hydrogen, cooling, air and humidification circuit.

The testing activities are the most important issue due to security and cost reasons in railway traction chain design. Therefore, having a representative simulation model of chain components is required. The simulated model of components such as PEMFC system has to be compatible with overall traction chain model. Thus, the simplicity and representativity of PEMFC system model are in priority. The model has to be suitable for system simulation and not design simulation. The PEMFC system model has to consist of the PEMFC stack model and the auxiliary circuits.

Several approaches are done to propose the static and dynamic model of PEMFC stacks [START_REF] Lee | Dynamic and Steady-State Performance of PEM Fuel Cells under Various Loading Conditions[END_REF]- [START_REF] Cultura | Dynamic Analysis of a Stand Alone Operation of PEM Fuel Cell System[END_REF] and [START_REF] Becherif | Electrical equivalent model of a proton exchange membrane fuel cell with[END_REF] which provide the polarization curve (voltage as a function of current). In last few years, the numerical modeling, computer simulation and PEMFC stacks performance have been proposed [START_REF] Dicks | [END_REF], [START_REF] Maher | Modelling of proton exchange membrane fuel cell performance based on semi-empirical equations[END_REF]- [START_REF] Khan | Dynamic Modelling and Simulation of a Fuel Cell Generator[END_REF]. State-space modelling of PEMFC is provided in [START_REF] Sachin | State-Space Modeling of Proton Exchange[END_REF]. These approaches allow us to establish equivalent electrical model of PEMFC. There are several approaches to propose the PEMFC stacks model using graphic way of power transfer in the system. In [START_REF] Bressel | Dynamical modeling of Proton Exchange Membrane Fuel Cell and parameters[END_REF]- [START_REF] Bressel | Bond Graph for modeling and diagnostics of Proton Exchange Membrane Fuel Cell[END_REF], PEMFC stack model is proposed based on bond graph methodology [START_REF] Gandanegara | Modeling and Multi-time Scale Analysis of Railway Traction Systems Using Bond Graphs[END_REF] which illustrates the energy exchange within overall system. Generally, the methodology of bond graph is applied to model the PEMFC stack in the following sections.

As mentioned in previous approaches, the polarization curve depends on numerous factors such as membrane thickness and state of hydration, catalyst loading, catalyst layer structure flow filed design and operating conditions (temperature, pressure, humidity, flowrates and etc.) [START_REF] Mahmah | Dynamic Performance of Fuel Cell Power Module for Mobility Applications[END_REF]. Experimentally, the polarization curve is also impacted by the transient load use, which causes a hysteresis effect on the dynamic characteristics of PEMFC stack [START_REF] Talj | Parameter analysis of PEM fuel cell hysteresis effects for transient load use[END_REF]. Based on several studies, the hysteresis effect occurs due to several reason such as double layer charging effect [START_REF] Fontes | Interactions Between Fuel Cells and Power Converters: Influence of Current Harmonics on a Fuel Cell Stack[END_REF] and its interaction with activation losses [START_REF] Gerteisen | Modeling the phenomena of dehydration and flooding of a polymer electrolyte membrane fuel cell[END_REF] and water management in the cells [START_REF] Le Canut | Impedance study of membrane dehydration and compression in proton exchange membrane fuel cells[END_REF]. Therefore, it would be necessary to take into account the PEMFC stack dynamic model, auxiliary circuits model, its control command system to realize the polarization curve of overall system in various operating modes. The detailed model of PEMFC system are presented in [START_REF] Li | Modeling and performance analysis of PEM fuel cell system[END_REF] and [START_REF] Eng | Modeling and Analysis of Renewable PEM Fuel Cell System[END_REF].

However, in railway application, a simplified model from user point of view is required to test the overall performance of system in different operating modes. Moreover, the detailed information and data about the PEMFC system are not often provided by the manufactures. Therefore, it is required to propose a model based on the overall profile mission of PEMFC system.

In this paper, a simplified model of PEMFC system is proposed using the experimental profile mission on train. The PEMFC stack model is obtained using bond graph methodology. To adapt the PEMFC stacks model parameters, the output polarization curves on experimental measurements affected by several auxiliary circuits is estimated by a specific average curve. This paper is organized as following. In section II, the general structure of PEMFC stack and PEMFC system are explained. The electrical model considering voltage drops of each cell under load condition is illustrated. The corresponding equations are provided. In section III, the parameter identification process is illustrated. The average polarization curve is obtained using the part of measured profile mission. The parameters of electrical circuit are then adapted using the average characteristics of PEMFC system. The railway experimental results in the other part of mission profile are compared with simulation results in section IV.

II. PEM FUEL CELL SYSTEM

A. PEMFC System Structure

A PEMFC system is constituted by several PEMFC stacks in series/parallel combination to achieve higher output voltage and power. Fig. 1 shows the example of basic PEMFC stack and PEMFC system configuration. As it is illustrated in Fig. 1 (a), the PEMFC stacks consists of a multitude of single cells stacked up. Each single cell is constituted by a membrane, electrodes which are the catalyst layer, a gas diffusion layer and a bipolar plate to separate the cells connected in series. The electrochemical reactions occur on anode (hydrogen) and cathode (oxygen). Therefore, in the PEMFC stack structure, the cathode of a single cell is electrically connected to the anode of the adjacent cell. In this paper, it is assumed that the parameters of individual cells can be lumped together to represent PEMFC stacks. Thus, it is assumed that the same current flows through each cell on the stack. Considering the overall current of PEMFC system as ܫ ௨௧ ሺݐሻ, the current for each single cell is as below:

ܫ ሺݐሻ ൌ ܫ ௨௧ ሺݐሻ ܰ ௦௧ǡ (1)
where ܰ ௦௧ǡ is number of series stacks in parallel. Supposing the output voltage response of each single cell as ܸ ሺݐሻ , the output voltage of PEMFC system can be calculated as [START_REF] Dicks | [END_REF].

ܸ ௨௧ ሺݐሻ ൌ ቀܰ ܰ ௦௧ǡ௦௦ ܸ ሺݐሻቁ െ ܸ ሺݐሻ െ ܸ ோ ಽ ሺݐሻ (2)
where ܰ and ܰ ௦௧ǡ௦௦ are numbers of single cell in PEMFC stacks, number of stacks in series, respectively. ܸ ሺݐሻ and ܸ ோ ಽ ሺݐሻ represent the voltage losses due to cables' impedance which can have an impact on dynamic response.

PEMFC stacks are obviously the heart of PEMFC system. However, it is required the supporting equipment to be useful and perform efficiently. Typically, the PEMFC systems involves several subsystems such as oxidant supply, fuel supply, heat management, water management, power conditioning and instrumentation and controls. Fig. 2 depicts a PEMFC stack with the involved auxiliary circuits.

Fig. 2 PEMFC stack and the involved auxiliary circuits

There is also a control system which sends or receives the commands depending on the operating modes. The control system structure is different depending on the designers.

The static and dynamic behavior of the system output are impacted by the subsystems and PEMFC stacks. However, as it is required to have a simple and representative model of PEMFC system, it is focused on the user point of view model. To do so, the impact of auxiliary circuits on polarization curve is considered in the proposed methodology in this paper considering the average polarization curve on train measurements. In the following sections, the PEMFC stack model is presented and the proposed methodology to predict the behavior of PEMFC system is discussed.

B. PEM Fuel cell stack model

The electrical model of a single cell of PEMFC stack is proposed in different literatures [START_REF] Wang | Dynamic models and model validation for PEM fuel cells using electrical circuits[END_REF]- [START_REF] Becherif | Electrical equivalent model of a proton exchange membrane fuel cell with[END_REF]. Fig. 3 depicts the PEMFC dynamic model considered. Once the electrical load is connected to the fuel cell, the internal potential will decrease due to the voltage drops caused by the kinetic electrochemical reactions, convection and diffusion process, resistance of electrons flow, and transferred ions through solid polymer membrane and electrode materials. These voltage drops are namely, activation, concentration and ohmic losses, respectively [START_REF] Wang | Dynamic models and model validation for PEM fuel cells using electrical circuits[END_REF]. There is also the double layer charging effect which is presented as a capacitor denoted C [START_REF] Dicks | [END_REF]. The voltage equations of internal potential and voltage losses in the presence of load are presented as ( 3)- [START_REF] Becherif | Electrical equivalent model of a proton exchange membrane fuel cell with[END_REF].

ܧ ሺݐሻ ൌ ܧ ǡ ݈݊ሺ ுଶ כ Ǥ ሺ ைଶ כ ሻ Ǥହ ሻ െ ܧ ௗǡ (3)
where ܧ ǡ ,  ுଶ כ ,  ைଶ כ and ܧ ௗǡ represent standard reference potential at standard state (298K and 1-atm pressure) which is 1.23 V, partial pressures (hydrogen and oxygen) and the fuel and oxidant delay of PEMFC, respectively. The fuel and oxidant delay are obtained as below [START_REF] Wang | Dynamic models and model validation for PEM fuel cells using electrical circuits[END_REF]:

ܧ ௗǡ ሺݐሻ ൌ ߣ ൬ܫ ሺݐሻ െ ܫ ሺݐሻ כ ݔ݁ ൬െ ݐ ߬ ൰൰ (4)
where ߣ and ߬ are constant value ( π ) and overall flow delay in seconds, respectively. Note that ( כ ) is the convolution product. ܧ ௗǡ ሺݐሻ is zero in steady state.

The ohmic voltage loss and state equation which represents the double layer charging effect are provided in ( 5) and ( 6) respectively.

ܸ ǡ ሺݐሻ ൌ ܫ ሺݐሻ ݈ ܵሺߪሻ (5) ܸ ሶ ǡ ሺݐሻ ൌ ͳ ܥ ܫ ሺݐሻ െ ͳ ൫ܴ ௧ǡ ܴ ǡ ൯ܥ ܸ ǡ ሺݐሻ ( 6
)
where ܵ is cell active area in cm 2 and ݈ is the thickness of membrane in cm (electrolyte of cell). ߪ is the membrane conductivity in mS/cm.

In [START_REF] Prasadsahu | Performance Study of PEM Fuel Cell under Different Loading Conditions[END_REF], ܥ represent double layer charging effect. ܴ ௧ǡ and ܴ ǡ are activation and concentration losses of each single cell of PEMFC, respectively. They can be expressed as [START_REF] Cultura | Dynamic Analysis of a Stand Alone Operation of PEM Fuel Cell System[END_REF] and [START_REF] Becherif | Electrical equivalent model of a proton exchange membrane fuel cell with[END_REF].

ܴ ௧ǡ ሺݐሻ ൌ ܴܶ ܫ ሺݐሻߙ݊ܨ ݈݊ ൬ ܫ ሺݐሻ ܫ ܫ ൰ (7)
where ݊ ൌ ʹǡ ߙǡ ܫ and ܫ are the number of transferred electrons in chemical reactions in PEMFC, the charge transfer coefficient, corresponding current of activation and the internal leakage current, respectively. ߙ, ܫ and ܫ are the parameters which should be found from datasheet and experimentally.

ܴ ǡ ሺݐሻ ൌ െ ܴܶ ܫ ሺݐሻߚ݊ܨ ݈݊ ൬ͳ െ ܫ ሺݐሻ ܫ ௧ ൰ (8)
where ܫ ௧ and ߚ are the limit diffusion current and a transfer coefficient of concentration loss, respectively. Note that ( 7) and ( 8) are valid while ܫ ሺݐሻ ് Ͳ.

Therefore, the dynamic output voltage of PEMFC ( 9) can be obtained using equations of ( 3)- [START_REF] Becherif | Electrical equivalent model of a proton exchange membrane fuel cell with[END_REF].

ܸ ሺݐሻ ൌ ܧ ሺݐሻ െ ܸ ǡ ሺݐሻ െ ܸ ǡ ሺݐሻ (9)
The output voltage of PEMFC single cell in steady state condition (without load transition) turns to [START_REF] Fmoun | Electrochemical Circuit Model of a PEM Fuel Cell[END_REF].

ܸ ሺݐሻ ൌ ܧ ሺݐሻ െ ܸ ௧ǡ ሺݐሻ െ ܸ ǡ ሺݐሻ െ ܸ ǡ ሺݐሻ (10)
where ܸ ௧ǡ ሺݐሻ and ܸ ǡ ሺݐሻ are the linked to ܴ ௧ǡ ሺݐሻ and ܴ ǡ ሺݐሻ respectively. The static polarization curve of PEMFC stack is as Fig. 4. Since the fuel cell stack developed for railway application are commercial products, authors have no access of detailed information about stack characteristics; and have taken some hypothesis for ܵ and ݈ constant based on the state of the art in the following sections. For the same reason, the control command of the auxiliary system is considered as a black box. Therefore, the dynamic behavior of PEMFC system is estimated using the profile mission measured on train. The estimated polarization curve is presented during the identification process in the following.

III. PARAMETER IDENTIFICATION

A. Experimental Profile Mission analysis

The parameter identification of a system is done by comparing the simulation results and system experimental measurement and tuning the model parameters. The optimal set of parameters is determined by minimizing the difference between the system response and simulation results.

To have the experimental results, the measurements are done on train (hybrid railway). Fig. 5 demonstrates the block diagram of the hybrid train including PEMFC system, battery, converter, motor, etc. It is also shown where the voltage and current measurements are done. The PEMFC system detailed block diagram is presented in Fig. 1 and2. Generally, the static polarization curve presented in the datasheet of PEMFC stacks are used to identify the parameters. Then, the dynamic response of PEMFC system which experimentally shows the nonlinear behavior, is obtained by a precise model of the overall PEMFC system. However, as it is aimed to have simplified model in user point of view, following approaches are proposed.

As depicted in Fig. 6, various voltage-current curves are observed based on the power demand to ensure the railway profile mission. All experimental polarization curves in profile mission are shown in Fig. 7 within 3h of measurement. Let us remind that the voltage and current shown in the following sections and figures are united by the maximum value of PEMFC voltage and current provided in the datasheet. The complete polarization curve can be obtained while progressive power range is demanded. Within 3h of measurement, 21 curves with the time constant of more than 5 seconds are registered. As illustrated in Fig. 7, it is distinguished two modes, denoted Mode1 and Mode 2. The Mode 1 defines the typical polarization curve impacted by the dynamic feature of system for different operating points considering the variation of temperature, level of the hydration etc. The voltage variation for the same demanded power is around 7%. This is due to the factors mentioned in section I which impact the polarization curve of PEMFC system. However, to identify the parameters, it is required to have unique polarization curve which represents a typical behavior of the PEMFC system during profile mission. This typical polarization curve allows us to identify the static parameters of system. Therefore, an average of polarization curve is obtained in next section.

The Mode 2 is due to the dynamic response of PEMFC system in fast power variation (<5 seconds response) which is out of the usual polarization curve. This mode occurs at maximum power while dramatically power reduction is demanded. At maximum power demand the hydration level of membrane is high and does not respond instantaneously. This mode allows us to identify the dynamic parameters of system.

B. Paramater identification of static model

Mainly, to identify the parameters, a part of measured voltage and current (last 2h of test procedure) is considered. Then, the identified parameters are inserted in the model and the simulation results are compared with the other part of measurement (1h in the beginning of test procedure) which is not the same as the part used in identification process to confirm the model.

The differential evolution (DE) method [START_REF] Aluffi-Pentini | Global Optimization and Stochastic Differential Equation[END_REF] is used here to identify the parameters. In this case the cost function for DE method is achieved by minimizing the difference between desired and calculated output voltage of PEMFC system. The cost function is as: [START_REF] Khan | Dynamic Modelling and Simulation of a Fuel Cell Generator[END_REF] where ߝ represents the sum of error, ሺܸ ݐݑ ሻ ǡௌ is the calculated output voltage of PEMFC system in steady state operation [START_REF] Fmoun | Electrochemical Circuit Model of a PEM Fuel Cell[END_REF]. ሺܸ ݐݑ ሻ ǡெ is the measured output voltage of PEMFC system in the last 2h of test procedure. To have a unique voltage response, the characteristics of PEMFC system in Mode 1 within the last 2h of measurements can be considered as an average polarization curve. The average of polarization curve is obtained using a window of 6% maximum measured current as Fig. 8. For each window, the average value voltage and current values are represented by (ܸ ௨௧ǡ , ܫ ௨௧ǡ ) operating points which are plotted by the red curve in Fig. 8.

Mode 1 Mode 2 ߝ ൌ ൫ሺܸ ௨௧ ሻ ǡெ െ ሺܸ ௨௧ ሻ ǡௌ ൯ (

C. Paramater identification of dynamic model

Once the desired parameters are obtained, the identification process has to be redone for the dynamic model. In this case one cycle of response voltage in Mode 2 is selected; the second identification process provides the values of ߣ , ߬ and ܥ . The cable impedance value is considered based on the its nominal value per meter.

IV. SIMULATION RESULTS AND COMPARISON WITH EXPERIMENTAL RESULTS

The PEMFC system model is simulated using (1)- [START_REF] Fmoun | Electrochemical Circuit Model of a PEM Fuel Cell[END_REF]. All the identified parameters in section III are inserted in the model. The measured power demand of PEMFC system in the first 1h of measurement is applied to the simulated model; then the results are compared with the measured output voltage as illustrated in Fig. 9. The average error of model during the profile mission in the first 1h of measurement is 1.85%. In mode 1 which consists the 90% of profile mission, the average value of error is about 1.7% and in mode 2 which consists the 10% of profile mission is about 5%. On the proposed model, the average of profile mission is considered to identify the parameters, therefore it is normal to have error on some specific and unusual point of profile mission (mode 2). The polarization curve obtained in simulation is compared with result on train in the first 1h of testing in Fig. 10. It is shown that simplified model of PEMFC system can represent two modes of experimental polarization curve. It is concluded that identifying the model parameter of PEMFC stack using profile mission allows to have a simple and representative model of PEMFC system with enough accuracy.

V. CONCLUSION

In this paper, PEMFC system is modeled using average experimental profile mission on train. This leads to have a simplified model of overall system in order to be easily simulate in railway application. The parameters of model are identified using the experimental average polarization curve. The two major characteristics of system in various operating modes are simulated using the dynamic model. Comparison between experimental measurements and simulation results depicts that the presented model sufficiently accurate to predict the behavior of PEMFC system.

In future work, the identification of PEMFC system model will be done considering the impact of control command circuits of PEMFC system.

ACKNOWLEDGMENT

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

… …

funding from the Shift2Rail Joint Undertaking (JU) under grant agreement No 826054. The JU receives support from the European Union's Horizon 2020 research and innovation program and the Shift2Rail JU members other than the Union.

  a) PEMFC stack (b) PEMFC system

Fig. 3

 3 Fig. 3 Dynamic model of a single cell of PEMFC In Fig.3, ܧ is the internal voltage of PEMFC provided by chemical and mass reactions in open circuit condition.Once the electrical load is connected to the fuel cell, the internal potential will decrease due to the voltage drops caused by the kinetic electrochemical reactions, convection and diffusion process, resistance of electrons flow, and transferred ions through solid polymer membrane and electrode materials. These voltage drops are namely, activation, concentration and ohmic losses, respectively[START_REF] Wang | Dynamic models and model validation for PEM fuel cells using electrical circuits[END_REF]. There is also the double layer charging effect which is presented as a capacitor denoted C[START_REF] Dicks | [END_REF]. The voltage equations of internal potential and voltage losses in the presence of load are presented as (3)-(8).

Fig. 4

 4 Fig. 4 Characteristics of PEMFC single cell in steady state

Fig. 5

 5 Fig. 5 Block diagram of train including fuel cell under measurement

Fig. 6

 6 Fig.6 depicts the measured current (a) and voltage of PEMFC system (b) on the part of profile mission on train united by the maximum value of voltage and current ( ೠ ೌೣ , ூ ೠ ூ ೌೣ ) of PEMFC system provided in datasheet. Note that the sampling time of measurement is 0.5 second.

Fig. 7

 7 Fig. 7 Measured voltage response versus the current extracted from profile mission for 3h of test procedure

Fig. 8

 8 Fig. 8 Average voltage-current points for profile mission in last 2h of measurements The average curve represents the simplified polarization curve of PEMFC system used to identify the parameters of the PEMFC stack steady state model. The objective is to reproduce the response voltage of PEMFC system which is measured on 1h in the beginning of the test procedure shown in Fig.6 (b). The optimization is started with 100 for number of population and 100 for number of evaluations. As a matter of facts, 10000 numbers of iterations are applied. The required parameters for steady state model of PEMFC stack are ܫ ǡ ܫ and ܫ ௧ . ߙǡ ߚ and ߪ are patent technological parameters given in PEMFC stack datasheet.

Fig. 9

 9 Fig. 9 Comparison of measured output voltage of PEMFC system and the simulated one using identified parameters in section III

Fig. 10

 10 Fig. 10 Comparison of measured voltage and current of PEMFC system on train in 1h in the beginning of test and the simulated one.


Roxane Mallouhy
email: reliasmallouhy@pmu.edu.sa

Chady Abou Jaoude
email: chady.aboujaoude@ua.edu.lb

Christophe Guyeux
email: christophe.guyeux@univ-fcomte.fr

Abdallah Makhoul
email: abdallah.makhoul@univ-fcomte.fr

Major earthquake event prediction using various machine learning algorithms

Keywords:

At least two basic categories of earthquake prediction exist: short-term predictions and forecast ones. Short term earthquake predictions are made hours or days in advance, while forecasts are predicted months to years in advance. The majority of studies are done on forecast, taking into consideration the history of earthquakes in specific countries and areas. In this context, the core idea of this work is to predict whereas an event is classified as negative or positive major earthquake by applying different machine learning algorithms. Eight different algorithms have been applied on a real earthquake dataset, namely: Random Forest, Naive Bayes, Logistic Regression, MultiLayer Perceptron, AdaBoost, K-nearest neighbors, Support Vector Machine, and Classification and Regression Trees. For each selected model, various hyperparameters have been selected, and obtained prediction results have been fairly compared using various metrics, leading to a reliable prediction of major events for 3 of them.

I. INTRODUCTION

Rapid global urbanization leads to an increase in earthquakes, a major catastrophe which has occupied 60% of deaths overall natural disasters [START_REF] Dong | A comprehensive review of earthquakeinduced building damage detection with remote sensing techniques[END_REF]. Obviously, human cannot stop natural disasters, while the application of machine learning is a powerful and invaluable method and technique, used by researchers as new area of study in geology, in order to reduce, as much as they can, the loss of life and billions of dollars in infrastructure and housing costs. Traditional classification methods rely on statistical assumptions for earthquakes that turn out to be unsatisfactory in danger state prediction. Therefore, the use of machine learning algorithms is becoming more widespread, that adapt and learn a problem by simulating a biological or natural system.

Introducing the machine learning techniques in the field of geology and earthquake beats the traditional and standard methods used in earlier years and provides the scientists a new fresh method for assessing seismic risk and triggering future earthquakes. Data mining applications have achieved great success in geophysics and geology such as ecology, weather prediction, modeling, etc. and finally predicting the earthquakes. Two basic categories of earthquake prediction exist: short-term predictions and forecast ones. Short term earthquake predictions are made hours or days in advance, while forecasts are predicted months to years in advance. Due to the chaotic and complex phenomena of earthquake process, it is considered that short-term predictions are very difficult to address. Thus, the majority of studies are done on forecast, taking into consideration the history of earthquakes in specific countries and areas.

The core idea of this work is to predict whereas an event is classified as negative or positive major earthquake by applying different machine learning algorithms. It is well known that there is no best algorithm or one solution that fits all the problems and datasets for machine learning since the performance of algorithms depends on many factors. Some algorithms are best for small data, while others perform better for a tons of data sample. Some algorithms require categorical inputs, while others need quantitative. Another important criterion while choosing the algorithm is the complexity of the dataset and how many features the model needs to learn and predict. This is why, in this work, eight different algorithms have been applied on an earthquake dataset, namely: Random Forest, Naive Bayes, Logistic Regression, MultiLayer Perceptron, AdaBoost, K-nearest neighbors, Support Vector Machine, and Classification and Regression Trees. For each selected model, various hyperparameters have been tested, and obtained prediction results have been fairly compared using various metrics, leading to a reliable prediction of major events for 3 of them.

The remainder of this paper is organized as follows. The next section contains the basic recalls. Section III discusses related work done by researchers in the field of earthquake machine learning. In Section IV, we apply experiments on the testing set after training the data using machine learning algorithms and for some models, parameters were modified, then results were compared. Finally, Section V concludes the paper.

II. BASIC RECALLS

We begin by recalling the various classical machine learning methods that have been proved useful in predicting earthquakes, before detailing, in the next section, some state-of-the-art articles that have used them in this context.

A. Naive Bayes

The Naive Bayes algorithm, a method based on the work of Thomas Bayes, is a heavily simplified probabilistic model that calculates probabilities by counting combinations of values and frequency in a data set. Bayesian classification assumes that the data belongs to a particular class, then the probability for the hypothesis to be true is calculated. Naive Bayes operates on a strong independence assumption, in other terms, the conditional probability of one attributes will not affect the probability of the other. However, the hypothesis can be updated each time new evidence is added.

Bayes Theorem calculates the probability of an event A given another event B has occurred:

P (A|B) = P (B|A).P (A) P (B)
where P(A) and P(B) are respectively probabilities of event A and event B, P(A|B) presents the probability of A given B, P(B|A) is the probability of B given A.

B. K-nearest neighbors

K-nearest neighbors (KNN) is an algorithm that can be used for regression and classification, but which is more frequently used for classification predictive problems. KNN is based on the principle that instances with similar properties exist generally in close proximity. The class value of the unclassified instance can be predicted by observing the class of its nearest neighbors. The KNN finds the only adjustable parameter in this algorithm, the K nearest neighbors to the query instance and classify it by determining the most frequent class label. The model can be made less or more flexible by modifying K.

C. Logistic Regression

Logistic regression known as logit model estimates the relationship between one dependent variable and independent variables. It predicts the probability of an event utilizing a logit function.

The logistic regression is a linear regression of the explanatory variables x i :

a 0 + a 1 * x 1 + a 2 * x 2 + a 3 * x 3 + ... + a k * x k ≈ y,
where the variable y could be (1) binary: only two possible outcomes such as good or bad. (2) Multinomial: three or more non-ordered categories such as single, married or divorced. Or (3) ordinal: three or more ordered categories such as food quality from 1 to 10.

D. MultiLayer Perceptron

Multilayer perceptron is a neural network [START_REF] Leya | Long short-term memory for predicting firemen interventions[END_REF] that consists at least of one hidden layer of nodes other than the input and output layers. This means, the layers of MLP should be a minimum of 3 layers including hidden, input and output. Each node of the layer, excluded the input, is called neuron. The neuron, a processing component, which is the main element of the neural network, collects information from a certain number of inputs, applies a weight, add a bias term and send the result to an active function which generates an output.

A multilayer perceptron model consists mainly of a linear activation function of all the neurons and backpropagation process for the training. The activation function maps the weighted inputs to the output of the neuron: it combines the input of the neuron with the weights, then adds a bias in order to generate the output. On the other hand, the backpropagation is supervised learning process that occurs with a continuous adjustment after each processing (based on the error in output) of the weights of the connections.

The backpropagation consists of two parts:

• Forward pass where the expected outputs that corresponds to given inputs are evaluated; • Backward pass where partial derivative of the cost function is propagated through the network.

E. Support Vector Machine

SVM is a supervised machine learning algorithm, which is based on the concept of classification hyperplane, created to isolate between positive and other cases allowing for maximization. It plots each instance as a point in a n-dimensional space where n is the number of features (if the number of features is only 2, the hyperplane will be a simple line).

Moreover, the classification is done by finding the hyper-plane that best segregate the two classes with a maximum distance between points of both classes. Support Vectors are the closer points for hyperplane used to maximize the margin between hyperplane and data points. Also, they affect the orientation and position of the hyperplane.

F. Classification and Regression trees

Decision trees are used to explain a value from a series of discrete or continuous variables. These are fairly efficient, non-parametric and non-linear methods of partitioning individuals, producing groups of individuals as homogeneous as possible from the point of view of the variable to be predicted, taking into account a hierarchy of the predictive capacity of the variables considered. This hierarchy makes it possible to visualize the results in a tree, and to constitute explicit explanatory rules.

Several iterations are necessary. To each of them: 1) individuals are divided into k (=2) classes, to explain the output variable; 2) the first division is obtained by choosing the explanatory variable that will best separate individuals; 3) this division defines sub-populations, represented by nodes of the tree; 4) each node is associated with a proportion measure, which makes it possible to explain the belonging to a class or the meaning of an output variable; 5) the operation is repeated for each sub-population until no further separation is possible. A decision tree can very quickly lead to overfitting, so it is necessary to prune the tree: stop at an adequate number of leaves when making the tree. In order to implement all this, three main algorithms exist: CART, C4.5, and CHAID [START_REF] Delen | Measuring firm performance using financial ratios: A decision tree approach[END_REF], [START_REF] Roger | An introduction to classification and regression tree (cart) analysis[END_REF], [START_REF] Ross | C4. 5: programs for machine learning[END_REF]. They proceed as follows:

• Choice of the decision variable:

measure of the χ 2 difference to independence for CHAID and the t of Chyprow, -Gini index (or split criterion) for CART entropy for C4.5.

• Adjusting the size of the tree:

post-pruning for CART and C4.5: the purest tree is made with all the segmentation, then a criterion is used to compare trees of different sizes pre-pruning for CHAID: a stop rule is set to stop the construction.

G. Random Forest

Decision trees have the following main flaws: performance too heavily dependent on the initial sample, and a tree topology that can change completely with the input of some additional observations.

To overcome these problems, several trees are used. And to avoid having equal trees, randomness is added: each tree has a fragmented vision of the problem, randomly drawn from the input observations, and from the explanatory variables. More precisely, the assembly of decision trees built on the basis of a random draw among the observations is the tree bagging algorithm. The random forests (RFs), proposed by Leo Breiman [START_REF] Breiman | Random forests[END_REF], [START_REF] Elghazel | Random forests for industrial device functioning diagnostics using wireless sensor networks[END_REF], add a feature sampling to the tree bagging [START_REF] Elghazel | Random forests for industrial device functioning diagnostics using wireless sensor networks[END_REF].

H. AdaBoost

While RFs build several trees in parallel, boosting also builds k trees (or other basic algorithms [START_REF] Couchot | Anonymously forecasting the number and nature of firefighting operations[END_REF]), but it does so in series. The k+1 tree will have access to its predecessor, or more precisely to the latter's error: it will concentrate its effort on correcting these errors. For a classification problem, prediction is no longer a majority vote, but a weighted sum of each of the weak algorithms.

The first implementation of boosting, proposed by Yoav Freund and Robert Shapire, is called AdaBoost (adaptive boosting [START_REF] Robert E Schapire | Explaining adaboost[END_REF]), which starts from the following idea: a meta-algorithm works successively weak algorithms, each having access to a different distribution of the problem, focusing on observations that are difficult to process, and thus forcing its successor to treat them correctly. The term boosting refers, in a broad sense, to methods operating on this principle of serial assembly of weak learners.

In the case of RF, weak learners are unit decision trees, built in a totally independent way. Each algorithm has the same importance for the final vote. Boosting is a little less democratic, and achieves a sum weighted by the final vote. The weighting coefficients α i , in AdaBoost, depend only on the errors ε i of each weak learner, as follows:

α i = 1 2 ln 1-εi i .
And at the AdaBoost level, the weak learner i + 1 receives a different distribution of the data than the i-th, the latter's errors having been overweighted.

III. RELATED WORK

Despite all the efforts made by researchers, it may never be possible to know the exact time of earthquakes because there is some random component in such a catastrophe. In this field of research, errors are obviously predominant, even if exceptions exist. T. Mastuwaza et al. [START_REF] Matsuzawa | Characteristic small-earthquake sequence off sanriku, northeastern honshu, japan[END_REF], researchers in the center for prediction of earthquakes and volcanic eruptions in Japan, compared for instance the recurrent seismograms every 5.3 ± 0.53 years and summed up their study by predicting that an earthquake will occur on November 2001 in Sanriku (Japan) with a probability larger than 0.99. And actually, as expected, on 13 November 2001 an earthquake of magnitude 4.8 happened.

Many machine learning based earthquakes research studies have been done with various motivations and goals such as determining the occurrence of earthquake, predicting the magnitude, the time, the location, detection of damage, and many other classification problems and methods. A technique done by A. Negarestani, 2002 [4] using layer neural network estimated the radon concentration in soil related to the environmental parameters. The change of soil radon is not only an earthquake precursor (anomaly phenomena in the earth) but is also controlled by the environmental parameters such as temperature, humidity, rainfall, etc. The data was obtained from a location in Thailand, processed using neural network of two hidden layers and analyzed. The result of this study indicates the ability to distinguish between time variation in radon concentration raised by earth anomaly phenomena such as earthquake and those caused by environmental parameters. This technique, in comparison with linear computational methods, can estimate better the radon variations related to environmental parameters.

Another investigation done by Rouet-Leduc et al. in 2017 [START_REF] Rouet-Leduc | Machine learning predicts laboratory earthquakes[END_REF] emphasized that predicting the magnitude and timing before an earthquake fails is a main goal for geoscientist. This was done by listening to acoustic hidden signals that precede lab quakes in a laboratory that resemble to the earth: they apply Random Forest model to a continuous time-series data in a laboratory where each tree will predict the time remaining before the next failure of lab quake based on statistical features derived from time windows. In comparison with Naive Bayes model that gives a performance of 0.3, the RF achieves an accuracy of 0.89 which is largely better.

A. Cooner et al. [START_REF] Cooner | Detection of urban damage using remote sensing and machine learning algorithms: Revisiting the 2010 haiti earthquake[END_REF] made a study after the earthquake happened on January 2010 near Port-au-Prince, Haiti, 7.0 moment magnitude with 316,000 deads [START_REF] Negarestani | Layered neural networks based analysis of radon concentration and environmental parameters in earthquake prediction[END_REF], in order to detect the damage caused by this devastating event by applying Random Forests (400 trees), Neural Networks (two hidden layers of 20 neurons each) and Radial Basis Function Neural Network (150 Gaussian functions) algorithms. The remote sensing data was sourced from DigitalGlobe Foundation and the goal was to compare between a pre-disaster image captured in 2009 and a post-disaster image after the earthquake. Damaged and undamaged pixels were used for training datasets, then preprocessed. The accuracy found was best in RBFNN with 77.26%. ANN and RF, for their part, had accuracy of 74.14% and 76.14% respectively.

Another study achieved by K. M. Asim et al. [START_REF] Km Asim | Earthquake magnitude prediction in hindukush region using machine learning techniques[END_REF] predicted the magnitude of earthquake in Hindukush region by using historic seismic activated with machine learning classifiers. They applied four machine learning algorithms on a dataset extracted from temporal distribution of past earthquake: Pattern recognition neural network, Recurrent Neural Network, Random Forest (50 trees) and linear programming boost ensemble. The accuracy to predict the earthquake occurrence in this study was best using Pattern recognition neural network or Linear Programming Boost Ensemble with 65% of accuracy. However, for Recurrent Neural Network and Random Forest, the scores were of 58% and 62% respectively.

J. Reyes et al. [START_REF] Reyes | Neural networks to predict earthquakes in chile[END_REF] presented in their study some neural networks to predict earthquakes in Chile in four different areas. The database was obtained from the Chiles National Seismological Service. They applied one neural network with 7 input neurons and 15 neurons in hidden layer, and then compared the prediction results with K-Nearest Neighbors (KNN), Support Vector Machine (SVM), and k-means. They finally found that none of the algorithms obtained a better result than their approach in more than one area.

A. Li and L. Kang [START_REF] Li | Knn-based modeling and its application in aftershock prediction[END_REF] proposed a method called PR-KNN that combined between Polynomial Regression and K Nearest Neighbors models to predict the aftershock with magnitude greater or equal to 4.0. They selected the values in training sample based on the KNN algorithm, then these attributes have been modeled by applying polynomial regression method. The experimental data was collected from Wenchuan website. Finally, PR-KNN was compared with Distance-Weighted KNN regression and traditional KNN regression algorithms, which shown that the maximum relative error was reduced by 7.751% and 6.012% respectively.

IV. EXPERIMENTS

Further introduction for the world of machine learning, and explanations of the mathematical concepts for algorithms used in this paper are in [START_REF] Bonaccorso | Machine Learning Algorithms: A Reference Guide to Popular Algorithms for Data Science and Machine Learning[END_REF].

A. Experimental protocol

In this study, a single time series set of data taken from an Earthquake data center in Nothern California has been considered, where the first reading was on 1967 and the last in 2003. Each data point is an average reading for 1 hour. Let us notice that there is no use of future or past information while making the prediction. Each prediction employs only the data onto one single time series of this earthquake center. Thus, by analyzing the huge number of data point, the events were classified between negative and positive major earthquake. It is known that major events or major shocked are followed by small earthquakes in the same geographical area called aftershocks, and only reading above 5 on the Rictor scale are considered major event. The objective of this dataset is not to detect the aftershocks, but to classify between negative or positive earthquake major event. Negative output are instances having readings less than 4 and preceded, in the past 512 hours, by at least 20 non-zero readings.

On the other hand, a positive case is defined by a major event which was preceded, for at least 512 hours, by another major event.

The experiment is based on machine learning algorithms or in other terms, statistical learning theory that involves a training test (size 322) with associated inputs and outputs. Then, the trained model is utilized to evaluate a testing set (size 139) which is totally independent of the training set.

As stated previously, we used the following algorithms, for the sake of comparison:

-Random Forest (RF); -Support Vector Machine (SVM); -Logistic Regression (LR); -Naive Bayes (NB); -KNN; -Multilayer Perceptron; -AdaBoost; -CART. All the methods were applied on the training and testing sets using Matlab release 2018b and Weka, on a 2.7 GHz core i7 processor with 8 GB RAM. In obtained results of classification, several evaluation criteria have been considered:

• TP or true positive: number of times the algorithm classifies the event as 1 and it is actually 1. • FP or false positive: number of times the algorithm classifies the event as 1 and it is actually 0. • TN or true negative: number of times the algorithm classifies the event as 0 and it is actually 0. • FN or false negative: number of times the algorithm classifies the event as 0 and it is actually 1. • Mean absolute error:

M AE = 1 n n j=1 |y j -x j |
where n is the number of all the values, y is the actual value, x is the predicted value. • Root mean squared error:

RM SE = 1 n n j=1 (y j -x j ) 2
• Prediction value: Percentage of correctly classified instances.

B. Obtained results

For all algorithms, we tried to change the batch size, but the score was not affected by this parameter.

1) Naive Bayes: By applying Naive Bayes algorithm, 66.9% of instances were classified correctly, however 33.1% were not accurate.

The prediction of the major event was the worst between all the algorithms tested with 62.58%. 52 instances were classified wrong out of 139.

2) K-nearest neighbors: We applied for the KNN algorithm different number of K (neighbors) on the earthquake dataset as shown in Table I (MAE stands for Mean Absolute Error, while RMSE is for Root Mean Squared Error). The convergence of percentage of prediction occurs after using 5 neighbors. The highest accuracy reached goes for 3 neighbors with 75.53% (Figure 1) of right instances classified. 3) Logistic Regression: The prediction of the major event was the worst between all the algorithms tested with 62.58% where 52 instances were classified wrong out of 139. However, after changing the ridge to 1, the accuracy increases to 68.34% and the number of wrongly classified instances becomes 44.

4) MultiLayer Perceptron: We use in this experiment two hidden layers with three different couple of nodes: (10,20), (20,30) and (40,50). The comparison between the prediction value, MAE and RMSE using these different couples of nodes is illustrated in Table II. After choosing (20,30) as the number of nodes for two hidden layers, we apply many tests on other parameters and realize that the best score for prediction and MAE was found by using 20 epochs (Figure 2), by setting learning rate to a value of 0.6 (Figure 3), the momentum to a value of 1 (Figure 4) and by using sigmoid activation function. After changing the values of the batch size, we found that this parameter does not have any impact on the results of the experiment. The algorithm that performs better in term of accuracy is Random Forest with 76.97% (see Figure 6), which is very close to KNN with 75.53% and both MLP and SVM with 74.82%. SVM, MLP, and AdaBoost perform between 72.66% and 74.82%, while 70.5% was the prediction for CART. On the other hand, MultiLayer Perceptron and Support Vector Machine are the algorithms having the minimum average of the errors in the prediction set with 0.25 as value of MAE (Figure 8). However, the algorithms giving the worst prediction percentage were Naive Bayes and Logistic Regression with 66.9% and 66.9% respectively. The comparison between TP, FP, TN, and FN is proposed in Figure 7.

V. CONCLUSION

Eight machines learning algorithms have been tested for our work to classify the major earthquake events between

Fig. 1 .

 1 Fig. 1. Percentage of prediction for different number of neighbors using KNN Algorithm

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Comparison between Mean Absolute Error using different number of epochs

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Percentage Prediction comparison between 8 algorithms

TABLE II

 II In our experiment, we used two different functions in SVM: PolyKernel and Normalized Poly Kernel with 64.02% and 74.82% accuracy respectively as shown in TableIII.7) Random Forest: Different number of trees are applied on the dataset. TableIVand Figure5show the statistical features for each number of trees. The best

	STATISTICAL FEATURES USING MLP
	5) Classification and Regression Trees: CART algo-
	rithm leads to 70.5% of accuracy while 28 instances were
	not classified correctly.
	6) Support Vector Machine:

TABLE IV STATISTICAL

 IV FEATURES FOR DIFFERENT NUMBER OF TRESS FORFor each algorithm, we collected the True Positive, False Positive, True Negative, and False Negative rates, together with Mean Absolute Error, Root Mean Squared Error and Percentage of accuracy in TableV.

		TP	FP	TN	FN	MAE	RMSE	Accuracy
	RF	5	2	102	30	0.312	0.4161	76.97
	NB	9	20	84	26	0.3257	0.5548	66.90
	LR	6	29	89	15	0.3851	0.3161	68.34
	MLP	0	35	104	0	0.2518	0.5018	74.82
	AdaBoost	11	14	90	24	0.2958	0.4177	72.66
	KNN	2	1	103	33	0.3002	0.4565	75.53
	SVM	0	0	104	35	0.2518	0.5018	74.82
	CART	7	13	91	28	0.3078	0.4607	70.50
				TABLE V
	STATISTICAL FEATURES COMPARISON BETWEEN 8 ALGORITHMS

RANDOM FOREST ALGORITHM Fig. 5. Percentage of prediction for different number of trees using Random Forest Algorithm C. Comparing results of all algorithms