Fatima Zahra

Michelle Schatzman

Computation of Lyapunov Exponents for Dynamical System with Impact

Keywords: Mathematics Subject Classification: Dynamical Systems Dynamical Systems, Impact problem, Lyapunov Exponents

Consider the dynamical system u + 2o:u + u = acoswt where the position u is constrained to remain above an obstacle of height Umin; when u reaches the obstacle, its velocity is reversed and multiplied by a restitution coefficient e E [0, 1]. For certain choices of parameters, the solutions are chaotic. We compute the Lyapunov exponents by three different methods, and we compare the results. The computation of these numbers is very sensitive to the method, and to the numerical parameters for a given method, even with a very accurate method.

Introduction

In this article, we consider a class of non smooth dynamical systems, which describe the motion of a mechanical system with one degree of freedom, subject to a unilateral constraint on the position. vVhen this constraint is saturated, the velocity is reversed and multiplied by a restitution coefficient.

More precisely, let f be a continuous function from IR 3 to IR, which is Lipschitz continuous with respect to its last two arguments. Assume that a real number Umin and a restitution coefficient e are given. The initial data are uo > Umin and Vo, such that Vo > 0 if uo = Umin• The solutions of the Cauchy problem, are defined as follows: u is a continuous function from [to, T] whose second derivative (in the sens of distributions) is a measure. This condition implies that the first derivative u of u is a function of bounded variation, and in particular, it is continuous almost everywhere, and it is continuous from the left and from the right at all points of (to, T] and [to, T) respectively. Thus, the difference uj(t, u, u) will be a measure 11; the function u and the measure 11 have to satisfy the following relations:

u=f(t,u,u)+/1 U > Umin suppf1 C {t: u(t) = Umin} if u(t) = Umin, then u(t + 0) = -eu(t-0) u(to) = Uo, u(to) = Vo (1) (2) (3) (4) (5)
A few words of comment on this definition are in order: the first derivative of u is expected to be discontinuous at impacts, for obvious geometric reasons; therefore, it makes sense to assume that the second derivative of u has Dirac masses, or more generally is a measure. The measure 11 is the reaction of the obstacle at impact, as can be seen from (1), which is basically Newton's law. Condition (2) means that u(t) remains inside the convex of constraints [umin, oo ). The reaction can be different from 0 only when there is a contact: this is relation [START_REF] Golub | Matrix Computations[END_REF]. Finally, condition (5) describes the constitutive law of the impact, with the help of the restitution coefficient e.

It has been proved in [START_REF] Paoli | Mouvement a nombre fini de degres de liberte avec contraintes unilaterales: cas avec perte d'energie[END_REF] that problem (1)-( 5) possesses a solution. However, uniqueness is not always true and an example has been given in [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF]. Nevertheless, generic uniqueness has been proved in a special case [START_REF] Carriero | Uniqueness of the one-dimensional bounce problem as a generic property in L 1 ([0, T]; IR)[END_REF] for n degrees of freedom; uniqueness has also been proved for one degree of freedom when f is analytic with respect to all its arguments [START_REF] Schatzman | Uniqueness and continuous dependence on data for onedimensional impact problems[END_REF].

A simple case of (1) -( 4) is a forced vibrating system with one degree of freedom defined by j(t, u, v) =a cos(wt) -2av-u. [START_REF] Paoli | A numerical scheme for a dynamical impact problem with loss of energy[END_REF] In this case, the uniqueness theorem of [START_REF] Schatzman | Uniqueness and continuous dependence on data for onedimensional impact problems[END_REF] applies.

For particular values of the parameters a, Umin and e, the problem (1 )- [START_REF] Paoli | Mouvement a nombre fini de degres de liberte avec contraintes unilaterales: cas avec perte d'energie[END_REF] with the choice (6) of function f shows typically non linear phenomena of sensitivity to initial data, and there are values of the parameter for which there exists a "strange" attractor.

In this article, we compute Lyapunov exponents for the system (1)-( 6), using two different numerical methods.

The first numerical method is the impact detection method, and it is very straightforward: suppose we are given initial data u(tk) = uk and u(tk+O) = vk at an impact time tk. vVe use the elementary and formulre for the free flight solution, and we seek the first time tk+l > tk for which u(tk+I)-Umin vanishes. vVe perform this search by sweeping forward in time so as to find a small interval where u-Umin changes sign; on this small interval, we find the impact time by Newton's method. The accuracy is limited only by the capability of the computer. At instant tk+l, we reverse the velocity according to rule (4), and we start the process again using the new values of the initial data.

The second numerical method is the [START_REF] Paoli | Schema N umerique pour un modele de vibrations avec contraintes unilaterales et perte d'energie aux impacts[END_REF][START_REF] Paoli | A numerical scheme for a dynamical impact problem with loss of energy[END_REF] numerical scheme defined in [START_REF] Paoli | Schema N umerique pour un modele de vibrations avec contraintes unilaterales et perte d'energie aux impacts[END_REF][START_REF] Paoli | A numerical scheme for a dynamical impact problem with loss of energy[END_REF].

Let us describe this scheme: denote PK the projection on the convex set K =

[(1 + e)umin, +oo); it is given by

PK(x) = max(x, (1 + e)umin), (7)
and a sequence Fn, which is defined by: (8)

The numerical scheme is given by the following relation:

It turns out that the computations using (7)-( 9) are much faster than the computations by the impact detection method. However, the scheme is not very accurate; it is not better than first order with respect to the position. It has been observed that in a chaotic case, the attractor of ( 1)-( 6) calculated by the impact detection method is well approximated by the attractor calculated by the numerical scheme (7)- [START_REF] Pollicott | Lectures on Ergodic Theory and Pesin Theory on Compact Manifolds[END_REF]. vVe would like to estimate how well the numerical scheme (7)-( 9) approximates more refined information, namely the Lyapunov exponents of (1)-( 6):

we reduce our continuous time dynamical system in three dimensions to a discrete time dynamical system, by using a Poincare map; this is easy because we cut the phase space by the planes t = t 0 + kT, where T = 21r jw is the period of the forcing. In the two dimensional case, the largest Lyapunov exponent, if it exists, describes the average rate of divergence of two infinitesimally close trajectories; the sum of the two Lyapunov exponents, if it exists, describes the average rate of evolution of infinetesimal volume in phase space.

There are classical methods for calculating numerically Lyapunov exponents; they are described for instance in [START_REF] Parker | Pratical Numerical Algorithms for Chaotic Systems[END_REF], and consist essentially in giving a reasonably stable numerical implementation of the definition of the Lyapunov exponents. The calculations of Lyapunov exponents are very close to the approximation of eigenvalues of a matrix by the power method, and it has been known for a long time that renormalizations are needed to make the computation possible.

The fact that our system is not smooth is not a serious problem, because the system can be linearized whenever the impact takes place with a strictly positive velocity, and we never observed tangential impacts for choices of the parameters leading to a chaotic behavior.

The article is organized as follows: in section 2, we recall the definition of the Lyapunov exponents, and we define a program which could lead to a theoretical treatment of the qualitative questions we consider, and we describe its state of advancement. In section 3, we differentiate the flow with respect to the initial data: we consider separately the flow in continuous time and the flow in discrete time, and we show that a crucial term disappears in the discrete time case. In section 4, we explain the three methods of computation of the Lyapunov exponents. In section 5, we describe and compare our numerical results. In section 6, we conclude.

Definition of Lyapunov exponents

Let us recall the definition of Lyapunov exponents in a smooth case: let snx be a trajectory of a discrete dynamical system; if 8 is very small, the largest Lyapunov exponent is the asymptotic rate of evolution of ISn(x + 8)-sn(x)l as time n tends to infinity, if it exists; similarly if 8 1 and 8 2 are independant, the sum of the first two Lyapunov exponents describes the asymptotic rate of evolution of the area of the paralellogram built on the vectors sn(x+8!)-Sn(x) and sn(x + 8 2 ) -sn(x). That such objects can be defined and exist is a non trivial fact, which is proved only in rather particular cases. Let us state just a few of the results which would enable one to prove that the Lyapunov exponents exist. Our sources for the following description have been [START_REF] Pollicott | Lectures on Ergodic Theory and Pesin Theory on Compact Manifolds[END_REF], [START_REF] Ruelle | Ergodic theory of differentiable dynamical systems[END_REF], [START_REF] Ruelle | Chaotic Evolution and Strange Attractors[END_REF] and [START_REF] Ya | Topics in Ergodic Theory[END_REF]. vVe give first some definitions which are classical in ergodic theory: consider a measurable space, i.e. a pair (X, x) where X is an abstract set, and x is a o--algebra of its subsets. A mapping S from X to itself is an endomorphism of (X, x) if for all C belonging to x, the set r-1 C belongs also to x. A measure 11 on (X, x) is invariant under the endomorphism S iff for all C E x,

J1(S-1 C) = J1(C).
The celebrated Birkhoff-Khinchin ergodic theorem states that if 11 is invariant under S, then for almost every x and for all f E L 1 (X, x, 11) the following limit

/(x) = lim 1 t f(SJx) n-+oo n + 1 j=O exists.
Moreover, the limiting function f is invariant, i.e;

/(Sx) = /(x), for f-L-almost all x E X.
LetS be an endomorphism of (X, x) and let fL be a probability measure that is invariant under S. The symmetric difference of two sets A and B is denoted A~B, and a set A E xis said to be invariant mod 0 iff f-L(A~(S-1 A)) = 0.

An endomorphism S is called ergodic if any invariant modulo 0 set is of measure 0 or 1. An equivalent formulation of this definition is that in the Birkhoff-Khinchin theorem, f is a constant:

f(x) =.I f(x) dfJ.(x).
Alternatively, it is possible to say that fL is ergodic iff the time and the space averages coincide. The existence of an invariant measure can be proved if X is a compact metric space, xis the Borel o--algebra and Sis an homeomorphism of X. This result is due to Bogoliubov and Krilov and can be found for instance in [START_REF] Pollicott | Lectures on Ergodic Theory and Pesin Theory on Compact Manifolds[END_REF], page 8, Lemma 1.2.

In order to prove that in our problem we have an invariant measure, we would have to follow the following program: first observe that with a periodic forcing, it makes sense to think of time as a periodic variable of period 21r / w = T; the mappingS maps the position and velocity at time t 0 to their image by the flow at time t 0 + T. The second step is to show that the system defined by (1)-( 6) has a bounded invariant set in (IR/T) x Ve; consider thew-limit set A= n{un~k{Snx}} k>O for some initial data x = ( u 0 , v 0 )• and show that A is compact. This step has been performed in a slightly different case by Angles, in his thesis [START_REF] Angles | Etude qualitative des vibrations en presence d'obstacle d'un systeme vibratoire[END_REF], and we believe that the methods of Angles apply with very little adjustement to the present case.

Then the really difficult task would be to prove that restricted to A, S is an homeomorphism; it is clear that it is continuous, thanks to the continuous dependence with respect to data; but its inverse could be very bad: S could even lack an inverse if (0, o)• belonged to A. Proving or disproving that (0, o)• belongs to A looks like a very hard problem. Even if this result were proved, our invariant measure on A could be non unique and very badly behaved. However, the existence of such an invariant measure would be an important fact, because it would build the foundations for the following question: consider a function of x, such as the average kinetic energy: it is defined by where ( u 0 , v 0 ) is any representative of the class x and <P 2 is the second component of <P. The existence of an average kinetic energy, i.e. of the limit Theorem 2.1 (Oseledec). Let S be a C 1 -diffeomorphism of a compact manifold X, and let 11 be an ergodic measure. Then there are two possibilities:

(i) there exists A E IR such for all v in the tangent space TxM to M at x, [START_REF] Ruelle | Ergodic theory of differentiable dynamical systems[END_REF] for almost all x in M, or (ii) there exists .\ 1 > .\ 2 and a splitting TxM = E~ EB E;: (with the maps x f------+ E~, E;: being measurable) such that for all v 1 in E~, and all v 2 in E;:: [START_REF] Ruelle | Chaotic Evolution and Strange Attractors[END_REF] and [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF] for almost all x in M In our case, Sis not a C 1 -diffeomorphism, and we have not proved that there exists an ergodic measure /1• Therefore, Oseledec's theorem does not enables us to establish the existence of Lyapunov exponents. However, our problem is not totally devoid of regularity: let us say that a strict impact time is an impact time for which the left limit of the velocity does not vanish. It is possible to calculate the differential of the mapping (u 0 ,u!) f------+ <P(t 1 -O,t 0 ,u 0 ,u!), provided that the interval (t 0 ,t!) contains only a finite number of impact times, and that they are strict impact times.

vVhen there is exactly one impact time, the differential is computed explicitly, and it is given by a 2 x 2 matrix; the most important term in this matrix comes from the differentiation of the impact time with respect to the initial data, which is possible only if the impact time is strict.

Thus, the Lyapunov exponents could be well defined though none of the conditions of Oseledec's theorem are satisfied; but their existence and their computation is a purely experimental matter. In what follows, it should be always understood that any object we consider should be complemented by the phrase "if it exists". But we shall not repeat it systematically.

In our case, the phase space is two dimensional, which leads to the computation of two Lyapunov exponents.

The larger one, A. 1 , measures the rate of evolution of a one dimensional infinitesimal element of the phase space. It can be given by the evolution of the distance between the reference trajectory and a neighbouring one.

The sum of the smaller one, A. Then the area of the parallelogram is multiplied by exp((A. 1 + A. 2 )t).

In particular, if A. 1 > 0 and A. 1 + A. 2 < 0, the areas are contracted, while the distances can be expanded: this is precisely a situation which can lead to chaotic behavior. For a dissipative dynamical systems the areas are contracted, as a rule.

3 Differentiation of the continuous and discrete flows

The continuous flow

In this section, we calculate the differential of the mapping x 0 f------+ <P ( t 1 , t 0 , x 0 ) when there is exactly one impact time, in the interval ( t 0 , t!), and the impact is strict, i.e. <P2(tc-O,to,x) < 0. vVe denoted Xo = (uo,vo), the initial condition in the phase space.

To do this, we express the flow on such an interval as the composition of the flow from time t 0 to time tC) the time of the impact, of a reflexion and of the flow from time tc to time t 1 . Since the reflexion law is given by [START_REF] Schatzman | Uniqueness and continuous dependence on data for onedimensional impact problems[END_REF] the phase after the impact is related to the phase before the impact by a linear transformation whose matrix is R = ( 1 0 ) 0 -e [START_REF] Ya | Topics in Ergodic Theory[END_REF] It is convenient to denote with respect to the second argument, and we obtain Therefore, if we observe that on intervals without impact the differential of <P with respect to the third argument is simply the fundamental solution of the corresponding constant coefficient homogeneous linear system, we can see that ( -u(tc + 0) ) D2<P(t1, tc + 0, Rxc) = D3<P(t1, tc + 0, Rxc) -¢(tc + O) • On the other hand Moreover, we apply the matrix identity to transform the expressions containing D 3 <P 1 into expressions containing D 3 <P. Therefore we may write (6.7) as with

U = R R ( u(tc-0) 0 ) 1 R ( u(tc + 0) 0 ) 1 ( 21 )
+ ¢(tc-0) 0 u(tc-0) + ¢(tc + 0) 0 u(tc-0)

A direct computation shows that the quantity in brackets is equal to which is precisely our claim.

The discrete flow

The discrete problem corresponding to (1)-( 4) is defined by ( 7), ( 8) and ( 9). Let the discrete velocity be un+l-un

vn=---h

The initial data for ( 6) and ( 7) are U 0 and ul-uo vo=---

h (22)

  2 , and the larger one measures the rate of evolution of infinitesimal volume elements in the phase space. If DS(x) had real distinct eigenvectors v 1 and v 2 at some point x, one could understand the definition of the second Lyapunov exponent in a more geometrical fashion: the parallelogram built on v 1 and v 2 is sent by the tangent mapping at x to S into another parallelogram built on DxSv 1 = exp(A. 1 t) and DxSv 2 = exp(A. 2 t)v 2 .

  Assume that the impact at tc is strict. Denote ¢( t) = f ( t, u( t), u( t Then the differential of x f------+ <P ( t 1 , t 0 , x) is equal to Proof. vVe differentiate the identity (16) with respect to x 0 and we obtain with W = { D1 <P(tc-0, to, Xo)Dtc + D3<P(tc-0, to, Xo)} vVe calculate the different quantities which appear in (19): thanks to (15), tc satisfies <P1 (tc -0, to, Xo) = Umin, which we differentiate with respect to x 0 ; as the impact is strict, we can see that D ( ) __ D3<P1(tc-O,to,xo) tc Xo -u( tc -0) . Let us calculate now D 2 <P: we differentiate the differential equation and the initial conditions satisfied by <P 1 (t, s, x)

	<P(tc-0, to, Xo) = Xc•	(15)
	Therefore, we may write now
		(16)
	Theorem 3.1 (17)

where

(18)

The perturbed solution (Jn corresponds to perturbed initial conditions U 0 and vo.

Let N (resp N) be the set of integers n such that the constraint is active on un (respectively on (Jn), i.e. and N = { n > 1 : 2 [Jn + ( e -1) [Jn-1 + h 2 ftm < ( 1 + e) Umin}.

Of course, for n EN, un+ 1 + eun-1 = (1 + e)umin• Lemma 3.2 For all integer M such that for all n E {1, ... M}, 2Un + ( e-1 )un-1 + h 2 pn is not equal to (1 + e )umin, there exists a p such that for IU 0 -U 0 1 + IV 0 -V 0 1 < p, then N n {0, ... 'M} coincides with N n {0, ... 'M}.

Proof. The mapping which assigns to ( u, v) the number w defined by w + eu = P(l+e)K(2v + (e-1)u + h 2 f(t, v, (w-u)/2h))

is Lipschitz continuous with a Lipschitz ratio estimated by Lh, where L is the Lipschitz constant off with respect to its last two arguments. Therefore, the principle of strict contractions implies that for Lh < 1, there exists a unique solution w of (23). If u' and v' are different data, we subtract (23) from w' + eu' = P(1+e)K(2v' + (e -1)u' + h 2 f(t,v', (w'-u')/2h)), and we obtain

which implies immediately our claim about continuity.

Therefore, for all n > 1' (Jn tends to un as I U 0 -U 0 I + I V 0 -V 0 I tends to 0. If 2Un + ( e-1 )un-1 + h 2 pn is not equal to (1 + e )umin, then for c small enough, 2Un + ( e-1 )(Jn-1 + h 2 pn is not equal to (1 + e )umi~, and therefore the constraints are saturated for un if they are saturated for un.

Let us define the discrete flow <Ph as the mapping which assigns to the discrete times t 0 and tn = t(O + nh) the solution Xn = (Un, Vn)T of the numerical scheme at time tn which satisfied the initial data X 0 = (U 0 , V 0 )T at the initial time t 0 .

It should be remarked here that the phase space for the discrete flow has not been studied in the same detail as the phase space for the flow in continuous time. Our present understanding is that this phase space is IR 2 . Theorem 3.3 Consider a discrete time interval { tk, ... , tl} on which the solution un of the numerical scheme satisfies n #p,p + 1 and tk < tn < tz ==::} 2Un + (e -1)Un-1 + h 2 Fn > (1 + e)umin, (24) and n = p or n = p + 1 ==::} 2Un + ( e -1) un-1 + h 2 pn < ( 1 + e) Umin.

( 25)

The derivative of the flow with respect to its spatial argument is given by (26) where

Proof. vVe can decompose the discrete flow on this interval as and

Then, the differentiation of xz with respect to Xk is performed by composition of differentiations. The assumption on the saturation of constraints at discrete times p and p + 1 implies so that the differential of (UP+ 1 ' VP+ 1 ) with respect to (UP-

The three computational methods

In this section, we compute numerically Lyapunov exponents for the problem (1)-( 6) using three different methods. The first two methods are based on the impact detection scheme. The third method is based on the scheme ( 7)- [START_REF] Pollicott | Lectures on Ergodic Theory and Pesin Theory on Compact Manifolds[END_REF]. vVe have to calculate ln(ID(Sk)(x)vol)/k. In order to obtain a reliable computation, we have to normalize at each step the iterates; otherwise, the norm of the vectors DSk(x)v 0 increases exponentially with k, and the computation leads to an overflow. This situation is reminiscent of the power method for the computation of eigenvalues of matrices ( [START_REF] Parlett | The Symmetric Eigenvalue Problem[END_REF], [START_REF] Golub | Matrix Computations[END_REF]), and the method employed is essentially the same: let Xj be the sequence

Then, by composition of differentials, (28)

Similarly, we have to calculate the second Lyapunov exponent: once again, the comparison with the methods used for the computation of matrix eigenvalues is illuminating: it is enough to compute ln I det D(Sk)(x)l, and to subtract from this number -\ 1 ,k• Computationnally, this procedure would be highly unreliable, since the image of a basis by D(Sk)(x) becomes extremely singular as k tends to infinity; thus, we modify the computation as follows:

(29)

The numerical differentiation of the continuous time and the discrete time flows is implemented in a very simple fashion: let c be a small positive number, and let Given X 0 and t 0 , we define the numerical differential Dc: S(Xo) of the Poincare map S as the matrix given by its column vectors

The same definition applies for the discrete time flow:

The three numerical schemes will be denoted • IDED Impact Detection, Exact Differentiation,

• IDND Impact Detection, Numerical Differentiation,

• NSND Numerical Scheme, Numerical Differentiation.

In IDED, we use compute an approximation of DS to machine precision of DS, which is given by formula (6.6) when t 1 = t 0 + T, and there is only one impact in the time interval (to, t 0 + T). In IDND, we replace DS by Dc: S, and in NSND, we use Dc: Sh. vVe did not try to compute exactly the derivative of Sh: it would have required the integration of a linear difference equation which admits an explicit solution because of its very simple nature. However, we claim that for c small with respect to h, the results of our computation must be bad; therefore, it would be of little interest to let c tend to 0: this is the reason why we did not perform the relevant calculation.

Numerical results

All our computations were performed in double precision.

All the results that are described in this section are obtained with e = 0.5 and the initial data u 0 = 0 and v 0 = 0.1. vVe start with a comparison of IDED and IDND: we performed the computation for 3000 periods, a time step h used as a sweeping parameter of h = T /2513, and c = lo-s in IDND; of course T is the period T = 27r /50. By convention, Aj,k is the numerical approximation of Aj obtained by a computation on k periods. The IDND give this results A1 3ooo = 3.05694 and A2 3ooo = -15.06065, ' '

and the IDED give the results A1 3ooo = 3.05695 and A2 3ooo = -15.06066

The results obtained by respectively by IDED and IDND have 5 common significative digits at least. vVe observe that the first Lyapunov exponent is positive, and the second is negative, and their sum is negative: this is in agreement with the chaotic character of the dissipative system we considered. Another computational observation is that IDED is much faster than IDND.

If we use values of c which are very small with respect to h in the NSND scheme, the sum of the two Lyapunov exponents is positive; for instance, when h = T /2513 and s = lo-s, we obtained This fact led us to compute the Lyapunov exponents with c somewhat larger than h: with the choice c = 3h, h = T /40000, the NSND method gave reasonable results

Conclusion.

The comparison of the different methods used for the computation of the Lyapunov exponents has showed a number of discrepancies and difficulties of numerical origin. It could be argued that NSND gives bad results because it is of low order; however, IDED does not perform too well. Thus, we observed that strictly numerical and apparently harmless parameters can seriously alter the result of a computation. Therefore, the conclusion is a caveat: any computations of a Lyapunov exponent should be subjected to a serious evaluation performed by applying several different numerical methods, and by appraising the effect of all the numerical parameters, including the ones which seem innocuous.


Thomas Dubot
email: thomas.dubot@onera.frjudicaëlbedouetjudicael.bedouet@onera.fr

Judicaël Bedouet

Stéphane Degrémont
email: stephane-y.degremont@aviation-civile.gouv.fr

Modelling, generating and evaluating sector configuration plans -methodology report of the SESAR VP-755 exercise

Thomas Dubot, Judicaël Bedouet, Stéphane

General Introduction

The European airspace is currently controlled by Air Traffic Management units called Area Control Centers (ACC), in charge of providing air traffic control services to controlled flights within their area of jurisdiction. European Area Control Centers (ACCs) are subdivided into volumes of airspace called elementary sectors that can be combined to form air traffic control sectors. This subdivision of the ACC into control sectors varies throughout the day, depending on the incoming traffic and the number of available controllers. At any moment, we call configuration the set of control sectors deployed to ensure the ACC's role. For instance, Fig. 1 shows one of the possible configurations with five control sectors deployed in the French ACC of Reims. The Flow Management Position (FMP), in charge of building these airspace configurations, compares predicted flight counts with sector capacities to assess the different airspace configurations available in its database and manually choose the best configuration for each time period of the day. We call sector configuration plan (or opening scheme), of a given ACC for a given day, the description of these different configurations throughout the day.

In order to adapt the airspace configurations to the traffic situation, sectors are typically split when controllers' workload increases and merged when it decreases. In case of hotspot (overloaded area), the sector is split in two control sectors, to share the ATC workload and resolve the demand and capacity imbalance, as illustrated by Fig. 2. If this process allows to build and modify dynamically sector configuration plans, some limitations remain:

MODELLING, GENERATING AND EVALUATING SECTOR CONFIGURATION PLANS -METHODOLOGY REPORT OF THE SESAR VP-755 EXERCISE

• Hotspot resolution systematically requires to increase the number of control positions;

• Some hotspots remain at the level of elementary sectors that cannot be split; • Only a small subset of predefined configurations is considered instead of exploring all the possible combinations of elementary sectors [START_REF] Gianazza | Airspace configuration using air traffic complexity metrics[END_REF]; • Some metrics currently used to assess the controller's workload often prove insufficient, if not irrelevant [START_REF] Gianazza | An efficient airspace configuration forecast[END_REF], which makes difficult for the FMP to assess and balance airspace configurations.

SESAR context and problem statement

To face these limitations, the European SESAR program [3] implements modular and flexible dynamic airspace configurations [START_REF]Dynamic Airspace Configuration Step2-V2 OSED[END_REF]. Large airspace blocks, such as ACCs, are hence decomposed into airspace building blocks, smaller than current elementary sectors, and delineating typical demand forecast patterns, e.g. traffic flows. These building blocks, which are not necessarily controllable, are grouped into control sectors named Controlled Airspace Blocks. In this way, control sectors are more adapted to traffic specificities, which enables to solve hotspots by reorganizing the frontiers of control sectors, with the same total number of control sectors, in order to balance the ATC workload, instead of splitting one of the existing control sector, as illustrated by Fig. 3. We present in this paper methodologies and tools developed within the SESAR VP-755 exercise of the SESAR 07.05.04 project (led by EUROCONTROL), which consists in a performance assessment of sectorization algorithms based on this new paradigm. Main objectives of this exercise are to assess the benefits that could emerge from this higher granularity and this dynamicity -in terms of time, but also in terms of flexibility in shapeand to introduce automated algorithms and tools that could support this evolution.

Modelling of building blocks

The focus of the VP-755 exercise was made on sector configuration processes. Nevertheless the sector design phase is a crucial step to ensure the relevance of elementary building blocks. The first data set used was realized from the Reims operational data of the 21 elementary sectors. The objective was to validate optimization algorithms described in §4.2 and to verify if such techniques could give relevant airspace configurations out of the classic FMP catalog.

Then a manual sector design was made to assess the benefits that could emerge from a higher granularity. 42 building blocks were created by the FMP from the current 21 elementary sectors, as shown by Fig. 4 hereafter. Basically, elementary sectors were divided according to the FMP expertise on specific traffic flows and complexity issues.

If the VP-755 relied on manual sector design, it has to be noted that the overall architecture has been realized to be compatible with any sector design. A next step could be to work with building blocks automatically generated, for instance by the EUROCONTROL ASTAAC algorithms [START_REF] Sergeeva | 3D Sectors Design by Genetic Algorithm Towards Automated Sectorisation[END_REF].

Once the building blocks have been defined and modelled with the data model presented in previous paragraph, interrelations between these building blocks were automatically given by an in-house algorithm based on the MATLAB program described in [START_REF] Sale | Development of a Computer Based Airspace Sector Occupancy Model[END_REF]. The result of this algorithm is an adjacency table between all building blocks, summarizing which building block is adjacent (a common surface in horizontal or vertical face) to another, as illustrated by Fig. 5. This sector configuration plan is composed of 24 airspace configurations associated to 24 variable time periods, and based on the deployment of 1 to 16 control positions. The sector configuration plans realized during the VP-755 exercise are based on these data. The division in time periods is the same as the reference one, and for each time period, the number of sectors is the same as the real ATC roster of the day considered, even if an optimization of this parameter would make sense in another context. For each run, the sector configuration plan is generated by the following process:

• Optimization algorithms provide the FMP with a set of good airspace configurations for each time period, and a default sector configuration plan, realized by selecting the best combination of these solutions to minimize the changes between time periods. • Then the FMP analyzes these airspace configurations, through a decision support tool presented in §4.3, and selects one airspace configuration for each time period, to build the final sector configuration plan.

Optimization principles

As described in our previous work [START_REF] Dubot | Building a smooth and dynamic opening scheme from graph partitioning-Exploring dynamic airspace configurations[END_REF], our optimization algorithms are based on the resolution of a combinatorial problem of graph partitioning, as illustrated by Fig. 7.

Fig. 7. Graph partition modelling

The graph G = (V, E) denotes the representation of the airspace, where V (the set of vertices) is the set of airspace building blocks and E (the set of edges), is such that (u, v) belongs to E only if there can be a direct trajectory from u to v.

The graph is valuated both on its vertices and edges as follows:

{ } k k S S P ,..., ) t ( 1 = ∆
. Each subset S i of the partition must be disjoint from the empty set and from the other subsets and the union of all the subsets must entirely cover the graph G and satisfies the connectivity constraint: the different elements of a subset, the aforementioned vertices, must be connected. To build the optimal partition of airspace into k sectors, we may consider different objectives to minimize:

• The workload imbalanced distribution measured by the sum of distances to the average of the density of each sector within each sector configuration period; • The total number of transfers measured by the sum of flights transiting from one sector to another within each sector configuration period; • The total number of overloads defined by the number of overloads of traffic (over a given threshold) in a sector during a given period of time; • The total number of reentries : a reentry corresponds to a flight that enters at least twice in the same sector; • The total number of short transits: a short transit corresponds to a flight that spends less than four minutes in a sector.

The complexity of such an optimization problem is considered as NP-complete [START_REF] Dubot | Building a smooth and dynamic opening scheme from graph partitioning-Exploring dynamic airspace configurations[END_REF]. If we don't consider the connectivity constraint, the number of possible solutions is given by the second Stirling number

∑ = -         - = k j n j k j j k k k n S 0 ) 1 ( ! 1 ) , (
where n is the number of building blocks k is the number of wanted partitions In reality, the controllers do not exploit much more than 80 sectors. In that case, the number of possible solutions is given by the binomial

        j k is the binomial coefficient )! ( ! ! j k j k - . For instance,
coefficient         k 80
. The number of possible solutions is still huge but if we consider the connectivity constraint, we have a limited set of possible sector configurations, as shown in Fig. 8. Those solutions are obtained by cutting in a sector tree the branches that definitively conduct to a configuration sector that does not satisfy the constraints. If the FMP catalog of sectors should be increased, it would be possible to reduce the exploration by eliminating configuration sectors that clearly conducts to an imbalance, independently of the traffic, e.g. a large sector with a small sector.

For each time period, we are then able to assess a limited number of sector configurations in respect of the previous minimization objectives and build successive Pareto fronts. A Pareto front contains all the solutions that cannot be dominated by any other solution of the same front. In other words, each front shows the optimal solution and compromises between the different objectives. So, we have different good solutions for each time period through these Pareto fronts. The second step is then to combine them to form the sector configuration plan of the day, as illustrated by Fig. 9. Between two time periods, we must minimize the distance between each sector configuration. It is a classic of the Bellman's Principle of Optimality [START_REF] Bellman | Dynamic programming and modern control theory[END_REF]. Given two partitions P and P' of the same graph G, we define the distance ( ) ' , P P D between these two partitions as the smallest sum of weights of any nodes of G whose removal causes the two induced partitions to be identical [START_REF] Gusfield | Partition-distance: A problem and class of perfect graphs arising in clustering[END_REF]. From an ATM stance, this is the total density workload of the building blocks which differ when switching from the first partition to the second one. To be more operational, the partition-distance has been revised to favor what we called the collapsing/de-collapsing operations. A collapsing operation corresponds to two sectors that are merged together to form a unique sector. A de-collapsing option corresponds to the inverse one.

At the end of this deterministic phase, we have at our disposal:

• For each time period, a set of very good airspace configurations; • A sector configuration plan built upon these solutions to be as smooth as possible.

Nevertheless, as these airspace configurations rely on the current FMP catalog, we complement this process with a stochastic approach to explore unknown solutions.

Our stochastic algorithms rely on the Simulated Annealing metaheuristic, as described in [START_REF] Dubot | Building a smooth and dynamic opening scheme from graph partitioning-Exploring dynamic airspace configurations[END_REF].

The objective is to explore new airspace configurations by exchanging building blocks, in order to check if criteria mentioned previously, e.g. the workload imbalanced distribution, could be improved without degrading the flow of airspace configurations throughout the day, and consequently their acceptance by ATCOs. Two main approaches have been considered. In the first one, the algorithm is initialized with one of the good solutions provided by the determinist algorithm described previously. For instance, the time period with largest number of control sectors is selected. Then we build the neighbor solutions (previous and next time period), by exploring solutions around this good solution, modified to match the requested number of control positions. The advantage of this technique is to introduce from the optimization phase a consideration on the stability of successive configurations. Nevertheless one of the main issues is that depending on the selected solutions, the distance to the selected reference can rapidly be irrelevant. In a second approach, we therefore take for each time period a reference configuration given by the determinist algorithm, i.e. the smoothest sector configuration plan found. We then explore through the same Simulated Annealing method if small exchanges around this "backbone" can improve significantly one or several objectives, while being close to operational situational awareness of ATCOs, and smooth by construction.

Decision support tool

Once optimization algorithms have provided a set of airspace configurations for each time period, the run consists in the human-in-theloop process of analyzing these configurations and selecting the most suitable one to form the final sector configuration plan. The FMP has at his disposal a set of tools to facilitate this process, as pictured in Fig. 10 hereafter.

Fig. 10. VP-755 FMP working position

The working position is composed of a 3D visualization tool, enabling the visualization of traffic and airspace data of each time period considered, and of an ad hoc HMI called Sector Configuration Plan creator, as pictured in Fig. 11 hereafter, synthesizing all airspace configurations obtained through optimization algorithms, and providing the FMP with a set of information to help him choose the best association of configurations.

Fig. 11. VP-755 Sector Configuration Plan creator

The interface displays dynamically the following information:

• For each airspace configuration, the differences with the previous and next configurations, and the associated distance computed between these configurations, as described in previous paragraph; • For each control sector of the airspace configuration, the density of flights, the number of transfers during the time period, the vertical compactness of the airspace volume, the number of overloads, the occupancy…

Once the FMP has selected, created or validated an airspace configuration for each time period, the final sector configuration plan of the run is recorded according to the data model described in §5.1.

Evaluation of sector configuration plans

The final step of our exercise is the evaluation of the sector configuration plan obtained through optimization algorithms and FMP expertise, as described previously.

Data model

The data model used in the VP-755 platform is based on the Aeronautical Information Exchange Model (AIXM) [11] and extended with new objects such as the building blocks.

The consistency of data used by different modules is ensured by the use of GAMME, an in-house meta-modelling tool, which enables to generate the different pieces of code that will be used to develop the simulation software and exchange data between the different components [START_REF] Bedouet | GAMME, a meta-model to unify data needs in simulation modeling (WIP)[END_REF].

Overall simulation architecture

Fig. 12 hereafter gives an overview of the VP-755 platform architecture. We describe in §3 and §4 the processes of sector design and airspace configuration generation through optimization algorithms. We will focus in this paragraph on the processes required to evaluate the sector configuration plan generated. The VP-755 platform is based on the data model previously described and the ONERA IESTA platform [START_REF] Bedouet | IESTA: evaluating air transport systems using a rich, modular, generic, managed platform[END_REF]. Fast-time simulation capabilities allow for instance to simulate the traffic with different levels of fidelity. One can choose:

• to use the IESTA Aircraft Simulation Module, based on EUROCONTROL BADA [14], in order to consider these 4D trajectories as new flight plans and to simulate the deviations to these flight plans; • to exactly follow them, only interpolating between these 4D positions.

Traffic generation

The day of traffic sample has been chosen according to local Reims ACC observation of recent highly loaded days. We used the EUROCONTROL Demand Data Repository (DDR2) [15] to build the trajectories on one specific day (26/06/15). More specifically, we used the M1 data (from Flight Plan) during the optimization process, and the M3 (Flight Plan updated with Radar Data from CFMU) during the evaluation phase. These DDR2 trajectories were then filtered, in order to keep only the flights with at least one point in the studied area. We then extrapolated these traffic data to build a Free Route traffic, as pictured by Fig. 13 hereafter.

Metrics

The Reims airspace is described as a set of building blocks composed of volumes. We computed the intersections between these volumes and the trajectories. Both are projected following a gnomonic projection centered on Reims. The advantage of this projection is to transform all great circles into straight lines. The trajectories are transformed into a set of 4D segments while the faces of the volumes are transformed into a set of oriented polygons. Then, we determine which oriented polygons are crossed by the different 4D trajectories in order to determine when and where a flight leaves a building block and enters another one. We can hence determine very quickly the different entries/exits inside the Reims airspace for a full day of operations, which is one of the basis of the different metrics described in next paragraph.

Topological metrics and constraints

The evaluation phase requires to define and implement relevant metrics, in order to evaluate at any moment the resulting configurations of control sectors obtained by aggregating airspace building blocks. The first metrics to be considered are geometrical ones, used in the airspace design phase to build the airspace building blocks [START_REF] Delahaye | Genetic algorithms for partitioning airspace[END_REF] and [START_REF] Trandac | Optimized sectorization of airspace with constraints[END_REF] identified the following constraints in the construction of sectors:

• Convexity constraint: an aircraft cannot enter the same sector twice; • Minimum distance constraint: the distance between a sector border and a network node must not be less than a given distance; • Minimum sector crossing time constraint: the aircraft must stay in each crossed sector at least a given amount of time, as illustrated by Fig. 15; • Connectivity constraint: the sector cannot be fragmented. As we focused in VP-755 on the airspace configuration phase, we consider that the airspace building blocks that we use have been defined according to such criteria. Nevertheless we still need to verify that the controlled airspace blocks formed by the gathering of these Airspace Building Blocks will not degrade any of these properties. For instance, when exchanging one building block from one sector to another, we always check that the resulting sectors are not fragmented (connectivity). Besides, as we consider both fixed route and free route operations, some measures such as the minimum distance constraint can be strongly different when analyzing historical traffic flows and real traffic flows on the day of operations, depending notably on weather conditions. Our evaluation module implements therefore the following metrics:

• Total number of short transits: number of flights transiting in a sector less than n seconds; • Total number of traffic nodes too close to sector borders: number of routes' intersections located at less than a given distance to the sector's frontiers; • Total number of re-entries: number of flights re-entering in the same sector within each sector configuration period.

Operational metrics

One of the objectives of the VP-755 evaluation was to enrich common operational metrics, such as the hourly entry count or the occupancy count [START_REF] Dalichampt | Hourly entry counts versus occupancy count -Relationship, definitions and indicators[END_REF], that prove to be insufficient to assess controllers' workload. We therefore introduced complexity metrics, such as the dynamic density [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], to assess in each timeframe the traffic density and the inherent complexity of tasks allocated to air traffic controllers. The use of this metric nevertheless requires allocating adequate weights to the different subparameters, as listed hereafter in Fig. 16. More generally, the evaluation module provides many metrics linked to the number of flights, such as those presented in §4.2: density, workload imbalanced distribution, total number of overloads, total number of transfers… Finally we implemented a metric linked to the stability of successive configurations in terms of geometrical shape. This metric is based on the Hausdorff distance described in [START_REF] Yousefi | Trigger metrics for dynamic airspace configuration[END_REF].

Conclusions

First results, based on both quantitative and qualitative analysis, seem promising, in terms of methodology and tools. The generation of sector configuration plans, through optimization algorithms and FMP expertise, gives for instance interesting results in terms of workload distribution, as illustrated by Fig. 17 hereafter. Besides, human factor analysis shows that decision support tools, such as those presented in VP-755 exercise, could facilitate the FMP tasks. Further studies should analyze how such tools could be integrated to the current FMP working tooling. Within the framework of the SESAR ATM system, such tools and methodology seem necessary to deal with the dynamicity required by unconstrained free route operations. Further studies should also assess the possibility to rapidly generate such opening schemes, with an important number of metrics and with building blocks becoming much smaller, as the pixel view illustrated in Fig. 18 hereafter. In any case, such techniques can only be seen as a way to enrich ATC centers' catalogs and as a decision support tool to operational experts. The FMP expertise will always be required to assess the overall relevance of the airspace configurations generated and their potential acceptation by air traffic controllers. Besides such a local optimization should be integrated within a bigger optimization loop, at a Functional Airspace Block (FAB) or European level, to consider side and network effects.

Fig. 1 .

 1 Fig. 1. Reims airspace configuration in 5 control sectors.

Fig. 2 .

 2 Fig. 2. Hotspot resolution by splitting one overloaded control sector (one additional sector).

Fig. 3 .

 3 Fig. 3. Hotspot resolution by reorganizing control sectors (same number of sectors).

Fig. 4 .

 4 Fig. 4. Manual sector design of Reims ACC in 42 building blocks

Fig. 5 . 4 Generation of sector configuration plans 4 . 1

 5441 Fig. 5. Adjacencies of the 21 Reims ACC building blocks

Fig. 6 .

 6 Fig. 6. Reference operational sector configuration plan

Fig. 8 .

 8 Fig. 8. Logarithm of the theoretical number of possible configurations in blue (right scale) and number of valid sector configurations in green (left scale).

Fig. 9 .

 9 Fig. 9. Smoothest combination of airspace configurations to form the sector configuration plan.

Fig. 12 .

 12 Fig. 12. Architecture of the SESAR VP-755 platform

Fig. 13 .

 13 Fig. 13. VP-755 Free Route trajectories Algorithms were developed to shift in time and space the different points of the trajectory, based on the great circle navigation from FL310. Fig. 14 shows for instance the modifications of the CES569 flight.

Fig. 14 .

 14 Fig. 14. Free Route trajectory extrapolation

Fig. 15 .

 15 Fig. 15. Minimum sector crossing time objective

Fig. 16 .

 16 Fig. 16. Parameters of the NASA Dynamic Density metric

Fig. 17 .

 17 Fig. 17. Significant gain in workload distribution with new sector configuration.

Fig. 18 :

 18 Fig. 18: Reims ACC airspace subdivided in 400 20x20NM cuboids

MODELLING, GENERATING AND EVALUATING SECTOR CONFIGURATION PLANS

Acknowledgments

Authors would like to express their gratitude to all members of the VP-755 exercise: Isabelle Luxembourg, Sylvain Mondou, Valérie Neyns from DSNA, Sébastien Aubry, Luis Basora, Thomas Chaboud, Antoine Joulia, Nicolas Huynh, Rémi Lafage, Xavier Olive from ONERA. This paper relies on their excellent work and their constructive suggestions during the reviewing of this paper. Authors are also grateful to all participants to VP-755 exercises: Marie-Laurence Bossy, Charlotte Chambelin, Yannick Migliorini, Bruce Palmouries from DSNA, Yevgen Pechenik and Leïla Zerrouki from EUROCONTROL. This study was co-financed by the European Union, through the SESAR Joint Undertaking. Authors would like to thank Kris Delcourte, project manager of the SESAR 07.05.04 project and all the European partners involved in this project: EUROCONTROL, DFS, DSNA, ENAIRE, ENAV, ERC, INDRA, INECO, NATS, NORACON, RAF, SICTA, Spanish AF, UK MOD.