N Wang

D Delahaye

M Mongeau

A Gondran

[EN-A-044] Optimal location of temporary military zones within civil aviation traffic

Keywords: optimization, military zone, location problem, simulated annealing algorithm

In this study, we focus on the problem of locating optimally temporary military zones with the aim of minimizing the number of civil flight trajectories impacted by the military activity, and the distance between the military zone and the military base. We model the military zones by 2D geometry shapes with a vertical extension associated to given flight levels during the temporary zone-activation time window. We propose a mathematical formulation of this problem as a constrained-optimization problem. We then introduce a global optimization methodology based on a simulated annealing algorithm featuring tailored neighborhood-search strategies. This is applied to one day of French traffic involving 8,836 civil flights. The results show that the proposed method is efficient to locate the military zone that is nearest from the military base, while minimizing the impact on civil flight trajectories.

Introduction

Along with the fast development of aviation technology, more and more aircraft are taking off and landing per hour; air traffic experienced its highest growth over the last five years [START_REF] Leahy | Airbus global market forecast for 2016-2035[END_REF]. The sky becomes more and more congested nowadays. Traditionally, prohibited airspace or no-fly zones are established for security reasons or military activities, where no civil flights are allowed, at all times. However, there are also needs for temporary prohibited zones to perform military training for short periods. The military zone, in this article, is a defined geographical surface above the earth. The location of this kind of military zone can be chosen dynamically according to the existing civil flights, so one minimizes the time spent by the civil aircraft in the military zone. The work in this article takes France as example, and proposes a method to locate a certain shape of military zone during a given period of time in the busy civil aviation traffic.

In France, there are more than 8,000 civil flights in one day; their trajectories cover almost the whole airspace above the French territory, as showed in Fig. 1. An example location of military zone is represented by the green rectangle in this figure, and the position of the military base is the center of the blue circle. Flights passing through this rectangle during the activation time of the military zone are affected and have to be modified. The military zone is to be relocated inside this circle at different time windows so as to minimize the impacts on the current civil flights. Moreover, from the fuel consumption saving point of view, minimizing the distance between the military zone and its military base is also desirable. For this type of problem, "quick and dirty" methods are recommended [START_REF] Moradi | Single facility location problem[END_REF]. For our problem, approximate models are used due to limited input data available for decision makers. Contrary to the classical facility problem, here the facility is considered "toxic", therefore one aims at minimizing the harmful effect (the impact on civil trajectories in our problem) inside the facility (the military zone) coverage range. Another specific feature of our problem is the fact that it deals with continuous candidate locations. The first survey on continous location problems was conducted by F. Plastria; in his work [START_REF] Plastria | Continuous location problems: Research, results and questions[END_REF], exhaustive topics on continous location theory are discussed. A recent survey on continous location-allocation problems was performed by J. Brimberg [START_REF] Brimberg | A survey of solution methods for the continuous location allocation problem[END_REF] in 2008; his work confirms that heuristic methods contributed significantly to improving the methods addressing this type of problem.

The remaining of this article is organized as follows: in Section 2, a mathematical model for the military zone optimal location problem is proposed. Section 3 introduces a simulated annealing algorithm to address this problem with dedicated neighborhood operator. Several test results on the French airspace are presented, and performance of the algorithm is analyzed in Section 4. Finally, conclusions are drawn along with future perspectives for military zone location problems at the end of this article.

Mathematical model

In this section, we model the military zone by a 2D geometry shape with vertical extension and then formulate our military zone location problem as a constrained optimization problem.

In this study, we aim to locate a military zone for one given military base, the position of this military base is given as input data. To simplify the exposition, the shape of this military zone is defined as a rectangle with a fixed vertical extension, i.e. a cuboid with fixe length L, width W and height H. It has also a fixed activation time duration D. As showed in Fig. 2, it can then be rotated and/or moved around the military base, where (ρ, θ) are the polar coordinates of this military zone (center of the rectangle), and φ is the rotation angle of the military zone. We also denote the altitude of the bottom side of this military zone as z, considering the military base as the origin, and its starting activation time as t.

To summarize, the decision variables are ρ, θ, φ, z and t. In theory, all five variables are continuous variables; In practice, aircraft follow Flight Levels (FL) in order to ensure safe vertical separation. Here we consider discrete elementary altitude shifts δz = 5 FL (1FL = 100 feet), and we consider the discrete set Z := {z min + δz, z min + 2δz, . . . , z min + Jδz} for the possible values of z, where J = (z max -z min )/δz ( z denotes the largest integer greater than or equal to z), z min and z max are respectively the lower and upper bounds of altitude constraints for the military zone. To simplify the time allocation, we also consider the discrete time shifts δt = 5 minutes, and we define the discrete set T := {t min + δt, t min + 2δt, . . . , z min + Kδz} for the possible values of t, where K = ( t max -t min )/δt , t min and t max are respectively the lower and upper bounds of activation time constraints for the military zone. Let us denote ρ max to be the maximum allowed distance between the center of the military zone and the center of the military base. Let M Z(x) denote the four-dimensional (time +space) military zone defined by an instantiation of our vector of decision variables x = (ρ, θ, φ, z, t). The given civil flight (4D) trajectories to be avoided are discretized into trajectory sample points P i ∈ R 4 , i = 1, 2, . . . , M at intervals of 15 seconds, where M is the total number of sample points for the whole trajectories. In order to minimize the impact of the military zone on the civil traffic, we aim at positioning (in 4D) the military zone so as to minimize the number of trajectory sample points intersecting it, we also wish to minimize the distance between the military zone and the military base (ρ). This is a bi-objective optimization problem that we propose to model as:

minimize x=(ρ,θ,φ,z,t) f (x) = M i=1 1 Pi∈M Z(x) + µρ (1) subject to 0 ≤ ρ ≤ ρ max (2) 0 ≤ θ ≤ 2π (3) 0 ≤ φ ≤ 2π (4) z ∈ Z (5) t ∈ T (6)
where

1 Pi∈M Z(x) is equal to 1 if P i ∈ M Z(x)
and is equal to zero otherwise. Hence

M i=1 1 Pi∈M Z(x)
corresponds to the number of sample trajectory points that lie inside the military zone, and µ is a weighting parameter whose value is to be set by the decider according to his priority between the two criteria: distance ρ vs. traffic impact. In our tests, we use µ = 1 so as to prioritize the impact on civil traffic. Remark that this model could easily be adapted so as to evaluate the traffic impact by counting rather the number of trajectories intersecting the military zone.

When we consider only minimizing traffic impact (case where the user sets µ = 0), this problem can also be considered as the dual problem of the classical maximum coverage problem involving only one set (one military zone in our case). This one-set situation removes the combinatorial nature of the maximum coverage problem. Nevertheless, instead of the typical discrete candidate sets or a planar (2-dimensional) continuous location problem, we deal here with a 5-dimensional location problem involving 3 continuous location parameters. The difficulty of this continuous optimization problem comes here from the fact that the objective function is non-differentiable and costly to evaluate (a black-box optimization problem), as counting the number of selecting points is time consuming and this computation must be performed at every iteration. This is why we propose here a heuristic resolution method instead of exact approaches.

Simulated annealing algorithm

Simulated Annealing (SA) is a typical Monte Carlo global optimization algorithm. It is a paradigm that simulates the relationship between the internal energy and the temperature T during the process of metal cooling [START_REF]Simulated Annealing: Theory and Applications[END_REF]. Metal is in liquid state at high temperature; its internal energy is at its highest level but this state is highly unstable because of the high disorder of its atomic structure. If a fast cooling procedure is applied, the temper phenomena appears: internal energy quickly drops down but does not converge towards a globally minimal value; the metal has then an amorphous atomic structure with high disorder. However, during a slow cooling process, at each temperature, the internal energy has some probability to reach a locally minimal value. When the temperature reaches room temperature, the internal energy of the metal stabilizes at a minimal value; metal has then a crystal atomic structure.

At the beginning of the cooling process, an initial state x 0 must be provided. This initial state can be randomly generated or can be particularly chosen to improve the SA performance; for instance, using prior knowledge on the problem, or using any heuristic. The temperature decrease at each iteration i is dealt with some cooling schedule such as: T i = α i T 0 , where T i is the temperature at iteration i, T 0 is the initial temperature chosen by user and α ∈ (0, 1) is a user-defined parameter. At each iteration/temperature, N t number of transitions (from one solution to another) are performed. At each transition, a neighbor state x of the current state x is generated by some neighborhood operator getN eighbor(x). This new state is then evaluated through the objective function f . If the new state is better (f (x ) ≤ f (x)), it is accepted. When the new state is worse than the current state, it is accepted with some probability proba i (x, x )related to the current temperature and the objective-function degradation. The most-used acceptance probability for a minimization problem is:

proba ( x, x ) = exp (f (x ) -f (x))
T The idea is that a bad transition is more likely to be accepted at high temperatures.

After acceptance or rejection of the new state, the algorithm proceeds to the next transition. The algorithm stops when T i goes below some pre-defined final temperature T f (for instance, we use T f = 0.0001T 0 in our tests), or some pre-defined target objective-function value is reached, or after performing N iterations, where N is the maximum number of iterations decided by user. SA is given in Algorithm 1:

Algorithm 1: Simulated annealing algorithm for a minimization problem

T ← T 0 , x ← x 0 while T > T f do for i = 1, . . . , N t transitions do x ← getN eighbor(x) if f (x ) ≤ f (x) then x ← x ; else x ← x with probability exp f (x )-f (x) Ti ; end T ← αT ; end return X;
The neighborhood operator getN eighbor(x) for our problem returns a neighbor state x (ρ , θ , φ , z , t ) of current state x(ρ, θ, φ, z, t) through the following five functions, chosen randomly with, in our tests, an identical probability, 1/5:

• T ranslate(x): return x , where ρ is chosen randomly in [0, ρ max ], and θ = θ, φ = φ, z = z, t = t.

• Rot θ (x): Rotate the military zone around the military base, return x where θ is chosen randomly in [0, 2π], and ρ = ρ, φ = φ, z = z, t = t.

• Rot φ (x): Rotate the military zone around its center, return x where φ is chosen randomly in [0, 2π], and ρ = ρ, θ = θ, z = z, t = t.

• Shif tAltitude(x): return x , where z is chosen randomly in Z, and ρ = ρ, θ = θ, φ = φ, t = t.

• Shif tActivationT ime(x): return x , where t is chosen randomly in T , and ρ = ρ, θ = θ, φ = φ, z = z.

Simulation results in the French airspace

In this section, we apply the simulated annealing algorithm on one day of French traffic involving 8,836 civil flight trajectories, discretized into 1,851,029 trajectory sample points at intervals of 15 seconds. We first present the input military zone parameters used in the numerical tests and the values we choose for the user-defined parameters involved in the SA algorithms. Then, numerical results obtained are presented and analyzed.

The input military zone parameters are listed in Table 1 The SA algorithm parameters are chosen as follows:

• Temperature decreasing parameter: α = 0.95

• Number of transitions for each temperature: N t = 200

• Initial temperature: T 0 calculated empirically at 80% acceptance rate We first apply the algorithm for three different values of the range ρ max : 100 NM, 200 NM, and 400 NM. We run the SA algorithms 10 times with the same parameter setting, average results will be computed over these 10 runs. In each case, the algorithm finds a military base with possible minimal impact on the civil flights: zero impact on the civil flights. The results obtained are presented in Table 2, where nbT rajs and nbP ts represent respectively the number of trajectories and sample points in the base range that can potentially be impacted, nbP ts * is the average minimal number of impacted trajectory sample points, ρ * is the average minimal distance of the military zone found by the SA algorithm, ∆nbP t * is the average deviation of nbP ts * for 10 runs, ∆ρ * is the average deviation of ρ * , CP U represents the average SA algorithm run time in seconds. In the objective function, the second criteria (distance from the military base is expressed in NM) is a small value in comparison with the first criterion (number of impacted points), since we chose to set the weighting parameter to µ = 1. As a consequence, we observe that the SA algorithm concentrates first on minimizing the first criteria (impact), and then pursues with minimizing the distance from the base while keeping the first criterion at its minimal value. The nearest military zone found is 51 NM away from the base for the range 100 NM. The objective function value converges faster with the increasing range limits, and the nearest military zone found is situated at 49 NM from the base for ranges 200 NM and 400 NM. The SA algorithm finds better locations when the range gets larger, because the probability of finding a location with zero impact on the civil flights increases with the range. To analyze whether the algorithm is sensitive to the shape and size of the military zone, we test different shapes and sizes of the military zone by modifying the input parameters: length, width, height and activation duration. The results obtained by SA for shape change and size reduction of the military zone are similar to those obtained above. However, if we increase the size of the military zone, SA may fail to find a military zone with zero impact on civil trajectories under small range limitations. Table 3 displays the results obtained for example when doubling the length of the military zone: ∆nbP t * and ∆ρ * are relatively large for the range of 100 NM since for 3 runs out of 10, the algorithm finds a military location near from the base but impacting hundreds of trajectory sample points.

• Final temperature: T f = 0.0001T 0

Conclusion and perspectives

As the air traffic becomes more and more crowded, it becomes more difficult to perform military activities in the air without impacting civil flights. In this article, we studied the location problem of a temporary military zone within civil aviation traffic dedicated for short-time military activities. The military zone was modeled by a 2D geometry shape, with vertical extension associated to given specified flight levels, that has the flexibility to move in space and time. A constrained-optimization formulation was proposed, it aims at minimizing both the impact of military activities on civil flights, and the distance between the military zone and the military base. A simulated annealing algorithm was designed to address this black-box optimization problem. Numerical experiments were conducted on instances involving one day of French traffic. The results indicates that the proposed methodology is a viable decision-aid tool to search for the nearest military zone having no impact on the traffic, as in our tests none of the existing 8,836 civil flights were impacted under different range limitations. Various tests on different input parameters show that the SA algorithm is robust and not sensitive to the shape and size changes of the military zone.

In future work, one may address the dynamic version of this location problem, in which the position of the center of the military zone may change with time during the given activation time window. This is relevant for instance when one deals with the military zone encompassing an in-flight supply plane. Another promising track of research could envisage using a deterministic black-box optimization methods such as those proposed in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF][START_REF] Audet | Derivative-Free and Blackbox Optimization[END_REF].

Finally, the military location problem may be combined with a real-time flight-trajectory conflict prediction and resolution problem. In this case, conflict resolution approaches will be applied on the trajectories impacted by the chosen military zone.

Figure 1 :

 1 Figure 1: Example of a military zone location (green rectangle) near the military base (center of blue circle) within civil flight trajectories (red lines) over France

Figure 2 :

 2 Figure 2: Military zone location model

Figure 3 :

 3 Figure 3: Evolution of the two criteria for ρ max = 100 NM

Figure 4 :Fig. 3 ,

 43 Figure 4: Evolution of the two criteria for ρ max = 200 NM

Figure 5 :

 5 Figure 5: Evolution of the two criteria for ρ max = 400 NM

Table 1 :

 1 : Input military zone parameters

	parameter	value	unit
	length (L)	40	NM
	width (W )	20	NM
	height (H)	5,000	feet
	zmin	10,000	feet
	zmax	30,000	feet
	δz	1,000	feet
	J	20	-
	D	120	minutes
	tmin	10 AM	-
	tmax	12 AM	-
	δt	5	minutes
	K	24	-

Table 2 :

 2 Simulation results with length L = 40NM with different values of ρ max

	ρmax (NM)	100	200	400
	nbT rajs	1339	1778	2113
	nbP ts	100,741 171,331 226,316
	nbP ts *	0	0	0
	∆nbP t *	0	0	0
	ρ * (NM)	53.2	51.1	50.8
	∆ρ * (NM)	2.3	3.8	1.4
	CP U (s)	114	299	285

Table 3 :

 3 Simulation results with length L = 40NM with different values of ρ max

	ρmax (NM)	100	200	400
	nbT raj	1,561	1,932	2,117
	nbP t	131,272 188,288 229,178
	nbP t *	62.4	0	0
	∆nbP t *	88	0	0
	ρ * (NM)	83.8	104	93.1
	∆ρ * (NM)	18.8	2.4	0.18
	CP U (s)	281	843	596

Acknowledgments

This work has been supported by Civil Aviation University of China and the program of China Scholarships Council.

Copyright

The authors of this article confirm that they and their institution, hold copyright of all original material included in their paper. They also confirm they have obtained permission, from the copyright holder of any third party material included in their paper, to publish it as part of their paper. The authors grant full permission for the publication and distribution of their paper as part of the EIWAC2017 proceedings or as individual off-prints from the proceedings.


Arianit Islami
email: arianit.islami@caa-ks.org

Supatcha Chaimatanan
email: supatcha@gistda.or.th

Moaxi Sun

Daniel Delahaye
email: daniel.delahaye@enac.fr

Optimization of military missions impact on civilian 4D trajectories

Keywords: Military area, optimization, mission trajectory, airspace

Flying is the fastest and one of the safest methods of transport. Air traffic is increasing constantly and according to statistics from ICAO and aviation industry, approximately every 15 years traffic is doubling and as the consequence the airspace becomes more and more congested. In order to cope with the increase and to improve safety, air traffic management system will rely on a so-called trajectory based operations concept that will increase air traffic capacity by reducing the controllers workload. In addition to the increase of traffic of civil aviation, we have to take into consideration requirements from military operations and their needs for missions in the airspace. In order to fulfill military mission requests, we shall minimize its impact into civil aviation whenever possible. This will be achieved by transferring tactical conflict detection and resolution tasks to the planning phase. In this future air traffic management paradigm context, this paper presents a methodology to address such military missions planning. The proposed methodology aims at minimizing the interaction of trajectories of military missions into civilian trajectories. Results that we gain through this methodology are by allocating alternative departure times, alternative horizontal flight paths, and alternative flight levels to the military missions involved in the interaction. This paper presents a mathematical formulation of this mission trajectory planning problem leading to a mixed-integer optimization problem, whose objective function relies on the new concept of interaction between trajectories. A computationally efficient algorithm to compute interaction between military missions trajectories with minimal impact into civilian trajectories for large-scale applications is presented. Resolution method based on simulated annealing algorithm have been developed to solve the above optimization problems. Finally, the methodology is implemented and tested with military missions planed to occur in real air traffic data taking into account uncertainty over the French airspace. Mission conflict-free and robust 4D trajectory planning are produced within computational time acceptable for the operation context, which shows the viability of the approach.

Historically, State agreements between military aviation units and Air Navigation Service Provider (ANSP) have focused on the needs of State defense, security and emergency procedures as well as military readiness and response requirements. There is now a clearly defined need to establish techniques that support the efficient integration of military and civil aviation in day-to-day operations. Aviation operations of all types contribute significantly to the economy of a State, and, as such, their growth needs to be protected and encouraged. The Flexible Use of Airspace (FUA) concept provides that airspace should no longer be des-ignated as either military or civil airspace but should be considered as one continuum and used flexibly on a day-to-day basis, in which all user requirements are accommodated to the greatest possible extent. [START_REF]Civil/military cooperation in air traffic management[END_REF] In light of this, a need for the real-time use of airspace allowing a safe separation between civil and military aircraft has been clearly identified. In this respect, European Commission has issued regulation in order to reiterate Flexible Use of the Airspace (FUA). To this end, Eurocontrol has provided Airspace Concept and Strategy for the ECAC Area, which provides targets for airspace developments. In this document, the concept of Dynamic Airspace Allocation (adapt- able airspace allocation moving in time and place) is required to respond in real-time to changing situation in traffic flows and to short-term changes of users intentions.

The European air traffic management modernization program SESAR also proposes the concept of Dynamic Mobile Area (DMA) that is described by integral part of the Mission Trajectory (MT) described by a 4D data set. The DMA defines the volume of airspace that satisfies specific requirements from different Airspace Users. Example of the DMA for aerial refueling mission is illustrated in Fig. 1.

In this paper, we are addressing the optimization of a mission trajectory and its corresponding DMA, in order to minimize its interaction with the civilian trajectories.

Air traffic regulations impose that aircraft must always be separated by some prescribed distance. Aircraft are considered to be in conflict when these minimum separation requirements are violated. As the global air traffic demand keeps on increasing, congestion problem becomes more and more critical. One of the key solutions is to balance the air traffic demand and the overall capacity of the civil/military Air Traffic Management (ATM) system. In order to cope with the increasing demand, the future ATM system will rely on the trajectory based operations concept. In this concept, aircraft will be required to follow a negotiated conflict-free trajectory, accurately defined in 4 dimensions (3 spatial dimensions and time) in order to reduce the need of controller's intervention during the tactical phase. In this perspective, the key factor to improve the ATM capacity and to accommodate military needs is an efficient strategic 4D trajectory planning methodology to compute a conflict-free 4D mission trajectory for each aircraft.

In this work, we propose a methodology to address such a strategic planning of trajectories at national scale. The goal of the proposed method is to minimize interaction (in both the three dimensional space and in the time domain) between the MT and civil trajectories, by allocating an alternative horizontal flight path, departure time, and flight level to the MT. An interaction between trajectories occurs when two or more trajectories have an effect on each other; for instance, when trajectories occupy the same space at the same period of time. Therefore, contrary to the concept of conflict, the measurement of interaction does not only refer to the violation of minimum separation requirements. It also allows us to take into account other separation criteria, such as separation between the civil trajectory and the DMA.

The following sections of this paper is organized as follows. Section 2 presents mathematical model of the proposed methodology. Section 3 describes resolution algorithm based on simulated annealing. Implementation of the resolution algorithm to our proposed mission trajectory planning is presented in Section 4. Section 5 reports the results from computational experiments. Finally, conclusion is discussed in Section 6.

Mathematical Modeling

In order to develop a resolution algorithm, the first stage in the optimization process consists of modeling the real problem using a mathematical abstraction, which is as faithful as possible. Using this abstraction, it becomes possible to develop solution algorithms which can be executed on a computer. The modeling stage, then, consists of characterizing the state space and the objective space.

Data

Before establishing the mathematical model, one must first consider the given data associated to the problem. We consider a given airspace and a mission trajectory (MT) that have to be designed between two points A and B. The shape and the size of the dynamic mobile area (DMA) that has to be routed between A and B is also given(see Fig. 2). In the following, this area will be called DMA. The center of the first one is located at 1 3 and the second one 1 3 .

In a given airspace we consider a set of civilian aircraft trajectories (red traffic on Fig. 2) that will be noted Γ. One element of this set will be noted γ i (t) (i ∈ [1, N] and N = |Γ|).

In this preliminary work, a direct route has been considered but the approach could easily extended to more complex mission trajectories. The current, the minimum and the maximum departure time of the mission will be noted ts 0 , ts min , ts max . We will consider a constant speed for the mission (v 0 ). Similarly, the current, the minimum and the maximum flight level of the mission will be noted FL 0 , FL min , FL max .

State Space

The mission trajectory will be noted γ DMA (t). We consider only en-route airspace and we propose to design this trajectory in terms of route (horizontal flight path) design, flight level and time of departure assignment. This represent the state space of our problem. In this problem, three decision variables will be considered :

1. Time of departure of the mission :

t s ∈ [ts min , ts max ].

Flight level of the mission :

FL ∈ [FL min , FL max ].

Shape of the mission trajectory

: γ DMA .
Common practice conducted us to rely on a discretization of the departure time t s and the flight level FL of the mission. Let δ t be the departure-time-shift step size such that t s = ts 0 + δ t . This yields N t = (ts maxts min + 1)/δ t possible departure time of the mission. Similarly, let δ FL be the flight level shift, this yields N FL = (FL max -FL min + 1)/δ FL possible flight level for the mission.

In our model, the shape of the mission trajectory is designed thanks to two waypoints (W 1 ,W 2 ) which are added along the route (see Fig. 3).

In order to locate such waypoints, we first build a unit segment (length = 1) for which we consider two rectangular areas (see figure 4). In each of those rectangles, we put one waypoint (w 1 and w 2 ) in order to design the initial shape. The size of the rectangle in the vertical dimension is fixed in order to limit the maximum deviation of the mission trajectory in terms of distance.

This initial shape is then scaled to the length of the segment [A, B] and rotated by an angle θ corresponding to the direction of the segment [A, B] (see figure 5).

For this rotation, we use the classical rotation matrix for computing the new coordinated of the waypoint in the airspace (W 1 ,W 2 ) :

x new y new = cos θ sin θ -sin θ cos θ x old y old A B θ Figure 5
Computed shape of the mission after scaling and rotation.

Objective

The first objective of the mission planing is to minimize the impact of the military area (DMA) on the civilian traffic.

Let γ DMA (t) be the position of the center of the DMA at time t, one must consider the civilian trajectories in intersection with the such DMA at that time. We note DMA(γ DMA (t)) the area for which the center is located at γ DMA (t). In this preliminary work, the DMA with is defined by a rectangular box with width (w DMA ), length (l DMA ), and height (h DMA ) to be set by the user.

If we note y(t) as a number of civilian trajectories in intersection with the DMA area, then :

y(t) = i = 1 i=N γ i (t) ∈ DMA(γ DMA (t))
If we note t e the time when the mission end (the DMA area reach the point B), the we can compute the overall objective function :

y = t=t e t=t s y(t)
This objective has to be minimize. Example of interaction between civilian traffic and the mission trajectory is illustrated in Fig. 6. Finally, the proposed mission trajectory optimization problem can be formulated as a mixed-integer optimization problem as follows: min t s ,FL,w y subject to

t s ∈ [ts min , ts ], FL ∈ [FL min , FL max ], w i ∈ W i ,
where w is the vector of waypoint w i for i = 1, . . . , M, where M is the number of waypoints that that are allowed to be used, and W i is the set that define feasible locations of the waypoints w i for i = 1, . . . , M.

3 Resolution algorithm: Simulated Annealing Simulated Annealing (SA) is one of the simplest and best-known meta-heuristic methods for addressing the difficult black box global optimization problems (those whose objective function is not explicitly given and can only be evaluated via some costly computer simulation). It is massively used on real-life applications.

In the early 1980s three IBM researchers, Kirkpatrick, Gelatt and Vecchi [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], introduced the concepts of annealing in combinatorial optimization.

These concepts are based on a strong analogy with the physical annealing of materials. This process involves bringing a solid to a low energy state after raising its temperature. It can be summarized by the following two steps:

• Bring the solid to a very high temperature until "melting" of the structure;

• Cooling the solid according to a very particular temperature decreasing scheme in order to reach a solid state of minimum energy.

In the liquid phase, the particles are distributed randomly. It is shown that the minimum-energy state is reached provided that the initial temperature is sufficiently high and the cooling time is sufficiently long. If this is not the case, the solid will be found in a metastable state with non-minimal energy; this is referred to as hardening, which consists in the sudden cooling of a solid.

In 1953, three American researchers (Metropolis, Rosenbluth, and Teller [START_REF] Metropolis | Equation of state calculation by fast computing machines[END_REF]) developed an algorithm to simulate the physical annealing. Their aim was to reproduce faithfully the evolution of the physical structure of a material undergoing annealing. This algorithm is based on Monte Carlo techniques which consist in generating a sequence of states of the solid in the following way.

Starting from an initial state i of energy E i , a new state j of energy E j is generated by modifying the position of one particle.

If the energy difference, E i -E j , is positive (the new state features lower energy), the state j becomes the new current state. If the energy difference is less than or equal to zero, then the probability that the state j becomes the current state is given by:

Pr{Current state = j} = exp E i -E j k b .T ,
where T represents the temperature of the solid and k B is the Boltzmann constant (k B = 1.38×10 -23 joule/Kelvin). The acceptance criterion of the new state is called the Metropolis criterion.

In the SA algorithm, the Metropolis algorithm is applied to generate a sequence of solutions in the state space S . To do this, an analogy is made between a multi-particle system and our optimization problem by using the following equivalences:

• The state-space points represent the possible states of the solid;

• The function to be minimized represents the energy of the solid.

A control parameter c, acting as a temperature, is then introduced. This parameter is homogeneous to the criterion that is optimized.

It is also assumed that the user provides for each point of the state space, a neighborhood and a mechanism for generating a solution in this neighborhood. We then define the acceptance principle:

The acceptance criterion for accepting solution j from the current solution i is given by the following probability:

Pr{ accept j} = 1 if f ( j) < f (i) exp f (i)-f ( j) c else.
Let i be the current state space, j be the neighbor of i defined by a given neighborhood function. And let k be the temperature step, L k be the number of iterations to be performed at each temperature step, and c k be the temperature at the k th step. The simulated annealing can be summarized the following way : One of the main features of simulated annealing is its ability to accept transitions that degrade the objective function.

At the beginning of the process, the value of the temperature c k is high, which makes it possible to accept transitions with high criterion degradation, and thereby to explore the state space thoroughly. As c k decreases, only the transitions improving the criterion, or with a low criterion deterioration, are accepted.

Finally, when c k tends to zero, no deterioration of the criterion is accepted.

Simulated annealing has been applied to many highly combinatorial problems coming from industry and operations; to mention a few:

• Airline crew scheduling [START_REF] Emden-Weiner | Best practice simulated annealing for the airline crew scheduling problem[END_REF] • Railway crew scheduling [START_REF] Hanafi | A hybrid constructive heuristic and simulated annealing for railway crew scheduling[END_REF] • Traveling salesman problem [START_REF] Bayram | A new simulated annealing approach for travelling salesman problem[END_REF] • Vehicle routing problem [START_REF] Mahmudy | Improved simulated annealing for optimization of vehicle routing problem with time windows (VRPTW)[END_REF] • Layout-routing of electronic circuits [?]

• Large scale aircraft trajectory planing [START_REF] Chaimatanan | A hybrid metaheuristic optimization algorithm for strategic planning of 4D aircraft trajectories at the continental scale[END_REF][START_REF] Islami | Large-scale 4D trajectory planning[END_REF] • Complex portfolio problem [START_REF] Crama | Simulated annealing for complex portfolio selection problems[END_REF] • Graph coloring problem [START_REF] Chams | Some experiments with simulated annealing for coloring graphs[END_REF] • High-dimensionality minimization problems [START_REF] Siarry | Enhanced simulated annealing for globally minimizing functions of many continuous variables[END_REF] Having introduced the simulated annealing, the next section presents the implementation of such algorithm to our optimization problem.

Algorithm Implementation

The overall structure of the resolution algorithm can be summarized by Fig. 7. It consists of two main components which are:

• Simulation environment module. Given input data and a candidate solution (state space), the simulation compute the corresponding mission trajectory and DMA. Then, place them in the airspace. After that it computes the value of objective function.

• Optimization module control the search for optimal solution. It evaluate the decision variables, and apply acceptation criteria. This process is iterated until the maximum number of interaction or the stopping criteria is fulfilled.

The algorithm proposes some decision variables (time of departure t s , flight level FL, and route of the mission γ) as a candidate solution. Then, it is used to build the mission trajectory that will be placed in the airspace in order to compute the associated objective function.

The process to simulate the mission trajectory and the DMA in the current airspace can be described as follows. Once the point of the state space ( X) is defined :

• the time of departure of the mission : ts . The algorithm locate the waypoint in the airspace by scaling and rotating the waypoint w 1 , w 2 in order to create W 1 , W 2 . Based on those airspace waypoint an can build the route of the mission. Having the route, the time of departure (ts) the flight level ( f l) and the speed of the mission (v 0 ), on can compute the mission trajectory (γ dma (t)), which is sampled every 15 seconds. This represents the trajectory on the center of the DMA area. This trajectory is then included in the airspace for evaluation.

In order to run the simulated annealing algorithm, a neighborhood operator is needed. For our application, the algorithm randomly select one of the feature of the state space (time of departure, flight level or route) and apply the following changes :

• Time of departure : ts new = ts old +δ; delta being a random variable; one must ensure that ts new ∈ [ts min , ts max ].

• Flight Level : f l new = f l old + u where is a random variable having values in -1; +1; one must ensure that f l new ∈ [ f l min , f l max ].

• Route : a waypoint w i to be modified is first selected. The position of this waypoint is randomly changed in the associated rectangle W i .

Having developed all the materials for implementing the simulated annealing algorithm to our problem, the next section presents the associated results to some realistic instances of the problem.

RESULTS

The simulated annealing algorithm adapted to solve the mission trajectory optimization problem is implemented in Java.

It is tested with two different missions and traffic sets over the French airspace. The overall methodology is implemented on a UNIX platform with 2.40GHz 4 processor and 20 GB DDR4 RAM (personal computer).

The algorithm has been first applied to the French airspace with simulated civilian traffic of 8,836 flights considering en-route traffic. Each trajectory being sampled every fifteen seconds, each sampled point is defined by four coordinates (x, y, z, t). This civilian traffic set is illustrated as blue lines on Fig. 8. The traffic is projected to Euclidean plane, and then normalized using separation standard. The initial mission trajectory and the DMA are set as follows:

• The mission is originated at coordinates (40,351) and destination is (13, 462) (in the normalized unit).

• The initial flight level, FL 0 , is set to FL350;

• The initial departure time, ts 0 is set to 04:10:00;

• The size of the DMA is as: l DMA = 45 NM, w DMA = 25 NM, and h DMA = 5,000 feet.

To solve this problem, we set the parameters that define the problem as given in Table 1, and empirically set the parameters that define the simulated annealing algorithm as given in Table 2.

The initial mission trajectory involves in interactions with 167 civilian flights (y init = 167). The resolution algorithm was able to find interaction-free solution within 48 seconds. The interaction-free mission trajectory is shown as green line in Fig. 8.

The evolution of the value of objective function for this data set is shown in Fig. 9. One can observe that in the beginning, the algorithm accepts solution that yields higher value of objective function. Then, this acceptation becomes less and less, until the algorithm is reached towards interaction-free solution.

Then, the resolution algorithm is applied to a larger set of simulated civilian traffic consists of 20,734 flights as illustrated in Fig. 10. The initial mission trajectory and the DMA are set as follows:

• The mission is originated from (44 • N,2 • E). The destination is set to (55 • N,2 • E).

• The initial flight level, FL 0 , is set to FL370;

• The initial departure time, ts 0 is set to 08:30:00;

• The size of the DMA is as: l DMA = 60 NM, w DMA = 40 NM, and h DMA = 4,000 feet.

The initial mission trajectory interacts with 732 civilian flights (y init = 732). To solve this problem, we use the same parameters that defined the problem and the resolution algorithm as given in Table 1 and 2, except the number of iteration L k is set to 1,000.

Again, the SA was able to find interaction-free solution for this larger set of civilian traffic within 289 seconds. The interaction-free mission trajectory is shown as thick yellow line in Fig. 11. And the evolution of the objective function value is shown in Fig. 12. Similar to the previous case,the algorithm accepts solution that yields higher value of objective function in the beginning, and then, the rate of acceptation becomes less and less, until the algorithm is reached towards interaction-free solution.

CONCLUSION

We have introduced a methodology to address mission trajectory planning in the framework of future airspace management concept that focus on increasing flexibility of the airspace usage by allocating airspace more dynamically, while satisfying the needs of both military and civilian airspace users.

The proposed mission trajectory planning methodology focuses on minimizing the interaction between military mission and civilian air traffic. The interaction is minimized by allocating alternative departure time, flight level, and horizontal flight path to the military mission.

The proposed model is formulated as a mixedinteger optimization problem, for which a resolution algorithm based on simulated annealing has been developed.

Computational experiments on a day instance of traffic over the French airspace, involving up to ≈ 20,000 flights, shows that the proposed methodology is able to find interaction-free mission trajectory.

[

  EN-A-074] Optimization of military missions impact on civilian 4D trajectories + A. Islami * M. Sun * + S. Chaimatanan * * D. Delahaye * 1 Introduction

Figure 1

 1 Figure 1 Example of DMA concept applied for tanker's mission trajectory.

Figure 2

 2 Figure 2 Initial mission trajectory connecting origin A to destination B in a given airspace.The moving area is represented by the green rectangle on the figure.

Figure 3

 3 Figure3The shape of the mission trajectory is changed by using two waypoints (w 1 ,w 2 ).

Figure 4

 4 Figure 4 Waypoints are located into two rectangular areas. The center of the first one is located at1 3 and the second one1 3 .

  Interaction between MT and civilian trajectories occurs when there are intersections between DMA and civilian trajectories.

Figure 6

 6 Figure 6 Example of interaction between two civilian trajectories with mission trajectory.

Simulated annealing 1 . 3 .; 4 . k := k + 1 ; 5 . 6 .

 134156 Initialization i := i start , k := 0, c k = c 0 , L k := L 0 For l = 0 to L k do • Generate a solution j from the neighborhood S i of the current solution i; • If f ( j) < f (i) then i := j ( j becomes the current solution); • Else, j becomes the current solution with probability e f (i)f ( j) c k Compute(L k , c k ); Until c k 0

•Figure 7

 7 Figure 7 Objective-function evaluation based on a simulation process

Figure 8

 8 Figure 8 Simulated France traffic consists of 8,836 flights, initial mission trajectory, and optimized mission trajectory.

Figure 9

 9 Figure 9 Evolution of the objective function value for data set consists of 8,836 civilian flights.

Figure 10 Simulated

 10 Figure 10 Simulated France traffic consists of 20,734 flights and initial mission trajectory (in thick yellow line).

Figure 11 Simulated

 11 Figure 11 Simulated France traffic consists of 20,734 flights and initial mission trajectory (in thick yellow line).

Figure 12

 12 Figure 12 Evolution of the objective function value for data set consists of 20,734 civilian flights.

Table 1

 1 Parameters that define the problem.

	Parameters	Value
	Minimum flight level FL min	FL 0 -20
	Maximum flight level FL max	FL 0 + 20
	Flight level shift δ FL	10
	Minimum departure time ts min	ts 0 -8 minutes
	Maximum departure time ts max	ts 0 + 8 minutes
	Departure time shift δ t	1 minutes
	Maximum route length extension 12%
	Parameters	Value
	Number of iterations L k	600
	Initial temperature c 0	3,284.2
	Temperature decreasing scheme c k+1 = 0.99 * c k
	Stopping criteria	c k < 0.0001 * c 0

Table 2

 2 Parameters that define the SA.


Ji Ma

Daniel Delahaye

Mohammed Sbihi

Paolo Maria Scala

Mujica Mota

+ J Ma

M Mujica Mota

A Study of Tradeoffs in Airport Coordinated Surface Operations

Keywords: Airport surface operations, Global optimization, Taxiways routing, Runways scheduling

Airports represent the major bottleneck in the air traffic management system with increasing traffic density. Enhanced levels of automation and coordination of surface operations are imperative to reduce congestion and to improve efficiency. This paper addresses the problem of comparing different control strategies on the airport surface to investigate their impacts and benefits. We propose an optimization approach to solve in a unified manner the coordinated surface operations problem on network models of an actual hub airport. Controlled pushback time, taxi reroutes and controlled holding time (waiting time at runway threshold for departures and time spent in runway crossing queues for arrivals) are considered as decisions to optimize the ground movement problem. Three major aspects are discussed: 1) benefits of incorporating taxi reroutes on the airport performance metrics; 2) priority of arrivals and departures in runway crossings; 3) tradeoffs between controlled pushback and controlled holding time for departures. A preliminary study case is conducted in a model based on operations of Paris Charles De-Gaulle airport under the most frequently used configuration. Airport is modeled using a node-link network structure. Alternate taxi routes are constructed based on surface surveillance records with respect to current procedural factors. A representative peak-hour traffic scenario is generated using historical data. The effectiveness of the proposed optimization methods is investigated.

INTRODUCTION

With the steady growth of air traffic, the current air network is facing capacity problems, leading to delays and congestions. One of the most critical parts is the airport and its surrounding airspaces. Increasing use of saturated airfield capacity will adversely impact predictability and punctuality. European SESAR (Single European Sky ATM Research) program [START_REF] Eurocontrol | European atm master plan[END_REF] and FAA's NextGen (Next Generation Air Transportation System) plan [START_REF]Nextgen implementation plan 2016[END_REF] aim to increase the network traffic throughput in order to accommodate all the forecast demand with a sufficient margin. To achieve this goal, new operational concepts and techniques need to be developed to support the increased traffic density. Efficient planning and optimization approaches of airport operations are critical to alleviate traffic congestions.

Airport operations involve ground movement [START_REF] Atkin | The airport ground movement problem: Past and current research and future directions[END_REF], runway sequencing and scheduling [START_REF] Bennell | Airport runway scheduling[END_REF], gate assignment [START_REF] Bouras | The airport gate assignment problem: a survey[END_REF], etc. Segregated researches on these do-mains have been conducted in the past years and have been proven to improve safety and efficiency. Recently, integrated study of these sub-problems are in trend since they are intimately linked and affected by one another. Holistic optimization can gain potential benefits and target a better synchronization. More and more, large-scale complex hub airports during peak hours are studied instead of limited toy example. Many works based on deterministic and stochastic optimization approaches were proposed. Mixed Integer Linear Programming (MILP) formulation is a deterministic approach usually used in this problem: In [START_REF] Balakrishnan | A framework for coordinated surface operations planning at dallas-fort worth international airport[END_REF], controlled pushback and taxi reroutes concepts are considered to minimize the total taxi time. [START_REF] Guépet | Optimization of airport operations: stand allocation, ground routing and runway sequencing[END_REF] addresses the integration of runway sequencing problem and ground routing problem including conflicts of the ramp area. As for stochastic approach, Genetic Algorithm (GA) has exclusively used. Gotteland et al. [START_REF] Gotteland | Optimisation du trafic au sol sur les grands aéroports[END_REF] presented a hybrid algorithm combining GA and branch and bound to solve ground movement prob-lem. Deau et al. [START_REF] Deau | Airport surface management and runways scheduling[END_REF] extended their previous work by considering runway sequencing problem and making it more consistent with ground movement. Recently, Weiszer et al. [START_REF] Weiszer | An integrated optimisation approach to airport ground operations to foster sustainability in the aviation sector[END_REF] introduced a multi-objective genetic algorithm to solve ground operation and departure runway scheduling problem. In general, combining the airside and ground problem and optimize together can gain more benefit. However, the complexity of the integrated problem would grow significantly as well.

In this work, we propose a methodology to address the problem of comparing different control strategies on the airport surface to investigate their impacts and benefits. An global optimization approach is used to solve in a unified manner the coordinated surface operations problem on network models of an actual hub airport. Controlled pushback time, taxi reroutes and controlled holding time are considered as main decisions to optimize the ground movement problem and the runway scheduling problem. Our first contribution is to analyze real alternate taxi routes used on a complex airport configuration. Then an integrated global model taking into account pushback, taxi routes, holding point, holding time, runway crossing and departure sequence at the same time is proposed. Several levers give the possibility to compare different scenarios.

The reminder of this paper is structured as follows. Section 2 describes the problem and analyzes the taxi routes using radar traffic data. Section 3 models the airport ground operation problem. Section 4 presents the solution approach. Section 5 performs tests and analyzes the results. Section 6 gives some conclusions and perspectives.

PROBLEM DESCRIPTION

Airport ground optimization aims at finding the best schedules and routes in order to minimize delays and maximize airfield capacity taking into account several constraints: taxiing separation, route choices, wake turbulence separation for landings and take-offs etc. A departure flight starts its taxi process with pushback at gate. Then it follows the assigned taxi route to reach the departure runway entry point and takes off. An arrival flight lands and exits the runway, then taxis to the assigned gate. Based on different configurations of airport, an arrival may cross a departure runway, and vice versa. In the following section, we give first an introduction about our airport model. Then we explain how to construct alternate taxi routes based on surface surveillance records with respect to current procedural factors. These alternate taxi routes will be used later in the optimization model.

CDG airport model

We choose to study Paris CDG airport because of its complexity and accessibility to the data. CDG air-port is one of the busiest passenger airports in Europe, composed of four parallel runways (two for landings and two for departures) and three terminals. In CDG, Ground controllers handle all intermediary taxiing routes. Local controllers and Apron controllers handle respectively the runway area and parking areas [START_REF] Chua | Simulating air traffic control ground operations: Preliminary results from project modern taxiing[END_REF]. Due to this different area classification and in order to simplify the problem, our model considers that taxiway starts with a defined meta-gate shown in Fig. 1, which is the exit point of the ramp area and the entry point of the taxiway area, and ends with runway entry point for departures. Ramp area is beyond the scope of this paper. For arrivals, taxiing path starts with runway exit point, and ends in meta-gate. We model CDG airport with a graph G = (N, L), where N and L represent the nodes set and links set respectively. Each node can be a runway entry/exit point, a holding point, an intersection or a meta-gate. Each link is composed of two nodes. We have in total 392 nodes and 617 links.

Ground traffic data preprocessing

In previous literature, three possible routing options are most used in the aircraft taxi problem: single path, alternate path and free path [START_REF] Guépet | Optimization of airport operations: stand allocation, ground routing and runway sequencing[END_REF]. In the first case, aircraft follow a predetermined taxi route, which is usually the standard route in the airport. In the second case, several routing options are proposed after applying, for instance, the k-shortest path algorithm [START_REF] Martins | A new implementation of yens ranking loopless paths algorithm[END_REF]. In the last case, any routes can be assigned to an aircraft.

In the operational point of view, most of the aircraft take standard taxi route with a preferential sens specified in the operations manual. Sometimes controllers deviate the predefined taxi route of one aircraft in order to avoid potential conflicts between flights or to forbid blocked or restricted areas. The alternate route choices obey some potential rules (e.g., not taking a long and unnecessary detour). Considering the taxiway configuration of Roissy airport, alternate path seems to be an appropriate option to formulate the problem in order to comply with reality. In the previous related work, alternate paths are generated by applying a classical k-shortest path algorithm in the graph network of airport [START_REF] Gotteland | Optimisation du trafic au sol sur les grands aéroports[END_REF]. Cost represents the travel time on one link, possibly augmented to avoid crossing some runway areas. Additional adjustments are applied to avoid passing two times the same node. Thereafter, feasible route options can not be simply defined as a set of validated shortest paths, because the distances from one to another are too small, which is algorithmically correct but not applicable in practice. Distinct alternate routes need to be found. Moreover, the value k should not be the same for each pair of origin and destination. A gate in the north side which is close to a north side runway will have less route options than a gate in the south side which is far from this runway. Based on these practical requirements, we decided to extract alternate routes sets by analyzing airport's flight radar records to find the operationally used potential routes set. To the best of our knowledge, this is the first time that alternate taxi routes are generated with radar data.

First, we define all the possible intersections of the taxi network, each with a circle centered on the intersection and with a radius to cover all the radar tracks. Remark that at the runway entry/exit point, aircraft moves faster, thereby our record radius must be large enough to capture all the aircraft plots. Thereafter, for each radar track, we record its route as a series of nodes, starting with runway exit point, ending with a meta-gate for arrivals and starting with a meta-gate, ending with runway entry point for departures. Fig. 2 shows an simplified example. We have one runway exit R1, one meat-gate M1 and 16 intersection nodes presented in Fig. 2a, the radar tracks are illustrated in red lines. The possible route options are { R1, 1, 4, 10, 14, 16, M1 } (Fig. 2b) with 10 aircraft passing, { R1, 1, 3, 5, 8, 15, 16, M1 } (Fig. 2c) with 5 aircraft passing and { R1, 1, 3, 5, 8, 15, 14, 10, 4, 5, 8, 15, 16, M1 } (Fig. 2d) with 1 aircraft passing with regard to the positions and time of radar tracks. Then we sort all the possible routes in a descending order by the number of flights going through the set of nodes. It has not escaped our notice that the route traversed by only one aircraft is usually in an abnormal case, see example on Fig. 2d. Therefore, it should not be selected in the potential route set. At last, for pair (R1, M1), we obtain two alternate routes: { R1, 1, 4, 10, 14, 16, M1 } and { R1, 1, 3, 5, 8, 15, 16, M1 }. To summarize, in order to generate alternate taxi routes, we proceed:

• Step 1: Generate routes of flights -For each flight, find the nodes for which the radar track passes inside the detection zone; -Obtain the route by chronologically sorting these nodes; -Count the number of flights using the same route, sorted in descending order. After analyzing 13 days of real traffic (February 2016), we generate all the feasible taxiing route sets for the west configuration in CDG. Fig. 3a illustrates an example of 9 possible route options from one north runway exit to a south meta-gate. In Fig. 3b, we found only one route option from the meta-gate to the runway entry, which can be explained by the short distance between origin and destination. However, in total 309 aircraft follow this route. We have in total 510 combinations of different pairs (runway meta-gate). In most cases we have only one standard route. Besides, other options exist between 2 and 9 routes. Few pairs possess more than 6 options. Table 1 lists the number of pairs admitting k routes option (k = 1, ..., 9).

MATHEMATICAL FORMULATION OF THE PROBLEM

In this section, we describe an integrated global optimization model for the airport ground operations problem. We first give flights input data. Next, decision variables are defined. Then, we clarify constraints. At last, an objective function is introduced.

Input data

We have a set of flights F = A D, where A denotes the set of arrival flights and D denotes the set of departure flights. For each f ∈ F, the following input data are given:

• C f : wake turbulence category;

• M f : meta-gate;

• E f : runway entry point for departure or runway exit point for arrival;

• P 0 f : initial off-block time for departure;

• L f : initial landing time for arrival;

• H f : initial holding point at runway threshold;

• R f : a set of alternate routes knowing the origin and the destination of f .

We have some assumptions in order to simplify the problem while keeping some level of reliability.

• Aircraft taxi with a constant speed for a given link. For each link we use the average speed value analyzed with the real data to take into account different taxiway types (e.g., taxiways near parking areas and near runways have significantly different speeds, as do straight taxiway segment and turning segment);

• Ramp area is beyond the scope of this work, instead we use the notion of meta-gate.

Decision variables

In order to optimize the ground movement, we now consider several potential control points as decisions. For each flight f ∈ F, the decision variables are defined as follows:

• r f ∈ R f : taxi-in or taxi-out route;

• t h f : holding time (waiting time at runway threshold for departures and time spent in runway crossing queues for arrivals);

• p f : pushback time;

• h f : holding point for arrival. CDG south-side runway layout shown in Fig. 4 motivates us to use arrival holding point as decision variable.

In reality, simultaneous flight crossings can enhance departure runway throughput.

Figure 4 CDG south side runway layout Furthermore, the following auxiliary variables are introduced:

• t E f f : final take-off time for departure or runway crossing time for arrival at E f . It is calculated based on the route chosen and the associated taxi speed;

• t f : completion time for flight f : t f = t E f
f for departures, t f is equal to in-block time for arrivals.

These decision variables are discretized considering that in practice, discretized time slots are assigned for the flights:

• t h f ∈ {0, ∆t, 2.∆t, ..., N h .∆t}, where ∆t is a time slot, N h is the maximum allowed number of holding time slots, T h = N h × ∆t is the maximum holding time. T h depends on the type of movement (arrival or departure). We define T a h and T d h as maximum holding time for arrival and for departure respectively.

• p f ∈ {P 0 f , P 0 f + ∆t, P 0 f + 2.∆t, ..., P 0 f + N p .∆t}, N p is the maximum allowed number of pushback delay time slots, T p = N p × ∆t is the maximum pushback delay.

Constraints

Airport operational constraints are taken into account:

• Minimum taxi separation of s = 60 meters [START_REF] Deau | Airport surface management and runways scheduling[END_REF] between two taxiing aircraft.

• Take-off runway wake turbulence separations shown in Table 2.

• Holding point capacity (the maximum number of lights waiting at holding point). For arrivals, it is usually one or two due to the fact that a landing flight can not hold too long time to vacate the position for the next landings. For departures, it's a parameter called runway pressure adjusted by controllers considering demand over the period. Based on the route network structure in Fig. 1, and in order to express the previous mentioned separation standards and capacity constraints, given an instantiation of decision variables, we define:

• C n -the total number of conflicts on nodes. For each node n, we record all the flights f passing this node with the time t n f , and sort according to t n f . For two successive aircraft f and g passing the node n, we must make sure that t n gt n f > t s , where t s is the minimum time separation calculated based on s and the taxi speed on node n. Otherwise we increase C n by 1.

• C l -the total number of conflicts on links. For each given link l, we record all the flights passing this link with the time at the link entry and exit. Then we sort into two lists and compare, if the orders of two aircraft are swapped, then a link conflict is detected and C l is increased by the rank difference between entry and exit. Besides, if two aircraft use the same link but come from opposite direction, the exit time of the previous aircraft must be earlier than the entry time of the latter one, otherwise C l is increased by 1.

The node-link conflict detection methodology is similar with previous work [START_REF] Ma | Merging flows in terminal maneuvering area using time decomposition approach[END_REF]. Moreover, we add the bi-directional link conflict detection in this work.

• C r -the total number of conflicts on runways.

For each departure runway r and for two successive take-off flights f and g, we have t r g -t r f ≥ s f g , where t r f and t r g are take-off times for flight f and g respectively, otherwise C r is increased by 1.

• C h -the total number of conflicts on holding points. For each holding point, we first make sure that the sequence of waiting flights remains the same. If not, C h is increased by the rank difference. Then by calculating the maximum number of aircraft simultaneously waiting in the queue, we compare it with the maximum holding capacity. If it exceeds the maximum holding capacity, we increase C h by the exceeded capacity.

• C = C n + C l + C r + C h -the

total number of conflicts

Then the previous separation and capacity constraint is transformed to C = 0 (conflict-free).

Objectives

Remark that one of our objectives in this paper is to investigate the impact of taxi reroutes on airport performance. One of the main roles of taxi reroutes is to avoid aircraft conflicts. Therefore, we decide to relax the conflict-free constraint and put C in our objective function.

The objective function that we want to minimize is:

C+α f ∈D (p f -P 0 f f ∈F h f +γ         f ∈D (t f -p f ) + f ∈A (t f -L f )         ,

SOLUTIONS APPROACHES

The benefits of integrated airport optimization, such as runway scheduling, taxiway routing and gate assignment are promising. However, the complexity of the integrated problem would grow, when in practice the computational time is critical. Heuristics and hybrid methods may have more potential than exact approaches for tackling this problem [START_REF] Atkin | The airport ground movement problem: Past and current research and future directions[END_REF]. Due to these high combinatorics, we propose a meta-heuristic algorithm -simulated annealing -to address the problem.

Simulated Annealing (SA) is a meta-heuristic that simulates the annealing of a metal, in which the metal is heated up and slowly cooled to move towards an optimal energy state. It can easily be adapted to largescale problems with continuous or discrete search spaces. In SA, the objective function to be minimized is analogous to the energy of the physical problem. A global parameter T is used to simulate the cooling process. A current solution may be replaced by a random "neighborhood" solution accepted with a probability e ∆E T , where ∆E is the difference between corresponding function values. We start cooling process from a high initial temperature T 0 (which can be determined by a heating process or defined by user), the current solution changes almost randomly at a higher temperature, thus the algorithm is able to trap out of local minima. The probability to accept a degrading solution become smaller and smaller when T decreases. Therefore, at the final stages of the annealing process, the system will converge to a near-global or global optimum. In order to adapt SA to the airport ground optimization problem, several parameters need to be defined: initial temperature and initial acceptance probabilities, cooling schedule, neighborhood function, equilibrium state and termination criterion. For our problem, some parameter values are listed in Table 3. Moreover, to generate a neighborhood solution, instead of simply choosing randomly a flight f , we proceed it by two steps: first resolve conflicts, then minimize time changes. In Algorithm 1 for conflicts resolution, for each aircraft, we use the number of conflicts as its performance indicator. For arrivals, the performance involves runway crossing conflicts and ground conflicts, denoted as crossingPerfo and groundPerfo respectively in Algorithm 1. For departures, we record their take-off conflicts (denoted as takeoffPerfo) and ground conflicts. The algorithm targets one flight involved in conflicts and changes its decision with regard to its sub-performances. For example, if a departure has ground conflict with another aircraft, it is clearly useless to change its holding time at the runway threshold in order to solve this type of conflict. Instead, the pushback time or the taxi route should be changed to generate new neighborhood solution. The fact that our neighborhood definition is based on the total number of conflicts, intensifies the neighborhood generation and accelerates conflicts resolution. Once a conflict-free solution is reached, we change our strategy, to target aircraft with time decision changes (pushback delay or holding time) and try to decrease this value by using the neighborhood function in Algorithm 2.

Algorithm 1 Neighborhood function for resolving conflicts

Require: For each flight, we record its takeoffPerfo, crossingPerfo and groundPerfo, the sum is denoted as totalPerfo; P c = crossingPerf/totalPerfo; P t = takeoffPerfo/totalPerfo; P g = groundPerfo/totalPerfo; 1: Choose one flight f involved in conflicts based on its performance; 2: Generate random number, ν = random(0,1);

3: if f ∈ A then 4:
if ν ≤ P c then choose with equal probability between holding point and holding time change; The SA terminates the execution either if the maximum number of transitions and the minimum tem-perature are achieved, or if an acceptable solution is obtained (for example, in the conflict-resolution case, SA stops when C = 0).

RESULTS

We test our methodology on a one-hour real data case at Paris CDG Airport. Numerical results with different settings of (user-defined) algorithm parameters are presented and discussed. The overall process is run on a 2.50 GHz core i7 CPU, under Linux operating system PC based on a Java code.

Real data analysis

We extracted one-hour dense traffic data from 9:00 a.m. to 10:00 a.m. on February 18, 2016. West configuration is activated over the course of the day. During this peak hour, we observed a long departure queue at runway 26R with radar tracks. A total of 31 arrivals and 69 departures were operated, including 65 Medium and 35 Heavy aircraft. Landing runway 26L and takeoff runway 26R in the south side are more charged with 22 arrivals and 38 departures respectively.

Three major aspects concerning airport ground performances are discussed in this section:

• Benefits of incorporating taxi reroutes on the airport performance metrics;

• Priority of arrivals and departures in runway crossing;

• Tradeoffs between controlled pushback and controlled holding time for departures.

Taxi reroute

In order to investigate the impact of taxi reroute on ground conflict resolution, we first set our objectives to be only C here. Remind that C includes runway conflicts, ground conflicts (link, node, bidirectional link) and holding conflicts. 30 random tests are launched and results are depicted in Table 4. In the case of "Without taxi reroute", we use standard route for each pair of runway and gate. At the end of algorithm running, we reached conflict-free solution for all the tests in taxi reroute case, while without taxi reroute, we have 2 times unsolved conflicts. Considering the average CPU time, taxi reroutes test is more than twice as fast compared to another case. Therefore, taxi reroutes can help reduce ground conflicts and reach a conflict-free solution faster.

Next, in order to test the taxi reroute influence on flight delays, we reset our objective function to with α = β = γ = 0.0001. Small values of weighting coefficients are chosen in order to ensure that conflict-resolution is the first priority. After reaching a conflict-free solution, the algorithm focuses on minimizing the delay. Even though we use relaxation on conflict constraints, by putting relatively small coefficients on other objectives, all the following results are conflict-free. Combining different objectives into one using a weighted sum approach is a common technique for multi-objective optimization problem. In future research, Pareto ranking-based algorithm [START_REF] Bandyopadhyay | A simulated annealing-based multiobjective optimization algorithm: Amosa[END_REF] can be a good option to target this problem.

C+α f ∈D (p f -P 0 f )+β f ∈F h f +γ         f ∈D (t f -p f ) + f ∈A (t f -L f )        
After running the algorithm, we observed that for landing runways 27R and 26L, the aircraft average holding times are similar in both situations. However, for take-off runways, a decrease of average holding time for runway 26R from 72 s to 63 s (12 %) and a decrease of average pushback delay for runway 27L from 62 s to 43 s (30 %) is reached after applying the taxi reroute strategy. The result is reasonable because in the case of Without Taxi Reroute, once a ground conflict is detected, the algorithm may modify the pushback time to solve conflicts, therefore more modifications are made compared to taxi reroute case.

Runway holding

To build a First-Come-First-Serve (FCFS) sequence, we allow a maximum pushback delay T p = 1 min with taxi reroute. In this way we obtain a ground conflictfree solution with initial time at runway threshold for both departures and arrivals. Thereafter, an arrival is obliged to cross the departure runway immediately as it reaches the holding point. Departures have to wait and take off in a FCFS order. In consequence, we build our FCFS sequence as shown at the top of Fig. 5.

Figure 6 Inter-arrival and inter-departure spacings for south-side runways To make a comparison, we still set T p = 1 min with taxi reroutes. Moreover, we choose T d h = 10 5 illustrates the holding time comparison of runway 26R in the south side. We can observe that much less time deviations are made in the optimized case compared to FCFS strategy. This is mainly thanks to the proper adjustment of arrival holding times. When two arrivals cross almost simultaneously the departure runway after a heavy take-off, it is beneficial to reduce runway waiting time and to increase runway throughput. Controllers give the same order facing high traffic demand in real cases. Blue lines represent the extra spacing between two consecutive aircraft. Even in the end of the second row a large extra spacing appears. This is due to the limit of aircraft earliest arrival time at the runway threshold.

Arrivals do have priority over departures. However, small adjustment of arrival holding time would reduce significantly departure holding time in high traffic demand period. Table 5 show that average departure holding time is reduced from 7 min to 2.6 min at the cost of increasing 0.5 min average arrival holding time. Moreover, the inter-departure spacing is reduced as well. Fig. 6 illustrates spacing reduction of the optimized case compared to the FCFS strategy.

Trade-offs between pushback time and holding time

We made several tests to see the trade-offs between pushback time and holding time for departures as shown in Fig. 7a. In case T p = 1 min, T d h = 15 min, aircraft are scheduled to start pushback as soon as possible, causing a large waiting time at the runway threshold with a high runway pressure. When we consider another extreme case T p = 10 min, T d h = 0, that is aircraft take off smoothly without any departure queue, the delay is transferred to the pushback at gate with engine-off. After testing several combinations of T p and T d h , one can observe that the departure runway waiting time can not be reduced without the cost of pushback delay. The average taxi time shown in Fig. 7b is impacted proportionally by the holding time as well. However, with a proper choice of T p and T d h , a compromise can be found to balance the two criteria. Therefore, controllers can choose the maximum holding time and maximum pushback delay as parameters according to their preference and current traffic demand.

In conclusion, our study shows that taxi reroutes can help reach a conflict-free solution faster and reduce delay for departures. The slight adjustment of arrival holding time can significantly reduce departure runway congestions in peak hour. The departure runway waiting time can not be reduced without the cost of pushback delay. However, the maximum holding time and maximum pushback delay can be used as decision support parameters for controllers to manage departure runway pressure.

CONCLUSIONS

This paper addresses the problem of comparing different control strategies (controlled pushback time, taxi reroutes, controlled holding time) on the airport surface to investigate their impacts and benefits. We propose an optimization approach to solve in a unified manner the ground movement problem and runway scheduling problem. A preliminary study case is conducted in a model based on operations of Paris CDG airport under the most frequently used configuration. Alternate taxi routes are constructed based on surface surveillance records with respect to current procedural factors. The preliminary analysis estimates that taxi reroutes can reduce ground conflicts and reach a conflict-free solution faster. The slight adjustment of arrival holding time can significantly reduce departure runway congestions during peak hour. The departure runway waiting time cannot be reduced without the cost of pushback delay. Future research will include uncertainty analysis considering pushback time, taxiing time and ramp area. In addition, multi-objective optimization of ground operations can be considered in a more holistic manner. More high-traffic demand scenarios need to be created for evaluation (e.g., evaluate one day traffic by time decomposition approach).

Integration of airside and airport movements can be investigated as well.

Figure 1

 1 Figure 1 CDG airport model in the west configuration

  (a) Define intersections on radar track (b) Route option 1 from R1 to M1 (c) Route option 2 from R1 to M1 (d) Route option 3 from R1 to M1

Figure 2

 2 Figure 2 Route generation based on radar tracks and intersections

(a) 9 Figure 3

 93 Figure 3 Taxiing route set example

where•

  C : Total number of conflicts; • f ∈D (p f -P 0 f ): Total pushback delay; • f ∈F h f : Total holding time; • f ∈D (t fp f ): Total taxi time for departures; • f ∈A (t f -L f ): Total taxi time for arrivals. and α, β and γ are weighting coefficients corresponding to pushback delays, holding time and taxi time respectively.

5 : 7 : 9 : 1 : 5 :

 57915 else choose with equal probability among taxi-in route, holding point and holding time change; else if f ∈ D then 8: if ν ≤ P g then choose with equal probability between pushback time change and taxi-out route change; else choose with equal probability among holding time, pushback time and taxi-out route change 10: end if 11: end if Algorithm 2 Neighborhood function for minimizing time changes Choose one flight f with time decision changes; 2: if f ∈ A then choose a new holding time between 0 and current one; 3: else if f ∈ D then 4: if pushback time changed then choose a new pushback time between 0 and current one; else if holding time changed then choose a new holding time between 0 and current one;

Figure 5

 5 Figure 5 Take-off time comparison between FCFS strategy and optimized case for runway 26R

Figure 7

 7 Figure 7 Airport performances

Table 1

 1 Route options count

	Number of route options k	1	2-5 6-9
	Number of pairs displaying k options	342 159	9

Table 2

 2 Single-runway separation requirements, s f g , in seconds.

	Category		Leading Aircraft, f
			Heavy Medium	Light
		Heavy	90	60	60
	Trailing Aircraft, g	Medium	120	60	60
		Light	120	60	60

Table 3

 3 Empirically-set parameter values of SA

	Parameter	Value
	Geometrical temperature reduction coefficient	0.99
	Number of iterations at each temperature step	100
	Initial rate of accepting degrading solutions	0.2
	Final temperature	0.0001*T 0

Table 4

 4 CPU comparison for T p = 10 min, T d h = 10 min, T a h = 3 min

	Decision Choice	With	Without
		Taxi Reroute	Taxi Reroute
	Average CPU (in s)	11	26
	Min CPU (in s)	4	4
	Max CPU (in s)	25	112
	Failed number	0	2

Table 5

 5 Average holding time comparison for southside runways

	Runway	FCFS average	Optimized average
		holding time	holding time
	26R	7.1 min	2.6 min
	26L	0	0.5 min

ACKNOWLEDGMENTS

This work has been partially supported by Civil Aviation University of China and by China Scholarship Council (CSC). The authors would like to thank Serge Roux for his assistance with data and technical support. The authors also would like to thank Jean-Baptiste Gotteland and Mathieu Cousy for their helpful discussions.