Dacheng Zhang

Piero Baraldi

Catherine Cadet

Nadia Yousfi-Steiner

Christophe Bérenguer
email: christophe.berenguer@grenoble-inp.fr

Enrico 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

Keywords: Proton Exchange Membrane Fuel Cell (PEMFC), Prognostics, Remaining Useful Life (RUL), Particle Filtering (PF), Ensemble, Dependent processes

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Proton Exchange Membrane Fuel Cell (PEMFC) has been considered as one of the most promising technologies for both stationary and transportation applications. However, it is not yet ready for large-scale industrial deployment and commercialization because of its limited durability [START_REF] Wang | Barriers of scaling-up fuel cells: Cost, durability and reliability[END_REF]. Prognostics and Health Management (PHM) approaches offer solutions to estimate the State Of Health (SOH) of fuel cell stacks and to predict their Remaining Useful Life (RUL). Prognostic results can be helpful in making decisions on maintenance scheduling and control strategy to properly operate components and systems [START_REF] Jouin | Prognostics and Health Management of PEMFC -State of the art and remaining challenges[END_REF][START_REF] Javed | State of the art and taxonomy of prognostics approaches, trends of prognostics applications and open issues towards maturity at different technology readiness levels[END_REF].

The degradation of a PEMFC stack can be assessed from different measurements. The stack voltage (or power) is the most commonly used indicator in literature [START_REF] Kimotho | PEM fuel cell prognostics using particle filter with model parameter adaptation[END_REF][START_REF] Jouin | Degradations analysis and aging modeling for health assessment and prognostics of PEMFC[END_REF][START_REF] Zhang | Some Improvements of Particle Filtering Based Prognosis for PEM Fuel Cells[END_REF][START_REF] Zhou | Degradation prediction of PEM fuel cell using a moving window based hybrid prognostic approach[END_REF]. Since this indicator can be continuously monitored in PEMFC applications, it is a good candidate for online prognostic purposes. However, the coexistence of reversible and irreversible degradation phenomena significantly limits its accuracy [START_REF] Zhang | Some Improvements of Particle Filtering Based Prognosis for PEM Fuel Cells[END_REF]. Degradation information can also be obtained by estimating the SOH using characterization measurements, such as Electrochemical Impedance Spectroscopy (EIS) [START_REF] Kim | A degenerated equivalent circuit model and hybrid prediction for state-of-health (SOH) of PEM fuel cell[END_REF] and polarization curves [START_REF] Lechartier | Proton exchange membrane fuel cell behavioral model suitable for prognostics[END_REF][START_REF] Bressel | Remaining useful life prediction and uncertainty quantification of proton exchange membrane fuel cell under variable load[END_REF]. Although the SOH typically provides an accurate estimation of the PEMFC degradation, it can only be measured parsimoniously given that this characterization is intrusive to the stack performance and may introduce additional complex degradation phenomena [START_REF] Zhang | Some Improvements of Particle Filtering Based Prognosis for PEM Fuel Cells[END_REF]. The contribution developed in this work aims at proposing a way of using jointly such two complementary deterioration indicators to elaborate a RUL prognosis result.

Prognostic approaches are typically distinguished into model-based, data-driven and hybrid approaches [START_REF] Kan | A review on prognostic techniques for non-stationary and non-linear rotating systems[END_REF][START_REF] Sutharssan | A review on prognostics and health monitoring of proton exchange membrane fuel cell[END_REF][START_REF] Tahan | Performance-based health monitoring, diagnostics and prognostics for condition-based maintenance of gas turbines: A review[END_REF]. Model-based approaches use mathematical equations to describe degradation phenomena and predict failure times. They can be applied only when knowledge about failure mechanisms, material properties and external loading is available [START_REF] Jouin | Degradations analysis and aging modeling for health assessment and prognostics of PEMFC[END_REF]. These approaches are used in specific applications, where an accurate analytic description of the system behavior has been developed. However, the necessary knowledge is not always available or mature, and for a real system, it is difficult (or even impossible) to obtain a degradation model in analytic form to integrate the degradation phenomena. Data-driven approaches such as those based on neural networks [START_REF] Javed | Prognostics of Proton Exchange Membrane Fuel Cells stack using an ensemble of constraints based connectionist networks[END_REF] and other meta-modeling techniques, do not require the availability of analytical models of the system behavior. They are often considered as "black boxes" because the behavior of the system is directly learned from historical data. They are usually easy to implement, as they do not explicitly model the links between internal phenomena and external observations. Data-driven approaches are, therefore, flexible to different problems but impose a high cost of data collection [START_REF] Si | Remaining useful life estimation -A review on the statistical data driven approaches[END_REF][START_REF] Tsui | Prognostics and health management: A review on data driven approaches[END_REF]. One of the main limitations of data-driven approaches lies in the requirement of training data, i.e. historical data used to infer correlations, establish patterns and evaluate data trends leading to failure. Hybrid approaches are the combination of the two previous types. They are based on physical equations, whose parameters change over time and are estimated by data-driven learning [START_REF] Zhou | Degradation prediction of PEM fuel cell using a moving window based hybrid prognostic approach[END_REF][START_REF] Liao | Review of hybrid prognostics approaches for remaining useful life prediction of engineered systems, and an application to battery life prediction[END_REF][START_REF] Liao | A hybrid framework combining data-driven and model-based methods for system remaining useful life prediction[END_REF]. The advantage of combining data-driven approaches (such as degradation trend regression or Artificial Neural Networks) with physical model-based approaches is that the data limitation and lack of knowledge can be mitigated. In [START_REF] Jha | Particle filter based hybrid prognostics for health monitoring of uncertain systems in bond graph framework[END_REF], RUL predictions are carried out using statistical degradation model obtained built on real degradation tests by a Bond-graph technique. The prognostics problem is formulated as the joint state-parameter estimation problem within a Particle Filtering framework where estimations of state of health (SOH) is obtained in probabilistic terms. In [START_REF] Zhou | Degradation prediction of PEM fuel cell using a moving window based hybrid prognostic approach[END_REF], an innovative robust prediction algorithm for performance degradation of PEMFC is proposed based on the combination of a degradation trend model and a Nonlinear Auto Regressive Neural Network (ANRNN) model. A recent review on PHM for PEMFC is presented in [START_REF] Sutharssan | A review on prognostics and health monitoring of proton exchange membrane fuel cell[END_REF]. Al-though several important research works have already been reported on PHM and RUL prognosis of PEMFC systems, this field is still at the early stage and needs further development. One of the main obstacles identified is the lack of available test and failure data. Another problem is that the aging and failure mechanisms of the PEMFC are not yet entirely clear, because of the influence of the different operating conditions.

In this work, we develop a hybrid approach based on an ensemble of models, which uses prediction outcomes provided by different degradation models fed by different deterioration measurements and properly combines them to provide the prognostic results. Ensemble of models have shown promising results for the prognostics of industrial systems [START_REF] Baraldi | Local fusion of an ensemble of models for the reconstruction of faulty signals[END_REF][START_REF] Baraldi | Ensemble neural network-based particle filtering for prognostics[END_REF][START_REF] Xing | An ensemble model for predicting the remaining useful performance of lithium-ion batteries[END_REF][START_REF] Al-Dahidi | Remaining useful life estimation in heterogeneous fleets working under variable operating conditions[END_REF][START_REF] Al-Dahidi | A locally adaptive ensemble approach for data-driven prognostics of heterogeneous fleets[END_REF][START_REF] Rigamonti | Ensemble of optimized echo state networks for remaining useful life prediction[END_REF]. For example in [START_REF] Al-Dahidi | Remaining useful life estimation in heterogeneous fleets working under variable operating conditions[END_REF], an ensemble approach based on a semi-Markov model and a fuzzy similarity model has been developed for the predictions of the RUL of a heterogeneous fleet of aluminum electrolytic capacitors used in electric vehicle power trains. The ensemble output is obtained by local aggregation of the outcomes of two prognostic models, assigning to each model a weight and a bias correction, which are estimated considering the models' local performances, i.e. the inaccuracies in predicting the RUL of validation patterns similar to the one under assessment. However, this kind of approach requires a sufficient amount of historical run-to-failure data, and it cannot be applied to PEMFC, due to the high cost of run-to-failure testing. Authors in [START_REF] Xing | An ensemble model for predicting the remaining useful performance of lithium-ion batteries[END_REF] presented an ensemble degradation model that integrates two regression models: an exponential model and a polynomial model for RUL prediction of Lithium-ion batteries. The ensemble model is a more accurate parametric model than each of the above two models because it takes into account on-board applications, and global and local regression characteristics. To be effective, ensemble-based systems and algorithms need diversity among the models, which can be achieved in different ways, [START_REF] Polikar | Ensemble Based Systems in Decision Making[END_REF]:

• Different training sets boosting, bagging, which is not the approach followed here, in particular because we work with time series;

• Different types of mathematical processing models;

• Same mathematical processing models but with different set of input features ; this is the case of the work presented in this paper.

Our work proposes an ensemble prognosis approach formed on the basis of different and dependent deterioration data sources, which makes it original when compared to existing ensemble prognosis approaches using rather a single data source processed through different models. More precisely, here, our proposition is to use an ensemble approach to take the best benefit from different measurements gathered at different levels in the monitored system, with a different quality and available with a different periodicity (namely here a stack voltage measurement and a SOH estimation), in order to improve the quality and the performance of the prognosis of the overall PEMFC system. Particle Filtering (PF), which has been developed recently among applications for PHM [START_REF] Jouin | Particle filter-based prognostics: Review, discussion and perspectives[END_REF], is applied so that the RUL uncertainties are also obtained.

The remainder of this paper is organized as follows. In Section 2, the problem formulation and the models for RUL prognosis are presented. Then, the Particle Filtering-based prognostic approach and the criteria for performance evaluation are presented in Section 3. The Ensembleof-models approach to RUL prediction is explained in Section 4. In Section 5, we pay particular attention to data simulation and the process dependencies, and two different procedures of data generation are explained. Finally, the results of a numerical experimentation are reported in Section 6 showing the good performance of the proposed ensemble approach for RUL prognosis.

Problem Formulation

Fuel Cell Degradation

For fuel cell systems, just like for any other technological system, performance degradation is unavoidable but can be minimized by proper operation and maintenance, based on a comprehensive understanding of degradation mechanisms. PEMFC degrades due to calendar aging, which can occur even under constant optimal conditions, start and stop cycles and inadequate operating conditions such as temperature, pressure and poor water management. In a fuel cell, the aging process reduces the component performance and modifies its material physical properties. V st (V)

Aging

Figure 1: Voltage degradation with aging (under constant current density) [START_REF] Gouriveau | IEEE PHM 2014 Data Challenge -Outline, Experiments, Scoring of results[END_REF].

During the lifetime of a PEMFC, its "health" and performance gradually deteriorate, due to irreversible physical and chemical changes, which take place with usage and with aging, until the moment the stack is no longer usable. The State Of Health (SOH) (e.g. from its Begin Of Life (BOL) status of full performance to its End Of Life (EOL) status of functional failure, i.e. its performance does not meet the desired operational standard) provides an indication (not an absolute measurement) of the performance which can be expected from the PEMFC in its current condition and of the amount of lifetime already spent by the component.

Any parameter significantly changing with age, such as cell impedance, can be used for indicating the SOH of the cell. These parameter changes are typically identified by performing characterization measurements such as polarization curves. The polarization curve describes the working performance of PEMFC. The variations of internal parameters, including physical and empirical ones, have great impact on the polarization characteristic. Figure 2 shows the variation of the polarization curves under aging. In this work, physical and empirical parameters are used to predict the performance of the fuel cell. We will consider two indicators of the PEMFC degradation: the stack voltage and the SOH. The stack voltage V st can be measured at a high frequency (≈ 0.6s) and the SOH are characterized every week in practice.

Models Description 2.2.1. Voltage Model

Being electrochemical cells, fuel cells obey to thermodynamic and kinetic laws. The static voltage of a fuel cell stack, depicted in Figure 2, is given by [START_REF] Lechartier | Proton exchange membrane fuel cell behavioral model suitable for prognostics[END_REF]:

V st = n • (E -V ohm -V act -V trans ) (1)
where V st is the stack voltage, n is the number of cells in the stack, E is the open circuit voltage (OCV), V act is the activation polarization, V ohm represents the ohmic losses (due to the electrical resistance of individual components and their contact), and V trans is the concentration polarization (due to mass transport limitation). For a stack operating at a current density j [START_REF] Larminie | Fuel cell systems explained[END_REF]:

V st = n • E -r • j -A • ln( j j 0 ) -m 1 • exp(m 2 • j) ( 2
)
where r is the internal resistance, j is the operating current density, A is the Tafel coefficient, j 0 is the exchange current density, m 1 and m 2 are the mass-transfer constants. Considering different current current density values, a static polarization curve is obtained.

SOH Degradation Model

A limitation of the stack voltage is that it does not allow separating the effect of the load variation, which causes current density variations, from that of the stack inner degradation, which influences the OCV [START_REF] Zhang | PHM-oriented Degradation Indicators for Batteries and Fuel Cells[END_REF] and the global resistance parameters [START_REF] Bressel | Remaining useful life prediction and uncertainty quantification of proton exchange membrane fuel cell under variable load[END_REF][START_REF] Jha | Particle filter based hybrid prognostics for health monitoring of uncertain systems in bond graph framework[END_REF]. Since physical laws describing the effects of the degradation on E and r are not known, in this work we adapt linear equations for simplicity of illustration and without loss of generality of the proposed approach. The changes in the two parameters are coupled by variable γ(t), which reflects the SOH degradation:

r(t) = r 0 (1 + γ(t)) E(t) = E 0 (1 -γ(t)) (3)
where r 0 and E 0 are the initial values of r and E. Since it has been proven in [START_REF] Bressel | Remaining useful life prediction and uncertainty quantification of proton exchange membrane fuel cell under variable load[END_REF][START_REF] Jha | Particle filter based hybrid prognostics for health monitoring of uncertain systems in bond graph framework[END_REF] that the SOH indicator γ(t) can be estimated from polarization curves, in this work we assume the availability of the procedure which returns the SOH degradation estimation γ(t) from characterization measurements of the PEMFC stack. Thus, γ(t) can be taken as an input for our prognostic procedure.

Prognostic Models for RUL Prediction

Two stochastic state transition models are used for describing the SOH deterioration γ(t) and the stack voltage degradation V st (t).

Problem Statement

When the fuel cell stack experiences variable operating conditions, a single degradation indicator is not able to provide a precise and robust RUL prediction. The stack voltage does not directly measure the component degradation but it is only related to degradation symptoms, which are significantly affected by operating conditions. The SOH provides aging information but it can only be measured at low frequency in industrial applications.

In this work, we consider prognostics based on two different measurements of the stack degradation:

• An external signal, such as the stack voltage, which is easily accessible and frequently measured, but of "poor quality", i.e. its measurement is affected by significant noise.

• A signal which provides an internal characterization of the component, such as the stack SOH, which is seldom measured due to the complexity and cost of the measurement procedure that requires to take the fuel cell stack out of service for the measurements.

The objective is to combine the predicted RUL outcomes based on the two signals.

Particle Filtering-based RUL Prognosis

RUL Prognosis

The RUL is the time remaining from the current moment and the moment when the system is considered failed. As been depicted in Figure 3, degradation measurements are used to train the prognostic model during the learning phase until the prediction time t λ . Then the learned behavior of this degradation path is used to predict the future evolution with time. The End Of Life (EOL) is the time when the estimated degradation state reaches the failure threshold where the RUL can be computed. The performance of the prognostic model is typically evaluated by comparing the RUL prediction with ground truth RUL (see Appendix 1).

Particle Filtering

Particle Filtering (PF) relies on state-space description of the system evolution and observation with possibly non-linear and non Gaussian features [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. It is a recursive state estimation techniques based on a Bayesian approach [START_REF] Patwardhan | Nonlinear Bayesian state estimation: A review of recent developments[END_REF].

For our purposes, the degradation dynamics and its observations are assumed to be governed by a discrete-time state transition model:

x k = f k (x k-1 , ω k-1 , Θ k-1 ) (4a)
z k = h k (x k , ν k ) ( 4b
)
where k is the time index, x is the system state, z is the measurement, f is the degradation model (state transition function), ω is the system noise, Θ is the vector of model parameters

(Θ = [θ 1 , θ 2 , ...]
), h is the measurement model and ν is the measurement noise. Both the process noise ω k and the observation noise ν k are assumed to be sampled from a zero-mean Gaussian distribution, i.e.

ω k ∼ N (0, σ ω 2 k ) and ν k ∼ N (0, σ ν 2 k
). The PF algorithm is summarized in Algorithm 1. The approximation of the probability distribution of the system state is based on sampled particles and associated weights. Bayesian updating is processed sequentially by propagating particles carrying probabilistic information on the unknown states and model parameters. The probabilistic model for the particles propagation relies on the state transition model (4a) and the probability distribution of the process noise ω k :

1. Propagate i = 1, . . . , n particles representing the system state probability density function (PDF) from x k-1 to x k by the state transition model described in Equation (4a) (Algorithm 1, line 5 ).

2. For each particle, estimate the associated weight by calculating its likelihood given an online measurement z k (Algorithm 1, line 6 ). This gives the corresponding weight of each particle (assuming the measurement noise

ν k ∼ N (0, σ ν 2 k ) is normally distributed): L(z k |x i k , σ ν i k ) = 1 √ 2πσ ν i k exp[- 1 2 ( z k -x i k σ ν i k ) 2 ] (5)
3. Perform resampling [START_REF] Li | Resampling Methods for Particle Filtering[END_REF] to remove the particles with small weights relative to a given weight limit and replicated those with large weights (Algorithm 1, line 10 to 17 ).

Algorithm 1 Particle Filtering

1: Initialize x i 0 , σ ω i 0 , σ ν i 0 and Θ i 0 // drawn from initial uniform distributions 2: Time step k = 1 3: while x i k > F T and k ≤ k p 4:
for i = 1, . . . , n // Importance sampling:

5: Draw particles x i k ∼ p(x i k |x i k-1 , σ ω i k-1 , Θ i k-1 ) using Equation (4a) 6:
Assign weight

w i k = L(z k |x i k , σ ν i k ) using Equation (5) 7:
end for

8:

Normalize weight

w i k = w i k / n i=1 w i k 9:
Calculate the cumulative sum of normalized weights:

{Q i k } n i=1 = Cumsum {w i k } n i=1 10:
for i = 1, . . . , n // Resampling (Multinomial):

11: j = 1 12:
Draw a random value u i ∼ U(0, 1]

13:

while Q j k < u i 14:

j = j + 1 15:
end while 16:

Update state

x i k = x j k Update noises σ ω i k = σ ω j k , σ ν i k = σ ν j k Update parameters Θ i k = Θ j k 17:
end for 18:

k = k + 1 19: end while
The process is performed until no measurement is available (prediction time t λ = k p • ∆t reached).

For the RUL prediction, the posterior PDFs of the state and model parameters, given the observation sequence up to time t λ , are used to estimate the future evolution of the particles. The RUL PDF can be obtained when the particles representing the system state reach the preset failure threshold, as illustrated in Figure 3. The prognostic procedure is summarized in Algorithm 2.

Algorithm 2 RUL prediction

1: k = k p // Start
while x i k > F T 5: k = k + 1 6:
Propagate particles

x i k = f (x i k-1 , σ ω i k-1 , Θ i k-1 ) using Equation (4a) 7:
end while 8:

Estimate RU L i k = (k -k p ) • ∆t 9: end for 4. Ensemble-based Prognostic Approach

Prognostic Models for RUL Prediction

For the PF-based estimation stage, the following simplified state models are used for the SOH deterioration γ(t) and the stack voltage degradation trend V st (t):

• γ(t) represents the SOH degradation. It assumes values from 0 (healthy) to 100% (failed), following a linear model:

γ(t + 1) = c (1) (t) • γ(t) (6)
where c (1) (t) is the time-dependent SOH degradation model parameter [START_REF] Bressel | Remaining useful life prediction and uncertainty quantification of proton exchange membrane fuel cell under variable load[END_REF].

• V st (t) represents a symptom of the stack degradation, which, according to Equation ( 6), follows a linear trend:

V st (t + 1) = c (2) (t) • V st (t) (7)
where c (2) (t) is the time-dependent voltage degradation parameter

The two linear degradation models of Equation ( 6) and ( 7) are used in two different particle filtering algorithms to provide the RUL predictions RU L (m) t (m=1, 2), respectively. Model 1 in Equation ( 6) uses measurements of good quality, but not frequently acquired, whereas Model 2 in Equation ( 7) uses measurements that are regularly available, whose quality can be poor due to higher measurements noise and lower correlation with the true health states. The objective is, then, to combine the individual estimates RU L [START_REF] Wang | Barriers of scaling-up fuel cells: Cost, durability and reliability[END_REF] and RU L [START_REF] Jouin | Prognostics and Health Management of PEMFC -State of the art and remaining challenges[END_REF] , taking into account their "local" qualities.

Ensemble of Models

Fusing the outputs of an ensemble of diverse prognostic models can improve overall prediction accuracy [START_REF] Bonissone | Fast meta-models for local fusion of multiple predictive models[END_REF]. Local aggregation dynamically assigns weights to each model according to its local performance, typically evaluated on the available historical patterns [START_REF] Rigamonti | Ensemble of optimized echo state networks for remaining useful life prediction[END_REF]. For prognostics, local aggregation requires the computation of the local performances of the individual models on a set of run-to-failure degradation trajectories.

RUL * t,k
LE (2) t,k LE (1) (2) t RUL (1) t RUL (1) t,k RUL (2) t,k RUL t LE (1) t,1 LE (2) t,N LE (1) t LE (2) t

t,k Ensemble SOH k Voltage k SOH Voltage Training set k N RUL
Step 1: Similar trajectories retrieval Step 2: Weights determination Step 3: RUL predictions 1 K } w (1) t B (1) t , B (2) t w (2) t , j 1 k LE (1) t,N } } } LE (2) t,1 } Figure 4: Scheme of the proposed prognostic approach.

Figure 4 presents the scheme of the ensemble-based prognostic approach proposed in this work. As mentioned previously, we assume the availability of the measurements of the signals x (1) = γ and x (2) = V st collected during the life of K identical fuel cell stacks:

{x train k } K k=1 = {(x (1),train k , x (2)
,train k

)} K k=1 (8)
These run-to-failure trajectories form a training set, which is also used within the ensemble approach for the aggregation of the individual model outcomes. The local fusion approach for the aggregation of the individual model outcomes is based on the following steps, illustrated in Figure 4:

Step 1 : Similar trajectories retrieval -For the test trajectory x test , evaluate the most similar trajectories in the training set by computing the point-wise Euclidean distance on a time window of length L. The lower the Euclidean distance, the greater is the similiarity between the trajectories. For each model m, identify among the training trajectories the most similar to the test trajectory x test , by finding the ones with the minimum Euclidean Distance (ED) [START_REF] Rigamonti | Ensemble of optimized echo state networks for remaining useful life prediction[END_REF]:

d (m) t = min{ED x (m),train (t-L):t,k , x (m),test (t-L):t } K k=1 (9)
where d

(m) t
is the minimum ED of the test trajectory for the m th model at time t, on a time window L. Then, select the N nearest (i.e. most similar) trajectories among the K training trajectories of each measurements for future analysis. on the RUL prediction defined at time step t for the m th model as:

LE (m) t = 1 N N k=1 |RU L * t,k -RU L (m) t,k | (10)
where RU L * t,k is the corresponding ground truth RUL and RU L (m) t,k is the estimated RUL of the k th trajectory predicted by the m th model. The larger the prediction error on the RUL predicted by a model, the poorer the performance of prediction of this model, at least locally (i.e. at the considered time step t). Thus, the local error LE (m) t provides information about the (poor) performance of the m th model in the reconstruction of the patterns of the training set which are closer to the test pattern. The weight associated to each model in the ensemble should then be a decreasing function of its local error, and we use the inverse of the local error to compute the weight [START_REF] Baraldi | Local fusion of an ensemble of models for the reconstruction of faulty signals[END_REF]:

w (m) t = 1/LE (m) t M m=1 1/LE (m) t (11)
where M is the number of models (in our case M = 2). The local weights w (m) are nonnegative and sum to 1. Note that the weights are "local" in the sense that the RUL estimation RU L (m) t is evaluated at different time steps dynamically. Before aggregating the RUL predictions with their corresponding weights, a bias correction B (m) t of the i th model is subtracted:

B (m) t = 1 N N k=1 RU L * t,k -RU L (m) t,k (12)
This quantity represents the accuracy of the RUL predictions obtained by each m th model on the N selected training trajectories. The reason of introducing the bias correction is that at the early prediction stage, due to insufficient available observations, the prognostic algorithm usually provides predictions characterized by large variability. Exploiting the historical data, the average variation can be learned from the training trajectories and used as an offset.

Step 3 : RUL predictions for the different models -Predict the RULs for the test trajectory using the PF method described in Section 3, based on the M models.

Step 4 : Ensemble RUL aggregation -Aggregate RUL predictions based on the individual models and weighted based on prognostic performances:

RU L t = M m=1 w (m) t • RU L (m) t -B (m) t (13)
where RU L

t , m = 1, 2, . . . , M , is the predicted RUL of the test trajectory x test and B (m) t

is the bias correction evaluated on all N training trajectories.

The ensemble approach allows obtaining the PDF density of the predicted RUL. Various mathematical methods and approaches for combining probability distributions are discussed in [START_REF] Clemen | Aggregating Probability Distributions[END_REF]. Among them, in this work, we consider the Linear Opinion Pool (LOP), which is a common method for weighted linear combination of the experts probabilities [START_REF] Berliner | A framework for multi-model ensembling[END_REF][START_REF] Bolger | Reliability updating in linear opinion pooling for multiple decision makers[END_REF][START_REF] Bolger | Deriving the probability of a linear opinion pooling method being superior to a set of alternatives[END_REF] and it is easily understood and calculated:

p( RU L t ) = M m=1 w (m) t • p RU L (m) t (14)
where p( RU L t ) represents the merged probability distribution, and p( RU L (m) t ), represents the RUL distributions predicted by the M particle filters.

Data Generation

Because of the lack of real data available to test and validate the proposed approach, we have to resort to simulated data. Under this situation, it is necessary to clearly show how the chosen procedure to generate the simulated data allows to mimic reasonably the reality and how their generation can be controlled to test different aspects, i.e. here the variability of the considered measurements, and their dependency. Hence, we aim at simulating a realistic evolution of the signals γ(t) and V st (t), properly accounting for temporal and stack-to-stack variability, and also the dependence between the two signals. By "realistic", we mean that both signals should be correlated, but not fully equivalent nor exchangeable with respect to the degradation information they carry. Given the unavailability of real data describing the degradation of a fleet of similar PEMFC stacks, the degradation trajectories are generated by applying the physics-based models of Equation ( 2) described in Section 2.2. This procedure allows obtaining the SOH and the voltage degradation paths of similar stacks, realistically taking into account their variaibility and dependence by resorting also to stochastic processes, and in particular here to Gamma processes. The simulated degradation trajectories are then divided into a training set made by K trajectories and a test set made by J trajectories.

Gamma Process

An homogeneous Gamma process is a stochastic process with independent, non-negative increments following a Gamma distribution. If (X t ) t≥0 is a Gamma process, then:

∆X t-s = X t -X s ∼ G(α.(t -s), β) for all 0 ≥ t ≥ s ( 15
)
where X 0 = 0, with probability equal to 1, ∆X t are independent, G(α, β) denotes the Gamma distribution with shape parameter α and scale parameter β, with the following probability density function

f (x; α, β) = β α x α-1 e -βx Γ(α) (16)
Over a time interval t, the average degradation rate (slope

) is x = α • β • t, the process variance V ar = α • β 2 • t.
The Gamma process is suitable to model gradual damage monotonically accumulating over time in a sequence of tiny increments, such as wear, fatigue, corrosion, crack growth, degrading health index, etc [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF]. Thus it is used here for simulating the irreversible degrading SOH of the PEMFC stack. The choice of α and β allows one to set different values for the deterioration average and variance, and hence to model various degradation behaviors, from almost-deterministic to very chaotic. Given the degradation measurements, both parameters can be estimated using classical statistical methods, such as maximum likelihood method, moment method, Bayesian statistics method, etc.

Using a stochastic process-based (specifically here a Gamma process-based) degradation models make it possible to take both the temporal variability and the item-to-item variability into account [START_REF] Pandey | The influence of temporal uncertainty of deterioration on life-cycle management of structures[END_REF]. Several kinds of stochastic processes could have been used ; one advantage of using a Gamma process for degradation modeling is that the required mathematical calculations are relatively straightforward. The RUL can be, thus, obtained in an analytic form if necessary.

Signal Simulation 5.2.1. SOH Simulation

The degradation path of γ(t) is generated by a Gamma process, which accounts for the randomness of the degradation process. The failure threshold F T γ , here set to the value of 0.15, is obtained by estimating the internal resistance from EIS characterization [START_REF] Kim | A degenerated equivalent circuit model and hybrid prediction for state-of-health (SOH) of PEM fuel cell[END_REF]. Figure 6 shows one simulated degradation path. The average Gamma process γ classifies the type of fuel cell stack, and the variation from stack to stack is represented by drawing different realizations from γ. The average End Of Life (EOL) can be found at the time point when γ crosses the threshold F T γ :

EOL = F T γ α • β (17)

Stack Voltage Simulation

The degradation pattern of PEM fuel cells degradation in stack power is not linear; the decreasing trends are not monotonic. According to Equation (2), the stack voltage is influenced by the loading current density j, which is here simulated by a Markov process [START_REF] Gamerman | Markov Chains[END_REF]. It is used here to simulate the operating conditions during stack usage (Figure 7). In order to be as realistic as possible, the current profile is simulated to its practical value around 0.7mA/cm with 5% random variation in a test bench [START_REF] Gouriveau | IEEE PHM 2014 Data Challenge -Outline, Experiments, Scoring of results[END_REF]. From a given γ, the degradation path of V st is simulated using Equation (2) where the failure threshold F T Vst is obtained by substituting t = EOL (Figure 8). Note that this failure threshold is deduced from the failure time EOL.

Observation

The two γ and V st trajectories simulated above are considered as the ground truth. Since measurements revealed by sensors are affected by noises, we randomly sample their values by adding to the ground truth states zero-mean Gaussian noises:

γ meas = γ + N (0, σ 2 γ ) (18)
V stmeas = V st + N (0, σ 2 Vst ) (19)
where γ and V st are the system true states, γ meas and V stmeas are the measurement readings, σ γ and σ Vst are the standard deviations of those two types of measurements, respectively. Note that σ γ < σ Vst given that the SOH measurements γ meas is more precise than the voltage measurements V st .

Data Availability

As mentioned in Section 2.2, the stack voltage can be measured more frequently than the SOH degradation. Thus, the measurement data of SOH degradation are constrained such that they are available only every 100 hours, whereas the measurement data of stack voltage are available every hour. The γ and V st measurements for one single stack are shown in Figure 9. For computational convenience, ten-time steps between two successive measurements are considered for the measurements of V st .

Degradation Simulation Procedure : Introducing Dependence between γ and V st

The algorithm developed in 5.2 allows the simulation of synthetic deterioration data with temporal and item-to-item variability and with measurements uncertainty. However, it does not take into account the dependence between γ and V st . It can be seen from Equations ( 2) and (3) that V st (t) is a symptom of the degradation γ(t), the deterioration levels of the two indicators are correlated. To properly produce realistic simulations of the degradation trajectories that model the different sources of variability, randomness, and dependence between the signals, it is necessary to introduce dependence in the proposed simulation scheme and to control the level of dependence between both deterioration indicators. Two approaches have been proposed to this aim :

• Approach 1 : the two indicators V st and γ are generated from the same realization of a Gamma degradation process, with different additive noises ;

• Approach 2 : the two indicators V st and γ are simulated from two different degradation processes, dependent by construction. To this aim, a bivariate dependent Gamma process is constructed by trivariate reduction in the case of bivariate Gamma random vectors [START_REF] Devroye | Complexity questions in non-uniform random variate generation[END_REF].

The technical details of these two approaches and the complete algorithms for synthetic data simulation are given in Appendix 2.

RUL Prognosis Results & Performance Evaluation

Considering that in real industrial applications we expect to have available a limited number of PEMFC stacks degradation trajectories. We simulate 100 trajectories of which we use each type of measurement: K = 50 for training and J = 50 for testing. By performing a sensitivity analysis (See Appendix 3) regarding the prediction accuracy and the computation time, we have set the number of nearest trajectories in the training set to N = 5 and the time window for the similarity calculation to L = 100 hours.

RUL Prognosis for Data Simulation Approach 1

The variance of the degradation process, αβ 2 , depends on the choice of the Gamma process parameters α and β. It stands for the similarity in degradation behavior of identical PEMFC stacks. As being discussed in the presentation of the simulation procedure (see Appendix 2), the objective of introducing the variance is to represent stack-to-stack variability around the average behavior. Figure 10 shows three examples of degradation paths with different levels of variance: 1) low variance (α = 0.6, β = 2.5e-4); 2) medium variance (α = 0.1, β = 1.5e-3); 3) high variance (α = 0.03, β = 5.0e-3). The RUL predictions for all the degradation trajectories are carried out by the Particle Filteringbased approach described in Section 3. For each trajectory in the test set (J = 50 trajectories), the RUL predictions are made every 100 time steps with Model 1 and every 10 time steps with Model 2. The RUL predictions based on Model 1 are less frequent than the ones based on Model 2, because the measurements that feed Model 2 are intermittently taken. Thus, to have a fair comparison between the two models, the missing predictions of Model 1 are reconstructed by linear interpolation. The simulation is carried out with the data dependence generation of Approach 1 and Gamma process with medium variance. Figure 11 shows the Local Error (LE) at different prediction time steps obtained for a single test trajectory (№40).

RUL Aggregation

Figure 12 shows the corresponding weights which are dynamically assigned to the two models according to their local error evaluated at each time step. Notice that:

• Model 1 weights are larger at the beginning of the component life compared to that of Model 2. This can be justified by the fact that Model 1 is fed by more precise SOH measurements and it is not influenced by loading current variations.

• Model 2 weights are larger than those of Model 1 after approximately 600 hours. This can be justified by the fact that Model 2 is trained by using more data. Thus, its prediction performance is improved much faster than the one of Model 1, especially near the end of life when Model 1 is no longer updated due to lack of new incoming measurements. The RUL predictions based on both models are aggregated according to Equation [START_REF] Tahan | Performance-based health monitoring, diagnostics and prognostics for condition-based maintenance of gas turbines: A review[END_REF]. Figure 13 shows the RUL predictions and the aggregation for one test trajectory. The ensemble RUL predictions take advantage of the complementary behaviors of individual models. Indeed, the analysis of Figure 13 suggests that:

• The predictions provided by the two models are comparable, even if Model 18 provides more accurate RUL predictions at the early life stages of the stack №40, Model 2 provides more accurate predictions when this stack approaches the EOL.

• The ensemble of the two models allows obtaining more accurate predictions throughout the RUL predictions of stack №40 than each individual model.

RUL Uncertainty Aggregation

Figure 16 shows the 25 th and 75 th percentiles of the RUL PDF provided by the Ensemble, which is obtained by merging the RUL PDFs of Model 1 and Model 2 according to Equation ( 14), for trajectory №40. By aggregating the two PDFs, we obtain not only the RUL but also the uncertainty of the predictions, which is very important for post-prognosis maintenance decisionmaking. As expected, the prediction becomes closer to the ground truth RUL and the uncertainties (PDFs) of the ensemble become smaller when approaching the end of life. Figure 17 depicts the RUL uncertainty at different life stages: at the early prediction time of 200 hours (Figure 17.1), at half-life of 500 hours (Figure 17.2) and near the EOL of 800 hours (Figure 17.3). All three models become less spread and centered to the true RUL accuracy zone when the prediction time approaches the end of life. The less spread distribution indicates that the RUL predictions become more accurate and more precise when more observations become available, which meets our expectation.

Prognostic Performance Evaluation

The quality of the RUL predictions of the individual models and the ensemble are evaluated using the prognostic performance metrics in Table 1, which reports the average performance over J = 50 test trajectories and all t λ time steps.

The metrics used for performance evaluation in this work consist of the accuracy index Acc, the α-λ accuracy αAc, the steadiness index Std, the risk index Rsk, the precision index P rc and the coverage index Cvg, as explained in Appendix 1. The values in Table 1 suggest that: • The Ensemble shows better performance than any individual model with respect to Acc, αAc, Std and Cvg indexes.

(1) (2) (3)
• Model 1 shows better performance in Rsk index, which means that the RUL predictions based on Model 1 are early notifications. This does not mean that all early predictions are good predictions: an early notification which is too far from the true failure time leads to unnecessary maintenance, which incurs extra cost. The Rsk performance needs to be considered jointly to the accuracy indexes (Acc and αAc). The Rsk of the Ensemble is between Model 1 and Model 2, with respect to both point values and uncertainty.

• The P rc index of the Ensemble is the weakest, whereas its Cvg index is the strongest. This is due to the fact that the PDFs of the Ensemble merges Models 1 and 2 PDFs. The spread of its distribution is, thus, broader than the individual models, but it provides a larger coverage.

Above all, in this example of Approach 1 with medium variance, we can conclude that the Ensemble-based approach globally provides the best prognostic performance.

RUL Prognosis for Data Simulation Approach 2

Similarly to what has been done for the data simulation Approach 1, three different levels of process variance are simulated. Furthermore, for each level of variance, seven different levels of processes dependence between γ(t) and V st are considered to represent the underneath correlation between the two signals.

Parameters Used for the Simulated Examples

The parameters used for the generation of the simulated examples of dependent Gamma processes are reported in Table 2. The correlation coefficient ρ (0 ≤ ρ ≤ ρ max = min(α 1 ,α 2 )

√ α 1 α 2 ) indicates the dependence level of the two final degradation processes after the trivariate reduction.

Acc(gain) = Acc Ensemble -max (Acc M odel1 , Acc M odel2 ) max (Acc M odel1 , Acc M odel2 ) Std(gain) = min (Std M odel1 , Std M odel2 ) -Std Ensemble min (Std M odel1 , Std M odel2 ) (20)
The values indicate the improvements of the Ensemble with respect to the best between Model 1 and Model 2. Gains above 0 indicate that the Ensemble performance is more satifactory than that of the individual models. With respect to the point values, notice that:

• Considering the Acc, αAc and Std metrics, the Ensemble always outperforms any of the individual models. Therefore, we can conclude that the Ensemble is more accurate than Models 1 and 2. The larger the process variance, the larger the gain in performance for the Ensemble approach. • Similar to the case of data simulation Approach 1, the Rsk index of the Ensemble tends to decrease, which means that the Ensemble provides RUL predictions exceeding the ground truth RUL, even though they are located in the accuracy zone. For the uncertainties aggregation, the analysis of Figure 19 indicates that: • Considering the Acc, αAc and Std metrics, the Ensemble model with any process dependence outperforms any of the individual models. The gains of Std are nearly the same as the one with point values. • The Rsk index for low variance processes is sometimes improved. It can also be noticed that this index is better than the one with point values aggregation, which indicates that with complete information (uncertainties), some "risky" predictions can be avoided. • Not surprisingly, the P rc index of the Ensemble is the weakest and, on the other hand, its

Cvg index is the strongest. This can be explained by the fact that the RUL PDF obtained by the ensemble approach results from the merging of of both RUL PDFs obtained with Model 1 and Model 2. The spread of the distribution is, thus, broader than those of the individual models, which provides for a larger coverage. Hence, the ensemble models approach can largely improve the prognostic performance for degradation processes with different variances and different dependencies.

Conclusion

In this work, the coexistence of two different sources of information on the degradation of a component, characterized by different levels of accuracy and acquisition rates, has been considered. We developed an ensemble approach that combines the RUL predictions from two different sources of information on the system deterioration, gathered at different levels. The RUL predictions of both models are dynamically aggregated according to their local weights estimated considering the prognostic performance evaluated on a set of historical data. The method has been applied to the prediction of the RUL of simulated PEMFC stack SOH and voltage degradation signals. The results show that the prediction accuracy is improved.

This work can be extended in several directions:

• In our approach, the weights are computed on the basis of offline historical run-to-fail data and a significant amount of historical data can thus be required, which makes it impossible to validate the proposed approach on a real PEFMC test-bench because of the cost of obtaining enough run-to-fail data. To improve the applicability of this approach it can be worthwhile to investigate other solutions for the weights definition, that are less dependent on the quantity of historical data.

• As already explained, ensemble approaches aims at taking advantage of the diversity in the data, processing methods and models. In this work, we have mainly investigated the interest of using deterioration data from different sources, with different characteristics, to improve the prognosis performance. It could be interesting to couple this use of different data sources with different ways of modeling the deterioration evolution (both physics-based and data-driven) and different ways of processing the data using these different models. Such an increased diversity could be a way to devise a more performant ensemble approach for prognosis.

Appendix 1. Prognostic Metrics

In order to evaluate the average performance of RUL predictions, the common way is to apply several RUL predictions at different time steps to obtain a sequence of predicted RULs [START_REF] Saxena | Evaluating prognostics performance for algorithms incorporating uncertainty estimates[END_REF][START_REF] Saxena | Metrics for Offline Evaluation of Prognostic Performance[END_REF].

To evaluate the quality of prognostic outcomes, a synthesis of the prognostic metrics is used [START_REF] Baraldi | Ensemble neural network-based particle filtering for prognostics[END_REF][START_REF] Saxena | Designing data-driven battery prognostic approaches for variable loading profiles: Some lessons learned[END_REF][START_REF] Hu | Online Performance Assessment Method for a Model-Based Prognostic Approach[END_REF][START_REF] Rigamonti | Particle Filter-Based Prognostics for an Electrolytic Capacitor Working in Variable Operating Conditions[END_REF][START_REF] Al-Dahidi | A locally adaptive ensemble approach for data-driven prognostics of heterogeneous fleets[END_REF]. Figure 20 illustrates the RUL predictions with uncertainties at different prediction time steps t λ . The uncertainties are represented by the Probability Density Function (PDF). The accuracy bounds of a width of 2α shrinks with the prediction time index t λ , which creates the α-λ accuracy zone covering the true residual life RU L * . The upper bounds and the lower bounds of the α-λ accuracy zone:

RUL* + RUL ^- RUL ^CI + CI _ α + α - RUL EOL Time t λ
α + = RU L * t • (1 + α) α -= RU L * t • (1 -α) (21)

RU L

+ and RU L are the upper and lower bounds of the predicted RUL uncertainties, whereas CI + and CI -are the bounds of the confidence interval.

Based on those characteristics, different metrics are described as the follows:

• The accuracy index Acc t directly reflects the prediction errors relative to the true RUL:

Acc t = 1 - |RU L * t -RU L t | RU L * t ( 22
)
where RU L * t the true RUL and RU L t the median value of predicted RULs at prediction time t t . A larger value of Acc t indicates a better accuracy.

• The α-λ metric considers whether the predicted RU L lies within the ±α interval stating whether the required accuracy is met at a given time t λ . As being illustrated in Figure 20, the probability of lying within the α-λ accuracy zone is described by Equation ( 23):

αAc t = p α - t ≤ RU L t ≤ α + t (23)
where α + t and α - t are the upper and lower bounds of the accuracy zone. A higher value represents a better performance.

• The coverage index Cvg t considers whether the true RUL lies within the RUL prediction interval at time index λ for each trajectory:

Cvg t = p RU L CI- t ≤ RU L * t ≤ RU L CI+ t (24)
The value of Cvg close to 80% indicates a good representation of the uncertainty [START_REF] Baraldi | Ensemble neural network-based particle filtering for prognostics[END_REF].

• The precision index P rc t computes the relative width of the prediction interval, which is defined by:

P rc t = RU L CI+ t -RU L CI- t RU L * t (25)
where RU L

where L is the length of a sliding time window filtering the variances of the predicted EOL. Smaller values of Std t indicate better performance.

• The risk index Rsk t is the probability of obtaining an estimated RUL larger than the true RUL:

Rsk t = p( RU L t > RU L * t ) (27)
This index indicates the probability of receiving a later notification of a failure such that scheduling a maintenance after the failure is risky. Lower values correspond to a lower risk, which means a better performance.

Approach 1

The two indicators V st and γ are generated from the same realization of a Gamma degradation process, with different additive noises, which introduces some dependence between both indicators. Different parameters α i and β i are used for the i th stack. The objective of the simulation procedure is also to represent stack-to-stack variability around the average behavior given by: γ(t) = γ(t -1) + ∆γ(t) [START_REF] Gouriveau | IEEE PHM 2014 Data Challenge -Outline, Experiments, Scoring of results[END_REF] where ∆γ(t) follows a Gamma law G(α∆t, β). The average values of the Gamma process parameters α and β are preset according to our knowledge of PEMFC stack degradation in the following way.

The failure threshold F T γ (degradation rate) is set to 0.15, the average End of Life EOL is set to 1000 hours, and the slope of the degradation path is fixed to the computed value of α • β = F Tγ EOL , whereas the degradation variance is the value of α • β 2 . Thus, for the i th PEMFC stack:

γ i (t) = γ i (t -1) + ∆γ i (t) (29)
where ∆γ i (t) follows a Gamma law G(α i ∆t, β i ), and α i (resp. β i ) is drawn from a normal distribution around α (resp. β) with 5% variation. The measurements data are simulated according to Algorithm 3. Generate realization γ i (t) of a Gamma process with parameters (α i , β i )

5:

Add noises to γ i (t) to obtain the SOH degradation indexes for building the signals for Model 1 and Model 2:

γ i 1 (t) = γ i (t) + ε 1 where ε 1 ∼ N (0, σ 2 1 (t)) γ i 2 (t) = γ i (t) + ε 2 where ε 2 ∼ N (0, σ 2 2 (t)) 6:
Generate V i st index via Equations ( 2) and (3) using γ i 2 (t)

7:

Add measurement noises:

γ i meas (t) = γ i 1 (t) + ε meas,1 where ε meas,1 ∼ N (0, σ 2 meas,1 (t)) V st i meas (t) = V i st (t) + ε meas,2 where ε meas,2 ∼ N (0, σ 2 meas,2 (t)) 8: end for Approach 2
The two indicators V st and γ are simulated from two different degradation processes, dependent by construction. To this aim, a bivariate dependent Gamma process is constructed by trivariate reduction in the case of bivariate Gamma random vectors [START_REF] Devroye | Complexity questions in non-uniform random variate generation[END_REF].

Let us first recall that an univariate Gamma process [START_REF] Kahle | Gamma Processes[END_REF] with parameters (α, β) (where α, β > 0) is a subordinator such that for every t ≥ 0, the random variable G(t) is Gamma-distributed (αt, β) with probability density function : f (x; α, β) = β α x α-1 e -βx Γ(α) f or x, α, β > 0 [START_REF] Zhang | PHM-oriented Degradation Indicators for Batteries and Fuel Cells[END_REF] The principle of this second approach is to generate a bivariate Gamma process γ t consisting of two dependent deterioration process γ 1,t and γ 2,t . γ 1,t is then used directly as the SOH deterioration index measurements, and γ 2,t to generate the stack voltage measurement V st , see Algorithm 4.

Starting from three independent univariate Gamma processes g j t with (a j , b j ) for j = 1, 2, 3, one can build two dependent Gamma processes (or a bivariate dependent Gamma process) by trivariate reduction:

γ 1,t = g 1,t + g 3,t γ 2,t = g 2,t + g 3,t (31)
The process γ t = (G 1,t , G 2,t ) is, then, a bivariate subordinator [START_REF] Mercier | A preventive maintenance policy for a continuously monitored system with correlated wear indicators[END_REF] with Gamma marginal processes and marginal parameters (α j , β j ) where α j = a j + a 3 for j = 1, 2. The linear correlation between the two random variables G 1,t and G 2,t is independent of time t and described by the Pearson's correlation coefficient [START_REF] Mercier | A preventive maintenance policy for a continuously monitored system with correlated wear indicators[END_REF][START_REF] Pham | An imperfect replacement policy for a periodically tested system with two dependent wear indicators[END_REF]:

ρ = a 3 √ α 1 α 2 ( 32
)
where ρ is the Pearson's correlation coefficient, α 1 and α 2 are the marginal gamma parameters. Consequently, we have the following link between the two parametrizations (a 1 , a 2 , a 3 ) and (α 1 , α 2 , ρ):

a 1 = α 1 -ρ √ α 1 α 2 a 2 = α 2 -ρ √ α 1 α 2 a 3 = ρ √ α 1 α 2 (33)
where 0 ≤ ρ ≤ min(α 1 ,α 2 ) √ α 1 α 2 . This link allows to choose a 1 , a 2 and a 3 so as to generate a bivariate Gamma process with desired α 1 , α 2 and ρ. Within the range 0 ≤ ρ ≤ min(α 1 ,α 2 )

√ α 1 α 2 , trivariate reduction leads to one of the fastest algorithms known to date for bivariate Gamma distributions [START_REF] Devroye | Complexity questions in non-uniform random variate generation[END_REF]. As in Approach 1, the desired parameters α 1 , α 2 , β are determined to be consistent with the available knowledge of PEMFC stack degradation : the failure threshold F T γ (on the SOH level) and the average End of Life EOL.

Appendix 3. Data Size Analysis

This appendix proposed the results of a sensitivity analysis applied to decide for the data size for the experiment presented in this work. Two parameters need to be taken into account: 1) the number of training trajectories K, and 2) the number of selected nearest neighbor trajectories N . The following analysis are carried out under data generation with medium variance (α = 0.1, β = 1.5e-3) and dependency ρ = 0.5ρ max .

Figure 21 shows the accuracy surface of the ensemble with different training size K from 2 to 50, and neighbor size N from 1 to 15. The optimal values are not strict. We can obtain best accuracy outcomes with a larger number of training size and a larger number of nearest neighbors.

During the simulation, it is noticed that the number of the nearest neighbors N strongly affects the simulation time (changes exponentially). Thus, we would prefer smaller N while maintaining a good accuracy. Generate the realizations of Gamma processes with parameters Generate V i st index via Equations ( 2) and (3) using γ i 2 (t)

a 1 = α 1 -ρ √ α 1 α 2 , a 2 = α 2 -ρ √ α 1 α 2 , a 3 = ρ √ α 1 α 2 , b = β: g i 1 (
8:
Add measurement noises: γ i meas (t) = γ i 1 (t) + ε meas,1 where ε meas,1 ∼ N (0, σ 2 meas,1 (t)) V st i meas (t) = V i st (t) + ε meas,2 where ε meas,2 ∼ N (0, σ 2 meas,2 (t)) 9: end for Figure 22 shows the impact of training data size K on the prognostic accuracy. It implies that larger size gives better accuracy.

Figure 23 shows the number of nearest neighbors K with the accuracy. It can be noticed that, after more than around 15 nearest trajectories, the accuracy does not vary significantly. With K = 50, a better accuracy can be obtained with a smaller N around 5.

As been mentioned before, a smaller N is preferred and compared with it, the training size K does not affect the simulation time too much. Therefore, we set the number of training trajectories K = 50 and the number of nearest neighbor trajectories N = 5 for the experiment presented in this work.

3

 3

  Figure 1 shows the evolution of the stack voltage degradation. The decreasing trend represents the irreversible degradation, whereas the voltage jumps represent the reversible behavior caused by operating conditions modification.

AgingFigure 2 :

 2 Figure 2: Polarization curves during aging [28].

Figure 3 :

 3 Figure 3: Degradation estimation and RUL prediction. t λ is the prediction time.

Figure 5 :

 5 Figure 5: Selection of similar trajectories for RUL prediction assessment.

Figure 5

 5 Figure5illustrates the similar training trajectories retrieval for one test trajectory. For example, five nearest trajectories in the training set are selected regarding their ED to the test trajectory, which is computed during a time window of 100 hours. Note that at different prediction time, different trajectories are considered as the nearest ones. Therefore, the similar trajectories are selected locally.Step 2 : Local weights determination -The local weight assigned to each model of the ensemble is computed based on its local performance in terms of RUL prediction accuracy on the N selected training trajectories[START_REF] Baraldi | Local fusion of an ensemble of models for the reconstruction of faulty signals[END_REF]. Consider the average local error LE (m) t

Figure 6 :

 6 Figure 6: Simulated average SOH degradation γ and one realization of signal γ representing one stack.

Figure 7 :

 7 Figure 7: Loading current density j of one stack.

Figure 8 :

 8 Figure 8: Voltage state Vst of one stack.

Figure 9 :

 9 Figure 9: SOH degradation γ and voltage Vst measurements of one stack.

Figure 10 :

 10 Figure 10: Data simulation with three levels of variance: (1) Low variance; (2) Medium variance; (3) High variance.

Figure 11 :

 11 Figure 11: Local error evaluated over 50 training trajectories for test trajectory №40.

Figure 12 :

 12 Figure 12: Weight assigned to each model for test trajectory №40.

Figure 13 :

 13 Figure 13: RUL predictions aggregation for trajectory №40.

Figure 14 Figure 14 :

 1414 Figure14provides a global view of the average local error for all 50 test trajectories. Since each trajectory (stack) has different EOL, we normalized the time index considering the EOL ratioλ j = t λ EOL j. Globally, Model 1 prediction errors are lower at earlier life stages, whereas Model 2 errors gradually decrease thanks to the updating by sufficient incoming measurements and finally becomes lower.

Figure 15 shows

 15 Figure15shows the average prediction error for each test trajectory. It is the average local error of each trajectory along the entire prediction horizon (from 100 hours to 1000 hours): we can see that the ensemble gives the smallest prediction error for almost all the test trajectories.

Figure 15 :

 15 Figure 15: Mean absolute error for 50 test trajectories.

Figure 16 :

 16 Figure 16: Aggregated RUL predictions with uncertainty for trajectory №40 with accuracy and confidence interval (CI) bounds.

Figure 17 :

 17 Figure 17: Histogram of aggregated RUL uncertainties for trajectory №40 at different prediction time steps: (1) t λ = 200 hours; (2) t λ = 500 hours; (3) t λ = 800 hours.

Figure 18 :

 18 Figure 18: Point values aggregation: prognostic performance gain vs. dependence.

Figure 19 :

 19 Figure 19: Uncertainties aggregation: prognostic performance gain vs. dependence.

Figure 20 :

 20 Figure 20: RUL predictions at different prediction times t λ .

Algorithm 3 1 : 3 :

 313 Data simulation Approach 1 Choose F T γ , EOL, α, β 2: for i = 1 : number of simulated stacks Draw α i , β i , from normal distributions with average values α and β 4:

Algorithm 4 1 :√ α 1 α 2 4:

 412 Data simulation Approach 2 Choose F T γ , EOL 2: Determine α 1 , α 2 and β for marginal Gamma distributions 3: Given ρ, 0 ≤ ρ ≤ min(α 1 ,α 2 ) for i=1 : number of simulated stacks 5:

Figure 21 :

 21 Figure 21: Accuracy surface with different training size K and neighbor size N .

Figure 22 :

 22 Figure 22: Accuracy Acc vs number of training trajectories K (with N = 5).

Figure 23 :

 23 Figure 23: Accuracy Acc vs number of nearest neighbors N (with K = 50).

  from the prediction time 2: for i = 1, . . . , n // For each particle

	3:	Use model parameters estimated at time t λ (from Algorithm 1) : Θ i k , σ ω	i k
	4:

Table 1 :

 1 Prognostic performance metrics (Approach 1, medium variance)

	Average	Model 1 Model 2	Ensemble
	Performance			Point PDF
	Acc	0.52	0.12	0.55 0.56
	αAc	0.31	0.24	0.54 0.49
	Std	0.16	0.14	0.07 0.07
	Rsk	0.34	0.49	0.44	0.36
	P rc	0.25	0.27	-	0.32
	Cvg	0.48	0.37	-	0.74

Table 2 :

 2 Parameters used for the simulated examples (ρmax = 0.9128) Figures 18 and 19 show the performance improvements of the Ensemble with respect to the two individual models, considering different dependence scenarios, for point values and uncertainty, respectively. The improvements of performance metrics are computed in terms of percentage increased in the metrics' values, for example:

	ρ	α 1	α 2	β	a 1	a 2	a 3
	0	0.60 0.50 4.00 0.60 0.50	0
	10%ρ max 0.60 0.50 4.00 0.55 0.45 0.05
	25%ρ max 0.60 0.50 4.00 0.47 0.38 0.13
	50%ρ max 0.60 0.50 4.00 0.35 0.25 0.25
	75%ρ max 0.60 0.50 4.00 0.22 0.13 0.38
	90%ρ max 0.60 0.50 4.00 0.15 0.05 0.45
	ρ max	0.60 0.50 4.00 0.10	0 0.50
	6.2.2. Prognostic Performance Evaluation

  are the upper and lower bounds of the Confidence Interval (CI) of the predicted RULs distribution (e.g. CI = 50%) while RU L * t is the corresponding true RUL. Smaller values of P rc t indicate more precise predictions.• The steadiness index Std t measures the variance of the estimated value of the End of Life (EOL) when new measurements become available. It is defined as:

	CI+ t	CI-and RU L t
		Std t =	var( EOL (t-L):t ) EOL *

  t) : Gamma process realization with parameters (a 1 , b) g i 2 (t) : Gamma process realization with parameters (a 2 , b) g i 3 (t) : Gamma process realization with parameters (a 3 , b)

	6:	Generate dependent SOH indexes by trivariate reduction:
		γ i 1 (t) = g i 1 (t) + g i 3 (t) γ i 2 (t) = g i 2 (t) + g i 3 (t)

7:

The posterior PDF built using resampling in step (3) is used as the prior for the following iteration.

Here "Model 1" stands for "the prognostic approach based on Model 1", "Model 2" for "the prognostic approach based on Model 2", and "Ensemble" for "the prognostic approach based on the ensemble of models". This simplification is to avoid the wordy expression, and is used in the rest of the paper.

Acknowledgment

This work is partly supported by InnoEnergy and Labex ACTION (ANR-11-LABX-0001-01). We thank Dr. Sameer Al-Dahidi for sharing the pearls of wisdom with us during the course of this research.

Appendix 2. Introducing and controlling dependence in simulated deterioration trajectories

This appendix gives technical details on the simulation of dependent trajectories of the two deterioration indicators γ and V st . Two approaches have been considered in this work.