<figure><img src="image_1.png" /></figure>

HAL
open science

# Remaining useful life prediction of PEMFC systems under dynamic operating conditions

Zhiguang Hua, Zhixue Zheng, Elodie Pahon, Marie-Cécile Péra, Fei Gao

<figure><img src="image_2.png" /></figure>

## To cite this version:

Zhiguang Hua, Zhixue Zheng, Elodie Pahon, Marie-Cécile Péra, Fei Gao. Remaining useful life prediction of PEMFC systems under dynamic operating conditions. Energy Conversion and Management, 2021, 231, pp.113825. 10.1016/j.enconman.2021.113825. hal-03252611

HAL Id: hal-03252611

https://hal.science/hal-03252611

Submitted on 19 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Remaining Useful Life Prediction of PEMFC Systems Under Dynamic Operating Conditions

Zhiguang Hua¹,², Zhixue Zheng³, Elodie Pahon¹,², Marie-Cécile Péra¹,², Fei Gao¹,²

¹FEMTO-ST Institute, Univ. Bourgogne Franche-Comté, UTBM, CNRS, Rue Ernest Thierry Mieg, F-90010 Belfort, France.

²FCLAB, Univ. Bourgogne Franche-Comté, UTBM, CNRS, Rue Ernest Thierry Mieg, F-90010 Belfort, France.

³ LMOPS lab, Université de Lorraine & CentraleSupélec, 2 Rue Edouard Belin, 57070 Metz, France.

zhiguang.hua@utbm.fr

*Abstract*—The Prognostic and Health Management (PHM) has been developed for more than two decades. It is capable to anticipate the impending failures and make decisions in advance to extend the lifespan of the target systems, such as Proton Exchange Membrane Fuel Cell (PEMFC) systems. Prognostic is a critical stage of PHM. Among various prognostic methods, the data-driven ones could predict the system lifespan based on the device's knowledge and historical data. In the Remaining Useful Life (RUL) prediction, the Health Indicators (HIs) should be able to reflect the health states of the PEMFC stack. Moreover, an effective HI could help to define an explicit degradation state and improve the prediction accuracy. The HIs of voltage and power are usually used under static conditions due to their monotonic decreasing characteristics. Besides, the measurements of voltage and current are implemented easily in practice. Nevertheless, the static HIs are unable to be directly used under the dynamic operating conditions because they are sensitive to the mission profiles. To overcome the weakness of static HIs, a convenient and practical HI named Relative Power-loss Rate (RPLR) is proposed herein. According to the polarization curve at the beginning of life, the initial power under different mission profiles can be identified. Then the actual power is obtained by monitoring the current and voltage continuously. Finally, the RPLR is calculated based on the initial power and actual power. Afterward, the RUL of PEMFC is predicted by some Artificial Intelligence (AI) prognostic algorithms. Among the various data-driven prognostic approaches, Echo State Network (ESN) has provided an efficient and promising solution for the RUL prediction of PEMFC systems. Compared with classical Recurrent Neural Network (RNN), it could accelerate the convergence rate and reduce the computational complexity. Nevertheless, the traditionally used single-input ESN structure is feeble to handle the varying mission profiles. As a scheduling variable, the current is an interesting parameter since it represents the working properties to some extent. Considering the system's dynamic characteristics, the stack current is regarded as another input of ESN, and the output matrix's dimension is increased at the same time. Therefore, a double-input ESN structure is proposed to enhance the
prediction performance. Based on the dynamic HI of RPLR, three dynamic micro-cogeneration (μ-CHP) durability tests of PEMFC systems are used to verify the improved ESN prediction structure.

Keywords—Fuel cell, Degradation, Health indicator, Remaining useful life, Dynamic operating condition, Data-driven prognostic

1. Introduction

Proton Exchange Membrane Fuel Cell (PEMFC) has many excellent properties such as rapid startup, high power density, environmentally friendly, and low working temperature [1-3]. It is considered a promising power conversion device in the areas of aerospace, vehicle transportations, portable devices, distributed generation, and so on [4-6]. Nevertheless, expensive cost and limited lifetime are two large obstacles for their extensive commercial application [7]. In practice, taking actions to enhance the controlling efficiency and prolong the lifespan of the PEMFC systems are indispensable to cope with the issue of short service life [8], [9]. The framework of Prognostic and Health Management (PHM) could identify the deteriorating trend of PEMFC systems and predict the Remaining Useful Life (RUL) in the early stage [10]. Then some optimization control methods are used in advance to extend the lifespan or remind the user to replace the PEMFC stack before its failure [11], [12]. Various prognostic methods are brought up and applied in the field of the PEMFC system's RUL prediction. Generally, these methods are divided into three categories: model-based, hybrid, and data-driven ones [13], [14]. Before implementing these prognostic methods, the Health Indicators (HIs) should be determined first to indicate the degradation state. An efficient Health Indicator (HI) could significantly improve the prediction accuracy and its selection depends on the actual working conditions.

Stack voltage and power are usually used as the HIs in the static and quasi-dynamic operating conditions as they can be applied directly [15]. Besides, the convenient installation of current and voltage sensors has provided a foundation for the measurements of these parameters. In 2014, the “IEEE PHM Data Challenge” (“Data Challenge” for short) was launched by the IEEE reliability society, FCLAB research federation, FEMTO-ST Institute, and the laboratory of excellence ACTION [16]. This challenge aimed at the RUL prediction of a PEMFC system and provided the data both in the static (FC1) and quasi-dynamic (FC2) operating conditions. In the “Data Challenge”, the RUL was defined as the time when the power losses reached a specific threshold (3.5%, 4.0%, 4.5%, 5.0%, and 5.5%). Lots of ensuing articles in lifespan prediction areas are based on this RUL definition and these two experimental databases. In [17] and [18], the stack power and stack voltage are regarded as the HIs separately and both of them are tested under the static and quasi-dynamic scenarios. In [17], the discrete wavelet transform is used to decrease the prediction length and four mathematical models are used to evaluate the RUL online. Similar to [17], the Long Short-term Memory (LSTM) method and the Auto-regressive Integrated Moving Average (ARIMA) model are combined in [18] to improve the prediction performance. Based on the Bayesian framework, the hybrid
prognostic methods of Particle Filtering (PF) [19], Adaptive Particle Filtering (APF) [20], Unscented Kalman Filter (UKF) [21], and Adaptive Unscented Kalman Filter (AUKF) [22] are used to evaluate the RUL of PEMFC systems. In these four articles, the stack voltage is regarded as the HI and all the data come from the "Data Challenge". In [19] and [20], prediction results of different models (linear, exponential, logarithmic, etc.) are compared and the RUL estimation is obtained by considering the successive probability distributions of degradation states. Based on the semi-empirical degradation models, the UKF and AUKF are used in [21] and [22] to extract the HIs, and then the prognostics of PEMFC systems are realized. Under the AUKF algorithm, the system process covariance can be adjusted automatically and the initial parameters setting problem of UKF can be solved as well. The HI of stack power which considers characterization disturbances and voltage recovery is also used in [23] and [24].

In the data-driven method, stack voltage is also the most popular static HI since it could reflect the overall behavior of the stack [25]. The data-driven method based on Artificial Intelligence (AI) attracts more and more attention in recent years because it does not depend on the system behavior models. In this method, a large amount of operation data (voltage, current, power, etc.) are required for future state prediction. Various data-driven technologies are proposed in the RUL prediction of PEMFC systems, such as the Relevance Vector Machine (RVM) [26], Adaptive Neuro-Fuzzy Inference Systems (ANFIS) [27], Summation Wavelet-Extreme Learning Machine (SW-ELM) [28], Gaussian Process State Space (GPSS) [29], Long Short-Term Memory (LSTM) [30], Echo State Network (ESN) [31], [32], etc. The RVM shows a satisfactory performance in the field of long-term ageing trend prediction, and the self-adaptive RVM is proposed in [26] to predict the degradation of the PEMFC systems. The ANFIS is used in [27] to predict the output voltage during a long-term operation condition and shows that the prediction accuracy is highly dependent on the data quality. In [28], the SW-ELM has a better performance at rapid learning and parameter initialization. Besides, dual activation functions in the neuronal node could improve the robustness of the prediction process. In [29], the GPSS framework is proposed for estimating the RUL and the degradation tendency is represented in the form of a probability distribution. To avoid gradient exploding and vanishing problems of traditional Recurrent Neural Network (RNN), LSTM is used to extend the lasted time of short-term memory [30]. At last, the aging phenomena of the PEMFC systems are explicitly predicted. To conquer the shortcomings of bifurcations, heavy computational burden, and tardy convergence speed of RNN, a novel training structure called Echo State Network (ESN) was put forward in [33], [34]. A randomly dynamic reservoir is adopted to replace the hidden layer of RNN. In ESN, the input weight and internal weight are random sparse matrices and need not be trained once their structures are decided. Thus, the computation time will be further reduced since only the output matrix is calculated by the linear regression. Recently, the influences of operating parameters on the prediction accuracy are analyzed in [35] and three types of ESN and their performance are compared in [36]. In the above-mentioned data-driven methods, the stack voltage is utilized as the HI. Based on the ESN framework, the cell voltage is also used as the HI in [31] and [32].
Except for the voltage and power, the HIs of Electro-Chemical Surface Area (ECSA) and Open Circuit Voltage (OCV) under the static conditions are also proposed in [37] and [38]. Compared with the voltage and power, measuring the ECSA needs to stop the power supply and change the gas feeding in the anode side. It is hardly achievable online. Measuring the OCV is reachable if no power is supplied by the fuel cell. At the same time, relaxation time is necessary to implement the relevant measurement, and this is a constraint towards the real application. To express the degradation trend of PEMFC precisely, the hybrid HIs are introduced in [39] and [40]. In [39], three feature parameters, i.e., inner resistance, stack voltage, and stack power are employed to calculate the HI of geodesic distance. Afterward, the RUL is evaluated by the Gaussian degradation model and the Unscented Particle Filter (UPF). In this method, the weight of these three feature parameters and the initialization of UPF should be adjusted carefully to ensure the prediction precision. To improve the adaptation capability of prognostic in PEMFC systems, two degradations HIs of stack voltage and State of Health (SoH) are lumped together in [40] to predict the lifetime. Two different physics-based models and different filters are utilized together to estimate the lifetime. However, adjusting the weight and bias correction between these two HIs is very difficult. The data of [39] and [40] come from the "Data Challenge", and the hybrid HIs are tested under static and quasi-dynamic operating conditions.

Therefore, finding a practical and convenient HI under the dynamic condition makes sense for the prognostic of the PEMFC systems. In practice, the current and voltage sensors are easy to be measured and they are always supervised for the control purpose. Based on the current and voltage, the dynamic HI of Relative Power-loss Rate (RPLR) is proposed in this paper to cope with the dynamic load conditions. The single-input ESN (SI-ESN) structure has been utilized in some articles to realize the lifespan prediction. With this structure, the historical RPLR is taken as the input of ESN and the predicted RPLR is regarded as the output. Moreover, the ESN can handle multiple input tasks and consider the parameters' influence which is closely related with HIs to improve its dynamic performance. In the μ-CHP applications, the scheduled current is an extremely meaningful operating parameter because it could reflect the system's dynamic characteristics. Besides, the degradation speed is also related to the load current. To enhance the prediction accuracy of RUL, a novel structure of double-input ESN (DI-ESN) is proposed for the PEMFC systems. As the double-input parameters, the scheduled stack current and the RPLR are combined to predict the degradation tendency of RPLR. And then the RUL of PEMFC is estimated under different dynamic operating conditions. This paper is structured as follows. Three long-term duration tests and health indicator extraction procedures are introduced in section 2. The mathematical background and the implementing principles of SI-ESN and DI-ESN are shown in section 3. The experimental results under the dynamic operating conditions and the sensitivity analysis of the parameters can be found in section 4. Finally, the conclusions are given in section 5.

## 2. Tests under dynamic operating conditions
In the framework of the French National Research Agency (ANR) project, PROPICE "Prognostic and Health Management of PEM Fuel Cell Systems" [41], three dynamic tests (382 h, 1000 h, and 405 h) under the micro-cogeneration (μ-CHP) conditions were performed. The main purpose of this project is to develop prognostic methods to estimate the RUL of PEMFC systems. Based on the calculated current power and the polarization curve at the Beginning of Life (BoL), the HI of Relative Power-loss Rate (RPLR) under different dynamic mission profiles is extracted. Then the SI-ESN and DI-ESN are implemented separately to predict the lifespan of PEMFC systems.

## 2.1 Three long-term duration tests

In this task, the experimental data come from the stationary PEMFC system which contains eight fuel cells and are obtained under dynamic mission profiles. The electrical power of the test bench is 1.0 kW and the active area of a single cell is 220 $cm^2$. Some of the nominal specifications and working parameters are shown in Table 1. For all the durability tests, the polarization curves and the Electrochemical Impedance Spectra (EIS) under different current densities are measured about once a week. They are used for analyzing the static and dynamic characteristics of PEMFC systems. The acquisition frequency of stack voltage and stack current is 1 Hz. Both the inlet and outlet reactant properties (flow rate, temperature, pressure, humidification, etc.) could be controlled during the tests and various sensors are installed to monitor the operating conditions [42]. The experimental platform is shown in Fig. 1.

<table><caption>Table 1. Nominal specifications and working parameters of the experimental platform.</caption><thead><tr><th>Parameters</th><th>Control range</th></tr></thead><tbody><tr><td>Dimensions</td><td>(220×160×186) mm</td></tr><tr><td>Temperature</td><td>80 °C</td></tr><tr><td>Current density</td><td>0.5 A/cm²</td></tr><tr><td>Anode/cathode stoichiometry</td><td>1.5/2</td></tr><tr><td>Anode/cathode inlet pressure</td><td>150/150 kPa</td></tr><tr><td>Cooling flow rate</td><td>2 L/min</td></tr><tr><td>Pressure drop</td><td>30 kPa</td></tr><tr><td>T<sub>inH2</sub>; T<sub>outH2</sub></td><td>Inlet/Outlet temperatures of H<sub>2</sub> (°C)</td></tr><tr><td>T<sub>inAir</sub>; T<sub>outAir</sub></td><td>Inlet/Outlet temperatures of air (°C)</td></tr><tr><td>T<sub>inWat</sub>; T<sub>outWat</sub></td><td>Inlet/Outlet temperatures of water (°C)</td></tr><tr><td>P<sub>inH2</sub>; P<sub>outH2</sub></td><td>Inlet/Outlet pressure of H<sub>2</sub> (mBar)</td></tr><tr><td>P<sub>inAir</sub>; P<sub>outAir</sub></td><td>Inlet/Outlet pressure of air (mBar)</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Fig. 1. The 10 kW experimental platform of the PEMFC system.</figcaption></figure>

The mission profiles of three long-term tests are shown in Fig. 2. The transient effects of stack current on the stack voltage are marked in red circles. The effects of start and stop operations on stack voltage are marked in green rectangles. Besides, the operating parameters (e.g., T<sub>in</sub>H<sub>2</sub> in this case) would also influence the stack voltage. The effects of T<sub>in</sub>H<sub>2</sub> changes on the stack voltage are marked in black triangles.

<figure><img src="image_5.png" /></figure>
Fig. 2. The mission profiles of three long-term tests: (a) stack current of Test-A, (b) stack voltage of Test-A, (c) stack current of Test-B,

(d) stack voltage of Test-B, (e) stack current of Test-C, (f) stack voltage of Test-C.

The 382 h duration data (Test-A) have three operating conditions and can be divided into 3 stages. The stack current and stack voltage are presented in Fig. 2 (a) and Fig. 2 (b).

(1) In stage 1 (0-200 h), the load current density is 0.77 A/cm² (maximum).

(2) In stage 2 (200 h-300 h), the load is cycling between 0.77 A/cm² and 0.38 A/cm².

(3) In stage 3 (300 h-382 h), the load current density is 0.38 A/cm².

The 1000 h duration data (Test-B) are divided into 4 stages and each stage lasts for 250 h. The stack current and stack voltage are presented in Fig. 2 (c) and Fig. 2 (d).

(1) In stage 1 (0-250 h), the load current density is maximum at 0.77 A/cm².

(2) In stage 2 (250 h-500 h), the mission profile changes dynamically between no load (0 A/cm²), 0.77 A/cm², and 0.38 A/cm² from 250 h to 380 h. The load current density alternates between 0.45 A/cm² and 0.23 A/cm² from 380 h to 500 h.

(3) In stage 3 (500 h-750 h), the load current density is 0.23 A/cm² from 500 h to 600 h. The load current density alternates between 0 A/cm² and 0.23 A/cm² from 600 h to 750 h.

(4) In stage 4 (750 h-1000 h), the load current density alternates between 0 A/cm² and 0.23 A/cm² from 750 h to 850 h. The load current density is 0.23 A/cm² in the rest time of stage 4.

The 405 h duration data (Test-C) are also divided into 4 stages. The stack current and stack voltage are presented in Fig. 2 (e) and Fig. 2 (f).

(1) In stage 1 (0-125 h), the load current density is 0.36 A/cm² (0-25 h) and it changes to 0.45 A/cm² during the rest of the time.

(2) In stage 2 (125 h-250 h), the load current density changes between 0.45 A/cm² and 0.23 A/cm² from 125 h to 225 h. It stays at 0.23 A/cm² during the rest of the time.

(3) In stage 3 (250 h-375 h), the load current density alternates between 0 A/cm² and 0.23 A/cm².

(4) In stage 4 (375 h-405 h), the load current density is 0.23 A/cm².

## 2.2 Dynamic health indicator extraction

In the static operating condition, the deviation of voltage and power are considered to be influenced only by the ageing degradation. Nevertheless, the voltage and power are also influenced by the mission profiles and it is improper to take them as the HIs in the dynamic or time-varying operating conditions. Finding the dynamic HI is difficult and just a few papers have presented the handling methods. One way is to transform the dynamic condition into the static representation [43], [44]. Another way is to extract some monotonic parameters during the degradation [45], [46]. The Linear Parameter Varying (LPV) equivalent models are built to fit the voltage signal segments in [43]. Based on the modeling identification method, the HI, i.e.,
virtual steady-state stack voltage is rebuilt from the model parameters. Later, an ensemble ESN structure is used in [44] to enhance the prognostic adaptability. Nevertheless, the extraction of HI is time-consuming, and the computing time is based on the number of LPV. It is worth mentioning that this method is more suitable for the conditions in which the load varies regularly. In [45], the overall resistance and limiting current density are extracted from the polarization curves, and then an Extended Kalman Filter (EKF) is adopted to evaluate the actual health state. Then, the inverse first-order reliability method is used in [46] to predict the RUL. Meanwhile, the quantification of uncertainty in the dynamic operating condition is also considered. Nevertheless, the HI comes from the periodic measurements of the polarization curves which are characterized every week in practice. Therefore, the extracted parameters are very limited to describe the complicated degradation phenomenon. Besides, the linearization hypothesis of overall resistance and limiting current density is inaccurate for the nonlinear PEMFC systems [47].

Finding a more general and convenient HI under the dynamic operating condition is necessary for practical use. In the dynamic tests, the measurements of polarization curves at the BoL and at the End of Life (EoL, which is a special case of the end-of-test) can be easily performed. Based on the polarization curves (voltage vs. current) at the BoL, the BoL power curve (power vs. current) can be obtained. Then a mathematical model is used to fit the BoL power curve. The power supplied by the fresh stack ($P_0$) under different loads can then be computed by the model. In the PEMFC systems, the $P_0$ is not a single value, but a curve: power (W) versus current (A). For each current, there is a corresponding $P_0$, measured on the fresh stack, before aging. In practice, during the lifetime of the stack (here the duration of the test), the current power ($P$) can be calculated from the stack current and stack voltage measurements. Based on the current power ($P$) and the BoL power ($P_0$) corresponding to the same current which is extracted from the models, the health indicator of Relative Power-loss Rate (RPLR, $\Delta P$) is defined as

$$ \Delta P = (P - P_0) / P_0 \quad (1) $$

So, the calculation of RPLR is conducted as follows:

Step1: measurements of BoL power ($P_0$) are carried out at different current values. A mathematical model of $P_0$ is built, which parameters are identified thanks to the experimental values.

Step2: at time step $t$, current power ($P_t$) can be calculated by the stack current ($I_t$) and stack voltage ($V_t$).

Step3: look at the table and find out the BoL power ($P_0-I_t$) under the stack current ($I_t$).

Step4: calculate the RPLR at time step $t$ based on the current power ($P_t$) and BoL power ($P_0-I_t$). With the same method, the RPLR during its lifetime can be calculated.
Comparing with the LPV models in the dynamic operating conditions [43], the computing complexity of RPLR is lower.

Comparing with the linearization technique in [45], the time interval of each two points is no longer limited by the interval of polarization curves. This provides the implementation foundation for data-driven methods. To validate the effectiveness of the RPLR, three dynamic tests with different durations are implemented under various operating conditions. For the tuning of the parameters, the Trust-Region optimization method is used in the BoL power modeling process. The parameters and the accuracy of these three models are shown in Table 2. The polarization and the BoL power curves are shown in Fig. 3, and the typical current points of each experimental test are marked in red circles.

<table><caption>Table 2. Parameters of three mathematical models.</caption><thead><tr><th>Test</th><th>Parameters</th><th>Accuracy</th></tr></thead><tbody><tr><td rowspan="2">A</td><td>a<sub>1</sub></td><td>b<sub>1</sub></td><td>c<sub>1</sub></td><td>R-square</td><td>RMSE</td></tr><tr><td>861.9</td><td>0.007438</td><td>0.003801</td><td>0.9999</td><td>2.588</td></tr><tr><td rowspan="2">B</td><td>a<sub>2</sub></td><td>b<sub>2</sub></td><td>c<sub>2</sub></td><td>R-square</td><td>RMSE</td></tr><tr><td>878.0</td><td>0.007241</td><td>0.004635</td><td>0.9999</td><td>2.717</td></tr><tr><td rowspan="2">C</td><td>a<sub>3</sub></td><td>b<sub>3</sub></td><td>c<sub>3</sub></td><td>R-square</td><td>RMSE</td></tr><tr><td>765.7</td><td>0.008437</td><td>0.004738</td><td>0.9999</td><td>2.206</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(f)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(g)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(h)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(i)</figcaption></figure>
<figure><img src="image_15.png" /></figure>

In particular, the dynamic behavior of Test-B is the harshest, and the degradation is huge after 300 h testing at maximal power. It is difficult to reach the maximal current value on the polarization curves, and then the limiting current value moves from 170 A to 100 A after 300 h testing due to the encountered degradation. The current and voltage are resampled at a regular interval of 30 mins. The assumption herein is that the degradation state in each time interval keeps the same. This is reasonable because the resampling time interval is short enough for the hundreds of hours' degradation. The RPLR and RUL definition of these three tests are shown in Fig. 4. Results show that the RPLR in each test has the monotone decreasing property and the current profile affects the degradation speed. A higher current level always leads to a higher degradation speed, e.g., in Test-B, $\Delta P_1$ during the first 300 h is bigger than $\Delta P_2$ during the rest 700 h. An efficient HI is useful for the prognostic methods to predict the SoH and proposing an index to indicate the End of Life (EoL) time is also an important work for PEMFC systems. In practice, the U.S. Department of Energy (DOE) has pointed out that the definition of a proper index depends on the actual needs of users. It is worth mentioning that the meaningful definition of EoL is the time when the PEMFC cannot transfer the required power under safe condition rather than the stack is out of use once it operates at a degradation state. According to the definition of the U.S. DOE, the Failure Threshold (FT) can be regarded as the PEMFC reaches a certain power decay (e.g., 10 % for the vehicle application and 20 % for $\mu$-CHP) which is the external manifestation of the degradation in component level [48]. During the testing, the power loss is tested under about 55%-65% rated stack current, and the transient operating, startup, and shutdown should also be included. Besides, the U.S. DOE metric does not represent the real EoL of original equipment manufacturers, and a higher or lower level of power degradation can be accepted in different applications [49]. Due to the limitations of the experimental conditions, the duration of these three dynamic tests is all within 1000 h. The failure points in these three dynamic tests do not mean that they have reached their real EoL. Nevertheless, EoL could be regarded as one special case of end-of-test (power decay reaches certain defined levels). Without loss of generality, the end-of-test in this paper is regarded as the time of EoL.
<figure><img src="image_16.png" /><figcaption>Fig. 4. The RPLR and RUL definition of three long-term tests: (a) Test-A, (b) Test-B, (c) Test-C.</figcaption></figure>

The first 50 h of data are abandoned due to the instability at startup. So, the training part is the data from 50 h to the time when the prediction starts ($t_{predict}$), and the prediction part is the data from $t_{predict}$ to the failure threshold ($t_{failure}$). The RUL ($t_{RUL}$) is usually represented by the time interval between $t_{predict}$ and $t_{failure}$. In general, the PEMFC stack always has a satisfactory performance at the BoL and the $t_{predict}$ usually starts when the system has been worked for certain hours. Besides, the prediction cannot be accurate at BoL because the historical data are inadequate for learning the degradation characteristics, especially for the data-driven methods. Based on the extracted HI of RPLR, the purpose of this paper is to compare the prediction performance of DI-ESN and SI-ESN. So, in each test, the precision improvement of DI-ESN can be verified when these two structures (SI-ESN and DI-ESN) have the same failure points. For comparison's purpose, the end-of-test in these three tests is regarded as the failure points. The failure time for Test-A, Test-B, and Test-C are 382 h, 1000 h, and 405 h respectively.
1

## 3. Echo State Network

2

### 3.1 Mathematical background

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

The Artificial Neural Network (ANN) has been applied in different areas and it is used to construct and simulate the working principle of the brain. There are two typical structures of ANN: Feed-Forward Neural Network (FFNN) whose information only flows forward and Recurrent Neural Network (RNN) which contains the feedback loop. By imitating the behavior of biological neurons, the ANN can solve the problem of complex engineering. There does not exist recurrent connections in FFNN, therefore, it lacks memory and cannot deal with temporal information. The RNN can deal with large dynamical memory and has high computational capabilities. It means that the current states of the neurons are not only affected by the inputs but also by the historical states of the network. The working principles of RNN are more likely to biological brains, the substrate of natural intelligence. Both the FFNN and RNN has been used in many fields. Nevertheless, the training complexity is a common problem. The error Back Propagation (BP) is one of the efficient training rules to decrease the complexity in ANN especially for the FFNN. There are still some weaknesses when applying the BP to RNN such as bifurcations, slow convergence, expensive computing, and local optimum. These shortcomings still hinder the large-scale deployment of RNN in practical applications. Recently, an improved training structure of the Echo State Network (ESN) is put forward to overcome these disadvantages. This proposal is based on some evidence that the RNN has a satisfactory performance even without adaptation to all the weights in the network. Compared to the RNN, the hidden layer is replaced by a large randomly generated reservoir which reflects the dynamic topologies of the neurons. The weights (input, internal, and feedback) of ESN are global scaled to reach a desired dynamic state. And all the weight matrices are randomly generated. Then the readout of the reservoir is calculated via linear regression. Echo State Property (ESP) is an important indicator to represent the dynamic behavior of the reservoir. The ESP should be considered carefully in designing the ESN, and the reservoir should wash out the neurons' initial states at a rate that is independent of the inputs [50], [51]. On one hand, the neurons in the reservoir should be dynamic enough to decrease the computation complexity of output weights. On the other hand, too dynamic neurons may move the network to an unstable boundary. An optimized reservoir means that the neurons have rich enough dynamics, and the fading memory should also be preserved. The basic representation of an ESN is shown in Fig. 5.
<figure><img src="image_17.png" /><figcaption>Fig. 5. The basic representation of Echo State Network.</figcaption></figure>

If the dimension of the input signal is $K$, the neuron number in the reservoir is $N$, and the dimension of the output signal is $L$, the external input signal $\mathbf{u}$ is a $K$-dimensional vector, the reservoir neuron activation signal $\mathbf{x}$ is an $N$-dimensional vector, and the output signal $\mathbf{y}$ is an $L$-dimensional vector. At time step $n$, the form of them are

$$ \left\{ \begin{aligned} \mathbf{u}(n) &= (\mathbf{u}_1(n), \dots, \mathbf{u}_K(n))^T, \quad n = 1, \dots, T \\ \mathbf{x}(n) &= (\mathbf{x}_1(n), \dots, \mathbf{x}_N(n))^T, \quad n = 1, \dots, T \\ \mathbf{y}(n) &= (\mathbf{y}_1(n), \dots, \mathbf{y}_L(n))^T, \quad n = 1, \dots, T \end{aligned} \right. \quad (2) $$

Where $T$ is the total data points in the training part. With the external input signal $\mathbf{u}$, the time series prediction task aims to find a model to minimize the error between the target signal $\mathbf{y}^{\text{target}}$ and the output signal $\mathbf{y}$. The error $E(\mathbf{y}, \mathbf{y}^{\text{target}})$ which is represented by the Root Mean Square Error (RMSE) is given as

$$ E(\mathbf{y}, \mathbf{y}^{\text{target}}) = \frac{1}{L} \sum_{i=1}^{L} \sqrt{\frac{1}{T} \sum_{n=1}^{T} (y_i(n) - y_i^{\text{target}}(n))^2} \quad (3) $$

Where $i \in [1, L]$ is the output dimension. The $\mathbf{W}_{in} \in \mathbb{R}^{N \times K}$, $\mathbf{W} \in \mathbb{R}^{N \times N}$, and $\mathbf{W}_{fb} \in \mathbb{R}^{N \times L}$ represent the input weight matrix, internal weight matrix, and feedback weight matrix respectively. The neurons' states in the reservoir can be represented in the time-continuous domain as

$$ \dot{\mathbf{x}} = \frac{1}{c} (-\alpha \mathbf{x} + f(\mathbf{W}_{in} \mathbf{u} + \mathbf{W} \mathbf{x} + \mathbf{W}_{fb} \mathbf{y})) \quad (4) $$

Where $c$ is the global scaling factor for the temporal dynamics, $\alpha$ is the leaking rate of reservoir neurons, and $f(\cdot)$ is a nonlinear transformation, and the sigmoid functions (tanh and Fermi) are commonly used. The neurons' behaviors are related
1 to the shape of the activation functions. During the selection of functions, a trade-off should be made between the dynamic
2 richness and the computing complexity based on the premise of differentiable and boundedness.

3 Based on the Euler method, the typical model of ESN in discretized time is given by

$$
\mathbf{x}((t+1)\delta) = \left(1 - \frac{\alpha\delta}{c}\right)\mathbf{x}(t\delta) + \frac{\delta}{c} f\left(\mathbf{W}_{in}\mathbf{u}(t\delta) + \mathbf{W}\mathbf{x}(t\delta) + \mathbf{W}_{fb}\mathbf{y}(t\delta)\right) \quad (6)
$$

4 Where $\delta$ is the discrete step size. Using the notation $\mathbf{x}((t+1)\delta) = \mathbf{x}(n+1)$, the model of ESN can be transferred into the
5 form

$$
\mathbf{x}(n+1) = \left(1 - \frac{\alpha\delta}{c}\right)\mathbf{x}(n) + \frac{\delta}{c} f\left(\mathbf{W}_{in}\mathbf{u}(n) + \mathbf{W}\mathbf{x}(n) + \mathbf{W}_{fb}\mathbf{y}(n)\right) \quad (7)
$$

6 The output vector $\mathbf{y}$ is given by

$$
\mathbf{y}((t+1)\delta) = g(\mathbf{W}_{\text{out}}[\mathbf{x}((t+1)\delta); \mathbf{u}(t\delta)]) \quad (8)
$$

7 Where $g(\cdot)$ is a nonlinear transformation, $\mathbf{W}_{out} \in \mathbb{R}^{L \times (N+K)}$ is the output weight matrix. In the Euler discretization, the
8 step size $\delta$ should be small enough to ensure the accuracy and stability of the approximation. The Echo State Property (ESP)
9 is the essential feature of ESN which represents the reservoir working at a stable state, and the ESP can be designed as follows

$$
\left\{
\begin{array}{ll}
|\lambda|_{\max}(\hat{\mathbf{W}}) < 1 & \text{; for zero inputs} \\
|\lambda|_{\max}(\hat{\mathbf{W}}) \ge 1 & \text{; for non-zero inputs}
\end{array}
\right.
\qquad (9)
$$

10 Where $|\lambda|_{\max}(\hat{\mathbf{W}})$ is the equivalent spectral radius of $\hat{\mathbf{W}}$ in the leaky integrator neuron, and the $\hat{\mathbf{W}}$ can be represented as

$$
\hat{\mathbf{W}} = \frac{\delta}{c} \mathbf{W} + \left(1 - \frac{\alpha \delta}{c}\right) \mathbf{I} \qquad (10)
$$

11 Where **I** is the identity matrix. Besides, the neuron in the reservoir in a single update step should less than its previous
12 excitation. So, all the neurons in the reservoir should under the rule of

$$
1 - \frac{\alpha \delta}{c} = 1 - \alpha \xi \le 0 \qquad (11)
$$

13 Where $\xi = \delta / c$ is the scale gain, then the model of ESN can be rewritten as
$$ \mathbf{x}(n+1) = (1 - \alpha \xi) \mathbf{x}(n) + \xi f(\mathbf{W}_{in} \mathbf{u}(n) + \mathbf{W} \mathbf{x}(n) + \mathbf{W}_{fb} \mathbf{y}(n)) \quad (12) $$

Assuming $\Gamma$ represents an ESN with the weights of $\mathbf{W}_{in}$, $\mathbf{W}$, and $\mathbf{W}_{fb}$, the leaking rate of $\alpha$, and the scale gain of $\xi$,

then the updated model of ESN can be expressed as

$$ \begin{aligned} \frac{1}{\xi} \mathbf{x}(n+1) &= \frac{1}{\xi}(1 - \alpha \xi) \mathbf{x}(n) + f(\mathbf{W}_{in} \mathbf{u}(n) + (\xi \mathbf{W}) \frac{1}{\xi} \mathbf{x}(n) + \mathbf{W}_{fb} \mathbf{y}(n)) \\ &= \frac{1}{\xi}(1 - \tilde{\alpha}) \mathbf{x}(n) + f(\mathbf{W}_{in} \mathbf{u}(n) + \tilde{\mathbf{W}} \frac{1}{\xi} \mathbf{x}(n) + \mathbf{W}_{fb} \mathbf{y}(n)) \end{aligned} \quad (13) $$

In the updated ESN of $\tilde{\Gamma}$, $\tilde{\alpha} = \xi\alpha$, $\tilde{\mathbf{W}} = \xi\mathbf{W}$. So $\tilde{\Gamma}$ has the same updated structure with $\Gamma$ except for the scaling factor

of $1/\xi$. The property of the updated internal weight matrix $\tilde{\mathbf{W}}$ can be scaled by its spectral radius, and the updated output

vector $\tilde{\mathbf{y}}$ can be presented as

$$ \tilde{\mathbf{y}}(n+1) = g(\tilde{\mathbf{W}}_{out}[\frac{1}{\xi}\mathbf{x}(n+1); \mathbf{u}(n)]) \quad (14) $$

It is identical to the form of Eq. (8). When $\xi = \tilde{\xi} = 1$, the updated and previous models are the alternative structures.

Without loss of generality, the dynamic state in the reservoir can be expressed as

$$ \mathbf{x}(n+1) = (1 - \alpha)\mathbf{x}(n) + f(\mathbf{W}_{in}\mathbf{u}(n) + \mathbf{W}\mathbf{x}(n) + \mathbf{W}_{fb}\mathbf{y}(n)) \quad (15) $$

Then the $\mathbf{W}_{out}$ can be calculated as

$$ \mathbf{W}_{out} = \arg \min \frac{1}{L} \sum_{i}^{L} \sqrt{\frac{1}{T} \sum_{n=1}^{T} (y_{i}(n) - y_{i}^{\text{target}}(n))^{2}} \quad (16) $$

Linear regression is a common technology to learn the value of $\mathbf{W}_{out}$.

$$ \mathbf{W}_{out} = \mathbf{y}^{\text{target}} \mathbf{x}^{\text{T}} (\mathbf{x} \mathbf{x}^{\text{T}} + \beta \mathbf{I})^{-1} \quad (17) $$

Where $\beta$ is the regularization parameter. The working principles of ESN can be divided into 4 steps and the implementation procedure is shown in Table 3.

Step 1: Generate a reservoir and build up the ESN. The number of reservoir neurons $N$, spectral radius $\rho$ of the internal matrix $\mathbf{W}$, and leaking rate $\alpha$ are determined. The input matrix $\mathbf{W}_{in}$ and internal matrix $\mathbf{W}$ are respectively assigned.
Step 2: Run the ESN and collect the reservoir activation states. The appropriate activation functions of $f(\cdot)$ and $g(\cdot)$ are selected, and the states of the neurons are initialized. The external input signal **u** is sent to the reservoir, and the neuron activation signal **x** and its updated state are calculated by Eq. (7).

Step 3: Calculate the output matrix. The training part aims to compute the output matrix $\mathbf{W}_{\text{out}}$ by Eq. (16) and Eq. (17). In the training dataset, the objective of the learning algorithm is minimizing the error between the target signal $\mathbf{y}^{\text{target}}$ the output signal of the ESN $\mathbf{y}$.

Step 4: The new output vector $\hat{\mathbf{y}}$ is predicted. In the prediction dataset, the output weight matrix $\mathbf{W}_{\text{out}}$ and the new input vector $\hat{\mathbf{u}}$ are used to calculate the new outputs $\hat{\mathbf{y}}$ by Eq. (8).

<table><caption>Table 3. The implementation procedure of ESN.</caption><thead><tr><th colspan="2">Algorithm: working principle of ESN</th></tr></thead><tbody><tr><td colspan="2"><b>Input:</b> <b>u</b>, <b>y</b><sup>target</sup>, <i>T</i>, <i>K</i>, <i>N</i>, <i>L</i>, <i>α</i>, <i>ρ</i>, <i>β</i>, <b>W</b><sub>in</sub>, <b>W</b>, <b>W</b><sub>fb</sub>, <i>f</i>(·), <i>g</i>(·), <b>u</b></td></tr><tr><td colspan="2"><b>Output:</b> <b>W</b><sub>out</sub>, <b>ŷ</b></td></tr><tr><td colspan="2"><b>Step1: Network initialization</b></td></tr><tr><td colspan="2">W<sub>in</sub> ∈ (-0.5, 0.5), W ∈ (-0.5, 0.5), W<sub>fb</sub> = 0,<br>α ∈ (0, 1), ρ ∈ (0, 1.5), N =400, x(0) = 0</td></tr><tr><td colspan="2"><b>Step 2: Training based on historical input u</b></td></tr><tr><td colspan="2">for n=1 to <i>T</i>, do</td></tr><tr><td colspan="2">update <b>x</b> by Eq. (7)</td></tr><tr><td colspan="2">collect the state of <b>x</b> by Eq. (7)</td></tr><tr><td colspan="2">collect <b>y</b> by Eq. (8)</td></tr><tr><td colspan="2">end</td></tr><tr><td colspan="2"><b>Step 3: Calculate W<sub>out</sub> by Eq. (16) and Eq. (17)</b></td></tr><tr><td colspan="2"><b>Step 4: Prediction based on new input ū</b></td></tr><tr><td colspan="2">for n=<i>T</i>+1 to <i>T</i>+<i>m</i>, do</td></tr><tr><td colspan="2">calculate the new output <b>ŷ</b> by Eq. (8)</td></tr><tr><td colspan="2">end</td></tr><tr><td colspan="2"><b>Performance evaluation</b> (<b>ŷ</b> vs. <b>ŷ</b><sup>target</sup>)</td></tr></tbody></table>
1

## 3.2 Implementation of the ESN

<figure><img src="image_18.png" /><figcaption>Fig. 6. The iterative process of SI-ESN and DI-ESN.</figcaption></figure>

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

Once the RPLR is extracted, the SI-ESN and DI-ESN both with an iterative prediction structure are used to evaluate the RUL of the PEMFC systems [52]. It should be noted that the stack current is usually taken as a planned variable. This double-input long-term prediction pattern is more applicable for applications where the stack current is schedulable or programmable, such as the μ-CHP application which load is related to the seasonal variations. For SI-ESN, the training part of RPLR can be expressed as $\{(t_1, y_1), \dots, (t_n, y_n), \dots, (t_T, y_T)\}$, where $t_n$ is the sampling time, $y_n$ is the RPLR at time $t_n$. The output matrix $W_{out}$ of the reservoir is computed by the training part and the value of it stays the same after the training. With the one-step ahead prediction structure, the multi-step ahead prediction is realized and the RPLR in the prediction part can be expressed as $\{(t_{T+1}, \hat{y}_{T+1}), \dots, (t_H, \hat{y}_H)\}$, where $H$ is the number of predicted points. For DI-ESN, the historical RPLR $\{(t_1, y_1), \dots, (t_n, y_n), \dots, (t_T, y_T)\}$ and corresponding stack current $\{(t_1, z_1), \dots, (t_n, z_n), \dots, (t_T, z_T)\}$ are regarded as the inputs, and the predicted RPLR $\{(t_{T+1}, \hat{y}_{T+1}), \dots, (t_H, \hat{y}_H)\}$ is regarded as the output, where $z_n$ represents the stack current at the time $t_n$. Also based on the one-step ahead prediction structure, the predicted RPLR, and the new scheduled stack current $\{(t_{T+1}, z_{T+1}), \dots, (t_H, z_H)\}$ are used to be the inputs of DI-ESN in the prediction part. The iterative process of SI-ESN and DI-ESN is shown in Fig. 6. The influences of stack current are marked in blue.

17

18

Where $\{\Gamma, [\mathcal{G}]\}$ represents the prediction of one-step ahead ESN and $[\mathcal{G}]$ represents the parameters in the ESN. And $p$ is the number of regressors, i.e., the quantity of previous data used for the prediction. In practice, the Prediction Horizon (PH) can
be used to evaluate the prediction ability, and it is defined as the time when the prediction is located in an allowable error bound (α = 5%) around the true RUL. Besides the PH, the other three criteria are also adopted to estimate the prediction performance: Root Mean Square Error (RMSE), Mean Average Percentage Error (MAPE), and the percentage error (%ErFT).

$$RMSE = \sqrt{\frac{1}{H} \sum_{i=1}^{H} (\mathbf{y}_i(n) - \mathbf{y}_i^{\text{target}}(n))^2} \quad (18)$$

$$MAPE = \frac{1}{H} \sum_{i=1}^{H} \left| \frac{\mathbf{y}_i(n) - \mathbf{y}_i^{\text{target}}(n)}{\mathbf{y}_i(n)} \right| \quad (19)$$

$$\%Er_{FT} = \frac{t_{RUL}^{act} - t_{RUL}^{pre}}{t_{RUL}^{act}} \times 100 \quad (20)$$

Where $\mathbf{y}_i(n)$ is the predicted signal, $\mathbf{y}_i^{\text{target}}(n)$ is the target signal, $H$ is the data points in the prediction part, $t_{RUL}^{act}$ is the actual RUL, $t_{RUL}^{pre}$ is the prediction RUL. In the practical application, the training part is always calculated offline, and the prediction part is implemented online. All the simulations are implemented in the Matlab 2018a version, and the calculating time in the prediction part is less than 20 s. This execution time is rather short for the hundreds of hours of lifetime prediction. The target statement of RUL and its prediction process is presented in Fig. 7.

<figure><img src="image_19.png" /><figcaption>Fig. 7. Target statement of RUL and its prediction process.</figcaption></figure>

# 4. Experimental Results

## 4.1 Results of dynamic tests

In dynamic tests, the results of SI-ESN are first presented. Afterward, the DI-ESN is implemented and compared with the SI-ESN. In terms of the Prediction Horizon (PH), a longer PH means that the prediction accuracy meets the requirement at an early stage [53]. Besides, criteria of RMSE, MAPE, and %ErFT are combined comprehensively to evaluate the prediction performances. In general, a small value of them means better performance. It is worth noting that there are two cases in the
RUL estimation: 1) underestimation (or early prediction) which means the predicted RUL is smaller than the actual RUL, 2) overestimation (or late prediction) which means that the predicted RUL is larger than actual RUL. In practice, underestimation is more meaningful and jeopardize less the operation of the system because it could warn the user to adjust the operating condition or replace the PEMFC in advance and then to prevent an incident. To analyze the influence of stack current ($I_s$) on the prediction effect, results of SI-ESN and DI-ESN with different training lengths are compared. All the predictions are implemented at about half of the whole data because premature prediction may have the problem of insufficient training data.

The parameters ($\alpha$, $\rho$, and $\beta$) of ESN in these three tests are determined based on the sensitivity analysis of them in section 4.2.

In Test-A, the RUL predictions with 60 %, 70 %, and 80 % of the whole data (382 h) for training are presented in Fig. 8. They show that the DI-ESN (with stack current $I_s$) has better performance than SI-ESN (without stack current $I_s$) especially during the load cycling period (200 h-300 h). On one hand, the DI-ESN could track the variation tendency of RPLR ($\Delta P$) and on the other hand, the dynamic effectiveness of stack current is also superimposed to the RPLR prediction. These superposition effects are more pronounced in Fig. 8 (b) and Fig. 8 (d). In Table 4, the prediction results of Test-A with different training lengths are given. With an increase of training length, the prediction error (%$Er_{FT}$, RMSE, and MAPE) has an overall downward trend both for the SI-ESN and the DI-ESN. The prediction results of DI-ESN are more accurate than those of SI-ESN as a whole. Nevertheless, two anomalies exist at 40 % and 50 % conditions from the view of RMSE and MAPE. Most probably, the prediction cannot be accurate if the training data are not enough. The two disturbances caused by the increase of $T_{in}H_2$ at 40 h and 100 h may also affect the prediction results. When considering the %$Er_{FT}$, the prediction RUL of DI-ESN is closer to the actual RUL both at 40% and 50% conditions. The RUL predictions of SI-ESN and DI-ESN in the whole lifespan of Test-A are shown in Fig. 9. The PH of SI-ESN is 117 h and the PH of DI-ESN is 157 h.
<figure><img src="image_20.png" /><figcaption>Fig. 8. The RUL prediction of Test-A: (a) SI-ESN with 60% training, (b) DI-ESN with 60% training, (c) SI-ESN with 70% training, (d) DI-ESN with 70% training, (e) SI-ESN with 80% training, (f) DI-ESN with 80% training.</figcaption></figure>

<table><caption>Table 4. Prediction results comparison of Test-A</caption><thead><tr><th>ESN type</th><th>Training length (%)</th><th>Actual RUL (h)</th><th>Prediction RUL (h)</th><th>%ErFT (%)</th><th>RMSE</th><th>MAPE</th></tr></thead><tbody><tr><td rowspan="6">SI-ESN</td><td>40</td><td>229</td><td>707</td><td>-208.6</td><td>0.01869</td><td>0.57772</td></tr><tr><td>50</td><td>191</td><td>329</td><td>-72.3</td><td>0.01067</td><td>0.24202</td></tr><tr><td>60</td><td>153</td><td>74</td><td>51.7</td><td>0.01867</td><td>0.26439</td></tr><tr><td>70</td><td>115</td><td>133</td><td>-15.7</td><td>0.00576</td><td>0.10435</td></tr><tr><td>80</td><td>76</td><td>67</td><td>11.8</td><td>0.00429</td><td>0.06791</td></tr><tr><td>90</td><td>38</td><td>41</td><td>-7.9</td><td>0.00360</td><td>0.05453</td></tr><tr><td rowspan="4">DI-ESN</td><td>40</td><td>229</td><td>47</td><td>79.4</td><td>0.02755</td><td>1.67360</td></tr><tr><td>50</td><td>191</td><td>349</td><td>-82.7</td><td>0.04238</td><td>1.88170</td></tr><tr><td>60</td><td>153</td><td>152</td><td>0.7</td><td>0.00650</td><td>0.11326</td></tr><tr><td>70</td><td>115</td><td>123</td><td>-7.0</td><td>0.00411</td><td>0.06186</td></tr><tr><td></td><td>80</td><td>76</td><td>94</td><td>-17.8</td><td>0.00373</td><td>0.06184</td></tr></tbody></table>
<figure><img src="image_21.png" /><figcaption>Fig. 9. The RUL results of SI-ESN and DI-ESN with 95% probability bounds in Test-A</figcaption></figure>

In Test-B, the RUL predictions with 40 %, 60 %, and 80 % of the whole data (1000 h) for training are presented in Fig. 10. Similar to Test-A, the DI-ESN has better performance than SI-ESN and the superposition effects of stack current (400 h-500 h, and 600 h-800 h) are more pronounced in Fig. 10 (b) and Fig. 10 (d). The prediction results of Test-B with different training lengths are shown in Table 5. The overall prediction performance of DI-ESN is better than SI-ESN. No matter for the SI-ESN or DI-ESN, the %Er<sub>FT</sub>, RMSE, and MAPE have a decreasing tendency with the increase of training length. Similar to Test-A, two anomaly stages exist also at 40 % and 70 % conditions from the view of RMSE and MAPE. The limiting data and the disturbance caused by the T<sub>inH2</sub> may affect the prediction result. This phenomenon weakened when more data are added to the training part (from 60%). The RUL predictions of SI-ESN and DI-ESN in the whole lifespan of Test-B are shown in Fig. 11. The PH of SI-ESN is 150 h and the PH of DI-ESN is 250 h. The results of Test-B also show that the system dynamic has a crucial effect on the prediction performance. In Fig. 11, an abnormal region exists both in SI-ESN (R1) and DI-ESN (R2) due to the cycling condition from 500 h to 800 h. This cycling condition would deteriorate the prediction results and increase the error. The abnormal region R1 of DI-ESN is smaller than the abnormal region R2 of SI-ESN. It means that the DI-ESN has stronger robustness in the dynamic load condition.

In Test-C, the RUL predictions with 50 %, 65 %, and 75 % of the whole data (405 h) for training are presented in Fig. 12. The superposition effects of stack current are more pronounced in Fig. 12 (b) and Fig. 12 (d). The prediction results of Test-C with different training lengths are shown in Table 6. The RUL predictions of SI-ESN and DI-ESN in the whole lifespan of Test-B are shown in Fig. 13. The results of Test-C are similar to Test-A and Test-B. They also prove that the DI-ESN structure performs better than SI-ESN. The PH of SI-ESN is 33 h and the PH of DI-ESN is 63 h. These three dynamic tests also show that the evaluation error is large at the BoL because of a lack of training data. This phenomenon is validated between 40% - 60% training length both in SI-ESN and DI-ESN structure. With the increasing of training data, the magnitude of the error decreases as the prediction time gets closer to EoL.
<figure><img src="image_22.png" /><figcaption>Fig. 10. The RUL prediction of Test-B: (a) SI-ESN with 40% training, (b) DI-ESN with 40% training, (c) SI-ESN with 60% training, (d)</figcaption></figure>

<caption>DI-ESN with 60% training, (e) SI-ESN with 80% training, (f) DI-ESN with 80% training.</caption>

<table><caption>Table 5. Prediction results comparison of Test-B</caption><thead><tr><th>ESN type</th><th>Training length (%)</th><th>Actual RUL (h)</th><th>Prediction RUL (h)</th><th>%Er<sub>FT</sub> (%)</th><th>RMSE</th><th>MAPE</th></tr></thead><tbody><tr><td rowspan="6">SI-ESN</td><td>40</td><td>600</td><td>468</td><td>22.0</td><td>0.01331</td><td>0.12065</td></tr><tr><td>50</td><td>500</td><td>497</td><td>0.6</td><td>0.00879</td><td>0.09849</td></tr><tr><td>60</td><td>400</td><td>174</td><td>56.5</td><td>0.02819</td><td>0.19113</td></tr><tr><td>70</td><td>300</td><td>110</td><td>63.3</td><td>0.02422</td><td>0.18926</td></tr><tr><td>80</td><td>200</td><td>122</td><td>39.0</td><td>0.01244</td><td>0.10711</td></tr><tr><td>90</td><td>100</td><td>90</td><td>10.0</td><td>0.00885</td><td>0.05098</td></tr><tr><td rowspan="6">DI-ESN</td><td>40</td><td>600</td><td>390</td><td>35.0</td><td>0.02065</td><td>0.14497</td></tr><tr><td>50</td><td>500</td><td>438</td><td>12.4</td><td>0.00788</td><td>0.07125</td></tr><tr><td>60</td><td>400</td><td>342</td><td>14.5</td><td>0.00810</td><td>0.06666</td></tr><tr><td>70</td><td>300</td><td>215</td><td>28.3</td><td>0.01197</td><td>0.10634</td></tr><tr><td>80</td><td>200</td><td>138</td><td>9.0</td><td>0.00631</td><td>0.03575</td></tr><tr><td>90</td><td>100</td><td>40</td><td>0.0</td><td>0.00440</td><td>0.02720</td></tr></tbody></table>
<figure><img src="image_23.png" /><figcaption>Fig. 11. The RUL results of SI-ESN and DI-ESN with 95% probability bounds in Test-B</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 12. The RUL prediction of Test-C: (a) SI-ESN with 50% training, (b) DI-ESN with 50% training, (c) SI-ESN with 65% training, (d) DI-ESN with 65% training, (e) SI-ESN with 75% training, (f) DI-ESN with 75% training.</figcaption></figure>

1

2

3

4

5
<table><caption>Table 6. Prediction results comparison of Test-C</caption><thead><tr><th>ESN type</th><th>Training length (%)</th><th>Actual RUL (h)</th><th>Prediction RUL (h)</th><th>%Er<sub>FT</sub> (%)</th><th>RMSE</th><th>MAPE</th></tr></thead><tbody><tr><td rowspan="6">SI-ESN</td><td>30</td><td>283.5</td><td>116.5</td><td>58.9</td><td>0.13765</td><td>0.54349</td></tr><tr><td>40</td><td>243</td><td>127</td><td>47.7</td><td>0.03882</td><td>0.27855</td></tr><tr><td>50</td><td>202.5</td><td>99.5</td><td>50.9</td><td>0.03396</td><td>0.26779</td></tr><tr><td>60</td><td>162</td><td>83</td><td>48.8</td><td>0.02254</td><td>0.20479</td></tr><tr><td>70</td><td>121.5</td><td>75</td><td>38.7</td><td>0.01407</td><td>0.15470</td></tr><tr><td>80</td><td>81</td><td>84</td><td>-3.7</td><td>0.00165</td><td>0.02045</td></tr><tr><td rowspan="6">DI-ESN</td><td>30</td><td>283.5</td><td>255</td><td>10.2</td><td>0.01100</td><td>0.26457</td></tr><tr><td>40</td><td>243</td><td>183</td><td>24.7</td><td>0.00953</td><td>0.08614</td></tr><tr><td>50</td><td>202.5</td><td>139</td><td>31.6</td><td>0.01199</td><td>0.11355</td></tr><tr><td>60</td><td>162</td><td>114</td><td>29.6</td><td>0.00937</td><td>0.09378</td></tr><tr><td>70</td><td>121.5</td><td>82</td><td>32.9</td><td>0.00976</td><td>0.11173</td></tr><tr><td>80</td><td>81</td><td>83</td><td>-2.5</td><td>0.00182</td><td>0.02273</td></tr></tbody></table>

<figure><img src="image_25.png" /><figcaption>Fig. 13. The RUL results of SI-ESN and DI-ESN with 95% probability bounds in Test-C</figcaption></figure>

2

3

4

## 4.2 Sensibility analysis of parameters

5

6

7

8

9

10

11

12

13

In industrial applications, parameter design is one of the bottlenecks of ESN's implementation and the sensitivity analysis is the basis of parameter design. There are three kinds of parameters in ESN: assigned parameters ($W_{in}$, $W$, $W_{fb}$, $K$, $L$), adjustable parameters ($N$, $\rho$, $\alpha$, $\beta$), and calculated parameter ($W_{out}$) [35]. Determining these parameters reasonably is significant to improve the dynamic performance of the reservoir, and a detailed parameter manually setting guide is presented in [54]. In general, the number of input $K$ and output $L$ is determined by the requirement of the user. The input weight matrix $W_{in}$, internal weight matrix $W$, and feedback weight matrix $W_{fb}$ are generated randomly (e.g., Gaussian distributions) and regarded to be fixed once they are assigned. In fact, the dynamic characteristics of ESN have slight differences when using different generation methods, and these assigned weight matrixes are optimized in [55]. A larger value of reservoir neurons $N$ is better for finding the linear combination between the input and output signals. The number of $N$ would increase with the task's
complexity, and using a reservoir containing more than $10^4$ neurons is a common phenomenon. The spectral radius $\rho$ is the maximal absolute eigenvalue of the $\mathbf{W}$, and it is always less than 1 to ensure the Echo State Property (ESP) for the zero-input case. Enough evidence also proves that smaller than 1 is not the necessary condition of ESP especially for the nonzero input and the reservoir could also have better properties when extending the value of $\rho$. The leaking rate $\alpha$ indicates the dynamic performance of the reservoir, and the echo property becomes more important when the $\alpha$ becomes less important. In other words, a small value of $\alpha$ represents that the value of the current step has a great influence on the next step. The parameters of $N$, $\rho$, and $\alpha$ are analyzed in [56]. Nevertheless, the inputs of ESN are the Relative Power-loss Rate (RPLR) and scheduled stack current ($I_s$) in the PEMFC lifespan prediction task, and the $\rho$ ($\rho \in (0,1)$) is better to be expended to a larger range ($\rho \in (0,1.5)$). The training part aims to calculate the output weight matrix $\mathbf{W}_{out}$, and the regularization parameter $\beta$ has an important direct effect on the calculation result. Compared with the $N$, analyzing the $\beta$ has a greater practical significance for the prediction results. So the leaking rate $\alpha$ ($\alpha \in (0,1)$), spectral radius $\rho$ ($\rho \in (0,1.5)$), and regularization parameter $\beta$ ($\beta \in (8 \times 10^{-3}, 8 \times 10^{-1})$) are selected as the key parameters in this RUL prediction task.

The objectives of parameter analysis are the quantification of the impact of the parameters on the model output and studying the interactions among the parameters. Analysis of Variance (ANOVA) is a collection of statistical models and procedures to compare the effects of different variables [57]. The ANOVA-representation of integrable function $f(x)$ in the $q$-dimensional space ($\mathbb{I}^q$) is shown as

$$ 
\begin{aligned}
f(x) &= f_0 + \sum_{s=1}^{n} \sum_{1 \le i_1 < \dots < i_s \le q} f_{i_1 \dots i_s}(x_{i_1}, \dots, x_{i_s}) \\
&= f_0 + \sum_i f_i(x_i) + \sum_{1 \le i < j \le q} f_{ij}(x_i, x_j) + \dots + f_{12 \dots q}(x_1, \dots, x_q)
\end{aligned} 
\quad (21) 
$$

Where $x$ are the input variables, $f(x)$ are the output variables, and $1 \le i < \dots < j \le q$. The premise of Eq. (21) is shown as

$$ \int_0^1 f_{i_1 \dots i_s}(x_{i_1}, \dots, x_{i_s}) dx_k = 0 \quad \text{for } k = i_1, \dots, i_s \quad (22) $$

Indeed, the members in Eq. (21) are orthogonal and the form of them can be expressed as
$$ egin{cases} 
\displaystyle \int f(x)dx = f_0 \\ 
\displaystyle \int f(x) \prod_{k \neq i} dx_k = f_0 + f_i(x_i) \\ 
\displaystyle \int f(x) \prod_{k \neq i, j} dx_k = f_0 + f_i(x_i) + f_j(x_j) + f_{i,j}(x_i, x_j) 
\end{cases} 
\quad (23) $$

The detailed description of ANOVA is shown in [58]. The Root Mean Square Error (RMSE) is used to estimate the prediction performance. A smaller RMSE indicates that the prediction has a preferable performance. The process of parameters sensibility analysis is as follows:

Step1: the data of Test-A and the SI-ESN structure are chosen as the specimen of parameter analysis, and 60% are used for training and 40% are used for prediction.

Step2: setting the key parameters ($\alpha$, $\rho$, $\beta$) by the trial-and-error method, then the ANOVA analysis is used to compute the effects and the contributions of these parameters on the prediction results.

Step3: based on the results of ANOVA analysis, the parameters of different training length of Test-A and the parameters in the other two tests (Test-A and Test-B) with different training structures can also be determined.

The parameters are divided into three different levels in step2: leaking rate $\alpha$ (0.3 for low level, 0.6 for middle level, and 0.9 for high level), spectral radius $\rho$ (0.5 for low level, 1.0 for middle level, and 1.5 for high level), and regularization parameter $\beta$ ($8 \times 10^{-3}$ for low level, $8 \times 10^{-2}$ for middle level, and $8 \times 10^{-1}$ for high level). Each experience is the result of a single simulation of a specific combination of parameters. There are thus $3^3=27$ experiment results. In Test-A, the results of different combinations are shown in Table 7 and the influence of different parameters is shown in Table 8.

<table><caption>Table 7. Prediction results of different combinations</caption><thead><tr><th>Parameters</th><th>RMSE</th><th>Parameters</th><th>RMSE</th><th>Parameters</th><th>RMSE</th></tr></thead><tbody><tr><td>α<sub>low</sub> + ρ<sub>low</sub> + β<sub>low</sub></td><td>0.00688</td><td>α<sub>mid</sub> + ρ<sub>low</sub> + β<sub>low</sub></td><td>0.00899</td><td>α<sub>hig</sub> + ρ<sub>low</sub> + β<sub>low</sub></td><td>0.01000</td></tr><tr><td>α<sub>low</sub> + ρ<sub>low</sub> + β<sub>mid</sub></td><td>0.01001</td><td>α<sub>mid</sub> + ρ<sub>low</sub> + β<sub>mid</sub></td><td>0.01098</td><td>α<sub>hig</sub> + ρ<sub>low</sub> + β<sub>mid</sub></td><td>0.01053</td></tr><tr><td>α<sub>low</sub> + ρ<sub>low</sub> + β<sub>hig</sub></td><td>0.00747</td><td>α<sub>mid</sub> + ρ<sub>low</sub> + β<sub>hig</sub></td><td>0.00770</td><td>α<sub>hig</sub> + ρ<sub>low</sub> + β<sub>hig</sub></td><td>0.00829</td></tr><tr><td>α<sub>low</sub> + ρ<sub>mid</sub> + β<sub>low</sub></td><td>0.00518</td><td>α<sub>mid</sub> + ρ<sub>mid</sub> + β<sub>low</sub></td><td>0.00525</td><td>α<sub>hig</sub> + ρ<sub>mid</sub> + β<sub>low</sub></td><td>0.00559</td></tr><tr><td>α<sub>low</sub> + ρ<sub>mid</sub> + β<sub>mid</sub></td><td>0.00691</td><td>α<sub>mid</sub> + ρ<sub>mid</sub> + β<sub>mid</sub></td><td>0.00539</td><td>α<sub>hig</sub> + ρ<sub>mid</sub> + β<sub>mid</sub></td><td><b>0.00498</b></td></tr><tr><td>α<sub>low</sub> + ρ<sub>mid</sub> + β<sub>hig</sub></td><td>0.00616</td><td>α<sub>mid</sub> + ρ<sub>mid</sub> + β<sub>hig</sub></td><td>0.00749</td><td>α<sub>hig</sub> + ρ<sub>mid</sub> + β<sub>hig</sub></td><td>0.00623</td></tr><tr><td>α<sub>low</sub> + ρ<sub>hig</sub> + β<sub>low</sub></td><td>0.01782</td><td>α<sub>mid</sub> + ρ<sub>hig</sub> + β<sub>low</sub></td><td>0.04008</td><td>α<sub>hig</sub> + ρ<sub>hig</sub> + β<sub>low</sub></td><td>0.06222</td></tr><tr><td>α<sub>low</sub> + ρ<sub>hig</sub> + β<sub>mid</sub></td><td>0.01017</td><td>α<sub>mid</sub> + ρ<sub>hig</sub> + β<sub>mid</sub></td><td>0.00950</td><td>α<sub>hig</sub> + ρ<sub>hig</sub> + β<sub>mid</sub></td><td>0.00835</td></tr><tr><td>α<sub>low</sub> + ρ<sub>hig</sub> + β<sub>hig</sub></td><td>0.01136</td><td>α<sub>mid</sub> + ρ<sub>hig</sub> + β<sub>hig</sub></td><td>0.00509</td><td>α<sub>hig</sub> + ρ<sub>hig</sub> + β<sub>hig</sub></td><td>0.00752</td></tr></tbody></table>

<caption>Table 8. Influence of different parameters</caption>
<table><thead><tr><td>Parameters</td><td>Influence (%)</td></tr></thead><tbody><tr><td>Leaking rate α</td><td>2.59</td></tr><tr><td>Spectral radius ρ</td><td>22.28</td></tr><tr><td>Regularization β</td><td>15.54</td></tr><tr><td>Interaction α * ρ</td><td>4.15</td></tr><tr><td>Interaction α * β</td><td>7.51</td></tr><tr><td>Interaction ρ * β</td><td>35.75</td></tr><tr><td>Error</td><td>12.18</td></tr><tr><td>Total</td><td>100</td></tr></tbody></table>

Results represent that setting the leaking rate $α$ at a high level, the spectral radius $ρ$ at a middle level, and the regularization parameter at a middle level could lead to a minimum RMSE (0.00498). They also show that the spectral radius $ρ$ is the most influential in these three key factors with a contribution equal to 22.28% of the total variance. The effect of the regularization parameter $β$ is the second most significant factor with a contribution equal to 15.54%. Any two of the three parameters have different interaction values. The most influential interaction parameters for this task are the spectral radius $ρ$ and regularization parameter $β$ with a contribution of 35.75%. Regarding the best combination in Table 7 as a benchmark could improve the adjusting efficiency of these three parameters in Test-B and Test-C. The key parameters in Test-A and Test-B are shown in Table 9. Thanks to the ANOVA analysis, influences of different parameters and their interactions are evaluated. And the trial-and-error research space of the parameters is reduced at the same time.

<table><caption>Table 9. Parameter levels and the typical values in different tests</caption><thead><tr><td></td><td>Leaking rate α</td><td>Spectral radius ρ</td><td>Regularization parameter β</td></tr></thead><tbody><tr><td>Low</td><td>0-0.3</td><td>0-0.5</td><td>&lt;8×10<sup>-3</sup></td></tr><tr><td>Middle</td><td>0.3-0.6</td><td>0.5-1.0</td><td>8×10<sup>-3</sup>-8×10<sup>-2</sup></td></tr><tr><td>High</td><td>0.6-0.9</td><td>1.0-1.5</td><td>8×10<sup>-2</sup>-8×10<sup>-1</sup></td></tr><tr><td>Test-B</td><td>0.9</td><td>0.7</td><td>8×10<sup>-2</sup></td></tr><tr><td>Test-C</td><td>0.9</td><td>0.6</td><td>8×10<sup>-2</sup></td></tr></tbody></table>

## 5. Conclusion

As a promising eco-friendly power converter device, the PEMFC stack suffers from a limited lifespan because of the degradation mechanisms. Based on the sufficient historical operating data, the data-driven prognostic methods could predict the degradation tendency without the PEMFC system behavior models. Then some actions, such as adjusting the fuel's supplying rate, changing the cooling temperature, and improving the control efficiency can be taken to extend the service life of the PEMFC systems. A dynamic health indicator named RPLR is proposed in this paper based on the polarization curve at the BoL and the parameters' continuous measurement. This dynamic health indicator can be used in practice due to its convenient extraction. The load current could reflect the system's dynamic and combining it with the RPLR could increase the dynamic property of ESN. Therefore, DI-ESN is proposed to enhance the prediction performance of lifespan under dynamic operating
conditions. In the prediction process, the output RPLR inherits historical data degradation tendency. Meanwhile, it takes into consideration the dynamics of stack current. The feasibility and effectiveness of the proposed dynamic health indicator and the improved multiple-input ESN structure are verified under three dynamic operating conditions. Different long-term prediction results with changing training length are carried out and compared. Results show that the RPLR is an efficient dynamic health indicator, and the ESN with double inputs could imitate the degradation properties more sufficiently and improve the prediction precision under the dynamic operating conditions. In general, the value of %Er<sub>FT</sub>, RMSE, and MAPE of DI-ESN are smaller than those of SI-ESN. The prediction horizon of DI-ESN is longer than that of SI-ESN in each dynamic test, and thus more time is given to the user to anticipate maintenance actions before failure. Besides, the effects and the contributions of the various parameters (leaking rate α, spectral radius ρ, regularization parameter β) are analyzed based on the ANOVA and multiple comparisons. In the next-step study, the prognostics method able to characterize the multi-time scale degradation of the PEMFC systems under various operating conditions will be investigated.

12

13

Acknowledgment

This work has been supported by the EIPHI Graduate school (contract "ANR-17-EURE-0002") and the Region Bourgogne

14

15

Franche-Comté, European Commission H2020 grant PANDA (H2020-LC-GV-2018), EU Grant agreement No: 824256, and

16

the ANR project BIPHOPROC_2 (ANR-14-OHRI-0018).

17

References

[1] Z. Hu, L. Xu, J. Li, M. Ouyang, Z. Song, and H. Huang, "A reconstructed fuel cell life-prediction model for a fuel cell hybrid city bus," Energy Conversion and Management, vol. 156, pp. 723-732, Jan. 2018.

[2] L. Vichard, R. Petrone, F. Harel, A. Ravey, P. Venet, and D. Hissel, "Long term durability test of open-cathode fuel cell system under actual operating conditions," Energy Conversion and Management, vol. 212, p. 112813, May. 2020.

[3] Q. Li, B. Su, Y. Pu, Y. Han, T. Wang, L. Yin, et al., "A state machine control based on equivalent consumption minimization for fuel cell/ supercapacitor hybrid tramway," IEEE Transactions on Transportation Electrification, vol. 5, pp. 552-564, May. 2019.

[4] Y. Zhou, A. Ravey, and M.-C. Péra, "Multi-objective energy management for fuel cell electric vehicles using online-learning enhanced Markov speed predictor," Energy Conversion and Management, vol. 213, p. 112821, Jun. 2020.

[5] H. Q. Nguyen and B. Shabani, "Proton exchange membrane fuel cells heat recovery opportunities for combined heating/cooling and power applications," Energy Conversion and Management, vol. 204, p. 112328, Jan. 2020.

[6] Q. Li, W. Yang, L. Yin, and W. Chen, "Real-time implementation of maximum net power strategy based on sliding mode variable structure control for proton-exchange membrane fuel cell system," IEEE Transactions on Transportation
1 Electrification, vol. 6, pp. 288-297, Feb. 2020.

2 [7] U.S. Department of Energy (DOE). Fuel Cells, 2016 (Update May 2017).

3 https://www.energy.gov/sites/prod/files/2017/05/f34/fcto_myrdd_fuel CELLS.pdf

4 [8] T. Wang, Q. Li, X. Wang, Y. Qiu, M. Liu, X. Meng, et al., "An optimized energy management strategy for fuel cell
5 hybrid power system based on maximum efficiency range identification," Journal of Power Sources, vol. 445, pp.
6 227333-227343, Jan. 2020.

7 [9] M. Ebrahimi and E. Derakhshan, "Design and evaluation of a micro combined cooling, heating, and power system based
8 on polymer exchange membrane fuel cell and thermoelectric cooler," Energy Conversion and Management, vol. 171, pp.
9 507-517, Sep. 2018.

10 [10] J. Lee, F. Wu, W. Zhao, M. Ghaffari, L. Liao, and D. Siegel, "Prognostics and health management design for rotary
11 machinery systems-reviews, methodology and applications," Mechanical Systems and Signal Processing, vol. 42, pp.
12 314-334, Jan. 2014.

13 [11] T. Wang, Q. Li, H. Yang, L. Yin, X. Wang, Y. Qiu, et al, "Adaptive current distribution method for parallel-connected
14 PEMFC generation system considering performance consistency," Energy Conversion and Management, vol. 196, pp.
15 866-877, Sep. 2019.

16 [12] X. Huang, Y. Li, A. Acharya, X. Sui, J. Meng, R. Teodorescu, et al., " A review of pulsed current technique for lithium-
17 ion batteries," Energies, vol. 13, p. 2458, May. 2020.

18 [13] M. Jouin, M. Bressel, S. Morando, R. Gouriveau, D. Hissel, M. C. Péra, et al., "Estimating the end-of-life of PEM fuel
19 cells: guidelines and metrics," Applied Energy, vol. 177, pp. 87-97, Sep. 2016.

20 [14] Z. Hua, Z. Zheng, F. Gao, and M. C. Péra, "Challenges of the remaining useful life prediction for proton exchange
21 membrane fuel cells " in IECON 2019 - 45th Annual Conference of the IEEE Industrial Electronics Society, pp. 6382-
22 6387, Oct. 2019.

23 [15] Z. Hua, Z. Zheng, E. Pahon, M. C. Péra, and F. Gao, "Health indicators for PEMFC systems life prediction under both
24 static and dynamic operating conditions," in IECON 2020 - the 46th Annual Conference of the IEEE Industrial
25 Electronics Society, pp. 3963-3968, Oct. 2020.

26 [16] R. Gouriveau, M. Hilairet, D. Hissel, S. Jemei, M. Jouin, E. Lechartier, et al., "IEEE phm 2014 data challenge: outline,
27 experiments, scoring of results, winners, " tech. rep., IEEE 2014 PHM Challenge, 2014.

28 [17] M. Ibrahim, N. Y. Steiner, S. Jemei, and D. Hissel, "Wavelet-based approach for online fuel cell remaining useful
29 lifetime prediction," IEEE Transactions on Industrial Electronics, vol. 63, pp. 5057-5068, Aug. 2016.

30 [18] R. Ma, Z. Li, E. Breaz, C. Liu, H. Bai, P. Briois, et al., "Data-fusion prognostics of proton exchange membrane fuel cell
1 degradation," IEEE Transactions on Industry Applications, Apr. 2019.

2 [19] M. Jouin, R. Gouriveau, D. Hissel, M. C. Péra, and N. Zerhouni, "Prognostics of PEM fuel cell in a particle filtering framework," International Journal of Hydrogen Energy, vol. 39, pp. 481-494, Jan. 2014.

3
4 [20] J. K. Kimotho, T. Meyer, and W. Sextro, "PEM fuel cell prognostics using particle filter with model parameter adaptation," in 2014 International Conference on Prognostics and Health Management, pp. 1-6, Feb, 2014.
5

6
7 [21] H. Liu, J. Chen, D. Hissel, and H. Su, "Remaining useful life estimation for proton exchange membrane fuel cells using a hybrid method," Applied Energy, vol. 237, pp. 910-919, Mar. 2019.

8
9 [22] H. Liu, J. Chen, C. Zhu, H. Su, and M. Hou, "Prognostics of proton exchange membrane fuel cells using a model-based method," IFAC-PapersOnLine, vol. 50, pp. 4757-4762, Jul. 2017.

10
11
12 [23] M. Jouin, R. Gouriveau, D. Hissel, M. Péra, and N. Zerhouni, "Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery," in 2014 International Conference on Prognostics and Health Management, pp. 1-6, Jun, 2014.

13
14
15 [24] M. Jouin, R. Gouriveau, D. Hissel, M. C. Péra, and N. Zerhouni, "Joint particle filters prognostics for proton exchange membrane fuel cell power prediction at constant current solicitation," IEEE Transactions on Reliability, vol. 65, pp. 336-349, Mar. 2016.

16
17
18
19 [25] K. Chen, S. Laghouche, and A. Djerdir, "Degradation prediction of proton exchange membrane fuel cell based on grey neural network model and particle swarm optimization," Energy Conversion and Management, vol. 195, pp. 810-818, May. 2019.

20
21 [26] Y. Wu, E. Breaz, F. Gao, D. Paire, and A. Miraoui, "Nonlinear performance degradation prediction of proton exchange membrane fuel cells using relevance vector machine," IEEE Transactions on Energy Conversion, vol. 31, pp. 1570-1582, Dec. 2016.

22
23
24 [27] R. E. Silva, R. Gouriveau, S. Jemeï, D. Hissel, L. Boulon, K. Agbossou, et al., "Proton exchange membrane fuel cell degradation prediction based on adaptive neuro-fuzzy inference systems," International Journal of Hydrogen Energy, vol. 39, pp. 11128-11144, Jul. 2014.

25
26
27 [28] K. Javed, R. Gouriveau, N. Zerhouni, D. Hissel, "Data-driven prognostics of proton exchange membrane fuel cell stack with constraint based summation-wavelet extreme learning machine," in the International Conference on Fundamentals & Development of Fuel Cells, pp. 1-8, Feb. 2015.

28
29 [29] L. Zhu and J. Chen, "Prognostics of PEM fuel cells based on Gaussian process state space models," Energy, vol. 149, pp. 63-73, Apr. 2018.

30 [30] R. Ma, T. Yang, E. Breaz, Z. Li, P. Briois, and F. Gao, "Data-driven proton exchange membrane fuel cell degradation
1 predication through deep learning method," Applied Energy, vol. 231, pp. 102-115, Dec. 2018.

2 [31] S. Morando, S. Jemei, R. Gouriveau, N. Zerhouni, and D. Hissel, "Fuel cells prognostics using echo state network," in IECON 2013 - 39th Annual Conference of the IEEE Industrial Electronics Society, pp. 1632-1637, Nov. 2013.

3 [32] S. Morando, S. Jemei, D. Hissel, R. Gouriveau, and N. Zerhouni, "Proton exchange membrane fuel cell ageing forecasting algorithm based on echo state network," International Journal of Hydrogen Energy, vol. 42, pp. 1472-1480, Jan. 2017.

4 [33] H. Jaeger, "The 'echo state' approach to analysing and training recurrent neural networks-with an erratum note," Fraunhofer Institute for Autonomous Intelligent Systems, 2010.

5 [34] H. Jaeger, "Tutorial on training recurrent neural networks, covering BPTT, RURL, EKF and the 'Echo State Network' approach," Technical report GMD-German National Research Center for Information Technology, 2002.

10 [35] Z. Hua, Z. Zheng, M. C. Péra, and F. Gao, "Remaining useful life prediction of PEMFC systems based on the multi-input echo state network," Applied Energy, vol. 265, p. 114791, May. 2020. https://doi.org/10.1016/j.apenergy.2020.114791.

15 [36] Z. Hua, Z. Zheng, M. C. Péra, and F. Gao, "Data-driven prognostics for PEMFC systems by different echo state network prediction structures" in ITEC 2020 - IEEE Transportation Electrification Conference and Exposition, Aug. 2020.

16 [37] X. Zhang, P. Pisu, "An unscented Kalman filter based approach for the health-monitoring and prognostics of a electrolyte membrane fuel cell polymer," in 2012 Annual Conference of Prognostics and Health Management Society, pp. 1-9, Jan. 2012.

20 [38] D. Zhang, C. Cadet, N. Y. Steiner, F. Druart, and C. Bérenguer, "PHM-oriented degradation indicators for batteries and fuel cells," Fuel Cells, Mar. 2017.

21 [39] J. Chen, D. Zhou, C. Lyu, and C. Lu, "A novel health indicator for PEMFC state of health estimation and remaining useful life prediction," International Journal of Hydrogen Energy, vol. 42, pp. 20230-20238, Aug. 2017.

23 [40] D. Zhang, P. Baraldi, C. Cadet, N. Y. Steiner, C. Bérenguer, and E. Zio, "An ensemble of models for integrating dependent sources of information for the prognosis of the remaining useful life of proton exchange membrane fuel cells," Mechanical Systems and Signal Processing, vol. 124, pp. 479-501, Jun. 2019.

26 [41] 2013-2016: ANR PROGELEC PROPICE: Prognostics et Health Management de systèmes Piles à Combustible de type PEMFC, available on http://www.propice.ens2m.fr/.

28 [42] E. Pahon, S. Morando, R. Petrone, M. C. Péra, D. Hissel, N. Y. Steiner, et al., "Long-term tests duration reduction for PEMFC µ-CHP application," International Journal of Hydrogen Energy, vol. 42, pp. 1527-1533, Jan. 2017.

30 [43] Z. Li, S. Jemei, R. Gouriveau, D. Hissel, and N. Zerhouni, "Remaining useful life estimation for PEMFC in dynamic
1 operating conditions," in 2016 IEEE Vehicle Power and Propulsion Conference (VPPC), pp. 1-6, Oct. 2016.

2 [44] Z. Li, Z. Zheng, and R. Outbib, "Adaptive prognostic of fuel cells by implementing ensemble Echo State Networks in
time varying model space," IEEE Transactions on Industrial Electronics, vol. 67, pp. 379-389, Jan. 2020.

4 [45] M. Bressel, M. Hilairet, D. Hissel, and B. O. Bouamama, "Extended Kalman filter for prognostic of proton exchange
membrane fuel cell," Applied Energy, vol. 164, pp. 220-227, Feb. 2016.

6 [46] M. Bressel, M. Hilairet, D. Hissel, and B. O. Bouamama, "Remaining useful life prediction and uncertainty
quantification of proton exchange membrane fuel cell under variable load," IEEE Transactions on Industrial Electronics,
vol. 63, pp. 2569-2577, Apr. 2016.

9 [47] B. Liu, H. Chen, T. Zhang, and P. Pei, "A vehicular proton exchange membrane fuel cell system co-simulation modeling
method based on the stack internal distribution parameters monitoring," Energy Conversion and Management, vol. 197,
p. 111898, Jan. 2019.

12 [48] J. Spendelow, J. Marcinkoski, and D. Papageorgopoulos, "Micro CHP fuel cell system targets," DOE Hydrogen and
Fuel Cells Program Record, Mar. 2012.

14 [49] J. Kurtz, S. Sprik, G. Saur, and S. Onorato, "Fuel cell electric vehicle durability and fuel cell performance," Technical
report NREL, Mar. 2019.

16 [50] H. Jaeger, M. Lukoševičius, D. Popovici, and U. Siewert, "Optimization and applications of echo state networks with
leaky-integrator neurons," Neural Networks, vol. 20, pp. 335-352, Apr. 2007.

18 [51] S. Zhong, X. Xie, L. Lin, and F. Wang, "Genetic algorithm optimized double-reservoir echo state network for multi-
regime time series prediction," Neurocomputing, vol. 238, pp. 191-204, May. 2017.

20 [52] R. Gouriveau and N. Zerhouni, "Connexionist-systems-based long term prediction approaches for prognostics," IEEE
Transactions on Reliability, vol. 61, pp. 909-920, Dec. 2012.

22 [53] A. Saxena, J. Celaya, B. Saha, S. Saha, and K. Goebel, "Metrics for offline evaluation of prognostic performance,"
International Journal of Prognostics and Health Management, vol. 1, pp. 2153-2648, Jan. 2010.

24 [54] M. Lukoševičius, "A practical guide to applying echo state networks," Neural Networks: Tricks of the Trade, pp. 659-
686, 2012.

26 [55] N. Chouikhi, B. Ammar, N. Rokbani, and A. M. Alimi, "PSO-based analysis of echo state network parameters for time
series forecasting," Applied Soft Computing, vol. 55, pp. 211-225, Jun. 2017.

28 [56] S. Morando, S. Jemei, D. Hissel, R. Gouriveau, and N. Zerhouni, "ANOVA method applied to proton exchange
membrane fuel cell ageing forecasting using an echo state network," Mathematics and Computers in Simulation, vol.
131, pp. 283-294, Jan. 2017.
1 [57] A. Azadeh, S. F. Ghaderi, and S. Sohrabkhani, "Annual electricity consumption forecasting by neural network in high energy consuming industrial sectors," Energy Conversion and Management, vol. 49, pp. 2272-2278, Aug. 2008.

3 [58] I. M. Sobol, "Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates," Mathematics and Computers in Simulation, vol. 55, pp. 271-280, Feb. 2001.