Loïc Shi-Garrier

Daniel Delahaye

Nidhal C Bouaynaya

Predicting Air Traffic Congested Areas with Long Short-Term Memory Networks

Keywords: Air traffic flow prediction, Complexity metrics, Long short-term memory (LSTM) networks

The Extended ATC Planning (EAP) function aims at bridging the gap between Air Traffic Flow & Capacity Management (ATFCM) and Air Traffic Control (ATC) by predicting air traffic congested areas tens of minutes before their formation, and by suggesting real-time and fine-tuning measures to alleviate airspace "complexity". Current Air Traffic Flow Prediction methods focus on crude aircraft count in a sector, and hence are unable to distinguish between low and high complexity situations for a similar aircraft count. Complexity indicators, on the other hand, aggregate air traffic measurements and workload to describe the perceived complexity. However, the evaluation of workload is a long-debated issue and an inherently ill-posed problem. In this work, we present an intrinsic complexity metric, independent of any traffic control system, and address the prediction task of the EAP function using a novel Encoder-Decoder Long Short-Term Memory (LSTM) neural network. The complexity is measured by the eigenvalues of a linear dynamical system fitted to the aircraft's speed vectors. The Encoder-Decoder LSTM network uses a sequence of (discrete) aircraft states to predict the complexity of the air traffic in all areas of an airspace, in a time horizon of 40 minutes. Simulated traffic corresponding to one day of traffic over the French airspace is used to train and validate the model. Our experiments show that the proposed model achieves a Mean Absolute Error of 0.08 in predicting the normalized complexity value 40 minutes in the future.

I. INTRODUCTION A. Background and Motivation

With the continuous rise in air transportation demand (a 4.3% annual growth of the worldwide passenger demand is expected until 2035 [START_REF]Long-Term Traffic Forecasts[END_REF]), the capacity of the Air Traffic Management (ATM) system is reaching its limits, leading to increased flight delays (expected to achieve 8.5 minutes per flight in 2035 with the current system [2]). The Covid-19 crisis has severely impacted the aviation industry, with a 90% decline of RPK (Revenue Passenger Kilometer) in July 2020 with respect to the same period in 2019. However, the traffic is expected to regain its growth rate once the crisis is over, with the worst-case scenario predicting a return to normality before 2025.

Congestion is defined as a situation where a set of trajectories strongly interact in a given area and time interval, leading to potential conflicts and hence requiring high monitoring from the controllers. Air traffic complexity is a measure of the difficulty that a particular traffic situation will present to air traffic control [START_REF] Prandini | Toward Air Traffic Complexity Assessment in New Generation Air Traffic Management Systems[END_REF], and was shown to negatively affect the controller's decision-making ability and increase the error rate [START_REF] Pfleiderer | Relationship of complexity factor ratings with operational errors[END_REF]. To address this complexity and improve the service provision to airspace users through reduced delays, better punctuality, less ATFCM regulations, and enhanced safety, the Extended ATC Planning (EAP) function was introduced by SESAR JU as the Solution #118 [START_REF]SESAR-Joint-Undertaking[END_REF]. The EAP function relies on automation tools that enable early measures to be taken by ATC before traffic enters overloaded sectors. The implemented tools aim at bridging the gap between Air Traffic Flow and Capacity Management (ATFCM) and Air Traffic Control (ATC) by facilitating the communication between the local flow management position (FMP) and the controllers' working positions, and providing information to help taking actions at tactical time (less than one hour) such as rerouting or sector configurations. The gap between ATFCM and ATC can be described as follows: although ATFCM measures have managed to balance the capacity with the demand, high complexity areas may appear at the control level. These "hot spots" or "traffic bursts" require an intense surveillance from air traffic controllers in order to decide if some trajectories should be modified to avoid any conflict. Two problems must be addressed to reduce this increased workload: 1) the prediction of hot spots tens of minutes before their formation, and 2) the mitigation of the hot spots by appropriate early actions. This paper focuses on the first problem of the prediction of congested areas tens of minutes ahead of formation.

B. Related Work

We can group the literature to tackle the general problem of congestion prediction into three main approaches. The first approach considers conflict detection tools in a time horizon lower than thirty minutes [START_REF] Tang | Conflict Detection and Resolution for Civil Aviation: A Literature Survey[END_REF], [START_REF] Kauppinen | European medium-term conflict detection field trials[END_REF]. The second approach focuses on predicting the air traffic flow [START_REF] Sridhar | Aggregate Flow Model for Air-Traffic Management[END_REF], [START_REF] Cao | Link Transmission Model for Air Traffic Flow Management[END_REF], [START_REF] Chen | A network-based dynamic air traffic flow model for short-term en route traffic prediction[END_REF], [START_REF] Lin | Deep learning based short-term air traffic flow prediction considering temporal-spatial correlation[END_REF], [START_REF] Liu | Research on the Air Traffic Flow Prediction Using a Deep Learning Approach[END_REF], [START_REF] Gui | Machine Learning Aided Air Traffic Flow Analysis Based on Aviation Big Data[END_REF]. The third approach computes various air traffic complexity metrics [START_REF] Laudeman | Dynamic density: An air traffic management metric[END_REF], [START_REF] Sridhar | Airspace Complexity and its Application in Air Traffic Management[END_REF], [START_REF] Hilburn | Cognitive Complexity in Air Traffic Control: A Literature Review[END_REF], [START_REF] Chatterji | Neural network based air traffic controller workload prediction[END_REF], [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], [START_REF] Kopardekar | Dynamic density: Measuring and predicting sector complexity[END_REF], [START_REF] Masalonis | Dynamic density and complexity metrics for realtime traffic flow management: Quantitative analysis of complexity indicators[END_REF], [START_REF] Flener | Air-traffic complexity resolution in multi-sector planning[END_REF], [START_REF] Ishutkina | Describing Air Traffic Complexity Using Mathematical Programming[END_REF], [START_REF] Delahaye | Air traffic complexity based on dynamical systems[END_REF], [START_REF] Prandini | Toward Air Traffic Complexity Assessment in New Generation Air Traffic Management Systems[END_REF], [START_REF] Nguyen | Classification in functional spaces using the BV norm with applications to ophthalmologic images and air traffic complexity[END_REF], [START_REF] Garcia | Air traffic complexity map based on linear dynamical systems[END_REF].

In the first approach of short-term conflict detection methods [START_REF] Tang | Conflict Detection and Resolution for Civil Aviation: A Literature Survey[END_REF], the Medium-Term Conflict Detection (MTCD) flight data processing is perhaps the most popular system. MTCD is designed to warn the controller of potential conflicts between flights in a time horizon extending up to thirty minutes ahead. The system integrates predictive tools that performs trajectory prediction, conflict detection, trajectory update, and trajectory edition using "what-if" scenarios and tools [START_REF] Tang | Conflict Detection and Resolution for Civil Aviation: A Literature Survey[END_REF]. The MTCD applications, currently implemented in operational context, e.g., the Eurocontrol MTCD project [START_REF] Kauppinen | European medium-term conflict detection field trials[END_REF], are mainly based on pairwise conflicts rather than on a global approach. Moreover, for longer time horizons, the uncertainties on the trajectories make it difficult to predict the exact trajectories that will be involved in conflicts.

Air traffic flow prediction (ATFP) focuses on aggregate models [START_REF] Sridhar | Aggregate Flow Model for Air-Traffic Management[END_REF] rather than simulating the trajectories of individual aircraft. Trajectory-based models result in a large number of states, which are susceptible to error and difficult to design and implement for air traffic flow management [START_REF] Chen | A network-based dynamic air traffic flow model for short-term en route traffic prediction[END_REF]. The ATFP approach develops models of the behavior of air traffic that can be used for the analysis of traffic flow management. Several methods were considered, ranging from probabilistic algorithms modeling the flow over a network [START_REF] Sridhar | Aggregate Flow Model for Air-Traffic Management[END_REF], [START_REF] Cao | Link Transmission Model for Air Traffic Flow Management[END_REF], [START_REF] Chen | A network-based dynamic air traffic flow model for short-term en route traffic prediction[END_REF] to more recent machine learning algorithms [START_REF] Lin | Deep learning based short-term air traffic flow prediction considering temporal-spatial correlation[END_REF], [START_REF] Liu | Research on the Air Traffic Flow Prediction Using a Deep Learning Approach[END_REF], [START_REF] Gui | Machine Learning Aided Air Traffic Flow Analysis Based on Aviation Big Data[END_REF].

In [START_REF] Sridhar | Aggregate Flow Model for Air-Traffic Management[END_REF], a linear time variant traffic flow model was developed based on historical data. The dimension of this model depends only on the number of considered control volumes and not on the number of individual aircraft. This model is able to forecast the aircraft count at an Area Control Center level with a time interval of ten minutes. Another flow-based model was later developed in [START_REF] Cao | Link Transmission Model for Air Traffic Flow Management[END_REF], called Link Transmission Model. In this model, a flight path is defined as a sequence of directed links passing through sectors. An aircraft is assumed to cross each of the links within an estimated crossing time such that the state of each link can be easily tracked. Aggregation of links in each sector yields a traffic forecast for that sector. This approach is extended in [START_REF] Chen | A network-based dynamic air traffic flow model for short-term en route traffic prediction[END_REF] where the authors introduce a dynamic network of air traffic flow characterizing both the static airspace topology and the dynamics of the traffic flow. Historical data is used to estimate the travel time corresponding to the weight of each edge of the network. Then, a probabilistic prediction method is implemented for the short-term prediction (fifteen minutes) of the air traffic flow.

Recently, several machine learning models were proposed for ATFP. In [START_REF] Lin | Deep learning based short-term air traffic flow prediction considering temporal-spatial correlation[END_REF], the authors define a 3D grid over a given airspace containing the number of flights in each cell. They rely on neural network models that combine convolutional operations to extract spatial features and recurrent operations to extract time-dependent features. The 3D grids of the last ten minutes are inputted into the model to predict the 3D grid of the next timestep. A similar task is addressed in [START_REF] Liu | Research on the Air Traffic Flow Prediction Using a Deep Learning Approach[END_REF] using several 3D convolutional neural networks to extract spatial features for one timestep then recurrent layers to extract temporal features. In [START_REF] Gui | Machine Learning Aided Air Traffic Flow Analysis Based on Aviation Big Data[END_REF], the authors used support vector machines (SVM) and recurrent neural networks to predict the hourly air flow in predefined routes from time information (such as the time of the day, or the season and holiday indexes). The main drawback of all the above-mentioned ATFP methods is that they do not define any concept of congestion since the predicted variable is the aircraft count whether in a sector, a route, or crossing a predefined waypoint. Hence, these models are unable to discriminate low complexity situations from high complexity situations for a similar aircraft count.

To address this limitation, the third approach redefines the forecast objective of the ATFP problem by considering complexity metrics. Research focusing on defining and quantifying air traffic complexity, and on analyzing its impact on air traffic controller workload, was extensively conducted during the last two decades [START_REF] Laudeman | Dynamic density: An air traffic management metric[END_REF], [START_REF] Sridhar | Airspace Complexity and its Application in Air Traffic Management[END_REF], [START_REF] Hilburn | Cognitive Complexity in Air Traffic Control: A Literature Review[END_REF], [START_REF] Chatterji | Neural network based air traffic controller workload prediction[END_REF], [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], [START_REF] Kopardekar | Dynamic density: Measuring and predicting sector complexity[END_REF], [START_REF] Masalonis | Dynamic density and complexity metrics for realtime traffic flow management: Quantitative analysis of complexity indicators[END_REF], [START_REF] Flener | Air-traffic complexity resolution in multi-sector planning[END_REF]. Extensive reviews on this topic can be found in [START_REF] Prandini | Toward Air Traffic Complexity Assessment in New Generation Air Traffic Management Systems[END_REF] and [START_REF] Nguyen | Classification in functional spaces using the BV norm with applications to ophthalmologic images and air traffic complexity[END_REF]. The traditional measure of air traffic complexity is the traffic density, defined as the number of aircraft crossing a given sector in a given period [START_REF] Hilburn | Cognitive Complexity in Air Traffic Control: A Literature Review[END_REF]. The traffic density is compared with the operational capacity, which is the acceptable number of aircraft allowed to cross the sector at the same time period. This crude metric does not take into account the traffic structure and the geometry of the airspace. Hence, a controller may continue to accept traffic beyond the operational capacity, or refuse aircraft even though the operational capacity has not been reached. The Dynamic Density introduced in [START_REF] Laudeman | Dynamic density: An air traffic management metric[END_REF] aims at producing an aggregate measure of complexity by combining complexity factors, including static air traffic characteristics (such as airway crossings) and dynamic air traffic characteristics (such as the number of aircraft and the closing rates). A list of these complexity factors can be found in [START_REF] Hilburn | Cognitive Complexity in Air Traffic Control: A Literature Review[END_REF]. The complexity factors can be combined either linearly as in [START_REF] Sridhar | Airspace Complexity and its Application in Air Traffic Management[END_REF] or through a neural network [START_REF] Chatterji | Neural network based air traffic controller workload prediction[END_REF], [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF]. Several versions of Dynamic Density were advanced [START_REF] Kopardekar | Dynamic density: Measuring and predicting sector complexity[END_REF], [START_REF] Masalonis | Dynamic density and complexity metrics for realtime traffic flow management: Quantitative analysis of complexity indicators[END_REF]. For example, the interval complexity introduced in [START_REF] Flener | Air-traffic complexity resolution in multi-sector planning[END_REF] is a variant of Dynamic Density where the chosen complexity factors are averaged over a time window.

As elaborated in [START_REF] Prandini | Toward Air Traffic Complexity Assessment in New Generation Air Traffic Management Systems[END_REF], these complexity indicators aggregate air traffic measurements and workload to describe the perceived complexity. However, the evaluation of workload is a long-debated issue and an inherently ill-posed problem. The difficulty of obtaining reliable and objective workload measures is the main motivation for investigating complexity metrics that are independent of the ATC workload. To address this issue, another approach has been developed in the literature: the intrinsic complexity metrics approach [START_REF] Ishutkina | Describing Air Traffic Complexity Using Mathematical Programming[END_REF], [START_REF] Delahaye | Air traffic complexity based on dynamical systems[END_REF]. In [START_REF] Nguyen | Classification in functional spaces using the BV norm with applications to ophthalmologic images and air traffic complexity[END_REF], a difference is made between the control workload and the traffic complexity. The control workload measures the difficulty for the traffic control system (human or not) to remedy a situation, whereas the traffic complexity is an intrinsic measurement of the complexity of the air traffic, independently of any traffic control system. In other words, intrinsic complexity metrics aim at modelling the level of disorder and the organization structure of the air traffic distribution, irrespective of its effect on the ATC workload. In [START_REF] Ishutkina | Describing Air Traffic Complexity Using Mathematical Programming[END_REF], the authors interpolate a velocity vector field satisfying certain constraints, e.g., the field shall be equal to the velocity of an aircraft if an aircraft is present at this point, the field shall be flyable by an aircraft. The problem is modeled as a mixed integer linear program. The complexity is measured as the number of constraints that must be relaxed to obtain a feasible problem. In our previous work [START_REF] Delahaye | Air traffic complexity based on dynamical systems[END_REF], we defined an intrinsic complexity metric using linear and nonlinear dynamical system models. The air traffic situation is modeled by an evolution equation (the aircraft trajectories being integral lines of the dynamical system). The complexity is measured by the Lyapunov exponents of the dynamical system. The Lyapunov exponents capture the sensitivity of the dynamical system to initial conditions: if a small variation in the current air traffic situation leads to a very different dynamical system, the complexity of the situation is considered to be high.

This paper uses state-of-the-art machine learning techniques to integrate ATFP with the intrinsic complexity metric developed in [START_REF] Delahaye | Air traffic complexity based on dynamical systems[END_REF]. The main objective of this work is to provide congestion predictions in an extended time horizon, between forty to sixty minutes. Due to the uncertainties related to such time scales, we will directly predict congested areas rather than address the classical problem of trajectory prediction or conflict detection and resolution. In contrast to classic ATFP approaches that predict the aircraft count, we introduce a set of aircraft states (e.g., latitude, longitude, altitude, ground speed, heading, and climb rate) and use a Long Short-Term Memory (LSTM) network to predict congestion. We use an Encoder-Decoder framework where the input is the state of the aircraft and the output is the spatiotemporal complexity map measuring the spatial complexity at every time instant.

The remainder of the paper is organized as follows. Section II surveys the mathematically-grounded intrinsic complexity metric considered in this paper. Section III introduces the machine learning model, namely the encoder-decoder LSTM architecture. In Section IV, the data and the pre-processing techniques that were performed to build the training dataset are presented. Implementation details are also discussed for reproducibility of the results. The Python code is made available in GitHub1 . Section V elaborates on the results of the model. Considering only the areas with a nonzero complexity value, the best model achieves a Mean Absolute Error of 8% in predicting the complexity value forty minutes in the future. Finally, Section VI summarizes the paper and suggests several improvements to system deployment as well as new applications that could benefit from this work.

II. COMPLEXITY METRIC

We first introduce the mathematical notations used in the paper.

A. Mathematical Notation

Scalars are represented by lower-case letters, e.g., x, x i , y ij . Vectors are represented by bold lower-case letters, e.g., x, b. All vectors are column vector. Matrices are represented by bold upper-case letters, e.g., A, W. • F is the Frobenius norm. Re is the real part function. λ(A) is the spectrum of the matrix A. The concatenation of two vectors x and y into a single column vector is denoted by [x T , y T ] T where T is the transpose operator. σ is the logistic activation function

σ : x → 1 1+exp(-x) .

B. Complexity Metric based on Linear Dynamical System

In this work, the intrinsic complexity metric is based on a linear dynamical system model that was recently introduced in [START_REF] Garcia | Air traffic complexity map based on linear dynamical systems[END_REF]. This metric computes a measure of complexity in the neighborhood of an aircraft at a given time (see Figure 2). A filter is applied to consider only the flights that may interact with the reference aircraft. For example, an aircraft that is vertically separated with the reference aircraft by 1000 ft (in RVSM) will not interact with the reference aircraft, and, thus, should not be considered in the metric computation. After this Figure 4: A Long Short-Term Memory (LSTM) layer. The current input x t is combined with the previous hidden state h t-1 to compute the three gates used to update the cell state c t-1 and compute a new hidden state h t .

filtering, the selected aircraft positions are extended by adding p forward positions and p backward positions as depicted on Figure 2 with p = 2. The aim of this extension is to take into account uncertainties on the true aircraft's positions.

Let n ac be the number of aircraft samples retained for the computation of the metric. Consider the matrices P and V, which denote, respectively, the positions and velocities of the n ac aircraft, i.e.,

P =   x 1 x 2 . . . x nac y 1 y 2 . . . y nac z 1 z 2 . . . z nac   ; V =   v x1 v x2 . . . v xn ac v y1 v y2 . . . v yn ac v z1 v z2 . . . v zn ac   ,
where x i , y i , z i are the coordinates and v xi , v yi , v zi the speed components of the i th sample. A linear dynamical system ẋ = Ax + b is fitted to the positions P and speeds V of these aircraft, such that V AP + b (here, b, representing the mean of the speed vector field, is broadcast to be added to the matrix AP). The matrix A and the vector b are computed as the solutions of the Least Mean Squares problem

min A,b V -(AP + b) 2 F . (1)
The complexity metric, c(A), is then defined as the absolute sum of the negative real part of the eigenvalues of the matrix

A c(A) = Re(λ(A))<0 |Re(λ(A))|. (2)
Recall that the evolution equation of a linear dynamical system has the form P(t) ∼ exp(At). By diagonalising the matrix A, we can see that the asymptotic behaviour of the system depends uniquely on the eigenvalues of A.

Positive real parts of the eigenvalues correspond to diverging behaviour along the associated eigendirections while negative real parts correspond to convergence to a critical point (Figure 1). The imaginary parts correspond to rotation behaviour. Hence, our metric c(A) measures the strength of the converging or shrinking behaviour of the dynamical system. Since our system is the closest linear dynamical system fitting the current positions and velocities of the aircraft, a strong converging behaviour corresponds to rapidly converging trajectories, which we associate with high complexity.

This metric can be computationally efficient using known techniques from Numerical Linear Algebra (NLA) [START_REF] Lin | Approximating Spectral Densities of Large Matrices[END_REF] and Randomized NLA (RandNLA) [START_REF] Drineas | Lectures on randomized numerical linear algebra[END_REF].

III. ENCODER-DECODER LONG SHORT-TERM MEMORY NETWORK

In this Section, we present the sequential machine learning models that will process the series of aircraft state vectors. Sequence models are first introduced, then integrated into an encoder-decoder network architecture. Finally, a onedimensional convolutional layer used in the initial processing of aircraft state vectors is described.

A. Sequential Models

Recurrent Neural Networks (RNN) are a class of neural networks designed to process sequential data such as text streams, audio clips, and time-series data [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF], [START_REF] Saon | English conversational telephone speech recognition by humans and machines[END_REF], [START_REF] Vinyals | Show and tell: A neural image caption generator[END_REF]. In these recurrent models, network weights are shared across time, building an internal memory that can capture sequential knowledge. RNNs have been extensively used in Natural Language Processing, including automatic translation [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] and speech-to-text [START_REF] Saon | English conversational telephone speech recognition by humans and machines[END_REF].

In a basic RNN layer (see Figure 3), the recurrent cell receives the input x t as well as its own output from the previous time step y t-1 . The output y t is computed as

y t = φ(W • [x T t , y T t-1 ] T + w 0 ), (3)
where W is the weights matrix, w 0 is the bias vector and φ is an activation function e.g., tanh or ReLU. Observe that the output y t depends not only on the current input x t but also on the entire sequence of previous inputs x 0 , . . . , x t-1 using the same number of learnable parameters regardless of the length of the input sequence. However, the basic RNN layer has difficulty learning long-term dependencies when the length of the sequence is large. This problem has been referred to as the vanishing gradient problem [START_REF] Pascanu | On the difficulty of training recurrent neural networks[END_REF]. To solve it, several architectures of RNN layers with long-term memory have been introduced such as Long Short-Term Memory (LSTM) [START_REF] Hochreiter | Long Short-Term Memory[END_REF] and Gated Recurrent Unit (GRU) [START_REF] Cho | Learning phrase representations using RNN encoder-decoder for statistical machine translation[END_REF] networks. Long Short-Term Memory networks (Figure 4) were first introduced in [START_REF] Hochreiter | Long Short-Term Memory[END_REF], and then improved over the years [START_REF] Sak | Long Short-Term Memory Based Recurrent Neural Network Architectures for Large Vocabulary Speech Recognition[END_REF], [START_REF] Zaremba | Recurrent Neural Network Regularization[END_REF]. LSTMs detect long-term dependencies by using three gates that control the gradient propagation in the recurrent network's memory, i.e., the forget gate f t , the input gate i t , and the output gate o t . In contrast to the basic RNN layer, the LSTM has two state vectors denoted h t (the hidden state) and c t (the cell state), which can be seen, respectively, as the short-term state and the long-term state. The LSTM's governing equations are given by

(input gate) i t = σ(W i • [x T t , h T t-1 ] T + w i ), (4)
(forget gate)

f t = σ(W f • [x T t , h T t-1 ] T + w f ), (5)
(output gate) o t = σ(W o • [x T t , h T t-1 ] T + w o ), (6)
(candidate cell state) ct = tanh(W c • [x T t , h T t-1 ] T + w c ), (7)
(new cell state) c t = f t × c t-1 + i t × ct , (8)
(new hidden state)

h t = o t × tanh(c t ), (9)
where ). These gate controllers are simply computed by a linear mapping followed by a logistic activation. In parallel, a candidate cell state ct is computed using the current input and the previous hidden state (Eq. 7). Eq. 7 can be seen as the LSTM equivalent of the basic RNN layer. In Eq. 8, the current cell state c t is then updated by adding the previous cell state c t-1 (filtered by the forget gate) and the candidate cell state ct (filtered by the input gate). The hidden state h t is finally computed by activating the cell state c t with the tanh function (or another activation function, such as ReLU), which is then passed though the output gate (Eq. 9). With the input gate i t , the LSTM is able to store information in the cell state c t , to extract it with the output gate o t and to discard it with the forget gate f t .

W i , w i , W f , w f , W o , w o , W c ,

B. Encoder-Decoder LSTMs

The first introduction of an encoder-decoder model was provided for an automatic translation task [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF]. This model has also been applied to various tasks, such as speech recognition [START_REF] Prabhavalkar | A Comparison of sequence-to-sequence models for speech recognition[END_REF] and video captioning [START_REF] Venugopalan | Sequence to sequence -Video to text[END_REF]. The Encoder-Decoder LSTM is a recurrent neural network designed to address sequence-to-sequence problems. Sequence-to-sequence prediction problems are challenging because the number of items in the input and output sequences can vary. The Encoder-Decoder architecture is comprised of two models: one for encoding the input sequence into a fixed-length vector, and a second for decoding the fixed-length vector and outputting the predicted sequence. The innovation of this architecture is the use of a fixed-sized internal representation in the heart of the model, which may be referred to as sequence embedding. Both the encoder and decoder models are defined as recurrent neural networks, i.e., LSTMs.

The input of the decoder consists of two elements: the output of the encoder and the previously predicted (i.e., decoded) output sequence term. During training, the previously decoded sequence is provided as the ground-truth sequence at the previous timestep. During inference, the true output sequence is obviously not known. Hence the input of the decoder network has to be redefined. A simple method consists in recursively using the prediction of the last timestep as the input for the next one. However, this method may result in accumulation of prediction errors and lead to poor performance, especially if the output sequence is long enough. A solution to this issue is to change the task of the decoder network such that the network does not predict an element of the output sequence at each timestep, but rather a distribution over the possible outputs. Then, using a beam search algorithm [START_REF] Graves | Sequence Transduction with Recurrent Neural Networks[END_REF], it is possible to maintain, at each timestep, a small number of possible output sequences using the partial probability of each sequence as the metric. At the last timestep, the model outputs the sequence with the highest probability.

This framework is more suited for classification tasks where the output distribution is discrete and finite. Even if it is possible to choose a family of continuous distributions to represent the set of possible output distributions of the decoder network for regression, a simpler approach has been retained here. The objective of the model is to find the output sequence Y that maximizes the conditional probability P (Y|X) given the input sequence X. More formally, we are solving the following optimization problem:

arg max Y P (Y|X) = arg max Y t P (y t |y 1 , . . . , y t-1 , X).
The network directly outputs a prediction y τ = arg max yτ P (y τ |y 1 , . . . , y τ -1 , X) for the τ th element of the output sequence. Hence, at time step τ , the current estimation of the probability of the partial sequence y 1 , . . . , y τ -1 is τ t max yt P (y t |y 1 , . . . , y t-1 , X). This greedy approach can be seen as approximating max Y P (Y|X) = max Y t P (y t |y 1 , . . . , y t-1 , X) by t max yt P (y t |y 1 , . . . , y t-1 , X). This approach is equivalent to the beam search algorithm with a beam width of 1.

C. One-dimensional (1-D) Convolutional Layer

The first layers of the encoder network are defined as 1-D convolutional layers to process the input sequence of aircraft states. A 1-D convolutional layer is composed of f different filters. A filter contains a kernel of learnable parameters with dimension d × m, where m is the dimension of the input sequence, and d, the width of the kernel, is a hyperparameter of the convolutional layer. Each filter aggregates a subsequence of the input sequence of length d such that each element of the output sequence of the layer contains information from several successive timesteps. More precisely, we have:

y ij = φ( d k=-d l w j kl x i+k,l ),
where y ij is the j th vector element of the i th sequence term outputted by the 1-D convolutional layer, and x kl is the l th vector element of the k th sequence term of the input; w j kl is the k th weight of the kernel of the j th filter associated with the l th dimension of the input sequence, and φ is an activation function. If i + k is lower than zero or greater than the input sequence's length, then x i+k,l = 0, for all l.

The convolution operation is performed only along the time dimension of the input sequence; this is why this convolutional layer is referred to as "one-dimensional". With the 1-D convolutional layers, we expect the encoder network to identify dependencies in a very short time window, such that the LSTM layers could process an intermediate sequence in which each timestep contain dynamical information from the previous and following timesteps. This may facilitate the computation of the embedding for the whole sequence. Without this convolutional layers, the LSTM layers would process one timestep after the other without any information from the next timesteps.

IV. METHODOLOGY AND EXPERIMENTS

We now elaborate on the implementation details of our model, starting with the description of the dataset and the construction of the training set. Then, we provide the architecture details of the encoder-decoder model for reproducibility of our results.

A. Data acquisition and preprocessing

In order to test the capacity of the model to predict the complexity of the traffic, we use a dataset of simulated trajectories. In our simulation, no deconflicting actions were applied, such that each flight follows its flight plan without interacting with the other flights. Hence, we obtain the true complexity before mitigation, as opposed with historical data where the controllers have already applied measures to lower the complexity. The dataset represents a regular day of traffic over the French airspace. A trajectory is defined as a sequence of aircraft states. Each aircraft state consists of the following 8 elements. The trajectories are created by inputting flight plans to a simulator, where the flight plans consist of a sequence of waypoints and a cruise flight level. The atmosphere is assumed to be in ISA conditions with no wind. The final dataset, plotted in Figure 5, is composed of 8, 011 simulated trajectories spread along 3, 025 timesteps (1 timestep corresponds to 15 seconds). A complexity value is then computed for each aircraft at each timestep as described in Section II.

The objective of our supervised learning model is a sequence-to-sequence regression task. To achieve this, a training set must be built using pairs of training inputs/outputs, composed of sequences of aircraft states paired with sequences of complexity values over the entire airspace. To build the training outputs, we defined, for each time step t, a n × n matrix H t that will be referred as a complexity matrix. This matrix is superimposed over the airspace, such that each element of the matrix covers a small rectangle area whose dimensions depend on the airspace size and on the parameter n. To avoid a too high output dimension, we do not consider the altitude axis. Hence, each area covers all flight levels and is only defined by the latitude and longitude of its center point. Then, for a given area, the corresponding element of the matrix takes the maximum value of the complexity metric of the aircraft inside this area at timestep t. The complexity values are then scaled by a logarithmic function x → log(1+x/T h), where T h is a fixed threshold, to obtain a more uniform distribution. The H t matrices are then flattened as n 2 vectors to be considered by the machine learning model.

At a given time t, the input vector x t is built by concatenating all the aircraft states currently in the airspace, ordered by increasing longitude (and by increasing latitude if the longitude is the same). This means that a given flight will not occupy the same position in x t for every t. However, a given position in x t will be occupied by aircraft flying at roughly the same geographical coordinates such that the positional information is preserved in the structure of the inputs. To get a fixed-size input vector, the dimension of x t is set to the maximum number of simultaneous aircraft states in the dataset multiplied by 8 (the length of an aircraft state). If the number of states is lower than the maximum, the remaining elements of x t are padded with zeros.

We can know describe the training set as pairs of input sequences X t0 and output sequences Y t0 where

X t0 = x t0 • • • x t0+tin-1 and Y t0 = H t0+tin • • • H t0+tin+t pred -1
, where t in is the length of the input sequence, while t pred is the time horizon of the prediction (length of the output sequence). In the following, the dataset will be randomly divided between training examples (90% of the dataset) and validation examples (10% of the dataset).

B. Architecture of the model

The proposed encoder-decoder LSTM model predicts a sequence of complexity matrices using a sequence of concatenated aircraft states as input (see Figure 6). Let us first consider the encoder network. The input sequence is fed to one or several 1-D convolutional layers. The output of the 1-D convolutional layers is fed to several LSTM layers. The hidden state of the last LSTM layer is the encoding sequence, which is the only encoder information that will be passed to the decoder network. The decoder network is composed of several LSTM layers (the hidden states of these layers are initialized with the encoding sequence) followed Another characteristic of the proposed model is that the decoder network is implementing a technique called "teacher forcing". This means that, during training, the decoder network takes as an additional input the true output sequence, but shifted by one timestep. In other words, the decoder network also receives as input the output it should have predicted at the previous timestep. During the inference process, the decoder network simply considers the prediction it has made at the previous timestep.

Several variations of the dataset and the network are tested in this work. First, we consider the dataset where the complexity matrices H t are smoothed using a Gaussian kernel. Gaussian smoothing was performed using a 3×3 approximate Gaussian kernel 1 16 1 2 1 2 4 2 1 2 1 to smooth the values of the complexity metric across the airspace. Smoothed data is then used to train the model. We then investigated the effect of the length of aircraft trajectory in the computation of the linear dynamical system metric. Specifically, we added p points corresponding to future positions and p points corresponding to past position of each aircraft. This experiment aims at modelling the uncertainty of the aircraft's positions along the curvilinear abscissa. In particular, we studied whether the choice of the parameter p affects the performances of the model. The entire framework is summarized in Figure 7.

The hyperparameters of the implemented model are as follows. The convolutional layer has a kernel of width d = 3 and f = 512 filters. The encoder's LSTM layer has a hidden dimension of 128. The number of time steps for the encoder network is fixed at t in = 160. Since one time step is 5 seconds, this corresponds to 40 minutes, i.e., the model has access to the last 40 minutes of traffic to compute its prediction. The decoder's LSTM layer has a hidden dimension of 128 and is followed by one dense layer with output dimension 10, 000 (since the complexity matrix is a 100 × 100). The dataset (input and output) is linearly normalized between 0 and 1. We use the ReLU activation function for the last dense layer. The sequence predicted by the decoder network is set to t pred = 160, which corresponds to the next 40 minutes of complexity values. During end-to-end model training, we used the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with a learning rate = 10 -3 , and hyperparameters ρ 1 = 0.9, and ρ 2 = 0.999. The loss function is the Mean Squared Error defined as

M SE( Ĥ, H) = 1 n 2 n 2 i=1 ( ĥi -h i ) 2 . ( 10
)
The mini-batch size is set to 128 and the model is trained over 100 epochs. The 8, 011 simulated trajectories of our dataset are spread across 3, 025 time steps. Since our model requires pairs of sequences of respective lengths t in and t pred , we built 3, 025 -t in -t pred = 2, 705 sequences. Among these 2, 705 sequences, 2, 434 were part of the training set, and 271 were used to validate the model.

V. RESULTS AND DISCUSSION

Figure 8 shows the true and predicted complexity matrices, as heat maps, for a sequence from the training set. The horizontal axis corresponds to the longitude and the vertical axis is the latitude. The color indicates the complexity value associated to each area of the airspace. The complexity values have been rescaled to their original distribution. We can visually check that the model is able to predict the future complexity of the airspace, 40 minutes ahead of time, with a reasonable accuracy. Figure 9 shows the true and predicted complexity heat maps for a validation sequence, i.e., that the network did not encounter during training. It can be seen that the model is less accurate with validation examples. This may be due to overfitting, although the validation loss curve seems to converge to the same value as the training loss. Several regularization methods were tested, including L2 regularization, dropout, and Batch Normalization, but they significantly decreased the overall performance. The small size of the dataset (2,705 sequences) may explain why the model struggles to generalize, as well as the absence of the vertical information in the target output.

To quantify the performance of the model, we used the absolute error between the predicted density and the true density over the validation set. Given that the validation set contains 271 examples, each example containing 10, 000 density values, we obtain a total of 2, 710, 000 values of absolute errors, among which 182, 982 errors are nonzero (which is expected since the airspace is mostly empty). After removing the points with zero error, we plot the histogram of the non-zero absolute errors for three variants of the model, along with the corresponding Mean Absolute Error (MAE). Figures 10a and10b show the histograms for the model with p = 10 (the same model used earlier in this Section) and p = 2, respectively. Figure 10c shows the distribution of the non-zero errors when the output matrices are smoothed with a Gaussian kernel. Comparing the two variants of the linear dynamical system metric (with p = 2 or p = 10), the two distributions feature two modes other than zero. The distributions have a decreasing overall trend but with a certain number of errors around or above 1. However, with p = 2, the errors are spread across a smaller range and the MAE is smaller. There is also a smaller number of non-zero errors overall. This result is coherent with the literature findings that smaller sequences are somewhat easier to predict that longer ones, notably due to known gradient vanishing or exploding issues. When the complexity measures are smoothed, the MAE is five times smaller than nonsmoothed output, reaching a value of 0.08. The total number of non-zero errors is higher, which is due to the fact that there are more non-empty areas among the smoothed outputs. The prediction task with a smoothed dataset seems to be easier for the encoder-decoder LSTM model.

VI. CONCLUSION

In this paper, we tackled the problem of predicting air traffic congested areas. This objective was achieved by combining Air Traffic Flow Prediction methods with a more rigorous definition of the congestion relying on an intrinsic complexity metric based on a linear dynamical system. The proposed model is an Encoder-Decoder LSTM network that uses a sequence of (discrete) aircraft states to predict the complexity of the air traffic in all areas of an airspace, in a time horizon of 40 minutes. The model achieves better performance when the output of the dataset is Gaussian smoothed prior to training, achieving a Mean Absolute Error of 0.08 on the validation set.

The impact of this work is significant. The Encoder-Decoder LSTM model can provide information to the Flow Management Position operator about the upcoming hot spots. Then, the FMP operator can decide to regulate flights involved in the hot spots or modify the configuration of the control room in order to avoid the formation of the predicted hot spots. The model can also be fed as an input to a mitigation system, which would suggest actions on trajectories (e.g., changing heading, speed, or altitude) to prevent the formation of congested areas. The ultimate goal is to bridge the gap between the Flow Management Position's strategic planning and the controller's tactical monitoring as part of the Extended ATC Planning concept. This work is an important step towards suggesting upstream strategies to remove the hot spots that are too local and unpredictable at the strategic time frame, but still generate complex situations at the tactical level.

Future research directions include the quantification of the uncertainty in the predictions. The air transportation system is a highly dynamic system with inherent uncertainty and stochasticity. In order to be accepted in an operational context, the model should be able to detect if its predictions are no longer reliable and inform the FMP operator. The model should also be further refined and tested with larger datasets as well as using historical data. Moreover, the model can be improved by incorporating weather information such as wind and temperature which constitute the main causes of trajectory prediction uncertainties. Another possible approach consists in predicting the probability of formation of hot spots above pre-determined waypoints rather than predicting the precise locations of these hot spots.

Figure 1 :

 1 Figure 1: Eigenvalues loci of the dynamical system matrix A for 4 typical situations: parallel, convergent, divergent, and round about trajectories. Observe that the real part of the eigenvalues are negative if and only if the trajectories are convergent.

Figure 2 :

 2 Figure 2: The neighborhood of a reference aircraft to compute the complexity metric. This neighborhood is defined as a 24.8 NM × 24.8 NM box and ±30 FL altitude, centered on the reference aircraft. Trajectories can be augmented to account for the uncertainties along the curvilinear abscissa (dotted points).

Figure 3 :

 3 Figure 3: Basic RNN layer (left) and the same layer unrolled through time (right). Σ represents linear mapping and φ represents the activation function.Figure4: A Long Short-Term Memory (LSTM) layer. The current input x t is combined with the previous hidden state h t-1 to compute the three gates used to update the cell state c t-1 and compute a new hidden state h t .

Figure 5 :

 5 Figure 5: Simulated trajectories of the dataset plotted in the horizontal plane.

  State = [timestamp, aircraft ID, latitude, longitude, altitude, ground speed, heading, rate of climb].

Figure 6 :

 6 Figure 6: Encoder-decoder LSTM model for complexity prediction. The encoder network (left) uses a sequence of aircraft states x t0 , • • • , x t0+tin-1 to compute an encoding sequence. The decoder network (right) uses the encoding sequence and the complexity matrix outputted at the last timestep Ĥt-1 to compute a prediction.

Figure 7 :

 7 Figure 7: Depiction of the entire framework. The dataset is composed of trajectories where x i t is the state vector of the i th trajectory at time t. The complexity values are computed as described in Section II. The dataset, composed of the aircraft states and the complexity matrices, is divided into training and validation sets. The encoder-decoder LSTM model is trained (validated) with the training (validation) sets.

Figure 8 :

 8 Figure 8: A heatmap of the complexity matrices of one sequence of trajectories from the training set. The longitude is represented horizontally and the latitude vertically. True values after 40 minutes (left) and predicted values (right).

Figure 9 :

 9 Figure 9: A heatmap of the complexity matrices of one sequence of trajectories from the validation set. The longitude is represented horizontally and the latitude vertically. True values after 40 minutes (left) and predicted values (right)

  and w c are the learnable weights and biases of the LSTM network. The current input x t and the previous hidden state h t-1 are used to compute the controllers of three gates: the input gate controller i t , the forget gate controller f t and the output gate controller o

t (Eqs.

[START_REF] Pfleiderer | Relationship of complexity factor ratings with operational errors[END_REF][START_REF]SESAR-Joint-Undertaking[END_REF][START_REF] Tang | Conflict Detection and Resolution for Civil Aviation: A Literature Survey[END_REF]

https://github.com/lshigarrier/PFE

ACKNOWLEDGMENT

The authors would like to thank Ian Levitt and Andrew Cheng from the US. Federal Aviation Administration (FAA) William J. Hughes Technical Center, and Elsa Birman from the CRNA Est (Reims ACC) for their valuable input concerning ATM research and operational expertise. Nidhal C.