key: cord-0481266-pvy3026q
authors: Dias, Fernando H. C.; Rey, David
title: Aircraft Conflict Resolution with Trajectory Recovery Using Mixed-Integer Programming
date: 2022-03-22
journal: nan
DOI: nan
sha: 225186ea6d633f49d903cba410ae6ed40f431091
doc_id: 481266
cord_uid: pvy3026q

To guarantee the safety of flight operations, decision-support systems for air traffic control must be able to improve the usage of airspace capacity and handle increasing demand. In this study, we address the aircraft conflict avoidance and trajectory recovery problem. The problem of finding least deviation conflict-free aircraft trajectories that guarantee the return to a target waypoint is highly complex due to the nature of the nonlinear trajectories that are sought. We present a two-stage iterative algorithm that first solves initial conflict by manipulating their speed and heading control and then identifying each aircraft's optimal time to recover its trajectory towards their nominal. The avoidance stage extends existing mixed-integer programming formulations, and for the recovery stage, we propose a novel mixed-integer formulation. We assume that speed and heading control are continuous variables for this approach while the recovery time is treated as a discrete variable. In this approach, it is shown that the trajectory recovery costs can be anticipated by inducing avoidance trajectories with higher deviation, therefore obtaining earlier recovery time within few iterations. Numerical results on benchmark conflict resolution problems show that this approach can solve instances with up to 30 aircraft within 10 minutes.

*Corresponding author(s). E-mail(s):

fernando.cunhadias@helsinki.fi; Contributing authors: david.rey@skema.edu;

To guarantee the safety of flight operations, decision-support systems for air traffic control must be able to improve the usage of airspace capacity and handle increasing demand. In this study, we address the aircraft conflict avoidance and trajectory recovery problem. The problem of finding least deviation conflict-free aircraft trajectories that guarantee the return to a target waypoint is highly complex due to the nature of the nonlinear trajectories that are sought. We present a two-stage iterative algorithm that first solves initial conflict by manipulating their speed and heading control and then identifying each aircraft's optimal time to recover its trajectory towards their nominal. The avoidance stage extends existing mixed-integer programming formulations, and for the recovery stage, we propose a novel mixed-integer formulation. We assume that speed and heading control are continuous variables for this approach while the recovery time is treated as a discrete variable. In this approach, it is shown that the trajectory recovery costs can be anticipated by inducing avoidance trajectories with higher deviation, therefore obtaining earlier recovery time within few iterations. Numerical results on benchmark conflict resolution problems show that this approach can solve instances with up to 30 aircraft within 10 minutes.

Airspace usage has seen an increase in demand throughout the last decades (except during specific periods such as post 09/11 attacks and the COVID-19 pandemic). Aligned with this trend, the limited airspace has put the current air traffic management (ATM) system under intense pressure. This is reflected in the increasing amount of control necessary to guarantee safety, which is a crucial aspect of air traffic control (ATC). With the advance of unmanned aerial systems and urban air mobility, denser and congested traffic configurations are also expected. This configuration may lead to impairment of aircraft safety. Nevertheless, state-of-the-art methods for aircraft traffic control are reaching their limits, and new approaches, including more automation, have recently received significant attention in the field [1, 2] . Introducing automation within ATC systems has the potential to reduce controller workload and improve airspace capacity [3] . Conflict detection and resolution (CDR) is vital for the workload model of air traffic controllers. Hence, this calls for advanced conflict algorithms capable of acting as efficient decision-support tools.

The aircraft conflict avoidance and resolution problem can be formulated as an optimization problem in which the goal is to find conflict-free trajectories while minimizing a cost function, e.g. the deviation to the original flight path. Many strategies have been proposed to address this problem based on the type of manoeuvres issued to aircraft: speed, heading or altitude control. These strategies can be applied separately or in combination. In this study, we refer to conflict avoidance all conflict resolution approaches that do not consider the aircraft trajectory recovery problem. Conflict resolution using global optimization has recently received growing attention due to its ability to provide optimal solutions that consider all traffic within an airspace region. Following, we review the literature on conflict avoidance in Section 1.1 before focusing on efforts to develop approaches for trajectory recovery in Section 1. 2 We outline our contributions relative to the field in Section 1.1.

One of the first global optimization approaches for aircraft conflict resolution was introduced by [4] which proposed two formulations: one focusing on speed control and another focusing on heading control, and both minimize overall flight time. In the proposed MIP (Mixed Integer Problem) formulation for conflict resolution with speed control, the authors derived linear pairwise aircraft separation constraints based on the geometric construction introduced by [5] . These separation conditions are obtained by projecting the shadow of an intruder aircraft onto the trajectory of a reference aircraft. [6] was the first to observe that this geometric construction provided a basis to characterize the set of aircraft pairwise conflict-free trajectories via linear half-planes in the relative velocity (speed and heading) plane. The authors introduced a nonconvex formulation for the conflict resolution problem with speed and heading control and proposed a convex relaxation based on semi-definite programming and a heuristic algorithm to find feasible solutions on problems with up to 10 aircraft.

Subsequent approaches proposed speed control and altitude levelassignment to minimize fuel consumption by metering aircraft at conflict points [2] . In [7] , the authors proposed a two-stage stochastic optimization model accounting for wind uncertainty and using speed control. Multiobjective optimization formulations attempting to balance flight deviation with the total number of manoeuvres (velocity, heading or altitude change), building on the work of [4] were proposed by [8, 9] . Subliminal speed control methods which focus on speed control only for conflict resolution has also proven to be efficient and with low impact in terms of deviation and fuel consumption, although they may fail to resolve all conflicts [3, 10] .

More recently, nonlinear global optimization approaches received increasing attention in the literature. [11] proposed a hybrid algorithm that uses the optimal solution of a MILP (Mixed Integer Linear Problem) as the starting point for solving a nonlinear formulation of the same problem. [10] proposed an MINLP (Mixed Integer Non-Linear Problem) approach for conflict resolution with speed control only which highlights that subliminal speed control alone may not be sufficient to resolve all conflicts in dense traffic scenarios. Using a similar framework, [12] presented a two-step approach where a maximum number of conflicts are first solved using speed control only, and outstanding conflicts are solved by heading control. [13] proposed a formulation based on bi-level optimization with multiple follower problems, each of which represents a two-aircraft separation problem. The authors presented two formulations, one using speed control and another using heading control. A cut generation algorithm is proposed to solve the corresponding bi-level optimization problems. More recently, [14] proposed a complex number formulation for speed and heading control without any form of discretization that was further expanded by [15] in which an exact constraint generation algorithm was proposed. A systematic review of mathematical programming methods in ATC can be found in [16] .

In mathematical programming approaches, aircraft trajectory is usually divided into two stages. Those stages usually are called avoidance (or action) which encapsulates the manoeuvres necessary to avoid any conflict, while the trajectory recovery corresponds to the manoeuvres necessary to restore the aircraft to its original trajectory. Different approaches (such as genetic programming and heuristics) might be considered trajectory recovery as part of conflict resolution. However, this is usually separated in mathematical programming due to its non-trivial aspects, bringing many non-convexity and non-linearity.

Despite their potential effectiveness, most efforts in conflict resolution have focused on ensuring conflict avoidance but have overlooked the costs and mechanisms for modelling aircraft's recovery to their original trajectory. This may be critical when conflict resolution is performed using heading control which may significantly deviate the aircraft from their initial trajectory, thus possibly increasing flight operating costs. Aircraft conflict resolution with trajectory recovery is the problem of finding conflict-free trajectories which ensure that aircraft recover their initial flight path upon completion of the manoeuvres performed. This problem has received very little attention in the literature due to its challenging nature. Meta-heuristics such as genetic algorithms exploited the fact that formulating an analytical trajectory is quite challenging (due to trigonometric functions) and expressed the aircraft path by discretizing it. In [17] , the aircraft path is divided into different segments based on possible conflict points and turning points, thus predefining a set of eligible trajectories. An ant colony algorithm was also proposed by [1] where the authors determined the target point for each aircraft trajectory and discretized its trajectory using a specific timestep and in [18] , where another formulation is attempted with discretized space and time. In order to generate a smooth aircraft path and derive conflict-free solutions, genetic algorithms were implemented. [19] proposed a model which uses an analogy with light propagation theory to create conflict-free aircraft trajectories with recovery. [20] proposed a B-splines model which uses way-points of a given trajectory to design conflictfree trajectories with recovery. In [21] , the author proposed a formulation providing parallel trajectory recovery while minimizing fuel consumption and delays. In this model, aircraft are assumed to perform a preventive manoeuvre before their intersection with other trajectories, and the formulation aims to find trajectories that are parallel to the aircraft initial's trajectory. Heading angles are discretized, and the optimization controls both aircraft heading and recovery time. Recently, [22] proposed a manoeuvre-discretized model in which predefined sets of manoeuvres are available for aircraft, and a cliquebased formulation is proposed to find the optimal combination of conflict-free manoeuvres.

This literature review highlights that despite recent improvements in computational optimization, there remain significant open challenges in designing scalable global optimization approaches for conflict resolution in air traffic control, especially on incorporating recovery in a scalable and effective way. Only a few methods for aircraft trajectory recovery have been proposed, and they typically assume a simplified trajectory design (as expressed beforehand). This can be explained by the complexity and nonlinear aspects present in the trajectory recovery formulation, especially using mathematical programming. Approaches that have jointly addressed the avoidance and recovery problems often assume the existence of complete manoeuvre sets before optimization or any other forms of variable discretization or simplification.

As highlighted in our literature review, most of the state-of-art models using mathematical programming focus on conflict avoidance only. Specifically, they combine speed and heading control manoeuvres in deterministic approaches. In this study, we present a new two-stage algorithm for aircraft conflict resolution with trajectory recovery. In this approach, the speed and heading of aircraft are first optimized to avoid conflicts while minimizing the deviation from their initial trajectories. Then, in a second stage, aircraft trajectories are modified to recover a target position on the aircraft's initial trajectories. These two stages are incorporated into an iterative algorithm, where the recovery cost is projected into the avoidance stage. Hence, this algorithm aims to obtain a nontrivial solution where the overall cost of aircraft manoeuvres is minimized. The main contributions of this study relative to the literature are i) the extension of state-of-the-art mixed-integer formulations for conflict avoidance to include on/off variables and constraints to select aircraft that are being controlled in the conflict resolution process; ii) the design of a penalty-based algorithm that iterates the conflict-avoidance and trajectory recovery stages to find non-trivial solutions to the problem at hand, and iii) extensive numerical experiments on benchmarking problems from the literature that demonstrate the benefits of the proposed approach.

The paper is organized as follows: in Section 2, the two-stage algorithm is detailed, starting with the conflict avoidance formulation and then the trajectory recovery portion. In Section 3, the experimental framework, the results and discussions are presented. In Section 4, the conclusion and future research directions are provided.

In this section, we present a novel approach for conflict resolution with trajectory recovery that is based on decomposing the problem into two stages. The proposed approach aims to repeatedly solve these two sub-problems, i.e. conflict avoidance and trajectory recovery using a penalty-based mechanism to anticipate recovery costs in the avoidance stage.

In aircraft trajectory formulation, we assume that aircraft current and target positions are known, and that aircraft are initially not in conflict. This sets the context of the optimization problem of interest: given a set of aircraft with known current and target positions, find least-deviating conflict-free trajectories for all aircraft, such that aircraft may safely reach their target destination. To address this problem, we propose decomposing the trajectory Aircraft conflict resolution with trajectory recovery using mixed-integer programming optimization problem into two stages: 1) conflict avoidance and 2) trajectory recovery. The first stage focuses on controlling aircraft heading and speed to avoid all conflicts, while the second stage focuses on calculating the optimal time for aircraft to start safely recovering towards their target position. We focus on the two-dimensional conflict resolution problem and only consider horizontal aircraft manoeuvres for brevity. The extension to the vertical case can be addressed by incorporating flight level change manoeuvres in the conflict avoidance stage [15] and ensuring safe recovery to aircraft target flight level and position. We leave this extension for future research.

The avoidance stage aims to determine the optimal variation in speed and angle for each aircraft to avoid conflicts. In the recovery stage, it is necessary to calculate and identify the manoeuvres required to recover aircraft initial trajectories. In order to impose airspace safety, it is necessary to impose the separation conditions in the recovery stage as well. In mathematical programming, recovery approaches are challenging for several reasons. Assuming that aircraft trajectories can be modelled as linear and that uniform motion laws hold, the conflict separation conditions are created via Euclidean distance. This leads to equations using quadratic and trigonometric elements. This creates difficulty in organizing and developing an optimization model using those equations. Thus, such formulations may not scale easily and be used in larger scenarios. Second, Euclidean distance creates quadratic components that are also nonlinear and nonconvex concerning decision variables and cannot be easily solved. As presented by [4] , [8] and [14] nonlinear constraints can be further simplified into a set of integer-linear components with regards to aircraft velocity. However, these formulations only focus on separating aircraft to solve existing conflicts and ignore the process of planning aircraft recovery to a target destination. As described Fig 1, trajectories such as the one highlighted in the blue and red segments are an example of conflict avoidance and trajectory recovery, respectively. At first, the aircraft manoeuvre as the blue segment represents. Another manoeuvre is performed to redirect the aircraft towards its target point (in yellow). In this study, we focus on developing a global optimization approach to construct piecewise linear conflict-free trajectories (avoidance and recovery) for a set of aircraft to minimize the total deviation from their initial original trajectories.

Consider a set of aircraft A sharing the same flight level. For each aircraft i ∈ A, assuming uniform motion laws apply, its position is:

which v i is the speed, x i and y i are the initial coordinates of i at the time of optimization, θ i is its initial heading angle, θ i is its deviation angle and q i is the speed deviation. The relative velocity vector of i and j, denoted v ij , can be expressed as: where:

Incorporating these elements into the equation of motion gives:

The relative position of aircraft i and j at time t can be represented as p ij (t) = p i (t) − p j (t). Let d be the horizontal separation norm, typically d = 5 NM. Two aircraft i, j ∈ A are horizontally separated if and only if: ||p ij (t)||≥ d, for all t ≥ 0.

Let P be the set by the pairs of aircraft, i.e. P = {i ∈ A, j ∈ A, i < j}. For each pair (i, j) ∈ P, the relative position vector p ij (t) is:

and the relative velocity

Imposing the separation condition, gives for each pair (i, j) ∈ P:

Let x ij = x i − x j and y ij = y i − y j . Squaring both sides in Eq. (7), we obtain:

The time instant t min ij represents the time of minimal separation of aircraft i and j. As noted in several studies [10, 12, 14] , if t min ij ≤ 0 then aircraft i and j are diverging and, assuming aircraft are separated at t = 0, they are thus separating for any t ≥ 0. Further, substituting t min ij in f ij (t) yields:

then aircraft i and j are separated. Writing the terms g ij and t min ij as functions of the relative velocity variables v x ij and v y ij , we obtain the following disjunctive pairwise aircraft separation conditions:

The separation condition (11) can be further linearized following the approach described by [14] and [15] . By alternatively fixing variables v x ij and v y ij and solving the resulting quadratic equations, we can obtain the solution for g ij (v x ij , v y ij ) = 0. By isolating each variable, we obtain the discriminants:

Assuming aircraft are initially separated, then x 2 ij + y 2 ij − d 2 ≥ 0 holds and thus the discriminants are positive, and the roots of equation g(v x ij , v y ij ) = 0 are the lines defined by the system of equations:

If all coefficients in Eq. (13) are non-zero, they define two lines, denoted

can be characterized based on the position of (v x ij , v y ij ) relative to these lines (see Figure ( 2a)). Recall that according to Eq. (9), the sign of the dot productp ij · v ij indicates aircraft convergence or divergence. Let (P ) be the equation of the line corresponding to the dot productp

The line defined by (P ) splits the plane {(v x ij , v y ij ) ∈ R 2 } in two half-planes, each of which representing converging and diverging trajectories, respectively. Consider the line normal to (P ), denoted (N ):

≥ 0 corresponds to a pair of conflict-free trajectories. Hence, the conflict-free region is nonconvex. According to the speed profile, we can divide the solution space into two regions: the region hashed in red corresponds to the conflict region C while the complement of this region represents conflict-free trajectories. We recall the formal definition of the conflict region. Definition 1 (Conflict region, [15] ) Consider a pair of aircraft (i, j) ∈ P. Let C be the subset of R 2 defined as:

C is the conflict region of (i, j) ∈ P.

The conflict region of pair of aircraft represents the set of relative velocity vectors (v x ij , v y ij ) which corresponds to conflicts. As depicted in Figure 2b , the feasible region is nonconvex, which can be solved via a disjunctive formulation. However, as presented in [14] and [15] , the set of linear equations described by

The inner box with black lines corresponds to the velocity bounds B in the deterministic scenario. The region is hashed in red corresponds to the conflict region C. If relative velocity B intersects with the conflict region B, then there exists a risk of conflict. 

The red lines represent the lines P and N . The dashed blue lines correspond to the linear equations R 1 and R 2 that are the roots of

is shown by the + and -pink symbols.

Eqs. 15 is equivalent to (10) as detailed by Theorem 1 in [15] . In each convex sub-region, the lines delineate the conflict-free region. The expressions of these lines depend on initial aircraft positions, i.e. x ij , y ij . Integer-linear separation conditions with regards to aircraft velocity components can be derived as follows, and we model this disjunction using the variable z ij ∈ {0, 1} defined as:

In conflict avoidance problems, a common objective is to minimize the combined deviations of all aircraft. This may lead several aircraft to perform minimal conflict avoidance manoeuvres, which may not be desirable from an operational perspective. As observed in [23] , even small deviations may result in costly recovery manoeuvres. In order to generate trajectories where fewer aircraft are controlled and those that are controlled have a reduced total cost, an additional binary variable is introduced. This variable determines whether an aircraft is controlled or not. In this formulation, the avoidance stage has as the objective to minimize the total deviation and the number of aircraft that are controlled. Let f i ∈ {0, 1} for i ∈ A represent the control variable: f i = 1 represents the case where aircraft i modifies its speed and/or heading; and f i = 0 represents the case where aircraft i does not perform any conflict avoidance manoeuvre.

Based on the binary variable f i , if aircraft i does not perform any manoeuvre, its speed and heading control should remain unchanged. This can be expressed as:

Considering this control variable, if a pair of aircraft is in conflict, at least one of these aircraft must perform manoeuvres. Let P 0 be the set the aircraft that are initially in conflict such as that

Recall that a pair (i, j) of aircraft is initially in conflict if and only if condition (11) is not satisfied. Therefore, for all pairs of aircraft that are initially in the conflict, the following cuts are valid inequalities:

(17) These cuts are added to strengthen the conflict resolution formulations by observing that at least one of them must perform an avoidance manoeuvre for any pair of aircraft initially in conflict.

To incorporate the impact of controlling aircraft in the objective function, we propose to minimize a weighted sum of two terms: a fixed cost linked to control variables (f i ) and a variable cost linked to trajectory deviation variables (q i , θ i ). We denote λ f the weight representing the fixed cost of controlling an aircraft, and we denote w in [0, 1] the weight used to capture the trade-off between heading and speed control deviations.

For each aircraft i ∈ A, we assume that the speed rate variable is lower bounded by q i and upper bounded by q i , i.e.:

We assume that the heading deviation is lower bounded by θ i and upper bounded by θ i , i.e.:

To derive lower and upper bounds on relative velocity components v x ij and v y ij , we re-arrange Eq. (6) using trigonometric identities:

Let v ij,x , v ij,x and v ij,y , v ij,y be the lower and upper bounds for v x ij and v y ij , respectively. These bounds can be determined using Eq. (21) and the bounds on speed and heading control provided in Eqs. (19) and (20) . The derived bounds on the relative velocity components can be used to define a box in the

Definition 2 (Relative velocity box) Consider a pair of aircraft (i, j) ∈ P. Let B be the subset of R 2 defined as

B is the relative velocity box of (i, j) ∈ P.

The relative velocity box B characterizes all possible trajectories for the pair (i, j) ∈ P based on the available 2D deconfliction resources, i.e. speed and heading controls. To characterize the set of conflict-free trajectories of a pair of aircraft (i, j) ∈ P, we compare the relative position of the relative velocity box B with the conflict region of this pair of aircraft.

The conflict avoidance stage can be formulated exactly as presented in [15] incorporating the penalty control variable f i as described in Eq. (16) and constraint (17) .

This formulation is nonconvex due to the velocity constraint, which is nonconvex quadratic (Eq. (6)). This results in a MINLP formulation which is challenging to solve and does not scale easily. Coefficients γ l ij , φ l ij and γ u ij , φ u ij (present in (15)) can be pre-processed based on the sign of x ij and y ij . Finally, B represents the bounds for the velocity variables (v x ij , v y ij ).

In this study, we are considering a simple trajectory recovery model in which deviated aircraft perform a second manoeuvre in order to change their trajectory towards their target point. Let (x i ,y i ) be the coordinate of the target point, those second manoeuvres are opposing to the deviations in the avoidance stage. The speed component q r is simply defined as the aircraft are returning their velocity profile to its nominal value, i.e. q r i = 1. For the heading changes, the deviation angle is based on the time each aircraft moves from its avoidance trajectory towards its target point. For a given aircraft i ∈ A, its recovery trajectory is defined as : p i (t) = [ x i + q ivi cos( θ i + θ i )t i +v i cos( θ i )t, y i + Aircraft conflict resolution with trajectory recovery using mixed-integer programming q ivi sin( θ i + θ i )t i +v i sin( θ i )t] T , where t i corresponds to the recovery time, i.e., the value of time in which each aircraft change from its avoidance trajectory to its trajectory recovery (see Figure 3) . Hence, the deviation angle θ r i can be calculated as:

where the d a i (t i ) (see Eq. (25)) corresponds to the distance flown during the conflict avoidance stage (from the initial position until the aircraft reach t i and changes to trajectory recovery) and d r i (t i ) (see Eq. (26) to the distance flown during the recovery stage (from the time t i where aircraft starts its trajectory until it reaches its destination point (x i ,y i )) (see Figure 3 ). If the aircraft did not change its heading angle, the recovery angle is not calculated. In the avoidance stage, the core idea is to compare the distance between each pair of aircraft and derive time-independent expressions that can be escalated to all instances and applied simultaneously. However, the same cannot be achieved when avoidance and recovery are calculated simultaneously. This is because the aircraft recovery trajectories are a function of the time when aircraft perform their recovery manoeuvre (t i ). It is evident that this expression is nonlinear concerning speed and heading variables. This solidifies that, by using mathematical programming, it is quite complex to reproduce the avoidance model to solve the recovery stage.

Given this context, it is expected that some simplification level is required to make any formulation using mathematical programming viable. One of these strategies were implemented by [23] . In their paper, a naive approach for trajectory recovery was implemented, and it is composed of a two-stage algorithm where avoidance and recovery are solved sequentially. In this formulation, the manoeuvres determined during the action stage are compensated in the recovery stage. In order to incorporate the cost of avoidance in the recovery stage, the total deviation on avoidance is passed on to the recovery stage. It determines how costly the deviation is in terms of speed and angle applied. In this way, the trajectory recovery can compensate for the cost during the avoidance stage. Another characteristic of this formulation is that the recovery time is a discrete variable. This reduces the feasible region and makes the options for trajectory recovery a limited, finite set. Based on those characteristics, this formulation can solve small to more significant instances (up to 30 aircraft) in a reasonable amount of time. More details can be found in [23] . An example of the solution obtained via this approach can be seen in Figure 4 .

The main drawback of this formulation is the lack of anticipation of the recovery costs at the avoidance stage. As both stages are solved separately, and each stage is solved only once, the solution obtained after solving both stages is myopic. This is not a significant issue for small instances due to the small amount of aircraft, but it becomes concerning in medium to large instances. As the number of aircraft and conflicts increase, the naive approach of [15] is expected to lead to weak solutions that do not anticipate recovery costs in the design of aircraft trajectories.

In order to explore the situation where the avoidance takes into account the recovery costs, we are proposing an alternative version of such an algorithm where the total cost is optimized throughout an iterative algorithm by projecting the cost of recovery operations into the avoidance stage to obtain non-trivial solutions. In the following section, we explain this algorithm in detail.

Each aircraft needs to perform opposing manoeuvres for trajectory recovery to cancel the deviation applied during avoidance. Similar to the avoidance model, our goal is to guarantee that all pairs of aircraft are separated throughout the recovery stage. Since the separation condition in Eq. (15) is based on linear motion, we need to distinguish the trajectory stage of each aircraft i ∈ A, i.e. before and after its recovery time t i . We denote A i the avoidance trajectory of aircraft i and R i its recovery trajectory. Given a pair (i, j) of aircraft, we need to ensure that aircraft are separated during all pairwise trajectory stages, denoted

Observe that separation for the stage A i A j is already ensured by the solution of Model 1. If aircraft i and j were to recover at the same time period, then aircraft will transition from A i A j to R i R j directly. Otherwise, if i (resp. j) recovers before j (resp. i),

To avoid trigonometric forms and to obtain a tractable formulation, we discretize aircraft recovery times. Let T be the set of time periods available for recovery, we require: .

where is the length of time periods. Abusing notation, we redefine the separation condition expressed in Eq. (11) as: g ij (m, n) ≥ 0 and t min ij (m, n) ≤ 0 where the pair (m, n) indicates the time period indices of recovery times t i and t j , respectively. Let Ω XiXj be the set of conflict-free pairs of recovery times for aircraft i, j ∈ A where X i represents the state of the trajectory of aircraft i, i.e. A i or R i ; and X j represents the state of the trajectory of aircraft j, i.e. A j or R j . This set can be specified into three different sets corresponding to the three different states during the recovery stage. The set Ω RiRj is defined as:

For the states A i R j and R i A j an extra condition is required. Consider the state A i R j : if the lines of motion corresponding to trajectories A i and R j are in conflict but aircraft i turns into recovery prior to the start of this conflict, then no conflict will occur. This illustrated in Figure 5 where g AiRj < 0 and t AiRj > 0. Let τ AiRj (t j ) be the smallest root of g AiRj = 0 if j recovers at time t j . If aircraft i recovers prior to τ AiRj (t j ), i.e. t i ≤ τ AiRj (t j ), then the conflict will be avoided. Accordingly, we define:

Let ρ im be a binary variable equal to 1 if aircraft i ∈ A recovers at time period m ∈ T and 0 otherwise. To track the states of aircraft pair (i, j) which are activated, we introduce two binary variables α ij and β ij . Those variables are used to identify whether t i < t j (α ij = 1) which activates state R i A j , or if t i > t j (β ij = 1) which activates state A i R j . Variables α ij and β ij are defined Aircraft conflict resolution with trajectory recovery using mixed-integer programming via the constraints:

We use the following constraints to exclude conflicting trajectories from the solution. Observe that states A i R j and R i A j are conditional on the recovery times of t i and t j and thus the corresponding constraints are only active if i and j do not recover at the time period.

Aircraft are assigned a recovery time via the constraint:

The second stage aims to identify the optimal time for aircraft to recover towards their target position. The avoidance cost function of aircraft i, denoted a i , is defined as:

For the recovery stage, we are proposing to minimize the total weighted recovery time, which accounts for manoeuvres applied in the avoidance stage as well as aircraft recovery times. i∈A m∈T

this expression comes from the approximated area created by the deviation during the avoidance stage, and it is quadratic in t 2 i . The trajectory recovery formulation is summarised in Model 2 which is a MIQP (Mixed Integer Quadratic Programming).

minimise i∈A m∈T

The main idea behind the proposed penalty-based approach is to capture the recovery cost during the avoidance stage. Therefore, its goal is to influence the behaviour of the avoidance stage to anticipate the cost of trajectory recovery and attempt to construct an efficient trajectory across both stages. This can be achieved by using the solution of the recovery stage as a preemptive cost in the avoidance stage. In this case, the algorithm aims to find a trade-off between the deviation costs incurred during the avoidance and recovery stages. Let T C i be the total cost per aircraft defined as the combined cost of avoidance and recovery.

where λ t is the weight for the recovery time component. Therefore, the objective function in the avoidance stage is modified to account for the anticipated cost of trajectory recovery. Let r i the recovery cost of each aircraft i ∈ A, where the avoidance cost is based on the deviation in the avoidance stage (q i and θ i being the optimal values for speed changes and heading angles).

By adding this expression into the avoidance stage, we aim to account for the impact of recovery costs within the avoidance stage. This approach aims to penalize the control of aircraft in the avoidance stage proportionally to their anticipated total trajectory deviation. Thus, the objective function in the avoidance stage can be rewritten as:

The proposed iterative approach is summarized in Algorithm (1). In each iteration, the avoidance stage (as described in Model (1)) is solved using the updated objective function described in Eq. (37).

Based on the solutions obtained from solving Model (1), the sets Ω AiRj , Ω AjRi and Ω RiRj are pre-calculated. The recovery stage is solved subsequently and the total cost and the cost variation ∆T C (based on Eqs. (35) and (38)) are calculated. The stop criteria is based on the overall total cost between two consecutive iterations, i.e.:

where n corresponds to the number of the iteration. If the value of ∆T C is below a predefined threshold, the algorithm converges and stops; otherwise, the recovery cost is updated and the algorithms proceeds. In the first iteration, we assume r i = 0; for the consecutive iterations, the anticipated recovery cost r i is calculated using the solution obtained in the previous iteration.

In order to solve Model (1), which is nonconvex and has nonlinear constraints that are notoriously challenging, we use the complex number formulation described by [14] and an adaptation of the solution method proposed by [15] .

The experimental framework used to test the proposed mixed-integer formulation for the trajectory recovery using continuous heading is introduced in Section 3.1. Then a detailed analysis of four groups of instances is presented in Section 3.2. The computational performance of the proposed approaches is thoroughly explored in Section 3.4, respectively.

We test the performance of the proposed mixed-integer formulations and algorithm using two benchmark problems from the literature: the Circle Problem (CP) and the Random Circle Problem (RCP). The two types of benchmarking instances are illustrated in Figure 6 . The CP consists of a set of aircraft uniformly positioned on the circle heading towards its centre. Aircraft speeds are assumed to be identical; hence the problem is highly symmetric (see Figure  6a ). The CP is notoriously difficult due to the geometry of initial aircraft configuration and has been widely used for benchmarking CD&R algorithms in the literature [1, 10, 12, 14, 24] . CP and RCP instances are named CP-N and RCP-N, respectively, where N is the total number of aircraft. For reproducibility purposes, all formulations and data used for testing are made available online at the public repository https://github.com/acrp-lib/acrp-lib.

Require: A,θ,v, q, q, θ, θ,  Ensure: q , θ , t i , Add cuts in Eq. (17) 14:

t ← Solve Stage 2 -Recovery using Model (2) 15:

if Infeasible then 16: return Infeasible 17: else 18: t ← t 19: q ← q 20: θ ← θ 21: t ← t 22: Calculate ∆T C 23: end if 24: if ∆T C ≤ then 

For all tests, a speed regulation range based on the subliminal speed control [−6%, +3%] is used, and the typical heading control range [−30 • , +30 • ] is also used. For the preference weight, we use as w = 0.5 and λ f = 1 in the objective in Eq. (33). This value was selected such that both heading and speed control terms were of a comparable order of magnitude, emphasizing penalizing heading control. For stage 2, a total of |T |= 15 time periods are used, with a step of = 2 minutes. To solve avoidance, the solution method used was an adaption of the algorithm implemented by [15] Disjunctive, and for the recovery stage, the model is described in 2. Both are solved with Cplex Python API and a time limit of 5 minutes per solving and 15 minutes per instance. For the first iteration, ρ i is assumed as 1, and it is calculated after the recovery in the following iterations. The proposed approach is compared to the algorithms presented by [23] . Those methods are respectively referred to as ExactNaive and GreedyNaive Recovery. The first method corresponds to a similar formulation of the method presented in this paper: it is the first iteration of Algorithm 1 where the avoidance and recovery stages are each solved only once. The second is a heuristic-based trajectory recovery procedure that iterates over all time steps and uses a priority list to decide which aircraft can be recovered at each time step. The priority list used is based on r i values (36). The algorithm first sorts aircraft accordingly and iterates overtime periods. At each time, the algorithm iterates over the sorted list of aircraft and check if each aircraft can be recovered at the current time. The process is repeated until no aircraft can recover at the current time. The proposed algorithm has a worst-case time complexity of O(|T ||A| 3 ). Details of our implementation of both benchmarking algorithms can be found in [23] and example of their solution can be found in Fig. (7) . Those models were slightly modified to be comparable to the algorithm presented in this paper. Specifically, the control variable f i is introduced in the avoidance stage. Since the solution methods presented in [23] only solve each stage a single time, they are comparable to the first iteration of Algorithm (1).

To illustrate the proposed two-stage iterative algorithm, the optimal solution obtained is plotted for CP instances with five aircraft and RCP instances with 10, 20 and 30 aircraft. In Figure (7) , dashed grey lines represent initial aircraft trajectories, green lines represent the avoidance trajectory of stage 1, and orange lines represent recovery trajectories of stage 2 using Model 2. In contrast, green lines represent the avoidance trajectory of stage 1, and the orange line represents the trajectory of stage 2 using Algorithm 1.

This behaviour is even more accentuated in RCP-10-1, RCP2-20-2 and RCP-30-2. In the first, it is shown that only three aircraft are necessary to be altered with relatively minor deviations, and that would be sufficient for guaranteeing separation conditions instead of deviating 9 out of 10 aircraft. The solutions obtained via the continuous angle causes the slightest disturbance in the network, which is a favourable trade for those formulations. For RCP-20 and RCP-30, there are more aircraft to be altered, but overall, they show the same behaviour: more aircraft with higher avoidance cost and reduced recovery costs, with overall cost improved.

To quantify the impact of the preference weights λ t and λ f in the proposed total cost function Eq. (35), we conduct numerical experiments on one instance of each of the two types of benchmarking instances for varying values of those parameters individually. This experiment focuses on the typical heading control range [−30 • , +30 • ]. Our goal is to show that by varying those preference weight in λ t ∈ ]0, 1[ and λ f ∈ ]0, 1[ the decision-maker can control the desired level of trade-off between the recovery time and total deviation. Recall that in the total cost function (37), λ f is the coefficient of f i which is minimal for f i = 0; therefore, one can expect that increasing λ f will tend to penalize the number of aircraft controlled and therefore decrease the recovery time. On the other hand, increasing the value of λ t will tend to penalize the recovery time and, therefore, lead to potentially a greater total deviation of the aircraft.

This behaviour is confirmed in our numerical experiments. Specifically, we solve the instances CP-4 and RCP-10-3 for λ t , λ f = 0.1, . . . , 0.9 in steps of size 0.1, i.e. for a total of 9 values per instance. Each parameter is changed separately. All instances are solved using Algorithm 1. The change in the total deviation Σ d and the total recovery time Σ t are summarized in Figures (8) and (9) . In Figure (8) , it shows the behaviour of CP-4 and RCP-10-3 when λ f is changed: using the proposed objective function, the decision-maker can control which manoeuvre is prioritized by scaling up or down the preference weights λ f accordingly. Higher values of λ f will minimize the total deviation, while smaller values of λ f will prioritize recovery time. We use λ f = 0.25 in the numerical experiments presented in the remaining of the paper. We find that increasing λ f monotonically decreases the total deviation and monotonically increases the recovery time for all four instances tested. For λ t (in Figure (9) ), higher values lead to the minimization of recovery time and bigger deviation while the increasing leads to smaller deviation in heading angle and speed while larger recovery time. We use λ t = 0.25 in the numerical experiments presented in the remaining of the paper. For the parameter w, we are using the value w = 0.5 based on the sensitivity analysis presented in [15] . 

The performance of Algorithm 1 is reported in Table 1 for CP instances with 4 to 15 aircraft. The performance of the two-stage iterative algorithm is examined. The following tables represent the results for CP and RCP instances. (1 − q i ) 2 + θ 2 i (in red) and Σ t represents the total recovery time defined as i∈A t 2 i (in green). All instances tested here are CP-4 and RCP-10.

In the header, there are four sections; the first is based on the instances containing |A| is the number of aircraft, and n c is the number of conflicts. The second group is based on the avoidance stage using the manoeuvre control variable: Obj. is the objective function, ∆T C is the optimality gap, Time(s) is the runtime in seconds, followed by total deviation in terms of |1 − q i |, |θ i | and f i , respectively. The third group contains the results from the recovery stage Table 1 summarises the results for CP instances and Table 2 for the RCP instances with 10,20 and 30 aircraft.

In the results, it is expected from those instances that all the aircraft are moving towards to centre. In terms of the objective function value the objective function is largely composed of i∈A f i . For those instances, it is consistently observed that one aircraft can remain in its initial condition in each instance. This differs considerably from the results using naive avoidance, where the global solution would imply that all aircraft needs to perform some manoeuvres. In instances with ten or more aircraft, there are only slight deviations. Similarly, the variation in the heading angle is considered negligible for smaller instances. However, it does not increase proportionally with the number of conflicts, revealing that this manoeuvre is used sporadically. Moreover, the current model results show that only a small amount of aircraft is required to be controlled, but the total deviation is more considerable. In terms of runtime, all instances with 12 aircraft can be solved within the time limit. In the recovery stage, the average recovery time does not increase with the number of aircraft (around 0.32 hours on average), with the minimum recovery time around 5 minutes while the maximum peaks at almost 1 hour. However, all the values of instances with 12 aircraft or more result in time out solutions. For recovery time, most instances can be solved within the time limit, but for the more prominent instances, the solution obtained in the avoidance stage timeout and the solution in the recovery stage. The algorithm only requires up to 3 iterations to obtain a solution within the convergence criterion used in terms of iterations. In comparison with the benchmark in Table 3 , we can conclude by the T C value that the speed variation is quite smaller and the angle deviation is overall larger. This means that the recovery time is also considerably larger. This is mainly due to the trivial solutions where the main focus is to solve the conflict as it comes without any consideration relating to the recovery (as it is done in the iterative algorithm). Another observation is that the recovery strategies presented in the ExactNaive and in the GreedyNaive methods always provided a strategy where the recovery time was slightly higher for the latter and considerably higher for the former. This solidifies that using the cost of recovery as a preemptive measure contributes to better performance and usage of the airs space configuration.

In Table 2 , the results for RCP instances are presented. For the avoidance stage, the objective function reflects the number of aircraft required to be separated. As indicated by f i , 2.16 for RCP-10, 6.66 for RCP-20 and 12.6 for RCP-30, which reflect the value obtained by However, compared to the heading deviation obtained using the complex number formulation, it is clear that heading deviation is considerably smaller and Aircraft conflict resolution with trajectory recovery using mixed-integer programming Table 3 : Summary of Comparison between the Penalty-Based Algorithm and the Exact-Naive and Greedy-Naive presented in [23] using the T C as comparison key for the CP Instances. Table 4 : Summary of Comparison between the Penalty-Based Algorithm and the Exact-Naive and Greedy-Naive presented in [23] using the T C as comparison key for the RCP Instances

ExactNaive GreedyNaive |A| i∈A T C i 1 |A| i∈A T C i Time i∈A T C i 1 |A| i∈A T C i Time i∈A T C i 1 |A| i∈A T C iTime

Exact Naive Greedy Naive Aircraft conflict resolution with trajectory recovery using mixed-integer programming that most of the total deviation is caused by speed changes. The opposite behaviour is observed here. The runtime for those instances is reasonably short, and instances with up to 30 aircraft can be solved in less than 10 s. Compared with the complex number formulation results, the Disjunctive presented gives solutions balanced with equal contribution from speed and heading control. These current results state that speed and heading control are kept to their minimum while most of the control is set around whether certain aircraft need to be altered. Finally, the optimality gap is negligible in all instances.

In the recovery stage, the runtime increases considerably with the number of aircraft given that the number of alternative routes to solve increases drastically: up to 2.14 s for RCP-10, 259 s for RCP-20 and time out for 100% of the instances under 5 minutes of the time limit. In terms of iteration, RCP-10 instances and RCP-20 instances require up to 4 iterations to achieve convergence while solving RCP-30 instances; only up to two iterations are executed due to the time limit. In comparison with the benchmark methods, the table 4, the results showed that via analyzing the T C values for the benchmarks, the recovery time in the naive recovery is significantly larger than the values obtained in the algorithm of this paper and in the greedy recovery is slightly bigger but only by a small margin. On the other hand, the runtime observed by the greedy recovery is relatively small, which highlights the computational trade-off at stake.

The findings are summarised in Section 4.1 and future research directions are discussed in Section 4.2.

A new mixed-integer formulation and a penalty-based algorithm for ACRP with trajectory recovery are proposed. The performance of this approach revealed that by echoing the cost associated with trajectory recovery, the avoidance could be manipulated to minimize pre-emptively the overall cost of ACRP. On average, in a couple of iterations, it showed that most instances could force the aircraft to have more significant deviations in the avoidance stage but with an earlier recovery time in the recovery stage. In this case, a trade-off between avoidance and recovery is observed. Another advantage of this addition is the concept of stability of the solutions throughout iterations. One of the main issues that happen by having an iterative algorithm is the variation in the profile of the solution. In each iteration, the algorithm provides a different group of aircraft as a possible candidate to recover, and each solution has a different deviation cost, which can be an improvement. As expected, the deviation cost is higher with a higher deviation angle, but the number of aircraft manoeuvring is small. In this situation, adding the f i variable introduces an extra hurdle in modifying the status of an aircraft. With multiple iterations, it can be expected that less deviation in the whole set of aircraft will be observed, although a more significant deviation in individual aircraft. This creates stability in the solution. In the point of view of air traffic controllers, solutions where there are many aircraft manoeuvring lead to higher workload and are not desirable or implementable in real applications. The comparison between the naive and the iterative approaches shows that the latter can outperform the former in a set of criteria. First, the naive approach shows that the total cost is relatively more significant because there is no compensation for the recovery costs while the avoidance is processed. In this way, the naive approach overcompensates the avoidance manoeuvre and ignores less trivial solutions where small avoidance could have been taken. On the other hand, using the iterative approach starts with a solution similar to the naive approach and iteratively reaches towards more balanced and nontrivial solutions where we observed that a smaller avoidance angle is obtained; hence a later recovery time also be applied. Second, the naive approach does not incorporate the concept of manoeuvre control. There is a lack of control of the minimum number of aircraft trajectories in that scenario. The naive approach follows the trivial optimization problem, whereas a flaw in the original formulation favours a symmetrical solution where all aircraft are required to perform some manoeuvre with a near similar deviation. Therefore, all aircraft in such instances will alter their trajectories. In the iterative approach, this assumption is optimized where only a subset of aircraft change their speed and angle. This solution is a bit more realistic, given that in actual airspace occasions, fewer aircraft will be impacted by eventual conflict. Another aspect of this decision choice is that the aircraft required to alter their trajectory has significant manoeuvres instead of the solution using the naive approach where most solutions are pretty small and even hard to be detected. Finally, the iterative approach has more time leverages as continuous variables, allowing for a broader range of possible trajectories as the avoidance angle. The main setback of the iterative approach is that the convergence of such a model is imposed via the reduction of the total cost, which might lead to many iterations if the convergence rate is low. Also, the optimality condition is not guaranteed in both approaches, although the iterative is the straightforward best approach given all the reasons explained throughout this paper.

Alternatively, an adaptation of the algorithm in [15] is proposed to solve the avoidance stage considering speed, heading control and manoeuvre control as decision variables. This is incorporated into an iterative two-stage algorithm that incorporates the projected cost of recovering the aircraft into the avoidance stage. A relatively small number of iterations showed that it could increase the overall deviation in the avoidance pre-emptively to reduce the recovery cost and, ultimately, the total cost. In the numerical experiments, the performance of the proposed algorithms shows that the number of aircraft required to be altered is reduced. Therefore, the deviation is more significant in heading angles for those required to be changed. However, the recovery time is shorter compared to the complex number formulation. In comparison with the Disjunctive results, it also shows that while the former provides solutions that are more balanced in terms of manoeuvres, they are also more invasive because it affects the whole set of aircraft. The algorithm present in this study presents an alternative solution where fewer aircraft are controlled will higher deviations, suggesting that this behaviour improves the total cost.

The most significant limitation in the iterative two-stage algorithm is conflict avoidance and trajectory recovery decomposition. Although this is a common practice in mathematical programming, there is no guarantee of the quality of the solution. Ideally, a global solution should be created by jointly optimizing both stages in a unified formulation. This is heavily challenged by the nonlinearity and complexity of such formulations. Therefore, further research is needed to model this problem more simply and efficiently. In addition, the cost of recovery is projected into the avoidance stage throughout the iterations. Alternatively, stochastic optimization methods can be used to determine the expected cost of recovery manoeuvres already incorporated in the avoidance stage. Finally, the discretization of any variable is a limitation and modelling recovery time as a continuous variable may help reduce the total cost of trajectories.

Further testing is required to fully assess the impact of the control variables on solution quality. Notably, quantifying the impact of first stage formulations with non-discretized heading changes and second stage formulation continuously is critical to evaluate the cost of manoeuvre discretization. Both formulations presented here are deterministic methods where no uncertainty is considered. Future research will focus on introducing uncertainty with the proposed formulation to incorporate trajectory prediction uncertainty capable of accounting for the expected cost of recovery at the conflict avoidance stage and develop robust conflict resolution algorithms to account for the uncertainty of aircraft trajectory prediction when planning.

Ant colony optimization for air traffic conflict resolution

Near real-time fuel-optimal en route conflict resolution

Subliminal speed control in air traffic management: Optimization and simulation

Conflict resolution problems for air traffic management systems solved with mixed integer programming

A geometric optimization approach to aircraft conflict resolution

Resolution of conflicts involving many aircraft via semidefinite programming

A two-stage stochastic optimization model for air traffic conflict resolution under wind uncertainty

Collision avoidance in air traffic management: A mixed-integer linear optimization approach

Exact and approximate solving of the aircraft collision resolution problem via turn changes

Maximizing the number of conflict-free aircraft using mixed-integer nonlinear programming

Hybridization of nonlinear and mixed-integer linear programming for aircraft separation with trajectory recovery

Mixed-integer nonlinear programming for aircraft conflict avoidance by sequentially applying velocity and heading angle changes

Detecting and solving aircraft conflicts using bilevel programming

Complex number formulation and convex relaxations for aircraft conflict resolution

Disjunctive linear separation conditions Aircraft conflict resolution with trajectory recovery using mixed-integer programming and mixed-integer formulations for aircraft conflict resolution

Airspace conflict resolution: A unifying mathematical framework and review

Optimal resolution of en route conflicts

Automatic aircraft conflict resolution using genetic algorithms

A lightpropagation model for aircraft trajectory planning

Solving air traffic conflict problems via local continuous optimization

A space-discretized mixed-integer linear model for air-conflict resolution with speed and heading maneuvers

Two decomposition algorithms for solving a minimum weight maximum clique model for the air conflict resolution problem

A two-stage algorithm for aircraft conflict resolution with trajectory recovery

Equity-oriented aircraft collision avoidance model