A Discrete-Continuous Algorithm for Free Flight Planning

We propose a hybrid discrete-continuous algorithm for ﬂight planning in free ﬂight airspaces. In a ﬁrst step, our DisCOptER method (discrete-continuous optimization for enhanced resolution) computes a globally optimal approximate ﬂight path on a discretization of the problem using the A ∗ method. This route initializes a Newton method that converges rapidly to the smooth optimum in a second step. The correctness, accuracy, and complexity of the method are goverened by the choice of the crossover point that determines the coarseness of the discretization. We analyze the optimal choice of the crossover point and demonstrate the asymtotic superority of DisCOptER over a purely discrete approach.


Introduction
optimum on a coarse, approximate, artificial network, and to use the resulting approximate solution to initialize a Newton method for the solution of a continuous optimal control problem. The correctness and effectiveness of this strategy depends on the choice of the crossover point between the discrete and the continuous part of the algorithm. The network for the discrete part must be coarse enough to be searched efficiently, and fine enough to guarantee sufficient proximity to the continuous optimum, such that a subsequent Newton iteration will converge to the latter. Clearly, the convergence radius of the continuous method strongly depends on the gradients of the wind field and affects, together with the approximation error associated with the graph, the computational complexity of both methods. We shall show that this idea is ideally suited for free flight settings.
Our aim in this paper is to demonstrate the potential of combining discrete and continuous optimization methods for the solution of problems that involve space or time discretizations. The goal is to achieve global optimality and rapid convergence at the same time. The model that we consider is simple, and the special characteristics of flight planning have left their imprint. The basic idea, however, should, with suitable modifications, be applicable to various problems of similar nature. In this vein, our paper intends to give a first indication of the usefulness of such approaches.
The paper is structured as follows. Sections 2.1-2.3 describe the free flight problem that we consider. The DisCOptER algorithm is introduced in Section 2.4 including error and complexity estimates. Finally, Section 3 provides a computational study and analyzes the choice of the switch over point.

Free flight planning
We consider in this paper an idealized version of the flight planning problem in 2D Euclidean space subject to a stationary wind field. We want to compute a flight path x : [0, T] → R 2 , τ → x(τ) that connects an origin and a destination x O , x D ∈ R 2 ; the path parameter τ measures the flight time. The path is influenced by a smooth field of stationary wind w : R 2 → R 2 of bounded magnitude w ≤w. Flying at an airspeed v : [0, T] → R 2 , τ → v(τ) of constant magnitude v =v >w, the aircraft arrives at time T ∈ R ≥0 , which we seek to minimize. Our setting is chosen for ease of exposition, but our method carries over to more complex 3D and/or time dependent versions, or other objectives, in particular, minimization of fuel consumption.

Continuous approach: optimal control
In free flight, the flight path is not restricted to a predefined airway network of waypoints and segments. Instead, any Lipschitz-continuous path x : [0, T] → R 2 , with x t − w =v almost everywhere, connecting origin x O and destination x D , is a valid trajectory. Among those, we shall find one of minimal flight duration T. This classic of optimal control is known as Zermelo's navigation problem [23].
In order to formulate the problem over a fixed interval [0, 1] independent of the actual flight duration, we scale time by T −1 as usual in free end time problems and arrive at the following optimal control problem for the flight duration T ∈ R, the flight path x ∈ H 1 ([0, 1]) 2 , and the airspeed v ∈ L 2 ([0, 1]) 2 : Here, the constraint c : Z → Λ maps from the primal domain Z :

Optimality conditions
Let us briefly recall the necessary and sufficient optimality conditions for the optimal control problem (1). With z = (T, x, v) ∈ Z and λ ∈ Λ, the Lagrangian is defined as If z is a (local) minimizer, the necessary first order or Karush-Kuhn-Tucker (KKT) condition and the necessary second order condition (3) are satisfied at z for some λ and the sufficient second order condition (SSC) holds for some α > 0, then z is a locally unique solution of (1) and stable under perturbations of the problem, e.g., due to sufficiently fine discretization. Let us point out that for wind fields with non-vanishing second derivative, in general there is no closed form solution of the necessary conditions (3), such that solutions must be approximated numerically. Approaches to computing solutions of (1) generally fall into two classes: indirect methods relying on Pontryagin's maximum principle [25,26] and direct methods based on discretization of the minimization problem (1) [27,28]. Indirect methods lead to a boundary value problem for state x and adjoint state λ together with a pointwise optimality condition for the control v, which can be solved by shooting type methods, collocation, or spectral discretization approaches [29,30]. The discretization of direct methods, usually by collocation or spectral methods, translates the optimal control problem in a finite dimensional nonlinear program to be solved by corresponding optimization problems [16,17,31]. While indirect methods using multiple shooting lead to smaller problem sizes than direct methods based on collocation in particular for high accuracy requirements, the latter are widely seen as being easier to implement and use, in particular in the presence of state constraints.
In both approaches, Newton-type methods for solving either discretized boundary value problems or nonlinear programming problems converge in general only locally towards a closeby minimizer. Thus, sufficiently good initial iterates need to be provided. The domain of convergence can be extended using line search methods or trust region methods, but without guarantee of global optimality. Special care has to be taken in case of non-convex problems, since the Newton direction need not be a descent direction if far away from a minimizer satisfying the sufficient second order conditions. Convexification, truncation of iterative solvers [32][33][34], or solvers for non-convex quadratic programs can be used to address this.

Discretization error
The discretization error introduced by the midpoint rule is well-known to be of second order. For given where δτ = max i=0,...,n−1 τ i+1 − τ i is the mesh width, see, e.g., [35]. This translates into a corresponding error in the objective, i.e., the flight time T. Different goal-oriented error concepts have been investigated [36,37] for equality constrained optimal control problems. The excess in actual flight time when the computed path is implemented depends quadratically on the path deviation x h − x L ∞ ([0,1]) , and is therefore of order O(δτ 4 ). In any case, the error depends directly on the mesh width δτ, which we therefore use as a coarse but simple qualitative measure of accuracy.

Newton-KKT solver
Analogous to the continuous Lagrangian (2), we may formulate its discretized counterpart and the corresponding necessary first order optimality condition Writing χ = [z h , λ h ] T , this can be solved using Newton's method by computing For smooth wind fields, there is some ω < ∞ related to an affine invariant Lipschitz constant of L h , such that holds [38]. Thus, Newton's method converges quadratically if started sufficiently close to a locally unique solution point χ * , i.e., if χ 0 − χ * < ω −1 . This convergence radius is in general mesh-independent [38,39], and does not depend on the final accuracy in terms of mesh width δτ, but only on problem parameters such as derivatives of the wind w.

Time complexity
The run time of the Newton-KKT solver is determined by the number of steps and the cost of each step. The computational effort of a Newton step is dominated by the cost of solving the linear equation (8). Due to the ODE structure, L h (χ) is an arrow-shaped matrix with band width independent of δτ. Assuming quasi-uniform meshes, i.e., there is a generic constant C > 0 such that τ i+1 − τ i ≥ Cδτ, this structure allows for an efficient solution in O(δτ −1 ) time using direct band solvers.
Starting sufficiently close to the solution, say ω χ 0 − χ * < 1, allows to bound the truncation error by linear convergence as Thus, the overall complexity in terms of mesh width δτ is Remark 1. If an inexact Newton method based on a geometrically refined sequence of meshes with mesh width δτ k = β k l for some l δτ and β < 1 is used, the number of iterations is determined by the linear convergence of the collocation discretization, while the truncation error is quickly diminished by the quadratic convergence of Newton's method [38]. Since the effort per Newton step grows geometrically to its final value, the complexity reduces to provided the initial error is sufficiently small, i.e. ω χ 0 − χ * ≤ β 2 .

Discrete approach: shortest paths in airway networks
If flight paths are restricted to a predefined airway network of waypoints and segments, flight planning becomes a special kind of shortest path problem on a digraph. Let V ⊂ R 2 be a finite set of waypoints including x O and x D , and We denote the set of all flight paths by P. The total flight duration T(p) for a path p = (x 0 , . . . , x n ) ∈ P is given in terms of the flight duration T(e i ) for an arc e i = (x i , x i+1 ) by The travel time on the arc e i can be computed from the local ground speed and the constant airspeed v =v by solving the quadratic equation This yields The discrete optimization problem to be solved is now min p∈P T(p).

Graph construction.
Discrete approaches to (flight) path planning fall into two classes. The first class are samplingbased algorithms. They construct the search graph by some kind of sampling during the execution of the shortest path algorithm, which is usually an A * -method. This class includes rapidly exploring random trees (RRT) and graph (RRG) algorithms, that are often used in robot path planning [8]. There are versions that guarantee convergence to a global optimum with probability one [40], however, although undoubtedly often extremely efficient in practice, these methods do not provide a priori error bounds or complexity estimates.
The second class are graph-based algorithms, that take the search graph as an input. Impressive super-fast performance in practice, and also theoretically for special classes of graphs, has been achieved by making use of preprocessing as well as sophisticated data structures. However, "proving better running time bounds than those of Dijkstra's algorithm is unlikely for general graphs" [6]. In many applications, including traditional flight plannning on airway networks, the search graph is canonical. In applications involving space and time discretization, like in free flight, graph construction is a degree of freedom. Using an appropriate discretization, a priori bounds on the runtime and the accuracy of the solutions can be derived. This is the approach that we take.
We cover free flight zones with "locally densely" connected digraphs, i.e., digraphs with a certain density of vertices and arcs.
if it satisfies the following conditions: We call h the vertex density and 2h + l the connectivity radius of an (h, l)-dense graph.
With this definition, |V| ∈ O(h −2 ) and |E| ∈ O((2h + l) 2 h −4 ) hold. Furthermore, it is easy to see that this graph structure implies that G is strongly connected and therefore P is nonempty, i.e. a path from x O to x D exists.

Discretization error
Similar to the collocation error (6), the error due to graph discretization depends on the vertex density h and the connectivity radius 2h + l. We restrict our exposition in this paper to a plausible argument for the error order in terms of h and l and refer the reader to [41] for rigorous error bounds of the same orders.
Spatial deviation due to vertex spacing is bounded by 1 = O(h), while the linear interpolation error depending on the curvature of the continuous optimal path is bounded by Fig. 1a. Together, the deviation measured as pointwise normal distance is bounded by O(l 2 + h). As mentioned in Section 2.2 above, this translates quadratically into a flight time error of order O(l 4 + h 2 ).
Assuming that the path actually contains arcs of length around l + 2h, this order of error is not only an upper bound, but also the expected error. The existence of arcs of length 2h + l in the optimal discrete path is likely, as the following argument based on Fig. 1b shows. With quasi-uniform vertex spacing, the number of adjacent vertices within a distance 2h +l ≤ 2h + l is of order O((2h +l) 2 /h 2 ), such that the average angle between arcs of length up to 2h +l, and thus the expected angular deviation α between the discrete path and the optimal continuous path is of order O((2 +l/h) −2 ). The geometric length difference between the discrete path and the continuous path is of . The expected length and therefore travel time error induced by angular deviation is smallest ifl is largest, i.e., arcs of maximum length 2h + l are to be preferred.
This analysis also implies that the total flight time error is of order O(l 4 + h 2 + (h/l) 4 ) and that h = o(l 2 ) ensures convergence for l → 0.

A * shortest path algorithm
The state-of-the-art for finding shortest paths in airway networks is the A * algorithm [10], which extends Dijkstra's algorithm by a (heuristic) lower bound for the distance of some vertex x i to the destination x D in order to prioritize the search. Depending on the tightness of the heuristic bound, A * can discard a substantial part of the graph and reduce the run time considerably. As a particularly fast and simple heuristic we employ the travel time along a straight line between x i and x D , assuming maximum tail wind, i.e., Tighter and more complex heuristics exploiting local bounds on the wind have been proposed for flight planning, also for time-dependent wind fields [10]. The travel time for each arc is again calculated via numerical integration of (12) using the same method as for the collocation with a fixed discretization. In order to draw up a fair comparison, we process arcs on the fly and thus eliminate any major preprocessing. This also makes the approach extendible to the time dependent case.

Time complexity
Using a Fibonacci heap for the priority queue, A * can be implemented to run in time where |A| and |V| are the numbers of arcs and vertices, respectively. As outlined above, we may assume l h. If we -conservatively -assume l 2 ≥ h 2 log h −1 , the total complexity becomes

Graph structure
Given the above considerations, we now address the question of which graph structure in terms of h and l is the most efficient, i.e., which structure achieves a minimal error for a given computational budget b or, equivalently, a desired accuracy with minimal effort. Models for computational work and accuracy are The constraint W(h, l) = b yields l 2 = h 4 b, which, when inserted into the objective, leads to the unconstrained problem The necessary optimality condition 4BH 3 + 1 − 2B −1 H −3 = 0 is a quadratic equation in H 3 with solution Inserting this into the constraint yields l 2 = h −1+ √ 33 8 . We conclude that graphs of optimal efficiency should follow the law such that the computational complexity of the discrete optimization can be expressed as with an associated error of O(l 4 ) in flight time and of O(l 2 ) in path approximation.

DisCOptER algorithm
Due to their superior angular resolution, arcs of length l will occur in the optimal discrete path, while the length of flight path segments in the collocation solution is around δτTv. Hence, we expect the accuracy of continuous and discrete optimization approaches to be comparable if l = δτTv. Obviously, for increasing accuracy l → 0, the collocation effort of order O(δτ −1 log δτ/ log(ω χ 0 − χ * )) grows much slower than the A * effort of complexity O(l −6 ). On the other hand, the collocation approach converges only locally.
We therefore propose a hybrid algorithm that combines the strengths of discrete optimization and optimal control, see Algorithm 1. First, a discrete shortest path of low accuracy is calculated using the A * algorithm as described in section 2.3. This is used as initial iterate for solving the KKT system (5) to the desired final accuracy using the ordinary Newton method. Employing linesearch or trust region globalization would enlarge the convergence domain of Newton's method and allow for coarser graphs to be used, but at the cost of less robust and efficient convergence. In favour of robustness and simplicity, we restrict the attention to the ordinary Newton method.
Algorithms that combine methods from discrete and continuous optimization have been proposed before. A reverse deterministic combination going from continuous to discrete has been proposed in [42] based on dynamic programming principles. Other combinations involve a stochastic discrete stage, like rapidly-exploring random trees (RRT, RRT*) (see, e.g., [43]) or Probabilistic Roadmaps (PRM) [44] in combination with a second NLP stage. Another approach was taken in [45], where the authors use a combination of A * and RRT* to find an optimal trajectory. These algorithms reveal remarkable performance when it comes to obstacle avoidance. Since our goal is, however, to develop an algorithm with a priori error estimates and bounded runtime, we do not make use of any stochastic approaches in the DisCOptER algorithm. (1)

Complexity
The runtime of this hybrid algorithm is comprised of the runtime of the A * algorithm (13) and the runtime (10) of the KKT-Newton solver for the collocation system. Provided that the initial iterate based on the discrete solution p is sufficiently close to the continuous optimum, i.e., ω χ 0 − χ * < 1, Newton's method converges locally as described by (10). For the path approximation error, we expect χ 0 − χ * ≈ Cl 2 , which leads to an overall complexity of or R H = O(l −6 + δτ −1 ) subject to ωCl 2 < 1 if an inexact Newton method is used. Both complexity bounds essentially suggest to choose a graph as sparse as possible (l → ∞), only restricted from above by the accuracy necessary for the Newton-KKT solver to converge locally, i.e. ωCl 2 < 1, which is independent of the final accuracy δτ.

Results
We validate in this section the effectiveness and the efficiency of our algorithm on a test set of four problems of increasing difficulty. Our aim is to demonstrate that our hybrid approach is asymptotically superior to a purely graph based alternative. We discuss the convergence properties of the method in relation to the choice of the crossover point and its dependence on the minimum required graph density. Based on these results, we evaluate the DisCOptER algorithm for varying graph densities, using the theoretically optimal graph structure of h = l 2 , cf. (14). We will see that -as expected -the best performance is achieved on very sparse graphs. The chapter concludes with a computational comparison of the hybrid algorithm with a purely discrete approach.

Test problems
We test our algorithm on a set of simple, but representative examples that are well suited to demonstrate our method, and not far from real world situations or arguably even more difficult. The wind fields are constructed in a such a way that the straight line from the origin to the destination is particularly unfavorable. The wind speed is limited tow = 1 2v , which is rather strong, but not unrealistic. Even though not formally required, the graph nodes are positioned on a uniform Cartesian grid, such that the diagonal distance of two adjacent nodes is 2h. For the sake of simplicity the graph structure will be described based on the node spacing in x-direction h x = √ 2h. This directly defines l = √ h by (14).
In the first test problem a), the wind field is a laminar flow of opposing parallel currents, namely, with H = 0.5, see Figure 2; a similar problem is discussed in [46]. Due to its simplicity, this problem has only one distinct minimum, a property that we make use of in the next section. In problems b)-d), the wind w is the sum of an increasing number of vortices w i , each of which is described by where s is the spin of the vortex (s=+1: counter-clockwise, s=−1: clockwise), r = x − z 2 is the distance from the vortex center z i , α is the angle with respect to the center and the x-axis with tan(α) = (x−z i ) 2 (x−z i ) 1 and the absolute vortex wind speedw is a function of r and the vortex radius R: Problem b) involves one large vortex with R = 1/2 at z = [0.5, −0.1] T , see Figure 3. This causes Newton's method to converge to a suboptimal trajectory (above the vortex) if initialized with the straight trajectory. That is the case if the DisCOptER algorithm is used with a minimally sparse graph with h x = 1 (see second subfigure). We observed, that the discrete shortest path passes below the vortex for any h x > 1. From there the globally optimal solution shown in the third subfigure was found for h x = 1/6. Problem c) involves 15 vortices with R = 1/8. One is centered at z = [0.5, −R/2] T , the others are regularly aligned as seen in Figure 4. Vortices with positive spin (clockwise) are colored green, vortices with negative spin (anti-clockwise) are colored red. Due to the turbulence of this wind field, a plain application of Newton's method is not guaranteed to converge. As an example, Subfigure 2 shows again the result with the trivial initialization. Further there are several local minima, Subfigures 3-4 show two of them. A relatively high graph density (h x ≤ 1/17) is required to find the globally optimal trajectory (see Subfigure 5). Blue dots: graph with h = l 2 and 1/h x = 1, 5, 6, and 16, respectively, Red: discrete optimal trajectory, Green: continuous optimal trajectory. Note that the straight trajectory is particularly unfavorable.
Problem d) involves 50 regularly aligned vortices of radius R = 1/16 ( Figure 5). This is clearly an exaggeration, and no commercial plane would ever try to traverse a wind field like this. We use the instance to show that the proposed algorithm outperforms existing methods even under the most adverse conditions. In fact, the high level of non-convexity exacerbates the situation regarding the convergence of Newton's method. Note that, the magnitude of derivatives of (18) is coupled directly to the vortex radius. With R = 1/16 the vortices here are half as large as in the previous example. In turn, the first and second order derivatives are 2 and 4 times larger, respectively. Consequently, a quite dense graph with h x ≤ 1/60 is required to make Newton's method converge reliably. Note that this wind field is point-symmetric with respect to [0.5, 0] T , which allows for two equivalent global optima (see Subfigures 2 and 3). It is a priori not obvious which one will be found. This depends on the discrete optimum.

Computational complexity
Before discussing the DisCOptER algorithm we validate that the optimal control methods are asymptotically more efficient than a purely graph-based approach. This can be investigated only for the first test problem. Due to the simplicity of this wind field, Newton's method will converge from a trivial initialization, i.e., the straight line from x O to x D . On the other hand, discrete flight paths were calculated with varying accuracy, that is, with varying l, which is, according to Section 2.3.2, a measure for the accuracy of the calculated path. For the sake of consistency, we indicate the accuracy of the continuous solution by l C := δτL, where L is the path length. Figure 6 clearly confirms the estimated time complexities of O(l −6 ) for the discrete approach (see eq. 13) and O(δτ −1 ) for the continuous approach (see eq. 10).  Figure 7 shows the runtime of the DisCOptER algorithm for various graph parameters (h, l). Some key observations can be made for all four test problems. Towards the top right of each figure the graph density increases, which comes with an increased runtime dictated by the graph searching part of the algorithm. The black dashed line represents the theoretically optimal graph structure h = l 2 , cf. (14); it is shown for orientation as the results in the following sections are computed with this setting.

Minimum graph requirements
From left to right the wind fields become increasingly more non-convex. As these plots reveal, this comes with generally higher runtimes, but more importantly with a region of low graph densities where the algorithm converges either to a local minimum or not at all (gray areas). As discussed before, this is not an issue in case a) since this problem is convex and Newton's method converges even from a trivial initialization. Even in case b) we found that a graph with h x < 1 is good enough such that we find the optimal flight path (only one additional column of nodes between start and destination is required). This outcome is not visible here due to the limited resolution of the plot. The convergence problem becomes all the more apparent with instances c) and d). In both cases a good number of local minima exists, each of which has a rather small radius of convergence such that Newton's method fails if not initialized with sufficient accuracy. Especially in case d) we see a patchwork of runs that converge presumably by chance in an unpredictable way. Some of this might be compensated by globalization techniques and an explicit treatment of non-convexity, which, however, would affect the efficiency of Newton's method. In order to provide reliable results, we need to use a graph density of at least l < 0.15 and 0.11 for case c) and d), respectively.
Interestingly, the node distance h appears to have a much stronger effect on the convergence than the connectedness width l. This might be explained to some extent in terms of the distance and angular error (cf. section 2.3.2). Low l decreases the angular resolution and thus induces an increased angular error. The discrete flight path then tends to zig-zag along the optimum. This can easily be smoothed by Newton's method. However, if the discrete flight path is parallely offset from the optimum (distance error due to large node distance h), it might quickly leave the convergence radius of Newton's method.

Optimal crossover point
In Section 2.3.2 we derived h = l 2 as the optimal graph structure, cf. (14). Using this setting and sampling over various graph densities for the test problems leads to the results presented in Figure 8. Obviously, an increased graph density comes with a computationally more expensive graph searching task (pink, top, cf. (15)). In turn the NLP part of the algorithm (green, bottom) becomes cheaper following (10), since the initial guess gets closer to the optimum. From test problem a) it can be seen that the best performance is achieved -independently of the overall accuracy δτ -with a relatively low graph density of 1/l ≈ 5, where the graph search is still negligibly cheap but the graph is already dense enough for Newton's method to converge with only one iteration. The exact numbers depend strongly on implementation details. We do not claim to have an optimal realization of either part of the algorithm. The trend towards low graph densities, however, is confirmed by the following examples.
As it turns out, the lower bound on the required graph density imposed by the non-convexity of the wind field is the more important criterion. Looking at examples c) and d), where we excluded the area that we identified as not trustworthy in the previous section, we conclude that we want the graph to be as sparse as possible and only as dense as necessary.

Computational complexity
We finally show that globally optimal shortest paths can be calculated more efficiently using the proposed DisCOptER algorithm than with a purely graph based approach, see Figure 9. In the previous section we showed that the algorithm performs best if the graph is chosen rather sparse while respecting the problem-specific minimum density. Consequently, the calculation of the discrete solution is relatively cheap and we can benefit from the computational efficiency of Newton's algorithm. We can also confirm that the time complexity of the DisCOptER Algorithm is O(δτ −1 ), see eq. 17, and that the time complexity of the purely graph-based approach is O(l −6 ), see (13).

Conclusion
In this paper we presented the novel DisCOptER algorithm to calculate flight paths in free flight airspaces utilizing a combination of discrete and continuous optimization. We demonstrated that the achieved efficiency is asymptotically much better than the conventional purely discrete alternative. Even though the algorithm was described for the static two-dimensional case it is strongly promising also for more complex cases, to which it can directly be transferred.
Our study also reveals a need for more theoretical analysis of the problem. In order to design a one-shot algorithm with theoretical efficiency guarantees, a priori error estimates allowing the determination of a minimum required crossover graph density is needed. This will of course depend mainly on the characteristics of the wind field including first and second order derivatives. On the other hand, a posteriori error estimates and adaptive coarse-to-fine graph refinement strategies will be necessary for robustness and efficiency in practice. We investigate some of these questions in [41].