UNIVERSITY OF BELGRADE FACULTY OF TRANSPORT AND TRAFFIC ENGINEERING Andrija D. Vidosavljević TRAJECTORY PLANNING ON PRE-TACTICAL AND TACTICAL LEVEL IN AIR TRAFFIC MANAGEMENT Doctoral Dissertation Belgrade, 2014 UNIVERZITET U BEOGRADU SAOBRAĆAJNI FAKULTET Andrija D. Vidosavljević UPRAVLJANJE PUTANJAMA VAZDUHOPLOVA U KONTROLI LETENJA NA PRE-TAKTIČKOM I TAKTIČKOM NIVOU doktorska disertacija Beograd, 2014 Thesis advisors: Dr Vojin Tošić Professor (retired), University of Belgrade – Faculty of Transport and Traffic Engineering, Belgrade, Serbia Dr Daniel Delahaye Professor, The French Civil Aviation University (ENAC), Toulouse, France Committee for evaluation and defence of the thesis: Dr Obrad Babić Professor, University of Belgrade – Faculty of Transport and Traffic Engineering, Belgrade, Serbia Dr Daniel Delahaye Professor, The French Civil Aviation University (ENAC), Toulouse, France Dr Feđa Netjasov Assistant Professor, University of Belgrade – Faculty of Transport and Traffic Engineering, Belgrade, Serbia Dr Radosav Jovanović Assistant Professor, University of Belgrade – Faculty of Transport and Traffic Engineering, Belgrade, Serbia Dr Mirjana Čangalović Professor, University of Belgrade – Faculty of Organizational Sciences, Belgrade, Serbia Public defence , 2014 Mentori doktorske disertacije: Dr Vojin Tošić Profesor u penziji, Univerzitet u Beogradu – Saobraćajni fakultet, Beograd, Srbija Dr Daniel Delahaye Redovni profesor, Francuski univerzitet za civilno vazduhoplovstvo (ENAC), Tuluz, Francuska Članovi komisije za ocenu i odbranu doktorske disertacije: Dr Obrad Babić Redovni profesor, Univerzitet u Beogradu – Saobraćajni fakultet, Beograd, Srbija Dr Daniel Delahaye Redovni profesor, Francuski univerzitet za civilno vazduhoplovstvo (ENAC), Tuluz, Francuska Dr Feđa Netjasov Docent, Univerzitet u Beogradu – Saobraćajni fakultet, Beograd, Srbija Dr Radosav Jovanović Docent, Univerzitet u Beogradu – Saobraćajni fakultet, Beograd, Srbija Dr Mirjana Čangalović Redovni profesor, Univerzitet u Beogradu – Fakultet organizacionih nauka, Beograd, Srbija Javna odbrana , 2014 This thesis has been supported by research grants of the Ministry of Science and Technological Development, Republic of Serbia, through projects of the Faculty of Transport and Traffic Engineering funded under the 2011- 2014 Research programs in technological development: Project TR36033, 2011-2014: “A support to sustainable development of the Republic of Serbia’s air transport system” To my father, who has always given me something to aspire. ACKNOWLEDGMENT First of all, I would like to express my special appreciation and thanks to my advisor, Prof. Vojin Tošić, for his guidance, inspiring discussions and valuable advice provided during the course of this research, as well as for his belief in me. He was a constant source of encouragement on both research and my career. I am also thankful to Prof. Daniel Delahaye, the thesis co-supervisor, who provides me valuable insight of application difficulties, optimization techniques and assured the necessary mathematic conscientiousness of this work. Without him, work on the thesis would be much harder. I would also like to thank my committee members, Ass. Prof. Fedja Netjasov, Ass. Prof. Radosav Jovanović and Prof. Mirjana Čangalović, for their insightful comments and helpful suggestions, especially Prof. Obrad Babić who made his valuable support available whenever I needed it. I am thankful to Chaimatanan Supatcha, whose programming code is used as a basis for the developed of optimization model. I owe my appreciation to other members of the Division of Airports and Air Traffic Safety: Tatjana Krstić Simić, Bojana Mirković, Nikola Ivanov and Marija Lazarević for their support and help in many occasions, brainstorming and friendship. I am also grateful to all members of Air Transport Department. I thank The French Civil Aviation University (ENAC) for providing traffic data and specially Cyril Allignol who assisted me manipulating and presenting data. Finally, I owe my deepest gratitude to my family, for their invaluable support during the whole course of my research. TRAJECTORY PLANNING ON PRE-TACTICAL AND TACTICAL LEVEL IN AIR TRAFFIC MANAGEMENT Summary: Global air traffic demand is continuously increasing, and it is predicted to be tripled by 2050. The need for increasing air traffic capacity motivates a shift of ATM towards Trajectory Based Operations (TBOs). This implies the possibility to design efficient congestion-free aircraft trajectories more in advance (pre-tactical, strategic level) reducing controller’s workload on tactical level. As consequence, controllers will be able to manage more flights. Current flow management practices in air traffic management (ATM) system shows that under the present system settings there are only timid demand management actions taken prior to the day of operation such as: slot allocation and strategic flow rerouting. But the choice of air route for a particular flight is seen as a commercial decision to be taken by airlines, given air traffic control constraints. This thesis investigates the potential of robust trajectory planning (considered as an additional demand management action) at pre-tactical level as a mean to alleviate the en-route congestion in airspace. Robust trajectory planning (RTP) involves generation of congestion-free trajectories with minimum operating cost taking into account uncertainty of trajectory prediction and unforeseen event. Although planned cost could be higher than of conventional models, adding robustness to schedules might reduce cost of disruptions and hopefully lead to reductions in operating cost. The most of existing trajectory planning models consider finding of conflict-free trajectories without taking into account uncertainty of trajectory prediction. It is shown in the thesis that in the case of traffic disturbances, it is better to have a robust solution otherwise newly generated congestion problems would be hard and costly to solve. This thesis introduces a novel approach for route generation (3D trajectory) based on homotopic feature of continuous functions. It is shown that this approach is capable of generating a large number of route shapes with a reasonable number of decision variables. Those shapes are then coupled with time dimension in order to create trajectories (4D). RTP problem is modeled as a mixed-variable optimization problem and it is solved using stochastic optimization technique, precisely Simulated Annealing (SA). The results indicate that, under certain conditions, with small increase of total planned costs, it is possible to increase robustness of the proposed solution providing a good alternative to the solutions given by an existing conflict-free trajectory planning models. Key words: 4D trajectory management, free-flight, robust trajectory planning, homotopy route generation model, flight de-confliction under uncertainty, flight interaction, decision cost, multi-objective, combinatorial optimization, simulated annealing Scientific field: Transport and Traffic Engineering Specific scientific field: Airports and Air Traffic Safety UDK: 656.7 (043.3) UPRAVLJANJE PUTANJAMA VAZDUHOPLOVA U KONTROLI LETENJA NA PRE-TAKTIČKOM I TAKTIČKOM NIVOU Rezime: Globalna potražnja za vazdušnim saobraćajem u stalnom je porastu i prognozira se da će broj letova biti utrostručen do 2050 godine. Potreba za povećanjem kapaciteta sistema vazdušnog saobraćaja motivisala je promene u sistemu upravljanja saobraćajnim tokovima u kome će u budućnosti centralnu ulogu imati putanje vazduhoplova tzv. “trajectory-based” koncept. Takav sistem omogućiće planiranje putanja vazduhoplova koje ne stvaraju zagušenja u sistemu na pre-taktičkom nivou i time smanjiti radno opterećenje kontrolora na taktičkom nivou. Kao posledica, kontrolor će moći da upravlja više letova nego u današnjem sistemu. Današnja praksa upravljanja saobraćajnim tokovima pokazuje da se mali broj upravljačkih akcija primenjuje pre dana obavljanja letova npr.: alokacija slotova poletanja i strateško upravljanje saobraćajnim tokovima. Međutim izbor putanje kojom će se odviti let posmatra se kao komercijalna odluka avio- kompanije (uz poštovanje postavljenih ograničenja od strane kontrole letenja) i stoga je ostavljen na izbor avio-kompaniji. Većina, do danas razvijenih, modela upravljanja putanjama vazduhoplova ima za cilj generisanje bez-konfliktnih putanja, ne uzimajući u obzir neizvesnost u poziciji vazduhoplova. U ovoj doktorskoj disertaciji ispitivano je planiranje robustnih putanja vazduhoplova (RTP) na pre-taktičkom nivou kao sredstvo ublažavanja zagušenja u vazdušnom prostoru . Robustno upravljanje putanjama vazduhoplova podrazumeva izbor putanja vazduhoplova sa minimalnim operativnim troškovima koje ne izazivaju zagušenja u vazdušnom prostoru u uslovima neizvesnosti buduđe pozicije vazduhoplova i nepredviđenih događaja. Iako predviđeni (planirani) operativni troškovi robustnih putanja mogu u startu biti veći od operativnih troškova bez-konfliktnih putanja, robusnost može uticati na smanjenje troškove poremećaja putanja jer ne zahteva dodatnu promenu putanja vazduhplova radi izbegavanja konfliktnih situacija na taktičkom nivou. To na kraju može dovesti i do smanjenja stvarnih operativnih troškova. U tezi je pokazano, da je u slučaju poremećaja saobraćaja bolje imati robustno rešenje (putanje), jer novo-nastali problem zagušenosti vazdušnog prostora je teško i skupo rešiti tj. zahteva značajnu promenu putanja vazduhoplova na taktičkom nivou. U ovoj tezi je uveden novi pristup za generisanje putanja vazduhoplova zasnovan na primeni homotopije neprekidnih funkcija. Ovim pristupom moguće je generisati veliki broj kriva, koje predstavljaju alternativne 3D trajektorije vazduhoplova, uz pomoć malog broja upravljačkih promenljiva (dve, tri). Potom se, primenom modela definisanog u tezi, svakoj tački dobijene 3D trajektorije pridodaje vreme u cilje dobijanja putanje vazduhoplova. RTP problem je modeliran kao problem matematičkog programiranja i rešen primenom tehnika stohastičke optimizacije, preciznije, simuliranog kaljenja. Rezultati pokazuju, da pod određenim uslovima, sa malim povećanjem ukupnih planiranih troškova moguće je znatno povećati robusnost predloženog rešenja, te razvijeni model predstavlja dobru alternativu konvencionalnim metodima za generisanje bez-konfliktnih putanja. Ključne reči: 4D upravljanje putanjama vazduhoplova, koncept slobodnog letenja, planiranje robustnih putanja vazduhoplova, model za generisanje putanja vazduhoplova primenom homotopije, bez-konfliktne putanje u uslovima neizvesnosti, interakcija između vazduhoplova, trošak upravljačkih akcija, više-kriterijumska funkcija cilja, kombinatrna optimizacija, simulirano kaljenje Naučna oblast: Saobraćajno inženjerstvo Uža naučna oblast: Aerodromi i bezbednost vazdušne plovidbe UDK broj: 656.7 (043.3) CONTENTS 1. INTRODUCTION ................................................................................................ 1 1.1. Motivation ....................................................................................................... 2 1.2. Current system and approach on solving congestion problem ............ 6 1.2.1. Air Traffic Flow Management ............................................................... 8 1.2.2. Review of literature in ATFM .............................................................. 12 1.2.2.1. Single airport ground holding problem - SAGHP ........................ 12 1.2.2.2. Multi-airport ground holding problem - MAGHP ....................... 13 1.2.2.3. Generalized traffic flow management problem - GTFMP ........... 15 1.2.2.4. Conclusion .......................................................................................... 20 1.2.3. Operational tools used in EU ATFM .................................................. 20 1.3. Future ATM ................................................................................................... 24 1.3.1. Trajectory modeling .............................................................................. 26 1.3.2. Previous research on Trajectory management .................................. 29 2. ROBUST TRAJECTORY PLANNING: PROBLEM AND METHOD DEFINITION ................................................................................................................ 33 2.1. Problem statement ....................................................................................... 33 2.2. Robustness and Flight interaction ............................................................ 35 2.2.1. Flight interaction .................................................................................... 36 2.2.2. Flight interaction metrics ...................................................................... 37 2.3. System boundaries and search space ....................................................... 41 2.4. Route generation model.............................................................................. 45 2.4.1. Homotopy ............................................................................................... 46 2.4.2. Mathematical formulation .................................................................... 48 2.4.3. Alternative route design ....................................................................... 50 2.5. Objective function ....................................................................................... 52 2.5.1. Cost of actions ........................................................................................ 54 2.6. Key model assumptions.............................................................................. 56 2.7. Mathematical formulation ......................................................................... 57 2.8. Complexity and Size of formulation ........................................................ 61 3. OPTIMIZATION TECHNIQUE ..................................................................... 63 3.1. Simulated annealing ................................................................................... 63 3.2. Application to our problem ....................................................................... 64 3.3. Flight interaction calculation ..................................................................... 66 4. RESULTS AND DISCUSSION ....................................................................... 68 4.1. General model settings and input data ................................................... 68 4.1.1. Fuel Burn Rate ........................................................................................ 68 4.1.2. Flight cost category ................................................................................ 69 4.1.3. Homotopy functions ............................................................................. 69 4.1.4. Flight interaction calculation ............................................................... 71 4.2. An academic test example .......................................................................... 71 4.2.1. Traffic data – Demand .......................................................................... 72 4.2.2. Test scenarios and numerical results of academic test example ..... 73 4.2.2.1. Base scenario ....................................................................................... 73 4.2.2.2. Low interaction cost scenarios ......................................................... 75 4.2.2.3. Constrained interaction scenarios ................................................... 78 4.2.2.4. Different aircraft type scenario ........................................................ 80 4.3. Real-world test example – French airspace ............................................. 81 4.3.1. Experimental setup ................................................................................ 81 4.3.1.1. Traffic sample ..................................................................................... 81 4.3.1.2. Traffic data source and format ......................................................... 82 4.3.1.3. Solution space ..................................................................................... 83 4.3.1.4. Interaction parameters settings ........................................................ 84 4.3.1.5. TMA modelling .................................................................................. 86 4.3.1.6. Scenarios settings ............................................................................... 87 4.3.2. Experimental results .............................................................................. 88 4.3.2.1. Conflict-free scenario ......................................................................... 88 4.3.2.2. Base scenario ....................................................................................... 90 4.3.3. Solution robustness testing .................................................................. 93 5. CONCLUSIONS ................................................................................................. 96 REFERENCES ............................................................................................................... 98 BIOGRAPHY .............................................................................................................. 108 LIST OF FIGURES Figure 1: Matching demand and capacity in European airspace ........................... 3 Figure 2: Departure delay by cause ............................................................................ 4 Figure 3: European wide horizontal en route flight efficiency ............................... 5 Figure 4: Breakdown of ATFM delay in September 2012 ........................................ 6 Figure 5: Operational Structure with CFMU Network Operations Unit ............ 22 Figure 6: ATFCM Activities ....................................................................................... 22 Figure 7: General functions in future ATM ............................................................. 25 Figure 8: The Business Trajectory Lifecycle ............................................................ 26 Figure 9: Lagrange interpolation result for a set of aligned points ...................... 27 Figure 10: Cubic Bezier curve .................................................................................... 28 Figure 11: Lunching rays – Light Propagation algorithm ..................................... 28 Figure 12: Principle of force field repulsive force ................................................... 31 Figure 13: Flights fully separated in time ................................................................. 37 Figure 14: Aircraft 3D position spaces ...................................................................... 37 Figure 15: Full 2D separation ..................................................................................... 37 Figure 16: Position time uncertainty ......................................................................... 38 Figure 17: Flight interaction functions ...................................................................... 39 Figure 18: Proximity/convergence referential for traffic complexity measure ... 40 Figure 19: Aircraft encounter on intersecting routes .............................................. 42 Figure 20: Time evolution of flights .......................................................................... 42 Figure 21: Flight delay ................................................................................................ 43 Figure 22: Flight rerouting.......................................................................................... 44 Figure 23: Speed management ................................................................................... 45 Figure 24: Example of two homotopic functions .................................................... 46 Figure 25: Route search space .................................................................................... 47 Figure 26: Vertical line test for a function ................................................................ 48 Figure 27: Symmetric homotopy with respect to reference function ............... 48 Figure 28: Alternative horizontal route .................................................................... 50 Figure 29: Horizontal route shape design and decoding ....................................... 51 Figure 30: Matching new route shape and given vertical profile ......................... 52 Figure 31: Specific Range relative to cruise speed .................................................. 55 Figure 32: Specific Range in function of cruise speed and altitude ...................... 56 Figure 33: Acceptance probability function ............................................................. 65 Figure 34: Classical trajectory representation .......................................................... 66 Figure 35: Four dimension space-time grid ............................................................. 67 Figure 36: Aircraft Fuel Burn Rates ........................................................................... 69 Figure 37: Boundary curves ....................................................................................... 70 Figure 38: One possible homotopy ........................................................................... 70 Figure 39: Resulting horizontal route ....................................................................... 70 Figure 40: Academic test example ............................................................................. 72 Figure 41: Roundabout test problem ........................................................................ 73 Figure 42: Base scenario - optimized trajectory in 2D ............................................ 74 Figure 43: Pareto frontier ............................................................................................ 76 Figure 44: Scenario 5 - optimized trajectory in 2D .................................................. 78 Figure 45: Set of solutions depending on required robustness ............................. 79 Figure 46: Distribution of flight cost categories ...................................................... 82 Figure 47: Input traffic data – French airspace ........................................................ 82 Figure 48: CATS file format ........................................................................................ 83 Figure 49: Alternative vertical profiles for one flight ............................................. 83 Figure 50: Flight interaction for initial (nominal) routes ....................................... 85 Figure 51: TMA airspace ............................................................................................. 86 Figure 52: Initial traffic demand at Paris TMA........................................................ 87 Figure 53: Conflict-free trajectories ........................................................................... 88 Figure 54: Traffic demand at Paris TMA for conflict-free trajectories ................. 88 Figure 55: Distribution of flight additional operating costs – Conflict-free scenario.......................................................................................................................... 89 Figure 56: Objective function convergence – Conflict-free scenario .................... 90 Figure 57: Robust flight trajectories .......................................................................... 91 Figure 58: Traffic demand at Paris TMA after optimization ................................. 91 Figure 59: Distribution of flight additional operating costs – Base scenario ...... 92 Figure 60: Objective function convergence – Base scenario .................................. 93 Figure 61: Distribution of conflicting points and flights ........................................ 95 Figure 62: Resulting congestion problem ................................................................. 95 LIST OF TABLES Table 1: Base scenario results ..................................................................................... 74 Table 2: Low interaction cost scenario results ......................................................... 75 Table 3: Low interaction cost – Scenario 5 results ................................................... 77 Table 4: Constrained interaction scenario results ................................................... 79 Table 5: Different aircraft type scenario results ...................................................... 80 Table 6: Test scenarios settings .................................................................................. 94 Table 7: Resulting congestion problem – conflicting points .................................. 94 Table 8: Resulting congestion problem – conflicting flights ................................. 94 1 1. INTRODUCTION Global air traffic demand is continuously increasing, and en-route congestion is cited as one of the principal restricting factor to future growth of the airline industry [Bertsimas and Stock Patterson, 2000]. The need for increasing air traffic capacity motivates a shift of ATM towards Trajectory Based Operations (TBOs). This implies the possibility to design efficient congestion-free aircraft trajectories more in advance (pre-tactical, strategic level). This thesis investigates the potential of robust trajectory planning (RTP) at pre-tactical level as a mean to alleviate en-route congestion in the airspace. A novel method for route generation based on homotopic feature of reference shapes is proposed. Centered on route generation model, a mathematical model is developed for robust trajectory planning that involves generation of congestion-free trajectories taking into account uncertainty of trajectory prediction and unforeseen events. In the remainder of this chapter the context of the research is set. It provides an overview of the problem of en-route congestion in European and USA airspace, in terms of magnitude and underlying causes. It then analyses the present way of dealing with this problem, indicating its principal limitations in this context. Operational systems currently in use at EUROCONTROL are also explained. A focused overview of literature in the field of Air Traffic Flow Management and Trajectory management is presented. In Chapter 2 the proposed model for robust trajectory planning is elaborated. After a statement of the proposed way forward, concept of flight interaction is presented as a measure of solution robustness. Several metrics of flight interaction are examined followed by system boundaries and solution space set-up. A technique for trajectory dimension space reduction using homotopic feature is investigated and proposed route generation model formulated. The specificities of the problem objectives are then analyzed, followed by the statement key model assumptions and the mathematical formulation itself. Finally, the complexity of the formulation is examined. Chapter 3 examines possible optimization techniques and explicates the choice of Simulated Annealing and specifics of built optimization model. To test the model, Chapter 4 describes a small hypothetical example that was developed first, illustrating the basic idea of the route generation model and 2 robust trajectory planning, and providing some valuable initial insights into model behavior. The model has been applied to a real life large-scale problem. Model results are examined and discussed. The advantage of RTP over conventional models is demonstrated by comparison of results of both models. Finally, Chapter 5 sums up the findings of this research, and cites the major contributions of the dissertation and point toward areas of future. 1.1. Motivation The continuous growth of the air transportation industry has placed an enormous strain on the Air Transportation System (ATS) and it has not been matched with an adequate increase of capacity (airport and airspace). In the late 80's, with a strong increase in air traffic, the total annual delay costs of European airlines due to congestion (including cost to passengers) were estimated to be $5 billion, and an average delay of all movement reached values of 5 minutes per movement (20 minutes per delayed movement) [Vranas et al., 1994]. Similar figures have been shown by United States airlines. This congestion crisis led to the decision in Europe to found the Central Flow Management Unit for coordinating air traffic flows across the EUROCONTROL’s Member States. European EUROCONTROL Central Flow Management Unit – CFMU (now Network Operation unit in the Network Management Directorate) and United States’ FAA Air Traffic Control Systems Command Center (ATCSCC former Central Flow Control Facility) have to assure that the traffic volume is compatible with the capacities declared by the appropriate ATS authority. Besides delays decrease in the first years of CFMU service (mainly between 1999-2003, red line Figure 1), delay is again rising with an air traffic increase [PRC, 2009, 2010, 2011]. Studies performed by the EUROCONTROL Experimental Centre showed that a 5% increase of traffic results in a 26% increase of delay [Bertsimas and Stock Patterson, 2000]. 3 Figure 1: Matching demand and capacity in European airspace [PRC, 2011] It is clear that the original forecast predicting three fold (lately reviewed to double) air traffic increase in 2020 compared to 2000 traffic [EC, 2001] won’t come true due to the global economic crisis and will be reached in 2050 instead, as shown by later studies [STATFOR, 2013]. However significant increases in demand expected to continue, and congestion is cited as one of the principal restricting factor to future growth of the airline industry. In 2012, about 17%1 of flights in Europe (9.55 million controlled flights being the annual total) arrived with more than 15 minutes delay compared with the schedule [PRC 2013]. IATA has estimated [IATA, 2013] that delays increased direct operating costs for European airlines by 4.5 billion Euros and an additional 6.7 billion for passengers. U.S. airlines have cited even larger numbers. Figure 2 shows departure punctuality and reveals that part of primary2 delay, not connected to turn around process, is almost equally caused by lack of airport and airspace capacity (including weather causes). One of the implications of both constraints at airport and in en-route airspace is that deriving good ATFM strategies has become a much more complicated task [Bertsimas et al., 2011]. 1 which is an improvement compared to 24% in 2010 with the same number of IFR flights 2 Primary delays are those which can be directly attributed to the reason of delay, whereas “reactionary” delays are the result of primary delays on earlier flight legs, which cannot be absorbed during the turnaround phase at the airport [PRC, 2005] 4 Figure 2: Departure delay by cause [PRC, 2013] EUROCONTROL’s initiative to enhance European flight efficiency through en route airspace design resulted in the implementation of free route airspace (FRA) starting from 2009. FRA refers to a specific portion of airspace where airspace users may freely plan their routes between an entry point and an exit point without reference to the fixed Air Traffic Services (ATS) route network [Eurocontrol, 2012a]. However flight remains subject to ATC at all times. Figure 3 represent the horizontal en-route flight inefficiency3 for the actual trajectory (red line) and the filed flight plan (blue line). The comparison of the annual values of flight inefficiency () shows a continuous improvement of flight efficiency due to more direct routing given by ATC in the tactical phase of flight. Savings are approximated to 1.3 million NM, representing the equivalent of 8 kiloton of fuel annually. Although the results are promising, local uncoordinated initiatives may not deliver the desired objective. With small airspaces, a large proportion of the observed inefficiency is due to the interface with adjacent airspace. This led to initiative for integration and enhanced cooperation between ACCs, forming Functional Airspace Block (FAB). FAB is defined [EC No 1070/2009] as airspace established regardless of State boundaries, where the provision of air navigation services and related functions is performed by an integrated provider. 3 An “inefficiency” is the difference between the length of the analyzed trajectory (great circle) and the achieved distance in percentage. Only portion between the departure and arrival terminal areas (radius of 40NM around airports) is considered [PRC, 2013]. 5 Figure 3: European wide horizontal en route flight efficiency [PRC, 2013] All these initiatives represent step forward to trajectory based operations and the Free Flight concept in the future. However, incremental changes of technology and procedures are no longer sufficient to keep up with the growth of traffic [Neal et al., 2011]. The International Civil Aviation Organization (ICAO) has developed the Global Air Navigation Operational Concept in response to this problem. This represents a fundamental change in the operating paradigm for air navigation services [ICAO 9854]. Major systems development programs are underway around the world including the Next Generation Air Transportation System (NextGen) program in the United States [FAA, 2012], and the Single European Sky ATM Research (SESAR) program in Europe [SESAR, 2013]. The ICAO Global Operational Concept envisages that the primary means by which the priorities such as safety, efficiency, and cost-effectiveness are to be achieved is via trajectory management. Also it proposes that airspace users “should retain primary responsibility for the provision of Conflict management (assuring separation between aircraft)”4. This implies the possibility to design efficient aircraft trajectories that are conflict-free at strategic level. Pre-tactical conflict management will reduce the need for separation provision (and collision avoidance) to a designated level. However, special events such as severe weather, volcanic ash, ATC strikes etc., represent large influencing factor to flight 4 However, the Operational Concept also states that the allocation of responsibilities is subject to the design of the ATM system. 6 inefficiency (Figure 4) today. This is not expected to change in the future, and therefore it is necessary to design robust trajectories that could alleviate the effects of an unforeseen event. Figure 4: Breakdown of ATFM delay in September 2012 [DNM, 2012] 1.2. Current system and approach on solving congestion problem Air Transportation System (ATS) has three main elements: infrastructure, vehicles (airplanes) and operations. Infrastructure refers to structure that allow an airplane to operate, and it consists of a routed network for which node are airports or beacons and links are airways. Airline companies are in charge of operating airplanes (public transport). General and Military aviation are another category of ATS users. Air Traffic Management (ATM) aims at ensuring safe, efficient, and expeditious aircraft movement and it is composed of the following services [Babic and Netjasov, 2011]:  Air Space Management (ASM) – includes airspace design/modeling which aims to better adapt ATC capacity (route/sector capacity, arrival/departure airport capacity) to continuously changing demand.  Air Traffic Flow Management (ATFM) – aims to optimally use available capacity in order to minimize delay and cost to airspace users. In Europe it has evolved into the new concept of Air Traffic Flow and Capacity Management (ATFCM).  Air Traffic Service – is service which regulate and assists aircraft once airborne and consists of: flight information service (FIS), air traffic advisory service (AS) and air traffic control (ATC). ATC is responsible for safe flow of air traffic. Air Traffic Controller Officer (ATCO) monitors the traffic, ensures minimum safety separation between aircraft at all time, and in case of separation violation take measures for conflict resolution (vectoring, altitude change, and speed regulation). Additional ATCO role 7 is to keep aircraft separated from ground and other aerial vehicles, out of restricted areas etc. When flying from origin to destination airport, an aircraft must follow ATS routes and beacons on the airways network as stated in its flight plan (FPL). When two adjacent aircraft routes converge to the same point, there exists a risk of collision. Therefore, the main task of ATCO is to ensure conflict free trajectories, as stated above. As human beings, ATCO has working limits, which is the main factor that determines the airspace capacity. The airspace is then partitioned into sectors, each of them being assigned to a team of controllers. On the other hand airport capacity is mainly related to runway capacity, defined as the number of aircraft that can land/take-off at a given period of time (arrival/departure capacity). Factors influencing airport runway capacity are: number and configuration of runways, aircraft mix and operating sequence, approach speed, etc. Both, airport and airspace capacities, are strongly influenced by weather. With increasing air traffic demand and/or capacity degradation, some parts of the air traffic system reach their operational limits and become congested generating delays. Solutions for traffic congestion problem vary according to considered time horizon:  Long-term approaches (years in advance) include infrastructure improvements such as: construction of additional runway, improving ATC technologies and procedures, airspace sectorisation, etc.  Medium-term approaches (several months up to weeks in advance) consider modification of traffic flows in order to reduce peak demands and adapt ATC capacity to demand. (Strategic Flow Management).  Short-term approaches (several days up-to day of operation). Pre-tactical Flow Management (days in advance) attempts to optimally use available capacity, minimizing costs to airspace users. To cope with congestion, or avoid bad-weather and restricted areas (no-fly zones) some aircraft trajectories are then recomputed and ATFM regulations are applied. Short-term solutions also include Tactical Flow Management occurring on the day of operations which aims to increase usability of available capacity as much as possible, and Tactical Control which solves the potential conflict in real time. Neither of last two solutions considers solution costs due to higher goals. 8 The level at which solution is applied dictates the uncertainty involved and the room for maneuver [Matos et al., 2001]. At the strategic level, there is scope for significant changes in the routing of flows, but information available on traffic demand a few months ahead are very poor. As the day of operations gets closer the information available becomes more accurate but the room for introducing changes to the routings diminishes, decreasing optimality of the solution. 1.2.1. Air Traffic Flow Management In the most general sense ATM can be viewed as a system supplying service with a finite number of elements [Tosic et al., 1995]. Each element (flight routes, sectors…) has limited capacity and is shared among several customers (flights). Therefore it is possible that total demand exceeds some element’s capacity at certain periods, and some customer may suffer delays or may receive control actions (assigned to different element). The flow control activity can be decomposed into two different phases [Bianco and Bielli, 1992]: congestion forecast and congestion prevention. The first phase requires accurate evaluation of air traffic capacity and demand, what is usually hard to forecast because of numerous influencing factors that vary with space and time. Once capacity and demand are evaluated, congestion prevention takes place and may involve control actions to keep delays and overall costs as low as possible. Such problem is called Air Traffic Flow Management Problem (ATFMP) in literature. Closely speaking of ATFMP the question is how to optimally control flights in the ATS network (assign route, speed and slot of departure) in order to minimize total costs (fuel, time, safety, passengers compensation costs, etc.) induced to users (airlines) with respect to the given system constraints. Mathematical formulation depends on the choice and meaning of the decision variable and the particular problem that model intends to solve (objective function and constraint definition). One general model, adapted from [Tosic et al., 1997], is here presented with the following notations used: Total number of capacitated elements Total number of time periods set of all flights 9 set of flight in that have alternative routes ( ) initial route for flight in set of initial routes for flights in , ⋃ set of alternative routes for flight in set of all alternative routes for flights from , ⋃ maximum allowed delay of an flight using route input binary variable equal to 1 if flight using route with delay demands service at element during period , and 0 otherwise ( , ̅̅ ̅̅ ̅̅ , ̅̅ ̅̅ ̅, t ̅̅ ̅̅ ̅ ) total capacity of element at period ( ̅̅ ̅̅ ̅, t ̅̅ ̅̅ ̅ ) ( ) total ground delay cost of flight using route delayed for permitted number of periods total additional constant delay cost of flight assigned to the alternative route . Equal to 0 if . decision binary variable equal to 1 if flight using route is delayed periods, and 0 otherwise ( , ̅̅ ̅̅ ̅̅ ) ∑ ∑( ( ) ) (1) subject to: ∑ ,  (2) ∑ ∑ ∑ , ( ) (3) ∑ ∑ ( ) ̅̅ ̅̅ ̅, t ̅̅ ̅̅ ̅ (4) { } ( ) , ̅̅ ̅̅ ̅̅ (5) The objective function (1) represents the sum of the total ground delay assigned to flights and additional costs of flights assigned to alternative routes. Constraints (2) and (3) ensure that each flight is assigned to one route only and is given one delay. It is an assignment constraint that guarantees that each flight occupies only one element during a given time period. The capacity constraint (4) expresses the requirement that total demand should not exceed the capacity at any 10 element and for all periods. Finally integrality constraint (5) ensures that decision variables are of binary type. In addition, if consecutive flights5 are modeled, a coupling constraint is necessary to ensure that takeoff time for a consecutive flight is scheduled after the landing time of a preceding flight. In the case of hub airport with banks of connected flights, coupling constrains ensure that all transfer passengers reach their connected flight. Tosic and Babic in [Tosic and Babic, 1995] discussed several metrics to be used as cost function ( ) (for delay measurements for their particular case) when solving the AFTMP:  Basic case where all flights that received controlled actions are treated equally without taking into account the magnitude of action. The objective is then to minimize the total number of flights which deviate from their initial-nominal route.  Using a single linear cost function for all flights will minimize the total amount of actions received (delayed minutes/additional flight miles).  In reality, unit costs differ from one flight to another because of: aircraft size (the bigger plane with more passengers carried, the greater operating costs are, etc.), flight type (feeder6 or non-feeder flight, cargo flight), etc. Using individual unit costs will minimize overall operating costs taking into account flight priority. On the other hand, use of individual unit costs especially in formulations that use dynamic landing priority rule, such as landing according to decreasing marginal cost, will raise equity (fairness) issue because of bias against “cheap” flights7 [Terrab, 1990], [Terrab and Paulose, 1993].  When the traffic assignment is controlled by a central authority (CFMU, ATCSCC), aiming to find a solution that is most desirable to the society (system- optimum solution), fairness is a very important issue and it could be incorporated into objective function using a non-linear cost function. In that way, larger control 5 Flights that are operated by the same aircraft on a given day 6 Flight that brings-in (feeds) passengers from small regional destinations to hub airports for onward journey (usually international or even intercontinental flights). Read more: http://www.businessdictionary.com/definition/feeder-airline.html#ixzz2vwqyZutz 7 Flights that have lower unit cost will be penalized more than others. 11 action will be penalized assuring balanced distribution of control actions among flights/airlines at the expense of higher total operating costs. Optimization models used to solve ATFMP may be classified according to different criteria:  Deterministic vs. Stochastic models – which are distinguished by whether the ATC capacity is assumed deterministic or probabilistic (presence of uncertainty). Although it is fairly common to consider airport and sector capacities as deterministic values, they are highly variable in reality and may change dramatically over time. When decisions are made under uncertainty one must consider the trade-off between possibilities to assign excessive inexpensive control action (e.g. ground holds) or to experience more expensive ones (e.g. holding at destination airport).  Static vs. Dynamic models – which are distinguished by whether or not the solutions are updated during time horizon. In static models control actions are chosen once at the beginning of the first time period. On the other hand with dynamic models, control actions are revised over time based on updated inputs (e.g. ATC capacities). Although dynamic models could yield better overall solution in the case of uncertainty, they will tend to penalize short-haul flights as shown in [Terrab, 1990], [Terrab and Paulose, 1993], [Bertsimas and Odoni, 1997].  Network optimization (Minimum cost flow [Terrab and Odoni, 1993], [Helme, 1992], Maximum flow problem [Bianco and Bielli, 1992], Assignment problem [Lindsay, et. al., 1993], Multi-commodity problem [Bianco and Bielli, 1992], [Lindsay, et. al., 1993], etc.), Dynamic programming (DP) [Andreatta and Romanin-Jacur, 1987], Linear (integer, mixed) programming (LP, IP, MILP) [Vranas et al, 1992], [Bertsimas and Stock Patterson, 1998], Non-linear programming (NLP) [Zenios, 1991], Stochastic programming (SLP) [Richetta and Odoni, 1994], Heuristic approach [Delahaye et al., 1994], [Delahaye and Odoni, 1997], [Tosic et al., 1997] - based on the optimization method (techniques) used.  Ground-holding problem (GHP) [Andreatta and Romanin-Jacur, 1987], [Terrab and Odoni, 1993], [Vranas et al, 1992], Generalized traffic flow management problem (GTFMP) [Lindsay, et. al., 1993], [Delahaye et al., 1994], [Tosic et al., 1995], [Bertsimas and Stock Patterson, 1998], Traffic flow management rerouting problem (TFMRP) [Tosic et al., 1997], [Delahaye and Odoni, 1997], 12 [Bertsimas et al., 2008] and their subclasses – according to the type of problem they address. 1.2.2. Review of literature in ATFM In the literature, problem of demand-capacity imbalance is addressed since 80's. Early efforts in the area of ATFMP research were made by the Italian National Research Council (CNR) [Tosic and Babic, 1995] and complete research results are summarized in [Bianco and Bielli, 1992]. Even though the literature has adopted the view that Andreatta in 1986 and later Odoni in 1987 first formalized and introduced ATFMP. Later paragraphs present some of the published research relevant to ATFMP in last 20 years. 1.2.2.1. Single airport ground holding problem - SAGHP Models and algorithms developed in 80's and early-90’s focus mainly on airport congestion. First models considered finding optimal aircraft release times (ground-holds assignment) in a very simplified single-airport configuration network - SAGHP. Andreatta and Romanin-Jacur in [Andreatta and Romanin- Jacur, 1987] consider a stochastic version of a single-airport problem for a single time period. They assumed that congestion may arise only at the arrival airport during a given single period and modeled airport capacity as a random variable. They used Dynamic Programming (DP) to find an optimal ground holding policy (ground-holds assignment) for a given landing priority rule. Terrab in his PhD thesis [Terrab, 1990] and [Terrab and Odoni, 1993] presented an extension of a single-airport problem for multiple periods. Deterministic version of the SAGHP (capacity of arrival airport is a deterministic function known in advance) is formulated as an Integer Program (IP). As the constraint matrix of this formulation is totally-unimodular, a simplex method is used to solve the relaxed version of the problem. Network-type algorithms were even more effective in this approach. A static stochastic extension of the problem was solved with DP whose time complexity highly depends on the number of time periods and requires enormous computational effort. A solution approaches (heuristic) to the dynamic version of the SAGHP, built on static algorithms, were also proposed. In order to take advantage of additional weather information, that becomes available as time passes, Richetta in [Richetta, 1991] and [Richetta and Odoni, 1994] developed an optimal stochastic and dynamic model of SAGHP. To reduce the size of the 13 problem, flights are classified into small number of different classes. Flights that have identical ground-holding delay cost and are scheduled to arrive at the destination airport at the same period of time belong to the same class. A stochastic linear programming is used to efficiently solve a real instance of the problem at the busiest airport in Massachusetts, Boston’s Logan. Opposed to Andreatta’s or Terrab’s formulation, Richetta used dynamic landing priority rule such as FCFS, landing according to the decreasing marginal cost, etc. Terrab and Paulose in [Terrab and Paulose, 1993] proposed a fully dynamic formulation (control actions are subject to revision in the latter stages) of the single-airport problem but only applied with two-stage decision making. In [Hoffman, 1997], [Ball et al., 1999] and [Ball et al., 2003] another IP formulation of static and stochastic SAGHP is proposed. Although delays are assigned to individual flight, decision variables are aggregated to some extent and represents the number of flights selected to land during a given period i.e. planned arrival capacity. This formulation gives several advantages: first, it treats individual flight that is consistent with Collaborative Decision-Making (CDM) approach; second, pre- and post-processing are possible to reduce the number of decision variables by an order of magnitude; third and the most important, the constraint matrix associated with this IP is totally unimodular, and the LP relaxation yields integer solutions (proved by the authors). A dynamic stochastic optimization based approach is presented in [Mukherjee, 2004] and [Mukherjee and Hansen, 2007] for SAGHP, where ground delays assigned to flights can be revised during different decision stages, based on weather forecasts. 1.2.2.2. Multi-airport ground holding problem - MAGHP Usually an aircraft performs more than one flight on a given day and when a specific aircraft is delayed, in many cases, the next flights (performed by the same aircraft) will also be delayed. In order to capture “network” effect, multi- airport ground holding problem has been formulated - MAGHP. Static deterministic formulation of the problem was addressed for the first time by Vranas, Bertsimas and Odoni in [Vranas et al, 1992] and [Vranas et al, 1994a]. They considered the network of the K-busiest airports and the associated flights departing from and arriving to an airport of K. They presented three models with 0-1 Integer formulation: (1) with departure and arrival capacity, (2) only with arrival capacity and (3) with arrival capacity and flight cancellation option. Flight 14 cancellation is introduced because the existence of delay’s upper bound could yield unfeasible solution. They proposed optimal solution as well as heuristic approach based on the LP relaxation of IP. Heuristics first collect a set of "critical" flights in LP (i.e. flights for which some integer constraints are violated) and gives a "rounding" scheme for those flights which leaves unchanged, as far as possible, the remaining flights. Same authors presented a dynamic extension of MAGHP in [Vranas et al, 1994b]. Dynamic and deterministic extension was also formulated as IP and efficiently solved. In dynamic and stochastic formulation, as airborne delay could not be totally avoided (depending on capacity scenario that is realized), ground delay could not be expressed in terms of landing assignment decision (equals 1 if flight land at given time). Those two decision variables are taken into account independently. Moreover, ground delay is not a binary variable. Because of the size of the formulation, heuristic approach was developed but was highly inefficient as stated by the authors. Andreatta and Brunetta in [Andreatta and Brunetta, 1998] proposed slightly different formulation of static and deterministic version of the MAGHP compared to Vranas et al. They model consecutive flights with one set of binary variables. Although the number of variables becomes very large (all pairs of flights) it is not dependent on the number of consecutive flights operated by an aircraft. In addition coupling constraints could be omitted in this formulation. Andreatta and Brunetta in [Andreatta and Brunetta, 1998] have performed a comparison of computational performance of three deterministic models for MAGHP: two models mentioned above and a model described in [Bertsimas and Stock Patterson, 1998]. The latter model is formulated as BIP with decision variables representing the arrival of a flight at a given point (could be airspace or airport) “by” the given time. As it is a general model of ATFMP it is explained in details in the following chapter. The authors concluded that Bertsimas and Stock model performs best in most cases especially due to strong formulation that yield integer solution of LP relaxation. Most of the models for MAGHP considered only arrival airport capacity as the main bottleneck taking into account that it takes longer for an aircraft to arrive than to depart. This is however not completely truth, as some runways are used for both operations at an airport. A way to extend models in order to incorporate the dependence between arrival and departure capacity was first introduced in [Vranas et al, 1994b]. It was shown that no additional variables are required, just 15 few more constraints. This was confirmed in [Bertsimas and Stock Patterson, 1998] using similar formulation. 1.2.2.3. Generalized traffic flow management problem - GTFMP In mid-90's, it has become increasingly evident that more delays are caused by airspace congestion especially in Europe. As noted in [Tosic and Babic, 1995] crucial difference between possible solution for ATFMP with respect to airport and airspace network is the fact that first could be solved using only temporal solution while second one requires both temporal and spatial solution. This is because the solution of airspace congestion requires enormous amount of ground- holds to be assigned which is very expensive and even not feasible (with constrained delay). Therefore, other options to resolve congestion by controlling a flight throughout its duration were introduced: speed management, rerouting, altitude change, etc. Zenios in [Zenios, 1991] propose a model to optimize the traffic flow and congestion (considering cost along with risk - dual objectives) at the target control sector by controlling ground-holds and flight level. Airspace is modeled as a multiple-period network with identical speeds assumed for all aircraft, one-class (commodity). The resulting nonlinear network problem was solved using numerical optimization technique. First, multi-period multi- commodity (multi-class) network type model was reported in [Bianco and Bielli, 1992] based on work of Bianco and Bielli in 1981 and 1982. The airspace structure was represented with direct multigraph (or pseudograph) where different arcs are associated to different flight levels along the same route. The optimal strategy of the traffic congestion problem was formalized as a maximum flow problem on the network: maximize source-generated flows so as to better match traffic demand with system capacity. The same authors proposed a second model to optimize flight plans of each aircraft by control actions such as altitude variation, speed control and departure delays. The objective was to minimize total delay and deviation from the nominal flight profile for all planes in the control region. Because of problem complexity and possible conflict with ATC procedures, flight planning was carried out according to FIFO (flight by flight) and the problem was decomposed and solved in sequential steps. Helme in [Helme, 1992] modeled airspace structure as the Space-Time Network very similarly to [Bianco and Bielli, 1992]. The problem was formulated as minimum cost flow with capacitated links. Hence airports were modeled with one additional node and artificial link between 16 them. Single destination problem (single sink in the network) is efficiently solved. Multi-destination problem, however, does not conform to the pure minimum cost flow structure. This is because every sink in the network is equally treated and flow demands between source-sink pairs cannot be preserved (“who goes where” is ignored). Therefore multi-destination problem is formulated as a multi- commodity minimum-cost flow on a network with commodity representing all flights going to a given airport. Implicit methods for handling capacity uncertainty as well as hubbing are proposed. While formulation of this model is straightforward, its computational performance were poor as referred in [Bertsimas and Odoni, 1997]. Time Assignment model of GTFMP presented in [Lindsay, et. al., 1993] is formulated as a Binary Integer Programming problem. The objective was to minimize the sum of ground and air delay costs for selected flights taking into account departure and arrival airport and airspace (sectors) capacities. Decision variables were times at which aircraft are to be at a selected fixes in their routes. Therefore it is possible to delay flights on the ground (by ground-holds) and in the air (by speed reduction, etc.). The problem was solved efficiently using custom developed software built upon the CPlex8 consisting of problem reduction techniques which are executed before CPlex optimization. Delahaye et al. in [Delahaye et al., 1994] considered controller workload in the traffic assignment problem rather than using capacitated airspace. The objective of the study was to minimize transportation costs and controller workload, where transportation cost included cost of distance travelled and congestion on the route and controller’s workload included monitoring, conflict and coordination workload. Opposite to previously mentioned studies on network flows, the authors assumed non-segmented origin-destination flows as equity constraint. They used genetic algorithm to solve this complex problem. Tosic et al. in [Tosic et al., 1995] proposed a discrete model of ATFMP with en-route capacities formulated as Binary Integer Programming with ground delays as decision variables. The model was classified as multi-element, multi-period, deterministic and static. Emphasized by authors, proposed model stays linear even when delay cost is not linear, because nonlinear dependence is expressed through the objective function. Several heuristic algorithms based on lexicographic approach have been developed. The idea was to set as many flights without delay, then sequentially 8 CPlex Mixed Integer Optimizer was developed by CPLEX Optimization Inc. and offered commercially starting in 1988. Today it is actively developed under IBM named IBM ILOG CPLEX Optimization Studio. 17 maximize the number of the remaining flights served with 1, 2 …Dmax periods of delay. Same authors proposed extension of the model to a spatial solution in [Tosic et al., 1997] using the set of alternative routes. This was one of the first studies that consider horizontal rerouting of flights. The objective was to minimize delay cost to the customers incurred both on the ground and in the air (alternative route). Delahaye and Odoni in [Delahaye and Odoni, 1997] proposed stochastic optimization model of general time-route assignment problem with the objective to significantly reduce the peak of workload in the most congested sectors and in the most congested airports. The model has ability to manage constraint of the airlines in a microscopic way (at the flight scale e.g. flight prioritization) and to take into account the connexion9 between flights. Promising results of hypothetical example were reported. Oussedik and Delahaye in [Oussedik and Delahaye, 1998] proposed similar GAs based model as previously reported by Delahaye and Odoni. The model was tested on real day traffic data across French airspace. The computation time was reported as the weak point of this model. First of the models reported by Bertsimas and Stock Patterson in [Bertsimas and Stock Patterson, 1998] was variant of Time Assignment model reported by [Lindsay, et. al., 1993]. Although decision variables in both models represent the times at which aircraft are to be at selected fixes in their routes, the difference lies in the formulation of decision variable. Binary decision variable in latter model equal to 1 if flight arrived in the sector by the given time i.e. at or before a given time. Such decision variable makes Binary Integer Programming formulation of the problem very efficient. Therefore the solution of the LP relaxation of this problem is almost always integer as reported by the authors10. This statement is not proved but it is confirmed by numerous experiments. As a result, the model is capable of solving large, realistic size problems in a reasonable amount of time. Most authors gave a credit that this is the model with the best performance among existing ones. Bertsimas and Stock Patterson examine several variations of the model taking into account: interdependency between arrival and departure capacities at airports, hubbing effect, multiple presequent flights11 and flight 9 Connected flights from the perspective of the passengers. Therefore one departure flight may be connected to many arrival flights. 10 Polyhedral structure of problem formulation was examine by the authors 11 At hub airports many airplanes are capable of performing any one of multiple consecutive flights. 18 rerouting. They show two possible approaches to model rerouting decision: “path approach” and “sector approach”. In the path approach, a set of possible routes that a flight may fly is defined. Decision variable is further extended and consider arrival of the flight at the sector by the time along the specified route. In the sector approach, route is defined on the run i.e. at each sector in its route it is decided which sector to enter next. Therefore, set of sectors that flight can enter after a given sector is defined and decision variable consider arrival time of the flight at the sector from specified sector. Although no extensive study was performed to determine which approach gives a better solution in terms of integer solution, it is clear that sector approach formulation result in more variables. In [Bertsimas et al., 2008] and [Bertsimas et al., 2011] sector approach is further studied with the scope of strengthening the formulation. Key achievement is that authors included rerouting decision without requiring any additional variable compared to [Bertsimas and Stock Patterson, 1998] model but include additional inequalities. Realistic problem instances of NAS level were solved in short computational time. Bertsimas and Stock Patterson in [Bertsimas and Stock Patterson, 2000] examine a problem of dynamical rerouting. They presented an aggregated model (dealing with flows) that was formulated as a dynamic, multi-commodity integer network flow problem. The network model is a mixture of the ones proposed in [Bianco and Bielli, 1992] and [Helme, 1992] however it is not transformed into static multi- period problem as the previous author did and therefore involves complicated side constraints. The commodities in the network are defined as origin-destination pairs of airports operated by an airline. With such definition issue of fairness become even more important as it is possible to assign all delay to a single airline. To find near-optimal solution, the authors proposed iterative algorithm were in a first step, Lagrangian relaxation of the problem is solved by LP giving objective lower bound, then randomized rounding heuristic is performed to round variables and decompose flows into routes, finally routes are packed into capacitated airspace system solving integer packing problem and upper bound is then updated. The procedure is repeated until upper and lower bound reaches desired accuracy. The computational results suggested that this approach is capable of efficiently solving real life problems. Leal de Matos et al. in [Matos et al., 2001] examine several optimization models for flight rerouting applicable to European flow management authority. According to current practice alternative routes are allocated on flight flows and not on an individual flight basis. The problem was solved in two stages: (1) route problem, where acceptable alternative routes for 19 each flow is identified, and (2) assignment problem, where the route is assigned to each flow so that the total cost of rerouting and congestion is minimized. They take congestion into consideration: using penalties whenever demand exceeds capacity of a sector, or constraining the demand. Three IP models are presented depending how congestion is modeled. The smallest and fastest model considers congestion penalties and only rerouting decision while two others consider limited sector capacities and ground-delay decision variable at different aggregation level. The work presents a nice overview of modeling approaches in relation to current practice and their usability but lacks of innovation. In [Dell’Olmo and Lulli, 2003] two-level hierarchical architecture for ATFMP is proposed: network airways system and single airway. At the system level dynamic multi-commodity (each flight represents commodity) network flow problem with side constraints was formulated and solved. Analyzing the network, the flight sequence for each way point is outputted. Then optimization of the lower level model is carried out on those arcs of the network where the flows exceed capacity. As a result of this model, the flight paths in the free flight scenario are obtained. Although authors considered ATM as a network, there were no fix routes i.e. a flight could follow any link of the network considering arrival time as the only constraint. The objective function was to minimize a cost function including: late and early arrival and deviation from its preferred speed. Due to nonlinearity, the computational time was very high, even when constraints such as sectors capacity, continued flights, etc. are not considered. The model is also validated using heuristic techniques. Lulli and Odoni in [Lulli and Odoni, 2007] presented an approximate, low-level-of detail deterministic model for the EU ATFM environment in order to highlight the critical characteristics of the system: (a) the intrinsic complexity of optimal ground and airborne holding strategies; (b) the fundamental conflict that may arise between the objectives of efficiency and equity. The model takes into consideration airport and en-route capacity constraints and assigns both ground and airborne delays (holding in terminal area) to flights with the objective to minimize total delay costs. The model can be viewed as macroscopic in that it omits certain details making assumptions such as: equal speed of travel, no rerouting, no consecutive flights, etc. The authors noted that even with perfect knowledge of future capacity (deterministic case), in many cases, the total delay cost can be reduced significantly by assigning more expensive airborne holding delays to selected flights (instead of ground delays). This is in sharp contrast to the situation in which only airport capacity constraints 20 exist. Richard, Constans and Fondacci in [Richard at al., 2011] proposed a MIP model of ATFMP. They consider dynamic search of routes during the sliding horizon loop process taking all flight characteristics into account. The problem was solved using branch-and-price and column generation techniques. On master level restricted problem with only a small subset of variables is solved and then this subset iteratively enlarged by adding some variables with negative reduced cost. In this context, columns to be generated correspond to feasible trajectories. Interesting columns are characterized by a negative reduced cost on the basis of the last optimal relaxed solution. Recently [Gupta and Bertsimas, 2011] and [Agustin et al., 2012] proposed the first models that deals with airspace uncertainty. 1.2.2.4. Conclusion Transportation network represents a well-adapted modeling framework for most processes that might be encountered in reality. Unfortunately, some of these problems are so complex that no transportation network algorithms are available to solve them [Delahaye, Puechmorel, 2013], neither they take account of equity between users. Dynamic and linear programming models have been proved to be efficient on simplified formulations but with more realistic formulation (route, slot and speed allocation, arrival, departure and space capacity, etc.) their use could be limited only to small-scale networks. Heuristic approaches are shown most suitable for solving large-scale problems. Moreover, the multi-objective approach could be used to identify several potential solutions, which can then be assessed by an expert. Nevertheless, none of the reviewed models are adapted for trajectory design (shape planning), that will be the corner-stone of the future ATM system as it will elaborated in the chapter 1.3. 1.2.3. Operational tools used in EU ATFM EUROCONTROL Central Flow Management Unit (CFMU) is mandated to provide Air Traffic Flow and Capacity Management services within the area of responsibility of participating European States. It was founded on behalf of Transport Ministers of the European Civil Aviation Conference (ECAC) States in 21 October 1988. The CFMU is based on the ICAO Centralised Traffic Management Organisation (CTMO) [ICAO 4444] concept and integrates Flow Management Positions (FMPs) in each Area Control Centre (ACC). From January 2011 with EUROCONTROL’s new organization, CFMU Directorate was integrated into a wider EUROCONTROL unit the Directorate Network Management (DNM) and has been renamed to CFMU Network Operations. When an imbalance between demand and available capacity is detected, CFMU NM in close cooperation with concerning ANSP's may develop a number of procedures to avoid congestion such as: the reconfiguration of some sectors, the activation of mandatory routes for certain trajectories, and the creation of slot allocation regulations [Matos and Omerod, 2000]. CFMU Network Operations has the following objectives:  To develop and maintain the highest level of quality of ATFCM service on behalf of both its ATS and AO users, as principal objective. It includes: - For ATS - The provision of flight plan data, the best utilization of available capacity, the smoothing of traffic flows and the assurance of protection against overloads. - For AOs - The provision of advice and assistance in flight planning and the minimization of penalties due to congestion.  Maintain and improve the cost effectiveness of its operations.  Adapt its procedures and systems to the evolution of its operational environment (Single European Sky (SES) initiative).  Provide reports and statistics on ATFCM operations and delay situation. The role of the CFMU Network Operations Unit can be illustrated by the main processes for exchanging operational information across the Network, as described in Figure 5. 22 Figure 5: Operational Structure with CFMU Network Operations Unit [Eurocontrol, 2011] Air Traffic Flow and Capacity Management is a service that is enhancing ATFM with the objective of managing the balance of demand and capacity by optimizing the use of available resources in order to enhance the quality of service and the performance of the ATM system. It consists of three activities as shown on Figure 6:  Strategic Flow and Capacity Planning;  Optimised Capacity Management;  Tactical Flow and Capacity Management. Figure 6: ATFCM Activities [Eurocontrol, 2004] In the strategic phase, based on traffic and capacity forecast, bottlenecks are identified and modifications of traffic flows are considered to adapt ATC capacity to demand. During pre-tactical phase, rerouting scenarios are examined in order 23 to achieve a global decrease of delays by spreading traffic. They can be in a form of advices or mandatory instructions when there is a risk of major imbalance between demand and capacity. On the day of operation when regulations are activated, the slot allocation procedures are applied. It includes the calculation of Calculated Take-Off Time (CTOT) for each flight passing thought the restricted area. Delays are carefully monitored and possible flights that would benefit from reroute are identified. Alternative route can be proposed automatically by Enhanced Tactical Flow Management System (ETFMS) or rerouting can be initiated by Aircraft Operators (AOs). Computer Assisted Slot Allocation (CASA) is the tool used at EUROCONTROL to calculate CTOT based on allocated slot. It is based on a heuristic model for solving GTFMP. Priority of flights is given according the “First Planned – First Served” principle, meaning that flights are sequenced in the order they arrive at a regulated point in the absence of any restriction. In the interest of fairness, CASA also reserves a portion of available capacity for short-haul flights and/or flights that may file a flight plan shortly before their intended departure. The slot allocation process can be described as follows:  Slot Allocation List (SAL) – First, for each regulated point (point on the route, sector or airport) the list of available slots are built based on the basic rate as a measure of capacity. For example, a basic rate of 30 flights per hour would result in the SAL with slots separated from one another by 2 minutes.  Pre-Allocation Stage – In this stage, traffic demand is evaluated based on most precise data available (Repetitive Flight Plan, Filed Flight Plan, Amendments, etc.). For each flight, provisional slot based on the order of their Estimated Time Over restricted location (ETO) is given. CASA automatically updates its solutions every few minutes, trying to improve slot allocation. When new flight data is received it is more likely that provisional slots would be changed. It’s not rare that initial change leads to chain reaction because flights whose slot has been taken tries to recover another slot by taking already reserved one, etc. Situation is reversed when a flight is canceled, which may improve slots given to other flights.  Allocation Stage – At fixed time before Estimated Off-Block Time (EOBT) of each pre-allocated flight, Slot Issue Time (SIT) is allocated to the flight and CTOT is calculated. It’s confirmed by sending Slot Allocation Massage (SAM) to 24 AO and ATC responsible of flight departure airport. Once confirmed, an allocated slot cannot be taken by another flight. If flight is subject to several regulations it is given delay of the most penalizing regulation.  Revision Process – Even after SAM is sent, CASA attempts to improve slot allocated depending on AO flight status. If the flight is in the RFI status (Request For Improvement), which is the default status, when an improvement is possible it immediately receives a Slot Revision Message (SRM). If flight status is set to SWM (SIP Wanted Message), then AO receives Slot Improvement Proposal (SIP) message which may be accepted or refused by sending appropriate response. There is no doubt that CASA is efficient tool capable of solving European level size problem within minutes. The unknown is the quality of the solution. Experimental results reported in [Vranas, 1996] shows that the total delays assigned by CASA were roughly 40% above optimum given by an exact IP algorithm. However CASA only needed half of minute to solve problems involving around 3000 flights and 25 sectors, while optimal solutions were found in approximately one hour. 1.3. Future ATM Historically, growth of air traffic has not been coupled with an adequate increase in capacity (airport and airspace). As consequence, rising presence of delays in the European ATM system has adverse effects on the economics of air transport [Jovanovic, 2011]. It is envisioned that incremental changes in technology and procedures are no longer sufficient to keep up with the growth in traffic [Neal et al., 2011]. ICAO Global Operational Concept proposes significant changes to the organization and delivery ATM services in the next 15 years [ICAO 9854]. Major systems development programs are underway around the world including the Next Generation Air Transportation System (NextGen) program in the United States [FAA, 2012], and the Single European Sky ATM Research (SESAR) program in Europe [SESAR, 2013]. One of the cornerstones of both programs is the Trajectory Based Operations (TBO) which implies moving from Airspace to 4D trajectory management. A major challenge in the design of future TFM systems is to design an adaptive system that can handle both variations in the magnitude and distribution of traffic. The system must be able: to handle two to three times current traffic, 25 new types of vehicles like very light jets and unmanned air vehicles, and to make good decisions in the presence of uncertainty. Elements of the future operational concept include changes to the organization and management of the airspace design to improve access and utilization, dynamic management of capacity to meet demand and respond to uncontrollable events (e.g., weather and emergencies), dynamic and flexible management of trajectories and traffic synchronization to improve safety and efficiency, and conflict management (Figure 7). Figure 7: General functions in future ATM [Neal et al., 2011] The primary means by which the priorities such as safety, efficiency, and cost-effectiveness are to be achieved is through trajectory management. Trajectory management and traffic synchronization, together with the other ATM components, will contribute to the efficient handling of traffic from gate to gate. There will be dynamic 4-D trajectory control and negotiated conflict-free trajectories. The ICAO proposes that airspace users should retain primary responsibility for the provision of ATM services. However, the Operational Concept also states that the allocation of responsibilities is subject to the design of the ATM system. Delegation of maintenance of spacing to the flight deck aims to increase traffic throughput while reducing ground system workload [ICAO 9854], [Dwyer and Landry, 2009]. Through trajectory management shown in Figure 8, 4D trajectory will evolve from preferred trajectory published by airspace user (Business Development Trajectory - BDT) to Reference Business Trajectory - RBT that airspace user agrees to fly and airspace service provider (ANSP, airport) agrees to facilitate. RBT constitutes an agreement between stakeholders for the whole airborne trajectory to destination [SESAR, 2013]. The trajectory management process starts with the collection of all BDTs and then is shared with ATM partners. These trajectories are then modified as necessary using a layered Collaborative Decision Making (CDM) planning process which takes account of identified constraints. The ATM planning process is one of continuous refinement 26 as better data become available. The output of CDM just before flight execution is the final trajectory called the Reference Business Trajectory (RBT) with related tolerance windows (time windows). Figure 8: The Business Trajectory Lifecycle [SESAR, 2013] TBO concept restore flexibility to an air transportation system by freeing aircraft from “the old highways in the sky” that are dependent on ground-based beacons [FAA, 2012]. TBO enables more direct, fuel efficient routes and provides alternatives for routing around disruptions or unexpected congestion. The introduction and implementation of TBO concept have to be followed by the development of new types of displays (cockpit and controller’s), decision support tools and management system (SWIM), communication, navigation and surveillance technologies, and evolution of meteorological service. Trajectory management concerns trajectory planning as an optimization problem equivalent to the traffic assignment problem with the current system. The objective of this type of planning is to determine an optimal trajectory shape (3D), departure time and speed profile that optimizes a given criterion [Delahaye and Puechmorel, 2013]. Methods used to manipulate trajectories, which are of infinite dimension, are examined in following section. 1.3.1. Trajectory modeling Trajectories are objects belonging to spaces with infinite dimensions. They are described as mappings from a time interval [ ] to a state space , having either or depending on whether speed is assumed to be part of the state. In order to manipulate such objects with algorithms, one must reduce the dimension of the search space. Survey of mathematical models for aircraft trajectory design is given in [Delahaye et al., 2013]. 27 Some dimension reduction tricks include:  Piecewise Linear Interpolation – Trajectory is characterized by way points (WP) and straight line segments between those points. The number of WP highly depends on trajectory curvature. If one wants to approximate trajectory with many shape turns, the number of way points needs to be increased in order to reduce the error between the model and the real trajectory. To improve concept performance, Lagrange and Hermite interpolation might be used to adjust polynomial function to a given set of way points. However, due to Runge’s phenomenom these interpolations induce oscillations between interpolation points (Figure 9). Figure 9: Lagrange interpolation result for a set of aligned points [Delahaye et al., 2013]  Piecewise Quadratic and Cubic Interpolation – Piecewise quadratic interpolation considers quadratic curves on the intervals [ ] ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ , each connecting two points is space. On each point, derivative of the previous curve has to be equal to the derivative of the next one, which cause the change of the interpolating curve whenever any point moves. The piecewise cubic interpolation avoids this drawback but do not insure that trajectory curvature is continuous.  Spline interpolation - also called natural spline because it represents the curve of a metal spline constrained to interpolate some given points. Spline interpolation is preferred over polynomial interpolation because the interpolation error can be made small even using low degree polynomials for the spline. Cubic spline represents a good approximation for aircraft trajectories [Delahaye et al., 2013].  Bezier Approximation Curve - A Bezier curve is defined by a set of n control points, with n represents curve order (n = 1 for linear, 2 for quadratic, etc.). The first and last control points are always the end points of the curve; however, the intermediate control points (if any) generally do not lie on the curve (Figure 10). When interpolation is not a hard constraint, one can use control points to 28 change the shape of a given trajectory without forcing trajectory to go through every single control point. Figure 10: Cubic Bezier curve [Delahaye et al., 2013]  Principal Component Analysis - When trajectories samples are available (for instance radar data), one can build a dedicated base which will minimize the number of coefficients (principal components) for trajectory reconstruction. Another approach for trajectory design is based on the wavefront propagation principle. It is a well-known fact in physics that waves of high frequencies tend to propagate along the minimum time trajectory. The principle of the wavefront propagation algorithms is to simulate such physical propagation model in order to find geodesic trajectory based on a criterion that has to be optimized (Figure 11). Some of existing methods are: Fast Marching algorithm [Sethian,1999], Ordered Upwind algorithm [Sethian and Vladimirsk, 2003], Light Propagation algorithm [Dougui et al., 2010], etc. Figure 11: Lunching rays – Light Propagation algorithm [Dougui et al., 2010] Path-planning is a term commonly used in the robotics to refer to the problem of generating an obstacle-free path to be followed by a vehicle in a two or three dimensional space containing obstacles [LaValle, 2006]. Because the vehicle 29 dynamics are not taken into account in these path-planning, resulting path cannot be followed exactly or even close by the vehicle. One way to ensure that the resulting paths correspond to feasible trajectories, satisfying the vehicle dynamics, is to use optimal control theory. The objective of optimal control theory is to determine the control input(s) that will cause a process (i.e., the response of a dynamical system) to satisfy the physical constraints, while, at the same time, minimize (or maximize) some performance criterion. First method that was possible to handle complicated differential equality constraints was proposed by Soviet mathematician L.S. Pontryagin (indirect method). Unfortunately, closed form solution to the previous problem is difficult [Delahaye et al., 2013]. In recent years, direct methods have become increasingly popular as they do not require closed form expression for the necessary conditions [Jain and Tsiotras, 2008]. The main idea behind direct methods is to discretize the states and controls of the original continuous-time optimal control problem in order to obtain a finite- dimensional nonlinear programming problem (NLP). A major issue with almost all current trajectory optimization solvers (direct or indirect) is the fact that their computational complexity is very high [Delahaye et al., 2013]. 1.3.2. Previous research on Trajectory management In the SESAR framework [SESAR, 2013] the need to increase capacity of air transportation system motivates 4D trajectory management. It concerns the generation of the trajectory from origin to destination taking into account several objectives: minimization of trajectory length and travel time while avoiding obstacles (fixed and moving). In literature such problem usually refers to finding conflict-free trajectories while minimizing given objective. Several classes of methods are used to address this problem. One class of models uses different optimization techniques for conflict resolution. One of the earliest approaches is based on genetic algorithms (GA). This optimization method is based on the evolution theory and uses basic operators: selection, mutation and crossover for generating new population of aircraft trajectories. The state space is a set of finite maneuvers in the horizontal plane (straight line, turning point and offset), as well as vertical maneuvers such as level-off. Those maneuvers are used currently by ATCo. GA generates trajectories with feasible operational maneuvers and maintains velocities within acceptable range [Durand and Alliot, 1997]. However, these techniques have not 30 been tested yet with curved trajectories [Delahaye et al. 2010]. The first approach for conflict resolution using stochastic optimization algorithms was done by Alliot et al. in [Alliot et al., 1993]. They consider horizontal airspace where every aircraft enters and leaves the airspace at the same time. All aircraft fly at the same speed and the only possible maneuvers are heading change (left or right by 30 degrees) at discretized time steps. The objective was to minimize the distance of the last point of the computed trajectory to the theoretic exit point. A second model proposed by the authors included speed changes and fuel cost in the objective function. The results of the model are compared to the classical graph search such as A*. A* needed a long time to compute the solution, while GA was able to calculate a good solution within a few seconds. In order to get closer to real ATC system Durand and Alliot in [Durand and Alliot, 1997] and [Durand et al., 1996] considered finding conflict free trajectories simultaneously by minimizing delay due to deviation imposed on aircraft, number of maneuvers and their duration. They take into account speed uncertainty and real aircraft performance model. The decision variable considers horizontal and vertical maneuvers. The model is tested on real problem over France airspace. Particle Swarm Optimization (PSO) techniques to perform two- and three- dimensional flight path optimizations using primitive maneuvers, compliant with operational constraints, is presented in [Blasi et al., 2011]. In this approach, flight paths are divided into a finite number of segments each associated with an elementary maneuver: straight flight, turn right, turn left, and alignment to the target. Sequence of maneuvers is controlled by the PSO. Prandini et al. in [Prandini et al., 1999] proposed probabilistic framework of conflict detection and resolution for a pair of aircraft flying at the same altitude thus allowing uncertainty in aircraft position to be explicitly taken into account. Another probabilistic framework (random wind disturbance) for the conflict resolution problem, where multiple aircraft in a specific region are required to reach a different target zone at the minimum expected time, while maintaining safe separation is presented in [Kantas et al., 2010]. Optimal policy, that automatically generates optimal and safe maneuvers for each aircraft, is computed using a sequential Monte Carlo approach. In [Eele and Richards, 2010] optimal trajectory for one aircraft constrained to avoid fixed obstacle is proposed. Novel rapid updating techniques for use with nonlinear branch and bound algorithm is 31 presented. The key feature of the model is the ability of finding a global optimal solution whilst retaining the full nonlinear dynamics model of the aircraft. Another class of methods uses a force field to generate a solution to a conflict. With this approach each aircraft is treated as a charged particle and electrostatic equations are used to generate resolution maneuvers. Each aircraft is attracted to its destination (positive charge) and the repulsive forces (negative charge) between aircraft or aircraft and obstacle provide maneuvers that avoid collisions. This enables automatic generation of conflict-free trajectories (Figure 12), with a mathematical proof. The major drawback of a force field method is a continuous aircraft maneuver in response to the changing force field [Kurchar and Yang, 2000]. Additionally, obtained solution does not respect ATM constraints such as speed limits or trajectory smoothness [Delahaye et al., 2010] and completely neglects objective optimization [Durand and Alliot, 1997]. However, several implementations of force field methods resolve these problems and shown that force field resolution can be effective when properly applied. A force field formulation based on the closest point of approach that follows the principle of increasing the minimum distance between mobiles (aircraft) is addressed in [Zeghal, 1998]. Control of multiple non-holonomic air vehicles using model predictive control and decentralized navigation functions was proposed in [Roussos et al., 2008] and [Roussos et al., 2009]. Figure 12: Principle of force field repulsive force [Delahaye et al., 2010] In [Delahaye et al., 2010] new methodology for trajectory planning using B- spline is presented. Aircraft conflict resolution is formulated as an optimization problem whose decision variables are the spline control points. B-splines are parameterized curves generalizing the Bezier curve concept and are an efficient approximation tool. This problem was solved using GA for which primary objective was to minimize the number of conflicts and the second one was to 32 minimize the total extra flown distance. Another trajectory dimension reduction technique using linear piecewise interpolation is presented in [Chaimatanan et al., 2012] and [Chaimatanan et al., 2013]. The problem addressed in this work is termed strategic de-confliction of aircraft trajectories where the objective is to minimize the number of conflicts by route and slot allocation. The proposed methodology for alternative trajectory design considers rerouting in the horizontal plane while keeping optimum speed and vertical flight profile. Therefore it does not only comply with aircraft performance but also gives trajectories that are close to optimal one chosen by the airline. Hybrid-metaheuristic algorithm based on Simulated Annealing (SA) has been developed to address this problem. Proposed optimization model is tested on air traffic data over the European airspace and is able to address more than 30.000 flights. Another method for generating conflict free 4D trajectory is presented in [Dougui et al., 2010]. Trajectory planning algorithm is based on the principle of least action: “The path of the light ray connecting two points is the one which time of transit, not the length, is minimum.”, and is able to compute smooth geodesic12 trajectories avoiding obstacles modeled by high-index areas. Light propagation algorithm (LPA) belongs to the wavefront propagation class algorithms, as any point on a wavefront can be considered as the source of tiny wavelets that propagate forward. Propagation is discretized in space and time and a branch- and-bound algorithm is used to solve the associated optimization problem. Two implementations with static and dynamic obstacles are considered and solved. Model extension to account for uncertainties in trajectory prediction is presented in [Dougui et al., 2012] by the same authors. In recent years, problem of trajectory planning becomes increasingly popular. Many new techniques for trajectory design have been proposed. However presented methods hardly take aircraft position uncertainty into account. In the next chapter an alternative approach, dealing with aircraft position uncertainty and unplanned situation, called robust trajectory planning is presented. 12 In geometry, a geodesic is a generalization of the notion of a “straight line” to “curved spaces”, and represent the shortest path between points in the space. 33 2. ROBUST TRAJECTORY PLANNING: PROBLEM AND METHOD DEFINITION 2.1. Problem statement The foundation of the future ATM Concept is a trajectory based operations – TBO [SESAR, 2013]. ATFM will shift into trajectory management, concept trough which 4D trajectory will evolve from preferred trajectory published by airspace users (Business Development Trajectory - BDT) to Reference Business Trajectory - RBT that airspace user agrees to fly and airspace service provided (ANSP, airport) agree to facilitate. RBT constitutes an agreement between stakeholders for the whole airborne trajectory to destination [SESAR, 2013]. The trajectory management process starts with the collection of all BDTs. These trajectories are then modified (if necessary) through CDM process taking into account identified system constraints. Outputs of planning process are RBTs and related tolerance windows (time windows). Inability to meet a tolerance window requires reassessment of the situation. When the tolerance window is not reachable or in the case of unplanned, crisis or contingency situation trajectory revision should be considered. This is the last option (back-up) as new trajectory might not be acceptable to the airspace user due to operational constraints and is generally a more expensive option. While there has been some work (chapter 1.2.2), which deal with ATFMP under capacity uncertainty, only few papers deals with aircraft position uncertainty. In strategic de-confliction model presented in [Chaimatanan et al., 2013] aircraft are separated taking into account “freedom margin” that is slightly higher than the separation norm (6 NM). At the tactical level, [Dougui, 2012] considered the longitudinal uncertainty of aircraft position and showed that the number of conflicts detected is 50% higher even using uncertainty reduction techniques such as Required Time of Arrival (RTA)13. Conflict solver, defined in [Durand and Alliot, 1997], assumed uncertainty in the aircraft future position (caused by error in ground speed and vertical speed prediction) but becomes saturated very fast. Finally, [Barnier and Allignol, 2012] take uncertainty with 13 As presented in the work, this mode enables to impose the aircraft to arrive at a given point of its trajectory at given time RTA with a tolerance of ±10 seconds. 34 respect to departure time. Takeoffs are randomly shifted by a bounded amount of time. It was shown, even small uncertainty of ±2 minutes generated a tremendous amount of ground delays as it was the only decision variable considered. It is clear that a novel approach to the problem is needed in order to cope with aircraft position uncertainty. An alternative approach, dealing with uncertainty in aircraft position (inability to cope with BT) and unplanned situation, is through building a more robust flight plan (robust flight trajectories) at a pre-tactical level. In this way, flight plan itself becomes unaffected by such disturbances reducing the need for tactical actions. The concept, proposed in this thesis, aims to find a set of robust 4D trajectories from origin to destination points minimizing total additional costs (including fuel, time...) caused by deviations from the user preferred trajectory (UPT) on pre-tactical level. Robustness includes both reducing the likelihood of disruption and increasing the number of options to recover from a disruption easily. This will affect the tactical controller workload alleviating conflict resolution task and traffic management at the tactical level as the primary target of future ATM [ICAO 9854]. Further, adding robustness to flight plan might result in higher airline planned costs, but disruption costs (cost of management actions taken to resolve conflict on tactical level due to disruptions) will be reduced due to reduced need for tactical interventions, finally leading to reduction in airline real operating costs. ICAO Global Operational Concept [ICAO 9854] envisioned that the control process for trajectory management needs to be adaptive, as it must respond to changes in anticipated demand and unforeseen events (weather, emergencies, etc.) One way to adapt is by adjusting decision criteria: safety, efficiency (from a single- flight perspective), and cost-effectiveness (of ATM system as whole). In current ATM system, the subjective importance of efficiency and cost-effectiveness to the air traffic controller will decrease as the anticipated level of workload increases [Leroux, 2000]. It is expected that a similar process will operate in the future ATM system, although the way that decision maker trade-off among decision criteria will vary. Therefore, problem termed Robust Trajectory Planning – RTP addresses how to optimally generate a set of 4D trajectories from origin to destination points in order to manage the two confronted objectives: 35  maximizing total robustness,  minimizing total additional costs caused by deviations from UPT. Generation of 4D trajectory assumes assignment of: route (3D trajectory shape), slot of departure and possibly speed profile. Origin and destination points could be ground fix (airport) or airborne fix (entering point of observed airspace). RTP aims to find system optimum (SO), a solution that is optimal from a viewpoint of a system as a whole, minimizing some aggregated indicators. On the other hand every user wants to use its optimal trajectory and thus issue of fairness, fair distribution of delay and/or route extension among airlines, flights, etc. has to be taken into account, as far as possible. Each individual user may have different optimality criteria: shortest distance, shortest time, minimum fuel burned, etc. RTP may be further constrained, by introducing ATS capacities e.g. arrival/departure airport capacity, including no-fly zones (severe weather cells, restricted or prohibited zones, etc.) that flights should avoid. Consecutive flight phenomena (daily sequences of flights operated by an aircraft) might also be taken into account, as a delay of the preceding flight usually cause a take-off delay of successor flights (delay propagation). The rest of the chapter is organized as follows. First, in the section 2.2, robustness is defined and some metrics that can be used to evaluate robustness of the solution are presented. Section 2.3 examines possible trajectory management actions and system search space. Description of the homotopy route generation model is presented in the section 2.4, followed by elaboration of the objective function in the section 2.5. Finally, key model assumptions, mathematical formulation of optimization model, and size of formulation are presented in the sections 2.6, 2.7 and 2.8 respectively. 2.2. Robustness and Flight interaction Future ATM might become more vulnerable to unforeseen disturbances because of the fact that the traffic level keeps growing and it is more difficult to control less structured networks. Building more robust flight schedule at pre- tactical level might help fight against such disturbances. To increase robustness, one must first define what “robustness” is? Which disturbances have to be taken into account? What are the indicators that quantify robustness? [Snelder et al., 2012] 36 Robustness definition is usually problem dependent [Lan, 2003]. Building robust (flight) plan that minimize some cost function (as RTP does) usually intends to find optimal solution that is not sensitive to small variation of the design variables or environmental parameters. However, ATM is a safety critical system and its main task is to provide a safe flow of air traffic before making it punctual and expeditious. Therefore, robustness is considered as "the ability of a system to resist changes without adapting its initial stable configuration" [Wieland and Wallenburg, 2012]. Increased robustness of the planning will reduce the probability of conflicts to occur and the need for tactical actions, finally leading in reduction of airline operating costs on average. Various indicators may be used to characterize the robustness in the context of trajectory planning. An alternative way of quantifying the solution robustness is to measure its vulnerability to environmental change as opposite indicator. While the robustness describes the strength of a solution, the vulnerability describes its weakness [Snelder et al., 2012]. In this thesis, flight interaction is chosen as a measure of solution vulnerability and indirect measure of solution robustness (higher flight interactions is, lower solution robustness). 2.2.1. Flight interaction Interaction between two flights may be defined as a situation where flights compete for the same point in 4D space. Interaction and conflict are not the same. The conflict is a real situation that might happen during flight evolution. Having fixed boundaries, the conflict occurs when two aircrafts are separated less than 5NM horizontally and 1000ft vertically (current standards). Conflicts are forbidden in normal operations as they jeopardize flight safety, thus many safety barriers exist to prevent them from happening [Netjasov et al., 2013]. On the other hand interactions are situations existing at the planning level and take into account the uncertainty propagation (deviations from the RBT). Maximum interaction is reached if 4D position (space + time) of the two aircraft equals on planning level, i.e. situation when both aircraft came into same point at the same time, and decrease with distance in 4D space. In the presence of uncertainty, minimizing interactions between flights at planning level will minimize conflict probability in reality and thus maximize robustness of the solution. 37 When flights’ 4D positions (including nominal and all alternative routes for a flight) are separated, there is no interaction between those flights no matter what action is taken by each individual flight (assumed that actions are bounded e.g. maximum rerouting, maximum delay, etc.). In such situations trajectory planning for each flight can be treated independently14. Two cases may happen:  Temporal separation – if the first flight (flight that departs earlier) arrives at the destination point before the second flight departs. This is the situation where aircraft operating those flights are not airborne at the same time (Figure 13). Figure 13: Flights fully separated in time  Spatial separation – if 3D position spaces of flights do not intersect. This means that aircraft operating those flights cannot occupy the same space (Figure 14). More specifically spatial separation may be viewed in horizontal plane and thus guarantee interaction free route not depending on vertical profile, departure time, etc. (Figure 15). Figure 14: Aircraft 3D position spaces Figure 15: Full 2D separation 2.2.2. Flight interaction metrics There are many uncertainty sources that can deflect an aircraft from its intended position (initial delay, wind, atmospheric temperature, actual aircraft 14 This is necessary condition to apply decomposition methods when solving trajectory planning problem. Height Latitude Longitude 38 weight, etc.). Difference between aircraft actual position and planned position may happen in space and/or time. Taking into account that advanced future avionics15 will enable full compliance with given route (as acknowledge by WP6 “Technological Enablers” of Episode 3 project in [Eurocontrol, 2009a]) time uncertainty remains the main cause of disruption (in general) and is therefore mainly considered in the thesis (Figure 16). Figure 16: Position time uncertainty Therefore in this thesis, flight interaction is defined as a function of time separation i.e. difference between times of arrival at intersecting point. Flight trajectories do not always intersect and might be close parallel. In that situation points of closest approach (CPAs) are determined representing points on flight trajectories where distance between two aircrafts reaches its minimum value [Netjasov et al., 2011, 2013]. If horizontal distance between CPAs is less than given norm, there is interaction between those two flights and flight interaction is calculated in the same manner as for intersecting trajectories. Taking into account that flight trajectories are curves in 3D space, it might be possible that two trajectories intersect more than once or have multiple CPAs. Then, interaction between two flights is computed as a sum of interaction between all intersection point and CPAs. Time interaction and time separation are inversely proportional as time interaction decrease when time separation increases. If represent time separation that gives no interaction (e.g. 10, 15 minutes) and is the exact time separation between points and of flights and respectively, then a measure of interaction at those points is defined as follows: { ( ) The function f may be constant, linear, polynomial, exponential, etc. Constant function represents interaction counter as any separation less than given threshold maxTS will be equally treated. To take into account magnitude of the time 15 Current advanced flight management system (FMS) gives accuracy equal to or better than +/- 1NM for about 95% of flight time (closed loop with feedback). Advance required navigation performance (A-RNP) concept is developed in order to enable further improvements. [Eurocontrol, 2010] 39 separation, linear function could be used. Finally non-linear function may be used in order to penalize larger flight interactions (Figure 17). Figure 17: Flight interaction functions Although space uncertainty is not considered in this work, several metrics that takes it into account and that could be used to compute flight interactions are elaborated in following paragraphs. Presented time interaction metric accounts for time uncertainty and might also implicitly account for small space uncertainty. To take aircraft position relation into account, traffic density metrics have to be used. One very simple traffic density metric is represented by the number of aircraft in airspace blocks (cells) during certain time period. By measuring number of flights that are in each cell during each time period, one can build a congestion map showing regions with high operational congestion. Although highly correlated with flight interaction this metric does not take account of the orientation of traffic and considers geometrically structured and disordered traffic in the same manner. To take density and disorder of traffic explicitly more sophisticated measures such as complexity metrics has to be used. In the literature, traffic complexity is used as a measure of controller’s workload as it represents traffic closeness in space and time. Although controller workload is not considered in the trajectory planning models, due to analogy, many factors identified affecting controller performance might be used. Based on interviews and survey techniques with ATCo performed at NASA ARC number of factors has been identified as significant contributors to the workload [Sridhar et al., 1998]. A proposed measure of complexity, dynamic density (DD), is calculated as a weighted sum of those factors. Some of the factors 40 are: number of aircraft with 3D Euclidean distance between 0-5NM, 5-10NM, and a number of aircraft with a lateral distance between 0-25NM, 25-40NM and 40- 70NM where vertical separation is less than standard norm, measured during a sample interval of time (one minute). In [Delahaye and Puechmorel, 2000] geometric approach is used to capture intrinsic traffic disorder that is highly related to the complexity. Metric is calculated as a combination of two indicators: proximity indicator that quantifies the level of aggregation in relation to a given volume, and convergence indicator that quantifies the geometric structure of the speed vectors of airplanes presented in a zone (Figure 18). Figure 18: Proximity/convergence referential for traffic complexity measure Another metric developed by the same authors based on dynamic systems is presented in [Puechmorel and Delahaye, 2009]. In this approach, air traffic is modelled by a dynamic system and complexity value in the airspace is assessed by the associated Lyapunov exponent. Lyapunov exponent is a measure of sensitivity to initial conditions of the underlying dynamical system. A small value of the Lyapunov exponent means that the future is highly predictable as the trajectory of a point of the dynamic system is not sensitive to small disturbance of the initial position (distance between the initial trajectory and the trajectory under disturbance is very small). As a result complexity at the point is low. On contrary high values of Lyapunov exponents represent zones of disordered traffic with 41 high complexity. Calculation of Lyapunov exponents is a difficult task mainly as it involves regression of the non-linear dynamic system, and its application in optimization models is very limited. 2.3. System boundaries and search space System boundaries are the physical boundaries of airspace: sector, region or whole national airspace, FAB or Europe as a whole, in which RTP problem is considered and solved within a chosen time horizon. Origin and destination points of each flight entering and exiting the system may be ground fixes (airports) or airborne fixes (entering points of observed airspace). This has an influence on search space such that flights entering the system at airborne fix cannot be assigned a slot of departure e.g. cannot be delayed. Both ground and airborne fixes are given as input and cannot be changed. To reduce interaction between flights, trajectories of involved flights might be separated in any of the four dimensions within the system boundaries. Possible management actions include change in:  slot of departure,  horizontal route,  vertical profile, and  speed. The effects of the management actions (one action is taken at the time) are illustrated on a very simple toy example with two flights with intersecting routes (Figure 19). Only cruise phase is considered, with constant speed. Flights operate at the same altitude. Departure time of the first flight is t1 and it arrives at the intersecting point at ti1, which is at the distance of d1 from departure point. Accordingly departure time, time of arrival at intersecting point and distance of the intersecting point of second flight are t2, ti2 and d2. 42 Figure 19: Aircraft encounter on intersecting routes Figure 20: Time evolution of flights Figure 20 illustrate the time evolution of flights along route. Assuming that maximum acceptable interaction is represented with minimum time separation at intersecting point tmin, the problem is to speed up the first flight or to delay the second flight in order to arrive at the intersecting point in ti1’ or ti2’ respectively. Alternative departure time Changing departure time of a flight result only in time shift at each position resulting in the so called temporal solution. This approach is based on the assumption that each flight will follow a single initially chosen optimal 3D trajectory and therefore not induce more fuel consumption while airborne. The range of time shift depends on the planning level (strategic, pre- tactical, and tactical) at which the action is conducted. At the strategic level, when desired trajectories are designed by airspace users, there is no precise plan what flights will be carried by what plane (Daily Operational Plans) and therefore there is more freedom in changing slots. Hence, the departure time of each flight can be shifted by a positive or a negative time shift. Although technically feasible, large shifts are not desirable as departure times is strongly influenced by traffic demand and these changes might reduce attractiveness of a flight and affect such demand. At pre-tactical level, room for maneuver is reduced as delay of one flight leads to delay of all subsequent flights performed by the same aircraft. Furthermore, for operational reasons (passengers that are scheduled for a flight) positive shifts are only allowable and known in the literature as ground holds (tactical level). Solution of the toy example with departure slot change is illustrated in , which depict situation where second flight is delayed for tmin - t time units (green line). Delay of departure cannot be greater than the given T denoted by the t1 t2 d2 d1 tmin ti1’ t1 d2 d1 t2 t tmin Time Distance along route ti2 t i 2’ t i 1 43 red line in the figure (constraint). As slot of departure cannot be moved earlier in time17 changing slot of the first flight would not solve the problem. Figure 21: Flight delay Alternative horizontal route With the increase of traffic volume, large amount of ground holds might be required to overcome the flight interaction problem. In such situation, flights have to be separated in spatial dimension, assigning alternative route to some flights. Although longer than the original one, this new route might be more acceptable as a result of less ground holds assigned to the flights and consequently less delay at destination. Deviation of alternative route compared to nominal one is limited by delay it produces at destination and additional cost due to more fuel burned. By changing flight route, one could try to manage intersection point, such that the time difference between arrivals of first and second flight at the point become greater than tmin, or even remove it (usually not possible or too expensive). Figure 22 shows solution of the toy example where second flight is assigned to the new route. Although the 4D trajectory of the first flight remains unchanged, arrival time at the intersecting point of both flights (ti1, ti2) is changing as a result of intersecting point relocation. Route change of the second flight has impact on its route length that is increased denoted by red extension line at the Figure 22. Longer route also results in late arrival at destination point of the second flight (without change of the speed). Another solution to the problem might be to change the route of the first flight in the same manner. 17 as explained, at pre-tactical and tactical planning level flight may only be delay for operational reasons (passenger may lose their flights) t1 d2 d1 t2 t tmin Time Distance along route ti2 t i 1 t i 2’ T 44 Figure 22: Flight rerouting Alternative vertical profile Cruise level, climb and descent profile are designed by flight engineers in order to best utilize aircraft performance. Although any deviation from the optimal vertical generates additional costs, in operations, controllers often require a change of flight level to resolve potential conflicts (mainly in TMA). Change of flight level results in more fuel burns than pre- planned profile at same cruise level due to deviation from optimum profile. Therefore several alternative vertical profiles might be planned for each flight consisting of a “primary” profile preferred by the carrier and “alternative” profiles. An altitude change of one flight in the toy example is the simplest solution as there are only two flights in the system. However this might be a more expensive solution than the other alternatives. The applicability of this solution depends on other traffic that is omitted in this example. Speed management Similarly to vertical profile, nominal speed profile for a flight is designed based on flight level and weight. Although increase of cruising speed will result in more fuel burn, often in operation it is used by airline to compensate flight delay. Opposite speed reduction may lead to flight delay but additional cost could be partially compensated by fuel savings. Small variance of the nominal cruising speed (-6%, +3%) are considered acceptable as impose no significant increase in perceived controller workload [Averty et al., 2007]. Currently speed management is mainly used as a tactical measure of resolving potential conflicts (ERASMUS project) [Weber et al., 2007], but it is suggested that it could be used as tool to manage flight delay costs [Cook et al., 2007]. Figure 23 shows a possible outcome of aircraft speed change applied to the toy example. Speed change limits are denoted by red lines in the figure (operating speed range). To solve the problem each flight could be either slowed down or t1 t2 d2 d1 d2’ d1’ t1 d2 d1 t2 t tmin Time Distance along route ti2 t i 1 t i 2’ d1’ d2’ ti1’ 45 speeded up. For instance, speed decrease of the second flight will result in late arrival at the intersecting point (ti2’) although its departure time is not changed (Figure 23a). Another possibility is to increase speed of the second flight so it arrives at the intersecting point tmin time unit before first aircraft. However required speed is out of operational range. Figure 23b shows another acceptable solution of the problem where speed of the first flight is increased in order to arrive earlier than planned at the intersecting point (ti1’). a) b) Figure 23: Speed management 2.4. Route generation model As it was shown in the review of the existing models for the ATFMP (chapter 1.2.2), current rerouting is viewed as allocation of a route from a pre- defined set of alternative routes. Given that each flight in the current ATM follows airways and beacons (routes usually given as a sequence of sectors crossed by flight), a small number of alternative routes (two-three) is completely reasonable and does not affect optimal solution (provided those alternatives are different enough from each other). Last statement does not hold for the future ATM that will enable more direct, fuel efficient routes using satellite positioning and advance navigation technologies. It will free airspace from the “old highways in the sky” [FAA, 2012]. As a matter of fact, each direct route corresponds to a link in a fully connected graph. If trajectory optimization takes into account other objectives (obstacles avoidance, congestion reduction) or considers wind effect then infinitive number of alternatives exists with curve shape trajectories. Consequently, the set of possible curve routes is infinitely large. Curves are objects belonging to spaces with infinite dimensions, and in order to manipulate such objects it is necessary to reduce the dimension of the search space [Delahaye et al., t1 d2 d1 t2 t tmin Time Distance along route ti2 t i 1 t i 2’ t1 d2 d1 t2 t tmin Time Distance along route ti2 t i 1 t i 1’ 46 2013]. In this thesis homotopy is used as dimension reduction technique. It was observed that homotopic feature of continuous functions may be used to easily map large continuous spaces using only small number of parameters. This feature is exploited in the route generation model to design the shape of the horizontal route (alternative horizontal route) that is matched with vertical profile to produce 3D route. Finally, 3D shapes are completed with time dimension in order to create trajectories. 2.4.1. Homotopy In topology, two continuous functions are called homotopic (Greek ὁμός (homós) = same, similar, and τόπος (tópos) = place) if one can be "continuously deformed" into the other. Such a deformation is called a homotopy between the two functions. Formally, a homotopy between two continuous functions f and g from a topological space X to a topological space Y is defined to be a continuous function   YXH  1,0: from the product of the space X with the unit interval [0,1] to Y such that, if Xx then    xfxH 0, and    xgxH 1, [Katok and Sossinsky, 2006]. Reference curves (two solid lines in Figure 24) given by functions  t0 and  t1 are homotopic relative to their endpoints because there is homotopy  ,tH that continuously deform  t0 into  t1 .        tttH 10 1,   Dashed lines in the Figure 24 represent iso-contours of this homotopy computed as a weighted sum (convex combination) of the reference functions for different values of the parameter  . Homotopy  ,tH maps space between reference functions using only one real number  ([ ] ). Figure 24: Example of two homotopic functions  47 In route generation model, homotopy is used as tool for horizontal route generation. The model performance heavily depends on the choice of reference functions as iso-contours of homotopy replicate their pattern. Further, for many nonholonomic18 systems, such aircraft [Roussos et al., 2008], the route geometry has critical influence on the feasibility and performance of route tracking. The following list represents general criteria that must be satisfied for a function to be a candidate for reference functions:  The reference function domain is enclosed by origin and destination point given on the axes passing through these two points (direct route). Function codomain (range) has a width of , where is length of the direct route. These two criteria ensure that an alternative route will not deviate too much from nominal direct route restricting route search space (Figure 25). It has to be noted that width of search space from the figure represent two times because two (symmetric) reference function are used for every homotopy. Figure 25: Route search space  Reference functions should be symmetric with a direct route as the axis of symmetry. This enables the direct route to be a result of homotopy for certain values of parameter ( ). This is required because the direct route is taken as nominal (initial) route;  Reference function should be a single-valued function, which means that, the image of an element of the domain is always a single element of the codomain. This restriction ensures that the aircraft is moving always towards the destination point, never away from it (Figure 26); 18 A system that is described by a set of parameters subject to differential constraints. Its state evolves along a path in its parameter space that is continuous. In robotics system is nonholonomic or anholonomic if the controllable degrees of freedom are less than the total degrees of freedom. [LaValle, 2006] D D/4 48 Figure 26: Vertical line test for a function  Reference function should belong to higher order differentiability class i.e. has a high order of continuity. Hence, function must be continuous and should have derivatives of k-first orders that are continuous too. A discontinuity in the derivative of the route may correspond to a sudden change of the speed vector, which would get the route infeasible. Besides, for two routes with the same length, and the same initial and final positions, better tracking performance can usually be achieved with the smoother route [Yiming, 2012]. 2.4.2. Mathematical formulation This sections examines computation of horizontal route as symmetric homotopy with respect to reference (primitive) function ( ), and as multiple homotopy of symmetric homotopies . Figure 27 shows one pair of symmetric reference functions ( ) and ( ), and symmetric homotopy with respect to those reference functions. Any points of function ( ) may be represented as vectors in Cartesian coordinate system, and for given value of , two point on reference functions are specified by ( ) and ( ), where ( ). Figure 27: Symmetric homotopy with respect to reference function ( ) 49 As previously defined, homotopy is a convex combination of reference functions and hence point at iso-contour for certain values of may be represented as a convex combination of points and . ( ) ( ) ( ) ( ) ( ) ( ( ) ( )) For any given symmetric homotopy with respect to reference functions ( ) ( ) etc. point are calculated in the same manner and expressed as function of and . ( ) ( ( ) ( )) ( ) ( ( ) ( )) Multiple homotopy of symmetric homotopies ( ), ( ), etc. with respect to reference functions ( ), ( ), etc. is calculated as a weighted arithmetic mean i.e. a weighted average of all points given by corresponding symmetric homotopies. The value 2 1 of the parameter  is taken as ‘identity element’ because the resulting iso-contour for any symmetric homotopy is direct (nominal) route. Weights are thus calculated as deviation of parameter  from 2 1 in absolute value 2 1 . Finally, multiple homotopy ( ) is computed as: ( ) | ⁄ | ( ) | ⁄ | ( ) | ⁄ | | ⁄ | ( ) ( ) ( | ⁄ | ( ) | ⁄ | | ⁄ | ( ) | ⁄ | ( ) | ⁄ | | ⁄ | ( )) ( ) | ⁄ | ( ⁄ ) ∑ | ⁄ | ( ) ( ∑ ( )) 50 2.4.3. Alternative route design Respecting given vertical profile, the route generation model is only focused on modifying the shape of a horizontal route. Taking into account that areas around airports – TMA will remain controlled19, horizontal route shape is only modified in the en-route segment (Figure 28). Figure 28: Alternative horizontal route To unify shape design process for all possible cases (all flights and all origin-destination pairs) and simplify route manipulation (length calculation, extension, etc.), new horizontal route shape is designed in horizontal rectangle [ ] [ ] (Figure 29) using previously defined model. Coordinates of the origin and destination point are ( ) and ( ) in this coordinate system respectively. Final route shape (horizontal route) is decoded by a scaling, rotation and translation for given real coordinates of origin and destination point. 19 Airports are physically constrained by the number of runways and their configurations, and it is assumed that traffic in the vicinity of airports must remain structured and to be operationally managed – sequenced. Today more mentioned eco-policies for noise and CO2 emission reduction further limit operations in terminal. Even SESAR and NextGEN foreseen that terminal area will remain controlled based on those reasons. TMA dependencies on ATM system could be modeled through its capacity. en-route segment TMA TMA 51 Figure 29: Horizontal route shape design and decoding Length of the alternative horizontal route is always larger than the length of nominal direct route. Therefore, direct mapping between given vertical profile and new horizontal route is not possible. One solution is to re-simulate flight using aircraft performance model, which is time consuming, or to extend the vertical profile to match the new route length. [Chaimatanan et al, 2013] has shown that a more efficient approach is to extend the vertical profile at the top of descent in the cruise phase of flight. The result is an acceptable approximation of the profile that respects optimal climb and descent profile as well as speed profile (Figure 30 bottom). Based on the original trajectory samples, the new shape given by homotopy route generation model is converted into trajectory by putting the former samples (produced by simulation) on this new shape (Figure 30). Procedure of alternative route design can be summarized in following steps:  step 1 – For given set of symmetric homotopies respect to reference functions ( ) and given set of parameters controlling homotopies , find new en-route horizontal route shape in rectangle [ ] [ ].  step 2 – Based on start and end point of en-route segment scale, rotate and translate given shape to decode it to real coordinate system.  step 3 – Calculate en-route extension and extend given speed and vertical profile by prolonging cruise phase of flight. This is done by multiplying (adding) flight segments at the top of descent as many times as necessary to make the length of such profile equal to the length of new horizontal shape.  step 4 – Match trajectory samples of the newly generated profiles (speed and vertical) with new horizontal shape (horizontal coordinates). This result in set of 4D points of new trajectory. 52  step 5 – for a given delay (if any), increase time dimension of each 4D points of new trajectory. Figure 30: Matching new route shape and given vertical profile 2.5. Objective function Opposite to the conventional models that minimize planned costs, objective of RTP model is the minimization of the total operating costs of the ATM users which is the sum of the planned costs and the cost of disruptions. Adding robustness to the schedule will certainly increase the planned costs but at the same time it will reduce costs of disruptions and finally the operating costs. Therefore, the aim of RTP model is to minimize total ATM user costs and total flight interactions as two confronted objectives. Minimum cost is obtained when each flight uses its nominal trajectory (nominal flight route, departure time, vertical and speed profile). Due to interaction problem, the nominal trajectory of some flights might not be available and alternative trajectories have to be used. Every alternative trajectory incurs additional costs. Choice of alternative trajectory represents a trade-off between efficiency (additional costs) and robustness (interaction reduction) of solution. Furthermore, control actions taken to solve interaction problem at one place may raise a new problem at other locations involving other flights (interdependent decisions). Consequently, it may be altitude distance along route climb cruise descent Prolonged cruise phase distance along route Route offset d1 d1 d2 d2 53 impossible to find single, ideal solution that is both optimal and robust and one of Pareto optimal solutions has to be chosen. RTP is represented as a multi-objective optimization (MOO) problem involving more than one objective function to be optimized simultaneously. The most common and straightforward approach to MOO is to sum different objectives into a single scalar value. This is so called Weighted Aggregation method. ( ) ∑ ( ) where are objectives and are weight coefficients. When using this method, only one Pareto optimal solution can be obtained. Although easy to use, it requires a priori knowledge to determine weights that is sometimes hard and tricky to get. Instead of using general weights, as a measure of the importance of each objective, it might be possible to express objectives in the same measurement unit. Since total costs as the first objective represent monetary value expressing total flight interaction, as second objective, into a monetary unit would solve the problem. The question to answer is: What is the financial benefit of flight interaction reduction i.e. what the additional cost is acceptable to achieve a certain level of robustness? If Ic stand for financial benefit to reduce total flight interaction by one unit, then aggregated objective function is defined as follows: CIcF I  wher I represents the total flight interactions and C the total flight costs. However, the value of robustness is hard to quantify. A further question is how to aggregate individual flight interaction values into a single measure. This is closely related to the definition of the robustness. If the general definition of robustness is considered: “the ability of a system to resist change without adapting its initial stable configuration”, then individual (flight in this case) with smallest robustness (poorest individual) is representative of the system. In the context of trajectory planning, however, two solutions that have equal values of the poorest individual are not identical as robustness of other individuals and quantity of poorest individuals matters too. Therefore, total interaction is calculated as the sum of all individual interactions where larger interactions are penalized using non-linear function. 54 In the academic test example (chapter 4) another approach is also examined where MOO problem is reduced to the SOO problem by choosing the single most important objective function. All other objective functions are used to form additional constraints setting lower and upper bounds for the objectives. This method is known as a Bounded objective function method. The idea is to set an upper limit for flight interaction while trying to minimize total cost. As a result, minimal cost required for certain level of robustness may be obtained. Choosing different levels of robustness will give a range of problem solutions as alternative scenarios for decision making process. 2.5.1. Cost of actions Management actions that alleviate interactions always incur additional costs. The cost may represent a number of modified flights, delay time, additional flight miles, etc. depending on possible actions (search space). When actions of different types are taken together, the more general cost metric must be considered. The most general metric (unit) is money. Therefore the cost function consists of:  delay cost – It is calculated as cost of flight delay at the destination point (time difference between new scheduled and nominal time of arrival). It consists of ground delay and airborne delay. Ground delay is generated at the origin airport by assignment of new slot of departure. Airborne delay is accumulated during flight evolution due to rerouting and speed decrease. On the other hand, the speed increase recovers already generated delay. Theoretically, it is possible that new scheduled time is before the nominal time of arrival, having then a negative delay. Negative delay is not taken into account i.e. it’ll not decrease total cost, or it could be even forbidden with the existence of the slot allocation policy at destination airport. Opposite to ground and airborne delay defined by EUROCONTROL in [Eurocontrol, 2012b], that have different costs20 as former includes fuel and other operating costs, in the context of the thesis both delays include only time cost and are hence treated jointly.  fuel cost – It is calculated as the cost of additional fuel burn due to longer route taken and/or speed increase/decrease and/or change of vertical profile 20 One minute of airborne delay cost three times more than one minute on the ground 55 (cruise altitude). If the direct route (as shortest one) is taken initially then a change of horizontal route will induce additional fuel burn due to longer route. It should be noted that shortest route is not always the one with the lowest operating costs considering other influencing factors: wind, route charges, etc. Since it is not possible to model all factors, direct routes are considered preferable by the airlines in this thesis. Speed change will always result in additional fuel burn if nominal speed is VMR21 (Figure 31). However, in real operation, positive CI22 is usually chosen denoted with speed V0 on Figure 31. Speed V0 that is higher than VMR gives lower specific range SR0. If V0 is nominal speed, speed increase will have the same effect as in ‘VMR case’ i.e. more fuel will be needed for the same distance. On the other hand, speed decrease will result in lower fuel burn but will increase delay cost. Figure 31: Specific Range relative to cruise speed [Delgado and Prats, 2009] Figure 32 shows SR evolution with cruising speed for different altitudes. One can note that SRmax and corresponding VMR increase with altitude. Maximum SRmax gives the optimal cruising altitude, and further altitude increase (until the ceiling) will result in lower SRmax. If optimal cruising altitude is initially chosen, then any altitude change will result in more fuel burned. 21 speed that maximize specific range SRmax. It is the optimal speed from consumption per NM point of view. 22 Cost Index - express the ratio between the cost of the fuel and the cost of the time 56 Figure 32: Specific Range in function of cruise speed and altitude [Delgado and Prats, 2009] 2.6. Key model assumptions Based upon previous discussion, this chapter lists the key model assumptions that have been used in the model.  Fixed demand. One of the first assumptions is that the number of flights which operate in the airspace at a given time period is fixed, and cancellation decision is not modeled. As flight cancellation is usually driven by excessive delay or rerouting, both are limited to ensure an acceptable solution.  Heterogeneous demand. The proposed model belongs to the multiclass- user traffic assignment model [Jovanovic, 2011]. It takes into account the heterogeneity of demand in terms of different aircraft types and consequently different users will have different cost functions.  Direct route is considered nominal. Beside route length, other factors such as: wind, ANS charges, etc. influence choice of the flight route in operations. However these factors are not considered in the model and shortest (direct) route is assumed to be the one with the lowest operating cost.  Finite set of alternative vertical profiles. This means that each flight can be assigned to different flight level if the nominal one is unattainable for any reason (congestion, etc.). Each alternative vertical profile is attached with initial cost due to additional fuel burn and/or late arrival at destination.  Nominal speed profile is preserved. Speed change has big influence on operating costs and is not considered in the model. Speed management is 57 mainly used as a tactical measure for solving conflicts, and by airline to compensate flight delay.  Airline operating costs are taken into account. The model only takes into account delay and fuel costs. Passenger compensation costs and other costs are not considered.  Cost functions are known. Fuel burn rate, fuel and delay cost are attached to each flight based on aircraft type that operate a given flight. Cost of flight interaction is also known and chosen by the decision maker.  No secondary flights. Departures can take off from origin airports as planned and no consideration is given to the delay propagation from earlier flights. This simplification is typically made in models at strategic level (until day of operation the set of flights for an aircraft is not known). This can be further justified by the fact that it is no longer obvious which aircraft will fly a subsequent flight in the hub network since many airplanes are capable of performing any one of multiple consecutive flights [Bertsimas and Stock Paterson, 1998].  En-route part of the flight is considered. As deregulation will only consider en-route airspace as perceived by [SESAR, 2013] takeoff and landing part of trajectories are truncated around airports within a given radius as the traffic is considered handled with specific procedures by the TMA control service in these zones [Barnier and Allignol, 2012] This is in accordance with the EUROCONTROL definition of en-route part of flight for flight efficiency calculation [PRC, 2013]. 2.7. Mathematical formulation Nomenclature. the length of considered time horizon; the set of arrival and departure points (airports, waypoints); total number of flights i.e. cardinality of a set ; the set of flights ; nominal departure time of flight , ; nominal arrival time of flight , ; 58 flight interaction between flights and , , ; flight interaction between points and , , , , ; steepness of the interaction exponential function; exact time separation between points and , , , , ; minimal time separation between points without interaction; flight interaction unit cost; delay unit cost of flight , ; fuel price per kilo; the number of symmetrical homotopies determining the shape of horizontal flight route; reference function of symmetric homotopy , { }, where [ ] [ ]; decision variable controlling homotopy for flight , { }, ; the set of indices of alternative vertical profiles for flight , ; jth alternative profile for flight , , ; initial costs of alternative vertical profile jth for flight , , ; [kg/min] fuel burn rate for flight using jth vertical profile, , ; 0-1 decision variable; equals to 1 if flight chooses alternative vertical profile j, and 0 otherwise, , ; decision variable representing (ground) delay of flight , ; maximum ground delay allowed;  set of possible ground delays; actual departure time of flight , ; actual arrival time of flight , ; 59 Decision variables. The decision variable for flight can be represented by a triple ( ). To separate trajectories in temporal space, a departure delay is associated to each flight . Actual departure time of flight is calculated as a sum of nominal departure time and assigned delay: . The nominal departure time of flight is respected when . Spatial separation of trajectories is done by association of ( ) for each flight in order to control its 3D trajectories. Vertical profile is selected from the set of alternative vertical profiles for flight . Vector ( ) controls the shape of the horizontal route, and consist of M parameter that controls multiple homotopy . Based on vector and set of symmetric homotopies respect to functions , point of the horizontal route shape is given by ( ) where: ( ) ∑ | ⁄ | ( ⁄ ) ( ) ∑ | ⁄ | Given horizontal route shape in rectangle [ ] [ ] is then decoded by scaling, rotating and translating based on coordinates of origin and destination point. Alternative profile is extended based on calculated extension of new horizontal route, and trajectory samples of newly generated profile are then matched with new horizontal route resulting in set of 4D point of new trajectory. Finally, time dimension of each 4D point is increased based on flight delay . Actual arrival time of flight is computed in route generation model. Based on computed value of actual arrival time, value of nominal arrival time and departure delay, airborne delay of a flight is calculated as: [( ) ]. Flight interaction between two flights is computed as a sum of interaction between all “conflicting” points (intersection point and CPAs), and , of those flights. Interaction between conflicting points is calculated as a function of time separation between them ( ), and exponentially increases with decrease of time separation. Given , that is input parameter, and , that is obtained by simulation as it depends on decision variables, interaction is expressed by following formula: { 60 Constraints. In practice, flow control deals with a time interval divided into a finite number of periods rather than with a continuous time variable. Therefore departure delay can be treated as discrete by dividing considered maximum allowed delay into periods of equal length  (  ). Then, the set of possible ground delays is:  {  ( - )   }. An important feature of the homotopy route generation model is that the shape and length of an alternative horizontal route are bounded by the reference functions. The only restriction is that the values of parameters controlling homotopy are between 0 and 1. Finally, the choice of vertical profile is limited to a pre-defined finite set of possible alternative profiles for each flight. Set also contains an index of the nominal vertical profile for which initial costs equals to zero ( ). Objective function. The Robust Trajectory Planning problem is formulated as an assignment of a vertical profile , a vector of parameters that controls homotopy ( ) and a delay for each flight , such that the objective function consisting of the total additional costs to network user and the total cost of flight interaction is minimized. The total flight interaction is calculated as the sum of flight interactions between all pairs of flights. Total additional costs to network users, due to deviation from UPTs, are calculated as sum of:  cost of alternative vertical profile – due to additional fuel burn and possible late arrival at destination. It is an input parameter for every given alternative profile, and equals zero when a nominal vertical profile is used,  cost of delay – due to late arrival at destination point, and  cost of fuel – due to a longer alternative horizontal route. It is calculated based on additional time that each flight spends in the air (the airborne delay). Then, the associated mathematical model is given as follows: ( ) F ∑ ∑ A F ∑ ( ) F ∑[( ) ] ∑ A F ∑ ∑ F F (1) subject to constraints: ∑ { } (2) 61  (3) { } (4) The first sum in the objective function (1) represents a total initial cost due to alternative vertical profiles. The second sum is the total delay cost due to late arrivals at destination. The third sum is the total cost of additional fuel burned due to a longer routes. Finally last sum represent total cost of flight interactions. The constraint (2) ensures that each flight can be only assigned to one vertical profile from its pre-defined set of profiles. The constraint (3) ensures that every flight delay takes a value from a given set of possible ground delays. Finally, the constraint (4) ensures that parameter controlling homotopy takes value between 0 and 1. 2.8. Complexity and Size of formulation It is shown in the literature of related problems (chapters 1.2.2 and 1.3.2) that RTP problem is NP_HARD. However the problem complexity will be briefly studied in the following paragraphs. The same nomenclature is considered. Let , continious decision variable controlling homotopy i of flight f, be discretized and let denote the set of possible values. The number of alternative horizontal routes R for flight f is then: R ∏| | M and are previously defined as a set of possible slots of departure and a set of alternative route for flight f respectively, then number of possible trajectories for flight f is: R | | | | where | |denote cardinality of the set. Then it is possible to express the number of decision variable combinations for all flights: | | ∏ N 62 If the number of trajectories is the same for every flight , then the cardinality of the state domain becomes: | | Moreover, those decision variables are not independent due to the connexion induced by flight interaction23, so decomposition methods cannot be applied [Delahaye and Odoni, 1997]. This formulation induces high combinatory in a huge space, and stochastic optimization has been used to solve it. Next chapter gives insight of optimization technique used in this thesis. 23 Although flight interaction is not constrained, since it is constituted in the objective function flights cannot be treated separately. 63 3. OPTIMIZATION TECHNIQUE The problem defined in the chapter 2 is known to be NP_HARD problem with non-separable state variables as shown in chapter 2.8. In addition, the value of the objective function is obtained by simulation. To solve real instances of problem involving high combinatory in huge state space, one must use stochastic optimization techniques. Genetic and other evolutionary algorithms (EA), ant colony optimization, simulated annealing etc. are known to be able to provide good solutions for this type of problem. Having in mind that RTP is basically multi-objective problem, using evolutionary algorithms to solve such problem would be rational choice. However one point in the state space needs a lot of memory. This is due to representation of flight trajectory as a sequence of route points each containing time, position and velocities vectors. On average, each flight contains 500 points (2hr flight sampled every 15 seconds). If one consider full day European traffic with 30.000 flights it is easy to calculate that each point in state space require about 1Gb of memory. This limit possible population size to the point that EA are no longer efficient. Then, Simulated Annealing as local search metaheuristic approach (new solutions are generated by modifying the only existing solution with local move) has been selected to address this problem based on the fact it needs less memory (only current solution is kept in the memory). 3.1. Simulated annealing Simulated annealing is a probabilistic method proposed by Kirkpatrick, Gelett and Vecchi in [Kirkpatrick et al., 1983] for finding the global minimum of a cost function that may possess several local minima. It is inspired by an annealing process in metallurgy. Solid is brought to a high energy state then it is slowly cooled-down so that eventually its structure is “frozen” when a minimum energy configuration is reached. In the optimization algorithm the objective function is equivalent to the energy of the physical system while control parameter T plays the role of the temperature. To escape local minima trap, optimization method is allowed to move to a worse solution with a given probability that decreases with parameter T. 64 The algorithm starts at random state ( ) for which neighbor states ( ) are generated. With probability 24 it select a neighbor ( ). Once is chosen, the next state ( ) is determined as follows:  If ( ) ( ) then ( ) .  If ( ) ( ) then: ( ) with probability ( ( ) ( ) ( ))25 ( ) otherwise. where is objective function or energy of given state. At each step, the temperature parameter T is decreased according to a pre- defined schedule and such process is repeated until termination criteria are satisfied (when minimum temperature or global minimum if known are reached, etc.). 3.2. Application to our problem The simulated annealing algorithm is very sensitive to the choice of parameters and most of them are problem specific. There are, however, some general rules that have to be followed when designing experiments. Thus, initial temperature has to be high enough for the final solution to be independent of initial one. In the optimization model for RTP it is determined through a heating procedure where an arbitrary number of states are randomly generated and initial temperature is then selected so it allows moves between generated states with desired probability. Cooling is done with a geometric scheme for which ( ) ( ). The cooling rate is constant in the experiment and set to 0.99 (low cooling rate) allowing better exploration of the search space. The number of iterations at each temperature level (i.e. temperature length) is controlled by temperature and it is increased with drop of temperature (allowing intensification of the search as system comes to slid state). As the objective value of global optimum is not known, the minimum temperature level is chosen as stopping criterion and set to four-hundredth part of initial temperature (empirical value). 24 In general algorithm each neighbor is chosen with equal probability (uniform distribution). However this probability might be controlled by heuristic method that is problem specific enabling better solution convergence and shorter computational time. 25 Acceptance probability function depends on energies of two states and on temperature. 65 Two different acceptance probability functions have been considered and tested. The first one is governed by the Metropolis probability [Metropolis et al., 1953]: ( ) ( ) ( ) , where probability to accept worse solutions decrease with temperature. The second one is similar to Metropolis probability (the only one that have mathematical proof of convergence) but incorporate additional parameters that are controllable and specific to the problem: . Base probability takes temperature influence into account and decreases with temperature. The multiplier is a function of energy difference and decreases when the energy difference increases. The idea behind is to increase acceptance probability of higher energy states at the high temperature to be able to effectively search rugged terrains of the objective function. Figure 33 presents an acceptance probability as a function of energy difference for three different temperature levels: high, medium and low. Figure 33: Acceptance probability function Set of neighborhood states and new state from this set are generated using heuristic approach. Heuristic favors neighbor with a higher interaction level when selecting a next state. At each iteration neighborhood is selected from those governed by completely random rule (pure simulated annealing) or local search heuristic, based on probability that changes with temperature (earlier stage – more random and vice versa). 66 3.3. Flight interaction calculation Although the shape of flight trajectory represents a curve in 3D space, it is not possible to use standard curve intersection methods26 to calculate the interaction between two trajectories. First, because time dimension has to be taken into account (however there is a workaround for this issue), and second, two trajectories can interact even if they do not intersect. Classical approach implies that trajectory is expressed as an order list of samples (position vectors) in 4D space (Figure 34). Then, a pairwise comparison of aircraft position could be used. However, this rough approach is time consuming and not suitable for a problem that involves a large number of trajectory samples. Calculation of the total interaction requires checking of all pairs of trajectories ( ) and therefore exponentially increasing with the number of trajectories ( ). Figure 34: Classical trajectory representation [Delahaye et al., 2013] In order to use position vector comparison on the large scale, a grid-based scheme is used [Chaimatanan et al., 2013]. First, the airspace is discretized using a four dimension grid as illustrated in Figure 35. Each cell in the grid has a unique location (address). Every position vector is associated with the address of the cell in the grid it belongs to. Now it is possible to calculate the interaction of a single trajectory position with the rest of the world only by checking grid cell that it belongs to and surrounding cells. If they are occupied by different trajectories, precise comparisons between pairs of aircraft position need to be performed. This is a very efficient method to compute the total interaction as its computational time depends linearly on the number of trajectories. 26 Methods that address the problem of computing the points at which two curves intersect. Predominant approaches are [Sederberg, 2009]: Bézier subdivision, interval subdivision, Bézier clipping. 67 Figure 35: Four dimension space-time grid [Chaimatanan et al., 2013] Using described grid-based scheme it is possible to identify all flight points that are separated less than given norms in space and time (e.g. 5 NM horizontally, 1000ft vertically and 10 minutes temporal) called “conflicting” points. Based on identified “conflicting” points, the total flight interactions could be measured using different metrics. 68 4. RESULTS AND DISCUSSION Method and model developed in previous chapters are first tested on a small hypothetical example, illustrating the basic idea of robust trajectory planning providing some valuable initial insights for the algorithm. The model will be then applied to a large scale real-life example. 4.1. General model settings and input data In addition to traffic data (flight plans) and specific flight data (operating an aircraft, priorities, etc.) given by the airline that differ from examples and scenarios, there are other general data necessary to make run the model. General model settings such as: reference function of symmetric homotopies that are used in these examples and methodology for flight interaction calculation are presented here. 4.1.1. Fuel Burn Rate Fuel Burn Rate (FBR) represents aircraft fuel consumption [kg] per time unit [min] in the cruise phase of flight. Although it depends on many factors such as engine type and model, number of engines, actual aircraft weight, speed and flight level, meteorological conditions, etc. it is not relevant to include them all in the model as this would increase its complexity, furthermore most of them are unknown at strategic level. In order to simplify the use of FBR, the following assumptions are made:  meteorology is completely disregarded and ISA atmospheric conditions are used,  wind effects are not taken into account,  aircraft of the same type are considered unique concerning power plants,  nominal aircraft weight27 and nominal speed28 are used. 27 BADA defined reference mass is roughly 70% of the way between the minimum (OEW) and maximum (MTOW) mass while corresponding to one of the mass values that are available in the reference trajectory data. [Eurocontrol, 2009b] 28 the most preferable speed for given flight level 69 Finally, FBR is given by a matrix which depends on aircraft type and flight level (Figure 36). FBR data were extracted from EUROCONTROL’s Advanced Emission Model – AEM that are based on BADA (Base of Aircraft Data) datasets. Figure 36: Aircraft Fuel Burn Rates 4.1.2. Flight cost category Cost category is assigned to each flight based on aircraft type. There are four aircraft categories considered: Low, Medium, High and Jumbo, with respect to operating costs. All operating costs: fleet, crew, maintenance, except fuel costs are considered. Although the two sometimes match, the above categorization is not directly related to ICAO Wake-turbulence categorization. The resulting cost category is later used to determine flight delay cost using unit cost29 for particular categories. Aircraft categorization is based on a study done by the Department of Transport Studies, University of Westminster in 2004 and updated in 2010 [Westminster, 2011] for PRU30, EUROCONTROL. The study considers cost estimate for twelve core aircraft: B733, B734, B735, B738, B752, A319, A320, A321, AT43, AT72, B744 and B763. These aircrafts were grouped into four cost categories. As RTP model contains wider aircraft database, aircraft not included in the study were given category based on similarity in weight, engines and number of passengers with identified (categorized) aircraft. 4.1.3. Homotopy functions In chapter 2.4 basic criteria for reference functions was defined and route generation model and its mathematical formulation were presented. It has been studied that trigonometric and power functions best fit criteria for reference functions: both are continuous and single-valued functions defined for all real 29 specific to each scenario 30 Performance Review Unit FL AC_TYPE 30 40 60 80 100 120 140 160 180 200 220 240 260 280 290 310 330 350 370 390 410 430 450 A320 29.8 29.9 33.9 34 34.2 34.4 47.1 47.3 47.5 47.7 47.8 48 48.2 47.9 46.5 43.7 41.4 39.3 37.7 36.5 - - - … AT45 7.4 7.5 12 12.1 10.3 10.4 10.5 10.6 10.7 10.3 9.7 9.2 8.9 - - - - - - - - - - … B744 122.9 123.3 128.6 129.5 130.3 131.2 183.3 184.1 184.8 185.6 186.3 187 187.6 184.1 179.4 170.9 163.7 158 153.8 151.3 150.2 150.5 152.2 70 numbers and belongs to a higher order differentiability class. In addition, trigonometric functions (sine function is used) are smooth function (they belong to differentiability class ); the range of the function is always [ ]make it very easy to implementation in the model. Based on those remarks, three pairs of symmetric reference functions have been selected (Figure 37,  1,0x ).   xxf I sin   xxf II 2sin   xxf III 2sin2 Figure 37: Boundary curves Each homotopy is controlled by independent parameter  (Figure 38). Weighted average of all homotopies gives resulting function – horizontal route (weighted line on Figure 39). Figure 38: One possible homotopy Figure 39: Resulting horizontal route Such route generation model offers a lot of freedom and flexibility in the choice of the route with a reasonable increase of problem complexity. Nevertheless, it is possible to reduce or increase the number of parameters, or to choose different reference functions. 1 = 0.65 2 = 0.60 3 = 0.52 71 4.1.4. Flight interaction calculation In the context of the following examples, flight interaction is calculated as a function of time separation. Interaction is computed for every pair of “conflicting” points identified using grid-based scheme described in chapter 3.3. (point that are separated in 3D less than given norms). There is no interaction between “conflicting” points if their time separation is larger than a given maximum separation and the interaction exponentially increases with decrease of time separation. The following function expresses the magnitude of interaction: { Parameter controls steepness of the exponential function (penalization of larger interaction) and bounds flight interaction. Both parameters might vary depending on the scenario. According to the definition, there could be multiple conflicting points between two flights and all are taken into account. This is reasonable to take into account as flight paths are curves that could intersect multiple times or be closely parallel. A number of conflicting points quantifies interactions between those two flights. Total flight interaction depends on both magnitude and quantity of interaction between flight points. 4.2. An academic test example In order to demonstrate the proposed RTP model, a simplistic academic test example was constructed. This model is solved using an exact optimization method (total enumeration search) which requires discretization of continuous search space. Due to high combinatoric, the model used in the hypothetical example is simplified from the one detailed in the chapter 2.7. The following are the main differences between the models:  no flight level change i.e. the model is tested in 2 space (horizontal flight); 72  high decision variable discretization steps (delay step to 2 minutes, homotopy coefficient to 0.1);  one symmetric homotopy in the route generation model – as the result, a set of possible trajectories is strongly reduced thus reducing model application only to the academic test example not real world problems;  linear cost function. 4.2.1. Traffic data – Demand Although simplified, model CPU running time still represents a limiting factor of the problem size. Therefore, test problem considers one air traffic control sector with four flights (Figure 40). Figure 40: Academic test example It is assumed that all flights have equal departure time, flight speed (450kts) and altitude level (FL350). Aircraft are initially located (departure point) on a circle of radius 150NM converging towards the circle center and flying to the diametrically opposed point (destination point). Length of preferable (initial) flight routes is 300NM. Although hardly ever expected, such challenging situation with high interaction between all flights is often seen in literature to test the capabilities of a model [Durand et al., 1996], [Dougui et al., 2013], [Delahaye et al., 2010], [Roussos et al., 2008]. This problem is known as a roundabout test problem because solution to this problem results in a roundabout like pattern on a circle center. To illustrate, space solution of the RTP model for eight flights test example is shown on Figure 41. 0 50 100 150 200 250 300 0 50 100 150 200 250 300 (3) (4) (2) (1) 73 Figure 41: Roundabout test problem Elementary time period of 15 seconds is assumed. Each flight is represented as a sequence of “route points” sampled every 15 seconds from origin to destination. Each point contains information about flight attributes including 4D coordinates (latitude, longitude, FL and time) of a point, current speeds, flight status (climb, descent or cruise) and flight intentions. Two successive points define flight segment as a direct portion of flight. 4.2.2. Test scenarios and numerical results of academic test example This section presents numerical results and discussion of several scenarios. Scenarios are carefully designed to illustrate insight model behavior depending on delay, fuel and interaction unit cost, aircraft type that operates certain flights, etc. 4.2.2.1. Base scenario In the base scenario all flights are operated by aircraft of the same type (A320). is set to 10 minutes. Delay unit cost of 30 EUR/min for the aircraft of “medium” cost category is assumed and total arrival delay at destination point is calculated. Fuel cost of 0.6 EUR/kg is assumed. For base scenario, a very high interaction unit cost of 500 EUR/unit is chosen resulting in no interaction and maximum robustness of the solution. It is important to note, although flight interaction is not constrained, for this specific example with lots of ‘free’ space for flight rerouting and due to high interaction unit cost it is cheaper to modify flight’s trajectory (change route and departure time) than have interaction between flights. 0 50 100 150 200 250 300 0 50 100 150 200 250 300 0 50 100 150 200 250 300 0 50 100 150 200 250 300 74 However, the value of interaction unit cost that gives no interaction will be different for different traffic demand. The results of model optimization for base scenario are given in Table 1. Table shows ground delay, alternative horizontal route extension compared to nominal one, total delay and total cost of actions applied for every flight. Since there were no interaction between flights, value of the objective function equals to the total cost of actions (EUR 1,100) taken to resolve initial flight interactions. A reference initial value of objective function, when all flights take nominal direct routes, is evaluated to EUR 167,000. This is due to high flight interaction and high flight interaction unit cost and does not represent real costs to airspace users. Table 1: Base scenario results Ground delay [min] Route extension Total delay [min] Total cost of actions [EUR] [NM] [% of route length] Flight 1 10 1.84 0.61 10.25 313.39 Flight 2 1 16.02 5.34 3.25 150.56 Flight 3 10 41.95 13.98 15.75 608.08 Flight 4 0 7.27 2.42 1.00 53.58 Total 21 67.08 - 30.25 1125.61 Average 5.25 16.77 5.59 7.56 281.40 Table 1 and Figure 42 represent one possible solution of the problem, although due to the symmetry of the problem and linear cost function several optimal solutions exist. Figure 42: Base scenario - optimized trajectory in 2D 0 50 100 150 200 250 300 0 50 100 150 200 250 300 (1) (2) (3) (4) 75 Another consequence of linear cost function is that some flights are penalized more than another. Flight 3 alone bears more than half of the total cost of actions applied to all flights. A more even distribution of cost among flights could be achieved using nonlinear cost function. 4.2.2.2. Low interaction cost scenarios Solution of base scenario represents only one Pareto optimum of this multi- objective problem. Change in weight coefficients of individual objective functions will possibly result in a different Pareto optimal solution. Since interaction cost is very hard to perceive and not possible to measure, providing decision makers with a set of different solutions rather than one solution is very encouraging. Therefore several optimizations of model are performed varying interaction unit cost. Except value of interaction unit cost, other parameters are not changed compared to base scenario. Table 2 shows optimal solution for all considered values of interaction unit cost. Every solution is represented by remaining flight interactions and associated costs; total delay, route extension and total costs of actions taken to resolve initial interaction problem. Last column represent value of the objective function for the solution and it is the sum of total cost of interaction and actions. Interaction is calculated as the number of “conflicting” pair of flight points divided into five interaction categories: (0-2], (2-4], (4-6], (6-8], (8-10) minutes of time separation. Table 2: Low interaction cost scenario results S ce n ar io Interaction unit cost Pareto optimum Interaction Total cost of interaction Delay [min] Route ext. [NM] Total cost of actions Objective function 0-2 2-4 4-6 6-8 8-10 Real Indicat. 1 50 0 0 0 0 0 0.00 0 18 81.26 1142.77 1142.77 2 40 0 0 0 0 0 0.00 0 18 81.26 1142.77 1142.77 3 30 0 0 0 0 10 81.2 130 20 61.59 1055.43 1136.63 4 20 0 0 0 0 21 113.69 273 20 52.84 988.45 1102.14 5 10 0 0 12 28 29 314.05 1437 10 38.66 581.30 895.35 6 5 0 0 37 23 14 206.39 2057 10 26.97 500.92 707.31 7 4 0 10 27 54 32 259.02 3213 7 21.54 400.74 659.76 8 3 0 34 21 36 34 245.89 4085 8 12.79 333.76 579.65 9 2 11 102 14 15 8 329.78 8274 0 18.22 133.95 463.73 10 1 69 52 28 2 23 175.22 11451 0 7.36 107.16 282.38 11 0.5 125 21 23 3 0 146.87 14698 0 1.84 13.39 160.26 12 0.25 167 0 0 0 0 83.50 16700 0 0.00 0.00 83.50 13 0.1 167 0 0 0 0 33.40 16700 0 0.00 0.00 33.40 76 As expected, high interaction unit cost results in no or little total interaction while interaction increases with decrease in unit cost. Although total flight interaction is increasing (see number of conflicting point in the table), the cost of total interaction does not, since unit cost is decreasing in parallel. To compare different scenarios and help in perceiving interaction level of given scenario, indicative value of total interaction is shown in the table, in addition to its real value. Indicative value of total interaction cost is calculated using unique unit cost for all scenarios no matter what real unit cost is used in optimization. Contrary, total cost of actions, taken to resolve flight interaction, decrease with the decrease in interaction unit cost. Figure 43 shows Pareto frontier of possible solutions of considered academic test example. Each point in the graph represents a single scenario from Table 2. Choice of “best” solution is left to decision maker based on information provided. From the figure, it is possible to conclude what is the minimal price needed to be paid if maximum robustness (minimum interaction) is required, what is the lower limit of action cost needed to be applied for desired robustness, as well as trade-off between solution robustness and its costs. Figure 43: Pareto frontier From the Figure 43, one can note that Pareto front is very steep for high values of solution robustness (left part of the graph). This implies that the total cost of actions is decreasing far faster than solution robustness, and it is possible to drastically reduce costs of the solution slightly sacrificing solution robustness. Scenario 5 that is emphasized in Table 2 and marked at the Figure 43 is a good 77 candidate for “best” solution. In this scenario, at the expense of the small interaction increase, it is possible to reduce cost of the solution by half. All flights in this scenario are separated more than 5 minutes, thus not significantly affecting the solution robustness. Table 3 shows resulting solution for scenario 5. Total delay at the destination airport (including ground and airborne delay) of all flights is 15.25 minutes in scenario 5 compared to base scenario that was 30.25 minutes (Table 1) giving 50% reduction. There is less but still significant decrease in route extension (38.66 NM compared to 67.08 NM). Table 3: Low interaction cost – Scenario 5 results Ground delay [min] Route extension Total delay [min] Total cost of actions [EUR] [NM] [% of route length] Flight 1 0 27.71 9.24 3.75 200.94 Flight 2 0 7.27 2.42 1.00 53.58 Flight 3 0 1.84 0.61 0.25 13.39 Flight 4 10 1.84 0.61 10.25 313.39 Total 10 38.66 - 15.25 581.30 Average 2.5 9.66 3.22 3.81 145.32 Information of total flight interaction in the Table 2 is given by a single aggregated indicator for all flights (total interaction cost). This value is used as an objective in the optimization. However, more detailed representation using a number of conflicting points per interaction category is easier perceived by humans, even if it is sometimes not sufficient to judge how robust some solutions are. Additional information such as: minimal time separation, geographical distribution of interactions, number of interaction zones and involved flights, etc. may help perceiving the whole picture. Figure 44 shows resulting solution of scenario 5. Interaction categories in the figure are graduated by line color and thickness. Red thick line depicts interaction category of (4-6] minutes separation, while lighter colors represent larger time separations. It is noticeable that most interactions involve flight 2 in combination with other flights. This gives very high solution robustness, since the most of the disturbances on the day of operation will mainly affect flight number 2 and small tactical actions on the very same flight could solve potential problems. 78 Figure 44: Scenario 5 - optimized trajectory in 2D 4.2.2.3. Constrained interaction scenarios This scenario represents another way to control tradeoff between solution cost and robustness. It considers changes of interaction limit, parameter , that controls the minimum time separation between two route points such that there is no interaction between them. The lower means lower time separation, resulting in a less robust solution in general. Opposite to previous scenarios, flight interaction is constrained to ensure that a given minimum time separation is maintained for all flights. This helps decision makers to perceive the level of robustness of each solution and to choose one that best suit the current situation. Table 4 shows optimum solutions of the problem for different values of . Scenario 0 is equivalent to Base scenario as the same value of is considered. When is decreasing the total interaction of all flight taking nominal trajectories is decreasing too, that requires less action to solve interaction problem. Hence, the total delay at destination airport, as well as the total cost of actions, is decreasing with . So, the largest robustness is given for scenario 0 when maximum value of is considered, resulting in the highest cost of the solution. On the other hand, scenario 10 gives no robustness but require no additional costs either. 0 50 100 150 200 250 300 0 50 100 150 200 250 300 (1) (2) (3) (4) 79 Table 4: Constrained interaction scenario results Scen. [min] Pareto optimum Ground delay [min] Route extension Total delay [min] Total cost of actions [EUR] [NM] [% of route length] 0 10 18 81.26 6.77 29.25 1142.77 1 9 18 52.84 4.41 25.25 928.45 2 8 18 33.88 2.82 22.75 794.51 3 7 12 45.25 3.52 17.75 668.08 4 6 8 45.63 3.80 14.25 574.87 5 5 6 34.98 2.92 10.75 434.50 6 4 6 23.29 1.94 9.25 354.13 7 3 4 17.86 1.49 6.50 253.95 8 2 0 23.65 1.97 3.25 174.13 9 1 0 7.36 0.61 1.00 53.58 10 0 0 0.00 0.00 0.00 0.00 Results of scenarios are presented in the Figure 45. Horizontal axis denotes solution robustness with the given value of . Vertical axis denotes solution costs. Each point in the graph represents one scenario. In the Figure 45 we see an almost linear relation between solution cost and robustness for this particular example. Although it is not by the nature, the graph may be interpreted as a Pareto frontier showing minimum cost of action for a certain level of required solution robustness. Figure 45: Set of solutions depending on required robustness 80 4.2.2.4. Different aircraft type scenario Here, we compare solutions of two scenarios: first where all flights are operated by aircraft of the same type (A320 as in the case of base scenario) and second where one flight (Flight 1 in the example) is operated by larger (jumbo-jet) aircraft Boeing 747-400 and all others use a base aircraft A320. Boeing 747-400 belongs to “jumbo” flight cost category with delay unit cost assumed to be 80 EUR/min. According to BADA (BADA code B744), FBR equals 158 kg/min. In comparison, delay unit cost of the A320 is assumed to be 30 EUR/min (“medium” cost category) and FBR is 39.3 kg/min. is set to 9 minutes, meaning that separation of 9 minutes between intersecting point of two flight is needed such that there is no interaction between them. The interaction unit cost is again set to high value (500 EUR/unit) resulting in no interaction between flights after optimization. In the discussion of base scenario it was concluded that, due to the use of linear cost function, actions applied and their costs are not evenly distributed across all flights. This implies that some flights are penalized more than others. When using different unit cost, one can expect that this difference will become even larger. Table 5 shows the results of these scenarios: Table 5a with all same aircraft, and Table 5b where flight number 1 is operated with B744. The first thing to note is that the total cost of actions taken to resolve the interaction problem is larger in the second scenario. This is expected, even if same actions are taken in both scenarios, because of the change in unit costs of one flight in the second scenario. Table 5: Different aircraft type scenario results Ground delay [min] Route extension [NM] Total delay [min] Total cost of actions [EUR] Flight 1 0 16.02 2.25 120.55 Flight 2 8 1.84 8.25 253.39 Flight 3 0 7.27 1.00 53.58 Flight 4 10 27.71 13.75 500.92 Total 18 52.84 25.25 928.45 Average 4.5 13.21 6.31 232.11 Ground delay [min] Route extension [NM] Total delay [min] Total cost of actions [EUR] Flight 1 0 1.84 0.25 43.70 Flight 2 8 1.84 8.25 253.39 Flight 3 0 27.71 3.75 200.92 Flight 4 10 27.71 13.75 500.92 Total 18 59.10 26.00 998.94 Average 4.5 14.77 6.50 249.73 a) same aircrafts A320 b) flight 1 – B744, others A320 81 More important is the fact that actions taken to resolve interaction are also changed. This is influenced by higher unit cost of the flight 1 in scenario 2 resulting in less action applied to that flight. Flight 1 in scenario 2 receives only small rerouting resulting in 1.84 NM route extension and 0.25 minute delay at the destination, compared to flight 3 in scenario 1 (flight that received the minimum actions) that has 7.27 NM route extension and 1 minute of ground delay. Reduction of actions applied to flight 1 in scenario 2 is at the expense of other flights. Therefore we have seen an increase in total route extension and total delay of all flights in scenario 2. Total route extension is 59.10 NM and total delay is 26 minutes in scenario 2 compared to 52.84 NM and 25.25 minutes respectively in scenario 1. 4.3. Real-world test example – French airspace Having illustrated basic idea of the proposed model for small hypothetical example, the next step is to test the model on large scale problem to see its capabilities. 4.3.1. Experimental setup 4.3.1.1. Traffic sample Traffic data included traffic in French Metropolitan airspace, more precisely flights that originate and/or are destined from/to an airport in French airspace and flights that overfly French airspace. August is chosen as one of the busiest months in the year and specific date (17.08.2008) has been chosen based on flight data availability. On given date there were 8845 flights in French airspace. To make the problem size manageable but still realistic three hour period (9h – 12h) from morning peak are selected. Traffic sample has 1755 flights in the French airspace of which: 172 (about 10%) were domestic flights, 700 (about 40%) were international flights (having either departure or arrival at French airports) while 883 or about 50% of all flights were overflights. Departures and arrivals from/to French airports amount to approximately 500 each, of which Paris Charles de Gaulle (LFPG) had the biggest share (in total about 300 operations). Flights were operated by 50 different aircraft types. Distribution of flight cost categories based on aircraft type is shown on Figure 46. 82 Figure 46: Distribution of flight cost categories Departure time (or the time when flight enters French airspace), departure and arrival points for each flight are extracted from real traffic data after flights were simulated using direct flight routes. Initial direct flight routes are shown on Figure 47. Figure 47: Input traffic data – French airspace 4.3.1.2. Traffic data source and format Data were obtained from the ENAC Traffic Simulator (CATS). The core of the CATS system is an en-route traffic simulation engine that uses a tabulated model of aircraft based on the BADA Eurocontrol performance model. Beside, CATS contains lots of different modules including conflict detection and resolution modules. It is lightweight and well suited for fast-time simulations as it takes only 15-30 minutes to simulate one day of European traffic [Alliot et al., 1997]. As an input it can use data from the French CAUTRA system, the CFMU European database (DDR) and many others. 83 The flight data were provided in “.cats” file format (Figure 48). As CATS is discrete model, each flight is represented as a sequence of “route points” sampled every 15 seconds from origin to destination. Each point contains the 4D position and velocities of aircraft with aircraft intention. Two successive points define flight segment as a direct portion of flight. Mini flight plan Number of point $ 1000 27:20:00 28:26:15 330 458 BIE2504 A321 LFPG GCFV 1000 NbPlots: 266 Flight points 27:20:00 -650 176169 -134 -329 125 2491 1 27:20:15 -686 176087 -135 -332 132 2445 1 27:20:30 -723 176005 -136 -335 138 2399 1 Figure 48: CATS file format Flight coordinates are given in Lambert conformal conic projection usually used for aeronautical charts because a straight line approximates an orthodrome (great-circle route). 4.3.1.3. Solution space Each flight is attached five vertical profiles corresponding to different flight levels. One nominal profile for nominal flight level and four alternative profiles: two at lower and two at higher flight level. For given flight levels flight is re- simulated and optimal vertical profiles are kept as inputs. Figure 49 shows five vertical profiles for one flight. Each profile has unique Top of climb (TOC), Top of descend (TOD) and different time of arrival to destination point that is same in all cases. The last is due to fact that climb and descent performance are respected and taken into account. Figure 49: Alternative vertical profiles for one flight 84 Although every alternative profile is optimal from the aircraft performance point of view (respect climbing and descending rates), only one profile can be optimal for specific origin-destination pair and given speed. Therefore other profiles incur additional cost to the aircraft operator. For each of the available alternative routes the cost was calculated as a sum of additional operating cost and additional fuel cost. Having in mind that aircraft operating cost is approximately a linear function of flight duration, for a given aircraft type [Swan and Adler, 2006], additional operating cost is computed by multiplying the additional flight duration (compared to flight length using the nominal profile) by the unit operating cost33 for the considered aircraft type. Additional fuel cost includes cost of additional fuel burned due to use of different flight level. Fuel burns rates were extracted from EUROCONTROL’s Advanced Emission Model – AEM that are based on BADA (Base of Aircraft Data) datasets as explained in chapter 4.1.1. Although both additional operating and fuel cost can be negative (not at the same time), the additional cost of the alternative route is always positive as it is assumed (and chosen in that manner) that nominal profile corresponds to the optimal flight level for that specific flight. The calculated initial cost of an alternative profile is stored as an input parameter attached to each profile. The horizontal route design is based on route generation model controlled by homotopy as explained in chapter 2.4 and using three symmetric homotopy functions respect to reference functions that are set in chapter 4.1.3. For this particular example maximum horizontal route offset from direct route is limited to 20 percent of direct route length. Ground delay is set to maximum 30 minutes. 4.3.1.4. Interaction parameters settings Due to the size of the problem, flight interaction is bounded to 3 minutes ( ), meaning that flights separated more than 3 minutes in time will not be in the interaction even if their 3D routes (not taking time into account) are not spatially separated. Flight interaction is calculated for the en-route phase only. Route points above 10.000ft or FL100 are considered to belong to en-route phase. This assumption is imposed for several reasons. First, as mentioned by [SESAR, 2013], deregulation will only consider en-route part of the airspace while traffic will remain structured in the vicinity of airports (terminal areas – TMA’s). Second, 33 Unit costs are taken from [Eurocontrol, 2012b] that are based on ICAO Base-line Aircraft Operating Costs 85 since all arriving flights for a single airport join at that airport runway or metering point and all departing flights disperse from it, reserving e.g. 3 minutes for each flight is operationally unfeasible neither preferred as it will significantly reduce airport throughput. Third, taking into account that departure and arrival points are static and given as input for all flights, spatial solution of the interaction problem by horizontal route and/or vertical profile change is only effective in en- route phase. Consequently interaction at departure airport and arrival airport and/or point (airborne fix) may only be resolved in time. Situation at departure airport is further restricted as interaction can be alleviated by controlling departure slot only, while there is no means to resolve interaction at departure point without collaboration with neighbor control centers. Figure 50: Flight interaction for initial (nominal) routes All “conflicting” points are further divided into six interaction categories depending on time separation (width of interaction equals to half minutes). In Figure 50 flight interaction categories are depicted using different colors grading from dark-red for highest interaction value to yellow for lowest interaction. Blue color represents part of the flight where no interaction with other flights is experienced. RTP is formulated as an unconstrained optimization problem since the static penalty method is used to control the interaction of the solution. If penalties are large enough, then the global minimum to the unconstrained problem is the solution to the original problem with constrained interaction. The use of large penalties, however, leads to rugged search terrains and deep local minima, making it difficult for global-search methods to escape from it [Bertsekas, 1982]. 86 Although RTP does not intend to constrain interaction but rather to find the best tradeoff between solution robustness and costs, selecting a suitable penalty imposes the same difficulties. A nonlinear interaction cost function is used as explained in chapter 4.1.4. Parameter controlling steepness of the exponential function is set to 0.9 far penalizing high interaction classes. The interaction unit cost is set to 1.000 EUR/unit. This value is experimentally chosen. If higher unit cost is selected one might expect greater robustness of the solution and higher operational cost of the solution. Using given parameters, nominal routes results in EUR 92.66 million of total initial interaction. This is just value of the objective function and it is not relate to the real cost of airspace users. 4.3.1.5. TMA modelling As explained, traffic remains structured in the TMAs and controllers will remain responsible for aircraft separation. Each controller can safely manage a limited number of flights and traffic in TMA should remain lower than TMA capacity at all time. In the RTP TMA is modelled as airspace with a given radius around the airport up to FL100 (Figure 51). Capacity of TMA is defined as the maximum number of flights that can simultaneously be in TMA in one hour. Although it is similar to the general definition of airspace capacity, the capacity of TMA in RTP takes into account airport runway system capacity. Figure 51: TMA airspace Capacity constraints are modeled through capacity violation costs, similar to approach for controlling flight interactions, and a static penalty is assigned to each violation of TMA’s capacity. For illustrative purpose, only one TMA is considered in this example, namely Paris TMA. It is centered on Paris Charles de Gaulle (LFPG) airport with 20NM in radius. TMA’s hourly capacity is considered constant over time and set to 110 flights. Capacity violation unit cost is set to EUR 10.000. Very high unit cost is chosen in order to constrain TMA capacity violations. ap R FL100 87 Figure 52 presents the initial demand for a Paris TMA when all flights use nominal routes. There is one period when traffic demand exceeds TMA capacity, and another one that is close to capacity level. This capacity violation accounts for an additional EUR 10.000 in costs. Objective function of initial (nominal) routes equals to EUR 92.67 million and it is the sum of total initial interaction and total capacity violation costs as no additional costs are experienced when nominal routes are used. Figure 52: Initial traffic demand at Paris TMA 4.3.1.6. Scenarios settings In this work two main scenarios are defined and tested on the given traffic over French airspace. First, the base scenario was tested for finding a robust solution as a balance between solution interaction and additional flight costs. Following unit costs are used in addition to unit cost of interaction and capacity violation that are set previously. Delay unit costs of light, medium, heavy and jumbo aircraft are assumed to be 15, 30, 60 and 80 EUR/min respectively; and fuel cost is set to 0.6 EUR/kg34. Then, one additional scenario is designed in order to evaluate robustness of the solution. Second scenario aims to find conflict-free solution not taking into account solution robustness i.e. interaction boundary is set to zero meaning that it is sufficient that flight are separated in space or in time. 34 This value is taken as long term average price. Data are collected from IATA jet fuel price monitor (http://www.iata.org/publications/economics/fuel-monitor/) that is sourced by leading energy information provider Platts. 88 4.3.2. Experimental results 4.3.2.1. Conflict-free scenario Considering the conflicts only, the general problem is relaxed and solution to the problem is easier to find. The initial number of conflicting points between nominal flight trajectories using direct routes was 18.268. In addition there was one capacity violation in the Paris TMA, resulting in total initial cost of EUR 18,278 million. The algorithm was able to find conflict-free solution and the best found solution had the value of the objective function equal to EUR 29.984. The best solution found without conflicts is presented in Figure 53. This solution resulted in no capacity violation in the Paris TMA as shown in Figure 54. Figure 53: Conflict-free trajectories Figure 54: Traffic demand at Paris TMA for conflict-free trajectories In the proposed solution 45% of all flights were modified. Most of the flights had their horizontal route changed (40%), about 15% received ground delays and around 35 flights or 2% are required to change their nominal flight level. The total route length extension was 0.25% with an average of 0.58% per modified flight. This results in 250 minutes of total en-route delay caused by flight rerouting. As a reference, today best value of annual flight inefficiency achieved was around 3.5% as shown in chapter 1.1. The total ground delay was 342 minutes with maximum ground delay of 22 minutes and average of 1.4 minutes per delayed flight. In addition some flights received a change of the flight level. Average flight level change received was 1.19 per affected flight meaning that the 89 majority of flights changed their nominal flight level by one level (±1 FL) and only few received change of ±2 FL (that was maximum level change possible). Total additional operating costs were EUR 29.984, giving an average of EUR 17 per flight. Most flights (1051) had not been modified and their initial nominal trajectories remained unchanged. Other flights received additional costs ranging up to EUR 680 with an average of EUR 41.5. Distribution of additional operating costs is shown on Figure 55. 90% of flights are in the range of ±50 EUR from the average cost. One might say that the rest 10% of flights are penalized as they receives an unfair portion of costs compared to other flights. Figure 55: Distribution of flight additional operating costs – Conflict-free scenario For this scenario the best solution of the problem is obtained in 18 hours of CPU time on the Intel Pentium Dual-Core 2.9GHz with 4Gb of RAM. As the number of iterations at each next step is increased, for better local search performance, the solution is slowly converging as shown at Figure 56. The solution with objective value in the range of 10% of the objective of the best known solution was found in less than 200 steps in 7 hours of CPU time. For the reference, if flight operating costs are excluded from the objective, and only number of conflict is counted, conflict-free solution is found in less than half of minute. 90 Figure 56: Objective function convergence – Conflict-free scenario 4.3.2.2. Base scenario In the search for the best robust solution, optimization algorithm was able to reduce the value of the objective function starting from initial 92.67 to 0.19 million of EUR. Not all flight interactions are solved meaning that some flight which are not separated in space37 were separated with less than 3 minutes (value of ). Remaining interactions result in EUR 4 thousand of costs in the objective value. Those “conflicting” points are located near Paris TMA where the densest traffic was experienced (marked in Figure 57) and include interaction between 20 flights in two locations. However, all except few flight points were separated more than 2 minutes. Change of unit costs will possibly result in different balance between cost of interaction and flight operational cost, and it is a way to adjust required solution robustness. 37 Separation minima are set to 5nm horizontal and 1000ft vertical distance as explained in 0. 91 Figure 57: Robust flight trajectories Figure 58: Traffic demand at Paris TMA after optimization No capacity violations are experienced in the Paris TMA after optimization as shown at Figure 58. By comparing initial traffic demand (Figure 52) and traffic demand after optimization (Figure 58) in the Paris TMA, one could note that the total traffic differs (337 compared to 338). This is due to the fact that TMA traffic includes overflights – flights that depart/arrive from/to nearby airports (Orly - LFPO in this case). With a choice of horizontal route, those flights could fly through or to be rerouted around TMA, increasing or reducing total TMA’s traffic. It is interesting to note that the total traffic in TMA was increased in this example. Taking into account that interactions close to Paris TMA and capacity violation of the TMA represent the greatest problem in the example this is opposite than one could expect. To find robust solution, as it was expected, many flights need to change their nominal trajectories (87% of all flights). Horizontal route change has been used far to alleviate flight interaction (more than 85% of all flights affected). Around 30% and 23% of flights receive ground delay and vertical profile change respectively. The total route length extension was 0.8% with maximum flight extension of 8%. This results in 781 minutes of additional flight time due to flight rerouting. In addition total ground delay was 2743 minutes with maximum ground delay of 27 minutes and an average of 5 minutes per delayed flight. Taking into account all sources of delay, average delay at the destination compared to nominal schedule (flight plan) is less than 2 minutes per flight. This is a very promising result considering that no additional disruptions are likely to happen because of solution robustness. Average flight level change was close to 92 1.5 (precisely 1.34) meaning that half of flights were altered to vertical profiles that are ±1 flight level from nominal level, and half to one that are ±2 flight levels. Total additional operating costs of the solution were EUR 192.417 with average of EUR 108 per flight. 363 flights remained unchanged and received no additional costs. More than half of all modified flights had a cost that was less than EUR 100. Average cost received was EUR 138 per modified flight with maximum of EUR 1288. Distribution of additional operating costs is shown on Figure 59. Almost 90% of flights had cost in the range (0-300] that is ±150 EUR from the average cost. The rest 10% of flights receive an “unfair” portion of costs compared to other flights. Figure 59: Distribution of flight additional operating costs – Base scenario Searching for robust solution was a very challenging task. Although in both scenarios (conflict-free and base) objective value of the optimal solution was not known, some initial estimation are possible in the conflict-free scenario exploiting the fact that optimal solution has to be conflict free and without terminal capacity violations. However the base scenario aims to find a balance between interaction and operating cost and as it depends on many factors, such as: traffic level and density, unit costs, state space, etc. it is not possible to estimate the optimal level of interaction in the solution. Once more stopping criterion for optimization process was temperature level. The best solution was obtained in 28 hours of CPU time, on the Intel Pentium Dual-Core 2.9GHz with 4Gb of RAM, in two stages. First stage (marked in black on Figure 60) that last 20 hours represents simulated annealing 93 optimization process. To further improve solution heuristics based on greedy search was run in the second stage and last additional 8 hours (marked in gray on Figure 60). Greedy algorithm only accepts better solution allowing intensification around local minimum. The solution with objective value in the range of 15% of the objective of the best known solution was found in less than 400 steps in 10 hours of CPU time. Excluding flight operating costs from the objective, it is possible to find a robust solution in one hour on average. Figure 60: Objective function convergence – Base scenario 4.3.3. Solution robustness testing To test robustness, solutions obtained in conflict-free and base scenarios are imposed for additional flight delays, and effects of these disturbances are measured and compared with congestion problems they produce to the system. The magnitude of the congestion problem is measured by the number of conflicting points it produces and the number of flights included in the conflicts. Numbers of scenario are defined taking conflict-free and robust trajectories as input. Number of delayed flights and delay magnitude varies from scenarios, and ranges from small to big disturbances. Taking into account that only portion of flights during a given day is considered, imposing long delays to some flights does not always yield worse congestion problem, as it might rescheduled these flights to the period with low traffic. Therefore, only small delays are considered (1.5 minutes) and size of the disturbance is controlled by the number of affected flights (Table 6). 94 Table 6: Test scenarios settings Scenario Number of affected flights Delay 1 20 1.5 minutes 2 50 1.5 minutes 3 100 1.5 minutes 4 500 1.5 minutes 5 (chaos) 1000 1.5 minutes The resulting congestion problem is very sensitive to the selected flights, the ones affected by the delay. Delaying a flight, not having interaction with other flights, will certainly not cause additional problems. For this reason, several groups of affected flights are identified for each scenario, then minimum, maximum and average value of the conflicting points and the conflicting flights are recorded. Table 7 shows the result of robustness tests based on the number of conflicting points, while Table 8 shows results based on the number of conflicting flight. Both tables clearly show that robust trajectories are less affected by the disturbances and therefore cause less congestion problems. Newly generated congestion problems are easy to solve at tactical level compared to one that resulting from the usage of conflict-free trajectories. Table 7: Resulting congestion problem – conflicting points Scenario 1 Scenario 2 Scenario 3 Scenario 4 Scenario 5 C o n fl ic t fr ee min. 120 302 588 3624 6000 max. 428 828 1328 4525 6966 avr. 217.2 473.2 1000.8 4064.4 6511.6 R o b u st min. 0 22 106 924 1830 max. 132 202 254 1556 2494 avr. 38 81.2 176 1230.8 2182 Table 8: Resulting congestion problem – conflicting flights Scenario 1 Scenario 2 Scenario 3 Scenario 4 Scenario 5 C o n fl ic t fr ee min. 20 52 85 388 571 max. 331 65 127 417 605 avr. 87 57.2 107.4 400.4 584.2 R o b u st min. 0 7 23 127 206 max. 12 21 37 155 256 avr. 5.4 12.2 28 138.2 231.8 95 Distribution of average number of conflicting points and flights over scenarios are illustrated in the Figure 61. Both figures show almost the same ratio between conflicting points and flights. Figure 61: Distribution of conflicting points and flights Figure 62 illustrates one congestion problem for the scenario 5. Left image shows congestion caused by disturbance of the conflict-free trajectories, and right of the robust trajectories. Conflicting points are marked by red color in the figure. Mentioned conclusions are more than obvious in the figure. Figure 62: Resulting congestion problem 96 5. CONCLUSIONS This thesis investigates the potential of robust trajectory planning – RTP (considered as an additional demand management action) at pre-tactical level as a mean to alleviate the en-route congestion in airspace. The most of existing trajectory planning models consider finding of conflict-free trajectories without taking into account uncertainty of trajectory prediction. It is shown in the thesis that in the case of traffic disturbances, it is better to have a robust solution otherwise newly generated congestion problems would be hard and costly to solve. RTP puts forward trajectory planning involving generation of congestion- free trajectories with minimum operating cost taking into account uncertainty of trajectory prediction and unforeseen event. Although planned cost could be higher than of conventional models, the underlying rationale is to reduce cost of disruptions adding robustness to schedules. As consequence, this might lead to reductions in operating cost. This thesis introduces a novel approach for route generation (3D trajectory) based on homotopic feature of continuous functions. This approach is capable of generating a large number of routes of different shape with a reasonable number of decision variables. Those routes are then coupled with time dimension in order to create 4D trajectories. RTP problem is modeled as a mixed-variable optimization problem and it is solved using stochastic optimization technique (Simulated Annealing). Application of developed optimization model is illustrated in the two examples: academic and real-life. Results show that the model is able to solve real instances of the problem, with computation time which corresponds to the intended use of the model (strategic, pre-tactical level). Further the results indicate that, under certain conditions, solution robustness could be considerably increased at the relative small expenses of solution costs providing a good alternative to the solutions developed by existing conflict-free trajectory planning models. A principal conclusion arising from the performed experiments is that the homotopy route generation model is very efficient because it provides large freedom in the route generation by controlling only few real value parameters. However, a performance of the model highly depends on the choice of reference (primitive) functions and additional research (improvement) could be done in the 97 future. Although it increases problem complexity, flight rerouting in the vertical dimension was found to be very effective in reducing flight delay and horizontal route extension. Finally, introduction of operating costs into objective function significantly increases problem complexity. The results of the optimization model are very sensitive to the choice of the static penalty assigned to traffic interaction. The use of large penalties leads to rugged search terrains and deep local minima, making it difficult for global-search methods to escape from it. However, it was shown that it is possible to find robust flight trajectories, reducing flight interaction, and keeping cost induced to the user to an acceptable level. Due to high combinatoric and the rugged shape of the objective function in the search space, slow exploration of the search space is inevitable. This influences computation time that should be further improved in order to address larger problem instances. Furthermore, a feature of the homotopy route generation model, that horizontal route is controlled by real value parameters, could be more exploited in future research using field congestion metrics. It is foreseen that in such a case, it would be possible to represent congestion metrics as a function of parameters that controls each homotopy, and not to calculate its value throughout the simulation. 98 REFERENCES - A. Agustin, A. Alonso-Ayuso, L.F. Escudero and C. Pizarro (2012), “On air traffic flow management with rerouting. Part II: Stochastic case”, European Journal of Operational Research, Vol. 219, 2012, pp. 167–177 - J. Alliot, H. Gruber, G. Joly and M. Schoenauer (1993), “Genetic Algorithms for solving Air Traffic Control conflicts”, CENA internal report, 1993 - J. Alliot, J. Bosc, N. Durand and L. Maugis (1997), “CATS: A Complete Air Traffic Simulator”, Toulouse, France, 1997 - G. Andreatta, G. Romanin-Jacur (1987), “Aircraft Flow Management under Congestion”, Transportation Science, Vol. 21, No. 4, November 1987, pp. 249-253 - G. Andreatta, L. Brunetta (1998), “Multiairport Ground holding problem: A computational evaluation of exact algorithms”, Operations Research, Vol. 46, No. 1, January-February 1998, pp. 57-64 - P. Averty, B. Johansson, J. Wise and C. Capsie (2007), “Could ERASMUS speed adjustments be identifiable by air traffic controllers?”, In Proceedings of the 7th USA/Europe ATM R&D Seminar, Barcelona, Spain, July 2007 - O. Babic and F. Netjasov (2011), “Air Traffic Control”, textbook (in Serbian language), Belgrade: Faculty of Transport and Traffic Engineering, 2011, ISBN 978- 86-7395-289-5 - M.O. Ball, R. Hoffman, A. Odoni and R. Rifkin (1999), “The static stochastic ground holding problem with aggregate demands”, Internal report, 1999 - M.O. Ball, R. Hoffman, A. Odoni and R. Rifkin (2000), “A stochastic Integer program with dual network structure and its application to the Ground holding problem”, Operations Research, Vol. 51, No. 1, February 2003, pp. 167-171 - N. Barnier and C. Allignol (2012), “Trajectory deconfliction with constraint programming”, The Knowledge Engineering Review, Vol. 27, No. 3, 2012, pp. 291- 307 - D. P. Bertsekas (1982), “Constrained Optimization and Lagrange Multiplier Methods”, Academic Press, 1982 - D. Bertsimas, A. Odoni (1997), “A critical survey of optimization models for tactical and strategic aspects of air traffic flow management”, Technical report, NASA, 1997 99 - D. Bertsimas, S. Stock Patterson (1998), “The Air Traffic Flow Management Problem with enroute capacities”, Operations Research, Vol. 46, No. 3, May-June 1998, pp. 406-42 - D. Bertsimas and S. Stock Patterson (2000), “The Traffic Flow Management Rerouting Problem in Air Traffic Control: A Dynamic Network Flow Approach”, Transportation Science, Vol. 34, No. 3, August 2000, pp. 239-255 - D. Bertsimas, G. Lulli, A. Odoni (2008), “Air Traffic Flow Management Problem: An Integer Optimization Approach”, Integer programming and Combinatorial Optimization, Volume 5035, 2008, pp. 34-46 - D. Bertsimas, G. Lulli and A. Odoni (2011), “An Integer Optimization Approach to Large-Scale Air Traffic Flow Management”, Operations Research, Vol. 59, No. 1, January-February 2011, pp. 211-227 - L. Bianco, M. Bielli (1992), “Air Traffic Management: Optimization Models and Algorithms”, Journal of Advanced Transportation, Vol. 26, No. 2, Summer 1992, pp. 133-167. - L. Blasi, S. Barbato and M. Mattei (2011), “Flight Path Optimization Using Primitive Manoeuvres: A Particle Swarm Approach” Automatic Control in Aerospace, No. 1, March 2011, ISSN 1974-5168 - S. Chaimatanan, D. Delahaye and M. Mongeau (2012), “A methodology for Strategic Planning of Aircraft Trajectories using Simulated Annealing”, ISIATM, 1st International Conference on Interdisciplinary Science for Air traffic Management, Daytona Beach, USA, 2012 - S. Chaimatanan, D. Delahaye and M. Mongeau (2013), “Strategic deconfliction of aircraft trajectories”, ISIATM, 2nd International Conference on Interdisciplinary Science for Air traffic Management, Toulouse, France, 2013 - A. Cook, G.Tanner, V. Williams and G. Meise (2007), “Dynamic Cost Indexing”, In Proceedings of 6th EUROCONTROL Innovative Research Workshops & Exibition, Bretigny sur Orge, France, December 2007 - D. Delahaye, J.M. Alliot, M. Schoenauer and J.L. Farges (1994), “Genetic algorithms for air traffic assignment”, In Proceedings of European Conference on Artificial Intelligence, 1994 100 - D. Delahaye and A. Odoni (1997), “Airspace Congestion Smoothing by Stochastic Optimization”, In Proceedings of the Sixth International Conference on evolutionary programming, 1997, pp. 163-176 - D. Delahaye, C. Peyronne, M. Mongeau and S. Puechmorel (2010), “Aircraft Conflict Resolution by Genetic Algorithm and B-Spline Approximation”, EIWAC ENRI International Workshop on ATM/CNS, Tokyo, Japan, November 2010 - D. Delahaye and S. Puechmorel (2000), “Air Traffic Complexity: Towards intrinsic metrics”, 3rd USA/Europe ATM R&D Seminar, Napoli, Italy, June 2000 - S. Puechmorel and D. Delahaye (2009), “Dynamical Systems Complexity with a view towards air traffic management applications”, Joint 48th IEEE Conference on Decision and Control and 28th Chinese Control Conference, Shanghai, P.R. China, December 16-18, 2009 - D. Delahaye and S. Puechmorel (2013), “Modeling and Optimization of Air Traffic”, Wiley-ISTE, ISBN: 978-1-84821-595-5, June 2013 - D. Delahaye, S. Puechmorel, P. Tsiotras and E. Feron (2013), “Mathematical Models for Aircraft Trajectory Design: A Survey”, EIWAC, 3rd ENRI International Workshop on ATM/CNS, Tokyo, Japan 2013 - L. Delgado and X. Prats (2009), “Fuel consumption assessment for speed variation concepts during the cruise phase”, GARS Workshop on ATM Economics, University of Belgrade, September 2009 - P. Dell’Olmo and G. Lulli (2003), “A new hierarchical architecture for Air Traffic Management: Optimisation of airway capacity in a Free Flight scenario”, European Journal of Operational Research, Vol. 144, 2003, pp. 179–193 - Directorate Network Management (2012), “Monthly Network Operations Report”, September 2012 - N. Dougui, D. Delahaye and M. Mongeau, (2010), “A New Method for Generating Optimal Conflict Free 4D Trajectory”, Procedia - 4th International Conference on Research in Air Transportation, Budapest, June 2010 - N. Dougui, D. Delahaye, M. Mongeau and S. Puechmorel, (2012), “Aircraft Trajectory Planning Under Uncertainty by Light Propagation”, Procedia - Social and Behavioral Sciences Volume 54, October 2012, pp 201-210 101 - N. Dougui, D. Delahaye, S. Puechmorel and M. Mongeau (2013), “A light- propagation model for aircraft trajectory planning”, Journal of Global Optimization, Volume 56, Issue 3, July 2013, pp 873-895 - N. Durand and J.M. Alliot (1997), “Optimal Reslution of En Route Conflicts”, CENA internal report, 1997 - N. Durand, J.M. Alliot and J. Noailles (1996), “Automatic aircraft conflict resolution using Genetic Algorithms”, Proceedings of the Symposium on Applied Computing, 1996 - J. P. Dwyer and S. Landry (2009) – “Separation Assurance and Collision Avoidance Concepts for the Next Generation Air Transportation System”, in M.J. Smith and G. Salvendy (Eds.): “Human Interface, Part II”, HCII 2009, LNCS 5618, pp. 748–757, Springer-Verlag, Berlin Heidelberg, 2009 - A. J. Eele and A. Richards (2009), “Path-Planning with Avoidance Using Nonlinear Branch-and-Bound Optimization”, Journal of Guidance, Control, and Dynamics, Vol. 32, No. 2, 2009, pp. 384-394. - EUROCONTROL (2004), “Air Traffic Flow and Capacity Management Strategy”, Edition 1.2, April 2004 - EUROCONTROL (2009a), “Overall description of the platform and its capabilities”, Episode 3 WP 6 “Technological Enablers”, V 2.0, August 2009 - EUROCONTROL (2009b), “Base of Aircraft Data (BADA) aircraft performance modelling report”, EEC Technical/Scientific Report No. 2009-009, February 2009 - EUROCONTROL (2010), “Advanced Required Navigation Performance (A- RNP) Real-Time Simulation”, Airspace Validation Unit, CRDS/SIM/RTS-10067- FEU, 2010 - EUROCONTROL (2011), “Air Traffic Flow and Capacity Management Operations – ATFCM Users Manual”, Edition 15.0, March 2011 - EUROCONTROL (2012a), “Free Route Developments in Europe”, Edition number 1.0, February 2012 102 - EUROCONTROL (2012b), “Standard Inputs for EUROCONTROL Cost Benefit Analyses”, Edition number 5.0, January 2012 - European Commission (2001) - Report of the group of personalities, “European Aeronautics: A Vision for 2020”, Office for Official Publications of the European Communities, Belgium, 2001 - European Commission (2009), “REGULATION (EC) No 1070/2009”, Official Journal of the European Union, October 2009 - FAA (2012), “NextGEN Implementation Plan”, March 2012 - S. Gupta and D. Bertsimas (2011), “Multistage Air Traffic Flow Management under Capacity Uncertainty: A Robust and Adaptive Optimization Approach”, TSL workshop, Asilomar, CA, June 2011 - M. Helme (1992), “Reducing Air Traffic Delay in a Space-Time Network”, Proceedings of the IEEE International Conference on Systems, Man and Cybernetics, Chicago, 1992, pp. 236-242 - R. Hoffman (1997), “Integer programming models for Ground-Holding in Air Traffic Flow Management”, PhD Thesis, University of Maryland, 1997 - IATA (2013), IATA Economic Briefing “Inefficiency in European Airspace”, December 2013 - ICAO 4444, “Air Traffic Management – PANS-ATM”, Doc 4444, Fifteenth Edition, 2007 - ICAO 9854, “Global Air Traffic Management Operational Concept”, Doc 9854, First Edition, 2005 - S. Jain and P. Tsiotras (2008), “Multiresolution-based direct trajectory optimization”, Journal of Guidance, Control, and Dynamics, Vol. 31, No. 5, 2008, pp. 1424–1436 - R. Jovanovic (2011), “Economic measures for air traffic demand management”, PhD Thesis, Faculty of Transport and Traffic Engineering, University of Belgrade, Serbia, July 2011 103 - N. Kantas, A. Lecchini-Visintini and J. M. Maciejowski (2010), “Simulation- based Bayesian optimal design of aircraft trajectories for air traffic management”, International Journal of Adaptive Control and Signal Processing, Vol. 24, August 2010 pp. 882–899 - A. Katok and A. Sossinsky (2006), “Introduction to Modern Topology and Geometry”, 2006 - S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi (1983) - “Optimization by Simulated Annealing”, Science, New Series, Vol. 220, No. 4598, May 1983, pp. 671- 680 - J. K. Kuchar and L. C. Yang (2000), “A Review of Conflict Detection and Resolution Modeling Methods”, IEEE Transactions n intelligent transportation systems, Vol. 1, No. 4, December 2000, pp. 179-189 - S. Lan (2003), “Planning for robust airline operation: Optimize aircraft routings and flight departure times to achieve minimum passenger disruptions”, PhD Thesis, Massachusetts Institute of Technology, June 2003 - S. LaValle (2006), “Planning algorithms”, Cambridge University Press, 2006 - M. Leroux (2000), “Cognitive aspects and automation”, in N. B. Sarter and R. Amalberti (Eds.), “Cognitive engineering in the aviation domain”, Lawrence Erlbaum Associates, 2000, pp. 99–133 - K. Lindsay, E. Boyd, R. Burlingame (1993), “Traffic Flow Management Modeling with the Time Assignment Model”, Air Traffic Control Quarterly, Vol. 1, No. 3, 1993, pp. 255-276 - G. Lulli and A. Odoni (2007), “The European Air Traffic Flow Management Problem”, Transportation Science, Vol. 41, No. 4, November 2007, pp. 431–443 - P. Leal de Matos, B. Chen, R. J. Ormerod (2001), “Optimisation Models for Re- Routing Air Traffic Flows in Europe”, The Journal of the Operational Research Society, Vol. 52, No. 12, December 2001, pp. 1338-1349 - P. Leal de Matos and R. J. Ormerod (2000), “The application of operational research to European air traffic flow management – understanding the context”, European Journal of the Operational Research, Vol. 123, 2000, pp. 125-144 - N. Metropolis, A. Resenbluth, M. Resenbluth, and A Teller (1953), “Equation of chemical physics”, Journal of Mathematics Computing and Modeling, Vol. 21, 1953, pp. 1087–1092 104 - A. Mukherjee (2004), “Dynamic Stochastic Optimization models for Air Traffic Flow Management”, PhD Thesis, University of California at Berkeley, Fall 2004 - A. Mukherjee and M. Hansen (2007), “A Dynamic Stochastic Model for the Single Airport Ground Holding Problem”, Transportation Science, Vol. 41, No. 4, November 2007, pp. 444–456 - A. Neal, J. Flach, M. Mooij, S. Lehmann, S. Stankovic and S. Hasenbch (2011), “Envisaging the Future Air Traffic Management System”, The International Journal of Aviation Psychology, Vol. 21, No. 1, January 2011, pp. 16-34 - F. Netjasov, A. Vidosavljevic, V. Tosic and H. Blom, “Systematic Validation of a Mathematical Model of ACAS Operations for Safety Assessment Purposes”, Proceedings of the 9th FAA/Europe Air Traffic Management Research and Development Seminar, Berlin, Germany, 2011, pp. 12 - F. Netjasov, A. Vidosavljevic, V. Tosic, M. Everdij and H. Blom, “Development, validation and application of stochastically and dynamically coloured Petri net model of ACAS operations for safety assessment purposes”, Transportation Research Part C, vol. 33, 2013, pp. 167-195. (ISSN: 0968-090X) - S. Oussedik and D. Delahaye (1998), “Reduction of Air Traffic Congestion by Genetic Algorithms”, In Proceedings of the Fifth International Conference on Parallel Problems Solving from Nature, 1998, pp. 885-864 - PRC, Performance Review Commission - EUROCONTROL (2005), “Performance Review Report – PRR 2004”, May 2005 - PRC, Performance Review Commission - EUROCONTROL (2009), “Performance Review Report Executive Summary – PRR 2008”, May 2009 - PRC, Performance Review Commission - EUROCONTROL (2010), “Performance Review Report Executive Summary – PRR 2009”, May 2010 - PRC, Performance Review Commission - EUROCONTROL (2011), “Performance Review Report Executive Summary – PRR 2010”, May 2011 - PRC, Performance Review Commission - EUROCONTROL (2012), “Performance Review Report Executive Summary – PRR 2011”, May 2012 105 - PRC, Performance Review Commission - EUROCONTROL (2013), “Performance Review Report Executive Summary – PRR 2012”, May 2013 - M. Prandini, J. Lygeros, A. Nilim and S. Sastry (1999), “A Probabilistic framework for aircraft conflict detection”, American Institute of Aeronautics and Astronautics, AIAA-99-4144, 1999 - O. Richard , S. Constans and R. Fondacci (2011), “Computing 4D near optimal trajectories for dynamic air traffic flow management with column generation and branchand-price”, Transportation Planning and Technology, Vol. 34, No. 5, June 2011, pp. 389-411 - O. Richetta (1991), “Ground holding strategies for Air Traffic Control under uncertainty”, PhD Thesis, Massachusetts Institute of Technology, Cambridge, May 1991 - O. Richetta, A. Odoni (1994), “Dynamic solution of the Ground-holding problem in Air Traffic Control”, Transportation Research, Vol. 28, No. 3, 1994, pp. 167-185 - G. P. Roussos, G. Chaloulos, K. J. Kyriakopoulos, J. Lygeros (2008), “Control of multiple non-holonomic air vehicles under wind uncertainty using Model Predictive Control and decentralized navigation functions”, 47th IEEE Conference on Decision and Control, Cancun, Mexico, December 2008, pp 1225-1230 - G. P. Roussos, D. V. Dimarogonas, K. J. Kyriakopoulos (2009), “3D Navigation and Collision Avoidance for nonholonomic aircraft-like vehicles”, International Journal of Adaptive Control and Signal Processing, 2009, pp 1-21 - T.W. Sederberg (2009), “Computer Aided Geometric Design” Course Notes, Brigham Young University, 2009 - SESAR Joint Undertaking (2013), “SESAR Concept of Operations Step 2”, Edition 1, July 2013 - J.A. Sethian (1999), “Fast marching methods”, SIAM review, Vol. 41, No. 2, 1999, pp. 199–235 - J.A. Sethian and A. Vladimirsky (2003), “Ordered upwind methods for static Hamilton-Jacobi equations: Theory and algorithms”, SIAM Journal on Numerical Analysis, Vol. 41, 2003, pp. 325–363 106 - M. Snelder, H.J. van Zuylen and L.H. Immers (2012), “A framework for robustness analysis of road networks for short term variations in supply”, Transportation Research Part A, Vol. 46, 2012, pp. 828–842 - B. Sridhar, K. S. Sheth and S. Grabbe (1998), “Airspace Complexity and its Application in Air Traffic Management”, 2nd USA/Europe ATM R&D Seminar, Orlando, USA, December 1998 - STATFOR, the EUROCONTROL Statistics and Forecast Service (2013), “Challenges of Growth 2013 - European Air Traffic in 2050”, June 2013 - W.M. Swan and N. Adle (2006), “Aircraft trip cost parameters: a function of stage length and seat capacity”, Transportation Research part E: Logistics and Transportation Review, Vol. 42, No. 2, pp. 105-115, 2006 - M. Terrab (1990), “Ground holding strategies for Air Traffic Control”, PhD Thesis, Massachusetts Institute of Technology, Cambridge, February 1990 - M. Terrab, A. Odoni (1993), “Strategic Flow Management for Air Traffic Control”, Operations Research, Vol. 41, No. 1, Special Issue on Stochastic and Dynamic Models in Transportation, January-February 1993, pp. 138-152 - M. Terrab, S. Paulose (1993), “Dynamic Strategic and Tactical Air Traffic Flow Control”, RPI Technical Report, 1993 - V. Tosic, O. Babic (1995), “Air Route Flow Management: Problems and Research Efforts”, Transportation Planning and Technology, Vol. 19, No. 1, 1995, pp. - V. Tosic, O. Babic, M. Cangalovic and D. Hohlacov (1995), “Some models and algorithms for en route Air Traffic Flow Management”, Transportation Planning and Technology, Vol. 19, No. 2, 1995, pp. 147-164 - V. Tosic, O. Babic, M. Cangalovic and D. Hohlacov (1997), “A model to solve En route Air Traffic Flow Management Problem: A temporal and spatial case”, Proceedings of the First Air Traffic Management Research and Development Seminar, Saclay, France, June 1997 - P. Vranas (1996) – “Optimal slot allocation for European air traffic flow management”, Air Traffic Control Quarterly, 1996 - P. Vranas, D. Bertsimas, A. Odoni (1992), “The Multi-Airport Ground- Holding Problem in Air Traffic Control”, Internal report OR 263-92, February 1992 107 - P. Vranas, D. Bertsimas, A. Odoni (1994a), “The Multi-Airport Ground- Holding Problem in Air Traffic Control”, Operations Research, Vol. 42, No. 2, March-April 1994, pp. 249-261 - P. Vranas, D. Bertsimas, A. Odoni (1994b), “Dynamic Ground-Holding Policies for a Network of Airports”, Transportation Science, Vol. 28, No. 4, November 1994, pp. 275-291 - R. Weber, JL. Garcia, G. Gawinowski, M. Brochard and S. Carotenuto (2007), “ERASMUS: Concept of Operations & Initial Modeling results”, In Proceedings of the 7th USA/Europe ATM R&D Seminar, Barcelona, Spain, July 2007 - Westminster (2011), “European airline delay cost reference values”, Final Report (Version 3.2), Department of Transport Studies, University of Westminster, London, March 2011 - A. Wieland and C. M. Wallenburg (2012), “Dealing with supply chain risks - Linking risk management practices and strategies to performance”, International Journal of Physical Distribution & Logistics Management, Vol. 42, No. 10, February 2012, pp. 887-905 - Z. Yiming (2012), “Efficient and robust aircraft landing trajectory optimization”, PhD Thesis, Georgia Institute of Technology, May 2012 - K. Zeghal (1998), “A review of different approaches based on force field for airborne conflict resolution”, American Institute of Aeronautics and Astronautics, AIAA-98-4240, 1998 - S. Zenios (1991), “Network based models for air-traffic control”, European Journal of Operational Research, Vol. 50, 1991, pp. 166-178 108 BIOGRAPHY Andrija Vidosavljević was born in Boljevac, Serbia, 08.11.1982, where he finished primary and secondary school. He graduated from Faculty of Transport and Traffic Engineering in 2007 and enrolled in the PhD program at the Department of Air Transport and Traffic in the following school year. After passing all required exams, the research topic for his PhD was submitted in June and accepted by the University of Belgrade in September 2011. He has been working as research assistant at Division of Airports and Air Traffic Safety since September 2007. He received a grant of the Ministry of Science and Technological Development, Republic of Serbia in 2011. Until December 2013 he has been employed at the Faculty of Transport and Traffic Engineering. During more than six years at the Faculty he has been engaged in labs for following undergraduate courses at the Department of Air Transport and Traffic: Air Traffic Control, Air Navigation and Communication Aids and Meteorology. He has interned in AENA, Spain (2007) in Air Navigation System Development Department, and at ENAC in the Laboratory of Applied Math (twice 2011 and 2013); and attended several certified aviation courses, seminars and workshops. He has been involved in many international and national research and professional projects.