ENRI International Workshop on ATM/CNS. Tokyo, Japan. (EIWAC 2009).

New Trends in Air Trafﬁc Complexity

Daniel Delahaye
Ecole Nationale Aviation Civile
Applied Math Laboratory
7 Avenue Edouard Belin
31055 TOULOUSE (FRANCE)
delahaye@recherche.enac.fr

Stéphane Puechmorel
Ecole Nationale Aviation Civile
Applied Math Laboratory
7 Avenue Edouard Belin
31055 TOULOUSE (FRANCE)
puechmor@recherche.enac.fr

Abstract

cept is ﬁxed on a per-sector basis and complexity is assessed by comparing the real number of aircraft with thesector capacity. It must be noted that under some circumstances controllers will accept aircraft beyond thecapacity threshold while rejecting trafﬁc at other times
although the number of aircraft is well below the maximum capacity. This simple fact clearly show that capacity as a crude complexity metric is not enough by itself
to fully account for the controller’s workload. In order
to better quantify the complexity, geometric features of
the trafﬁc have to be included.

This paper describes a novel approach for deﬁning an
air trafﬁc complexity indicator, aimed at improving overthe current operational deﬁnition that is only the number of aircraft present in a sector at a given time. The
key concept in this work is that a level of disorder or
equivalently unpredictibility is the right indicator. Our
approach is relevant for many applications in the air
trafﬁc management area, like sector design ([2]), trafﬁc assignment ([1, 3]) and will be very well suited forpredicting complex situations in a future autonomous
aircraft environment. Moreover, since this modellingtakes explicitely into account the trajectories of the aircraft, either observed or planned, it ﬁts perfectly with
4D based ATM like SESAR or NextGen.

1.1. Complexity vs workload
As previously stated, depending on the trafﬁc structure, ATC controllers will perceive differently situations, even if the number of aircraft present in the sector
is the same. Furthermore, exogenous paramaters likethe workload history can be inﬂuential on the perceived
complexity at a given time (a long period of heavy load
will tend to reduce the efﬁciency of a controller). Some
reviews of complexity in ATC have been completed,mainly from the controller’s workload point of view
[5, 11], and have recognized that complexity is related
to both the structure of the trafﬁc and the geometry of
the airspace. This tends to prove that controller’s workload has two facets :1. Introduction to complexity
The ATM system has to cope with an increasing
number of ﬂights, pushing the capacity to its limits. As
an example, the average daily trafﬁc above Europe was
26286 ﬂights/day, with a peak trafﬁc demand in excess
of 31000 ﬂights [4]. Although delays are kept low, it isexpected from the same reference that capacity has to beextended in the furture. Basically, two strategies can bedevised : adapt the demand to capacity (slot-route allocation, collaborative decision making, . . . ) or adapt thecapacity to the demand (Airspace design, 4D trajectory
planning, autonomous aircraft, . . . ). The ﬁrst approach
can be used in the context of current ATM system, while
innovative future designs will mainly follow the second
strategy.
Currently, complexity of the trafﬁc is measured
only as an operational capacity : the maximum number of aircraft that ATC controllers are willing to ac-• An intrinsic complexity related to trafﬁc structure.
• A human factor aspect related to the controller itself.
While most complexity metrics tend to capture those effects within a single aggregate indicator, the purpose ofthis work is to design a measure of intrinsic complexity
only since it is the most relevant metric for an highly
automated ATC system (no human factors).

55

with ti , xi , vi respectively the measurement time, the position and the velocity of the i-th aircraft, the interpolation problem is to ﬁnd a smooth enough vector ﬁeld :X : R × R3 → R3 such that :1.2. Dynamic density
The ﬁrst complexity indicator incorporating structural considerations along with the simple number ofaircraft is the “Dynamic Density” of Laudeman etal. from NASA [10]. The “Dynamic Density” is aweighted sum of the trafﬁc density (number of aircraft),the number of heading changes (> 15 degrees), the number of speed changes (>0.02 Mach), the number of altitude changes (>750 ft), the number of aircraft with 3D Euclidean distance between 0-25 nautical miles, thenumber of conﬂicts predicted in 25-40 nautical miles.
These factors are summed together using weighting factors that were determined by showing different trafﬁc
scenarios to several controllers. B.Sridhar from NASA[12], has developed a model to predict the evolution ofsuch a metric in the near future. Efforts to deﬁne “Dynamic Density” have identiﬁed the importance of a wide
range of potential complexity factors, including structural considerations. However, the instantaneous position and speeds of the trafﬁc itself does not appear to beenough to describe the total complexity associated withan airspace. A few previous studies have attempted toinclude structural consideration in complexity metrics,
but have done so only to a restricted degree. For example, the Wyndemere Corporation proposed a metricthat included a term based on the relationship between
aircraft headings and dominant geometric axis in a sector [8]. The importance of including structural consideration has been explicitly identiﬁed in work at Eurocontrol. In a study to identify complexity factors using
judgment analysis, Airspace Design was identiﬁed as
the second most important factor behind trafﬁc volume[9]. Histon, et, al. [6, 7] investigated how this structure
can be used to support structure-based abstractions that
controllers appear to use to simplify trafﬁc situations.

∀i = 1 . . . N, X(ti , xi ) = vi
Unfortunately, this problem is ill-posed since this requirement will not uniquely deﬁne X. A classical way
to deal with this problem is to add an extra requirement
on X, namely that it realizes the minimum of a so-called
energy functional. Most of the time, this functional is
obtained as a spatial L2 -norm of the form :R3

LX(t, x)2 dx

with L a differential operator. It has to be noted thatvariable t is not taken into account, so that the solution
is expected to be a stationary vector ﬁeld. A lot of works
have been devoted to this kind of problem that falls in
the category of interpolating vector splines ([14],[13]).
It is well known that the optimal X satisfying the interpolation condition is of the form :NX(t, x) = ∑ λit G(x, xi )i=1

Where G(x, y) is the Green’s function associated withthe differential operator Lt L and λi is a 3-dimensional
vector coefﬁcient. A convenient L for many applications in physics is the div-curl operator :L = α∇div + β ∇curlwith α, β positive real numbers tuning the respective
contributions of the div and the curl part : a high α
will yield to a nearly constant divergence while a highβ will enforce constant curl. It has to be noted that if
α = β , the criterion degenerates to a square laplacian
operator, which has the nice property of being diagonal,
so that the vector problem reduces to 3 one dimensional
problems. While this particular approach has been exploited in some of our past works, it is more relevant
for our application to allow time-varying interpolating
vector ﬁelds : in the following X : R × R3 will be a timedependent vector ﬁeld with time as ﬁrst coordinate. Theenergy functional is chosen to be :


  
2
∂X 
E(X) =
+
R R3  ∂t 

2. Dynamical System Modeling
The main idea behind this new complexity indicator is to consider that observed trafﬁc is a sample ofan underlying ﬂow. Within this framework, aircraft canbe thought as particles moving in a stream, the geometry of the corresponding system being the most inﬂuencial factor for complexity. While this modelling is
not directly related to a physical reality it can be noted
however that in an highly automated context, the manoeuvres generated by an automatic conﬂict solver will
adhere to this principle quite accurately. The ﬁrst stepfor complexity computation is thus to ﬁnd an interpolating vector ﬁeld based on the past aircraft positions andintents. Given a set of N observations :μ ∇(div + curl)X2 dtdx
where μ > 0 tunes the relative part of time variation anddiv-curl criterion. Intuitively, X will be the most stationary ﬁeld with minimal div-curl variation. The differential operator associated with this variational problem is(ti , xi , vi )Ni=1

56

Like the heat kernel, the fundamental solution obtained
is singular for x = 0,t = 0. Before introducing a way ofobtaining computable kernels, it is interesting to look atsome properties of p(t, x) :given by :∂2
P = 2 − μΔ2∂t
An elementary solution can be found in S  as :√1exp(−|t| μ|ξ |2 )pf √μξ 2• Error function is rapidly converging to 1 when its
argument goes to +∞. In fact, er f (4) is equal to 1
at machine precision when computing with single
precision ﬂoat numbers. It thus means that contributions in the interpolating ﬁeld coming from aircraft far from the evaluation point are very close
to being proportional to a standard 1/x potentialfunction.

where the previous expression is the fourier transformin S  , the space a rapidly decreasing smooth functions,with respect to the spatial coordinates of the elementary
solution. It turns out that this distribution is indeed ob1 (integrable on compact sets) mapping
tained from a Lloc
so that the inverse Fourier transform can be obtained
readily by integration as :1p(t, x) = 3
8π• Unlike classical div − curl splines obtained fromthe ﬁrst non time-varying functional, these new
splines are decreasing at inﬁnity, allowing to drop
contributions of far enough aircraft (strictly speaking, this is the case even for the div-curl splines,
but decrease is due to cancellations between contributions and is slower).

1√μξ 2R3√exp(−|t| μξ 2 )exp(ix, ξ )dξby Fubini’s theorem and polar change of variables it becomes :p(t, x) =


8π

S2

1
3√



μ

R

√
exp(−|t| μr2 )The main issue in using the kernel p is that the reconstructed ﬁeld :exp(ix, rs)dσ (s)drwith dσ the solid angle measure.
parametrization of the unit sphere :NX(t, x) = ∑ λi p(t − ti , x − xi )Using a polar

i=1

is singular at observation points (ti , xi , vi ). A simple wayof avoiding this is to change the interpolation criterion :instead of enforcing that X(ti , xi ) = vi , i = 1 . . . N, we impose that the mean value of the ﬁeld in a neighboroodof xi has to be equal to vi . In order to allow simple computation of the interpolating splines in closed-form, it
has been chosen to compute the mean of X with respect
to a gaussian density, so that the interpolation condition
becomes :p(t, x) =√1exp(−|t| μr2 )√38π μ R 2π  π00exp(ixr cos θ ) sin θ dθ dφ drand ﬁnally :1√
2π 2 μ√
sin xr
exp(−|t| μr2 )
dr
xr
R

p(t, x) =1(2πσ )3/2= vi

by Parseval equality :−x−xi 22σ 2X(t, x)dxp(t, x) =

σ
x

er f
x
σ 2 + |t|so in terms of error function erf :1
erf
4π 2 xR3e

σ is the standard deviation for the gaussian density andis a tuning parameter for the shape of the ﬁeld. Solving the functional equation with this new interpolation
criterion yields a non singular kernel :


π
1
√
2
4xπ
|t| μ


 x
2
2
2π
π ω
dω
exp − √
x
|t| μ
− 2π

p(t, x) =p(t, x) =



x
 √
2 |t| μHere again, it worth notice that the asymptotic behaviour of p is to revert to a classical 1/x potential.

57

3. Lyapunov exponents

Lyapunov exponents over the time (assuming samplingtimes (t1 , . . . ,tn ) :The metric chosen for complexity computation relies on a measure of sensitivity to initial conditions of
the underlying dynamical system called Lyapunov exponents. In order to ﬁgure out what Lyapunov exponents are, let consider a point and look at its evolution
when transported by the dynamical system.
Let x be ﬁxed (initial point) and let φ be a point trajectory of the dynamical system associated to the vector
ﬁeld X given by :φ (t, x) = x + t0X(s, φ (s, x))ds

κ(x) =

1 i=n∑ Dx X(φ (ti , x))2
n i=1

(3.3)

The ﬁgure 1 shows an example of Lyapunov exponents map for which full organized miles in trail trajectories (from south west to north east) cross two random
trafﬁc situations. This ﬁgure shows clearly a complexity(3.1)

Assume now that trajectory is disturbed by a small
perturbation ε, we have :φ (t, x + ε) = φ (t, x)+Dx φ (t, x) · ε + o(ε)
where Dx φ (t, x) is the differential of the vector ﬁeld at
x that satisﬁes :Dx φ (t, x)=dt
Dx X(φ (t, x)) · Dx φ (t, x)

Figure 1. Miles in trail traﬁc between disordered areas(3.2)

valley on the miles in trail direction. This organization
may have been detected even if the miles in trail trajectories would have been structured on a curve trajectory.
That is the strong point of this metric: Lyapunov exponents are able to identify any kind of trajectory organization. The next ﬁgure is an example of a complexity
map computed on real trafﬁc over France : Major air-

The Lyapunov exponents are closely related to the singular values of the matrix Dx φ (t, x) and can be thought
as local shear values for the dynamical system.
When Lyapunov exponents are high, the trajectoryof a point under the action of the dynamical system isvery sensitive to initial conditions (or to the parameterson which the vector ﬁeld may depend), so that situation in the future is unpredictable. On the other hand,small values of the Lyapunov exponents mean that the
future is highly predictable (expected to be comfortable
for a controller). So, the Lyapunov exponent map determines the area where the underlying dynamical system
is organized. It identiﬁes the places where the relative
distances between aircraft do not change with time (lowreal value) and the ones where such distance change alot (hight real value).
Let us now describe the practical procedure forcomputing complexity maps.
First of all, the optimal dynamic div-curl approximation for the observed trajectories is computed, basedon the deﬁning equations. That step requires a linear
system solving.
The second step computes the second derivatives
matrix Dx φ at each point of the grid for φ trajectory
starting at x. This is done by solving the differential
equation 3.2 with a Runge-Kutta integrator. The complexity value at point x is then obtained by averaging

Figure 2. Complexity over France
ports are clearly visible as high complexity areas. It has
to be noted that on this picture, organized trafﬁc is classiﬁed as low complexity as expected.

58

4. Computational issues

largest Lyapunov exponents : the matrix A(t, x) is factored as A(t, x) = Q(t)R(t) with Q(t) orthogonal andR(t) upper triangular, both matrices beeing smooth int. The Lyapunov exponents are directely related to the
diagonal elements of R(t). The problem is then to update continuously a QR factorisation (or a polar decomposition which is very similar except that the right handside matrix is symetric instead of being upper triangular). The differentiel equation satisﬁed by Q can be established ﬁrst by noticing that the matrix :Computing Lyapunov exponents amounts to integrate a differential equation. Under standard smoothness assumptions on X :dDx φ (t, x) = Dx X(t, φ (t, x))Dx φ (t, x)
dt
giving a linear cocyle equation. Putting
A(t, x) = Dx φ (t, x)

M(t) = Qh (t)A(t)Q(t) − Qh (t)Q̇(t)the cocyle equation is :is upper triangular so adding its conjugate yields a symetric matrix with the same coefﬁcients except on the
diagonal where there are doubled. However :dA(t, x) = Dx X(t, φ (t, x))A(t, x)
dtwhich is a linear differential equation in matrix function A(t, x). Nearly all standard algorithms can be used
for that purpose but one quickly realize that the problem is far from being well conditionned. In fact, mostof the time it is easy to obtain the ﬂow φ (t, x) with good
accuracy but by construction A(t, x) tends to grow exponentially fast in some directions (corresponding to positive Lyapunov exponents) and to decay exponentiallyfast in others (corresponding to negative Lyapunov exponents) : the condition number is thus increasing again
exponentially. Then it will be extremely inaccurate tocompute Lyapunov exponent merely by integrating thelinear cocyle equation : one needs a kind of rescaling
to recover good numerical properties. There is a abundant literature on the subject, however all methods fall
into one of two categories : spatial integration or temporal integration. Direct application of the deﬁnition
gives the second approach while ergodic theorem gives
the ﬁrst one. Both have advantages and drawbacks :M h (t) + M(t) = Qh (t)(A(t) + Ah (t))Q(t)since Qh (t)Q̇(t) + Q̇h (t)Q(t) = 0 (Q is unitary so thederivative of Qh (t)Q(t) = Id vanishes). So M(t) can beotained without Q̇(t). Now, simply use the relation :Q̇(t) = A(t)Q(t) − Q(t)M(t)This differential equation can be solved readily by standard runge-kutta integrator. Orthogonality of Q is nevertheless hard to preserve, so that efﬁcient implementations require specﬁc algorithms. A natural approachis to use integrators working on the Stiefel manifold on
which Q lives. Another one is to use a parametrizationof Q using elementary rotations (we need 6 such rotations for a 3x3 matrix). Regardless of the method used,
computing a complexity map on large airspaces is stilla challenge. Many improvements have been made thatallow now to produce country-sized complexity timein about 10 minutes with several thousands aircraft involved. Investigation on paralell computing and more
speciﬁcally on GPU computing is ongoing : the expected increase in computation speed will probably allow real-time computations even for relatively large airspaces.

• Spatial integration is efﬁcient and free from slow
convergence phenomenon occuring sometimes in
temporal integration. However, in order to use theergodic theorem, one must ﬁnd an invariant measure. Most of the time, it has to be done by covering algorithms.
• Temporal integration can suffer from slow convergence. Moreover, it requires periodic rescaling to
avoid numerical problems.

References[1] D Delahaye, J.L Farges, M Schoenauer, and J.M Alliot,
Genetic algotithms for air trafﬁc assignment, Proceedings of the European Conference on Artiﬁcial Intelligence, ECAI, 1994.
, Genetic algotithms for partitioning airspace,[2]
Proceedings of the 10th Conference on Artiﬁcial Intelligence for Applications, IEEE, 1994.
[3] D Delahaye and A Odoni, Aispace congestion smoothing by stochastic optimization, Proceedings of the Evolutionnary Programming Conference, IEEE, 1997.

In our complexity application, only the temporal approach has been tested.

4.1. Algorithms based on differential geometry
All the machinery used in this part is borrowed
from a paper from T.J. Bridges and S.Reich. Basically,
the trick is to use the Perron lemma to compute the k59[4] Eurocontrol, Atfcm and capacity report 2006, Eurocontrol, 2006.
[5] B. Hilburn, Cognitive complexity in air trafﬁc control :A literature review, Tech. Report 04/04, EEC, 2004.
[6] J. Histon, G. aigoin, D. Delahaye, R.J. Hansman,and S. Puechmorel, Introducing structural considerations into complexity metrics, Proceedings of the Fourth
USA/EUROPE ATM R&D Seminar (Eurocontrol-FAA,ed.), 2001.
[7] J. Histon, R.J. Hansman ans G. Gottlieb, H. Kleinwaks,
S. Yenson, D. Delahaye, and S. Puechmorel, Structural considerations and cognitive complexity in air trafﬁc control, Proceedings of the 21th Air Trafﬁc Management for Commercial and Military Systems (IEEEAIAA, ed.), 2002.
[8] Windemere inc, An evaluation of air trafﬁc control complexity, Tech. Report 2-14284, NASA, 1996.
[9] B. Kriwan, R. Scaife, and R. Kennedy, Investigatingcomplexity factors in uk air trafﬁc management, Human
Factors and Aerospace Safety, 2001.
[10] I.V. Laudeman, S.G. Shelden, R. Branstrom, and C.L.
Brasil, Dynamic density : an air trafﬁc management
metric, Tech. Report TM-1998-112226, NASA, 1998.
[11] R.H. Mogford, J.A. Guttman, and S.L. Morrowand P.
Kopardekar, The complexity construct in air trafﬁc control : A review and synthesis of the literature, Tech. Report DOT/FAA/CT-95/22, FAA, 1995.
[12] B. Shridar, K.S. Seth, and S. Grabbe, Airspace complexity and its application in air trafﬁc management, Proceedings of the Fourth USA/EUROPE ATM R&D Seminar (Eurocontrol-FAA, ed.), 2001.
[13] D. Suter, Motion estimation and vector splines.
[14] G. Wahba, Spline models for observational data, SIAM
(1990).

Copyright“Copyright Statement The authors conﬁrm that
they, and/or their company or institution, hold copyrightof all original material included in their paper. They
also conﬁrm they have obtained permission, from the
copyright holder of any third party material included intheir paper, to publish it as part of their paper. The authors grant full permission for the publication and distribution of their paper as part of the EIWAC 2009 proceedings or as individual off-prints from the proceedings.”

60