Dr Sameer Alam

Professor Frédéric Barbaresco

Professor Sonia Cafieri

Professor Fedja Netjasov

Valentin Polishchuk

Catya Zuniga

Serge Roux

François Huchet

Mohammed Sbihi

Gilles Baroin

Andrija Vidosavljevic

Arianit Islami

Clément Bouttier

Imen Dhief

Jun Zhou

Ji Ma

Ning Wang

Olga Rodionova

Romaric Breil

Marina Sergeeva

Tambet Treimuth

Vincent Courjault-Radé

Yi-Ping Jiang

Sabine Wen-Juan Wang

Antoin

Charles

Damian Pascal

Tihoni

Keywords: Air traffic management, Trajectory based operation, Merging and sequencing, Multiple parallel runway

Congestion in Terminal Manoeuvring Area (TMA) at hub airports is the main problem in Chinese air transportation system. At most of the hub airports, the capacity is near saturated or even overloaded. Civil Aviation Administration of China (CAAC) reported that Air Traffic Management (ATM) is the main cause of delays. Despite the already overloaded ATM system, the Chinese airplanes fleet is continuing to expand. China will become the largest traffic flow in the world before the end of 2035. There is an urgent need to develop a new more efficient method for sequencing and merging arrival flows in TMA, so that airports can maximise the benefits from the emerging Communication, Navigation and Surveillance (CNS) techniques, and consequently increasing capacity.

Automation can be highly efficient in ATM, however, due to safety considerations, full automation in ATM is still a challenge. Facing extremely dense operations in complex TMA, we can consider reducing traffic complexity by solving all potential conflicts in advance with a feasible trajectory control for controllers, or automating a large proportion of routine operations, such as sequencing, merging and spacing. As parallel runways are a common structure of Chinese hub airports, in this thesis, we propose a novel system to integrated sequencing and merging aircraft to parallel runways. Our methodology integrates a Area Navigation (RNAV)-based 3D Multi-Level and Multi-Point Merge System (MLMPMS), a hybrid heuristic optimization algorithm and a simulation module to find good, systematic, operationally-acceptable solutions. First, a Receding Horizon Control (RHC) technique is applied to divide 24-hour traffic optimization problem into several sub-problems. Then, in each sub-problem, a tailored Simulated Annealing (SA) algorithm and a trajectory generation module worn together to find a near-optimal solution. Our primary objective is to rapidly generate conflict-free and economical trajectories with easy, flexible and feasible control methods. Based on an initial solution, we continuously explore possible good solutions with less delay and shorter landing interval on runway.

Taking Beijing Capital International Airport (BCIA) as a case to study, numerical results show that our optimization system performs well. First, it has very stable de-conflict performance to handle continuously dense traffic flows. Compared with Hill Climbing (HC), the tailored SA algorithm can always guarantee a conflict-free solution not only for the mixed iii or segregated parallel approach (arrivals only) pattern, but also for the independent parallel operation (integrated departures and arrivals) pattern. Second, with its unique Multi-Level Point Merge (ML-PM) route network, it can provide a good trajectory control solution to efficiently and economically handle different kinds of arrival flows. It can realize a shorter flying time and a near-Continuous Descent Approach (CDA) descent for arrival aircraft, compared with baseline. For normal traffic, with near-equal traffic demand for two runways, with landing balancing function ON, the average flying time on different Standard Terminal ARrival (STAR) routes can be reduced up to 8 minutes compared with the baseline. It also realizes an easier re-sequencing of aircraft with more relaxed position shifting as well, compared with conventional sequencing method. Theoretically, the Maximum Position Shifting (MPS) can be up to 6 positions, overcoming the hard constraint of 3 position shifts (MPS <= 3). For asymmetric traffic, with big difference on traffic demand for two runways, with runway balancing function ON, it is more likely to find a conflict-free solution compared with the runway balancing function OFF, and again reduces the average flying time. Third, it is efficient for the segregated parallel approach patterns. Compared with hard constrained position shifting, which is often used in current Arrival Manager (AMAN) system and controller's manual-control First Come First Served (FCFS) method, it can reduce the average delay, average additional transit time in super dense arrival situations. The average time flown level per flight is less than 12% of total transit time in TMA. Fourth, in independent parallel patterns, it can provide a range of useful information concerning the associated objective value, the average flying time, crossing trajectories in hot spots between arrivals and departures, the efficiency of using different designed sequencing legs in ML-PM route network. Thus, it helps the Air Navigation Service Provider (ANSP) to find the best configuration of ML-PM route network to efficiently satisfy the traffic demand. Last but not least, the computation time of our system is reasonable. It generally needs around 290s-350s for 2 hours of heavy traffic demand with the mixed parallel approach.

In conclusion, theoretically, our system realizes good trajectory planning of dense flows at busy airports. It can guarantee a conflict-free solution, increase runway throughput, and minimize delay. At the same time, it can simplify merging, re-sequencing, and improve the economical descent profile with advanced ML-PM route design. Although the methodology defined here is illustrated using the BCIA airport, it could be easily applied to airports worldwide.

Résumé

La congestion dans les Terminal Manoeuvring Area (TMA) des aéroports en hub est le principal problème dans le transport aérien chinois. La plupart des aéroports en hub sont surchargés voire presque saturés. La Civil Aviation Administration of China (CAAC) a signalé que la mauvaise performance de l'Air Traffic Management (ATM) est la principale cause de retard. Bien que le système ATM soit déjà surchargé, la flotte d'avions chinoise continue de se développer. La Chine deviendra le flux de trafic le plus important dans le monde avant la fin 2035. Face à la forte demande dans un proche avenir, pour bénéficier des nouvelles techniques de Communication, Navigation and Surveillance (CNS), il est urgent de développer de nouvelles méthodes pour un séquencement efficace et une bonne fusion des flux d'arrivée dans les TMA, afin d'augmenter la capacité de l'aéroport. L'automatisation permet d'atteindre des performances élevées, mais en raison des considérations de sécurité, l'automatisation complète de l'ATM est encore un défi. Face au trafic extrêmement dense dans les TMA, nous pouvons envisager d'automatiser une grande partie des opérations de routine, comprenant la planification, le séquencement et la séparation. Nous proposons dans cette thèse un nouveau système automatisé de séquencement des avions et de fusion des flux vers des pistes parallèles, qui sont utilisées dans la plupart des aéroports chinois. Notre méthodologie intègre un réseau de route 3D nommé Multi-Level Point Merge (ML-PM) basé sur le concept de l'Area Navigation (RNAV) et un algorithme d'optimisation heuristique hybride pour trouver une solution correcte, opérationnellement acceptable. Un algorithme de Simulated Annealing (SA) spécifique et un module de génération de trajectoire collaborent pour rechercher la solution quasi optimale. Notre objectif est de générer en temps réel des trajectoires sans conflit, minimisant la consommation de carburant et permettant des méthodes de contrôle faciles et flexibles. Dans ce but, nous explorons en permanence les solutions avec le moins de retard et assuront l'atterrissage le plus rapide. Nous déterminons quatre variables de décision pour contrôler chaque vol : l'heure et la vitesse d'entrée dans la TMA, le temps de vol sur l'arc de séquencement et le choix de la piste utilisée. La simulation de trajectoire dans les différentes phases de vol est basée sur le modèle de performances BADA.

Dans le cas de l'aéroport de Beijing Capital International Airport (BCIA), les résultats v numériques montrent que notre système d'optimisation de bonnes performances sur le séquencement et la fusion des trajectoires. Tout d'abord, il permet d'assurer des performances de résolution conflit très stables pour gérer les flux de trafic continuellement denses. Par rapport à l'algorithme Hill Climbing (HC), le SA peut toujours trouver une solution sans conflit, non seulement pour l'approche parallèle mixte ou séparée (pour les arrivées), mais aussi pour les configurations parallèles indépendantes (départs et arrivées intégrés). Ensuite, avec un réseau d'itinéraires ML-PM unique, il peut fournir une bonne solution de contrôle de la trajectoire pour traiter efficacement et économiquement différents types de flux d'arrivée. Il peut réaliser un temps de vol plus court et une descente vers le bas en Continuous Descent Approach (CDA) pour l'avion d'arrivée. Pour un trafic normal bien réparti sur les deux pistes, avec la fonction d'équilibrage d'atterissage ON, le temps de vol moyen sur les différentes routes Standard Terminal ARrival (STAR) peut présenter jusqu'à 8 minutes de moins que la référence. Il peut réaliser un re-séquencement plus facile des avions avec un déplacement de position plus relâché. Théoriquement, les Maximum Position Shifting (MPS) peuvent atteindre 6 positions, surpassant la contrainte difficile de 3 positions (MPS <= 3). Pour une demande de trafic très asymétrique sur les deux pistes, avec la fonction d'équilibrage des pistes ON, il est de plus possible de trouver la solution sans conflit qui permet d'équilibrer les pistes avec un temps de vol moyen moins élevé, par rapport à la fonction d'équilibrage des pistes OFF. Troisièmement, l'algorithme montre son efficacité dans un modèle d'approche parallèle séparé avec une capacité de séquencement plus relâché. Par rapport au décalage de position forcé dur, qui est souvent utilisé dans le système actuel Arrival Manager (AMAN) et la méthode First Come First Served (FCFS) utilisé par les contrôleurs, il peut réduire le délai et le temps de transit moyens dans une situation d'arrivée très dense. Le palier par vol est inférieur à 12% du temps de transit total dans la TMA. Quatrièmement, en configuration parallèle indépendant, il peut fournir des informations différentes concernant la valeur objectif associée, le temps de vol moyen, les trajectoires de croisement en point chaud entre les arrivées et les départs, l'efficacité avec différents arcs de séquencement conçus dans le réseau de route ML-PM etc.. Ainsi, il aide les Air Navigation Service Provider (ANSP) à trouver la meilleure configuration du réseau de route ML-PM pour satisfaire efficacement la demande de trafic. Enfin, le temps de calcul de notre système est raisonnable. Il nécessite généralement environ 290s-350s pour 2 heures de forte demande de trafic avec une approche parallèle mixte.

En conclusion, théoriquement, notre système réalise une bonne planification de trajectoire en flux denses. Il peut assurer une solution sans conflit, augmenter le débit, minimiser le délai. Dans le même temps, il peut simplifier la fusion, le re-séquencement, l'amélioration du profil de descente économique associé avec la conception de l'itinéraire avancé ML-PM. Bien que la méthodologie définie ici soit appliquée à l'aéroport BCIA, elle pourrait également être vi appliquée à d'autres aéroports dans le monde.

Mots

Motivation

Due to rapid economic growth, the demand for air services in China has significantly increased in the last 10 years. The number of aircraft movements grew at an average rate of 9.9% per year between 2006 and 2015. However, the on-time performance of flights in this period dropped from 81.48% to 68.33% (CAAC, 2016), see Fig. 1.1. Flight delays have several negative impacts. Firstly, from an economic point of view, given the uncertainty of delay occurrence, passengers tend to plan to travel many hours before their appointments to ensure they arrive on time, thus increasing their trip costs. Airlines incur penalties, fines and additional operational costs, such as crew and aircraft retentions in airports [START_REF] Ball U. A ; Ball | Total delay impact study[END_REF]. Secondly, from a sustainability point of view, delays may also cause environmental damage by increasing fuel consumption and gas emissions (Ryerson u. a., 2014). Thirdly, frequent and long delays generate passenger strong discomfort, which may lead to bad behaviour towards airlines and airport staff, threatening air transportation safety. CAAC reported that ATM is now the most important cause of delay, producing 30.68% of flight delays in 2015, ahead of all the other factors like weather, airport, passenger and airlines etc. The key cause of delays in ATM is the disequilibrium between the high traffic flow demand and the low capacity. Only around 30% of the airspace is reserved for civil aviation operations, compared to about 80% of the airspace in the US. In the en-route part of the civil airspace, the average traffic flow in the 13 busiest airways in China was more than 500 flights per day per route in year 2014. At hub airports, capacity is close to saturation or even overloaded. For example, the seven busiest airports in China exceeded their capacities by more than 30% during their peak traffic periods. Average landing delay at the ten busiest airports is around 26 minutes, and average flight holding time on-ramp, from closing doors to pushing back, is around 12 minutes [START_REF] Caac | [END_REF]. With the successful implementation of some projects, for example the Reduced Vertical Separation Minimum (RVSM) between 8900 meters (29100 feet) and 12500 meters (41100 feet), the en-route traffic capacity has been significantly enhanced. Nowadays, the air traffic bottleneck is shifting from en-route segments to terminal airspace around the busy airports.

Despite the already overloaded ATM system, the Chinese fleet is continuing to expand. Boeing (2015) forecasts that over the next 20 years, China's commercial airplane fleet will nearly triple: from 2570 airplanes in 2014 to 7210 airplanes in 2034. Airbus (2016) forecasts that domestic China will become the largest traffic volume before the end of 2035, supplanting domestic US. However, the airspace for civil aviation operation in China is still expected to be very limited due to state security concerns. Parallel runways are the main structure of Chinese hub airports. More and more parallel runways are built in existing or new hub airports, such as Beijing DaXing international airport which will have 5 parallel runways. Facing high demand in the near future, current conventional control methods are no longer suitable. There is an urgent need to develop a novel approach to efficiently manage arrival flows to parallel runways, to significantly increase capacity.

In the United States, the Super Density Operations (SDO) project has been proposed as a part of the NextGen research program. Its purpose is to enable significantly increased and robust throughput at the most congested metropolis airports while minimizing environmental impact. The concept provides for a transition from current operations to a terminal system that relies on automation for a large portion of routine operations (including scheduling, sequencing, spacing) and suggests leveraging the complex problem solving abilities of humans to manage recovery from off-normal events [START_REF] Isaacson U. A | [END_REF]. In Europe, under the concept of Trajectory Based Operation (TBO) in the SESAR program, a planning tool named Medium Term Conflict Detection and Resolution (MT-CD&R) systems has been designed to help controllers manage the synchronization of 4D arrival trajectories at a tactical level [START_REF] Ruiz U | [END_REF]. Reducing traffic complexity by planning conflict-free trajectories in advance is an important component of increasing capacity and reducing delay. The system developed in this thesis will be applied to BCIA, aiming to alleviate the heavy delay in China.

Background

An Air Navigation Service Provider (ANSP) is a public or a private legal entity providing Air Navigation Service (ANS). Depending on the specific mandate, an ANSP provides one or more of the following services to airspace users:

1. Air Traffic Management (ATM), 2. Communication, Navigation and Surveillance (CNS), 3. Meteorological service, 4. Search and rescue, 5. Aeronautical information services.

ATM is the primary service provided by ANSs. ATM is dynamic, integrated management of air traffic and airspace, safely, economically and efficiently, through the provision of facilities and seamless services in collaboration with all parties and involving airborne and groundbased functions (ICAO Doc4444, 2007). It includes Air Traffic Services (ATS), Airspace Management (ASM), and Air Traffic Flow Management (ATFM), see Fig. 1.2. , 2011). Airspace organization will establish airspace structures in order to accommodate the different types of air activity, volume of traffic and differing levels of service. ICAO Doc.9713 also mentions that ASM is a planning function with the primary objective of maximizing the utilization of available airspace by dynamic time-sharing and, at times, the segregation of airspace among various categories of users based on short-term needs. 3. ATFM is a service established with the objective of contributing to a safe, orderly and expeditious flow of air traffic by ensuring that ATC capacity is utilized to the maximum extent possible, and that the traffic volume is compatible with the capacities declared by the appropriate ANSP. ATFM provided by ANSP is to balance air traffic demand with a ATM system capacity. ATFM activities can be divided into three phases: strategic phase, pre-tactical phase and tactical phase (detailed information is described in A.2). The main mission of a flow manager or network manager is to optimize traffic flows according to air traffic control capacity, while enabling airlines to make correct decisions that ensure safe and efficient flights. ATFM does not solve all potential conflicts, however, it can provide strategies to avoid the occurrence of long delays, such as the ground holding procedure for aircraft. The relationships among flight phases, ANSP, ATS, and controlled airspace are illustrated in Tab. 1.1. In this thesis, we will focus on the terminal airspace around the airport, which is named TMA, see Fig. 1.3. TMA is a control area usually established at the confluence of ATS routes in the vicinity of one or more major aerodromes. Inside TMA, most aircraft are in the flight phases of departure, arrival or approach. Thus, the associated traffic complexity is very high. Our contribution targets optimizing the current traffic management system in TMA, and then at increasing the capacity in TMA.

Contributions

How to increase capacity, and then to reduce the serious delay in main airports has been an important research issue for a long time. Runways are the crucial resource of the airport. To increase runway throughput, we have to optimize landing and take-off. Because landing usually takes priority over taking off, there is a significant body of work focused on the airport landing problem.

Significant advancements in optimization models and algorithms have been achieved in a range of disciplines. For example, automated arrival scheduling and sequencing technologies have been heavily investigated in the past 30 years. [START_REF] Erzberger Und Nedell ; Erzberger | Design of automated system for management of arrival traffic[END_REF] ;[START_REF][END_REF]; Erzberger (1995) developed an automated sequencing and merging system based on a hierarchy of advisory tools for controllers, referred to as the Center-TRACON Automation System (CTAS), which was installed and evaluated at Denver and Dallas/Ft.Worth air traffic control facilities [START_REF] Denery | [END_REF]. Later, the AMAN system was devel-oped in Europe and US as a decision tool to support approach controllers to handle the tasks of aircraft sequencing and metering (Fairclough, 1999). Nowadays, AMAN plays an important role in dense airport operation management, especially when there is bad weather or a runway closure (Hasevoets und Conroy, 2010). However, the AMAN system can not provide a conflict detection and resolution strategy for controllers, so in dense TMA operations, there is still a lot of pressure on controllers to ensure safety. In CTAS, controllers are only capable of using some automation tools. Therefore, full automation of CTAS is still a challenge.

New techniques and procedures in CNS have greatly improved the accuracy of aircraft positioning and trajectory prediction, making the modernization of the ATM system more feasible. More precisely, the application of RNAV2 techniques in all flight phases has directly contributed to improve airspace optimization and flight efficiency. The Automatic Dependent Surveillance Broadcast (ADS-B)3 system benefits both pilots and controllers. It improves both flight safety and flight efficiency. Recently, a number of researchers have considered the benefits of combination of advanced avionics capabilities with ATM. [START_REF] Erzberger | [END_REF]; [START_REF] Erzberger | Transforming the NAS: The next generation air traffic control system[END_REF]; [START_REF] Ky Und Miaillier ; Ky | SESAR: towards the new generation of air traffic management systems in Europe[END_REF] stated that TBO with sharing of the trajectory information via ADS-B between air and ground and between air and air, will enable a safer and more efficient handling of flights in the next generation of air traffic control system. Prevot u. a. ( 2005) presented a concept of cooperative air traffic management to address the integration of trajectory based operations and Airborne Separation Assistance Systems (ASAS). [START_REF] Oberheid Und Söffker ; Oberheid | Cooperative arrival management in air traffic control-a coloured Petri net model of sequence planning[END_REF] demonstrated that cooperation between air and ground enables the arrival planning system to establish an optimized sequence to the runway. [START_REF][END_REF] analyzed the effects of automated arrival management, airborne spacing, controller tools, and data link, while Haraldsdottir u. a. ( 2009) presented arrival management architecture and performance analysis with advanced automation and avionics capabilities. Erzberger u. a. (2010) described the concept of automated conflict resolution, arrival management, and weather avoidance for air traffic management.

However, several challenges remain before we have a good, systematic, operationallyacceptable solution for efficient managing of complex flows in TMA. Potts u. a. ( 2009) reviewed the research studies on the Airport Landing Problem (ALP) and the Airport Take-Off Problem (ATP), and clearly stated that we have to pay attention to the following issues: 2. A quick and good (near-optimal) solution, instead of a slow and optimal solution, is more interesting in the view of a controller. 3. The control problem should be considered, not just the decision problem. 4. More precise information in advance will help the controllers, airports and airlines to do Collaborative Decision Making (CDM). 5. Developing an integrated model including runway assignment, scheduling take-off and landing together, and gate assignment, etc. is a major challenge for airport operation optimization.

Here, we contribute to resolve the following challenges. The first challenge is to integrate aircraft sequencing, merging and runway assignment together to support real-time traffic operations. The second challenge is to make an automated, cooperative, and economic control of trajectories. Here, automated control means our system can automatically provide a feasible solution to an aircraft according to the prescribed objective. Cooperative control means that different stakeholders can make a collaborated decision for the management of trajectories. Economic control means less fuel consumption, and less CO 2 emissions can be realized. The last challenge is to overcome the constraint of "MPS less than three" (explained later, see section 2.1.1), which has limited the optimization of the Aircraft Sequencing Problem (ASP) for a long time. A more dynamic re-sequencing proposed by our system can gain more benefits on capacity without increasing the workload of controllers.

The three key goals of our research are:

1. To propose a novel system to integrated sequencing and merging aircraft to parallel runways in terms of multi-objectives: under hard constraints on collision avoidance, soft constraints on the deviation from ETA, optimized fuel consumption with continuous descent approach, and balanced runway landing rate. The system should quickly generate a conflict-free, least-delay, and operationally-acceptable good solution to efficiently schedule dense traffic in busy TMA. This system should reduce the traffic complexity in TMA, which in turn should greatly increase capacity. 2. To design an advanced RNP-based 3D ML-PM route network to support our concept.

Using this network, the proposed system should realize efficient and economical trajectory control. More precisely, first, we should realize a good decision and control support for runway allocation, balancing the asymmetric traffic on different runways. Second, we should realize an easy and more relaxed re-sequencing for arrival flights. Third, we should realize an economical descent profile with less level-off flying time during the approach phase. 3. To derive a robust and optimal solution both at the pre-tactical and tactical levels in a dynamical fashion. The proposed system should well solve not only the strategic and pre-tactical decision problem of sequencing, but also the tactical control problem of merging. In addition, it should provide rich information concerning different operational patterns on parallel runways. Such information includes: conflict resolution, average square delay, re-sequencing, average additional transit time in TMA, and the relationship between capacity and efficiency, which will help different actors (ANSPs, airlines and airports) to make good decisions.

Outline of the dissertation

The dissertation is organized as follows: First, state of the art on air traffic management in TMA is introduced in Chapter 2. The chapter consists of a literature review, recent developments and trends, a brief introduction to local traffic regulations and the characteristics of traffic flows at BCIA airport. In Chapter 3 and Chapter 4, we present our system for optimizing the traffic for multiple runways in busy TMA. Chapter 3 focuses on how we solve the optimization problem. It includes an overview of methodology, analysis of the topological design, the mathematical formulation of the optimization problem, and the choice of the optimization method. Chapter 4 explains in detail simulation modules, which are used for objective evaluation. It consists of the computation of 4D trajectories and the conflict detection. The experiments and the numerical results are discussed in Chapter 5. Three kinds of operational patterns are analyzed separately: mixed parallel approach operation, segregated parallel approach operations and independent parallel operations with integrated arrivals and departures. Finally, conclusions and perspectives are presented in Chapter 6.

Chapter 2

State of the art

Literature review

Arrivals are generally given priority at airports over departures. In this section, we focus on the arrival management problem. The activities of managing the arrival flows can be divided into three levels: strategic level, pre-tactical level and tactical level. At strategic level, an initial plan is developed that ensures all arrival aircraft can land on the runway, or all departure aircraft can join different airways; the pre-tactical level coordinates the plan after a collaborative decision making process involving operational partners, such as ATC units and aircraft operators; the tactical level is to update the plan according to real time traffic demand. Our efforts focus on problems associated with optimizing arrivals at the pretactical and tactical levels, including: sequencing, merging, and runway assignment. Among them, sequencing and runway assignment are decision problems, and merging is a control problem.

Sequencing optimization

Sequencing is the queue management of the arrival flows over a time window of 30-45 minutes in TMA. The most commonly used sequencing strategy in the world is the First Come First Served (FCFS). FCFS is simple to implement by controllers. However, it is likely to produce excessive delays and is not suitable for high-density operations, because both the required minimum landing interval between two successive aircraft on the runway and the required minimum separation between two successive aircraft in flight vary with aircraft types and their relative positions, theoretically, it is interesting to shift the position of aircraft in the FCFS sequence to produce an optimal sequence with the shortest makespan (the fastest landing). This optimization problem is denoted as ASP [START_REF] Dear | The dynamic scheduling of aircraft in the near terminal area / Cambridge[END_REF].

ASP has been studied for a long time, in static or dynamic approaches, from theoreti-cal or operational points of view, for single runway or multi-runway operations. In the late of 90's, researchers mainly considered ASP using a static approach, i.e all the aircraft are considered to be in a holding stack, and they can land at any time; thus we can shift their position in the landing queue without limitation. This method does not consider the uncertainty of aircraft operations, aircraft performance or availability of control techniques; it is not realistic from the operational point of view. In 1976, [START_REF] Dear | The dynamic scheduling of aircraft in the near terminal area / Cambridge[END_REF] proposed a decision methodology termed Constrained Position Shifting (CPS), which considered that the composition of the aircraft mix could be changed over time. CPS prohibits an aircraft from being shifted more than a Maximum Position Shifting (MPS) from its FCFS position, and takes the operational constraints on the re-arrangements of the sequence into consideration. This is more realistic than previous approach. In practice, the first generation of air traffic decision support systems for arrival scheduling, such as CTAS in US, COMPAS in Germany, MAESTRO in France, widely used Branch-and-Bound algorithms to schedule arrival aircraft [START_REF] Volckers | Arrival Planning and Sequencing with COMPAS-OP at the Frankfurt ATC-Center[END_REF]Garcia, 1990;Brinton, 1992). Although these scheduling systems could produce a near-optimal delay minimization sequences, and could advise controllers on whether to accelerate, maintain or slow down aircraft to keep this strategic plan, they were found to be difficult for controllers to use in busy TMAs [START_REF] Robinson | [END_REF]. Meanwhile, the MPS is very limited in the mentioned operational systems, normally less than three (Balakrishnan und Chandran, 2006). Later, with the development of Artificial Intelligence (AI) theory, researchers at various labs worldwide experimented with a wide range of AI algorithms for optimizing sequencing in TMA. For example, [START_REF][END_REF] proposed models and heuristic algorithms for real-time control of TMA. For both static and dynamic cases, heuristic algorithms were proposed and computational results were discussed. [START_REF] Robinson | [END_REF] proposed a fuzzy reasoning-based sequencing of arrival aircraft in TMA. This method considered both performance criteria and workload criteria. Their operational test results were very positive. Beasley u. a. ( 2001) developed a population based heuristic algorithm to optimize aircraft landings at London Heathrow. They demonstrated that heuristic algorithms could be successfully applied to CPS-based ASP, as they could rapidly search the near-optimal sequences for controllers.

All the research efforts mentioned above followed a centralized approach in which all aircraft are managed in one step for the entire time horizon of traffic prediction, and this time horizon is normally a snapshot of the entire arrival flows period, see Fig. 2.1. Referring to the initial arrival queue based on ETA, we can apply three optimization methods. In FCFS queue, the make-span is longest, and it is easy for controller to realize. In optimal queue, the make-span is shortest, and we can gain a large margin, compared with the FCFS queue, to accept more coming aircraft. However, it is the hardest for controller to realize. In the optimized queue with MPS <= 3, the make-span is between the two methods mentioned above. The possible gained margin in this case is smaller than the optimal queue, but it is relatively easier for controller to realize. In the centralized approach, we did not consider the dynamics of traffic flows and its uncertainty. The ATC operational environment requires a new methodology which can robustly and directly model the dynamic and uncertain features of the ASP problem over a long period of time with dynamic approach. At the beginning of the 2000's, Hu und Chen (2005b) applied the concept of Receding Horizon Control (RHC) to the problem of arrival scheduling and sequencing in a dynamic environment. The rolling horizon approach divides the problem into multiple periods, enabling the dynamic management of aircraft for large time horizons, see Fig. 2.2. In the decompression approach, we only consider the aircraft queuing problem in a limited window, we apply the optimization method to those aircraft. Hu und Chen (2005a) also integrated the RHC strategy into a Genetic Algorithm (GA) for real-time implementations in a dynamic environment of air traffic control. Simulation results show that a RHC-based GA achieves much better performance than a pure GA without RHC. Later, many studies improved on the performance of meta-heuristic sequencing algorithms, such as Zhan u. a. (2010) ;[START_REF] Samà U. A ; Samà | Rolling horizon approach for aircraft scheduling in the terminal control area of busy airports[END_REF]. Meanwhile, in the domain of deterministic method, a Dynamic Programming (DP) based approach that scales linearly the number of aircraft was presented by [START_REF] Balakrishnan Und Chandran | [END_REF]. They presented an approach for minimizing makespan in the presence of CPS for a single runway, and have demonstrated its effectiveness in a real-world setting at Denver International Airport. The most important contribution of this work is that the approach they presented can handle precedence constraints that can arise from operational constraints or airline preferences, and take into account restrictions on possible arrival times of aircraft. They envisioned that their procedure could replace existing heuristic techniques for computing CPS sequences in decision support tools such as CTAS.

Nowadays, with the complexity of runway configurations in modern airports, research efforts are focused on the study of multi-runway operations. Several studies look at the mathematical model of ASP for multiple runways. For example, Beasley u. a. (2000) modelled the multi-runway sequencing problem in the static case using a Mixed Integer (MI) zero-one formulation, and tested the approach on a problem involving up to 50 aircraft and 4 runways. [START_REF] Bojanowski U | [END_REF] investigated the problem of multi-runway aircraft sequencing at congested airports in a dynamic way. They used an inductive algorithm for scheduling N aircraft from P classes onto M runways. The approach was tested in both static and dynamic frameworks. The algorithm is factorial in the number of runways, exponential in the number of classes, and polynomial in the number of aircraft in each class. Other studies have considered about how to handle the ASP with different multiple runway operational conditions. Kupfer (2009) did a study on the scheduling problem for closely spaced parallel approaches. [START_REF] Farrahi | [END_REF] studied the pair-scheduling problem for landing air-craft in very closely spaced parallel approaches to facilitate pairing of aircraft while meeting a schedule. [START_REF] Samà U. A ; Samà | Rolling horizon approach for aircraft scheduling in the terminal control area of busy airports[END_REF]Samà u. a. ( , 2014) ) studied the RHC-based approach for aircraft scheduling in the terminal control area of busy airports with two runways. Lieder und Stolletz (2016) analyzed the aircraft scheduling problem for heterogeneous and interdependent runways, and defined the problem as a MI problem and provided an efficient DP approach to optimally solve the ASP.

In summary, up to now, MPS is one of the most stringent constraints for the optimization of landing sequence, which generally allows no more than three (MPS <= 3) [START_REF] Balakrishnan Und Chandran | [END_REF]. This constraint limits the potential benefits gained from a more relaxed re-sequencing approach. In addition, the re-sequence problem must consider its compatibility with the dynamic arrival flows in the real world. Current AMAN tools for arrival planning controllers mainly modify the aircraft who have conflict to land on the same runway, but they generally can not provide an appropriate method, such as vectoring, path stretching, speed changes or holding, for the aircraft to meet its time or position in the sequence. For multiple runway operations, it is even more stringent. Thus, in order to gain real benefit from the re-sequencing, control problems should be considered, in addition to the decision problem of sequencing.

Runway assignment

unnecessary flight time, minimize delays and maximize airport capacity.

Several papers have emphasized the importance of the runway re-assignment problem. Isaacson u. a. (1997) studied knowledge-based runway assignment for arrival aircraft in TMA at the strategic level. In this work, the knowledge base for runway assignment uses a set of hierarchical rules and decision logic that evaluates both performance and workload criteria. [START_REF][END_REF] presented a tool named Multiple Runway Planner, developed by Boeing for systematic analysis of arrival sequencing, scheduling and runway assignment with alternative performance objectives from airlines and ANSP. Kim u. a. ( 2014) presented an optimization model for simultaneously assigning aircraft to runways and scheduling the arrival and departure operations on these runways such that the total emissions produced in the terminal area and on the airport surface are minimized. [START_REF] Vela U. A | [END_REF] studied the problem of strategically balancing departure demand at runways in order to reduce departure delays at airports with multi-runway configurations. Delsen (2016) focused on research in flexible arrival and departure runway allocation using Mixed-Integer Linear Programming (MILP) to optimize fuel and noise.

These tools can provide a strategic schedule of runway assignment for controllers, but they do not specify how to realize the real-time runway re-assignment at the tactical level. Changing the runway-in-use is easier for the departing aircraft than for the arrival aircraft. For departing aircraft, controllers can change its runway-in-use by issuing another taxiing route. However, for the arrival aircraft, controllers have to think about changing the runwayin-use as early as possible. This is because with the conventional route network, it is difficult for controllers to deviate aircraft to runway other than the scheduled one, as there are very few spaces to insert aircraft in a landing queue on middle or final approach segments. Thus, to realize a successful runway re-assignment decision, the control problem needs to be considered as well.

Merging control

Merging is used to handle real-time traffic flows at the tactical level. It is a control problem. TMA is a transition airspace between airports and the network of airways, and it is set up in the vicinity of one or more major airports. Due to a large number of ascending and descending aircraft, the traffic complexity in TMA is very high. The more aircraft approaching the runway, the more difficult it is for controllers to change the FCFS schedule of arrivals. The way of merging the aircraft at each approach segment has to consider not only adapting to the strategic schedule, but also operation streamlining around merging points. Recent conventional methods of merging aircraft flows mostly rely on radar vectors, such as heading change, speed control, and flight level change, issued by controllers. This method is flexible, however, controllers and pilots have to keep a high frequency of radio communication, which can easily produce high workloads for both. Nowadays, based on new emerging technologies in CNS, the controller can sometimes choose CTA for the aircraft; the crew then uses the aircraft Flight Management System (FMS) to fly the aircraft towards the Required Time of Arrival (RTA). Avionics-based ATM operations and cooperation between air and ground by sharing flight trajectory information produces significantly more efficient merging in TMA.

On the side of investigating merging topology design, a successful example is the Point Merge System (PMS). PMS is a systemized method for sequencing arrival flows developed by the EUROCONTROL Experimental Center in 2006. It is a RNAV-based route topology design, using predefined legs at iso-distance to the center merging point for path shortening or stretching. It creates a linear holding pattern instead of traditional holding stacks. Meanwhile, it creates a novel kind of merging and sequencing topology, which will been introduced in detail in Chapter 3. In addition to conventional tree-based topologies, other RNAV-based topology designs have also been considered. For example, Zúñiga u. a. (2013); Chida u. a. (2016) studied the automated merging problem with a fish-bone shaped merging topology, as an extension of a tree-merging model. Polishchuk (2016) studied an open merging topology, instead of merging at one unique way-point, aircraft merge close to an area near the final approach. The different kinds of topology designs are shown in Fig. 2.3. These topologies can merge n flows into 1 main flow for single runway, but can not be used to merge n flows into 2 flows for two parallel runways. Meanwhile, other researches have considered a range of topics around merging based on RNAV and RNP routes. For example, Becker u. a. (2004) ;Becher u. a. (2005) conceived the Spacing of Performance-based Arrivals on Converging Routes (SPACR) concept for allowing aircraft to remain on RNAV routes while managing the flow of the aircraft to the final approach segment. Alam u. a. (2010) proposed a methodology to generate aircraft-specific dynamic CDA routes in 3D. This methodology involves the discretization of the terminal airspace into concentric cylinders with artificial way-points and uses enumeration and elimination (based on aircraft performance envelopes) from one way-point to another to identify all possible routes. However, the route for each aircraft is different, and the set of routes is too complex to be implemented. The RNAV-based route structures provide more flexibility for merging control, and will play an important role in future ATM systems.

Concerning automated merging, Niedringhaus (1995) proposed a Stream Option Manager (SOM) concept using linear programming techniques. SOM can find modified flight paths for each aircraft that satisfy all requirements to stay as close as possible to the pilot-preferred paths, however, it has only be implemented in 2D. With two main ATM modernization programs, SESAR and NextGen, some studies have demonstrated the potential benefits of integrating automated arrival management with conflict-free and airborne spacing. For example, [START_REF][END_REF] analysed the effects of automated arrival management, airborne 2009) presented arrival management architecture and performance analysis with advanced automation and avionics capabilities. The benefits of noise abatement, fuel consumption, and CO 2 emissions have also been analyzed (Korn u. a., 2006;Haraldsdottir u. a., 2007;Scharl u. a., 2008;Coppenbarger u. a., 2009;[START_REF][END_REF]Toratani u. a., 2015b,a). Again, automation can achieve high efficiencies in ATM, however, due to safety considerations, full automation in ATM is still a challenge. Facing extremely dense operations in complex TMA in the near future, we can consider automating a large portion of routine operations including scheduling, sequencing, and spacing [START_REF] Isaacson U. A | [END_REF]. We have to consider that when humans participate in managing recovery from abnormal events, they prefer traffic complexity as lower as possible.

Recent trends in transition airspace TMA

CDA and CCO

Continuous Descent Approach (CDA) and Continuous Climb Operation (CCO) are aircraft operating techniques in TMA enabled by airspace design, instrument procedure design and facilitated by ATC, see Fig. 2.4. The ideal CDA starts at the top of a descent and ends when the aircraft starts the final approach and follows the glide slope to the runway. Typically, in a conventional, non-CDA approach, the aircraft descends stepwise, with portions of level flight in-between. By performing a CDA the aircraft remains higher for longer and operates at lower engine thrust. Both of these elements induce a reduction in fuel use, emissions and noise along the descent profile prior to the point at which the aircraft is established on the final approach path. Similarly, the ideal CCO is a continuously fuel optimal climbing path with an optimal fuel-conserving rate from the runway to the top of climb. The fuel used in climbing to the most fuel efficient cruise level can be a significant part of the overall fuel used for the flight. CCO allows the aircraft to reach the initial cruise flight level at optimum air A single CDA and CCO compared to a non-optimised climb or descent profile can result in fuel savings of 50 -200 kilograms of fuel per flight (Eurocontrol, 2009b). ICAO estimated that savings from the planned implementation of CDA and CCO in Europe could save as much as 500 kilotons of fuel per year. In addition, using CDA can reduce noise by 1-5dB compared to a non-CDA operation. However, for many airports, the opportunity to implement a CDA is very limited because of the high volume of air traffic on approach and in the vicinity of the airport, especially during busy daytime periods. When approaching traffic is heavy, a pilot may need to adjust throttles, flap settings, and extend landing gear to maintain safe and consistent spacing with other aircraft in the terminal airspace. Extending flaps, and landing gear increases drag, which requires the application of additional thrust to keep the aircraft flying at the same speed. Nonetheless, in Europe, more and more airports are to use CDA as much as possible and to gradually increase the percentage of CDA flights. In addition, when CCOs or CDAs are used, appropriate airspace design and ATC procedures should be used to avoid the necessity of resolving potential conflicts between the arriving and departing traffic flows through ATC level or speed constraints.

PBN and PMS

PBN is a new operational concept presented by the ICAO in 2006, for the purpose of integrating the operational practices and technical standards of RNAV and RNP in various parts of the world. It enables airspace designers to develop and implement new automated flight paths that increase airspace efficiency and optimize airspace use. Thus, it enables a safer and more accurate flight model and a more efficient ATM operation, including enhanced safety, increased efficiency, reduced carbon emissions, and reduced fuel costs.

As shown in Fig. 2.5, conventional routes use ground-based navaids, such as Very-High and navigation in which the flight path had no longer to be tied directly to over flying the ground navigation stations. RNP is built on RNAV, it allows an aircraft to fly a specific path between two 3D-defined points in space. Compared with RNAV, RNP systems require onboard performance monitoring and alerting, and can provide better use of airspace, such as a seamless vertical path, a Direct-to-a-Fix (DF), a curved path etc., which can be successfully supported by modern FMS [START_REF] Herndon U. A | [END_REF][START_REF] Herndon U. A | [END_REF].

PBN, combining advanced on-board equipment with satellite-based navigation and other state-of-the-art technologies, covers all phases of flight from en route and terminal area to approach and landing (the detailed PBN navigation specification defined for different flight phases is introduced in A.3). Curved paths are a key PBN capability enabling precise departure, arrival and approach procedures. The curved path capability is accomplished through Constant-Radius-to-a-Fix (RF) legs that enable procedure designers to adapt SID, STAR, and approach segments in ways not possible with straight segments. The RF leg provides the procedure designer with the flexibility to avoid specific terrain, abate noise or protect environmentally sensitive areas, avoid restricted airspaces, make arrivals more efficient. Thanks PMS is a systemized method for sequencing arrival flows. It is a RNAV-based route topology design, using predefined legs at iso-distance to the center merging point for path shortening or stretching. It creates a linear holding pattern instead of traditional holding stacks. In the PMS system, aircraft remain in lateral mode, after they enter the sequencing leg, they fly at an economical speed, and when Direct to instruction is issued at the specific time from controller to pilot, the aircraft perform a turn towards the merge point, at the same time it can perform a CDA descent. The airborne separation between preceding and trailing aircraft is maintained only by speed adjustments ordered by controllers. The benefits of PMS are clear:

1. It reduces controllers' workload under high traffic demand in TMA. First, because heading instructions needed in this system are much less than conventional radar control method, controllers' workloads are dramatically reduced, especially for the feeder controller who integrates the flow in the final approach phrase. Boursier u. a. ( 2007) presented small-scale experimental results, and found the PMS method to be efficient, safe and accurate. Ivanescu u. a. ( 2009) compared the PMS method with a radar vectoring method 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. 2. It is an innovative linear holding pattern, instead of a conventional stack holding mode.

Consequently, it is able to easily adapt to the different volumes of traffic in TMA. In low capacity periods, its sequencing legs can be shorter, and in high capacity periods, its sequencing legs can be longer. From this point of view, its route structure is very suitable for handling dynamic traffic flow.

3. The new arrival flow integration technique enables ATM to successfully manage more complex scenarios and realize an advanced continuous descent [START_REF][END_REF]. Favennec u. a. ( 2009) performed an investigation on how PMS could be adapted to typical terminal area configurations with more complex environments. It is now one of the ICAO Aviation System Block Upgrades and is referenced as a technique to support CDA (ICAO-Doc.9931, 2008). Up to 2016, PMS has been successfully implemented in Oslo, three Norwegian regional airports, Dublin, Seoul, Paris ACC, Kuala Lumpur, Lagos, Canary Islands, Hannover, London City and Biggin Hill. 4. PMS also supports: a better pilot situational awareness; more orderly flows of traffic with a better view of arrival sequences; a better trajectory prediction, allowing for improved flight efficiency; a standardisation of operations; and better airspace management.

Local traffic in Beijing TMA

Known as "China's No.1 Gateway", BCIA is the most important, largest and busiest international aviation hub in China. As the growing international airline network connects into Beijing, BCIA is becoming one of the busiest airports in the world. Everyday, 1700 flights from 94 airlines connects Beijing with 244 cities in 54 countries around the world. In this thesis, we use BCIA as a case study. Referring to the official Aeronautical Information Publication (AIP), the declared operational capacity at BCIA is 88 movements per hour, including departures and arrivals. As shown in Fig. 2.7, it operates under saturation from 6:00 to 24:00. From 6:00 to 9:00, there is a heavy demand for departures. From 21:00 to 24:00, there is a heavy demand for arrivals. Between 9:00 to 21:00, the number of arrivals and the departures almost equal each other.

Flight procedures and APP services

As shown in Tab. 1.1, APP services and some parts of ACC are provided to air traffic in TMA. There are designated SID routes for guiding the departing aircraft into the airways, and STAR routes for helping the arrival aircraft approach the destination airport. As shown in Fig. 2.8, there are six entry points in Beijing TMA for arrival flights. Among these, there are four entry points in the South: JB, BOBAK, VYK and DOGAR, and two entry points in the North: KM and GITUM. The arrival flows are very heavy in the South. For departure flights, there are seven exit points, including YV, CDY, TONIL, LADIX, RENOB, SOSDI, and KM. The point KM is special, it serves both for departure and arrival flights.

Radar control has been implemented within Beijing TRACON. The minimum horizontal radar separation is 6km for aircraft within 50km of Aerodrome Reference Point (ARP) and 10km for aircraft beyond 50km of the ARP. The minimum vertical radar separation is 300m (1000ft). Normally, arrival aircraft are vectored and sequenced from VYK, KM and JB, BOBAK, GITUM, DOGAR or some specific transfer of control points to the appropriate final approach track, or to the time when runway is in sight. Instructions about radar vectors, ascent (descent) altitudes or speed adjustment will be issued for spacing and separating the aircraft, so that stipulated radar intervals and wake intervals are maintained, taking into account aircraft characteristics or control regulations. During rush hour in a day, arrival aircraft will be vectored. Radar vectoring tracks will be different to the published STARs. The landing aircraft of medium type and below should fully vacate the runway within 50 seconds after flying over the runway threshold, and landing aircraft of heavy type and above should fully vacate the runway within 70 seconds after flying over the runway threshold. Departing aircraft operate according to SID procedures or are vectored to join the SID routes by a radar controller. Published STAR and SID charts can be found on the following website: http://www.eaipchina.cn/Version/201513/.

Parallel runways operations

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.

As shown in Fig. 2.9, 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 departures and arrivals, runway 36R-18L is mainly used for departures, and runway 01-19 is mainly used for arrivals. During departure rush hour, the three parallel runways will be used for departure. During arrival rush hour, the three parallel runways will be used for arrivals. The distance between runway 18R-38L and runway 18L-36R is 1960m, the distance between runway 18L-36R and runway 01-19 is 1525m. Theoretically, independent instrument parallel approaches and parallel departures are feasible at BCIA. However, dependent parallel approaches are operationally used between adjacent runways. This means that: (i), a minimum of 1000 ft (300 m) vertical separation, or a minimum of 3Nm (5.6 km) radar separation, are provided between aircraft during turn-on to parallel ILS localizer courses; (ii), a minimum of 3Nm (5.6 km) radar separation is required between successive aircraft on the same ILS localizer course, unless increased longitudinal separation is required due to wake turbulence; and (iii), a minimum of 2Nm (3.7 km) radar separation is provided between successive aircraft on adjacent ILS localizer courses, see Fig. 2.10. Here, we assume that independent instrument parallel approaches are used between runway 18R-36L and runway 01-19, and that segregated instrument parallel approaches and/or departures are used between runway 18L-36R and runway 01-19, see Fig. 2.11.

Local traffic characteristic

We collected several days of real flight data in June, 2017. We observed different kinds of flight operations in the Beijing TMA, for example helicopter, general aviation, over flights, departing flights, and arrival flights etc. We selected trajectories that either departed from or arrived at BCIA. In a given 24-hour time period, the operation pattern of operations at BCIA airport can vary from day to day. Further, there are South-inbound and outbound operations (Magnetic number of the runway-in-use (QFU)18) and North-inbound and outbound operations (QFU 36). Therefore, there was a change of runway-in-use for some flights. After separating the traffic according to the QFU, Fig. 2.12 shows the characteristics of local traffic at BCIA. We observed that: first, QFU 36 is more frequently used, since most arrival traffic is from South, especially via the entry point DOGAR (red colour). Second, a large proportion of the trajectories are manoeuvring flights used for absorbing delays. At the same time, the available manoeuvring airspace is very limited. For example, in QFU 18, in the north of airport, some trajectories from DOGAR are deviated by controllers over the airport, so as to balance the landing rate on runways. Third, departure trajectories are less Chaotic than arrivals. However, departing flights have to share the TMA with arrivals. Hence, since there are many arrival deviations, there are a large number of departure deviations as well. For example, in QFU 18, many of the departing flights to the North are deviated (dark blue colour).

We continuously selected one day traffic to study the South-inbound traffic characteristic at BCIA. In order to reduce the noise of trajectories for each arrival flow, we clustered them with a preprocessing algorithm. In the end, the South-inbound traffic characteristic Figure 2.12: Four kinds of traffic pattern with real traffic data Figure 2.13: South-inbound daily real traffic and clustered traffic is demonstrated in Fig. 2.13. We use this clustered traffic as a baseline, against which to compare our numerical simulation results. In Fig. 2.13, we see that flights from JB will usually join to BOBAK first, 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 little manoeuvring airspace. As a result, they execute frequently the holding procedure. However, flights from KM normally land directly with very few manoeuvring actions. Further, there is a large manoeuvring space near KM, which is used to absorb traffic flows from JB, BOBAK, and VYK. There is another small manoeuvring space near runway 01-19, which is used to absorb traffic flows from GITUM, VYK and DOGAR. All three runways are used for the South-inbound landing in order to accommodate the dense traffic flow.

Conclusion

In this chapter, we first review the literature with the numerous models and optimization methods, which have been proposed to solve the problem of arrival flow management in TMA. Especially, our efforts focus on three optimization problems at pre-tactical and tactical levels: sequencing, merging and runway assignment. For landing sequencing optimization problem, FCFS queue, optimal queue, and MPS based optimized queue managements are analyzed in static and dynamic way. It is pointed out that a feasible optimization solution for resequencing has to consider the constraint from controllers' workload, and M P S <= 3 is a strict constraint for further improving the efficiency of sequencing in real world. For runway assignment problem, with the conventional route network, runway re-assignment will lead to heavy workload for controllers to deviate aircraft. Its bottleneck is how to realize an efficient real time runway re-assignment at tactical level. For merging control, different merging topology and automated merging control are discussed. Emerging techniques in CNS provides a more efficient approach to address the merging control problem.

Then, we introduce the recent develop trend in transition airspace TMA, including the concepts of CDA and CCO, PBN and PMS. CDA and CCO can significantly improve the climbing and descending performance of aircraft in TMA, consequently reduce the fuel consumption and noise. However, they are very limited on busy airports with high volume of air traffic, due to re-sequencing and merging control of aircraft. Efforts in Europe are now paid on using CDA and CCO at more and more airports to the extent possible, and on gradually increasing the percentage of CDA and CCO. On the other side, PBN and PMS are two new concepts to improve the use of airspace. Curved path is a keyPBN capability enabling precise departure, arrival and approach procedures with more flexibility. It allows optimized use of airspace. PMS is a systemized method for sequencing arrival flows. It shows many benefits on the sequencing and merging of air traffic flows in TMA in terms of controllers' workload, use of airspace, linear holding, and delay etc.

In the end, we conclude the detailed operational information about local traffic at BCIA, including flight procedures, APP services, parallel runway operation modes. In addition, we analyze the local traffic characteristics at BCIA. The available airspace in Beijing TMA is very limited, the main arrival flows come from and go to South. Its hourly operational capacity is 88 movements per hour. Total movement per hour from 6:00 to 23:00 is over its capacity. It is necessary to find an optimization approach to manage the traffic flows in TMA, increase its capacity and reduce the potential delay. The proposed approach should be based on the emerging techniques in ATM/CNS. It should be able to improve the performance of air traffic management system. It should be feasible for the future implementation.

In the next chapter, we will introduce the optimization approaches used to address this problem.

Chapter 3

Optimization methods and algorithms 3.1 Methodology overview Fig. 3.1 gives an overview of the methods discussed in this chapter. We divide the methods into three categories:

1. Route network optimization: aims to design a novel RNP-based route network to support efficient control of trajectories, to meet management requirements for arrival flows to parallel runways. 2. Optimization module: is a decision optimization solver. The designed optimization algorithm searches the decision variables to meet the constraint requirements, to reach some optimized results, such as conflict-free, minimizing delay. 3. Simulation module: which is split into several components concerning all necessary modelling approaches to compute a realistic 4D trajectory. It is able to calculate airspeed, aerodynamic forces, motion of aircraft, as well as fuel consumption, flying distance, flying time (or transit time), etc.. which are the operational performance indices used to support decisions made by different actors.

Note that, over the course of our research period, we built three versions of the simulation models. The first version consists of three decision variables and a constant speed profile to simulate the segregated parallel approach operation. The second version has four decision variables, a step-by-step variable speed profile, and a near CDA descent design to simulate the mixed parallel runway operation. The third version improves on the second one by applying the BADA1 -based flight performance model, resulting in more precise trajectory and fuel consumption. In this thesis, we use the third version to simulate the flight trajectory in

Route network design for parallel runways

Horizontal profile design: strategy of sequencing and merging

The topology of arrival and departure route networks plays an important role in commercial aircraft operations. In the horizontal route network design, we have to consider how to perform sequencing and merging, and how to implement the optimized sequencing decision by tactical merging control.

It is well known that the most important requirement for high density operations is to maximize capacity. This means that in a fixed time interval, an airport can safely accept as many aircraft as possible to land. As previously discussed, given an initial arrival queue based on Estimated Time of Arrival (ETA), and the required minimum separation between aircraft, three possible sequences that can be proposed for controllers: the FCFS, the optimal queue, and the CPS-based optimized queue. The FCFS queue has the longest make-span, however, it is easily implemented by controllers. The optimal queue has the shortest makespan, and can provide more available space to accept additional coming aircraft in a fixed time interval. Although, it can maximize the landing throughput, it is the hardest for controllers to implement. The optimized queue with MPS ≤ 3" is a compromise between the FCFS queue and the optimal queue. It combines the problem of optimization and implementation. The resulting capacity depends on the MPS. Generally, the larger the MPS, the more benefit will Figure 3.2: Radar vector techniques be gained. But, in a busy TMA, one small re-sequencing action to aircraft can increase heavy workload for controllers. Controllers' workload is the key limitation for the augmentation of MPS. As a result, to improve the performance of CPS, we need to increase the flexibility of sequencing, so that we can increase the upper limit of MPS. In this context, we have to study the nature of manoeuvring control.

From a microscopic point of view, conventional re-sequencing control techniques require aircraft to change their headings, so as to deviate from their planned route. Different deviation techniques impose different workloads to controllers and pilots. Fig. 3.2 shows six radar vector deviation techniques: "Dog-leg", "Detour", "Parallel offset", "Short-cut", "Holding", and "Direct-to". For each action, a minimum number of heading changes are required. For example, for "Dog-leg", "Detour", "Short-cut", and "Holding", a minimum of two heading changes are required. For "Parallel offset", a minimum of three heading changes are required. For"Direct-to", only one heading change is needed.

From a macroscopic point of view, there exists a conventional way of handling the complex traffic in TMA, it is based on the experience of radar controllers. Its topology is shown in Fig. 3.3 and can be summarized by a rectangular shape consisting of four sequencing legs. Traffic coming from different directions will join the associated sequencing leg in order to integrate with the main landing flow. There are: crosswind leg, downwind leg, base leg, and final leg. If we want to re-sequence the landing queue, the six types of radar vector techniques can be applied. If traffic flow increases, then we can double the "downwind leg" to absorb more traffic. With this conventional way of sequencing and merging, the re-sequencing actions will dramatically increase the workload of controllers. Moreover, the trajectories will be chaotic under heavy traffic loads, such as the recorded trajectory sample in Beijing TMA.

In order to overcome the constraints of MPS ≤ 3, while avoiding overly heavy workloads for the controllers, we propose a novel route topology based on the existing PMS, named the Multi-Level and Multi-Point Merge System (MLMPMS), or ML-PM for short. First, considering the diversity of operational modes for parallel runways, we design two kinds of horizontal topologies. One is designed for segregated parallel approach operations, the other is designed for mixed parallel approach operations, see Fig. 3.4. The shaded areas blue and pink are the possible merging areas for aircraft. Each coloured area indicates the possible manoeuvring airspace for some specific STAR route. We can see that in the conventional topology, different coloured areas are well separated from each other, while in the proposed ML-PM topology, the different coloured areas have a common part. With a limited airspace, the proposed ML-PM topology may provide a greater manoeuvring area for aircraft compared to the conventional topology. Moreover, in the conventional topology, controllers issue at least two radar vector instructions, usually three instructions in method of "parallel offset", to guide the aircraft to the nearest runway. In the ML-PM topology, controllers issue only one "Direct to" instruction to guide aircraft to the merge point. In addition, in the mixed ML-PM topology, aircraft can be guided to change to an alternative runway by holding its flight level on the sequencing leg. This is very convenient for dynamic runway allocation. Conversely, runway re-allocation at tactical level is very hard or sometimes impossible to do using conventional topology and segregated ML-PM topology, due to safety considerations. Taking the runway reallocation in Beijing TMA as an example, controllers can only change the allocated runway for those flights arriving in the TMA with the same orientation as Figure 3.4: Sequencing and merging topologies for parallel runways QFU. Otherwise, they vector the aircraft to join another runway by deviating them to fly over the airport. Second, because there is an overlapping area in the merging zones for the ML-PM topology, we have to design the vertical plane, so as to reduce traffic complexity. The idea is shown in Fig. 3.5. Aircraft with different wake turbulence categories and arriving from different directions, fly on the sequencing leg with different flight levels2 . There are three parallel flight levels for a group of aircraft coming from the same entry point, "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 a unique projection onto the horizontal level. After joining the sequencing leg, aircraft in the same category need to keep their CAS and flight level until turning into the merge point.

Economical descent profile design

As previously mentioned, CDA can reduce fuel consumption and noise compared with the conventional step descents method. However, it is hard to be totally implemented. In CDA, the aircraft descend from en-route TOD (Top of descent) to touchdown with low engine thrust settings and a low drag configuration, thereby reducing fuel burn and emissions during descent. Generally, the optimum vertical profile takes the form of a continuously descending path (near 3 degrees). However, in reality aircraft normally maintain a part of level flight to be transferred from one sector to another sector, and there is a great difficulty of implementing CDA in the congestion circumstance. As a result, in the ML-PM system, we design a target altitude at each way-point to guide aircraft to execute a near-CDA descent. In this kind

Mathematical formulations

Assumptions

We assume that there is a set of aircraft F = {1, 2, ..., n} planned to land at BCIA in a no wind condition. Aircraft fly in an ISA environment (Appendix A.4). For each aircraft i ∈ F, the following six kinds of data are given:

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

t e

i -the ETA at the entry point.

v e

i -the initial CAS of aircraft i at the entry point. 4. r e i -the initial landing runway. For aircraft coming from KM, JB and BOBAK, r e i = 1. For aircraft coming from VYK, DOGAR and GITUM, r e i = 0. 5. ET A L i -the estimated time of landing. 6. cat i -the wake turbulence category.

State space

For each arrival aircraft i, there are four decision variables: All these variables are modelled by some discrete variables. The time t E i is adjusted by a number of slots denoted by j, and the duration of each slot is ∆ = 2s,

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

v E i = v e i (1 + g), g = 0, ±1%, ±2%, ..., ±15%. (3.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.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. For the runway assignment, r i is defined as below:

r i =    1 if i expected to merge at M1, 0 if i expected to merge at M2. (3.4)
Here, M1 and M2 are merge points in ML-PM. For BCIA case, we can decide that M1 is the merge point connecting to runway 18R-36L, and M2 is another merge point connecting to runway 01-19. The coordinates of M1 and M2 will be designed later.

Constraints

Some operational constraints must be carefully considered. They are crucial to the fairness and safety of aircraft. First, t E 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 Und Chandran | [END_REF][START_REF] Carr U. A | [END_REF][START_REF] Lee Und Balakrishnan ; Lee | Hamsa: A study of tradeoffs in scheduling terminal-area operations[END_REF]. In this paper, we suppose that the earliest time of arrival is limited to 5 minutes before ETA, and the latest time of arrival time is 15 minutes after ETA, then we have

t e i -5min ≤ t E i ≤ t e i + 15min. (3.5)
Second, speed change is limited by the performance of commercial aircraft in descent Figure 3.6: Minimum separation profile. In TMA, when aircraft fly below 10000ft, their airspeed must not be higher than 250kt due to bird ingestion damage, and at the same time it must not be less than the minimum clean configuration speed due to low speed stall, thus v E i is subject to:

v e i (1 -15%) ≤ v E i ≤ v e i (1 + 15%), (3.6
)

v e i (1 -15%) ≥   
230kt for Heavy aircraft, 220kt for Medium aircraft.

(3.7) Third, t T i must comply with geographic constraints on the sequencing leg:

t T imin ≤ t T i ≤ t T imax . (3.8)
Fourth, the required minimum radar separation and wake turbulence separation must be considered to avoid conflicts. According to ICAO regulations, two aircraft are considered to be in conflict if their horizontal separation is less than the minimum separation standard and their vertical separation is less than 1000ft, see Fig. 3.6 . In the ML-PM system, vertical separation is partially assured by vertical profile design, especially at the entry points of sequencing legs. If aircraft can not match the vertical separation, then a horizontal separation must be assured.

In the approach airspace, we have to consider two kinds of minimum separation: wake turbulence minimum separation and approach radar separation. The required aircraft minimum separations in the TMA, denoted by s min i,j are listed in Table .3.1a. Since time-based system is much more convenient for detecting conflicts and metering flows, we convert the distance-based separation to a time-based separation (Nikoleris u. a., 2014). The reference velocity used for this computation is based on the average final approach speed of commercial Table 3.1: ICAO minimum separation aircraft in the specified weight category and assuming no wind. This is also the minimum speed during all the approaching phases. Thus, based on this assumption, a safety margin have been considered. The minimum time-based wake turbulence separations are showed in Tab.3.1b and are used for computing conflicts. It should be noted that commercial aircraft normally follow an airline CAS speed schedule, thus their final approach CAS is usually relatively steady. For other special cases, such as when a trailing aircraft is significantly faster than the leading aircraft on the final approach, or when there is a heavy crosswind at airport, the time separation must be adjusted as necessary to ensure that the minimum required separation distance is not violated anywhere on the final approach. In this thesis, we do not explicitly consider such cases.

Objectives

There are a number of different stakeholders in the air transportation system, including ANSP, airlines, airports, and government. Each stakeholder has their own set of objectives. For example, ANSP aims to ensure safety and efficiency of aircraft. Airlines focus mainly on minimizing fuel costs, maximizing flights punctuality, and keeping fairness between airlines. The airports goal is to minimize delay for departure and arrival flights. Government's preference is to minimize the environmental effects, such as noise and pollution. As a result, there is a trade-off between various and sometimes conflicting objectives. Consequently, our optimization problem is a multi-objective optimization problem.

A multi-objective optimization problem is an optimization problem that involves multiple objective functions. In mathematical terms, a multi-objective optimization problem can be formulated as

min(f 1 (u), f 2 (u), • • • , f k (u)) (3.9) s.t. u ∈ U, (3.10)
where the integer k ≥ 2 is the number of objectives and the set U is the set of feasible decision vectors. An element u ∈ U is called a feasible solution or a feasible decision. A vector z := f (u) ∈ R k for a feasible solution u is called an objective vector or an outcome. In multi-objective optimization, there does not typically exist a feasible solution that minimizes all objective functions simultaneously. Therefore, we need to consider compromised solutions.

A priori method is applied to convert the multi-objective problem into a single-objective problem (scalarizing). By varying the parameters used in the scalarization, different optimal solutions are obtained. The most well-known scalarization is the linear scalarization,

min u∈U k i=1 α i f i (u), (3.11)
where the weights of the objectives α i > 0 are the parameters of the scalarization.

In the vertical profile design of ML-PM, we have already improved the descent profile by using a set of target levels on the way-points. Consequently, we mainly focus on the four following objectives: conflict-free, minimum delay, fast-landing makespan, and constrained position shifting. For a flight i, we gather all decision variables into a vector

u i = (t E i , v E i , t T i , r i ).
We define the position of aircraft i as p i = (x i , y i , h i ), then the objective function of our optimization problem is:

z = min C + α 1 D + α 2 S + α 3 P, where (3.12) C = n i=1 C i ( u i , p i ), (3.13) D = 1 n n i=1 D i , (3.14) D i =      t L i ( u i ) -ET A L i 2 if t L i ( u i ) > ET A L i , 0 otherwise. (3.15) S = 1 n (t L last -t L f irst ) (3.16) P = n i=1 P i ( u i ) (3.17) P i = P Actual -P F CF S . (3.18)
Here, in Equation 3.12, C is the total number of conflicts, D is the average square delay from ETA, S is the average landing interval, and P is the total position shift. The first term aims to minimize the number of conflicts. The second term searches for a near-optimal resolution with less deviation from the initial landing time. The third term aims to speed up the whole landing process. The fourth term aims to minimizing the total positions of shift. In addition, α 1 , α 2 and α 3 are weighting parameters. These parameters are chosen to fit with the values of the magnitude of D, S, P . In our case, D ∈ [0, 100], S ∈ [0, 50], then, the default setting are α 1 = 0.01 and α 2 = 0.02. We don't care about P in mixed parallel operations, so in this case α 3 = 0. In segregated parallel approach operations, α 3 = 1.

Then in Equation 3.13, n is the number of flights, C i is the number of conflict encountered by aircraft i (it depends on u i ), in Equation 3.14, D i is the square delay of aircraft i, in Equation 3.15, t L i is the actual landing time of flight i on the runway (it depends on u i ), and in Equation 3.16 t L last and t L f irst are the actual landing times of the first landing aircraft and the last landing aircraft at the airport. In Equation 3.17 and 3.18, for aircraft i, P i is the number of position shifts between the aircraft's position in the actual queue and the FCFS queue. P i is only applicable to the segregated parallel approach operation, because when an aircraft change their initial landing runway, then their position in the new queue is not the same as in the previous queue.

Optimization method

Choice of optimization algorithm

The choice of the optimization method depends on the problem considered and the choice of mathematical modelling techniques (discrete variables, continuous or random modelling uncertainties or not, etc.). It also depends on the nature and properties of the objective function and constraints: linearity, convexity, presence or absence of constraints, the existence of an analytical formulation of functions and derivatives. Further, it depends on the goal as well, to find a local optimum in a certain neighbourhood, or to search for a global optimum of the objective function. Therefore, we have to think about local or global optimization, continuous or discrete (combinatorial or integer) or mixed (combining continuous and integer) variables, with or without constraints, deterministic or stochastic, convex or non-convex. In terms of implementation, optimization methods usually use iterative algorithms. From a starting point chosen in the search space, or more points for the population algorithms, they try to iteratively improve the objective function. These algorithms are themselves either deterministic or stochastic. In the latter case, they use a random walk guided by a heuristic. These are called meta-heuristics, and include such as evolutionary algorithms, SA, differential evolution, and particle swarm.

The ATM system is complex and dynamic. In real situations, controllers require an algorithm which can quickly find a good solution (near-optimal) rather than an optimal solution achieved after a lengthy computation. In addition, controllers prefer that the optimization program can be periodically updated, which can enforce their observation of dynamic traf-fic. Exact methods, such as Mixed Integer Programming (MIP), Integer Programming (IP), Constraint Programming (CP), and the Satisfiability Problem (SAT), can always find an optimal solution if one exists, but on larger instances their running time increases dramatically. Further, to use a deterministic method for our optimization problem, we would have to decide the sequence of aircraft first. In the merging zone of the designed ML-PM topology, aircraft on the sequencing leg have the same opportunity to join the merge point, because they have the same distance from the merge point. They may dynamically change runway allocation as well. In these circumstances, a deterministic method will not be suitable to handle the mixed parallel runway operational pattern, because of high induced combinatorics. Stochastic methods, such as the Hill Climbing (HC) algorithm and SA, can quickly find near-optimal solution, especially on large instances. For the aircraft sequencing problem, HC is faster than SA, however SA produces solutions with better quality. Both are very flexible in modelling additional requirements [START_REF] Fahle U | [END_REF]. Consequently, in addition to be a suitable approach for solving the current optimization problem for parallel runway operations, we can easily add new decision variables, such as route selection, to model more complex traffic situation in the future.

Moreover, in our optimization problem, the objective function is computed thanks to a simulation process, which considers not only the decision problem but also the control problem. 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 E i , v E i , t T i , r i ), the number of possible control solutions for each flight is N S i = 600 × 30 × 100 × 2. For n flights, at time t, the total number of the possible trajectories is N T raj all |t = n × N S i . To guarantee a conflict-free solution is the most important objective. The position of the aircraft at every moment (x i , y i , h i , t i ) has to be considered for conflict detection, i.e. the computation of C i . If we use a pairwise conflict evaluation method, then for the duration of scenario T scenario , then the number of conflict evaluation positions for an aircraft is N C i all = N T raj all |t × T scnario have to be checked. Consequently, the combinatorics associated with 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 time window, we have decided to address such global optimization problem by using a stochastic heuristic approach.

RHC-based system dynamics and control

Our optimization problem needs to consider the space-time trajectories of each aircraft, as well as the control problem of each aircraft. Recall ( p i ) t is the position of aircraft i at time t, we consider a discrete-time linear system of the form,

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

n i=1 ( p i ) t , n i=1 ( u i ) t ∈ C t i , i ∈ F (3.20)
where C t i is the constraint set. The instantaneous objective value of the 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 using the objective function,

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

Based on system dynamics and control theory, we choose to solve a 24 hour traffic optimization problem by applying the RHC technique, which decomposes the original problem into several sub-problems. The RHC policy works as follows. At stating time t INIT , we consider a time interval, called a sliding window, extending k steps towards the final 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.

We then carry out the following steps:

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. Minimize 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.

z t N i=1 ( p i ) t , N i=1 ( u i ) t , i ∈ N ⊂ F (3.22) subject to N i=1 ( p i ) t = f N i=1 ( p i ) t-1 , N i=1 ( u i ) t-1 (3.23) N i=1 ( p i ) t , N i=1 ( u i ) t ∈ C t i (3.24)
where, N ⊂ F is the number of aircraft in the current sliding window.

Fig. 3.7 illustrates how the sliding window operates along the time axis of 24 hours. The first sliding window begins at t INIT , and a SA optimization algorithm is applied to the corresponding time interval [T s (0), T e (0)]. Next, the sliding window moves toward in the future by S, and the current optimization interval becomes [T s [START_REF] Lee Und Balakrishnan ; Lee | Hamsa: A study of tradeoffs in scheduling terminal-area operations[END_REF], T e [START_REF] Lee Und Balakrishnan ; Lee | Hamsa: A study of tradeoffs in scheduling terminal-area operations[END_REF]]. We repeat this process until we reach the k th sliding window with T e (k) = t FINAL . According to the relative relationship between the aircraft arrivals [t E i , t L i ] and the active time window [T s (k), T e (k)], we classify the status of aircraft into four categories: Completed, On-going, Active and Planned. Completed means that the aircraft has already landed. On-going means the aircraft has not landed yet, but the decisions on changing its trajectory have already been taken, its descent trajectory is frozen. Active means the decisions for this aircraft can still be changed, thus its trajectory is not frozen. Planned means the aircraft is not in the current window. Only Active and On-going aircraft will be included in the optimization process step 2. on-going aircraft are considered as constraints for the active aircraft. Because the entry time of aircraft in TMA is a decision variable, we have to take into account the uncertainty of t E i and t L i of each aircraft, to dynamically define the status of aircraft.

Tailored SA algorithm

SA originated in the domain of thermodynamics (Kirkpatrick u. a., 1983;[START_REF] Černỳ | Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm[END_REF]. It is well known for its ability to escape from the local minima by allowing random neighbourhood changes. Moreover, it can easily be adapted to different kinds of problems with continuous or discrete space states. This method stems from an analogy with the physical phenomenon of slow cooling found in metal in a state of fusion which leads to a solid, low-energy state. The temperature must be reduced slowly, with steps which are sufficiently long for thermodynamic equilibrium to be attained at each temperature level. For materials, this low energy results in a regular, crystal-like atomic structure.

The annealing process thus consists of bringing a solid into a low energy state after raising its temperature, a process which may be summarized in the following two steps: 1) raise the solid to a very high temperature in order to reach the point of "fusion"; 2) cool the solid, following a specific temperature reduction plan in order to attain a solid state with minimal energy.

In 1953, Metropolis developed an algorithm to simulate the physical process of annealing on a computer (Metropolis u. a., 1953). Given a current state i of energy E i , a state is generated by applying a disturbance which transforms the current state into a new state.

1. If E j -E i ≤ 0 the state j is accepted as the new current state. 2. If E j -E i > 0 the state j is accepted as the new current state with probability P a :

P a = e E i -E j k b T .
with T : temperature and k b : Boltzmann constant.

The temperature influences the probability of accepting of a higher energy state. For a high temperature, the probability of accepting any given movement tends towards 1: all changes will be accepted. If the cooling process is sufficiently slow, the solid attains a state of equilibrium at each temperature. In the Metropolis algorithm, this equilibrium is attained by generating a large number of transitions at each temperature. Thermal equilibrium is characterized by the Boltzmann statistical distribution. This distribution gives the probability that the solid will be in a state i of energy E i at temperature T :

P r {X = i} = 1 Z(T ) e - E i k b T ,
where X is random variable associated with the current state of the solid; and Z(T ) is the distribution function of X allowing normalization:

Z(T ) = j∈S e - E j k b T .
In the SA algorithm, we can apply the Metropolis algorithm to generate a sequence of solutions in the state space S. To do this, we create an analogy between a multi-particular system and our optimization problem using the following equivalences: 1) the admissible solutions represent the possible states of the solid; 2) the function to be optimized represents the energy of the solid. We then introduce a control parameter, Cl, which plays the role of the temperature. Let Cl k be the value of this parameter and G k the number of transitions generated at iteration k. Using this notation, we can summarize the principle of SA in Algorithm 1. At the beginning of the process, the values of Cl k are high, allowing us to accept transitions with major degradations in the objective function, thus exploring the state space in a homogeneous manner. As Cl k decreases, only transitions which improve or barely damage the criterion are accepted. Finally, as Cl k tends towards zero, no deterioration of the criterion will be accepted, and the SA algorithm behaves in the same way as a local search algorithm.

Algorithm 1 Simulated Annealing

Require: x i , Cl 0 , G 0 , k = 0, Cl EN D 1: repeat 2: for g = G 0 → G k do 3:
Generate a solution x j from the neighbourhood S( x i ) of the current solution x i ;

4:

If f ( x j ) < f ( x i ) then x j becomes the current solution;

5:

Otherwise x j becomes the current solution with probability p = e

f ( x i )-f ( x j ) Cl k 6:
end for 7:

k = k + 1 8: Calculate (G k , Cl k ) 9: until Cl k Cl EN D

Application of the optimization method

For our optimization problem, at a given time t in the process of RHC, we first identify the on-going and active aircraft by the sliding window management algorithm 2. We then put them in the decision vector ( u i ) t , see Tab.3.3. We add for each aircraft a performance value which represent how the performance of decision on such aircraft. The performance value will be served for guiding the optimization algorithm to search the near-optimal solution. Such values are gathered together into a vector

Q t = (Q t 1 , Q t 2 , ..., Q t N )
. At the beginning of the optimization process, this vector Q t is set to 0. The decision vector ( u i ) t is then given to the simulation module which generates the trajectories of aircraft, computes the objective z t , and updates the performance Q t i for each trajectory. z t and Q t i are used by the SA algorithm to find a near-optimal solution for the objective defined by function 3.22. The default values of the parameters in the RHC-SA algorithm can be found in Tab. 3.2.

Algorithm 2 Sliding Window Management

1: procedure SlidingWindow 2: k ← 0; 3: T s (k) ← t INIT ; 4: T e (k) ← T s (k) + W; 5:
Determine each flight status relative to sub-window; end while 17: end procedure More precisely, the neighbourhood operation helps the SA to explore new areas in the state space domain. In order to speed up the algorithm, this operation selects more often the decisions with a bad performance Q t i , and it is only executed among the active aircraft in a given window. Once a decision vector has been selected, we have to change its associated decision variables in order to obtain a new decision vector for generating a new trajectory. The default strategy for choosing the combination of j, g, h and r follows the rule of Roulettewheel selection using weighted probabilities. The highest desirable decision variable has the largest share of the roulette wheel, while the lowest desirable decision variable has the smallest share. The process is described in Algorithm 3. The default weights for the four decision variables t E i , t T i , v i , r i are: 25, 25, 25, 25. The control parameter β is designed to control the probability of choosing runway 18R-36L as a landing runway, i.e the value of r i . For example, we randomly generate a number which is in the range [0,[START_REF] Lee Und Balakrishnan ; Lee | Hamsa: A study of tradeoffs in scheduling terminal-area operations[END_REF]. If this number is less than β, then r i = 1, and the aircraft chooses runway 18R-36L to land. Otherwise, r i = 0, and the aircraft chooses runway 01-19 to land. The default value of β is 0.5. With Algorithm 3, we generate the new decision variable. Then, the former decision vector is replaced by the new decision

O A O A A A A Decision ( u i ) t ( u 1 ) t ( u 2 ) t ( u 3 ) t ... ... ... ( u N ) t Performance indicator Q t i Q t 1 Q t 2 Q t 3 ... ... ... Q t N Algorithm 3
The strategy for choosing the combination of j, g, h and r Require: weights (w 1 , w 2 , w 3 , w 4 ) for j, g, h and r w 1 , w 2 , w 3 , w 4 ∈ [0, 100] and w 1 + w 2 + w 3 + w 4 = 100 1: for i ← 1, 4 do 2:

Randomly select an integer I I ∈ [0, 100]

3: if I <= w 1 then 4:
Change the value of j Update the entry time end if 18: end for vector, consequently the new objective value z t and performance field (Q t 1 , Q t 2 , ..., Q t N ) are computed thanks to the simulation module. The neighbourhood selection process is repeated again until the temperature in SA reaches its minimum value.

Conclusion

In this chapter, we first present an overview of optimization methodologies. Four main components are listed: route network, mathematical formulation, optimization module, and simulation module.

Then, we build up a MLMPMS in the part of route network design, after making a systematic analysis of the strategy of sequencing, merging, descending in the TMA. This route network is the innovation of the PMS system for parallel runway operation. It can conveniently realize a near-CDA operation.

After that, we analyze mathematical formulations for our optimization problem, including assumptions, state space, constraints, objectives. For each aircraft, six kinds of data are given. We define a multi-objective optimization objective function with four decision variables and its associated constraints on operation and safety. RHC-based approach is used to address our dynamic optimization problem. Tailored SA algorithm is designed to solve the optimization problem. It searches decision variables to match the objective requirements, such as conflictfree, minimizing delay, etc.

In the next chapter, we will build up the simulation module, which compute a realistic 4D trajectory and the performances of decision variables. Its results will be used to support optimization module to find a good solution.

Chapter 4

Dynamic system based simulation modelling

The optimization module generates the decision variables which then are used by the simulation module to generate the trajectory, and the associated performance indicators(C, D, S, and P ). This chapter explains modelling method used in the simulation module, including continuous dynamics modelling, discrete dynamics modelling, and conflict detection modelling, see Fig. 4.1. The part concerning continuous dynamics consists of an equation of motion for an aircraft, airlines scheduled CAS speed computations, and the fuel consumption computation. The discrete dynamics modelling consists of flight procedure modelling and multi-phase trajectory generation. All the computations in the simulation module are supported by ISA, BADA, and topology design.

Discrete dynamics

Flight procedure modelling

The flight procedure model is closely related to the topology design, since it depends on the operational model at the airport. The main goal of flight procedure modelling is to support the robustness of ML-PM system, while at the same time reducing the complexity of the traffic.

Taking Fig. 4.2 as an example, we model the 3D flight procedures L by a 3D network of

V = w i , i ∈ N, (4.1)
E = e i , i ∈ N, (4.2)
H = h i , i ∈ N, (4.3) R = (V, E) = R i , i ∈ N, (4.4)
L = (R, H) (4.5)
where V is a set of way-points, including both fixed way-points and the dynamic way-points which corresponds to the turning points on the sequencing leg. E is a set of edges, including the fixed links and the dynamic links. Dynamic links correspond to the following links: the sequencing leg and the link between the turning point and the merge point. All the other links are considered as fixed links. R is the available routes for aircraft, and H is a set of target altitudes at the associated way-points.

The dynamics in the merging zone in the ML-PM topology is the most difficult part of flight procedure design. For the equal-distance segregated ML-PM topology, the solution is simple, because it consists of two independent single-runway PMSs. Each merging zone is independent. There is no crossing traffic between them. While, for mixed parallel operation patterns, the ML-PM topology is more complex. Consider Fig. 4.2 as an example, after being transferred by adjacent sectors, aircraft will follow different arrival routes to corresponding sequencing legs, the blue and red legs in Fig. 4.2. On these sequencing legs, aircraft maintain level-off and keep lateral separation on different legs, then they are sequenced to turn toward the appropriate merge point one by one, according to operational requirements. The advantage of this topology design is that aircraft can easily change landing runway if some constraints arise, such as unavailability of the initial landing runway, weather, or military control. The challenges are: first, compared with traditional PMS, here the distances from each point on the sequencing legs to the corresponding merge points (M1 or M2) are not equal. And second, we have to avoid intersections of trajectories in the merging areas.

Finally, the flight operation for the merging zones in the mixed parallel operation pattern is designed as shown in Fig. 4.3. First, different segregated flight levels on the sequencing legs are provided for incoming aircraft according to their wake turbulence categories. "Heavy" aircraft will use the higher level, "Medium" aircraft will use the middle level and "Light" aircraft will enter the lower level. Once aircraft enter the sequencing leg, they fly at the same constant CAS, then they will perform a "Direct to" turn towards the merge point. Aircraft are only allowed to take a turn in the allowable range on the sequencing leg for M1 or M2. More precisely, if aircraft a want to merge to M1, it can only turn between P s and P m (the midpoint of the whole sequencing leg). If it wants to merge to M2, first it has to follow P s to P m , then it can make a turn to M2 between P m and P e . After they have turned, aircraft perform a near-optimal continuous descent during the merging process.

Multi-phase trajectory control

Commercial aircraft follow a flight plan from the departure airport to the destination airport. The flight phases from descent to landing could be further sub-divided into four segments: arrival descent, initial approach segment, intermediate approach segment, and final approach segment. In the ML-PM system design, aircraft first enter TMA at an entry point with a target flight level and CAS, they then follow an arrival route to the entry point of the sequencing leg at another target flight level. They keep a constant CAS and a level-off flight on the sequencing leg until the time they turn toward the merging point. They then descend to another target flight level at the merging point, reach Final Approach Fix (FAF) and finally land on the runway. In summary, there are five different designed phases during this whole process of approach and landing. At each phase, an aircraft has a specific configuration to control its flight status. We model this multiphase flight trajectory control problem as a deterministic finite-state machine M . Definition 1 A finite state machine M is defined by a 5-tuple (Σ, Q, q 0 , F, δ), where 1. Σ is the set of symbols representing input to M,

2. Q is the set of states of M, 3. q 0 ∈ Q is the start state of M , 4. F ⊆ Q is a set of final states of M, 5. δ : Q × Σ → Q is the transition function.
F = {q 9 , q 10 } is the final state. Each state q n has an associated operational procedure (OP), aerodynamic configuration (AC) and dynamic motion (AM). The transition function is given by: δ(q 0 , l 1 ) = q 9 (4.6) δ(q 0 , l 2 ) = q 10 (4.7) For each state of M, we calculate the aircraft trajectory based on the BADA 3.13 aircraft performance model. The Total Energy Model (TEM) is used to predict aircraft trajectories in BADA. In the TEM, any two of the three variables of thrust, speed (CAS), or Rate of Descent (ROD) are controlled, the third variable is then determined as a function of the first two. In our simulation module, both ROD and CAS are controlled, then the necessary thrust and fuel flow are calculated. Because the ROD and descent angle in each route segment are fixed, and CAS is constant, the associated True Air Speed (TAS) changes linearly with flight altitude. In this thesis, to simplify our research problem, we take the A320 BADA performance parameters as representative of all aircraft with "Medium" wake turbulence category, the B747-800 for "Heavy", and the CESSNA for "Light". Algorithm 4 describes the entire trajectory computation process.

l 1 = {x ∈

Algorithm 4 Multiphase trajectory generation

1: procedure Multiphase trajectory control 2: get the decision vector x i for aircraft i; 3: set the initial state q 0 of aircraft i; for each segment u i ∈ E do 6: verify the status of aircraft i in δ(q 0 , l 1 ) or δ(q 0 , l 2 ) 7:

get the altitudes required at the beginning and ending way-points;

Continuous dynamics 4.2.1 Aircraft equations of motion

For each state of M, the motion of an aircraft follows a continuous dynamics rule. The model, illustrated in Fig. 4.5, describes the forces applied to the centre of gravity of a descending aircraft. We assume that aircraft operate in ISA and in no wind conditions, that earth is a sphere of radius R T , and the g is constant. An aircraft's state is denoted by X = [x, y, h, V a , Ψ a , m], where (x, y, h) is the 3D position of the aircraft, V a is the true airspeed, Ψ a is the horizontal path angle, and m is the mass of the aircraft. The command vector is denoted by U = [T hr, γ a , Φ], where T hr is the total thrust, γ a is the vertical flight path angle, and Φ is the bank angle. The motion of a descending aircraft satisfies the following set of equations [START_REF] Nuic | [END_REF]: Where D is the aerodynamic drag, and L is the aerodynamic lift. To solve this problem, we have to calculate X.

ẋ = V a cos γ a cos Ψ a , (4.10) ẏ = V a cos γ a sin Ψ a , (4.11) ḣ = -V a sin γ a , ( 4
First, V a is modelled as a function of CAS, because it is common for commercial aircraft to descend with a constant CAS value, see Section 4.2.2. With a given CAS, which will be further explained in Section 4.2.2, we can calculate V a using Equation 4.16 [START_REF] Nuic | [END_REF],

V a =   2 µ p ρ      1 + p 0 p   1 + µ 2 ρ 0 p 0 CAS 2 1 µ -1     µ -1      1 2 , (4.16)
where, p is the air pressure, p 0 is the air pressure at mean sea level, ρ is air density, ρ 0 is the air density at mean sea level, µ = 1 3.5 . Lift and drag are calculated using the following equations [START_REF] Nuic | [END_REF]: (4.18) where C L is the lift coefficient, C D is the drag coefficient, and S is the wing reference area. Under normal conditions, C D is specified as a function of C L as follows [START_REF] Nuic | [END_REF]:

L = 1 2 ρSV 2 a C L , (4.17) D = 1 2 ρSV 2 a C D ,
C D = C D0 + C D2 × (C L ) 2 . (4.19)
With different AC values, i.e cruise, approach and landing configuration, the values of C D0 and C D2 change. In the BADA Operations Performance File (OPF), we can get the values of m, g, S, C L , C D0 and C D2 . If we fix the values of γ a , and Ψ a in each route segment, referring to the route network and flight altitude requirements in Section 4.1.1, then we can calculate L, and D using Equations 4.14,4.15,4.17,4.18,and 4.19,and calculate (x,y,h) using Equations 4.10,4.11,and 4.12. Last, we need to find a way to compute T hr. Equation 4.13 is also called TEM in BADA, which is used to predict aircraft trajectories. Recall that in TEM, any two of the three variables of thrust, speed, or ROD are controlled, the third variable is the function of the first two. Here, both ROD and CAS are controlled, and the necessary thrust is calculated. Equation 4.13 can be rewritten as [START_REF] Nuic | [END_REF]:

ROD = dh dt = (T hr -D)V a mg 1 + V a g dV a dh -1 . (4.20)
We introduce an energy share factor as a function of Mach number, labelled by f {M } [START_REF] Nuic | [END_REF]:

f {M } = 1 + V a g dV a dh -1 , (4.21)
and it exits the following function in a condition of constant CAS descent below the (4.23) where k = 1.4, R is the real gas constant for air, β is the ISA temperature gradient with altitude below the tropopause, and T is the air temperature. The descent thrust T hr can then be calculated by:

f{M } =    1 + kRβ 2g M 2 + 1 + k -1 2 M 2 -1 k-1    1 + k -1 2 M 2 k k-1 -1       -1 , (4.22) M = V a / √ kRT ,
T hr = ROD f{M } × mg V a + D. (4.24)

Airline CAS schedule (descent)

The CAS schedule for jet and turboprop aircraft, is given in Tab.4.1. For different flight level ranges, CAS is calculated using a different function, with some specific parameters and the landing stall speed.

The meaning of the variables are as follows: C V min is the minimum speed coefficient, (V stall ) LD is the landing stall speed, V des,1 is the standard descent CAS [knots] between 3000/6000 and 10000 ft, and V des,2 is standard descent CAS [knots] between 10000 ft and Mach transition altitude. M des is the standard descent Mach number above the Mach transition altitude. All the parameters mentioned here can be found in the BADA PTF, APF and OPF files. Note that, V des,2 equals v E i in the simulation module of this thesis.

Fuel consumption computation

For the jet and turboprop engines, the thrust specific fuel consumption, η [kg/(min.kN)], is specified as a function of the TAS, V a [kt] [START_REF] Nuic | [END_REF]:

η = C f 1 × (1 + V a C f 2 ). (4.25)
Hence, the nominal fuel flow, f nom [kg/min], can be calculated as a function of the thrust, T hr [START_REF] Nuic | [END_REF]:

f nom = η × T hr (4.26)
These expressions are used in all flight phases except during idle descent and cruise, where the following expressions are used. The minimum fuel flow, f min [kg/min], corresponding to idle thrust descent conditions for jet engines, is specified as a function of the geopotential pressure altitude, H p [ft] (here, H p = h), that is [START_REF] Nuic | [END_REF]:

f min = C f 3 × (1 - H p C f 4 ) (4.27)
Note that for jet engines, the idle thrust part of the descent stops when the aircraft switches to approach and landing configuration, at which point thrust is generally increasing. Hence, the calculation of fuel flow during the approach and landing phases shall be based on the nominal fuel flow and limited to the minimum fuel flow if necessary [START_REF] Nuic | [END_REF]:

f ap/ld = M AX(f nom , f min ) (4.28)
For different aircraft configuration (cruising, approach and landing), given the values of T hr, h, V a and the coefficient parameters

C f 1 , C f 2 , C f 3 , C f 4
found in the BADA APF and OPF files, we can calculate the associated fuel flow.

Conflict detection modelling

M i contains all merge conflicts within a merging zone. There are conflicts between aircraft approaching the same merge point M1 or M2.

If we consider the flight plan of an aircraft as shown in Fig. 4.6, the conflict C i of aircraft i, i ∈ F is calculated as below: (4.29) (4.30) (4.31)

C i = L i + N i + M i ,
L i = L1 + L2,
N i = N 1 + N 2 + N 3 + N 5,
M i = L3 + L4 + N 4. (4.32)
The link conflict L i is determined based on the time difference at the beginning and the end of the arc. Assuming that leading aircraft, i, and trailing aircraft, j, will enter into the same arc u ∈ E defined by way-points w a and w b , we define ∆t u i,j as the entry time difference and ∆t u i,j as the exit time difference, see Fig. The node conflict N i is determined based on the time difference when passing a common way-point. Assuming that the leading aircraft, i, and the following aircraft, j, will pass a common way point w c ∈ V , see Fig. 4.8, then if the following constraint could not be satisfied, then N i increases by one:

∆t wc i,j := t wc j -t wc i ≥ s min i,j | j, i ∈ F. (4.36)
The calculation of M i is more complicated, see Fig. 4.9. It depends on the type of route network topology and the flight procedure. For convenience, we transfer the route network topology of both merging zones into a virtual time-based network to detect the merge conflicts. All the links are time-based segments. In the following context, we will explain the calculation of M i in details.

Turning points equidistant from a merge point

For single-runway and segregated ML-PM topologies, the calculation of M i is simple, because there is just one available merge point for each sequencing leg, and the distances between each point on the sequencing leg and merge point are the same, see Fig. Here, b and c are not fixed, because t T i for each aircraft is different. Based on this concept, if aircraft can not comply with the following constraints, then M i increases by one:

∆t T i,j := t T j -t T i ≥ s min i,j . (4.37)
where t T j and t T i are the decision variables controlled by the optimization module. The radius d 1 and d 2 are constant values defined in the ML-PM topology design. Therefore, the passing time over the merge point M can be calculated. And the times over the merge point M for two successive aircraft i and j (supposing that i is turning from point b, j is turning from point c ) should satisfy the following constraints:

t M j = t T j + 2 × d 1 v T j + v M j , ( 4.38
)

t M i = t T i + 2 × d 2 v T i + v M i , (4.39) ∆t M i,j := t M j -t M i ≥ s min i,j . (4.40)
where v T i and v T j are speeds on the turning points, v M i and v M j are speeds on the merging points. They are computed according to the airline CAS speed schedule and BADA performance files. We suppose that the acceleration is constant, because the rate of descent (or climb) is constant in our case.

Turning points not-equidistant from a merge point

For the mixed ML-PM topology, the calculation of M i is more difficult, because there are two possible merge points for each sequencing leg, and the distances between each point on the sequencing leg and the chosen merge point are different. As shown in Fig. 4.11, for an aircraft entering from West, if its r i = 1 it will dynamically turn at some position on the sequencing leg (the red color), then approach M 1 to land. If its r i = 0, it will travel over mp2 first then it will dynamically turn at some position on the sequencing leg toward M 2 to land. The same procedure holds for the aircraft entering from East. Based on this, if aircraft can not comply with the following constraints at the merging points, then M i increases by one:

∆t T i,j := t T j -t T i ≥ s min i,j (r i × r j + (1 -r i )(1 -r j )). (4.41)
Next, we have to compute M i on M 1 or M 2 respectively. To calculate the times over M 1 or M 2, we have to calculate the distance between the dynamical turning point and the merging point. We assume that there are two aircraft i and j arriving on different sequencing legs.

Aircraft i flies on the outer sequencing leg, and aircraft j flies on the inner sequencing leg. Given the coordinate (x a , y a ) of the centre point a, see Fig. 4.11, the coordinates of b (x 3 , y 3 ) and c (x 4 , y 4 ) can be calculated by:

θ 2 = - v i (t T i -t T imin ) d 2 + θ imin (r i ) (4.42) θ 1 = v j (t T j -t T jmin ) d 1
+ θ jmin (r j ) (4.43)

x 3 = d 2 × cos θ 2 + x a (4.44)
y 3 = d 2 × sin θ 2 + y a (4.45) x 4 = d 1 × cos θ 1 + x a (4.46) y 4 = d 1 × sin θ 1 + y a (4.47)
where θ imin corresponds to t T imin , and θ jmin corresponds to t T jmin . Note that θ imin and θ jmin depend on r i and r j respectively. If an aircraft does not need to change runway, then its earliest turning time is the entry of the sequencing leg, and its latest turning time is the midpoint of the sequencing leg; if it needs to change runway, then its earliest turning time is the midpoint of its sequencing leg, and its latest turning time is the end of the sequencing leg. Consequently, l 1 , l 2 , l 3 , l 4 can be calculated based on b(x 3 , y 3 ), c(x 4 , y 4 ). The time difference of aircraft i and j passing the associated merge point can also be calculated.

t M j = t T j + 2 × l 1 r j + l 2 (1 -r j ) v T j + v M j , (4.48) t M i = t T i + 2 × l 3 r i + l 4 (1 -r i ) v T i + v M i , (4.49) ∆t M i,j := t M j -t M i ≥ s min i,j . (4.50)
where v T i and v T j are speeds on the turning points, v M i and v M j are speeds on the merging points. We suppose that the acceleration is constant, because the rate of descent (or climb) is constant in our case.

Conclusion

In this chapter, we explain all necessary modelling approaches to compute a more realistic 4D trajectory in the simulation module. The computations concern airspeed, aerodynamic forces, motion of aircraft, as well as fuel consumption, flying distance, flying time (or transit time), etc.. The approach to compute the trajectory is based on discrete dynamic system control theory. We model the flight procedure in 3D with a set of nodes, links and altitudes. Then we model the multiphase flight trajectory control problem as a deterministic finite state machine. The state of aircraft is progressively switched from entry point to land on the runway. The state of aircraft in each phase is computed based on BADA performance model.

After aircraft trajectories generation, we detect the conflicts between aircraft, so as to help our optimization module to find a conflict-free solution. The conflict detection modelling is divided into three parts: link conflict, node conflict and merge conflict. The total conflict is the sum of them. Geometric analysis for each type of conflict detection is made. Specially, for the merge conflict, two more cases are further analyzed, they are: turning points equidistant from merge point and turning points not-equidistant from merging point. In the next chapter, we will present the application of our algorithm coding in JAVA to realistic instances of the problem. West, after take-off, use the outer sequencing leg to join one of the exit points KM, SOSDI and RENOB. Flights departing to North-est, after take-off, use the inner sequencing leg to join one of the exit points TONIL, CDY, and YV. Flights departing to South, after take-off, join the middle point of sequencing leg, then fly towards the exit point LADIX. For all the available departing routes see Appendix A.5.

If we combine the departure ML-PM with the arrival ML-PM, see Fig. 5.3, we can see clearly that there are three hot spots of intersecting traffic. For safety, in these hot spots we design a strict flight level constraint for aircraft passing through these areas. For example, in South-bound operations, in hot spot 1, arrival aircraft can not fly lower than 4500 m. In hot spot 2, arrival aircraft can not be lower than 4200 m. In hot spot 3, departing aircraft can not be higher than 3900 m. The same design principle is applied for the North-bound operation; in hot spot 1, arriving aircraft from KM maintain 4500 m until they have passed the boundary. With this kind of topology design, we can simplify the crossing trajectories, and ensure safety.

Taking QFU 18 as an example, we design the total optimized descent profile for BCIA using the following steps. Firstly, the flight altitudes on the five entry points of TMA are kept as real operational values in reference to the control transfer agreement between adjacent sectors. Once an aircraft enters the TMA, we keep it as high as possible before initiating descent. Then, according to the distance from the runway threshold and the descent perfor- ). The points with the same colour are designed for a specific route. In total there are six different arrival routes to guide an aircraft to land. We can see that most of the target altitudes are below, but not far from the optimal 3 • descent profile. Thus, they will not exceed aircraft engine limitations. In particular, we design two levels to separate "Heavy" and "Medium" aircraft. At W4, the entry point of outer sequencing legs, 2400 m is for "Medium", 2700 m is for "Heavy". At W6, the entry point of inner sequencing legs, 3000 m is for "Medium", 3600 m is for "Heavy". Remark that there are no Light aircraft operating at BCIA in the historical flight data, hence we do not need to design a flight level for them. Similarly, there are few A380 landing at BCIA, hence it is not necessary to design a separate level for them. Instead, they can join the level for "Heavy" aircraft, if necessary. Furthermore, the altitudes at M1 and M2 are 1500 m and 1200 m respectively. Finally, with these target altitudes at different significant way-points, aircraft from different arrival routes (JB-A, BOBAK-A, KM-A, GITUM-A, VYK-A, DOGAR-A) can maintain a higher altitude as long as possible, then they can execute a near-optimal descent profile to land on the allocated runway. For all the target altitudes at the way-points see Appendix A.5.

Mixed parallel approach operations

Data preparation

The first dataset Data set 1 is built up based on the flight plan of flights on BCIA on a specific operational day in December, 2015. As shown in Tab. 5.1, there are 670 flights planned to land at BCIA in 24 hours, 76.27% of which are "Medium", and 23.58% of which are "Heavy". The geographic distribution of the traffic is as follows: 14.03% of traffic comes from KM, 15.97% from JB, 12.24% from BOBAK, 18.66% from VYK, 18.66% from DOGAR, and 20.45% from GITUM. In total, 57.76% of flights are planned to land on runway 01-19, and 42.24% of flights are planned to land on runway 18L-36R. There is no planned switch of runway-in-use.The hourly distribution of flight data in the 24 hour period is shown in Fig. 5.5. We focus on the period between 6:00 to 20:00 at BCIA, because outside of this period the traffic density is too low.

Performance of the automated conflict resolution

We chose the peak period (6:00-8:00) to evaluate the performance of the proposed algorithm in regards to automated conflict resolution. In Fig. 5.6, the X-axis is the passing time with 6 sliding windows, the primary Y-axis is the number of conflicts and the secondary Y-axis is the number of "Active" aircraft in the current window. The blue line represents the "Link conflicts", the orange line represents "Node conflicts", and the grey line represents "Total conflicts". The yellow column represents the number of aircraft (A/C for short). We found that within each sliding window the algorithm can quickly find a conflict-free solution, and that 100% of conflicts can be successfully resolved. In the third window, there are 48 "Active" aircraft and the initial number of conflicts is 600. In the sixth window, there are 10 "Active" aircraft and the initial number of conflicts is close to 70. We also found that the running time of the algorithm increased with the number of the "Active" aircraft in the current window and the number of initial conflicts.

Second, we compared the de-conflict performance of the four decision variables, and ranked them in order of influence. This also provided insight into how the interactions between decision variables influence the de-conflict performance of our system. Different decision variables have different impacts on the performance of our algorithm. Moreover, different actors, such as airport, controllers, airlines, would prefer to take different decisions to solve the conflicts based on their requirements and real traffic situations. If an aircraft is still far from the TMA, then it may be more interesting to change the entry slot time or the planned runway than turning control on the sequencing leg or speed control on the entry point. If an aircraft is close to the TMA, then turning control on the sequencing leg and speed control on the entry point may be preferable. This reference is realized by assigning different combinations of weights to the decision variables. Taking 2 hours (6:00-8:00) of Data set 1 with the highest demand as input, we defined 8 sets of weights for x i (t E i , t T i , v i , r i ): G1 (45-25-15-15), G2(35-35-25-5), G3(10-30-20-40), G4(40-30-20-10), G5(0-0-0-100), G6(0-0-100-0), G7(0-100-0-0), G8(100-0-0-0) to study the impacts of different decision preferences. The results are shown in Fig. 5.8. The de-conflict performance of each individual variable is separately demonstrated by G5, G6, G7, and G8. We found that the runway assignment r i has the weakest de-conflict performance (Test G5), which is normal due to the fact that in this case the state space dimension is strongly reduced. The next weakest performance was the turning time control on the sequencing leg t T i (Test G7), followed by the entry speed v i (Test G6). The strongest de-conflict variable is the entry time t E i (Test G8). It alone can solve a large number of the conflicts in one sliding window. Tests G1, G2, G3, and G4 are the mixed combination of decision variables. They all have a good conflict resolution performance, 100% of conflicts can be resolved.

Third, we further test the algorithm de-conflict performance by increasing the scenario Figure 5.6: De-conflict performance in period (6:00-8:00) Figure 5.7: De-conflict performance in period (6:00-12:00) time period from 2 hours (6:00-8:00) to 6 hours (6:00-12:00), and by augmenting the size of the sliding window W from 3600 s to 7200 s. Note, the time S shifts from 1800 s to 3600 s. There are 286 aircraft in this 6-hour period. The results are shown in Fig. 5.7. We can see that the total number of conflicts in one sliding window increases with the size of the sliding window. For example, in the first sliding window, there are more than 1000 initial conflicts. The more "Active" aircraft in one sliding window, the more initial conflicts to be solved. The de-conflict performance of our algorithm is stable, 100% conflicts can still be successfully solved with the same parameter settings as in the SA. However, it took around 700 s to optimize 6 hours traffic, using substantially more time to find a solution.

Fourth, we compared the solution quality of SA with the solution quality of HC. As shown in Tab. 5.2, the quality of SA is better than the quality of HC. SA is always able to compute the conflict-free solution. HC can not guarantee a conflict-free solution, and Figure 5.8: De-conflict performance with different combinations of weights on decision variables produces solutions with higher variance in quality than SA. However, HC is faster than SA.

Trajectory control

We studied the aircraft trajectory under the default weighting scheme (25-25-25-25) and input data from the period (6:00-8:00). First, we looked at the flight performance of individual aircraft. As an example, we consider the Medium aircraft labelled A3245 coming from JB with entry CAS 300 kts. It follows the route "JB-W1-W2-W3-W4-Sequencing leg-M1-Rwy(18R-36L)". As shown in Fig. 5.9a, this aircraft makes a CAS constant descent, and the associated TAS is linearly changing with the flight altitude. Its fuel flow profile is shown in Fig. 5.9b. At high altitude, the rate of fuel flow is relatively low, then, as the aircraft is approaching the threshold of runway with a fixed angle of descent, it increases rapidly. The descent profile of A3245 with distance is shown in Fig. 5.9c. After entering the TMA, it maintains the entry altitude as long as possible, and then executes a continuous descent. Because there is no level-off on the sequencing leg, the total descent profile has no disruption. The distance travelled profile can also be seen in Fig. 5.9d.

Next, let us look at the results for all trajectories. Fig. 5.10a is the vertical view of all the trajectories. Aircraft travel time is calculated from entry into the TMA to the runway threshold (in seconds), the altitude of aircraft is measured in 100 feet. A group of trajectories with the same colour are aircraft from the same entry point into the TMA. In total, there are 6 groups: KM, JB, BOBAK, VYK, DOGAR and GITUM. We can see that there are four level-off lines between 8000 ft to 12000 ft. The two higher lines are for the inner sequencing legs, and the two lower lines are for the outer sequencing legs. These four level-off lines mean that aircraft are on the sequencing legs. The level-off time on the sequencing leg is different for each aircraft, and depends on the decision variable t T i . The maximum level-off time is 5 minutes, the minimum level-off time is 0 minute. The lowest level-off sequencing legs for aircraft (FL80) is still above the designed flight level at which the full approach or landing configuration is triggered (FL30 for A3245, see Fig. 5.9b). The fuel flow rate on the sequencing legs is still relatively low. Therefore, the additional fuel consumption due to re-sequencing and runway allocation is still acceptable. Fig. 5.10d shows the distribution of final decisions taken to compute all the trajectories in this scenario. We can find the range of deviations for three decision variables: the entry speed, the entry time (slot change), and the turning time. The percentage of aircraft landing on runway 18R-36L is the positive value at the area of "Runway allocation", and the percentage of aircraft landing on runway 01-19 is the negative value. The average speed change is 2%, the average entry time change is -4%, and the average holding time on the sequencing leg is 28%. Fig. 5.10b is the static horizontal view of all trajectories in 2D. The left part of Fig. 5.10b is the real radar-based trajectories, the right part of Fig. 5.10b is the simulated trajectories proposed by our system. We observe that the proposed merging trajectories (right half) in the two merging zones are very structured, compared with real traffic (left half). Aircraft can change their initial landing runway, however, there are no crossing trajectories between two merging zones, thanks to the tailored procedure designed for runway re-assignment operation. Fig. 5.10c is the static view of all trajectories in 3D, which again are very structured. Last, the average fuel consumption and travel time of the aircraft with different categories can be also evaluated. The results are listed in Table 5.3. Based on the different entry waypoints, we can get the information on the average fuel consumption [kg] and the average flying time in the TMA [sec(min)] for different categories of aircraft, "Heavy" or "Medium". For example, a "Heavy" aircraft entering from KM has an average fuel consumption of 948 kg, and an average travel time of 13 minutes. While, a "Medium" aircraft entering from the same entry point consumes an average fuel of 136 kg and spends an average flying time of 14 minutes to land. KM and GITUM are less far from the runways, hence the average flying time and fuel consumption on the routes from these two entry points are relatively lower than from other entry points. This information related to the optimized trajectories of aircraft can be shown directly to the flow managers to assist them in decision making. These information related to the optimized trajectories of aircraft can be shown directly to the flow managers to assist them in decision making. Compared to the real traffic in the left half of Fig. 5.10b, the advantage of our system is clear. Further, up to 8 minutes of flying time can be saved when our system is used to design the flight trajectories.

Re-sequence ability

In order to study the re-sequencing ability of the proposed ML-PM system, we set the weighting scheme for decision variables to be (10-50-40-0). We take input data from the period (6:00-8:00). We focus on the re-sequencing induced mainly by controlling the speed change and the turning time on sequencing legs, in the case when there is no runway changes for aircraft and very few allowable deviations on the entry time to the TMA. We executed ten tests, all of which have conflict-free results. The associated MPS and CPS are shown in Tab.5.4. We find that most aircraft keep less than or equal to 3 position shifts, with an average value of 76.23% on runway 01-19 and 80.45% on runway 18R-36L. On average, 23.77% of aircraft execute a CPS > 3 deviation on runway 01-19, 19.55% of aircraft on runway 18R-36L respectively. The MPS of aircraft to land on runway 01-19 is around 6, and on runway 18R-36L it is 6. We also find that in some individual tests the value of MPS can reach up to 7 for runway 01-19 and 8 for 18R-36L.

An example of the sequence results can be seen in Fig. 5.11 in the part of "MPS > 3". This kind of landing queue position shifting needs controller to issue only one instruction "Direct to" or "speed regulation", compared with traditional methods of radar vectoring with several heading changes, it requires less radio communications between the controller and the pilot. Compared with the conventional way "MPS <= 3 (less deviation from ETA)" in Fig. 5.11, the makespan of the case "MPS > 3" is shorter than the case "MPS <= 3", and the landing times of the last aircraft on the runway are generally earlier. However, the associated average fuel consumption per flight and the average flying time in the TMA per flight in the case "MPS > 3" is a little higher than that in the case "MPS <= 3". These results are listed in Tab. 5.5.

Fuel saving with runway landing balancing

We simulated four scenarios to investigate the advantage of balancing landing rates on asymmetric arrival flows on two runways. In each test, there were 128 aircraft in a fixed time window of 2 hours. As shown in Tab. 5.6, S1 and S2 have the same input traffic data, they represent the extra unequal initial arrivals on the two runways with 69% demand for runway 18R-36L and 31% demand for runway 01-19. S3 and S4 have the same input traffic data. They represent the normal initial arrivals on the two runways with 52% demand for runway 18R-36L and 48% demand for runway 01-19. In S1 and S3, aircraft can not change runways, and can not change the initial turning time on the sequencing leg. They have the same combination of weights (50-0-50-0), which give more chance to speed change and entry time change to resolve conflicts. They simulate the conventional way of de-conflict used by controller. In S2 and S4, there are no limitation on the decision variables. They use the weighting scheme (25-25-25-25), and represent the optimized way proposed by our algorithm. In S1, due to super dense arrival flows at one runway, a lot of aircraft are delayed and hence the value of D is out of the default range [0,100]. As a result, our system could not find a conflict-free resolution. In contrast, S2 can find a conflict-free solution. Compared with S1, S2 achieved an average reduction in fuel consumption of 1.0% per aircraft, and saved, on average, 1 minute of flying time in the TMA. Further, the average "D" in S2 is much lower (less delay), and "S" is shorter (landing is faster). S3 found a conflict-free solution without changing the landing ratios on the two runways. Compared with S3, S4 not only optimized the values of "D" and "S", but also saved an average 3.4% on fuel consumption per flight, and reduced the average flying time by 1 minute in the TMA. The comparison of horizontal trajectories for the different tests is shown in Fig. 5.12.

Computation time

We tested the computation time of our system with five different traffic scenarios: (6:00-8:00), (8:00-10:00), (10:00-12:00), (6:00-12:00), and (6:00-20:00). The parameter settings in the Simulated Annealing Algorithm were the same as in the Tab. 3.2. 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.5.7. 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 311 seconds to find a conflict-free solution for heavy traffic, 189 seconds for medium traffic, and only 75 seconds for light traffic. When we increased the size of the window and the shifting interval, i.e. W = 7200s, and S = 3600s, the running time of our system did not change significantly. The computation time also increased with the scenario time period. It took approximately 664 s for optimizing a 6 hour scenario, and 1450 s for a 12 hour scenario.

Segregated parallel approach operations

In order to ensure continuous traffic demand at runways and maximize runway usage, a minimum level of queuing is required. 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. As discussed in the results of re-sequencing ability in the section of Mixed parallel approach operation, MPS > 3 can be realized with our designed ML-PM route network. Consequently, higher runway landing throughput is possible. In this section, 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 maximally 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, see Fig. 5.13.

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 (Eurocontrol, 2009a).

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 E i , v E i , r E 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. 5.14. We have:

∆ ini i = ETA L i -t e i , ( 5.1)
∆ Act i = t L i -t E i , (5.2)
ε in i = ∆ Act i -∆ ini i , (5.3) ε out i =        t E i -t e i if t E i > t e i , 0 otherwise.
(5.4)

ε i = ε in i + ε out i , ( 5.5
)

ε average = n i=1 ε i i ∈ F.
(5.6) Fig. 5.15 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 can not be reduced to zero, however, with the P -OFF setting the total number of unsolved conflicts can be reduced to zero. This means that more dynamic position shifting helps to resolve conflicts. Fig. 5.15f 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. 5.15c 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, see Fig. 5.16.

Second, for constraints, we added a wake turbulence separation constraint between the departing aircraft and arriving aircraft successively operating on the same runway. Here, only the runway 18R-36L is a mixed operation. In addition, we also have to consider the runway occupancy time of the arrival aircraft as well. Given the maximum time to evacuate a runway at BCIA is 50 seconds, we add the following safety constraints:

t L i -t dep j ≥ s min i,j , i ∈ F arr , j ∈ F dep , ( 5

Data preparation

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 three departures-only traffic samples, each lasting 60 minutes (with 60, 70, 72, 74, 76 flights respectively). We use the QFU-18 SID/STAR network in Fig. 5.3 to study the independent parallel operations with integrated arrivals and departures.

Automated de-conflict performance

For departure-only scenarios, the results are shown in Fig. 5.18 and Tab. 5.8. Our system can solve all the potential conflicts for the different traffic demands with 60, 70, 72, 74 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 colour in Fig. 5.18. Fig. 5.18a 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. 5.18b, 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.

For the mixed arrival and departure scenarios, the results are shown in Fig. 5.19. There are four different traffic demands with 100 flights, 110 lights, 120 flights, and 130 flights respectively, of which 50% are arrivals and 50% are departures. Our system can still successfully resolve all the conflicts, see Fig. 5.19a. The results from iteration 0 to iteration 2800 are shown, and there are three window shifts. In the first shift, the number of initial

Climbing and descending performance

For the departures only case, the climbing performances with 60, 70, 72, 74 and 76 flights are shown in Fig. 5.20. For all j ∈ F dep , and with different decisions for t dep j , t T j and r j , the vertical profiles of all departures almost keep a similar shape, see Fig. 5.20a,Fig. 5.20b,Fig. 5.20c,Fig. 5.20d. In each figure, the lines with violet colour 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. 5.20 shows that the holding time on the sequencing legs increases with the number of departure For mixed arrival and departure traffic, we have to pay more attention to the hot spot areas, see Fig. 5.3, because some of the arrival and departure trajectories cross inside these hotspot areas (here, we only consider QFU-18). Taking the traffic sample with 130 flights as an example, Fig. 5.21a and Fig. 5.21b show the integrated departures and arrivals trajectories. We then zoom in the traffic around the hot spot areas. Fig. 5.21c shows that arrivals are all above the departures in hot spot area 1. Fig. 5.21d shows that departures going to VM may have intersection with the arrivals coming from DOGAR in spot area 2. Fig. 5.21e 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. 5.21f shows the result of this modification. We found that there is sufficient vertical separation between departures and arrivals now, because the departures

Capacity and efficiency

We compared our system performance regarding capacity with the maximum hourly movements observed at BCIA on the 07/09/2015 (Fig. 2.7, (baseline)). 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 of 72 arrivals and 77 departures. The declared total movement capacity is 88 aircraft per hour. By applying our system, we can reach a maximum of 99 arrivals without changing the runway-in-use, and a maximum of 120 movements per hour with integrated arrivals and departures. This corresponds to a 30% increase in maximum landing movements per hour, and a 36% increase in maximum total movements per hour, see Fig. 5.22, however, for the departure only case, there is not significant improvement.

As discussed in Chapter 2, PMS is able to easily absorb additional volume of traffic demands in the TMA. Let us look at further 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. 5.23a and Fig. 5.23b.

In Fig. 5.23a, 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 We can see clearly that as the traffic demand increases from 100 flights to 120 flights per hour, the gradient of the trend line also increases. In Fig. 5.23b, 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 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.

Conclusion

In this chapter, we have designed the SID and STAR for BCIA, which are then used for the evaluation of our optimization model. These advanced ML-PM route networks include Northinbound and outbound (QFU-36), South-inbound and outbound (QFU-18). Our experiments have been tested referring to QFU-18, which has higher complexity.

For the mode of mixed parallel approach operation, we have tested several scenarios to merge the arrival aircraft to land on two parallel runways. With the numerical results, it is found that: 1) our proposed system could achieve a near CDA profile (except a reasonable level-off on the sequencing leg) for aircraft to approach the runway, an easy re-assignment of landing runway, an automated conflict-free merging, and a more flexible sequencing. 2) This optimization system is dynamical and has stable de-conflict performances to handle the routine traffic in busy TMA. The proposed optimization solutions concern the control of entry time, entry speed, turning time on sequencing leg and runway assignment, which make the theoretical near-optimal solution easily implemented in the real world.

3) The user-defined parameter settings in our algorithm could be flexibly changed according to the different requirements. They can provide real time runway allocation, the predicted average flight time, and fuel consumption information to assist airport and airlines to make good decisions on reducing the operating cost.

For the mode of segregated parallel approach operation, we have tested the approach efficiency of our ML-PM route structure in terms of sequencing and merging efficiency and vertical flight efficiency. Two different sequencing techniques have been designed: a constrained sequencing technique which minimizes the value of total position shifts P -ON, and tries to maintain the FCFS sequencing, and another relaxed sequencing technique with P -OFF, which is our proposed system. The numerical results show that our proposed sequencing technique shows advantages on handling the much denser traffic situation. With more dynamic position shifting, it can find conflict-free solution for the super-dense case with 99 arrivals. The average delay and average flying time in TMA, average additional flying time are all shorter than the conventional sequencing technique with P -ON. In addition, the percentage of time flown level during approach is less than 12% of the total flown time in TMA with our ML-PM route structure, either with P -ON or P -OFF sequencing technique. The percentage of average level-off time maintains relatively stable between 7% and 11%.

For the mode of independent parallel operation with integrated arrivals and departures, we adjust the optimization algorithm for departure, including the new decision variables for departure aircraft and the required minimum separation between the landing aircraft and take-off aircraft on the same runway. The numerical results show that our system can provide a stable de-conflict performance. The mixed arrival and departure aircraft show a good continuous descending and climbing performance. Meanwhile, we also analyze the hot spots between departures and arrivals. In addition, we study the relationship between the efficiency and capacity. Compared with radar-based baseline, our proposed system could significantly increase the capacity at BCIA, from maximum 72 landing movements per hour to 99 landing movements per hour, from maximum 88 total movements per hour to 120 total movements per hour. We also study the relationship between the performance of our system and the different available length of sequencing legs.

Chapter 6 Conclusion and future work 6.1 Conclusion

In this dissertation, novel integrated sequencing and merging techniques have been developed to devise more efficient aircraft trajectory planning tools for parallel runway operations at busy airports.

Since aircraft arrival management in the TMA is very complex, sequencing and merging arrivals involves not only a complex decision problem but also an efficient control problem. Our first contribution was to design a good route network system, ML-PM, capable of supporting our approach. The proposed network system is based on advanced avionics techniques. With this in place, we then contributed a complete modelling framework capable of dynamically handling large traffic demand in routine dense operations. It realizes automated conflict resolution and efficient trajectory planning for multiple parallel runway. The proposed modelling approach is based on a hybrid RHC and SA algorithms, together with a suitable mathematical optimization formulation. Finally, we used the proposed modelling framework to solve a large number of trajectory planning problems. We suggest several conclusions based on results obtained:

1. The ML-PM route network is very good at supporting an economical 3D descent to the runway, an easy re-assignment of landing runway for the mixed parallel approach, and more flexible sequencing for the segregated parallel approach in the TMA. 2. Our system has a stable conflict-resolving performance, and is able to handle dense routine traffic in busy TMAs. The proposed optimization solutions taking the control of entry time, entry speed, turning time on sequencing leg and runway assignment as control variables, makes the theoretical near-optimal solution easily to implement in the real world.

3. The user-defined parameter settings in our algorithm can be changed according to different user preferences. Our system can provide real time runway allocations, the predicted average flight time and fuel consumption information etc. to assist airport and airlines to make good decisions. 4. Our system can find good solutions to flight planning problems considering issues: (i) with segregated approach to parallel runways, (ii) mixed approach to parallel runways, and (iii) integrated arrivals and departures with multi-parallel runways. It can also increase the efficiency of trajectory planning by dynamically controlling the length of sequencing legs in the ML-PM according to the traffic demand.

Overall, we conclude that the ML-PM-based integrated sequencing and merging techniques studied in this dissertation have a strong potentiality to handle dense routine arrivals for multiple parallel runway operations, resulting an increase in airport capacity, and supporting robust and efficient arrival management.

Future work

We see three main directions for future research that could improve the proposed optimization approach:

1. The optimization of arrival management, creating an efficient sequencing and conflictfree merging, is based on accurate trajectory generation. In our model, we did not consider the effect of wind. In future work, wind forecast modelling could be added to improve the predictability of the trajectory.

2. There are a range of heuristic algorithms that can solve the multi-objective optimization problems. For example, the GA algorithm, and the Particle Swarm Optimization algorithm. It will be interesting to test these algorithms on our problem, and compare their performance to the RHC-SA algorithm employed. Furthermore, intelligent learning algorithm could be applied to dynamically control the length of sequencing legs according to the traffic demand in the TMA. This could lead to more precisely and economically trajectory control, resulting in fuel savings and reduced emissions.

3. Route. Stage in which most of the flight takes place, flying at the established level, with practically constant configuration and speed. This phase takes place in the CTA, which includes airways (AWY).

4. Arrival. After the Route stage, the aircraft begins to descend, following the STAR, going to specific points where, if necessary, it will wait until it is authorized to continue with the approach stage. Like the departure stage, this takes place within the TMA.

5. Approach. This begins at the way point where the arrival route ends. These operations take place in the TMA, until they are transferred to the CTR just before landing begins.

The approach is determined by the approach charts, which can be instrumental or visual, depending on the type of approach the aircraft is carrying out. If once the decision speed is reached (DA/H) visual contact has not been established with the runway, the approach must be abandoned.

6. Landing. Flight stage which begins after the approach and ends when the aircraft has come to a complete stop. Like take-off and approach, it takes place within the CTR. In the case of failures during the landing stage, if the circumstances and the configuration of the aircraft allow it, the landing must be aborted.

A.2 ATFM phases

The objective of ATFM is to optimize traffic flows according to air traffic control capacity while enabling airlines to operate safe and efficient flights. In Europe, the ATFM activities are divided into three phases:

1. Strategic phase. About one year before the flight takes place until one week before real time operations. During this phase, the Network Manager Operations Center (NMOC) helps the ANSPs to predict what capacity they will need to provide in each of their air traffic control centres. In addition, a routing scheme is prepared -a structure of air routes across Europe designed to balance the air traffic flows and maximize capacity. This also includes avoiding imbalances between capacity and demand for events taking place a week or more in the future (large-scale military exercises, major sports events, etc.).

4

 4 -clés: Gestion du trafic aérien, Trajectoire basée sur l'opération, Fusion et séquençage, Climbing and descending performance . . . . . . . . . . . . . . . . . 5.4.5 Capacity and efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Conclusion and future work 6.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. 1

 1 CAS speed schedule (descent) . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Geographical flight distribution in Data set 1 (Mixed parallel approach) . . . 5.2 Solution quality of SA compared with HC . . . . . . . . . . . . . . . . . . . 5.3 Average fuel consumption and flight time in different routes (Mixed parallel approach) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 The position shifts on the arrival sequencing (Mixed parallel approach) . . . 5.5 Comparison of makespan, fuel consumption, and flying time between "MPS > 3" and "MPS <= 3" (Mixed parallel approach) . . . . . . . . . . . . . . . . 5.6 Test results with or without runway landing balancing (Mixed parallel approach) 5.7 Comparison of computation times (Mixed parallel approach) . . . . . . . . . 5.8 Performance with departures only case . . . . . . . . . . . . . . . . . . . . . A.1 PBN Navigation specifications for different flight phases . . . . . . . . . . . . A.2 ISA (mean sea level conditions) . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Way-points, coordinates, target altitude . . . . . . . . . . . . . . . . . . . . . A.4 Designators of links . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Designators of routes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 Set of routes for North-inbound or South-inbound . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Total runway movement and on-time performance in China between 2006-2015

1 .

 1 ATS provided by ANSP is the most important part of ATM. It aims to provides safe, secure, and efficient management for the assigned airspaces, as well as facilitate the smooth operation and punctuality of flight schedules. ATS includes Flight Information Service (FIS), Alerting Service (AS), Air Traffic Advisory Service (ATAS), and Air Traffic Control (ATC) service. The ATC service is the most important element of ATS. It is to prevent collisions between aircraft, between aircraft and obstructions on the manoeuvring areas at airports, and expedite and maintain an orderly flow of air traffic. ATS is further divided into three sub-services: Area Control (ACC) service, Approach Control (APP) service, and Aerodrome Control (TWR) service, according to different flight phases (more information about flight phases refer to A.1). ATSs are provided by ANSPs to air traffic during all phases of flights. Each ANSP has branch control units that provide ATS according to the flight phase. Generally, there are Air Traffic Control Towers (ATCTs), Terminal Radar Approach Control Centers (TRACONs), and Air Route Traffic Control Centers (ARTCCs) 1 . In the controlled airspace aircraft will be under the control of a single control unit at all times, and different control units cooperate together to ensure the safety of flight.

Figure 1 . 2 :

 12 Figure 1.2: Key elements of ATM system

Figure 1 .

 1 Figure 1.3: Airspace system

Figure 2 . 1 :

 21 Figure 2.1: A centralized approach: sequencing problem with different solutions

Figure 2 . 3 :

 23 Figure 2.3: Examples of merge topology

Figure 2 . 4 :

 24 Figure 2.4: Concepts of CDA and CCO

Figure 2 . 5 :

 25 Figure 2.5: Conventional routes compared to PBN-based routes

Figure 2

 2 Figure 2.6: Basic PMS topology for single runway

Figure 2

 2 Figure 2.7: Hourly movements at BCIA on 07/09/2015

Figure

  Figure 2.8: TMA of BCIA

Figure 2 . 9 :

 29 Figure 2.9: Layout of runways at BCIA

Figure 3 . 1 :

 31 Figure 3.1: Optimization framework

Figure 3

 3 Figure 3.3: Conventional sequencing and merging in TMA

Figure 3 . 5 :

 35 Figure 3.5: Segregated vertical levels on sequencing leg for specific category of aircraft

Figure 3

 3 Figure 3.7: RHC approach to decomposes the original problem into several sub-problems

F

  OPT ← Active and on-going flights; 7: while T e (k) < t FINAL do 8: if at least one active flight in F OPT then 9:Sub-problem: optimize considering F OPT ;

Figure 4 . 1 :

 41 Figure 4.1: Relationship between simulation module and optimization module

Figure 4 .

 4 Figure 4.3: Flight procedure in the merging zones for mixed parallel approach

  {Entry point, IAF, Turning point, M1, FAF1, runway18R-36L}} (4.8) l 2 = {x ∈ {Entry point, IAF, Turning point, M2, FAF2, runway01-19}} (4.9)

4 :

 4 load the route network for aircraft i; 5:

  , ROD, Drag, Thrust, fuel flow rate with different values of AC and flight altitude; 11: calculate the motion of aircraft based on the value of AM and OP ; 12: update current altitude, travel time and travel distance of flight i 13:end for 14: end procedure

Figure 4 . 5 :

 45 Figure 4.5: Motion of aircraft

m

  Va = T hr -D -mg sin γ a ,

Figure 4 . 6 :

 46 Figure 4.6: Example of conflict calculation

Figure 4 . 7 :

 47 Figure 4.7: Example of link conflict detection

  4.10.The node "Turning point" is a time-based point. For aircraft flying on the outer sequencing leg, it is denoted by b. For aircraft flying on the inner sequencing leg, it is denoted by c.

Figure 4 .

 4 Figure 4.9: Example of conflict detection in merging zone

Figure 4 .

 4 Figure 4.10: Sequencing legs equidistant from single merge point

Figure 4 .

 4 Figure 4.11: Merge conflict detection in mixed parallel approach

Figure 5 . 1 :

 51 Figure 5.1: ML-PM based STAR design for BCIA

Figure 5 .Figure 5 . 4 :

 554 Figure 5.2: ML-PM based SID design for BCIA

Figure 5 . 5 :

 55 Figure 5.5: Hourly landing demand in Data 1 (Mixed parallel approach)

Figure 5 . 9 :

 59 Figure 5.9: An example of "Medium" aircraft descent performance: A3245

Figure 5 .

 5 Figure 5.10: Total trajectories in BCIA TMA (Mixed parallel approach)

Figure 5 .

 5 Figure 5.13: Geographical distribution of flights (Segregated parallel approach)

Figure 5 .

 5 Figure 5.14: Definition of additional time

( a )

 a Total number of unsolved conflicts (b) Average square delay (c) Average flying (transit) time (d) Average fuel consumption (e) Value of ε average (f) Maximum number of position shifts

Figure 5 .Figure 5 .

 55 Figure 5.15: Results with minimizing position shift P ON/OFF

Figure 5 .

 5 Figure 5.18: Automated de-conflict performance with only departures

Figure 5 .

 5 Figure 5.19: Automated de-conflict performance with integrated arrivals and departures

Figure 5 .

 5 Figure 5.20: Climbing performance in departures only cases

( a )Figure 5 .

 a5 Figure 5.21: Hot spots areas with integrated arrivals and departures

Figure 5 .

 5 Figure 5.23: Performance with different traffic demand and available sequencing length

  . . . . . 6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Flight phases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Total runway movement and on-time performance in China between 2006-2015 1.2 Key elements of ATM system . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Airspace system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Radar vector techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Conventional sequencing and merging in TMA . . . . . . . . . . . . . . . . . 3.4 Sequencing and merging topologies for parallel runways . . . . . . . . . . . . 3.5 Segregated vertical levels on sequencing leg for specific category of aircraft . 3.6 Minimum separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 RHC approach to decomposes the original problem into several sub-problems 4.1 Relationship between simulation module and optimization module . . . . . . .1 ANSP, Air Traffic Services (ATS) and controlled airspace in different flight phases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 International Civil Aviation Organization (ICAO) minimum separation . . . 3.2 Experience-based parameters setting . . . . . . . . . . . . . . . . . . . . . . 3.3 Performance indicator in an active window at time t . . . . . . . . . . . . . .

	List of Figures
	Bibliography
	1.1
	A Appendix
	A.1 xiii xiv

A.2 Air Traffic Flow Management (ATFM) phases . . . . . . . . . . . . . . . . . A.3 PBN navigation specifications . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 International Standard Atmosphere (ISA) . . . . . . . . . . . . . . . . . . . A.5 Route network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 A centralized approach: sequencing problem with different solutions . . . . . 2.2 A decompression approach: RHC-based sequencing problem . . . . . . . . . 2.3 Examples of merge topology . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Concepts of CDA and CCO . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Conventional routes compared to PBN-based routes . . . . . . . . . . . . . . 2.6 Basic PMS topology for single runway . . . . . . . . . . . . . . . . . . . . . 2.7 Hourly movements at BCIA on 07/09/2015 . . . . . . . . . . . . . . . . . . . 2.8 TMA of BCIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Layout of runways at BCIA . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Independent and dependent parallel approaches . . . . . . . . . . . . . . . . 2.11 Parallel runway operation pattern in our study case . . . . . . . . . . . . . . 2.12 Four kinds of traffic pattern with real traffic data . . . . . . . . . . . . . . . 2.13 South-inbound daily real traffic and clustered traffic . . . . . . . . . . . . . . 3.1 Optimization framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 4.2 An example of ML-PM route network . . . . . . . . . . . . . . . . . . . . . . 4.3 Flight procedure in the merging zones for mixed parallel approach . . . . . . 4.4 Multiphase dynamical control process . . . . . . . . . . . . . . . . . . . . . . xv List of Tables 1

Table 1 .

 1

1: ANSP, ATS and controlled airspace in different flight phases Phases ANSP ATS Controlled airspace Take-off Air traffic control tower TWR Control Zone (CTR) Departure Terminal radar approach control centre APP CTA or TMA Route Air route traffic control centre ACC Control Area (CTA) Arrival Air route traffic control centre ACC CTA or TMA Approach Terminal radar approach control centre APP CTA or TMA Landing Air traffic control tower TWR CTR 2. ASM is the process by which airspace options are selected and applied to needs of the airspace users (ICAO

  E i -the actual entry speed at TMA. 3. t T i -the actual turning time on the sequencing legs. 4. r i -the actual landing runway.

1. t E

i -the actual entry time at TMA.

2. v

Table 3

 3

	.2: Experience-based parameters setting
	Sliding Window
	Duration of window			3600 seconds
	Window shifting interval			1800 seconds
	Simulated Annealing Algorithm
	Initial Temperature for heating			0.01
	Heating rate			1.1
	Number of transition for heating or cooling 200
	Cooling rate			0.99
	Cooling stopping criterion			T < 0.0001 × T init
	Table 3.3: Performance indicator in an active window at time t
	No. Flights	1	2	3	... ... ... N
	Status of flights (On-going or Active )

  + w 2 + w 3 ) < I <= (w 1 + w 2 + w 3 + w 4 ) then

	5:	else if w 1 < I <= (w 1 + w 2 ) then
	6:	Change the value of g	Update the entry speed
	7:	else if (w 1 + w 2 ) < I <= (w 1 + w 2 + w 3 ) then
	8:	Change the value of h	Update the turning time
	9: else if (w 1 10: Change the value of r	Update the landing runway
	11:	Generate randomly a value S runway	S runway ∈ [0, 1]
	12:	if S runway <= β then
	13:	r ← 1	Choose runway 18R/36L
	14:	else
	15:	r ← 0	Choose runway 01-19
	16:	end if
	17:

Table 4 .

 4 1: CAS speed schedule (descent)

	Altitude	CAS (knot)
	from 0 to 999 ft	C V min × (V stall ) LD + 5
	from 1,000 to 1,499 ft	C V min × (V stall ) LD + 10
	from 1,500 to 1,999 ft	C V min × (V stall ) LD + 20
	from 2,000 to 2,999 ft	C V min × (V stall ) LD + 50
	from 3,000 to 5,999 ft	min(V des,1 , 220)
	from 6,000 to 9,999 ft	min(V des,1 , 250)
	from 10,000 ft to Mach transition altitude V des,2
	above Mach transition altitude	M des
	tropopause:

  4.7. Consequently, if the following constraints can not be satisfied, then L i increases by one:

	∆t u i,j := t wa j -t wa i ≥ s min i,j ,	(4.33)
	∆t u i,j := t w b j -t w b i ≥ s min i,j ,	(4.34)
	∆t u i,j × ∆t u i,j ≥ 0.	(4.35)

Table 5 .

 5

1: Geographical flight distribution in Data set 1 (Mixed parallel approach)

Table 5 .

 5 2: Solution quality of SA compared with HC

	Instance	No. Active A/C	Remaining conflicts	Total running time[s]
			HC SA	HC SA
	Small size window (W=3600s, S=1800s)
	1th window	22	0	0	21	371
	2th window	38	2	0
	3th window	53	8	0
	4th window	38	10	0
	5th window	17	2	0
	6th window	10	0	0
	Big size window (W=7200s, S=3600s)
	1th window	40	6	0	58	385
	2th window	82	28	0
	3th window	23	8	0
	Super big size window (W=10800s, S=7200s)
	1th window	128	32	0	23	224

Table 5 .

 5 3: Average fuel consumption and flight time in different routes (Mixed parallel approach)

	Entry way-point	Average Fuel Consumption [kg]	Average Flying Time [sec(min)]	Baseline [sec(min)]
		Medium	Heavy	Medium	Heavy
	KM	136	948	833(14)	804 (13)	880(15)
	JB	193	1140	1286 (21)	1171 (20)	1698(28)
	BOBAK	186	1118	1254 (21)	1135 (19)	1663(28)
	VYK	259	1161	1448 (24)	1188 (20)	1663(28)
	DOGAR	258	1167	1455(24)	1246 (21)	1236(21)
	GITUM	182	855	799 (13)	756 (13)	1116(19)

Table 5 .

 5 4: The position shifts on the arrival sequencing (Mixed parallel approach)

	Test No.	Landing runway	CPS>3 (Perc. of A/C)	CPS<=3 (Perc. of A/C)	MPS
	Test 1	01-19	21.31%	78.69%
		18-36	20.90%	79.10%
	Test 2	01-19	14.75%	85.25%
		18-36	17.91%	82.09%
	Test 3	01-19	18.03%	81.97%
		18-36	19.40%	80.60%
	Test 4	01-19	18.03%	81.97%
		18-36	20.90%	79.10%
	Test 5	01-19	24.59%	75.41%
		18-36	19.40%	80.60%
	Test 6	01-19	24.59%	75.41%
		18-36	17.91%	82.09%
	Test 7	01-19	27.87%	72.13%
		18-36	25.37%	74.63%
	Test 8	01-19	31.15%	68.85%
		18-36	17.91%	82.09%
	Test 9	01-19	29.51%	70.49%
		18-36	19.40%	80.60%
	Test 10	01-19	27.87%	72.13%
		18-36	16.42%	83.58%
	Average	01-19	23.77%	76.23%
		18-36	19.55%	80.45%

Table 5 .

 5 6: Test results with or without runway landing balancing (Mixed parallel approach)

	Tests	Arrival Flow	Landing	Initial ratio	Optimized ratio	Criteria of objectives	Average transit
			balancing	18-36 01-19	18-36	01-19	C	D	S	Fuel[kg] Time[min]
	S1	Asymmetric	OFF	69%	31%	69%	31%	2	176.00	1.42	497	21
	S2	Asymmetric	ON	69%	31%	62%	38%	0	36.80	1.39	492	20
	S3	Normal	OFF	52%	48%	52%	48%	0	30.80	1.44	500	21
	S4	Normal	ON	52%	48%	54%	46%	0	13.68	1.41	483	20

Table 5 .

 5 7: Comparison of computation times (Mixed parallel approach)

	Running times

Table 5 .

 5 8: Performance with departures only case

	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
				88

  Table A.6: Set of routes for North-inbound or South-inbound

				continued from previous page continued from previous page continued from previous page continued from previous page continued from previous page continued from previous page
	No.waypoint X(NM) No.waypoint X(NM) No.waypoint X(NM) No.link beginning No. waypoint ending No. waypoint flagOfArcs flagPhase Y(NM) Z(ft) Z(m) Name of waypoint Y(NM) Z(ft) Z(m) Name of waypoint Y(NM) Z(ft) Z(m) Name of waypoint No.link beginning No. waypoint ending No. waypoint flagOfArcs flagPhase No. Route for ARR or DEP A set of links
	12 47 82	65.5879 69.9514 88.4994 flagOfArcs: 0-normal link,1-sequencing leg, 2-merging link 8.0810 18696 5700 BOBAK 71.7727 5904 1800 centre of circle 2 18.7266 19680 6000 VM flagOfArcs: 0-normal link,1-sequencing leg, 2-merging link R11 10 8 9 91 14 16 17 18
	13 48 83	22.7754 110.2918 91.7575 flagPhase: 1-departure, 2-arrival 83.3019 55.4109 64.9627 flagPhase: 1-departure, 2-arrival R12	14760 4500 17712 5400 11808 3600 11 12 14 16 17 18 KM TONIL AMVIK
	14 49 84	90.5271 85.0334 101.5193 9 93 DER 10 76		105.6154 95.6326 95.8680		11808 3600 12792 3900 12792 3900 8 83	GITUM w9 CDY	0 0	1
	15 50 85	68.5145 92.9876 74.3708 R13 11 77	14 92		75.3850 8.8073 104.6043		3936 19680 6000 1200 12792 3900 5 95 43 35 36 25 M1 LADIX YV	0 0	1
	16 51 86	70.8306 22.7758 74.2610 R14 12 78	5 19		75.8830 83.3016 50.5136		2952 14760 4500 900 4920 1500 6 86 43 35 36 26 M2 KM D1	0 0	1
	17 52 87	69.9493 56.0169 76.7234 R15 13 79	4 78		71.7710 2.9647 50.7242		2952 17712 5400 900 5904 1800 21 87 43 35 36 27 FAF1 JB D2	1 0	1
	18 53 88	72.2372 65.6008 93.9106 R16 14 80	6 86		70.8648 8.0980 53.3919		1968 18696 5700 600 11808 3600 27 97 43 35 37 38 28 29 30 FAF2 BOBAK beginning of outer sequencing leg for departure 1 2 1
	19 54 89	72.2572 73.3121 63.9266 R17 15 81	21 97		66.6838 12.3209 46.7727 ARR	1	114.8 17712 5400 35 9840 3000 15 92 43 35 37 38 28 31 Threshold of 18R VYK beginning of inner sequencing leg for departure 2 1 1 1 21 23 16 17 18;1 13 15 19 20
	20 55 90	74.2823 83.2220 80.6335 R18 16 82	27 97		66.2397 10.6040 33.7764 ARR	2	114.8 13776 4200 35 9840 3000 16 91 43 35 37 38 32 33 Threshold of 19 DOGAR midpoint of outer sequencing leg 2 1 1 2 3 4 5 21 23 16 17 18;2 3 4 5 13 15 19 20
						ARR	3	6 3 4 5 21 23 16 17 18;6 3 4 5 13 15 19 20
	21 56 91	61.6694 90.4758 79.7342 R19 17 83	16 91		90.1337 105.6310 37.6740 ARR	4	7872 11808 3600 2400 11808 3600 18 93 43 35 37 38 32 34 virtual point for merging to M1 GITUM midpoint of inner sequencing leg 0 1 1 7 8 9 91 14 16 17 18;7 8 9 91 22 24 15 19 20
	22 57 92	60.0168 50.8892 60.1047 R20 18 84	18 87		93.7968 43.6883 45.5925 ARR	5	7872 18696 5700 2400 9840 3000 20 98 44 39 41 42 25 midpoint of outer sequencing leg waypoint1 ending point of outer sequencing leg 0 2 1 10 8 9 91 14 16 17 18;10 8 9 91 22 24 15 19 20
	23 58 93	61.6694 81.9506 89.9579 R21 19 85	15 98		90.1337 98.9842 52.7784 ARR	6	9840 11808 3600 3000 11808 3600 17 93 44 39 41 42 26 midpoint of inner sequencing leg waypoint2 ending point of inner sequencing leg 0 1 1 11 12 14 16 17 18;11 12 22 24 15 19 20
	24 59 94	82.2957 77.4844 76.7234 R22 20 86 South-inbound 17 98	92.5434 80.5667 50.7242 DEP	7	7872 11808 3600 2400 5904 1800 19 90 44 39 41 42 27 end of outer sequencing leg waypoint3 centre of circle 2 0 1 1 79 84 86 87 77 68;78 80 81 77 68
	25 60 95	50.1618 83.1354 50.8892 R23 21 87	4 90		75.5202 53.1289 43.6883 DEP	8	9840 11808 3600 3000 13776 4200 22 92 44 39 40 28 29 30 end of inner sequencing leg waypoint4 wp1-2 1 1 1 79 84 86 87 77 69;78 80 81 77 69
						DEP	9	79 84 86 87 77 70;78 80 81 77 70
	26 61 96	69.9514 53.4531 110.2918 6 90 R24 22 88		71.7727 23.3265 55.4109 DEP	10	7872 18696 5700 2400 17712 5400 23 82 44 39 40 28 31 centre of circle waypoint5 TONIL 1 0 79 84 86 88 72;78 80 82 71 72	1
	27 62 97	61.6694 78.7523 76.2610 R25 23 89	22 54		90.1337 37.0458 50.5136 DEP	11	7872 7872 9840 27 71 44 39 40 32 33 2400 virtual point for merging to M2 2400 M1 3000 virtual point from D1 1 0 2 79 84 85 73;78 80 82 83 73
	28 63 98	73.7243 80.9297 76.7234 R26 24 90	23 55		63.9325 36.4526 50.7242 DEP	12	114.8 8856 11808 3600 35 2700 21 71 44 39 40 32 34 Threshold of 36R M2 virtual point from D2 1 0 79 84 85 76 74;78 80 82 83 76 74	2
	29 64 99	72.5171 74.4680 92.9876 ARR 25 91	45 100		64.9957 55.8343 8.8073 DEP	13	114.8 3936 19680 6000 35 1200 30 6 79 84 85 76 75;78 80 82 83 76 75 Threshold of 36L FAF1 LADIX 0 0	2
	30 65 100	22.7772 76.6385 86.56471756 66.88731909 13776 4200 83.2998 13776 4200 55.0524 4920 1500 45 31 9 7 R27 26 92 45 46 56 57 64 65 KM2 FAF2 wp100 0 0 ARR 14 45 46 58 59 60 66 67;45 46 56 57 64 65	2
						ARR	15	47 49 58 59 60 66 67;47 49 56 57 64 65
	31 66	22.5597 72.5171 R28 27	45		66.7890 64.9957 ARR	16	13776 4200 114.8 35 32 47 49 56 57 64 65 SOSDI Threshold of 36L 0 48 49 58 59 60 66 67;48 49 56 57 64 65
	32 67	20.8248 74.5291 R29 28	46		45.3136 64.1645 ARR	17	21320 6500 114.8 35 34 48 49 56 57 64 65 RENOB Threshold of 01 0 50 59 60 66 67;89 63 57 64 65
	33 68	88.4994 64.3110 R30 29	34		18.7266 31.3209 ARR	18	19680 6000 7872 2400 33 89 63 57 64 65 VM Beginning point of outer sequencing leg 0 51 59 60 66 67;90 63 57 64 65
	34 69	91.7575 94.1642 R31 30	33		64.9627 38.5068 Table A.4: Designators of links 11808 3600 AMVIK 9840 3000 Beginning point of inner sequencing leg 50 0 Table A.5: Designators of routes 90 63 57 64 65 ARR 19 52 53 54 55 61 60 66 67;52 53 54 55 62 63 57 64 65
	35 70	101.5193 84.8398 34 R32 31 North-inbound	95.8680 19.5048 DEP	20	12792 3900 7872 2400 48 52 53 54 55 62 63 57 64 65 CDY midpoint of outer sequencing leg 0 44 39 41 42 25;43 35 36 25
	36 71	74.3708 83.9405 No.link beginning No. waypoint ending No. waypoint flagOfArcs flagPhase 104.6043 12792 3900 YV 23.4024 9840 3000 midpoint of inner sequencing leg 32 46 49 0 DEP 21 44 39 41 42 26;43 35 36 26 No. Route for ARR or DEP A set of links R33 45 46 58 59 60 66 67 DEP 22 44 39 41 42 27;43 35 37 38 28 29 30
	37 72	68.5145 98.1169 flagOfArcs: 0-normal link, 1-sequencing leg, 2-merging link 75.3850 5904 1800 D1 39.1203 7872 2400 end of outer sequencing leg 33 49 35 0 ARR R34 47 49 58 59 60 66 67 DEP 23 44 39 40 28 29 30;43 35 37 38 28 29 30
	38 73	70.8306 68.1329 flagPhase: 1-departure, 2-arrival 75.8830 32.5011 34 49 R1 R35 DEP 24	7872 9840 36 1 13 15 19 20 2400 D2 3000 end of inner sequencing leg 0 48 49 58 59 60 66 67 44 39 40 28 31;43 35 37 38 28 31
	39 74	82.2957 80.9297 R2 35 1 R36	13 37		92.5434 36.4526 DEP	25	9840 7872 4 41 2 3 50 59 60 66 67 3000 end of inner sequencing leg for departure 2400 centre of circle 0 2 2 4 5 13 15 19 20 44 39 40 32 33;43 35 37 38 32 33
	40 75	50.1618 83.9405 R3 36 2 R37	11 41		75.5202 23.4024 DEP	26	11808 3600 7872 2400 1 45 6 3 4 51 59 60 66 67 end of outer sequencing leg for departure virtual point for merging to M1 0 2 1 5 13 15 19 20 44 39 40 32 34;43 35 37 38 32 34
	41 76	82.2957 84.4319 R4 37 3 R38	1 41		92.5434 23.5256		9840 7872 2 44 7 8 52 53 54 55 61 60 66 67 3000 virtual point from D1 2400 virtual point for merging to M2 0 2 1 9 91 22 24 15 19 20
	42 77	82.2957 72.2572 R5 38 4 DER	2 44		92.5434 66.6838		11808 3600 114.8 35 3 46 10 8 9	virtual point from D2 Threshold of 18R 0 1 91 22 24 15 19 20 2
	43 78	60.0168 73.7696 R6 39 5 R39	3 38		93.7968 63.9813		9840 114.8 4 42 11 12 22 24 15 19 20 3000 midpoint of outer sequencing leg 35 Threshold of 18L 0 2 2 78 80 81 77 68
	44 79	61.6694 22.7772 R7 40 6 R40	12 42		90.1337 83.2998		11808 3600 13776 4200 1 46 1 21 23 16 17 18 midpoint of inner sequencing leg KM2 0 2 1 78 80 81 77 69
	45 80	46.1998 22.5597 R8 41 7 R41	10 42		76.1903 66.7890		9840 13776 4200 3000 8 43 2 3 4 78 80 81 77 70 beginning point of outer sequencing leg SOSDI 0 2 1 5 21 23 16 17 18
	46 81	80.3438 20.8248 R9 42 8 R42	8 43		89.0282 45.3136		11808 3600 21320 6500 7 45 6 3 4 78 80 82 71 72 beginning point of inner sequencing leg RENOB 0 2 1 5 21 23 16 17 18
		9	R10	7	continued on next page continued on next page 100 continued on next page 7 8 9 91 14 16 17 18 0 continued on next page	2
					continued on next page continued on next page

These items are US terms, not international standard

RNAV allows an aircraft to choose any course within a network of navigation beacons, rather than navigate directly to and from beacons. This can reduce flight distance, congestion, and allow aircraft to fly towards airports without beacons.

ADS-B is a surveillance technology in which an aircraft determines its position via satellite navigation and periodically broadcasts it, enabling it to be tracked. The information can be received by air traffic control ground stations as a replacement for secondary radar. It can also be received by other aircraft to provide situational awareness and allow self separation.

The runway systems at major airports are highly constrained resources. Runway assignments have to be considered in multi-runway operations, because the effective use of airport capacity depends on the optimal use of all available runways[START_REF][END_REF]. Runway assignment is typically dependent on the airport configuration, wind direction and speed, the direction of arriving aircraft, departure routes and the gate assigned to the aircraft(Brinton, 1992). For an aircraft, its runway assignment is generally specified by the flight plan according to predicted wind situation at airport. While an aircraft is approaching the runways, adjustments can be made to the flight plan by assigning the aircraft to an alternative runway, which is known as runway allocation. For controllers, the arrival flows into and out of the TMA are asymmetric due to different airlines schedules. Generally, aircraft are scheduled to land on the closest runway to their entry side of TMA, which can reduce the crossing of trajectories, thus reducing the complexity of traffic. But it may create heavy demands on one runway, while another one is underutilized. It also causes unnecessary delays. Airlines prefer their aircraft to land on or depart from runways close to their ramps or base terminal, which is convenient for passengers and also economical for airlines. For airports, if a lot of aircraft land on one runway which is far from their assigned gate, congestion on the connecting taxiway often occurs in peak traffic periods. With runway re-assignment (runway allocation), we could balance the runway landings and departures at all available runways, so as to reduce

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

BADA (Base of Aircraft Data) is the world's leading aircraft performance model, managed by EURO-CONTROL for use by the aviation community. The main application of BADA is trajectory simulation and prediction.

ICAO mandates separation minima based upon wake turbulence categories. These minima are typically categorized 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.

The methodology defined here was only applied to Beijing airport. Given the complexity of traffic in busy TMAs, studies for simultaneously scheduling arrivals for a group of busy airports in the same TMA using a ML-PM-based route structure, or simultaneously scheduling arrivals including helicopter missions, may be interesting.

Pre-tactical phase. Six days before real time operations. The task of the NMOC staff is to: coordinate the definition of a daily plan aimed at optimizing the overall ATM network performance and minimizing delay and cost, after a collaborative decision making process involving operational partners -such as ATC units and aircraft operators;

Acknowledgements

Definition 3 A transition function is a mapping of a set of states and a set of symbols onto the original set of states. The interpretation of δ(q, x) = p is that for a finite state machine M in state q with input x will be in state p after processing the input string x.

Definition 4 OP is a set of operational procedure variables. It represents different operational requirements associated with state q n , OP ∈ {Angle, CAS, F L}, which correspond to a 3 degree path angle descent, constant CAS, maintain a constant Flight level. Consider the mixed ML-PM system as an example, see Fig. 4.4. Here, Σ = {w i , i ∈ N} is a list of way-points, Q = {q n , n = 0, 1, ..., 10} is the state of aircraft in some specific phase, and Conflict detection modelling enables the computation of conflicts among aircraft. Let us introduce several definitions before presenting the approach.

Definition 7

The total number of conflicts C i of each aircraft i ∈ F is the union of the following three subsets:

1. L i contains all link conflicts, referring to catch-up conflicts and overtake conflicts between two successive aircraft flying on the same route.

Chapter 5

Experiments and numerical results

In this chapter, we first introduce our SID and STAR route network design for BCIA. Then, based on the default configuration QFU 18, we present the experiments and numerical results for three different operational patterns: the mixed parallel approach, the segregated parallel approach and the integrated arrival management with departure.

Route network design

In respect of real operational environment, referring to the published RNAV STAR chart and the instrument approach chart at BCIA, the horizontal STAR route network of ML-PM system is designed, see Fig. 5.1. The left part is for North-inbound operations, and the right part is for South-inbound operations. Each part can model either segregated ML-PM topology without considering runway allocation, or the mixed ML-PM topology considering runway allocation. In QFU-18 arrivals, M1/M2 are much closer to FAFs, because the airspace in the North of Beijing TMA is very limited, see Fig. 2.8. We use the South-inbound layout as the default layout. Flights entering from West (coming from entry points KM, JB and BOBAK) follow a fixed set of routes, and are initially planned to merge at M1. Flights entering from East (coming from entry points GITUM, DOGAR, and VYK) follow another set of routes, and are initially planned to merge at M2. For all the available arrival routes see Appendix A.5. PMS is generally designed for managing arrival flows, here we plan to apply the ML-PM route structure for the departure flows as well, which will be helpful for implementation of CCO at the busy airport. Therefore, referring to the published RNAV SID chart at BCIA, the horizontal SID route network of ML-PM system was designed, see Fig. 5.2. The left part is for North-outbound operations, the right part is for South-outbound operations. Each part models either segregated ML-PM departing topology or the mixed ML-PM departing topology. We use the South-outbound layout as the default layout. Flights departing to

Vertical flight efficiency

Here, we only consider the level-off on the sequencing leg in the descent. As shown in Fig. 5.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 more than 12% of the total flown time in the TMA.

Independent parallel operations with integrated arrivals and departures

Adjustment of optimization algorithm

First, for decision variables, we design three new decision variables for dynamically modifying departure trajectory. For each aircraft j ∈ F dep , we have the departure time t dep j , the turning time on the sequencing t T j , and the departing runway allocation r j . For the departure runway allocation, aircraft going to KM, SOSDI, RENOB will initially use runway 18R-36L to take off, and aircraft going to VY, CDY, TONIL, LADIX will initially use runway 18L-36R. The associated constraints for the variation of these decision variables are:

(5.7) 0% ≤ δt T j ≤ 100%, (5.8)

1 if j to KM, SOSDI, RENOB, 0 otherwise .

(5.9)

Appendix A

Appendix

A.1 Flight phases

There are six flight phases: Take-off, Departure, Route, Arrival, Approach, and Landing, see 1. Take-off. The aircraft needs to be guided as it taxis along the runway and ascends until completing the second segment of take-off, 400 feet up. This takes place in the CTR, a controlled airspace extending from the ground to a certain upper limit. If a failure happens before the decision speed is reached, the take-off will be aborted. If the failure is detected after the speed is reached, the take-off must continue.

2. Departure. After the second segment of take-off is complete (about 400 feet above the runway), the aircraft continues to climb, following the SID, until it reaches the flight level at which the route begins. This route takes place within the airspace structure called the CTA, and especially, in the TMA, where control services are provided. inform them about the ATFM measures that will be in force in European airspace on the following day via the publication of the agreed plan for the day of operations.

3. Tactical phase. The day of operations. We monitor and update the Daily Plan made the day before based on current reality. We continue working on capacity optimisation according to real time traffic demand, and where aircraft are affected by a regulation, we offer alternative solutions to minimize delays. Flights taking place on that day receive the benefit of the flow management service, which includes inter alia the allocation of individual aircraft departure slots, re-routings to avoid bottlenecks and alternative flight profiles in an attempt to maximize flight efficiency and make the best use of the available capacity.

A.3 PBN navigation specifications

Referring to ICAO-Doc.9613 (2008), PBN specifies that aircraft RNP and RNAV systems performance requirements be defined in terms of accuracy, integrity, availability, continuity, and functionality required for the proposed operations in the context of a particular airspace, when supported by the appropriate navigation infrastructure, see Tab. A.1.

A RNAV 1 means aircraft must maintain a total system error of not more than 1 nautical mile for 95% of the total flight time. A RNP 1 means that a navigation system must be able to calculate its position to within a square with a lateral dimension of 1 nautical miles. The key difference between them is the requirement for on-board performance monitoring and alerting, see Fig.A.2. A navigation specification that includes a requirement for on-board navigation performance monitoring and alerting is referred to as an RNP specification. One not having such a requirement is referred to as an RNAV specification. Therefore, if ATC radar monitoring is not provided, safe navigation in respect to terrain shall be self-monitored by the pilot and RNP shall be used instead of RNAV.

A.4 ISA

The ISA is an atmospheric model of how the pressure, temperature, density, and viscosity of the Earth's atmosphere change over a wide range of altitudes or elevations. It has been established to provide a common reference for temperature and pressure and consists of tables of values at various altitudes, plus some formulas by which those values were derived. The ISA is defined in ICAO Document 7488/2. The ISA assumes the mean sea level (MSL) conditions as given in Tab. A.2. Temperature decreases with altitude at a constant rate up to the tropopause:

The temperature remains at a constant value of -56.5 • C (216.65 • K) from the tropopause up to 20000m (65600f t).

Standard pressure p at a given altitude h is calculated by the function:

here the unit of T 0 is • K, and h is in meters.

A.5 Route network

The route network includes four tables. They are the list of nodes, see Tab. A.3, the list of links, see Tab. A.4, the list of routes, see Tab. A.5, and the route sets, see Tab. A.6.