Areski Hadjaz
email: areski.hadjaz@thalesgroup.com

Gaétan Marceau
email: gaetan.marceaucaron@thalesgroup.com

Pierre Savéant
email: pierre.saveant@thalesgroup.com

Marc Schoenauer
email: marc.schoenauer@inria.fr

Online Learning for Ground Trajectory Prediction

Keywords: Air Traffic Control, Trajectory Prediction, Hybrid System, Total-Energy Model, Black-Box Optimization, Machine Learning

This paper presents a model based on an hybrid system to numerically simulate the climbing phase of an aircraft. This model is then used within a trajectory prediction tool. Finally, the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) optimization algorithm is used to tune five selected parameters, and thus improve the accuracy of the model. Incorporated within a trajectory prediction tool, this model can be used to derive the order of magnitude of the prediction error over time, and thus the domain of validity of the trajectory prediction. A first validation experiment of the proposed model is based on the errors along time for a one-time trajectory prediction at the take off of the flight with respect to the default values of the theoretical BADA model. This experiment, assuming complete information, also shows the limit of the model. A second experiment part presents an on-line trajectory prediction, in which the prediction is continuously updated based on the current aircraft position. This approach raises several issues, for which improvements of the basic model are proposed, and the resulting trajectory prediction tool shows statistically significantly more accurate results than those of the default model.

I. INTRODUCTION

Trajectory Prediction (TP) is the core component of automated systems in Air Traffic Control (ATC). Many functionalities of the Decision Support Tool directly rely on an accurate TP: controller posting, workload estimation, arrival sequencing, loss of separation detection, and conflict resolution, to name the most prominent ones. As a consequence, TP is the weakness of the current automation ATC systems and a major issue in the ATC research community, even more with the new paradigmatic shift toward 4D trajectories in both the SESAR Joint Undertaking and the NextGen project. The challenge is to reduce the uncertainty of the prediction of the aircraft states on a temporal horizon of at least 20 minutes. To achieve this, the information of the current state of the aircraft and its environment shall be reliable. As a matter of fact, the Flight Management System (FMS) has access to the measurements from the sensors of the aircraft and creates its own TP, which is updated frequently. Therefore, we should expect that this TP is the most accurate one. Thanks to the Data Link, it should be possible for the ground control to receive data from this onboard TP. Unfortunately, today, this promising technology is not implemented in all aircrafts. Moreover, a decision support tool requires the capacity to efficiently generate 'what-if' scenarios. Today, it is unrealistic to receive many trajectories from the on-board system of every aircraft in the airspace and merge them into different airspace scenarios in real-time. Therefore, the ground TP is still an essential component of the future air traffic control systems.

A trajectory can roughly be separated in three phases: climbing, level flight, and descent. Phases with altitude changes are the most difficult to handle due to the differences of aircraft performances and the effect of the weather conditions. Operationally, the controllers of approach control centers, who deal with climbing and descent phases, must have a good representation of the vertical evolution of the aircraft to ensure their separations. In most airspace, this separation shall be at least 10 flight levels. For these reasons, this study will focus mainly on the prediction in the vertical axis.

The paper is organized as follows: Section II surveys some recent works in the domain of trajectory prediction. Section III details the proposed model. The validation experiments are described in Section IV, while the on-line prediction method is experimented with in Section V. Section VI summarizes the paper and gives some hints for further research in the area of trajectory prediction.

II. RELATED WORK

Research on TP includes a set of methodologies for its specification, implementation and evaluation. First, the Eurocontrol Specification for Trajectory Prediction gives the requirements for an operational TP and a validation methodology useful for the Air Navigation Service Providers in their choice of a new ATM system compatible with SESAR. This document covers many functionalities like flight plan and clearance processing, airspace constraints and real-time monitoring. Unfortunately, it does not specify the accuracy requirements. Besides, [START_REF] Musialek | Literature Survey of Trajectory Predictor Technology[END_REF] provides an important literature survey of trajectory prediction technology with 282 reviewed documents and 20 selected for further studies. From the selected set, many implement the point-mass model where the rotational moments are not modeled. This is an acceptable assumption for airline aircraft. Also, [START_REF] Coppenbarger | Climb trajectory prediction enhancement using airline flight-planning information[END_REF] and [START_REF] Swierstra | Common Trajectory Prediction Capability for Decision Support Tools[END_REF] enumerate the principal difficulties inherent to TP, i.e. the uncertainty on the input data, the controller and pilot intents, and quantify the errors accordingly. The input data mainly refers to aircraft characteristics which are given in the Base of Aircraft DAta (BADA) Aircraft Performance Model. Still, these values are only nominal and can be different from the real situation. Therefore, the SESAR project 5.5.2 [START_REF] Chistien | Final report of project 05[END_REF] concludes that sharing airline operational control data, like the mass and the speed schedule, could provide quick improvement to ground trajectory prediction, with limited investment. As a matter of fact, these parameters are determinants for computing the TP with a point-mass model.

More generally, a parametric approach refers to a model based on flight equations and aircraft characteristics. The point-mass model is an example of a parametric approach. As we will see, many parameters can be used to tune the model to reality. To generate trajectories, [START_REF] Gallo | Trajectory Computation Infrastructure Based on BADA Aircraft Performance Model[END_REF] uses six ordinary differential integrators depending on the longitudinal motion instruction. Similarly, [START_REF] Glover | A Stochastic Hybrid Model for Air Traffic Control Simulation[END_REF] and [START_REF] Kamgarpour | Hybrid Optimal Control for Aircraft Trajectory Design with a Variable Sequence of Modes[END_REF] use hybrid systems to model the change of differential equations according to the control law and the aircraft states. A discrete space, e.g. an automaton, is defined to represent the modes of the system. Every mode defines the differential equations and creates trajectories in a continuous space. As seen in [START_REF] Glover | A Stochastic Hybrid Model for Air Traffic Control Simulation[END_REF], this model is well-suited for implementing BADA. [START_REF] Alligier | Energy Rate Prediction Using an Equivalent Thrust Setting Profile[END_REF] exposes a technique to find a general thrust settings, i.e. a control law, that could be used in such framework. The idea of fitting the mass parameter of BADA on a few past points is used. Results on the accuracy of this generic control law in the vertical plane should be expected in a near future.

A common flaw in parametric approaches concerns the nominal values used for every parameters. As shown on figure 1, the effect of the mass parameter clearly shows the importance of tuning these parameters. As an example, 400 seconds after the takeoff, there is already a difference of approximately 250 flight levels between the minimal and the maximum mass which is enormous for an application like the TP. This can easily become a burden when the model is rather complex. To overcome this inherent difficulty, non-parametric approaches are studied in order to obtain an aircraft model from the past trajectories. Non-parametric approaches rely on machine learning and statistical inference: [START_REF] Fablec | Using Neural Networks to predict aircraft trajectories[END_REF] uses neural networks, [START_REF] Alligier | Programmation évolutionnaire appliquée à la prévision de trajectoire[END_REF] uses genetic programming in order to learn the structure of the variables of a linear regression and [START_REF] Ghasemi | Prédiction de trajectoires d'avions à l'aide de la régression floue[END_REF] uses fuzzy regression with k-nearest neighbor. The main drawback of non-parametric approaches is that it requires lots of historical data and the model is learned for a specific context and can hardly be generalized because of the airspace constraints, e.g. aircraft following a Standard Instrument Departure.

Combining both approaches is an interesting research question addressed by [START_REF] Alligier | Energy Rate Prediction Using an Equivalent Thrust Setting Profile[END_REF] and [START_REF] Crisostomi | Combining Monte Carlo and worst-case methods for trajectory prediction in air traffic control: a case study[END_REF]. The latter combines Monte-Carlo Simulation and worst-case scenario for modifying the parameters of BADA while integrating a wind model. However, this work is limited to the descent phase and the experiments are performed on trajectories obtained by simulation. The main contribution of our work on this question is to show that tuning the parameters of BADA during the progress of a flight is a way to overcome the difficulties of parametric approach and avoiding to learn a model from the historical data. To the best of our knowledge, no work has answered to this question with results of the vertical plane.

III. MODEL

The basic idea is to create an hybrid system to generate a trajectory in the vertical plane using BADA. The system describes the transitions between the different modes in function of the states and the control laws. Following the mode, a given differential equation is integrated to obtain the trajectory. This section states the assumptions, defines the model and gives an integration schema for generating the trajectories.

A. Assumptions

Depending on the requirements of the simulation, one may choose different levels of complexity for the coordinate systems and the aircraft control system. Since we consider only the climbing phase, the system will evolve on a short period of time on a limited geographical area. From these considerations, we assume that the flat earth model is a reasonable approximation (cf. [START_REF] David | Fundamentals of Airplane Flight Mechanics[END_REF]). This model assumes that the earth is flat, non-rotating and can be defined as an inertial reference frame. Also, the gravitational force is constant and perpendicular to the ground. The atmosphere is at rest relative to the earth and atmospheric properties solely depend on the altitude. With these assumptions, one can define the differential equations of kinematics and dynamics. For the kinematics, the equations imply solely that the displacement is proportional to the speed projected on the vertical or horizontal plan. For the dynamics, we need to model the thrust, the drag, the lift and the weight. To our knowledge, BADA is the most complete model describing the capabilities of many aircraft types.

B. Total-Energy Model

This section relies on the BADA 3.10 User Manual [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF]. BADA is based on the total-energy model where the rate of work is equal to the rate of potential and kinetic energy. Here, the rate of climb is obtained by controlling the speed and the throttle. From [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF] p.14:eq.3.2-7 and by transforming the quantities into functions of altitude h, speed V and mode q, the rate of climb is defined by:

ḣ(V, q) = T (h) -∆T T (h) (T hr(h) -D(V, h)) • V m • g •f (V, h, q)
(1) where the function T (h) is defined in BADA 3.10 User Manual [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF] p. 10, T hr(h) at p.22, D(V, h) at p.20 and f (V, h, q) at p.15-16. As a first approximation, the mass m, the gravitational acceleration g and the temperature differential ∆T are constant during the climb phase. According to [START_REF] David | Fundamentals of Airplane Flight Mechanics[END_REF], this is a reasonable assumption since the climb phase lasts around 20 minutes. The terms h, V and q are variables that evolve with the system and one must specify the evolution of the two remaining independent variables. However, the acceleration in function of the aircraft dynamic is not specified in BADA. From [START_REF] David | Fundamentals of Airplane Flight Mechanics[END_REF], it is given by:

V (h, ḣ) = 1 m (T hr(h) -D(V, h) -m • g • sin(γ)) (2)
where sin(γ) = ḣ V . So, from eq. 2, we can see that the acceleration evolves independently of the mode given the rate of climb. Also, when ḣ is high, sin(γ) is high and the acceleration V is low. This goes along with the hypothesis the total-energy model. Finally, the next section will present the evolution of the variable q.

C. Mode Definition

In this section, we define the modes and their transitions based on BADA 3.10 User Manual [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF]. First, let Q = {CAS, M ACH} × {LOW, HIGH} × {DEC, CST, ACC} be the mode space. The first two modes, transition altitude and tropopause, depend solely on the altitude:

q 1 (h) = CAS if h ≤ H trans MACH otherwise q 2 (h) = LOW if h ≤ H trop HIGH otherwise
where H trans is the transition altitude and H trop is the tropopause geopotential pressure altitude. Finally, the last feature q 3 is the mode of acceleration. The simplest controller of this mode is given by:

q 3 (V, h) =      ACC if V ≤ V * (h) - DEC if V ≥ V * (h) + CST otherwise
where ACC is acceleration, DEC is deceleration and CST is constant speed and V * be a target speed at altitude h where V tends to converge. V * can be chosen according to the speed schedule defined in the Airline Procedures Model of BADA 3.10 User Manual [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF] p.29, where three target speeds

(V 1 , V 2 , M
) are required as parameters. The nominal values can be found in the airline procedure files of BADA. Finally, ∈ R is a threshold value to avoid jitter. Next, the energy share factor function f takes its values according to the following flight conditions: [START_REF] Fortin | Analyse Numerique pour ingenieurs[END_REF] Deceleration in climb where V M is the Mach speed and V CAS is the calibrated speed. These flight conditions can be defined in terms of conjunctions of q 1 , q 2 and q 3 . From [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF] p.15-16, f is discontinuous when q jumps from one mode to another. Also, f is bounded and therefore, with the flight envelope, one can also bound ḣ and V . One must pay attention to the Zeno behavior (cf. [START_REF] Sastry Lygeros | Hybrid Systems: Foundations, Advanced Topics and Applications[END_REF]) where a system can make an infinite number of jumps in a finite amount of time. With these considerations in minds, one can use a numerical procedure to integrate eq.1 and eq.2.

1) Constant V M above tropopause, 2) Constant V M below tropopause, 3) Constant V CAS above tropopause, 4) Constant V CAS below tropopause, 5) Acceleration in climb,

D. Trajectory Generation

In order to generate the trajectory, one must specify eq.1 and eq.2 as functions of time. With respect to BADA, let T hr(h(t), t) = T hr(h(t)), D(h(t), V (t), t) = D(h(t), V (t)) and T (h(t), t) = T (h(t)), that is the evolution of the thrust, the temperature and the drag are time-invariant. Moreover, the aircraft dynamic functions shall be specified with respect to the flight envelope constraints. In this study, we choose a nominal thrust function. In eq.1, T hr is replaced by the maximum thrust T hr max and the whole equation is multiplied by a reduced climb power coefficient C red , which is supposed to give realistic profiles (cf. [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF] p.24). To simulate the system, a common fourth-order Runge-Kutta method is used. Let f 1 (t, h, V, q) = ḣ(V, q) and f 2 (t, V, h, ḣ) = V (h, ḣ). Then, one obtains the following integration scheme:

dh 1 = f 1 (t n , h n , V n , q(V n , h n )) dv 1 = f 2 (t n , V n , h n , ḣn ) dh 2 = f 1 (t n + ∆t 2 , h n + dh 1 ∆t 2 , V n + dv 1 ∆t 2 , q(V n + dv 1 ∆t 2 , h n + dh 1 ∆t 2
))

dv 2 = f 2 (t n + ∆t 2 , V n + dv 1 ∆t 2 , h n + dh 1 ∆t 2 , dh 2 ) dh 3 = f 1 (t n + ∆t 2 , h n + dh 2 ∆t 2 , V n + dv 2 ∆t 2 , q(V n + dv 2 ∆t 2 , h n + dh 2 ∆t 2
))

dv 3 = f 2 (t n + ∆t 2 , V n + dv 2 ∆t 2 , h n + dh 2 ∆t 2 , dh 3 ) dh 4 = f 1 (t n + ∆t, h n + dh 3 ∆t, V n + dv 3 ∆t, q(V n + dv 3 ∆t, h n + dh 3 ∆t)) dv 4 = f 2 (t n + ∆t, V n + dv 3 ∆t, h n + dh 3 ∆t, dh 4 ) h n+1 = h n + ∆t 6 (dh 1 + 2dh 2 + 2dh 3 + dh 4 ) V n+1 = V n + ∆t 6 (dv 1 + 2dv 2 + 2dv 3 + dv 4 ) t n+1 = t n + ∆t
where the initial conditions t 0 , q 0 , h(t 0 , q 0 ) = h 0 , V (t 0 , h 0 ) = V 0 are given. The choice of the Runge-kutta method is justified by the jumps in the function f . As a matter of fact, this method will minimize the impact of a jump during a timestep by averaging the variations at the beginning, twice at the middle point and at the end of the timestep (cf. [START_REF] Fortin | Analyse Numerique pour ingenieurs[END_REF]). Nevertheless, stability questions shall be addressed in the later in order to validate the approach.

E. Parameter Tuning

Now that the trajectory generator is defined, we would like to tune the model parameters according to observations. First, there should be a trade-off between the number of parameters and the capacity of the model to approximate real trajectories. As a bad example, instead of using the default thrust controller defined in [START_REF] Nuic | User Manual for the Base of Aircraft Data(BADA) Revision 3.10[END_REF], we tried to find the optimal controller for the acceleration mode for a given trajectory. The resulting controller was a Bang-Bang control function where we have observed that the mode switches between accelerating and decelerating at every timestep. Clearly, this function can approximate any real trajectories, but does not generalize from one trajectory to another and does not reflect the real behavior of the aircraft. So, from eq.1 and eq.2 and the speed schedule V * , the parameters m, V 1 , V 2 ,M and ∆T are good candidates since they are time-invariant, contrary to a discrete control law. Moreover, their values are fixed to nominal values and BADA explicitly suggests to tune them.

One way to tune the parameters is to minimize the position errors between the TP and the real trajectory. This leads to an optimization problem where the function to minimize is expressed by the trajectory generation scheme. Here, we decide to use a black-box optimization algorithm that will perform parameter estimation on the hybrid system. Due to the relations between the parameters and the complexity of the BADA models, we assume that parameter estimation is a non-trivial optimization problem.

F. Black-Box Optimization

As already mentioned, parameter tuning pertains to nonconvex black-box optimization, and several methods could be used to tackle it. Furthermore, no information whatsoever is available regarding the modality of the objective function with five parameters. As a matter of fact, it is unimodal when the mass parameter or the differential temperature are the only ones to be tuned. These have an effect on the whole trajectory. But, the speed parameters have a local effect depending on the speed schedule. Moreover, the differences between the speed parameters will induce acceleration phases that will transform the trajectory. Finding the best value for a speed parameter to fit locally the trajectory will create a local optimum that could be worse than finding the two speed parameters that will avoid an acceleration phase that is not undertaken in reality. Even if we did not prove that the objective function is multimodal, we can think that it can possibly be the case. Finally, the objective is non-differentiable (or at least the analytical derivative is out of reach). Hence general-purpose derivative-free optimization method is required.

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF] is today a state-of-the-art derivative-free optimization method that has demonstrated outstanding performances for problems up to a few hundred variables, in several official comparisons (see, among others, the CEC 2005 challenge [START_REF] Hansen | Comparison of Evolutionary Algorithms on a Benchmark Function Set[END_REF], and both Black Box Optimization Benchmark workshops at ACM-GECCO 2009 [START_REF] Hansen | Comparing Results of 31 Algorithms from the Black-Box Optimization Benchmarking BBOB[END_REF] and 2010 [START_REF]BBOB workshop[END_REF]), as well as on a large number of real-world applications [START_REF] Hansen | References to CMA-ES Applications[END_REF]. CMA-ES is an Evolution Strategy [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technisher Systeme nach Prinzipien des Biologischen Evolution[END_REF], [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF] that uses Gaussian mutation with adaptive parameter setting. A Gaussian mutation is defined here by its step-size and its covariance matrix. The step-size is increased (resp. decreased) if the cumulated path of the current best solution is smaller (resp. larger) than that of a random walk, and in the original version [START_REF] Hansen | Adapting Arbitrary Normal Mutation Distributions in Evolution Strategies: The Covariance Matrix Adaption[END_REF], [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF], the covariance matrix was updated by adding a rank-one matrix with eigenvector the direction of progress. An improved version with rank-µ update was later proposed [START_REF] Hansen | Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES)[END_REF], and several additional variant made it more and more powerful. The most recent version is the so-called bi-pop-CMA-ES [START_REF] Hansen | Benchmarking a BI-Population CMA-ES on the BBOB-2009 Function Testbed[END_REF], that evolves both a large and a small population and outperforms previous versions in case of multimodal functions. All source code is available on the author's web page (http://www.lri.fr/ ∼ hansen/index.html), in different programming languages, including the bi-pop version. Using CMA-ES for parameter estimation is then straightforward, and amounts to interfacing the objective function for BADA TP with the core CMA-ES program, after eventually normalizing its parameters.

IV. MODEL VALIDATION

This section deals both with the empirical validation of the nominal model and the estimation algorithm used to fit the parameters on the whole trajectory.

A. Dataset

To validate our model, we will use a dataset composed of 262 real departure trajectories of A320. These trajectories have been recorded via a radar systems during one month. For one trajectory, there is a data vector at every 5 seconds composed of the aircraft position, the rate of climb and the true airspeed. Then, the top of climb is calculated as the first highest point of the trajectory. Also, these trajectories do not have a rate of climb that is equal to zero for more than 30 seconds, they shall reach a cruise level at least of 300 FL and their durations are at least of 1100 seconds. These filtering conditions permits to keep trajectories that are not affected by a level flight clearance.

B. Methodology

To assess the model, we evaluate the position error between the TP and the real trajectory. Let H(θ; s 0 ) be the sequence of altitudes generated by the simulation of the hybrid system with the parameters θ

= [m ∆T V 1 V 2 V M ]
T and the initial condition vector s 0 = [t 0 q 0 h 0 V 0 ] T . Then, the measure is simply the sum of absolute errors. In the case of parameter estimation, we search in the feasible space of parameters Θ, the point θ ∈ Θ that minimizes this measure. So, we have the optimization problem:

θ * i,j = argmin θ∈Θ j n=i |H(θ; s 0 ) n -T n | (3)
where H(θ, s 0 ) is the sequence of ordered altitudes obtained for given parameters θ and initial conditions s 0 and T is the sequence of observed altitudes from a trajectory of the dataset. We suppose that (j -i) ≤ |H(θ, s 0 )| ≤ |T | and that the timestamp associated to the point H(θ, s 0 ) n is the same than T n . When fitting the whole climbing phase, some difficulties with eq.3 may arise. First, we need to give a termination criterion when simulating the hybrid system that is to stop when the TP reaches the level flight. Depending on the parameters, the cardinal of the resulting sequence of points will be different. So, it might be necessary to add points to the real trajectory T since the TP can reach the top of climb after T, e.g. a high value for the mass parameter. To compute this error, we simply add points at the level flight until the TP reaches it. Inversely, we add these points to the TP if it reaches the top of climb before the real trajectory. Besides, i < j are bounds to restrain the optimization on any subset of contiguous points of the trajectory. Thereafter, this notation will be useful for the online predictor.

C. Results

Table I shows the evolution of the mean errors and the standard deviation with time for both models: with nominal values and with tuned values. The first important evidence is the inaccuracy of the Total-Energy Model to model the positions at the beginning of the trajectory. As a matter of fact, for the tuned values, the errors at 2 minutes are the highest. An explanation consists in the fact that the aircraft states (position, speed, heading) change rapidly at the beginning of the trajectory and the selected five parameters are not sufficient to capture this complexity. Furthermore, the optimization of the equation eq.3 has a global scope and so, the errors generated by local behaviors of the aircraft are ignored in favor of the common behavior. This depends directly on the selected values of i,j in θ * i,j , which in this case are i = 0 and j = t toc where t toc is the time at top of climb. From our dataset, this common behavior happens around 10 minutes where the errors are the smallest after an acceleration phase which happens around 5 minutes. As an example, on the figure 2, we can distinguish three main behaviors: the initial climb from 0s. to approximately 200s, a short acceleration phase from 200s. to 300s. and the common behavior after 300s. From figure 3, we can see that the initial climb is characterized with high variability in the rate of climb and an acceleration phase between 50s. and 100s., shown by a huge decrease in the rate of climb. During this phase, the predicted rate of climb is far from reality for both parameter sets. Thereafter, both models capture the acceleration phase and finally, they average the rate of climb during the common behavior, which fits well the positions as shown on figure 2.

Another interesting result is the evolution of the standard deviation, which is in parentheses, for the model with nominal values. It increases with time from 2 min. to 15 min. and afterwards, it seems to stabilize around 9 FL. This can be interpreted as an uncertainty cone, which is often used in the air traffic community, but here, we can see that the cone stops to grow at 15 minutes and becomes a corridor of uncertainty. First, we understand that the cone grows with the flight envelope, but afterwards, the uncertainty is bounded by the fact that the trajectories, as functions of time, are strictly increasing because of the filtering conditions of the dataset (cf. Section IV-A) and the upper bound that is the cruise level.

V. ONLINE TRAJECTORY PREDICTOR

In this section, we present an online trajectory prediction which uses the observed positions of a current flight to tune at the same time the parameters of the model used to predict the rest of the trajectory. Different kind of algorithms can be used to undertake this task. Traditionally, a probabilistic approach is used in this online configuration, e.g. Kalman Filter, where the uncertainty is explicitly modeled. These performs very well on short periods of time but, for longer period like in this context, the linearity of the model could be too limited and should be subject to future experiments. In our approach, we use the BADA model in conjunction with the optimization algorithm CMA-ES in order to fit the parameters of the model with an objective function that defines a distance between the observed positions and the modeled one. As we will see, the problem of overfitting, which is studied by statistics and machine learning, arises in our context. The current section presents the difficulties and the solutions chosen with the associated results.

A. Design of Experiment

In order to verify that the idea of an online predictor, as the one mentioned previously, is valid, we must do an empirical evaluation of the chosen algorithm by replicating the same context. The main hypothesis is that from the observation set, we can determine the values of parameters that will be fitted to the current flight. So the trajectory is separated in two subsets: the observed altitudes from the beginning to the present and the future altitudes from the present to the top of climb. As in the subsection IV-C, we use metering points to evaluate the quality of the prediction by computing the errors between the predicted altitudes and the real ones. To distinguish if a set of errors is statistically greater than the other one, we use a Wilcoxon signed-rank test. The null hypothesis of this test is that two related paired samples come from the same distribution. In our case, this test is a relevant choice because two approaches are tested on the same trajectory dataset. As usual, we reject the null hypothesis if the p-value is lower than 0.05.

B. Methodology

The most naive way learn the parameters of the model from the observed altitudes is to directly apply Eq.3 from subsection IV-B and to apply them to generate the rest of the trajectory. By doing so, the default model is always better than the fitted one with a significant p-value. The reason behind this result is simply that the fitted model does not generalize over all the behaviors of the aircraft. In other words, it is fitted only to the behavior captured in the observations. This problem is referred to overfitting in Machine Learning litterature. To circumvent this problem, we must think of the trajectory as a time serie where the observations arrive with an determined order, i.e. the temporal order. So, at first, we will always observe the initial climb where we know from Table I that the inaccuracy of the model is the greatest. Furthermore, we are more interested by the parameter values that fit better the positions near the present time than at the beginning of the trajectory. To this end, we will add a weighting vector α that will penalize more the errors that are near the present time. But still, this is not sufficient because, depending on the present time, some parameters will not have any effect on the trajectory. As a matter of fact, from 0 to FL60, a predefined schedule is applied and only the mass parameter has an effect in BADA. The scope of the parameter V 1 is from FL60 to FL100, the scope of V 2 is from FL100 to the transition altitude (around FL277) and finally, the scope of V m is over the transition altitude. Furthermore, we add the constraint that V 2 is greater than V 1 to the optimization problem. To avoid that the optimization algorithm assigns them some arbitrary values resulting in unrealistic trajectories, we use a regulation method that penalizes any deviation from the default parameters. A meta-parameter λ is associated to the weight of the penalty, which controls the tradeoff between exploration and exploitation. Consequently, the resulting objective function is:

θ * 0,t = argmin θ∈Θ   t i=0 α i |H(θ; s 0 ) i -O i | + λ |θ| i=0 θ i -θ d i   (4
) where θ d is the default parameters of BADA. Notice that the penalty occurs on the parameters space, but could also be applied on the trajectory space because the mapping from the parameter space to the trajectory space is not linear in the sum of differences of altitudes.

Finally, in order to set the value of λ, we use a crossvalidation approach where we partition the observation set in two: the learning set and the validation set. We choose the validation set to be just before the current altitude. One must notice that the samples are not independent and identically distributed and that we create a bias in favor of the points located just after the current altitude. Because of our extrapolation context, a bias is inevitable and this one seems the most justifiable one in order to gain accuracy in predicting the future positions. Figure 4 shows the partition of the trajectory. The cross-validation technique used in this study consists in learning the parameters of the model on the learning set with the objective function and to use these parameters for generating the points on the validation intervals. Then, we compare these points with the real ones. We do it for multiple values of λ and we choose the parameter values where the validation error is the lowest to generate the rest of the trajectory.

C. Results

The approach is validated on the same dataset than subsection IV-A. We choose three different time slices in order to represent the online aspect of the method. The validation set size is fixed to 36 points or 180 seconds. This choice must do the trade-off between the validation purpose of avoiding overfitting which could lead to a large validation set size and the learning purpose of finding the best parameter values which could lead to a large learning set size and therefore, a small validation set size. At least, the validation set size must be higher than the acceleration phase, where the local behavior is the most different from the global one. In the learning objective function (cf. Eq.4), we choose a linear weight function where α i = i t-1 . The initial lambda value is arbitrarily set to 100 and are doubled until the penalty is high enough so that the fitted values equal the default ones. Then, the parameter values generating the lowest validation error are chosen.

Also, to avoid that the algorithm changes the parameter values based on a poor learning performance, we arbritrarily set a threshold error at 5 FL, which is higher than the results at subsection IV-C. When the threshold is exceeded, the BADA default values are chosen.

Table II, III and IV show the results of the proposed methodology for each time slice. At P = 400s, the two model performances are not significantly different as shown with the high p-values. We can see that our approach increases the accuracy by 1 FL in average for the metering points at 2 minutes and 5 minutes after the current time slice with a pvalue significantly under 0.05. For the metering points at 10 minutes, the p-value is higher than 0.05 and the difference is not statistically significant due to the high standard deviation values. This can be explained by the fact that the model is not very accurate during initial climb (cf. model validation section) and the validation set covers the acceleration phase. Also, because the learning error is too high, the algorithm can choose the BADA default values. So, the default choice ratio is around 20% which is rather high.

At P = 500s, the fitted model performs better at 2 and 5 minutes after the current position with small p-values. For 10 minutes, the two models are not significantly different because of the high value of the standard deviation. In fact, this can be interpreted as the two models are equally affected by the uncertainty around the possible maneuvers of the aircraft. In order to perform better, more information on the flight intents are required to reduce the variability in the trajectories. Here, the ratio of the default choice is 16%.

At P = 600s, the results are similar to P = 500s. The reason is that the aircraft keeps the same behavior between 500s. and 600s. which is different from the behavior at 400s.

There is some kind of regularity that explains the fact that the prediction is enhanced up to 5 minutes. This regularity is captured more easily by the learning algorithm when the behavior is stable during the validation interval. In this case, the ratio of the default choice is 14%.

VI. CONCLUSIONS AND FUTURE WORK

In conclusion, this article presents a flight model for the climbing phase defined as a hybrid system based on BADA. An integration scheme is defined in order to generate the trajectories from this system. Then, tuning parameters are identified in order to be used in the online context. The method is validated through measuring the vertical errors between real trajectories and generated ones: both for default and fitted parameter values. In the validation context, fitting is done on the entire trajectory i.e. with total information. The measured errors are considered as the accuracy limit of this five parameters model. Afterwards, the model is applied in the online context, which evolves with time. Known altitudes are used to fit the parameters and then, these are used to predict the remaining points. To avoid overfitting, the known points are partitioned in a learning set and a validation set. The validation set is used to fit the regulation parameter, which penalizes the deviation from the default values. Results shows that the initial climb, which is before the main acceleration phase, is not modeled accurately in order to fit the parameters. On the contrary, when the flights adopt a common behavior after this acceleration phase, the online learning method increases the accuracy of the trajectory prediction. The gain is about 1 FL for 2 minutes and 5 minutes after the current time. After that, the two models are not significantly different because of the huge uncertainty on the trajectory. This study shows that the uncertainty becomes too important between 5 and 10 minutes with minimal information. Consequently, without further information on the flight intents and airspace constraints, ground trajectory prediction is not accurate enough for automated tasks such as conflict resolution. Furthermore, this study confirms the need to use the aircraft derived data to feed the BADA model in order to build a ground trajectory prediction as the foundation stone of automated Air Traffic control systems. The next step should be the identification of the relevant onboard data improving significantly the trajectory predictability on ground.

Figure 1 .

 1 Figure 1. The Effect of the mass parameter

Figure 3 .

 3 Figure 2. Trajectory Fitting

Figure 4 .

 4 Figure 4. Online Prediction

TABLE II COMPARISON

 II OF ONLINE MODELS AT P = 400S

	Time after takeoff	Nominal (FL)	Tuned (FL)	p-value
	2min.	3.3029 (2.6698)	3.1699 (2.6740)	0.3401
	5min.	6.7553 (5.6084)	6.5518 (5.6578)	0.6726
	10min.	8.7851 (7.0757)	9.1846 (7.5687)	0.4541
		TABLE III
	COMPARISON OF ONLINE MODELS AT P = 500S
	Time after takeoff	Nominal (FL)	Tuned (FL)	p-value
	2min.	4.0406 (3.2758)	3.2834 (2.7237)	5.612e-4
	5min.	8.1290 (6.0885)	7.0567 (4.8281)	0.02049
	10min.	9.0872 (7.0085)	9.4205 (6.7658)	0.7939
		TABLE IV
	COMPARISON OF ONLINE MODELS AT P = 600S
	Time after takeoff	Nominal (FL)	Tuned (FL)	p-value
	2min.	4.5110 (3.4354)	3.5912 (2.4845)	1.289e-05
	5min.	6.7936 (4.9209)	5.7231 (4.0456)	1.289e-03
	10min.	8.5131 (6.6410)	9.4992 (7.8805)	0.09098

VII. ACKNOWLEDGMENTS

This author Gaétan Marceau is funded by the scholarship CIFRE 710/2012 established between Thales Air Systems and INRIA-Saclay, and the scholarship 141138/2011 from the Fonds de Recherche du Québec -Nature et Technologies.


R Al-Nazer

V Cattin

M Montaru

BROADBAND IDENTIFICATION OF BATTERY ELECTRICAL IMPEDANCE FOR HEV

Keywords: battery impedance, spectroscopy, broadband signals, Li-ion battery, coherence, identification

In recent years, Li-ion batteries have been proposed as an essential element for hybrid electrical vehicles (HEV) and electrical vehicles (EV). In such applications, the most possible accurate estimation of the battery states is needed to optimize its operation. Accordingly, battery electrical impedance is known to be able to provide useful states information. Though that electrical impedance spectroscopy has firmly established itself as one of the most informative investigation method especially because of its accuracy, it cannot be easily implemented in embedded systems. In this paper, broadband excitation signals, frequently used in system identification applications, are proposed to perform impedance measurements on a battery cell. Moreover, spectral coherence is an advanced parameter estimated in order to determine the frequency bands where the transfer function of the system is accurately identified. We propose in this study to test and compare the identification performances of such signals for the broadband monitoring of a battery.

I. INTRODUCTION

Batteries are an integral part and critical backup system of EV and HEV. Li-ion battery technology is believed to be the most attractive for these applications. It ensures higher power and energy densities, long cycle life, low cost of raw materials and superior safety characteristics [START_REF] Wakihara | Lithium ion batteries: fundamentals and performance[END_REF] [START_REF] Matsuki | General Concepts[END_REF]. In order to ensure safety in vehicles and improve the use of batteries, a Battery Management System (BMS) is involved [START_REF]Battery management systems for large Lithium-Ion battery packs[END_REF]. It should provide an efficient way to monitor battery performances and assessment of its condition in order to increase the reliability of EV and HEV systems. Several studies [START_REF] Montaru | Frequency and temporal identification of a Li-ion Polymer battery model using fractional impedance[END_REF][5] [START_REF] Kozlowski | Electrochemical Cell Prognostics using Online Impedance Measurements and Model-Based Data Fusion Techniques[END_REF] point out the usefulness of cell impedance measurements. For this reason, electrochemical impedance spectroscopy (EIS) [START_REF] Buller | Impedance-Based Simulation Models for Energy Storage Devices in Advanced Automotive Power Systems[END_REF] [START_REF]Frontmatter, in Impedance Spectroscopy: Theory, Experiment, and Applications[END_REF] is frequently used to better investigate the states of a battery. Electrochemical impedance spectroscopy can be performed either in a galvanostatic or in a potentiostatic mode. Following the first approach, a small AC current flows through the storage device under investigation and its AC voltage response is measured. The electrical impedance is estimated from the single frequency AC current superimposed a DC current (charging or discharging and used to define the overall working point of the cell) and the measured AC voltage response [START_REF] Carbonini | Degradation behaviour of 6013-T6, 2024-T3 alloys and pure aluminium in different aqueous media[END_REF][10] [START_REF] Uruchurtu | Electrochemical investigations of the activation mechanism of aluminium[END_REF]. Though its accuracy, an important drawback appeals research to seek for new ways of operating. Indeed, EIS is still a laboratory technique that cannot be easily implemented in embedded systems. Accordingly, broadband impedance identification using different types of excitation pattern is proposed. The concept consists in measuring the system response at multiple frequencies at the same time.

This paper focuses on the test of these broadband excitation signals and the comparison of their identification performances for the estimation of battery electrical impedance. After a brief review of nonparametric broadband identification basics for linear and time-invariant systems, the method is applied using simulated data. Spectral coherence is computed to select the frequency band where the system deals with LTI hypothesis. Finally, experimental results that validate the relevance of this approach are presented.

II. NON PARAMETRIC IDENTIFICATION THEORY

A linear and time invariant (LTI) single input single output (SI/SO) system H is completely characterized by its impulse response h n [ ] or its frequency response function H λ ( ), which are related by a Fourier transform:

H λ ( ) = h n [ ] n =-∞ +∞ ∑ e -j 2πλn (1)
In this equation, j = -1 and λ ∈ -

1 2 ; 1 2      
is the normalized frequency, leading to the frequency f in Hertz when multiplied by the sampling frequency.

Non-parametric identification of LTI systems aims to estimate the frequency response function H λ ( )

from input x n
[ ] and noisy output measurements z n [ ] (Figure 1) without the use of any model. Thanks to this last assumption, one can write:

S zx λ ( ) = S yx λ ( ) = H λ ( )S xx λ ( ). (2)
Therefore, on the frequency bands where the input PSD S xx λ ( ) ≠ 0 , the unknown frequency response function H λ ( ) can be calculated through:

H λ ( ) = S zx λ ( ) S xx λ ( ) if S xx λ ( ) ≠ 0 . (3)
This finally leads to the frequency domain identification of the unknown system H . Eq. ( 3) clearly shows that it is advantageous to use broadband input signals x n [ ] since they allow the computation of H λ ( ) on a wide frequency band as a whole.

An essential quantity in such a method is the spectral coherence between measured signals x n [ ] and z n

[ ] [10][11]:

c zx λ ( ) 2 = S zx λ ( ) 2 S xx λ ( )S zz λ ( ) . ( 4
)
This statistical quantity is bounded by 0 and 1, and measures the linear dependency or correlation between x n [ ] and z n [ ] at each frequency λ [12][13] and is used in the following to compute confidence limits for different spectral estimators.

The PSD and CPSD used in Eq. ( 3) and ( 4) can be easily estimated through the well-known Welch modified periodogram [START_REF] Bendat | Random data: analysis and measurements procedures[END_REF]. In this method, measured signals are split-up into L data segments of length N and as an example, the corresponding estimator of the CPSD between x n [ ] and z n [ ] is

given by:

ˆ S zx λ ( ) = A L Z k λ ( )X k * λ ( ) k = 0 L -1 ∑ , (5)
where:

-A is a normalization factor,

-Z k λ ( ) (resp. X k λ ( )) is the Fourier transform of the k th windowed segment of z n [ ] (resp. x n [ ]),
-* denotes the complex conjugate.

It yields that the spectral coherence c zx λ ( )

2 and the frequency response function H λ ( ) can be estimated by using Eq. ( 5) in Eq. ( 3) and ( 4):

ˆ c zx λ ( ) 2 = ˆ S zx λ ( ) 2 ˆ S xx λ ( ) ˆ S zz λ ( ) , (6)
ˆ

H λ ( ) = ˆ S zx λ ( ) ˆ S xx λ ( ) if ˆ S xx λ ( ) ≠ 0 . (7)
Halliday [START_REF] Halliday | A framework for the analysis of mixed time series/point process data-theory and application to the study of physiological tremor, single motor unit discharges and electromyograms[END_REF] includes the line

1 1 %) 95 1 ( 1 - - - L
on the coherence plot as an estimate of the upper 95% confidence limit under the hypothesis of independence. Estimated values of coherence lying below this line can be taken as evidence for the lack of a linear association between input and output. Halliday [START_REF] Halliday | A framework for the analysis of mixed time series/point process data-theory and application to the study of physiological tremor, single motor unit discharges and electromyograms[END_REF] also used the estimated coherence to compute upper and lower 95% confidence limits for the gain and

phase estimates ˆ G λ ( ) = ˆ H λ ( ) (Eq. (8)) and ˆ P λ ( ) = arg ˆ H λ ( )
{ } (Eq. ( 9)).

log 10 ˆ G λ ( ) { } ± 1.96 log 10 e ( ) ( ) 2 2L 1 -ˆ c zx λ ( ) 2 ˆ c zx λ ( ) 2 (8)
ˆ P λ ( ) ± 1.96 1 2L 1 -ˆ c zx λ ( ) 2 ˆ c zx λ ( ) 2 (9)
Finally, Eq. ( 5) and ( 7) constitute the "identification algorithm" appearing in Figure 1 and used to estimate the frequency response function H λ ( ) of an unknown LTI system through its input x n [ ] and noisy output z n [ ]. Eq. ( 6), ( 8) and ( 9) are used to evaluate the performance of the algorithm by computing 95% confidence limits of the previous estimators.

In what follows, batteries are modelled as electrical LTI systems and are stimulated in galvanostatic mode (the current conrresponds to the input and the voltage response to the output). The corresponding frequency response function is then the electrical impedance of the battery, and this quantity is estimated thanks to the previous set of equations.

III. BROADBAND IDENTIFICATION OF ELECTRICAL IMPEDANCE FOR LI-ION BATTERIES

A. Battery modelling and simulation

Equivalent electrical circuits (EEC) are commonly used to reproduce the battery electrical behavior.

They consist on passive (resistors, capacitors, inductors, constant phase elements) and active (voltage and current sources) elements [START_REF]Modélisation électrique et énergétique des accumulateurs Lithium-Ion[END_REF]. For electrical engineers, such models are able to characterize electrochemical phenomena, and lead to perform a quick analysis and prediction of the battery behavior in frequency as in time domains [START_REF] Dong | Dynamic Modeling of Li-Ion Batteries Using an Equivalent Electrical Circuit[END_REF]. Adapted Randles model (Figure 2) is the EEC adopted in our study. It was developed, used and implemented using Simulink in [START_REF] Halliday | A framework for the analysis of mixed time series/point process data-theory and application to the study of physiological tremor, single motor unit discharges and electromyograms[END_REF]. It includes the modelling of theoretical impedance is considered to be the analytical expression of the EEC impedance. In order that the battery meets the requirements of a LTI system described in section II, the whole measurement time should be chosen so that only a little variation of the battery SOC occurs during the identification step.

We focused on the frequency band from 1 to 7 Hz. The simulations are undertaken under the same operating conditions of SOC (60%) and DC current (1A). We consider a time duration of 90 seconds, and the sampling frequency s f is 700 Hz. These signals are cut into L=10 disjoint segments of length N=6300 samples. A white Gaussian noise with zero mean and a variance such that the signal to noise ratio (SNR) is 0 dB is added on the output voltage signal to simulate measurement noise. Spectral coherences and electrical impedances are then estimated from this noisy voltage and the excitation signals through Eq. ( 6) and [START_REF] Buller | Impedance-Based Simulation Models for Energy Storage Devices in Advanced Automotive Power Systems[END_REF]. Estimated electrical impedances are finally compared to the theoretical impedance value in order to evaluate the quality of the broaband identification process.

B. Broadband excitation signals

Broadband identification is valid when the excitation signal shows an almost flat power spectral density in the frequency band of interest. Pseudo random binary sequences (PRBS), swept square (square pattern with a period varying continuously) and a square wave (square pattern with constant period) are broadband signals frequently used in system identification applications [START_REF] Pintelon | System identification, a frequency domain approach[END_REF]. The latter signal is not on itself a broadband signal but can be considered so if its harmonics are in the frequency band of interest. Such signals allow the estimation of the frequency response function of LTI systems, in particular the battery impedance, over a large bandwidth from a single set of measurements. In this work, these broadband signals are introduced as identification patterns that can be used in embedded systems.

Estimate PSDs are computed using Eq. ( 5) and the same values as those given in section III.A are chosen (L=10, N=6300,

C. Coherence results

We consider a reasonable measurement noise of SNR=0 dB.

D. Confidence limits results

Confidence limits upon the gain and the phase of the estimated impedance quantify the estimation performance reached by this identification method and can be easily computed by using Eq. ( 8) and ( 9).

In this section, the results concerning the gain factor for PRBS, swept square and square are given as an illustrative example, and they are exclusively represented in the selected frequency band. Phase factor results are very similar.

The different figures of Figure 5 show that the PRBS (Figure 5(a)) has large confidence limits near the upper limit of the band max f , which confirms the coherence results of Figure 4. Moreover, tight confidence limits are observed upon the whole selected frequency band using a swept square (Figure 5(b)). The square wave (Figure 5(c)) provides, as expected, tight confidence limits only around its odd harmonics frequencies. We infer from those results that broadband impedance can be accurately identified with signals composed of square patterns. Such signals are easy to apply to a battery from simple electronic components, for example by using electronic switches.

IV. EXPERIMENTAL RESULTS

A. Hardware and implementation

This work was realized on a graphite/LiFePO4 cell with a nominal capacity of 2.3 Ah (ANR26650m1 battery from A123 Systems Company Ltd). Experiments have been carried out at room temperature. As a first step, we consider the same operating conditions that those taken in simulations (SOC of 60% and DC current of 1A, number of blocks L=10, T=9s). We also consider the same frequency band [1Hz -7Hz] for different SOCs. The goal is to verify if broadband impedance identification can recover the known distortions Nyquist plots with the evolution of the battery SOC. A specific electronic circuit was designed to perform the experiments and to allow application of input current with squared patterns, particularly swept square, PRBS and square.

B. Experimental results

1)

PRBS, swept square and square for a frequency band of 1 -7Hz

The signal as predicted This shows that the battery can be considered as a LTI system under the operating conditions used, and that its electrical impedance will be correctly estimated. As shown by simulations of section III. it can also be noticed that the use of a PRBS current induces a decrease in the coherence near the upper limit frequency max f . Accordingly, the variability of the estimated impedance with the different signals is very small all over the frequency band, unless near the upper frequency for PRBS. Based on the high coherence values obtained during experiments, this battery can be considered as a LTI system within the chosen frequency band and under the chosen operating conditions. Moreover, confidence limits of impedance estimators are sufficiently small to affirm that the electrical impedance is accurately identified all over the frequency band, and that the swept square input current leads to better results than the PRBS.

2)

PRBS for a frequency band of 1 -7Hz at different SOCs Several studies reveal that the SOC distorts the battery electrical impedance at lowest frequencies. In this context, we design a PRBS that can excite the frequency band of [1 -7] Hz. The choice of PRBS is based on the previous results that showed the strong performances of this excitation signal near the lower limit frequency band. Exciting low frequencies yields that the number of blocks L used for averaging cannot be higher than done previously (because of LTI assumptions). In that case it will be limited to L=10 blocks of duration 9 seconds each. The sampling frequency is unchanged. The battery is stimulated under several SOC (45, 60, 90 %) and for each SOC a broadband identification step is performed.

Figure 7 shows that the PRBS is able to reflect the SOC effect on the battery electrical impedance. With a small time measurement (90 s), this method is thus able to accurately estimate the battery impedance for low frequencies. We can ensure the quality of the estimation thanks to the coherence values which are once again very close to 1 .

Figure 1 :

 1 Figure 1: Non-parametric identification of a LTI system H in the frequency domain. The unknown additive measurement noise b n [ ] is supposed uncorrelated with x n[ ] and therefore with

  connectors and electrolyte (R,L), passivation film (Rf // CPEf), charge transfer (Rtc ) and double layer phenomena (CPEdl). The open circuit voltage (OCV) is given in a look up table with respect to the current intensity and the battery state of charge. This circuit introduces constant phase elements (CPE) to accurately reflect the behaviour of the battery observed on impedance spectroscopy measurements [17][18].

Figure 2 .

 2 Figure 2. Equivalent electrical circuit (adapted Randles model) of a graphite/LiFePO4 battery Circuit parameters are function of the state of charge (SOC) and the current intensity. In simulations,

sf

  =700Hz, T=9s). The resulting PSDs are shown in Figure 3. Under assumption of constant power for the excitation signals, we note that the spectral information is different in levels and frequency bandwidths.

Figure 3 .

 3 Figure 3. PSD of the three excitation signals.

  Figure 4 shows the coherences estimated with the three excitation signals previously defined, and it can be clearly noticed that the results are coherent with the PSD. Higher coherence values are obtained within the considered frequency band.

Figure 4 .

 4 Figure 4. Coherence plots for the three excitation signals with a measurement SNR=0 dB.

Figure 5 .

 5 Figure 5. 95% confidence limits results using: (a) PRBS, (b) swept square, (c) square as excitation signals with a measurement SNR=0 dB.

  broadband signals based on squared patterns are applied to the battery. The corresponding estimated coherences are plotted in Figure 6 (d), (e) and (f) , while Figure 6 (a), (b) and (c) show the (b) (c) estimated electrical impedance. The coherence is clearly close to 1 all over the specified frequency band for PRBS and swept square signals while it has only strong values at odd harmonics of the square

Figure 6 .

 6 Figure 6. Experimental results: Nyquist plots -(a) PRBS, (b) square, (c) swept square coherence plots -(d) PRBS, (e) square,(f) swept square.

Figure 7 .

 7 Figure 7. Experimental results with PRBS signal and under several SOCs:-(a) Nyquist plots, (b) coherence plots.V. CONCLUSIONThis paper focuses on the usefulness of broadband excitation signals for the identification of a Li-ion battery electrical impedance. After a review of non-parametric identification theory and advanced parameters such as coherence function and confidence limits, the simulation results ensure that this method can be applied to battery systems. Experimental tests performed at a low frequency band reveal that signals based on square pattern like swept square, PRBS and square lead to correct broadband identification results. They are particularly well suited for electronic implementation. Experimental results also obtained at a low frequency band for different SOCs confirm the possibility to apply PRBS to batteries to monitor their SOC evolution. These results are promising to improve and enhance actual BMS.

The authors thank Marco Ranieri for its contribution to this work.


Hascoët, L Claire Lauvernet
email: claire.lauvernet@irstea.fr

Le Dimet

Laurent Hascoët
email: laurent.hascoet@inria.fr

Frédéric Baret
email: baret@avignon.inra.fr

Using Automatic Differentiation to study the sensitivity of a crop model

Keywords: adjoint mode, agronomic crop model, sensitivity analysis 1

Automatic Differentiation (AD) is often applied to codes that solve partial differential equations, e.g. in geophysical sciences or Computational Fluid Dynamics. In agronomy, the differentiation of crop models has never been performed since these models are more empirical than fully mecanistic, derived from equations. This study shows the feasability of constructing the adjoint model of a crop model referent in the agronomic community (STICS) with the TAPENADE tool, and the use of this accurate adjoint to perform some sensitivity analysis. This paper reports on the experience from AD users of the environmental domain, in which AD usage is not very widespread.

terms of quantity and quality, environment in terms of drainage and nitrate leaching, and to soil characteristics evolution under cropping system 1 .

The two key output variables simulated by STICS that we will need in this paper are the Leaf Area Index (LAI) and the biomass. The LAI is the total one-sided area of leaf tissue per area of ground surface (unitless). This is a canopy parameter that directly quantifies green vegetation biomass. As the leaves are considered to be the main interfaces with the atmosphere for the transfer of mass and energy [START_REF] Rosenberg | Microclimate: the biological environment[END_REF], the LAI indirectly describes properties such as potential of photosynthesis available for primary production, plant respiration, evapotranspiration and carbon flux between the biosphere and the atmosphere, and gives evidence of severely affected areas (fires, parasites. . . ). Because it is the most observable canopy parameter by remote sensing, the LAI is very commonly used e.g., in crop performance prediction [START_REF] Guérif | Data assimilation and parameter estimation for precision agriculture using the crop model STICS[END_REF], in models of soil-vegetation-atmosphere [START_REF] Olioso | Future directions for advanced evapotranspiration modeling: Assimilation of remote sensing data into crop simulation models and SVAT models[END_REF], in crop models [START_REF] Brisson | Author-produced version of the article[END_REF][START_REF] Brisson | STICS: a generic model for simulating crops and their water and nitrogen balances. II: model validation for wheat and maize[END_REF], in radiative transfer models [START_REF] Verhoef | Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model[END_REF]. Its values can range from 0 for bare soil to 6-7 for a crop during its life cycle, and up to 15 in extreme cases (tropical forests). STICS simulates the crop growth from sowing to harvest, focusing on the evolution of the LAI at a few selected [START_REF] Brisson | Author-produced version of the article[END_REF] vegetative stages shown on Fig. 1. These stages involve process thresholds, accounting for some of the differentiation problems described in Sect. 3.2. For a wheat crop, the main phenological stages are known as ear at 1 cm, heading, flowering, and maturity. In this work we do not simulate grain yield but only the total biomass. As we focus on the LAI, we only consider the vegetative stages namely: LEV (emergence or budding), AMF (maximum acceleration of leaf area index, equivalent to ear at 1 cm), LAX (maximum LAI i.e. end of leaf growth), and SEN (start of net senescence).

Sensitivity Analysis

A model is a more or less realistic or biased simplification of the state variables it simulates. This is especially true for agronomic models, since the functioning of vegetation is not a priori described by exact equations: agronomic models attempt to predict the behavior of the crop by incremental improvements of the simulation code, based on observations made on the field and then published by specialists. Thus, in some parts of the model, this empirical approach is not based on the equations of some underlying physics or chemistry. Sensitivity analysis, which studies the impact of perturbing the control parameters on the model output, gives insights useful to improve or even simplify the model. Sensitivity analysis requires two essential ingredients:

• A model: F(X, K) = 0, where X is the state variable (LAI, biomass . . . ) and K the control variables (parameters, forcing variables . . . ). F is a differential operator a priori non-linear finite-dimensional, that describes implicitly the evolution of X for a given K. We assume that the system has a unique solution X(K). In this study, what we call the model is exactly the STICS computer program. • A response function G which combines one or more elements of X into a scalar value, e.g. the final value or the integral over time of an output.

The problem is to evaluate the sensitivity of G with respect to K or in other words the gradient of G with respect to K. With the help of the adjoint model, computing the gradient takes only 2 steps: run the direct model once for the given K, then solve the adjoint model once [START_REF] Dimet | Sensitivity analysis in variational data assimilation[END_REF]. The classical justification is:

∇G = dG dK t = dG dX . dX dK t = dX dK t .
dG dX t where we observe that dG dX is easily computed from the definition of G alone and the product of dX dK t with a vector is achieved by feeding this vector to the adjoint code of STICS, produced by the adjoint mode of Automatic Differentiation.

Sensitivity analysis using an adjoint model is the only way to calculate formally the gradient of the response function at a cost that does not depend on the size of K. It is particularly suitable when the number of entries K is large compared to the size of the response function G [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF].

One can also compute the gradient accurately with tangent-linear differentiation, at a cost that is proportional to the size of K. The other sensitivity methods only approximate the gradient: finite difference approximation of the gradient require extensive direct model computations [START_REF] Castaings | Sensitivity analysis and parameter estimation for distributed hydrological modeling: potential of variational methods[END_REF]. Stochastic sampling techniques require less mathematical insight as they consist (roughly speaking) in exploring the space of control to determine an overall global sensitivity [START_REF] Saltelli | Sensitivity Analysis[END_REF][START_REF] Ionescu-Bujor | A comparative review of sensitivity and uncertainty analysis of large-scale systems. I: deterministic methods[END_REF]. Their cost grows rapidly with the dimension of K. These methods have been widely applied to the agronomic models and in particular on STICS [START_REF] Houlès | Evaluation of the ability of the crop model stics to recommend nitrogen fertilisation rates according to agro-environmental criteria[END_REF][START_REF] Ruget | Sensitivity analysis of a crop simulation model, STICS, in order to choose the main parameters to be estimated[END_REF][START_REF] Varella | Global sensitivity analysis measures the quality of parameter estimation: The case of soil parameters and a crop model[END_REF].

If in many cases, the response function G is a differentiable function of K, it can happen that the model is driven by thresholds e.g., the code uses a lot of branches. Theoretically, a piecewise continuous function is not continuously differentiable, but it has right-and left-derivatives. Differentiation of such a code can only return a sub-gradient. Actually, the methods that do not rely on derivatives (divided differences, stochastic,. . . ) behave better in these cases, although they remain expensive. In practice, this problem is not considered serious as long as the local sensitivity is valid in a neighborhood of the current K.

Automatic Differentiation of STICS

The TAPENADE Automatic Differentiaton tool

TAPENADE [START_REF] Hascoët | TAPENADE 2.1 user's guide[END_REF] is an Automatic Differentiaton (AD) tool based on source transformation. Given a source program written in FORTRAN, TAPENADE builds a new source program that computes some of its derivatives. In "tangent" mode, TAPE-NADE builds the program that computes directional derivatives. In "adjoint" mode, TAPENADE builds the program that computes the gradient of the output with respect to all input parameters.

Considering the complete set of derivatives of each output with respect to each input, i.e. the Jacobian matrix of the program's function, the tangent mode yields a column of the Jacobian whereas the adjoint mode yields a row of the Jacobian. Therefore in our particular case where the output is a scalar G, one run of the adjoint code will return the complete gradient. In contrast, it takes one run of the tangent mode per input to obtain the same gradient. Although we will experiment with the two modes, the adjoint mode fits our needs better.

However, the adjoint mode evaluates the derivatives in the inverse of the original program's execution order. This is a major difficulty for large programs such as STICS. The AD model copes with this difficulty by a combination of storage of intermediate values and duplicated evaluation of the original program, at a cost in memory and execution time. In TAPENADE, the strategy is mostly based on storage of intermediate values, combined with the storage/recompute tradeoff known as checkpointing, applied automatically at each procedure call.

STICS adjoint : the pains and sufferings of an AD end-user

The STICS model being written in FORTRAN 77, TAPENADE can in theory build its adjoint. However, there were shortcomings with the early versions of TAPE-NADE, before 2005. Later versions brought notable improvements but we believe it is worth describing the main problems that we encountered at these early stages.

AD allows for instructions which the symbolic differentiation systems cannot process. It also provides a real gain in computational time. However, a few good programming practices are recommended: the input parameters involved in derivatives must be clearly identified and if possible separate from the other variables. The same holds for the outputs to be differentiated. The precision level of all floating point variables must be coherent, especially for validation purposes: if the chain of computation is not completely "double precision", then the divided difference that is used to validate the analytic derivatives will have a poor accuracy, validation will be dubious and may even fail to detect small errors in the differentiated code.

Validation helped us detect small portability problems in STICS. As divided differences requires to call STICS twice, we discovered that two successive calls to STICS apparently with the same inputs gave different results. In fact the first call was different from all the others, which pointed us to a problem of hidden uninitialized remanent global. Fixing this gave us correct divided differences, and a more portable STICS code.

More specifically to this agronomy application, we had problems with the high number of tests and other conditional jumps in an average run. In more classical situations of Scientific Computing, programs are derived from mathematical equations, typically a set of ODE's or PDE's. This forces some regularity into the code that discretizes and solves these equations: even if branches do occur, they rarely introduce discontinuity and the derivative itself often remains continuous. In our application, the program itself basically is the equation. The model evolves by introducing by hand new subcases and subdivisions, i.e. more tests. If this evolution is not made with differentiation in mind, it may introduce sharp discontinuities that do not harm the original code but make it non-differentiable. It took time to replace faulty branches with a cleaner, differentiable implementation. On the other hand, users agreed that this resulted in a better code.

Still, the number of branches in the STICS model is very large: thresholds, conditions, loops, and other control, all are tests that the adjoint code must remember to run backwards. STICS consumes an unusually large memory for that. Until recently, TAPENADE did not store this control efficiently, using a full INTEGER value to store only a boolean in general.

Checkpointing the time stepping was difficult. Before binomial checkpointing [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation[END_REF] was implemented in TAPENADE, we had to split the main time loop of 400 iterations into two nested loops of 20 iterations each, and place these two loops into two new subroutines to force checkpointing. This tedious manipulations are now spared with the new TAPENADE directives for binomial checkpointing.

More than five years after this sensitivity study, both STICS and TAPENADE have evolved. The latest version 6 of STICS is more readily differentiable than before. TAPENADE 3.6 had several bugs fixed and more importantly provides a set of user directives to control checkpointing better. These checkpointing directives are also the answer to the serious performance problem discussed in Sect. 3.3.

Validation of the adjoint model

Validation was performed in two steps as usual, and for several directions of perturbation. First, the tangent derivatives were compared with divided differences, and they agreed up to the 8 th decimal for an increment of 10 -8 in the one-sided divided difference. Second, the adjoint derivatives were compared with the tangent derivatives ("dot-product" test [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF]) and they agreed up to the 14 th decimal. At the time of the study, the run times were: Direct model : 0.21s Tangent model : 0.39s Adjoint model : 30.96s

The run time of the adjoint code is much higher than the custommary 5-to 10-fold slowdown. The problem was left to the TAPENADE developers to go on with the sensitivity study. Identifying its causes was hard, and pointed to the need for specific profiling tools for adjoint codes. Profiling instructions must be inserted by the AD tool itself, and tools are missing to help interpret the profiling results. Eventually, the problem was found to come from the systematic checkpointing on procedure calls on a chain of 4 nested procedure calls, each of them doing little else than calling the next nested call, cf Fig. 2. Checkpointing [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF] one call to P reduces the peak memory used by the adjoint. This reduction is roughly proportional to the run-time of P. On the other hand, it costs one extra run of P, plus some memory (a "snapshot") to restore the execution state. Checkpointing nested calls causes an increasing number of extra runs. This is inherent to the approach and beneficial in general, but is a waste for procedures that are nearly empty shells around a deeper call. In our case, the problem was amplified by the size of a very big work array that was restored at each checkpoint. The answer is to deactivate checkpointing on the calls to the "empty shell" procedures. This is known as the "split" mode of adjoint AD [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF], and is sketched on the right of Fig. 2. This required development in TAPENADE, plus new directives ($AD NOCHECKPOINT) to let the user trigger this split mode on selected procedure calls. Conversely in other cases, it is useful to trigger checkpointing on pieces of a procedure, and TAPENADE new directives ($AD NOCHECKPOINT) and ($AD NOCHECKPOINT) let the user do that. This results in the following times obtained with TAPENADE 3.6:

Direct model : 0.22s Tangent model : 0.52s Adjoint model : 0.86s

Results : sensitivity analysis of STICS

We decided to compute the gradients of two response functions G: LAI and biomass, and more precisely their integrals over the simulation time from sowing to harvest. These response functions capture well the growth dynamics.

Author-produced version of the article published in Springer's Lecture Notes in Computational Science and Engineering, vol. Recent Advances in Algorithmic Differentiation, vol. 87, p. 59-69

The original publication is available at http://www.springer.com/series/3527

G LAI = T ∑ i=1 LAI(t i ) G biomass = T ∑ i=1 biomass(t i )

Selection of input parameters for sensitivity analysis of output variables

For this feasibility study, the control variables correspond to wheat crops from the Danube's plain in Romania in 2000-20012 [START_REF] Baret | The adam database and its potential to investigate high temporal sampling acquisition at high spatial resolution for the monitoring of agricultural crops[END_REF]. The gradient was calculated with respect to the following input parameters3 : for LAI, we chose the varietal parameters acting on the dynamics of LAI, and dlaimaxbrut that strongly characterizes the aerial growth. Parameters were adapted to the ADAM database, including the variety of wheat (Flamura) used here for its particular cold resistance. For biomass, efficiencies at three important phases of the cycle of wheat (juvenile phases, vegetative and grain filling) and vmax2 were chosen following the experience accumulated by users of the crop model. Table 1 describes the role of these parameters, and their values for this sensitivity study.

Sensitivity results of LAI and biomass

One goal of this sensitivity study was to establish the hierarchy of influent parameters. Therefore Fig. 3 shows the 10 influences normalized as percentages, totalling 100%. Among the 10 selected, the most influential parameters on the LAI are adens (47%), dlaimaxbrut (21%), stlevamf (17%), jvc (10%), and finally stamflax (2%). adens represents the ability of a plant to withstand increasing densities, and since it depends on the species and varieties, its influence may be particularly strong for this type of wheat and less for other crops. For biomass, we observe that the hierarchy is modified by the strong influence of the efficiency efcroiveg (maximum growth efficiency during vegetative phase) which is similar to that of adens (27%). This means that we can ignore the estimate of efcroiveg if we only want to assimilate LAI data, but absolutely not if we need to simulate biomass. stlevmaf and dlaimaxbrut are of similar importance (14 and 12%). Finally, there is a relatively low sensitivity (5% and 3%) of biomass integrated over the life cycle to the other two parameters of efficiency efcroirepro and efcroijuv, meaning that the biomass is not so dependant on the juvenile and the grain filling phases but essentially on the vegetative phase.

The fact that only the integral over the entire cycle was studied involves a very small influence of the parameters efcroirepro and efcroijuv, as opposed to efcroiveg. These efficiencies with a small influence matter only during short phenological stages: only a sensitivity study restricted to these stages can modify the hierarchy of influent parameters, opening the way to estimation of these low-influence parameters [START_REF] Ruget | Sensitivity analysis of a crop simulation model, STICS, in order to choose the main parameters to be estimated[END_REF]. LAI is actually dependant on 4 parameters and biomass on 5 on the 10 tested, which will help the user concentrate on these and estimate them better. Uncertainty on the other parameters is of relatively smaller importance.

Conclusion and outlook

This case study illustrates the interest of AD for sensitivity analysis of agronomic models. Coupled with other models, for example radiative transfer model [START_REF] Lauvernet | Multitemporal-patch ensemble inversion of coupled surface-atmosphere radiative transfer models for land surface characterization[END_REF], it will allow to assimilate remote sensing data into crop models by using the adjoint to minimize the discrepancy cost function. This work shows the feasability of applying and developing variationnal methods in agronomy, in the same way as in oceanography or meteorology. Author-produced version of the article published in Springer's Lecture Notes in Computational Science and Engineering, vol. Recent Advances in Algorithmic Differentiation, vol. 87, p. 59-69

The original publication is available at http://www.springer.com/series/3527

For the agronomic community, the adjoint model of STICS is an interesting tool to perform sensitivity analysis since it requires the calculation only once for each agro-pedo-climatic situation. The most difficult work is the differentiation of the model, which must be done only once, and with the help of AD tools that keep improving. However, the local sensitivity analysis is valid only in a small neighborhood and the hierarchy of sensitivities may vary under different conditions. These results are only a first step. Following work could concentrate on:

1. a "multi-local" sensitivity analysis, keeping the crop management and climate of the base ADAM, but letting the parameters vary in a given range. This would require many runs of the adjoint modes on a representative sample of possible parameter values. This would return a parameter hierarchy with a more general validity. 2. an application of this analysis to other conditions (climate, soil...) to see whether the hierarchy is preserved in general. Extending to other varieties is also important. Actually, it seems unlikely that this hierarchy is preserved since the change of climate and soil conditions may rapidly hit limiting factors (stress for the plant) and thus modify the parameters influence. 3. a study of the sensitivity at selected phenological stages of the cycle to study the effect of variables temporally valid (especially efficiency) on the general hierarchy.

The adjoint code is able to compute the sensitivities of one response function to all parameters in just one run. There are more parameters to STICS than the 10 we have selected for this sensitivity study. Looking at the influence of all parameters will guide the attention of STICS users on some parameters and modules, according to the users' objectives. Sensitivity study is a preliminary to parameter estimation: many of these agronomic parameters (yield, balance . . . ) are not directly observable by remote sensing. On the other hand the outputs (biomass) can be measured. The adjoint of the model, by returning the gradient of any discrepancy cost function, is the key to estimate these hidden agronomic parameters from the ones we can measure.

Fig. 1

 1 Fig. 1 Simplistic scheme of the stages simulated by the STICS model on dynamics of LAI.

Fig. 2

 2 Fig.2The cost of checkpointing long chains of nested calls

Fig. 3

 3 Fig.[START_REF] Brisson | STICS: a generic model for simulating crops and their water and nitrogen balances. II: model validation for wheat and maize[END_REF] Relative sensitivity (%) to selected STICS parameters of ouput variables LAI (left) and biomass (right) computed by the adjoint.

Table 1

 1 Parameter role and values for the ADAM conditions

	parameter definition	value
	dlaimaxbrut maximum rate of gross leaf surface area production	0.00044
	stlevamf	cumulated development units between the LEV and AMF stages 208.298
	stamflax	cumulated development units between AMF and LAX stages	181.688
	jvc	days of vernalisation (cold days needed to lift)	35
	durvieF	lifespan of a cm of adult leaf	160
	adens	compensation between number of stems and plants density	-0.6
	efcroijuv	maximum growth efficiency during juvenile phase (LEV-AMF)	2.2
	efcroiveg	maximum growth efficiency during vegetative phase (AMF-DRP) 2.2
	efcroirepro maximum growth effiicency during grain filling phase (DRP-MAT) 4.25
	vmax2	maximum rate of nitrate absorption by the roots	0.05

ADAM experiment (Data Assimilation through Agro-Modelling). Project and database at http://kalideos.cnes.fr/spip.php?article68

All the parameters of STICS are described in http://www.avignon.inra.fr/agroclim stics eng/notices d utilisation

Acknowledgements

This study was conducted thanks to a grant provided by CNES within the ADAM project (http://kalideos.cnes.fr/spip.php?article68), during the Phd of the first author at INRA Avignon and the University of Grenoble.