Thibault Lehouillier

Moncef Ilies Nasri

Jérémy Omer

Guy Desaulniers

Solving the Air Conflict Resolution Problem under Uncertainty as an Iterative Bi-Objective Mixed Integer Linear Program *

Keywords: Air Traffic Management, Air Traffic Control, Mixed Integer Linear Programming, Uncertainties, Multi-Objective Programming

In this paper, we tackle the aircraft conflict resolution problem under uncertainties. We consider errors due to the wind effect, the imprecision on the aircraft speed prediction, and the delay in the execution of maneuvers. Using a geometrical approach, we derive an analytical expression for the minimum distance between aircraft, along with the corresponding probability of conflict. These expressions are incorporated into an existing deterministic model for conflict resolution. This model solves the problem as a maximum clique of minimum weight in a graph whose vertices represent possible maneuvers and where edges link conflict-free maneuvers of different aircraft. We then present a solution procedure focusing on two criteria, namely fuel efficiency and the probability of re-issuing maneuvers in the future: we iteratively generate solutions of the Pareto front to provide the controller with a set of possible solutions where he/she can choose the one corresponding the most to his/her preferences. Intensive Monte-Carlo simulations validate the expressions derived for the minimum distance and the probability of conflict. Computational results highlight that up to 10 different solutions for instances involving up to 35 aircraft are generated within three minutes.

Introduction

Automating air traffic control

In the current air traffic management (ATM) organization, the air traffic control (ATC) is in charge of maintaining safety. To this end, controllers monitor the traffic to ensure the separation between all aircraft at all times. A projected loss of separation between two aircraft is called a conflict and must be solved by the controller. To this end, avoidance maneuvers are issued to the pilots of the involved aircraft to prevent the loss of separation. Maintaining safety in the

The air conflict resolution problem

One complex and central problem encountered in ATC is the air conflict resolution problem (CR). A conflict occurs when two aircraft are too close to each other regarding predefined horizontal and vertical separation distances of respectively 5NM and 1000ft. To solve a conflict, the controllers issue maneuvers that can consist of speed, heading or altitude changes. Given the current position, speed, acceleration and the predicted trajectory of a set of aircraft, the CR problem corresponds to identifying the conflict-free maneuvers that minimize a given cost function. The CR problem can be tackled following two different settings, namely deterministic and stochastic. The first one assumes that aircraft follow exact trajectory predictions, along with maneuvers applied without any errors. However, uncertainties are one of the reasons why ATC is a complicated task. The weather conditions, along with the incomplete knowledge of the physical characteristics of the aircraft and the imprecision during the communication and maneuver execution processes represent the main factors of uncertainty in ATC [START_REF] Erzberger | Conflict detection and resolution in the presence of prediction error[END_REF]). In this context, the uncertainties cause a perturbation of the trajectory, inducing cross and along-track errors in the prediction of the trajectory. The along-track error (or longitudinal error) is the distance between the predicted aircraft position and the projection of the actual aircraft position on the predicted trajectory. The cross-track error (or lateral error) corresponds to the distance between the actual aircraft position and the predicted trajectory. Figure 1 illustrates these errors. They can jeopardize the conflict resolution process. To tackle this issue, stochastic resolution methods aim at solving the CR problem while taking into account these perturbations.

Literature review on the CR problem

The CR problem is one of the most widely studied problems in ATM. We provide a synthetic analysis of the studies that were most influential to our work, both in a deterministic and a stochastic setting.

A complete coverage of the existing literature on the deterministic CR problem may be found in the review performed in [START_REF] Martín-Campo | The collision avoidance problem: methods and algorithms[END_REF]. Mixed integer linear and nonlinear programming are powerful theoretical frameworks for the study of CR. With the realistic restriction that the aircraft perform at most one maneuver at the initial time, [START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF] exploit the geometry of the separation constraints to develop two mixed integer linear programs (MILPs) that allow either speed changes with constant headings or heading changes with constant speeds. Alonso-Ayuso et al. ( 2012) extend the model of [START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF] by introducing continuous instead of instantaneous speed changes. More recently, [START_REF] Omer | A space-discretized mixed-integer linear model for air-conflict resolution with speed and heading maneuvers[END_REF] develops a MILP with a space discretization using only the points of interest for the conflict resolution.

Uncertainties can be gathered and modeled as having a global impact on the trajectory prediction. [START_REF] Ballin | An analysis of landing rates and separations at the Dallas/Fort Worth International Airport[END_REF] quantify the along-track error by comparing prediction and actual data for the Dallas Fort Worth Airport. Results highlight that for a time horizon shorter than 20 minutes, the error follows a normal distribution. [START_REF] Irvine | A geometrical approach to conflict probability estimation[END_REF] develops an expression of the minimum distance and the corresponding probability of conflict using a geometrical approach. The author models cumulative cross and along track errors that are affected to each aircraft at the beginning of the observation. After applying this initial perturbation, aircraft are assumed to evolve in a deterministic environment. Uncertainties can also be divided into different categories than can be modeled more specifically. For instance, [START_REF] Lygeros | Aircraft and weather models for probabilistic collision avoidance in air traffic control[END_REF] model the effect of the wind and the resulting FMS correction. [START_REF] Cole | Assessment of the Rapid Update Cycle (RUC) with near real-time aircraft reports[END_REF] and [START_REF] Schwartz | Accuracy of ruc-1 and ruc-2 wind and aircraft trajectory forecasts by comparison with acars observations[END_REF] conduct statistical studies comparing predictions to actual data in the Denver area to derive the correlation structure of the wind. [START_REF] Chaloulos | Effect of wind correlation on aircraft conflict probability[END_REF] study the perturbations due to imprecisions in the speed and air temperature measures. The authors model the error as a normal distribution.

When the uncertainties become too complex to derive exact probability expressions, Monte-Carlo simulations are often performed. [START_REF] Prandini | A probabilistic approach to aircraft conflict detection[END_REF] use Monte-Carlo simulations to develop a model where the wind correlates the cross and along track errors.

Critical analysis and contribution statement

The literature review highlights that a lot of results have been established to solve the CR problem, both in the deterministic and the stochastic setting. Nevertheless, some features still need to be addressed. More specifically, we formulate three observations that we consider important when designing a resolution tool for the CR problem. The first one relies on the fact that robustness is critical in ATC. A large span of factors can have a dramatic impact on the conflict resolution. As a consequence, it is necessary to provide the controller with a tool as robust as possible. In other words, the controller needs to be ready to handle every possible situation. To this end, the mathematical framework in the developed decision tools needs to remain valid, whatever the hypotheses followed. Unfortunately, a lot of models lack of consistency when it comes to the modification of hypotheses, like the introduction of uncertainties, or other modeling features concerning the aircraft dynamics. For instance, the constraints in [START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF] are linear when aircraft perform either a heading change or a speed change, but become nonlinear when both are performed. The second observation is related to the multi-objective nature of the CR problem. Indeed, focusing on only one objective, like the fuel consumption, or the delays, does not necessarily reflect all the aspects of the problem, nor does it respect the users' preferences. Several multi-objective approaches of the CR problem have been performed [START_REF] Menon | Optimal strategies for free-flight air traffic conflict resolution[END_REF]; [START_REF] Tomlin | Conflict resolution for air traffic management: A study in multiagent hybrid systems[END_REF]; Alonso-Ayuso et al. ( 2016)), and research needs to be conducted in this direction. The last observation we formulate is that in conflict resolution the notion of optimality is subjective. Indeed, depending on the objective to optimize, the optimal solution is not necessarily far better than other good solutions. As a consequence, providing the controller with only one solution can be restrictive, depending on the context and the controller's preferences. Few work has been done on methods generating a set of solutions instead of a single solution. For instance, satisfycing game theory (see [START_REF] Stirling | Satisficing games[END_REF] for a description of the theory) allows to generate a set of satisfycing solutions regarding two criteria representing the preferences of the players in terms of efficiency or resource consumption. Applications to ATC have been considered (see, e.g., [START_REF] Archibald | A satisficing approach to aircraft conflict resolution[END_REF]), but the hypotheses are quite restrictive, and the model suffers from a lack of computational power.

Our main contributions in our effort to provide an answer to the aforementioned remarks are twofold. First, we provide an analytical expression of the minimum distance and the probability of conflict in a context allowing complex uncertainties: the error in wind predictions is considered, along with the error on the aircraft speed prediction. We introduce the uncertainty on the delay in the execution of maneuvers, which to our knowledge has not been studied yet in the literature, although it is a reality in ATC. With this approach, we are able to cover a large span of uncertainties involved in ATC. Besides, these computations are fast compared to a simulation-based approach that can be more time consuming. Second, we model the CR problem as a bi-objective problem minimizing fuel consumption and the probability that the controller has to reissue maneuvers: we sequentially solve a mono-objective MILP. With this approach, we benefit from the powerful results yielded by MILPs, namely the guarantee of finding an optimal solution (if existing) in a short time, even for large and complex instances. Each iteration generates a solution that is immediately available to the user. The set of generated solutions is a tight approximation of the Pareto front of the solution. This method allows the user to choose which solution to apply within the generated set, depending on his/her preferences or other factors. The MILP used is taken from a preliminary study performed by Lehouillier et al. (2015b,a). It was chosen because it fully separates the modeling of the aircraft dynamics, maneuvers and cost function from the resolution process. As a consequence, the hypotheses considered do not jeopardize the validity of the proposed mathematical framework, and in particular the introduction of uncertainties. Besides, the fact that we are able to introduce uncertainties in the model from Lehouillier et al. (2015b) validates its robustness.

To evaluate the model, we first validate the computations derived for the probability of conflict by running Monte-Carlo simulations. We use several test beds generating 2000 random scenarios to verify the correctness of the developed theory. After the validation of the computations, we test our iterative resolution procedure by conducting intensive simulations on a benchmark of structured and random instances that are complex to solve. The aim of the experiments is to verify that our algorithm is able to provide the user with a set of solutions in a short period of time, while ensuring that separation is maintained in complex situations.

The organization of the paper will be as follows. We formulate the problem in Section 2. We describe the mathematical model to be adapted in Section 3. We detail the iterative optimization procedure used to generate the set of solutions in Section 4. The method is then tested and analyzed through intensive experiments described in Section 5.

Problem Formulation

Aircraft dynamics

As in the majority of the literature, we use a three-dimensional point-mass model for aircraft dynamics. This model establishes relationships between the different physical parameters of each aircraft.

dp x dt = V cos γ cos χ (1) dp y dt = V cos γ sin χ (2) dp z dt = V sin γ (3) dγ dt = g 0 V (n cos φ -cos γ) (4) dχ dt = g 0 V n sin φ cos γ ( 5
)
dV dt = F T -F D m -g 0 sin γ (6)
The position of the aircraft is given by the coordinates (p x , p y , p z ) of its center of gravity in a local coordinate system, (p x ,p y ) being its coordinates in a horizontal plane and p z its altitude. The aircraft flies at speed V and the angles χ, φ and γ correspond respectively to its heading, roll and pitch. F T and F D denote the norm of the thrust and drag forces respectively, m is the aircraft mass, n is the load factor and g 0 corresponds to the gravitational acceleration.

In this article, we make the assumption that aircraft are stabilized and follow a planar motion in a single flight level. Aircraft follow their trajectory with a stepwise constant acceleration. This assumption is realistic since it respects the time-continuity of speed, and it corresponds to a setting where maneuvers are performed smoothly.

Aircraft maneuvers

The maneuvers are horizontal maneuvers consisting in heading and speed changes. These maneuvers are performed dynamically in order to avoid a significant error in separation distance. Aircraft execute a speed or a heading change with a constant acceleration and turn angle, respectively, according to values extracted from [START_REF] Paielli | Modeling maneuver dynamics in air traffic conflict resolution[END_REF]. Other types of maneuvers, i.e., flight level changes, could be considered without changing the validity of the mathematical resolution.

Aircraft trajectory recovery

We consider that aircraft follow a 4D contractual trajectory, which represents a compromise between the user's preferences and the capacity constraints of the network. The trajectories of the aircraft then have to meet time and space requirements over a sequence of 4D points. Noncompliance with this contract induces penalty fees to companies. As a consequence, it is important to make sure that, after resolving every conflict, every aircraft recovers its initial 4D trajectory. Ensuring a strict velocity control can be very costly and almost impossible in practice. Physicial recovery is required, whereas time recovery is optional, but it is favored by giving a penalty on the time shift between the 4D contract and the 4D trajectory after the maneuvers are performed.

Maneuver cost

The cost of a maneuver corresponds to the additional burnt to perform the maneuver, along with a time shift penalty. This measure serves as an indicator of the perturbation of the 4D trajectory induced by the executed maneuvers.

For a jet commercial aircraft f with constant altitude, the fuel consumption by time and distance unit is given by ( 7) and ( 8):

C t, f t, V f (t) = c 1, f 1 + V f (t) c 2, f F T, f (t) (7) C d, f t, V f (t) = C t, f t, V f (t) V f (t) (8)
where c 1, f and c 2, f are numerical constants depending on the type of aircraft f that are extracted from the BADA performance tables EUROCONTROL (2011).

The time shift penalty is computed according to the method found in [START_REF] Omer | Hybridization of nonlinear and mixed-integer linear programming for aircraft separation with trajectory recovery[END_REF]. The penalty corresponds to the extra fuel burnt to make up for the time shift.

Modeling the uncertainties

In this subsection, we detail the models used to describe the different uncertainties.

Error on wind prediction.

The aircraft are considered as flying within a wind field. Control commands are issued to reach the desired airspeed v a , while the control units monitoring the aircraft speed are ground-based. As a consequence, the groundspeed v g can be linked to the airspeed. Let w(p, t) denote the windspeed at position p at time t. We have that:

v g (t) = v a (t) + w(p(t), t) (9)
The wind vector is decomposed in a nominal part corresponding to weather forecasts, and a random part describing the difference between the actual wind and its nominal part. The impact of the nominal wind of the aircraft dynamics is quite complex and was briefly studied in the literature. Most publications focus on the random part of the wind, and do not consider the nominal part. In this paper, we focus solely on the random wind.

The wind field is a set of random vectors W(p, t) depending on the time and the point of space considered. Taking the wind into account complexifies the conflict resolution. Indeed, aircraft that are close from each other undergo highly correlated winds that will impact the conflict resolution. In this case, the error of prediction for the different aircraft become correlated. We follow the models presented in [START_REF] Lymperopoulos | Sequential Monte Carlo methods in air traffic management[END_REF]. The authors simplify the computations performed in [START_REF] Cole | Assessment of the Rapid Update Cycle (RUC) with near real-time aircraft reports[END_REF] and [START_REF] Schwartz | Accuracy of ruc-1 and ruc-2 wind and aircraft trajectory forecasts by comparison with acars observations[END_REF] in order to save execution time. The wind is stationary and isotropic, and each random vector W(p, t) follows a zero-mean normal distribution such that the following conditions hold:

E[W(p 1 , t 1 )] = 0, ∀t 1 ∈ R + , ∀p 1 ∈ R 2 (10) E [ W(p 1 , t 1 )|W(p 2 , t 2 ) ] = 2 f (t 1 , p 1 , t 2 , p 2 ), ∀(t 1 , t 2 ) ∈ R 2 + , ∀(p 1 , p 2 ) ∈ R 4 ( 11
)
where f is the correlation function associated with the random wind developed in [START_REF] Cole | Assessment of the Rapid Update Cycle (RUC) with near real-time aircraft reports[END_REF].

We assume that the flight management system (FMS)compensates for the lateral errors, but does not correct the along-track errors. Indeed, the majority of commercial aircraft are equipped with 3D FMS which track only the cross-track errors.

Error on aircraft speed measures.

We consider the uncertainties due to the imprecision of speed and air temperature measures presented in [START_REF] Chaloulos | Effect of wind correlation on aircraft conflict probability[END_REF]. These errors have an impact on the along-track speed of the aircraft which is modeled as a zero-mean normal variable independent from the other aircraft. Since these two uncertainties are highly time-correlated, the authors assumed they were constant over time.

Delays in the execution of maneuvers.

We model uncertainties induced by delays in the execution of maneuvers, which to our knowledge has not been studied yet.

In the literature, models always assume that the performance of the maneuvers is instantaneous. However, there are several actions required before the maneuver can actually be performed. First, the automated decision tool has to provide the controller with a feasible solution. Then, the controller has to process the solution and then communicate the corresponding instructions to the pilots, before they can execute the maneuvers.

More formally, let T i denote the maneuver delay for aircraft i. T i is decomposed as follows:

• the time required for the resolution tool to provide the controller with a feasible solution, denoted T s ;

• the time during which the controller analyses the solution and communicates it to the different aircraft, denoted T c ;

• the time required for the pilot of aircraft i to execute the communicated maneuver, denoted

T p i .
In other words T i is the sum of a term shared by all aircraft including the solution process and the controller's communication, and a term depending on the pilot of i. Figure 2 summarizes the whole process resulting in the delay. In this paragraph, we detail the work presented by [START_REF] Irvine | A geometrical approach to conflict probability estimation[END_REF] which serves as the foundation for the method we use to derive the expression of the probability of conflict. In his article, Irvine models the global impact of the uncertainties and the resulting cross and along-track errors, instead of modeling each source of error differently. Let A i and A j be two aircraft flying at a stabilized altitude at speed v i and v j , respectively. Their trajectories intersect in O with a crossing angle θ ij . Let x i (t) and x j (t) denote the curvilinear abscissa at time t of A i and A j in a coordinate system centered on O. The distance between the two aircraft can be computed as follows:

d(t) 2 = x i (t) 2 + x j (t) 2 -2x i (t)x j (t) cos θ ij (12)
If d(t) is replaced by the separation distance required between A i and A j , denoted d sep , Equation ( 12) defines an ellipse in the coordinate system (O, x i , x j ). The aircraft follow straight trajectories at constant speed, hence the set of points x i (t), x j (t) t≥0 defines a straight line of slope

dx j dx i = dx j d t dx i d t = v j v i = m ( 13
)
where m is the speed ratio between the two aircraft. If this line intersects the ellipse, then the aircraft are said to be in conflict. Figure 3 illustrates this condition.

Figure 3: Ellipse in the coordinate system (O, x i , x j )

Tangent of equation

x j = mx i + d sep λ Tangent of equation x j = mx i -d sep λ x i x j O Conflict-free trajectory

Conflicting trajectory

To derive an analytical expression of this condition, the author uses the two tangents of the ellipse that are parallel to the parametric line x i (t), x j (t) t≥0 . Their equations are given as follows:

x j = mx i ± d sep λ where λ = sin θ ij m 2 -2m cos θ ij + 1
The minimum distance between A i and A j in the deterministic case, denoted d min , can then be expressed as a function of the initial curvilinear abscissa of the two aircraft, x 0 i and x 0 j .

d min = |λ(x 0 j -mx 0 i )| (14)
A i and A j are in conflict if and only if the minimum distance d min is strictly less than the minimum separation distance allowed d sep :

-d sep < λ(x 0 j -mx 0 i ) < d sep
Irvine then considers along-track errors and makes the assumption that within the range of along-track distances for which conflict is possible, the along-track error is approximately constant and that the aircraft flies with its predicted speed. This assumption is used to quantify the cumulative along-track error between t = 0 and the instant where the two aircraft are the closest from each other in the deterministic case, denoted t = τ ij , given by

τ ij = (x 0 j + mx 0 i ) cos θ ij -(x 0 i + mx 0 j ) ||v i ||(1 -2m cos θ ij + m 2 ) (15)
This value is chosen for τ ij because the computation of the instant where the two aircraft are the closest from each other in the stochastic case is hard in practice. Consequently, deriving a handy formula of the probability of conflict would not be possible. Besides, considering τ ij is a realistic assumption, since for the time intervals considered, the difference due to the approximation would be negligible.

The cumulative along-track error, denoted ∆L(τ ij ), follows a normal distribution N (0, α σ τ ij ), where α σ is a constant. This error is applied to the initial position of the aircraft who then evolves in a entirely deterministic environment. This yields a new expression of the minimum distance in an uncertain setting, denoted D min .

D min = |λ(x 0 j + ∆L j (τ ij ) -mx 0 i -m∆L i (τ ij ))| ( 16
)
D min is the sum of a deterministic term with the sum of independent random variables. The sum of independent, normally distributed random variables is also normally distributed, with a mean equal to the sum of the means of the individual distributions, and variance equal to the sum of the variances of the individual distributions. As a consequence, we have that D min follows a normal distribution of mean µ d and variance σ 2 d where

µ d = λ(x 0 j -mx 0 i ) (17) σ 2 d = (α σ τ) 2 (1 + m) 2 (18)
Irvine applies a similar reasoning for the impact of cross-track errors, but since we assume that the FMS compensates for these errors, in our article, we do not give any details about it.

The probability of conflict P c is given by

P c = P |λ(x 0 j + ∆L j (τ) -mx 0 i -m∆L i (τ))| < d sep (19) = 1 σ d √ 2π d sep -µ d -d sep -µ d exp - u 2 2σ 2 d du (20) = Φ d sep -µ d σ d -Φ -d sep -µ d σ d ( 21
)
where Φ is the cumulative distribution function of the standard normal distribution.

Enriching the formula.

In this subsection, we modify the formula derived by Irvine by introducing the errors on the wind prediction, the speed measures, and the delay in the execution of maneuvers. These errors are independent. We note these cumulative errors ∆X i (τ ij ) and ∆X j (τ ij ), respectively. They can be decomposed as follows:

∆X i (τ ij ) = ∆W i (τ ij ) + ∆S i (τ ij ) + ∆D i (τ ij ) (22) ∆X j (τ ij ) = ∆W j (τ ij ) + ∆S j (τ ij ) + ∆D j (τ ij ) (23)
where ∆W . , ∆S . and ∆D . denote the cumulative error due to the wind, the speed prediction and the maneuver delay, respectively. Subsections 2.5.1 and 2.5.2 yield the following expressions for ∆W . and ∆S . :

∆W i (τ ij ) = W|v i ||v i || τ ij ∆S i (τ ij ) = Υ i τ ij (24) ∆W j (τ ij ) = W|v j ||v j || τ ij ∆S j (τ ij ) = Υ j τ ij ( 25
)
where Υ i and Υ j denote the error due to speed measures for A i and A j , respectively. The cumulative along-track error due to the maneuver delay is slightly more complex to determine. For the sake of clarity, we give an illustrative example in Figure 4 where two aircraft A i and A j flying with a speed v 0 i and v 0 j have to perform a heading change of value θ i and θ j , respectively. They perform these maneuvers with a delay corresponding to random variables denoted T i and T j , respectively. i and x0 j and the ones in the deterministic setting x 0 i and x 0 j . This difference, denoted ∆D 0 i and ∆D 0 j , is computed as follows:

∆D 0 i = T i ||v 0 i || -T i ||v 0 i || cos θ i - T j ||v 0 j || sin θ j sin θ ij - T i ||v 0 i || sin θ i cos θ ij sin θ ij ( 26
)
∆D 0 j = T j ||v 0 j || -T j ||v 0 j || cos θ j - T i ||v 0 i || sin θ i sin θ ij - T j ||v 0 j || sin θ j sin θ ij sin θ ij ( 27
)
If t < T i then A i has not started its maneuver yet and flies at speed v 0 i . If t ≥ T i then A i has flown during T i at speed v 0 i before changing its speed to v i . Since after T i , aircraft A i flies at speed v i like it is supposed to, the cumulative along-track error due to the delay T i until τ ij is in fact cumulated on the interval ]0, T i ]. The value of this error is derived by

∆D i (T i ) = T i (||v 0 i || -v i ) (28)
Variables ∆D i (τ ij ) and ∆D j (τ ij ) are derived with

∆D i (τ ij ) = ∆D i (T i ) + ∆D 0 i (29) ∆D j (τ ij ) = ∆D j (T j ) + ∆D 0 j (30)
yielding the following expressions:

∆D i (τ ij ) = 2T i ||v 0 i || -T i v i -T i ||v 0 i || cos θ i - T j ||v 0 j || sin θ j sin θ ij - T i ||v 0 i || sin θ i cos θ ij sin θ ij (31) ∆D j (τ ij ) = 2T j ||v 0 j || -T j v j -T j ||v 0 j || cos θ j - T i ||v 0 i || sin θ i sin θ ij - T j ||v 0 j || sin θ j cos θ ij sin θ ij (32)
To derive the new expression of the minimum distance between aircraft A i and A j , we aggregate the errors ∆X i (τ ij ) and ∆X j (τ ij ) into

D min = λ x 0 j + ∆X j (τ ij ) -m x 0 i + ∆X i (τ ij (33) = λ x 0 j + ∆W j (τ ij ) + ∆S j (τ ij ) + ∆D j (τ ij ) -m x 0 i + ∆W i (τ ij ) + ∆S i (τ ij ) + ∆D i (τ ij ) (34)
We simplify Equation (34) in order to determine an approximation of the distribution of variable D min .

In order to derive an analytical expression of the probability, we perform the computation with the hypothesis of a constant wind, which modifies the along-track speed of an aircraft flying at speed v by a factor w|v ||v|| . This assumption seems reasonable since the considered intervals of detection and resolution are quite small and the wind is highly time-correlated. The terms related to the wind then correspond to

∆W j (τ ij ) -m∆W i (τ ij ) = T j W|v 0 j ||v 0 j || + (τ ij -T j ) W|v j v j -m T i W|v 0 i ||v 0 i || + (τ ij -T i ) W|v i ||v i || (35)
We approximate T i and T j by the mean of their distribution µ T i and µ T j in the quadratic terms, in order to find the analytical expression of the probability of conflict. This approximation is acceptable, since in Section 5 we use a distribution of T i where the standard deviation is small compared to the mean. Equation ( 35) can be rewritten:

∆W j (τ ij ) -m∆W i (τ ij ) = W µ T j v 0 j ||v 0 j || + (τ ij -µ T j ) v j ||v j || -mµ T i v 0 i ||v 0 i || -m(τ ij -µ T i ) v i ||v i || (36) = W|u (37)
The terms involving the error on speed prediction can be simplified into

∆S j (τ ij ) -m∆S i (τ ij ) = Υ j τ ij -mΥ i τ ij (38) = τ ij (Υ j -mΥ i ) (39)
The terms related to the maneuver delay correspond to

∆D j (τ ij ) -m∆D i (τ ij ) = r j T j -mr i T i (40)
where

r j = ||v j || -||v 0 j || cos θ j + ||v 0 j || sin θ j tan θ ij + m||v 0 j || sin θ j sin θ ij r i = ||v i || -||v 0 i || cos θ i + ||v 0 i || sin θ i tan θ ij + ||v 0 i || sin θ i m sin θ ij
Equations ( 35), ( 39) and ( 40) yield a simplified expression for D min :

D min = λ(x 0 j -mx 0 i ) + λ (Υ j -mΥ i )τ ij + r j T j -mr i T i + W|u (41)
Equation ( 41) expresses D min as a deterministic term λ x 0 jmx 0 i , added with the sum of the following independent random variables

• (Υ j -mΥ i )τ ij ∼ N 0, (1 + m)σ Υ τ ij ; • r j T j ∼ N (r j µ T j , r j σ T j ) ; • -mr i T i ∼ N (-mr i µ T i , mr i σ T i ); • W|u ∼ N (0, σ W ||u||).
D min follows a normal distribution of mean µ D and variance σ 2 D given by

µ D = λ (x 0 j + r j µ T j ) -m(x 0 i + mr i µ T i ) (42) σ 2 D = λ 2 σ 2 Υ (1 + m) 2 τ 2 ij + (mr i σ T i ) 2 + (r j σ T j ) 2 + (σ W ||u||) 2 (43)
The probability of A i and A j being in conflict corresponds to the probability of the event |D min | < d sep :

P(|D min | < d sep ) = P(-d sep < D min < d sep ) = Φ d sep -µ D σ D -Φ -d sep -µ D σ D

Deterministic Model

In this section, we describe the resolution method developed by Lehouillier et al. (2015b,a) that will serve as a foundation for the optimization procedure presented in Section 4. The main idea is to model the CR problem as a maximum clique of minimum weight problem. To this end, we build a graph whose vertices represent maneuvers for the different aircraft, and where edges link conflict-free maneuvers of different aircraft. A maximum clique of minimum weight yields a conflict-free situation of minimal cost. The advantage of this process is that it fully separates the modeling of aircraft dynamics, the separation verification and the costs computations from the resolution: whatever the hypotheses considered, and in particular taking into account uncertainties, the proposed mathematical framework will remain valid. The remainder of this section highlights the key elements of modeling and resolution of the model.

Graph construction

In this subsection, we introduce the conflict graph G = (V, E ) used to model the CR problem.

The set of vertices is defined as V = 1; |M| , where M denotes the set of possible maneuvers for all aircraft. We denote V f the set of vertices corresponding to aircraft f . Let (i, j) ∈ V × V be a pair of vertices representing maneuvers (m i , m j ) ∈ M × M of aircraft ( f i , f j ) ∈ F × F . For i = j, we write m i m j when no conflict occurs if aircraft f i follows maneuver m i while aircraft f j performs maneuver m j . The set of edges E corresponds to the pairs of maneuvers performed by two different aircraft without creating conflicts: For a more synthetic presentation, we consider in this subsection that maneuvers and vertices are equivalent without loss of generality. As explained in the previous subsection, the cost of a maneuver depends on its execution, which itself varies with the maneuvers of the other aircraft. As a consequence, we need to define the cost of the edges before the cost of the vertices.

E = (i, j) ∈ V × V, i = j : m i m j (
Again, for ease of presentation, an edge e = (i, j) is considered as a pair of maneuvers. We compute the cost of an edge e = (i, j) as a pair constituted of the cost of maneuvers i and j, denoted C (i,j) i and C

(i,j) j . Let us consider a maneuver i. The cost of each edge linking i to another maneuver j corresponds to an execution time t j i which is the minimum time during which i and j have to be executed before a safe return can be performed by at least one of the corresponding aircraft.

To determine the cost of i, denoted c i , we need to compute the time t i during which it is actually applied. If i is not in the optimal solution, then t i = 0. Otherwise, t i is given by

t i = max j∈V ∩C t j i (45)
Equation ( 45) states that maneuver i has to be applied long enough in order to be conflict-free with every other chosen maneuver. As a consequence, we can determine c i :

c i =    max j∈V ∩C C (i,j) i if i ∈ C 0 otherwise

MILP formulation

In our model the costs of the vertices are not determined a priori, since they depend on which vertices are in the clique. As a consequence, the dedicated algorithms of existing graph theory libraries cannot be used in this study. To address this issue, we formulate our problem as a MILP that can be solved with any generic MILP solver.

The decision variables of the model all relate to the vertices of the graph. They correspond to the choice of the vertices in the clique and the cost of each vertex:

• x i = 1 if vertex i is part of the maximum clique 0 otherwise • c i ∈ R + is the cost of vertex i.
The clique search can then be modeled as the following MILP, denoted MIP:

minimize ∑ i∈V c i (46) subject to x i + x j ≤ 1, ∀(i, j) ∈ V × V \ E (47) ∑ i∈V x i = |F | (48) c i ≥ C (i,j) i (x i + x j -1), ∀(i, j) ∈ E (49) x i ∈ {0; 1}, ∀i ∈ V (50) c i ≥ 0, ∀i ∈ V (51)
The objective function ( 46) minimizes the cost of the maneuvers. Constraints ( 47) are clique constraints stating that two nonadjacent vertices must not be part of the clique. In terms of conflict resolution, it means that two maneuvers in conflict must not be part of the solution. Constraint (48) defines the cardinality of the maximum clique. Constraints ( 49) are used to compute the cost of the vertices: if a vertex is in the maximum clique, then its cost must be greater than the cost on every edge connecting it to other vertices in the clique. Otherwise, no particular constraint is imposed on the vertex cost. Constraints ( 50)-( 51) are binarity and nonnegativity constraints, respectively.

Inserting uncertainties into the deterministic model

In this subsection, we explain how the expression of the probability of conflict between two aircraft derived in Subsection 2.6 is used to modify the deterministic model presented in this section.

An edge exists between two maneuvers if they are conflict-free. In other words, if the probability of conflict associated with these maneuvers is 0. If they are in conflict (i.e if the probability of conflict was 1), then no edge is drawn between the corresponding vertices. To take into account the uncertainties, we change the necessary condition to build an edge.

The set of edges E is defined by

E = {(i, j) ∈ V × V, i = j : P c (i, j) < δ s } (52)
where δ s is a security threshold restricting the set of possible maneuvers. δ s represents an upper bound on the probability that a conflict remains after the maneuvers are issued. We remind here that a remaining conflict will always be solved: the controller will issue another set of maneuvers. In other words, δ s can be regarded as an upper bound on the probability that the controller uses a recourse to solve the problem once again. This adaptation of the deterministic setting makes a good pairing with the expression computed in Subsection 2.6. Indeed, when the set of possible maneuvers becomes very large, the number of probabilities to compute would require a huge computational effort if they were determined through simulation, whereas with our approach, we determine these values instantaneously.

Bi-Objective Optimization Procedure

In this section, we detail the bi-objective approach designed to solve the problem. This method optimizes the CR problem according to two criteria depicting the efficiency and the probability of having recourse related to a solution. It iteratively solves the model MIP according to the first criterion, while imposing a certain improvement on the second criterion between two consecutive resolutions. Each resolution results in a solution approximating the Pareto front of the problem. In the end, the method provides the air traffic controller with a set of solutions corresponding to different trade-offs between efficiency and probability of using a recourse.

Optimization Criteria

The first criterion of optimization corresponds to the objective function of the mathematical program MIP presented in Section 3:

z e = ∑ i∈V c i ( 53
)
This objective is the total amount of additional fuel burnt induced by the chosen maneuvers, and represents the aspect of a solution related to its economical efficiency. Indeed, it is an indicator on the perturbation of the planned trajectories and gives an insight into the effort required to catch up with the initial flight plan after the maneuver is performed. In addition to the perturbation of the set of aircraft itself, this criterion also illustrates the perturbation on surrounding traffic.

The second criterion is given by

z s = ∑ i∈V ∑ j∈V j =i P ij x i x j (54)
where P ij is the probability of conflict of maneuvers i and j, and x i and x j are the decision variables corresponding to whether or not maneuvers i and j are chosen. Value z s corresponds to the expected number of conflicts potentially remaining after the solution is applied. This is a relevant measure for the controller as it gives an idea of the potential additional effort required to solve the problem once again in the close future. The higher z s is, the higher the probability of re-issuing avoidance maneuvers will be. Variable z s represents an indicator of the additional workload and cognitive charge that will potentially be required in order to definitely solve the problem.

To keep the expression of z s linear in the decision variables, we apply Fortet's linearization [START_REF] Fortet | Applications de lalgebre de Boole en recherche opérationelle[END_REF]. We introduce a new set of binary variables y ij respecting the following constraints:

y ij ≤ x i , ∀i ∈ V, ∀j ∈ V (55) y ij ≤ x j , ∀i ∈ V, ∀j ∈ V (56)
y ij ≥ x i + x j -1, ∀i ∈ V, ∀j ∈ V (57) (58)
yielding a new expression of z s .

z s = ∑ i∈V ∑ j∈V P ij y ij (59)
Algorithm 1 describes the mechanics of the iterative procedure. The user-defined parameters are the security threshold δ s used to build the conflict graph, and an improvement thresholds δ 2 for the second criterion, respectively. The algorithm starts by solving the program MIP: it finds the optimal solution for the first criterion of value z e . We compute z s the value of this solution for the second criterion. The point (z e , z s ) is a Pareto-optimal point, since it is globally optimal for the first criterion. The value of p is then used to add the constraint (60) to MIP:

∑ i∈V ∑ j∈V P ij y ij ≤ p -δ i (60)
Constraint (60) simply reflects the minimum improvement required on the second criterion.

The value of the parameter δ i can be considered as a factor of granularity of the Pareto front.

MIP is then solved once again, and the values of the two criteria are updated. The algorithm continues until the value of the second criteria becomes smaller than the threshold p f .

Algorithm 1 Iterative bi-objective optimization procedure (IBIOP)

1: procedure IBIOP(δ s , p f , δ i ) 2:
Input: Set of aircraft F , set of maneuvers M 3:

Parameters: security thresholds δ s , p f , improvement threshold δ i 4:

Build the conflict graph according to F , M and δ s 5:

z e ← +∞, z s ← +∞ 6:
Solve MIP 7:

z e ← optimal value of MIP 8:

z s ← value of second criterion for optimal solution of MIP 9:

while z s ≥ p f do 10: Add constraint ∑ i∈V ∑ j∈V P ij y ij ≤ p -δ i to MIP 11:
Solve MIP 12:

z e ← optimal value of MIP 13:

z s ← value of second criterion for optimal solution of MIP

Results

This section is organized as follows. Subsection 5.1 describes the values assigned to the different parameters for the experiments. Subsection 5.2 tests the validity of the assumptions made in the computations in Subsection 2.6. Computational results are detailed in Subsection 5.3. Subsection 5.4 provides a quantitative analysis of a Pareto front for a particular example.

Parameter values and simulations of the uncertainties distributions

In this subsection, we define the values assigned to the different parameters of the random variables distributions, and we describe the methods used to generate the random samples used for the Monte-Carlo simulations.

Parameter values of the uncertainties distributions.

We give the values assigned to the parameters of the distributions of the uncertainties.

Probability distribution of the wind:

The simulated wind follows a zero-mean normal distribution of standard deviation σ W = 5.4kt, according to the model described in [START_REF] Chaloulos | Effect of wind correlation on aircraft conflict probability[END_REF].

Probability distribution of the error on speed prediction: We follow the model presented in [START_REF] Chaloulos | Effect of wind correlation on aircraft conflict probability[END_REF], where the error on speed prediction is a zero-mean normal variable of standard deviation σ Υ = 7.9kt.

Probability distribution of the maneuver delays:

As no data on these delays exist to our knowledge, we interviewed an experienced air traffic controller to obtain an insight into what those values could be. As a result, we decided to use the following values:

• T c + T s ∼ N (µ T c,s , σ T c,s
) where µ T c,s = 30 seconds and σ T c,s = 10 seconds;

• ∀i ∈ F , T p i ∼ N (µ T p i , σ T p i
) where µ T p i is a random variable uniformly distributed between 20 and 40 seconds, and σ T p i = 10 seconds;

Monte-Carlo simulations.

For the simulation of aircraft trajectories, we generate random values for the wind according to the method developed by [START_REF] Lymperopoulos | Sequential Monte Carlo methods in air traffic management[END_REF]. The author performs a time and space discretization of the wind field, and iteratively computes at each time step the wind values at point of the grid according to the wind values computed at the previous time step, using correlation functions. The values of the normal distribution for the error on speed prediction and on the maneuver delay are generated according to the Box-Muller method described in [START_REF] Rubinstein | Simulation and the Monte Carlo method[END_REF], which simulates centered normal random variables using uniformly distributed random variables.

Validating the calculus through simulations

In this subsection, we validate the assumptions made in Subsection 2.6, by checking the validity of the approximation of the probability of conflict derived in Subsection 2.6 by comparing it to the probability obtained through simulations.

To this end, we study a test case representing a conflict situation with two aircraft i and j crossing each other with an angle θ (the set of values for theta is {60 • ; 90 • ; 120 • }). To avoid the conflict, they start their maneuver at 100NM from the crossing point of their trajectories. We designed different scenarios corresponding to a couple of maneuvers (M i , M j ) consisting in either speed or heading changes. The speed maneuvers range from -6% to 6% with a 1% step and the heading changes range from -10 • to 10 • with a 1 • step. In total, we have 545 scenarios. For each scenario, 2000 independent random samples are generated.

To validate the approximation of the probability of conflict derived in Subsection 2.6, for each scenario we computed the value of the approximated probability of conflict, and we simulated 2000 random scenario samples to estimate the probability of conflict. The same process was performed for the minimum distance between the two aircraft.

Configuration

Difference 1 highlights that the average simulated minimum distance is close to the computed one. Depending on the type of maneuver and the crossing angle, it ranges from 0.05 NM to 0.16 NM. The variance is really small, meaning that the difference between the two distances is usually close to the mean. Probabilities of conflict computed and simulated are almost identical, with a difference always under 1% for the considered instances, meaning that the expression derived in Subsection 2.6 is a really good approximation of the simulated probability.

Figure 5 highlights the results of Table 1 graphically, and focuses on two aircraft crossing with an angle of 90 • and performing heading changes to avoid the conflict. Figure 5a compares the average of the computed distribution of the minimum distance with the average minimum distance simulated. Figure 5b compares the computed and simulated probability of conflict. Scenarios are sorted on the horizontal axis following decreasing values of the compared quantities. The plotted functions are stepwise, because the difference in maneuver magnitudes from one scenario to the other are discrete, which has a significant impact on the minimum distance and the resulting probability of conflict. Random benchmark. This benchmark consists of random instances, where aircraft are uniformly distributed within a square sector with side length 50NM. To avoid generating infeasible instances, we perform a preprocessing before solving the problem: for each pair of aircraft that First, we observe that all the instances met the stopping criteria, yielding solutions with less than 0.01 expected conflicts. Results highlight that for instances with less than 10 aircraft, the solution time stays shorter than 30 seconds while on average 5 different solutions are generated. This achievement is meaningful for the air traffic controller, since he/she is able to access a small set of different solutions within a short period of time. For instances with more than 10 aircraft, solution times tend to slightly increase to reach up to two and a half minutes to generate 9 solutions (instance R 16 ). Even though it seems far from real-time, the advantage of our procedure is that it generates solutions on the fly, meaning that the controllers has at least one or two possible solutions within the first seconds of the resolution if he/she needs to act quickly. Those solutions have a higher probability of using a recourse, but by concept it gives to the controller a certificate that at least each pair of aircraft has less than 5% chances of needing another avoidance maneuver in the future. An observation worth mentioning is that the total solution time is closely linked to the value of the security threshold δ s and the stopping criterion. Indeed, their values will mostly influence the number of generated Pareto-front solutions. Moreover, for the experiments we choose to have a large set of possible maneuvers for each aircraft, which also has an impact on the solution time. A smaller set of maneuvers would result in a shorter execution time, especially if we delete the largest maneuvers, which represent the vertices with higher degrees. Indeed, those nodes have an impact on the number of cliques and on the search of a maximal clique of minimum weight.

Quantitative analysis of a solution set

We now focus on the solutions generated by our procedure applied to the instance depicted on Figure 7. One aircraft intersects a train of two aircraft separated by 20 NM. The algorithm parameters were identical to the ones for the other tests.

In order to get an insight into the effect of considering uncertainties on the chosen maneuvers, we remind that the deterministic solution is worth 5.54 kilograms of fuel. The generated solutions during the resolution are displayed on Figure 8. The first solution (6.66; 0.35), was computed in 0.29 seconds. Every half second, a new solution is generated until the eighth and last one (47.47; 0.009) which was computed in 6.22 seconds.

We can divide the generated solutions into two clusters, depending on their geometrical The interpretations of these results are two-fold. First, it shows that there is a discontinuity in the geometry of the generated solutions, instead of having the same pattern repeated with a different magnitude. Second, it provides the controller with a visual outlook of the solutions, and he/she can easily identify which ones will be easier to communicate, or which ones will be the more robust.

Depending on the preferences of the controller, the quality of the solution will differ regarding one criterion or the other. For instance, if the controller aims at efficiency, he/she will apply the solution costing 6.66 kilograms of fuel, but where in 35% of the scenarios he/she will need to issue new maneuvers in order to ensure separation. If he/she aims at saving potential workload, he/she will choose a solution with less than 1% chances of having to re-issue maneuvers, but costing 47.47 kilograms of fuel. If his/her preferences are more mixed,

Conclusions

In this article, we tackled the air conflict resolution problem under uncertainty. With our work, we provide the controller with a decision analysis tool generating a set of solutions representing different trade-offs between several criteria. The model is robust, since we consider a large span of uncertainties, including errors due to the wind, the imprecision of speed prediction. We also presented a new type of uncertainty: the delay in the execution of maneuvers. As a consequence, we cover a large part of the possible uncertainties that can arise during a conflict resolution. We computed the conflict probabilities and integrated them within an optimization model which is flexible, since it fully separates the modeling from the resolution process. Hence, the underlying mathematical framework remains valid, whatever the hypotheses considered. We focused on two objective that are relevant for conflict resolution, namely the extra fuel consumption, which serves as a performance index, and the expected number of conflicts, which serves as an indicator of potential additional workload required to re-issue avoidance maneuvers. We solve the problem by iteratively solving the aforementioned model, hence taking advantage of its power. Each resolution generates a solution of the Pareto front of the problem. At the end of the simulation run, the controller has a set of solutions where he can choose the one to apply, depending on the context and its preferences. Monte-Carlo simulations validated the theory, and intensive simulations highlighted interesting results. Complex instances with up to 20 aircraft are solved within seconds, and an average of 5 different solutions are generated within two minutes.

Further research will focus on the extension of the method to aircraft with changing altitudes, in order to consider a bigger variety of problems. A rolling-horizon procedure will also be the center of new research, allowing us to run continuous simulations of real-life data sets.

Figure 1 :

 1 Figure 1: Cross and along-track errors on an aircraft trajectory predicted position

Figure

  Figure 2: Structure of the maneuver delay

Figure 4 :

 4 Figure 4: Illustration of a maneuver delay for two aircraft performing heading changes

  et al. (2015b)) links the cliques in G to the CR problem: Proposition 3.1 Let C be a clique in graph G. Then C represents a set of conflict-free maneuvers for a subset of F of cardinality |C|.

Figure 5 :

 5 Figure 5: Comparison of calculus and simulations for two aircraft performing heading changes to avoid at the intersection of their planned trajectories making a 50 • angle

Figure 7 :

 7 Figure 7: Example for the analysis of the generated solutions

Figure 9 :

 9 Figure 9: Two different geometrical solutions for the example

Table 1 :

 1 Comparison of the simulation results and the calculus for the minimum distance and the probability of conflict Table

				of	Difference of
			separation distance	probabilities
	Maneuvers	Crossing	Absolute	Variance	Mean	Variance
		angle ( • )	Mean (NM)	(NM)	(%)
	H/H	60	0.09	0.01	0.68	0.02
	H/H	90	0.10	0.02	0.70	0.01
	H/H	120	0.13	0.01	0.61	0.04
	S/H	60	0.16	0.01	0.80	0.03
	S/H	90	0.12	0.01	0.72	0.12
	S/H	120	0.11	0.01	0.50	0.14
	S/S	60	0.05	0.01	0.62	0.02
	S/S	90	0.16	0.03	0.66	0.04
	S/S	120	0.14	0.01	0.66	0.11

Table 2 :

 2 Computational results

			Instance size	Solutions explored	Resolution
		|F | |V |	|E |	1 sol	D sol	T r	Nb p
	R 4	4	60	571	(12.26 ; 0.053)	(15.19 ; 0.005)	0.67	2
	R 6	6	90	1296	(32.12 ; 0.107)	(44.46 ; 0.009)	6.69	7
	R 8	8	120 2256	(77.04 ; 0.0917) (100.7 ; 0.008)	8.58	8
	R 10	10 110 2647	(151.9 ; 0.072)	(252.3 ; 0.008)	11.69	5
	R 12	12 156 6215	(337.6 ; 0.184)	(370.3 ; 0.007)	24.17	4
	R 16	16 208 10924	(562.2 ; 0.137)	(914.8 ; 0.008)	139.7	11
	F 1,60,10	2	30	83	(5.85 ; 0.035)	(9.13 ; 0.001)	0.09	2
	F 2,60,10	4	60	717	(20.73 ; 0.037)	(21.79 ; 0.002)	1.26	2
	F 4,60,10	8	120 3543	(55.90 ; 0.045)	(59.48 ; 0.009)	7.46	2
	F 6,60,10 12 156 6027	(134.7 ; 0.023)	(137.6 ; 0.008)	42.92	2
	G 2,10	8	120 3690	(73.15 ; 0.078)	(90.20 ; 0.008)	30.9	5
	G 3,10	12 156 11034	(382.8 ; 0.020)	(480.1 ; 0.007)	82.89	4
	U 15	15 195 6940	(14.42 ; 0.0512) (22.15 ; 0.002)	45.12	4
	U 25	25 325 12313	(65.07 ; 0.014)	(89.24 ; 0.004)	7.46	2
	U 35	35 455 33127	(89.15 ; 0.021) (120.05 ; 0.009)	130.02	2

Acknowledgement

This work was carried out under the project OPR-601 funded by the CRIAQ (Consortium on Research and Innovation in Aerospace in Quebec).

will loose separation within the first 30 seconds of observation, we randomly delete one of the aircraft. For a desired number of aircraft, we generate 15% more aircraft to anticipate the effect of the preprocessing. If more aircraft than desired remain after the preprocessing, extra aircraft are then randomly removed until the number is reached.

Computational results.

All tests were performed on a computer equipped with the following hardware: Intel Core i7-3770 processor, 3.4 GHz, 8-GB RAM. The algorithms were implemented in C++ and relies on CPLEX 12.5.1.0 CPL (2014) with default options to solve every instance. For instances with up to 10 aircraft, heading changes range from -10 • to 10 • with a 2 • step. For instances with more than 10 aircraft, the possible heading changes are ±5 • , ±10 • , ±15 • and ±20 • . In addition, aircraft can also perform speed changes of ±3% and ±6%. The parameters of the IBIOP procedure were assigned the following values:

• security threshold for a maneuver δ s = 5%;

• improvement threshold for safety δ i = 1%;

• the stopping criterion p f ≤ 1%.

Table 2 gathers information about the instance dimensions, the generated solutions and computational results. The headings are given as follows:

• |F |: number of aircraft;

• |V |: number of vertices;

• |E |: number of edges;

• 1 sol : first generated solution of the Pareto front, expressed as the pair (c 1 , c 2 ) of the values of the two criteria, where c 1 is expressed in kilograms of fuel and c 2 is the expected number of conflicts;

• L sol : last generated solution of the Pareto front, expressed as the pair (c 1 , c 2 ) of the values of the two criteria, where c 1 is expressed in kilograms of fuel and c 2 is the expected number of conflicts;

• T r : resolution time (in seconds);

• Nb p : number of generated solutions of the Pareto front.


peripheral indole groups stabilize benzene-1,3,5-tricarboxamide helical assemblies

Gaëtan Basuyaux, a Alaric Desmarchelier, a Geoffrey Gontard, a Nicolas Vanthuyne, b Jamal Moussa,* a Hani Amouri, a Matthieu Raynal,* a and Laurent Bouteiller The bottom-up assembly of small synthetic molecules provides an efficient strategy for the construction of artificial helices and their implementation in chirality-based applications. 1 In this context benzene-1,3,5-tricarboxamide (BTA) 2 molecules constitute an iconic class of monomers given their ability to predictably form helical assemblies with a preferred handedness, and their utility as asymmetric catalysts, 3 circularly-polarized light emitters 4 and chiral templates. 5 Stacks of BTA molecules in the form of a threefold hydrogen bonded helix connecting the amide units has been characterized in the crystalline state 6 in the liquid-crystalline state, 7 in gel 8 and in solution. 2a,9 An important challenge for imparting these supramolecular polymers with further functionality is to introduce functional groups in their side chains without destabilizing the helices. Comparatively to BTAs bearing alkyl chains (alkyl BTAs), a significant destabilization of the assemblies has been reported for BTAs endowed with short ethylene glycol chains 10 or derived from -amino esters (ester BTAs).

11

In the case of ester BTAs, the hydrogen bond competing nature of the ester group is reflected by the possible formation of a dimeric structure in which the N-H moieties are bound to the ester C=O, instead of amide C=O in the case of helical BTA stacks.

11c-e

Moreover, the nature of the substituent attached to the aminoester -carbon significantly influences the assembly properties of ester BTAs. Previous studies revealed that stereoelectronic 11d and conformational effects 11e can strongly shift the assembly of ester BTAs from dimers to stacks. Also, weak non-covalent interactions involving the side chains of the BTAs derived from Lphenylalanine 11a, 11c, 11e, 11f or L-methionine 6e, 11c, 11e help stabilize the helical assemblies in apolar solvents. Nevertheless, these helices are less stable, by at least one order of magnitude, than alkyl BTAs.

11e

Much effort has been devoted towards the stabilization of -helical peptides and mimetics 12 and of the collagen triple helix, 13 which can be achieved by chemical mutation of a particular amino-acid, amongst other possible strategies. Moreover, the helical conformation of synthetic covalent polymers such as poly(isocyanide)s 14 and poly(acetylene)s 15 can be stabilized by means of hydrogen bonding interactions involving the polymer side chains. Such a toolbox for the stabilization of artificial supramolecular helices by strong interactions between peripheral groups is currently lacking. [START_REF] Masuda | For covalent fixation of helical BTA assemblies see[END_REF] Taking into consideration the possible array of interactions involving the indole moiety, [START_REF] Cochran | Proc. Nat. Acad. Sci[END_REF] we became interested in the self-assembly properties of BTA monomers derived from tryptophan alkyl esters (BTA Trp). We report herein that BTA Trp monomers, thanks to additional hydrogen bonding interactions provided by the indole groups, self-assemble into remarkably stable supramolecular helices in apolar solvents. Please do not adjust margins

Please do not adjust margins Examining the packing arrangement in the crystal shows that (S)-BTA Trp i-Pr molecules stack upon each other and form supramolecular helices by means of a double hydrogen bond network: (i) an inner helical network involving the amide functions only (see the red dotted lines in Fig. 1a and the strands colored with different shades of grey in Fig. 1b) similar to the one observed in the crystalline structures of other alkyl and ester BTAs. 6 The distance (N1 S1). The refinement of the Flack x parameter, equal to 0.00 (17), ascertains the absolute configuration of the structure, i.e. (S)-BTA Trp i-Pr assembles into right-handed helices which is in agreement with the handedness of helices found for other ester BTAs having (S)-configuration. 6d, 6e (ii) an outer hydrogen bond network involving the indole N-H groups and the amide C=O (see the green dotted lines in Fig. 1a and the red, green and blue colored strands in Fig. 1b). The N•••O distance (2.84(1) Å) is consistent with that found for indole N-H group of Trp residues involved in hydrogen bonding interactions in peptides and proteins. 17b The NH•••O angle of only 139(1)° is likely due to the bifurcated nature of the hydrogen bond since the acceptor, the amide C=O, is already engaged in the inner hydrogen bond network. The two threefold hydrogen bonded helices, which adopt opposite handedness, are emphasized in Fig. 1b. Importantly, the additional interaction provided by the amino ester side chain allows each consecutive monomer inside the stacks to be maintained by six hydrogen bonds (Fig. 1a), compared to three in previously reported BTA structures, thus potentially strengthening the overall assembly (vide infra). , 293 K), zoom on the N-H region (see the C=O region in Fig. S1b). The FT-IR spectra of the models (indole in acetone, 5×10 -3 mol.L -1 , and BTA Aib in MCH, 10 -2 mol.L -1 ) are shown with an offset of 600 L.mol -1

.cm -1

. The simulated spectrum is obtained by addition of the FT-IR spectra of indole (×3) and of BTA Aib.

A set of experiments has been carried out in order to probe whether this unusual type of helices with a double hydrogen bond network also exists in solution. (S)-BTA Trp C12 has been selected for the solution studies given its good solubility in apolar solvents. Firstly, the FT-IR spectra of (S)-BTA Trp C12 in methylcyclohexane (MCH) and in solid state have been compared to that of solid (S)-BTA Trp

i-Pr

. The three FT-IR spectra exhibit the same bands in the N-H and C=O regions which suggests that (S)-BTA Trp C12 and (S)-BTA Trp i-Pr form the same type of self-assemblies (Fig. S1). Then, the nature of these bands have been definitively established by using BTA Aib 11e in MCH and indole in acetone as models of the inner and outer hydrogen bonding interactions, respectively, found in the structure of (S)-BTA Trp i-Pr . Thus, the bands with absorption maxima at = 3361 cm -1 and = 3235 cm -1 correspond to indole N-H and amide N-H groups, respectively, both hydrogen bonded to amide C=O (Fig. 2). As a matter of comparison, free indole and amide N-H bands both exhibit maxima at ca. = 3440 cm -1 (Fig. S2a). The simulated FT-IR spectrum, obtained by combining those of the models, resembles the experimental one except that the indole N-H band is broader which is likely due to a less directional hydrogen bonding interaction for indole in acetone relatively to the same interaction taking place in the helical structure of (S)-BTA Trp C12 . The absorption maxima corresponding to ester band (= 1742 cm -1 , free) and amide I band (= 1633 cm in the crystalline state (Fig. S2b).

The Small-Angle Neutron Scattering (SANS) analysis of (S)-BTA Trp C12 in C 7 D 8 (4.4×10 -3 mol.L -1 ) shows a q -1 dependency of the scattered intensity at low q values which is characteristic of the presence of rigid one-dimensional objects (Fig. 3). ‡ The data extracted from the fit indicate that these objects are very long cylinders (L>1000 Å, DP w >275) with a radius of 12.9 Å which is close to the one of the isolated molecule in its fully extended conformation (17 Å). § This analysis thus

Please do not adjust margins

Please do not adjust margins , 293 K). The curve was fitted ( ) according the form factor for rigid rods of infinite length with a circular cross section and a uniform scattering length density (L =length, r =radius). , MCH, Fig. 4). § § In contrast to monomers, assemblies of (S)-BTA Trp C12 display a set of intense CD signals with a negative Cotton effect centered at ca. 280 nm and two successive bisignated CD signals (¯= 237 nm/⁺= 227 nm and ¯= 218 nm/⁺= 207 nm). The wellstructured and intense CD band at 280 nm is related to the longer wavelength absorption maxima of the indole unit. CD spectra of previously studied BTA assemblies usually display a single (broad) CD signal in the 200-250 nm region, as it is the case for the BTA derived from the dodecyl ester of L-Phenylalanine (see the overlaid spectra in Fig. S5). 11c This corroborates that, in (S)-BTA Trp C12 assemblies, the indole side chains are located in a well-defined chiral environment. As expected, (S)-BTA Trp C12 and (R)-BTA Trp C12 feature mirrorimage CD spectra reflecting the opposite handedness of their helical assemblies. This full set of analyses (FT-IR, SANS and CD) demonstrates the ability of (S)-BTA Trp C12 monomers to assemble into long, onedimensional helical stacks in solution similar to the ones reported for other BTA assemblies but with additional hydrogen bonding interactions involving the indole moieties of the side chain. FT-IR (Fig. S6) and CD (Fig. S7) analyses in MCH indicate that at room temperature the helical assemblies are identical over several order of magnitudes of concentrations (2×10 -2 mol.L -1 -2.0×10 -6 mol.L -1

). Variable-temperature CD and UV-Vis analyses (Fig. S8) have then been performed to probe the thermal stability of the helical assemblies under very diluted conditions (2.0×10 -6 mol.L -1 ). The transition occurs between helical stacks (with a double hydrogen bond network) and monomers at a remarkably elevated temperature (≈ 360 K). [START_REF]In this highly diluted condition, the error bar in the measurement of the ellipticity (< 5 mdeg) is too[END_REF] It infers that these assemblies are more stable, by at least one order of magnitude, than those formed by the other ester BTA assemblies 11e and are at least as stable as alkyl BTAs investigated to date. 7b, 19 The importance of the extra hydrogen bonds in stabilizing the helical assemblies of (S)-BTA Trp C12 is further supported by the fact that (S)-BTA NMe Trp, which lacks the indole N-H group, does not form long assemblies in MCH (Fig. S9). Finally, we have examined whether the high stability of these helical assemblies, explained at the molecular scale by the six hydrogen bonds connecting adjacent BTA monomers within the stacks, translated into an improved property at the macromolecular level. With this aim in mind, we have prepared 2×10 -2 mol.L -1 toluene solutions of (S)-BTA Trp

C12

, of (S)-BTA NMe Trp and of a set of ester BTAs and alkyl BTAs previously-found to form helical assemblies in solution. From all of these samples, only (S)-BTA Trp C12 forms a self-standing gel for more than 2 minutes after inverting the tube (Fig. means that the overlap concentration is below 0.1 wt%, which yields a lower estimate for degree of polymerization of the stacks: DPn>290 at a concentration of 0.1 wt%. 20 In this study, we brought to light the ability of the peripheral indole N-H groups to stabilize helical BTA assemblies through hydrogen bonding interactions with the amide C=O. Without these additional interactions, the presence of ester C=O as a hydrogen bond acceptor would actually have destabilized the formation of supramolecular polymers under the same conditions. As a consequence of this extra stabilization, helical assemblies of BTA Trp are very long and remain stable at low concentrations and high temperatures. In BTA Trp, the supramolecular chirality is expressed at two levels, the inner and outer hydrogen bond networks, and this may offer the opportunity to amplify the chirality to a greater extent compared to classical BTAs. Studies towards this goal are underway in our laboratory as well as the implementation of these assemblies as scaffolds for asymmetric catalysis. This work was supported by the Labex MiChem, part of French state funds managed by the Agence Nationale de la Recherche (ANR) within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02, and by the ANR (project ANR-13-BS07-0021 SupraCatal). The GDR 3712 Chirafun is acknowledged for stimulating discussion. Jacques Jestin (LLB, Saclay) is acknowledged for assistance with SANS experiment.

  Trp) self-assemble into helices with an inner threefold hydrogen bond network surrounded by a second network involving the indole N-H groups. As a consequence of this extra stabilization of its helical assemblies, BTA Trp forms more viscous solutions than a range of ester and alkyl BTAs.

Chart 1

 1 (a) Chemical structures of the BTA monomers investigated in this study. (b) ORTEP representation of an isolated molecule of (S)-BTA Trp i-Pr together with heteroatom labels. BTA Trp monomers (Chart 1) with dodecyl (BTA Trp C12 ) or isopropyl (BTA Trp i-Pr ) ester chains have been prepared via conventional synthetic procedures and have been isolated in pure form according High-Resolution Mass Spectrometry and Fourier-Transform Infrared (FT-IR) analyses (ESI †). Chiral HPLC analyses certify that both (S)-BTA Trp C12 and (S)-BTA Trp i-Pr are enantiopure (de> 99%, ee> 99%). (S)-BTA NMe Trp has been synthesized from N-methyl (S)-tryptophan with the purpose of determining the influence of the indole N-H group on the selfassembly properties of BTA Trp monomers. Crystals of (S)-BTA Trp i-Pr have been obtained by cooling of a concentrated solution in CH 3 CN. (S)-BTA Trp i-Pr crystallizes in the hexagonal P6 3 space group, a view of the isolated molecule is shown in Chart 1.

Fig. 1

 1 Fig. 1 X-ray structure of (S)-BTA Trp i-Pr . (a) ORTEP representation of a trimer of (S)-BTA Trp i-Pr with the two types of hydrogen bonds present in the helical polymer shown in green (indole/amide) and red (amide/amide). The ester side chains (CO2i-Pr) and hydrogen atoms are omitted for clarity. The six hydrogen bonds connecting the lower monomer pair are shown with a thicker dotted line. Selected distances (Å): N1•••O1: 2.89(1), N2•••O1: 2.84(1). Selected angles (°): N1-H1A•••O1: 162(1), N2-H2•••O2: 139(1). (b) Diagram of a hexamer of (S)-BTA Trp i-Pr highlighting the inner helical threefold hydrogen bond network (O, C, NH helical strands shown with different shades of grey, BTA ring in black) surrounded by the second helical network involving the indole moiety (CH, CH2, N ring corresponding to different strands are shown in blue, red and green). The amide C=O are common to both hydrogen bond networks.

Fig. 2 FT

 2 Fig. 2 FT-IR spectrum of (S)-BTA Trp C12 in MCH (10 -2 mol.L -1

Fig. 3

 3 Fig. 3 SANS analysis ( ) of (S)-BTA Trp C12 in C7D8 (5.61 g.cm -3 , 4.4×10 -3 mol.L -1

Fig. 4

 4 Fig. 4 CD spectra at 293 K of (S)-BTA Trp C12 (10 -4 mol.L -1 , MCH, stacks), (S)-BTA Trp C12

Fig. 5

 5 Fig. 5 Top: Pictures of 2×10 -2 mol.L -1 toluene solutions of various ester and alkyl BTAs 2 minutes after inverting the vials. From left to right: (S)-BTA Trp C12 , (S)-BTA NMe Trp, (S)-BTA Met, BTA Aib, BTA C8 and BTA C8*. (S)-BTA Met: BTA derived from the dodecyl ester of L-methionine, BTA C8: BTA endowed with thee octyl side chains. BTA C8*: BTA endowed with two octyl side chains and one (S)-1-methylheptyl side chain. See Chart S1 for the chemical structures of these BTAs. Bottom: Relative viscosity of (S)-BTA Trp C12 and BTA C8* solutions in toluene versus concentration (293 K).

Conflicts of interest

There are no conflicts to declare.

Notes and references ‡ (S)-BTA Trp

C12 adopts a similar supramolecular structure in toluene and MCH as indicated by the virtually identical FT-IR spectra obtained in these solvents (Fig. S3). § As obtained by geometry optimization using molecular mechanic method (MM3) with the software Scigress (Fujitsu).

§ § (S)-BTA Trp

C12 is fully assembled down to 2.0×10 -6 mol.L -1 at 293 K in MCH as demonstrated by FT-IR and CD analyses performed at various concentrations (Figs. S6 and S7, respectively).