D Delahaye
email: delahaye@recherche.enac.fr

S Puechmorel
email: puechmor@recherche.enac.fr

P Tsiotras
email: tsiotras@gatech.edu

E Feron
email: feron@gatech.edu

Mathematical Models for Aircraft Trajectory Design : A Survey

Air traffic management ensures the safety of flight by optimizing flows and maintaining separation between aircraft. After giving some definitions, some typical feature of aircraft trajectories are presented. Trajectories are objects belonging to spaces with infinite dimensions. The naive way to address such problem is to sample trajectories at some regular points and to create a big vector of positions (and or speeds). In order to manipulate such objects with algorithms, one must reduce the dimension of the search space by using more efficient representations. Some dimension reduction tricks are then presented for which advantages and drawbacks are presented. Then, front propagation approaches are introduced with a focus on Fast Marching Algorithms and Ordered upwind algorithms. An example of application of such algorithm to a real instance of air traffic control problem is also given. When aircraft dynamics have to be included in the model, optimal control approaches are really efficient. We present also some application to aircraft trajectory design. Finally, we introduce some path planning techniques via natural language processing and mathematical programming.

Introduction

Aircraft trajectory is one of the most fundamental objects within the frame of ATM. However, partly due to the fact that aircraft positions are most of the time represented as radar plots, the time dependence is generally overlooked so that many trajectory statistics conducted in ATM are spatial only. Even in the most favorable setting, with time explicitly taken into account, trajectory data is expressed as an ordered list of plots labeled with a time stamp, forgetting the underlying aircraft dynamics. Furthermore, the collection of radar plots describing the same trajectory can have tenths more samples, nearly all of them redundant. From the trajectory design point of view, this redundancy is real handicap for the optimization process. In this survey, alternative trajectory representations are presented with a description of their advantages and limits. Such new approaches may be applied in many areas :

Aircraft trajectories data compression. As it has been previously mentioned, ATM system manage aircraft trajectories and control them in order to guarantee safety and airspace capacity. Currently those trajectories are represented by the mean of plot lists which are manipulated by ATM software. Every day, all aircraft trajectories are registered into large database for which huge capacity is needed. Based on this new trajectory representation for which redundancy has been removed, the trajectories database may be strongly improved from the capacity point of view. This compressed trajectory format may also be used for improving the trajectories transmission between ATM entities.

Aircraft trajectories Distance Computation. Although trajectories are well understood and studied, relatively little investigation on the precise comparison of trajectories is presented in the literature. A key issue in performance evaluation of ATM decision support tools (DST) is the distance metric that determines the similarity of trajectories. Some proposed representation may be used to enhance trajectory distance computation.

Aircraft model Inference. All aircraft models are based on ODEs(Ordinary Differential Equation), including tabular ones. Control input includes condition and model parameters. The model refinement (and computational complexity) ranges from tabular to many degrees of freedom. The aircraft model inference consists in answering the following question: Given a parametrized model and a goal trajectory, can we infer the best parameter values? A model can be viewed as a mapping from the control space into the trajectory space. The way to answer the previous question is then given by the closest model to the goal trajectory.

Trajectory prediction. Air traffic management research and development has provided a substantial collection of decision support tools that provide automated conflict detection and resolution [START_REF] Erzberger | Conflict detection in the presence of prediction error[END_REF][START_REF] Brudnicki | User request evaluation tool (uret) conflict probe performance and benefits assessment[END_REF][START_REF] Vink | Eatchip medium term conflict detection: Part 1 eatchip context[END_REF],trial planning [START_REF] Mcnally | Field test evaluation of the ctas conflict prediction and trial planning capability[END_REF], controller advisories for metering and sequencing [START_REF] Swensen | Design and operational evaluation of the traffic management advisor at the forth worth air route traffic control center[END_REF][START_REF] Coppenbarger | Design and development of the enroute descent advisor (eda) for conflict-free arrival metering[END_REF],traffic load forecasting [START_REF] Meckiff | The tactical load smoother for multi-sector planning[END_REF][START_REF] Masalonis | Using probabilistic demand prediction for traffic flow management decision support[END_REF], weather impact assessment [START_REF] Kirk | Problem analysis resolution and ranking (parr) development and assessment[END_REF][START_REF] Sud | Air traffic flow management collaborative routing coordination tools[END_REF][START_REF] Evans | Reducing severe weather delays in congested airspace with weather support for air traffic management[END_REF]. The ability to properly forecast future aircraft trajectories is central in many of those decision support tools. As a result, trajectory prediction (TP) and the treatment of trajectory prediction uncertainty continue as active areas of research and development (eg [START_REF] Slattery | Trajectory synthesis for air traffic automation[END_REF][START_REF] Swierstra | Common trajectory prediction capability for decision support tools[END_REF][START_REF] Mondoloni | Trajectory modeling accuracy for air traffic management decision support tools[END_REF][START_REF] Ryan | Review of trajectory accuracy methodology and comparison of error measure metrics[END_REF][START_REF] Mondoloni | Impact of factors, conditions and metrics on trajectory prediction accuracy[END_REF]). Accuracy of TP is generally defined as point spatial accuracy (goal attainment) or as trajectory fol-lowing accuracy. The last one can be rigorously defined by the mean of trajectory space. The first one is a limit case of the second by adding a weight function in the energy functional. Since we may prescribe smoothness accuracy of a simplified model relative to a finer one, may be computed.

Major flows definition. When radar tracks are observed over a long period of time in a dense area, it is very easy to identify major flows connecting major airports. The expression "major flows" is often used but never rigorously defined. Based on an exact trajectory distance and a learning classifier, it is possible to answer the following questions: Given a set of observed trajectories, can we split it into "similar" trajectory classes? If yes, classes with highest number of elements will rigorously define the major flows. Given those classes and a new trajectory, can we tell if it belongs to a major flow and which one? The principle of the major flows definition is to use shape space to represent trajectory shapes as points and to use a shape distance (the shape of a trajectory is the path followed by an aircraft, that is the projection in the 3D space of its 4D trajectory. The speed on the path has no impact).

Trajectory planning. To improve Air Traffic Management, projects have been initialized in order to compel the aircraft in position and in time (4D trajectory) so as to avoid potential conflict and allow for some optimality with respect to a given user cost index, environmental criteria (noise abatement, pollutant emission . . . ). Depending on the time horizon, several kind of plannings can be designed:

• At a strategical level, only macroscopic indicators like congestion, mean traffic complexity, delays can be taken into account, considering the high level of uncertainty; • at a pre-tactical level, the accuracy of previous indicators, specially congestion and complexity increases while at the same time early conflict detection can be performed; • finally, at the tactical level, conflict resolution is the major concern and optimality of the trajectories is only marginally interesting.

As we can see, there are many areas of ATM where trajectories are the main objects that have to be manipulate.

The second part of this survey presents some relevant features of aircraft trajectories. The third part, presents dimension reduction tricks for optimization approaches. The fourth part describes approaches based on wave front propagation in isotropic and anisotropic environments. The fifth part presents automatic control approaches with some application to air traffic control. Finally, the sixth part introduces some path planning techniques via natural language processing and mathematical programming. Finally, the seventh part gives some air traffic management applications of such trajectory design approaches.

Some trajectories features

In the following all aircraft trajectories will be described as mappings from a time interval [a, b] to a state space E with E either R 3 or R 6 depending on whether speed is assumed to be part of aircraft state or not. Extension to trajectories on a sphere (typically long haul flights) will be sketched only.

Notations and terminology

The reference for this section is [START_REF] Berger | Differential geometry: Manifolds, curves and surfaces[END_REF]. Let γ[a, b] → E be a trajectory. The origin of the trajectory is γ(a) and the destination is γ(b). Those two points are called the endpoints of the trajectory. All trajectories are assumed to be at least continuously differential (class C 1 ) so that the length of a trajectory γ[a, b] → E is well defined as:

l(γ) = b a γ ′ (t) dt (1) If γ ′ (t) = 0 for some t ∈ (a
s(t) = t a γ ′ (t) dt (2)
We see that by the non-singularity assumption on γ, s ′ (t) = γ ′ (t) > 0 for any t ∈ (a, b), so that s is an invertible mapping. Now, γ • s -1 is a mapping from the open interval (0, l(γ)) to E satisfying:

(γ • s -1 ) ′ = (γ ′ • s -1 ) • (s -1 ) ′ = 1 (3)
In the following, we will simply write γ(s), s ∈ (0, l(γ)) for a curve parametrized by arc length, dropping the variable t.

Remark 1. One must be careful with the respective definitions of trajectories and curves: a curve is defined on an open interval and thus has no endpoints. Nevertheless, any trajectory γ[a, b] → E has an associated curve, namely γ(a, b) → E. It is generally more convenient to deal with curves to avoid special treatment of the endpoints.

Remark 2. The non singularity assumption on the underlying curve is very natural when dealing with aircraft trajectories in R 3 since it is not possible for an aircraft to stop except at the endpoints of the trajectory.

Remark 3. While the case E = R 3 is very natural and intuitive, care must be taken when E = R 6 since all the preceding definitions apply in a completely different setting: for example, the non singularity assumption does not implies nowhere zero speed, but only that speed and acceleration cannot both vanish at the same time. The arc length parametrization allows to define very important geometrical quantities when E = R 3 .

Definition 1. Let:

γ(0, l) → R 3 (4)
be a C 1 curve parametrized by arc length. The unit tangent vector to γ at s ∈ (0, l) is :

τ(s) = γ ′ (s) (5)
It is clear from the definition of parametrization by arc length that τ(s) is a unit vector.

Definition 2. Let:

γ(0, l) → R 3 (6)
be a C 2 curve parametrized by arc length. The curvature of γ at s ∈ (0, l) is :

K(s) = γ ′′ (s) (7)
The curvature can be explicitly computed even if the curve γ is not parametrized by arc length. The general formula is:

K(t) = γ ′ (t) ∧ γ ′′ (t) γ ′ (t) 3 (8)
with ∧ the vector cross product. Curvature is of primary importance for ATM related studies since as mentioned before aircraft trajectories are mainly made of straight lines and arcs of circle and so have piecewise constant curvature. If at point t the curvature is not zero, the curve is said to be biregular at t. For a curve γ parametrized by arc length, the unit normal vector ν(s) is defined at all biregular points by :

ν(s) = γ ′′ (s) K(s) (9)
Remark 4. A straight line has everywhere zero curvature. However, it is clearly possible to define a unit normal vector. At a biregular point, τ(s) and ν(s) are well defined. Taking their cross product gives a new vector β (s) = τ(s) ∧ ν(s). If the curve γ is assumed to be C 3 , it can be shown that β (s) and ν(s) are collinear :

β (s) = T (s).ν(s) (10)
The real number T (s) is called the torsion of the curve at s and represents an obstruction for the curve to be planar. As for the curvature, it is possible to compute the torsion even if the curve is not parametrized by arc length :

T (t) = - det(γ ′ (t), γ ′′ (t), γ ′′′ (t)) γ ′ (t) ∧ γ ′′ (t) 2 (11)
Torsion is not so useful as curvature for en-route data analysis since only a few number of trajectories have non zero torsion. However, it is very relevant in terminal areas.

Remark 5. The E = R 6 case is again very different, since the geometric meaning of curvature and torsion is not obvious in this setting. Furthermore, the extra degrees of freedom will impose using higher order derivatives in order to build up an equivalent description. A complete treatment goes beyond the scope of the present paper and has little interest for our purpose (in practical applications, the speed information, when available, is used to improve estimates of curvature and torsion and not to study a trajectory in R 6 ).

Trajectory Models for Optimization

This section presents some dimension reduction tricks in order to reduce the dimension of the state space for which an optimization process is searching for an optimal vector of parameter.

This approach is summarized by the figure 1.

The optimization process controls the parameter vector which is then used to build the trajectory γ for evaluation.

Each coordinate can be considered separately in order to build a given trajectory:γ(t) = [x(t), y(t), z(t)] T .

In this section, several trajectory models are presented and compared. Simple models are first presented.

Straight line segments

One of the easiest way to design trajectory is to use way points connected by straight lines (see figure 2. This easy principle ensures continuity for the trajectory but not for its derivatives. If one want to approximate trajectory with many shape turns, one have to increase the number of way points in order to reduce the error between the model of the real trajectory.

In order to improve concept Lagrange interpolation process adjust a polynomial function to a given set of way points.

Lagrange interpolation

Given n + 1 real numbers y i ,0 ≤ i ≤ n, and n + 1 distinct real numbers x 0 < x 1 < ... < x n , Lagrange polynomial [START_REF] Jeffreys | Methods of mathematical physics[END_REF] of degree n (L n (x)) associated with {x i } and {y i } is a polynomial of degree n solving the interpolation problem :

L n (x i ) = y i , 0 ≤ i ≤ n (12) L n (x) = n ∑ i=0 y i .l i (x) (13)
where

l i (x) = ∏ j =i (x -x j ) (x i -x j ) (14)
An example of Lagrange interpolation is given on figure 3 for which four points are interpolated by the black curve which represents L 4 (x). The four polynomial functions {l 0 (x), l 1 (x), l 2 (x), l 3 (x)} are also given by the red, blue, green and yellow curves.

When derivatives have also to be interpolated, Hermite interpolation has to be used.

Hermite interpolation

Hermite interpolation [START_REF] Bartels | An introduction to splines for use in computer graphics and geometric modeling[END_REF] generalizes Lagrange interpolation by fitting a polynomial (H(x)) to a function f that not only interpolates f at each knot but also interpolates a given number of consecutive derivatives of f at each knot. This means that the first derivative of the polynomial H(x) have to fit the first derivatives of the function f (x) :

∂ j H(x) ∂ x j x=x i = ∂ j f (x) ∂ x j x=x i ( 15
)
for all j = 0, 1, ..., m and i = 1, 2, ..., k This means that n(m + 1) values (x 0 , y 0 ), (x 1 , y 1 ), . . . , (x n-1 , y n-1 ), (x 0 , y ′ 0 ), (x 1 , y ′ 1 ), . . . , (x n-1 , y ′ n-1 ), . . . . . . . . .

(x 0 , y (m) 0 ), (x 1 , y (m) 1 ), . . . , (x n-1 , y (m) n-1 ) (16)
must be known, rather than just the first n values required for Lagrange interpolation. The resulting polynomial may have degree at most n(m + 1) -1, whereas the Lagrange polynomial has maximum degree n -1. These interpolation polynomials seem attractive but they both induce oscillations between interpolation points (Runge 's phenomenon. Runge's phenomenon is a problem of oscillation at the edges of an interval that occurs when using polynomial interpolation with polynomials of high degree (which is the case for Lagrange and Hermite interpolation). An example of such Runge's phenomenon is given on figure 4 for which Lagrange interpolation has been used. We can conclude that interpolation with high degree polynomial is risky. In order to avoid this drawback of high degree polynomial interpolation one must use piecewise interpolation.

Piecewise Linear Interpolation

This is the simplest piecewise interpolation method.

Given n + 1 real numbers y i ,0 ≤ i ≤ n, and n + 1 distinct real numbers x 0 < x 1 < ... < x n , we consider the n linear curves lin i (x) = a i x + b i on the intervals [x i , x i+1 ] for i = 0, ...n -1 (lin i (x) represent linear functions for which a i is the slope and b i a constant).

Each l i (x) has to connect two points ((x i , y i ),(x i+1 , y i+1 )

y i = a i x i + b i x i and y i+1 = a i x i+1 + b i x i+1 (17)
In order to associate a piecewise formulation of this interpolation method, the following "tent" functions are defined :

ψ i (x) =      x-x i-1 x i -x i-1 i f x ∈ [x i-1 , x i ] x i+1 -x x i+1 -x i i f x ∈ [x i , x i+1 ] 0 elsewhere (18)
Then,

f (x) = i=n ∑ i=0 y i .ψ i (x) (19)
An example of such a linear piecewise interpolation is given on figure 5 The derivative of the resulting curve is not continuous. In order to fix this drawback, one can use piecewise quadratic interpolation x 0 x i x i+1 x i-1

x n

Piecewise Quadratic Interpolation

We consider the n quadratic curves ψ i (x) = q i (x) = a i x 2 + b i x + c i on the intervals [x i , x i+1 ] for i = 0, ...n -1.Each q i (x) has to connect two points ((x i , y i ),(x i+1 , y i+1 ); ⇒ y i = a i x 2 i + b i x i + c i and y i+1 = a i x 2 i+1 + b i x i+1 + c i . Furthermore, on each point, the derivative of the previous quadratic has to be equal to the derivative of the next one; ⇒ 2a i +b i = 2a i-1 +b i-1 . For the first segment the term 2a i-1 +b i-1 is arbitrarily chosen (this will affects the rest of the curve). An example of piecewise quadratic interpolation is given on figure 6. The main drawback of piecewise quadratic interpolation is linked to the effect induced on the curve by moving on point. As a matter of fact moving one point may totally change the shape of the interpolating curve. The piecewise cubic interpolation avoid this drawback.

Piecewise cubic interpolation

This interpolation is also called Hermite cubic interpolation [START_REF] Hearth | Scientific computing, an introductory survey[END_REF]. For this interpolation :

ψ i (x) = C i (x) = a i x 3 + b i x 2 + c i x + d i ( 20
)
and we have the following constraints :

C i (x i ) = y i C i (x i+1 ) = y i+1 ( 21
)
C ′ i (x i ) = y ′ i = y i+1 -y i-1 x i+1 -x i-1 C ′ i (x i+1 ) = y ′ i+1 = y i+2 -y i x i+2 -x i (22)
An example of piecewise cubic interpolation is given on figure 7.

Fig. 6 Piecewise quadratic interpolation. The shape of the entire curve depend of the choice of the initial slope.Between two points, a quadratic polynomial is fitted.

x 0 x i x i+1 x i-1

x n

Initial slope

Fig. 7 Piecewise cubic interpolation. The derivative at point x i is given by line joining the point (x i-1 , y i-1 ) and (x i+1 , y i+1 ).Between two points, a cubic polynomial is fitted. The term h represents the distance two consecutive points.

x i x i+1 x i-1 i y y i-1 y i+1 y i+2 x i+2 h slope in i slope in i+1
Moving a point do not affect all the curve which is the main advantage of this interpolation.The resulting curve is C 1 but not C 2 (the second derivative is not continuous). The curvature radius of a curve may be expressed by the following expression :

R = 1 + d f (x) dx 3 2 | d 2 f (x) dx 2 | (23)
The piecewise cubic interpolation do not insure that trajectory curvature is continuous which is not adapted for aircraft trajectory mainly in TMA 1 areas and cubic spline interpolation has to be used.

Cubic Spline Interpolation

This method has been developed by General Motor in 1964 [START_REF] Birkhoff | Piecewise polynomial interpolation and approximation[END_REF]. For this piecewise interpolation ψ i (x) = S i (x) with the following constraints :

S i (x i ) = y i S i (x i+1 ) = y i+1 S ′ i (x i ) = S ′ i-1 (x i+1 ) S ′ i (x i+1 ) = S ′ i+1 (x i+1 ) S ′′ i (x i ) = S ′′ i-1 (x i+1 ) S ′′ i (x i+1 ) = S ′′ i+1 (x i+1 ) (24)
One can show that S i (x) for x ∈ [x i , x i+1 ] is given by :

S i (x) = σ i 6 . (x i+1 -x) 3 x i+1 -x i + σ i+1 6 . (x-x i ) 3 x i+1 -x i + y i . x i+1 -x x i+1 -x i -σ i 6 .(x i+1 -x i )(x i+1 -x) + y i+1 . x-x i x i+1 -x i - σ i+1 6 .(x i+1 -x i )(x -x i ) (25)
where

σ i = d 2 S i (x) dx 2 (26)
An example of such interpolation is given on figure 8. Such spline is also called natural spline because it represents the curve of a metal spline constrained to interpolate some given points. x i+2

1 TMA : "Terminal Maneuvering Area"

When interpolation is not a hard constraint, one can use some control points which change the shape of a given trajectory without forcing this trajectory to go through such control point; such approach is called approximation for which one of the famous methods is the Bézier curve.

Bézier Approximation Curve

Bézier curves [START_REF] Farin | The essentials of cagd[END_REF] were widely publicized in 1962 by the French engineer Pierre Bézier, who used them to design automobile bodies. But the study of these curves was first developed in 1959 by mathematician Paul de Casteljau using de Casteljau's algorithm [START_REF] Farin | Curves and surfaces for computer aided geometric design[END_REF], a numerically stable method to evaluate Bézier curves. A Bézier curve is defined by a set of control points P 0 through P n , where n is called its order (n = 1 for linear, 2 for quadratic, etc.). The first and last control points are always the end points of the curve; however, the intermediate control points (if any) generally do not lie on the curve. Given points P 0 and P 1 , a linear Bézier curve B(t) is simply a straight line between those two points (see figure 9). The curve is given by :

B(t) = P 0 + t(P 1 -P 0 ) = (1 -t)P 0 + tP 1 , t ∈ [0, 1] (27)
With four points (P 0 , P 1 , P 2 , P 3 ), a Bézier curve of degree three can be built. The curve starts at P 0 going towards P 1 and arrives at P 3 coming from the direction of P 2 . Usually, it will not pass through P 1 or P 2 ; these points are only there to provide directional information (see figure 10). Properties • The polygon formed by connecting the Bézier points with lines, starting with P 0 and finishing with P n , is called the Bézier polygon (or control polygon). • The convex hull2 of the Bézier polygon contains the Bézier curve.

• The start (end) of the curve is tangent to the first (last) section of the Bézier polygon.

The explicit form of the curve is given by :

B(t) = (1 -t) 3 P 0 + 3(1 -t) 2 tP 1 + 3(1 -t)t 2 P 2 + t 3 P 3 , t ∈ [0, 1]. ( 28
)
B(t) = n ∑ i=0 b i,n (t)P i , t ∈ [0, 1] (29)
where the polynomials

b i,n (t) = n i t i (1 -t) n-i , i = 0, . . . n ( 30
)
are known as Bernstein basis polynomials of degree n. So, if there are many points, one has to manipulate polynomials with high degree. In order to circumvent this weak point one must use Basis-Splines.

Basis Splines

A B-spline [START_REF] De | A practical guide to splines[END_REF] is a spline function that has minimal support with respect to a given degree, smoothness, and domain partition. B-splines were investigated as early as the nineteenth century by Nikolai Lobachevsky. A fundamental theorem states that every spline function of a given degree, smoothness, and domain partition, can be uniquely represented as a linear combination of B-splines of that same degree and smoothness, and over that same partition. It is a powerful tool for generating curves with many control points, B stands for basis. A single B-spline can specify a long complicated curve and B-splines can be designed with sharp bends and even "corners". B-Spline interpolation is preferred over polynomial interpolation because the interpolation error can be made small even when using low degree polynomials for the spline. Furthermore, spline interpolation avoids the problem of Runge's phenomenon which occurs when interpolating between equidistant points with high degree polynomials.

Uniform B-Splines of Degree Zero

We consider a node vector T = {t 0 ,t 1 , ...,t n } with t 0 ≤ t 1 ≤, ..., ≤ t n and n points P i . One want to build a curve X 0 (t) such that :

X 0 (t i ) = P i (31) ⇒ X 0 (t) = P i ∀t ∈ [t i ,t i+1 ]. X 0 (t) = ∑ i B i,0 (t).P i ( 32
)
where

B i,0 (t) = 1 if t ∈ [t i ,t i+1 ] 0 elsewhere (33)
The shape of the X 0 (t) function in one dimension is given on figure 11.

Uniform B-Splines of Degree One

We are searching for a piecewise linear approximation X 1 (t) for which :

X 1 (t) = 1 - t -t i t i+1 -t i P i-1 + 1 - t -t i t i+1 -t i P i ∀t ∈ [t i ,t i+1 ] (34)
One can write X 1 (t) :

X 1 (t) = ∑ i B i,1 (t).P i ( 35
)
where

i t =4 t i+1 t i B (t)
i,1 (t) =      t-t i-1 t i -t i-1 if t ∈ [t i-1 ,t i ] t i+1 -t t i+1 -t i if t ∈ [t i ,t i+1 ] 0 elsewhere (36)
The shape of the X 1 (t) function in one dimension is given on figure 12.

Uniform B-Splines of Degree Three

Those B-Splines have been developed at Boeing in the 70s and represent one of the simplest and most useful cases of B-splines. Degree 3 B-Spline with n + 1 control points is given by :

X 3 (t) = n ∑ i=0 B i,3 (t).P i 3 ≤ t ≤ n + 1 ( 37
)
where

B i,3 (t) = 0 if t ≤ t i or t ≥ t i+4 . X 3 (t) = j ∑ i= j-3 P i .B i,3 (t) t ∈ [ j, j + 1], 3 ≤ j ≤ n ( 38
)
When a single control point P i is moved, only the portion of the curve X 3 (t) is changed (with t i < t < t i+4 ) insuring local controlproperty. The basis functions have the following properties :

• They are translates of each other i.e B i,3 (t) = B 0,3 (t -i) • They are piecewise degree three polynomial • Partition of unity ∑ i B i,3 (t) = 1 for 3 ≤ t ≤ n + 1 • The functions X i (t)
t i-1 i-1,1 B (t) 1 Fig. 12 Uniform B-Splines of Degree One B i-2,3 (t) = 1 h                    (t -t i-2 ) 3 if t ∈ [t i-2 ,t i-1 ] h 3 + 3h 2 (t -t i-1 ) + 3h(t -t i-1 ) 2 -3(t -t i-1 ) 3 if t ∈ [t i-1 ,t i ] h 3 + 3h 2 (t i+1 -t) + 3h(t i+1 -t) 2 -3(t i+1 -t) 3 if t ∈ [t i ,t i+1 ] (t i+2 -t) 3 if t ∈ [t i+1 ,t i+2 ] 0 otherwise ( 39
)
where h is the distance between two consecutive points.

Those basis functions are shown on figure 13.

Principal Component Analysis

When trajectories samples are available (from radar for instance), one can build a dedicated bases which will minimize the number of coefficient for trajectory reconstruction. Principal component analysis (PCA) is a mathematical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components. The number of principal components is less than or equal to the number of original variables. In the example presented in figure 14 a set of trajectories γ i (t), i = 1 . . . n are used to build K = 4 principal components (ψ k (t)) which can be used to reconstruct the initial trajectories.

B (t) 2, 3 3,3 B (t
γ i (t) = k=K ∑ k=1 a ik ψ k (t) (40)
When probability density functions of the coefficient a ik can be identified, one can use this trick to plug a stochastic optimization process which generates random coefficients in order to produce relevant random trajectory. More information about Functionnal PCA can be found in the reference [START_REF] Ramsay | Functional data analysis[END_REF].

Homotopy trajectory design

An easy way to build trajectory is to used reference trajectory (regular trajectories used by aircraft) and to compute a weighted sum of such reference trajectories to build a new one. If we consider two (or more) references trajectories joining the same origin destination pair (see figure 15) (past flown trajectories may be considered) :

γ 1 , γ 2 (41)
One can create a new trajectory γ α by using an homotopy :

γ α = (1 -α)γ 1 + αγ 2 (42)
In this example only one coefficient α has been used but one can extend this principle to several parameters.

All the models described in this section, may be used in an optimization process for which dimension reduction is needed. Depending on the targeted properties of the designed trajectories one must select the most adapted representations for the underlying application.

After having review algorithms based on optimization, the next section presents wave front propagation approaches.

Generalities

It is a well-known fact in physics that waves of high frequencies tend to propagate along the path whose traveling (or propagation) time is shortest (such minimum time path are called geodesic). The knowledge of the wave velocity at each point of space allows for the computation of wave fronts that are the set of all points reached by the wave at a given time, assuming it has started at a point source. The principle of the wavefront propagation algorithms is to simulate such physical propagation model in order to find optimal path based on a criterion that has to be optimized. As an example, congestion will be taken into account by velocity reduction, so that paths crossing congested areas will be penalized. Depending on the fact that the metric is or not isotropic (the later case being the one to be investigated when the wind is used into the criterion and has a non negligible velocity), two classes of algorithms are used.

Fast Marching Algorithms

The Fast Marching method, presented by Sethian [START_REF] Sethian | Fast marching methods[END_REF] is a part of the more general methods called Level Set [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. These techniques are designed to track the evolution of interfaces. The evolution of the wavefront can be compared to deform a curve or a surface from a partial differential equation. The Fast Marching method is used in the particular case where the wavefront speed is isotropic. It can still be applied if the anisotropy is low enough. For air traffic applications, this last assumption is valid in some areas of the airspace (of course it is not the case in the vicinity of jet streams).

In the particular case where the Fast Marching method is applicable, the calculus of the minimum time T to reach any points of the environment from the initial point is equivalent to solve the Eikonal equation of the form :

|∇T (x)| = 1 F(x) , F(x) > 0, and T (x initial ) = 0 (43)
where x ∈ R 2 represents the position in space, T ∈ R the minimum time and F : R → R the speed of propagation.

In free space, the wave speed F is equivalent to the aircraft speed. When we have forbidden areas, we force the propagation speed at zero in order to get a barrier value since the time to reach this point will be equal to infinity. Thereby, we have guaranteed avoidance property for those areas. For the congestion, the method is different, we want to penalize some areas where the congestion is high but we do not want to ban aircraft from driving through these areas. We just need to reduce the propagation speed. Thus, the time is increased proportionally to the congestion value, penalizing the crossing.

To design the optimal path between the arrival point and the departure point, we can then perform a gradient descent using the calculated values of T on the space, from the arrival point to the initial point. There is no risk to get stuck on a local minimum since the function T has only one optimum which is global.

The numerical resolution is like the graph search algorithms. However, in opposition to these graph search algorithms, the Fast Marching method is consistent since when the grid is refined, the obtained solution converges on the exact solution of the Eikonal equation [START_REF] Sethian | Fast marching methods[END_REF] that is a geodesic curve.

Ordered upwind algorithm

When wind is to be taken into account and has a non negligible speed with respect to the one of the aircraft, the propagation is no longer isotropic. The speed of the wavefront depends on the position and the directions of wind. A specific algorithm, called Ordered Upwind, has been developed to overcome this problem in [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: Theory and algorithms[END_REF], at the expense of a higher algorithmic complexity. Basically, an extra parameter, the anisotropy ratio, is considered: it is the ratio of the fastest to slowest propagation speed for each points. Given a point in space, the algorithm first considers the points on the current wavefront that are closest to it: it gives a time to travel when taking as propagation speed the slowest. Now, the other points to be considered are located no farther than the anistropy ratio times the minimal distance. By maintaining a list of potential points contributing to the information at the point of interest, the ordered upwind algorithm can still be implemented in a single pass.

In order to keep the efficiency of the Fast Marching algorithm, Petres proposed an extension of the algorithm Fast Marching in [START_REF] Ptrs | Path planning for autonomous underwater vehicles[END_REF], he assumed the field is smooth. He applied this extension to plan a path for autonomous underwater vehicles taking underwater currents into account. We propose here a similar extension to Petres's method of the Fast Marching method, our extension is specific to aircraft trajectories.

In the next section, we present an example of application of such wave propagation algorithm for air traffic management problem.

Light Propagation Algorithm

In geometric optics, light behavior is modeled using rays. Light emitted from a point is assumed to travel along such a ray through space. In an effort to explain the motion through space taken by rays as they pass through various media, Fermat (1601-1665) developed his principle of least action [START_REF] Douglas | Physics for scientists and engineers with modern physics[END_REF]:

The path of a light ray connecting two points is the one for which the time of transit, not the length, is a minimum.

In the framework of geometrical optics, light propagates through space respecting Descartes' Law:

Let n 1 be the index of refraction of medium 1 in which the incident ray propagates and n 2 that of the medium in which the refracted ray propagates (medium 2). We obtain:

n 1 . sin(θ 1 ) = n 2 . sin(θ 2 ) v 1 = c n 1 , v 2 = c n 2 , (44)
where v 1 ,v 2 are the speed values in media 1 and 2, θ 1 the angle of incidence, θ 2 the angle of the refracted ray and c the speed of light in a vacuum.

We can make several observations as a result of Fermat's principle :

• In a homogeneous medium, light rays are rectilinear. That is, within any medium where the index of refraction is constant, light travels in a straight line. • In an in-homogeneous medium, light rays follow smooth geodesic curves with minimum transit time.

Light therefore tends to avoid high index areas where rays are slowed down. Light reaches lowest speed for the highest encountered index.

Based on this principle of least action, we introduce an optimal path planning algorithm which computes smooth geodesic trajectories in environments with static or dynamic obstacles.

For a given mobile, we wish to find the shortest path between two points of R n , taking account of a given metric (time, distance etc.) whilst avoiding obstacles and respecting speed constraints. The trajectories produced must also respect a regularity criterion characterized by a maximum curvature value.

This algorithm mimics light propagation between a starting point towards a destination point, with obstacles modeled by high-index areas. By controlling the index landscape, it is possible to ensure that the computed trajectories meet the speed constraints and remain at a specified minimum distance from obstacles.

We begin by positioning a light source at the departure point (origin) which emits rays in a hemisphere oriented towards the destination point (this restriction prevents the generation of unrealistic trajectories which begin by turning 180 before turning back towards the destination). The path followed by the first beam to reach the arrival point corresponds to the geodesic trajectory we wish to obtain(see figure 16).

Optimal Control for Trajectory Generation

Optimal Trajectory Generation

In the physical space, a trajectory is occasionally represented as a four-dimensional flight path, following the tradition of air traffic control [START_REF] Chakravarty | Four-dimensional fuel-optimal guidance in the presence of winds[END_REF], with time as the fourth dimension, in addition to the normally used three-dimensional representation of a Destination Curent node Origine Fig. 16 Launching rays from the departure point. Geodesic computation (A * like algorithm or Triangle mesh algorithm).

path. Generating time-parameterized paths necessitates the incorporation of the aircraft dynamics and/or kinematics, which makes the problem much more difficult than simply finding a path that avoids obstacles in the physical three-dimensional space. Path-planning is a term commonly used in the robotics and artificial intelligence communities to refer to the problem of generating an obstacle-free path to be followed by a vehicle (robot, aircraft, vehicle, etc) in a two or three dimensional space containing obstacles [START_REF] Lavalle | Planning algorithms[END_REF].

Because the vehicle dynamics are not taken into account in these path-planning methods (the solution of which only considers the geometric constraints of the problem) it is often the case that the resulting path is infeasible, that is, it cannot be followed exactly or even closely by the vehicle. One way to ensure that the resulting paths correspond to feasible trajectories satisfying the vehicle dynamics, is to use optimal control theory. The objective of optimal control theory is to determine the control input(s) that will cause a process (i.e., the response of a dynamical system) to satisfy the physical constraints, while, at the same time, minimize (or maximize) some performance criterion. Feasibility of the trajectories is automatically ensured using this approach. The typical optimal control problem (OCP) can be stated as follows:

Given initial conditions x 0 , final conditions x f ∈ X , and an initial time t 0 ≥ 0, determine the final time t f > t 0 , the control input u(t) ∈ U and the corresponding state history x(t) for t ∈ [t 0 ,t f ] which minimize the cost function

J(x, u) = t f t 0 L(x(t), u(t)) dt, (45)
where x(t) and u(t) satisfy, for all t ∈ [t 0 ,t f ] the differential and algebraic constraints :

ẋ(t) -f (x(t), u(t)) = 0, (46)
C(x(t), u(t)) ≤ 0. ( 47
)
Optimal control has its roots in the theory of calculus of variations, which originated in the 17th century by Fermat, Newton, Liebniz, and the Bernoullis, and was subsequently further developed by Lagrange, Weirstrass, Legendre, Clebsch and Jacobi and others in the 18th and 19th centuries [START_REF] Ewing | Calculus of variations with applications[END_REF]. Calculus of variations deals with the problem of minimizing [START_REF] Jeffreys | Methods of mathematical physics[END_REF] subject to the simple differential equality constraint of the form ẋ(t)u(t) = 0, and is not able to handle more complicated differential equality constraints such as [START_REF] Lavalle | Planning algorithms[END_REF] or algebraic constraints such as [START_REF] Liu | Probabilistic aircraft mid-air conflict resolution using stochastic optimal control[END_REF]. It was not until the middle of the 20th century when the Soviet mathematician L. S. Pontryagin developed a complete theory that could handle constraints such as [START_REF] Lavalle | Planning algorithms[END_REF] and [START_REF] Liu | Probabilistic aircraft mid-air conflict resolution using stochastic optimal control[END_REF]. Simply put, Pontryagin's celebrated Maximum Principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF] states that the optimal control for the solution of the problem ( 45)-( 47) is given as the pointwise minimum of the so-called Hamiltonian function, that is,

u opt = argmin u∈U H(t, x, λ , u), ( 48
)
where

H(t, x, λ , u) = L(x, u) + λ T f (x, u)
is the Hamiltonian, and λ are the co-states, computed from

λ (t) = - ∂ H ∂ x (x(t), λ (t), u(t)). ( 49
)
subject to certain boundary (transversality) conditions on λ (t f ). Unfortunately, an analytic solution to the previous problem is difficult. The optimal control formulation of a trajectory optimization problem using Pontryagin's Maximum Principle (PMP) leads to a Two-point Boundary Value Problem (TBVP), or a Multi-point Boundary Value Problem (MBVP) when the optimal trajectory is composed of multiple phases. Numerical techniques such as shooting and multiple shooting methods can be applied to solve accurately these TBVP and MBVP, but their convergence is very sensitive to the choice of an initial guess for the solution. A software that solves the optimal control problem using this approach is BNDSCO [START_REF] Oberle | BNDSCO -A program for the numerical solution of optimal control problems[END_REF]. BNDSCO is an example of a class of numerical optimization methods which are often referred to as indirect methods. The term indirect reflects the fact that in these methods a solution is sought not by maximizing (or minimizing) the cost (45) but, rather, by computing potential optimizers by solving the corresponding necessary optimality conditions ( 48)- [START_REF] Berger | Differential geometry: Manifolds, curves and surfaces[END_REF].

In recent years, direct methods have become increasingly popular for solving trajectory optimization problems, the major reason being that direct methods do not require an analytic expression for the necessary conditions, which for complicated nonlinear dynamics can be intimidating. Moreover, direct methods do not need an initial guess for the co-states whose time histories are difficult to predict a priori. As mentioned earlier, direct methods, do not try to satisfy the necessary conditions of optimality from PMP, instead, they minimize directly [START_REF] Jeffreys | Methods of mathematical physics[END_REF] subject to ( 46)- [START_REF] Liu | Probabilistic aircraft mid-air conflict resolution using stochastic optimal control[END_REF].

The main idea behind direct methods is to discretize the states and controls of the original continuous-time optimal control problem in order to obtain a finitedimensional nonlinear programming problem (NLP). The solution of this NLP, which consists of discrete variables, is used to approximate the continuous control and state time histories. Typical direct methods are collocation methods, which discretize the ordinary differential equations (ODEs) of the problem using collocation or interpolation schemes [START_REF] Russell | A collocation method for boundary value problems[END_REF][START_REF] Weiss | The application of implicit Runge-Kutta and collocation methods to boundary value problems[END_REF][START_REF] Enright | Discrete approximation to optimal trajectories using direct transcription and nonlinear programming[END_REF][START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF]. They introduce the collocation conditions as NLP constraints together with the initial and terminal conditions. The so-called "pseudospectral" methods use orthogonal polynomials to choose the collocation points and are very efficient, exhibiting superlinear convergence when the solution is smooth. Numerical optimal control software packages that implement direct methods for the solution of OCPs include SOCS [START_REF] Betts | Sparse optimal control software SOCS[END_REF], RIOTS [START_REF] Schwartz | Theory and implementation of numerical methods based on runge-kutta integration for solving optimal control problems[END_REF], DIDO [START_REF] Ross | user's manual for DIDO: A MATLAB application package for solving optimal control problems[END_REF], PSOPT [START_REF] Becerra | Psopt optimal control solver user manual[END_REF], GPOPS [START_REF] Rao | User's manual for GPOPS version 4.x: A matlab package for software for solving multiple-phase optimal control problems using hp-adaptive pseudospectral methods[END_REF], MTOA [START_REF] Jain | Multiresolution-based direct trajectory optimization[END_REF] and DENMRA [START_REF] Zhao | Density functions for mesh refinement in numerical optimal control[END_REF] among many others. A recent survey of numerical optimal control techniques for trajectory optimization can be found in [START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF].

In all these direct methods, the convergence rate and the quality of solution depends on the grid used to discretize the equations, the cost and the problem constraints. Uniform or fixed grid methods tend to perform poorly, especially when the problem has several discontinuities or irregularities. Not surprisingly, adaptive grid methods have been developed to accurately capture any discontinuities or switchings in the state or control variables. The main idea behind all these adaptive grid methods is to use a high resolution (dense) grid only in the vicinity of control switches, constraint boundaries etc, and a coarse grid elsewhere. Examples of such adaptive gridding techniques for the solution of optimal control problems are [START_REF] Betts | Compensating for order variation in mesh refinement for direct transcription methods[END_REF][START_REF] Binder | Grid refinement in multiscale dynamic optimization[END_REF][START_REF] Schwartz | Theory and implementation of numerical methods based on runge-kutta integration for solving optimal control problems[END_REF][START_REF] Betts | Mesh refinement in direct transcription methods for optimal control[END_REF][START_REF] Jain | Trajectory optimization using multiresolution techniques[END_REF][START_REF] Gong | Spectral algorithm for pseudospectral methods in optimal control[END_REF].

A major issue with almost all current trajectory optimization solvers (direct or indirect) is the fact that their computational complexity is high and their convergence dependents strongly on the initial conditions, unless certain rather stringent convexity conditions hold. As a result, the solution of trajectory optimization problem in real-time is still elusive. A common line of attack for solving trajectory optimization problems in real time (or near real time) is to divide the problem into two phases: an offline phase and an online phase [START_REF] Seywald | Neighboring optimal control based feedback law for the advanced launch system[END_REF][START_REF] Lu | Regulation about time-varying trajectories: Precision entry guidance illustrated[END_REF][START_REF] Jardin | Neighboring optimal aircraft guidance in winds[END_REF][START_REF] Yan | Real-time computation of neighboring optimal control laws[END_REF]. The offline phase consists of solving the optimal control problem for various reference trajectories and storing these reference trajectories onboard for later online use. These reference trajectories are used to compute the actual trajectory online via a neighboring optimal feedback control strategy typically based on the linearized dynamics. Another strategy for computing near-optimal trajectories in real-time is to use a receding horizon (RH) approach [START_REF] Ohtsuka | Quasi-Newton-type continuation method for nonlinear receding horizon control[END_REF][START_REF] Bellingham | Stable receding horizon trajectory control for complex environment[END_REF][START_REF] Williams | Application of pseudospectral methods for receding horizon control[END_REF]. In a receding horizon approach a trajectory that optimizes the cost function over a period of time, called the planning horizon, is designed first. The trajectory is implemented over the shorter execution time and the optimization is performed again starting from the state that is reached at the end of the execution time. A third approach is to use a two-layer architecture, where first an acceptable (in terms of length, safety, etc) path is computed using common path-planning techniques, and then an optimal time-parameterization is imposed on this path to yield a feasible trajectory. As mentioned earlier such an approach needs to be carefully designed to ensure compatibility of the resulting path with the vehicle dynamics. However, when successful, such an approach is numerically very efficient and can be implemented in real-time with current computer hardware. Even if the resulting trajectory is not exactly feasible, it is often close to a feasible trajectory, or it can be made as such using smoothing techniques [START_REF]Stable receding horizon trajectory control for complex environment[END_REF]. As a result, alternatively, the final trajectory can be used as a good initial guess for a follow-up optimal trajectory generator. The next section summarizes this approach for applications related to aircraft maneuvering under strict time and fuel constraints. For more details the interested reader is referred to [START_REF]Time-optimal parameterization of geometric path for fixed-wing aircraft[END_REF][START_REF] Bakolas | Initial guess generation for aircraft landing trajectory optimization[END_REF][START_REF] Zhao | Efficient and robust aircraft landing trajectory optimization[END_REF].

Emergency Aircraft Trajectory Design

Aircraft maneuvering was one of the first areas where optimal control theory was used to generate optimal trajectories. Not surprisingly, traditionally, most work has been focused on military aircraft. Relevant references on this subject are too many to enumerate here. We just mention the work on fuel and range optimization studied in [START_REF] Grimm | Periodic control for minimum-fuel aircraft trajectories[END_REF][START_REF] Seywald | Range optimal trajectories for an aircraft flying in the vertical plane[END_REF][START_REF] Seywald | Long flight-time range-optimal aircraft trajectories[END_REF], and the minimum-time, three-dimensional aircraft trajectory optimization problem considered in [START_REF] Schultz | Three-dimensional trajectory optimization for aircraft[END_REF]. In the latter, an approximation of the aircraft dynamics using an energy state was used to reduce the dimension of the problem for better convergence. This type of model reduction technique is commonly used for aircraft trajectory optimization [START_REF]Control and dynamic systems: Advances in theory and applications[END_REF].

Some of these results have been extended to commercial airline operations. For example, the work of [START_REF] Burrows | Fuel-optimal aircraft trajectories with fixed arrival times[END_REF] deals with the problem of minimum-fuel trajectories with fixed time-of-arrival (TOA) for several civil aviation aircraft including B737, B747 and B767. Trajectory planning problems have also been studied in the context of air traffic management (ATM) and automation. Reference [START_REF] Jackson | Sensitivity of trajectory prediction in air traffic management[END_REF] performed a sensitivity analysis of trajectory prediction for ATM. The aircraft trajectory synthesis problem is studied in [START_REF] Slattery | Trajectory synthesis for air traffic automation[END_REF] to provide some basic tools for air traffic automation. Somewhat related is the recent work of Sridhar [START_REF] Sridhar | Aircraft trajectory optimization and contrails avoidance in the presence of winds[END_REF], in which he considered the generation of wind-optimal trajectories for cruising aircraft while avoiding the regions of airspace that facilitate persistent contrails formation. A shooting method was employed to solve the associated optimal control problem by minimizing a weighted sum of flight time, fuel consumption, and a term penalizing the contrail formation. The airspace avoidance problem has also been considered in Ref. [START_REF] Jacobson | Airspace constraints in aircraft emission trajectory optimization[END_REF]. In that reference, the avoidance of restricted airspace is formulated as a non-convex constrained trajectory optimization problem; it is claimed that with a feasible starting guess, the efficiency of the optimization algorithm is not degraded too much by the non-convex airspace constraints. Finally, some researchers have used ideas from stochastic optimal control to deal with issues related to the unpredictability of future trajectory, wind effects, presence of additional aircraft in the ATM airspace etc [START_REF] Liu | Probabilistic aircraft mid-air conflict resolution using stochastic optimal control[END_REF][START_REF] Tomlin | Conflict resolution for air traffic management: A study in multi-agent hybrid systems[END_REF].

In this work we deal with the problem of efficiently generating minimum-time and minimum-fuel (with fixed TOA) landing trajectories for commercial aircraft.

The former problem is of relevance in case of an on-board emergency where the pilot has to land the airplane quickly and safely to the closest airport or airfield; the latter problem is of interest for typical terminal ATC phase applications. Prior work in emergency landing includes the abort landing problem in the presence of windshear [START_REF] Bulirsch | Abort landing in the presence of windshear as a minimax optimal control problem Part I: Necessary conditions[END_REF][START_REF]Abort landing in the presence of windshear as a minimax optimal control problem Part II: Multiple shooting and homotopy[END_REF][START_REF] Miele | Optimal trajectories and guidance trajectories for aircraft flight through windshears[END_REF], and emergency landing during loss of thrust [START_REF] Strube | Dynamic flight guidance recalibration after actuator failure[END_REF]. The latter references generate feasible trajectories using segments of trajectories corresponding to selected trim condition maneuvers. The search results are however limited to those that can be generated by connecting trim state trajectory segments with stable transitions. Because the unstable flight conditions are not considered in the search, the algorithm cannot identify any feasible trajectories containing unstable flight modes. Furthermore, the path length is used as the search criterion, which is less appropriate when compared to flight time for emergency landing, or fuel consumption for normal flight. Related work includes the investigation of Atkins et al [START_REF] Atkins | Emergency flight planning applied to total loss of thrust[END_REF], where the problem of emergency landing due to the loss-of-thrust was studied using a hybrid approach. A two-step landing-site selection/trajectory generation process was adopted to generate safe emergency plans in real time under situations that require landing at an alternate airport. In the trajectory generation routine, a heuristic path planner was used to generate a three-dimensional trajectory connecting the current position of the aircraft to the runway, which consists of straight lines and circular arcs. This method is fast and simple. However, it is limited to conservative aircraft maneuvers (typically Dubins paths) in order to reduce the chance of obtaining an infeasible trajectory. As a result, the optimality of the generated trajectory could be unacceptable for emergency landing, and further research is necessary to reduce such a conservatism.

In our approach, we start with the assumption that the path to be followed by the aircraft is given. Note that this does not mean that the trajectory to be followed is given. A trajectory requires a time-parameterized path and it is, indeed, the main goal of this approach to provide such a time parameterization so as to meet feasibility along with certain optimality specifications. The assumption that the path is given is not as unusual or atypical as one may initially think. Commercial airliners during the terminal landing phase, are required to follow strict Air Traffic Control (ATC) rules, which guide the airplanes to follow "virtual" three-dimensional corridors all the way to the landing strip. Furthermore, since our approach leads to very fast computation of feasible trajectories, one can use the approach over new, locally modified paths repeatedly till a satisfactory path is found.

To this end, let a path in the three-dimensional space, parameterized by the path coordinate s, be given as follows:

x = x(s), y = y(s), z = z(s), where s ∈ [s 0 , s f ].
The main objective is to find a time-parameterization along the path, i.e., a function s(t), where t ∈ [0,t f ] such that the corresponding time-parameterized trajec- tory (x(s(t)),y(s(t)),z(s(t))) minimizes either the flight time t f (emergency landing case) or fuel (terminal landing operation). As shown in [START_REF] Zhao | Efficient and robust aircraft landing trajectory optimization[END_REF] all control (thrust, angle of attack, load factor, etc) constraints can be mapped into constraints involving the specific kinetic energy of the aircraft, E = v 2 /2 where v is the aircraft velocity of the form g w (s) ≤ E(s) ≤ g w (s),

for some path-dependent functions g w (s) and g w (s). The original problems therefore reduce to the following simplified problems:

For the Minimum-Time Problem we have min

T s f s 0 ds 2E(s) (50a) subject to E ′ (s) = T (s) m -D(E(s), s) -g sin γ(s), (50b)
g w (s) ≤ E(s) ≤ g w (s), (50c)
T min ≤ T (s) ≤ T max , (50d)
and for the Minimum-fuel Problem with fixed TOA we have min

T s f s 0 T (s) ds, (51a)
subject to

E ′ (s) = T (s) m -D(E(s), s) -g sin γ(s), (51b)
t ′ (s) = 1 2E(s) , (51c)
g w (s) ≤ E(s) ≤ g w (s) (51d)
T min (s) ≤ T (s) ≤ T max (s), (51e)
where D(E(s), s) is the drag, T is the thrust, γ is the flight-path angle, and where prime denotes differentiation with respect to path length s. The main advantage of these problem formulations is the dimensionality reduction of the problem that can be leveraged to solve both of these problems very efficiently and reliably. In fact, these OCPs are simple enough so that the optimal switching structure of the optimal solution can be unraveled using the necessary conditions from PMP. For the minimum-fuel problem the switching structure varies depending on the given TOA. However, for a given path and a fixed TOA, the structure is uniquely determined. This helps tremendously the convergence properties of the algorithm.

The overall architecture for optimal on-line trajectory generation is shown in Fig. 17. As shown in this figure, the method first generates a trajectory by assigning an optimal time parameterization along the path given by the geometric path planner via the solution of one of the previous two optimization problems. If the trajectory is feasible, then it is used as an initial guess for the numerical optimal control solver. Meanwhile, such a feasible trajectory is also stored as a back-up plan in case of the failure of the NLP solver. If the trajectory generated by the time-optimal path tracking method is not feasible, then the path is revised using the path smoothing method described in [START_REF]Stable receding horizon trajectory control for complex environment[END_REF], and the optimization is applied again to the smoothed path. The process is repeated until either the trajectory is feasible, or the maximum number of iterations is reached. If no feasible trajectory can be obtained after reaching the iteration limit, the infeasible trajectory is passed to the numerical optimal control algorithm, which makes a last attempt to produce a feasible trajectory. If this last attempt is not successful, then it does not exist a feasible trajectory that solves the problem.

6 Path planning techniques via natural language processing and mathematical programming: A paradigm for future aircraft trajectory management

Aircraft trajectory planning has reached enough maturity to shift the trajectory planning problem from the mathematical optimization of the aircraft trajectory to the automated parsing and understanding of desired trajectory goals, followed by their re-formulation in terms of a mathematical optimization program. To a large extent, we propose a possible evolution of the Flight Management System currently used on all transport aircraft to become a full-fledged autonomous logic that may also be used on unmanned aerial vehicles as an alternative to the current -and bulky-Remotely-Piloted Vehicle (RPV) paradigm. What follows is a direct extension of work originally presented in [START_REF] Schouwenaars | Linear programming and language processing for human/unmanned-aerial-vehicle team missions[END_REF].

The overall proposed architecture is shown in Fig. 18. It consists of several nested feedback loops. At the operator level, the information is presented to the operator in the form of sentences expressed in natural language (eg that used by air traffic control phraseology). At the level of trajectory planning automation, the information is presented as a mix of continuous parameters (aircraft position and speed), and discrete parameters describing mission status (completed tasks, tasks remaining to be completed). The operator formulates the vehicles's goals (eg flight plans) using Fig. 18 Basic Autonomous Mission Management Loop for future Aircraft and UAVs natural language. A natural language interpreter and task scheduler transforms the operator's requirements into tractable mathematical optimization programs that may be executed by the the vehicle through its flight control computer. The vehicle's innermost dynamics (that consist of raw vehicle dynamics and stability augmentation system), although critical to vehicle stability, are not shown.

Natural language parsing and generation

The main goal of using a Natural Language Interface (NLI) for interacting with a computer-based system is to minimize the workload on the operator. Using normal English sentence commands and reports indeed allows an operator to communicate efficiently and effectively with an aircraft, as if it were a human pilot. The NLI module that we have developed for demonstration purposes consists of two major components. The first one takes sentence commands from the operator (presumably an air traffic controller) and turns them into a formally coded command that looks like the formulation of an optimization problem. The second component takes a coded command set from the aircraft and generates natural language responses for the air traffic controller to interpret. A sample dialog between the human operator and the machine could be as follows: Controller: Flight AA1234, this is Air Traffic Control. Aircraft: Go ahead, Air Traffic Control. Controller: Add new waypoint. Proceed to waypoint Echo-Charlie 5 in minimum time. and wait for further instructions after the task is completed Aircraft: Roger. Acknowledge task information -proceeding to waypoint Echo-Charlie 5. Controller: AA1234, out.

The NLI module analyzes the natural sentences produced by the air traffic controller using parsing, which is the process of converting an input sentence, e.g. "Proceed to waypoint Echo-Charlie 5 in minimum time," into a formal representation. The latter is typically a tree structure which can in turn be translated into an explicit formal command. In our system, parsing consists of first applying entity extraction to all the individual concepts (e.g. "Flight AA1234" or "Echo-Charlie 5") and then combining these concepts through cascades of finite-state transducers using techniques derived from those described in [START_REF] Roche | Parsing with finite-state transducers, Finite-State Language Processing[END_REF]. While natural language processing represents our progress so far, it is easy to imagine that it could now be completed by a voice recognition device to further ease the level of communication between controller (or operator) and aircraft.

Task Scheduling and Communications Interfacing

The task scheduling and communications processing components are designed to centralize all of the aircraft mission processing in one module. Together with the Natural Language Interface, it provides flexibility for an operator to insert and change flight plan during the flight. The aircraft software keeps track of the flight tasks, waypoint locations and known obstacles to pass on to the guidance algorithm. The communications processing component provides the air trafic controller or operator with the authority to send commands and receive status updates, threat or obstacle avoidance information and acknowledgement messages. It also provides remote pilots monitoring the flight with the ability to override the guidance system in the event of an emergency or error. The system sends threat and override information to the air traffic controller before any status or update information in an effort to send the most important data relevant to the demonstration before any auxiliary information. Input/Output data are processed every 1 Hz frame before the task planner and guidance step to ensure that the most up-to-date information is used by the aircraft trajectory planner. The task scheduling component operates like a Flight Management System and allows the aircraft operator or the air traffic manager to enter a flight plan using a pre-defined list or as programmed during a mission. Many additional features may be added to such a task scheduler, such as orders to follow loiter patterns, 'take me home' or low-emission approaches functionalities. In addition, he or she has the option of providing (in real-time via the NLI) the optimization metric used by the trajectory generation algorithm (i.e. minimum time, minimum fuel, or the amount of time to finish the flight). Next, the operator can either give the aircraft a new plan or change the current plan it is performing. A "New Plan" command is added to the end of the aircraft task list and is executed after all of the tasks currently scheduled have been completed. A "Change Plan" command, on the other hand, modifies the current task performed by the aircraft. Once a task is completed, it is removed from the list. After each of these actions, an acknowledgement is sent to the air traffic controller and the updated task information is included in the data sent to the Trajectory Generation Module.

Trajectory planning

After the Natural Language Interface and Flight Planning and Scheduling components have converted the flight plan into a series of tasks for the aircraft to perform, the Trajectory Generation Module guides the vehicle from one task to the next, i.e. from an initial state to a desired one, through an obstacle field while optimizing a certain objective. The latter can be to minimize time, fuel or a combination of both. Much of the functionality described below becomes increasingly available in today's avionics systems, and also include such real-world factors as weather and wind conditions. For our purposes, 2D scenarios were considered in which specialuse airspace and other no-fly zones are viewed as obstacles and detected while the flight proceeds. The environment is always fully characterized inside a certain detection region D around the aircraft. The resulting formulation can, however, be easily generalized to account for any detection shape, such as a radar cone, and for unknown areas within that shape. Since trajectories must be dynamically feasible, the aircraft dynamics and kinematics should be accounted for in the planning problem. For optimization purposes, the vehicle is characterized by a discrete time, linear state space model (A, B) in an inertial 2D coordinate frame (east-north). As such, the state vector x consists of the east-north position (x, y) and corresponding inertial velocity ( ẋ, ẏ). Depending on the particular model, the input vector u is an inertial acceleration or reference velocity vector. In both cases, however, combined with additional linear inequalities in x and u, the state space model must capture the closed-loop dynamics that result from augmenting the aircraft with a waypoint tracking controller. Since the environment is only partially-known and further explored in real-time, a receding horizon planning strategy is used to guide the vehicle towards the desired destination.

The destination is denoted by x f and is an ingress/egress state of a waypoint with a corresponding inertial velocity vector. At each time step, a partial trajectory from the current state towards the goal is computed by solving the trajectory optimization problem over a limited horizon of length T. Because of the computation delay, the initial state x 0 = (x 0 , y 0 , ẋ0 , ẏ0 ) in the optimization problem should be an estimate x estim of the position and inertial velocity of the aircraft when the plan is actually implemented. The solution to the optimization problem provides a sequence of waypoints (x i , y i ) and corresponding inertial reference velocities ( ẋi , ẏi ) to the aircraft for the next T time steps. Typically, however, only the first waypoint and reference velocity of this sequence are given to the waypoint follower, and the process is repeated at the next time step. As such, new information about the state of the vehicle and the environment can be taken into account at each time step. By introducing a cost function J T over the T time steps, the general trajectory optimization problem can be formulated as to

Minimize x i ,u i J T = T ∑ i=1 f i (x i , u i , x f ) + f T (x T , x f ) Subject to                x i+1 = A x i + Bu i , i = 0, . . . , T -1 x 0 = x estim x i ∈ X 0 , i = 1, . . . , T u i ∈ U 0 , i = 0, . . . , T -1 (x i , y i ) ∈ D 0 , i = 1, . . . , T (x i , y i )\ ∈ O 0 , i = 1, . . . , T
The objective function consists of stage costs f i (x i , u i , x f ) corresponding to each time step i, and a terminal cost term f T (x T , x f ) that accounts for an estimate of the cost-to-go from the last state x T in the planning horizon to the goal state x f . The sets X 0 and U 0 represent the (possibly non-convex) constraints on the vehicle dynamics and kinematics,such as bounds on velocity, acceleration and turn rate. Here, the 0subscript denotes the fact that these constraints can be dependent on the initial state. Lastly, the expressions (x i , y i ) ∈ D 0 and (x i , y i )\ ∈ O 0 capture the requirement that the planned trajectory points should lie inside the known region D 0 , but outside the obstacles O 0 as given at the current time step i = 0. Note that they are assumed to hold for x 0 ; if not, the trajectory optimization problem would be infeasible from the start. As demonstrated in [START_REF] Schouwenaars | Receding horizon path planning with implicit safety guarantees[END_REF], however, despite the detection region and avoidance constraints, the above receding horizon strategy has no safety guarantees regarding avoidance of obstacles in the future. Namely, the algorithm may fail to provide a solution in future time steps due to obstacles that are located beyond the surveillance and planning radius of the vehicle. For instance, when the planning horizon is too short and the maximum turn rate relatively small, the aircraft might approach a no-fly zone too closely before accounting for it in the trajectory planning problem. As a result, it might not be able to turn away in time, which translates into the optimization problem becoming infeasible at a future receding horizon iteration. In [START_REF] Schouwenaars | Receding horizon path planning with implicit safety guarantees[END_REF], a safe receding horizon scheme was therefore proposed based on maintaining a known feasible trajectory from the final state x T in the current planning horizon towards an obstacle-free holding pattern. The latter must lie in the region D 0 of the environment that is fully characterized at the current time step, and is computed and updated online. Assuming that the planned trajectories can be accurately tracked, at each time step, the remaining part of the previous plan together with the holding pattern can then always serve as an a priori safe backup or "rescue" plan. In practice, we have found that formulating the problem of finding the nominal and rescue trajectories optimization as mixed-integer programs (MIP) works very well in practice, though such choices are not mandatory and may be replaced by the other techniques discussed in this paper.

Application to Air Traffic Management

Aircraft trajectory design can be applied in many areas of air traffic managements like strategic planning, pre-tactical planning, tactical planning, SID-STAR design, emergency trajectory design, etc... In the following three Air Traffic Management applications are presented for which trajectory design is critical. The first one presents the strategic aircraft trajectory planning (optimization approach), the second one the pre-tactical and tactical planning by light propagation algorithm (wave propagation approach) finally the last one describes Emergency trajectory design by optimal control (optimal control approach).

Strategic aircraft trajectory planning

The aim of strategic planning is to reduce airspace congestion by modifying takeoff slots and routes for a set of flights on a country-wide scale.We thus considered a sectorized airspace and a set of flight plans for a given day. For each flight, we defined a set of alternative routes and a set of possible takeoff slots. Our goal was to produce an optimal assignment to minimize airspace congestion. Our problem can be presented as follows:

• We consider all of the flight plans associated with the airspace of a country.

• For each airplane k, we suppose that the following elements are known:

-a set of possible routes (+ associated costs) -a set of possible takeoff times (+ associated costs) -the set of flights connected with flight k at departure and arrival points (hub phenomenon)

We wish to obtain the configuration which allows us to reduce airspace congestion, minimize cost and respect connection plans.

The architecture of the chosen approach (see figure 19) is made up of a traffic simulator and a genetic algorithm (GA).

The simulator is used in pre-treatment to obtain the data needed for optimization. It receives flight plans (airplane identification, airplane type, origin-destination, departure time, flight profile, route as a set of beacons) and constructs airplane trajectories, while saving the necessary data (sector input and output, flight position and direction every minute, etc.).

Genetic algorithms simulate the process of natural selection in a hostile environment linked to the problem under consideration [START_REF] Davis | Handbook of genetic algorithms[END_REF]. They use a vocabulary similar to that found in natural genetics, without neglecting the fact that the underlying principles involved in the two domains are considerably more complex in their natural context. We thus refer to individuals in a population, and often individuals will be reduced to a single chromosome. These chromosomes themselves are made up of genes which contain the hereditary characteristics of the individual. We also use principles of selection, crossing (crossover), mutation etc.

In the context of optimization, each individual represents a point in the state space to which we associate the value of the criterion to optimize. We then randomly generate a population of individuals from which the genetic algorithm aims to select the best specimens while ensuring efficient exploration of the state space.

The optimization algorithm then manipulates possible plannings, represented by disturbances, in terms of routes and slots, to the initial flight plans.

The optimization generates plannings which are increasingly efficient in terms of minimizing congestion. The best planning is then put forward as the result of the optimization process.

This approach has been successfully applied for removing conflict in European airspace with 32000 aircraft. The purpose of this research is to design a gate-to-gate More information about this work may be found in [START_REF] Chaimatanan | Conflict free strategic planning[END_REF]. The objective of the proposed method, based on an analogy with optics, is to find an optimal 4D trajectory for each aircraft which avoids congestion or conflicts and minimizes a criterion based on a local metric. Congestion and the protection zone (volume surrounding the aircraft where no other aircraft may enter) of other aircraft will be modeled as high-index areas. Our light propagation algorithm (LPA) is designed from a particular aircraft point of view. It is assumed that the aircraft knows the surrounding aircraft trajectories (the set of trajectories of the other aircraft is a given input of the algorithm). Some examples of geodesic trajectories computation are given on the following figure 22 for which the red areas contain high index values. This algorithm has also been used at tactical phase in order to produce conflict free trajectories. To generate a trajectory, we use a wavefront propagation algorithm in 2D + time (in order to conform to operational practice, we carry out resolution in terms of heading, but with a detection method which takes account of altitude) with temporal sampling (the wave is propagated with a time step dt).

Real flight plans of 12 August 2008 with about 8000 flights have been used to produce reference trajectories. The initial trajectories (before conflict resolution) induce a total number of around 4000 conflicts. Based on a sliding time window, trajectory segments are extracted in order to built conflict cluster for which LPA is applied. To address the conflict resolution, LPA is sequentially applied to aircraft. We assign a trajectory to the first aircraft disregarding the other aircraft (without considering any constraints). Then, LPA looks for a trajectory for the subsequent aircraft by considering the trajectory of the first aircraft as a constraint, and so on, up to the m th aircraft which considers the m -1 previous aircraft trajectories as constraints. In practice, some operational criteria may also be used in order to select a specific sequence (for instance: first-come first-served rule, some aircraft may have higher priority, trajectory length, etc.). The time window is then shifted and the process is applied again.

The algorithm nearly solves all conflicts, with only 28 situations for which conflict-free trajectories have not been found. However, these situations correspond to some aircraft being already in conflict at the beginning of the simulation, for instance at their starting point. Only 1501 trajectories have been modified to reach such a conflict-free planning. In many cases, the new computed trajectories are shorter than the initial ones (those that follow waypoints), due to the fact that LPA is searching for the shortest path trajectories and proposes direct routes when possible. More information about this work may be found in [START_REF] Dougui | A light-propagation model for aircraft trajectory planning[END_REF].

Emergency trajectory design by optimal control

This work presents an application of optimal control for computing backup trajectories in case of emergency. In such a situation, one must be able to design backup landing trajectories for aircraft that have lost having loose their thrust. The proposed algorithm begin to find a flyable path to avoid obstacles (Dubins paths generation with continuous descent), then find a feasible trajectory to follow along this path. This requires the solution of optimal time parametrization (or velocity generation) problem. The latter is a one-dimensional optimal control problem that can be solved very efficiently.

Two cases are presented in the following, the first one is link to the Swissair 111 crash and the second one to the US Air 1549 crash.

On Wednesday, 2 September 1998, the Swiss Air 111 flight crashed into the Atlantic Ocean southwest of Halifax International Airport (due to fire on board).

The Swissair 111 trajectory is given on figure 23 The algorithm has been applied to this case and has produced several backup trajectories bringing back the aircraft at Halifax airport as it can be shown on figure 24. This algorithm has also been applied to the US Air 1549 case. On January 15, 2009, the US Air 1549 flight struck a flock of Canada Geese during its initial climb out, lost engine power, and ditched in the Hudson River off midtown Manhattan. Again, the algorithm found 4 backup trajectories bringing back safely the aircraft the Laguardia airport (see figure 25).

Conclusion

This survey has shown several approaches for trajectory modeling. As it has been mentioned, trajectories are belonging to infinite dimension space for which dimension reduction has to be implemented. Using trajectory samples vectors is really re- dundant and inefficient. After having presented some definitions and some features of aircraft trajectories, some dimension reductions methods have been presented in order to be included in an optimization process. Interpolation and approximation technique have been presented and some information have also been given about Principal Component Analysis and Homotopy. The next section has presented wave front propagation approaches and gives some details on Fast Marching Algorithms, Ordered upwind algorithms and the light propagation algorithm (LPA). The fourth part has focused on methods coming automatic control for which vehicle dynamics are explicitly included in the models. Some applications to air traffic management have also been presented. The sixth part has presented an original path planning

  , b) the point t is said to be singular. A parametrized curve of class C p (or more concisely a C p curve) will be a C p mapping from an open time interval (a, b) to the state space E with no singular points. Any C 1 curve can be parametrized by arc length. Let γ(a, b) → E be such a curve. Defining the mapping s(a, b) → (0, l(γ)) by:

Fig. 1 .

 1 Fig.1.The optimization process control the X vector in order to build a trajectory γ for evaluation.

Fig. 2

 2 Fig. 2 Trajectory defined by four way points connected by straight lines.

Fig. 3

 3 Fig. 3 L n (x) is represented by the black curve. The others curves are the polynomials l i (x).

Fig. 4

 4 Fig. 4 Lagrange interpolation result for a set of aligned points.

Fig. 5

 5 Fig. 5 Piecewise linear interpolation.

Fig. 8

 8 Fig. 8 Cubic Spline Interpolation.

Fig. 9 Fig. 10 Cubic

 910 Fig. 9 Bézier Curve with 2 points.

Fig. 11

 11 Fig. 11 Uniform B-Splines of Degree Zero

B

Fig. 13

 13 Fig. 13 Order 3 basis function

Fig. 15

 15 Fig.[START_REF] Bulirsch | Abort landing in the presence of windshear as a minimax optimal control problem Part I: Necessary conditions[END_REF] One new trajectory is built by using a weighted sum of two reference trajectories.

Fig. 17

 17 Fig. 17 Schematic of landing trajectory optimization.

Fig. 19

 19 Fig. 19 Optimization structure

  Figure 20 shows the European traffic with interaction in red.

Fig. 20

 20 Fig. 20 Direct route traffic over Europe

Figure 21

 21 Figure 21 represents the conflict free traffic after optimization.More information about this work may be found in[START_REF] Chaimatanan | Conflict free strategic planning[END_REF].

Fig. 21

 21 Fig. 21 Traffic picture after optimization

Fig. 22

 22 Fig. 22 Example of optimal geodesic curves. The red color indicate high index areas.

Fig. 23

 23 Fig. 23 Swissair 111 trajectory reconstruction

Fig. 24

 24 Fig. 24 Swissair 111 case backup trajectories generation in green. The blue trajectory represent the real flown trajectory.

Fig. 25 US

 25 Fig. 25 US Air 1549 backup trajectories in green . The blue trajectory represent the real flown trajectory ending in Hudson river.

  are of degree 3 for any set of control points

				X (t) 1
					Pi
	0	1	2	3	i t =4	5	6	7	8	9
					t i	t i+1

The convex hull or convex envelope of a set X of points in the Euclidean plane or Euclidean space is the smallest convex set that contains X.

techniques mixing natural language processing and mathematical programming. Finally, some air traffic management applications have been presented for realistic cases.