Florence Nicol
email: nicol@recherche.enac.fr

Stephane Puechmorel
email: stephane.puechmorel@enac.fr

Information Geometry for Safety Data Analysis

Keywords: curve clustering, probability distribution estimation, functional statistics, safety data analysis

Air traffic system generates huge amounts of data, most of them being only partly exploited. One important issue is the functional nature of the samples that consist mainly in curves. Analyzing this kind of highly structured data requires dedicated tools. Most of the time, the functions considered are expanded on a truncated Hilbert basis then usual multivariate statistics tools are applied on them. When additional constraints are put on the curves, like in applications related to air traffic where operational considerations are to be taken into account, this approach is no longer applicable. The present article surveys some recent results obtained by the authors using a new framework where curves are represented as regularized currents. The question of high dimensional data with intrinsic redundancy is also discussed and a possible extension of the already obtained results is sketched.

I. INTRODUCTION

The air traffic system outputs huge volumes of data coming from surveillance systems, like Radar or ADS-B, and from monitoring and maintenance sources. Although very different in sampling rates and samples dimensions, they share in common the fact that almost all the time, data can be represented as functions.

As an example, clustering aircraft trajectories is a central problem in Air Traffic Management (ATM). It arises naturally when designing arrival and departure procedures and when some special constraints are to be considered: noise exposure, pollutants emission. For this application, one must first identify the major flows and then assess the statistical structure of the discrepancy of the real flight paths around them.

In the en-route phase, clustering is applied in order to extract mean flight paths, which are further used to optimally design the airspace elements (sectors and airways). A special instance of this problem is the automatic generation of safe and efficient trajectories, but in such a way that the resulting flight paths are still manageable by human operators. Several algorithms can automatically produce conflict free trajectories but generally fail to issue planning schemes that may be used by air traffic controller in the event of a system failure. Clustering may be a solution to this problem using a two phase procedure: in a first step, a set of flight plans or intended reference business trajectories is submitted to the optimizer which outputs a conflict free planning. In a second step, a clustering algorithms extracts medial lines from this trajectories set, which form a dynamic route system. If it is simple enough, a controller may use it in a manual control context.

Finally, a very important applicative field of trajectory analysis is related to safety data. This subject covers in fact many different applications, depending on the data source. When only radar tracks are used, only a small subset of the aircraft states can be inferred. However, it is still possible to get extract valuable information: runway adherence condition is a worked example of that. On-board data logging facilities give access to an almost complete picture of the aircraft state: the main issue when trying to analyze it is the volume and the internal redundancy of the samples. Tools dealing with such a situation are yet to be designed. A work in progress tries to pave the way to future systems, where almost independent features will be identified.

In all the above situations, a common denominator is the need for algorithms that take into account the functional nature of the samples. After a brief description of previous related works in general statistics and on the targeted air traffic application, a new framework based on densities associated with curve systems will be introduced. It relies on currents that can be seen as an extension of the delta distribution. An application to major flows identification will be presented. In a second part, a lie group modeling will allow to take into account flows orientations and velocities. Finally, this work will be further extended to flag manifolds, that may be viewed as an abstract model for independent subspaces analysis. While still in early stage of research, it represents the first attempt to treat the issue of functional data with values in a structured, high dimensional state space.

II. PREVIOUS RELATED WORK

The field of functional data is quite recent and active. The results collected in [START_REF] Ramsay | Functional Data Analysis[END_REF], [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF] give a picture of the multivariate algorithms that have been transposed to the functional framework.

The most obvious way of dealing with such data is simply to sample at evenly spaced times and to collect the values obtained in a finite dimensional vector. This procedure was used in a study conducted by the Mitre corporation on behalf of the Federal Aviation Authority (FAA) [START_REF] Enriquez | Identifying temporally persistent flows in the terminal airspace via spectral clustering[END_REF]. Due to the very high dimension of the vectors considered, clever numerical workarounds must be found. In the Mitre study, random projections reduce the dimensionality, then a classical PCA is performed. The huge computational cost of the required singular values decomposition is thus alleviated, allowing use on real recorded traffic over several months.

The most important limitation of the plain sampling approach is that the shape of the trajectories is no taken into account when applying the clustering procedure. Furthermore, there is no simple mean to put a constraint on the mean trajectory produced in each cluster: curvature may be quite arbitrary even if samples individually comply with flight dynamics.

Another approach is taken in [START_REF] Mahrsi | Graph-based approaches to clustering network-constrained trajectory data[END_REF], where an underlying graph structure is assumed. It is a variation of the original work described in [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. It is well adapted to road traffic as vehicles are bound to follow predetermined road segments. For air traffic applications, it may be of interest for investigating present situations, using the airways and beacons as a structure graph, but will mis-classify aircraft following direct routes which is a quite common situation, and is unable to work on an unknown airspace organization. This last limitation prevents its use for airspace redesign applications or in the context of trajectory based operations. Varying the similarity measure will put the emphasis on different aspects of the trajectories [START_REF] Kim | Spatial and temporal characterization of travel patterns in a traffic network using vehicle trajectories[END_REF],without removing the aforementioned drawbacks An interesting vector field based algorithm is presented in [START_REF]Vector field k -means: Clustering trajectories by fitting multiple vector fields[END_REF]. An salient feature is the ability to distinguish between close trajectories with opposite orientations. It falls within the frame of landmarks based algorithms, first introduced for shape analysis by Kendall [START_REF] Kendall | Shape manifolds, Procrustean metrics, and complex projective spaces[END_REF].

Due to the functional nature of trajectories, that are basically mappings defined on a time interval, some authors resort to techniques based on times series, as surveyed in [START_REF] Liao | Clustering of time series data -a survey[END_REF], [START_REF] Rani | Recent techniques of clustering of time series data: A survey[END_REF]. The algorithms pertaining to this category are based on sequences, possibly in conjunction with dynamic time warping [START_REF] Meesrikamolkul | Shape-based clustering for time series data[END_REF].

Finally, as mentioned at the beginning of the section, one can assume that data come from an unknown underlying stochastic process whose sample paths belong to a given Hilbert space with countable basis. Using truncated expansions of this basis, functional data revert for implementation to usual multivariate algorithms. For the same reason, modelbased clustering may be used in the context of functional data even if no notion of probability density exists in the original infinite dimensional Hilbert space [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF]. A nice example of a model-based approach working on functional data is funHDDC [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF].

III. CURVE SYSTEMS

When dealing with aircraft trajectories, some specific characteristics must be taken into account. First of all, flight paths consist mainly in straight segments connected by arcs of circles, with transitions that may be assumed smooth up to at least the second derivative. This last property comes from the fact that pilot's actions result in changes on aerodynamic forces and torques and a straightforward application of the equations of motion. When dealing with sampled trajectories, this induces a huge level of redundancy within the data, the relevant information being concentrated around the transitions. It is well known from numerical analysis that it renders covariance and related matrices illconditioned, inducing instabilities in the algorithms. From a practical standpoint, random pre-conditioners will greatly improve the situation, but render the interpretation of the final results difficult. At the same time, it is highly desirable that trajectories produced by algorithms be of smallest possible curvature, so that they may represent potential flight paths.

A second issue comes from the fact that trajectories of interest are defined as mapping from a time interval [a, b] to R 3 which is not the usual setting for functional data statistics: most of the work is dedicated to real valued mappings and not to vector ones. A simple approach will be to assume independence between coordinates, so that the problem falls within the standard case. However, even with this simplifying hypothesis, vertical dimension must be treated in a special way as both the separation norms and the aircraft maneuverability are different from those in the horizontal plane. This issues becomes even more important when dealing with on-board recorded flight data as the dimension of the state space may be in the order of several hundreds, with a strong correlation between coordinates.

Within the frame of functional data statistics, a clever choice of the hilbert state space may yield a solution to both problems. As a matter of fact, one can consider Sobolev spaces of the form:

H L = {f | Lf (x) 2 dx < +∞} (1)
where L is a differential operator with 0 kernel. Taking for example L = Id + AD xx , with A > 0 a tuning parameter, allows the acceleration of the curve to be considered. Suitable hilbert basis in such a case can be obtained using green functions of the L * L operator, with L * the adjoint of L, yielding the so-called spline basis [START_REF] Wahba | Spline Models for Observational Data[END_REF]. Such Sobolev spaces are reproducing kernel hilbert spaces (RKHS), with kernel the green function, which allows computation saves when evaluating inner products or distances. This relationship is exposed in detail in the above reference. The spline basis has been proved asymptotically optimal under special choices of L, fully justifying its use when working on truncated expansion.

When no obvious choice of L can be made or when the green function cannot be computed, the method is less straightforward to apply. Furthermore, the nice optimality property is generally lost when working with general basis: the truncation threshold may become so high that the method will no compete with plain sampling. Finally, geometric features of curves, like curvature, cannot be easily formulated as a Sobolev norm, and may not be preserved in the algorithms outputs.

A different approach is considered in this work, building upon geometry of shapes [START_REF] Michor | Riemannian geometries on spaces of plane curves[END_REF] and information geometry [START_REF] Puechmorel | Entropy minimizing curves with application to automated flight path design[END_REF]. Curves will be modeled as smooth currents [START_REF] Griffiths | Principles of Algebraic Geometry[END_REF], then regularized using a kernel functions. Given a complete set of trajectories, it will give rise to a density, that may enter an optimization process. As one of the most important requirement about flight paths is that they must be as straight as possible, an entropy criterion was selected. In the next section, the contributions coming from [START_REF] Puechmorel | Entropy minimizing curves with application to automated flight path design[END_REF], [START_REF] Puechmorel | Geometry of curves with application to aircraft trajectory analysis[END_REF], [START_REF] Iaria | Unsupervised Aircraft Trajectories Clustering: a Minimum Entropy Approach[END_REF] are summarized.

IV. THE ENTROPY OF A SYSTEM OF CURVES

Let a set γ 1 , . . . γ N of smooth curves with value in a bounded open subset Ω of finite dimensional vector space E be given. If U k , k = 1 . . . P is a partition of Ω, then one can construct a density estimator for the curve system as:

d k = λ -1 N i=1 1 0 1 U k (γ i (t)) dt (2)
where the normalizing constant λ is obtained as:

λ = P k=1 N i=1 1 0 1 U k (γ i (t)) dt (3) = N i=1 1 0 P k=1 1 U k (γ i (t)) dt (4)
and since U k , k = 1 . . . P is a partition:

λ = N i=1 1 0 dt = N (5)
Density can be viewed as an empirical probability distribution with the U k considered as bins in an histogram.

Proceeding the same way as in non-parametric density estimation [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], a kernel smoother is used in place of the characteristic functions 1 U k :

d : x → N i=1 1 0 K ( x -γ i (t) ) dt N i=1 Ω 1 0 K ( x -γ i (t) ) dtdx (6)
it comes:

Ω K ( x -γ i (t) ) dx = R 2 K ( x ) dx provided that Ω contains the set: {x ∈ R 2 , inf i=1...N,t∈[0,1] x -γ i (t) ≤ A}
where the interval [-A, A] contains the support of K. When the kernel has unit integral, the expression of the density simplifies to:

d : x → N -1 N i=1 1 0 K ( x -γ i (t) ) dt (7)
A common choice for K is the Epanechnikov function:

K : x → 1 -x 2 1 [-1,1] (x)
that is presented here in its un-normalized version. The density map of one day of traffic over France is considered and pictured on figure 2: Unfortunately, this simple procedure suffers a severe flaw: it tends to overemphasis areas when the aircraft are slow, as they will represent a longer time occupancy for the same size. A change is made so as to overcome this issue, yielding:

d : x → N i=1 1 0 K ( x -γ i (t) ) γ i (t) dt N i=1 l i (8)
with l i the length of curve i. A unit integral kernel is assumed in the above expression. The resulting density map is more satisfying, as shown in figure 3: It turns out that the expression can be interpreted as a sum of regularized smooth currents [START_REF] Griffiths | Principles of Algebraic Geometry[END_REF], with the curves γ i considered as 0currents and K the smoothing kernel. Other choices of K may be more natural in this context: Gaussian kernels will be a good starting point, but are not compact (this is not in issue in real-world applications as the truncation error is generally negligible). Varying the bandwidth of K will change the amount of smoothing applied to the original curves: very low ones tend to reproduce accurately the trajectories while larger ones will spread the density over Ω. All the theoretical properties of smooth currents apply verbatim to the density of a curve system. It will be used in a future work to compute limiting objects.

The density defined that way is obviously positive and of unit integral: its entropy is thus well defined and can be computed as: Since the entropy is minimal for concentrated distributions, it is quite intuitive to figure out that seeking for a curve system (γ 1 , . . . , γ N ) giving a minimum value for E(γ 1 , . . . , γ N ) will induce the following properties:

E(γ 1 , . . . , γ N ) = - Ω d(x) log d(x) dx (9)
• The images of the curves tend to get close on to another.

• The individual lengths will be minimized: it is a direct consequence of the fact that the density has a term in γ within the integral that will favor short trajectories.

Using a standard gradient descent algorithm on the entropy produces an optimally concentrated curve system, suitable for use in an unsupervised clustering algorithm. The displacement field for trajectory j is oriented at each point along the normal vector to the trajectory, with norm given by:

Ω γ j (t) -x γ j (t) -x N K ( γ j (t) -x ) log d(x)dx γ j (t) (10)
-

Ω K ( γ j (t) -x ) log d(x))dx γ j (t) γ j (t) N (11)
+ Ω d(x) log( d(x))dx γ j (t) γ j (t) N , (12)
where the notation v |N stands for the projection of the vector v onto the normal vector to the trajectory. An overall normalizing constant:

1 N i=1 l i ,
has to put in front of the expression to get the true gradient of the entropy.

Applying the gradient based algorithm to the traffic situation studied above, yields after convergence the picture presented in figure 4. The clustering procedure was also used in the context of automated trajectory planners. In the example presented in figure 5, 30 flights were simulated on the depicted route system (the figure is a snapshot of the situation, so that not all aircraft are present). A multi-agent system was applied in order to solve all conflicts, but with no constraints put on trajectories: the resulting planning scheme makes each aircraft move in a free flight fashion, and may not be controlled by a human ATC controller. Applying the minimum entropy clustering scheme allows the emergence of a route system. Since conflicts are introduced by this procedure, an iterative application of the multi-agent solver and clustering must be used. After convergence, the conflictfree planning of figure 6 was obtained.

V. EXTENSION TO TRAJECTORIES IN A LIE GROUP

While satisfactory in terms of traffic flows, the previous approach suffers from a severe flaw when one considers flight paths that are very similar in shape but are oriented in opposite directions. Since the density is insensitive to direction reversal, flight paths will tend to aggregate while the correct behavior will be to ensure a sufficient separation in order to prevent hazardous encounters. Taking aircraft headings into account in the clustering process is then mandatory when such situations have to be considered.

Solving this issue is not as straightforward as one may think at a first glance. It is always possible to add a penalty term related to the inner product of speed vectors for each pair of curves, but tuning its relative importance is quite difficult and may impair the overall convergence of the algorithm. A classical workaround that several authors have used in usual multivariate PCA is to enrich the state space of the curves. One can consider for example a state vector of the form:

γ(t) γ (t) . ( 13
)
Instead of using the speed vector to enrich the state, successive positions may be used:

γ(t) γ(t + τ ) . ( 14
)
where the time shift τ has to be chosen to capture instantaneous motion of the aircraft. While satisfactory, two points are not addressed by this formulation:

• Just comparing speed vectors (or couples) using standard euclidean distance is not scale invariant: if velocities are multiplied by a fixed positive scalar, then the inner product will scale with the square of this amount. This is not relevant in most applications as speed vectors comparison is generally assumed not to change under translations/rotations/scalings. • Salient features of the speed vector in an operational context are magnitude and heading. There are not treated as independent parameters in the aforementioned representation.

Keeping the operational requirements in mind, it is more natural to describe the aircraft state as a centroid position and attitude. Since only observed trajectories are considered here, attitude will be described by a unit vector pointing in the direction of the motion and a positive real number giving the velocity. From a mathematical point of view, such a situation may be described by a group action on a reference vector, chosen here to be:

0 d e 1 ,
with e 1 the first basis vector, and 0 d the origin in R d . It is equivalent to model the state in 13 as a linear transformation:

0 d ⊗ e 1 → T (t) ⊗ A(t)(0 d ⊗ e 1 ) = γ(t) ⊗ γ (t),
where T (t) is the translation mapping 0 d to γ(t) and A(t) is the composite of a scaling and a rotation mapping e 1 to γ (t). Considering the vector (γ(t), 1) instead of γ(t) allows a matrix representation of the translation T (t):

γ(t) 1 = Id γ(t) 0 1 0 d 1 .
From now, all points will be implicitly considered as having an extra last coordinate with value 1, so that translations are expressed using matrices. The origin 0 d will thus stand for the vector (0, . . . , 0, 1) in R d+1 . Gathering things yields:

γ(t) γ (t) = T (t) 0 0 A(t) 0 d e 1 . ( 15
)
The previous expression makes it possible to represent a trajectory as a mapping from a time interval to the matrix Lie group G = R d × Σ × SO(d), where Σ is the group of multiples of the identity, SO(d) the group of rotations and R d the group of translations. Please note that all the products are direct. The A(t) term in the expression [START_REF] Michor | Riemannian geometries on spaces of plane curves[END_REF] can be written as an element of Σ ⊗ SO(d). Starting with the defining property A(t)e 1 = γ (t), one can write A(t) = γ (t) U (t) with U (t) a rotation mapping e 1 ∈ S d-1 to the unit vector γ (t)/ γ (t) ∈ S d-1 . For arbitrary dimension d, U (t) is not uniquely defined, as it can be written as a rotation in the plane P = span(e 1 , γ (t)) and a rotation in its orthogonal complement P ⊥ . A common choice is to let U (t) be the identity in P ⊥ which corresponds in fact to a move along a geodesic (great circle) in S d-1 . This will be assumed implicitly in the sequel, so that the representation A(t) = Λ(t)U (t) with Λ(t) = γ (t) Id becomes unique.

The Lie algebra g of G is easily seen to be R d × R × Asym(d) with Asym(d) is the space of skew-symmetric d × d matrices. An element from g is a triple (u, λ, A) with an associated matrix form:

M (u, λ, A) =   0 u 0 0 0 0 λId + A   . (16)
The exponential mapping from g to G can be obtained in a straightforward manner using the usual matrix exponential:

exp((u, λ, A)) = exp(M (u, λ, A)).
The matrix representation of g may be used to derive a metric:

(u, λ, A), (v, µ, B) g = Tr M (u, λ, A) t M (v, µ, B) .
Using routine matrix computations and the fact that A, B being skew-symetric have vanishing trace, it can be expressed as:

(u, λ, A), (v, µ, B) g = nλµ + u, v + Tr A t B . (17)
A left invariant metric on the tangent space T g G at g ∈ G is derived from [START_REF] Griffiths | Principles of Algebraic Geometry[END_REF] as:

X, Y, g = g -1 X, g -1 Y g , with X, Y ∈ T g G.
Please note that G is a matrix group acting linearly so that the mapping g -1 is well defined from T g G to g. Using the fact that the metric (17) splits, one can check that geodesics in the group are given by straight segments in g: if g 1 , g 2 are two elements from G, then the geodesic connecting them is:

t ∈ [0, 1] → g 1 exp t log g -1 1 g 2 .
where log is a determination of the matrix logarithm. Finally, the geodesic length is used to compute the distance d(g 1 , g 2 ) between two elements g 1 , g 2 in G. Assuming that the translation parts of g 1 , g 2 are respectively u 1 , u 2 , the rotations U 1 , U 2 and the scalings exp(λ 1 ), exp(λ 2 ) then:

d(g 1 , g 2 ) 2 = (λ 1 -λ 2 ) 2 + ( 18
)
Tr log U t 1 U 2 log U t 1 U 2 t + u 1 -u 2 2 . ( 19
)
An important point to note is that the scaling part of an element g ∈ G will contribute by its logarithm: this ensures the scale invariance of the metric as any common scaling factor will cancel out in the distance computation. The part related to unit matrices can be interpreted as finding an optimal sequence of rotations that will make the eigenvectors coincident and compute the total angular variation needed. It has to be noted that when dealing with planar curve, the value obtained that way is exactly the angle between the speed vectors. Finally, the metric differs from the affine invariant one in that orientation is explicitly taken into account.

Based on the above derivation, a flight path γ with state vector (γ(t), γ (t)) will be modeled in the sequel as a curve with values in the Lie group G:

Γ : t ∈ [0, 1] → Γ(t) ∈ G, with: Γ(t).(0 d , e 1 ) = (γ(t), γ (t)).
In order to make the Lie group representation amenable to statistical thinking, it is needed to define probability densities on the translation, scaling and rotation components that are invariant under the action of the corresponding factor of G.

VI. DENSITY ESTIMATION ON G

The kernel density estimator presented in the first part of the paper cannot be used directly as it relies explicitly on a vector space structure. It has to be adapted for the Lie group case, so as to ensure invariance by the group action.

Since the Lie group G is a product, the three components can be treated separately.

The translation part is the additive group R d and its action is given by a vector translation. As a consequence, its group structure is identical to the one of R d considered as a real vector space: the density estimation on it will not differ from the one presented in IV. The same estimator will be used, with either the epanechnihov or the Gaussian kernel.

Concerning the SO(d) component, it well known that an arbitrary rotation U can be represented as a sequence of points on spheres of decreasing dimension: the procedure was already detailed in the section V. Estimating a density on SO(d) can thus be done using kernels on spheres: it becomes a classical problem in directional statistics.

A commonly used choice is the von Mises-Fisher (vMF) distribution on S d-1 which is denoted M(m, κ) and given by the following density expression [START_REF] Mardia | Directional Statistics, ser. Wiley Series in Probability and Statistics[END_REF]:

K V M F (x; m, κ) = c d (κ) e κm T x , κ > 0, x ∈ S d-1 , ( 20
)
where

c d (κ) = κ d/2-1 (2π) d/2 I d/2-1 (κ) (21)
is a normalization constant with I r (κ) denoting the modified Bessel function of the first kind at order r. The vMF kernel function is an unimodal p.d.f. parametrized by the unit mean-direction vector µ and the concentration parameter κ . κ is a smoothing parameter that plays the role of the inverse of the bandwidth parameter as defined in the linear kernel density estimation. Large values of κ imply greater a concentration around the mean direction and lead to undersmoothed estimators whereas small values provide oversmoothed circular densities [START_REF] García-Portugués | Kernel density estimation for directional-linear data[END_REF]. Indeed, if κ = 0, the vMF kernel function reduces to the uniform circular distribution on the hypersphere. Note that the vMF kernel function is by design invariant by the action of SO d . The VNf density enjoys properties analogous to those of multivariate Gaussian distribution: a central limit theorem exists in the context of spheres, with limiting distribution the VMf.

The vMF distribution may be expressed by means of the spherical polar coordinates of x ∈ S d-1 [START_REF] Mardia | Statistics of directional data[END_REF].

Given the random vectors X i , i = 1, . . . , n, in S d-1 , the estimator of the spherical distribution is given by:

f (x) = 1 n n i=1 K V M F (x; X i ) = c d (κ) n n i=1 e κX T i x , κ > 0, x ∈ S d-1 .
The density estimator on S d is constructed componentwise, using the previous expression on each of the spheres involved in its decomposition.

As for the scaling component of G, the usual kernel functions such as the Gaussian and the Epanechnikov kernel functions are not suitable due to a bias introduced in close to the domain boundary at 0. An asymmetrical kernel function on R + such as the log-normal kernel function is a more convenient choice. Moreover, this p.d.f. is invariant by the action of the multiplicative group of positive real numbers. Let R 1 , . . . , R n be univariate random variables from a p.d.f. which has bounded support on [0; +∞[. The radial density estimator may be defined by means of a sum of log-normal kernel functions as follows:

g(r) = 1 n n i=1 K LN (r; ln R i , h), r ≥ 0, h > 0,
where

K LN (x; µ, σ) = 1 √ 2πσx e -(ln x-µ) 2
is the log-normal kernel function and h is the bandwidth parameter. The resulting estimate is the sum of bumps defined by log-normal kernels with medians R i and variances (e h 2 -1)e h 2 R 2 i . Note that the log-normal (asymmetric) kernel density estimation is similar to the kernel density estimation based on a log-transformation of the data with the Gaussian kernel function: in practice, the estimator IV is used on the logarithms of the scale components.

VII. UNSUPERVISED ENTROPY CLUSTERING

The derivation of the algorithm used for clustering curves with value in the Lie group is similar to the vector space case [START_REF] Rani | Recent techniques of clustering of time series data: A survey[END_REF]. The kernel has to be replaced with the respective VMf and Log-normal ones on the rotation and scaling components. The term γ (t) in the original expression of the density, that is required to ensure invariance under reparametrization of the curve, has to be changed according to the metric in G and is replaced by γ (t), γ (t)

1/2 γ(t)
. The density at x ∈ G is thus:

d G (x)) = N i=1 1 0 K (x, γ i (t)) γ i (t), γ i (t) 1/2 γi(t) dt N i=1 l i (22)
where l i is the length of the curve in G, that is:

l i = 1 0 γ i (t), γ i (t) 1/2 γi(t) dt (23)
The expression of the kernel evaluation K (x, γ i (t)) is split into three terms. In order to ease the writing, a point x in G will be split into x r , x s , x o components where the exponent r, s, t stands respectively for translation, scaling and rotation. Given the fact that K is a product of component-wise independent kernels it comes:

K (x, γ i (t)) = K t x t , γ t i (t) K s (x s , γ s i (t)) K o (x o , γ o i (t))
where:

K t (x t , γ t i (t)) = ep x t -γ t i (t) (24)
K s (x s , γ s i (t)) = 1 x s σ √ 2π exp - (log x s -log γ s i (t)) 2 2σ 2 (25) K o (x o , γ o i (t)) = C(κ) exp κTr x ot γ o i (t) (26)
with C(κ) the normalizing constant making the kernel of unit integral. Please note that the expression given here is valid for arbitrary rotations, but for aircraft trajectory the heading is the only observable parameter unless onboard data is downlinked. As a consequence it is enough to consider it only a standard von-mises distribution on S d-1 :

K o (x o , γ o i (t)) = C(κ) exp κx ot γ o i (t)
with normalizing constant as given in [START_REF] Mardia | Directional Statistics, ser. Wiley Series in Probability and Statistics[END_REF]. The entropy of the system of curves is obtained from the density in G:

E(d G ) = - G d G (x) log d G (x)dµ G (x) (27)
with dµ G the left Haar measure. Using again the fact that G is a direct product group, dµ is easily seen to be a product measure, with dx t , the usual Lebesgue measure on the translation part, dx s /x s on the scaling part and the lebesgue measure dx o on S d-1 for the rotation part. It turns out that the 1/x s term in the expression of dx s /x s is already taken into account in the kernel definition, due to the fact that it is expressed in logarithmic coordinates. The same is true for the Von-Mises kernel, so that in the sequel only the (product) lebesgue measure will appear in the integrals.

Finding the system of curves with minimum entropy requires a displacement field computation as detailed in [START_REF] Puechmorel | Entropy minimizing curves with application to automated flight path design[END_REF]. For each curve γ i , such a field is a mapping η i : [0, 1] → T G where at each t ∈ [0, 1], η i (t) ∈ T G γi(t) .Compare to the original situation where only spatial density was considered, the computation must now be conducted in the tangent space to G. Even for small problems, the effort needed becomes prohibitive. The structure of the kernel involved in the density can help in cutting the overall computations needed. Since it is a product, and the translation part is compactly supported, being an epanechnikov kernel, one can restrict the evaluation to points belonging to its support. Density computation will thus be made only in tubes around the trajectories.

Second, for the target application that is to cluster the flight paths into a route network and is of pure spatial nature, there is no point in updating the rotation and scaling part when performing the moves: only the translation part must change, the other two being computed from the trajectory. The initial optimization problem in G may thus be greatly simplified.

Let be an admissible variation of curve γ i , that is a smooth mapping from [0, 1] to T G with (t) ∈ T γi(t) G and (0) = (1) = 0. We assume furthermore that has only a translation component. The derivative of the entropy E(d G ) the t curve γ i is obtained from the first order term when γ i is replaced by γ i + . First of all, it has to be noted that d G is a density and thus has unit integral regardless of the curve system. When computing the derivative of E(d G ), the term

- G d G (x) ∂ γi d G (x) d G (x) dµ G (x) = - G ∂ γi d G (x)dµ G (x)
will thus vanish. It remains:

- G ∂ γi d G (x) log d G (x)dµ G (x)
The density d G is a sum on the curves, and only the i-th term has to be considered. Starting with the expression from ( 22), one term in the derivative will come from the denominator. It computes the same way as in [START_REF] Puechmorel | Entropy minimizing curves with application to automated flight path design[END_REF] to yield:

γ t i (t) γ i (t), γ i (t) G N E(d G ) (28)
Please note that the second derivative of γ i is considered only on its translation component, but the first derivative makes use of the complete expression. As before, the notation | N stands for the projection onto the normal component to the curve.

The second term comes from the variation of the numerator. Using the fact that the kernel is a product K t K s K o and that all individual terms have a unit integral on their respective components, the expression becomes very similar to the case of spatial density only and is:

- G K (x, γ i (t)) log d G (x)dµ G(x) γ t i (t) γ i (t), γ i (t) 1/2 G N (29) + R d e(t)K t x t , γ t i (t) log d G (x) γ i (t), γ i (t) 1/2 G dx t (30)
with:

e(t) = γ t i (t) -x t γ t i (t) -x t N VIII. RESULTS
Although the original Lie group formulation can be simplified in the framework of aircraft trajectories, the computation cost is still too high to be used in practical applications. A preprocessing phase was thus applied to the dataset in order to limit the number of terms used in the expression of the entropy or its gradient. The information about a trajectory can be summarized by the horizontal component of its centroid and its mean heading, yielding a three dimensional point. Using it within an octree allows the density computation to be limited not only to spatial areas close enough to existing trajectories but also to curves with similar orientation. It becomes possible to save a lot of useless kernel evaluations in estimating the density on the rotation part.

The arrivals and departures at Toulouse Blagnac airport were analyzed. The algorithm performs well as indicated on Figure 7. Four clusters are identified, with mean lines represented through a spline smoothing between landmarks. It is quite remarkable that all density based algorithms were unable to separate the two clusters located at the right side of the picture, while the present one clearly show a standard approach procedure and a short departure one. For the case of Figure 7, that represents 1784 trajectories, the computation time is 5 minutes on a XEON 3Ghz machine and with a pure java implementation.

IX. TOWARD ANALYSIS OF HIGH DIMENSION DATA

When dealing with data coming from quick access recorder (QARs), the size of the samples arises as a new issue. It is no longer possible to use the techniques presented in the previous sections as the estimation of the density will become prohibitive in terms of computational complexity. Furthermore, the various components appearing in the samples may be of very different nature: it is thus required to perform first a segmentation into homogeneous groups. A promising approach is the one described in [START_REF] Marrinan | Numerical Mathematics and Advanced Applications -ENUMATH[END_REF] where the data is represented as point in a flag of invariant subspaces. It is in some sense a generalization of the classical PCA decomposition, but with principal subspaces used instead of vectors: the obvious benefit is that less computation is needed and that all vectors pertaining to the same subspace may be considered equivalent. This gives a clue for a adaptation of the entropy based algorithms to very high dimensional data:

• Compute the flag decompositions;

• For a sample point, get its distance to each of the subspaces involved in the flag. Please note that it can be viewed as a kind of generalized coordinates, the usual expansion on a orthonormal basis being recovered by letting the flag by composed only of 1-dimensional subspaces; • Use the previous expansion as the input data for an entropy based algorithm.

This work is still in early research stage, however it seems to be a very promising way of addressing the issue of QARs data clustering. Two important points are yet to be solved:

• Find an optimal way to construct the flag of invariant subspaces;

• Adapt the entropy based algorithm so as to be able to deal with samples with dimension in the order of 10, typical of what is provided by the flag expansion.

X. CONCLUSION AND FUTURE WORK

The entropy associated with a system of curves has proved itself efficient in unsupervised clustering application where the information coming from geometry is of primary importance. Extending this setting to the case of Lie group valued curves adds the ability to take into account features of the data that are not of spatial nature, like heading or scale. On several test problems, the clusters produced extract relevant aeronautical information, yielding one of the first mathematically sound algorithm able to deal with such data. The next major step will be the introduction of invariant subspaces decompositions, that will naturally produce a summary of very high dimensional samples like produced by the QARs. It is anticipated to pave the way for a new generation of algorithms able to operate in this context.

Fig. 1 .

 1 Fig. 1. Traffic over France the 12th February 2013

Fig. 3 .

 3 Fig. 3. Geometric density for the 12th February 2013

Fig. 4 .

 4 Fig. 4. Bundling one day of traffic

Fig. 5 .

 5 Fig. 5. Initial flight plan.

Fig. 6 .

 6 Fig. 6. Entropy minimal curve system from the initial flight plan.

Fig. 7 .

 7 Fig. 7. Bundling trajectories at Toulouse airport


Nicolas Peteilh

Towards a robust multidisciplinary design optimization model for the airplane in the air transport system

Keywords:

In an airplane preliminary design phase, the optimization for costs and fuel consumption is based on anticipating the future operations of the airplane by defining one representative mission. How effective is this process? In this study, based on data exploration and statistical analysis, we identify and represent the operational parameters impacting the variability of the actual missions of an airplane. These parameters include the airline route network and flight optimization, Air Traffic Management (ATM) considerations, and the atmospheric and meteorological phenomena. Our results show how significant some impacts are, mainly on actual flight distances but also on flying conditions. A conclusion of our study is that airplanes are operated on missions different from the reference mission considered during the preliminary design phase and this has a significant impact on operating costs and fuel burn efficiency.

INTRODUCTION

In its Strategic Research and Innovation Agenda [1] (SRIA), the Advisory Council for Aeronautics Research in Europe (ACARE) has identified the protection of the environment and the energy supply as a challenge for the future and the definition of the air vehicles of the future as a key action area. It also emphasizes the need for developing an integrated resilient air transport system, where the airplane, together with the other actors of the civil aviation world such as the airports, the airlines, the air navigation service providers, the authorities, etc. interact one with each other and influence each other's design. In this perspective, considering the airplane within its real operational environment appears to be a way of improvement, even in its early preliminary design phase, when decisions have a great impact on future operational efficiency and on costs.

Focusing on the early design phase of an airplane, the benefits of multidisciplinary design process and the limits of airplane design optimization based on a reduced number of reference flight conditions have already been described in the literature. See for instance [START_REF] Liem | Multimission Aircarft Fuel-Burn Minimization via Multipoint Aerostructural Optimization[END_REF], where ways forward are proposed through multi-mission and multipoint approaches. Reference [START_REF] Liem | Multimission Aircarft Fuel-Burn Minimization via Multipoint Aerostructural Optimization[END_REF] describes a method in order to obtain a set of representative missions from inservice data and demonstrates the interest of considering real operational conditions in improving the results of multidisciplinary design optimization (MDO). However, it relies mainly on the BTS flight database [3], which features aggregation of operational data, the atmosphere is considered as standard, and the mission profile is determined by performing theoretical-mission analysis.

In this paper, we shall demonstrate the benefits of taking into account more details in describing the conditions in which the airplane is operated in the real world, in order to capture the potential effects of airlines routes network, operational flight optimization, Air Traffic Management (ATM) considerations, and the atmospheric and meteorological conditions.

One key challenge here is to find reliable and detailed operational data that enable one to have a deeper insight of the real conditions in which airplanes operate. Another objective is to extract from these data the relevant information. This paper is organized as follows. Section 2 describes the main identified operational uncertainty parameters. Section 3 addresses the data processing from the selection of sources to the methods used. Section 4 presents the analysis of the statistical results. Section 5 concludes the paper.

MAIN UNCERTAINTY PARAMETERS

In order to capture the potential effects of airline route network, operational flight optimization, Air Traffic Management (ATM) considerations and atmospheric and meteorological conditions, we first investigated how each of these conditions affects the airplane operations. This investigation is made by engineering judgement, at the best knowledge of the actual operational world. The parameters identified are those the most exposed to uncertainties.

The payload

The very first parameter that affects the airplane is the payload. It changes at every flight the mass of passengers and fret.

Airlines aim at optimizing their revenue and improving their load factor. Nevertheless, they must sometimes operate empty or half empty airplanes. The first data of interest is therefore the actual weight of the airplane at take-off.

Route network

The route network has a direct impact on the ranges the airplanes fly. Indeed, even though the use of a new type of airplane may allow an airline to modify its route network, the localization of the airports to be operated and the flight schedule to meet passenger demand will determine the ranges flown by airplanes. Therefore, the second set of data needed includes the departures, destinations and frequencies of flights performed the type of airplane under study.

Figure 1. Short-range flight profile

Operational flight optimization

When preparing their flight, airlines choose the ground track to be followed as well as the flight path profile. This analysis depends on the payload and ensures that all operational limitations of the airplane are met, the safety requirements are addressed, and the operational costs are minimized. As a result, the chosen path from one airport to another might well not be the direct route.

The third data set needed is the actual ground track chosen by the airlines, which corresponds to the latitude and longitude coordinates all along the airplane trajectory.

Remark also that airlines might choose longer ground tracks in order to obtain shorter air distances. The fourth parameter that matters is therefore either the wind data or the airplane true airspeed.

ATM considerations

Once the airline has defined its optimized mission, it submits the flight plan for approval. After being processed together with the other airline flight plans, for instance in Europe by the NMOC (Network Manager Operations Centre), the final approved flight plan may be different from the one requested. The modifications may impact departure time, flight levels, or even the ground track. Furthermore, the ground trajectory may also be impacted by ATC (Air Traffic Control).

The fifth parameter of importance is therefore the altitude of the flight. As shown in Fig. 1 and Fig. 2, range and ATC can have an impact on altitude. As airplanes consider barometric altitude, it is even more important to get the atmospheric pressure all along the flight track.

The atmospheric and meteorological conditions

The atmospheric conditions are the last effect to be investigated. As atmospheric pressure is already considered, the sixth parameter that is needed is the air temperature all along the flight path.

There are also local specific phenomena, like storms, that require real-time flight path modification by the pilot, and will impact the flight conditions. However, studying such phenomena goes beyond the scope of this paper.

Synthesis

The main identified uncertainty parameters considered in this study are the following: The data investigation was done in order to observe and quantify the variability encountered in service.

DATA PROCESSING

Researching data, we identify two sources of interest. We then collect data and after exploring them, we can represent the operational variability with respect to the main identified uncertainty parameters. We focus our analysis on Airbus long range airplanes and more precisely on A340 types.

Data sources

It was not possible to find any set of data including all the relevant parameters we wanted to explore. As a result, our approach is to get sufficient data to explore independently each parameter. Here are the sources on which we base our analysis.

Reference [3] includes data related to payload (passengers and freight) carried per month, airline, airplane type, and route (departure and destination). Direct route distances are also included. This database contains only U.S. domestic and international air traffic data. The main drawback of this first database is that it does not contains the trajectory detail (latitude, etc.) described at the end of section 2.

Reference [START_REF]MOZAIC IAGOS Database[END_REF] provides the all-along-the-flight-path data that we need. Nevertheless, payload and mass information is missing.

Data processing

The approach used for data processing followed five steps. Data are extracted from [START_REF]MOZAIC IAGOS Database[END_REF], for one specific airplane during one year. More than 400 trajectories were analysed. Fig. 1 and Fig. 2 show two examples of vertical flight path.

Figure 7. Total air distance per flight altitude

We then focus on the variability related to air temperature. Fig. 8 shows the difference between the air temperature actually measured on the flight path, and the International Standard Atmosphere temperature at the altitude considered.

ANALYSIS OF THE STATISTICAL RESULTS

The results presented illustrate the variability related to the main identified uncertainty parameters.

Variability and operational parameters

The direct route distance distribution (Fig. 3) shows that even long-range type airplanes can be operated on very short-range trips. This is not specifically related to the airplane analysed in this study but is rather a more general observation, as shown in Fig. 9. This situation cannot be improved by the airlines as there is no conventional distance between the airports to be flown into.

The distance that the airplane actually flies is the air distance. From the analysis of Fig. 4, Fig. 5 and Fig. 6, we identify the operational conditions and parameters that affect the air distance flown.

When optimizing the ground trajectory, the airlines take into account costs, landscape, payload, atmospheric conditions, ATM, airports constraints, airplanes limits and safety considerations. This study shows that these parameters can result in significant increase of the actual ground distances flown by the airplane under consideration in this paper. Fig. 4 shows that, even for long-haul flights, the actual ground distance can be more that 20% longer than the direct route distance. Fig. 5 and Fig. 6 illustrate the differences between ground and air distances. The main parameter that explains these differences is the wind. The distribution presented appears mainly symmetrical, which means that the wind has both positive and negative impacts and almost in the same proportion. It appears consistent with operational conditions. However, one might have expected to observe more drawbacks than benefits from the wind. Indeed, if head wind could explain this symmetry, side wind would induce a shift towards longer air distances. The balance between positive and negative impacts observed in Fig. 5 can however be explained by the strategy of the airlines to modify the flight track to reduce the negative impact and maximise the benefits of the wind.

The variability related to altitude (Fig. 7) is mainly related to ATM. The cruise altitudes are aligned with the available flight levels defined to ensure vertical separation between the airplanes. However, there are some cruise altitudes that are different from the ATM predefined flight levels. We have not identified the operational parameter related to this variability yet.

The temperature shows variability too (Fig. 8). The airplane can operate on colder or warmer environment. The specific airplane considered in this study mainly operated in the northern hemisphere. The shape of the distribution shows two peaks. Our first hypothesis to explain this is the seasonal variation of the temperature. The left peak would correspond to winter and the right one to summer. Finally, Fig. 9 reveals the variability related to the payload is very large. Our analysis is that the uncertainty on payload depends on the routes, the day of the week, the period of the year and the efficiency of the airline to sell flight tickets. It could also be explained by operational limitations at the airport for instance. The correlation has not been investigated yet.

Variability for the airplane design

During the preliminary phases of an airplane design, this one is usually optimized with respect to one referenced mission. Aerodynamics, structures, systems and engines are then optimized to minimize the weight, the fuel consumption and the operating costs or to maximize the range and the payload. Based on these analyses, decisions are made to define the specific features of the future airplane. Reference [START_REF] Liem | Multimission Aircarft Fuel-Burn Minimization via Multipoint Aerostructural Optimization[END_REF] underlines the risks related to optimising an airplane based on one mission only.

The representations we obtained of the operational variability during the work presented here can be the basis for developing models that will be adequate for airplane preliminary design and for proposing robust methods to cope with the described uncertainty parameters.

CONCLUSION

The main identified uncertainty parameters related to airplane operations and real civil aviation world activity have been listed and described. Using databases and statistical data processing, the related variability was plotted and analysed to make a link between its statistical distribution and various operational parameters.

The benefit demonstrated by this study is that it seems possible to correlate the observed variability to operational parameters, such as airlines routes network, operational flight optimization, ATM, and the atmospheric and meteorological conditions. Based on this work, the next step will be to elaborate models that will capture and represent this variability so that it can be taken into account in the preliminary design process of an airplane, where the decisions taken are of great importance for the future efficiency of the designed airplane.

Figure 2 .

 2 Figure 2. Long-range flight profile

1. Airplane type 2 .

 2 Route information (departure, destination and frequency) 3. All along the trajectory: latitude and longitude coordinates, wind data (or true airspeed), barometric altitude, and air temperature 4. Payload

Fig. 3

 3 Fig.3 displays the direct route distance distribution over one year for one specific airplane. It shows the variability of the airlines' route networks. It also sets a reference for the remaining of the study.

Figure 3 .

 3 Figure 3. Direct route distance distribution Fig.4 compares the actual ground path with the direct route distance.

Figure 4 .

 4 Figure 4. Relative increase of ground distanceIn order to take into account the effect of the wind, we plot the difference between air distance and ground distance on Fig.5.

Figure 5 .

 5 Figure 5. Air-to-ground distance distribution Fig.6 illustrates difference between air and ground distances.

Figure 6 .

 6 Figure 6. Relative increase of air vs. ground distancesFig.7illustrates the variability related to cruise altitude by plotting the total air distance flown by the airplane during one year (Fig.7) for each cruising altitude. As it is calculated based on atmospheric pressure along the flight path, this plot provides also information on the pressure variability.

Figure 8 .

 8 Figure 8. Delta ISA Temperature (in °C)

Figure 9 .

 9 Figure 9. Payload-Range distribution


Thierry Miquel
email: thierry.miquel@enac.fr

DEVELOPMENT AND ILLUSTRATIVE EXAMPLES OF AN AIRCRAFT LONGITUDINAL REFERENCE PROFILE TO COMBINE CONTINUOUS DESCENT OPERATIONS AND TAILORED ARRIVALS

Keywords: Tailored arrival, Continuous Descent Operation, reference speed profile, reference vertical profile

This paper addresses a specific aspect of air traffic control services, namely the achievement of an orderly and expeditious flow of air traffic under time constrained continuous descent approach. The task of establishing properly spaced landing sequences is quite demanding for air traffic controllers, especially under heavy traffic conditions. Indeed, terminal control area (TMA) air traffic controllers have to merge two or more streams into a single stream by means of radar vectoring and speed instructions. This high level task of sequencing aircraft is not currently communicated to the pilot. In this paper, the task of merging an aircraft over a specified meter fix is addressed through a novel clearance in which air traffic control clears an aircraft to track an ad-hoc computed reference trajectory. This enables the aircraft to merge a specified meter fix at a given time. This kind of application may be envisioned as an enhancement of the use of the ground based Arrival MANager (AMAN), which is a tactical controller assistance system enabling the computation of rendez-vous time at meter fix to meet the runway capacity and absorb the traffic. The time computed by the AMAN is envisioned to be used by the air traffic controller to request the aircrew to overfly the meter fix at the desired time. In comparison with current operations, the change is that the controller would communicate decisions on traffic flow organization at a higher level to the pilot rather than requiring the controller to calculate and communicate lower-level guidance instructions.

While keeping the controller responsible for making the traffic flow decisions, the envisioned application involves new avionics capabilities for merging operations, including the tracking of height and speed. Thus, air traffic control improvement is achieved through a greater involvement of pilots in cooperation with air traffic controllers. This type of application is clearly in the scope of the 4D trajectory concept promoted by the European programme SESAR and by the US programme NextGen.

In this paper a futuristic 4D trajectory application where the air traffic controller will ask an aircraft to overfly a meter fix at a specific time is addressed. The overflying time is assumed to be specified by the air traffic controller to settle properly the arrival sequence and could be given for example through an Arrival Manager (AMAN). In addition, this clearance is assumed to be given after the Top Of Descent (TOD) of the aircraft. This paper focuses on the issue of computing reference airspeed and height profiles for time constrained descent. The proposed approach is illustrated through some examples.

INTRODUCTION

Nowadays, environmental impact and efficiency have become the two very important aspects in aviation industry after safety. New operations such as Continuous Descent Operation (CDO) can significantly reduce the noise impact of landing aircraft by keeping them longer at higher altitudes and by avoiding steps during descent. Such new operations need a collaborative work between airlines and Air Navigation Service Providers (ANSPs) to be defined and operated. CDO is defined by ICAO as "an aircraft operating technique aided by appropriate airspace and procedure design and appropriate ATC clearances enabling the execution of a flight profile optimized to the operating capability of the aircraft, with low engine thrust settings and, where possible, a low drag configuration, thereby reducing fuel burn and emissions during descent. The optimum vertical profile takes the form of a continuously descending path, with a minimum of level flight segments only as needed to decelerate and configure the aircraft or to establish on a landing guidance system (e.g. ILS)" [START_REF] Miquel | Continuous Descent Operations (CDO) Manual[END_REF]. [7] aims to facilitate the merging of traffic from a number of Area Navigation (RNAV) arrival routes. The technique is based upon aircraft flying a quasi-arc, up to 30NM long, with a radius of more than 20NM from the designated merging point. Each arc has a published height that the aircraft must have reached before establishing on the arc and a predefined speed to fly it. In general the arc nearest to the merging point has the highest height while the other has the lowest height so that the external sequencing leg is free from traffic from the internal sequencing leg during descent of the traffic. On April 2011 Oslo became the world's first airport to implement a PMS in their airspace. Other countries and airports will follow.

In this paper, the task of merging an aircraft over a specified meter fix is addressed through a novel clearance in which air traffic control clears an aircraft to merge a specified meter fix at a given time. This kind of application may be envisioned as an enhancement of the use of the ground based Arrival MANager (AMAN), which is a tactical controller assistance system enabling the computation of rendez-vous time at meter fix to meet the runway capacity and absorb the traffic [8]. The time computed by the AMAN is envisioned to be used by the air traffic controller to request the aircrew to overfly the meter fix at the desired time. In comparison with current operations, the change is that the controller would communicate decisions on traffic flow organization at a higher level to the pilot rather than requiring the controller to calculate and communicate lower-level guidance instructions.

While keeping the controller responsible for making the traffic flow decisions, the envisioned application involves new avionics capabilities for merging operations including the tracking of height and speed. Thus air traffic control improvement is achieved through a greater involvement of pilots in cooperation with air traffic controllers. This type of application is clearly in the scope of the 4D trajectory concept promoted by the European programme SESAR [9] and by the US programme NextGen [10].

The main benefit expected from this application [11] is to improve flight efficiency by more precise maneuvering resulting from onboard capabilities, and also noise abatement and fuel saving. More precise maneuvers are also expected to increase sector capacity. Indeed moving from radar vectoring to monitoring pre-computed trajectories would contribute to decrease controller's workload, and therefore to increase sector capacity. This paper presents a new methodology to compute a reference trajectory for time based continuous descent operations. The aircraft is considered as a point mass model and focuses on aircraft longitudinal motion including known wind. As far as time constrained operations are assumed, final time as well as final altitude and along track distance to be flown are imposed. We propose a new methodology to compute a reference indicated airspeed (IAS) and a reference vertical profile to achieve imposed final position and altitude at a prescribed time. The computed trajectory is a time parametrized trajectory which will be used as a reference trajectory by some envisioned tracking controller installed on board the aircraft. Nevertheless the design of the tracking controller is out of the scope of this paper.

The structure of the paper is the following: first the equations of motion which will be used are presented as well as the constraints to be satisfied. Then the Two-point Boundary Value Problem (TPBVP) to be solved is defined. The core of the paper consists in describing the proposed methodology to compute the reference indicated airspeed (IAS) and reference vertical profile to solve the problem. Then a scenario is proposed and simulation results are provided to illustrate the effectiveness of the proposed approach.

EQUATIONS OF MOTION

In the following, the aircraft is modeled as a point mass flying on a non-rotating and flat earth. The aircraft is assumed to fly in an time varying atmospheric wind field comprising of both horizontal and vertical components. In addition the side-slip angle is assumed to be zero so that the flight is coordinated as depicted on Figure 1.

Figure. 1. Reference frame and notations

The notations which will be used are the following:



V is the true airspeed (TAS) of the aircraft  x and y denote the position of the aircraft in the horizontal plane.

 Wh is the horizontal component of the wind  ψw is the direction from where the wind is blowing  ψ is the heading of the aircraft Kinematic equations in the horizontal plane reads as follows:

{ ẋ =V p cos(ψ)+W h cos( ψ w +π) ẏ=V p sin(ψ)+W h sin(ψ w +π ) (1)
Denoting by γ the flight path angle the relationship between Vp and V is the following:

V p =V cos  (2)
 s is the horizontal path length flown by the aircraft. Infinitesimal path length increment ds is linked to infinitesimal positions increments dx and dy by the following relationship

ds 2 =dx 2 +dy 2 (3)
 Vg is the ground speed of the aircraft

ṡ=V g =√ ẋ2 + ẏ2 (4)
 χ is the track angle of the aircraft

{ ẋ=V p cos W h cos w  =V g cos X  ẏ=V p sin  W h sin  w  =V g sin X  (5)
Rearranging (5) leads to the following relationships:

{ ṡ=V g =  V p 2 W h 2 t -2 V p W h cos  - w  = X -arcsin  W h V p sin X - w   (6)
 z is the vertical position (altitude) of the aircraft  Wz is the vertical component of the wind

The kinematic equation in the vertical plane reads as follows:

ż=V sin W z t  (7)
Pilots usually use the indicated airspeed on board aircraft. Conversion between true airspeed V and indicated airspeed Vi can be done using a nonlinear relationship provided for example within BADA documentation [12]:

V =g (V i , z) (8) 3 y W h V g V cos(γ)
x π+ψ w 0

When the indicated airspeed Vi is approximated by the equivalent airspeed Ve we get:

V =V e √ ρ(0) ρ(z) ≈V i √ ρ(0) ρ(z ) (9)
Where ρ(z) is the air density which depends on altitude z and ρ(0) is the air density at sea level.

PROBLEM FORMULATION

In aviation, a standard terminal arrival route or standard terminal arrival (STAR) is a published procedure followed by aircraft on an IFR flight plan just before reaching a destination airport [13].

Figure. 2. Example of standard terminal arrival route

On usual approach chart the flight path to be followed by the aircraft is known by the aircrew. As a consequence the track χ of the aircraft can be easily computed as a function of the horizontal path length s meaning that χ(s) is known during the all the path to be followed by the aircraft. Thus the path length s will be envisioned as the natural independent variable of the aircraft path rather than time t.

Taking the path length s as the new independent variable, we will denote by t' the derivative of time t with respect path length s and by z' the derivative of altitude z with respect path length s. Equations ( 6) and ( 7) are then rewritten as follows:

{ t ' = d t ds = 1 V g z ' = dz dt d t ds (10)
That is:

{ t ' = 1  V p 2 W h 2 s-2 V p W h cos   - w  z ' = V sin W z  V p 2 W h 2 s-2V p W h cos - w  (11)
Heading ψ(s) can be computed as a function of track χ(s), direction ψw(s) from where the wind is blowing, horizontal component Wh(s) of the wind and speed Vp:

= X -arcsin  W h V p sin X - w   (12)
Where speeds Vp is linked to the indicated airspeed Vi and flight path angle γ by ( 2) and (8).

V p =V cos(γ)=g (V i , z )cos(γ) (13)
Airliners are equipped with autopilot systems. Some autopilot modes control the indicated airspeed Vi as well as the flight path angle γ. So the two variables which are chosen to form the control vector u of the aircraft are the indicated airspeed Vi and the flight path angle γ:

u=[ V i  ] T (14)
Equations ( 11) can be rewritten as follows:

x ' =f ( x , u , p )

Where:

 u is the control vector defined in (14) 

x is the state vector

x=[t z ] T ( 16
)
 p is the vector of exogenous parameters depending on path length s

p= [ W h (s) W z (s) ψ w (s) Χ (s) ] T (17)
In order to satisfy the time constrained maneuver the following constraints are imposed on initial and final values of time t with respect to path length s:

{ t (0)=0 t (s f )=t f (18)
The preceding constraints indicate that the aircraft starts the maneuver at initial time 0 and finishes it at the prescribed time of arrival tf after having flown a distance equal to sf.

Similarly the constraints on initial and final values of the vertical components of state vector x are the following:

{ z 0=z 0 z s f = z f (19)
The preceding constraints indicate that the aircraft starts the maneuver at altitude z0 and finishes it at altitude zf at the prescribed arrival path length sf.

In terms of vertical speed we impose the following initial and final values:

{ ż 0=V z0 ż t f =V zf (20)
The preceding constraints indicate that the vertical speed at the beginning and at the end of the maneuver are imposed to be Vz0 and Vzf respectively. Moreover the indicated airspeed Vi is set to be Vi0 at the beginning of maneuver and Vif at the prescribed time of arrival tf:

{ V i 0=V i0 V i t f =V if (21)
Then the problem to be solved consists in finding and expression of controls Vi(s) and γ(s) such that differential equations (11) hold under constraints (18), ( 19), ( 20) and (21). Constraints (18), ( 19), ( 20) and (21) associated with the set of coupled ordinary differential equations (11) constitute a Two-Point Boundary Value Problem (TPBVP).

The key step of our approach to solve this Two-Point Boundary Value Problem is the parametrization of the indicated airspeed Vi and flight path angle γ, that is the control vector u, along the path.

CONTROL PARAMETRIZATION

We will denote by τ the ratio between the actual path length s and the total distance to be flown during the maneuver sf:

τ= s s f (22)
As a consequence when path length s varies from 0 to sf parameter τ varies from 0 to 1. When identifying r ' (s) with r ' (s f τ) we get:

d d s r (s)=r ' (s)≡r ' (s f τ)= 1 s f d d τ r (s f τ ) (23)
To control the set of coupled ordinary differential equations (11) we shall provide the expression of the commanded indicated airspeed Vi and commanded flight path angle γ. Quadratic Bézier curves will be used to build both control. The expression of the commanded indicated airspeed Vi is chosen as follows:

V i =1- 2 V i0 2 1-u V  2 V if (24)
It is clear that whatever constant parameter uV is the proposed expression for the indicated airspeed Vi satisfied the constraints (21) on initial and final values:

{ V i (0)=V i 0 V i (1)=V if (25)
Similarly the vertical speed z' is obtained thanks to a Quadratic Bézier curves.

z ' =1- 2 z ' 0 2 1-u z  2 z ' f (26)
The values of constants z'0 and z'f are chosen in order to satisfy the constraints on initial and final values in equation (20). From ( 7) and (11) we get:

{ z ' 0= z' 0 = V z0 V g 0 z ' 1=z ' f = V zf V g s f  (27)
Where Vg(s) is the ground speed::

V g s=  V p 2 W h 2 s-2 V p W h cos   - w  (28)
In the preceding equation heading ψ is computed thanks to (12) whereas flight path angle γ is obtained thanks to (11):

  =arcsin  z ' V g   -W z   V p    (29)
Integrating (28) leads to the expression of altitude z(τ):

z = -1- 3 3 z ' 0  2 1- 2 3 u z   3 3 z ' f (30)
We choose vertical speed z' as the control variable rather than flight path angle γ because the value of z is needed to compute the true airspeed V.

Indeed as expressed in (8) true airspeed V is a function of indicated airspeed Vi and altitude z.

Thus the proposed approach consists in computing constant parameters uV and uz such that differential equation (11) holds under constraints (18) and ( 19):

find [u z u V ] T st { ∫ 0 1 t ' d = t f s f ∫ 0 1 z'  d = z f -z 0 s f (31)
This is a classical nonlinear equation solving problem which can easily be handled with numerical computational tools such as Scilab. Scilab is an open source platform for numerical computation [14].

Once the reference indicated airspeed Vi and the reference vertical profile z have been computed the feasibility of those reference profiles shall be checked in order to know if those profiles are flyable. This feasibility check is done thanks to the dynamics equations of a point mass fixed wing aircraft:

{ V = T -D m -g sin   = L-mg cos   mV (32) Where  L is the lift force L= 1 2 z S V 2 C L (33)  ρ(z)
is the air density which depends on altitude z  S is the reference wing surface area

 CL is the lift coefficient.  D is the drag force D= 1 2  z S V 2 C D ( 34
)
 CD is the drag coefficient. We will assume a parabolic drag coefficient which is specified as a function of the lift coefficient CL

C D =C D0 C D2 C L 2 (35)



F is the thrust, which is assumed to be in the opposite direction of the drag D  g is the gravitational acceleration  m is the aircraft mass Thanks to the knowledge of the reference indicated airspeed Vi and reference vertical profile z true airspeed V is known and the values of the required thrust T and required lift coefficient CL can be computed from (32) and (33):

{ T =m V m g sin   D C L = m V mg cos   0.5 h S V 2 (36)
The last step consists in checking that the required thrust T and the required lift coefficient CL are compatible with aircraft performance.

SIMULATION RESULTS

The scenario which has been simulated is the following: the aircraft starts the descent at FL 140 with an indicated airspeed (IAS) of 220 kts. The final approach fix (FAF) is assumed to be 35 NM ahead and shall be overflown at 2500 feet with an indicated airspeed of 170 kts. The maneuver duration is set to T = 540 sec (9 min). The aircraft is assumed to descent in clean configuration. The hypothetical approach path is depicted on Figure 3: the first leg is east oriented and its length is 15 NM; the second leg is north-east oriented and its length is 5 NM and the last leg is north oriented and its length is 10 NM.

Figure. 3. Hypothetical path to be followed

As far as wind is concerned we will assume no vertical wind. The horizontal wind is set to 5 kts blowing from the north in the first scenario and to 30 kts still blowing from the north in the second scenario. We have also checked that the distance flown is 30 NM for both scenario, as expected; this is a question of integrating the ground speed Vg.

CONCLUSION

This paper addresses a specific aspect of air traffic control services, namely the achievement of an orderly and expeditious flow of air traffic under time constrained continuous descent approach. This paper presents a new methodology to compute a reference indicated airspeed (IAS) and a reference vertical profile for time constrained descent which satisfy length and endpoint constraints.

The proposed approach is based on the shaping of Quadratic Bézier curves which are applied to the controlled indicated airspeed (IAS) of the aircraft and its vertical profile. An appropriate choice of the available degree of freedom when computing the reference indicated airspeed (IAS) and the vertical speed enables to solve the Two-point Boundary Value Problem (TPBVP).

Simulation results which include wind illustrate the efficiency of the proposed approach.

Future developments include the test of the robustness of the proposed design with respect to unexpected wind. This may be addressed through the periodic update of the computation of reference indicated airspeed and vertical speed as well as the design of some tracking controller installed on board the aircraft.

The proposed approach can be extended to the case where constraints over multiple fixes are imposed. Indeed, the way to compute the reference speed takes explicitly into account the initial and final values of aircraft's position and speed and can thus be used to accommodate the reference speed to each segment of flight.

Figure. 4 .

 4 Figure. 4. Reference altitude (FL) for the 5 kts North wind scenarioAs expected, the aircraft descents smoothly from FL 140 to 2500 feet following a continuous descent.The following figure presents the reference indicated airspeed (IAS in kts) with respect to path length (in NM) for the first scenario:

Figure. 5 .

 5 Figure. 5. Reference indicated airspeed (kts) for the 5 kts North wind scenarioAs expected, the indicated airspeed changes from 220 kts to 170 kts. We can see on the figure that there is a slight increase of the indicated airspeed just before the achieving 170 kts which may be disturbing for the pilot during the approach phase.The next figure depicts the ground speed Vg (in kts) with respect to path length (in NM) for the first scenario:

Figure. 9 .

 9 Figure. 9. Groundspeed (kts) for the 30 kts North wind scenario We can see on the preceding figure two discontinuities; those are obviously the effect of the change in the wind direction.