Business Research Ethics Paper

European Journal of Operational Research 241 (2015) 502–512 Contents lists available atScienceDirect European Journal of Operational Research journal homepage:www.elsevier.com/locate/ejor Innovative Applications of O.R.

Solving air traﬃc conﬂict problems via local continuous optimization Clément Peyronne a,∗, Andrew R. Conn b, Marcel Mongeau c,d , Daniel Delahaye c,d aCapgemini, 15 av. du Dr Maurice Grynfogel, 31000 Toulouse, FrancebIBM, T.J. Watson Research Center, P.O. Box 218, Yorktown Heights, NY 10598, USAcENAC, MAIAA, F-31055 Toulouse, FrancedUniversité de Toulouse, IMT, F-31400 Toulouse, France article info Article history:

Received 29 August 2012 Accepted 31 August 2014 Available online 28 September 2014 Keywords:

Air traﬃc conﬂict problem B-splines Continuous optimization Genetic algorithms Semi-inﬁnite programming abstract This paper ﬁrst introduces an original trajectory model using B-splines and a new semi-inﬁnite program- ming formulation of the separation constraint involved in air traﬃc conﬂict problems. A new continuous optimization formulation of the tactical conﬂict-resolution problem is then proposed. It involves very few optimization variables in that one needs only one optimization variable to determine each aircraft trajec- tory. Encouraging numerical experiments show that this approach is viable on realistic test problems. Not only does one not need to rely on the traditional, discretized, combinatorial optimization approaches to this problem, but, moreover, local continuous optimization methods, which require relatively fewer iterations and thereby fewer costly function evaluations, are shown to improve the performance of the overall global optimization of this non-convex problem. ©2014 Elsevier B.V. All rights reserved. 1. Introduction This introduction describes the importance of the conﬂict- resolution problem, followed by some details on related work. 1.1. Context Air traﬃc management (ATM) aims at ensuring smooth running of the transportation system under safety and schedule alignment con- straints while keeping ﬂights on schedule. In order to reach this goal, air traﬃc is planned at different time frames.Strategicalplanning is done several months before take-off and consists of assigning ﬂight plans for a whole day of traﬃc with an emphasis on an even distribu- tion of aircraft density in space and time.Pre-tacticalplanning then updates the strategical planning trajectories using information such as weather or airspace congestion. It takes place two hours before the aircraft reaches the considered airspace.Tacticalplanning, the subject of this study, is performed within a 20-minute time horizon and con- sists mainly of conﬂict detection and resolution. This tactical planning work has always been done by air traﬃc controllers who are in charge of the most critical aspect of ATM, namely, ensuring suﬃcient sep- aration between airplanes. Air traﬃc controllers are responsible for ∗Corresponding author. Tel.: +33 627025703.

E-mail addresses:[email protected](C. Peyronne), [email protected](A. R. Conn),[email protected](M. Mongeau), [email protected](D. Delahaye). ensuring the respect of regulatory separation rules that are currently 5 nautical miles horizontally and 1000 feet vertically.

Aconﬂictsituation happens when an aircraft enters the standard safety zone of another aircraft (a situation where a regulatory sepa- ration rule is not respected). Our work focuses on tactical planning, which is currently handled by air traﬃc controllers on each airspace sector. To solve a conﬂict, controllers can use a set of maneuvers: off- set, turning point, speed change, and ﬂight-level change. For instance, Fig. 1shows maneuvers implying direction changes.

As a consequence of increasing traﬃc, controllers in charge of an airspace sector must handle more and more ﬂights (seeFig. 2). The current approach is to decrease the size of the control sectors in order to compensate for the growth of traﬃc. However, the traﬃc is reach- ing the point where a decrease of the size of sectors is no longer eﬃ- cient. In reality, ATM has already used every available resource in an attempt to increase airspace capacity. However, from now to 2030, air traﬃc is expected to increase by a factor of two or three (SESAR Joint Undertaking, 2009). Consequently, ATM will have to deal with this overload while ensuring at least similar standards of safety (SESAR Joint Undertaking, 2009). As is illustrated inFig. 2, the difference of traﬃc capacity between 1970 and 2010 is rather dramatic.

In this context, the long-term vision aims at lowering the workload of the air traﬃc controllers by reducing the conﬂict-resolution task.

An example is the 4D-trajectory concept, which consists of deﬁning precisely a trajectory in space and time. One option is to create an automatic conﬂict-resolution tool to provide advisory solutions to the controller. Some previous work has been done in the direction http://dx.doi.org/10.1016/j.ejor.2014.08.045 0377-2217/ ©2014 Elsevier B.V. All rights reserved. C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512503 Initial trajectories Turning point Oﬀset Fig. 1.Two maneuvers implying direction changes. Fig. 2.Evolution of traﬃc, number of sector, controllers and the number of ﬂights per controller.

of automatic-conﬂict resolution (seeBosc, Durand and Maugis 1997) but to our knowledge, none has been used operationally yet.

1.2. Previous related work One approach for air traﬃc automatic conﬂict resolution uses nav- igation functions (Dimarogonas & Kyriakopoulos, 2005; Roussos & Kyriakopoulos, 2009). What follows is a presentation of the repul- sive force technique that inspired the navigation function methods (the repulsive force technique is much more intuitive and easier to comprehend). A repulsive force technique considers the airspace as a potential ﬁeld, and aircraft as particles navigating in it. Negative charges represent obstacles for the aircraft (other aircraft, congested areas). The destination is associated with a positive charge. As a result, each aircraft is attracted by its destination while being repulsed by obstacles, as is illustrated byFig. 3. Each trajectory is determined in its own functional space so that an aircraft cannot be attracted by the destination of another aircraft. This enables the automatic generation of conﬂict-free trajectories.

Navigation function methods were shown to ensure collision avoidance while connecting the departure and destination points. The major drawback of navigation functions is that the obtained solution does not necessarily respect ATM constraints such as the particular bounded speed (an aircraft cannot ﬂy below or above a certain range of speed), or trajectory smoothness. Furthermore, they can lead to major delays and overcosts as they tolerate large deviations from the direct route, i.e. the straight line between the departure and arrival points.

Optimization methods have also been used in air traﬃc automatic conﬂict resolution. InPallottino, Feron and Bicchi (2002), two ap- proaches based on a local optimization method are presented. One approach considers using only speed changes, and the other relies solely on direction changes. However, the speed-change approach cannot solve every conﬂict situation (simply consider a face-to-face situation for example), and their direction-change approach is re- stricted to straight-line maneuvers.

Further developments are presented inAlonso-Ayuso, Escudero, and Martin-Campo (2011)combining speed and altitude changes and using a mixed-integer linear optimization approach. This method pro- vides very interesting results both from the point of view of computa- tional time and the quality of resolution on conﬂict situation involving up to 50 aircraft (which are not all involved in the same conﬂict how- ever). InAlonso-Ayuso et al. (2011)a method is presented that relies on exact optimization, which is a signiﬁcant advantage. However, the use of altitude changes is a drawback as it induces costly maneuvers that are avoided by air traﬃc controllers due to their high costs in fuel that are unacceptable for airline companies.

Remaining optimization methods for this problem rely on heuris- tics. The authors ofDurand and Alliot (1995)andMédioni, Durand and Alliot (1994)obtain relatively good results on real traﬃc using Genetic algorithms (GA). However, their approach is restricted to offset and turning-point maneuvers, i.e. piecewise-linear trajectories. Similarly, using ant colony optimization and modelling trajectories as a path in a graph,Durand and Alliot (2009)andOlive (2006)are conﬁned ObstaclePotential line AircraftRepulsive force Sliding force Aircraft 1Aircraft 2 Sliding forceSliding force Potential line Potential line Relative speed 2/1 Relative speed 1/2 Fig. 3.Repulsive force between an aircraft and an obstacle (left), between two aircraft (right). 504C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512 to produce piecewise-linear trajectories. Finally, another recent so- lution for tactical conﬂict resolution, presented inDougui, Delahaye, Puechmorel and Mongeau (2012), called LPA (Light Propagation Al- gorithm), is inspired by an analogy with light propagation. LPA uses a Branch-and-Bound technique to build trajectories and obtains good results on a full day of traﬃc over France. However, the trajectories are built sequentially, which can cause signiﬁcant deviations for some aircraft. For more detail on air traﬃc conﬂict detection and resolution, see the reviews byDelahaye and Puechmorel (2013)andKuchar and Yang (2000).

1.3. Limitation of scope There are three means to solve an air traﬃc conﬂict: speed changes, direction changes and altitude changes (of course, reduc- ing the traﬃc also naturally decrease the number of conﬂict). This paper focuses on tactical en-route conﬂict resolution, which in our case means we consider a 20-minute time horizon. Thus, this study considers only direction changes, as speed changes are not as eﬃ- cient on such a short time frame, and as altitude changes are only used as a last resort by air traﬃc control (because of their cost and the passenger discomfort they engender). That is why, in this paper, each aircraft is also assumed to preserve its imposed vertical proﬁle.

Consequently, as the study concentrates on en-route traﬃc, aircraft are considered to ﬂy at a stable altitude, which means solving the conﬂicts solely intwodimensions.

This paper presents a tractable and practical method to solve tac- tical conﬂicts on the considered time horizon, involving only a small part of the trajectory (corresponding to the considered 20-minute window). One of the contributions of this paper is to model trajec- tories with B-splines. The overall trajectory is to be managed by an algorithm that aggregates the separate 20-minute windows. To treat a large air traﬃc instance (for example the French en-route air traﬃc), the conﬂict solver presented in the paper is to be applied at regular time intervals using a moving time-window. Furthermore, the solver can be applied on different geographical zones at the same time. See Dougui et al. (2012)for more details on this moving time-window process. For these reasons, computational eﬃciency is an important aspect of conﬂict resolution. The concept of a moving time-window (decreasing the number of considered aircraft) and the need for com- putational eﬃciency are the main reasons why local optimization has been tested in this study. 1.4. Problem description and overview As already mentioned, the context of this study is the tactical plan- ning phase. The aim here is to obtain, from a conﬂicted situation, an optimal conﬂict-free solution by deviating the trajectories smoothly andaslittleaspossible.

This paper introduces an original trajectory model using B-splines and a new semi-inﬁnite programming formulation of the separation constraint involved in air traﬃc conﬂict problems. Another related contribution consists of a new continuous-optimization formulation of the tactical conﬂict-resolution problem that involves very few op- timization variables: only one real optimization variable per trajec- tory. Finally, encouraging numerical experiments are obtained us- ing genetic algorithms, a ﬁnite-difference interior-point method, and derivative-free optimization.

The input data of the problem consist of :

1)N: the number of aircraft involved in the conﬂict 2) αi start and αi end : the initial and ﬁnal positions of the aircrafti(i= 1,2,...,N) respectively, as illustrated inFig. 4.

3)v i: the speed of aircrafti, which is assumed constant.

Note that each piece of trajectory considered is ﬁrst to be modelled by a B-spline, deﬁned by a single continuous parameter. The output α4 start α1 start α2 start α3 start α4 end α3 end α1 end α2 end Aircraft 3 Aircraft 2 Aircraft 4 Aircraft 1 Fig. 4.A typical input of the problem: one conﬂict involving trajectory segments (in bold).

should consist of smooth conﬂict-free trajectories, with trajectoryi going from αi start toαi end at constant speedv i. 2. B-spline trajectory model and decision variables The use of smooth trajectories is not possible in today’s operational context as air traﬃc control is currently restricted to tra- jectories involving linear segments connected with small constant- curvature turns. However, the challenging research projects SESAR and NextGen consider smooth trajectories as an option for the future Flight-Management Systems. The ability of future FMS (Flight Man- agement System) to ﬂy such trajectories opens new opportunities for optimization with respect to environmental criterion. In this context, the use of cubic B-splines (smooth piecewise-cubic polynomials) to design aircraft trajectories allows one to describe an aircraft trajectory deviation between the points αstart and αend with a single continuous parameter. These parameters (one per aircraft) will be the optimiza- tion variables in this study. Furthermore, it should be clear that it is more desirable in practice to use smooth trajectories. This section describes the main ingredients of B-spline theory in a more general context. 2.1. Elements of B-spline theory B-splines are parameterized curves determined by a set of points calledcontrol points. One considers here a set of control points (Xk,Yk)∈R 2. One can deﬁne γ(s), the B-spline curve determined by the control points (Xk,Yk)as follows: γ(s)=(γx(s),γy(s)),s∈[0,N cp−1], wheresis thenatural parameterof the B-spline, andN cpis the number of control points. The B-spline curve is obtained as a linear combina- tion of the B-spline function basis (Bk)k=−1,...,N cp. Each element,B k,of the function basis is a cubic-polynomial function. The basis function B 0is the interpolation natural cubic spline of the following points: (−2,0 ),(−1, 16),(0,23),(1,16),(2,0 )centered on 0, and the remaining elements of the basis are obtained by simple translations ofB 0. C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512505 Fig. 5.The four B-spline basis functions relevant for the interval [1,2]. The B-spline basis functions used to calculate the B-spline curve γ(s)fors∈[1,2] are the ones with a non-null value on this interval [1,2] (in bold inFig. 5). Thus, the unique 2D B-spline curve that is determined by theN cp control points (Xk,Yk)along with three so- calledphantom pointsis given by the following linear combination of theB k’s : γx(s)= Ncp+1 k=−1XkBk(s)and for γy(s)= Ncp+1 k=−1YkBk(s).

The B-spline ﬁtting curve is a piecewise cubic polynomial function ofs,wheretheknots,s k, are the points where the pieces join. Denote γk x,dγk xds,d2γk x ds2and d3γk x ds3respectively the B-spline desired values: the B-spline desired ﬁrst, second and third derivative values fors=s k.

The following requirements yield standard B-spline properties for unit knot intervals: γk x=X k+1 +4X k+X k−1 6, d γk x ds=X k+1 −X k−1 2, d 2γk x ds2 =X k+1 −2X k+X k−1 , d 3γk x ds3 =X k+2 −3X k+1 +3X k−X k−1 .(1) Since the Taylor’s expansion is exact for polynomials (providing enough terms are taken) one then obtains the B-spline value for all s∈[s k,sk+1 ]: γx(s)=γk x+(s−s k)dγk x ds+ (s−s k)2 2d 2γk x ds2 +(s−s k)3 6d 3γk x ds3.(2) Higher-degree derivatives are zeros, as γx(s)is a cubic B-spline. As we use three control points ( αstart and αend and one movable, α(u)), there are three knots, inducing two sub-intervals. The trajectory γ(s) is represented on the interval [0,2] and on the two sub-intervals [0,1] and [1,2]. For more details, the B-spline theory is described, for example, inDuncan (2005).

In the following, the indicesiandjwill be used to determine the considered aircraft, andkwill be used to determine the considered control point. Thus, γistands for the trajectory of theith aircraft, and γkfor the value of γat the knots kfor the considered trajectory. Start PointEnd Point Moveable Point Fig. 6.A B-spline trajectory determined by three control points. B-spline ﬁtting is a very eﬃcient tool for trajectory modelling in terms of both ﬁtting quality and computational time. Moreover, B-splines feature interesting properties such asC 2-continuity, which is crucial for modelling ﬂyable smooth aircraft trajectories. In addi- tion, a very attractive property is the fact that, by construction, the B-splines minimize the quantities: (γ x(s))2dsand (γ y(s))2ds. This is important in an operational context as it induces low energy con- sumption and passenger comfort.

B-splines have already been used in trajectory optimization and one can ﬁnd an example inMilam, Mushambi and Murray (2000). 2.2. Optimization variables deﬁning the trajectories In this study, one models a trajectory with a cubic B-spline de- ﬁned by three control points: the start and end points (given input data), plus a middle control point (Fig. 6). This middle control point is called themovable control point, as it will be used to deviate the trajectory. The exact position of this control point will be monitored by a single real-valued parameter. The vector of decision variables of the optimization problem we are about to deﬁne will be made of these real-valued parameters.

One wants a compromise between allowing the trajectory to deviate freely (in any direction) from the direct route in order to avoid conﬂicts, and staying as close as possible to the direct route.

For that purpose, one deﬁnes a ﬁxed maximal bandwidth (interval 506C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512 Dmax Dinit pi mid d αi start αi(ui) α i end Fig. 7.Direct route, maximal deviation bandwidth and modiﬁed trajectory. [−D max ,Dmax ]) that depends on the direct route length (D init). The trajectory will be imposed to stay within the maximal bandwidth (seeFig. 7).

The middle control point coordinates corresponding to the trajectory of aircrafti(i=1,2,...,N) is denoted αi(ui).Itsexactposi- tion will be determined using the parameteru iwhich will represent a bandwidth proportion. It is deﬁned as:u i=± dDmax (percent), where D max =λDinit,dis the distance between the control point αi(ui)and the direct route for aircrafti,and λis a user-deﬁned proportion. In the numerical tests, λis set to 0.3, so as to expect a maximal distance increase of about 10 percent from the original trajectory (this was observed empirically).

The optimization variable of the problem is therefore the vector uwhoseith component is that percentageu icorresponding to tra- jectoryi=1,2,...,N. The search space is therefore the hypercube [−100,100] N.

In order to calculate the deviation of each aircraft (which will be required for computing the objective function of the optimization formulation of the next section), one must ﬁrst build the trajectories driven by the above-deﬁned control points. From a given deviation percentage vector,u, one calculates the control point locations, and the B-splines driven by them. In order to illustrate this, consider the trajectory of a single aircraft,i. First, let us deﬁne thedirectionvector δi:= αi end −αi start of theith direct route. The middle of the original trajectory direct route is called thesiteand is denoted byp i mid .The middle control point αiis placed on the line orthogonal to the direct route and intersecting the direct route atp i mid (seeFig. 7). In other words, pi mid =αi start +12δi.(3) The vector of normD i max (maximal deviation bandwidth of theith aircraft) that is orthogonal to the direct route, and, by convention, oriented toward the left of the aircraft trajectory, is denoted byq i∈ R 2. More precisely, it is obtained by solving the system:

(δi)Tqi=0 ||q i|| =D i max .(4) The ﬁrst line of this system ensures the orthogonality betweenqand δ, while the second line deﬁnes the length ofq. Solving this system leads to two possible vectorsq. To choose between those two possibilities, we deﬁne a convention stating thatq ishould be oriented toward the left of the original trajectory. Consequently, ifu i>0, theith aircraft will turn left. Hence, one can write the following relation between the control point αi(u)and theith optimization variableu icorresponding to aircrafti: αi(ui)=p i mid +u iqi.(5) Thus, givenu i, one knows the three required control points ( αi start , αi(ui)and αi end ), and can compute the B-spline deﬁning theith aircraft trajectory. 3. Optimization problem formulation The problem consists in designing, if possible, an optimal conﬂict- free situation. This section proposes an optimization formulation of this problem whose constraints will require that there is no conﬂict.

The objective function of the problem can be deﬁned in several ways:

one can consider minimizing the total deviation distance or total fuel consumption for instance. This paper concentrates on ﬁnding conﬂict- free trajectories that minimize the average total deviation distance (with respect to the direct routes). The formulation presented can easily be adapted to other objective functions. 3.1. Formulation A new idea introduced in this paper is to consider the continuous separation distance as a constraint, making the approach more direct than previous methods that rely on a discretization of the trajectories, not only for numerical purposes but in the actual formulation of the problem to be solved. However, this constraint must be respectedat all times, leading to an inﬁnite number of constraints. As described above, one usesu, the vector containing all the middle control point locations,u i,i=1,2,...,Nas the optimization variable.

Letf (u)be the objective function representing the average trav- elled distance for every aircraft (details are given inSection 3.2), which is written as follows: f(u)=1 N N i=1 fi(ui),(6) wheref i(ui)is related to the travelled distance of the aircrafti following its modiﬁed trajectory. This function is further explained inSection 3.2.

In order to guarantee a conﬂict-free situation, one must ensure that all aircraft are, at least, τ:=5 nautical miles away from each other (in practice, τ=5.001 is taken to give a margin tolerance).

The separation-norm constraint can be expressed explicitly by im- posing, for any pair of aircraftiandjwhose trajectories are denoted respectively γi(ui;s)and γj(uj;s): cij(u;t ):= γi(ui;s(t))−γj(uj;s(t)) 2 2≥τ2, ∀t∈[0,t ij min ],(7) wheret ij min :=min (ti end ,tj end ),witht i end is the arrival time of aircrafti to its end point αend i . The natural spline parameter,s, and the time,t, are related by the following bijective relation betweensandt: s(t)=θ−1(ui;t)andt= θ(ui;s),(8) where θ(ui;s)= s 0 dγx(ui,ξ) dξ 2+ dγx(ui,ξ) dξ 2dξrelates time and the natural spline parameter,s, through the arc-length closed-form ex- pression.

Therefore, the problem can be formulated as: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩min u f(u)=1 N N i=1 fi(ui) s.t.c ij(u;t )≥τ2 ∀t∈[0,t ij min ];i=1,...,N−1; j=i+1,...,N.(9) In order to computec ij(u;t ), one needs to evaluate γi(ui;s(t))and γj(uj;s(t))at timet. This is achieved in practice by sampling the trajectory in time, noting that one has for eachs, the values (xi,yi,di,t), whered iis the distance travelled by aircraftiwhose co-ordinates are (xi,yi)at timet. These values are sampled with respect to time, using linear interpolations (d iandtare exact; andx iandy iare approximated via sampling). C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512507 3.2. The objective function and its derivatives In this sub-section, the objective function is deﬁned and closed- form expressions of the objective-function derivatives with respect to the optimization variables are obtained. Recall that the objective- function is expressed in (6), and that the optimization variables are meant to describe the middle control-point locations ( αi(ui)= p i mid +u iqi).

Sincefi(u)depends only on theithucomponent, theith partial derivatives of the objective-function is: ∂f ∂ui=1 N ∂fi ∂ui. As described in the previous section, each trajectory can be rep- resented by a parametric curve inR 2:γi(ui;s)=:

γi x(ui;s) γ i y(ui;s) ,where s∈[0,2] (=[0,N cp−1]), and it is deﬁned by a B-spline driven by three control points corresponding tos∈ {0,1,2 }.

To deﬁne the objective function, let us deﬁne: ˆ Ti(ui):= 2 0 d γi x(ui;s) ds 2 + d γi y(ui;s) ds 2 ds. and, Ti(ui):= 2 0 d γi x(ui;s) ds 2 + d γi y(ui;s) ds 2 ds.(10) The former gives the total distance travelled by the aircraftiwhile the latter represents an energy. Minimizing the energyT ialso minimizes ˆ T i, therefore, one can rather minimizeT iin order to avoid the square root (see, for exampledo Carmo 1992, pp. 190–200). We choose the total energy rather than the travelled distance as a criterion in order to ease the derivatives calculation.

In order to have a normalized objective functionf (u), one deﬁnes the functionsf i(ui)∈[0,1] as follows: fi(ui)=T i(ui)−T i(0) Ti(1)−T i(0).(11) A tedious calculation, detailed inAppendix A, yields the closed-form derivatives, fori=1,...,N:

∂ ∂uifi=1 Ti(1)−T i(0) 32 15u i q i x 2+ q i y 2 . 3.3. The constraint functions Inspired byConn and Gould (1987)andVisweswariah, Haring and Conn (2000), one can reformulate the constraint functionc ijas: Cij(u)= tij min 0 max {τ2−c ij(u;t );0}dt=0. Doing so, one transforms thesemi-inﬁnite constraints(each con- straints is deﬁned over a time interval [0,t ij min ]) into a single equal- ity constraint (see alsoStein 2012for a more theoretical survey of semi-inﬁnite optimization). This constraint ensures that the separa- tion norm is respected between aircraftiandj. Indeed, if, for some pair of aircraft (i,j), there is a conﬂict, then τ−c ij(u;t )>0 during some time interval, leading to a non-zero integral (constraint violation).

One can rewriteC ij(u)as: Cij(u)= ij τ2−c ij(u;t ) dt,(12) with ij:= κ[tij in κ,tij out κ ], where the union is over each time interval [t ij in κ,tij out κ ] during which the aircraftiandjare in conﬂict (violat- ing the separation constraint), and where κis the number of such >5Nm >5Nm Fig. 8.Different trajectory conﬁgurations leading at most to: zero conﬂict (left), one conﬂict (center), and two conﬂicts (right).

time intervals. Since there is only one movable control point for each trajectory (which means an aircraft cannot do more than one turn), there can only be at most κ=2 violating time intervals for each pair of aircraft (seeFig. 8). Consequently, without loss of generality, one can assume there are at most three possible conﬁgurations corresponding to κ=0,1,2 depending on the situation. A situation with κconﬂicts corresponds to κviolating intervals.

Standard optimization methods commonly require providing the objective-function and constraint derivatives. No satisfying results have been obtained for the computation of the constraint derivatives although we are working on such an improvement (and have suc- cessfully provided derivatives for the objective function). Thus, in this study, when constraint derivatives will be required by an optimiza- tion method, we shall be content with ﬁnite-difference gradients. 4. Optimization methods This section details the different optimization methods applied to the problem. First, a genetic optimization method (GA) is used to han- dle the combinatorial aspect of the problem, as this optimization ap- proach is the most used for conﬂict resolution problems. One can next take advantage of the fact that the B-spline trajectory model allows one to apply a standard local continuous optimization method such as interior-point methods. Finally, because of the relatively large num- ber of function evaluations required for the use of ﬁnite-difference approximations of the gradients, it is natural to try one of the modern derivative-free optimization method. 4.1. Genetic algorithm and problem-speciﬁc genetic operators The ﬁrst optimization method considered to solve the problem is genetic algorithms (GA;Goldberg, 1989). We emphasize that this is also a contribution of this paper since we use the formulation (9).

Our implementation of GA selects the best individuals of the pop- ulation at each iteration using a deterministic (λ,μ)-tournament se- lection which randomly selects λindividuals and keeps the μbest elements, where λ>μare user-deﬁned parameters. This step is re- peated until a new intermediate population is completed. Genetic operators (crossover, mutation or nothing) are then applied with speciﬁed user-deﬁned probabilities (noted respectivelyp c,pmand, 1−p c−p m). Ultimately, one obtains the next generation of chromo- somes. This generational process is repeated until some user-deﬁned termination condition has been reached. In this study, a maximal number of generations is imposed.

The chromosome encoding used here represents the trajectories ofNaircraft using a vector ofNreal numbers (theu i’s, seeFig. 9). 508C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512 Aircraft 1 u1 Aircraft 2 u2 ... ... Aircraft 19 u19 Aircraft 20 u20 Fig. 9.Chromosome encoding for a 20-aircraft instance. It respects thelocality principle: two individuals close to each other in the search space represent close solutions. The crossover operator is used to mix the features of two good individuals (good from the point of view of candidates for improving upon the current solution), called parents, from the previous generation. It consists of picking the most conﬂicted aircraft of each parent (more precisely, we seek for the highest value of j =iCij(u)among all aircrafti=1,...,N), say aircraftifor one parent and aircraftjfor the other, and of modifying its B-spline trajectory over [s 0,s2] using a barycentric transformation of the two parents’ corresponding movable control point (Fig. 10).

The crossover results in two children that are expected to yield local improvement.

The mutation operator, used to diversify the genes in the popu- lation in order to explore widely the search space, consists, in this application context, of choosing randomly one movable control point and to assign to it a new value chosen randomly (using a uniform distribution) in the interval (percent) [−100; 100] (seeFig. 11).

Finally, the ﬁtness, which quantiﬁes the ability of an individual to solve the problem is deﬁned as follows: f(u)=1 N N i=1 fi(ui)+ω N i=1N j=1 Cij(u),(13) where ωis a penalization parameter, weighting the importance of feasibility with respect to optimality, set empirically by the user.

This penalty function is also invoked in order to enable the use of the derivative-free optimization method BOBYQA described in Section 4.2.1. 4.2. Local optimization method This section details brieﬂy the local optimization methods applied and the reasons why they were chosen.

4.2.1. Derivative-free optimization The lack of closed-form expressions of the constraint deriva- tives naturally leads to derivative-free optimization methods (Conn, Scheinberg & Vicente, 2009). The method chosen here is Powell’s BOBYQA (Bound Optimization BY Quadratic Approximation;Powell, 2009), one of the most effective derivative-free optimization methods available. It is based on a trust-region model described in Chapter 10 ofConn et al. (2009). Note that this class of methods does not han- dle constraints directly (except for simple bounds). Consequently, the penalized objective function deﬁned in (13)isused.

For our numerical experiments, we use all the default values proposed inPowell (2009)for the various parameters deﬁning the method.

4.2.2. Local differentiable optimization In fact, a closed-form expressions of the derivatives of the ob- jective function can be obtained here (this calculation is detailed in Appendix A). One can therefore consider applying a standard local differentiable optimization method. We choose constrained interior- point methods because they are state-of-the-art methods for non- linear programming. The numerical experiments are conducted us- ing the Matlab routine fmincon (Byrd, Gilbert & Nocedal, 2000). The CROSSOVER Parent 2 Parent 2 Child 1 Child 2 CROSSOVER Fig. 10.Parent chromosomes and children obtained via a barycentric transformation (averaging the chosen bandwidth percentage) and their corresponding trajectories. MUTATION MUTATION Fig. 11.Initial and mutated chromosome and their corresponding trajectories. C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512509 gradient of the objective function is supplied to this routine. How- ever, the derivatives of the constraint functions (12) (which are not trivial to obtain) are approximated automatically by fmincon using ﬁnite differences. Current work aims at obtaining a closed form of the gradient of the constraint functions.

Here again, we use in our tests the default values provided by Matlab for the various parameters involved in the algorithm.

5. Numerical results This section presents comparative numerical results obtained with the different optimization methods proposed above. First, an aca- demic test problem, called theroundabout, is presented. Then, a more realistic problem (so-calledoperational-like test problem)isproposed.

Although rather artiﬁcial, both test problems feature some opera- tional aspects of real-life problems and thereby allows one to test the viability of our methodology.

The results we obtain with GA are used as a reference. Indeed, we expect GA, considering its ergodicity property, to be able to reach the neighborhood of any desired point of the search space within the (large) number of iterations allowed. Consequently, with a signiﬁcant number of function evaluations, we expect GA to ﬁnd a solution whose value is relatively close to the optimal value. Moreover, GA is known to be ﬁt to handle the air traﬃc conﬂict problem (Médioni et al., 1994), which makes it a good reference to compare with the results of our local optimization approach. To summarize, we shall compare the results we obtain with local optimization method with (fmincon) and without derivatives (BOBYQA), to the ones we obtain with our GA implementation.

The proposed genetic algorithm is implemented in Java. As mentioned above, the differentiable local optimization method used is the routine fmincon from the Matlab Optimization toolbox. The derivative-free optimization method BOBYQA is coded in fortran 77 and is called from Matlab via a mex interface. As a consequence of this difference in programming language, computational time cannot be considered as an objective comparison criterion. We rather rely here on comparing the number of function evaluations, which rep- resents most of the computational time spent by the optimization methods. This corresponds to common practice in black-box opti- mization where the objective and/or constraint functions are costly to evaluate. However, in order to give an idea of the order of computa- tion time involved, let us simply present the calculation time for one evaluation of both the objective and constraint functions for a same point, and on a same traﬃc situation. In both cases (Java and Matlab), we use a 2.53 GHz processor Intel Core 2 Duo on a Ubuntu 12.04 LTS operating system. Using Java, one evaluation requires 7 milliseconds while with Matlab, it needs 26 milliseconds.

Here are the parameter values used to implement the different above-mentioned algorithms:

1) Population size: 100; number of generations: 100 (hence GA will evaluate 10,000 times the objective and constraint functions) 2) Mutation probability:p m=0.3 3) Crossover probability:p c=0.6 4) Constraint penalization parameter: ω=0.01 5) Stopping criterion of the local optimization methods: ||uk+1 − u k||<10 −6 6) Approximate global optimization value tolerance for GAf −f< f:=10 −4, for feasible solutions, wheref is the best solution value found by GA. Thisf value will be used as a comparative quality criterion in the numerical results.

The local optimization methods fmincon and BOBYQA start with the best feasible point from 100 randomly-generated points from the search space ([−1,1] N). These 100 extra function evaluations are taken into account in the function-call counts for these methods Fig. 12.Roundabout test problem conﬁguration forN=6aircraft. (in a suitable parallel environment, choosing many different starting points could be relatively inexpensive).

5.1. Roundabout test problem First, a simple but diﬃcult academic test with conﬂict situations, that is widely used for air traﬃc conﬂict problems, is considered. Each instance of the roundabout problem involvesNaircraft uniformly distributed on a circle of radius 100 nautical miles. Each of theN aircraft ﬂies to the diametrically opposed point at a common speed (forNeven, each point on the circle has an outgoing and an incoming trajectory, seeFig. 12).

This problem is convenient to test conﬂict-resolution methods be- cause, due to its symmetry, there is an expected solution. Indeed, air traﬃc controller common practice would then make all the aircraft involved in the conﬂict turn to the same side. In the sequel, this obvi- ous solution is referred to as a global solution of the problem. This test problem is widely used because automatic air traﬃc conﬂict resolu- tion algorithms are not designed to take advantage of this particular symmetry.

The roundabout test problems considered here involve succes- sivelyN=2, 4, 6, 8, and 16 aircraft. Recall that the chosen start point of the local optimization methods is the best point among a 100 randomly-generated points of [−1,1] N. GA initial population will also be randomly generated in [−1,1] N.

The computational results obtained with the three optimization approaches (GA, fmincon, and BOBYQA) on this academic problem reveals that:

1) GA and BOBYQA always (forN=2,4,6,8, and 16) ﬁnds a global minimum (all aircraft turn to the same side).

2) ForN=2, fmincon always ﬁnds a global optimum.

3) ForN=4, fmincon ﬁnds a local minimum for which objective- function value vary in [f ,3×f ] depending on the chosen starting point.

4) ForN=6 fmincon ﬁnds a local minimum for which objective- function value vary in [f ,8×f ].

5) ForN=8 and 16 aircraft, a local optimum is not always found with fmincon, depending on the chosen starting point. 510C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512 Fig. 13.Operational-like test problem conﬁguration forN=6aircraft. No further detail is given as this is an academic problem that does not correspond to real-world problems.

5.2. Operational-like test problem In order to create test problems that are more realistic, an operational-like test problem can be created by introducing some perturbations in the roundabout test problem. This study is limited to conﬂicted situations involvingN=6 aircraft, as according to air traﬃc controller experience, conﬂicts involving more than four air- craft are not typical in en-route traﬃc (only one four-aircraft conﬂict happens over France each year). For each aircraft, one randomly dis- places (or not) the position of the start and end points on the circle by 20 ◦to obtain a conﬂict situation involving various crossing angles, losing thereby the special symmetry of the roundabout test problem (seeFig. 13). A randomly-generated speedv i∈[550; 650] knots (nau- tical miles per hour) is assigned to each aircrafti,i=1,...,N, during the generation of the test-cases (for the academic roundabout test problem, the speed of each aircraft wasv i=600 knots). Following this methodology, 100 different instances were generated. Each op- timization method is applied to these 100 instances in order to make an elementary statistically-valid comparative study. To compare the results of the different methods on these 100 instances, GA results are taken as a reference, as, with 10,000 func- tion evaluations, it is expected to ﬁnd a (nearly) globally optimal value f for each of these instances. Several criteria are taken into account.

First, the percentage of success of the local optimization method is reported (whether it converges to a locally-optimal solution or whether it fails). The objective-value statistics (mean, median, worst, best) are taken over the 100 instance runs. For each instance, the dif- ference in the objective-function values is calculated with respect to f when the local method converges to a feasible solution. The best (respectively worst) objective-function value differences with respect tof presented in the table correspond to instances for which fmin- con and BOBYQA performed the best (respectively worst). Finally, the statistics on the number of function evaluations required to converge are displayed. The results are summarized inTable 1.

The percentage of success for the local optimization methods is high. This shows that local optimization methods are a viable alternative for automatic conﬂict resolution. The objective-function difference reveals that a local method can even ﬁnd better feasible objective function values than genetic algorithms. More precisely, fmincon obtains better results (in terms of the objective-function value) than GA in 48 instances (over the 95 solutions fmincon ﬁnds), which is very promising. BOBYQA obtains better results (again, in terms of the objective-function value) than GA for 23 instances, which is still signiﬁcant. A discussion on the corresponding num- ber of function evaluations is given below. Of course, in the remain- ing cases, the local methods converged to high (i.e. bad) objective- function values (local minima) when compared with GA, which is why there is a discrepancy between the median and the average of the objective-function value difference (a gap of 0.30 can be considered as signiﬁcantly bad).

As mentioned previously, GA performs 10,000 function evaluations here, although it ﬁnds a feasible (conﬂict-free) solution before. After how many function evaluations does GA ﬁnd such a feasible solution? The ﬁrst conﬂict-free solution is discovered by GA quickly (110 evaluations on average, roughly equivalent to the pure random-search starting-point strategy of fmincon and BOBYQA).

However, such feasible solutions are of poor quality. In any case, the global optimum found by GA requires much more time to be found with an fprecision (6 times the average number of function evaluations needed by fmincon, and 15 times for BOBYQA). The lo- cal optimization methods are therefore more eﬃcient than GA for ﬁnding a feasible solution in terms of computational time for compa- rable quality. We remark however that GA seems to be more sta- ble than fmincon, as it ﬁnds feasible solutions for 98 out of the 100 test cases, and the two cases where GA does not ﬁnd a fea- sible solution, the GA solution is close to being feasible (the con- straint violation is of order 10 −2, when, for fmincon infeasible so- lution, the violation is signiﬁcantly larger). However, BOBYQA is the most stable methods among these three as italwaysﬁnd a feasible solution. Table 1 Comparative numerical results over the 100 instances. GA fmincon BOBYQA Feasibility success rate (percent) 98 95100 Average objective-function value0.0440.049 0.091 Median objective-function value0.0440.070 0.11 Objective-function value difference with respect tof Mean−0.0290.067 Median−0.0060.047 Worst−0.340.29 Best−−0.054−0.069 Number of function evaluations Mean 10000 1135400 Median 10000 820393 Worst 10000 68251033 Best 10000178194 C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512511 Finally, note that the number of functions reported for fmincon is strongly impacted by the fact that it relies on ﬁnite-difference con- straint gradients (inNdimensions, each ﬁnite-difference gradient approximation requiresN+1 function evaluations, and gradients are needed at least once at every iteration). Indeed, the average num- ber of iterations is only 54. Once the constraint gradients will be computed (current research work), one can expect a much lower number of function evaluations for local differentiable optimization methods. 6. Conclusion Standard approaches to air traﬃc conﬂict resolution problems rely on discretization of the search space yielding a combinatorial optimization formulation of the problem and, thereby, computation- intensive optimization methods must be used. This paper shows that the use of an original B-spline trajectory model involving only one continuous variable per aircraft, and of a semi-inﬁnite programming formulation of the separation constraints, permits one to obtain good results via local optimization methods. A major advantage of a local optimization methods is that it requires much fewer (costly) function evaluations. Moreover, in a safety operational context, competent authorities may prefer decision-aid tools based on a deterministic optimization methods with an underlying convergence theory and a more eﬃcient use of function evaluations.

Note that the optimization formulation introduced here involves only differentiable functions. Also, one expects that, ultimately, it will not be necessary to rely either on derivative-free optimiza- tion methods, or on ﬁnite-difference variants of classical differen- tiable local optimization methods. Current work concentrates on computing closed-form constraint derivatives in order to decrease signiﬁcantly the computational time required by local optimization methods. Finally, the elementary random-search start-point selec- tion step combined here with the local optimization method could be replaced with a more adapted hybrid two-phase (global/local) method. Acknowledgments This work has been supported by French National Research Agency (ANR) through COSINUS program (project ID4CS no. ANR-09-COSI- 005) and through JCJC program (project ATOMIC no. ANR 12-JS02- 009-01). The authors are indebted to Professor Stéphane Puechmorel for his help on the calculation and computation of the derivatives, and to Laurent Lapasset for interesting discussions on the operational aspects of the problem. The authors are grateful to two anonymous reviewers for many detailed and helpful comments. Appendix A. Objective-function derivatives Recall that the expression for the objective function is given by (6), (10)and(11), and that given a control point list, αi start ,αi(ui) and αi end ,thexcomponents of the B-spline approximation associated with movable control point α(ui)for aircraftiis given by γi x(ui;s).

For notational simplicity, in the sequel, the aircraft indexiis dropped in the notation, except on the optimization variableu i,withrespect to which we differentiate. Moreover, the factor1 Ti(1)−Ti(0), being a normalization constant, it is not included in the following calculation.

Each term in the integral (in the expression ofT i(ui)) can be considered separately. The calculation is only detailed for the term involving γi x.

Using the Taylor’s expansion expressed in Eq. (2), one can obtain γi x(ui;s)for eachsin the interval [s k,sk+1 ]. This Taylor’s expansion is exact as the B-spline curve γi x(ui;s)is a third-degree polynomial.

The trajectory model involves three control points and two phan- tom control points. Thus, the B-spline is to be expressed on two intervals, [s 0,s1]and[s 1,s2]. Let us begin with the case for whichs∈[s 0,s1]:=[0,1]. Using (2) withk=0, one obtains: γi x(ui;s)=γ0 x+(s−s 0)dγ0 x ds+ (s−s 0)2 2d 2γ0 x ds2 +(s−s 0)3 6d 3γ0 x ds3. To calculate our B-spline fors∈[0,2], ﬁve control points are used: the phantom pointX −1 :=2X 0−X 1=2X 0−(pi mid x+u iqi x), the starting pointX 0:= αstart x, the movable control pointX 1:= αi x(ui)=p i mid x+ u iqi x, the ending pointX 2:= αend xand the second phantom point X 3:=2 αend x−(pi mid x+u iqi x). To compute on each interval [s k,sk+1 ], the B-spline method uses four control pointsX k−1 ,Xk,Xk+1 ,Xk+2 .In our case, whens∈[s 0,s1]=[0,1] (i.e. fork=0), the B-spline com- putation usesX −1,X0,X1,andX 2,andwhens∈[s 1,s2]=[1,2] (i.e.

fork=1), the computation usesX 0,X1,X2,andX 3. The dependen- cies inu iare determined accordingly. In the following, for notational simplicity, we will useX 0for αstart x,andX 2for αend x.Thus,forall s∈[s 0,s1]=[0,1], one has: γi x(ui;s)=X 0+ p i mid x+u iqi x−X 0 s + X 2−2p i mid x−2u iqi x+X 0 s3 6. Then, one can obviously deduce, by differentiating the expression above, that, for alls∈[s 0,s1]: dγi x(ui;s) ds= p i mid x+u iqi x−X 0 + X 2−2p i mid x−2u iqi x+X 0 s2 2. Consequently, for alls∈[s 0,s1]: d γi x(ui;s) ds 2 = p i mid x+u iqi x−X 0 2 +s 4 4 X 2−2p i mid x−2u iqi x+X 0 2 +s 2 pi mid x+u iqi x−X 0 × X 2−2p i mid x−2u iqi x+X 0 . It is necessary to calculate the integral smax0 dγi x(ui;s) ds 2+ dγi y(ui;s) ds 2ds in order to obtainT i. As mentioned before, only the terms in γi xare detailed, and, in this part, the calculation focuses on the interval where s∈[0,1]: 1 0 d γi x(ui;s) ds 2 ds= 1 0 p i mid x+u iqi x−X 0 2 +s 4 4 X 2−2p i mid x−2u iqi x+X 0 2 +s 2 pi mid x+u iqi x−X 0 × X 2−2p i mid x−2u iqi x+X 0 ds = p i mid x+u iqi x−X 0 2 +15 X 2 2−p i mid x−u iqi x+X 0 2 2 +13 p i mid x+u iqi x−X 0 × X 2−2p i mid x−2u iqi x+X 0 . Then, differentiating the above expression with respect tou i, one obtains the closed-form expression of the integral for s∈[s 0,s1]:=[0,1] : ∂ ∂ui 1 0 d γi x(ui;s) ds 2 ds=q i x 16 15 p i mid x+u iqi x −6 5X 0+215X 2 . 512C. Peyronne et al. / European Journal of Operational Research 241 (2015) 502–512 The calculation now concentrates on the interval wheres∈ [s 1,s2](k=1). Using Taylor’s expansion (2) on this interval, one obtains: γi x(u;s )=γ1 x(s1)+(s−s 1)dγ1 x(s1) ds+ (s−s 1)2 2d 2γ1 x(s1) ds2 +(s−s 1)3 6d 3γ1 x(s1) ds3 . Now, this value has to be calculated fors∈[s 1,s2]=[1,2], where the optimized control point isX k=p mid x+u iqx. Consequently, due to the change ofkfrom 0 to 1 for this interval, one has:X k−1 :=X 0=αstart x; X k+1 :=X 2=αend x;Xk+2 :=2X 2−p mid x−u iqx. Thus, we obtain for all s∈[1,2]: γi x(ui;s)=1 6 X 2+4p i mid x+4u iqi x+X 0 + (s−1 ) 2 (X2−X 0) +(s−1 )2 2 X 2−2p i mid x−2u iqi x+X 0 + (s−1 )3 6 2p i mid x+2u iqi x−X 2−X 0 . Differentiating with respect tos, one obtains for alls∈[1,2]:

dγi x(ui;s) ds=1 2(X2−X 0)+(s−1 ) X2−2p i mid x−2u iqi x+X 0 + (s−1 )2 2 2p i mid x+2u iqi x−X 2−X 0 ,(A.1) and therefore, d γi x(ui;s) ds 2 =1 4(X2−X 0)2+(s−1 )2 X2−2p i mid x−2u iqi x+X 0 2 +(s−1 )4 4 2p i mid x+2u iqi x−X 2−X 0 2 +(s−1 )(X2−X 0) X2−2p i mid x−2u iqi x+X 0 + (s−1 )2 2 (X2−X 0) 2p i mid x+2u iqi x−X 2−X 0 + (s−1 )3 2p i mid x+2u iqi x−X 2−X 0 × X 2−2p i mid x−2u iqi x+X 0 . Integrating this value over [1,2] yields: 2 1 dγi x(ui;s) ds 2ds=1 4(X2−X 0)2+13 X 2−2p i mid x−2u iqi x+X 0 2 +120 2p i mid x+2u iqi x−X 2−X 0 2 +12(X2−X 0) X2−2p i mid x−2u iqi x+X 0 +1 6(X2−X 0) 2p i mid x+2u iqi x−X 2−X 0 +1 4 2p i mid x+2u iqi x−X 2−X 0 × X 2−2p i mid x−2u iqi x+X 0 . One then differentiates with respect tou i: ∂ ∂ui 2 1 d γi x(ui;s) ds 2 ds=q i x 16 15 (pi mid x+u iqi x)+215X 0−65X 2 . As mentioned earlier, similar formulas can be obtained for the terms in γi y. Thus, the exact formula for the objective-function deriva- tive is written as follows: ∂ ∂uifi(ui)=q i x 16 15 p i mid x+u iqi x −6 5X 0+215X 2 +q i x 16 15 p i mid x+u iqi x +2 15X 0−65X 2 +q i y 16 15 p i mid y+u iqi y −6 5Y 0+215Y 2 +q i y 16 15 p i mid y+u iqi y +2 15Y 0−65Y 2 =32 15q i x pi mid x+u iqi x +32 15q i y(pi mid y+u iqi y) −1615q i x(X2+X 0)−1615q i y(Y2+Y 0) Finally, since by construction, we havep i mid x=X2+X02 ,andp i mid y= Y2+Y02 ,theith component of the gradient vector ∇if(u)is: ∂ ∂uif(u)=1 N(Ti(1)−T i(0)) 32 15u i q i x 2+ q i y 2 . References Alliot, J. M., Bosc, J. F., Durand, N., & Maugis, L. (1997). CATS: A complete air traﬃc sim- ulator. In:IEEE Digital Avionics Systems Conference, American Institute of Aeronautics and Astrophysics(Vol. 2, pp. 8.2-30–8.2-37). IEEE: USA.

Alonso-Ayuso, A., Escudero, L. F., & Martin-Campo, F. J. (2011). Collision avoidance in the air traﬃc management: A mixed integer linear optimization approach.IEEE Transactions on Intelligent Transportation Systems,12, 47–57.

Byrd, R. H., Gilbert, J. C., & Nocedal, J. (2000). A trust region method based on inte- rior point techniques for nonlinear programming.Mathematical Programming,89, 149–185.

Conn, A. R., & Gould, N. I. M. (1987). An exact penalty function for semi-inﬁnite pro- gramming.Mathematical Programming,37, 19–40.

Conn, A. R., Scheinberg, K., & Vicente, L. N. (2009).Introduction to derivative-free opti- mization. MPS-SIAM series on optimization. Society for Industrial and Applied Math- ematics Philadelphia, PA: USA. 173–205.

Delahaye, D., & Puechmorel, S. (2013).Modeling and optimization of air traﬃc.Wiley- ISTE: USA.

Dimarogonas, D. V., & Kyriakopoulos, K. J. (2005). A feedback stabilization and collision avoidance scheme for multiple independent nonholonomic non-point agents. In:

Proceedings of the IEEE international symposium on robotics and automation(pp. 820– 825). IEEE: USA.

do Carmo, M. P. (1992).Riemannian geometry. Birkhäuser: Switzerland.

Dougui, N., Delahaye, D., Puechmorel, S., & Mongeau, M. (2013). A light-propagation model for aircraft trajectory planning.Journal of Global Optimization,56(3), 873–895.

Duncan, M. (2005).Applied geometry for computer graphics and computer-aided design (2nd ed.) Springer: UK.

Durand, N., & Alliot, J. M. (2009). Ant colony optimization for air traﬃc conﬂict res- olution. In:Proceedings of the 8th USA/Europe air traﬃc management research and development seminar. Napa, USA.

Durand, N., Alliot, J. M., & Noailles, J. (1996). Automatic aircraft conﬂict resolution using genetic algorithms. In:11th annual ACM symposium on applied computing (ACM/SAC 96).ACM:USA.

Goldberg, D. E. (1989). Genetic algorithms in search optimization and machine learning.

Addison-Wesley: USA.

Kuchar, J. K., & Yang, L. C. (2000). A review of conﬂict detection and resolution modeling methods.IEEE Transactions on Intelligent Transportation Systems,1, 179–189.

Médioni, F., Durand, N., & Alliot, J. M. (1995). Air traﬃc conﬂict resolution by ge- netic algorithms. Lecture Notes in Computer Science, Vol. 1063.Artiﬁcial evolution (pp. 370–383). Springer: USA.

Milam, M. B., Mushambi, K., & Murray, R. M. (2000). A new computational approach to real-time trajectory generation for constrained mechanical systems. In:Proceedings of the 39th IEEE conference on decision and control. IEEE: USA.

Olive, X. (2006). Résolution de conﬂits par algorithmes stochastiques parallèles (Con- ﬂict resolution via parallel stochastic algorithms) (Ph.D. thesis). Toulouse, France:

École Nationale Supérieure de l’ Aéronautique et de l’ Espace.

Pallottino, L., Feron, E., & Bicchi, A. (2002). Conﬂict resolution problems for air traﬃc management systems solved with mixed integer programming.IEEE Transactions on Intelligent Transportation Systems,3, 3–11.

Powell, M. J. D. (2009). The Bound Optimization BY Quadratic Approximation (BOBYQA) algorithm for bound constrained optimization without derivatives. Technical Report NA2009/06. Cambridge,England.

Roussos, G., & Kyriakopoulos, K. J. (2009). Towards constant velocity navigation and collision avoidance for autonomous nonholonomic aircraft-like vehicles. In:Pro- ceedings of the Joint 48th IEEE conference on decision and control and 28th Chinese control conference. IEEE: USA.

SESAR Joint Undertaking (2009). SESAR Joint Undertaking Brochure: Today’s partner for tomorrow’s aviation. Technical Report. Brussels, Belgium.

Stein, O. (2012). How to solve a semi-inﬁnite optimization problem.European Journal of Operational Research 233, 312–320.

Visweswariah, C., Haring, R., & Conn, A. R. (2000). Noise considerations in circuit op- timization.IEEE Transactions on Computer Aided Design of Integrated Circuits and Systems,19, 679–690.