An integrated rolling horizon and adaptive-refinement approach for disjoint trajectories optimization

Planning for multiple commodities simultaneously is a challenging task arising in divers applications, including robot motion or various forms of traffic management. Separation constraints between commodities frequently have to be considered to ensure safe trajectories, i.e., paths over time. Discrete decisions to ensure at least one of often multiple possible separation conditions renders planning of best possible continuous trajectories even more complex. Hence, the resulting disjoint trajectories optimization problems are mostly solved sequentially or with restricted planning space, potentially leading to losses in the usage of sparse resources and system capacities. To tackle these drawbacks, we develop a graph-based model for disjoint trajectories optimization with general separation requirements. We present a novel technique to derive a discretization for the full available space of motion. This can depict arbitrary, potentially non-convex, restricted areas. This necessitates solving an integer linear optimization program whose size scales with the number of discretization points. Thus, even for moderately sized instances a sufficiently detailed representation of space and time leads to models too large for state of the art hard- and software. To overcome this issue, we develop an adaptive-refinement algorithm: Starting from an optimal solution to the integer program in a coarse discretization, the algorithm re-optimizes trajectories in an adaptively-refined discretized neighborhood of the current solution. This is further integrated into a rolling horizon approach. We apply our approach to the integrated trajectory optimization and runway scheduling in the surrounding of airports. Computational experiments with realistic instances demonstrate the efficiency of the method.


Introduction
Safely separated and overall optimal trajectories for multiple commodities are crucial in many applications, including air traffic management, packet routing, or multi robot motion. Thereby, the challenging task of determining optimal trajectories for the individual commodities is paired with having to maintain some kind of separation, which could be temporal or spatial. Often only one of several possible separation conditions has to be fulfilled, which adds discrete decisions to the continuous trajectory planning problem. In real world settings, such as arrival planning at airports, the inherent complexity of such discrete-continuous optimization problems leads to highly restrictive and heuristic procedures to maintain a manageable workload for authorities. However, limiting the considered planning space to highly restricted networks can reduce the system's theoretical capacity and increase costs.
In this work we present an approach for the challenging task of disjoint trajectory optimization. This general modeling approach we develop for the disjoint trajectories problem is capable of depicting arbitrary separation requirements for commodities and uses a graph-based representation for space and time. Hereby, trajectories are paths in a corresponding time-expanded network. As the resulting problem is NP-hard even in very restricted settings (Hoch et al. (2020)), we derive an integer linear program. Using (mixed) integer linear programming (MILP) techniques, this can be solved to global optimality -with respect to the selected representation.
To overcome the use of highly restricted, predefined networks, we further develop a discretization technique that enables trajectory optimization in the full available space of motion. Such an approach is particularly beneficial in, for example, approach planning to airports, where broadly one-dimensional arrival routes are currently used. Considering the whole three-dimensional airspace in a free-flight spirit as introduced by Force (1995) offers the prospect of not only expanding airspace capacity, but also improving efficiency of runway usage. An additional application is traffic planning of ships in canals, as considered in Lübbecke et al. (2019). Here, capacities could be further increased by allowing speed and route adjustments using the full geometry of the canal rather than following predefined tracks.
Our approach to discretizing three-dimensional space and time allows planning in an arbitrary, potentially non-convex, environment. It comes at no surprise that -depending on the desired level of detail -the presented model quickly exceeds the capabilities of current hardware. Therefore, we use advanced algorithm engineering techniques to keep instance sizes manageable: We develop an adaptiverefinement approach that starts with a coarse discretization of space and time and computes globally optimal trajectories with respect to this discretization. Later, in an adapted search space around said trajectories the discretization is refined and re-optimized considering those search spaces. To allow for even larger instances, we further embed this adaptive-refinement in a rolling-horizon approach. Here, we follow the trajectories over time and consider refinement in segments.
We concretize this general approach for an application in air-traffic management, namely approach planning at airports. Approach planning requires runway 1 3 An integrated rolling horizon and adaptive-refinement approach… scheduling and aircraft trajectories to be optimized simultaneously -a task well within our model's capabilities as it can depict all necessary en-route distances and ensure minimal temporal distances between aircraft. For optimization we use an objective function to fairly mediate between fuel consumption and deviation from the scheduled time.
In the computational study conducted for this application, we consider realistic instances with up to ten aircraft. Results illustrate not only the algorithmic properties, but also the efficiency, of this approach.
The remainder of the paper is organized as follows: A review of related literature is followed by the development of the general disjoint trajectories model in Sect. 2. In Sect. 3, a meta-algorithm for iterative reoptimization of trajectories in smaller sub-problems is presented. The discretization method for three-dimensional space and time is given in Sect. 4, while in Sect. 5 the adaptive-refinement and the rolling horizon approach are introduced. Section 6 presents the application of approach planning at airports and the corresponding computational experiments are detailed in Sect. 7. We close with the conclusion and outlook in Sect. 8.
Related literature: The applications for multi-commodity trajectory optimization are vast and so too is the existing literature. It is therefore not possible to provide an exhaustive survey, but we will highlight the studies most relevant to the general problem and those with applications in air-traffic management. Where appropriate, we provide pointers towards related fields.
In a space-time discretized setting, the disjoint trajectories planning problem becomes a variant of the disjoint paths problem as introduced by Karp (1975). This problem has been widely studied, with an overview of complexity results in Fortune et al. (1980); Kobayashi and Sommer (2010). Algorithmic approaches to deal with the classic versions can be found in Yu and LaValle (2016); van Den Berg et al. (2009); Sharon et al. (2015), to name but a few. While all approaches are promising for the described settings, they are not easily transferable to large scale timeexpanded networks. Such networks are necessary to properly depict free movement in the considered space. Time-expanded networks allow to interpret safety separations as arc dependencies. In the context of network design, arc dependencies are studied by Oellrich (2008). For scheduling applications the alternative graph model by Pacciarelli (2002) is often used to impose proper separation and precedence.
Optimization for single trajectories has been treated both in discrete as well as continuous settings. Continuous methods use the equations of motion and apply optimal control approaches to solve the resulting problems. Although globalization methods do exist, these methods can generally only provide locally optimal solutions. von Stryk and Bulirsch (1992); Pecsvaradi (1972); Hagelauer and Mora-Camino (1998), or more recently Khardi (2012), provide exhaustive overviews on those methods.
In a discrete setting, global optimal solutions can be obtained, despite limitations in describing the commodities dynamics of motion. In Blanco et al. (2016), so called super optimal winds are used in an A * algorithm for trajectory planning for one aircraft. A graph-based approach is also used in Blanco et al. (2017) to minimize the crossing cost of a trajectory. Fischer and Helmberg (2014) propose an exact method to compute shortest paths in a time-expanded network with preferences for early arrivals, where the graph is built dynamically.
When it comes to safe planning of multiple trajectories, numerous conflict detection and resolution approaches aim to de-conflict given trajectories: Kuenz (2015) describes an algorithm for large scale conflict detection, and trial and error conflict resolution based on sub-dividing the search space. In Dias et al. (2020), a two stage algorithm for conflict resolution with recovery to the given initial trajectories is proposed. A common feature of most of the conflict resolution approaches is consideration of a pre-defined set of actions to avoid the conflict. These usually include a set of possible heading or speed changes. Pelegrín et al. (2021) introduce fairness in deconfliction for given trajectories in a track network. For general reviews on conflict detection and resolution approaches, especially in aviation, see Kuchar and Yang (2000); Martín-Campo (2010); Ribeiro et al. (2020).
Direct optimization for safe trajectories of multiple commodities is also frequently considered: In Richards and How (2002), a strongly simplified continuous approach for combined trajectory planning and collision avoidance is presented. Li et al. (2020) study unmanned aerial vehicle (UAV) trajectory planning in a two dimensional grid using an ant-colony algorithm and safely separated routes for ship traffic are for example considered in Lübbecke et al. (2019). A model for safety separation due to headway constraints in railway traffic is also presented in Caimi et al. (2018).
The substantial body of work on coordinated multi-robot path planning should also be mentioned in this context. Gawrilow et al. (2008) perform optimization of trajectories for multiple commodities with a strictly restricted network and conservative reservation of space to avoid conflicts. Yan et al. (2013); Wagner and Choset (2015) give a broad overview over related literature in this field. In a time-expanded setting, Hoch et al. (2020) study the theoretical complexity of the non-stop disjoint trajectories problem.
Two very recent publications from the field of aircraft trajectory optimization are closely related to the approach and the application presented here: First, Borndörfer et al. (2021) develop a discrete-continuous approach for single trajectory optimization. They provide criteria on graph 'density' that allow to determine tight error estimates for the obtained discrete reference trajectory if post-processed in a continuous optimization step. Second, in the context of unmanned aerial vehicles (UAVs), Schmidt and Fügenschuh (2021) present a mission planning and trajectory optimization approach for given wind fields and convex no-fly zones. They use a linearized discretization of the equations of motion to determine trajectories for the UAVs.
Contribution: We present a novel concept to ensure disjointness of trajectories on time-expanded graphs. This is used in an integer programming model for disjoint trajectories optimization that can be solved to global optimality with respect to the underlying graph. Further, we develop a new method to generate arbitrarily fine representations of complex, possibly non-convex, three-dimensional space environments. Combined with the presented general adaptive-refinement framework, this modeling enables free planning of trajectories. Advanced algorithmic engineering is applied to limit the fast growing instance sizes, which arise from a finer 1 3 An integrated rolling horizon and adaptive-refinement approach… discretization of space and time. We thus integrate the adaptive-refinement in a rolling horizon approach.
We concretize this general framework for the real world application of combined trajectory optimization and runway scheduling in the surroundings at airports. Computational experiments with realistic scales of up to ten aircraft demonstrate efficiency of the algorithmic approach.

Disjoint trajectories -general formulation
In this section, we describe a general graph-based approach to determine disjoint trajectories for multiple commodities. We also introduce the necessary basic notation and concepts used throughout the remainder of the manuscript. An overview of the notation used can be found in Table 1.
First, we give a formal description of the problem as an extension of the disjoint trajectories problem introduced in Hoch et al. (2020). An integer programming formulation is subsequently derived.
To determine trajectories for commodities, we consider a base-network to be given by a directed graph G = (V, A) without multiple edges. Trajectories are paths in G followed over time. Therefore, consider a planning horizon T and a set of time-steps T∶={0 = < < < ... < = T} to be given. Consider a function ∶ A → I ∪ � , with I being the set of closed intervals of ℝ + 0 . This Table 1 Notation for the disjoint trajectories problem as used in Sect. (2) Symbol Description Indicator variable whether commodity 1 ≤ i ≤ k uses arc ∈ A T traversal-time function gives the possible traversal times for each a ∈ A . We can thus explicitly state the time-expanded graph G T = (V T , A T ) with Depending on the situation, we use (v, t) ∈ V T or a bold print ∈ V T to refer to vertices in the time-expanded graph. Similarly we use ∈ A T for arcs. As a trajectory we consider a path W = (u 1 , t 1 ), (u 2 , t 2 ), ..., (u , t ) in the timeexpanded graph G T without direct returns, that is u i ≠ u i+2 for all 1 ≤ i ≤ − 2 . In a slight abuse of notation we also write W ∈ G T to imply that the set of nodes of W are a path in G T .
We consider k ∈ ℕ (heterogeneous) commodities with individual connection requests. For each commodity 1 ≤ i ≤ k the corresponding connection request is given by a source-destination pair s i , d i ∈ V × V together with release time r i ∈ T , where the corresponding trajectory has to start at s i . Further, for each 1 ≤ i ≤ k earliest resp. latest times of arrival i , i ∈ ℝ for the commodity at its destination are given. As commodities are heterogeneous they are assigned individual arc costs c i ( ) ∈ ℝ ∪ {∞} for each ∈ A T . We say that a trajectory W satisfies the request of commodity the associated costs for commodity i are finite.
To determine disjointness, we introduce a generic definition: For each pair of commodities 1 ≤ i, j ≤ k, i ≠ j and each arc ∈ A T , consider a (possibly empty) set of separation-arcs S(i, j, ) ⊆ A T . This set of arcs must not be used by commodity j if commodity i uses .
Thus, we call trajectories W(i), W(j) of commodities 1 ≤ i, j ≤ k, i ≠ j disjoint if, for all arcs i implied by W(i) , no arc of W(j) is in S(i, j, i ) and vice versa.
Definition 1 (The disjoint trajectories problem (DTP)) Let a time-expanded network G T = (V T , A T ) , k commodities with connection requests (s i , r i ), (d i , i , i ) , individual arc costs functions c i and separation-arcs S(i, j, ) be given.
The disjoint trajectories problem (DTP) is to determine trajectories W(i) ∈ G T for each 1 ≤ i ≤ k that satisfy the corresponding connection request and are pairwise disjoint with respect to the separation-arcs S.
Defining disjointness by using separation-arcs enables a high flexibility and enormous modeling power. In the spirit of alternative arcs as in Pacciarelli (2002), it allows modeling of precedence and disjunction constraints that are only base-network dependent. In addition, time dependent links between resources of the basenetwork can be depicted. These cover, e.g., traversal speed dependent constraints, as well as day-time dependencies. Further, headway constraints on arcs as modeled in e.g. Caimi et al. (2018) are a special case of separation arcs.

Integer programming formulation:
To solve DTP, which in general is NP-complete ( Hoch et al. (2020)), we develop an integer programming formulation in the following.
We introduce binary variables x i ( ) for each 1 ≤ i ≤ k and ∈ G T that indicate whether commodity i uses arc ( x i ( ) = 1 ) or not ( x i ( ) = 0 ). Each individual commodity 1 ≤ i ≤ k has to fulfill an extended flow-property on G T ((1a)-(1b)). This differs from classical flow formulations in that not a single sink is possible, but one of the feasible time-steps to arrive at the given destination must be reached. Hence, in the time-expanded network the possible destinations are all vertices (d i , t i ) with i ≤ t i ≤ i . This is formalized in (1c).
Avoidance of direct returns of commodities is assured by Disjointness for two commodities 1 ≤ i ≠ j ≤ k can be represented by the following set of constraints: The linear objective function is given by Now, for any feasible solution to (1a)-(1f) with finite objective value (1g), the variables x i directly indicate a trajectory in the considered sense. As disjointness with respect to the separation-arcs is also assured by constraint (1f), the model fully describes the DTP as stated in Definition 1.
Note: There is an alternative formulation to ensure disjointness: (2) where M is a sufficiently large real number that is chosen so as not to impose any restrictions on the trajectory for j if x i ( i ) = 0 . Depending on the value of M this formulation in general leads to weaker relaxations and should be used carefully, see e.g. Bonami et al. (2015).
Remark 1 Waiting at vertices is covered by the presented formulation. This is possible if there are loops in the corresponding base-network. Respecting the maximal waiting time at a vertex is directly enforced by the condition that avoids direct returns: Choosing successive waiting arcs of the same vertex violates constraint (1e).

Remark 2
We observe that the presented model can also cover the use-case where a commodity may have multiple possible sources or destinations. We give a short explanation for multiple destinations with the well known concept of adding an artificial depot to G. Here is considered as new destination for all requests and artificial arcs from all originally possible destination vertices to that depot are added. As expanding the artificial depot is technically not necessary, it can also be added to the time-expanded graph instead of the base graph with appropriate arc costs. In that case, constraint (1c) has to be changed to: Note: By adding an artificial depot to the time-expanded graph constraints (1c) can generally be reformulated in that way.

An iterative meta-algorithm
Increasing the degree of detail in the base-network G or the number of considered time-steps can drastically increase the problem size for Problem (1) on the timeexpanded graph G T . Hence, for many applications, consideration of the complete underlying base-network in a sufficiently fine time-resolution exceeds the current technical capabilities. It is thus often beneficial to consider a restricted graph and to adaptively refine only when necessary.
We will present a meta-algorithm in the spirit of adaptive-refinement methods. Such methods have successfully been applied to solving general linear programs, as by Gleixner et al. (2016), and even more complex settings such as non-convex optimization, see e.g. Kuchlbauer et al. (2020);Geißler et al. (2012).
A plethora of studies for adaptive refinement, or coarse to fine approaches exist that are more tailored for specific applications: A prominent example of the successful application of this concept is shortest path planning in road networks via contraction hierarchies. For applications in air-traffic management, authors in Blanco et al. (2016) showed that for even one trajectory A * algorithms outperform contraction hierarchies. Fischer and Helmberg (2014)

3
An integrated rolling horizon and adaptive-refinement approach… instance sizes resulting from time-expanded networks by developing a dynamic approach for shortest path computations. They embed this into train timetabling for conflict-free schedules. One key assumption for their shortest path computation -that is not fulfilled in our setting in general -is a preference for early arrival.
A graph aggregation approach for maximum flow computations is also successfully considered in Bärmann et al. (2015). For general flows, Liers and Pardella (2011) show that reducing the original graph and finding good initial solutions has great positive effects. A column generation approach to determine refinement actions based on heuristic reduced cost computations is successfully applied by Borndörfer et al. (2014) and subsequent works that build on that approach.
For the DTP, none of these approaches are easily extendable to deal with multiple commodities and consider any kind of separation constraints, especially if full information about the finest level is not available or not desired in storage. This also holds true for other advanced path-planning algorithms such as RRTs or D * algorithms (see e.g., Yang et al. (2016) for a broad survey) that work well for single commodities.
Hence, to deal with the necessary separations and arc dependencies, we introduce the following iterative meta-heuristics: For each commodity 1 ≤ i ≤ k and its optimal trajectory W(i) , this function returns a new sub-graph G � i T ← N(W(i)) to be considered for re-optimization. 4. Solve the reduced Problem (1) resulting from the new sub-graphs G ′ i T obtained in Step 3.
Step 3 can be seen in the spirit of general neighborhoods in graphs, as they are introduced in Mladenović and Hansen (1997). In Section 5 we give an explicit construction that refines space and time resolution in each iteration. This is performed in a geometrical neighborhood of the current trajectory.
The idea is formalized in Algorithm 1. Its running time is determined by the number N ∈ ℕ of iterations, and solving the restricted versions of Problem (1). This depends on the size of the initial graph, as well as the selected neighborhood graphs. As all the restricted problems are also generalized instances of the NP−hard k disjoint paths problem, we cannot expect to find polynomial time algorithms to solve the restricted problems. In the worst case, we have to expect exponential running times with respect to the number of arcs in the neighborhood graphs and the number of commodities considered.
Remark 3 Algorithm 1, although an heuristic, provides an important global point of view: The initial restricted problem is solved to global optimality with respect to the considered sub-graphs. Further, a feasible solution of the initial iteration will also be feasible for the full instance. For subsequent iterations and neighborhood functions, the problems will again be solved to global optimality with respect to the considered sub-graphs, and feasible solutions also will be feasible in G T .
While in general no optimality guarantees with respect to the full instance on G T can be given, there is the following that can be ensured: We say that an algorithm has an improvement property if, starting from a feasible solution, the objective value is not worsened in each iteration.
Lemma 1 Let a number N of iterations for Algorithm 1 and a set of neighborhood functions N , 1 ≤ ≤ N on a time-expanded graph G T be given. If for each 1 ≤ ≤ N and any valid trajectory W in G T it holds that Algorithm 1 has the aforementioned improvement property, if the initial restricted problem was feasible.
The proof of Lemma 1 follows directly from the condition that the previously optimal solution is still a valid solution.
Note: An algorithm engineering issue that has to be considered when applying Algorithm 1 is the following: A careful selection of the initial restricted subgraphs is necessary to avoid infeasible sub-problems where the overall problem would be feasible. Whilst there are no general rules on how to select initial subgraphs and neighborhood functions for arbitrary instances, in Section 4 and Section 5 we describe a technique that led to promising results in the computational experiments conducted in Section 7. If feasibility should be granted, dummy arcs with prohibitively high costs can be added to the original time-expanded graph and considered in each sub-graph.

Modeling for free commodity movement
The model described in Sect. 2 can be used for general trajectory planning for multiple commodities. Thereby, it can not only be used in classical road-like networks, but also in environments where free movement in the Euclidean plane or threedimensional space is allowed.
In this section we present a discretization method for trajectory optimization in a free, three-dimensional environment. In the process, we introduce an approach that allows the necessary information about arcs in the base-network to be derived properly, without having the need to explicitly store the base-network completely. We also explain how the possible arc traversal-times for different commodities are determined to enable time expansion. This approach can be used for several applications, including robot motion planning, trajectory optimization in air traffic management or for naval traffic. Table 2 gives an overview of additional notation introduced and used throughout the following sections.

Spatial discretization -the base-network G
We consider a cubic portion of three dimensional space to be represented by a grid of space-vertices. Further, we consider a rectangular coordinate system and neglect earth curvature. The vertices representing space are evenly spaced grid points with uniform grid-step distances p ∈ ℝ in the horizontal plane and h ∈ ℝ for the vertical height levels. Arcs represent straight line connections between two vertices. In this way, the base-network can directly carry information about arbitrarily shaped forbidden zones, such as topography, building or no-fly zones in air traffic management.
As the base-network itself can become intractable for fine resolutions of p, h and a large portion of space considered, we use an implicit representation. We Maximal climbing reps. sinking angle therefore consider the respective volume to be divided into regular three-dimensional cubes. For each of these cubes we store the information on whether or not it is a forbidden zone in a three-dimensional environment matrix M(G) . In this way, for each point in the considered volume, the corresponding cube and thus the environment information of M(G) can be determined. If desired, the entries of M(G) can carry further information such as weather or noise exposition. For any straight line connection between two vertices v, w ∈ V in the base-network, the necessary arc parameters can now be derived from M(G) . To achieve this, we have to determine which cubes are affected by that line segment and check the corresponding information in M(G) . If the information implies that the arc does not cross restricted zones, it is considered in the base-network. Figure 1 illustrates an example for crossed cubes. The filled circles mark the start and end position of the segment under consideration, while the gray-filled squares mark the cubes that are crossed by the segment.

Remark 4
In applications, the set of cubes crossed by a line segment is not necessarily the set of cubes affected by that segment: -If commodities are not treated as points, the dimensions of the commodities may actually be considered to obtain the affected cuboids. Therefore, the amount of affected tiles may increase compared to Fig. 1. -If, for a particular application, the resolution of the environment matrix is substantially finer than the base-network resolution, an inexact approach may be sufficient. See Sect. 7, Fig. 7 for an example.
For this general modeling section, we therefore do not provide an explicit procedure for determining the affected cuboids of a line segment in M(G).
For practical reasons we also do not consider every possible space-vertex pair in the base-network for connections, but rather consider a maximal Euclidean arc-length D max and a maximal climb resp. sinking angle ∈ [0, 2 ] . Figure 2 illustrates these rules.
Overall, this method of modeling leads to a highly flexible approach in several ways: -The degree of detail in the environment description can be controlled by the resolution of the environment matrix. With M(G) provided, the arc-information for base-networks with arbitrary resolutions can be properly derived. -The number and shape of restricted zones has no influence on the modeling and no direct influence on the running time of the algorithm. This is particularly useful for e.g. modeling airspace where mountains as well as cities imply an enormous variety of classically non-convex no-flight zones. -A full a priory storing of G with all arcs will consume considerably more storage than considering environment in M(G) . This allows an even finer representation of space in M(G) being used than that needed for the desired degree of detail in the base-network G.
Remark 5 This representation of the base-network will be especially useful for the iterative approaches later detailed. It ensures that arcs that are only relevant in a finer discretization only are considered and evaluated when they are really needed in a sub-graph, and not for the entire space, that may be visited.

Time expansion
For time expansion we consider, as before, a given horizon T and an additional timestep t . This is used to construct the set of considered time-steps as As arc-costs in the time-expanded graph are commodity dependent, we also determine the possible traversal-times (a) for a = (u, v) ∈ A of the base-network commodity dependent. Therefore, for each commodity 1 ≤ i ≤ k , we use a black box cost function i ((u, t 1 ), (v, t 2 )) → ℝ ∪ ∞ . This function shall be built such that Note: If in addition to general environment information there is also timedependent information, this can be depicted by an (additional) four-dimensional environment matrix. In that case, the cost function can be designed to incorporate necessary checks analogously to the three dimensional case described above.

Separation-Arcs
The determination of separation requirements between commodities are highly application dependent. As the separation-arc sets do not affect the general discretization we leave it open for the reader to determine those according to their application. If any vertex related restrictions have to be fulfilled, they can be ensured by considering corresponding couplings between the in-and out-going arcs.
For an explicit example of the construction of separation arcs, see Sects. 6 and 7.

Iterative adaptive-refinement algorithms
In this section we present two iterative adaptive-refinement algorithms for the DTP in the setting described above. First, we present a static version where all trajectories are fully considered for refinement at once. This directly follows the idea of Algorithm 1 and is formalized by stating the neighborhood functions used. Second, we introduce a rolling horizon approach for the refinement iterations. This follows the same neighborhoods as introduced for the static version, but moves the part of the trajectory to be refined through time.
For both algorithms, the idea is to start with a coarse discretization of both space and time to obtain initial trajectories fulfilling all separation constraints. Afterwards, a four-dimensional refined neighborhood 'tube' is considered around each trajectory. In this tube, space and time are finer discretized. In this way more detailed trajectories can be obtained while the problem size in each iteration remains tractable.
Both algorithms start from the same initial sub-graph, which is why we explain how to determine this first. Afterwards, we introduce the general approach for the iterative refinement along complete trajectories, and the rolling horizon algorithm. This can enable 'larger' neighborhoods, whereby larger can either imply consideration of larger deviations from the original trajectory or a finer discretization. We conclude the section by stating some properties of the described neighborhood search and providing some notes on how to avoid undesired side-effects.
1 3 An integrated rolling horizon and adaptive-refinement approach…

Initial sub-graph
To determine the initial sub-graph G ′T i for each aircraft 1 ≤ i ≤ k , we assume that the environment matrix M(G) is given. We start at the corresponding source s i as 'anchor' and build a three dimensional grid with step length p 0 in the plane and height-level-step of h 0 . Thereby, p 0 has to be a multiple of the desired finest resolution p , and h 0 a multiple of h.
We consider all vertices whose positions are in one of the volumes covered by M(G) . This does not exclude the possibility that the destination d i is not part of the built grid. We therefore add d i to the resulting space grid. For time expansion, we also start at r i and use increment t 0 , which also has to be a multiple of the desired finest resolution t.

Complete trajectory refinement algorithm
Having obtained an optimal solution to the previous problem, we refine the search space in proximity to the resulting optimal trajectories by defining appropriate neighborhood functions. To do so, in the following we present those neighborhood functions N for each refinement iteration 1 ≤ ≤ N . They can be described by a set of parameters: -Refined step lengths p , h , t . These have to be selected such that the previous step lengths are multiples of the new step length, and the current step lengths remain multiples of the finest desired one. That is, * −1 ∕ * ∈ ℕ has to hold for all three step lengths, as well as * ∕ * ∈ ℕ . Again it is not necessary fo all dimensions to be refined by the same factor. -Maximal deviations p ∈ ℝ in the plane, h ∈ ℝ in height and time deviation of t ∈ ℝ. -A maximal arc-length D max, ∈ ℝ in space.
For now, we first consider one arc and explain how the parameters are used to construct the vertex set of the search-neighborhood. The neighborhood function uses this concept along all arcs of a given trajectory and adds arcs between them following the rules described before.
First, we construct the vertex set of the time-expanded graph to be considered. Therefore, let an arc ∈ A T be given. To determine the nodes in the finer discretization that best follow the given line segment, we use a four-dimensional extension of Bresenhams line drawing algorithm Bresenham (1965). We call the resulting node set L ( ) ⊂ V T line nodes, see Figure 3 for an illustrative example in two dimensions.
These nodes are now used to determine the set of all 'refinement' nodes to be considered in the neighborhood function. To the set R ( ) ⊂ V T of so called refinement nodes for each ∈ L ( ) we add all nodes with b ∈ ℤ 4 such that Figure 4 illustrates for each line node the set of added nodes in the resulting box.
We use this constructive approach for the vertices in the neighborhood graph instead of determining nodes in the finer resolution that would have a maximal given distance to the arc. This avoids unnecessary checks for very remote vertices.
To obtain a proper neighborhood function N for each arc in a trajectory W , we use the procedure described above to determine the full set of refinement nodes R (W) = ⋃ ∈W R( ) . This yields the node set of the resulting sub-graph N (W) . For two 1 , 2 ∈ R (W) , an arc = ( 1 , 2 ) is added if ∈ A T and the Euclidean distance in three-dimensional space is not larger than the maximal arc length of D max, ≤ D max . Finally, we add the arcs of W to the arc set. This ensures that the neighborhood functions have the desired improvement property as stated in Lemma 1.

Rolling trajectory refinement algorithm
We now integrate the approach into a rolling horizon framework. This can allow for larger deviations or finer resolutions from one refinement iteration to the next compared to the static version. Rolling horizon has successfully been applied for many problems, where problems get to large to deal with the whole instance, see e.g. Bean and Smith (1984); Chand et al. (1990) or more recently Fattahi and Govindan (2022); Glomb et al. (2022).
In addition to the parameters given in Table 3 to determine the refinement neighborhood graph, for each iteration 1 ≤ ≤ N we now consider a planning horizon H and a horizon step H . We are then able iterate over the overall planning horizon T in steps of H and consider a time interval of length H for refinement. Only line-nodes (v, t) ∈ L with t within the refinement interval are considered to determine the refinement nodes. For all time-wise earlier parts of the trajectory, the trajectory computed thus far has to be maintained. For all time-wise later parts, the three-dimensional trajectory in space must remain the same, but we allow to accommodate possible time deviations from the previous trajectory. Therefore, a maximal temporal deviation of t on when to arrive at the corresponding vertices of the base-network is allowed. Figure 5 gives an impression of how the refinement period is determined for a trajectory projected onto the x-axis and time, and indicates how to refine this trajectory.
The exact procedure, and especially how to connect the 'fixed' parts earlier and later than the refinement horizon with the refined search graph, is detailed in the following.
Let t s be the starting time of the current refinement period. This implies that the end of the refinement period is t e = t s + H . Further, as in the previous section, let p , h , t , p , h , t and D max, be given. Any given trajectory t 1 ), ..., (v , t ) is now split into up to three parts. If t 1 < t s < t e < t holds, those are: p 0 ∈ ℝ Grid-step length in the plane for initial sub-graphs h 0 ∈ ℝ Height-level-step length for initial sub-graphs t 0 ∈ ℝ Time-step increment for initial sub-graphs Step lengths in iteration p ∈ ℝ Maximal deviation in the plane in iteration h ∈ ℝ Maximal height deviation in iteration t ∈ ℝ Maximal time deviation in iteration D max, ∈ ℝ Maximal three dimensional arc length in iteration H Refinement period length in the rolling adaptive-refinement algorithm H Step width in the rolling adaptive-refinement algorithm L ( ) ⊂ V T line-nodes for neighborhood function R (W) ⊂ V T refined node set for a trajectory W in iteration ( ) ⊂ V T time-deviation nodes for a node ∈ V T in iteration t s , t e ∈ ℝ Start, resp. end time of the refinement period in the rolling horizon approach W A , W , W Trajectory segments that are fixed ( W A ), fully refined ( W ) resp. with temporal adaptions ( W ) in the rolling horizon approach 2. W = (v p 1 , t p 1 ), ..., (v p 2 , t p 2 ) such that t p 2 −1 ≤ t e ≤ t p 2 3. W = (v p 2 , t p 2 ), ..., (v , t ).  An integrated rolling horizon and adaptive-refinement approach… If t s < t 1 ≤ t e holds, W A is empty and W starts with (v 1 , t 1 ) . If t 1 > t e , both W A and W are empty and W = W. If t s < t < t e holds, W ends with (v , t ) and W is empty and if t < t s , both W and W are empty and W A = W.
Ad W A ∶ This part of the trajectory is no longer considered in the planning horizon; nodes and corresponding arcs are directly added to the search graph.
Ad W ∶ We construct the set of refinement nodes for this sub-trajectory similarly to the previously described construction. For all p 1 + 1 ≤ n ≤ p 2 − 2 we use the same procedure as in the previous section to determine the set of refinement nodes corresponding to arcs (v n , t n ), (v n+1 , t n+1 ) . If p 1 = 1 , we do the same for (v 1 , t 1 ), (v 2 , t 2 ) and, if p 2 = , we do the same for (v −1 , t −1 ), (v , t ).
Special attention has to be paid to the first and last arcs of W if they are not the beginning or the end of W: If p 1 ≠ 1 , determine the set of line nodes W.l.o.g., we assume the line nodes to be indexed such that t 1 from the set of line nodes. For those remaining, we add nodes to the current set of refinement nodes as described in (3). Further, we add (v p 1 , t p 1 ) to this set.
For the last trajectory segment in W , we proceed analogously if t p 2 < t . Finally, we consider a set of so-called time-deviation nodes for (v p 2 , t p 2 ) . These connect the refinement graph to the last part of the new sub-graph, where the three-dimensional routing is fixed but adjustments in the traversal times are possible. The time-deviation nodes are given by Finally, we add ((v p 2 , t p 2 )) to the set of refinement nodes R . We add all arcs ( 1 , 2 ) ∈ A with 1 , 2 ∈ R with an Euclidean distance no greater than D max, in space to the refinement graph.
Ad W ∶ For each arc in W , we iteratively consider the time-deviation nodes for each arc. That is, for n = p 2 , ..., − 1 we consider ((v n , t n )) and ((v n+1 , t n+1 )) and add all possible arcs between them to the refinement graph. Again, to ensure the improvement property of Lemma 1, we add the arcs of W to the arc set. This procedure fully describes how the new restricted sub-graphs are determined by the neighborhood function N ,t s . This function now not only depends on the refinement step, but also on the current start time for refinement.
Algorithm 2 describes the full algorithmic procedure. After solving the restricted problem on the initial sub-graph, for each refinement step 1 ≤ ≤ N a total number of ⌈ T H ⌉ iterations have to be performed. Hence, in total restricted Problems (1) have to be solved. The computation times for the restricted problems again depends on the size of the neighborhood graph, implied by the parameters of Table 3.
Note: By construction, the improvement property of Lemma 1 holds for Algorithm 2, as well as for the static version. Still, for arbitrary costs we cannot give any optimality guarantees of the obtained trajectories with respect to the finer discretization: Example 1 Consider a two-dimensional base grid as given in Figure 6. Assume the black vertices to be in the initial sub-graph, unit traversal times for all arcs and cost as follows: Solid black lines: M , snaked black lines M − 1 and dashed gray lines 0. For the initial graph, sum the cost of the arcs of the more finely resolved graph.
The number of arcs in the path, and hence the arrival time, shall have no influence on the cost.
An optimal trajectory in the initial graph would now go straight to the right and then up to d, with costs of 2M − 2 . Now, for any neighborhood function that has p ≤ 3 , this will remain the optimal trajectory, whilst for the complete instance the optimal solution would have value 0 by first moving left unit the second to last vertex, then up and finally right again.
Despite this, there are a few pitfalls to this approach that can be avoided: -Initial discretization: Ensure that D max ≥ p 0 holds. Further, time-step t 0 has to be chosen such that arcs of the (restricted) base-network can actually be traversed by commodities in times that are implied by the time discretization. This ensures a non-empty time-expanded graph. The finer that t 0 is with respect to the maximal arc length in the base-network and the maximal commodity speed, the more speeds to traverse the arc are in the model.
-Base-Network: Ensure that the initial discretization is chosen sufficiently fine to avoid having a false infeasible instance. Avoidance of restricted areas can also become prohibitively conservative.

Free-flight trajectories and runway scheduling
The model and algorithms presented so far are designed for general trajectory optimization and a three-dimensional space with free movement. In this section we apply this to approach planning at airports, where trajectories and runway scheduling have to be considered jointly. We first provide a short overview of the specifically related literature. After, we explain how we can naturally include runway scheduling in the graph-based setting. Finally, we discuss the objective chosen to explicitly and jointly consider both fuel-costs and delay, while ensuring a fair distribution of costs among the involved aircraft. The additional notation used in this section can be found in Table 4.

Related literature
A vast body of literature focuses on runway scheduling as an isolated problem. Beasley et al. (2000) introduce a basic scheduling model, that has been used as a basis for many more recent approaches. An exhaustive survey on runway scheduling can be found in Bennell et al. (2011). Even robust approaches to runway scheduling under uncertainties have been intensively studied, see e.g., Heidt et al. (2016). Some approaches have already been developed for the combined optimization of arrival sequence and trajectories at airports: Grüter et al. (2016) present a bi-level approach to solve the sequencing and trajectory optimization problem by combining gradient based methods and genetic algorithms. Toratani (2016) also treats Fig. 6 Instance for Example 1; dotted path (red) with costs 2M − 2 optimal in initial graph; dash dotted path (blue) with costs 0 globally optimal trajectory and sequence optimization simultaneously, ensuring only safe arrival times, but not en-route safety distances. Building on works of ; Grüter et al. (2016); Bittner et al. (2015)), Semenov and Kostina (2020) provide another hierarchical approach to the problem, which shows good performance. All these approaches to combined optimization of trajectories and runway scheduling include local or hierarchical optimization approaches, leading to possibly fast, but not necessarily globally optimal solutions. In the following we explain how runway scheduling can be included in the DTP by using the modeling trick explained in Remark 2.

Runway scheduling
If runway scheduling has to be considered in addition to trajectory planning, the general model given in the previous section can be adapted to naturally fit the framework of the model, as described in Sect. 2. Thereby the methodology is flexible enough to not be fixed for the final approach fixes, but if desired any fixed point in the approach procedure that may allow for sufficient manual control.
Consider a set F ⊂ V of possible final approaches. We add an artificial depot to the set of time-expanded vertices V T of the base-network and connect each possible To model the necessary safety separations between final approaches, let d 1 ,d 2 (i, j) denote the minimal separation time between a leading aircraft 1 ≤ i ≤ k and trailing aircraft 1 ≤ j ≤ k if aircraft i uses approach d 1 ∈ F and j approach d 2 ∈ F.
For any pair of aircraft 1 ≤ i, j ≤ k , each d ∈ R and any time-step t ∈ T , the following constraint has to hold: Horizontal safety distance for commodities v ∈ ℝ Vertical safety distance for commodities Weighting factor for fuel-costs and time-deviation costs ∈ {0, 1} Indicator whether continuous descent shall be enforced by the model or not 1 3 An integrated rolling horizon and adaptive-refinement approach… Note: Constraint (5) is of the form of constraint (2), with M = 0 . In contrast to the general spatial separation constraints, we add the runway scheduling constraints directly to the model, and not as feasibility cuts as for the general en-route safety distances. This is implied by Nikoleris and Erzberger (2014), where computational tests show that 80 − 90% of en-route conflicts in the terminal maneuvering area are avoided by proper arrival spacing.

Remark 6
From an algorithm-engineering point of view, the possible separation times d 1 ,d 2 (i, j) should be considered when deciding on an initial time-step increment t 0 : For an aircraft-runway combination (i, j), is large, this can lead to solutions 'far away' from the optimal solution of a finer resolution.

En-route separations
For en-route separations, we assume that either a minimal horizontal separation of h ∈ ℝ or a minimal vertical separation of v ∈ ℝ have to be maintained between any two commodities at any given point in time. To determine commodity position on the trajectory, we assume that arcs are traversed with constant speed. Directly considering the full set of separation-arcs and the implied constraints would add many unnecessary constraints (trajectories may be in different areas of the considered space). Hence, we do not pre-compute these, but add them dynamically as feasibility cuts in the solution process. For more details see Sect. 7.

Objective
In air traffic management, stakeholders have different objectives for optimization that have to be considered in a fair way. This ensures everybody's needs are satisfied and high acceptance rates are reached among multiple stake holders (e.g., airport, airlines, citizens). We consider the following values explicitly: 1. Cost induced by fuel consumption 2. Deviation from the scheduled time If further values such as noise exposition to citizens shall be included, two possible approaches are available: The value can be incorporated as a direct cost, following an approach as described in the following. Alternatively, an indirect approach can be adopted, excluding arcs from the network if using them violates a certain threshold. The necessary data could again be stored in the environment matrix. Fuel Consumption : To determine fuel consumption for a trajectory, we consider the cost function i ( ) as a cost estimator for the optimal fuel cost for this arc. For this work we assume that i ( ) depends on the initial vertex, the heading angle, length of the segment and the time to traverse the segment, i.e. assuming fixed and uniform weather conditions. If locally or temporally changing weather should be included this can be achieved by using fuel cost estimators that consider further edge information that may be derived from information stored in M(G) . For the given application of arrival planning in the terminal maneuvering area, we neglect changes in the aircraft mass due to fuel burn along the trajectory, as the relative change in mass will only be approximately 0.5 %.
Deviation from planned arrival time : For each aircraft 1 ≤ i ≤ k approaching an airport, a scheduled arrival time t ST i is given. At this scheduled time, the aircraft is expected to arrive at any of its possible approaches d ∈ F i and deviating from this time incurs direct and indirect cost. We assume that being delayed is more costly, and thus we use the following assessment of cost: We set a linear penalty for being early, whilst quadratically penalizing delay.This ensures that one heavily delayed aircraft is more expensive than several slightly delayed ones. For each d ∈ F i and t ∈ T with t ≤ i , we use the following cost function to price deviation from arrival time: for some scaling factor ∈ ℝ.
This choice to assign cost based on derivation from arrival time ensures fairness in the sense that delaying one aircraft for 30 time units is more costly than delaying tree aircraft by ten time units each. Thus, the optimal solution avoids assigning a long delay to only one aircraft, where possible.
Similarly, as time expansion of space arcs was individually performed for each aircraft, for any approach d ∈ F�{F i } as well as for times t > i we do not consider the arcs ((d, t), ) for aircraft 1 ≤ i ≤ k . That is, we assume i ((d, t), ) = ∞.
Overall, the objective shall give the possibility to depict a weighted sum of total used fuel costs and time deviation cost . Therefore, we introduce the weighting factor ∈ [0, 1].
Arc costs in the time-expanded network for arc = ( , ) ∈ A T for aircraft 1 ≤ i ≤ k are given by: Continuous descent: Climbing represents unfavorable behavior for approaching aircraft, both due to fuel consumption and noise exposition. Thus, in the construction

3
An integrated rolling horizon and adaptive-refinement approach… of arcs in the base-network continuous descent operation can be enforced by only considering descending or height-preserving arcs. Continuous descent is induced by the parameter . Further, enforcing continuous descent in the terminal area significantly reduces the size not only of the base-network, but also of the resulting optimization problem. Note: Preliminary studies empirically showed that if only one aircraft is considered, the optimal trajectories obtained through the approach follow the continuous descent paradigm, even tough it was not enforced in those experiments.
Remark 7 (Further modelling assumptions) The described modeling implicitly includes some further modeling assumptions. First, we do not assume a given earliest time of arrival for flights. This time is implied by the underlying time-expanded network and the traversal times. Further, the model allows for arbitrary speed changes, and only direct returns are forbidden. If speed changes or turning should be restricted, this could be incorporated in the presented formulation by introducing 'separation'-sets S(i, i, ) that include all arcs forbidden to aircraft 1 ≤ i ≤ k if ∈ A T is used.
For the considered instances in Section 7, the turning behavior of aircraft was always reasonable with respect to the given network and direct returns to the previous space position never occurred.

Computational experience
In this section, we first describe some further algorithmic details of the implementation before presenting results for the class of circle instances and for some more realistic instances with individual release times of aircraft.

Implementation details
The model described in the previous sections is implemented in python 3.6. For the solution of MIPs, we use the solver GUROBI Gurobi Optimization, LLC (2022) version 9.1.0.

Feasibility cuts for disjoint trajectories:
Using GUROBI allows the use of so-called lazy constraints, i.e., constraints that are added as feasibility cuts on demand, as mentioned in Sect. 6.
The solution process is started using only the flow and scheduling constraints as described earlier. Constraints of type (1e) and (1f), that refer to direct returns and en-route separation are omitted in the beginning. Now, if an optimal integral feasible solution (with respect to the current set of constraints) is found, direct returns and conflict freeness are checked. While it is straightforward to check direct returns to a previous space vertex, the check for conflicts is a done as follows: We compute aircraft positions in steps of one second and pairwise check whether the minimal safety distances are maintained. If this is the case and no direct returns are observed, the solution found is globally optimal with respect to all possible separation constraints. If any direct return or conflict of two aircraft is detected, the corresponding constraints of type (1e) or (1f) are added and the model re-optimized.
In the worst case, this will add all constraints of type (1e) and (1f), slowing computations. However, but experiments imply that the number of added en-route safety separation constraints is indeed marginal if a feasible runway scheduling is ensured.
Note: If a conflict is shorter than one second this is not detected that way, but for the considered application this is negligible. If a fine time resolution for determination of conflicts is necessary, further advanced methods can be applied to identify conflicting arc pairs.

Fuel cost model and environment matrix
To compute the fuel costs, we use a blackbox fuel cost model kindly provided to us by the Institute of Flight System Dynamics at TU Munich ( Schweighofer et al. (2022)). The model is an analytical approximation to the response surface of an Optimal Control Problem formulated with a 3-DoF aircraft model based on BADA 3. It represents the cross-couplings between distance, flight path angle and time, as well as the effects of initial altitude and aircraft mass, and is restricted to an approximate/conservative flight envelope. In its current form, the model captures only the basic characteristics of aircraft behavior, and cannot be relied upon for quantitative analysis of flight plans.
The environment representation used in Subsection 7.3 is also kindly provided by the Institute of Flight System Dynamics at TU Munich. It considers a medium-sized German airport and the no-fly zones are derived from publicly available information on topography and restricted air-spaces.
To determine the affected cubes, for this application we use generalization of Bresenhams line drawing algorithm Bresenham (1965). For a two-dimensional illustration see Fig. 7. Note: Using this approach, we cannot exclude the possibility of slight inexactness, in the sense that small segments of the arc may traverse a cuboid that would be classified a no-fly zone in M(G) . However, for the resolutions of the base-network and the environment matrix later used, this slight inexactness is negligible.

Circle instances
Circle instances where commodities have their start and destination position on a circle are often used in the literature for two-dimensional conflict resolution and avoidance Rey et al. (2015); Dias et al. (2020); Schmidt and Fügenschuh (2021). In these studies, homogeneous aircraft are uniformly placed on a circle and all head to the opposite side through the center.
We adapt this concept for our three-dimensional problem, including runway separation as follows: On flight level (FL) 230 we consider a circle with radius of 90 km and place a number of n aircraft evenly on the circumference of this circle. The final approach is on FL 50 in the center of this circle. The underlying grid has a step size of 30 km in the plane and 10 Fl in height. The maximal arc length is set to D max = 59 km and the maximal sinking angle to 10 • . Climbing operations are not allowed. Time discretization for the initial run is set to 100 sec. Both release time and scheduled time is 0 for all aircraft and the weighting factor for the delay is set to 0. Hence, only minimal fuel consumption and disjointness are considered. Safety distances are the classical 5 nm in the plane and 10 Fl in height. All aircraft are considered to have medium weight, so runway separation is set to 75 sec for each possible combination of leading and trailing aircraft. No no-fly zones are considered for this setup.
The instances are solved on machines with 2x Intel Xeon E5-2643 v4 CPUs. They have 12 cores/ 24 threads with 3.4 GHz and a RAM of 256 GB. Figure 8 shows time for model building, as well as the time for optimization for the circle instances with 2 − 10 aircraft involved. Also the time spend in lazy constraints is depicted. The model generation covers individual graph generation for each aircraft, as well as deriving the optimization model. With the generic environment, the time for model building grows almost linearly in the number of aircraft. It can clearly be seen that the time for optimization of the model and the time spend in lazy constraints grow non-linearly with the number of considered aircraft. For increasing aircraft number the time to sovle the model becomes the determining factor of the overall running time in this initial iteration.
In Fig. 8 the long running time for nine aircraft is particularly striking. Moreover, the instance of eleven aircraft is not depicted in the figure, where after a runtime of three hours not even a feasible solution had been found and over 250 GB of the available RAM had already been used. Hence, while it is evident, that running-times increase drastically with an increasing number of aircraft, it is implied to grow faster if extrapolating from the nine aircraft instance and slower if extrapolating from the ten aircraft instance. What we can be sure of is that for aircraft number nine there is a large amount of time ( ∼ 1400sec) spent in the GUROBI lazy constraints, i.e. checking for en-route separation and adding those constraints. That the time spend for this increases is to be expected as: -With an increasing number of aircraft appearing in the airspace at the same time, spatial separation along the trajectory becomes increasingly challenging, especially with a coarse space-time discretization. -Similar behavior can be seen for example in Schmidt and Fügenschuh (2021), where circle instances where also considered for planning of conflict-free trajectories for unmanned aerial vehicles. They also reported an exponential increase of running times for an increasing number of aircraft. In their approach, running times of 3600 seconds are also quickly exceeded for similar instance sizes.
However, as the algorithm was designed in the expectation of successively arriving aircraft, we do not further focus on those classical deconfliction instances. Instead, we move to instances where aircraft are released in the considered airspace with temporal separation as would be expected in the terminal maneuvering area around an airport.

Successive release instances
For a set of more realistic instances, we use the matrix M(G) mentioned in Subsection 7.1 on a 200 × 200 km square grid considering a cube size of 1000m × 1000m × 100m . Figure 9 shows the considered airspace schematically. We assume a single final approach on FL 50 ( = 1500 m) in the center of the considered grid in the plane. This final approach can be reached with an unrestricted approach angle. For aircraft numbers from two to ten we generate ten instances each, where aircraft arrive in the considered airspace at any of the four vertices of the plane grid on flight level 200. The selection of corners, release times and originally scheduled times are random, but such that in each initial approach there is a minimal separation of 50 seconds between successively released aircraft. Other than that, the temporal distance in seconds between the release of an aircraft in any initial approach is independently identically distributed (iid.) in [0, 250] . The scheduled time for a flight is computed by iid. selecting a value in [−1800, 1800] (seconds) and adding it to the release time. To avoid unnecessarily large instances, the latest arrival time for any flight is 1800 seconds after its release. For an aircraft number of two to ten we randomly generated ten instances each. We assume all aircraft to be medium weight, which implies a runway separation of 75 seconds.
In the following we analyze the running times for the initial optimization, the static-in comparison to the rolling refinement and multiple refinement iterations. Afterwards, we look at how the deviation from the scheduled times behaves with the different approaches. Finally, effects of refining dimensions sequentially are examined. As the algorithmic concept remains our focus, we refrain from going into further detail on the involved fuel costs model.
All computations shown in this subsection are run on machines with Intel Xeon E3-1240 v5 or Intel Xeon E3-1240 v6 CPUs, respectively. Each of these has four cores with 3.5 GHz each and a RAM of 32 GB.

Running-time analysis
For a running-time analysis of the different refinement approaches we, used the data displayed in Table 5 for the initial discretization.
Initial optimization: In Fig. 10 the running times for the initial optimization of the instances are shown in blue, and the geometric mean over the running times in red. The geometric mean increases only slightly with an increasing number of aircraft and is significantly lower than those presented for the circle instances. The outliers that can be seen for example at an aircraft number of eight, are instances where the temporal separation at the source of the aircraft is close to the minimally enforced 50 seconds. This leads to individually unfavorable maneuvers for aircraft to maintain the necessary safety distances. In a further experiment we looked at the following: What are the solution times for a the problem, when all lazy constraints, that are added through our presented approach, are already included in the model. For the instances with eight aircraft, those solution times where all between 21.4 and 29.95 seconds. This indicates, that it can be beneficial for applications to add the constraints that are violated by the current solution, and further to determine tailored heuristics to add several cuts ad once. This has the potential to minimize the number of re-optimization rounds in the solution process.
Static vs. rolling refinement: For the instances with six aircraft, we compare the computation times for one refinement iteration using the static and the rolling approach. The neighborhood specification of Table 6 for the second iteration, i.e. the first refinement, is used. For the rolling approach we sum the times for model generation and optimization needed to complete one full refinement run. A side by side comparison of the running times for graph generation and optimization for both the rolling and the static adaptive-refinement approach can be seen in Fig. 11.
What may be surprising at first glance is, that although the rolling refinement approach requires multiple sub-graphs to be built, the overall time needed for generating the models is generally lower than in the static approach. This can be explained by the selected rolling horizon parameters, where H and H match. Hence, every

Fig. 10
Running times for optimization in initial iteration segment of the trajectory is essentially only considered once for refinement. Thus, some longer arcs that are in the sub-graph of the static approach never have to be considered in the rolling approach under these conditions. Of course, if H < H , parts of the trajectory are considered multiple times for the computation of the refinement nodes, leading to a rise in overall computation times for the rolling approach.
Multiple refinement iterations: Figure 12 shows the computation times for three iteration rounds and six aircraft. The rolling horizon length and horizon step equal the data given in Table 6.
For the initial step, the model generation is the time spent to build the whole graph and model. For the refinement rounds, both the time for model generation and solution of the optimization model of each step in the rolling horizon approach are summed. We do not depict the average graph generation or optimization time for the refinement steps in the rolling horizon approach, as these times vary considerably. This is depending on the length of the trajectory segments actually considered for refinement.
It is evident,that on average the most time is spent in the model generation steps of the refinement iterations. This is significantly slower than generating the original models. This can be explained by the restrictions to the base-network and the time expanded nodes of the sub-graphs,restrictions to the base-network Fig. 11 Running times for one refinement iteration with static and rolling approach Table 6 Instance specification leading to running-times of Fig. 10. Iteration one refers to the initial problem, while iteration two and three are the first, resp. second, refinement iteration and the time expanded nodes of the sub-graphs, which enforces more checks on whether arcs can be added or not. Note: As the sub-graphs for each aircraft are computed individually, there is significant potential for speeding up this part of the computations by parallelization. However, this is not the focus of our study and therefore left to potential users of the algorithm according their specific context.

Punctuality
In addition to running time, we shall also examine the punctuality increase for instances. Table 7 gives an overview of the above described instances with two to ten aircraft, with ten instances each. We give the minimal and maximal individual delay decrease for an aircraft over all instances as well as the minimal, maximal, and mean total delay decrease for instances in seconds. It is to be noted that negative  individual delay decreases exist, implying an additional delay in the favor of some further overall objective decrease. Note that a maximal individual delay decrease higher than t as given in Table 6 is possible as we are considering the rolling horizon approach. Hence, for every rolling step we are allowed the maximal deviation t . This constitutes a further advantage of the rolling approach over a full trajectory refinement.
In Fig. 13 trajectories are represented in x and time for an instance of ten aircraft. After one full refinement iteration in the rolling adaptive-refinement approach, we see that most aircraft arrive earlier than in the initial solution, and even the arrival sequence is changed for the turquoise and violet trajectory. In Fig. 14 we see the trajectories of the same instance projected on the plane. It is apparent, that some aircraft follow the same route when projected to the plane. We would like to highlight the detours of the green and one dark blue trajectory: These are considerably more severe after the initial iteration and can be flattened out through the rolling adaptiverefinement approach. Table 8 provides a comparison of the delay reductions with respect to the initial optimal solution for different numbers of refinement iterations on the six aircraft instances. Included are: One refinement iteration of the static refinement algorithm and the rolling refinement approach with one and two refinement iterations. The specification of the neighborhoods follows the data displayed in Table 6. Presented are the same key-values as in Table 7.
As expected, the total delay decrease is the highest for two rolling refinement iterations, while it is the lowest for the static refinement. The previously described effect regarding the potentially higher delay reduction through the Fig. 13 Trajectories over time for ten aircraft. Left: after the initial iteration; right: after one full refinement in the rolling adaptive approach algorithm design of the rolling approach becomes even more noticeable in Table 8. In fact, for none of the considered instances the total delay decrease of the static refinement was as high as of the rolling refinement after one full iteration.

Sequential refinement of dimensions
The results presented so far considered a simultaneous refinement of space and time to highlight the general properties of the presented approach. Another question that arises for applications is, whether a sequential refinement of dimensions can be beneficial. We thus compare the joint refinement as described in Table 6, iterations 1  and 2, with two sequential approaches as given in Table 9. The first scheme represent a time-first, space-second refinement, while the second represents space-first, time-second refinement. We tested these iteration schemes on the six-aircraft instances: The overall running times for the joint-discretization are between 1034 and 1977 seconds. For the time-first scheme they are between 628 and 1060 seconds. For the space-first scheme the total running times are between 696 and 999 seconds. This general decrease can be explained by a significantly smaller amount of nodes and arcs in the refinement sub-graphs. Even though three full iterations have to be carried out instead of two, the time for model generation is lower. An interesting effect can be observed for the obtained costs: In comparison to the joint discretization, for the time-first scheme there is a mean increase of 5.3% for the costs for deviations from the scheduled time.The fuel costs, however, decrease by 2.0% (mean value). In contrast, for the space-first refinement, the deviation costs are decreased by 9.4% and the fuel costs increased by 6.8% (mean values). Hence, depending on the desired application, consideration of different iteration schemes and evaluation of effects on the resulting solutions can be beneficial.

Conclusion and outlook
In this work we develop a new model for disjoint trajectories optimization for multiple commodities. The integer programming model is based on a graph and is able to consider heterogeneous commodities. It can be solved to global optimality with respect to the underlying discretization by integer programming techniques. Further, we present a discretization technique for three-dimensional space and time that allows trajectory optimization in a full dimensional environment and respecting arbitrarily shaped, potentially non-convex, restricted areas. To deal with huge instances that can result from considering sufficiently fine discetizations of space and time, a rolling horizon approach is integrated into an adaptive-refinement algorithm. In addition to the model presented, the developed algorithmic framework is evaluated for the real world application of approach planning and runway scheduling at an airport under the free-flight paradigm. An exhaustive computational study illustrates the properties of the approach and proves successful applicability of the methodology.
Interesting research questions remain for future work, including investigations to further improve the algorithm. Here, especially gaining and exploiting application specific insights to further reduce instance sizes in the refinement approach could be rewarding. E.g., are there general rules on resolution sequences, or can the set of refinement vertices that are actually relevant be characterized. Another research direction would be to consider prerequisites to provide optimality guarantees of the solutions obtained by the integrated rolling horizon and adaptive-refinement algorithm in comparison to the global optimal solution when considering the whole possible space and time with the same discretization. Uncertainties in arc costs, release dates or traversal times should also be considered to enable robust optimization of trajectories and even widen the field of applicability.