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

HAL
open science

# Functional interior point programming applied to the aircraft path planning problem

Stéphane Puechmorel, Daniel Delahaye

► To cite this version:

Stéphane Puechmorel, Daniel Delahaye. Functional interior point programming applied to the aircraft path planning problem. EuroGNC, European Aerospace Guidance, Navigation and Control Conference 2015, Apr 2015, Toulouse, France. 2015, Apr 2015, Toulouse, France. 10.1007/978-3-319-17518-8_30. hal-01113889

HAL Id: hal-01113889

https://hal.science/hal-01113889

Submitted on 6 Feb 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
Functional interior point programming applied
to the aircraft path planning problem

Puechmorel, S. and Delahaye, D.

**Abstract** Multiple aircraft trajectory planning is a central problem in future air traf-
fic management concepts where some part of the separation task, currently assumed
by human controllers, will be delegated to on-board automated systems. Several
approaches have been taken to address it and fall within two categories: meta-
heuristic algorithms or deterministic methods. The framework proposed here mod-
els the planning problem as a optimization program in a space of functions with
constraints obtained by semi-infinite programming. A specially designed innovative
interior point algorithm is used to solve it.

# 1 Introduction

The path planning problem for multiple robots evolving in a possibly dynamic en-
vironment with obstacles is a very active area of research among the automatics
and robotics communities, and is addressed in many classical references [6, 7]. For
application to aircraft trajectories, nonholonomic constraints have to be considered
since bounds on velocity and curvature are dictated by flight dynamics and aircraft
operations. In such a context, even for two dimensional motion, it is known that an
exact shortest path computation for a single mobile avoiding polyhedral obstacles
is NP-Hard [3]. In order to make the problem tractable, only approximate solutions
are to be sought after. One possible way of dealing with the intrinsic complexity
is to resort to solvers based on metaheuristics. Since the state space exploration
is made on a random fashion, there is no insurance on the quality of the solution
obtained after a finite number of iterations and it is not even possible for some al-
gorithms to prove convergence in a probabilistic sense. However, due to the nature

<footnote>Delahaye, Daniel.
ENAC, 7 avenue Edouard Belin, 31055 Toulouse, e-mail: delahaye@recherche.enac.fr</footnote>

<footnote>Puechmorel, Stephane.
ENAC, 7 avenue Edouard Belin, 31055 Toulouse e-mail: stephane.puechmorel@enac.fr</footnote>
of the metaheuristic approach, a wide class of constraints and optimality criteria
can be taken into account. In the field of air traffic management, one can impose
that the solution is built using only maneuvers used by human controllers. Suc-
cessful applications of stochastic optimization algorithms can be found in [4, 5].
The second way of addressing the issue of computational complexity is to allow
sub-optimality of the solution. Turning the original problem into a continuous op-
timization program under constraints allows traditional algorithms to be used and
provides a locally optimal solution. As the NP-hardness of the original formulation
cannot be avoided, global optimality cannot be reached generally, unless some kind
of exhaustive search among all local optima is performed. Hybridization between
the two approaches can be done in order to randomly sample the set of local op-
tima, yielding in turn an improved admissible solution. The algorithm that will be
presented here falls in the second category, and is built using a special geometrical
property of complex configuration spaces. The path planning program is designed in
a functional way, with criterion based on a cost associated to individual trajectories
and constraints expressed using semi-infinite programming. The paper is organized
as follows: in a first section, general results on complex configuration spaces will
be briefly recapped, as they will provide the theoretical framework on which the
interior point formulation will be built. The second section will introduce the op-
timization program associated to the path planning problem, and a relation with a
recently introduced complex harmonic navigation function will be pointed out. Fi-
nally, a possible algorithmic implementation of the solver yielding an approximate
solution to the original problem will be detailed.

## 2 Complex configuration spaces

Given a set of $N$ mobiles with planar motion, the complex configuration space $\mathcal{C}^N$ is defined to be the Cartesian product of $n$ copies of the complex plane $\mathbb{C}$ with the set $\Delta_N = \{(z_1, \dots, z_N), \exists i \neq j, z_i = z_j\}$ of simultaneous positions removed. Any path connecting two points in $\mathcal{C}^N$ is an admissible collision-free planning of $N$ trajectories. It is easy to show that $\mathcal{C}^N$ is a path connected space using a sequential planning argument: let $(z_1^0, \dots, z_N^0)$ be any initial configuration and let $(z_1^1, \dots, z_N^1)$ be the desired final situation. Let $z_1, \dots, z_N$ be fixed. Then it exists a path joining $z_1^0$ and $z_1^1$ in $\mathbb{C} - \{z_1, \dots, z_N\}$ since this space is path-connected. The same procedure can then be applied iteratively to the remaining points. The resulting path is a sequence of moves along coordinates axis in $\mathcal{C}^n$ and assumes that only one mobile moves at a time.

Geometrical insights about $\mathcal{C}^N$ can be gained from the knowledge of its coho-
mology group, obtained in [2]. It turns out that the generators of this group will be
used to define the constraints in the optimization program constructed later. They
consist of degree one forms:
$$ \omega_{k,l} = \frac{dz_k - dz_l}{z_k - z_l}, \quad 1 \le k \le N, 1 \le l \le N, k \ne l $$

Please note that in the original paper [2] a factor $(i2\pi)^{-1}$ appears in the expression: it is removed here as the real part of $\omega_{k,l}$ is of primary interest. Let a $C^1$ path $\Gamma: [0, 1] \to \mathcal{C}^N$ be given. Then:

$$ \int_{\Gamma} d\omega_{k,l} = \int_{[0,1]} \frac{\Gamma'_{k}(t) - \Gamma'_{l}(t)}{\Gamma_{k}(t) - \Gamma_{l}(t)} dt $$

with $\Gamma_j, j = 1...N$ denoting the j-th component of $\Gamma$. Rewriting the right hand side, it comes:

$$ \int_{\Gamma} d\omega_{k,l} = \int_{[0,1]} \frac{(\Gamma'_{k}(t) - \Gamma'_{l}(t)) (\overline{\Gamma_{k}(t)} - \overline{\Gamma_{l}(t)})}{|\Gamma_{k}(t) - \Gamma_{l}(t)|^2} dt $$

Given any two complex numbers $z_1 = x_1 + iy_1, z_2 = x_2 + iy_2$, the real part of the product $\overline{z_1}z_2$ is the inner product of the vectors $(x_1, y_1)$ and $(x_2, y_2)$ and imaginary part their determinant. The expression:

$$ \frac{(\Gamma'_{k}(t) - \Gamma'_{l}(t)) (\overline{\Gamma_{k}(t)} - \overline{\Gamma_{l}(t)})}{|\Gamma_{k}(t) - \Gamma_{l}(t)|^2} $$

thus admits an interpretation as an expansion rate (resp. rotation rate) for the vector $d_{kl}(t) = \Gamma_k(t) - \Gamma_l(t)$ and in turn the integral:

$$ \int_{\Gamma} d\omega_{k,l} $$

will give $\log(||d_{kl}(1)||) - \log(||d_{kl}(0)||)$ as its real part and $2\pi\theta$ for its imaginary part, with $\theta$ the winding number of the path $d_{kl}$. Considering the plane $z_k = z_l$ that is a subset of $\Delta_N$, one can interpret $||d_{kl}(t)||$, for $t \in [0, 1]$, as twice the distance of the path $t \mapsto \Gamma_k(t)$ (resp. $t \mapsto \Gamma_l(t)$) to the constraint $z_k = z_l$. In the spirit of interior point algorithms, it is natural to consider $-\log(||d_{kl}(t)||)$ as a barrier function for the constraint $z_k \ne z_l$, with the major difference that it gives rise to a mapping instead of a single value.

## 2.1 Basic assumptions for the path planning problem

As the path planning problem is targeted towards air traffic applications, some restrictions on the maneuvers that an aircraft can do are coming from flight dynamics and passengers comfort. First of all, velocity has to be bounded below and from above, with a quite narrow interval of efficiency dictated by engines performance. Second, curvature cannot be made arbitrary high and it is advisable to limit its total integrated value so as to minimize its detrimental effect on passengers comfort.
It is not intended in this work to consider aircraft in terminal maneuvering areas, nor the climb and descent phases: as a consequence, the path planning algorithm will not perform any change in altitude, yielding a problem that conforms to the complex configuration space modelling. Furthermore, only collision avoidance is considered, the compliance with separation norms will be addressed in a future work.

Finally, it is assumed that the planner is used in a free flight context where the aircraft are not bound to predefined routes and at a tactical level with a time horizon not exceeding 20 minutes. In such a case, one can assume that the level of uncertainty is low enough to allow a deterministic approach to be taken. The effect of wind will not be included in the model. This may be unrealistic at first glance since it is one of the most influential factor on aircraft trajectories, but given the ability of future FMS systems to infer and broadcast the wind experienced along the flight path, it seems reasonable that the wind field will be known with a sufficient degree of accuracy to adjust the initial planning.

Perfect knowledge of the aircraft positions within the airspace of interest is assumed, no communication issues are considered.

## 3 Path planning as a penalized optimization program

To turn the path planning problem into an optimization program amenable to interior point algorithms, it is needed to define first a usable criterion and second a mean of getting a tractable set of barrier functions that will be used to penalize the criterion. A major concern for the last point is the functional nature of the state space: a mean of turning the continuous time constraints into a vector of real value must be sought after.

### 3.1 Applicable criterion on a set of trajectories

For en-route traffic, overall flight cost is the dominant factor that airlines want to optimize. Since most of the time aircraft are flying near their efficient altitude and at constant velocity, it can be related to trajectory length or flight time. A second consideration that can be taken into account is the search for a flight path minimizing total curvature, as every turns induces an increase in fuel consumption and has a detrimental effect on passengers comfort.

Based on the previous remarks, a natural choice for the optimization criterion is the cumulative length. Let a smooth path $\Gamma: [0, 1] \to \mathcal{C}^N$ in the configuration space for $N$ planes be given. Its cumulative length is defined as:

$$L(\Gamma) = \sum_{j=1}^{N} \int_{0}^{1} \|\Gamma'_{j}(t)\| dt$$
In order to take into account possible priorities between aircraft, a weight may
added:

$$
L(\Gamma) = \sum_{j=1}^{N} w_j \int_{0}^{1} \| \Gamma'_{j}(t) \| dt
$$

with $w_j \in [0, 1]$ and $\sum_{j=1}^N w_j = 1$. Let $\epsilon$ be a positive value and $H: ]-\epsilon, \epsilon[ \times [0, 1] \to \mathcal{C}^N$ a smooth mapping such that for all $t \in [0, 1]$, $H(0, t) = \Gamma(t)$ and for all $s \in ]-\epsilon, \epsilon[$, $H(s, 0) = \Gamma(0)$, $H(s, 1) = \Gamma(1)$. The derivative:

$$
\frac{\partial H}{\partial s}(0, \cdot): [0, 1] \to \mathbb{C}^n
$$

may be interpreted as a tangent vector at base point Γ in the manifold of smooth
paths of Cⁿ [8]. The derivative of L applied to it can be computed from the defini-
tion:

$$
\frac{\partial L(H(s, \cdot))}{\partial s}\Big|_{s=0} = \sum_{j=1}^{N} w_j \int_{0}^{1} \frac{\left\langle \frac{\partial^2 H_j}{\partial s \partial t}(0,t), \frac{\partial H_j}{\partial t}(0,t) \right\rangle}{\left\| \frac{\partial H_j}{\partial t}(0,t) \right\|} dt
$$

using an integration by parts and, since the endpoints of H are fixed:

$$
\frac{\partial H}{\partial s}(0,0) = \frac{\partial H}{\partial s}(0,1) = 0
$$

it comes:

$$
\frac{\partial L(H(s, \cdot))}{\partial s}\Big|_{s=0} = - \sum_{j=1}^{N} w_j \int_{0}^{1} \left\langle \frac{\partial H_j}{\partial s}(0,t), \kappa_i(t) N_i(t) \right\rangle dt
$$

where $\kappa_i$ is the curvature of the path $\Gamma_i$ and $N_i$ its unit normal vector. The expression obtained is the classical first order variation for the length of a smooth path [8], summed over all components of the path $\Gamma$. In the language of manifolds,

$$
\frac{\partial L(H(s, \cdot))}{\partial s}\Big|_{s=0}
$$

gives the value obtained by applying the form $dL$ to the tangent vector

$$
\left( \Gamma, \frac{\partial H}{\partial s} (0, \cdot) \right)
$$

3.2 Turning collision avoidance constraints into semi-infinite
barrier functions

The results on complex configuration spaces presented in the first section make
natural the choice of an integral form of the barrier function [9]. For the constraint
$z_i \neq z_j$ it will be defined as:

$$ - \int_{0}^{1} \log \| \Gamma_{i}(t) - \Gamma_{j}(t) \| dt $$

Please note that such an integral does not necessary take an infinite value if $\Gamma_i(t_0) = \Gamma_j(t_0)$ and thus differs from the finite dimensional setting. However, it turns out that there is a close relationship between curvatures of the component paths $\Gamma_i$ and the norms $\|\Gamma_i(t) - \Gamma_j(t)\|$ that can be used to make the barrier functions effective.

Let $H: [-\varepsilon, \varepsilon] \times [0, 1] \to \mathbb{C}^N$ be a smooth mapping satisfying the same properties than in 3.1. The penalized criterion is:

$$ \sum_{j=1}^{N} w_j \int_{0}^{1} \| \Gamma'_{j}(t) \| dt - \sum_{j \neq k} \lambda_{jk} \int_{0}^{1} \log \| \Gamma_{j}(t) - \Gamma_{k}(t) \| dt $$

where the $\lambda_{jk}$ are strictly positive penalty weights. Let:

$$ \Theta(\Gamma)_{jk} = \int_{0}^{1} \log \| \Gamma_{j}(t) - \Gamma_{k}(t) \| dt $$

be the j, k barrier function. The first order variation formula for $\Theta(\Gamma)_{jk}$ gives:

$$ \frac{\partial \Theta_{jk}(H(s, \cdot))}{\partial s} \Big|_{s=0} = \int_{0}^{1} \frac{\left\langle \frac{\partial H_j}{\partial s}(0,t) - \frac{\partial H_k}{\partial s}(0,t), \Gamma_j(t) - \Gamma_k(t) \right\rangle}{\| \Gamma_j(t) - \Gamma_k(t) \|^2} dt $$

Reordering terms, the variation formula for the complete criterion is:

$$ - \sum_{j=1}^{N} w_j \int_{0}^{1} \left\langle \frac{\partial H_j}{\partial s}(0,t), \kappa_i(t) N_i(t) \right\rangle dt - \sum_{j=1}^{N} \sum_{k \neq j} \lambda_{jk} \int_{0}^{1} \frac{\left\langle \frac{\partial H_j}{\partial s}(0,t) - \frac{\partial H_k}{\partial s}(0,t), \Gamma_j(t) - \Gamma_k(t) \right\rangle}{\| \Gamma_j(t) - \Gamma_k(t) \|^2} dt $$

Assuming that the path $\Gamma$ is minimal with respect to penalized criterion, the first order necessary optimality conditions yields, for all $j = 1...N$:

$$ w_j \kappa_j(t) N_j(t) = \sum_{k \neq j} \lambda_{jk} \frac{\Gamma_k(t) - \Gamma_j(t))}{\| \Gamma_j(t) - \Gamma_k(t) \|^2} $$

Some important facts may be derived from the expression above. First of all, since $N_i(t)$ is of unit norm, it comes:

$$ w_j |\kappa_j(t)| \leq \sum_{k \neq j} \frac{\lambda_{jk}}{\| \Gamma_j(t) - \Gamma_k(t) \|} $$
So that there is a control on the curvature given by the penalty weights. In the special case of a path planning problem with only two mobiles, that turns out to be associated in the frame of air traffic applications with two planes encounters, the relation becomes an equality:

$$w_j|\kappa_j(t)| = \frac{\lambda_{jk}}{\|\Gamma_j(t) - \Gamma_k(t)\|}, \quad j = 1 \dots 2, k \neq j$$

In such a case, bounding the curvature of the paths will in turn bound the minimal distance separating the two mobiles: it improves over the simple collision avoidance guarantee by allowing separation norms to be enforced.

The general situation with an arbitrary number of mobiles is not so simple to deal with. However, from the expression linking $\kappa_j(t)N_j(t)$ and the sum of terms coming from the barrier functions, it is clear that a collision is not possible with two mobiles if the curvature is bounded (otherwise, the corresponding term will go to infinity while the remainder will be bounded). Only cases involving encounters with 3 or more mobiles and a special symmetry may violate this fact.

Addressing this issue can be done by noticing that the barrier functions were derived from the complex mapping $\log(z_i - z_j)$. It is in fact the complex potential generated by a simple sheet uniform distribution on the plane $z_j = z_j$. When encounters involving $p > 2$ mobiles are considered, the forbidden area in $\mathcal{C}^N$ is an intersection of $p-1$ hyperplanes and has an expression:

$$\Delta_{i_1, \dots, i_k} = \{z : z_{i_1} = z_{i_2} = \dots = z_{i_k}\}$$

where the sequence $i_1, \neq i_2, \dots, \neq i_p$ is extracted from the complete set of indices $1 \dots N$ and denotes the mobiles in interaction. The complex potential generated by a simple sheet distribution on $\Delta_{i_1, \dots, i_k}$ can be obtained using the procedure described below. First of all, the projection of a point $z$ in the configuration space onto $\Delta_{i_1, \dots, i_k}$ is given by:

$$P_{i_1, \dots, i_k}(z) = (z_1, \dots, z_{i_1-1}, h, z_{i_1+1}, \dots, z_{i_k-1}, h, z_{i_k+1}, \dots, z_n)$$

where $h = k^{-1} \sum_{j=1}^k z_{ij}$ is the mean value of the components belonging to the forbidden set. The complex potential generated by a simple sheet distribution is then:

$$\left(\sum_{j=1}^{k} z_{ij} - h\right)^{2-k} = \|z - Pz\|^{2(2-k)}$$

This potential remains harmonic, but not pluriharmonic as the one based on the log.

Taking the integral of $\|z - Pz\|^{2(2-k)}$ along a path $\Gamma$ in $\mathcal{C}^N$ yields the additional barrier function needed to deal with the encounter situation described by $\Delta_{i_1, \dots, i_k}$. It is clear from the expression of this set that $2^{N-1}$ terms will have to be taken into account in the final penalty term. It can be viewed as a consequence of the NP-Hardness of the original problem and shows how the intrinsic combinatorics
of the optimization program will appear indirectly. For practical applications, is it extremely uncommon to consider encounters involving more than 4 mobiles: on a the French airspace, it happens only a few times a year. Exceeding 5 occurs only on simulated traffic. The choice made was thus not to include the extra barrier terms and to postpone the complete formulation to a future work, where a conflict detection algorithm will be launched in a pre-processing phase in order to keep only the really needed high order barrier functions.

## 4 Interior point solver

The solver is implemented using a very simple procedure that consists in discretizing the trajectories at points regularly located in the interval [0, 1]. The flight path of aircraft i, i = 1...N is then described as a sequence $x_{i,j} = \Gamma_i(j/m)$ where m is the number of samples on each trajectory. Following the derivations made in [1], a maximum curvature $\kappa_{max}$ can be imposed using the relation:

$$ \frac{\|x_{i,j+1} - x_{i,j-1}\|}{\|x_{i,j} - x_{i,j-1}\|} \geq \sqrt{4 - \kappa_{\max}^2 \|x_{i,j} - x_{i,j-1}\|} $$

Under the near constant velocity assumption:

$$ \|x_{i,j} - x_{i,j-1}\| \approx \frac{L_i}{m} $$

with $L_i$ the total length of $\Gamma_i$, and the expression reduces to a lower bound condition on $\|x_{i,j+1} - x_{i,j-1}\|$. The integrals involved in the expression of the functional criterion turn to finite sums, in a way amenable to standard non-linear optimization algorithms. The expression of the gradient of the penalized criterion at a single vertex $x_{i,j}$ may be obtained using the following approximation of the curvature and normal vector, again assuming constant velocity and putting $l_i = L_i/m$:

$$ \bullet |\kappa_i(j/m)| = \frac{\|x_{i,j+1} - 2x_{i,j} + x_{i,j-1}\|}{l_i} $$

$$ \bullet N_i(j/m) = \frac{x_{i,j+1} - 2x_{i,j} + x_{i,j-1}}{\|x_{i,j+1} - 2x_{i,j} + x_{i,j-1}\|} $$

The gradient is then, with respect to the i-th trajectory:

$$ w_i \kappa_i(j/m) N_i(j/m) + \sum_{k \neq i} \lambda_{i,k} \frac{x_{i,j} - x_{k,j}}{\|x_{i,j} - x_{k,j}\|^2} $$

and summing up all the contribution gives rises to the complete gradient of the penalized criterion, that is used within a classical finite dimensional optimizer.
# 5 Conclusion and future work

A framework based on a functional description of the path planning problem has been described. Specially designed barrier functions inspired both by complex potential theory and semi-infinite programming allows to turn the original problem with an infinite number of constraints into a more tractable one, involving only a finite number of integrals. Bounding the curvature below allows to ensure collision avoidance. A possible implementation of the algorithms based on sampling trajectories in an even fashion was described. The complete procedure is in an early stage of development and a large amount of work still has to be performed:

* Pre-process the input traffic so as to identify encounters geometries in order to be able to use higher order barrier functions without incurring too much computational complexity.

* Investigate alternative representations of trajectories, especially those based on expansions on functional basis.

* Perform an extensive benchmark test on both synthetic and real traffic so as to assess the performance of the algorithm. A special attention must be paid to sensitivity analysis and robustness assessment against random perturbations.

* Integrate the algorithm within a realistic traffic simulator so as to quantify the effects of flight path tracking errors.

# References

1. Path planning for multiple robots: an alternative duality approach, author = Motee, N. and Jabdabaie, A. and Pappas, G., booktitle = 2010 American Control Conference, year = 2010

2. Arnol'd, V.: The cohomology ring of the colored braid group. Mathematical notes of the Academy of Sciences of the USSR 5(2), 138-140 (1969). DOI 10.1007/BF01098313. URL http://dx.doi.org/10.1007/BF01098313

3. Canny, J., Reif, J.: New lower bound techniques for robot motion planning problems. In: Proceedings of the 28th Annual Symposium on Foundations of Computer Science, SFCS '87, pp. 49-60. IEEE Computer Society, Washington, DC, USA (1987). DOI 10.1109/SFCS.1987.42. URL http://dx.doi.org/10.1109/SFCS.1987.42

4. Delahaye, D., Durand, N., Alliot, J.M., Schoenauer, M.: Genetic algorithms for air traffic control system. In: 14th IFORS Triennial Conference (1996)

5. Durand, N., Alliot, J.M.: Ant colony optimization for air traffic conflict resolution. In: 8th USA/Europe Air Traffic Management Research and Development Seminar (2009)

6. Latombe, J.C.: Robot Motion Planning: Edition en anglais. The Springer International Series in Engineering and Computer Science. Springer (1991). URL http://books.google.com.au/books?id=Mbo_p4-46-cC

7. LaValle, S.: Planning Algorithms. Cambridge University Press (2006). URL http://books.google.com.au/books?id=-PwLBAAAQBAJ

8. Michor, P., Mumford, D.B.: Riemannian geometries on spaces of plane curves. Journal of the European Mathematical Society 8(1), 1-48 (2006). DOI 10.4171/JEMS/37. URL http://dx.doi.org/10.4171/JEMS/37

9. Shapiro, A.: Semi-infinite programming, duality, discretization and optimality conditions. Journal of Optimization 58(2), 133-161 (2009). DOI 10.1080/023319309027300070