<figure><img src="image_1.png" /></figure>

HAL
open science

# Entropy minimizing curves with application to automated flight path design

Stéphane Puechmorel, Florence Nicol

► To cite this version:

Stéphane Puechmorel, Florence Nicol. Entropy minimizing curves with application to automated flight path design. GSI 2015, Oct 2015, Paris, France. hal-01162820

HAL Id: hal-01162820

https://hal.science/hal-01162820

Submitted on 15 Jun 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Entropy minimizing curves with application to automated flight path design

S. Puechmorel¹ stephane.puechmorel@enac.fr and F. Nicol nicol@recherche.enac.fr²

Ecole Nationale de l'Aviation Civile (ENAC)
Dept. SINA/MAIAA
7, Avenue Edouard Belin
F31500 TOULOUSE FRANCE

**Abstract.** Air traffic management (ATM) aims at providing companies with a safe and ideally optimal aircraft trajectory planning. Air traffic controllers act on flight paths in such a way that no pair of aircraft come closer than the regulatory separation norm. With the increase of traffic, it is expected that the system will reach its limits in a near future: a paradigm change in ATM is planned with the introduction of trajectory based operations. This paper investigate a mean of producing realistic air routes from the output of an automated trajectory design tool. For that purpose, an entropy associated with a system of curves is defined and a mean of iteratively minimizing it is presented. The network produced is suitable for use in a semi-automated ATM system with human in the loop.

## 1 Introduction

Based on recent studies [1], traffic in Europe is expected to grow on an average yearly rate of 2.6%, yielding a net increase of 2 million flights per year at the 2020 horizon. Long term forecast gives a two fold increase in 2050 over the current traffic, pointing out the need for a paradigm change in the way flights are managed. Two major framework programs, SESAR (Single European Sky Air traffic management Research) in Europe and Nextgen in the US have been launched in order to first investigate potential solutions and to deploy them in a second phase. One of the main changes that the air traffic management (ATM) system will undergo is a switch from airspace based to trajectory based operations with a delegation of the separation task to the crews. Within this framework, trajectories become the basic object of ATM, changing the way air traffic controllers will be working. In order to alleviate the workload of controllers, trajectories will be planned several weeks in advance in such a way that close encounters are minimized and ideally removed. For that purpose, several tools are currently being developed, most of them coming from the field of robotics [6]. Unfortunately, flight path issued by these algorithms are not tractable for a human controller and need to be simplified. The purpose of the present work is to introduce an automated procedure that takes as input a set of trajectories and outputs a
simplified one that can be used in an operational context. Using an entropy associated with a curves system, a gradient descent is performed in order to reduce it so as to straighten trajectories while avoiding areas with low air craft density, thus enforcing route-like behavior. This effect is related to the fact that entropy minimizing distributions favor concentration.

# 2 Entropy minimizing curves

## 2.1 Motivation

As previously mentioned, air traffic management of the future will make an intensive use of 4D trajectories as a basic object. Full automation is a far reaching concept that will probably not implemented before 2040-2050 and even in such a situation, it will be needed to keep humans in the loop so as to gain a wide societal acceptance of the concept. Starting from SESAR or Nextgen initial deployment and aiming towards this ultimate objective, a transition phase with human-system cooperation will take place. Since ATC controllers are used to a well structured network of routes, it is advisable to post-process the 4D trajectories issued by automated systems in order to make them as close as possible to line segments connecting beacons. To perform this task, in an automatic way, flight paths will be moved Add "iteratively" to dictionary so as to minimize an entropy criterion, that enforces avoidance of low density area and at the same time penalize length. Compared to already available bundling algorithms [3] that tend to move curves to high density areas, this new procedure generates geometrically correct curves, without excess curvature.

## 2.2 Entropy of a system of curves

Let a set $\gamma_1, \dots, \gamma_N$ of smooth curves be given, that will be aircraft flight paths for the air traffic application. It will be assumed in the sequel that all curves are smooth mappings from $[0, 1]$ to a domain $\Omega$ of $\mathbb{R}^2$ with everywhere non vanishing derivatives in $]0, 1[$. This last condition allows to view trajectories as smooth immersions with boundaries and is sound from the application point of view as aircraft velocities cannot be 0 expect at the endpoints. A classical performance indicator used in ATM is the aircraft density [2], obtained from sampled positions $\gamma_i(t_j), j = 1 \dots n_i$. It is constructed from a partition $U_k, k = 1 \dots P$ of $\Omega$ by counting the number of samples occurring a given $U_k$ then dividing out by the total number of samples $n = \sum_{i=1}^N n_i$. More formally, the density $d_k$ in the subset $U_k$ of $\Omega$ is:

$$d_k = n^{-1} \sum_{i=1}^{N} \sum_{j=1}^{n_i} 1_{U_k}(\gamma_i(t_j)) \quad (1)$$

with $1_{U_k}$ the characteristic function of the set $U_k$. It seems natural to extend the density obtained from samples to another one based on the trajectories
themselves using an integral form:

$$d_k = \lambda^{-1} \sum_{i=1}^{N} \int_{0}^{1} 1_{U_k}(\gamma_i(t)) dt \qquad (2)$$

where the normalizing constant $\lambda$ is obtained as:

$$\lambda = \sum_{k=1}^{P} \sum_{i=1}^{N} \int_{0}^{1} 1_{U_k}(\gamma_i(t)) dt = \sum_{i=1}^{N} \int_{0}^{1} \sum_{k=1}^{P} 1_{U_k}(\gamma_i(t)) dt$$

and since $U_k, k = 1...P$ is a partition:

$$\lambda = \sum_{i=1}^{N} \int_{0}^{1} \gamma_{i}(t) dt \qquad (3)$$

Density can be viewed as an empirical probability distribution with the $U_k$ considered as bins in an histogram. It is thus natural to extend the above computation so as to give rise to a continuous distribution on $\Omega$. For that purpose, a kernel function $K: \mathbb{R} \to \mathbb{R}^+$ is selected and a smooth version of the density [5] is defined as a mapping $d$ from $\Omega$ to $[0, 1]$:

$$d: x \mapsto \frac{\sum_{i=1}^{N} \int_{0}^{1} K (\|x - \gamma_i(t)\|) dt}{\sum_{i=1}^{N} \int_{\Omega} \int_{0}^{1} K (\|x - \gamma_i(t)\|) dt dx} \qquad (4)$$

Standard choices for the $K$ function are the ones used for non-parametric kernel estimation like the Epanechnikov function:

$$K: x \mapsto (1 - x^2) 1_{[-1,1]}(x)$$

When $K$ is compactly supported, which is the case of the Epanechnikov function and all its relatives, it comes:

$$\int_{\Omega} K (\|x - \gamma_i(t)\|) dx = \int_{\mathbb{R}^2} K (\|x\|) dx$$

provided that $\Omega$ contains the set:

$$\{x \in \mathbb{R}^2, \inf_{i=1...N, t \in [0,1]} \|x - \gamma_i(t)\| \le A\}$$

where the interval $[-A, A]$ contains the support of $K$. The case of kernels with unbounded support, like Gaussian functions, may be dealt with provided $\Omega = \mathbb{R}^2$. In the application considered, only compactly supported kernels are used, mainly to allow fast machine implementation of the density computation. Normalizing the kernel is not mandatory as the normalization occurs with the definition of $d$. It is nevertheless easier to consider only kernels satisfying:

$$\int_{\mathbb{R}^2} K (\|x\|) dx = 1$$
Using the polar coordinates $(\rho, \theta)$ and the rotation invariance of the integrand, the relation becomes:

$$2\pi \int_{\mathbb{R}^+} K(\rho)\rho d\rho = 1$$

Which yields a normalizing constant of $2/\pi$ for the Epanechnikov function, instead of the usual 3/4 in the real case. When the normalization condition is fulfilled, the expression of the density simplifies to:

$$d: x \mapsto N^{-1} \sum_{i=1}^{N} \int_{0}^{1} K(\|x - \gamma_{i}(t)\|) dt \quad (5)$$

As an example, one day of traffic over France is considered and pictured on figure 2.2 with the corresponding density map, computed on a evenly spaced grid with a normalized Epanechnikov kernel:

<figure><img src="image_3.png" /><figcaption>Fig. 1. Traffic over France and associated density</figcaption></figure>

Unfortunately, density computed this way suffers a severe flaw for the ATM application: it is not related to the shape of trajectories but also to the time behavior. Formally, it is defined on the set **Imm**([0, 1], $\mathbb{R}^2$) of smooth immersions from [0, 1] to $\mathbb{R}^2$ while the space of primary interest will be the quotient by smooth diffeomorphisms of the interval [0, 1], **Imm**([0, 1], $\mathbb{R}^2$) / **Diff**([0, 1]). Invariance of the density under the action of **Diff**([0, 1]) is obtained as in [4] by adding a term related to velocity in the integrals. The new definition of d becomes:

$$\tilde{d}: x \mapsto \frac{\sum_{i=1}^{N} \int_{0}^{1} K(\|x - \gamma_{i}(t)\|) \|\gamma'_{i}(t)\| dt}{\sum_{i=1}^{N} \int_{\Omega} \int_{0}^{1} K(\|x - \gamma_{i}(t)\|) \|\gamma'_{i}(t)\| dtdx} \quad (6)$$

Assuming again a normalized kernel and letting $l_i$ be the length of the curve $\gamma_i$, is simplified to:

$$\tilde{d}: x \mapsto \frac{\sum_{i=1}^{N} \int_{0}^{1} K(\|x - \gamma_{i}(t)\|) \|\gamma'_{i}(t)\| dt}{\sum_{i=1}^{N} l_{i}} \quad (7)$$
The new **Diff**-invariant density is pictured on figure 2.2 along with the standard density. While the overall aspect of the plot is similar, one can observe that routes are more apparent the right picture and that the density peak located above Paris area is of less importance and less symmetric is due to the fact that near airports, aircraft are slowing down and this effect exaggerates the density with the non-invariant definition.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Density (left) and **Diff** invariant density (right) for the 12th February 2013 traffic</figcaption></figure>

Having a density at hand, the entropy of the system of curves $γ_1, ..., γ_N$ is defined the usual way as:

$$E(γ_1, ..., γ_N) = - \int_{Ω} \tilde{d}(x) \log (\tilde{d}(x)) dx$$

## 2.3 Minimizing the entropy

In order to fulfill the initial requirement of finding bundles of curve segments as straight as possible, one seeks after the system of curves minimizing the entropy $E(γ_1, ..., γ_N)$, or equivalently maximizing:

$$\int_{Ω} \tilde{d}(x) \log (\tilde{d}(x)) dx$$

The reason why this criterion gives the expected behavior will become more apparent after derivation of its gradient at the end of this part. Nevertheless, when considering a single trajectory its is intuitive that the most concentrated density distribution is obtained with a straight segment connecting the endpoints.

Letting $ϵ$ be a perturbation of the curve $γ_j$ such that $ϵ(0) = ϵ(1) = 0$, the first order expansion of $-E(γ_1, ..., γ_N)$ will be computed in order to get a maximizing displacement field, analogous to a gradient ascent<fnref n="1" /> in the finite

<footnote n="1">¹ Choice has been made to maximize the opposite of the entropy, so that the algorithm will be a gradient ascent one</footnote>
dimensional setting. The notation:

$$ \frac{\partial F}{\partial \gamma_j} $$

will be used in the sequel to denote the derivative of a function F of the curve $\gamma_j$ in the sense that for a perturbation $\epsilon$:

$$ F(\gamma_j + \epsilon) = F(\gamma_j) + \frac{\partial F}{\partial \gamma_j}(\epsilon) + o(\|\epsilon\|_2) $$

First of all, please note that since $\tilde{d}$ has integral 1 over the domain $\Omega$:

$$ \int_{\Omega} \frac{\partial \tilde{d}(x)}{\partial \gamma_j}(\epsilon) dx = 0 $$

so that:

$$ -\frac{\partial}{\partial \gamma_j} E(\gamma_1, \dots, \gamma_N)(\epsilon) = \int_{\Omega} \frac{\partial \tilde{d}(x)}{\partial \gamma_j}(\epsilon) \log(\tilde{d}(x)) dx \quad (8) $$

Starting from the expression of $\tilde{d}$ given in equation 7, the first order expansion of $\tilde{d}$ with respect to the perturbation $\epsilon$ of $\gamma_j$ is obtained as a sum of a term coming from the numerator:

$$ \int_{0}^{1} K (\|x - \gamma_j(t)\|) \|\gamma_j'(t)\| dt \quad (9) $$

and a second one coming from the length of $\gamma_j$ in the denominator. This last term is obtained from the usual first order variation formula of a curve length:

$$ \int_{[0,1]} \|\gamma_j'(t) + \epsilon'(t)\| dt = \int_{[0,1]} \|\gamma_j'(t)\| dt + \int_{[0,1]} \left\langle \frac{\gamma_j'(t)}{\|\gamma_j'(t)\|}, \epsilon'(t) \right\rangle dt + o(\|\epsilon\|_2) $$

Using an integration by part, the first order term can be written as:

$$ \int_{[0,1]} \left\langle \frac{\gamma_j'(t)}{\|\gamma_j'(t)\|}, \epsilon'(t) \right\rangle dt = \quad (10) $$

$$ - \int_{[0,1]} \left\langle \left( \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right)_N, \epsilon \right\rangle dt \quad (11) $$

with:

$$ \left( \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right)_N = \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} - \frac{\gamma_j'(t)}{\|\gamma_j'(t)\|} \left\langle \frac{\gamma_j'(t)}{\|\gamma_j'(t)\|}, \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right\rangle $$

the normal component of:

$$ \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} $$
for a curve in $\mathbb{R}^2$.

The integral in 9 can be expanded in a similar fashion. Using as above the notation $()_N$ for normal components, the first order term is obtained as:

$$ \int_{[0,1]} \left\langle \left( \frac{\gamma_j(t) - x}{\|\gamma_j(t) - x\|} \right)_N, \epsilon \right\rangle K' (\|\gamma_j(t) - x\|) \|\gamma'_j(t)\| dt \quad (12) $$

$$ - \int_{[0,1]} \left\langle \left( \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right)_N, \epsilon \right\rangle K (\|\gamma_j(t) - x\|) dt \quad (13) $$

From the expressions in 12 and 10, the first order variation of the entropy is:

$$ \frac{1}{\sum_{i=1}^{N} l_i} \left( \int_{[0,1]} \left\langle \int_{\Omega} \left( \frac{\gamma_j(t) - x}{\|\gamma_j(t) - x\|} \right)_N K' (\|\gamma_j(t) - x\|) \log \tilde{d}(x) dx, \epsilon \right\rangle \|\gamma'_j(t)\| dt \right) \quad (14) $$

$$ \int_{[0,1]} \left\langle \int_{\Omega} \left( \frac{\gamma_j(t) - x}{\|\gamma_j(t) - x\|} \right)_N K' (\|\gamma_j(t) - x\|) \log \tilde{d}(x) dx, \epsilon \right\rangle \|\gamma'_j(t)\| dt \quad (15) $$

$$ - \int_{[0,1]} \left( \int_{\Omega} K (\|\gamma_j(t) - x\|) \log \tilde{d}(x) dx \right) \left\langle \left( \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right)_N, \epsilon \right\rangle dt \quad (16) $$

$$ + \left( \int_{\Omega} \tilde{d}(x) \log(\tilde{d}(x)) dx \right) \int_{[0,1]} \left\langle \left( \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right)_N, \epsilon \right\rangle dt \quad (17) $$

As expected, only moves normal to the trajectory will change at first order the value of the criterion: the displacement of the curve $\gamma_j$ will thus be performed at $t$ along the normal vector $N_{\gamma_j}(t)$ and is given, up to the $(\sum_{i=1}^{N} l_i)^{-1}$ term by:

$$ \int_{\Omega} \left( \frac{\gamma_j(t) - x}{\|\gamma_j(t) - x\|} \right)_N K' (\|\gamma_j(t) - x\|) \log \tilde{d}(x) dx \|\gamma'_j(t)\| \quad (18) $$

$$ - \left( \int_{\Omega} K (\|\gamma_j(t) - x\|) \log \tilde{d}(x) dx \right) \left( \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right)_N \quad (19) $$

$$ + \left( \int_{\Omega} \tilde{d}(x) \log(\tilde{d}(x)) dx \right) \left( \frac{\gamma_j''(t)}{\|\gamma_j'(t)\|} \right)_N \quad (20) $$

The first term in the expression will favor moves towards areas of high density, while the second and third one are moving along normal vector and will straighten trajectory.

In practical implementation, the scaling factor in front of the whole expression is dropped and moves are made proportionally to the given vector. As usual with gradient ascent algorithms, one must carefully select the step taking in the maximizing direction in order to avoid divergence. A simple fixed step strategy was selected here and gives satisfactory results.

In the example of figure 3, the problem of automatic conflict solving is addressed. In the initial situation, aircraft are converging to a single point, which is unsafe. Air traffic controllers will proceed in such a case by diverting aircraft
from their initial flight path so as to avoid each other, but only using very simple maneuvers. An automated tool will make a full use of the available airspace and the resulting set of trajectories may fail to be manageable by a human: in the event of a system failure, no backup can be provided by controllers. The entropy minimization procedure was added to an automated conflict solver in order to end up with flight paths still tractable by humans. The final result is shown on the right part of figure 3, where encounters no longer exists but aircraft are bound to simple trajectories, with a merging and a splitting point. Note that since the automated planner acts on velocity, all aircraft are separated in time on the inner part.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Initial flight plans and final ones</figcaption></figure>

# 3 Conclusion and future work

Algorithms coming from the field of shape spaces emerge as a valuable tool for applications in ATM. In this work, the foundations of a post processing procedure that may be applied after an automated flight path planner are presented. Entropy minimization makes straight segments bundle emerge, which fulfills the operational requirements. Computational efficiency has to be improved in order to release an usable building block for future ATM systems. One way to address this issue is to compute kernel density estimators using GPUs which excel in this kind of task, very similar texture manipulations. Furthermore, theoretical insights have to be gained in the next step of the work.

# References

1. EUROCONTROL/NMD/STATFOR. Eurocontrol seven-year forecast. https://www.eurocontrol.int/sites/default/files/content/documents/official-documents/forecasts/seven-year-flights-service-units-forecast-2014-2020-sep2014.pdf, 2014.
2. W.H. Harman. Air traffic density and distribution measurements. Lincoln Laboratory, MIT, Report ATC-80, May 1979.

3. C. Hurter, O. Ersoy, and A. Telea. Smooth bundling of large streaming and sequence graphs. In *Visualization Symposium (PacificVis), 2013 IEEE Pacific*, pages 41–48, Feb 2013.

4. Peter W Michor and David Mumford. Riemannian geometries on spaces of plane curves. *J. Eur. Math. Soc. (JEMS)*, 8:1–48, 2006.

5. Emanuel Parzen. On estimation of a probability density function and mode. *Ann. Math. Statist.*, 33(3):1065–1076, 09 1962.

6. G.P. Roussos, D.V. Dimarogonas, and K.J. Kyriakopoulos. Distributed 3d navigation and collision avoidance for nonholonomic aircraft-like vehicles. In *Proceedings of the 2009 European Control Conference*, 2009.