Sameer Alam
email: sameeralam@ntu.edu.sg

Supatcha Chaimatanan
email: supatcha@gistda.or.th

Daniel Delahaye
email: lahaye@recherche.enac.fr

Eric Feron
email: feron@gatech.edu

A Meta-Heuristic Approach for Distributed Trajectory Planning for European Functional Airspace Blocks

Functional airspace blocks (FAB) concept is adopted by the European airspace to allows cooperation between airspace users to manage the air trac ow, while ensuring eciency, safety, and fairness without the constraints of geographical boundaries. This integration of airspaces allow for exibility in airspace management and aircraft trajectory planning. This paper proposes a distributed air-trac ow management model to address four-dimensional (4D) trajectory planning over the European FAB. The proposed method is based on a meta-heuristic approach which uses a hybrid algorithm of simulated annealing and hill-climbing local search to separate a given set of aircraft trajectories in space and time domain (we term this as ight interaction), by allocating an alternative ight plan (route and departure time) to each ight. An innovative data structure, termed as FAB-Flight interaction matrix, captures the ight interaction in-

I. Introduction

The main objective of an Air Trac Flow Management (ATFM) system is to ensure safety, reduce delays, and balance demand and capacity among dierent components of the air transportation system [1]. In European airspace, ATFM activities are carried out by Eurocontrol's Network Manager Operations Centre (NMOC, previously called CFMU) which acts as a coordinator between air navigation service providers (ANSPs) and key stakeholder, such as airlines, airports, and the military.

However, fragmentation of airspace along national boundaries creates structural ineciencies, sub-optimal ight routing and makes implementation of ATFM strategies quite challenging as highlighted by [START_REF] Button | Single European Sky and the functional airspace blocks: Will they improve economic eciency[END_REF] [START_REF] Button | Single European Sky and the functional airspace blocks: Will they improve economic eciency[END_REF] and Dubot et al. 2015 [START_REF] Dubot | Building a smooth and dynamic opening scheme from graph partitioning-Exploring dynamic airspace congurations[END_REF]. To address these challenges, a concept of Functional Airspace Blocks (FAB) is being developed by Eurocontrol. A FAB is dened in the Single European Sky legislative package as an airspace block based on operational requirements and established regardless of state boundaries, where the air navigation services and related functions is performance-driven and optimised through enhanced cooperation among ANSPs or an integrated provider" [START_REF]Single European Sky legislative package, Tech. Rep. Regulation (EC)[END_REF].

As illustrated in Figure 1, nine FAB initiatives have been taken: UK-Ireland FAB, Baltic FAB, BLUE MED FAB, Danube FAB, Denmark-Sweden FAB, FAB CE, FABEC, North European FAB and South West FAB, two of these have already been implemented, namely the UK-Ireland and Denmark-Sweden FABs and rest are in progress [START_REF]Eurocontrol Functional Airspace Block[END_REF].

Establishment of FABs will have consequences for NMOC operations. Angenendt, 2005 [6] has argued that one of the key challenges will be to implement the centralized ATFM strategies in the individual FAB and to develop the basis for the cooperation among the FABs. Delegation of ATFM to individual FABs can lead to demand-capacity imbalances as each FAB may try to optimize trac ow in its own airspace without considering other FAB requirements.

As per regulation (EC) No 551/2004 of the European Parliament on the organization and use of the airspace in the single European sky by European Commission, one of the specic requirements for the implementation of FABs is that there implementation should be supported by a safety case i.e. implementing of FAB must not result in increased instances of loss of separation. Subsequent evaluation of Functional Airspace Block initiatives and their contribution to performance improvement by Eurocontrol Performance Review Report 2008 [START_REF] Eurocontrol | Performance review report[END_REF], Adler et al., 2014 [8] and Luppo et al., 2016 [9] have concluded that without increased coordination and information exchange amongst Airlines, ANSPs, Military and Airports FAB may not be able to deliver its objectives. Arijana et al., 2014 [START_REF] Modi¢ | Performance Scheme Implementation in Functional Airspace Block Central Europe[END_REF] in their analysis of expected ATM changes in central Europe have concluded that higher degree of collaboration amongst the ANSPs and that the states must ensure that ATM can provide services in conjunction with their neighbors with greater deal of cooperation.

Button and Neiva, 2013 [START_REF] Button | Single European Sky and the functional airspace blocks: Will they improve economic eciency[END_REF] examined the economic eciency of functional airspace blocks and found that greater local coordination of ANSPs may play in the longer-term policies of establishing a full Single European Sky. Alvarez et al., 2013 [START_REF] Ichoua | A scenario-based approach for the air trac ow management problem with stochastic capacities[END_REF] suggested that one of the indicators of FAB success is that ATM performances at network level are achieved in a collaborative distributed decisions making context.

The above research motivated us to establish a distributed ATFM model which forms the basis of interaction and information exchange among FABs to implement ATFM strategies. The goal of distributed ATFM model will be to generate strategic 4D trajectories which minimize ight interaction, a concept which considers both space and time elements of overlapping trajectories.

Thus we attempt to address both, the safety and coordination, issues in FAB implementation.

The proposed method is based on a meta-heuristic approach which uses a hybrid algorithm of simulated annealing and hill-climbing local search to separate a given set of interacting aircraft trajectories by allocating an alternative ight plan (route and departure time) to each ight. The distributed decision making for ight interaction resolution stems from an innovative data structure, termed as FAB-Flight interaction matrix, which captures the ight interaction information between and within the FABs. Simulated Annealing belongs to a special class of search methods, known as evolutionary algorithms which were introduced by Holland in 1992 [START_REF] Holland | Adaptation in natural and articial systems: an introductory analysis with applications to biology, control, and articial intelligence[END_REF]. Since, air trac ow management problems are dicult to solve by other traditional methods because of non-convex, non-linear, or non-analytic characteristics, evolutionary search methods oers an eective approach of addressing open research problems in this domain [START_REF] Hansen | Genetic search methods in air trac control[END_REF]. Further, seminal work of Cheng et al., [START_REF] Cheng | Air trac control using genetic search techniques, Control Applications[END_REF], Du et al. [START_REF] Du | Identifying vital edges in Chinese air route network via memetic algorithm[END_REF] and Delahaye and en-route conict resolution has added empirical support for evolutionary algorithms as an optimization methodology that is both eective and ecient in analyzing complex air trac ow management problems.

The paper also extends a trajectory planning method, based on route and departure-time allocation technique, developed by co-authors in Chaimatanan, 2014 [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] by including a local search strategy for better exploitation of search knowledge and a neighborhood structure for eective identication of FABs and associated ights which causes high interaction. This enables interaction resolution of ights from FABs which causes higher interaction in the overall airspace. This leads to faster convergence of algorithm especially when number of ights is large.

The paper is organized as follows: Section II describes the concept of the proposed distributed ATFM model for Functional Airspace Blocks. Section III presents the proposed model and methodology in a mathematical framework. Then, a method to compute interaction between trajectories is presented in Section IV. Resolution algorithm to the problem is explained in Section V. Finally, numerical results are presented and discussed in Section VI.

II. Concept Description

Centralized ATFM though oers a fair distribution and demand-capacity balance at a global level but stakeholder can only provide inputs and are not part of decision making. It may also lead to large number of pairwise reversals, i.e., the resulting sequence of ight arrivals can be quite dierent as compared to the original published ight schedules [START_REF] Bertsimas | Fairness in air trac ow management[END_REF]. Because of this deviation from the original ight ordering, it may becomes dicult to implement such a solution locally.

In a distributed ATFM the decision-making responsibilities are shared between a number of airspace users (airlines, ANSPs, Military, Airports). Some examples of distributed ATFM are Ration by Schedule (RBS) [START_REF] Barnhart | Demand and Capacity Management in Air Transportation[END_REF] and Ground Delay Program (GDP) [START_REF] Glover | Stochastic optimization models for ground delay program planning with equityeciency tradeos[END_REF][START_REF] Mukherjee | Predicting Ground Delay Program at an airport based on meteorological conditions[END_REF]. However the current stateof-the art is limited mostly to strategic planning, and the users' participation in planning reduces as the planning interval becomes smaller. With the gradual implementation of FAB concept in the European airspace, any distributed ATFM systems must also take into consideration multiple FAB Figure 2 and 3 illustrates the dierence in the centralized and distributed ATFM approach in terms of trajectory planning and optimization. In centralized ATFM the optimization algorithm controls the vector of decision variables, X, which are used by the simulation process in order to compute the performance (quality), y, of such decisions. However, In this approach, evolutionary algorithms may not be adapted to address such problems, mainly when the simulation environment requires huge amount of memory space as is often the case in nowadays real-life complex systems. As a matter of fact, in the case of a evolutionary algorithm based approach, the simulation environment has to be duplicated for each individual of the population of solutions, which may require an excessive amount of memory. In order to avoid this drawback, one may think about having distributed approach where only one simulation environment which could be used each time a point in the population has to be evaluated as follows. In order to evaluate one population, one rst consider the rst individual. Then, the simulation environment is initiated and the simulation associated with the rst individual is run. The associated performance is then transferred to the optimization algorithm. After that, the second individual is evaluated, but the simulation environment must be rst cleared from the events of the rst simulation. The simulation is then run for the second individual, and so on until the last individual of the population is evaluated. In this case the memory space is not an issue anymore, but the evaluation time may be excessive and the overall process too slow, due to the fact that the simulation environment is reset at each evaluation.

In a FAB scenario, the trac ow management will be highly interdependent and will demand a signicant cooperation with other FABs. One way to achieve this is by having distributed 4D trajectory planning amongst FABs which can ensures conict-free trajectory for each aircraft. This will then translates into trajectory based operations where aircraft are required to y a negotiated conict-free trajectory through respective FABs.

In long term planning, airline operation center (AOC) generates ight schedules and ight plans and takes into consideration seasonal demand/trends, passenger projection, spoke-hub congurations, maintenance Schedule, pilots/Technical sta roster and available slots at airports. While in medium term planning Air Navigation Service Providers project trac loads based on airline schedules and ight plans for 2 to 12 Hours and takes in to consideration weather patterns. In Europe, this is done centrally by Eurocontrol Central Flow Management Unit (CFMU), now known as NMOC, for all the member states. CFMU also ensures that in the medium term all ight plans are conict and congestion free. In short term planning (Airspace Flow Planning) ight coordination is done on a regional scale with 30 Minutes to 2 Hour Timescale. Flow restrictions are used for overcrowding which includes Miles-in-Trail and Time-in-Trail. While in Tactical planning (Sector Trac Planning) Sector Air Trac Controller Executive (R-Side) is responsible for managing aircraft through the use of radar (5 minutes interval) and Sector Air Trac Controller Planning (D-Side) is responsible for coordinating hand-os with surrounding sectors (5 to 20 Minute Timescale).

The focus of this paper is on medium term planning where we envision that the demand and capacity balancing and initial ight planning will be taken care by NMOC and the responsibilities of ensuring that the ight plans are conict and congestion free will rests on individual FABs where the ight is originating (i.e. the controlling FAB). This is a loose form of distributed air trac ow management where the decision of de-coniction/de-congestion (collectively, we call it ight interaction) rests with FAB and not with a centralized body.

A. Proposed Model

The main idea proposed in the paper is a distributed decision making model which can enable eective information sharing among the FABs for 4D trajectory planning which are interaction free.

The goal of the proposed method is to separate a given set of aircraft trajectories in space and time domain by allocating an alternative ght plan (route and departure time) to each ght in a given FAB. Since the focus of this paper is on medium-term planning, we assume that the problem of demand-capacity balancing are addressed in the strategic long-term planning. In medium term planning the key issues is minimizing the interaction between trajectories. An interaction between trajectories occurs when two or more trajectories have an eect on each other; for instance, when trajectories occupy the same space at the same period of time. Therefore, contrary to the concept of conict, 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 minimum separation time between aircraft crossing at the same point.

The proposed concept is developed as follows: as illustrated in Figure 4, airspace users (airlines, airports, ANSPs etc.) submits relevant information (ight plan, slots, capacity) to NMOC, which then applies centralized trac ow management strategies to match demand with capacity and other airspace constraints and generate revised ight plans. These ight plans are then used as an input to a FAB-Flight interaction Matrix. It is a 2D matrix which captures the ight interaction information between and within FABs. One dimension of the matrix is termed Controlling FABs and the other dimension is termed Intermediate FABs. A Controlling FAB is dened as a FAB where a given ight is originated or activated (in case of an enroute ight entering European airspace), whereas an intermediate FAB is dened as a FAB through which a given ight traverse (over y), terminates (lands) or exits the European airspace.

B. FAB-Flight Interaction Matrix

The FAB-Flight Interaction Matrix is developed as follows: as illustrated in Figure 6, for N FABs in a given airspace A, a 2D matrix of N rows and N columns is developed. The row vector of the matrix represents the number of ight interactions caused by a Controlling FAB C j in the Intermediate FABs I i for i = 1 to N . The column vector of the matrix represents the number of ight interactions caused by the Controlling FABs C j for j = 1 to N in an Intermediate FAB I i Fig. 6 FAB-Flight Interaction Matrix

I i = IN T (C j , I A ) IN T (C j , I B ) ... IN T (C j , I i ) (1) C j =             IN T (C A , I i ) IN T (C B , I i ) .... IN T (C N , I i )             (2)
For example, as illustrated in Figure 6, ight interactions due to ights controlled by FAB Therefore, the total number of ight interactions U , in a given FAB i, can be given by summing the column vector:

U i = N j=1 IN T (C j , I i ) (3)
The total ight interaction V in a given airspace A (which comprises of N FABs) can be given by

V = N j=1 U j (4)
The average contribution of each FAB, in overall ight interactions, for a given airspace A can be given by:

U j rel = U j /V
(5) for j = 1 to N.

C. Trac Flow Management Strategy

The distributed Trac Flow Management strategy is developed as follows: rst the Intermediate FAB with highest number of ight interaction is identied as a candidate FAB (equation 6).

F AB I i = M AX(U A , U B ..., U N ) (6)
Then, for the FAB I i , the Controlling FAB C j which generated highest number of ight interaction is identied (equation 7).

F AB C j = M AX(IN T (C A , I i ), IN T (C B , I i ), ..., IN T (C N , I i )) (7)
The ATFM strategies (Space-Time separation) are then applied on randomly selected (tnessproportional selection) ights in Controlling FAB C j . The tness proportional method implies that the probability of a ight to be selected is simply proportionate to the number of interactions it causes. Flight interactions are recomputed given the revised ight plans and the FAB-Flight Interaction Matrix is updated. This process is repeated until the FAB-Flight Interaction Matrix is interaction free. Figure 7 illustrates the updated process, where the decision made by each FABs are evaluated by the optimization process. Then, the information of interaction based on FABs-Flight interaction matrix is feed back to each FABs, which then make new decision and repeat the process until a solution that leads to minimum overall interaction is reached.

III. Mathematical Modelling

This section sets the mathematical framework of the distributed air trac ow management methodology we are proposing. First, a denition of interaction between trajectories is given. Then, the route / departure-time allocation techniques adapted for the distributed ATFM model is presented.

A. Interaction between trajectories

Interaction between trajectories indicates, when two or more trajectories occupy the same space at the same period of time. It is dierent from the conict situation, which corresponds simply to a violation of the minimum separation (i.e. 5 NM horizontally and 1,000 ft vertically). In this paper, the concept of interaction takes into consideration the duration of conict. However, additional separation conditions, such as time separation, topology of trajectory intersection, distance between trajectories, etc. can also be taken into account.

Consider a given set of N discretized 4D trajectories, where each trajectory i is a time sequence of 4D coordinates, P i,k (x i,k , y i,k , z i,k , t i,k ), specifying that aircraft must arrive at a given point (x i,k , y i,k , z i,k ) at time t i,k , for k = 1, . . . , K i , and K i is the number of sampling points of trajectory i.

Consider a point k of trajectory i, interactions at point P i,k , denoted Φ i,k , may be dened as the total number of times that the protection volume around point P i,k is violated. Figure 8 illustrates interaction in the horizontal plane between N = 3 trajectories measured at point P i,k . The interaction associated with trajectory i, denoted Φ i , is therefore dened to be:

Φ i = Ki k=1 Φ i,k . (8)
Finally, the total interaction between trajectories, Φ tot , for a whole trac situation is simply dened as:

Φ tot = N i=1 Φ i = N i=1 Ki k=1 Φ i,k . (9)
One can observe that the measurement of the interaction between trajectories implicitly take into account the duration of conict between trajectories. A practical methodology to compute the value of the interaction between trajectories in a large-scale context is presented in Section IV.

B. Route/Departure-time allocation

In order to separate the trajectories in 3D space and time domain, we formulate a route/departure-time allocation techniques based on [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF]. The objective is to nd alternative 4D trajectory for each ight, so as to minimize the total interactions between trajectories.

Given data. A problem instance is given by:

• A set of initial N discretized 4D trajectories with associated controlling FAB;

• The discretization time step, ∆t;

• The number of allowed virtual waypoints, M ;

• The maximum allowed advance departure time shift of each ight i, δ i a < 0;

• The departure time shift step size, δ s ;

• The maximum allowed delay departure time shift of each ight i, δ i d > 0;

• The maximum allowed route length extension coecient of each ight i, 0 ≤ d i ≤ 1;

• The length of the initial en-route segment of each ight i, L i,0 .

The alternative departure time and the alternative route to be allocated to each ight are modelled as follows.

Alternative departure time. The departure time of each ight can be shifted by a positive (delay) or a negative (advance) time shift. Let δ i ∈ ∆ i be a departure time shift attributed to ight i, where ∆ i is a set of acceptable time shifts for ight i. The departure time t i of ight i is therefore t i = t i,0 + δ i , where t i,0 is the initially-planned departure time of ight i. The departure time shift

∆ i :={-N i a .δ s , -(N i a -1).δ s , . . . , -δ s , 0, δ s , . . . , (N i d -1).δ s , N i d .δ s }. (10)
Alternative trajectory design. In this work, an alternative trajectory is constructed by placing a set of virtual waypoints, denoted

w i = {w m i |w m i = (w m ix , w m iy )} M m=1 , (11)
near the initial en-route segment and then by reconnecting the successive waypoints with straightline segments as illustrated in Fig. 9. To limit the route length extension, the alternative en-route prole of ight i must satisfy:

L i (w i ) ≤ (1 + d i ), (12)
where L i (w i ) is the length of the alternative en-route prole determined by w i . Fig. 9 illustrated initial and alternative trajectories, constructed with M = 2 waypoints, where the location of each waypoint is constrained to be in a rectangular-shape possible location. Let W m ix be a set of all possible normalized longitudinal locations of the m th virtual waypoint on trajectory i. For each trajectory i, the normalized longitudinal component, w m ix , is set to lie in the interval:

W m ix := m 1 + M -b i , m 1 + M + b i , (13)
where b i is a (user-dened) parameter that denes the range of possible normalized longitudinal component of the m th virtual waypoint on trajectory i. To obtain a regular trajectory, the normalized longitudinal component of two adjacent waypoints must not overlap, i.e.

m 1 + M + b i < m + 1 1 + M -b i (14)
and hence the user should choose b i so that

b i < 1 2(M + 1) . ( 15
)
Let W m iy be a set of all possible normalized lateral locations of the m th virtual waypoint on trajectory i. Similarly, the normalized lateral component, w m iy , is restricted to lie in the interval:

W m iy := [-a i , a i ], (16)
where 0 ≤ a i ≤ 1 is a (user-dened) model parameter that denes the range of possible normalize lateral location of the m th virtual waypoint on trajectory i, chosen a priori so as to satisfy [START_REF] Holland | Adaptation in natural and articial systems: an introductory analysis with applications to biology, control, and articial intelligence[END_REF].

Let us set the compact vector notation: δ := (δ 1 , δ 2 , . . . , δ N ), and w := (w 1 , w 2 , . . . , w N ),

We shall denote by u i the components of u. It is a vector whose components are related to the modication of the i th trajectory, thereby our decision variable is:

u := (δ, w).
Finally, the interaction minimization problem can be formulated as a mixed-integer optimization problem, as follows:

min u=(δ,w) Φ tot (u)
subject to

δ i ∈ ∆ i , w m ix ∈ W m ix , w m iy ∈ W m iy , for all i = 1, . . . , N, m = 1, . . . , M, (P1)
where W m ix , and W m iy are dened by [START_REF] Hansen | Genetic search methods in air trac control[END_REF], and ( 16) respectively.

IV. INTERACTION DETECTION

In order to evaluate the objective function, at a candidate solution, (u), one needs to compute interaction between the N aircraft trajectories. To avoid the N (N -1) 2 time-consuming pair-wise comparisons, which are prohibitive in our continental-scale application context, we use a grid-based interaction detection scheme [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF] which is implemented in a so-called hash table.

First, the airspace is discretized using a four-dimensional grid (3D space + time) for conict detection, as illustrated in Figure 10, as a time series of 3D grids which is sampled with discretization time step ∆t = t n -t n-1 . The size of each cell in the 3D grids is dened by the minimum separation requirements (N h and N v ).. The size of each cell in the 4D grid is dened by the minimum separation requirement and the discretization time step, ∆t (see below). Then, for each given 4D coordinate

P i,k (x i,k , y i,k , z i,k , t i,k
) of each trajectory i, we identify which cell, says C i,j,k,t , of the 4D grid contains

P i,k (x i,k , y i,k , z i,k , t i,k ).
Next, we consider each such cell C i,j,k,t and we successively check its surrounding cells (there are 3 3 = 27 such neighbouring cells, including cell C i,j,k,t itself). If one cell is occupied by an aircraft other than aircraft i itself, the horizontal distance (d h ) and the vertical distance (d v ) between the corresponding aircraft coordinates are measured. A violation of the protection volume is identied

when both d h < N h and d v < N v .
In order not to underestimate interaction, and to avoid using small value to ∆t which leads to large number of trajectory samples and long computational time, we propose an inner-loop Fig. 10 Four-dimension (3D space + time) grid for conict detection, illustrated as a time series of 3D grids.

algorithm, detecting interaction between two sampling times, t and t + ∆t, by interpolating aircraft positions with a suciently small step size, t interp . The interpolation is performed only if no interaction is detected at time t. Then, one checks each pair of these interpolated points. The algorithm stops when an interaction is identied or when every pair of the interpolated points has been checked.

V. RESOLUTION ALGORITHMS

The 4D trajectory planning methodology for a distributed ATFM model, presented in this paper, relies on the interaction minimization problem introduced in Section III whose objective function values are obtained by simulation through the interaction detection scheme developed in Section IV. To solve the problem, a hybrid metaheuristic approach adapted to handle an air-trac assignment problem at the continent scale is used. It relies on a classical simulated annealing (SA) algorithm and two dierent local-search (LS) modules. The LS allows the system to intensify the search around a potential candidate solution while the SA allows the system to escape from a local trap and thereby ensuring the exploration of the solution space.

The concept of Simulated annealing (SA) in combinatorial optimization was introduced by S. Kirkpatrick et al. in 1983 [22]. It is a metaheuristic stochastic method of optimization that is well known for its ability to escape from local minima by allowing occasional moves that deteriorate the value of the objective function, such deteriorating moves being less and less as the number of iterations grows.

The concept of SA is 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;

• Cool 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 minimumenergy state is reached provided that the initial temperature is suciently high and the cooling time is suciently 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.

Temperature in a physical system is analogues to Control parameters in the optimization problem, where the objective function of the problem is analogous to the energy state of the system. As, from an initial solution, Simulated Annealing algorithm can proceeds in several iterations and at each iteration, a random neighbor is generated. Moves that improve the cost function are always accepted. Thus, at a given temperature, the lower the increase of the objective function, the more signicant the probability of accepting the move and higher the temperature, the more signicant the probability of accepting a worst move.

The proposed hybrid algorithm combines the SA and the local search algorithm such that the LS is considered as an inner-loop of the SA, which will be performed when a pre-dened condition is satised. The structure of the proposed hybrid algorithm of SA and LS methods is illustrated in Figure 11.

Previously, hybrid optimization algorithm have been applied to solve 4D trajectory planning at continent scale in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF], however the air trac ow were managed based on a centralized decision making. In order to apply the proposed hybrid SA-LS algorithm to distributed ATFM model, the LS search strategy and the neighborhood structure are modied according to proposed ATFM strategy, as described in Subsection II C.

For our problem, the simulated annealing proceeds as follows. , where ∆Φ = Φ N -Φ C is the dierence of energy between current state C and new state N . When the maximum number of iterations, n T , at a given temperature is reached, the temperature is decreased according to the user-provided pre-dened schedule, and the process is repeated until the pre-dened nal temperature, T f inal , is reached. More detail on simulated annealing can be found, for instance, in [START_REF] Dreo | Metaheuristics for hard optimization[END_REF].

Local search modules. The local search modules we use are heuristic methods that accepts a new solution only if it yields a decrease of the objective function. The process repeats until no further improvement can be found or until the maximum number of iterations n T LOC is reached.

The two local-search modules correspond to the two following strategies:

• Intensication of the search on one Particular Trajectory (PT). Given a ight i, this state-exploitation step focuses on improving the current solution by applying a local change from the neighborhood structure only to ight i (only the decision variables (w i , δ i ) are aected).

• Intensication of the search on the Interacting Trajectories (IT). Given a ight i, this state-exploitation step applies a local change, from the neighborhood structure, to every ight that is both subjected to the same controlling FABs as ight i, and currently interacting with ight i.

Neighborhood structure The proposed hybrid algorithm relies on a neighborhood structure to determine the next move. First, the controlling FAB, which generates the highest proportion of interaction, is identied. Then, a ight from the identied FAB is chosen such that the interaction associated with trajectory i, Φ i ≥ τ • Φ avg , where τ is a user-dened parameter and

Φ avg = Φ tot /N
is the average value of interaction.

In order to generate a neighborhood solution for a given ight, i, from the current conguration (w i , δ i ) C , one has to determine whether to modify the location of waypoints or to modify the departure time in the next move. In general, searching for the solution in the time domain would be more preferable since it does not induce extra fuel consumption. However, empirical tests show that limiting the search to only that degree of freedom results in prohibitive computational time.

Therefore, we introduce a user-dened parameter P w to control the probability to modify the location of the waypoints w i and such that the probability to modify rather the departure time is 1 -P w .

For a given ight i, the neighborhood operator generates a new set of virtual waypoints or a new alternative departure time according to this probability P w .

Hybrid algorithm (SA and LS). Here is how the above-mentioned methods are combined.

The methods are carried out according to pre-dened probabilities, which are proportional to the control temperature, T . The probability to carry out simulated annealing step, P SA , is:

P SA (T ) = P SA,min + (P SA,max -P SA,min ) • T 0 -T T 0 , (17)
where P SA,max and P SA,min are the maximum and minimum probabilities to perform the SA (predened by the user). The probability of running the LS module, P Loc , is given by:

P Loc (T ) = P Loc,min + (P Loc,max -P Loc,min ) • T 0 -T T 0 , (18)
where P Loc,max and P Loc,min are the maximum and minimum probabilities to perform the local search (dened analogously). And, nally the probability of carrying out both SA and the local search (successively), P SL , is:

P SL (T ) = 1 -(P SA (T ) + P Loc (T )). (19)
A key factor in tuning this hybrid algorithm is to reach a good trade o between exploration (diversication) and exploitation (intensication) of the solution space.

VI. Numerical experiments

The proposed distributed ATFM model is implemented in Java. The overall methodology is tested with air trac data involving ights over the European FABs, consisting of nine FABs listed in Table 1. First, it is tested with a set of trac consisting of 4,000 ights over the European FABs on a UNIX platform with 1.7 GHz processor and 8 GB memory. The parameter values chosen to specify the optimization problem are given in Table 2. The parameter values that specify the resolution algorithm are given in Table 3. Then, it is tested with a full day of en-route air trac over the European FABs consisting of 26,122 ights on a UNIX platform with 2.4 GHz processor and 32 GB memory, using the same parameters as given in Table 2 and The parameters in Table 2 which were domain dependent, were chosen based on the feedback from domain experts. For example, the number of ight level changes is kept at maximum of two ight levels (2000 feet) to minimize major changes in originally led cruise ight level. Similarly, the maximum departure time shift is kept to two hours keeping in line with airline practices. The parameters in Tables 2 and 3 are adopted from the centralized AFTM methodology developed by co-authors in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF]. These parameters were determined by initial set of experiments and ne tuning.

These parameters are kept same in this paper for comparison purpose. We have used following matrices for comparison:

• Initial and nal interaction between trajectories,

• The number of iterations performed resolves the interactions,

• The computation time to resolve the interactions,

• How the solution converged (Color Map) Time to converge and resolution of interactions for continental scale trac were the two key metrics for the algorithm evaluation. Since this paper also introduces a new data structure-FAB-Flight Interaction matrix, which forms the basis for information exchange and identication of the ights that need interaction resolution, one of the key hypothesis is how the distributed algorithm will scale with increase in trac from country level (4,000 ights) to continent level (26,000 ights) trac sample.

The initial and nal total interaction between trajectories, the computation time, and the number of iterations performed to solve the distributed problems compared to the one solved based on centralized decision making methodology (proposed in [START_REF] Chaimatanan | Large Scale 4D Trajectory Planning[END_REF]) are reported in Table 4. The FAB-Flight interaction matrix is presented as color map showing the number of interaction caused by controlling FAB into the Intermediate FAB. To better visualize the evolution of interaction resolution, the FAB-Flight interaction matrix color maps are reported at four key stages of experimentation: at the start, at 30%, 70% and at the end of the process.

A. Small Flight Sample: 4000 ights Results indicate that, in terms of iterations the distributed model takes 8,306 iterations as compared to 5,035 iterations in centralized model to solve 48,272 interactions with 4000 ights.

However in terms of CPU time the dierence between the two approaches is not signicant.

Figure 12 and Figure 14 presents the color map showing how interactions are resolved at dierent stages of the optimization process. Initially, the FAB that has the highest level of initial interaction is FAB 5 (FAB Europe Central). The controlling FABs which generates the highest number of interactions in FAB 5 are FAB 1 (Baltic FAB),4 (Danube FAB), and 5 (FAB Europe Central). In Distributed algorithm demonstrate approx 20% improvement in performance over centralized model with large trac sample. In terms of CPU time, the dierence between the two approaches is signicant at 194 minutes (over 3 hours). This is a signicant improvement over centralized approach as in medium term planning Air Navigation Service Providers project trac loads based on airline schedules and ight plans for 2 to 12 Hours and takes in to consideration weather patterns.

Figure 16 and Figure 18 shows the color map of trajectory interaction and their subsequent resolution at the start of process, 30%, 70% and at the end of process. The start of the process shows that both the algorithms start with the same set of interactions. At the end of 30% of the process, it can be observed that the distributed model rst resolves the high interactions caused

VII. Conclusion

In this paper, we have presented a distributed air trac ow management model, aiming at minimize interactions between aircraft 4D trajectories at continent-scale. The objective was to develop a basis of information exchange and interaction among FABs for implementing distributed AFTM with a population based search method based on Simulated annealing and local search.

The overall methodology is implemented and tested on a continent-size air trac over the European FABs, and then compared with the one obtained using centralized ATFM model. Although the trajectories can be separated only by modifying the horizontal ight prole and the departure time of each ight, the resolution algorithm nds an interaction-free solution in both the cases.

In distributed model, the resolution algorithm tries to minimize the total interaction by modifying ight plans of ights associated with the FABs that generate high interaction in the overall air trac. While, in the centralized model case, the resolution algorithm modies ight plans of ights that involves high interaction, without taking into account the FAB-Flight interaction information.

The convergence rate of the distributed model is much better than centralized model and is viable for strategic planning as well as pre-tactical planning purpose.

In the case of small data set, the computation time to reach interaction-free solution is not significantly dierent in the distributed and centralized model. However, when the trac size increases, the distributed model converges to interaction-free solution 3 hours faster than the centralized model. This is a signicant improvement over centralized model as in medium term planning Air Navigation Service Providers project trac loads based on airline schedules and ight plans for 2 to 12 Hours and takes in to consideration weather patterns.

The proposed method of distributed air trac ow management using FAB-Flight interaction matrix demonstrate that a trade-o between exploration and exfoliation in a population based search method for interaction resolution yields better results that the centralized trac ow management approach. This trade-o is achieved using a heuristic algorithm consisting of a population based search method augmented with a local search and a neighbourhood search function. This approach dynamically assess trajectory interactions from origination and controlling FABs while resolving them. This ensures that the convergence time, specically when it comes to plan continental scale trac, is small. However, in practice, interfaces and dierent implementations, time delays and so forth need to be managed which present diculties. These diculties may not be compensated by these theoretical advantages.

Fig. 1

 1 Fig. 1 European airspace is divided into nine Functional Airspace Blocks viz. UK-Ireland FAB, Baltic FAB (Lithuania, Poland), BLUE MED FAB (Albania, Cyprus, Greece, Italy and Malta), Danube FAB (Bulgaria, Romania), Denmark-Sweden FAB, FAB CE (Austria, Bosnia & Herzegovina, Croatia, Czech Republic, Hungary, Slovak Republic, Slovenia), FABEC (Belgium, France, Germany, Luxembourg, the Netherlands and Switzerland), North European FAB (Estonia, Finland, Latvia, and Norway) and South West FAB (Portugal, Spain).

Fig. 2 Fig. 3

 23 Fig. 2 Trajectory planning and optimization for FABs in Centralized ATFM.

Fig. 4

 4 Fig. 4 Proposed Concept for information sharing among the FABs for interaction-free trajectory planning.

Fig. 5

 5 Fig. 5 An airspace divided into three FABs (FAB A, FAB B and FAB C) with four ight ( A, B, C and D) scenario.

C

  A in the Intermediate FAB I A is given by row FAB C A and column FAB I A and denoted by IN T (C A , I A ). Similarly, the number of ight interactions due to ights controlled by FAB C C in the same intermediate FAB i.e. FAB I c is given by row FAB C c and column FAB I c .

Fig. 7 FAB

 7 Fig. 7 FAB-Flight Interaction Matrix update process.

Fig. 8

 8 Fig. 8 Interactions, Φ i,k , at sampling point P i,k of trajectory i.

Fig. 9

 9 Fig. 9 Initial and alternative trajectories with rectangular-shape possible location of M = 2 virtual waypoints.

  δ i will be limited to lie in the interval ∆ i := [δ i a , δ i d ]. Common practice in airports conducted us to rely on a discretization of this time interval using time-shift step size δ s . This yields N i a := -δ i a δs possible advance slots and N i d := δ i d δs possible delay slots of ight i. Therefore, we dene the set, ∆ i , of all possible departure time shifts of ight i by

  Fig. 11 Structure of the proposed hybrid algorithm of simulated annealing and hill-climbing local search methods.

Fig. 12

 12 Fig. 12 Evolution of FAB-ight interaction matrix for 4,000 trajectories using distributed ATFM model.

Fig. 13

 13 Fig.[START_REF] Hansen | Genetic search methods in air trac control[END_REF] Evolution of interaction in each FAB for 4,000 trajectories using distributed ATFM model.

Fig. 15

 15 Fig. 15 Evolution of interaction in each FAB for 4,000 trajectories using centralized ATFM model.

Figure 13 ,

 13 Figure 13, and Figure 15 presents the resolution of interaction over number of iterations in each FAB during the optimization process for distributed and centralized models respectively. It can be seen that for with small trac sample (4,000 ights) the centralized model converges faster to

  by the controlling FAB (for example FAB 5) on other intermediate FABs. The color map also demonstrate the advantage of incorporating the Local Search (LS) strategy with neighbourhood search to facilitate the exploration of the problem search space initially (till 70% of process) and then exploitation of knowledge gained as demonstrated in quick resolution of most of the interaction during last stage of optimization process.

Figure 17 and

 17 Figure17and Figure19shows the resolution of interaction over number of iterations in each FAB

Fig. 17

 17 Fig. 17 Evolution of the interaction in each FAB for 26,122 trajectories using distributed ATFM model.

Fig. 19

 19 Fig. 19 Evolution of the interaction in each FAB for 26,122 trajectories using centralized ATFM model.

  3, except that this time the number of iteration N I is set to 2,700.

	Table 1 European FABs.
	No. FAB name
	1	Baltic FAB
	2	Blue Med
	3	FAB Central Europe
	4	Danube FAB
	5	FAB Europe Central
	6	NEFAB
	7	NUAC FAB
	8	SW Portugal-Spain FAB
	9	FAB UK Ireland

Table 2

 2 Chosen (user-dened) parameter values for the optimization problem.

	Parameter	Value
	Discretization time step, ∆t	20 seconds
	Discretization time step for possible departure-time shift, δs 20 seconds
	Maximum departure time shift, delta i a = δ i d := δ	120 minutes
	Maximum allowed route length extension coecient, di	0.20
	Maximum allowed ight level shifts, li,max := lmax	2
	Maximum number of virtual waypoints, M	3

Table 3

 3 Empirically-set (user-dened) parameter values of the resolution methodology.

	Parameter	Value
	Minimum probability to perform SA step, PSA,min	0.8
	Maximum probability to perform SA step, PSA,max	0.9
	Minimum probability to perform local search step, PLoc,min	0.4
	Maximum probability to perform local search step, PLoc,max	0.6
	Number of iterations at each temperature step, NI	400
	Number of iterations of the inner-loop local search step, nT LOC 5
	Geometrical temperature reduction coecient, β	0.99
	Final temperature, Tf	(1/500).T0
	Inner-loop interpolation sampling time step, tinterp	5 seconds
	Probability to modify horizontal ight prole, Pw	1/3
	Probability to modify ight level, Pl	1/3
	Threshold value, Φτ	0.5 Φavg

Table 4

 4 Initial and nal interaction between 4,000 trajectories.

	case	N	ATFM	initial nal	solved	no. of	cpu time
			strategy	Φ D tot	Φ D tot interactions iterations (mins)
			Distributed		0	100%	8,306	2.47
	1	4,000		48,272
			Centralized		0	100%	5,035	2.08
			Distributed		0	100 %	509,924	369.67
	2	26,122		266,318
			Centralized		0	100%	632,002	563.53

Acknowledgement

Authors would like to thank Leilla Zerrouki from Eurocontrol for useful discussions on FAB design that greatly improved the manuscript.

VIII. References