Functional Interior Point Programming Applied to the Aircraft Path Planning Problem

Multiple aircraft trajectory planning is a central problem in future air traffic 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 models 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.


Introduction
The path planning problem for multiple robots evolving in a possibly dynamic environment 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 algorithms to prove convergence in a probabilistic sense. However, due to the nature 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. Successful 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 optimization 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 optima, 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 optimization program associated to the path planning problem, and a relation with a recently introduced complex harmonic navigation function will be pointed out. Finally, a possible algorithmic implementation of the solver yielding an approximate solution to the original problem will be detailed.

Complex configuration spaces
Given a set of N mobiles with planar motion, the complex configuration space C N is defined to be the Cartesian product of n copies of the complex plane C with the set ∆ N = (z 1 , . . . , z N ) , ∃i = j, z i = z j of simultaneous positions removed. Any path connecting two points in C N is an admissible collision-free planning of N trajectories. It is easy to show that C N is a path connected space using a sequential planning argument: let z 0 1 , . . . , z 0 N be any initial configuration and let z 1 1 , . . . , z 1 N be the desired final situation. Let z 1 , . . . z N be fixed. Then it exists a path joining z 0 1 and z 1 1 in C−{z 1 , . . . 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 C n and assumes that only one mobile moves at a time.
Geometrical insights about C N can be gained from the knowledge of its cohomology 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: Please note that in the original paper [2] a factor (i2π) −1 appears in the expression: it is removed here as the real part of ω k,l is of primary interest. Let a C 1 path Γ : [0, 1] → C N be given. Then: Rewriting the right hand side, it comes: Given any two complex numbers z 1 = x 1 + iy 1 , z 2 = x 2 + iy 2 , the real part of the product 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: thus admits an interpretation as an expansion rate (resp. rotation rate) for the vector d kl (t) = Γ k (t) − Γ l (t) and in turn the integral: Γ dω k,l will give log ( d kl (1) ) − log ( d kl (0) ) as its real part and 2πθ for its imaginary part, with θ the winding number of the path d kl . Considering the plane z k = z l that is a subset of ∆ N , one can interpret d kl (t) , for t ∈ [0, 1], as twice the distance of the path t → Γ k (t) (resp. t → Γ 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 = z l , with the major difference that it gives rise to a mapping instead of a single value.

Basic assumptions for the path planning problem
As the path planning problem is targeted towards air traffic applications, some restrictions on the manoeuvers 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 manoeuvering 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.

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.

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 the optimization criterion is the cumulative length. Let a smooth path Γ : [0, 1] → C N in the configuration space for N planes be given. Its cumulative length is defined as: In order to take into account possible priorities between aircraft, a weight may added: using an integration by parts and, since the endpoints of H are fixed: where κ i is the curvature of the path Γ 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 Γ . In the language of manifolds, gives the value obtained by applying the form dL to the tangent vector

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 = z j it will be defined as: Please note that such an integral does not necessary take an infinite value if Γ i (t 0 ) = Γ 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 Γ i and the norms Γ i (t) − Γ j (t) that can be used to make the barrier functions effective. Let H : ]−ε, ε[×[0, 1] → C N be a smooth mapping satisfying the same properties than in 3.1. The penalized criterion is: where the λ jk are strictly positive penalty weights. Let: be the j, k barrier function. The first order variation formula for Θ (Γ ) jk gives: Reordering terms, the variation formula for the complete criterion is: Assuming that the path Γ is minimal with respect to penalized criterion, the first order necessary optimality conditions yields, for all j = 1 . . . N: Some important facts may be derived from the expression above. First of all, since N i (t) is of unit norm, it comes: 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: 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 κ 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 C N is an intersection of p − 1 hyperplanes and has an expression: where the sequence i 1 , = i 2 . . . , = i p is extracted from the complete set of indices 1 . . . N and denotes the mobiles in interaction. The complex potential generated by a simple sheet distribution on ∆ i 1 ,...,i k can be obtained using the procedure described below. First of all, the projection of a point z in the configuration space onto ∆ i 1 ,...,i k is given by: P i 1 ,...,i k (z) = z 1 , .., z i 1 −1 , h, z i 1 +1 , .., z i k −1 , h, z i k +1 , .., z n where h = k −1 ∑ k j=1 z i j is the mean value of the components belonging to the forbidden set. The complex potential generated by a simple sheet distribution is then: 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 Γ in C N yields the additional barrier function needed to deal with the encounter situation described by ∆ i 1 ,...,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.

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 = Γ i ( j/m) where m is the number of samples on each trajectory. Following the derivations made in [1], a maximum curvature κ max can be imposed using the relation: Under the near constant velocity assumption: with L i the total length of Γ 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: The gradient is then, with respect to the i-th trajectory: 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.

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 trafic simulator so as to quantify the effects of flight path tracking errors.