Joint particle filters prognostics for PEMFC power prediction at constant current solicitation Marine Jouin, Rafael Gouriveau, Member, IEEE, Daniel Hissel, Senior Member, IEEE, Marie-Cécile Péra, Member, IEEE, Noureddine Zerhouni Member, IEEE

Abstract-Proton Exchange Membrane Fuel Cells (PEMFC) are promising energy converters but still suffer from a too short life duration. Applying Prognostics and Health Management seems to be a great solution to overcome that issue. More precisely, developing prognostics to anticipate and try to avoid failures is a critical challenge. To tackle this problem, a hybrid prognostics is proposed. It aims at predicting the power aging of a PEMFC stack working at constant operating condition and constant current solicitation. The main difficulties to overcome are the lack of adapted modeling of the aging for prognostics and the occurrence of disturbances creating recovery phenomena through the aging. Consequently, this work propose a new empirical modeling for power aging that takes into account these recoveries based on different features extracted from the data. These models are used in a joint particle filters framework directly initialized by an automatic parameter estimate process. When sufficient data is available, the prognostics can give very accurate behavior predictions compared to experimentation. Remaining useful life estimates can be given with an error smaller than 5% for a horizon of 500 hours on a life duration of 1750 hours which is clearly long enough for decision making.

Index Terms-Prognostics, Particle Filter, Proton Exchange Membrane Fuel Cell (PEMFC), Aging, Remaining useful life. A power source failure is very often dramatic as all the systems relying on it become inoperative. Anticipating such a failure thanks to Prognostics and Health Management (PHM), and more particularly to its prognostics part, can be of great concern. Especially when the considered power source is a fuel cell and that the technology considered still suffers from a too short life duration. Fuel cells, and more precisely among all the existing technologies Proton Exchange Membrane Fuel Cells (PEMFC), are considered as promising electrical energy sources as they convert chemical energy coming from hydrogen and air into electricity, water and heat. Moreover they offer a great range of applications from transportation to micro-CHP (Combined Heat and Power generation) for building or power supply for portable devices [START_REF] Sharaf | An overview of fuel cell technology: Fundamentals and applications[END_REF]. Current PEMFC technology hardly meets the requirements for a large scale deployment: 2000 hours to 3000 hours of functioning can be achieved when 5000 hours are required for transportation applications and 100 000 hours for stationary ones. In [START_REF] Jouin | Prognostics and health management of PEMFC state of the art and remaining challenges[END_REF], the PHM solution is proposed to overcome the life duration issue of PEMFC. Various challenges are highlighted on all the layers of PHM, but the prognostics one is particularly empty. Indeed, although prognostics is a key process in PHM [START_REF] Gouriveau | Connexionist-systems-based long term prediction approaches for prognostics[END_REF], [START_REF] Ramasso | Remaining useful life estimation by classification of predictions based on a neuro-fuzzy system and theory of belief functions[END_REF], it just starts developing for PEMFCs. Prognostics approaches are classified in three categories: model-based, data-driven and hybrid [START_REF] Gouriveau | Connexionist-systems-based long term prediction approaches for prognostics[END_REF]. Among these approaches, two are present in prognostics of PEMFC applications. Four different types of works can be found in the literature: two datadriven approaches using echo-state networks [START_REF] Morando | Fuel cells prognostics using echo state network[END_REF] or adaptive neuro fuzzy inference systems [START_REF] Silva | Proton exchange membrane fuel cell degradation prediction based on adaptive neuro fuzzy inference systems[END_REF] and two hybrid prognostic approaches with Unscented Kalman Filter (UKF) in [START_REF] Zhang | An unscented kalman filter based approach for the health-monitoring and prognostics of a polymer electrolyte membrane fuel cell[END_REF] and with Particle Filter (PF) in [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF]. In [START_REF] Zhang | An unscented kalman filter based approach for the health-monitoring and prognostics of a polymer electrolyte membrane fuel cell[END_REF], the authors propose a prognostic model focusing on the degradation of the electrochemical surface area (ECSA) of a 1-cell PEMFC. The UKF-based framework is used to predict the decreasing of the ECSA and the Remaining Useful Life (RUL) of the cell. The results are very interesting when analyzed with the αmetric defined in [START_REF] Saxena | On applying the prognostic performance metrics[END_REF], even if the RUL estimates are all late predictions. However, this work suffers from some drawbacks. First, it is performed on a single cell with a very short period of functioning (≈ 300 hours). By looking at the lifetime targets listed earlier, the prediction horizon remains very short and in real case applications the PEMFC is never composed of a single cell. Moreover, although tracking the ECSA evolution is a good idea as the output power of the PEMFC is closely linked to its size, this surface is hard to measure without disturbing the operation of the PEMFC. Additionnaly, this kind of in-situ measurements tends to aggravate the aging of the PEMFC. So this prognostics could be hardly used for real PEMFC applications. On their side, the authors in [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF] use PF to predict the voltage aging of a 5-cell PEMFC stack responding to a constant current solicitation. Three empirical models are tested to set a basis for future prognostics works. The RUL predictions are promising and allowed to validate the choice of PF as a prognostic tool for such application. Nevertheless, the use of very simple models does not allow taking into account disturbances encountered during the aging. This work is pursued in [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF] in which a starting point for a solution to that limit is proposed. New models are introduced. The prognostics results start to be better, but the filters have to be initialized manually when the data or the horizon of prediction change. So although it gives interesting results, the models proposed are still incomplete and an automatic prognostics initialization has to be proposed. This paper pursues the work published in [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF]. The objective here is to set a prognostics able to take into account disturbances introduced by planned characterization while estimating with a minimum error the behavior of the PEMFC. The main contributions rely on a literature review of PEMFC aging, precise models for power degradation with time varying parameters, the introduction of an automatic parameter estimate process to initialize particles filter for prognostics. To support this, the paper is organized as follows. The first section is dedicated to PEMFC detailed presentation: the functioning principle is explained as well as the different components of the system. This allows introducing the ongoing issue of degradation understanding and modeling related to that particular system. Section III presents the modeling methodology adopted. The choice of using empirical models for prognostics is justified and the different steps for building an adapted set of models are explained. This set of models implies the computing of a parameter estimate process but also models' transformations to recursive forms. On this basis the association of different particle filters to perform prognostics is presented. This is the purpose of Section IV. It is directly followed by the application of the whole proposition in Section V. A comparison with our previous work [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF] is made to highlight the improvement of the prognostics estimates. The results are discussed to show the strengths as well as the weaknesses that should be corrected. Finally, some concluding remarks are given.

ACRONYMS

II. BACKGROUNDS: AGING OF PEMFC

Before presenting the new framework for prognostics of PEMFC, all the key information needed to understand what will follow regarding PEMFC stacks and their aging is presented in the section.

A. System overview

A PEMFC stack is a power device that converts chemical energy into electricity, water and heat (Figure 1). To do so, it uses the oxidation of hydrogen at a first electrode, the anode, to obtain electrons that flow into an external electrical circuit and protons that go through a proton exchange membrane. These protons and electrons meet again at the second electrode, the cathode, to reduce the oxygen [START_REF] Larminie | Fuel cell systems explained[END_REF]. The electrochemical reactions involved are:

2H 2 → 4H + + 4e - (1)
O 2 + 4H + + 4e -→ 2H 2 O (2)
(1) occurs at the anode while (2) occurs at the cathode. The general process can be given by:

2H 2 + O 2 → 2H 2 O + electricity + heat (3)
These reactions are not directly happening at the stack level but at the cell level. Indeed a stack is composed of a certain number of cells associated in series. The number of cells depends of the maximum power required from the stack. A cell is basically composed of:

• 2 bipolar plates to bring the reactant and collect electrons;

• 2 gas diffusion layers (GDL) to diffuse reactants toward the electrodes; • 2 electrodes: anode and cathode where occur the oxidation and reduction reactions; • 1 proton exchange membrane;

• sealing gaskets to ensure the impermeability of the cell. As the stack is only an energy converter it cannot work on its own and different ancillaries surround it for providing reactants, controlling the operating conditions or collecting the produced electricity. For our prognostics purpose, we consider that all the ancillaries work properly with no malfunctioning and no failure. That allows us to focus only on the stack. Operating conditions and current mission profiles are very important as far as the lifetime of a stack is concerned. Indeed, to respond in the best conditions to the current demand, operating conditions, namely temperatures, pressures, reactants flows or products evacuations, must be properly controlled. If not, it can lead to early failures preventing the stack to provide sufficient power. The value of the power is defined by the current setpoint imposed to the stack. In this study, the current profile used in experiments remains constant and the operating conditions are kept optimal. In that way, a power drop during aging will only be the result of inner degradations of the stack. This is a major hypothesis for all the following work. It is also important to mention that the experiments are conducted in a continuous manner. The stack was never remove from the test bench for storage, nor stopped and never experienced any current variations in an other occasion that the characterizations presented in the next paragraphs.

B. Reversibility and Irreversibility

During its lifetime, the stack looses gradually its ability to deliver its initial power. This is due to what are known in the fuel cell community as reversible and irreversible degradations. Reversible degradation refers to phenomena that lead to power losses but that can be partially or totally canceled by changing the operating conditions during a short period before returning back to nominal ones. Irreversible degradation refers to degradation as defined by the International standard [START_REF] Staff | European Standard EN 13306[END_REF]: "An irreversible process in one or more characteristics of an item with either time, use or an external cause". To avoid any ambiguity, the terms are redefined as follows:

• reversible degradation, that is in contradiction with the standard, is left aside in favor of "reversible" or "recoverable" phenomena, • irreversible degradation is now simply called degradation. The vocabulary set, some further explanations can be provided.

1) Reversible phenomena: These phenomena have been observed in a lot of works [START_REF] Cleghorn | A polymer electrolyte fuel cell life test: 3 years of continuous operation[END_REF], [START_REF] Kundu | Reversible and irreversible degradation in fuel cells during open circuit voltage durability testing[END_REF], [START_REF] Wu | Proton exchange membrane fuel cell degradation under close to opencircuit conditions: Part i: In situ diagnosis[END_REF], [START_REF] Prasanna | Effects of MEA fabrication method on durability of polymer electrolyte membrane fuel cells[END_REF] but are still not well understood. They can be caused by the interruption of the continuous testing for resting periods or for characterization with in-situ methods. Main reasons for reversibility of these phenomena can be water and thermal management which may change the water content and distribution in the different components in the cells [START_REF] Wang | Investigation of the recoverable degradation of PEM fuel cell operated under drive cycle and different humidities[END_REF].

2) Degradation: Degradation of a PEMFC stack is quite complex due to the different levels (stack-cell-components) and the numerous phenomena involved (chemical, electrochemical, mechanical, thermodynamics). Indeed degradations can come from:

• individual components in one cell (electrodes, membrane, etc.); • interfaces between components of a same cell;

• interactions between two consecutive cells;

• position of a cell within the stack (in the middle or near the edge of the stack);

Degradation type

• gradients (temperature, water) creating heterogeneities within the stack and affecting the aging; • varying operating conditions accelerating some specific degradation phenomena; All of these degradations have not been studied, and precise modeling of the ones that have been studied is very often not available. As writing a deep literature review of PEMFC degradation phenomena is not the point here, the reader may refers to [START_REF] Borup | Scientific aspects of polymer electrolyte fuel cell durability and degradation[END_REF], [START_REF] De Bruijn | Review: durability and degradation issues of pem fuel cell components[END_REF], [START_REF] Schmittinger | A review of the main parameters influencing long-term performance and durability of PEM fuel cells[END_REF], [21] for global reviews and to the references given in Table I for component-specific explanations. An important fact regarding degradation is that all the consequent dramatic failures that can prevent the stacks from providing an expected power have time constants in hours. As regards the failure definition, again it refers to the International standard [START_REF] Staff | European Standard EN 13306[END_REF]: "termination of the ability of an item to perform a required function". As examples, failures in a PEMFC can be the impossibility for the stack to deliver a minimum power or the impossibility of the membrane to perform its hydrogen-oxygen separator by letting the hydrogen crossing to the cathode side. Having time constants in hours is interesting for prognostics as it means that one data point per hour can be enough to use prognostic algorithms.

C. Existing aging models at the stack level

A complete aging model of a PEMFC stack should be able to take into account all previously listed: the aging of all the components within the stack, the interactions of the aging phenomena one with each others and the way they impact the efficiency of the stack, etc. Such a modeling is proposed in [START_REF] Franco | Modelling and analysis of degradation phenomena in pemfc[END_REF], but all the parts of the model are not explicitly given and reaching nano-scales does not seem appropriate for prognostics. Other models focusing on voltage or power modeling seem closer to prognostics expectations. In [START_REF] Lu | A semi-empirical voltage degradation model for a low-pressure proton exchange membrane fuel cell stack under bus city driving cycles[END_REF], the authors propose a semi-empirical modeling of the voltage degradation at low pressures to explain and model the data observed from a PEMFC stack used in a bus. The main idea is to model the different losses that can be observed on polarization curves (i.e. curves representing the voltage in function of the current) in which the model is not time dependent in its classic form and to introduce a time dependency. The same idea is also used in [START_REF] Xinfeng | Statistic analysis on voltage degradation of pemfc under road environment[END_REF] where a hundred stacks aged in small demonstration vehicles are modeled. But in both cases, most of the parameters of the models are empirical, so closely linked to their data. Moreover, in the applications considered, the stacks functioned with dynamic load profiles which is not our study case.

For prognostics, modeling of the power seems the most appropriate. Indeed, the standards for PEMFC stacks end of life (EoL) are all expressed in terms of power, as in [START_REF]The department of energy hydrogen and fuel cells program plan[END_REF] where the standard EoL is set at 10% of power loss during the lifetime. The context set, the modeling choice made for prognostics can be now explained in details.

III. MODELING OF PEMFC POWER DEGRADATION A. Physic-based vs empirical models

To use a particle-filtering-based prognostics, an aging model of the stack is needed. Two main criteria to define this model are its complexity and the possibility to access and measure its parameters on the stack. As described before, modeling all the aging process occurring within the stack, also the reversible phenomena that cannot be ignored, can lead to very complex models. In this case leaving aside physic-based modeling and using empirical models based on parameters monitored during the experiments offers a great alternative. Moreover, as we are working with constant operating conditions and constant current demand, the influence of parameters such as temperature, reactant pressures or current has not to appear explicitly in the model. The only mandatory parameter in this application is time.

B. Data

As this can be helpful to have visual examples to illustrate the modeling process, the data used for the experiments are described. Only one dataset is available, it comes from a commercial 5-cell PEMFC stack with an active area of 100cm 2 . The measurements performed on this stack are made thanks to the test bench and procedures described in [START_REF] Hissel | Characterization of polymer electrolyte fuel cell for embedded generators. test bench design and methodology[END_REF]. The experiment, realized on another FCLAB project, matches perfectly our hypotheses with a constant current profile and constant operating conditions and can be used for this prognostics purpose. For more clarity, the description of the dataset is given in Table II. This table provides the name that will be used in this paper for the dataset (D1), the mission profile in terms of current, the duration in hours of the experiments and also the number of characterizations performed on the stack during its lifetime.

For the experimentations, different measurements are performed:

• voltage measurements in continuous;

• current measurements in continuous (to compare with the input command); • punctual characterizations: polarization curves + electrochemical impedance spectroscopy. In this study, we propose to focus only on the power delivered by the stack. Consequently, as P = U.I only the voltage and current measurements are used. However, as their effects cannot be ignored, characterizations have to be mentioned. Power measurements can be seen on Figure 2. Let's now focus on how to model the power during the aging of the stack.

C. Power global modeling 1) Formalisation:

The main idea is to create a set of models that is able to describe accurately the behavior of the stack during its aging. Although it will lead to purely empirical models, it can be interesting to chose trends or descriptors that could be linked to real phenomena occurring within the fuel cell stack. First, the data set can be divided into two types of events (see Figure 2):

1) parts of continuous aging 2) interruptions due to characterization processes As the operating conditions and current are varying during the characterizations, these phases will not be modeled. However, assumptions will be made a little further regarding the consequences induced. Let's focus on the different parts of continuous aging. D1 counts thirteen parts of continuous aging of different durations. Based on existing works [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF], [START_REF] Hou | An improved dynamic voltage model of PEM fuel cell stack[END_REF] and different model fittings, it appears that each part follows the equation:

P (t) = -a.ln(t) -b.t + c ( 4
)
with t representing time but with different values of coefficients a, b and c. For comparison purpose, t is set to zero in the beginning of each segment. This empirical modeling is now justified by actual phenomena occurring within the stack.

2) Justification: It can be observed that as time evolves, the aging effects occur faster and faster leading to more serious power drop. Moreover the impact of the logarithmic part of the model seems to accentuate and last longer. It is illustrated in the bottom part of Figure 2, which shows the normalized powers of all segments to compare them visually. It is resulted from dividing the power by its initial value in each segment. This can be explained if we give a physical meaning to the different parts of the model. Due to the diffusion phenomena of gas and water and also to reversible phenomena, the stack goes through a transient phase after the operating conditions are set to their nominal values. It is represented by the logarithmic part of the model. Once that transient phase is over, the stack enters a steady power decay represented by the linear part. This distinction between the transient and steady parts is further explained in [START_REF] Kundu | Reversible and irreversible degradation in fuel cells during open circuit voltage durability testing[END_REF]. As the stack ages, its components might have some difficulties to fulfill their functions leading, for example, to slower gas diffusion through the gas diffusion layers, to presence of contaminants at the membrane and electrodes or to greater water accumulation at the cathode by loosing hydrophobic properties. Combined together, these phenomena, coming from both degradation and reversible events, may lead to longer transient phases with a greater power drop. Degradation, namely of conductivity and activity of the membrane and the electrodes, also lead to more severe power drop during the steady state, reflected by greater coefficients in the linear part of the model. It is now obvious that our model coefficients, a and b (c left apart because not used afterward), will evolve with time. It can be interesting to figure out if they will follow particular predictable trends.

D. Parameter trend modeling

To extract trends for a and b from (4), these coefficients are determined by fitting the equation to the data on the different aging segments. The whole data set is used to that goal. Equation ( 4) is fitted successively to all the aging segments and the different values of a and b are determined. This is achieved by using a least squares fitting. These values, which are constant in each segment, are shown on Figure 3 by the blue dots. Both coefficients seem to follow trends that can be modeled easily with exponential functions (Figure 3, red curves). Indeed the evolution of the transient phase (coefficient a) can be easily modeled by:

a(t) = a 1 .exp(a 2 .t) + a 3 .exp(a 4 .t) (5)
whereas the acceleration of the linear part is harder to model precisely:

b(t) = b 1 .exp(b 2 .t) + b 3 (6)
However, although this last seems less obvious when the global trend is compared to the data, it might not be so far from the reality. Indeed, the authors in [START_REF] Zhang | The characteristics of voltage degradation of a proton exchange membrane fuel cell under a road operating environment[END_REF] showed that the degradation rate for constant current follows a bathtub-like curve. Here, if we consider that the linear part of the model ( 4) can be assimilated to a degradation rate in a steady state and that the stack has already left its early life (manufacturers usually let the stack working around 100 hours to leave the early life stage before delivering the stack to the customers), the equation ( 6) can find a justification. Now that we are able to model the different continuous aging part, we have to wonder what the initial power at the starting point of each part is .

E. Recovery modeling

After each characterization phase, the stack seems to recover some power (see Section II-B). However, it is not clearly explained. It may be linked to change in gas and liquids' repartitions within the stack when the operating conditions return to their nominal values after variations during the characterizations. It can be seen on Figure 4 that this recovery does not remain constant during the aging. In terms of difference of power between and after the characterization phase, it increases. But by looking at the values reached after all the characterizations, the maximum recovery power decreases with time. Indeed, as the stack ages, the recovery is limited by the degradation of the stack components. For our prognostics purpose, we choose to follow the evolution of the maximum recovery power through time. As for the previous parameters, the recovery is extracted from the data and a global trend of it as aging is found. It follows an exponential form:

Rec(t) = r 1 .exp(r 2 .t) + r 3 .exp(r 4 .t) (7)
This exponential trend is not a surprise, it is coherent with the fact that the main degradations impacting the power drop follow exponential trends. Among these degradations, the loss of active surface area [START_REF] Zhang | An unscented kalman filter based approach for the health-monitoring and prognostics of a polymer electrolyte membrane fuel cell[END_REF] and the increase of hydrogen crossover from the anode to the cathode [START_REF] Schmittinger | A review of the main parameters influencing long-term performance and durability of PEM fuel cells[END_REF], [START_REF] Baik | Effects of operating parameters on hydrogen crossover rate through nafion membranes in polymer electrolyte membrane fuel cells[END_REF] can be In this part, we have selected different parameters and built models that will help follow the aging of the stack. To stick to the reality, all the parameters chosen can have physical interpretations. This has now to be integrated in a prognostics framework to help predicting the remaining useful life of PEMFC stacks.

IV. PROGNOSTICS OF PEMFC

To achieve power behavior prediction and RUL estimation, our prognostics framework is divided in two main part: 1) a parameter extraction part; 2) a particle filter-based prognostics. This is illustrated by Figure 5. The main reason why parameter extraction and prognostics are computed together is that the parameter estimate part directly creates the initial distributions needed by the particle filters. When the length of the learning set varies, the parameter estimate results change giving different initializations to the prognostic part. As more data become available, the parameter estimate part gives better results and this directly impacts the prognostics predictions. Before detailing the whole framework, a short reminder on particle filters theory is given here.

A. Particle filters

Particle filters are great tools to solve nonlinear Bayesian tracking problems. Indeed, the models defined above meets completely the hypotheses of this kind of problems as they are nonlinear, non-exact (with unknown coefficients), nonstationary and they might contain non Gaussian noise. A Bayesian tracking problem is defined by two equations [START_REF] Doucet | Sequential Monte Carlo Methods in Practice, ser. Information Science and Statistics[END_REF], [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF]:

• the state equation

x k = f (x k-1 , ϑ k , ν k ) (8)
• and the observation model

z k = h(x k , µ k ) (9)
In the state equation, {x k , k ∈ N} represents the state that has to be followed and predicted, ϑ k parameters of the model and ν k an independent identically distributed (i.i.d.) noise.

These variables change at each step and the passage from step k -1 to k is made thanks to the transition function f . Regarding the observation model, {z k , k ∈ N} represents the measurements, µ k an i.i.d. noise and h the observation function. As measurements are directly available from the data, no observation model has to be built in this work.

To obtain a distribution of the possible states of x at time k, the probability density function p(x k |z 1:k ) has to be built. The starting point is given by the initial state distribution p(x 0 |z 0 ) ≡ p(x 0 ). It is assumed that this initial pdf is known. The optimal Bayesian solution is obtained by repeating a prediction and an update stages: 1) prediction

p(x k |z 1:k-1 ) = p(x k |x k-1 )p(x k-1 |z k-1 )dx k-1 (10) 2) update p(x k |z 1:k ) = p(z k |x k )p(x k |z 1:k ) p(z k |z 1:k-1 ) (11)
However, this optimal solution cannot be obtained analytically in most of the cases. Here comes the particle filter to give an approximate solution.

Particle filters are from the family of Monte Carlo-based tools using the Baye's theorem as basis. The first stage consists in splitting the initial state distribution p(x 0 ) into n samples named particles. Three steps are then repeated until the desired results is obtain (Figure 6): 1) Prediction: particles are propagated thanks to the state model from step k -1 to k giving a new pdf. 2) Update: the measurement z k is used to calculate the likelihood function p(z k |x k ) and give weights to the particles. Particles representing states close the last measurement have the higher likelihood and represent the most probable states. 3) Re-sampling: particles with the lower weights are eliminated and particles with the higher weights are manifolded. This procedure allows avoiding a degeneracy of the filter, namely the increasing number of particles with low weights that would lead to bad results in the previous stages. This process is applied all along the learning of the prognostics. In fact, as soon as no data become available, no measurements z k can be used to calculate the likelihood function and the prognostics enters its prediction phase. In that step, the particles are simply propagated by the state model. Once the failure threshold is reached, the final distribution of the particles give the most probable state, i.e., the state represented by the weighted average of the particles, and its uncertainty distribution, i.e., the distribution of the particles.

B. Parameter extraction 1) Parameter extraction on raw data:

The raw data at the entry of the parameter extraction box is composed of the power signal and the characterization calendar. Thanks to the calendar, the different parts of continuous aging available for the learning are identified. The fitting procedure uses the data to identify the coefficients a and b from equation (4). A robust least squares algorithm is used for that purpose. The recovery values after each characterization are also identified at that stage. These values are stored in tables.

2) Model identification and filter initialization: In the previous section, we have defined four models and three of them ( 5) to [START_REF] Zhang | An unscented kalman filter based approach for the health-monitoring and prognostics of a polymer electrolyte membrane fuel cell[END_REF] have unknown coefficients that should be identified. Thanks to the data stored in the previous step the parameters a 1 to a 4 , b 1 , b 2 and r 1 to r 4 are estimated. As their identification relies on values estimated on the previous step, these estimates may come with uncertainties. As our procedure does not estimate this uncertainty, we decide to associate uncertainty ranging from ±1% to ±5% depending on the order of magnitude of the parameter. Consequently, our model identification directly gives us distributions for intializing the particle filters used just after. These distributions are uniform distributions centered on the value obtained by identification with a width given by the associated uncertainty. This step imposes a first limitation to the prognostics framework. Indeed, as equations ( 5) and ( 7) have each one four coefficients to identify, it implies that four characterizations have already been performed in the data used for the learning. Now that the extraction of parameters is done, the information extracted is used to perform prognostics.

C. Joint particle-filtering-based prognostics

For the record, the final objective of that work is to predict the behavior of the power delivered by the stack through time and to estimate the RUL at various points in time. To achieve this goal, we have a global power model ( 4) and three other models dedicated to specific parameters from model ( 4). As all the future evolutions of these models have to be predicted, a framework with four particle filters is set. Of course, the framework could have been reduced to two by including equations ( 5) and ( 6) in the first one. But, in order to evaluate the quality of the different models, keeping them separate is more interesting for further investigation. To integrate them in the particle filters as state equations, all the equations have to be written in a recursive manner. However, one important assumption is made regarding the global power model: coefficients a and b remain constant on the full length of a continuous aging part. This assumption is made to remain consistent with the way these parameters were extracted in the first place. The continuous aging parts are numbered and this is represented by the index i in the next equation. Consequently, although the parameters a and b are estimated continuously, only the first value estimated at the beginning of the aging part i is stored as a i and b i . Then when a characterization happens and a new aging part starts, the last value of a and b obtained at time (k + length of segment i) is used to define a i+1 and b i+1 . So, the characteristics of the different filters are as follows:

• Filter 1: Global power -State equation:

P k = a i .log(k -1/k) -b i + P k-1 (12)
-Measurements: power obtained from voltage and current measurements

• Filter 2: Coefficient of the logarithmic part a -State equation:

a k = a 1 .exp(a 2 .k).(1 -exp(-a 2 )) + a 3 .exp(a 4 .k).(1 -exp(-a 4 )) + a k-1 (
Rec k = r 1 .exp(r 2 .k).(1 -exp(-r 2 )) + r 3 .exp(r 4 .k).(1 -exp(-r 4 )) + Rec k-1 (15)
-Measurements: recovery extracted during the parameter extraction It can be noticed that no noise has been introduced in the different state equations. As general trends and not precise models are expected, we consider that this noise can be ignored for prognostics. As far as the measurement noises are concerned, as the raw data are used, the noise is already included in them. So for each filter, the z k used are the noisy measurement. However, the noise contained in z k should be taken into account when comparing the filter estimation to the actual data. As the characteristics of the noise are unknown in that application, it is assumed that it is a Gaussian noise with an unknown standard deviation σ. It allows to use the method proposed in [START_REF] An | A tutorial for model-based prognostics algorithms based on matlab code[END_REF]. It states that the standard deviation of the measurement noise can be considered as a variable to be estimated and that it becomes part of the state vector. Consequently, this standard deviation is estimated and propagated thanks to the particles. It then intervenes when the likelihood of each particle is calculated which becomes a function of the state, the model parameters and the noise standard deviation estimates: L(z k |x i k , param i k , σ i k ). The four filters are working in parallel, synchronized on the same time step. During the period of continuous aging, they are all predicting their attributed state model with no interaction with each other. When the date of a characterization is met, the particles distributions for parameters a and b in Filter 1 are updated thanks to the last distributions predicted by Filter 2 and Filter 3. Moreover, the last distribution of the recovery Rec k is used to replace the state particles of P k . This process is represented on Figure 7 for more clarity.

D. Uncertainty propagation

The uncertainty coming all along the prognostic process is one of the major issues of prognostics applications. This uncertainty may come from the data, the parameter estimate and the prognostic part. Consequently, all the framework was built to take into account the uncertainty and to propagate it from one step to another. Regarding the extraction part, the treatment of uncertainty has already been explained: distributions centered on the parameters identified are drawn. For the prognostics part, the uncertainty propagation is ensured by Filter 4 before the last characterization and after by Filter 1. Indeed, as particle filters are working with particle distributions, they naturally give the uncertainty coming with the prediction (see Section IV-A). As the power given by Filter 1 is regularly updated by external

Continuous aging part i

Continuous aging part

2 -a(i) b(i) x ln(k-1/k) - + P(k-1) P(k) = Filter 1 Filter 3 Rec(k) Filter 4 b(k) a(k) Filter 2 -a(i+1) b(i+1) x ln(k'/k'+1) - + P(k') P(k'+1) = Filter 1 a(k'+1) Characterization Time step k Time step k' = k+ duration charact. + 1 P(k') Filter 1 Rec(k') b(i+1) a(i+1) b(k') a(k') Time step k' + 1 …… Fig. 7.

Filters interactions between and just after characterizations

Distributions coming from R estimations propagate the uncertainty until the last characterization P propagates the uncertainty that will be kept at the end Fig. 8. Models role in uncertainty propagation parameters, the uncertainty coming with power behavior predictions changes after each characterization. When the power is concerned, only the recovery prediction given by Filter 4 propagates the same uncertainty from the beginning to the last characterization. After this last one, the last distribution given by Rec is used by Filter 1 to predict the power until the failure threshold (Figure 8). By proceeding that way, the RUL distribution, which is the final goal of our prognostics framework, takes into account all the possible sources of uncertainty of the framework. To demonstrate the ability of this framework to give good predictions, series of tests are performed on the dataset presented earlier.

V. EXPERIMENTS AND DISCUSSION

A. Experiment settings

The main parameters to define for these experiments are the EoL threshold and the length of the learning. Regarding the threshold, by referring to what is explained earlier, it should be set at 10% of the initial power value. Nevertheless, if this threshold is respected, only a small part of the data would be used: power would be around 184 W corresponding to less than 1300 hours. To benefit from all the data available, the EoL threshold is chosen to be located at 16.3% of the initial power. This value is calculated thanks to the difference between the value of the power at t=0 and the value at the end of the experiments.

Then, the learning length has to be defined. In order to estimate the largest prediction horizon that can be reached with a reasonable prediction error, different length of learning are tested from 500 hours to 1700 hours with a step of 100 hours. The minimum value of 500 hours cannot be reduced as at least four points are needed to identify coefficients in equations ( 5) and [START_REF] Zhang | An unscented kalman filter based approach for the health-monitoring and prognostics of a polymer electrolyte membrane fuel cell[END_REF].

It is also important to mention that the four particle filters are SIR (Sampling Importance Resampling) type [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF]. The filters are programmed in Matlab language and inspired from [START_REF] An | A tutorial for model-based prognostics algorithms based on matlab code[END_REF]. To avoid any degeneracy problem during the testing, they are all working with 10 000 particles. This can take a long computing time for obtaining results, but this duration is not yet a performance criterion in this work. Before showing the results, it should be mentioned that in Figures 9, 11 and 15, vertical lines appear at the time of characterizations as already observed in Figure 2. As the current varies during the characterization breaking the hypothesis of constant current solicitation for prognostics, no prediction is made during that period. For convenience of use, the state vector is filled with zeros at the dates of these events creating the vertical lines.

B. Power aging predictions

1) General comments: For each learning the prognostics returns the curves shown on Figure 9. The left upper part represents the prediction of the coefficient a by Filter 2, the right upper part shows b predicted by Filter 3, the left lower part shows the recovery from Filter 4 and finally the last part presents the power predictions by Filter 1. For illustration purpose, the confidence interval given by the filters appears on the power and b graphs. A zoom of the confidence interval coming with power aging can be found on Figure 10. The results obtained for power aging predictions can be divided in two parts according to the duration of the learning: bad predictions with a learning between 500 and 1100 hours and good predictions for 1200 hours and more. Indeed, for small learnings, the models seems to fail catching the behavior of the stack (Figure 11). It indicates that one or more parameters of the models cannot be predicted accurately.

2) Parameter estimate vs learning duration: The results given by the prognostics are strongly determined by the parameter estimate part. Indeed, when very few points are available identifying parameters a, b and Rec trends by fitting their models can lead to false approximations. As more data become available, the trends start to draw closer to the reality. This can clearly be seen on Figure 12. The recovery Rec seems to be the major limitation when less than six points are available, corresponding to a learning of For 700 hours, another part of the explanation concerns the coefficient a driving the logarithmic part of the model. For short learnings, this coefficient is not well estimated. Some of the extreme data points are gathered between 500 and 700 hours preventing the parameter estimate from approximating the good trend. However, after 1000 hours the feature extracted starts to converge toward the real trend.

Regarding the coefficient b which drives the linear part of the model, the parameters extracted are more dispersed and most of the time overestimates the trend. The impact of this dispersion is noteworthy on the results for short training, particularly for 900 and 1000 hours where a negative sign appears. However, it seems that the overestimation of this coefficient is more or less counterbalanced in the power aging prediction by the error on a approximation.

3) Evaluation of the error: To further evaluate the predictions, it can be interesting to take a closer look a the RMSE and the MAPE observed during learning and prediction for the different learning duration (Figure 13 and Table III). For the learning the MAPE is bounded between 0.4% and 0.9% with a mean of 0.574%. As regards predictions, the MAPE seems to decrease exponentially with the increasing of the learning duration. This is more visible by using the RMSE values available in Table III. The maximum MAPE value is 5.06% obtained for a learning of 600 hours to reach its minimum of 0.33% for 1700 hours. The decreasing is logical, as when more data become available and the EoL is closer, the predictions become more accurate. These values also reflect what happens on Figure 11 and the ability of the models to predict the behavior as discussed above.

C. RUL estimates

Calculating and representing the remaining useful life evolution through time can give us an idea of the prognostics horizon that can be expected from the framework. Figure 14 shows the RUL estimates with the lower and upper bounds of their distribution (red lines) and also the relative error between actual RUL and estimates (blue lines). Only the results for a learning equal and greater than 1100 hours are shown on this figure, as all the previous estimates exhibit errors greater than 600 hours. The reasons for that have already been discussed earlier.

RUL estimates are very promising. Indeed, for a learning greater than 1300 hours, the maximum relative error is 57 hours for 1400 hours. It means that for a prediction horizon of 450 hours, the RUL can be given with an error smaller than 5%. To define more precisely what could be the larger horizon with a maximum error of 5%, new tests are launched. The step between two prediction is reduced from 100 hours to 10 hours. The results are summarized in Table IV. Although the error is twice superior to 5% after that date, the maximum horizon seems to be located around 518 hours for a learning. Assuming that the stack is working continuously, it gives 21 days to react and to take the right mitigation actions which is fairly correct.

It can be seen thanks to RUL bounds on Figure 14 and Table V that RUL distributions are very small. Indeed, the larger RUL distribution is obtained for a learning of 1500 hours and it is 21 hours long but the mean value is 9.3 hours. This demonstrates a very small uncertainty coming with RUL estimates. In a real case application, it means that the date when the stack is going to reach its EoL could be known to the precision of a single day. Considering that almost all the prediction are early predictions, the prognostics will be very interesting for deciding mitigation strategies or system maintenance.

D. Comparison with previous work

As a comparison between models that take into account disturbances (in [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF]) and those that don't (in [START_REF] Jouin | Prognostics of PEM fuel cell in a particle filtering framework[END_REF]), was already done in [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF], showing the improvement in estimating the state and reducing the prognostics uncertainty, this comparison only focuses on comparing the results presented here and those presented in [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF]. The EoL threshold was exactly the same as in this paper. As stated in the introduction, the main differences rely in the modeling and in the automatic initialization of the filters. In [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF], no modeling for a existed and the only variations of the parameter was due to the estimation of the particle filter in charge of power prediction. Without the model, the linear part was dominating leading to poor behavior estimations on the parts where the logarithmic effect was accentuating. The RUL estimates seemed to be globally better mainly thanks to 2 factors:

1) the behavior prediction luckily meets the EoL threshold at the right place, if another threshold was defined the RUL estimates would not be so good; 2) the particle filters were manually initialized, different tries were required to define the initial distributions that would lead to acceptable results: in practical applications such a practice would not be usable. So, the RUL estimates are left apart and only the behavior prediction is considered. Figure 15 shows both the old and new models behavior estimations when using all the data available for training. The new model catches better the behavior of the power during the aging of the stack. The adding of the logarithmic part modeling improves clearly the results. Moreover, it shows that the automatic initialization of the filters can lead to results as good as for manual initialization, with the exception of the learning duration limitation.

E. Discussion

The analysis of the results shows some good points as well as some limitations to the prognostics developed here. First, the occurrence of the characterizations is defined by a calendar. This calendar is known to perform prognostics. The framework should also be able to deal with unplanned characterization. Then, modeling independently the different phenomena related to the occurrence of disturbances seems to be a good solution. Indeed, it allows modeling precisely the aging of the different parameters with time. The empirical expressions defined for each of them are simple but sufficient to obtain behavior Comparison between the old model from [START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF] and the new model predictions with a small MAPE. But using different models with unknown coefficient also creates one major limitation of that work: the amount of data needed to have good trends from the parameter estimate based on the least squares method covers more than 1100 hours which is around 60% of the stack lifetime in this case study. Using a different parameter extraction method could be a solution to that problem. Trying to find a deeper link with physics to have fix values of some coefficients might be another one.

VI. CONCLUSION

This work presents a new modeling for prognostics of PEMFC under constant current solicitation but facing events creating disturbance in the natural aging behavior. A global model for power aging is set including time-variant parameters. These parameters can be modeled independently as they show specific trend as the stack degrades. Moreover, they allow highlighting different physical phenomena as transient effects, steady state aging and recoveries. These models are parts of a global prognostics framework composed by a parameter extraction part and a particle filtering part. The automatic parameter estimate based on least square fitting allows identifying the unknown coefficients of all the models and thereby initializing four particle filters. These lasts interact together to allow predicting the power aging of the stack and estimating the remaining useful life. The whole proposition is quite convincing as it can predict the RUL with less than 5% of error with a horizon of more than 500 hours. However, a major limitation relies in the length of the data needed for the learning as more than 1100 hours are needed, covering more than the half of the stack life. Consequently, the parameter estimate part should be improved to use less data while the prognostics keeps offering good results. To ensure that this prognostics is sufficiently generic, new experiments are conducted to provide new aging data coming from different fuel cell stacks. Once this method is definitely validated, the next step of this work will be to include the current in the modeling and so being able to manage variable current profiles. This will lead to prognostics closer to real automotive or stationary applications.

Fig. 1 .

 1 Fig. 1. Stack, cell and functioning

Fig. 2 .

 2 Fig. 2. Top: Dataset D1 -Bottom: Normalized voltages of the aging parts

Fig. 3 .

 3 Fig. 3. Trends extracted for coefficients a and b

Fig. 4 .Fig. 5 .

 45 Fig. 4. Trend extracted for the recovery

Fig. 6 .

 6 Fig. 6. Particle filter principle

13 )-• Filter 3 :-•

 133 Measurements: values of a extracted during the parameter extraction Coefficient of the linear part b State equation: b k = b 1 .exp(b 2 .k).(1 -exp(-b 2 )) + b k-1 (14) -Measurements: values of b extracted during the parameter extraction Filter 4: Power recovery -State equation:

Filter

Fig. 9 .

 9 Fig. 9. Example of results for a learning of 1300 hours

Fig. 10 .

 10 Fig.10. Zoom on the confidence interval for power prediction at 1300 hours

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Parameter estimates with different learnings

Fig. 14 .

 14 Fig. 14. RUL evolution through time

  Fig. 15. Comparison between the old model from[START_REF] Jouin | Prognostics of proton exchange membrane fuel cell stack in a particle filtering framework including characterization disturbances and voltage recovery[END_REF] and the new model

TABLE I PROPOSED

 I LECTURES FOR DEGRADATION OF COMPONENTS

TABLE II CHARACTERISTICS

 II OF THE DATASETS

	Name	Current profile	Duration	Nb charact.
	D1		1750 hours	12

ACKNOWLEDGMENT

The authors would like to thank the ANR project PROP-ICE (ANR-12-PRGE-0001) and the Labex ACTION project (contract "ANR-11-LABX-01-01") both funded by the French National Research Agency for their support.

Systems. In 2008, he joined the FEMTO-ST (CNRS) Institute and became Head of the "Energy systems modelling" research team. Since 2012, he is the Head of the "Hybrid & Fuel Cell Systems" research team in the same Institute and also currently the Director of the FCLAB (Fuel Cell Lab) Research Federation (CNRS), devoted to Fuel Cell Systems Research and Technology. His main research activities are concerning fuel cell systems dedicated to automotive and stationary applications, modelling, non linear control and energy optimization of these systems and fuel cell system diagnosis. He was Associate Editor of IEEE Transactions on Industrial Electronics between 2004 and 2012 and is Associate Editor of ASME Fuel Cell Science and Technology. He is also the Chair of the IEEE VTS French Chapter and member of the advisory board of the MEGEVH network, the French national network on EV and HEV. He has published more than 290 scientific papers in peer-reviewed international journals and/or international conferences.

Prof. Marie-Cécile Péra was born in 1968 in Paris, France. She received a PhD in electrical engineering from the Institut National Polytechnique de Grenoble, in 1993. From 1994 to 1999, she was associate professor at the University of Reims Champagne Ardennes, where she studied non-linear dynamics of electrical systems, based on chaos theory. Since 1999, she has joined the University of Franche Comté (UFC). In September 2008, she became a full Professor. She's the deputy director of the FEMTO-ST Institute. She's a researcher of the Energy Department of FEMTO-ST and a member of the FCLAB Research Federation. She develops models based on the Energetic Macroscopic Representation (EMR) of power generation system (fuel cells, PEMFC and SOFC, supercapacities, batteries), their control and the energy management. She works on diagnosis and prognosis of fuel cells. She has contributed to more than 200 publications in international scientific journals and international conferences.

Prof. Noureddine Zerhouni received his engineering degree from National Engineers and Technicians School of Algiers (ENITA) in 1985. He received his Ph.D. Degree in Automatic Control from the Grenoble National Polytechnic Institute in 1991. In September 1991, he joined the National Engineering School of Belfort (ENIB) as Associate Professor. Since September 1999, Noureddine Zerhouni has been a Professor at the national high school of mechanics and microtechniques of Besanc ¸on. As for investigation, he is with AS2M department of FEMTO-ST Institute. His main research activities are concerned with intelligent maintenance systems, and PHM.


João Marcelo

Castro Monteiro

Jean-Pierre Daniel

A PROPOSAL FOR COLLABORATIVE NAVIGATION IN CORRIDORS

The current study presents the proposal of a collaborative method to take advantage of the available information from surrounding airplanes in a corridor and provide a back-up navigation aid with a performance equivalent to GPS. Here it is considered that an aircraft may use the information of distance between itself and other airplanes provided by TCAS as well as the absolute position information of aircraft made available by ADS-B, creating a dynamic constellation similar to GPS where the satellites are replaced by the surrounding airplanes. A simulation study considering two types of traffic scenarios is developed to access the resulting accuracy of the proposed system.

Introduction

With the expected huge increase in air traffic over the next decades, it is believed that a large proportion of traffic will be organized in main air streams operating under advanced CNS procedures. Many studies and proposals of solution for a future congested ATM are developed and presented in many scientific events and magazines, and one of the most considered is known as "Flow Corridors". In order to enable this type of traffic flow, with minimal or no interference at all of controllers, the airplanes must follow certain traffic specifications of the corridor and be able to maintain self separation and surveillance of the inside traffic. In such environment, self-controlled dense traffic, the navigation should be based on accurate and reliable systems throughout the entire corridor.

The current study aims to present a collaborative method to take advantage of the many measurements from surrounding airplanes to improve their navigation performance. For the purpose of this study, it is understood that collaborative navigation is characterized by the utilization of measurements performed by surrounding aircraft and exchanged to improve the navigation performance of the each of them. Said in another way, an airplane surrounded by other airplanes should experiment a better navigation performance than an aircraft flying alone.

The main objective of this study is to assess the possibility of using TCAS and ADS-B signals, already made available in commercial aircraft, as means to improve navigation performance. To improve navigation performance, an aircraft may use the information of distance between itself and other airplanes provided by TCAS, since the estimation of the absolute position of each airplane is available (ADS-B), and to create a dynamic constellation similar to GPS where the satellites are replaced by the surrounding airplanes. In this situation, instead of having the ephemeris with the satellite position we have an ADS-B message with the aircraft position and instead of measuring the distance aircraftsatellite using the satellite clock and received message time, we use the TCAS to measure the distance between aircraft. The challenges of this approach are not only the level of uncertainty of each measurement system, but also the lack of synchronization between measurements.

The mathematical model used in this study is similar to the one used for GNSS positioning applications, with a slight difference adopted to fit the proposed method. The concept used by GNSS is the Time of Arrival (TOA) where the time that a signal sent by the transmitter (satellite) to reach the user receiver is taken into account and its multiplication by the speed of light gives the distance between them. This supposes that the information about distance satellite-user and satellite position are synchronized. However in the considered case, all surrounding airplanes send their ADS-B message at different times and even the distance measurement and position signals of a given aircraft are not synchronized. So, to apply a least square technique to provide an additional estimation of position for each aircraft, it is necessary first to synchronize all this information with respect to a unique time line, a reference aircraft time line.

Collaborative Navigation in Corridors

Here the main ATM concepts used in this study are introduced.

Corridors

Many studies and proposals of solution for a future congested ATM are developed and presented in many scientific events and magazines, and one of the most considered is known as "Flow Corridors". Corridors are like long tubes with parallel trajectories in the vertical and horizontal planes, where the aircraft could fly for long distances with minimal interference from the conventional surrounding traffic as Youse et al. describe in [START_REF] Youse | Nextgen Flow Corridors Initial Design, Pro-cedures and Display Functionalities[END_REF]. Figure 1 shows a simplified corridor view with two flight levels and two lanes for each level (2 x 2).

Figure 1. Basic flow corridor

Several characteristics permit to differentiate the corridors from the common route besides their layout:

-The parallel trajectories inside a corridor have all the same direction; -Separation between airplanes inside the corridor is done by the aircraft itself, without interference of controllers; -Corridors have entrance and exit points while in conventional routes aircraft can intercept the route where they want.

Flow corridors are not restricted to high flight levels, but their implementation seems to be more natural in this airspace since the corridors appear to be more advantageous for long distance flights.

In order to enable this type of traffic flow, with minimal or no interference at all of controllers, the airplanes must follow certain traffic specifications of the corridor and be able to maintain self separation and surveillance of the inside traffic. In such environment, dense traffic and self-controlled, the navigation should be based on accurate and reliable systems throughout the entire corridor.

Collaborative Navigation

The study performed by MAIAA Lab about the air streams theory tries to evaluate many aspects of this kind of traffic flow, such as surveillance, corridors rules, transitions between the corridor and the outside airspace and also the navigation performance inside the corridor. Professionals and students are involved in this study led by MAIAA, each one developing a different aspect of the research, and the part intended for my internship is related to collaborative navigation. Definition: For the purpose of this study, we should understand collaborative navigation as the utilization of measurements performed by surrounding aircraft and exchanged to improve the navigation performance of the each of them. Said in another way, an airplane surrounded by other airplanes should experiment a better navigation performance than an aircraft flying alone.

There are different ways to improve the navigation performance with the available navigation signals. The direction of this study is to assess the possibility of using TCAS and ADS-B signals, already made available in commercial aircraft, as means to improve navigation performance. Some studies using TCAS for improvement of positioning accuracy only for surveillance purposes have been already published [START_REF] Abeloos | Potential co-operations between the TCAS and the ASAS[END_REF][START_REF] Wu | TCAS-aided multilateration for terminal surveillance with improved accuracy[END_REF].

Figure 2. Triangulation with airplanes constellation

To improve navigation performance, an aircraft could use the information of distance between itself and other airplanes provided by TCAS, since the estimation of the absolute position of each airplane is available (ADS-B), and to create a dynamic constellation similar to GPS where the satellites are replaced by the surrounding airplanes (Fig. 2).

In the configuration above, instead of having the ephemeris with the satellite position we have an ADS-B message with the aircraft position and instead of measuring the distance aircraft-satellite using the satellite clock and received message time, we use the TCAS to measure the distance between aircraft. The challenges of this approach are not only the level of uncertainty of each measurement system, but also the lack of synchronization between measurements.

Position Determination Principle Triangulation with airplanes constellation

Cartesian coordinates are adopted here considering that distances between involved aircraft are very small with respect to the size of Earth. The Euclidean distance between two points in Cartesian coordinates can be expressed as follows:

A a z z y y x x r a r a r a a        2 2 2 ) ( ) ( ) (  ( 1
)
where A is the set of involved aircraft with n A  . Here, scalar a  is the TCAS measurement between a surrounding aircraft a and the reference aircraft. Here (x a , y a , z a ) are the Cartesian coordinates of the surrounding aircraft a and (xr,, yr, zr) are the Cartesian coordinates of the reference aircraft which have to be estimated.

There are different nonlinear methods [START_REF] Bancroft | An algebraic solution of the GPS equations[END_REF][START_REF] Oszczak | New algorithm for GNSS positioning using system of linear equations[END_REF][START_REF] Lim | Alternative Algorithms for Least Square Solution of Pseudo-range Measurements[END_REF] to solve the pseudo-range equation (1) for multiple surrounding airplanes (at least 3 of them) and to find the reference aircraft position. Among these methods, the linear approximation method appears with particular interest since the linear method is not only easy from the computational side but also brings the into view the concept of Dilution of Precision (DOP). This method produces a mean square error solution by solving iteratively linear approximations of the set of equations [START_REF] Youse | Nextgen Flow Corridors Initial Design, Pro-cedures and Display Functionalities[END_REF]. The coordinates that present the minimum mean square error value are considered to be the best estimation.

So, let at time t the variations at iteration k be: are the new estimates of the position coordinates. We will write:

       t k r P , t k r t k r t k r t k r t k r t k r t k r t k r t k r z z z y y y x x x , 1 , , , 1 , , , 1 , , ˆˆ            (2)
 , 2 , 1 ˆ, , , ,             k z y x P t k r t k r t k r t k r (3)
Let the separation coordinates considered at iteration k with aircraft a be given by:

t k r a t k a r t k r a t k a r t k r a t k a r z t z Z A a y t y Y x t x X , 1 , , , 1 , , , 1 , , ) ( ) ( )
(           (4) Let  , 2 , 1 ), ( ), ( ),
                        )) 2 ( ) 1 ( ( ) 1 ( ˆ)) 2 ( ) 1 ( ( ) 1 ( ˆ)) 2 ( ) 1 ( ( ) 1 ( , 0 , 0 , 0 (5)
The first order approximation of equations (1) at time t for point

1 1 1 , , ˆ   k r k r k r z y x is given by: A a Z Y x X t t k a r t k a r t k a r t k a r t k r t k a r t k a r t k a r a           , 1 , , , , 1 , , , , , 1 , , , , 1 ,
) (

     (6)
where

A a z t z y t y x t x t k r a t k r a t k r a t k a r            2 , 1 2 , 1 2 , 1 , 1 , ) ) ( ( ) ) ( ( ) ) ( (  (7)
Considering the approximate equations ( 6) for the n surrounding airplanes and writing them in matrix format, we can write:

t k r t k r k r P H D , ,    (8)
with: -If the number of surrounding aircraft is less than 3, this system has no particular solution (infinite number of solutions).

                      t k n r n t k r t k r k r t t D , 1 , , 1 2 , 2 , 1 1 , 1 ) ( )
(        (9) and                                           t k n r t k n r t k r t k r t k r t k r t k n r t k n r t k r t k r t k r t k r t k n r t k n r t k r t k r t k r t k r t k r Z Z Z Y Y Y X X X H , 1 , , 1 , , 1 2 , , 1 2 , , 1 1 , , 1 1 , , 1 , , 1 , , 1 2 , , 1 2 , , 1 1 , , 1 1 , , 1 , , 1 , , 1 2 , , 1 2 , , 1 1 , , 1 1 ,
-If there are 3 surrounding aircraft, there is a unique solution when H is invertible.

-If there are more than 3 surrounding aircraft, there is in general no solution. Then, minimizing

  k r T t k r t k r T t k r t k r D H H H P      , 1 , , , ˆ (11)
So the updated position at iteration k will be:

t k r t k r t k r P P P , , 1 ,     (12)
Iterations will be stopped at k*, either when k*=kmax or when condition:

   t k r P *, ˆ (13)
is satisfied, where  is a very small positive number.

Theoretical performance evaluation

In the previous subsection we described the iterative method to estimate the position of the reference aircraft, and as in the classical GNSS linear method, it is possible to determine the covariance of the position error to assess the accuracy of the considered navigation process. are also Gaussian variables with zero mean. Then:

  T P P E P      ) cov( ( 14
) or   1 1 ) ( ) ( ) cov(        H H H D D E H H H P T T T T (15)
When considering that the components of D  are statistically independent and are identically distributed, it is possible to write: where  2 is a pseudo-range error, then:

1 2 ) ( ) cov(     H H P T  (16)
Writing i h the i th column of H, we get:

    ij j T i T d h h H H     1 1 . ) ( ( 17
) Then: ] [ ) cov( 2 ij d P    (18)
and a position dilution of precision (PDOP) is given by:

33 22 11 d d d PDOP     (19)

Time frames and offsets

An important adjustment to be performed to the GPS model to fit with the current study is the removal of the clock offset parameters since here there is no such offset: the transmitter clock is also the receiver clock, since the TCAS generates the transmission time of interrogation and also the received time of reply.

The time frame of data reception is probably the biggest issue to be solved in the present approach. All the measurements, ADS-B and TCAS, are received by the transponder (XPDR) but not at the same time because the XPDR is not able to receive different sources simultaneously. So, the data received from other airplanes is not synchronized with the TCAS time frame of the reference aircraft. Therefore, the time frame has to comply with the unsynchronized issue and also with the update interval chosen equal to one second.

Figure 3. Time frame of received data from surrounding aircraft

All data is received with the same interval but at different times, while to deploy the GNSS model it is necessary to have synchronized information. For GPS systems it is not a problem because the receiver is able to receive many satellites at the same time and the distance measurement is made with the same signal used to transmit the satellite position. This implies that the information about distance satelliteuser and satellite position are synchronized. In our case we do not have the same advantage. All surrounding airplanes send their ADS-B message at different times and even the distance measurement and position of a given aircraft are not synchronized. So, to apply the proposed method we have to synchronize all this information with respect to a unique time line, the reference aircraft time line.

TCAS synchronization

Initial simulation trials were made just repeating the distance measurement in the reference time but it was easy to notice a delay in the final position results. Therefore, it was decided to estimate a correction for TCAS measurements considering the delay between TCAS time and reference time ( ttcas). The method proposed was to use the information of the previous TCAS measurement and estimate the speed that the distance increases or decreases. The result is a scalar value that when multiplied by t tcas and added to the last TCAS information it produces an estimated distance at reference time:

tcas k k k sync k t TCAS TCAS TCAS TCAS / ) ( 1     (20)
Figure 4 displays a situation where the TCAS measures the distance between the reference aircraft and the sur-rounding airplane #7 for the first time (k=1) at 0.27 seconds and the second measurement (k=2) at 1.27 seconds. The distance values are 4828 m and 4880 m respectively. Thus, to synchronize with the reference time at 2 seconds we need to calculate the TCAS delay ( t tcas = 1 -0.27 = 0.73 s), to multiply it by the difference between both measurements (4880 -4828 = 52 m) and to sum to the last TCAS data to find the estimated distance (520.73 + 4880 = 4917.96 m).

Figure 4. TCAS data synchronization

ADS-B synchronization

For the position data received from ADS-B message we have the same problem of synchronization but the extra information about the velocity vector makes easier the estimation of surrounding airplane positions at reference time. The ADS-B message contains the velocity vector (V x ; V y and V z ) of the transmitter aircraft which contains also some uncertainty but still helpful to estimate the synchronized position. Unlike the TCAS estimation, since we have now this velocity vector it is not necessary to use previous reference aircraft position information to estimate the current position, only the last position and velocity data from the transmitter aircraft. So, to estimate the synchronized position of the transmitter air-craft we have just to multiply the velocity component by the delay time and to add to the respective coordinate:

ads x k k t V x x    ˆ (21)
Example -The ADS-B message received brings the information that the x coordinate of the airplane #7 is 8945 m and the x component of its velocity is 142 m/s. If tads is equal to 0.61 seconds, so the estimated x coordinate for synchronized time is 9031.62 m (8945 + 0.61 142 = 9031.62 m). We shall apply the same for y and z coordinates.

Simulation Study

To simulate these situations, we adopted the following definitions and simplifications:

The reference time line is the clock of the reference aircraft, which starts on zero and has interval update of t = 1 second (the main line in figures 3, 4 and 5).

The interval update of all data has the same value, t = 1 second also.

For each data, the first received message time is zero plus a random value in seconds. Each data has its own random time start, so there is no order in receiving messages. For example, the TCAS data from airplane #2 could be the first message to be received, the next could be the ADS-B data from airplane #7 and so on.

Two types of scenario were specified: A corridor scenario, the main theory considered in the study, and a scenario named "free route", where all the trajectories do not follow straight lines or same direction.

The main objective to have both scenarios was to be able to compare the performance of the proposed approach in both situations and assess the possible effect caused by the organization of aircraft flows in the corridor scenario.

Figure 5. ADS-B data synchronization

Corridor Scenario

To set up a configuration, considering the reference aircraft performing a flight inside a corridor in cruise phase, initially are introduced the adopted minimum allowed separations and spacing between the involved aircraft:

 Minimum longitudinal spacing of 5 NM as performed in dense traffic airspace;  Vertical separation of 1000 ft as specified by RVSM airspace;  Horizontal separation of 4 NM adopted in RNP 1. Using the minimum separation parameters defined above, the next step is to define the corridor layout to be applied in the simulation for different scenarios. The adopted corridor layout has 8 parallel lanes in 3 different flight levels (FL), totaling 24 trajectories inside the corridor, as shown in figure 6 below.

Figure 6. Chosen corridor structure

For codification in the simulation algorithm, the lanes have been identified by numbers taken between 1 and 24. Indexes begin at the right upper corner and going left and from upper to lower lanes. The reference aircraft, which is supposed to perform collaborative navigation, is assigned to lane 12 (the red aircraft in the figure above).

Figure 7. Corridor airspace configuration with 23 surrounding aircraft

Simulations were performed for 8 different traffic corridor configurations, using respectively 23, 20, 17, 14, 12, 11, 8 and 5 surrounding airplanes (see figure 7 for a 23 aircraft scenario). The scenario with 12 airplanes was used to compare with the free route scenario, but it was not taken into account in the corridor configuration tendency study.

Free Route Scenario

This scenario has two main differences comparing to the corridor one:

the aircraft trajectories are not straight lines but follow arcs with different radius;

not all airplanes fly the same direction Since this scenario has been introduced with a comparative purpose, a single configuration with 12 surrounding aircraft has been considered (see figure 8) in the present paper.

Figure 9 shows a comparative picture between the free route and corridor scenario, both with 12 surrounding aircraft.

Simulation results

A post processing of the points estimated in relation to true position of the reference aircraft has been performed. Figures 10 and11 show these distributions for the free route and the corridor scenarios respectively.

Conclusions

In this communication, we proposed an additional navigation device totally independent of actual onboard navigation systems in a modern transportation aircraft and based on collaboration with surrounding traffic. Then, the design of a backup navigation system to face a failure of a classical navigation system appears feasible. This is very important in high traffic density situations. The accuracy performance of the new system, obtained through collaborative navigation, seems to be comparable to current GNSS systems, when there is enough surrounding traffic.

  be computed from equation (7) and the TCAS signals. The solution w.r.t.

  depends on the number of surrounding aircraft:

For

Figure 8 .Figure 9 .

 89 Figure 8. Free Route airspace configuration with 12 surrounding aircraft

Figure 10 Free route horizontal error distributionFigure 11 .

 1011 Figure 10 Free route horizontal error distribution

Figure 12 .

 12 Figure 12. Evolution of error with the number of aircraft


Xavier Delen
email: xavier.delen@institutoptique.fr

Loïc Deyra

Aurélien Benoit

Marc Hanna

François Balembois

Benjamin Cocquelin

Damien Sangla

François Salin

Julien Didierjean

Patrick Georges

Xavier Délen

Hybrid master

come

There is an increasing need for narrow linewidth, high average power ultraviolet (UV) laser sources for highresolution spectroscopy, remote environmental sensing, laser-induced fluorescence, semiconductor inspection, and LIDAR applications [1,2]. A standard way of designing efficient sources in the UV is to frequency convert near-infrared diode-pumped solid-state lasers [3,[START_REF] Diening | Proc. SPIE[END_REF]. In this case, obtaining high nonlinear conversion efficiency requires a large peak power in the fundamental wave. At a fixed repetition rate and average power, this results in a trade-off between attainable linewidth, related to the pulse duration, and the final nonlinearly converted average power.

Ytterbium-doped fiber amplifiers are particularly attractive to design high power sources at 1030 nm, providing compactness, efficiency, and high beam quality. In narrow linewidth regime (<100 MHz), the physical effect that limits the achievable power is stimulated Brillouin scattering (SBS). The use of large mode area (LMA) fibers and various SBS mitigating techniques has allowed scaling up of the peak power to around 1 kW in a 10 MHz bandwidth [START_REF] Jeong | [END_REF]. However, this peak power level is still too low for efficient frequency quadrupling to reach the wavelength of 257 nm. In this case, resonant cavities might be used to generate UV light efficiently [6], at the price of an increased complexity.

Recently, single-crystal fiber (SCF) amplifiers have demonstrated a great potential to scale fiber systems beyond their limits fixed by nonlinear effects, while retaining the possibility to generate high average powers [7]. This power scaling has been recently reported in a chirped-pulse femtosecond amplifier [8], allowing the generation of 1 mJ 380 fs pulses with a moderate stretching ratio impossible to attain with glass fibers due to self-phase modulation.

In this Letter, we demonstrate that Yb:YAG SCF amplifiers can be used to scale the power of SBS-limited 15-150 ns pulse ytterbium-doped fiber sources. Subsequent single-pass nonlinear conversion in bulk crystals permits the generation of Fourier transform-limited 15 ns pulses with 3.2 W of average power, or 150 ns pulses with 430 mW of average power at a repetition rate of 30 kHz at 257 nm.

The source is built following a master oscillator power amplifier (MOPA) architecture similar to [9], with an experimental setup depicted in Fig. 1. It starts with a CW single-frequency distributed-feedback (DFB) laser diode emitting 5 mW at 1030 nm with spectral linewidth below 1 MHz. First, this signal is amplified to 50 mW in a singlemode fiber preamplifier with 6 μm core diameter. A fiberpigtailed acousto-optic modulator (AOM) is then used to carve nanosecond pulses with a variable duration of 15-150 ns at a high repetition rate of 330 kHz in order to limit the amplified spontaneous emission (ASE). This signal is amplified in two stages of single-mode 6 and 15 μm core diameter Yb-doped fiber amplifiers to an average power of 400 mW. A second AOM further decreases the repetition rate to a few tens of kilohertz before coupling to the final LMA fiber amplifier. This amplifier is made of a 2 m long polarizing double-clad 40∕200 μm Yb-doped photonic crystal fiber pumped by a 25 W fiber-coupled laser diode system at 976 nm. The beam is then further amplified in an Yb:YAG SCF amplifier and converted into UV through a fourth-harmonic generation stage.

To study the SBS threshold in the last fiber amplification stage and characterize the amplified signal, we measured the output spectral linewidth. Part of the CW signal from the DFB diode is sent into a 2 km long optical fiber and used to establish a beating signal with the amplified output after the LMA fiber or the SCF amplification stage. It allows for heterodyne measurement of the narrow linewidth spectrum after different amplification stages [10]. The initial optical frequency is shifted twice by the AOMs through the amplification process. The beating frequency is observed at 90 MHz corresponding to the difference between the two shifts of 110 and 200 MHz of the AOMs. For example, Fig. 2 shows the beating signal and its Fourier transform for 150 ns pulses after the LMA fiber. In this case the linewidth at FWHM is 7 MHz, which is close to the Fourier transform limit.

The final repetition rate is first set to 10 kHz to study the SBS threshold as a function of linewidth at the output of the LMA fiber amplifier. The SBS threshold is detected by looking at the temporal profile of the output pulses, which is distorted when significant SBS occurs. Since the output pulse energy level is over 100 μJ, significant saturation is observed in the fiber amplifier, leading to possible modifications of the temporal profile both in intensity and phase in the time domain, as shown in Fig. 3. (right) The heterodyne measurement is therefore used to measure accurately the optical linewidth for each selected pulse duration. Figure 3 (left) shows that the large core diameter and short length of the fiber power amplifier allow to reach peak powers of 7 kW for 15 ns pulses (35 MHz linewidth at FWHM), with an SBS threshold decreasing to about 3 kW for pulse durations above 75 ns corresponding to linewidth below 13 MHz. Although at the state of the art for narrow linewidth fiber amplifier, the peak powers obtained are not sufficient for efficient frequency-quadrupling.

To further scale the power, an Yb:YAG SCF doublepass amplifier is added at the output of the LMA fiber amplifier. The crystal has an Yb doping concentration of 1 at. %, a diameter of 1 mm, and a length of 4 cm. It is set in a water-cooled mount at 16°C and pumped in single pass by a 75 W high brightness laser diode system at 940 nm coupled to a 105 μm fiber with NA 0.15. The free-propagating signal waist diameter inside the SCF is 400 μm, and the highly multimode pump beam is focused on the same diameter inside the crystal rod. In the SCF, the pulse energy is not limited by SBS but by damage threshold on the AR coated facets, occurring at 5 to 7 J∕cm 2 for 500 ps pulses. The repetition rate was therefore adjusted to 30 kHz to remain well below this limit for reliable operation of the system. Figure 4 (left) shows the output power as a function of pump power, for input pulse duration of 15 ns and input average powers ranging from 100 mW to 2.5 W. Between 95% and 97% of the pump power is absorbed by the SCF depending on the seed power. A maximum average power of 22.5 W is reached, corresponding to an optical-to-optical efficiency of 27%, and pulse energy of 750 μJ. It corresponds to a peak power higher than 50 kW, which is seven times higher than the peak power at SBS threshold in glass fiber for this pulse duration.

For 22 W average power, the optical spectrum obtained using an optical spectrum analyzer is displayed in Fig. 4 (right). The spectrum is essentially free of ASE with over 97% of the power in the main peak. Heterodyne linewidth measurements for pulses amplified in the SCF confirm that they remain Fourier-transform-limited, resulting in linewidths of 5-35 MHz at 1030 nm depending on the pulse duration. As shown in Fig. 4 (inset), the beam quality at the output of the SCF is diffractionlimited with M 2 < 1.1 in both transverse directions.

Frequency conversion to 257 nm is performed using two successive second-harmonic generation (SHG) stages. The first stage uses a 3 × 3 × 20 mm LBO crystal cut for noncritically phase-matched type I SHG at 1030 nm. It is AR coated at 1030 and 515 nm, and mounted in an oven at 187.8°C. The amplified beam is focused on a diameter of 120 μm inside the crystal. The power obtained at 515 nm and conversion efficiency are plotted as a function of fundamental power in Fig. 5 (left) for a pulse duration of 15 ns. The efficiency reaches 58% for a maximum power of 12.4 W at 515 nm, showing that the peak power is sufficient for efficient frequency conversion. The green beam quality remains excellent, with M 2 values lower than 1.1 in both directions.

The second SHG stage is based on an uncoated 4 × 9.5 × 9.5 mm BBO crystal cut for critically phasematched type I SHG at 515 nm. No temperature regulation of the crystal was implemented. To maximize the conversion efficiency, we used elliptical focusing down to 40 × 150 μm waists to remain within the angular acceptance of the crystal. The power obtained at 257 nm and conversion efficiency of the second SHG stage are displayed in Fig. 5 (right). A maximum power of 3.2 W is obtained, an unprecedented value to our knowledge for a single frequency laser at 257 nm without using nonlinear conversion in resonant cavities. Although the UV average power is stable for hours below 2 W, above this level, power instabilities related to absorption and thermal effects in the crystal start to appear. An optimized thermal management of the BBO crystal should allow increasing this 2 W threshold for appearance of thermal instabilities.

To highlight the trade-off between linewidth and available average power, Table 1 summarizes the results obtained for pulse durations ranging from 15 to 150 ns at 1030 nm after amplification. The generated UV power remains above 400 mW in all cases, making this source a versatile tool for LIDAR and remote sensing applications.

In conclusion, we demonstrate the use of a Yb:YAG SCF amplifier to scale the peak power of SBS-limited Yb-doped single-frequency fiber-based MOPA sources. Furthermore, we show subsequent efficient nonlinear conversion to the fourth harmonic at 257 nm. The performance of the presented hybrid source is currently essentially limited by the available pump power and damage thresholds of coatings of the SCF and the LBO crystal. Higher pump power should also allow power scaling to the 50 W level at 1030 nm without any further modification of the system architecture. Frequency conversion efficiencies could also potentially be improved with higher damage thresholds of the coatings. We believe that this type of architecture could be the basis for next generation LIDAR systems in the UV spectral range.

Fig. 1 .

 1 Fig. 1. (Color online) Experimental setup.

Fig. 2 .

 2 Fig. 2. Left: beating signal obtained after the LMA fiber amplifier for a 150 ns pulse. Right: Fourier transform of the beating signal.

Fig. 3 .

 3 Fig. 3. (Color online) Left: SBS threshold peak power at 10 kHz as a function of the output linewidth in the LMA fiber amplifier. Right: pulse shapes after the LMA fiber amplifier for different initial durations.

Fig. 4 .

 4 Fig. 4. (Color online) Left: output power as a function of pump power for the SCF amplifier (inset: beam profile at full pump power for a seed power of 2.5 W). Right: optical spectrum at maximum power and pulse width of 150 ns.

Fig. 5 .

 5 Fig. 5. (Color online) Left: power at 515 nm and conversion efficiency as a function of input power at 1030 nm. Right: power at 257 nm and conversion efficiency as a function of input power at 515 nm. The fundamental pulse duration is 15 ns.

Table 1 .

 1 Summary of the Hybrid Source Performances for Various Initial Pulse Widths at a Repetition Rate of 30 kHz

	Δt (ns)	15	50	100	150
	Δν (MHz)	40	21	11	6.2
	Power at 515 nm (W)	12.4 10.4	7.8	6.3
	UV power at 257 nm (W)	3.2	1.8	0.85	0.43

The authors acknowledge the DGA for the funding of this project.