J Zhou
email: junzhou@recherche.enac.fr.

S Cafieri
email: sonia.cafieri@enac.fr.c

D Delahaye
email: daniel.delahaye@enac.fr.

M Sbihi

Optimization-Based Design of Departure and Arrival Routes in Terminal Maneuvering Area

The efficient design of departure and arrival routes in the airspace surrounding airports, called the Terminal Maneuvering Area (TMA), is crucial for increasing the capacity of such areas, and thus alleviating congestion around airports caused by worldwide air traffic growth. In this paper, an efficient method of designing departure and arrival routes in TMA is proposed, taking into account the configuration of the airport and nearby environment, as well as related operational constraints, such as obstacle avoidance and route separation. Each route is modeled in three dimensions (3D), and consists of two components: a curve in the horizontal plane and a cone in the vertical plane. A Branch and Bound (B&B)-based approach is developed, whose branching strategies are tailored to how the obstacles are avoided. Routes are generated sequentially, and each route is initially built in isolation. If the route is found to be in conflict with previously generated routes, it is perturbed locally around the conflict zones.

Numerical tests, performed on artificially generated problems and the TMA of Paris Charles-de-Gaulle (CDG) airport, demonstrate that the proposed method is efficient and could be embedded in a decision-aid tool for procedure design.

I. Introduction

A Terminal Maneuvering Area (TMA) is an area surrounding one or more neighboring airports, that is designed to handle aircraft arriving at and departing from such airports. Most airports

have pre-designed routes to guide departures and arrivals, called Standard Instrument Departure (SID) routes and Standard Terminal Arrival Routes (STAR). A SID is a flight route followed by an aircraft from its take-off phase to the en-route phase. A STAR is a route which connects the last en-route waypoint to the Initial Approach Fix (IAF). The route segment connecting the IAF to the Final Approach Fix (FAF) is called the Initial Approach. After the FAF, the flight path is usually specified by ground-based instruments. In the present work, a STAR refers to the section of the route connecting a TMA entry point to a FAF. Currently, SIDs and STARs are designed manually according to operational requirements [START_REF]Aircraft Operations (Doc 8168)[END_REF] and airport layout. This kind of design is generally not very efficient and is not designed to optimize any specific objective. Optimizing departure and arrival procedures in TMA is crucial for regulating air traffic flows, whose continual growth [START_REF] Eurocontrol | EUROCONTROL Long-Term Forecast: IFR Flight Movements 2010-2030[END_REF] can affect air traffic operations. The aim of this work is to automatically design SIDs and STARs in three dimensions (3D), with the objective of minimizing the total route length while satisfying various constraints, including obstacle avoidance and route separation. The routes are designed to be suitable for a range of aircraft, since the specifications and capabilities of aircraft using the routes can vary significantly. The proposed design is at a strategic level, thus only static obstacles are taken into account.

Performance Based Navigation (PBN) is one of the key navigation concepts to enable future airspace and traffic flow design [START_REF]Performance Based Navigation (PBN) Operational Approval Handbook[END_REF]. It offers operational benefits such as enhanced safety and increased efficiency. Required Navigation Performance (RNP), a typical way of navigating within the PBN concept, is especially useful in complex airspaces such as TMAs, since it provides higher design flexibility. Currently the RNP-1, a type of RNP with a specified performance level, enables the Radius-to-Fix (RF) functionality, which can be applied in SID/STAR design. A RF leg is defined as an arc with specified radius between two defined waypoints [START_REF] Eurocontrol | European Airspace Concept Handbook for PBN Implementation[END_REF]. Route construction in our proposed method is based on the RNP concept, in particular the RF functionality in SIDs/STARs.

The problem considered fits into the framework of path planning. More specifically, it is a route design problem. In contrast to trajectory design, the routes to be designed are not associated with any notion of time. The problem of path planning has been studied since the 1980s, especially in the domain of robotics, e.g., [START_REF] Gallina | A Technique to Analytically Formulate and to Solve the 2-Dimensional Constrained Trajectory Planning Problem for a Mobile Robot[END_REF][START_REF] Souissi | Path planning: A 2013 survey[END_REF]. Recently, optimal flight path design has received increased attention, and a range of solution approaches have been proposed (nicely summarized in [START_REF] Delahaye | Mathematical Models for Aircraft Trajectory Design: A Survey[END_REF]). In the following, several methods for optimal route design proposed in the literature are reviewed, according to different points of view: designing in 2D vs. in 3D, one route design vs. multiple routes design, exact design approach vs. heuristic design approach.

For 2D route design, a typical approach for computing the shortest path between two points that avoids a fixed set of polygonal obstacles, is to use a Visibility Graph, introduced by Lozano-Pérez and Wesley in 1979 [START_REF] Lozano-Pérez | An Algorithm for Planning Collision-free Paths Among Polyhedral Obstacles[END_REF]. It has been proven that the 2D shortest path avoiding polygonal obstacles is the shortest path on the corresponding Visibility Graph [START_REF] Lozano-Pérez | An Algorithm for Planning Collision-free Paths Among Polyhedral Obstacles[END_REF]. In the case when both polygonal and curved obstacles are present, the notion of a Tangent Graph is introduced by Liu and Arimoto in 1992 [START_REF] Liu | Path Planning Using a Tangent Graph for Mobile Robots Among Polygonal and Curved Obstacles[END_REF]. A Tangent Graph allows one to take into account obstacles in the form of circular discs, such as in [START_REF] Pocchiola | Minimal tangent visibility graphs[END_REF][START_REF] Kim | Shortest paths for disc obstacles[END_REF]. In [START_REF] Kim | Shortest paths for disc obstacles[END_REF], the authors prove that the shortest path lies in a convex hull of a few circular obstacles around the line segment connecting the start and end points. This significantly reduces the size of the search space. In addition to graph-based methods, several other approaches have been proposed for solving 2D optimal route design problems. For example, a Branch and Bound (B&B) method that builds an optimal route avoiding circular obstacles is proposed in [START_REF] Eele | Path planning with avoidance using nonlinear branch and bound optimization[END_REF]. The branching strategy is based on the decision to avoid an obstacle by turning clockwise or counter-clockwise. In [START_REF] Pierre | Aircraft Trajectory Planning with Dynamical Obstacles by Artificial Evolution and Convex Hull Generations[END_REF], the path planning problem in a TMA is solved using convex hull generation combined with a Genetic Algorithm (GA). Again, the obstacles are avoided by turning counter-clockwise or clockwise around their boundaries. In [START_REF] Jimmy | Turn-Constrained Route Planning for Avoiding Hazardous Weather[END_REF], a dynamic programming approach, based on a variant of the Bellman-Ford method, is applied to compute an optimal turn-constrained route in the presence of weather hazards. The generated route is in a polygonal form in 2D.

Optimal route design in 3D is more difficult than in 2D, since the vertical profile has to be taken into account. One approach proposed in the literature is to extend the 2D visibility graph to 3D with polyhedral obstacles [START_REF] Sharir | On Shortest Paths in Polyhedral Spaces[END_REF][START_REF] Jiang | Finding the 3D shortest path with visibility graph and minimum potential energy[END_REF]. However, the problem becomes NP-hard in 3D [START_REF] Canny | New lower bound techniques for robot motion planning problems[END_REF], primarily because the shortest path does not necessarily pass through the vertices of polyhedrons. Other works have contributed to the 3D-route design in the aeronautical domain. In [START_REF] Gianazza | Allocating 3D-trajectories to air traffic flows using A* and genetic algorithms[END_REF], two different approaches are considered for designing routes connecting a pair of airports: an A* algorithm and a GA. The horizontal route is only allowed to vary at the boundary of a TMA, thus it is polygonal in shape. The associated vertical profile is specified by connecting fixed level flights using linear climb or descent paths. In [START_REF] Pfeil | Optimization of airport terminal-area air traffic operations under uncertain weather conditions[END_REF], the author designs airport terminal routes that avoid obstacles using a modified A* algorithm. Other operational constraints are also considered (e.g., the need for an aircraft to perform smooth turns). The vertical profile is in the form of a cone, bounded by the minimum and maximum aircraft slopes. In [START_REF] Zhou | Optimization of Arrival and Departure Routes in Terminal Maneuvering Area[END_REF], the departure and arrival routes are designed using the Fast Marching Method with an obstacle avoidance constraint. However, the proposed algorithm does not specifically deal with the take-off or descent slopes of aircraft. In [START_REF] Zhou | Optimizing the design of a route in Terminal Maneuvering Area using Branch and Bound[END_REF], a B&B method is used to design a 3D-route avoiding obstacles. The branching strategy is tailored to the way the obstacles are avoided, namely turning clockwise or counter-clockwise, or imposing a level flight below an obstacle.

The majority of the works cited above concern the design of a single route ( [START_REF] Lozano-Pérez | An Algorithm for Planning Collision-free Paths Among Polyhedral Obstacles[END_REF][START_REF] Liu | Path Planning Using a Tangent Graph for Mobile Robots Among Polygonal and Curved Obstacles[END_REF][START_REF] Pocchiola | Minimal tangent visibility graphs[END_REF][START_REF] Kim | Shortest paths for disc obstacles[END_REF][START_REF] Eele | Path planning with avoidance using nonlinear branch and bound optimization[END_REF][START_REF] Pierre | Aircraft Trajectory Planning with Dynamical Obstacles by Artificial Evolution and Convex Hull Generations[END_REF][START_REF] Jimmy | Turn-Constrained Route Planning for Avoiding Hazardous Weather[END_REF][START_REF] Sharir | On Shortest Paths in Polyhedral Spaces[END_REF][START_REF] Jiang | Finding the 3D shortest path with visibility graph and minimum potential energy[END_REF][START_REF] Canny | New lower bound techniques for robot motion planning problems[END_REF][START_REF] Zhou | Optimizing the design of a route in Terminal Maneuvering Area using Branch and Bound[END_REF]). The main difficulty when designing multiple routes is ensuring that distinct routes meet separation constraints, which is particularly important from the Air Traffic Management (ATM) point of view. In order to design multiple routes satisfying the separation criterion, some authors apply a sequential 1-againstn strategy where the routes are generated one after another according to a user-defined priority order (for example according to the traffic load on each route). The previously constructed routes become obstacles for the route under construction. For the sequential strategy, the quality of a solution depends on the order in which routes are generated. In [START_REF] Gianazza | Allocating 3D-trajectories to air traffic flows using A* and genetic algorithms[END_REF] and following works [START_REF] Gianazza | Separating air traffic flows by allocating 3D-trajectories[END_REF][START_REF] Gianazza | Assessment of the 3D-separation of Air Traffic Flows[END_REF], routes are generated sequentially according to traffic load, and the route separation is realized through an A* algorithm. The modified A* method proposed in [START_REF] Pfeil | Optimization of airport terminal-area air traffic operations under uncertain weather conditions[END_REF] deals with route separation in a similar way.

In [START_REF] Zhou | Optimization of Arrival and Departure Routes in Terminal Maneuvering Area[END_REF], a Simulated Annealing (SA) method is applied to optimize the order in which the routes are generated. In [START_REF] Polishchuk | Generating Arrival Routes with Radius-to-Fix Functionalities[END_REF], the design of STARs based on the RF functionality is considered, where routes are generated sequentially and merged gradually to the runway. Some other works apply a global strategy, where all routes are generated simultaneously in order to minimize a global cost associated with the set of routes. For example in [START_REF] Gianazza | Allocating 3D-trajectories to air traffic flows using A* and genetic algorithms[END_REF], a GA is applied to generate routes simultaneously and the route sections involved in conflicts are penalized in the global cost. Similarly, in [START_REF] Zhou | Optimal Design of SIDs/STARs in TMA Using Simulated Annealing[END_REF], a SA algorithm is used to simultaneously generate multiple 3D-routes in TMA.

As the design of routes is generally achieved by minimizing or maximizing a given objective, the problem is naturally modeled as an optimization problem. Exact or heuristic methods can be applied to find a solution. An exact approach guarantees an optimal solution, if such a solution exists, but may be time-consuming for difficult problems. Two typical exact approaches based on graph models, Dijkstra's algorithm [START_REF] Dijkstra | A Note on Two Problems in Connexion with Graphs[END_REF] and Bellman-Ford algorithm [START_REF] Bellman | On a Routing Problem[END_REF], are applied in the previously cited works [START_REF] Kim | Shortest paths for disc obstacles[END_REF][START_REF] Jimmy | Turn-Constrained Route Planning for Avoiding Hazardous Weather[END_REF]. Another widely used exact approach is the B&B method [START_REF] Lawler | Branch-And-Bound Methods: A Survey[END_REF], which is used to solve discrete and combinatorial optimization problems. Two applications of the B&B method in the considered context are presented in [START_REF] Eele | Path planning with avoidance using nonlinear branch and bound optimization[END_REF][START_REF] Zhou | Optimizing the design of a route in Terminal Maneuvering Area using Branch and Bound[END_REF]. These exact methods are efficient for the design of an individual route. However, the optimal design of multiple 3D-routes satisfying the numerous constraints of TMAs (route separation, obstacle avoidance, noise abatement, etc.) is a much more complex problem. One of the main difficulties is finding suitable analytic mathematical expressions for some of the problem constraints, making exact methods difficult to apply. In contrast to exact approaches, a heuristic approach seeks to produce good-quality but not necessarily optimal solutions in reasonable computation time. Thus, heuristics are more commonly used for solving 3D or multiple route design problems. The A* algorithm [START_REF] Hart | A Formal Basis for the Heuristic Determination of Minimum Cost Paths[END_REF] is a widely used heuristic approach in graph-based shortest path design problems. Several previously mentioned works apply the A* (or modified A*) method, such as [START_REF] Liu | Path Planning Using a Tangent Graph for Mobile Robots Among Polygonal and Curved Obstacles[END_REF][START_REF] Gianazza | Allocating 3D-trajectories to air traffic flows using A* and genetic algorithms[END_REF][START_REF] Pfeil | Optimization of airport terminal-area air traffic operations under uncertain weather conditions[END_REF]. Meta-heuristics [START_REF] Dreo | Metaheuristics for hard optimization[END_REF], including GA and SA, are also often applied in the ATM domain, such as in [START_REF] Pierre | Aircraft Trajectory Planning with Dynamical Obstacles by Artificial Evolution and Convex Hull Generations[END_REF][START_REF] Gianazza | Allocating 3D-trajectories to air traffic flows using A* and genetic algorithms[END_REF][START_REF] Zhou | Optimization of Arrival and Departure Routes in Terminal Maneuvering Area[END_REF][START_REF] Zhou | Optimal Design of SIDs/STARs in TMA Using Simulated Annealing[END_REF].

In the present paper, a method based on an exact approach is proposed to solve the problem of designing multiple routes in a TMA. In contrast to previous work, where routes are usually considered in 2D, our methodology deals with routes in 3D by taking into account the specific aircraft take-off and landing slopes. Moreover, the most commonly used method to separate two routes in the literature is to modify routes in the horizontal plane. The method proposed in this paper extends previous approaches by also allowing modifications in the vertical plane. The obtained routes are in compliance with the RNP concept, and can be regarded as preliminary guidelines for manual design. In the proposed methodology, obstacle avoidance and route separation are considered as the primary constraints. An obstacle is modeled as a cylinder in 3D. Each 3D-route consists of a horizontal curve associated with a vertical cone. The horizontal curve is composed of a series of arcs, used to avoid obstacles, and straight line segments that tangentially connect two arcs. In the vertical plane, for each route a cone is designed and the possibility of imposing a level flight below an obstacle is allowed. This is a potentially effective way of avoiding other flight paths, as it enriches the space of possible maneuvers and corresponds to what is done in reality (for instance, it is used on one STAR in Paris Roissy airport, see Sec .V). The problem of designing multiple SIDs/STARs is modeled as an integer optimization problem. A simpler subproblem of designing a single route avoiding obstacles is solved by B&B [START_REF] Zhou | Optimizing the design of a route in Terminal Maneuvering Area using Branch and Bound[END_REF], where the branching strategies correspond to the different ways of avoiding an obstacle. To address the design of multiple routes, a separation constraint (3NM in the horizontal plane, or 1000ft in the vertical plane) between pairwise routes must also be satisfied. Two routes are defined as in conflict when a loss of separation occurs. A conflict detection technique is proposed to check pairwise conflicting routes. Moreover, a B&B-based method is developed to sequentially build minimal-length routes that satisfy obstacle avoidance and route separation constraints. This paper is organized as follows. Section II introduces the routes and obstacles modeling.

Section III explains the approach used for generating a single route. Section IV presents the method for solving the multiple route design problem. Section V presents some simulation results and, finally, Sec. VI draws conclusions.

II. Problem statement and mathematical model

In this study, the SIDs/STARs design in the TMA surrounding a single airport is considered.

The input data are presented in Tables 1 and2. In the following, the route and obstacle modeling as well as the constraints taken into account are presented.

Table 1 Input data related to airport configuration the total number of obstacles to be avoided, M ∈ N the total number of the available runways the usage of each runway (take-off or landing) the direction of take-off or landing on each runway the coordinates and altitudes of thresholds and center for each runway the coordinate and altitude of the corresponding FAF for a runway used for landing The obstacles (together with their protection areas) are modeled as cylinders in 3D, whose projections onto the horizontal plane are disks. Let M ∈ N be the number of cylinders. Each cylinder Ω j , j = 1, • • • , M is defined by (C j (x j , y j ), r j , z j inf , z jsup ), where C j (x j , y j ) is the center and r j is the radius of the top and the bottom of the cylinder; z j inf and z jsup are the altitude of the lower and upper bases respectively.

Let N ∈ N be the total number of routes to be built. A 3D-route γ i , i = 1, • • • , N is defined by two elements: a curve γ i H in the horizontal plane and a cone γ i V in the vertical plane. The horizontal curve γ i H is designed in the form of a series of segments and arcs, motivated by the fact that the shortest path among circular obstacles in the horizontal plane consists of segments (tangentially connecting two obstacles) and arcs (lying on the border of obstacles) [START_REF] Liu | Path Planning Using a Tangent Graph for Mobile Robots Among Polygonal and Curved Obstacles[END_REF]. An example is shown in Fig. 1a). The form of the horizontal curve both simplifies the route modeling, and corresponds to the shape of SIDs/STARs under RNP. In particular, a segment defined by two tangent points corresponds to a standard point-to-point leg in ATM and an arc defined by start and end tangent points, turn center and radius, corresponds to a RF leg in ATM. The vertical cone γ i V is designed to contain all ascent (or descent) profiles of the aircraft flying on this route. The vertical of each segment on the function corresponding to the lower bound is either α min (on a climbing or descending route section) or 0 (on a route section where the flight level is maintained). Similarly, the slope of each segment on the function corresponding to the upper bound is either α max or 0.

By abuse of language, this band is called a cone. The idea of taking a cone that contains all vertical profiles is inspired by the behavior illustrated in Fig. 2, which shows some real take-off and landing data at Paris Charles-de-Gaulle (CDG) airport. From the figure it can be seen clearly that in both cases the vertical profiles are contained in a region defined by two straight lines. The variation in vertical profile is mainly due to different aircraft mass and performance, and to the effect of the wind.

In the horizontal plane, γ i H connects a start point A i (x Ai , y Ai ) to an end point B i (x Bi , y Bi ). In a SID case, the start point is at the midpoint of a runway threshold and the end point is an exit point of a TMA. In a STAR case, for the simplicity of implementation, the route is built by considering a FAF as start point and an entry point of TMA as end point. The FAF is usually aligned to the central line of the corresponding runway, and aircraft follow specific guidance (Instrument Landing System (ILS) as example) from the FAF to the runway threshold, thus the route section between FAF and the runway threshold is not considered in this study. The horizontal route γ i H is a smooth map defined as:

γ i H : [0, 1] → R 2 , (1)
where γ i H (0) = (x Ai , y Ai ) and γ i H (1) = (x Bi , y Bi ). In the vertical plane, the start point

A i (x Ai , y Ai )
is associated with an altitude H Ai , corresponding to the airport elevation in a SID case, and to the altitude of the FAF in a STAR case. The vertical profile γ i V is defined as:

γ i V : [0, 1] → I R t → [h i inf (d i (t)), h isup (d i (t))] (2)
where I R defines the set of intervals of R,

d i (t) = t 0 γ i H (s) 2 ds is the flown distance until t in the horizontal plane, [h i inf (d i ), h isup (d i )]
is the interval defined by the cross section of the cone at d i , and

γ i V (0) = [H Ai , H Ai ]. Figure 1 illustrates how γ i H is associated with γ i V in the case of a
SID, where H Ai = 0 and α min,T O (respectively, α max,T O ) is the minimum (respectively, maximum) take-off slope of aircraft on this route.

B. Decision variables

An obstacle is defined as active when it is touched by a route and must be avoided according to one of the following maneuvers: turn counter-clockwise, turn clockwise or impose a level flight. For a given route γ i , each cylinder Ω j is associated with two decision variables s ij and t ij . The variable s ij defines whether Ω j is active or not with respect to the route γ i :

s ij =        0, if Ω j not active 1, if Ω j active ∀i = 1, • • • , N ; ∀j = 1, • • • , M, (3)
while t ij defines the way an active obstacle Ω j is avoided on route γ i :

t ij =                0, if turn counter-clockwise 1, if turn clockwise 2, if impose a level flight below obstacle Ω j ∀i = 1, • • • , N ; ∀j = 1, • • • , M. (4)
Note that when an obstacle is not active (s ij = 0), the value of t ij is not applicable in the design.

The proposed decision variables are based on the available obstacle avoidance strategies and take into account the considered route shape.

C. Constraints

TMA is one of the most complex types of airspace; many constraints have to be satisfied. These constraints fall into two main categories: operational constraints related to air traffic operations (such as obstacle avoidance and route separation), and environmental constraints (such as noise abatement). SIDs/STARs are designed to satisfy these constraints and to deal with the dense traffic converging to and diverging from airports. The following constraints are considered.

Obstacle avoidance

This constraint is considered and handled through the choice of the decision variables.

Route separation

The standard separation norm between pairwise aircraft in a TMA is 3NM in the horizontal plane or 1000ft in the vertical plane. Routes are designed to respect these separation criteria, as a consequence, aircraft following two different routes are automatically separated.

Runway alignment

A SID must join a route section smoothly after a take-off leg and a STAR must go directly to the corresponding runway before the FAF. To deal with the runway alignment constraint, the notion of buffer obstacles is proposed, as shown in Fig. 3. Each route γ i is associated with one buffer obstacle Ω bi in a cylinder shape, defined by (C bi (x bi , y bi ), r bi , z bi inf , z bi sup , t Ω b i ), where t Ω bi is the turn orientation on Ω bi and the other parameters are as defined in Sec. II A. The values of these parameters are input data and ensure that the tangent segment connecting the buffer obstacle and the start point is parallel to the corresponding runway, so that the designed route joins the runway directly. A buffer obstacle is always active and the choice of its avoidance strategy is either counter-clockwise turn or clockwise turn.

RF related constraints

The first constraint related to RF is the minimum turn radius R min . The turn radius depends on the bank angle Φ and ground speed V GS [START_REF]The United States Standard for Area Navigation (RNAV) (Order 8260[END_REF][START_REF]United States Standard for Performance Based Navigation (PBN) Instrument Procedure Design[END_REF], and is calculated using the following equation:

R = V 2 GS g • tan(Φ) , ( 5
)
where g is the acceleration of gravity. According to Eq. ( 5), for a fixed value of V GS , the lowest radius of a flyable RF leg is given by

R min (V GS ) = V 2 GS g•tan(Φmax)
, where Φ max is the maximum value of the bank angle. Consequently, in order to ensure that a RF leg is flyable for all aircraft following it, the radius of the RF leg has to be greater than

R min = V 2 GS,max g•tan(Φmax)
, where V GS,max is the maximum ground speed in a TMA. Taking V GS,max = 400kt and Φ max = 25 • (according to [START_REF] Eurocontrol | User Manual for the Base of Aircraft Data (BADA) Revision 3.9[END_REF]), the corresponding minimum radius R min is equal to 5NM. To satisfy the constraint on minimum turn radius, a pre-processing is applied to the obstacles: for an obstacle with radius lower than R min , its radius is increased to R min . Note that, as aircraft ground speed is very low near the runway, the radius of a runway buffer obstacle may be lower than R min . In order to avoid a very big turn, an upper bound on the RF leg radius is also imposed. Let R max be the largest radius of a RF leg. According to [START_REF] Garmin | Garmin Radius to Fix Leg Project Report[END_REF], the value of R max is set to be 13NM. In our pre-processing, all input obstacles are modeled as cylinders whose radius of lower and upper bases is lower than R max .

Another constraint related to RF is that RF legs must terminate at least 2NM prior to the FAF [START_REF]The United States Standard for Area Navigation (RNAV) (Order 8260[END_REF].

Monotonicity of the vertical profile

The vertical profile for a SID should be monotonically ascending and for a STAR should be monotonously descending. That is, for t 1 , t 2 ∈ [0, 1], if t 1 ≤ t 2 , Eq. ( 6) must be satisfied:

h i inf (t 1 ) ≤ h i inf (t 2 ), i = 1, • • • , N h isup (t 1 ) ≤ h isup (t 2 ), i = 1, • • • , N. (6)

Level flight related constraints

Since the changes in flight levels have an influence on the fuel consumption and CO2 emissions, the number of level flights on each route is bounded above by N max , usually set to 2. For each route

γ i : m j=1 max(t ij -1, 0) ≤ N max , i = 1, • • • , N. (7)
Moreover, as the altitudes of imposed level flights have a direct impact on the noise pollution, a minimum altitude H min for each level flight is defined. In practice, the following constraints are imposed: for an obstacle Ω j , if z j inf < H min , then no level flight is allowed below it, therefore ∀i, t ij ∈ {0, 1}. In particular, for the obstacles connecting to the ground (e.g., mountains, buildings), no level flight is allowed below them. They are avoided only by horizontal turns. The obstacles that can be avoided by level flight can be military areas and conflict zones between routes (to be presented in Sec. IV). Finally, to ensure the passenger comfort, the length of each level flight should not be too short, therefore a minimum length L min for each level flight is required.

D. Objective function

To define our objective function, let us first define, for each route γ i , a weighted sum L γi of two terms: the length of the horizontal curve γ i H and the total length of level flights. More precisely:

L γi = c 1 1 0 γ i H (t) 2 dt + c 2 i LF , (8)
where i LF is the length of route sections corresponding to level flights projected onto the horizontal plane. The coefficients c 1 and c 2 are penalty parameters, whose values are user-defined and depend on the weight associated with the corresponding term. The objective is to minimize the sum of

L γi , i = 1, • • • , N : L = N i=1 L γi (9)
(P) [START_REF] Pocchiola | Minimal tangent visibility graphs[END_REF] This is a combinatorial optimization problem as the decision variables take integer values. The numerous constraints to be considered make the SIDs/STARs design a very complex problem.

                                               min L = N i=1 L γi s.t.
Therefore, in Sec. III the simpler subproblem of designing a single optimal route is considered. The proposed solution approach for the design of multiple routes is presented in Sec. IV.

III. Designing one route

The strategy proposed in [START_REF] Zhou | Optimizing the design of a route in Terminal Maneuvering Area using Branch and Bound[END_REF] is used to design an individual route γ i avoiding obstacles in TMA in an optimal way. Note that all of the constraints presented in Sec. II C, except route separation, are considered. This is the basis of our proposed solution approach to build multiple routes, where the constraint of route separation is considered additionally. To design a single optimal route, a B&B method is applied. This is inspired by the approach presented in [START_REF] Eele | Path planning with avoidance using nonlinear branch and bound optimization[END_REF], where a path planning problem avoiding circular obstacles is studied in 2D. A branching strategy is proposed, where, for each obstacle, two branches are created depending on whether an obstacle is avoided by turning clockwise or counter-clockwise. This branching strategy is extended to also allow obstacles to be avoided by imposing a level flight below the obstacle. Our branching strategy is illustrated in Fig. 4.

More precisely, in order to branch on an obstacle Ω ij , it is first set as active (s ij = 1) or not (s ij = 0); when it is active, three branches are developed in order to account for the 3 possibilities of avoiding it: counter-clockwise (t ij = 0), clockwise (t ij = 1) or imposing a level flight (t ij = 2), as illustrated in Fig. 5. More precisely, the counter-clockwise and clockwise turn avoid the obstacle along its border (Figs 5a), 5b)), and the level flight is imposed under the lower basis of the obstacle (Fig 5c )).

Let SP be a subproblem in B&B for the generation of γ i , where only some of the values of

a) (s ij , t ij ) = (1, 0), 2D View b) (s ij , t ij ) = (1, 1), 2D View c) (s ij , t ij ) = (1, 2), 3D View
LB SP = c 1 1 0 γ SP H (t) 2 dt + c 2 L min M j=1 max(t ij -1, 0) (11)
where γSP H is the horizontal profile of γSP , and t ij , j = 1, • • • , M are related to the choices of the obstacle avoidance strategies in SP . The obtained lower bound is used to identify whether a branch requires further exploration. The pseudo-code of the B&B method for generating a single optimal

route γ i , i ∈ {1, • • • , N } is presented in Algorithm 1.
Algorithm 1 Generation of a single route using B&B Require: start and end points (Ai, Bi), the altitude of the start point HA i , buffer obstacle Ω bi , obstacles

Ωj, j = 1, • • • , M 1: Initialize: subproblems list SP L = ∅, current best length LCB = +∞ 2:
Build the direct route connecting Ai and Bi, denoted as γ 0 i 3: if γ 0 i satisfies all constraints then 4:

γ 0 i is the optimal route, γi := γ 0 i , LCB := L γ 0 i 5: else 6:
Develop new branches according to Fig. 4, add them to SP L Build the corresponding vertical cones γSP V

12:

Compute the lower bound LBSP according to Eq. ( 11)

13:

if LBSP < LCB then 14:

if γSP satisfies all constraints and γSP is touched by all active obstacles then 15:

Compute Lγ SP according to Eq. ( 8)

16:

if Lγ SP < LCB then else if SP has remaining obstacle to be branched on then else if SP has remaining obstacle to be branched on then The route is built in the following way. First, the active obstacles (excluding the buffer obstacle)

avoided by a turn (t ij = 0 or 1), are numbered in increasing order of length(A i , P roj (AiBi) C j ),

where P roj (AiBi) C j is the projection of C j (the center of the two bases of obstacle Ω j ) onto the line

(A i B i ).
Next, the horizontal curve γ i H is computed by tangentially connecting, the buffer obstacle and successive active (s ij = 1) obstacles which are associated with counter-clockwise (t ij = 0)

or clockwise (t ij = 1) turns, following the order determined above. The alternative scenarios for connecting tangent points are illustrated in Fig. 6.

After obtaining γ i H , the associated vertical profile γ i V is in the form of a cone bounded by two straight line segments whose slopes are α min,T O , α max,T O in a SID case (respectively, α min,LD , α max,LD in a STAR case) (Fig. 1). For an active obstacle Ω j with t ij = 2, before imposing a level flight, we check whether the cone intersects Ω j both in the horizontal and vertical plane. We first check the intersection in the horizontal plane, as shown in Fig. 7. The length of route sections from A i to the two intersection points are denoted as d 1 and d 2 . When an intersection exists in the horizontal plane, we continue to check in the vertical plane whether the cone intersects Ω j . If neither of the two conditions:

h i inf (d 1 )
≥ z jsup (Fig. 8a))

h isup (d 2 ) ≤ z j inf (Fig. 8b)) (12)
is satisfied, then the vertical intersection exists. If intersection occurs in both the horizontal and the vertical planes, then the corresponding level flight is imposed under the obstacle. Moreover, the vertical profiles of the lower and the upper bounds are computed separately. For each bound, the level flight starts at the point where the altitude reaches z j inf , and ends at the point where the corresponding length in the horizontal plane is d 2 . Afterwards, the vertical profiles of the upper and the lower bounds continue to increase with the slopes α max and α min respectively. Note that if an active obstacle with t ij = 2 is not intersected by the cone associated with the horizontal route, then the corresponding level flight cannot be imposed. In such a case, the route is infeasible regarding our definition of "active obstacle".

IV. Designing multiple routes

To design N conflict-free routes avoiding obstacles, where a conflict between two routes is defined as a violation of the minimum separation (3NM in the horizontal plane, or 1000ft in the vertical plane), a sequential 1-against-n approach is proposed, based on the B&B method presented in Sec. III. Routes are generated sequentially in decreasing order of their traffic load. Previously built routes become obstacles for new routes. Note that, ordering the routes by traffic load provides better chance of finding near-optimal solutions, but does not guarantee optimality.

Each new route is initially generated using the B&B method proposed in Sec. III. Pairwise conflicts are then detected between the new and the previously generated routes. The method used for conflict detection is presented in Sec. IV A. The detected conflicts are clustered into groups.

Each group of conflicts is associated with a fictitious obstacle, modeled as a cylinder as in Sec.

II A. The method used to cluster conflicts and create the fictitious obstacles is described in Sec.

IV B. In the case when more than one group of conflicts exist, the corresponding fictitious obstacles are ordered according to their projection lengths onto the line connecting the start and end points of the route, as explained in Sec. III. In order to solve the conflict corresponding to the first

fictitious obstacle Ω f (C f (x f , y f ), r f , z f inf , z fsup )
, the natural solution is to perturb the considered route locally around this conflict zone. Again, three possible perturbations are proposed: turning counter-clockwise, turning clockwise, or imposing a level flight. Afterwards, the remaining route section is built by re-applying B&B, and then applying the conflict detection technique and route perturbation method to the newly generated section of the route. Thus, the new route is built progressively, and each of its sections is perturbed based on the optimal route structure given by the B&B method. Note that a conflict is not guaranteed to be solved by one of the available route perturbation strategies. It may happen that there are not enough degrees of freedom available to be able to solve the conflict. In this case, the proposed algorithm provides information on the unsolved conflict and on the length of involved route sections. A pseudo-code of the proposed B&Bbased sequential approach is presented in Algorithm 2. The following subsections describe the key intermediate steps.

A. Conflict detection

In steps 6 and 15 of Algorithm 2, a conflict detection technique is applied on pairs of routes.

Detecting conflicts between 3D-routes, especially when the vertical profiles are cones instead of curves, is not an easy problem. A two-steps scheme is proposed to deal with this problem. First, a check is performed in the horizontal plane, to see if any pairwise routes violate the 3NM separation.

Then, for the route sections found to be in violation of horizontal separation constraint, a check is performed in the vertical plane, to see if they violate the 1000ft separation. If a route section violates both horizontal and vertical separation constraints, then it is in conflict.

Algorithm 2 Generation of multiple routes using the B&B-based approach Require: routes start and end points, altitudes of the corresponding start points, and the associated buffer obstacles {(Ai, Bi), HA i , Ω bi |i = 1, • • • , N } (classified in an increasing order of traffic load), obstacles

{Ωj|j = 1, • • • , M } 1:
Initialize: list of optimal routes list = ∅ 2: Apply B&B to generate the optimal route γ1 connecting (A1, B1) call Algorithm 1((A1,B1), HA i , Apply B&B to generate an initial route γ 0 i connecting (Ai, Bi) call Algorithm 1(A1,B1),

Ω b1 , {Ωj|j = 1, • • • , M })
HA i Ω bi , {Ωj|j = 1, • • • , M }) 6:
Detect pairwise conflicts between γ 0 i and routes in list

7:
if no conflict detected then 8:

Set γi := γ 0 i , add γi to list while conflict detected do 12:

Cluster conflicts and associate a fictitious obstacle with each cluster 13:

Order the fictitious obstacles with respect to (Ai, Bi) in the same way as in Sec. III

14:

Perturb γ cur i locally around the conflict zone corresponding to the first fictitious obstacle

Ω f call Algorithm 3(γ cur i , Ω f , {Ωj|j = 1, • • • , M }, list) 15:
Detect pairwise conflicts between the remaining route section after the last perturbation on γ cur i and routes in list 2D-grid is used to cover all the horizontal curves. The idea of using a 2D-grid for conflict detection has been applied in [START_REF] Chaimatanan | A Hybrid Metaheuristic Optimization Algorithm for Strategic Planning of 4D Aircraft Trajectories at the Continental Scale[END_REF]. The dimension of a cell on the 2D-grid is 3NM×3NM, as defined by the horizontal separation norm. Each cell is identified by two indices (I x , I y ), where I x is the index on the x-axis and I y is the index on the y-axis, I x , I y ∈ N.

Each horizontal curve is discretized by using a discretization step δt. In order to guarantee that in each cell occupied by a curve there is at least one discretization point taken on that curve in that cell, a post-processing to add discretization points is applied. Let

0 = t i 0 < t i 1 < t i 2 < • • • < t i Ni = 1 be a subdivision of [0, 1]. Each t i k corresponds to a discretization point P i,k on curve γ i H , where P i,k = γ i H (t i k ).
A curve section is defined as the section of a curve between two successive discretization points. Let P i,k and P i,k+1 be two successive discretization points on curve γ i H , then the curve section between P i,k and P i,k+1 is denoted by γ i k,k+1 . An illustration of curve discretization on a 2D-grid is presented in Fig. 9a).

A horizontal conflict occurs when the minimum distance between two curve sections (in the form of arc or segment) on two different curves is less than 3NM. This minimum distance is computed analytically. For γ i H , the horizontal conflict is detected by successively considering the cells occupied by its curve sections. For each occupied cell, instead of checking the whole grid, only the neighboring non-empty cells are checked, since the horizontal conflict with other curves only occurs in the same or neighboring cells. In such a way, the efficiency of conflict detection is improved. More precisely, suppose that P i,k (respectively, P i,k+1 ) is located in cell (I x i,k , I y i,k ) (respectively, (I x i,k+1 , I y i,k+1 )). Then the neighboring cells where potential horizontal conflict may occur are

{(I x , I y )| min(I x i,k , I x i,k+1 ) -1 ≤ I x ≤ max(I x i,k , I x i,k+1 ) + 1, min(I y i,k , I y i,k+1 ) -1 ≤ I y ≤ max(I y i,k , I y i,k+1 ) + 1}.
For example, consider the curve γ 1 H in Fig. 9a). The discretization points on this curve are

P 1,k , k = 1, • • • , 9.

Vertical detection

If a horizontal conflict is detected, a further check in the vertical plane is needed. Denote H i,k inf (H i,k sup ) the lower (upper) bound of the cross section in the vertical plane at discretization point P i,k on route γ i , where

H i,k inf = h i inf (t i k ) and H i,k sup = h isup (t i k ). Suppose that two curve sections, γ i k,k+1
and γ i l,l+1 , on curves γ i H and γ i H respectively (i = i ) have a horizontal conflict. These two curve sections are in vertical conflict when neither of the following conditions is satisfied:

min H i,k inf , H i,k+1 inf -max H i ,l sup , H i ,l+1 sup > 1000ft min H i ,l inf , H i ,l+1 inf -max H i,k sup , H i,k+1 sup > 1000ft (13)
If one of these conditions is satisfied, which implies that two curve sections are separated vertically at their extremities, then they are separated along the sections. This is because the vertical profiles are monotonically increasing. Even though this detection method brings an additional margin in the separation, it has the advantage of being simple to implement. Finally, if two curve sections violate both horizontal and vertical separation constraints, they are in conflict.

B. Clustering conflicts and creating fictitious obstacles

After detecting pairwise conflicts between the current route and previously generated routes, each conflict zone is associated with a cylinder-shaped fictitious obstacle (step 12 in Algorithm 2), modeled using the method presented in Sec. II A. In order to generate fictitious obstacles, a conflict cell is defined as a cell containing at least one conflicting curve section on the current route. The basis of the cylinder is the minimum (maximum) altitude associated with the conflict cells in the corresponding group, as defined in Eq. ( 14) and Eq. ( 15). An example is illustrated in Fig. 10. if γ lf i solves the corresponding conflict then

z inf (I x , I y ) = min{H i,k inf |k ∈ K i , i ∈ I} (14) z sup (I x , I y ) = max{H i,k sup |k ∈ K i , i ∈ I} (
Ω f (C f (x f , y f ), r f , z f inf , z fsup ) (step 14 in
Algorithm 2). A fictitious obstacle is associated with two decision variables, s ij and t ij . When a fictitious obstacle is active (s ij = 1), it can be avoided by three strategies: turn counter-clockwise (t ij = 0), turn clockwise (t ij = 1) and impose a level flight (t ij = 2). Since the Continuous Climb Operations (CCO) and Continuous Descent Operations (CDO) are preferred in TMAs, in our route perturbation strategies horizontal perturbation has a priority over a vertical perturbation. If the horizontal perturbation is able to solve the considered conflict, then there is no need to try the vertical perturbation. A pseudo-code of the route perturbation is presented in Algorithm 3. In the following, the different route perturbation strategies are presented in more detail.

Perturbing a route by counter-clockwise and clockwise turns

In the case of avoiding a fictitious obstacle by a turn (t ij = 0 or 1), the perturbed route must be both smooth and flyable, as well as close to the current route γ cur i . Indeed, in the route perturbation process, the aim is not only to reduce conflicts, but also to keep the perturbed route as close as possible to the optimal route generated by B&B. To meet these objectives, for each fictitious obstacle, four smoothing obstacles tangent to both the considered fictitious obstacle and γ cur i are generated. The radius of the smoothing obstacles r s (r s ≥ R min ) is a user-defined parameter and it is the same for all smoothing obstacles. Figure 11 • First, set A i and B i as the start and end points of γ ccw,par i respectively.

• Second, create a list of obstacles, containing all the active obstacles when building γ cur i

• Third, add the fictitious obstacle Ω f to the list of obstacles and order the obstacles in the list according to their projection length to the line (A i B i ), as described in Sec. III.

• Then, for obstacles whose order index is lower than that of Ω f , keep them as active and keep their avoidance strategies; for obstacles whose order index is higher than that of Ω f , set them as not-considered.

• If an active real obstacle is avoided by a turn and it is intersected with the fictitious obstacle, also set it as not-considered. • Add smoothing obstacle 1 and 3 before and after the fictitious obstacle respectively.

• Finally, build the route in the same way as presented in Sec. III.

After the partial route has been constructed, pairwise conflicts are detected between γ ccw,par i and previously generated routes. If any residual conflict is detected, the radius of the fictitious obstacle is iteratively increased by a user-defined value δ R , until the conflict is solved or the radius reaches the maximum radius R max . In the case when the radius of the fictitious obstacle reaches R max and a residual conflict still exists, γ ccw,par i is set to be the perturbed route portion satisfying all constraints except route separation and with the smallest length involved in conflicts obtained in the various iterations.

To complete the perturbed route γ ccw i (i.e. build the route section between B i and B i ), a list of obstacles is created, including all the real obstacles as well as the fictitious and smoothing obstacles created when perturbing the route. Then, in order to keep the perturbed route section between A i and B i unchanged in B&B, the indices of the obstacles as well as the values of the decision variables considered when building γ ccw,par i are not changed. Afterwards, the other obstacles not-considered when building γ ccw,par i are ordered according to their projection length on the line (A i B i ). To build the remaining part of the route, the values of the decision variables for the non-active obstacles are determined in B&B. In the case when the clockwise avoidance strategy is chosen, smoothing obstacles 2 and 4 are taken to build γ cw,par i and their turn direction is counter-clockwise, as presented in Fig. 12b). The approach for building γ cw i is similar to the one for building γ ccw i .

Perturbing a route by imposing a level flight

In the case of avoiding fictitious obstacles by imposing a level flight (t ij = 2), there are two different scenarios, depending on the relative position between the fictitious obstacle and γ cur i . Each fictitious obstacle has two intersections with γ cur i in the horizontal plane, since the corresponding conflict zone is located around γ cur i . Denote the intersection closer to B i as B i .

The first case occurs when B i is located on an arc of γ cur i , as shown in Fig. 13a). Then the partially perturbed route γ lf,par i connecting A i and B i is built in the following way:

• First, set A i and B i as the start and end points of γ lf,par i respectively.

• Second, create a list of obstacles, consisting of all the active obstacles when building γ cur i

• Then, add the fictitious obstacle to the list of obstacles and order them according to their projection length to line (A i B i ).

• For the obstacles in the list whose order index is lower than that of the obstacle where the second intersection point is lying on (including the obstacle itself), keep their avoidance strategies;

for the other obstacles, set them as not-considered.

• Finally, build the route in the same way as presented in Sec. III.

The second case occurs when B i is located on a segment of γ cur i . To maintain the horizontal form of the perturbed route, two smoothing obstacles are added, which are tangent to both the fictitious obstacle and the line that the segment is lying on, as shown in Fig. 13b). Then the partially perturbed route γ lf,par i connecting A i and B i is built in the following way:

• First, set A i and B i as the start and end points of γ lf,par i .

• Second, create a list of obstacles, consisting of all the active obstacles when building γ cur i

• Then, add the fictitious obstacle to the list of obstacles and order them according to their projection length to line (A i B i ).

• For the obstacles whose order index is lower than that of the fictitious obstacle, keep their avoidance strategies; for the other obstacles, set them as not-considered.

• Finally, build the route in the same way as presented in Sec. III.

After building γ lf,par i , pairwise conflicts are detected between γ lf,par i and the previous routes. If any residual conflict is detected, the altitude of the lower basis of the fictitious obstacle is iteratively decreased by a user-defined value δ H , until the conflict is solved or the altitude of the lower basis reaches H min . In the case when the altitude of the lower basis of the fictitious obstacle reaches H min and a residual conflict still exists, γ lf,par i is set to be the perturbed route portion satisfying all constraints except route separation, and with the smallest length involved in conflicts obtained in the iterations.

To complete the route section between B i and B i , a list of obstacles is created, including all the real obstacles as well as the fictitious and smoothing obstacles created when perturbing the route.

These obstacles are then ordered according to their projection length on the line (A i B i ). In order to keep the form of γ lf,par i unchanged, the values of the decision variables for the obstacles whose order index is lower than that of the fictitious obstacle (including the fictitious obstacle itself) are kept. Afterwards, the methods for completing the remaining route section in the first case (the second intersection point is located on an arc) and in the second case (the second intersection point is located on a segment) are slightly different.

• In the first case, the values of the decision variables of the obstacle where the second intersection point is located is kept. The values of the decision variables for the other obstacles are determined in the modified-B&B method.

• In the second case, the values of the decision variables for the other obstacles (including the two smoothing obstacles) are determined in the modified-B&B method. Meanwhile, only counter-clockwise turn is allowed on smoothing obstacle 1 when it is active and only clockwise turn is allowed on smoothing obstacle 2 when it is active.

D. Post-processing technique (step 17 in Algorithm 2)

In the route perturbation process, each horizontal perturbation is realized by sequentially avoiding a smoothing obstacle, a fictitious obstacle and a second smoothing obstacle, as shown in Fig. 12.

Thus the perturbed route section is composed by three successive arcs lying on the three obstacles.

This kind of curve is flyable in reality, since a RF leg can join to another RF leg with an opposite turn direction and a different radius [START_REF]Performance Based Navigation (PBN) Operational Approval Handbook[END_REF]. However, in order to simplify the operation for pilots, a post-processing technique on γ cur i is proposed. For a real SID/STAR, a horizontal perturbation occurs to avoid an obstacle or to avoid another route, afterwards the route usually connects directly to the next route section. Thus, in a SID (respectively, STAR) case, for each horizontal route perturbation, the second (respectively, first) smoothing obstacle is set as non-active and the values of the decision variables for the other obstacles are preserved and the route γ cur i is re-built. If no more conflict with the previously generated routes is induced by the post-processing, then the new γ cur i is accepted.

segment. Next, γ 2 is generated. It is perturbed clockwise to pass above γ 1 (Fig. 14b)). Afterwards, γ 3 is computed. It is perturbed counter-clockwise to avoid the conflict with γ 1 (it consequently passes below γ 1 (Fig. 14c))). Next, route γ 4 is generated. As it is in conflict with γ 3 , and no perturbation in the horizontal plane is able to solve this conflict, a level flight is imposed to let γ 4 pass below γ 3 . An illustration of the level flight on γ 4 in the vertical plane is shown in Fig. 14e).

Another conflict occurs when γ 4 meets γ 1 , thus a clockwise perturbation is executed on the route section that follows the level flight on γ 4 (Fig. 14d)). Finally, γ 5 is generated and two successive perturbations are performed to avoid the conflict with γ 2 and with γ 3 respectively. The proposed B&B-based approach solves all conflicts with a total simulation time of 2s.

Case 2: 9 obstacles (N = 5, M = 9)

In this case, 9 obstacles are added to the previous configuration (case 1), to test our approach in a more complex situation. With the presence of obstacles, the space for route perturbation is reduced, thus the complexity of the problem is increased. The simulation result is presented in Fig. 15. The layout of the obstacles, as well as the optimal routes obtained are presented in Fig. 15a). Figure 15b) shows that both γ 1 and γ 2 circumnavigate the obstacles before they meet each other due to obstacle avoidance. When their paths cross, their altitudes are not the same as in case 1. Thus, route γ 2 is perturbed clockwise with a bigger turn radius compared to case 1. Due to the presence of obstacles, the γ 3 generated is also different to that generated in case 1. As the shapes of γ 1 and γ 2 are also changed, no conflict occurs between γ 3 and γ 1 , nor between γ 3 and γ 2 .

Next, route γ 4 requires a longer level flight compared to case 1, so that it passes below both γ 3 and γ 1 . The profile of γ 4 in the vertical plane is illustrated in Fig. 15e). Note that, in case 1, γ 4 avoids γ 1 by a clockwise turn, however this is not possible here, since the clockwise perturbation intersects an obstacle. Finally, γ 5 is generated and consists of two successive horizontal perturbations, as shown in Fig. 15f). Its shape is different from the one in case 1 (Fig. 14f)), since the post-processing in case 2 does not provide a better solution. The proposed B&B-based approach solves all conflicts with a total computing time of 1.1s. The simulation result is then compared with the standard published SIDs and STARs.

Input data

CDG airport has four parallel runways, as shown in Fig. 16a), where the numbering of the 8 thresholds (two thresholds per runway) is also illustrated. In general, two runways are used for take-off and the other two are used for landing. The radar data used for testing correspond to a day when the runway thresholds used for take-off (respectively, landing) are 09R and 08L (respectively, 27R and 26L). The same thresholds are used in our simulation.

The radar data used are illustrated in Fig. 17 for jet aircraft are selected from the published Jeppesen charts [START_REF] Jeppesen | LFPG[END_REF]. These standard routes, as well as the related waypoints, are also presented in Fig. 17, where the gray segments represent SIDs and black ones represent STARs.

In order to compute the traffic load on each principal route, an exit TMA window is defined for each exit point on a SID. Similarly, an entry TMA window is defined for each entry point on a STAR. The traffic load on each principal route is calculated by summing the number of flights passing through the corresponding exit or entry window. In Fig. 17, the selected routes are numbered according to decreasing traffic load. The information on the selected routes, including the threshold used and traffic load (in percentage) are presented in Table 4.

In our simulation, the design of 15 routes (denoted as γ i , i = 1, • • • , 15) is considered. Their start and end points, as well as traffic loads correspond to that of the 15 selected standard principal routes. More precisely, in a SID case, the start and end points of γ i are the same as that of the standard route i, while in a STAR case, the end point of γ i corresponds to the TMA entry point on the standard route i, and the start point is the FAF of the corresponding runway threshold. As described in Sec. II, each route to be designed is associated with a buffer obstacle. The characteristics of the buffer obstacles are presented in Table 5. It is assumed that the routes using the same runway use the same buffer obstacle. The TMA of CDG airport is located in a relatively simple geographical environment, the only area to be avoided is Paris city, for the sake of noise abatement. This area is modeled as an obstacle. Thus M = 1 in our simulation, and the characteristics of this obstacle are (x 1 , y 1 , r 1 , z 1 inf , z 1sup ) = (92.43NM, 113.16NM, 5NM, 0, 50000ft). The buffer obstacles as well as the obstacle representing Paris city are illustrated in Fig. 16b).

Simulation results

Routes are built sequentially according to decreasing order of traffic load. The simulation result is presented in Fig. 18. The presented results show that two main conflict zones occur in the simulation. One conflict zone is between departure routes γ 9 , γ 12 and arrival routes γ 4 , γ 15 , the other is between departure routes γ 3 , γ 5 and arrival route γ 6 . This is coherent with our analysis of the real radar data. For departure routes γ 9 , γ 12 and arrival routes γ 4 , γ 15 in the simulation, route γ 4 is the first of the four to be generated, since it has a higher traffic load. Routes γ 9 and γ 12 are generated next. In order to avoid the conflict with γ 4 , route γ 9 is perturbed clockwise and route γ 12 is perturbed counter-clockwise. Route γ 15 is the final route to be generated, and is initially in conflict with both γ 12 and γ 9 . To avoid γ 12 , a level flight is imposed, and to avoid γ 9 , a counter-clockwise turn is applied. These four routes are illustrated in 3D in Fig. 18c), and the vertical profile of γ 15 is illustrated in Fig. 18d). For departure routes γ 3 , γ 5 and arrival route γ 6 , route γ 6 is generated after γ 3 , γ 5 , and a level flight is imposed on γ 6 , so as to avoid conflict with previously generated routes. These three routes are illustrated in 3D in Fig. 18e), and the vertical profile of γ 6 is illustrated in Fig. 18f). The computing time for obtaining the 15 optimal routes is about 9.8s. the dashed sections represent level flights. The total length of the 15 selected standard routes is 1358.26NM, and it is reduced by around 45NM using the proposed B&B-based approach. Except routes γ 9 , γ 12 and γ 15 , which are perturbed to avoid conflicts, the other routes obtained by the B&B-based approach are all shorter than the standard routes. This is especially true for the routes with relatively heavy traffic load, since they have a priority order in the design and therefore the generated routes are straight line segments. The length of level flight on γ 6 is 25.7NM, longer than the corresponding level flight on the standard route, which is 20NM (according to [START_REF] Jeppesen | LFPG[END_REF]). The length of level flight on γ 15 is 30Nm, again longer than the corresponding level flight on the standard route, which is about 25NM (according to the radar data). On average, CDG airport accommodates a total of 1300 flights per day (including both departures and arrivals). By using the optimal routes obtained in our simulation, the total flown distance is reduced by 5889NM (10906km) per day, and the total length of level flights is increased by 640NM (1185km) per day.

From an economic point of view, this can lead to considerable savings. Consider an A320 aircraft as an example. The average fuel consumption of one engine in the take-off and landing phases is about 6kg/km. Thus, the total fuel consumption is reduced by about 131ton per day, leading to about $65,000 in the reduction of jet fuel costs (assuming a jet fuel price of 61.3$/bbl (497.5$/ton) [START_REF]Fuel Price Analysis[END_REF]). Since the length of level flights is only slightly increased, and the affected traffic load is only 8.8%, the impact of the increase in the total length of level flights on fuel consumption is negligible.

In addition, because aircraft flying on different routes are automatically separated, the workload for the controllers is reduced. This may lead to a reduction in the expenses related to ground services.

Readers interested by the detailed input data (including the coordinates of the start and end points, characteristics of the buffer obstacles, and traffic load of each route), or detailed numerical results (including the initial length of route sections involved in conflicts, the optimal route length, the length of level flight on each optimal route) of tests 1 and 2, may refer to [START_REF] Zhou | [END_REF] for more information. More simulation results on the design of a single or multiple routes are presented in [START_REF] Zhou | Optimisation des procédures de départ et d'arrivée dans une zone terminale[END_REF].

VI. Conclusions

In this paper, a decision support methodology is proposed for the design of multiple Standard Instrument Departure routes (SIDs) and Standard Terminal Arrival Routes (STARs) in a Terminal

Maneuvering Area (TMA) at a strategic level. The objective of the proposed approach is to design SIDs/STARs that minimize total route length, while taking into account operational and environmental constraints; in particular, obstacle avoidance and route separation constraints. A route is modeled as a curve formed by segments and arcs in the horizontal plane, associated with a cone enveloping all ascent (or descent) aircraft profiles in the vertical plane. Three different ways to avoid an obstacle are proposed: turning clockwise, turning counter-clockwise, and imposing a level flight under the obstacle. These maneuvers are compatible with operational requirements and correspond to what is currently implemented.

The problem is modeled as an Integer Programming optimization problem. To design a single optimal route, a Branch and Bound (B&B) method is developed, whose branching strategies correspond to the three ways obstacles can be avoided. To design multiple routes, the proposed B&B method is used to sequentially generate routes, and a route-modification step is introduced to deal with conflicts between routes. The proposed methodology is first tested on an artificially generated problem, and then applied to the design of 15 routes in the TMA of Paris Charles-de-Gaulle (CDG) airport. The routes obtained are continuous and smooth. These routes result in a reduction in total route length compared with the published standard charts in the case of Paris CDG airport, which in turn results in a sizable reduction in jet fuel consumption.

  the total number of routes to be built, N ∈ N the start point Ai and starting altitude HA i for route i, i = 1, • • • , N the end point Bi for route i, i = 1, • • • , N traffic load on each route minimum (maximum) take-off slope αmin,T O (αmax,T O ) minimum (maximum) landing slope αmin,LD (αmax,LD) minimum (maximum) radius on a RF leg Rmin (Rmax) maximum number of level flights on each route Nmax minimum altitude of each level flight Hmin minimum length of each level flight Lmin A. Obstacle and route modeling

Fig. 1

 1 Fig. 1 Examples of γ iH and γ iV .

Fig. 2

 2 Fig. 2 Take-off and landing profiles in Paris CDG airport.

Fig. 3

 3 Fig. 3 Runway buffer obstacle.

Fig. 4

 4 Fig. 4 Branch and Bound Branching Strategy.

Fig. 5

 5 Fig. 5 The routes associated with different values of the decision variables.

17 :

 17 γi := γSP , LCB := Lγ SP 18:

19 :

 19 Develop new branches according to Fig.4, add them to SP L

22 : 26 :Fig. 6

 22266 Fig. 6 Different scenarios for connecting tangent points. In the notations [t i,k ] and [t i,k , t i,k+1 ], t i,k indicates the turn direction to avoid obstacle Ω k on route γ i .

Fig. 7 Fig. 8

 78 Fig. 7 Checking intersection between γ iH and Ω j in the horizontal plane.

3: Add γ1 to list 4 :

 4 for i=2 to N do 5:

Apply the post-processing technique on γ cur i 18 :, add γi to list 19

 1819 Set γi := γ cur i Route projections in the horizontal plane γ i H , i = 1, • • • , N are in the form of 2D-curves. A

The first curve section on γ 1 HFig. 9

 19 Fig. 9 An example of horizontal conflict detection.

Fig. 10

 10 Fig. 10 An example of clustering conflict cells and creating fictitious obstacles.

  minimum and maximum altitudes associated with a conflict cell (I x , I y ) are denoted as z inf (I x , I y ) and z sup (I x , I y ) respectively. Let I ⊂ {1, • • • , N } be the subset denoting the indices of the previously generated routes, and let {γ i k,k+1 |k ∈ K i , K i ⊂ {0, • • • , N i }, i ∈ I} be the set of curve sections in conflict with the current route in cell (I x , I y ). Then z inf (I x , I y ) and z sup (I x , I y ) can be computed:

15 )

 15 Hence, the geometries related to the conflicting route sections are transferred to the geometries of the corresponding conflict cells. Afterwards, the conflict cells are clustered into groups in such a way that two adjacent conflict cells are in the same group. Two cells (I k x , I k y ) and (I l x , I l y ) are adjacent when |I k x -I l x | ≤ 1 and |I k y -I l y | ≤ 1. Each group of conflict cells is associated with a fictitious obstacle, which is the smallest cylinder containing them. The altitude of the lower (upper)

C. Perturbing γ cur i Algorithm 3 1 :clockwise ⇒ obtain γ cw i and γ ccw i 2 : if both γ ccw i and γ cw i solve the corresponding conflict then 3 :γ cur i equals to the one with the shorter length 4 :and γ cw i solves the corresponding conflict then 5 :: else 7 :

 3123457 Route perturbationRequire: the current route γ cur i , the first fictitious obstacle in the ordered list Ω f , obstacles {Ωj|j = 1, • • • , M }, the list of previously generated routes list Re-generate the portion of route γ cur i involved in the conflict by perturbing it clockwise and counter-Set else if only one of γ ccw i Set γ cur i equals to the one that solves the corresponding conflict 6Re-generate the portion of route γ cur i involved in conflict by imposing a level flight ⇒ obtain γ lf i 8:

  i which has the minimum length involved in conflict on the perturbed portion of route

  illustrates the four smoothing obstacles given γ cur i and a fictitious obstacle. Let γ ccw i and γ cw i denote the new routes obtained by perturbing the portion of route γ cur i involved in a conflict counter-clockwise and clockwise respectively. When the counter-clockwise avoidance strategy is chosen, smoothing obstacles 1 and 3 are taken to build γ ccw i and the turn direction on both smoothing obstacles is clockwise, as shown in Fig. 12a). The point B i is the tangent point of γ cur i and smoothing obstacle 3. The partially perturbed route γ ccw,par i between the start point A i and point B i is built in the following way:

Fig. 11

 11 Fig. 11 Configuration of the smoothing obstacles

a) Case 1 :

 1 Fig. 13 Route perturbation in the vertical plane.

a)

  Fig. 14 Test 1: case 1, N = 5, M = 0, simulation results (runway buffer obstacles are

a)

  Fig. 15 Test 1: case 2, N = 5, M = 9, simulation results (runway buffer obstacles are

  , where the real traffic arriving at and departing from Paris CDG airport during one day are presented. The light gray areas represent departure a) Four parallel runways b) Buffer obstacles and a real obstacle

Fig. 16

 16 Fig. 16 Paris CDG airport configuration.

Fig. 17

 17 Fig. 17 One day radar data of Paris CDG airport, and principal SIDs and STARs.

Figure

  Figure. 18b) represents the comparison between the selected standard routes and our simulation

Fig. 18

 18 Fig. 18 Test 2: simulation results.

Table 2

 2 Input data related to route design

Table 4

 4 Test 2: traffic load and waypoints of the selected standard routes

	No.	1	2	3	4	5	6	7	8	9	10	11 12 13 14	15
	route type STAR SID SID STAR SID STAR STAR STAR SID STAR SID SID SID SID STAR
	threshold 27R 09R 08L 27R 08L 26L	26L	27R 09R 26L 09R 09R 09R 08L 27R
	traffic (%) 12.77 10.82 10.72 8.7 7.44 7.33 7.16 6.97 5.63 5.61 4.9 4.76 3.38 2.33 1.48

Table 5

 5 Test 2: runway thresholds and the corresponding buffer obstacles bi , y b i , r bi ) (NM) (z bi inf , z bi sup ) (ft)

			buffer obstacle
	threshold
		(x tΩ b i
	09R	(96.07, 124.59, 2)	(0, 50000)	1
	08L	(93.29, 117.61, 3)	(0, 50000)	0
	27R	(112.55, 129.61, 4.9)	(0, 50000)	0
	26L	(115.07, 120.4, 2.5)	(0, 50000)	1

Acknowledgments This work has been partially supported by Civil Aviation University of China (CAUC), by National Natural Science Foundation of China (NNSFC) through grant NNSFC 61201085, by the Fundamental Research Funds for the Central University of China through grant number ZXH2012H005, and by French National Research Agency (ANR) through grant ANR 12-JS02-009-01 "ATOMIC".

The authors would like to thank Serge Roux and Alana Moore for interesting discussions and feedback on the manuscript. The authors would also like to thank Gilles Baroin for his help in formatting the figures.

V. Simulation results

In this section some numerical results on the design of multiple routes are presented. The proposed methodology is first tested on an artificially generated problem. Then the design of 15 routes in the Paris CDG airport is presented. A comparison of the designed routes with the standard routes published in Jeppesen charts [START_REF] Jeppesen | LFPG[END_REF] is also presented. The coordinates are given in a Cartesian Coordinate System where the center is (47 • N, 0 • E). Tests were run on a Linux platform with a 2.4 GHz processor and 8 GB RAM. The input data, common to all the presented tests, are given in Table 3. The results of tests carried out on an artificially generated problem are presented first. The routes to be built are denoted γ i , i = 1 • • • , 5, where γ 1 and γ 2 are SIDs, and γ 3 , γ 4 and γ 5 are STARs. The routes are generated sequentially according to the increasing order of their indices i.

Table 3 Common input data

Case 1: 0 obstacles (N = 5, M = 0) First, consider the case when there are no obstacles. The simulation result is presented in Fig. 14.

Since there are no obstacles, route γ 1 , which is generated first, takes the form of a straight line


Jun Zhou
email: junzhou@recherche.enac.fr

Sonia Cafieri
email: sonia.cafieri@enac.fr

Daniel Delahaye
email: daniel.delahaye@enac.fr

Mohammed Sbihi
email: mohammed.sbihi@enac.fr

Optimizing the Design of a Route in Terminal Maneuvering Area Using Branch and Bound

Keywords: TMA, SID/STAR design, obstacle avoidance, modeling, Branch and Bound

The sharp increase in air traffic flow causes traffic congestion in airspaces near airports, called Terminal Maneuvering Areas (TMA). The departure and arrival traffic of airports follow pre-designed routes named Standard Instrument Departure (SID) routes and Standard Terminal Arrival Routes (STAR). Optimizing these routes is crucial to regulate air traffic. Currently, SIDs and STARs are designed manually, based on the airport layout and nearby constraints. The objective of this research is to propose a methodology for designing an arrival/departure route in TMA, taking into account some constraints including obstacle avoidance. The shape of a route in horizontal plan is a succession of arcs of circles and segments. The originality of our study is, on the one hand, that the horizontal route is associated with a cone in vertical plan enveloping all ascent (or descent) aircraft profiles, and on the other hand, a branching strategy in a Branch and Bound (B&B) framework tailored on the problem is proposed.

Introduction

The continuously increase in air traffic flow density causes traffic congestion in the areas surrounding airports, thus affecting the normal operation of air traffic. An area surrounding one or more neighboring airports is called Terminal Maneuvering Area (TMA), and it is designed to handle aircraft arriving to and departing from the airports. Optimizing departure and arrival procedures in TMA is therefore crucial to regulate air traffic flows. Most of the airports have pre-designed procedures indicating how aircraft depart from or arrive to airports. These procedures are called Standard Instrument Departure (SID) routes and Standard Terminal Arrival Route (STAR). A SID is a flight route followed by aircraft after take-off from an airport until the start of en-route phase. A STAR is a route which connects the last en-route way-point to the Initial Approach Fix. Currently, SID/STARs are designed manually according to operational requirements (ICAO Doc 8168), taking into account airport layout and nearby constraints. However, this kind of design is generally not very efficient and not expected to optimize any specific criterion. The objective of this work is to automatically design SID/STARs in 3D with respect to certain optimization criteria. Being this study at a strategic level, only static obstacles are taken into account. Rather than considering an individual flight, we deal with flows of flights, that is to say the flights following the same SID (or STAR) belong to the same flow.

The considered problem is in the framework of path planning. Specifically, it is a route design problem: as contrarily to trajectory design, we aim at designing routes that are not associated to the notion of time. The problem of path planning has been studied since 1980s in the robotic domain [START_REF] Gallina | A technique to analytically formulate and to solve the 2dimensional constrained trajectory planning problem for a mobile robot[END_REF] [START_REF] Souissi | Path planning: A 2013 survey[END_REF]. Nowadays planning optimal aircraft paths becomes a rich and dynamic research area, some approaches have been summarized in [START_REF] Delahaye | Mathematical models for aircraft trajectory design: A survey[END_REF]. In particular, computing the shortest path between two points, given a number of obstacles, is one of the most extensively studied topic [START_REF] Lozano-Pérez | An algorithm for planning collison-free paths among polyhedral obstacles[END_REF] [8] [10] [7] [START_REF] Kim | Shortest paths for disc obstacles[END_REF]. Most of these works search for the shortest path in a horizontal plan, while in our study we consider searching for the shortest path in 3D. Moreover, in several works the obstacle is modeled as polygon [START_REF] Lozano-Pérez | An algorithm for planning collison-free paths among polyhedral obstacles[END_REF] [8] [START_REF] Storer | Shortest paths in the plane with polygonal obstacles[END_REF], and in a few works as circle [START_REF] Pocchiola | Minimal tangent visibility graphs[END_REF] [START_REF] Kim | Shortest paths for disc obstacles[END_REF]. In the present study we model the obstacle as cylinder in 3D and the projection to the horizontal plan is in form of circle.

Specifically in aircraft paths designing domain, even though there is a large number of researches, routes design in TMA is a particular problem for which to our knowledge there is not a rich literature. In [START_REF] Pfeil | Optimization of airport terminal-area air traffic operations under uncertain weather conditions[END_REF] the author designs terminal routes getting around obstacles with a modified A* algorithm. In our study, we design a route not only getting around obstacles, but also allowing level flights. Indeed, imposing a level flight in vertical plan is also an effective way to avoid obstacles, as it enriches the possible maneuvers and corresponds to what is done in the reality in a TMA. The different ways to avoid obstacles allow us to define specific branching strategies in a Branch and Bound tailored on the problem, that is another contribution of this work.

This paper is organized as follows. Section 2 introduces the route and obstacles modeling. Section 3 presents the proposed approach to solve the problem. Section 4 gives some preliminary simulation results. Finally, Section 5 draws conclusions and proposes future directions.

Problem Modeling

TMA is one of the most complex types of airspace. Many constraints have to be satisfied, falling into two categories: operational constraints related to air traffic operations (such as obstacle avoidance and flyable routes), and environmental constraints (such as noise abatement). SID/STARs are designed to satisfy these constraints and to deal with the dense traffic converging to and diverging from airports. The constraints in TMA make the SID/STARs design a very complex problem. Therefore in this study we consider the simpler subproblem of designing a single route avoiding obstacles and satisfying some other operational constraints. The obstacles in TMA could be mountains, cities, military area, etc. In the following, we present the way we model routes and obstacles.

A 3D route γ is defined by two elements: a curve γ H in a horizontal plan which is composed by a succession of arcs of circles (to bypass obstacles) and segments (to connect tangentially two arcs); a cone γ V in a vertical plan that contains all ascent (or descent) profiles of the aircraft flying on this route. The cone is defined by two straight lines whose slopes are the minimum and maximum values of the take-off (or landing) rate of the aircraft on this route. The idea of taking a cone that contains all vertical profiles is inspired by the behavior illustrated in Fig. 1, which shows some real take-off data in runway 08L of Paris CDG airport. From the figure we can see clearly that the vertical profiles are contained in a cone defined by two straight lines. The vertical profiles for landing are similar. This behavior is mainly due to the different aircraft mass and performance and to the effect of the wind.

In a horizontal plan, we define a starting point A (x A , y A ) and an ending point B (x B , y B ). In a SID case, the starting point is at the runway threshold and the ending point is an exit point of a TMA. In a STAR case, the starting point is an entry point of a TMA and the ending point is the Final Approach Fix (FAF). The horizontal route γ H is a smooth mapping defined as:

γ H : [0, 1] → R 2 (1)
where γ H (0) = (x A , y A ) and γ H (1) = (x B , y B ). In a vertical plan, γ V is defined as:

γ V : [0, 1] → I R t → [z in f (d(t)), z sup (d(t))] (2)
where I R defines the set of intervals of R, and

d(t) = t 0 γ H (s) 2 ds is the flown distance until t in horizontal plan, [z in f (d), z sup (d)]
is the interval defined by the cross section of the cone at d. The vertical profile γ V is associated to the horizontal curve γ H through the flown distance d. Figure 2 illustrates an example of how γ H is associated with γ V in the case of a SID, where α min and α max are the minimum and maximum values of take-off rate of aircraft on this route.

In the case of a SID, the starting point A is associated with an altitude z A ; the ending point B is associated with an altitude interval [z B , z B ] which indicates the altitude interval to exit TMA. Therefore the boundary conditions are γ

V (0) = [z A , z A ] and γ V (1) ⊂ [z B , z B ].
Similarly, in the case of a STAR, the altitude interval of the starting point A is [z A , z A ] which indicates the interval to enter TMA; the altitude of the ending point B is z B . The boundaries conditions are:

γ V (0) = [z A , z A ] and γ V (1) = [z B , z B ]. Note that, in order to guarantee that any SID reaches B within [z B , z B ], we set z B = z A + α min • d(A, B)
, where α min is the minimum taking-off slope, and d(A, B) is the Euclidean distance between the starting and ending points A and B. Similarly, to guarantee that any STAR reaches z B at FAF, we set

z A = z B + α min • d(A, B) and z A = z B + α max • d(A, B)
, where α min (α max ) is the minimum (maximum) landing slope. A level flight can be imposed to ensure that z B (respectively, z B ) is not exceeded in the case of a SID (respectively, STAR).

The obstacles, together with their protection areas, in number of m ∈ N, are modeled as cylinders in 3D, whose bases are parallel to the horizontal plan as presented in Fig. 3. Each cylinder

Ω i , i = 1, • • • , m is defined by (C i (x i , y i ), r i , z i in f , z i sup )
, where C i (x i , y i ) and r i are the center and the radius of the two bases respectively; z i in f and z i sup are the altitude of the lower and upper bases. These obstacles are numbered in Fig. 3 Obstacle Modelization Fig. 4 Obstacles Numbering an increasing order of length(A, Pro j (AB) C i ) where Pro j (AB) C i is the projection of C i onto the line (AB). An illustration is presented in Fig. 4.

Let us define an active obstacle as an obstacle that is touched by a route and has to be avoided according to one of the following maneuvers: turn counter-clockwise, turn clockwise or impose a level flight. Each cylinder Ω i is associated with two decision variables s i and t i : s i defines whether Ω i is active or not:

s i = 0, if Ω i not active 1, if Ω i active (3)
while t i defines the ways an active obstacle Ω i is avoided:

t i =    0, if turn counter-clockwise 1, if turn clockwise 2, if impose level flight (4)
Once the values of decision variables are chosen, the horizontal route is computed by connecting tangentially the successive active (s i = 1) obstacles with t i = 0 or t i = 1 in the increasing order of their numbering. Active obstacles with t i = 0 are bypassed counter-clockwise and those with t i = 1 are bypassed clockwise. This horizontal route is hence built piecewise: it is composed by (∑ m i=1 s i -∑ m i=1 max(t i -1, 0) + 1) straight line segments and (∑ m i=1 s i -∑ m i=1 max(t i -1, 0)) arcs of circles. Routes are therefore constrained to lie on the border of obstacles. In order to ensure flyable routes, the radius of the arcs of circles is imposed to be at least equal to 3 Nm (FAA Orders 8260.54A and 8260.58). Note that these arcs can be followed using a type of Performance Based Navigation (PBN) named Required Navigation Performance (RNP). Then, a vertical profile is associated to the horizontal route, taking into account α min , α max , and imposing a level flight below (respectively above) the active obstacle Ω i in a SID (respectively STAR) case when t i = 2. If some active obstacles with t i = 2 is not intersected by the cone associated with the horizontal

(s i ,t i ) = (1, 0), 2D View. c (s i ,t i ) = (1, 1), 2D View. d (s i ,t i ) = (1, 2), 3D View
route, then the route is unfeasible regarding to our definition of "active obstacle". Note that the way of building a horizontal route simplifies the computation but it does not necessarily lead to the shortest horizontal route between A and B.

An illustration of different values of the decision variables for an example of a SID in a TMA with one obstacle is presented in Fig. 5. In Fig. 5a, the obstacle is not active, so s i = 0. The horizontal route is a straight line segment connecting A and B. It is associated with a cone in the vertical plan. This route in the considered example is not a feasible one, because it intersects the obstacle. Then when the obstacle is active (s i = 1), 3 possibilities are considered to avoid it (Figs. 5b,5c,5d): turn counter-clockwise, turn clockwise and impose a level flight under the obstacle at altitude z i in f , corresponding to t i = 0, 1, 2 respectively.

Two examples of SIDs are presented in Fig. 6 showing how a route is computed, given the values of the decision variables, in the case with more than one obstacle. In example 1 (Figs. 6a,6b), given (s 1 ,t 1 ) = (1, 0), (s 2 ,t 2 ) = (1, 1), the horizontal route is composed by five parts: three segments and two arcs of circles. The three segments are used to connect tangentially the starting point to Ω 1 , Ω 1 to Ω 2 and Ω 2 to the ending point respectively. The two arcs are used to bypass Ω 1 counter-clockwise and Ω 2 clockwise respectively. In example 2 (Figs. 6c,6d), given (s 1 ,t 1 ) = (1, 2), (s 2 ,t 2 ) = (1, 1), the horizontal route is constructed by only bypassing Ω 2 , thus it is composed by two segments and one arc of circle. In vertical plan, when the route reaches the altitude of the lower basis of Ω 1 , a level flight is im-

∑ i=1 max(t i -1, 0) ≤ N max (5)
Second, as the altitudes of imposed level flights have a direct impact on the noise pollution, a minimum altitude H min for each level flight is defined. In practice, we impose the following constraints: in a SID case (respectively a STAR case), for an obstacle Ω i , if z i in f < H min (respectively z i sup < H min ), then no level flight is imposed below (respectively above) it, therefore t i ∈ {0, 1}.

Third, as to take into account the passengers comfort, the length of each level flight should not be too short, a minimum length L min for each level flight is imposed.

We minimize a weighted sum L γ of the length of the route γ in the horizontal plan and the length related to the level flights. More precisely:

L γ = c 1 × 1 0 γ H (t) 2 dt + c 2 × L min × m ∑ i=1 max(t i -1, 0) (6)
The coefficients c 1 and c 2 are two penalty parameters, their values depend on the importance of the corresponding term. The obtained problem (denoted by (P)) is a combinatorial optimization problem. In the next section we explain the proposed solution approach for this problem.

Solution Approach: Branch and Bound

To solve the problem, we apply a Branch and Bound method. In [START_REF] Eele | Path planning with avoidance using nonlinear branch and bound optimization[END_REF], a path planning problem avoiding circular obstacles is studied in 2D. A branching strategy is proposed, where, for each obstacle, two branches are created depending on the clockwise or counter-clockwise obstacle bypassing. We extend this branching strategy to take into account the specificity of our problem, where obstacles can be avoided also by imposing a level flight below (SID case) or above (STAR case) the obstacle. Our branching strategy is illustrated in Fig. 7. We start by setting Ω i as active (s i = 1) or not (s i = 0); when it is active, we develop three branches in order to account for the 3 possibilities of avoiding it: counter-clockwise (t i = 0), clockwise (t i = 1) or imposing a level flight (t i = 2). The lower bound for each subproblem is computed by generating the route which bypasses only the active obstacles and by calculating its length according to the objective function [START_REF] Pfeil | Optimization of airport terminal-area air traffic operations under uncertain weather conditions[END_REF]. The obtained lower bound is then used to identify whether a branch requires further subdivisions.

We present a step-by-step illustration (Fig. 8) to show how the Branch and Bound method works. The starting and ending points as well as two obstacles Ω 1 , Ω 2 are presented in Fig. 8a. We take c 1 = 1, c 2 = 1, that is to say we penalize the length of level flights in the objective function.

Step 1: We develop 4 branches on Ω 1 . We start by deviating the route counterclockwise, and we obtain a route that does not intersect Ω 2 . The lower bound in this case is 100320 m. Besides, the value of the objective function associated with solution (s 1 ,t 1 ) = (1, 0), s 2 = 0 is equal to the lower bound. Therefore, no further exploration is needed.

Step 2: Another branch on Ω 1 is developed by deviating the direct route clockwise around Ω 1 . The length of this horizontal route is the lower bound of this subproblem, the value is greater than the current best value. There is no possibility to get a better solution by further branching on Ω 2 , therefore we cut this branch.

Step 3: The third branch on Ω 1 is obtained by imposing a level flight, as shown in Figs. 8d,8e. The lower bound is greater than the current best value, so the branch is cut.

Step 4: The last branch on Ω 1 is with s 1 = 0. The route intersects Ω 2 and its lower bound, which corresponds to the length of the direct route, is less than the current best length, thus 4 branches on Ω 2 are developed. In the case s 2 = 0, the route is not feasible thus not accepted.

Step 5: By branching counter-clockwise around Ω 2 , the obtained route is still intersecting Ω 1 , thus it is not accepted.

Step 6: By branching clockwise around Ω 2 , we obtain a feasible route with length greater than the current best value, so it is not accepted.

Step 7: The last branch is obtained by imposing a level flight under Ω 2 as shown in Figs. 8i,8j. This route is still encountered by Ω 1 , so it is not accepted.

All the possible branches are considered. The best distance is 100320 m by taking

(s 1 ,t 1 ) = (1, 0), s 2 = 0.
The maximum size of the search space of our problem is 4 m , where m is the number of obstacles. Indeed we have 4 possibilities to deal with each obstacle. The size of the search space can be reduced significantly by applying the pre-processing techniques developed in [START_REF] Kim | Shortest paths for disc obstacles[END_REF]. In fact, the authors prove firstly that the shortest path connecting two points and avoiding circular obstacles must lie into an ellipse containing these two points. Moreover, they prove that the shortest path lies in a convex hull of a few circular obstacles around the line segment connecting the starting and ending points. We apply these pre-processing techniques to reduce the number of considered obstacles before applying our approach. These filters still hold when considering cylindrical obstacles. A 3-dimensional feasible route can be indeed built based on the 2D shortest path avoiding the obstacle projections on the plan. Some simulation results with and without the pre-processing techniques are presented in Sect. 4.

Simulation Results

In this section we present two simulation examples. The first one is the case of a SID, the second one is the case of a STAR. Tests were run on a Linux platform with a 2.4 GHz processor and 8 GB RAM.

Different strategies are applied to choose the next sub-problem to solve and the next obstacle to branch on. The simulation results presented in this section use the Best Lower Bound (BLB) for selecting the next sub-problem, and use the First In- There are nine obstacles. Table 1 gives the center (x i , y i ), radius r i and altitudes of the two bases

(z i in f , z i sup ) of obstacle Ω i , i = 1, • • • , 9.
The unit is in meter. The simulation results are presented in Table 2 and Fig. 9. We carried out two tests with different values of c 1 , c 2 : the first one with c 1 = 1, c 2 = 0; the second one with c 1 = 1, c 2 = 1.

When c 1 = 1, c 2 = 0, the length of level flights is not penalized in the objective function, thus Figs. 9a, 9b show that the optimal route is obtained with a level flight to avoid the obstacle with the center coordinates (29632, 25928). However, when c 1 = 1, c 2 = 1 the length of level flight is penalized in the objective function, and as a result a counter-clockwise turn is made to avoid the mentioned obstacle instead of using a level flight, as shown in Figs. 9c,9d. Moreover, thanks to the ellipse and convex hull filters, which reduces the number of the potential obstacles from 9 to 5, the computation time and the number of iterations are reduced effectively.

In the second example, the input data are (unit in meter):

• Starting point A, (x A , y A , [z A , z A ]) = (0, 0, [4454, 7722]) • Ending point B, (x B , y B , z B ) = (96304, 96304 , 914)
, c 1 = 1, c 2 = 0. b γ V , c 1 = 1, c 2 = 0. c γ H , c 1 = 1, c 2 = 1. d γ V , c 1 = 1, c 2 = 1
• α min = 2.6%, α max = 5%

• N max = 2, L min = 9260, H min = 914 There are seven obstacles defined by the parameters presented in Table 3. The simulation results are presented in Table 4 and Fig. 10. Note that the optimal route in the case with c 1 = 1, c 2 = 0 is the same as the one in the case with c 1 = 1, c 2 = 1. Even though by taking c 2 = 0, level flights are not penalized in the objective function, no level flight is imposed in the optimal solution. The reason is that the routes with level flights do not offer a better value of the objective function or do not satisfy the constraints that mentioned before. By applying the pre-processing techniques, the number of the potential obstacles is reduced from 7 to 3, therefore the computation time and the number of iterations are reduced significantly.

Conclusion and Perspectives

In this paper, we introduce a methodology for generating a 3D route in TMA at strategic level, performed by a deterministic global optimization approach. The route is represented by a horizontal curve associated to a cone in the vertical plan. We develop three ways to avoid an obstacle: turn clockwise or counter-clockwise (2D) and level flight (3D) which correspond to branching strategies in a Branch and Bound tailored on the problem. By setting appropriately the penalization coefficients, it is possible to obtain continuous and smooth routes which are available for Continuous Climb Operations (CCO) and Continuous Descent Operations (CDO).

To summarize, this approach can be regarded as a decision support tool for the designing of SID/STARs. In future work, we will consider the problem of designing several routes taking into account the avoidance of obstacles and the separation between them. One possible approach is to decompose the problem in three parts: (a) generating each route separately by applying the Branch and Bound method presented in this paper; (b) detecting the conflicts between the generated routes; (c) eliminating the conflicts by adding fictitious cylinder obstacles at the position of the conflicts.

Fig. 1 Fig. 2

 12 Fig. 1 Take-off Profiles in CDG Airport

Fig. 5

 5 Fig. 5 The Routes Associated to Different Values of the Decision Variables. as i = 0, 3D View. b (s i ,t i ) = (1, 0), 2D View. c (s i ,t i ) = (1, 1), 2D View. d (s i ,t i ) = (1, 2), 3D View

Fig. 6

 6 Fig. 6 Routes Construction. a Example 1: Horizontal Plan. b Example 1: Vertical Plan. c Example 2: Horizontal Plan. d Example 2: Vertical Plan. Note that in b and d, slopes appear discontinuous as an effect of a projection of a 3D image on a plan

  m

Fig. 7

 7 Fig. 7 Branch and Bound Branching Strategy

Fig. 8 •

 8 Fig. 8 Branch and Bound Illustration. a Step 0. b Step 1. c Step 2. d Step 3, Horizontal Plan. e Step 3, Vertical Plan. f Step 4. g Step 5. h Step 6. i Step 7, Horizontal Plan. j Step 7, Vertical Plan

Fig. 10

 10 Fig. 10 Example 2: Simulation Results. a γ H, c 1 = 1, c 2 = 0 or c 1 = 1, c 2 = 1. b γ V , c 1 = 1, c 2 = 0 or c 1 = 1, c 2 = 1

Table 1

 1 Example 1: Characteristics of Nine Obstacles

	(x i , y i ) (m)	r i (m)	(z i in f , z isup ) (m)
	(7408,14816)	3704	(701,2377)
	(9260,35100)	11112	(884,1707)
	(9260,57412)	4630	(2286,4389)
	(29632,25928)	9260	(2682,3109)
	(33336,48152)	5556	(2835,5425)
	(42596,14816)	5556	(2286,4389)
	(48152,51856)	7408	(3536,6828)
	(61116,29632)	9260	(3048,6706)
	(70376,51826)	9260	(4154,8022)

Table 2

 2 Example 1: Simulation Results

	(c 1 , c 2 )	no pre-processing Time (s) Iterations	with pre-processing Time (s) Iterations	L γ (m)
	(1, 0)	0.47	254	0.09	46	100082
	(1, 1)	0.57	307	0.20	115	102478

Fig. 9 Example 1: Simulation Results. a γ H

Table 3

 3 Example 2: Characteristics of Seven Obstacles

	(x i , y i ) (m)	r i (m)	(z i in f , z isup ) (m)
	(22224,22224)	9260	(0,5486)
	(27780,74080)	14816	(0,3048)
	(48152,96304)	7408	(0,4572)
	(59264,51856)	12964	(0,4572)
	(72228,24076)	9260	(0,4572)
	(77784,81488)	14816	(0,3048)
	(96304,48152)	12964	(0,4572)

Table 4

 4 Example 2: Simulation Results (c 1 , c 2 )

		no pre-processing Time (s) Iterations	with pre-processing Time (s) Iterations	L γ (m)
	(1, 0)	0.64	349	0.09	45	142230
	(1, 1)	0.36	188	0.09	45	142230
			a		b

Acknowledgements This work has been partially supported by Civil Aviation University of China, by National Natural Science Foundation of China (NNSFC) through grant NNSFC 61201085 and by French National Research Agency (ANR) through grant ANR 12-JS02-009-01 ATOMIC. The authors would like to thank Serge Roux from the French Civil Aviation University (ENAC), for the interesting discussions and useful comments.