Ahmed Khassiba
email: ahmed.khassiba@capgemini.com

Daniel Delahaye
email: delahaye@recherche.enac.fr

Simulated-Annealing Hyper-Heuristic for Demand-Capacity Balancing in Air Traffic Flow Management

Keywords: demand-capacity balancing, hyper-heuristics, simulated annealing

Very critical convective weather leads to sharp drop in air-traffic en-route sector capacity keen on creating severe demand-capacity imbalances, called also hotspots. Following the knock-on effect, these imbalances spread across the network, causing the so-called netspots. The problem of demand-capacity hotspot/netspot mitigation with demand-side measures aims at solving these imbalances using ground delays and reroutings, while minimizing the total delay at arrival. A simulated-annealing hyper-heuristic relying on efficient neighborhood operators is proposed. A study on the combination of different delay-based and rerouting-based neighborhood operators is conducted. The best combination is used to solve a large-scale and challenging instance, in a relatively short computation time.

I. INTRODUCTION

Up to 2019 just before Covid-19 crisis, "ATFM en-route" has been the second most important cause of delay (after airlines operations) in air transport in Europe [START_REF] Walker | All-causes delay and cancellations to air transport in Europe -annual report for 2020[END_REF], with the highest delays being recorded in summer 2018. The upturn in air traffic, more recently, paves the way unfortunately for delay figures to soar again. Delays due to "ATFM en-route" are typically those impacting flights take-off time, when enroute sectors cannot accommodate the high air traffic demand. A sector capacity can fall during some time intervals due to local critical convective weather conditions (e.g., storms), which consequently, creates an imbalance between demand and capacity. An imbalanced sector during a given period of time, e.g., a sector-hour, is called a demand-capacity hotspot.

In the current operational process, on the day of operations, Flow Manager Positions (FMPs) use Computer-Assisted Slot Allocation (CASA) system to mitigate such hotspots, by defining regulations for the imbalanced sectors. CASA issues new (delayed) take-off times for flights that are impacted by regulations. Following the knock-on effect, regulated flights propagate delay through the network, by creating new demandcapacity hotspots, hundreds of nautical miles far from the convective weather zone, as observed in summer 2018. The aggregation of demand-capacity hotspots spread across the network, e.g., the European airspace, created by the knock-on effect, while their origin is a local convective weather, is called a netspot. The use of conventional CASA regulations reveals inefficient to solve netspots as excessively large delays are generated, especially for flights that are captured by several regulations. Several improvements of the current operating concepts in air traffic flow management have been identified to "hedge" against the impact of very-critical weather conditions on the air traffic situation at the network level. In this context, SESAR-H2020-funded project ISOBAR (for Artificial Intelligence Solutions to Meteo-Based DCB Imbalances for Network Operations Planning) seeks to leverage Artificial-Intelligence techniques to better predict high-risk convective weather and sector capacity drop, to identify and to efficiently mitigate hotspots at the network level, i.e., netspots, both at pre-tactical (one day before operations) and tactical (the day of operations) stages.

This paper focuses on hotspot mitigation at the pre-tactical stage, defined as the problem of accommodating air traffic through the airspace, by delaying take-off times and changing planned routes, in order to meet reduced sectors capacity, while minimizing deviations with respect to the original flight plans. A hyper-heuristic algorithm using an efficient combination of neighborhood operators is proposed as a solution method. In the following Subsection I-A, the general principles of hyper-heuristics are presented. The scope of the paper is presented in Subsection I-B.

A. Hyper-heuristics

Search algorithms called metaheuristics are approximate solution methods that iteratively search for a good-quality solution to a given optimization problem, while trying to escape local optima. Metaheuristics operate directly on the solution space, and keep track of either: a single solution (e.g., simulated annealing, tabu search, etc), or a population of solutions (e.g., genetic algorithm, etc). On the other hand, hyper-heuristics are approximate optimization algorithms with a higher level of automation than metaheuristics. They operate on a space of low-level heuristics (each of which operates directly on the solution space), and use a strategy for either selecting or generating low-level heuristics, throughout the iterations. For an extensive literature review on hyper-heuristics, the reader is redirected to [START_REF] Burke | A classification of hyper-heuristic approaches: revisited[END_REF], [START_REF] Burke | Hyper-heuristics: A survey of the state of the art[END_REF]. A selection-based hyperheuristic relies on: (1) a pool of low-level heuristics, (2) a strategy to select the low-level heuristic for the next iteration, and (3) a strategy to update low-level heuristics scores. Lowlevel heuristics (LLH), called also neighborhood operators in the metaheurisitcs vocabulary, are procedures that operate on the solution space in order to generate a new solution from a current one. Low-level heuristics are problem-specific, and can be classified either as diversification or as intensification operators. Very often, an initial score is given to each lowlevel heuristic, then is updated throughout the iterations based on each LLH performance. There are several strategies to select the LLH to apply for the next iteration, among them: the max strategy, and the roulette-wheel strategy. The max strategy selects always the LLH with the best score, while the roulette-wheel randomly selects LLHs with probabilities that are proportional to their scores. Hyper-heuristics have been applied successfully to many optimization problems; to name a few recent applications in air traffic management: [START_REF] Juntama | Hyperheuristic approach based on reinforcement learning for air traffic complexity mitigation[END_REF], [START_REF] Dalmau | Comparison of various temporal air traffic flow management models in critical scenarios[END_REF].

B. Scope of the paper

In this paper, we consider an optimization problem similar to the one presented in [START_REF] Bolić | Reducing ATFM delays through strategic flight planning[END_REF], where demand-capacity hotpsots are to be mitigated using two types of demand-side measures: ground delay, and alternative trajectories selection. Alternative trajectories are assumed to be agreed on between airspace users and the Network Manager at a strategic level. We focus on situations of very-critical convective weather, leading to sharp en-route sector capacity drop with respect to nominal values. We opt for entry count as the capacity metric, since finer-grained metrics, such as occupancy count [START_REF] Dalichampt | Hourly entry counts versus occupancy count-relationship, definitions and indicators[END_REF] or geometrical complexity [START_REF] Delahaye | Air traffic complexity map based on non linear dynamical systems[END_REF], cannot be used with a lookahead time of 24 hours before operations, due to high uncertainty. We propose a simulated-annealing-based hyper-heuristic using efficient low-level heuristics relying on the concept of flight overload factor. Recently, a similar concept has been applied successfully to airspace congestion mitigation using departure slots assignment [START_REF] Lavandier | Selective simulated annealing for large scale airspace congestion mitigation[END_REF]. We conduct a numerical study to identify the best combination of neighborhood operators, the best strategy to select them, and the best strategy to update their score, at every iteration.

The remainder of this paper is organized as follows. Section II presents more formally the demand-capacity-balancing problem, with demand-side mitigation measures. Mathematical notations are summarized in Table I. The solution method based on a simulated-annealing hyper-heuristic is described in Section III. The numerical study is reported in Section IV. Section V concludes this paper.

II. PROBLEM STATEMENT

Consider a time horizon decomposed into P equal-length, possibly overlapping, time periods. Each time period p ∈ P = {1, 2, . . . , P } has a start time t p and an end time t p > t p , such that: t p -t p = l > 0, where l is the time-period duration or

T = T f f ∈ F
Vector of the selected trajectories decision variables,

t 0 = t f,0 f ∈ F
Vector of the take-off times decision variables, w 1

Weight for the total arrival delay term w 2

Weight for the total unsustainable overload term w ∈ (0, 1]

User-defined coefficient in roulette-wheel selection α ∈ (0, 1)

Scale factor of the cooling schedule β Minimum acceptance ratio threshold during cooling length; the same for all time periods, e.g., l = 60 minutes.

We have also: t p+1 -t p = m > 0 and m ≤ l, where m is the time shift, or move, between two consecutive time periods such that p ∈ {1, 2, . . . , P -1}. Typically, m = 20 minutes. Consider a set of S sectors indexed by s ∈ S. Let C p s ≥ 0 denote the capacity of sector s during time period p ∈ P, expressed as the maximum number of flight entries to sector s during time interval [t p , t p ), that can be handled by air traffic controllers (ATCOs). During a given time period p ∈ P, a sector s ∈ S can be either open or closed. An open sector is controlled by a couple of ATCOs, hence it has a positive capacity C p s > 0. When sector s is grouped with other sectors into a larger collapsed sector, it is no more monitored on its own, while the larger sector is. In this case, sector s is said to be closed during time period p, and its capacity is set to

C p s = ∞.
Consider a set of flights indexed by f ∈ F. Each flight f has an expected take-off time (commonly abbreviated ETOT) t org f,0 , and an original trajectory T org f . A trajectory is represented as a sequence of entry times to a subset of sectors. Assume a flight f crossing a subset of sectors S f ⊂ S following a trajectory T f , then:

T f = {(t f,s , s), s ∈ S f } (1)
where t f,s is the entry time of flight f to sector s. Let F p s denote the set of flight indices that enter sector s during period p:

F p s = {f ∈ F | s ∈ S f and t f,s ∈ [t p , t p )} (2)
Let D p s denote the demand of sector s during time period p ∈ P, expressed as the number of flight entries to sector s during the time interval [t p , t p ). We establish that:

D p s = |F p s | if sector s is open during period p 0 if sector s is closed during period p (3)
When demand exceeds capacity in a given sector s during a time period p, i.e., D p s > C p s , we say that sector s is overloaded during period p. The overload of open sector s during period p is expressed as:

O p s =    D p s -C p s C p s if sector s is open during period p 0 if sector s is closed during period p (4)
In practice, small overloads can still be managed by ATCOs. Hence, we define v ∈ [0, 1) an overload-tolerance threshold, set by the user. A typical value for overload tolerance is v = 10%. A demand-capacity hotspot occurs in a sector s during period p when:

D p s > C p s × (1 + v), i.e., O p s > v.
In order to mitigate hotspots, that is to ensure that demand and capacity are balanced, flight departure times can be delayed and original trajectories can be changed.

We consider F ⊂ F a subset of flights that can be modified, i.e., delayed or rerouted, while flights in F\ F are considered as fixed. For example, in Europe, flights eligible for modification are those departing from European airports, i.e., typically their airport code starts with L, E, GC, or GM. However, flights planned to depart from America or Asia and to use the European airspace cannot be modified by the Network Manager Operating Center. We assume that every flight in F has a set of alternative trajectories T f ∋ T org f . The Demand-Capacity Balancing (DCB) problem in AT-FCM, defined for a set of counting periods P and set of sector indices S, consists in finding for every flight f ∈ F, an adjusted take-off time t f,0 ≥ t org f,0 , and a trajectory T f ∈ T f , such that all demand-capacity hotspots are mitigated, i.e. O p s ≤ v for all (s, p) ∈ S ×P, while minimizing the total delay at arrival, incurred by all flights. Moreover, the take-off delay per flight may be bounded by d max f for every flight f ∈ F. The complete model of the demand-Capacity Balancing (DCB) problem in ATFCM reads:

min t0,T f ∈ F (t f, * -t org f, * ) (5)
s.t. O p s ≤ v (s, p) ∈ S × P (6) t f,0 -t org f,0 ∈ 0 , d max f f ∈ F (7) T f ∈ T f f ∈ F ( 8
)
where t org f, * , and t f, * denote respectively the original, and the adjusted, arrival time of flight f . These arrival times can be computed directly from take-off times, t f,0 and t org f,0 , using the flight duration of the corresponding trajectories.

Relaxed DCB problem

Finding a feasible solution to the DCB problem is not straightforward, especially when the initial overload situation is too critical. For this reason, we consider the relaxed problem where capacity constraints [START_REF] Bolić | Reducing ATFM delays through strategic flight planning[END_REF] are integrated into the objective function, as follows:

min t0,T w 1 f ∈ F (t f, * -t org f, * ) + w 2 (s,p)∈S×P max + (O p s -v)
s.t. Constraints ( 7) and ( 8)

(Relaxed DCB)
where w 1 and w 2 are user-defined weights, and max + (•) = max (0, •). Note that for the relaxed problem, the solution with no delay at departure and no change in trajectories, is a feasible solution. We call the term

(s,p)∈S×P max + (O p s -v) the total unsustainable overload.
A feasible solution of the relaxed problem (Relaxed DCB) that has a null total unsustainable overload is a feasible solution of the original problem, given by equations ( 5) - [START_REF] Delahaye | Air traffic complexity map based on non linear dynamical systems[END_REF]. We define the best solution as the one that has, most importantly, the lowest total unsustainable overload (ideally no unsustainable overload). Then, for two solutions with equal total unsustainable overload, the one with the smallest total arrival delay is the better.

III. SIMULATED-ANNEALING HYPER-HEURISTIC

Simulated annealing (SA) is a well-known metaheuristic that mimics the annealing process from metallurgy [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. While classical SA relies on a single neighborhood operator, we propose a simulated-annealing-based hyper-heuristic to solve the relaxed version of the DCB problem (Relaxed DCB), using a pool of neighborhood operators, designed using the idea of flight overload factor. Each neighborhood operator is given a performance score, that is updated along the iterations according to predefined rules. At every iteration, a neighborhood operator is selected according to specific strategies. In our study, we evaluate different rules to update the performance scores of neighborhood operator, and different strategies to select the neighborhood operator at every iteration. In Subsection III-A, we briefly review SA applied to our problem. The hyper-heuristic procedure is presented in Subsection III-B. The pool of neighborhood operators is presented in Subsection III-C. The studied rules to update the operators' score and the strategies to select them are presented in Subsection III-D.

A. Simulated annealing

SA starts with an initial solution, represented by a vector of decisions. In our problem, each decision specifies, for each flight, the on-ground delay and the trajectory, i.e.

t f,0 -t org f,0 , T f . Using a predefined neighborhood operator, the current solution is evolved into a new one, and evaluated in terms of objective function. The acceptance of the new solution is driven by the Metropolis acceptance criterion [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], that stipulates that an improving solution is always accepted, while a non-improving solution is accepted with a probability controlled by an internal parameter, called the temperature. Assuming obj new and obj curr are the objective function values of the new and the current solutions respectively, and Temp the temperature parameter, the acceptance probability, AP , according to the Metropolis criterion, in the minimization case, reads:

AP = 1 if obj new < obj curr , e -(obj new -obj curr )/Temp otherwise.
In SA, the temperature parameter is initialized, either statically or dynamically through a heat-up procedure, then it decreases following a cooling schedule. The Metropolis criterion ensures that as the temperature decreases, the probability of acceptance of non-improving solutions converges to zero, which allows for exploration at the early search stages, and exploitation at the late ones. Note that a too fast cooling schedule could lead to pre-mature convergence of the search to a local optimum, while a too slow cooling may extend the computing time. Reheating or temperature reset are also common strategies to help the search escape from local optima. Several strategies to set the initial temperature, as well as to conduct cooling, and reheating or temperature reset, are proposed in [START_REF] Franzin | Revisiting simulated annealing: A component-based analysis[END_REF].

B. Hyper-heuristic procedure

Algorithm 1 Hyper-heuristic procedure end for 13:

Temperature cooling or reheating 14: end while

We propose a hyper-heuristic algorithm based on simulated annealing, as sketched in Algorithm 1. In the initialization phase (lines 1-3), we assign and initialize a performance score to each neighborhood operator, we compute the initial temperature using a heat-up procedure (see Subsection III-B1), and we set the current solution to the "no delay, no rerouting" solution. At the beginning of each round of iterations, a neighborhood operator is selected according to a predefined strategy (line 6), then used to generate a new solution (line 7). The Metropolis criterion is applied to decide whether to accept the new solution or not (lines 8-10). Subsequently, the performance score of the selected neighborhood operator is updated according to predefined rules (line 11). After every round of iterations, the temperature parameter is revised, either cooled or re-heated (line 13), as described in Subsection III-B2. The algorithm stops either when the temperature is too low, or when a predefined time limit is exceeded (line 4).

1) Heat-up procedure: We set the temperature to a first value, for example, computed as the maximum difference between two neighbor solutions [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. We generate a random solution. Then, it is evolved during nb iter iterations using neighborhood operators, where nb iter is a user-defined parameter. The heat-up procedure is stopped if the acceptance ratio over the last nb iter iterations is greater than or equal to a predefined target value (typically 99%). Otherwise, the temperature is increased by a heat-up factor (e.g., by a factor of 2), and a new round of neighboring solution generation is restarted for another nb iter iterations.

2) Cooling and re-heating schedule: We opt for a geometric progression with a scale factor α ∈ (0, 1) for temperature cooling; Temp i+1 = Temp i × α, where Temp i and Temp i+1 denote the temperature values for two successive iterations. We define β the minimum acceptance ratio threshold over nb iter iterations, typically β = 90%. During the cooling process, if the acceptance ratio over the last nb iter iterations is below β, then we suspend the cooling schedule. The temperature is kept at the same value for additional nb iter add iterations. The cooling schedule is resumed when either the acceptance ratio is greater than or equal to β, or a predefined maximum number of iterations at the same temperature, denoted max nb iter, is reached. This procedure has the merit to avoid too early convergence without causing divergence, unlike re-heating schemes where the temperature is simply increased.

C. Neighborhood operators based on flight overload factor

In order to design an efficient hyper-heuristic, one needs to design well-performing neighborhood operators tailored to the problem at hand. For the DCB problem with demand-side decisions, one key idea is to modify flights that are very much involved in hotspots. To measure the degree of involvement of individual flights in the overall demand-capacity overload situation across the region of interest, we define the flight overload factor o f ≥ 0, for a given flight f ∈ F, as the sum of all fractions of unsustainable overload of the sector-periods crossed by f . It reads:

o f = s∈S f p∈P |t f,s ∈[t p ,tp) max + (O p s -v) (9)
The definition of flight overload factor not only takes into account the number of unsustainably overloaded sector-periods crossed by a given flight but also the amplitude by which the overload tolerance is exceeded in each sector-period. Note that a flight f with a null overload factor, i.e., o f = 0, is a one that does not cross any hotspot. We shall call a flight with a positive overload factor, i.e., o f > 0, an overloaded flight. The neighborhood operators should focus on flights with high overload factor, by either delaying them, to overstep the overloaded time periods, or by rerouting them, to hopefully deviate from some overloaded sectors. To allow a flexible selection among highly overloaded flights, we rely on roulettewheel selection based on flights overload factor, as follows. First, the largest flight overload factor, o max , is found. Then, all flights with an overload factor within [w × o max , o max ] are shortlisted, where w ∈ (0, 1] is a user-defined coefficient. Finally, a roulette-wheel selection is applied on the shortlisted flights where each flight's selection probability is proportional to its own overload factor, as described in Algorithm 2. The proposed neighborhood operators are summarized in Table II. The operators from DO-0 through DO-2 apply delay, while operators from RO-1 and RO-2 apply rerouting. The neighborhood operator LS is a local search procedure, that can be applied directly after any DO-x operator, so that the delay assigned to the flight selected in DO-x can be locally optimized. When applied in a given iteration, DO-x operators assign a new delay to a selected flight, regardless its last delay value, from the previous iterations. Hence, DO-x can be considered as diversification operators, while LS can be seen as an intensification operator. DO-2 and RO-2 are further described in Algorithms 3 and 4 respectively.

D. Strategies for neighborhood operators management

We study two strategies to update neighborhood operators scores:

• Additive: if the new solution improves the objective function's value, the score is incremented by 1. Otherwise, the score is decremented by 1, while ensuring that the minimum score is 1. • Customized: if the new solution improves the objective function's value, the score is incremented by 1. Otherwise, the score is square-rooted. We study two strategies to select the next neighborhood operator:

• Max strategy: the neighborhood operator with the highest score is systematically selected. • Roulette-wheel strategy: neighborhood operators are selected with probabilities proportional to their scores.

IV. NUMERICAL STUDY

We conduct a numerical study to identify the best combination of neighborhood operators, the best rule to update neighborhood operators' score, and the best strategy to select the next operator. First, we test the three delay-based neighborhood operators, with and without the local search procedure, introduced in Subsection III-C. The best delay-based operator is the one that solves (ideally) all unsustainably-overloaded sector-periods, while generating the smallest total delay. After selecting the best delay-based neighborhood operator, combinations with rerouting operators are tested, assuming that at every iteration of the hyper-heuristic, the two neighborhood operators are selected with a 50% probability each. The best combination between delay-based and rerouting-based operators is selected to test the rules to update neighborhood operators' score, and the best strategy to select the next operators, among those introduced in Subsection III-D.

Details about the studied instance and the parameter values is given Subsection IV-A, and Tables III and IV therein. Results of delay-based neighborhood operators are presented in Subsection IV-B. Results of the combination of the best delay-based neighborhood operator and rerouting operators are given in Subsection IV-C. Rules to update neighborhood operators' score and strategies to select neighborhood operators are studied in Subsection IV-D. Finally, more insights are given in Subsection IV-E about the best solution obtained in terms of the overload situation, the issued delay and rerouting. The proposed hyper-heursitic was developed in JAVA programming language, and all results are obtained on a Windows 11 platform, with an Intel(R) Core(TM) i5 CPU @ 1.60GHz and 12 GB RAM. The maximum time limit is set to 5 minutes.

A. Instances and parameter values

We focus on July 27th 2019, where very high delay was recorded due to very convective weather in Europe. All elementary and collapsed sectors in ACCs: Karlruhe (EDU-UUAC), Zagreb (LDZOACC), Barcelona (LECBACC), Marseille (LFMMACC), and Vienna (LOVVACC) are considered. Eurocontrol DDR2 traffic files, so6 and t5, are used. Entry times to sectors are retrieved from t5 files, while meta data on flights (callsigns, aircraft operators, etc) are taken from so6 files. Predicted capacity data was generated by machine learning, and provided as an input from ISOBAR consortium partners. Details on the instance size and the initial overload situation are provided in Table III. The hyper-heuristic parameter values are given in Table IV.

Alternative trajectories: To retrieve a set of alternative trajectories for each flight, we analyse the so6 and t5 files for the day of interest, and we group all trajectories by four criteria: origin (O), destination (D), aircraft type (AC) and airspace user (AU). Then, for every set of trajectories sharing the same values for the four criteria, we eliminate equivalent trajectories. Two trajectories are assumed to be equivalent if they cross the same sectors in the same order. The final sets of non-equivalent trajectories grouped by (O, D, AC, AU) are used to fill in the set of alternative trajectories for every flight. It is noteworthy that the merit of the retrieved alternative trajectories is twofold: they are realistic, and they integrate the airspace user's dimension.

B. Results of delay-based neighborhood operators

Results of the three delay-based neighborhood operators, with and without the local search procedure, are given in Table V and Figure 1. Figure 1 gives the evolution of the best solution's total unsustainable overload (upper figure) and total delay (lower figure) over iterations, for the 6 delaybased neighborhood operators. Note that total unsustainable overload of the best solution decreases over the iterations, for all neighborhood operators. DO-2 and DO-2+LS exhibit the fastest decrease, while DO-1 and DO-1+LS are stuck, very early, at a large value of the total unsustainable overload. The best solution with random delay operators, DO-0 and DO-0+LS, has a slow decrease in total unsustainable overload, but still ends at a low value. In terms of total delay, the best solutions using DO-1 and DO-1+LS present the lowest delay. DO-0 and DO-0+LS have a bell-shaped total delay evolution. For DO-0, after a dramatic increase, the total delay decreases significantly to reach competitive values with DO-2 and DO-2+LS. However, DO-0+LS reveals to be much slower (due to the additional local search at every iteration) than DO-0, hence, it has not enough time to decrease its total delay figure. Delaying systematically the most overloaded flight (DO-1) at every iteration reveals to be inefficient, since the total unsustainable overload is only slightly decreased. Clearly, highly overloaded flights are the best candidates to mitigate hotspots. The first expected effect of the additional local search procedure is to extend the time to accomplish one iteration, and hence to reduce the number of iterations made within the time limit of 5 minutes. Nevertheless, a better performance of DO-1+LS and DO-2+LS was achieved compared to their simple counterparts, in less iterations. All in all, the best delaybased neighborhood operator is DO-2+LS, since it mitigates all unsustainable overloads, while generating the smallest total delay. In the next subsection, we present the results of the combination of DO-2+LS with the two proposed rerouting operators RO-1 and RO-2.

C. Results of the combination of delay and rerouting operators

Table VI shows the total residual unsustainable overload, as well as information on the total delay, the number of delayed flights, and the number of rerouted flights (column "# rerouted flights"), for the two combinations of DO-2+LS with RO-1 and RO-2. Since both tests reached the time limit of 5 minutes, the "CPU" column was omitted. We notice that using rerouting operators decreases significantly the total arrival delay, as expected. Moreover, using a roulette wheel over flights with a high overload factor performs better than automatically selecting the flight with the highest factor at every iteration. The best combination with rerouting operators is DO-2+LS+RO-2, where delay and rerouting operators are selected with a 50% fixed probability each, along the iterations. In the following subsection, we study the effect of the two proposed score update rules (additive, customised) and the two neighborhood selection strategies (max, roulette wheel).

D. Results for strategies to select neighborhood operators and to update their score

Table VII shows the results for the retained neighborhood operators DO-2+LS+RO-2 with the four possible combinations of score update rules and neighborhood operator selection strategies, in addition to the baseline case (denoted "No update, Roulette-Wheel") retained from Subsection IV-C. All tests reached the time limit of 5 minutes. Notice that the max rule to select neighborhood strategies lead to not calling the rerouting operator. On the other hand, updating the neighborhood scores, either with the additive or the square-root rule, improves on the case with no update. This reveals that at different search stages, delay and rerouting operators do not have the same performance, and the hyper-heuristic calls them with different probabilities. The percentage of cumulative calls to the two operators is plotted in Figures (2a) and (2b), for the additive and the custom score update rules respectively. The best combination as reported in Table VII is to update scores using the additive rule and to use the roulette-wheel for neighborhood operator selection. In the next subsection, we show more details on the solution returned by this best combination.

E. Best solution details

The best combination of neighborhood operators (DO2-LS-RO2) with an additive score-update scheme, and a roulettewheel selection over operators was applied again on the studied instance with a time limit of 1 minute. The same solution was returned, showing that the best solution is found in the early stages of the search and is not updated later. Also, this highlights the performance of the designed hyperheuristic to find very-good quality solutions in a relatively short computation time. Some details of the best solution are given in Figures 3, 4a, and4b. The distribution of overloads over sector-periods, before and after optimization, is given in Figure 3. All unsustainably overloaded sector-periods are mitigated, at the price of an increase in the overload ranges of -10 -0% and 0 -10%. The distribution of arrival delay, and of extra mileage are given in Figures 4a and4b respectively. Flights with a negative arrival delay, correspond to rerouted flights, with alternative trajectories that are shorter than their planned trajectories. Although, such reroutings may be interesting to solve the DCB problem, in real life, airlines may not be keen on accepting such reroutings. Indeed, there are situations where the shortest trajectory yields a high operating cost for the airline, due to high overflight fees of some countries. Moreover, 24.73% of the rerouted flights do not cross the sectors considered for the DCB problem anymore. In this sense, they do not contribute to the load of the studied sectors, however, their new trajectories must increase the load of sectors beyond the region of interest. This suggests that neighbouring sectors should also be monitored for demand-

V. CONCLUSION

The problem of demand-capacity hotspot mitigation using ground delay and rerouting was considered in the case of very critical convective weather, leading to a sharp decrease in en-route sector capacity. A simulated-annealing-based hyperheuristic using efficient low-level heuristics relying on the concept of flight overload factor was proposed. A numerical study was conducted to identify the best combination of lowlevel heuristics. Applying a roulette-wheel on a subset of flights with a high overload factor reveals to be the best flight selection rule along the algorithm iterations, either to apply delay or rerouting. Moreover, following the delay heuristic by a local search procedure decreased significantly the total arrival delay. Future work will focus on integrating uncertainty on capacity data, as predictions 24 hours before operations cannot have high accuracy. Also, a two-stage procedure can be put in place in order to update the mitigation plan, once a more accurate capacity drop prediction is provided on the day of operations.

Although a DCB solver based on the optimization of individual flights, called cherry-picking measures, can reduce the total delay compared to CASA, its applicability is still very questionable. On the one hand, FMPs think usually in terms of regulations applied to sector-periods or traffic volumes. Also, airspace users may not be keen on delaying their flights, just for the sake of "global optimality", while they do not cross hotspots. Consequently, such a paradigm change may need time, and hybrid DCB mitigation technologies as well as, more collaborative operations, may be required.

Algorithm 2

 2 Roulette-wheel selection based on flight overload factor Require: w ∈ (0, 1] a user-defined coefficient 1: Find the largest flight overload factor, omax 2: Select all flights with an overload factor within [w × omax , omax], and store them in F max . 3: Let: prob f = o f k∈ F max o k , for every f ∈ F max 4: Apply a roulette-wheel selection among the selected flights where prob f is the selection probability of flight f ∈ F max .

Require: d min and d max 1 : 2 :Algorithm 4 RO- 2 : 2 :

 12422 -overloaded flight LS Search for the best delay of a given flight, around a given delay value RO-1 Reroute the most overloaded flight RO-2 Reroute a highly-overloaded flight Algorithm 3 DO-2 : Delay a highly-overloaded flight Select a flight using the roulette-wheel procedure defined in Algorithm 2 Apply a random delay between d min and d max to the selected flight. Reroute a highly-overloaded flight 1: Select a flight using the roulette-wheel procedure defined in Algorithm 2 Randomly select an alternative trajectory for the selected flight.

Figure 1 .

 1 Figure 1. Evolution of the best solution's total unsustainable overload (upper figure) and total delay (lower figure) over iterations for the 6 single neighborhood operators based on delay. The same time limit of 5 minutes was set.

  (a) Additive rule. (b) Square-root-based rule.

Figure 2 .

 2 Figure 2. Percentage of cumulative number of calls to delay and rerouting operators, with a roulette-wheel selection strategy, and different score-update rule.

Figure 3 .

 3 Figure 3. Sector-period overload distribution: comparison between the initial and the final situations.

  (a) Arrival delay distribution. (b) Extra-mileage distribution among rerouted flights.

Figure 4 .

 4 Figure 4. Best solution characteristics.

TABLE I .

 I NOTATIONS SUMMARY

	Notation	Definition
	F	Set of flight indices
	F ⊂ F	Set of indices of flights eligible for hotspot
		mitigation measures
	F max	Set of indices of highly-overload flights
	P	Set of counting periods
	S	Set of sectors of interest
	S f ⊂ S T f	Set of sectors crossed by flight f ∈ F Set of alternative trajectories for flight f ∈	F
	P	Number of counting periods in P
	t p	Start time of counting period p
	tp	End time of counting period p
	m	Time shift between two consecutive time periods
	l C p s D p s O p s T f ∈ T f T org f t f,0 t org f,0	Length of one time period Capacity of sector s during period p Demand of sector s during period p Overload of sector s during period p Selected trajectory decision variable for flight f ∈ Original trajectory for flight f ∈ F F Take-off time decision variable for flight f ∈ Original take-off time for flight f ∈ F	F
	t f,s t f, * t org f, * d max f	Entry time of flight f ∈ F to sector s ∈ S F Arrival time for flight f ∈ Original arrival take-off time for flight f ∈ F Maximum possible take-off delay for flight f ∈	F
	o f	Flight overload factor for flight f ∈ F
	omax	Maximum flight overload factor over all flights in F
	v	Overload tolerance threshold

TABLE III .

 III INSTANCE DETAILS

	Day of traffic:
	Nb counting periods (l = 60 min, m = 20 min)	72
	Nb sectors of interest		175
	Nb flights crossing airspace of interest		12 224
	Nb flights eligible for hotspot mitigation measures	10 924
	(departing from airports: L* E* GC* GM*)
	Predicted overload situation:
	Maximum overload per sector-period		56%
	Total unsustainable overload		48.84
	Nb overloaded sector-periods		522
	Nb unsustainably-overloaded sector-periods	330
		TABLE IV. PARAMETER VALUES
	Parameter	Value	Parameter	Value
	α	0.99	w 1	1.0
	β	0.90	w 2	2000.0
	w	0.2	v	10%
	nb iter	200	d max f	180 min
	nb iter add	50
	max nb iter	400

TABLE V .

 V RESULTS FOR THE DELAY-BASED NEIGHBORHOOD OPERATORS. TILIM STANDS FOR TIME LIMIT = 5 MINUTES.

	Neigh. Op.	CPU Res. unsus. Total delay # delayed Avg. delay Max delay
		(sec)	overload	(min)	flights	(min)	(min)
	DO-0	215.7	1.03	142 346	9980	14.3	178
	DO-1	TiLim	37.56	10 775	104	103.6	177
	DO-2	TiLim	0	168 980	1850	91.3	180
	DO-0+LS	TiLim	4.83	736 058	9633	76.4	180
	DO-1+LS	TiLim	36.82	9 172	104	88.2	178
	DO-2+LS	TiLim	0	122 205	1403	87.1	180

TABLE VI .

 VI RESULTS FOR THE COMBINATION OF THE TWO REROUTING NEIGHBORHOOD OPERATORS WITH THE RETAINED DELAY-BASED OPERATOR.

	Neigh. Op.	Res. unsus. Total Arrival # delayed Avg. delay # rerouted
		overload	delay (min)	flights	(min)	flights
	DO-2+LS+RO-1	0	102 950	1185	85.4	75
	DO-2+LS+RO-2	0	88 761	1177	75.4	181

TABLE VII .

 VII RESULTS FOR THE STRATEGIES TO SELECT NEIGHBORHOOD OPERATORS AND TO UPDATE THEIR SCORES.

	Score-update Neigh. Op.	Res. unsus. Total Arrival # delayed Avg. delay # rerouted
	rule	selection	overload	delay (min)	flights	(min)	flights
	No update	Roulette-Wheel	0	88 761	1116	75.4	181
	Additive	Roulette-Wheel	0	73 893	985	68.5	279
		Max	0	115 543	1329	86.9	0
	Custom	Roulette-Wheel	0	85 096	1086	73.2	244
		Max	0	115 543	1329	86.9	0
	capacity balancing, and new DCB problems, with a larger
	geographical scope, may need to be solved.

ACKNOWLEDGMENT ISOBAR project has received funding from the SESAR Joint Undertaking (JU) under grant agreement No 891965. The JU receives support from the European Union's Horizon 2020 research and innovation programme and the SESAR JU members other than the Union.


Chunyao Ma

Sameer Alam

Qing Cai
email: qcai@ntu.edu.sg

Daniel Delahaye
email: daniel.delahaye@enac.fr

Air Traffic Flow Representation and Prediction using Transformer in Flow-centric Airspace

Keywords: air traffic management, flow-centric operation, air traffic flow prediction, transformer neural networks

The air traffic control paradigm is shifting from sector-based operations to cross-border flow-centric approaches to overcome sectors' geographical limits. Under the flow-centric paradigm, prediction of the traffic flow at major flow intersections, defined as flow coordination points in this paper, may assist controllers in coordinating intersecting traffic flows which is the main challenge for implementing flow-centric concepts. This paper proposes to predict the flow at coordination points through a transformer neural network model. Firstly, the flow coordination points, i.e., the major flow intersections, are identified by hierarchical clustering of flight trajectory intersections whose location and connectivity characterize daily traffic flow patterns as a graph. The number of coordination points is optimized through graph analysis of the daily flow pattern evolution. Secondly, air traffic flow features in the airspace during a period are described as a "paragraph" whose "sentences" consist of the time and callsign sequences of flights transiting through the identified coordination points. Finally, a transformer neural network model is adopted to learn the sequential flow features and predict the future number of flights passing the coordination points. The proposed method is applied to French airspace based on one-month ADS-B data (from Dec 1, 2019, to Dec 31, 2019), including 158,856 flights. Results show that the proposed prediction model can approximate the actual flow values with a coefficient of determination (R 2 ) between 0.909 to 0.99 and a mean absolute percentage error (MAPE) varying from 27.4% to 11.7% with respect to a 15-minute to 2-hour prediction window. The sustainability of the prediction accuracy under an increasing prediction window demonstrates the potential of the proposed model for longer-term flow prediction.

I. INTRODUCTION

The traditional Air Traffic Control (ATC) service is established in geographically partitioned sectors, organized and managed by the national Air Navigation Service Provider (ANSP) under expected traffic demand and airspace availability [START_REF] Gianazza | Learning air traffic controller workload from past sector operations[END_REF]. However, such local sector-based approaches prevent the regional level management and optimization of air traffic as it only manages individual flight paths within a sector. Facing the growth of air traffic demand and difficulties in subdividing heavily loaded sectors, researchers have started examining and testing the concept of sectorless ATC, which views the airspace as a whole instead of the current practice of dividing the airspace into small sectors [START_REF] Korn | Validating sectorless atm in the hungarian airspace: Results of human in the loop simulations[END_REF].

One primary practice of sectorless airspace is flow-centric operation [START_REF] Zeki | Business models for flight-centric air traffic control[END_REF], which relies on controlling and monitoring flowbased formation and evolution of air traffic instead of a fixed number of geographical sectors. The flow-centric approach is based on the management of dynamic flow corridors. It opens the opportunity to distribute air traffic more efficiently in the airspace without being constrained by sector boundaries.

Despite the benefits of the flow-centric concept, its implementation has been rather limited. One primary challenge is the coordination between traffic flows to avoid potential conflicts [START_REF] Gerdes | From free-route air traffic to an adapted dynamic main-flow system[END_REF]. The highly dynamic nature of air traffic flow brings uncertainties into the flight routes, traffic pattern, and airspace constraints, adversely affecting the effective, efficient, and safe coordination and regulation of air traffic. An accurate short-term prediction of air traffic flows, i.e., one to two hours look-ahead time [START_REF] Lin | Deep learning based short-term air traffic flow prediction considering temporal-spatial correlation[END_REF], in the airspace can provide network-wide flow information to ANSPs and airspace users. It can be a vital support to flow coordination beforehand to avoid conflict situations involving a large number of flights. For instance, flow regulation measures, such as traffic re-routing, can be executed in advance to provide a smooth and efficient flow of traffic when the predicted flow exceeds available capacity [START_REF] Sridhar | Aggregate flow model for air-traffic management[END_REF].

In the literature, most research modeled air traffic flows as the number of flights at specific geographical locations in the airspace, such as airways, waypoints, or sectors [START_REF] Zhang | A hierarchical heuristic approach for solving air traffic scheduling and routing problem with a novel air traffic model[END_REF]. one way of making air traffic flow prediction is to predict the individual flight trajectories over time and count the future number of flights in the airspace [START_REF] Fernández | Dart: a machine-learning approach to trajectory prediction and demandcapacity balancing[END_REF]. The prediction error of such trajectory-based methods increase exponentially when the forecast time horizon increases beyond 20 minutes, which may be caused by its sensitivity to the prediction accuracy of various flights under consideration [START_REF] Lymperopoulos | Model based aircraft trajectory prediction during takeoff[END_REF]. Furthermore, the dimension of the models depend on the number of flights in scope, making such models untenable in real world practice due to limited computational resources.

An alternation to the trajectory-based approach is the aggregate air traffic flow prediction model focusing on predicting the overall distribution of the air traffic flows [START_REF] Wei | Total unimodularity and decomposition method for large-scale air traffic cell transmission model[END_REF]. The aggregated prediction methods are less sensitive to the uncertainty factors related to individual flights, which can achieve a longer forecast time horizon with less prediction errors and requires significantly less computational cost comparing to the trajectory prediction based approaches.

One way to make aggregate flow prediction is to describe the number of future flights handled by a control center or a sector as a linear function of a collection of influencing observations. Examples include the Linear Dynamic System Models (LDSM) based on state transition matrices [START_REF] Sridhar | Modeling and optimization in traffic flow management[END_REF], [START_REF] Chen | A network based dynamic air traffic flow model for en route airspace system traffic flow optimization[END_REF] and the linear regression model [START_REF] Gilbo | A new model to improve aggregate air traffic demand predictions[END_REF]. The major shortcoming of those linear models is their incapability in describing complex relations other than linearity among the air traffic flow and the contributing variables. Considering the air traffic flow in the airspace is the consolidated operation of a variety of stakeholders over space and time, e.g., controllers, airspace users and ANSPs, simple linear models may not be able to capture the complex and dynamic behaviors of air traffic flow to produce accurate short-term predictions. Another widely discussed approach for aggregated flow prediction is the machine learning based model [START_REF] Polson | Deep learning for short-term traffic flow prediction[END_REF], such as time-series extraction and prediction of the number of flights passing a given waypoint using Extreme Learning Machine based models [START_REF] Zhang | Research on air traffic flow forecast based on elm non-iterative algorithm[END_REF], prediction of air traffic demand between two airports using long Short-Term Memory (LSTM) and Support Vector Machine (SVM) [START_REF] Gui | Machine learning aided air traffic flow analysis based on aviation big data[END_REF], and end-to-end spatial-temporal time-series prediction based on the concatenation of Convolutional Neural Networks or Graph Convolutional Networks (GCN) with recurrent neural networks, such as LSTM, to respectively capture spatial and temporal features [START_REF] Lin | Deep learning based short-term air traffic flow prediction considering temporal-spatial correlation[END_REF], [START_REF] Ma | Sector entry flow prediction based on graph convolutional networks[END_REF]. Most of the time-series models use the number of flights passing areas or points in the airspace as the only input feature. Those models mainly analyze the temporal or spatialtemporal variation of the number of flights [START_REF] Dalmau-Codina | A machine learning approach to predict the evolution of air traffic flow management delay[END_REF], which may fail to consider more detailed air traffic features, such as the flight sequences transiting through different flow coordination points, and therefore, limits the prediction accuracy.

In view of the above analysis, this paper proposes a transformer base model for predicting traffic flow at the major flow intersections.First, the flow coordination points are identified by clustering the flight trajectory intersections using a hierarchical clustering algorithm. The number of clusters is determined through graph analysis of the flow pattern evolution. Then, the complex and dynamic airspace flow features are described as a "paragraph" with "sentences" representing sequences of the time and callsigns of flights transiting through the coordination points. Finally, These sequential flow details are used as the input to a transformer neural network model to extract the contextual relations of elements in long sequences and predict future flow on the coordination points.

II. METHODOLOGY

A. Methodology Overview

Focusing on the flow-centric concept, this paper investigates air traffic flow prediction in two main steps. The first step is to develop a mechanism to identify the flow coordination points where major air traffic flows intersect from air traffic data. The second step is developing a flow prediction model to predict the number of flights passing the identified coordination points during different look-ahead periods. Fig. 1 presents a concept diagram of the proposed flow prediction method.

In the flow-coordination-points identification step, the daily flow patterns are represented as a graph whose nodes are extracted from clusters of daily flight trajectory intersection points and whose edges are constructed based on the connectivity between the nodes. To determine the optimal number of clusters for flow representation, this paper proposes a graph analysis approach that evaluates the homogeneity of nodes and edges across graphs constructed from different days of air traffic. By maximizing the homogeneity between the graph representations of the daily air traffic flow, the cluster number of trajectory intersection points is determined, and the centers of the clusters are identified as the flow coordination points.

In the air traffic flow prediction step, each flight trajectory is represented as a sequence of the identified coordination points through trajectory registration which minimizes the dissimilarity between the original trajectory and the trajectory formed from the coordination points. Then, the flow feature in the airspace during a period is modeled as sequences of flights transiting through the set of coordination points. It is forwarded to a transformer based flow prediction model to learn the contextual description of the complex dynamics of the air traffic flow and predict the behavioral dynamics of the flow components, i.e., the number of flights transiting through the coordination points during a future period.

B. Flow Regulation Points Identification

ADS-B data for aircraft surveillance can give highly accurate aircraft position and velocity information, providing a valuable source for analytical solutions to effective and efficient airspace usage [START_REF] Ma | Airway network management using braess's paradox[END_REF]. This paper identifies the flow coordination points based on the flight trajectories analysis using ADS-B data, including intersection points clustering, daily flow representation, and graph analysis.

1) Intersection Points clustering: As mentioned in Section II-A, the flow coordination points in this paper are defined as the positions in the airspace where air traffic flows intersect. Therefore, this paper proposes to determine the coordination points clustering the intersection points of flight trajectories.

The intersection points between flight trajectories in this paper are calculated as the intersections of the mapped trajectories onto the earth's surface determined by latitudes and longitudes. Then clustering algorithms are applied to extract the natural groupings of the trajectory intersection points to identify the flow coordination points.

Hierarchical clustering relies on the hierarchical decomposition of the data based on group similarities to find a multilevel hierarchy of clusters, where clusters at one level are joined as clusters at the next level [START_REF] Cai | Hierarchical clustering of bipartite networks based on multiobjective optimization[END_REF]. Air transportation networks are commonly a nodal hierarchy that follows the spoke-hub structure in which traffic flows are often concentrated around the traffic hubs, and flight routes are organized as a series of "spokes" that connect outlying areas to a hub area. There is a hierarchy of air traffic flows ranging from regional feeders to international hubs [START_REF] Seabra | Determinants of brazilian international flights: The role of hub-and-spoke and infrastructure variables[END_REF]. Therefore, this paper adopts the hierarchical clustering algorithm to discover the representative flow coordination points. Fig. 2 shows an example of the clusters of intersection points. The white curves depict the flight trajectories, and the solid red circles denote the major flow intersections, which are the centers of the identified trajectory intersection clusters. The size of each circle is proportional to the number of intersection points in the cluster.

Prediction Model

Identified flow Coordination points

Intersection points clustering

Graph analysis

Flow context extraction

The number of clusters to be identified by the hierarchical clustering algorithm is an important parameter, as the clustering outcome should be able to discover and represent the natural organization of traffic flow intersections in the airspace.

On the one hand, if the number of clusters is too large, the clustering result will be susceptible to small fluctuations in air traffic flow, which may be unable to capture major flow patterns to produce credible flow prediction results. On the other hand, when the number of clusters is too small, the identified flow intersections can largely deviate from real flow paths, which dilutes the implementation significance as a supporting tool to controllers for flight coordination.

2) Daily Flow Pattern Representation: To determine the optimal number of clusters that should be identified by the clustering algorithm, this paper proposes a graph [START_REF] Ma | Critical links detection in spatial-temporal airway networks using complex network theories[END_REF] based approach which represent the daily air traffic flow pattern as a weighted graph G = (V, E), where V is the set of nodes denoting the air traffic flow components, i.e., the flow coordination points, which are represented as the centers of the identified clusters. The flow connections between the nodes can be described by the weighted edges E. If two nodes are passed by a flight consecutively, they will be considered connected, and there will be an edge between them. The weight of each edge is proportional to the air traffic volume transiting through it, i.e., the number of flights whose trajectories consecutively cross the two nodes connected by the edge.

C. Graph Analysis

Considering the necessity of consistent and dependable performance of airspace users for improvement of the ATM system predictability, [START_REF] Murc ¸a | Flight trajectory data analytics for characterization of air traffic flows: A comparative analysis of terminal area operations between new york, hong kong and sao paulo[END_REF], the representation of the traffic flow patterns should be able to describe the behavioral consistency of air traffic flow in the airspace as well as the daily alternations in the geographical positions of the intersection points and the variations in the flow structure. Thus, this section proposes to determine the optimal number of clusters by modeling the consistency of daily flow intersection patterns versus the changes in the number of clusters.

The flow patterns can evolve temporally following the daily alternations in the air traffic flow, geographically and structurally. Therefore, the flow pattern consistency is evaluated from two perspectives: 1) geographical consistency in node positions; 2) structural consistency in node connections.

a) Nearest neighbour analysis: To measure the geographical consistency in the daily node positions, a nearest neighbour based analysis is conducted to match the nodes on the temporal horizon. Concretely, given the traffic pattern graphs for two consecutive days, for each node in a graph, this analysis first searches in the other graph to find its nearest node based on the geographical coordinates. Considering the daily flow consistency, this pair of nodes is expected to represent the same flow coordination point in the airspace despite some minor positional shifts. Therefore, under a proper graphical flow pattern representation, the resulting nodes pair, a node with its identified nearest neighboring node, should be identical for the nodes in both graphs. Conversely, an unqualified graph representation may result in significantly different matches between the nodes. Therefore, the geographical consistency in node positions is quantified as the ratio of the mutually matched nearest node pairs to the total number of nodes in the graphs.

Given a number of nodes n in the graph G k , let as (φ i k+1 , λ i k+1 ), the nearest neighbour of node v i k in V k+1 is identified according the great circle distance:

V k = {v 1 k , v 2 k , ..., v i k , ..., v n k } represent
a i = arg min {j} j=1,...,n D(ij) = R • 2sin -1 (min(1, √ h)) (1)
where

h = sin 2 ∆φ 2 + cos φ i k cos φ j k+1 sin 2 ∆λ 2 ∆φ = φ i k -φ j k+1 ∆λ = λ i k -λ j k+1 (2)
R is earth's radius (mean radius = 6,371km), and h is the haversine formula which determines the great-circle distance between two points on a sphere given the longitudes and latitudes.

By identifying the nearest neighbouring node v ai k+1 of node v i k for i = 1, ..., n and identifying the nearest neighbouring node v bj k of node v j k+1 through the same calculation, two sets of matched node pairs between the two graphs can be obtained:

S 1 : (v 1 k , v a1 k+1 ), (v 2 k , v a2 k+1 ), ..., (v n k , v an k+1 ) and S 2 : (v b1 k , v 1 k+1 ), (v b2 k , v 2 k+1 ), ..., (v bn k , v n k+1
). Thus, the geographical consistency gc in node positions is quantified as:

gc 1 = |S 1 ∪ S 2 | n (3)
where |S 1 ∪ S 2 | represents the number of node pairs in the union of S 1 and S 2 , which will be denoted as l in the paper. b) Graph structure analysis: Let S = S 1 ∪ S 2 represent the set of paired nodes of graph G k and graph G k+1 , while

C k = {c 1 k , c 2 k , ..., c l k } and C k+1 = {c 1 k+1 , c 2
k+1 , ..., c l k+1 } denote the nodes from graph G k and G k+1 in S respectively. Let e i-j c-k represent the edge connecting nodes c i k and c j k . Let w i-j c-k represent the weight on edge e i-j c-k . The structural consistency in the air traffic flow pattern is measure by the divergence in the weighted edges between graphs. More specifically, with the set of paired nodes between graphs, the structural consistency is evaluated by the mutual flow connections of nodes in C k and C k+1 versus all the flow connections in graph G k and G k+1 , which is formulated as:

gc 2 = l-1 i=1 l j=i+1 min(w i-j c-k , w i-j c-k+1 ) 0.5 × ( n-1 i=1 n j=i+1 w i-j k + n-1 i=1 n j=i+1 w i-j k+1 ) (4
) By calculating the geographical consistency gc 1 and the structural consistency gc 2 versus varying number of clusters, a "saddle point" on the curve can be adopted as the optimal number of clusters for traffic flow representation.

D. Air Traffic Flow Prediction

Most traffic flow prediction model adopts recurrent neural networks, e.g., LSTM and GRU [START_REF] Islam | Foreign exchange currency rate prediction using a gru-lstm hybrid network[END_REF], to capture the temporal relations in a time series data. However, data in such sequential models need to be processed in order, i.e., the model architecture requires the output of the previous state to calculate the next state, which makes it difficult to capture long-term relations in the sequences, reducing the model's prediction power when the prediction window increases.

Transformers are the state-of-the-art neural networks that learn context by tracking relationships in sequential data like the words in the sentence. A transformer combines the attention mechanism (a mechanism that enables to highlight of relevant features of the input sequence dynamically) based representations and a parallel architecture in processing elements in a sequential input. Without the use of recurrent connections, the transformer model allows the ingestion of an entire sequence in parallel and has the potential to understand the relationship between sequential elements that are far from each other.

The seminal transformer uses an encoder-decoder architecture. The encoder extracts features from an input sequence, and the decoder uses the features as input to produce an output sequence, such as translating an English sentence into a French sentence. Therefore, the air traffic flow prediction model in this paper proposes to adopt a transformer neural network architecture to extract the air traffic flow sequence features and predict future air traffic flow.

Moreover, most existing air traffic flow prediction models are flow-density-based models that only analyze the temporal or spatial-temporal variations of the number of flights. They may not be a potent representation of the airspace's complex air traffic flow features. In this consideration, air traffic flow information in this paper is described by a detailed "paragraph" whose "sentences" are the sequences of flights passing the flow coordination points in a time sequential order.

In view of the above analysis, this section first registers the flight trajectories onto a sequence representation of the flow coordination points. The traffic flow context in the airspace during a period is then represented as a "paragraph" whose "sentences" are the sequences of flights (characterized by the flight callsigns) in the airspace. Then this contextual description of air traffic flow is used as input to a prediction model to learn from the flow context and predict the future air traffic flow. Details of each step will be illustrated in the following paragraphs.

1) Trajectory Registration: The objective of trajectory registration in this paper is to search for a sequence of flow coordination points that can optimally approximate the original trajectory. This objective is formulated as finding the minimum dissimilarity between the original trajectory and a representative trajectory constituted by a subset of the coordination points. Let T o represent a flight trajectory which is depicted by a sequence of points {p

1 o → p 2 o , → ..., → p M o } obtained from ADS-B data. Let F RP = {p 1 r , p 2 r , ..., p k r , ..., p N r } rep- resent the set of N identified flow coordination points. Let T r = {p r1 r → p r2 r →, ..., → p rm
r } represent a sequence of m flow coordination points with r 1 being a numerical index referring to the r 1 -th points in the original coordination points set. Therefore, the objective of trajectory registration is to find a representative trajectory T * r , which has the largest similarity to the original trajectory T o , from all the possible T r . The trajectory similarity in this paper is measured by the F re´chet distance [START_REF] Lyu | A partial-fréchet-distance-based framework for bus route identification[END_REF], a measure of similarity between curves by taking into account the location and ordering of the points along the curves. Thus, the objective function of the trajectory registration is formulated as follows:

obj = arg min Tr δ F (T o , T r ) T r = {p r1 r → p r2 r →, ..., → p rm r } {p r1 r , p r2 r , ..., p rm r } ∈ F RP (5)
The function δ F (T o , T r ) denotes the F re´chet distance between T o and T r . After trajectory registration, the trajectory T o will be represented as the sequence T * r = {p 2) Flow Context Extraction: As mentioned before, the traffic flow context in the airspace during a period is then represented as a "paragraph" whose "sentences" are the sequences of flights (characterized by the flight callsign) in the airspace. Use t 0 and t 1 to denote the start and end times of a period. Let F P k r denote the traffic flow at coordination point p k r during t 0 and t 1 . F P k r can be described as a sequence of flights (callsigns) according the time flights passing p k r , i.e., F P k r :

f 1 k , f 2 k , ..., f i k , ..., f m k k
with f i k denoting the i-th flight in the sequence and m k denoting the total number of flights passing p k r . Analogizing F P k r as a "sentence" depicting the flow context at p k r , the combinations of "sentences" during t 0 and t 1 , i.e., "t 0 -t 1 :

f 1 1 f 2 1 ... f m1 k ; f 1 2 f 2 2 ... f m2 2 ; ...; f 1 N f 2 N ... f m N N .
", constitute a "paragraph" description of the flow context in the entire airspace. The paragraphs for various periods will be used as inputs to the transformer based prediction model, introduced in the next section, to learn the contextual relations between air traffic flow and make predictions about the future air traffic flow.

3) Transformer based Prediction Model: The structure of the transformer based prediction model is shown in Fig. 3. The main components are tokenization, embedding, transformer encoder, and a fully-connected layer to produce the prediction outcome.

a) Tokenization: Given an input sequence, i.e., the traffic flow context including the time and the flight callsign sequences at different flow coordination points, tokenization [START_REF] Sun | Spectral-spatial feature tokenization transformer for hyperspectral image classification[END_REF] is applied to convert the input into a list of integers that can be embedded into a vector space. The flight callsigns in the input sequence refer to specific flights in the airspace, which differ from the natural word composition, whose prefixes, suffixes, and infixes can change their inherent meaning. Therefore the tokenization in this study adopts word-tokenization, which splits the data based on natural breaks and meaning, such as time (number), callsigns (words), and sequence separations (delimiters).

b) Embedding: After tokenization, the token embedding layer converts the list of integers, i.e., tokenized flow sequence, into a list of vectors, as is the case in language processing applications in general. For the model to use the order of the elements in the sequence, positional embeddings, which contain information on the relative or absolute position of the elements, are added to the token embeddings as the input to the transformer encoders. The embeddings are trained jointly with the rest of neural network. Back-propagation is carried through all the network layers up to the embeddings that are updated as other parameters. After embedding the elements in the input sequence, each of them flows through the transformer encoders which will be described in the follow section.

c) Transformer encoder: The transformer encoder processes each element in the input sequence and compiles the information it captures into a context tensor, i.e., an array of numbers. The encoder sends the context tensor to the fullyconnected layer to project the context to the predicted flows. As shown in Fig. 3, the transformer encoder is composed of a stack of N e encoder blocks. Each encoder block consists of two sub-layers. The first layer is a multi-head self-attention mechanism, and the second layer is a position-wise fully connected feed-forward network.

The self-attention mechanism allows the encoder to look at other elements in the input sequence when encoding a specific element. The self-attention creates three vectors from each of the encoder's inputs, Query vector, Key vector, and Value vector, by multiplying the embedding by three matrices trained during the training process. The three vectors are used to score the relevance of other elements in the sequence against the specific element. The intuition of self-attention is to keep the values of relevant elements intact and drown out irrelevant elements. Multi-headed attention refined the self-attention layer by adding a mechanism called "multi-headed" attention that runs through an attention mechanism several times in parallel. Intuitively, multi-head attention allows the model to jointly attend to information from different representation subspaces at different positions [START_REF] Vaswani | Attention is all you need[END_REF].

The output from the encoder is then forwarded to a fully connected layer to obtain the flow prediction results for different looking-ahead windows, i.e., future period 1, period 2,..., and period N p .

III. EXPERIMENTAL STUDY

To verify the efficacy of the proposed method, an experimental study has been carried out on the French airspace using one-month ADS-B data from December 1 to December 31, 2019, including a number of 158856 flights. This study focus on the en-route air traffic above 10,000 ft. The prediction target in this experimental study is set as: using the past one-hour traffic information in the airspace to predict the future number of flights that will transit through a flow coordination point in the coming 15 minutes, 30 minutes, 45 minutes, 1 hour, 1.25 hours, 1.5 hours, 1.75 hours, and 2 hours.

A. Flow Regulation Points Identification

By finding intersections of flight trajectories in the French airspace and clustering the intersection points to identify the major flow intersections on a daily basis, a graph representation of the daily air traffic flow pattern can be obtained. The nodes in the graph show the major flow intersections, and the edges describe the flow volume between these intersections. Fig. 4 shows the constructed graphs for the air traffic flow in the French airspace on December 1, 2019. The sub-figures 4a, 4b, 4c, and 4d display the constructed graph when the number of clusters is set as 300, 600, 1000, and 2000 respectively. Each green node in the graphs denotes a flow intersection point, and the node's size is proportional to the number of trajectory intersections in the corresponding cluster. Each red line represents an edge in the graph, while the line thickness is proportional to the number of flights traveled between the edge's start/end nodes. It can be observed from Fig. 4 that, on the one hand, when the number of clusters is small, the graph depicts air traffic flow features between city pairs, as most of the identified flow intersections are featured by the geographical locations of the metropolis and small cities. Such a graph may fail to represent important intersecting flows in the en-route phase of flights. On the other hand, when the number of clusters is large, the graph structure becomes complex as it may be affected by trivial interactions between the air traffic flow. In accordance with a set of different cluster numbers ranging from 100 to 1500, this paper calculated the geographical consistency gc 1 and the structural consistency gc 2 . Fig. 5 shows the changes of gc 1 and gc 2 which are scaled between 0 an 1 versus the varying number of clusters. A "saddle point" is observed for gc 1 and gc 2 around the cluster number 605. Therefore, this paper takes the value 605 as the number of clusters to be identified by the hierarchical clustering algorithm. The centers of identified clusters are determined as the flow coordination points.

C. Result Analysis

Given that there are over 600 flow coordination points identified, this section uses the coordination point, which is over the Paris area control center and handles the highest traffic volume, as an example to present the prediction result on the test dataset. hours respectively. The blue solid shows the actual number of flights passing the coordination point, while the red line shows the predicted value with the proposed transformer neural network prediction model. We can observe from Fig. 7 that the proposed method sustainably gives forecasts in close proximity to the actual flow value as the prediction window increases from 15 minutes to 2 hours. Furthermore, the model performs better in approximating the true flow value under a larger prediction window, as the prediction curve is more laminated on the actual curve for the 2-hour prediction window than the 15-minute window.

The proposed model has been compared with the canonical Multi-layer Perceptron (MLP) neural network. Blue dashes in Fig. 7 show the predicted values using MLP neural network. We can observe that the proposed model gives better predictions for all prediction windows.

Table I shows the quantified prediction performance of the proposed method, including four metrics: Mean Absolute I that the MAE, MSE, and R 2 values of the proposed method increase with the prediction window, while the MAPE value decreases with the prediction window. Fig. 8 visualizes the changes in the prediction performance of the proposed method versus MLP for varying prediction windows. It can be observed that the proposed method outperforms MLP in all of the metrics. The MSE, MAE, and MAPE values of the proposed method are lower than MLP, indicating more minor prediction errors. At the same time, the R 2 value is higher than MLP, showing a better approximation to the actual flow curve. The growth in the MSE and MAE values shows that the gap between the predicted value and the true value increases as the prediction window increases. The reason for this may be the increase in the number of flights passing the coordination points in a 2-hour duration compared to 15minute duration. As can be seen from Fig. 7, the maximum number of flights in a 15-minute duration is 24, while in a 2hour duration is 120. In contrast, the decreasing MAPE value, as a percentage, and the increasing R 2 , as a regression performance metric, indicate that the model's prediction accuracy improves as the prediction window increases. The reason may be that the number of flights transiting in the air space during the next 2 hours has less variation compared to a future 15minute duration, which is more stable and predictable.

IV. CONCLUSIONS

Aimed at contributing to the future flow-centric ATC paradigm, this paper proposed a transformer neural network approach for air traffic flow prediction at flow coordination points. Firstly, the flow coordination points were identified through hierarchical clustering of the flight trajectory intersections and the number of clusters was identified by analysing the consistency of daily traffic flow patterns modelled as graphs. Secondly, this paper described the air traffic flow as a "paragraph" of the flights sequences passing the coordination points and a transformer neural network model was adopted to learn from the contextual descriptions of air traffic and predict future air traffic flow at the coordination points. An experimental study on French airspace using the flight data from December 1-31, 2019. Results showed that the proposed prediction model approximated the true flow values with a R 2 between 0.909 to 0.99 and a MAPE varying from 27.4% to 11.7% for a 15-minute to 2-hour prediction window. The sustainability of the prediction performance under the increasing prediction window demonstrated the potential of the proposed model for longer-term flow prediction.

The proposed air traffic flow prediction method provides a graph-based approach to modelling air traffic flow patterns under the flow-centric paradigm. Predicting the future flow on the coordination points can inform flow-centric operations of the anticipated intersecting flows so that flow coordination can be applied in advance to avoid complex and conflicting situations involving a large number of flights. This paper's prediction of air traffic flow only considered flight numbers as the prediction target. However, the mixture of the intersecting flows can also affect the complexity and difficulty in flow management. A flow scenario with more crossing traffic is more challenging than a scenario with more parallel traffic, although the same number of flights are involved. Therefore, in the future, this traffic flow prediction method can be improved by considering the prediction of flow mixture to deliver a more meaningful prediction result.

Figure 1 :

 1 Figure 1: Conceptual diagram of the proposed method for coordination points flow prediction using transformer neural networks.

Figure 2 :

 2 Figure 2: An example of the intersection point clustering.

  The flight crossing time at point p r * k r is approximated as the flight crossing time at the closest point to p r * k r in T o .

Figure 3 :

 3 Figure 3: Proposed flow prediction model based on transformer neural networks.

Figure 4 :

 4 Figure 4: The constructed graphs for the air traffic flow in the French airspace on Dec 1, 2019 with different nodes.

Figure 5 :

 5 Figure 5: cluster Number determination.

Figure 6 :

 6 Figure 6: Comparison of the identified flow coordination points for two different days: 1st (red dots) and 2nd (green dots), December 2019.

Fig. 6

 6 Fig.6shows the identified flow coordination points for December 1 (red dots) and 2 (green dots), 2019. The sizes of the dots are proportional to the traffic volume passing the corresponding coordination points. It can be observed that the identified coordination points show consistent patterns in geographical distributions and traffic volumes, although there are some variations in their geographical locations due to the differences in daily air traffic flow organizations.

Figure 7 :

 7 Figure 7: Flow prediction results for Dec 30, 2019 over different time horizons as shown by the y-axis.

Fig. 7

 7 Fig.7presents the traffic flow prediction result from 00:00 to 23:59 on Dec 30, 2019. The eight panels plot the prediction results for the future 15 minutes, 30 minutes, 45 minutes, 1 hour, 1.25 hours, 1.5 hours, 1.75 hours, and 2 hours respectively. The blue solid shows the actual number of flights passing the coordination point, while the red line shows the predicted value with the proposed transformer neural network prediction model. We can observe from Fig.7that the proposed method sustainably gives forecasts in close proximity to the actual flow value as the prediction window increases from 15 minutes to 2 hours. Furthermore, the model performs better in approximating the true flow value under a larger prediction window, as the prediction curve is more laminated on the actual curve for the 2-hour prediction window than the 15-minute window.The proposed model has been compared with the canonical Multi-layer Perceptron (MLP) neural network. Blue dashes in Fig.7show the predicted values using MLP neural network. We can observe that the proposed model gives better predictions for all prediction windows.TableIshows the quantified prediction performance of the proposed method, including four metrics: Mean Absolute

Figure 8 :

 8 Figure 8: Prediction performance (MAE, MSE, MAPE, and R 2 ) of the proposed method versus MLP as the prediction window increases.

  the set of nodes in the graph constructed for day D k , and letV k+1 = v 1 k+1 , v 2k+1 , ..., v j k+1 , ..., v n k+1 represent the set of nodes in the graph G k+1 constructed for day D k+1 . For each node v i k in V k , find its nearest node v ai k+1 in V k+1 . By representing the (latitude, longitude) of node v i

k as (φ i k , λ i k ) and node v i k+1

TABLE I :

 I Prediction performance of the proposed method under different prediction windows in terms of MAE, MSE, MAPE, and R 2 .

		MAE	MSE	MAPE	R 2
	15 minutes 1.683 6.082	0.274	0.909
	30 minutes 2.339 12.184 0.217	0.951
	45 minutes 2.755 17.092 0.179	0.968
	1 hour	3.102 21.819 0.155	0.977
	1.25 hours	3.378 28.561 0.139	0.980
	1.5 hours	3.794 34.200 0.129	0.983
	1.75 hours	3.958 34.934 0.122	0.988
	2 hours	4.129 35.837 0.117	0.990

Error (MAE), Mean Squared Error (MSE), Mean Absolute Percentage Error (MAPE), and R-squared (R 2 ). It can be seen from Table

ACKNOWLEDGMENT

This research is supported by the National Research Foundation, Singapore, and the Civil Aviation Authority of Singapore, under the Aviation Transformation Programme. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not reflect the views of National Research Foundation, Singapore and the Civil Aviation Authority of Singapore.