Improving recognition of complex aerial scenes using a deep weakly supervised learning paradigm Praveer Singh , Nikos Komodakis Abstract-Categorizing highly complex aerial scenes is quite strenuous due to the presence of detailed information with large number of distinctive objects. Recognition happens by first deriving a joint relationship within all these distinguishing objects, distilling finally to some meaningful knowledge that is subsequently employed to label the scene. However, something intriguing is whether all this captured information is actually relevant to classify such a complex scene ? What if some objects just create uncertainty with respect to the target label, thereby causing ambiguity in the decision making ? In this paper, we investigate these questions and analyze as to which regions in an aerial scene are the most relevant and which are inhibiting in determining the image label accurately. However, for such Aerial Scene Classification (ASC) task, employing supervised knowledge of experts to annotate these discriminative regions is quite costly and laborious; especially when the dataset is huge. To this end, we propose a Deep Weakly Supervised Learning (DWSL) technique. Our classification-trained Convolutional Neural Network (CNN) learns to identify discriminative region localizations in an aerial scene solely by utilizing image labels. Using the DWSL model, we significantly improve the recognition accuracies of highly complex scenes, thus validating that extra information causes uncertainty in decision making. Moreover, our DWSL methodology can also be leveraged as a novel tool for concrete visualization of the most informative regions relevant to accurately classify an aerial scene. Lastly our proposed framework yields state-of-the-art performance on existing ASC datasets.

Index Terms-Deep Learning, Aerial Scene Classification, Scene Complexity, Weakly Supervised, Multi-Instance Learning.

I. INTRODUCTION

Scene Classification is of paramount importance in remote sensing community to automatically categorize images for further scrutiny. These images are later on utilized by expert annotators for varied roles, including detection or segmentation of objects of interests.

One characteristic feature that holds Remote Sensing (RS) scenes distinctively apart from natural scenes is the widespread scale at which the area of interest is captured. Consequently, the extent of distinctive objects captured in a RS scene are also quite large. One of the challenges encompassing RS scene recognition is the added complexity when these discriminative objects are present simultaneously in multiple scenes. For e.g., in Fig. 1, a scene of Mountain (middle) has distinctive patches of grasslands (similar to Meadow on left) however the most representative parts of the scene are the white snowy tracts as would be illustrated in section IV.

Recently deep learning models have shown significant improvement in performance for varied remote sensing tasks such

The authors are affiliated to École des Ponts ParisTech & Université Paris Est, France (e-mail: {firstname.lastname}@enpc.fr).

Low Complexity

Mid Complexity High Complexity Fig. 1: Scenes representing various levels of complexity.

as aerial scene segmentation [START_REF] Volpi | Dense semantic labeling of subdecimeter resolution images with convolutional neural networks[END_REF][START_REF] Audebert | Beyond rgb: Very high resolution urban remote sensing with multimodal deep networks[END_REF], hyper-spectral classification [START_REF] Li | Classification of hyperspectral imagery using a new fully convolutional neural network[END_REF] or change detection [START_REF] Lyu | Learning a transferable change rule from a recurrent neural network for land cover change detection[END_REF]. Inspired by architectures of the primate visual cortex [START_REF] David | Receptive fields, binocular interaction and functional architecture in the cat's visual cortex[END_REF], these models interpret simultaneously various complex concepts which were lacking in conventional hand crafted methods [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF]. [START_REF] Li | What do we learn by semantic scene understanding for remote sensing imagery in CNN framework?[END_REF] defined scene complexity, based upon how much attention a person devotes to understand a particular scene. For e.g., in Fig. 1, a low complexity scene of Meadow is easy to interpret while comparatively more complex scenes of Mountain or Square have much more detailed information that needs to be distilled effectively and hence requires fairly larger period of attention.

However, instead of pivoting simply on greater attention for more complex scenes [START_REF] Li | What do we learn by semantic scene understanding for remote sensing imagery in CNN framework?[END_REF], we rather propose an alternate strategy. We argue that by limiting the amount of information gathered from a scene, we tend to minimize the ambiguity at the time of knowledge distillation and thus yield higher performance. Earlier we had seen that the amount of information to be processed differs from scene to scene, with highly complex aerial scenes exhibiting more number of distinctive regions. We postulate that only some of these regions are relevant for characterizing an aerial scene while the other regions are either irrelevant or inhibit in the recognition performance. This is specifically pertinent in the case of more complex scenes where the number of distinctive regions are significantly large and often lead to confusion in the final decision of a network. We, therefore, propose to remove this ambiguity using a novel methodology which allows the network to learn to choose the most and least relevant regions in a scene that aids in improving the overall recognition accuracy.

Nevertheless, selecting important discriminative regions in aerial scenes is a tedious task. This is mainly because images are of very high resolution and require an expert annotator. We choose to automatize the task of selecting relevant regions by introducing an end-to-end deep weakly supervised learning model in the context of aerial scene classification. Precisely, we make the network learn to select the most relevant regions in a scene that can predict the scene label with higher accuracy. Since we leverage only the class labels for localizing important regions in our aerial scene without object bounding box annotations, we call it weakly supervised. With further investigation, we also empirically showcase that RS scene recognition using this weakly supervised paradigm is much more beneficial for more complex scenes, thus substantiating our previous hypothesis of ambiguity removal.

Instances % OA Resnet-DWSL (without object bank) k=3,m=0 92.56 ± 0.22 k=3,m=3

93.02 ± 0.57 In a way, our proposed methodology is a relaxation of the prominent Multiple Instance Learning (MIL) [START_REF] Dietterich | Solving the multiple instance problem with axis-parallel rectangles[END_REF] technique where an image is characterized as a bag of instances (region). These bags or images are assigned a label: positive if it contains at least one positive instance and negative if all the instances are negative (Negative Instances in Negative bag or NiN). NiN is a fairly strong assumption simply because for e.g., absence of tennis court label doesn't imply that it is absent in the actual scene as seen in Square scene (rightmost in Fig. 1). It simply implies that the person who labeled the scene based his judgment upon the most predominant region in the scene. Thus, relaxing on the NiN presumption, our network tries to maximize the prediction of the correct class labels by utilizing both the top-K regions (positive instances) and bottom-M regions (negative instances).

Inspired from [START_REF] Hu | Transferring deep convolutional neural networks for the scene classification of high-resolution remote sensing imagery[END_REF], we showcase our own mechanism to visualize the prediction scores by highlighting the most discriminative regions which help to correctly classify any aerial scene. In a way, it furnishes us with a nice interactive tool to visualize the most important regions in a scene as simple as in a single forward pass through our network.

Lastly, we test our proposed methodology on all the existing aerial scene classification datasets [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF][START_REF] Zou | Deep Learning Based Feature Selection for Remote Sensing Scene Classification[END_REF][START_REF] Xia | Structural high-resolution satellite image indexing[END_REF] as well as on a fairly new and challenging one [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF], utilizing two prominent deep learning architectures namely Vgg-16 and Resnet-101. Experimental results clearly exhibit our Deep Weakly Supervised Learning (DWSL) method giving state-of-the-art (SoA) results both with Resnet-101 and Vgg-16 architecture. We also note that Resnet performs notably better than Vgg primarily In a nutshell, 1) We propose the first deep weakly supervised learning technique for aerial scene classification that automatically localizes most prominent regions using scene labels. 2) Our model removes ambiguity in complex aerial scenes by selecting only the most important objects out of a large number of distinctive objects that usually confuse the network in overall decision making. 3) We outperform SoA scene classification methods using only few relevant regions thus eradicating the need of utilizing entire scene as done in the past. 4) We present a compact tool for visualizing the most discriminative regions, thus giving us a better understanding of where network focuses upon while classifying a scene.

II. PROPOSED FRAMEWORK AND METHODOLOGY

As illustrated in Fig. 2a, our proposed Deep Weakly Supervised Learning (DWSL) pipeline consists of 3 major steps:

1) Feature Extraction using a VGG16 or Resnet-101 deep learning architecture. 2) Weakly Supervised Learning to select the most discriminative regions and jointly pooling them. 3) Object Bank Strategy to concatenate features from multiple scales and training an SVM classifier on top for ASC. We discuss each of the steps in more detail in following subsections.

A. Feature Extraction Module

Feature extraction module is employed to compute a fixedsize feature representation from the input image. We first perform a random scaling of the input image to a size between Motivated from [START_REF] Durand | WELDON: Weakly Supervised Learning of Deep Convolutional Neural Networks[END_REF], our WSL Module consists of two major components: Ws1 and K-max M-min Pooling (represented as s(.) in Fig. 2b).The Ws1 layer is similar to the fully connected layers and is often witnessed as terminal layers in most deep learning architectures. The only difference is that here we have used them in the form of multiple fully convolutional layers each of size 1 × 1. Further, these layers are applied individually to each of the 𝑝 × 𝑝 cells (as shown by the curved arrow) of the previous block. This operation yields us a 𝑝 × 𝑝 block but with a depth size equal to the class number 𝑐. Now, we treat the new output block as stacked heatmaps for each individual class with dimension of 𝑐 × 7 × 7 (corresponding to an input image of size 224). These class heatmaps are then fed to a K-max M-min Pooling layer 𝑠(.) which (a) selects the top-K scoring activations in a particular class heatmap (highlighted by green cells in Fig. 2b,(b) caters to lowest-M scoring activations in the heatmap (red cells).

B. Weakly Supervised Learning Module

Both top and lowest scoring regions are projected onto the original image and visualized as green and red bounding boxes in Fig. 7. Finally, these cells are aggregated individually for each of the 𝑐 class heatmaps using 𝑠(𝑐) which is given as: 𝑠(𝑐) = 𝑠 𝑡𝑜𝑝 (𝑐) + 𝑠 𝑏𝑜𝑡𝑡𝑜𝑚 (𝑐).

C. Object Bank Strategy

As a final step, we adopt an object bank strategy as highlighted in Fig. 2a. Similar to [START_REF] Durand | WELDON: Weakly Supervised Learning of Deep Convolutional Neural Networks[END_REF], we concatenate features computed from multiple scales of the input followed by training a support vector machine (SVM) classifier on top. This trained SVM model is then used to classify scenes from test dataset thus yielding final aerial scene classification score.

Hence, at the time of training, the K-max M-min pooling trains the network weights to accurately localize both the most discriminative regions (positive instances) that can correctly classify the scene plus those regions which have no correlation with the class (negative evidences). The reason for choosing negative instances is that some classes which have small interclass variability (Playground and stadium class in Fig. 3) might result in high classification scores for both by only using positive instances, while classifying a scene lets say of stadium class. This is due to the fact that there is presence of positive instances ( playground field and stadium roof tops respectively) for both the classes. In such a scenario, negative instances provide with complementary information (stadium roof giving clear evidence of absence of playground class) thus resulting in correct classification of the scene as stadium. This is elucidated empirically for with and without negative instances for our Resnet-DWSL model without object bank strategy (table I). Boost in overall classification score clearly highlights the necessity of using negative instances.

Method

Overall Accuracy (%) Pretrained GoogLeNet [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] 85.84 ±0.92 Hierarchical Coding [START_REF] Wu | Hierarchical coding vectors for scene level land-use classification[END_REF] 86.4 ±0.7 Pretrained VGG-VD-16 [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] 87.18 ±0.94 Pretrained CaffeNet [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] 88.25 ±0.62 Deep filter banks [START_REF] Wu | Deep filter banks for land-use scene classification[END_REF] 90.4 ±0.6 Two-stage deep feature fusion [START_REF] Liu | Scene classification based on two-stage deep feature fusion[END_REF] 92.37 ±0.72 Our DWSL method 96.19 ± 0.50

Method

Overall Accuracy (%) Pretrained GoogLeNet [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] 86.39 ±0.55 Pretrained CaffeNet [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] 89.53 ±0.31 Pretrained VGG-VD-16 [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] 89.64 ±0.36 𝑠𝑎𝑙𝑀 3 𝐿𝐵𝑃 -𝐶𝐿𝑀 [START_REF] Bian | Fusing local and global features for high-resolution scene classification[END_REF] 89.76 ±0.45 Combining 2 FC Layers [START_REF] Chaib | Deep feature fusion for vhr remote sensing scene classification[END_REF] 91.87 ±0.36 Two-stage deep feature fusion [START_REF] Liu | Scene classification based on two-stage deep feature fusion[END_REF] 94.65 ±0.33 Our DWSL method 95.06 ± 0.35

III. TRAINING AND IMPLEMENTATION DETAILS

A. Datasets

We conduct our experiments on 4 different datasets namely, AID [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF], UCMerced [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF], WHURS19 [START_REF] Xia | Structural high-resolution satellite image indexing[END_REF], RSSCN7 [START_REF] Zou | Deep Learning Based Feature Selection for Remote Sensing Scene Classification[END_REF]. AID is a fairly recent large scale dataset composed of 10000 images coming from 30 different classes, drawn from Google imagery captured using multiple imaging source. An additional level of complexity arises from the larger intra-class variations as each sample is collected from different locations over varied time period and seasons. The other 3 datasets are limited in number and somewhat saturated in terms of performance.

B. Training

Our experimental setup is somewhat similar to [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF]. We randomly draw 50% of our dataset as training set and rest is kept for testing. We repeat this step thrice to avoid the influence of randomness and compute precision accuracy for each run. We report overall mean and standard deviation over all runs in table II. We perform training on VGG-16 and Resnet-101 architectures as mentioned in section II.

For the baseline models (VGG-16 and Resnet-101), we first clip the final softmax layer + fully connected layer. Then, we add a fully connected layer with outputs equal to the number of classes present in the target dataset, finally appended with a softmax layer. Additionally, we pre-train these models on ImageNet and then, fine-tune them on target dataset.

For proposed DWSL models, we clip the last max pooling layer (just after convolution layers) + all the fully connected layers including soft max layer. Moreover, we simply add the Ws1 layer followed by K-max M-min pooling for Resnet-101 model. In case of VGG-16, an additional convolution layer of 7 × 7 × 2048 convolutions is added before the Ws1 layer. We add a softmax layer at the end for both the models. Our training loss is a log loss function which depends on the probabilities computed by softmax layer.

To effectively fine tune our networks, we utilize a dampening factor layer (with dt=10) to divide the back-propagating Fig. 6: Scenes from Center and Viaduct categories gradient just before it reaches the feature extraction module. This is mainly to assure that the pre-trained feature extraction weights are not tampered much. We train our models using Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] which computes adaptive learning rates for each parameter over the course of iteration. We set the initial learning rate to 1𝑒 -4. Each run consists of 50 epochs and we choose the best epoch based upon the test accuracy.

IV. RESULTS AND DISCUSSIONS A. Quantitative Results

Firstly, we compare the performance of our proposed model (DWSL) with the baseline model through overall mean accuracy over all runs as reported in Table II. Our weaklysupervised learning technique outperforms the current baselines by a considerable margin over all datasets. Especially in the case of VGG architecture, we witness a significant performance gain over the much recent AID [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] dataset.

Additionally for completeness, we also show comparisons with other state-of-the-art methods on AID and RSSCN7 benchmark datasets where we outperform all other handcrafted as well as deep-learning methods by considerable margin. Our results are summarized in tables III and IV. Our method is superior to second best results by a margin of 3.8% and 0.4% on RSSCN7 and AID datasets respectively.

The results clearly emphasize that learning to pay attention to few discriminative regions in a scene and recognizing using only these regions is much more beneficial than recognition by visualizing the entire scene. Another conclusion that can be drawn here is since the design of Resnet architecture naturally preserves spatial information throughout the network, it tends to learn more meaningful discriminative regions as compared to VGG network thus contributing to better performance.

Since Resnet architecture gives state of art performance on [START_REF] Xia | AID: A benchmark dataset for performance evaluation of aerial scene classification[END_REF] which is a fairly challenging dataset, therefore, we stick to Resnet to demonstrate the effectiveness of our method for all other well known datasets.

Next, we study the recognition accuracies from the perspective of scene complexities as depicted in Fig. 4. We observe that enhancement in performance is significantly higher in case of medium (5.66 point) and high complexity (3.22 point) regions than in case of low-complexity (2.88 point) regions. Delving deeper into class-wise predictions (in Fig. 5), we witness huge improvements in recognizing mid-level or high level complexity scenes such as Bridges (6.66 points), Parking (6.33 point), Mountain (12 point), Stadium (9 point), Airport (12.66 point) and Dense Residential (5.66 point). Thus we can fairly conclude that our weakly supervised methodology assists in minimizing ambiguity while recognizing more complex scenes where there is the presence of much more distinctive regions.

Point to note here is for the accuracy test results for high complexity scenes, we find that most of them have, on an average gain of at least 1.3 -2 percentage points which is quite significant considering that for many of them the baseline performance is already quite high. Only in two of the cases, i.e., for center and viaduct the performance becomes comparable. In case of Center class (Fig. 6), we observe that most of the scenes in this category have one big roof structure in the center of the image instead of large number of scattered objects throughout the scene. Thus it is more straightforward to classify a scene from the complete image rather than selecting few relevant regions. Similarly in case of viaduct, there are large-scale loopy structures located in the center of image which again makes both our method and baseline that uses full image competing enough.

B. Qualitative Results

Fig. 7 demonstrates qualitative results for our methodology by depicting aerial scenes (column 1), most and least distinctive regions marked by green and red boxes respectively (column 2) and the overlaid resulting heatmaps (column3). It is quite evident that the network selects only those relevant regions which are characteristic feature of a particular class. At the same time it also picks up insignificant regions which best demonstrate the absence of a class.

For eg. in case of the bridge scene, a neural network might confuse it with a river if overlooking at the entire scene. However by providing a weak localization supervision using image labels in our method, we tend to make the network decide on which are the most representative regions in a scene that best describe a scene. Similarly in case of mountain scene, the network, by looking at the entire scene, might interpret it to be a Bareland or Forests or Meadows, however by focusing on just the snow covered regions it tends to correctly classify it as Mountains. Thus we can conclude that making the network to localize distinctive regions and only focusing upon them for recognizing the scene, considerably removes the ambiguity caused while visually inspecting the complete scene.

V. CONCLUSIONS

In this paper, a deep weakly supervised technique is proposed in the perspective of aerial scene classification. We illustrate how a network can learn to localize the most predominant regions in an image simply from the scene labels. We also conclude that using these discriminative regions instead of the entire scenes results in state-of-the-art performance while at the same time providing meaningful interpretable information that allows one to elucidate better the decision of a network. We substantiate empirically as to how, by adopting our weakly supervised paradigm, we tend to remove ambiguity in more complex aerial scenes, resulting in a boost in their performance. In addition to this, our experimental results also underline that Resnet architectures suits much better to our task due to its characteristic feature of spatial information preservation. Finally we showcase a nice visualization tool to highlight most relevant regions in aerial scenes.

Fig. 2 :

 2 Fig. 2: Overall Pipeline and Network Architecture

Fig. 3 :

 3 Fig. 3: Scenes from playground and stadium categories

Fig. 5 :

 5 Fig.5: Class-wise accuracies for VGG-16 (baseline) and VGG-DWSL (our method) for the AID dataset.

Fig. 7 :

 7 Fig. 7: Qualitative Results for Bridge and Mountain.

  This is followed by cropping the image to a fixed size of 224×224. We perform several other data-augmentation strategies like horizontal and vertical flip or rotation by a small angle theta to avoid over-fitting. The resulting image is then passed through the convolution layers of either Resnet-101 or Vgg-16. In Resnet-101 model, the final convolution layer results in an output of size 2048 × 7 × 7 which can be treated as the first block of the Weakly Supervised Learning (WSL) Module (see Fig.2b). In case of Vgg-16, the output from the final convolutional layer is of size 512 × 14 × 14 which is then fed to WSL module as input.

			Levels of Complexities
	% OA	90 95
	low-complexity	mid-complexity high-complexity
		VGG-16			VGG-DWSL
	Fig. 4: Overall Accuracies (OA) for different levels of complexities
	[224, 244]. Dataset	Arch.		Baseline	DWSL	% Gains
		AID [12]	VGG		91.33 ±0.46	95.06 ±0.35	4.11
		AID [12]	Resnet-101	95.44 ±0.26	96.96 ±0.34	1.52
	WHURS19 [11]	Resnet-101	96.94 ±0.92	98.67 ±0.61	1.73
	RSSCN7 [10]	Resnet-101	94.48 ±0.33	96.19 ±0.50	1.71
	UCMercred [6]	Resnet-101	97.24 ±0.35	97.81 ±0.26	0.57

TABLE II :

 II Results comparing baseline with our DWSL method

TABLE III :

 III Performance Comparison Over RSSCN7

TABLE IV :

 IV Performance Comparison Over AID dataset


Man Liang

Daniel Delahaye

Pierre Marechal

Conflict-free Arrival and Departure Trajectory Planning for Parallel Runway with Advanced Point-Merge System

Keywords:

In this paper, an efficient trajectory planning system is proposed to solve the integration of arrivals and departures on parallel runways with a novel route network system. Our first effort is made in designing an advanced Point Merge (PM) route network named Multi-Level Point Merge (ML-PM) to meet the requirements of parallel runway operations. Then, more efforts are paid on finding a complete and efficient framework capable of dynamically modelling the integration of arrival and departure trajectories on parallel runways, modelling the conflict detection and resolution in presence of curved trajectory and radius-to-fix merging process. After that, a suitable mathematical optimization formulation is built up. Receding Horizon Control (RHC) and Simulated Annealing (SA) algorithms are proposed to search the near-optimal solution for the large scale trajectories in routine dense operations. Taking Beijing Capital International Airport (BCIA) as a study case, the experimental results show that our system shows good performances on the management of arrivals and departures. It can automatically solve all the potential conflicts in presence of dense traffic flows. With its unique ML-PM route network, it can realize a shorter flying time and a near-Continuous Descent Approach (CDA) descent for arrival aircraft, an economical climbing for departure aircraft, an easier runway allocation together with trajectory control solutions. It shows a good and dynamic sequencing efficiency in Terminal Manoeuvring Area (TMA). In mixed ML-PM mode, under tested conditions, our proposed system can increase throughput at BCIA around 26%, compared with baseline. The methodology defined here could be easily applied to airports worldwide.

Introduction

Motivation

Flight delays lead to many negative impacts, such as increasing passengers trip costs, airlines operation costs, environmental damage due to more fuel consumption and gas emissions [START_REF] Ball | Total delay impact study[END_REF][START_REF] Ryerson | Time to burn: Flight delay, terminal efficiency, and fuel consumption in the national airspace system[END_REF]. The key cause of delays is the disequilibrium between the high traffic flow demand and the low capacity. In 2015, China's total civil aviation turnover, passenger turnover, and cargo-mail turnover reached 85.165 billion ton-kilometers, 728.255 billion passenger-kilometers, and 20.807 billion ton-kilometers, a year-on-year increase of 13.8%, 15.0%, and 10.8%, respectively [START_REF] Caac | China national civil aviation efficiency report[END_REF]. However, less than 30% of the airspace is reserved for civil aviation operations in China, compared to about 80% of the airspace in the US. Meanwhile, a large amount of radar vectorings are frequently used by approach controllers to manage dense arrivals and departures in China. The current Standard Instrument Departures (SIDs) and Standard Instrument ARrivals (STARs) route design method does not provide a good use of airspace. Area Navigation (RNAV)/Required Navigation Performance (RNP) routes are not being fully exploited in today's baseline operation. As a result, there is an inefficient conflict resolution with a number of low-level radar vectoring, a low runway utilization, and a low flight efficiency in Terminal Manoeuvring Area (TMA). Consequently, at the first seven busiest airports in China, their peak traffic volumes exceeded their capacities by more than 30%. Average landing delay at the ten busiest airports is around 26 minutes, and average flight holding time on-ramp from closing door to pushing back is around 12 minutes [START_REF] Caac | China national civil aviation efficiency report[END_REF]. The air traffic bottleneck in the terminal airspace around the busy airport is an important issue nowadays.

On the other hand, parallel runways are the main structure of Chinese Top 16 airports, for example Beijing Capital International Airport (BCIA) with 3 parallel runways, Shanghai Pudong airport with 4 parallel runways, Guangzhou Baiyun airport with 3 parallel runways, Chongqing airport with 3 parallel runways. More and more parallel runways are being built in existing or new hub airports, such as Beijing DaXing airport with 5 parallel runways. Facing high demand in the near future, there is an urgent need to develop a novel approach to efficiently manage departing and arriving flows at busy airports with parallel runways, to significantly improve flight efficiency and increase capacity.

Furthermore, the total volume of airspace for civil aviation in China will not be changed in the near future, due to state security reasons. The volume of TMA in China is usually very small, but with huge traffic demand. To solve this problem, firstly, it is necessary to design an efficient and integrated RNAV/RNP-based SID/STAR route network. It should not only reduce the traffic complexity, but also adapt the constraint of small airspace, and support an optimized descent and climb in TMA. Secondly, a more dynamic re-sequencing optimization should be studied to achieve an integrated arrival and departure optimization. It should improve the predictability at busy airports on the strategic level, further increasing capacity. Thirdly, in dense traffic cases, automated conflict resolution by planning conflictfree trajectories in advance could potentially be a good solution. It should provide a good way to continuously reduce the traffic complexity and increase capacity. Especially, in a complex traffic situation, a good and feasible automation solution for controllers and pilots without requiring specialized equipage is crucial from an operation point of view. The system developed in this paper will be applied to BCIA, aiming to alleviate the heavy delay in China in the near future.

Traffic at BCIA

BCIA is the most important, largest and busiest international aviation hub in China. Everyday, around 1700 flights from 94 airline companies operate at BCIA, which makes Beijing in a close bond with 244 cities in 54 countries around the world. There are three parallel runways at BCIA: runway 01-19, runway 18L-36R, and runway 18R-36L. According to published Aeronautical Information Publication (AIP), runway 36L-18R is used for departure and arrival, runway 36R-18L is mainly used for departure, runway 01-19 is mainly used for arrival. Three runways could be all used for landing in case of dense traffic demand. The declared operation capacity at BCIA is 88 movements per hour, including departures and arrivals. As shown in Fig. 1, its peak traffic period is from 6:00 to 24:00. From 6:00 to 9:00, there is a heavy demand for departure. From 21:00 to 24:00, there is a heavy demand for arrival. Between 9:00 to 21:00, the arrival and the departure almost equals each other.

place Fig. 1 about here In Beijing TMA, there are designed SID routes for guiding the departing aircraft into the airways, and STAR routes for helping the arrivals aircraft to approach the airport. There are six entry points in Beijing TMA for arrival flights. Among which, there are four entry points in the South: JB, BOBAK, VYK and DOGAR, two entry points in the North: KM and GITUM. The traffic flows are very heavy in the South. Meanwhile, there are seven exit points for departure flights, including YV, CDY, TONIL, LADIX, RENOB, SOSDI, and KM. The significant point KM is special, it serves both for departure and arrival flights.

place Fig. 2 about here To understand baseline traffic at BCIA, several days of real flight data in the June, 2017 (see Fig. 2) were collected. We chose South-inbound arrivals and departures to study. Furthermore, we filtered out the 10th day traffic to study the South-inbound traffic characteristic at BCIA. In order to reduce the noise trajectories for each arrival flow, the initial data were clustered with a preprocessing algorithm. In the end, the South-inbound traffic characteristic is demonstrated in Fig. 2 "South-inbound arrivals on 10th day". It is found that flights from JB will normally join BOBAK firstly, then proceed together with the flows over JB to land. Flights from VYK have two possibilities to land, either join the flows from JB or join the other flows from DOGAR. Flights from GITUM have very few manoeuvring airspace. As a result, they execute frequently the holding procedure. However, flights from KM normally go direct to land with very few manoeuvring actions. Besides, there is a big manoeuvring space near the KM, which is used to absorb the traffic flows from JB, BOBAK, and VYK. There is another small manoeuvring space near the runway 01-19, which is used to absorb the traffic flow from GITUM, VYK and DOGAR.

State of the art

Arrival and Departure Manager

In industry, Arrival Manager (AMAN) system has been developed to support controllers' daily work for many years. In Europe, arrival management is a general term given to the process of safely and effectively arranging arrivals into a smooth efficient flow for landing at a destination airport [START_REF] Eurocontrol | Arrival Manager, Implementation GUIDELINES and Lessons Learned[END_REF]. Basic AMAN system provides simple Time to Lose/Time to Gain solutions to controllers, rather than more complex direct trajectory management solutions, such as airspeed to be flown. High level aims of AMAN are to assist the controller to optimise the runway capacity (sequence) and/or to regulate/manage (meter) the flow of aircraft entering the airspace, such as a TMA. It also aims to provide predictability for its users (both ground and air) and at the same time minimise the impact on the environment, by reduced holding and low-level vectoring [START_REF] Eurocontrol | Arrival Manager, Implementation GUIDELINES and Lessons Learned[END_REF].

Point Merge System (PMS) developed by the EUROCONTROL Experimental Center is a systemized method for sequencing arrival flows. It is a RNAV-based route topology design, using predefined sequencing legs at iso-distance to the center merging point for path shortening or stretching. [START_REF] Boursier | Merging arrival flows without heading instructions[END_REF] proved that this conventional PMS with only one merging point is efficient, safe and accurate in a small-scale experiment. Then, [START_REF] Ivanescu | Models of air traffic merging techniques: evaluating performance of point merge[END_REF] compared the PMS method with a radar vectoring method by using fast time simulation. The results showed that the PMS model reduced the mean controller task load by 20%, and the number of instructions to pilots by 30%, compared with vectoring. Nowadays, PMS has been successfully implemented in different airports in the world, such as Oslo, Dublin, Seoul, Paris, and London City. The topology of PMS reduces the traffic complexity in TMA, requires smaller manoeuvring airspace compared with conventional radar vectoring method, less workload for controllers and pilots. It is very suitable for handling dense traffic in a small TMA.

Nowadays, AMAN/Departure Manager (DMAN) are frequently used together at busy airports to support a good planning and implementation of arrival/departure sequence for controllers. They aim to increase throughput and predictability at an airport (Eurocontrol, 2017). Because, AMAN and DMAN have too often been developed separately, the integrated AMAN/DMAN is an important research topic in SESAR ATM Master Plan. Based on it, some key constraints in current AMAN and DMAN systems need to get more attentions from researchers. Firstly, it is found that the sequence proposed by AMAN is constrained by Maximum Position Shifting (MPS), which generally allows no more than three (MPS ≤ 3) [START_REF] Balakrishnan | Scheduling aircraft landings under constrained position shifting[END_REF]. AMAN mainly reproduces the sequence for aircraft who have conflict to land on the same runway, and indicates the Time to Lose/Time to gain information for controllers, but it can not provide an appropriate and feasible action for controllers/pilots to meet the target time or position in the sequence. Moreover, the re-sequencing between two arrivals on the same route that are already in-trail of one another generally is not allowed in current AMAN. Secondly, the departure re-sequencing is realized only before take-off. Large number of re-sequencing of aircraft in the take-off order will increase workload for the pilots and controllers by frequently modifying the routine taxi route. Clearly, a more relaxed re-sequencing will shorten the landing and departing make span, however, it is hard to realize in reality. Integrated AMAN/DMAN with PMS will support a more efficient traffic synchronization in TMA. It provides the ability to improve the overall arrival and departure management in high density traffic environment.

Parallel Runway Operations

In TMA, the runway systems at major airports are highly constrained resources. The effective use of airport capacity depends on the best use of all available runways [START_REF] Berge | The multiple runway planner: Modelling and analysis for arrival planning[END_REF]. The main objective of multiple runway operations is to increase runway capacity and aerodrome flexibility. The largest increase in overall capacity requires the use of independent approaches to parallel or near-parallel runways1 . According to ICAO-Doc.9643 (2004), there are four different modes of operational concepts relating to simultaneous operations on parallel or near-parallel instrument runways:

1. Independent parallel approaches are simultaneous approaches to parallel runways where radar separation minima are not prescribed between aircraft using adjacent Instrument Landing Systems (ILSs). 2. Dependent parallel approaches are simultaneous approaches to parallel runways where radar separation minima between aircraft using adjacent ILS are prescribed. 3. Independent parallel departures are independent instrument departures for aircraft departing in the same direction from parallel runways. 4. Segregated parallel approaches and/or departures are simultaneous operations on parallel or near-parallel instrument runways in which: one runway is used exclusively for approaches, and the other runway is used exclusively for departures.

Controllers can use simultaneous operations on parallel runways following specific rules and depending on the airport layout. Where parallel instrument runways are intended for simultaneous use, the minimum distance between their centre lines should be:

1. 1035m (3400ft) for independent parallel approaches 2. 915m (3000ft) for dependent parallel approaches 3. 760m (2500ft) for independent parallel departures 4. 760m (2500ft) for segregated parallel operations If the minimum distance between each runway centre line is below 760 m, the simultaneous operation for take-off and landing is not possible. The runways are considered as a single runway with regard to vortex wake separation. Specially, in segregated parallel operations, if the arrival runway is staggered toward the arrival aircraft, then the minimum distance between the centre lines may be decreased by 30 m (98 ft) for each 150 m (500 ft) to a minimum of 300 m (984 ft), if the arrival runway is staggered away from the arriving aircraft, the minimum distance between the centre lines should be increased by 30 m (98 ft) for each 150 m (500 ft).

For airport with multiple parallel runways, such as BCIA, the runway assignment is typically dependent on the airport configuration, the wind direction, the direction of arriving aircraft, and departure routes of the aircraft [START_REF] Brinton | An implicit enumeration algorithm for arrival aircraft[END_REF]. For arrivals, adjustments can be made to the flight plan by assigning the aircraft to an alternative runway; for departures, adjustments can also be made by issuing another taxi route for aircraft to an alternative runway. It is known as runway allocation problem. For controllers, the arrival flows into and out of the TMA are asymmetric due to different airlines scheduling. With runway re-assignment (runway allocation), the runway landing or departing could be balanced at all available runways, so as to reduce the unnecessary flight time and taxiing time, minimize the delays and maximize airport capacity. For instance, [START_REF] Kim | Runway assignments that minimize terminal airspace and airport surface emissions[END_REF] analysed runway assignments that minimize terminal airspace and airport surface emissions problem. [START_REF] Vela | Evaluation of strategic and tactical runway balancing[END_REF] studied the problem of strategically balancing departure demand at runways in order to reduce departure delays at airports with multi-runway configuration. [START_REF] Delsen | Flexible arrival and departure runway allocation using mixed-integer linear programming: A schiphol airport case study[END_REF] focused on research in flexible arrival and departure runway allocation by using Mixed-Integer Linear Programming (MILP) to optimize fuel and noise. These tools could provide an optimized strategic schedule of runway allocation, however they could not indicate the solution to realize the real-time runway re-assignment at the tactical level. To realize a successful runway re-assignment decision, the flight procedures on each STARs and SIDs need to be considered [START_REF] Samà | Rolling horizon approach for aircraft scheduling in the terminal control area of busy airports[END_REF][START_REF] Samà | Optimal aircraft scheduling and routing at a terminal control area during disturbances[END_REF]. A good route network should be established to support the aircraft to land on the assigned runway with short flying time in TMA.

Tactical Trajectory Control

As TMA is a transition airspace between airports and the network of airways, there are a large number of climbing and descending aircraft, the traffic complexity in TMA is very high. The way of trajectory controlling at each STAR and SID route segments has to consider not only adapting well the strategic schedule, but also operation streamlining around the critical points, such as merging points and intersection points. Conventional method of trajectory controlling in the TMA mostly relies on radar vector, such as heading change, speed control, flight level change. The more aircraft are operating in the TMA, the higher workloads are for controllers and pilots. Some innovative approach procedures, the Staggered Approach Procedure (SGAP) and the Steeper Approach Procedure (SEAP), were proposed in combination with the baseline conventional approach procedure under Instrument Meteorological Conditions (IMC) and Instrumental Flight Rules (IFR) at two large airports, one in Europe and other in the US [START_REF] Janic | Modelling the capacity of closely-spaced parallel runways using innovative approach procedures[END_REF]. An analytical model for calculating the ultimate arrival, departure, and mixed operation capacity of closely-spaced parallel runways was successfully built up.

Nowadays, based on new emerging technologies of Communication, Navigation and Surveillance (CNS), sometimes the controller can choose control time of arrival for the aircraft, the aircrew then uses the aircraft Flight Management System (FMS) to fly an optimized trajectory towards the target Required Time of Arrival (RTA). Avionics-based Air Traffic Management (ATM) operation produces a significant benefit for more efficiently 4D trajectory controlling in TMA [START_REF] Haraldsdottir | Arrival management with required navigation performance and 3d paths[END_REF][START_REF] Scharl | A fast-time required time of arrival (rta) model for analysis of 4d arrival management concepts[END_REF][START_REF] Korn | 4 dimension trajectory management in the extended terminal manoeuvring area: coupling arrival manager and 4 dimension flight management system for optimized approach trajectories[END_REF]. However, these benefits of avionics-based operation produces are limited, they can be mainly achieved in a relatively orderly and steady traffic flow. They are hard to be achieved in a dynamic and heavy congestion traffic situation. Because, the mainly de-conflict solution provided by FMS system is speed regulation. In targeting a single control time via speed control, the FMS solution only ensures that the conflict-free at the control point itself and not along the entire trajectory leading to that point. In a relatively small time horizon, the effect of speed regulation is limited. Other airborne-spacing solutions, such as heading change or flight level change, are not frequently proposed in a congestion situation, because of the traffic complexity.

Automation can achieve high efficiency in ATM. For example, [START_REF] Prevot | Effects of automated arrival management, airborne spacing, controller tools, and data link[END_REF] analyzed the effects of automated arrival management, airborne spacing, controller tools, and data link. Its benefits on noise abatement, fuel consumption, CO 2 emission have been demonstrated as well [START_REF] Coppenbarger | Field evaluation of the tailored arrivals concept for datalink-enabled continuous descent approach[END_REF][START_REF] Erzberger | Automated conflict resolution, arrival management and weather avoidance for atm[END_REF]Toratani et al., 2015b,a). But, due to safety consideration, full automation in ATM is still a challenge. Instead, it was considered by realizing an automated merging with conflict free for a large portion of routine traffic at congested airport [START_REF] Isaacson | A concept for robust, high density terminal air traffic operations[END_REF]. Therefore, for modernization of current ATM system with high level of automation, the routine trajectory Degree Of Freedom (DOF) employed in solutions must be carefully chosen. For example, the commercial aircraft are flying in TMA with a typical speed profile predefined in FMS. Below 10000 ft, the application of speed regulation to solve conflict is very limited. A Calibrated Air Speed (CAS) constant descent is widely used by FMS between Top of Descent (TOD) in en route airspace and 10000 ft in TMA. A feasible speed regulation in TMA for automated conflict resolution should be applied to CAS above 10000 ft. Another example, in a fuel optimal CDA and Continuous Climb Operation (CCO), there is no level-off, but in actual operation, a short level-off segment could be useful in reducing uncertainty in CDA trajectories proposed fully by airborne FMS. If there is some level-off in the automated conflict resolution, in the view of fuel consumption, it will be better to choose high altitude level-off, and the time of level-off should be as shorter as possible. Therefore, choosing the right control variables is crucial for a successful automation in trajectory control.

Contributions

To address the problem of 4D trajectory planning for parallel runway at busy airports, the authors began with solving the arrival management optimization problem. Firstly, designing an advanced topology of route network is important for reducing the airspace complexity and solving abilities of humans to manage recovery from off-normal events. [START_REF] Liang | A framework of point merge-based autonomous system for optimizing aircraft scheduling in busy TMA[END_REF] introduced the framework of an autonomous system with ML-PM for BCIA. Based on the classic PMS, which is used only to arrival management for single runway operation, the proposed ML-PM route network aimed at adapting to the segregated parallel runway operation. Hybrid algorithm with Receding Horizon Control (RHC) and Simulated Annealing (SA) is applied to find the solution, however the runway re-assignment problem was not considered. Then, [START_REF] Liang | Potential operational benefits of multi-layer point merge system on dense TMA operation[END_REF] made a brief evaluation of this ML-PM potential operational benefits in terms of flight efficiency and runway throughput, but the airspeed model and the fuel consumption model needs to be further improved. Later, [START_REF] Liang | Integrated sequencing and merging aircraft to parallel runways with automated conflict resolution and advanced avionics capabilities[END_REF] improved the previous model and proposed an advanced modelling to handle high-density arrivals to mixed parallel runways operation, meanwhile with more accurate airspeed and fuel consumption models.

Up to now, PMS is only used for arrival management, none of literatures has considered to use it for departure management. In fact, for departure aircraft, how to realize an easy transition in airspace from runway to different exit points of TMA is an important research issue as well. In addition, [START_REF] Liang | Integrated sequencing and merging aircraft to parallel runways with automated conflict resolution and advanced avionics capabilities[END_REF] proved that the application of ML-PM system could realize a more relaxed re-sequencing for arrivals, easily relax the hard constraint of MPS ≤ 3. Following the principle, to design a ML-PM for departures, it may provide an alternative choice for re-sequencing departure aircraft, consequently, reduces some unnecessary taxi time on the ground. Based on previous researches, the contribution of this paper are the followings:

1. Application of PMS to departure route network design, instead of only use it for arrivals. It could help to realize a feasible, dynamic and efficient departure trajectory management for parallel runway operation. Together with ML-PM arrival route network, the whole network will not only reduce the traffic complexity in the TMA, improve the continuous descending and climbing performance, but also provide a good situation awareness for controllers/pilots. 2. Integration of departures and arrivals with automated conflict-free solution. Taking into account the separation between successive departures and arrivals on the same runway, This system could support a more efficient and realistic trajectory management in the TMA. 3. Improvement of optimization model for more complicated operation situations. Scenario concerning mixed ML-PM with integrated arrival and departure, and scenario concerning segregated ML-PM with only arrivals, will be both studied in details to demonstrate our proposed model performances in terms of capacity, flight efficiency, and sequencing efficiency.

Route Network Design and Mathematical Optimization

The methodology applied to our problem consists of three main components:

1. Route network optimization. It aims to design a novel RNP-based route network to support an efficient control of trajectory, to maximumly use the available airspace for handling a large traffic from or to parallel runways. 2. Optimization module. It is a decision making solver. The designed optimization algorithm searches the decision variables to match the objective requirements, such as conflict-free, less delay etc.. 3. Simulation module. It is split up into several parts illustrating all necessary modelling approaches to compute a more realistic 4D trajectory. It solves the problem of computation about airspeed, aerodynamic forces, motion of aircraft, as well as fuel consumption, flying distance, flying time (or transit time), etc. [START_REF] Liang | Integrated sequencing and merging aircraft to parallel runways with automated conflict resolution and advanced avionics capabilities[END_REF] has introduced simulation module in details. In the following parts, the route network design and the optimization module will be present in details.

Route Network Design

The topology of arrival and departure route networks plays an important role in the commercial aircraft operation. A novel route topology Multi-Level and Multi-Point Merge System (MLMPMS), or ML-PM for short, has been developed for arrival management [START_REF] Liang | Integrated sequencing and merging aircraft to parallel runways with automated conflict resolution and advanced avionics capabilities[END_REF]. It could overcome the constraints of MPS ≤ 3, meanwhile avoid imposing too much workload to the controllers, adapt for two parallel runway operation. For departure management, the application of advanced PM is similar to the arrival ML-PM.

place Fig. 3 about here Firstly, considering the diversity of operation modes for parallel runways, we design two kinds of horizontal topology for departures, see Fig. 3. One is named segregated ML-PM, the other is named mixed ML-PM. The difference between them is whether there is runway re-allocation for aircraft. Mixed ML-PM has it, while segregated ML-PM does not have it. The shadowed areas with two different colors are the possible diverging areas for aircraft. Each colored area indicates the possible manoeuvring airspace for some specific SID route. In conventional topology, departing aircraft generally choose the runway close to its exit way-point of TMA, so as to maximumly reduce the potential conflicts with other traffic. The manoeuvring zones for different runways are well separated with each other. While, in the proposed ML-PM topology, the different colored areas (Diverging zones) overlap together. Aircraft are guided to join different sequencing leg after take-off, then follow the SID to the exit way-point. In a limited airspace, the proposed ML-PM topology may provide more possible manoeuvring area for aircraft instead of that in the conventional topology. For different ML-PM topologies, there are two sequencing legs: the inner leg and the outer leg, with segregated flight altitudes for categories of aircraft. Each sequencing leg has a single operation direction. Aircraft departing from runway either join the inner sequencing leg for East, or join the outer sequencing leg for West. In mixed ML-PM topology, the departing aircraft can take off from different runways, and easily and safely cross each other by holding on the sequencing leg, and then to go different directions. However, in segregated ML-PM topology, aircraft need to take off on the nearest runway to join the corresponding exit point in TMA.

Secondly, because there is an overlapped area in diverging zone of ML-PM topology, then we have to well design the vertical plan to reduce the traffic complexity. The idea is that aircraft with different wake turbulence category fly with different flight altitudes on the sequencing leg. All the available flight altitudes on the inner sequencing leg are higher than that on the outer sequencing leg. Then, there are three parallel flight altitudes for a group of aircraft coming to the same sequencing leg, "Heavy" aircraft will choose the higher level, "Medium" aircraft will use the middle level and "Light" aircraft will enter the lower level, all of the three layers have an unique projection on horizontal level. The reason is that for the commercial aircraft, lower altitude decreases the fuel economy, "Heavy" aircraft need more thrust to counter drag than other categories of aircraft. Taking BCIA as an example, because there is no "Light" aircraft operating at BCIA, thus there is only two segregated flight levels for "Heavy" and "Medium" aircraft. As shown in Fig. 4, a1&a2 are flying on the inner sequencing leg with different higher altitudes than a3&a4 on the outer sequencing leg.

place Fig. 4 about here Thirdly, CCO can reduce fuel consumption and noise compared with the conventional step climb method. However, it is hard to be totally implemented. In an ideal optimal CCO, the aircraft climb from threshold of runway to Top of Climb along an optimized flight path. In reality, aircraft normally maintain a part of level flight. Therefore, we consider to design an economical climb profile, which has the target altitudes at each way-point, to guide aircraft to conveniently complete a near-CCO climb in the ML-PM system. Some reasonable level-off will be accepted for the purpose of well managing the departure aircraft to different direction.

In the end, based on current RNAV SID in reality and the ML-PM based STAR in the previous study, the route network for BCIA South inbound and outbound operation is designed, see Fig. 5. The solid lines are the STAR routes, and the dashed lines are SID routes. The orange area in the North of airport is the merging zone for arrivals, the blue area in the South of airport is the diverging zone for departures. In this ML-PM system, there are different target altitudes on the sequencing legs, for example, in the merging zone for arrivals, 2400 m and 2700 m are two target altitudes on the outer sequencing leg. While, in the diverging zone for departures, 3600 m and 3900 m are two target altitudes on the inner sequencing leg. The higher altitude is for "Heavy" aircraft and the lower altitude is for "Medium" aircraft. In addition, there is target altitude assigned for each way-point. Especially, for way-point W8 and W7, there are two target altitudes, 5400 m is for flows coming from VYK and 4200 m is for flows coming from DOGAR. Remark that controllers in China use meters as unit for flight level.

place Fig. 5 about here

Mathematical Formulation 3.2.1. Assumptions

There is a set of aircraft F arr = {1, 2, ..., n} planned to land at BCIA and another set of aircraft F dep = {1, 2, ..., m} planned to depart from BCIA in a no wind condition (due to lack of wind field data). All the aircraft fly in an International Standard Atmosphere (ISA) environment. Then, the distance between runway 18R-38L and runway 18L-36R is 1960m, the distance between runway 18L-36R and runway 01-19 is 1525m. It is assumed that the independent instrument approaches between runway 18R-36L and runway 01-19, the segregated parallel approaches and/or departures between runway 18L-36R and runway 01-19 are used at BCIA in our study case, see Fig. 6. Runway configuration changes are not discussed in this paper as well.

place Fig. 6 about here For each arrival aircraft i ∈ F arr , the following six kinds of data are given:

1. p e i -the entry point in TMA .

t e

i -the Estimated Time of Arrival (ETA) at the entry point.

v e

i -the initial CAS at the entry point. 4. r e i -the initial runway-in-use.

ET A L

i -the estimated time of landing. 6. cat i -the wake turbulence category.

For each departure aircraft j ∈ F dep , the following six kinds of data are given:

1. p e j -the exiting point in TMA. 2. t e j -the initial departure time. 3. v e j -the initial take-off climbing CAS. 4. r e j -the initial runway-in-use. 5. ET A O j -the estimated time of exiting TMA. 6. cat j -the wake turbulence category.

State space

For each arrival aircraft i ∈ F arr , there are four decision variables, which can be used to control the trajectory of aircraft:

1. t a i -the actual entry time at TMA. 2. v a i -the actual entry CAS at TMA. 3. t T i -the actual turning time on the sequencing legs. 4. r a i -the actual runway-in-use. All these variables are modelled by some discrete variable, considering the operational requirement and commercial aircraft performance. The time t a i is adjusted by a number of slots denoted by j, and the duration of each slot is ∆ = 2s,

t a i = t e i + j∆, j ∈ Z. (1)
The speed v a i is changed in a discrete way as follows,

v a i = v e i (1 + g), g = 0, ±1%, ±2%, ..., ±15%. (2)
The time t T i is controlled by the formulation mentioned below,

t T i = t T imin + h(t T imax -t T imin ), h = 0, 1%, 2%, ..., 100%, (3)
where t T imin is the earliest turning time for aircraft i, and t T imax is the latest turning time. Note that t T i depends on the length of the sequencing leg. As shown in Fig. 7, the aircraft is flying on the outer sequencing leg with a constant CAS. If it turns at the entry point with t T imin , it will have the shortest holding time on the sequencing leg; if it turns at the exit point with t T imax , it will have the longest holding time on the sequencing leg. The aircraft may also choose to turn at specific moment t T i on the sequencing leg, which is defined by the formulation 3. 10 place Fig. 7 about here For the runway assignment, r a i is defined as below:

r a i = 1 if i expected to merge at M1, 0 if i expected to merge at M2. (4)
Here, M1 and M2 are merge points in ML-PM. For BCIA case, M1 is merge point connecting to runway 18R-36L, and M2 is another merge point connecting to runway 01-19.

Similarly, for each departing aircraft j ∈ F dep , there are three decision variables for dynamically modifying the departure trajectory, they are: 1. the actual departure time t d j , 2. the actual turning time on the sequencing t T j , 3. the actual departing runway allocation r d j . The method of modelling the departure decision variables is similar to the arrival decision variables. Note that there is:

r d j =
1 if j to KM, SOSDI, RENOB, 0 otherwise .

(5)

Constraints

Some operational constraints must be carefully considered. They are crucial to the fairness and safety of aircraft. Firstly, for arrivals, t a i must vary in a reasonable range. If aircraft arrive too early before t e i , they need to fly at a higher speed before entering the TMA, which induces a high fuel consumption. If aircraft arrive too late after t e i , they will produce propagation of delay at the destination airport [START_REF] Balakrishnan | Scheduling aircraft landings under constrained position shifting[END_REF][START_REF] Carr | Fast-time study of airline-influenced arrival sequencing and scheduling[END_REF][START_REF] Lee | A study of tradeoffs in scheduling terminal-area operations[END_REF]. In this paper, it is supposed that the earliest time of arrival on the entry point of TMA is limited to 5 minutes before ETA, and the latest time of arrival time is 15 minutes after ETA, then there is:

-5min ≤ δt a i ≤ +15min. (6)
Similarly, for departing aircraft, it exists:

-5min ≤ δt d j ≤ +10min. (7)
Secondly, speed change is limited by the performance of commercial aircraft while in descent profile. In TMA, aircraft airspeed must not be less than the minimum clean configuration speed due to low speed stall, thus v a i is subject to:

v e i (1 -15%) ≤ v a i ≤ v e i (1 + 15%), (8)
v a i ≥
230kt for Heavy aircraft, 220kt for Medium aircraft.

For departing aircraft, there is no speed regulation, aircraft follow its take-off climbing performance. Thirdly, t T i for arrival and t T j for departure on different sequencing legs must comply with geographic constraints:

t T imin ≤ t T i ≤ t T imax , i ∈ F arr (10) t T jmin ≤ t T j ≤ t T jmax , j ∈ F dep (11)
Fourthly, radar separation and wake turbulence separation minimum 2 must be considered to avoid aircraft conflict. According to International Civil Aviation Organization (ICAO) regulations, two aircraft are considered to be in conflict if their horizontal separation is less than the minimum aircraft separation and their vertical separation is less than 1000 ft. If aircraft can not match the vertical separation, then a horizontal separation must be assured. In approach airspace, two kinds of minimum separation have to be considered: wake turbulence minimum separation and approach radar separation. The required aircraft minimum separation in TMA, denoted by s min k,g , k, g ∈ F . Since the time-based system is much more convenient for detecting conflict and metering flows, we consider to transfer the distance-based separation to time-based separation [START_REF] Nikoleris | Autonomous system for air traffic control in terminal airspace[END_REF]. The reference velocity used for this computation is based on the average final approach speed of commercial aircraft by aircraft weight category in no wind condition, which is also the minimum speed of aircraft operating in TMA. It should be noted that for commercial aircraft, they normally follow an airline CAS speed schedule, thus their final approach CAS is relatively steady in normal case. The used reference velocity for computation in our case is True Air Speed (TAS), and it is computed based the corresponding final approach CAS in ISA condition on mean sea level. For other special cases, such as the trailing aircraft was significantly faster than the leading aircraft on final approach, or there is a heavy crosswind condition at airport etc., the time separation must be reconsidered as necessary to ensure that minimum required separation distance was not violated anywhere on final approach. In this paper, the general situation is only considered. Then, the uncertainty in atmospherics is not considered in this paper, because it is assumed that aircraft operate in ISA condition. The uncertainty of pilot execution varies with different culture and individual human behaviours. This value has to be evaluated by human factor experts. In this paper, it is assumed that pilot executes the de-conflict action without delay. If the exact value of uncertainty is given in some case, it can be easily added into current time-based separation. Remark that accounting for these uncertainties will likely require added separation buffers that will limit capacity improvements and may introduce additional delay. In the end, the minimum time-based wake turbulence separations are shown in Tab.1b, it will be used for the conflict computation later.

In the air, link and node conflict computation method is applied. In this paper, all the conflict detections are based on time. Note that the proposed speed regulation in this paper is in terms of CAS, thus there is a need to convert CAS to TAS. The effect of altitude dependencies from CAS to TAS is considered with ISA condition, and the GS (Ground speed) is equal to TAS in this paper, because of zero wind assumption. Then, the link conflict L k is determined based on the time difference at the beginning and the end of arc. Assuming that preceding aircraft k and trailing aircraft g will enter into the same arc u ∈ E defined by two way-points w a and w b , we define ∆t u k,g as the entry time difference and ∆t u k,g as the exit time difference. Consequently, if the following constraints could not be satisfied, then L k increases by one:

∆t u k,g := t w a g -t w a k ≥ s min k,g , (12)
∆t u k,g := t w b g -t w b k ≥ s min k,g , (13)
∆t u k,g • ∆t u k,g ≥ 0. ( 14
)
The node conflict N k is determined based on the time difference on passing a common way-point. Assuming that the preceding aircraft k and the following aircraft g will pass a common way point w c ∈ V, if the following constraint could not be satisfied, then N k increases by one:

∆t w c k,g := t w c g -t w c k ≥ s min k,g . (15)
On the ground, the wake turbulence separation between the departing aircraft and arrival aircraft successively operating on the same runway is considered. In addition, the evacuating runway time of the arrival aircraft is considered as well. Given the maximum evacuating runway at BCIA is 50 seconds, assuming that t L g is the actual landing time of flight g on the runway, and t d k is the actual departure time of flight k, then, the safety constraints are defined as following:

t d k -t L g ≥ s min i, j , k ∈ F arr , g ∈ F dep , (16)
t d g -t d g+1 ≥ s min g,g+1 , g, g + 1 ∈ F dep (17)

Objectives

There are a number of different stakeholders in air transportation system, including Air Navigation Service Provider (ANSP), airlines, airports, and government. Each of them has its own set of objectives. For example, ANSP aims to ensure safety and efficiency of the aircraft. Airlines prefer minimizing the operating fuel costs, maximizing the punctuality of flights, and fairness between airlines. Airports prefer maximizing the throughput. Government prefers minimizing the environmental effects, such as noise and pollution. As a result, there is a trade-off between various and not necessarily aligned objectives in the formulation of our optimization problem. Obviously, our optimization problem is a multi-objective optimization problem. We will apply the linear scalarization (the weighted sum approach) method to convert a multi-objective optimization problem into a single-objective optimization problem.

In our case, we plan to improve the descent and climbing performance by designing a novel route network, then for the objective function, we mainly focus on the four following objectives: conflict-free, less delay, higher throughput, and constrained position shifting. We gather all decision variables for the arrival flight i ∈ F arr to a decision vector u i , for the departure flight j ∈ F dep to a decision vector u j , define the position of aircraft i as p i = (x i , y i , h i ), the position of aircraft j as p j = (x j , y j , h j ), then the objective of our optimization problem is formulated as follows:

z = min C + α 1 D + α 2 S + α 3 P, ( 18
) C = i∈F arr C i ( u i , p i ) + j∈F dep C j ( u j , p j ), (19)
D = 1 n + m ( n i=1 D i + m j=1 D j ), ( 20
)
D i =        t L i ( u i ) -ET A L i 2 if t L i ( u i ) > ET A L i , 0 otherwise. ( 21
) D j =        t O j ( u j ) -ET A O j 2 if t O j ( u j ) > ET A O j , 0 otherwise. ( 22
) S = 1 n + m (t last -t f irst ), (23)
P = n i=1 P i ( u i ), i ∈ F arr ( 24
)
P i = P Actual -P FCFS . ( 25
)
Here, C is the total number of conflicts, D is the average square delay, S is the average operation interval, equals to the total elapsing time divided by number of aircraft in operation, and P is the total position shift for all arrivals. The first term aims at minimizing the conflicts, the second term aims at searching a near-optimal resolution with less deviation from initial planning, the third one aims at speeding up the whole operation process, the fourth aims at sequencing control. In addition, α 1 , α 2 and α 3 are weighting parameters. These parameters fit with the values of the magnitude of D, S , P. In our case, D ∈ [0, 100], S ∈ [0, 50], then, the default setting α 1 = 0.01 and α 2 = 0.02. If arrival sequencing constraint is considered, then α 3 = 0, otherwise, α 3 = 1. Then, n is the number of arrival flights, m is the number of departure aircraft. C i (or C j ) is the conflict of aircraft i (or j), it depends on u i (or u j ). D i (or D j ) is the square delay of aircraft i (or j), it measures the delay propagating to the next phase. t L i is the actual landing time of flight i on the runway, t O j is the actual overflying time of flight j at the exit point of TMA, and t last and t f irst are the actual ending time of the first aircraft and the last aircraft in the instance. P i is the number of position shift between its position in the actual queue and the First Come First Served (FCFS) queue at the runway. P i is only applicable to the segregated parallel approach operation. Because when an aircraft changes its initial landing runway, then its position in the new queue is not the same as in the previous queue.

Optimization Algorithm

In our optimization problem, we have to generate conflict-free trajectories. Taking arrival sequencing and merging optimization problem as an example, firstly, our objective function directly depends on the value of C, D, S , and P, and indirectly depends on the u i and time. According to the state spaces of the four decision variables (t a i , v a i , t T i , r a i ), the number of possible control solutions for each arrival flight N S i = 600 × 30 × 100 × 2. For n flights, at the time t, then the total number of the possible trajectories N T ra j all |t = n × N S i . In order to guarantee a conflict-free solution, the positions of aircraft at every moment (x i , y i , h i , t i ) have be considered, e.i. the computation of C i . If we use a pairwise conflict evaluation method, then for a duration of scenario T scenario , all the positions for an aircraft conflict evaluation N C i all = N T ra j all |t × T scnario . Obviously, the combinatorics associated to our problem is very high and no separation can be directly identified in the objective function. Having several hundreds of aircraft to optimize in a large given time window, we have decided to address such global optimization problem by using a stochastic heuristic approach.

Secondly, our optimization problem needs to consider not only the space-time trajectories for each aircraft, but also the control problem of each aircraft. In this paper, we consider a discrete-time linear system of the form with time t to address this problem,

n i=1 ( p i ) t+1 = n i=1 f ( p i ) t , ( u i ) t , i ∈ F (26)
where n i=1 ( p i ) t ∈ R n is the system state, the positions of all the aircraft. n i=1 ( u i ) t ∈ R m×n is the control input. The state and input must satisfy some constraints, expressed abstractly as

( p i ) t , ( u i ) t ∈ C t i , i ∈ F (27)
where C t i is the constraint set. The instantaneous objective of system depends on both the current state and control action, and is denoted z t n i=1 ( p i ) t , n i=1 ( u i ) t . We judge the quality of control by the objective,

z = T t=0 z t        n i=1 ( p i ) t , n i=1 ( u i ) t        , i ∈ F (28)
where, T is the terminal time. As with the dynamics data, we subscript the constraint set and objective computation function with the time t, to handle the case when they vary with time.

Based on the discrete-time system, we choose to solve a 24 hour traffic optimization problem by application of the RHC technique, which decomposes the original problem into several sub-problems. The RHC policy works as follows, see Fig. 1. Form a predictive model. Compute the performance of system by using the estimate of control input available at time t = T s (k). 2. Optimize. Solve the problem of minimizing the objective, subject to the dynamics and constraints. Here, the objective, dynamics and constraints are estimates, based on information available at time t = T s (k). 3. Execute. Choose the control input to be the value obtained in the optimization problem of step 2.

In step 2. The sub-optimization problem takes the form as following:

minimize

T e (k) t=T s (k) z t         N i=1 ( p i ) t , N i=1 ( u i ) t         , i ∈ N ⊂ F (29) subject to N i=1 ( p i ) t = N i=1 f ( p i ) t-1 , ( u i ) t-1 (30) ( p i ) t-1 , ( u i ) t-1 ∈ C t-1 i (31)
where, N ⊂ F is the number of associated aircraft in current sliding window. Four traffic samples for mixed arrival and departure traffic have been built, each lasting 60 minutes (with 100, 110, 120 and 130 flights respectively). Half of the flights are arrivals, and the other half are departing flights. For arrival flights, there are 14% from KM, 16% from JB, 15% from BOBAK, 25% from VYK, 15% from DOGAR, 15% from GITUM. For departing flights, there are 14% to KM, 16% to SOSDI, 15% to RENOB, 20% to LADIX, 15% to TONIL, 10% to CDY, and 10% YV. There are also five departures-only traffic samples, each lasting 60 minutes with 60, 70, 72, 74, 76 flights respectively.

Automated de-conflict performance

For departure-only scenarios, the results are shown in Fig. 9 and Tab. 2. Our system can solve all the potential conflicts for the different traffic demands with 60, 70, 72, 74, 76 departures respectively, however, it cannot solve all the potential conflict for the demand equal and more than 76 departures. The associated performance of conflict resolution in a sliding window are illustrated by color in Fig. 9. Fig. 9a shows the results from iteration 0 to iteration 2800. There are three window shifts. We zoom in the results in the range of iteration 935 to iteration 1135, see Fig. 9b, to see the process of conflict resolution in the second shift. We found that the number of initial conflict is in proportion with the number of departure flights. The average square delay D and average flying time in the TMA both increase as the volume of departing traffic increases.

place Fig. 9 about here For the mixed arrival and departure scenarios, the results are shown in Fig. 10. There are four different traffic demands with 100 flights, 110 flights, 120 flights, and 130 flights respectively, of which 50% are arrivals and 50% are departures. Our system can successfully resolve all the conflicts, see Fig. 10a. The results from iteration 0 to iteration 2800 are shown, and there are three window shifts. In the first shift, the number of initial conflict is the biggest. We zoom in the results in the range of iteration 0 to iteration 240, see Fig. 10b. Compared with the designed operational capacity of BCIA 88 aircraft per hour in BCIA, our proposed system can handle the demand with up to 120 aircraft per hour at BCIA. place Fig. 10 about here

Climbing and descending performance

For the departures-only case, the climbing performances with 70, 72, 74 and 76 flights are shown in Fig. 11. For all j ∈ F dep , and with different decisions for t d j , t T j and r d j , the vertical profiles of all departures almost keep a similar shape, see Fig. 11a,Fig. 11b,Fig. 11c,Fig. 11d. In each figure, the lines with violet color refer to the flights with wake turbulence category "Medium", the lines with green colour refer to the flights with wake turbulence category "Heavy". The shape of the departure is structured. Fig. 11 shows that the holding time on the sequencing legs increases with the number of departure flights, because of conflict resolution. The maximum level-off time on the sequencing leg is less than 6 minutes.

place Fig. 11 about here For mixed arrival and departure traffic, more attentions should be paid to the hot spot areas, see Fig. 5, because some of the arrival and departure trajectories cross inside these hotspot areas. The route design should not gain a capacity increase in compromise of safety. Even though, current ML-PM system could propose an automated conflict-free solution for managing arrivals and departures, the system crash situation must be considered. Thus, it is important that the departure routes and the arrival routes are well separated in the hot spot areas, so that controllers do not give clearance unless separation is assured.

Taking the traffic sample with 130 flights as an example, Fig. 12a and Fig. 12b show the integrated departures and arrivals trajectories. The green trajectories represent departures, the blue trajectories represent arrivals. We then zoom in the traffic around the hot spot areas. Fig. 12c shows that arrivals are all above the departures in hot spot area 1. Fig. 12d shows that departures going to VM may have intersection with the arrivals coming from DOGAR in spot area 2. Fig. 12e shows that the arrivals are all above departures in hot spot 3. Consequently, there is sufficient separations between departure and arrival aircraft in the hot spot area 1 and area 2, however, we have to modify the arrival procedure from DOGAR for safety. We modify the first segment of the route for aircraft from DOGAR. These aircraft will join the way-point W7 directly instead of passing way-point W8. Fig. 12f shows the result of this modification. We found that there is sufficient horizontal separation between departures and arrivals now, because the departures going to VM will always pass over the arrivals from DOGAR.

place Fig. 12 about here 4.1.4. Throughput We compared our system performance regarding throughput with the maximum hourly movements observed at BCIA on the 07/09/2015 (baseline in Fig. 1). In the baseline, the three runways can be used simultaneously for both departures and arrivals, while in our system, only runway 18R-36L can be used for both departures and arrivals; runway 18L-36R is only for departures, and 01-19 is only used for arrivals. In the baseline, there is a maximum throughput of 72 movements per hour for arrivals, 74 movements per hour for departures, and 95 movements per hour for total. The declared capacity at BCIA for total movements with departure plus arrivals is 88 aircraft per hour. Remark that throughput is the airports actual air traffic movement on the day of operations. Capacity is the theoretical air traffic movement capability of an airport [START_REF] Eurocontrol | Arrival Manager, Implementation GUIDELINES and Lessons Learned[END_REF].

In order to get different maximum throughputs of our system, in addition to the four mixed-traffic and five departures-only traffic scenarios, four arrival-only traffic scenarios were added. These included: four 90 minutes scenarios with 66, 77, 88, and 99 flights respectively. By applying our system under the tested conditions, we can reach a maximum throughput of 74 movements per hour with only departures, 99 movements per hour with only arrivals, and 120 movements per hour with integrated arrivals and departures. This corresponds to a 30.6% increase in maximum landing movements per hour, and a 26.3% increase in maximum total movements per hour, see Fig. 13. Note that for the departure only case, there is no significant improvement. Because there are only three decision variables for departure case, compared with four decision variables in the arrival case. The speed regulation in departure case is not available. In addition, it is remarked that the operational capacity generally depends on the workload of controllers, the human-in-loop validation work need to be planned, so as to evaluate the real benefit of capacity increase by our proposed system. place Fig. 13 about here 4.1.5. Traffic demand and efficiency Let us further look at the relationship between traffic demand and efficiency. We consider three different arrival traffic samples: 100, 110, and 120 per hour. For each sample, we performed 10 runs with different available lengths of sequencing leg respectively. Here, the available length of sequencing leg is defined as a percentage of the total designed sequencing leg. The results are illustrated in Fig. 14a and Fig. 14b.

In Fig. 14a, the X-axis is the percentage of total designed sequencing leg, and the Y-axis is the value of objective function. For each traffic sample, we first plot the associated values on the graph, and then we build the trend line using a second-order polynomial function. Note that, all the objective values are less than 1. This means that all the potential conflicts were successfully solved. However, the efficiency of trajectory planning with different available sequencing legs vary. This efficiency is evaluated using two factors: the average square delay and the average flying time in the TMA. Since the value of C is always 0, thus the objective function mainly represents the average square delay, see Fig. 14a. We can see that with different available percentages of the total designed sequencing leg, the square delays are different. The objective value for each kind of traffic sample has both a min and max value. We can see clearly that as the traffic demand increases from 100 flights to 120 flights per hour, the gradient of polynomial trend line (the dashed line) also increases. In Fig. 14b, the average flying time in the TMA slightly increases when the percentage of the total designed sequencing leg increases. However, it does not increase with the number of traffic demand. For 100 and 110 flights per hour, the minimum objective value is reached within 0% of the total designed sequencing leg; for 120 flights per hour, it is within 20%. This means that the minimum objective value with the average square delay may not correspond to the minimum available sequencing leg.

place Fig. 14 about here

Computation time

We tested the computation time of our system with five different traffic scenarios. For each type of traffic, we ran 10 tests and collected the Max, Min, and Mean running time. The associated results are shown in Tab.3. We found that the running time increased with the density of traffic demand. For example, for a 2 hour scenario, with W = 3600s, and S = 1800s, it took an average of 508 seconds to find a conflict-free solution for heavy traffic, 362 seconds for medium traffic, and only 220 seconds for light traffic. The computation time also increased with the scenario time period. With 100 movements per hour, it took approximately 362 seconds for optimizing a 2 hour scenario, and 1137 seconds for a 6 hour scenario and 2440 seconds for a 12 hour scenario. To simulate the real situation, we can change the parameter S , W and t FINAL to realize a non-stop rolling process. But, the computation time may increase significantly when using real wind, weather and non ISA data.

Segregated ML-PM with only arrivals

In order to ensure continuous traffic demand at runways and maximize runway usage, a minimum level of queuing is required to keep pressure on the runways. However, additional time in holding is detrimental to operational efficiency, fuel consumption and the environment. Therefore, there exists a trade-off between approach efficiency and runway throughput. In this part, we continue to investigate the approach efficiency with different sequencing techniques. Two different sequencing techniques are considered:

1. minimizing position shift P ON. It tries to maximumly maintain the FCFS queue. It is a constrained sequencing, α 3 = 1. 2. minimizing position shift P OFF. There is no constraints on the sequencing. It is a relaxed sequencing, α 3 = 0.

We will study the difference between the two techniques in terms of sequencing and merging efficiency, and vertical flight efficiency.

Data preparation

Four traffic samples are built. Each of them lasts 90 minutes (with 66, 77, 88, and 99 flights respectively). The samples are generated according to the geographical distribution of flights at BCIA: 14% from KM, 16% from JB, 15% from BOBAK, 25% from VYK, 15% from DOGAR, 15% from GITUM.

Sequencing and merging efficiency

The Unimpeded transit time and the additional transit time were introduced. They are developed by the Performance Review Unit of EUROCONTROL to characterize the performance of arrival management process. The unimpeded transit time is the transit time in the area without congestion. The additional transit time is the difference between the actual transit time and the unimpeded transit time. It represents the extra time generated by the arrival management and is a reference for the level of inefficiency (holding and sequencing) of the inbound traffic flow during times when the airport is congested. For simplicity, here we define the unimpeded transit time ∆ ini i as the flying time of aircraft between the entry-point to the TMA and the landing runway with initial states (t e i , v e i , r e i ), following the defined STAR route without level-off on the sequencing leg. We define the actual transit time ∆ Act i as the flying time of aircraft between the entry-point in the TMA and the landing runway with actual status (t a i , v a i , r a i , t T i ). The additional transit time ε in i is the difference between ∆ Act i and ∆ ini i . It is the additional transit time in the TMA airspace, which represents the cost of sequencing, merging and spacing an aircraft i in TMA. We define ε out i as the cost of sequencing and spacing on aircraft i at the entry point of the TMA. It represents the cost of delay to the adjacent sector due to our de-conflict strategy "entry slot time change". The total cost ε i is the sum of ε in i and ε out i . It represents the total effect of sequencing, merging and spacing in the TMA. The relationship between these variables is illustrated in Fig. 15. We have:

∆ ini i = ETA L i -t e i , (32)
∆ Act i = t L i -t a i , (33)
ε in i = ∆ Act i -∆ ini i , (34)
ε out i = t a i -t e i if t a i > t e i , 0 otherwise. ( 35
)
ε i = ε in i + ε out i , (36)
ε average = 1 n n i=1 ε i i ∈ F . ( 37
)
place Fig. 15 about here place Fig. 16 about here Fig. 16 shows the comparison of different performance indicators with minimizing position shift P either ON or OFF. We found that when the number of flights in this scenario reached 99, the total number of unresolved conflicts with the P-ON setting, which allows no more than three position shifts, can not be reduced to zero, however, with the P-OFF setting, which allows more than three position shifts, the total number of unsolved conflicts can be reduced to zero. This means that more dynamic position shifting helps to resolve conflicts. Fig. 16f shows that the MPS increases with the number of flights. With the P-ON, when there are less than 88 flights, the MPS is less than 3. However, when there are 99 flights, it is hard to maintain the FCFS queue, meanwhile it is hard to resolve all the conflicts. As a result, more dynamic position shifts are required.

The average square delay increases with the number of flights. However, the relaxed sequencing performs better than the constrained sequencing. Fig. 16c shows that the relaxed sequencing also performs slightly better than the constrained sequencing in term of average flying time. However, in terms of average fuel consumption, there is no obvious difference between them. ε average is the cost of sequencing, merging and spacing in the TMA. It is found that ε average increases with the number of flights. The relaxed sequencing performed better than the constrained sequencing.

We also found that the scenario with 66 flights is a non-congested condition, and it is easy to solve with our segregated ML-PM system. And the scenario with 99 flights is a super-dense condition, and Runway 01-19 is the first runway to be saturated.

Vertical flight efficiency

Here, we only consider the level-off on the sequencing leg in the descent. As shown in Fig. 17, the average times flown level per flight for the different scenarios are almost the same in the relaxed sequencing situation, about 90 seconds, 1.5 minutes, about 9.0% of total flying time in the TMA. In the constrained sequencing situation, it slightly increases with the number of flights, but is never no more than 12% of the total flown time in the TMA.

place Fig. 17 about here

Conclusion

Integration of aircraft arrivals and departures in a busy TMA is a problem in a complex network. An efficient and feasible sequencing and managing of arrivals and departures involves not only the strategic decision making but also tactical control. In this paper, a novel system has been developed to devise an integrated aircraft trajectory planning concept for parallel runway operation. A ML-PM route network is designed to support more efficiently arrival merging and departing, and it is built up based on the current operational route network. A systematic approach is proposed to dynamically managing the large trajectories in routine dense operation with automated conflict-free solutions.

According to the numerical results, we can conclude that our system shows good performances. It can automatically solve all the potential conflicts in departure-only cases and mixed departure/arrival cases. It can realize a shorter flying time and a near-CDA descent for arrival aircraft, an economical climbing for departure aircraft, an easier runway allocation together with trajectory control solution. It shows a good sequencing efficiency in TMA as well. Compared with baseline, in mixed ML-PM mode, our proposed system can increase throughput around 26% under tested conditions with ISA environment. Although the methodology defined here is illustrated using the BCIA airport, it could be easily applied to airports worldwide.

However, due to lack of wind data, we did not consider the effect of wind in this paper, but wind effect mathematical model has been integrated in the simulation model. In future work, historical wind data will be collected and then used in the simulations. In addition, weather will be taken into account for improving the robustness of our model.

Figure 1 :

 1 Figure 1: Hourly movements at BCIA on 07/09/2015

Figure 3 :

 3 Figure 3: Advanced Point Merge network for Departures

Figure 5 :

 5 Figure 5: Advanced Point Merge network for BCIA

  8. At starting time T INIT , we consider a time interval, named sliding window, extending k steps into the ending time T FINAL . Four parameters are introduced: 1. W: the time length of the sliding window; 2. S: the time shift of the sliding window at each iteration; 3. T s (k): the starting time of the k th sliding window, T s (k) = T INIT + kS; 4. T e (k): the ending time of the k th sliding window, T e (k) = T INIT + kS + W. place Fig.8 about here We then carry out the following steps:

Figure 8 :Figure 11 :

 811 Figure 8: Sliding window in RHC technique

  Figure 12: Hot spots areas with integrated arrivals and departures

Figure 13 :

 13 Figure 13: Hourly movements comparison between baseline and our system

Figure 15 :

 15 Figure 15: Definition of additional time

Table 2 :

 2 Performance of ML-PM based departure

	Traffic volume C D [min2] S [min] Average flying time [min]
	60	0 4.70	0.99	15
	70	0 27.58	0.97	15
	72	0 44.67	0.95	16
	74	0 77.79	1.03	16
	76	2 108.21	1.10	16
	4. Experimental Results
	4.1. Mixed ML-PM with integrated arrivals and departures
	4.1.1. Data preparation

Table 3 :

 3 Computation times (Mixed ML-PM with integrated arrivals and departures)

	Running times

Near-parallel runways are non-intersecting runways whose extended centre lines have an angle of convergence or divergence of 15 degrees or less.

ICAO mandates separation minima based upon wake turbulence categories. These minima are typically categorised as follows: Light (L) Maximum Take Off Weight (MTOW) of 7000 kilograms or less; Medium (M) MTOW of greater than 7000 kilograms, but less than 136000 kilograms; Heavy (H) MTOW of 136000 kilograms or greater; Super (J) Refers only to the Airbus A380.

Acknowledgement

The authors would like to thank Serge Roux for the collection and the treatment of ADS-B flight data, and thank the anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions.