Optimal trajectory management for aircraft descent operations subject to time constraints
Abstract
The growth in traffic increased the pressure on the environmental sustainability of air transport. In this context, many research effort has been devoted to minimise the environmental impact in the different phases of flight. Continuous descent operations, ideally performed with the engines at idle from the cruise altitude to right before landing, have shown to reduce fuel, noise nuisance and gaseous emissions if compared to conventional descents. However, this type of operations suffer from a well known drawback: the loss of predictability from the air traffic control point of view in terms of overfly times at the different waypoints of the route. Due to this loss of predictability, air traffic controllers require large separation buffers, thus reducing the capacity of the airport. Previous works investigating this issue showed that the ability to meet a controlled time of arrival at a metering fix could enable continuous descent operations while simultaneously maintaining airport throughput. In this context, the planning and guidance functions of state-of-the-art flight management systems need to be modernised. On-board trajectory planners capable to generate an optimal trajectory plan satisfying time constraints introduced during the flight are seldom, mainly because the real-time optimisation of aircraft trajectories is still elusive. Furthermore, the time scale and spatial resolution of the wind forecasts used by these trajectory planners are far from being adequate to generate accurate flight time predictions. Finally, there exist guidance strategies capable to accurately comply with time constraints enforced at a certain fix in the trajectory plan, yet they are not specifically designed to minimise the environmental impact. This PhD thesis aims at investigating fast optimisation techniques to enable real-time updates of the optimal trajectory plan subject to time constraints during the descent; wind networking concepts to generate more accurate and up-to-date wind forecasts and, consequently, time predictions; and more robust an efficient guidance strategies to reduce the environmental impact at the maximum extent while complying with the time constraints of the trajectory plan. First, the feasible time window at a metering fix that could be achieved during a descent requiring neither thrust nor speed brakes usage is quantified as a function of the aircraft states (altitude, distance to the metering fix and airspeed), aiming to assess the feasibility of guidance strategies that take advantage of time and energy management concepts. Then, the performance of four of these guidance strategies is compared in terms of environmental impact mitigation and ability to satisfy operational constraints. Results from the comparison reveal that model predictive control, a strategy based on a frequent re-calculation of the optimal trajectory plan during the execution of the descent, is the most robust in terms of energy and time deviation at the metering fix, providing at the same time excellent environmental impact mitigation figures. However, the execution time required to solve a rigorous trajectory optimisation problem at each re-calculation instant remains a critical limitation for practical applications. In order to address this issue, a variant of the model predictive control strategy that allows for fast updates of the optimal trajectory plan based on parametric sensitivities is proposed, which shows analogous results yet halving the time needed to update the optimal trajectory plan. Finally, the potential benefits of using wind observations broadcast by nearby aircraft to reconstruct the wind profile downstream right before updating the optimal trajectory plan when using model predictive control is also investigated. Promising results show that the combination of model predictive control with wind networking concepts could enable optimal descents without degrading the capacity of the airport.
Full text
Optimal trajectory management for aircraft descent operations subject to time constraints RAMON DALMAU CODINA Aeronautical Engineer Advisor DR. XAVIER PRATS I MEN ´ ENDEZ Doctorate program in Aerospace Science and Technology Department of Physics – Aeronautics Division Technical University of Catalonia – BarcelonaTech A dissertation submitted for the degree of International Doctor of Philosophy May 2019
Optimal trajectory management for aircraft descent operations subject to time constraints Author Ramon Dalmau Codina Advisor Dr. Xavier Prats i Men´ endez Reviewers Dr. Daniel Delahaye Dr. Erik-Jan van Kampen Thesis committee Dr. Daniel Delahaye Dr. Erik-Jan van Kampen Dr. Manuel Soler Doctorate program in Aerospace Science and Technology Technical University of Catalonia – BarcelonaTech May 2019 This dissertation is available on-line at the Theses and Dissertations On-line (TDX) repository, which is managed by the Consortium of University Libraries of Catalonia (CBUC) and the Supercomputing Centre of Catalonia (CESCA), and sponsored by the Generalitat (government) of Catalonia. The TDX repository is a member of the Networked Digital Library of Theses and Dissertations (NDLTD) which is an international organisation dedicated to promoting the adoption, creation, use, dissemination and preservation of electronic analogues to the traditional paper-based theses and dissertations http://www.tdx.cat This is an electronic version of the original document and has been re-edited in order to fit an A4 paper. PhD. Thesis made in: Department of Physics – Aeronautics Division Esteve Terradas, 5. 08860 Castelldefels Catalonia (Spain) This work is licensed under the Creative Commons Attribution-Non-commercial-No Derivative Work 3.0 Spain License. To view a copy of this license, visit http://creativecommons.org/ licenses/by-nc-nd/3.0/es/deed.en_GB or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco, California, 94105, USA.
Als meus pares i al meu germà Olga, Ramon i Ferran
Contents List of Figures ................................................. vii List of Tables ................................................. ix List of Publications .............................................. xi Agra¨ ıments .................................................. xiii Resum ..................................................... xv Abstract .................................................... xvii Notation .................................................... xix List of Acronyms ............................................... xxv CHAPTER I Introduction ......................................... 1 I.1 Environmental impact of descents ................................. 2 I.2 Predictability of descents ...................................... 4 I.3 New concepts to enable CDOs in busy TMAs .......................... 6 I.4 Motivation of this PhD ........................................ 9 I.5 Objectives of this PhD thesis .................................... 12 I.6 Scope and limitations of this PhD thesis ............................. 14 I.7 Outline of this PhD thesis ...................................... 15 CHAPTER II Framework on trajectory management ......................... 17 II.1 Models needed for trajectory management ............................ 18 II.2 Optimal trajectory planning ..................................... 22 II.3 Trajectory guidance ......................................... 34 CHAPTER III The energy-neutral time window ............................ 41 III.1 State of the art ............................................ 42 III.2 Setup of the experiment ...................................... 43 III.3 Results ................................................ 49 III.4 Discussion .............................................. 56 v
CHAPTER IV Comparison of guidance strategies to meet time constraints in optimal descents ............................................. 57 IV.1 State of the art ............................................ 58 IV.2 Setup of the experiment ...................................... 59 IV.3 Results ................................................ 66 IV.4 Discussion .............................................. 73 CHAPTER V Model predictive control to meet time constraints in optimal descents ..... 75 V.1 State of the art ............................................ 76 V.2 Setup of the experiment ...................................... 76 V.3 Results ................................................ 83 V.4 Discussion .............................................. 94 CHAPTER VI Wind-networking to improve time predictability and fuel efficiency of descents 97 VI.1 State of the art ............................................ 98 VI.2 Setup of the experiment ...................................... 98 VI.3 Results ................................................100 VI.4 Discussion ..............................................107 CHAPTER VII Concluding remarks ....................................109 VII.1 Summary of contributions .....................................109 VII.2 Future Research ...........................................111 vi
List of Figures I-1 Comparison between a CDO and a current descent operation . . . . . . . . . . . . . . . . . . 2 I-2 Comparison between two CDOs with different mass and CI . . . . . . . . . . . . . . . . . . 4 I-3 Vectoring instructions used by the ATC at Barcelona-El Prat airport . . . . . . . . . . . . . . 5 I-4 Scheme of the three-degrees deceleration approach concept . . . . . . . . . . . . . . . . . . . 7 I-5 Scheme of the 3D-path arrival management concept . . . . . . . . . . . . . . . . . . . . . . . 8 II-1 Theoretical and empirical weather models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 II-2 Methods to solve optimal control problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 II-3 Discretisation process for a generic optimal control problem . . . . . . . . . . . . . . . . . . 27 II-4 Comparison between interpolating and smoothing spline approximation for the maximum thrustofanA320 ........................................... 31 II-5 Jacobian matrix for a generic NLP optimisation problem . . . . . . . . . . . . . . . . . . . . 33 II-6 Tactical and strategic behaviour in a 180 degree turn . . . . . . . . . . . . . . . . . . . . . . . 35 III-1 Scenario to investigate the feasible time window of energy-neutral descents . . . . . . . . . 43 III-2 Smoothing spline approximating aircraft performance data obtained from PEP . . . . . . . 45 III-3 Definition of the phases and constraints for an energy-neutral descent . . . . . . . . . . . . 47 III-4 Examples of earliest and latest energy-neutral trajectories . . . . . . . . . . . . . . . . . . . . 50 III-5 Baseline scenario: earliest and latest times of arrival at the metering fix . . . . . . . . . . . . 51 III-6 Baseline scenario: feasible time window at the metering fix . . . . . . . . . . . . . . . . . . . 51 III-7 Baseline scenario: fuel consumption at the metering fix . . . . . . . . . . . . . . . . . . . . . 52 III-8 Baseline scenario: extra fuel consumption with respect to the minimum fuel . . . . . . . . . 53 III-9 Effects of aircraft mass on the feasible time window at the metering fix . . . . . . . . . . . . 53 III-10 Effects of longitudinal wind on the feasible time window at the metering fix . . . . . . . . . 54 III-11 Effects of the position of the metering fix on the feasible time window . . . . . . . . . . . . . 55 III-12 Effects of the initial speed on the feasible time window at the metering fix . . . . . . . . . . 56 IV-1 Setup of the experiment to compare guidance strategies . . . . . . . . . . . . . . . . . . . . . 59 IV-2 Routes at Barcelona-El Prat chosen to compare guidance strategies . . . . . . . . . . . . . . 60 IV-3 Drag polars approximating aerodynamic data for the A320 . . . . . . . . . . . . . . . . . . . 61 IV-4 Polynomials approximating propulsive data for the A320 . . . . . . . . . . . . . . . . . . . . 62 IV-5 Definition of the phases and constraints for the comparison of guidance strategies . . . . . 64 IV-6 Specific energy and time deviations per guidance strategy . . . . . . . . . . . . . . . . . . . 67 IV-7 Planned and executed states per guidance strategy . . . . . . . . . . . . . . . . . . . . . . . . 68 IV-8 Planned and executed controls per guidance strategy . . . . . . . . . . . . . . . . . . . . . . 69 IV-9 Mean metrics per guidance strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 IV-10 Metrics distribution per guidance strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 vii
IV-11 Correlation between metrics and longitudinal wind speed error . . . . . . . . . . . . . . . . 72 V-1 BOSSS TWO standard arrival procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 V-2 BADA 4 propulsive model for the A320-214 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 V-3 BADA 4 aerodynamic model for the A320-214 . . . . . . . . . . . . . . . . . . . . . . . . . . 79 V-4 Definition of the phases and constraints for the assessment of NMPC strategies . . . . . . . 80 V-5 RAP wind forecast RMSE distribution from June 2017 to June 2018 . . . . . . . . . . . . . . 81 V-6 Setup of the experiment to compare NMPC guidance strategies . . . . . . . . . . . . . . . . 82 V-7 RAPwindforecastandanalysis................................... 83 V-8 Optimal trajectory updates for perturbations in the parameters vector . . . . . . . . . . . . . 84 V-9 Planned and executed trajectories by guidance strategy . . . . . . . . . . . . . . . . . . . . . 86 V-10 Time deviation at the metering fix for different NMPC variants . . . . . . . . . . . . . . . . . 87 V-11 Specific energy deviation at the metering fix for different NMPC variants . . . . . . . . . . . 88 V-12 Additional specific energy difference with respect to the open-loop execution for the differentNMPCvariants.......................................... 90 V-13 Fuel consumption difference with respect to the open-loop execution for the different NMPCvariants ............................................ 91 V-14 Normalised execution time for the NLP and QP algorithms . . . . . . . . . . . . . . . . . . . 92 V-15 Time and specific energy error at the metering fix as a function of the mean longitudinal winderror ............................................... 93 V-16 Specific energy added/removed with respect to the initial plan as a function of the mean longitudinal wind error for the INMPC and OL . . . . . . . . . . . . . . . . . . . . . . . . . . 93 V-17 Fuel consumption difference with respect to initial plan as a function of the mean longitudinalwinderror............................................ 93 VI-1 Poissondistribution.......................................... 99 VI-2 Planned and executed trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 VI-3 RAPwindforecastandanalysis................................... 102 VI-4 Time deviation at the metering fix for different values fo µ................... 103 VI-5 Specific energy deviation at the metering fix for different values of µ............. 104 VI-6 Specific energy added by means of Twith respect to the initial plan for different values of µ105 VI-7 Fuel consumption difference with respect to the initial plan for different values of µ.... 106 viii
List of Tables I-1 Altitude and overfly time of two CDOs with different mass and CI at various distances to go 4 III-1 Definition of the phases and constraints for an energy-neutral descent . . . . . . . . . . . . 46 III-2 List of case studies for the energy-neutral feasible time window assessment . . . . . . . . . 48 IV-1 Definition of the phases and constraints for the comparison of guidance strategies . . . . . 63 IV-2 List of case studies per each of the four guidance strategies assessed . . . . . . . . . . . . . . 65 V-1 Definition of the phases and constraints for the comparison of NMPC strategies . . . . . . . 80 V-2 Metrics of energy removed by using speed brakes for the different NMPC variants and windforecastlookaheadtimes ................................... 88 V-3 Percentage of energy neutral trajectories for the different NMPC variants and forecast lookaheadtimes .............................................. 90 VI-1 Metrics of energy removed by using speed brakes for different values of µand wind forecastlookaheadtimes ......................................... 104 VI-2 Percentage of energy neutral trajectories for different values of µand forecast look-ahead times .................................................. 105 ix
l’exactitud de la traject` oria calculada. Els resultats d’aquest ´ ultim experiment demostren que si s’implement´ es model predictiu com a estrat` egia de guiatge i les aeronaus cooperessin per compartir observacions de vent, es reduiria l’impacte ambiental sense degradar la capacitat de l’aeroport. xvi
Abstract The growth in traffic increased the pressure on the environmental sustainability of air transport. In this context, many research effort has been devoted to minimise the environmental impact in the different phases of flight. Continuous descent operations, ideally performed with the engines at idle from the cruise altitude to right before landing, have shown to reduce fuel, noise nuisance and gaseous emissions if compared to conventional descents. However, this type of operations suffer from a well known drawback: the loss of predictability from the air traffic control point of view in terms of overfly times at the different waypoints of the route. Due to this loss of predictability, air traffic controllers require large separation buffers, thus reducing the capacity of the airport. Previous works investigating this issue showed that the ability to meet a controlled time of arrival at a metering fix could enable continuous descent operations while simultaneously maintaining airport throughput. In this context, the planning and guidance functions of state-of-the-art flight management systems need to be modernised. On-board trajectory planners capable to generate an optimal trajectory plan satisfying time constraints introduced during the flight are seldom, mainly because the real-time optimisation of aircraft trajectories is still elusive. Furthermore, the time scale and spatial resolution of the wind forecasts used by these trajectory planners are far from being adequate to generate accurate flight time predictions. Finally, there exist guidance strategies capable to accurately comply with time constraints enforced at a certain fix in the trajectory plan, yet they are not specifically designed to minimise the environmental impact. This PhD thesis aims at investigating fast optimisation techniques to enable real-time updates of the optimal trajectory plan subject to time constraints during the descent; wind networking concepts to generate more accurate and up-to-date wind forecasts and, consequently, time predictions; and more robust an efficient guidance strategies to reduce the environmental impact at the maximum extent while complying with the time constraints of the trajectory plan. First, the feasible time window at a metering fix that could be achieved during a descent requiring neither thrust nor speed brakes usage is quantified as a function of the aircraft states (altitude, distance to the metering fix and airspeed), aiming to assess the feasibility of guidance strategies that take advantage of time and energy management concepts. Then, the performance of four of these guidance strategies is compared in terms of environmental impact mitigation and ability to satisfy operational constraints. Results from the comparison reveal that model predictive control, a strategy based on a frequent re-calculation of the optimal trajectory plan during the execution of the descent, is the most robust in terms of energy and time deviation at the metering fix, providing at the same time excellent environmental impact mitigation figures. However, the execution time required to solve a rigorous trajectory optimisation problem at each re-calculation instant remains a critical limitation for practical applications. In order to address this issue, a variant of the model predictive control strategy that allows for fast updates of the optimal trajectory plan based on parametric sensitivities is proposed, which shows analogous results yet halving the time needed to update the optimal trajectory plan. Finally, the potential benefits of using wind observations broadcast xvii
by nearby aircraft to reconstruct the wind profile downstream right before updating the optimal trajectory plan when using model predictive control is also investigated. Promising results show that the combination of model predictive control with wind networking concepts could enable optimal descents without degrading the capacity of the airport. xviii
Notation Throughout this PhD thesis and as a general rule, scalars and vectors are denoted either with lower or upper case letters. Vectors are noted with the conventional overhead arrow, like for example ~a or ~ ψ. Sets are denoted using caligraphic fonts, like for example A,Bor X, while matrices use the same font but in bold series, like R. The derivative of magnitude a(t)with respect to tis expressed by da(t) dt= ˙a. Furthermore, if not otherwise noted, all vectors are column vectors and a transposed vector is denoted by [·]T. For an ordered set or a vector, for example A, the element at the ith position is denoted by Ai. Similarly, Ai:jis the subset of Athat includes the elements whose position ranges from ito j(inclusive). Being consistent with this notation, Ai:is the subset of Athat includes the elements with position greater or equal than i. Note that A1: =A. Furthermore, (·)∗indicates optimal. Finally, for a scalar function (·)+=max((·),0), and for a vector-valued function of nelements, (·)+=h(·)+ 1,(·)+ 2,...,(·)+ ni. Next, the principal symbols that are used throughout this dissertation are shown along with their meaning. The reader should note that this list is not exhaustive. BB-Spline basis function CDdrag coefficient CLlift coefficient CD,BADA drag coefficient of the BADA v4 model CD,sim drag coefficient of the simulator C0 D,sim drag coefficient of the simulator increased by a factor of KCD CD0parasite drag coefficient CD2induced drag coefficient CDβ coefficient that represents the increase in drag coefficient when speed brakes are deployed Daerodynamic drag Ekkinetic energy of the aircraft Eppotential energy of the aircraft Esspecific energy of the aircraft Ettotal energy of the aircraft Jcost function of the optimal control problem KCDfactor to increase or decrease the nominal drag coefficient KTidle factor to increase or decrease the nominal thrust idle Kβweighting parameter that determines howmuch the use of speed brakes is penalised Laerodynamic lift force MMach number Nnumber of discretisation intervals in the whole time horizon N1revolutions of the engine fan N1idle minimum revolutions of the engine fan xix
N1max maximum revolutions of the engine fan Pnumber of phases of the optimal control problem Rperfect gas constant Swing surface area Tengine thrust Tidle,sim idle thrust of the simulator T0 idle,sim idle thrust of the simulator increased by a factor of KTidle Tidle idle thrust Tmax maximum thrust ∆Esdifference of specific energy ∆τdiscretisation step Λforgetting factor Πquadrature function of the running cost <set of real numbers αHellman exponent βspeed brakes deflection χgground track angle δnormalised pressure of the air δ11 normalised pressure of the air at the tropopause inf infeasibility criteria opt optimality criteria ηlanding gear status γaerodynamic flight path angle γaa specific heat ratio of the air γmin minimum aerodynamic flight path angle πthrottle setting λhstandard temperature lapse rate Weq weighting matrix for the slack variables associated with equality constraints Win weighting matrix for the slack variables associated with inequality constraints Eset associating the index of the last time sample to its corresponding phase Gac active set of inequality constraints Gin inactive set of inequality constraints Iset associating the index of the first and interior time samples to its corresponding phase Llagrangian of the nonlinear programming optimisation problem Nnormal distribution Oset of wind observations Palgorithm that solves the nonlinear programming optimisation problem Qalgorithm that solves the quadratic programming optimisation problem Tset associating the index of each time sample to its corresponding phase Uuniform distribution µparameter describing the expected number of occurrences in a Poisson distribution µexpectation of the simulated wind speed error ωweight associated to each wind observation φend cost or Mayer term πrunning cost or Lagrange term ρSSL SSL density of the air at standard sea level σnormalised density of the air σstandard deviation of the simulated wind speed error τdiscretised time sample τSSL temperature of the air at standard sea level τinf threshold for the infeasibility criteria τopt threshold for the optimality criteria Pos Poisson distribution θnormalised temperature of the air θ11 normalised temperature of the air at the tropopause ~ Fevolution function of the state vector ~ Xcurrent states vector ˆ ~ Xpredicted states vector xx
~ λdual variables associated with inequality constraints ~µ dual variables associated with equality constraints ~ ψterminal constraints ~ ϑeq equality interior-point constraints ~ ϑin inequality interior-point constraints ~ beq equality path constraints ~ bin inequality path constraints ~c control points of the spline ~ dfixed parameters of the model ˆ ~ dpredicted fixed parameters of the model ~eeq vector of slack variables for the equality constraints ~ein vector of slack variables for the inequality constraints ~ fdynamics of the state vector ~g inequality constraints of the nonlinear programming optimisation ~ hequality constraints of the nonlinear programming optimisation problem ~p vector of parameters of the nonlinear programming optimisation problem ˆ ~p predicted vector of parameters of the nonlinear programming optimisation problem ~u controls vector ~w wind vector ~x states vector ~z primal variables of the nonlinear programming optimisation problem ζaerodynamic configuration of the aircraft aspeed of sound fcost function of the nonlinear programming optimisation problem ggravity acceleration hgeometric altitude hrreference altitude of the Hellman model h11 standard altitude of the tropopause hCO cross-over altitude hF AP altitude at which the FAP has to be intercepted hIF altitude at which the IF has to be intercepted kmax maximum number of QP iterations before triggering a NLP trajectory update mmass of the aircraft nzload factor pprobability pSSL pressure of the air at standard sea level qnominal fuel flow qidle idle fuel flow sdistance to go sF AP distance to go of the FAP sIAF distance to go of the IAF sIF distance to go of the IF ttime tFfinal time of the time horizon tIinitial time of the time horizon vtrue airspeed vCAS calibrated airspeed vF AP calibrated airspeed at which the FAP has to be intercepted vIF calibrated airspeed at which the IF has to be intercepted vmax maximum true airspeed vmin minimum true airspeed weeast wind component wnnorth wind component wrreference speed of the Hellman model wslongitudinal wind component wxcross wind component simulated wind speed error xxi
xxii
List of Acronyms 3D-PAM 3D-path arrival management ACARS aircraft communication and reporting system ADS-B automatic dependent surveillance-broadcast AGL above ground level AMAN arrival manager AMDAR aircraft meteorological data relay ANSP air navigation service providers APERO avionics prototyping environment for research and operations APM aircraft performance model AsNMPC advanced-step NMPC ATC air traffic control ATFM air traffic flow management ATM air traffic management BADA base of aircraft data CAS callibrated airspeed CDA-MP continuous descent approach for maximum predictability CDO continuous descent operation CI cost index CONCORDE flight operations for novel continuous descents CTA controlled time of arrival DAE differential algebraic equation EDA efficient descent advisor FAA Federal Aviation Administration FAP final approach point FASTOP fast optimizer for continuous descent approaches FMS flight management system GAMS general algebraic modelling system GD green dot speed xxiii
GFS global forecast system GRACE generic research aircraft cockpit environment IAF initial approach fix ICAO International Civil Aviation Organization IF intermediate fix ILS instrumental landing system INMPC ideal NMPC ISA international standard atmosphere KKT Karush-Kuhn-Tucker conditions MLM maximum landing mass MMO maximum operative Mach MPC model predictive control NASA National Aeronautics and Space Administration NextGen next generation air transportation system NLP nonlinear programming NLR Netherlands Aerospace Centre NMPC nonlinear model predictive control NOAA National Oceanic and Atmospheric Administration NWP numerical weather prediction ODE ordinary differential equation OL open-loop guidance OPD optimised profile descent OPTIMAL optimised procedures and techniques for improvement of approach and landing PARTNER partnership for air transportation noise and emission reduction PEP Airbus performance engineering program PM particulate matter QP quadratic programming RAP rapid refresh RMSE root-mean-square error RNAV area navigation RTA required time of arrival SbNMPC sensitivity-based NMPC SESAR Single European Sky Air traffic management Research SLQP sequential linear-quadratic programming SNOPT sparse nonlinear optimiser SQP sequential quadratic programming SSL standard sea level STAR standard arrival procedure STATFOR statistics and forecast service SUNDIALS SUite of Nonlinear and DIfferential/ALgebraic Equation Solvers TA tailored arrival TAAM total airport and airspace model TAS true airspeed TDDA three-degrees deceleration approach xxiv
TEMO time and energy managed operations TMA terminal maneouvering area TOD top of descent TRL technology readiness level UPC Technical University of Catalonia UPS United Parcel Service VMO maximum operative CAS VNAV vertical navigation WAFC world area forecast model xxv
6Chapter I - Introduction Since field tests of CDOs in high traffic scenarios are difficult to conduct, Wilson & Hafner (2005) presented a simulation-based approach to quantify the benefits of CDOs in busy TMAs taking into consideration spacing and sequencing issues. Fast-time simulations at ATL airport were performed with the total airport and airspace model (TAAM) by removing altitude constraints associated to the waypoints of the original procedures, as well as allowing aircraft to directly fly from the initial approach fix (IAF) to the base-leg of the landing runway. Promising results showed significant fuel and time savings with minimum increase on ATC workload. Aiming at assessing safety aspects of CDOs, Park et al. (2016) proposed a trajectory-based methodology to estimate the encounter rate of CDOs in high traffic scenarios. The proposed method was validated by comparing to Monte-Carlo simulation results. Khan et al. (2009) also performed fast-time simulations in high traffic scenarios to compare current ATC procedures with ATC assisted with various systems of ground automation, in terms of CDO success rate and throughput. Statistical analysis from the fast-time experiment performed by Robinson & Kamgarpour (2010) revealed that fuel savings of CDOs are marginal for high traffic scenarios due to the tactical instructions required to provide separation assurance, which have negative effects on the fuel-efficiency of CDOs. In particular, fuel savings due to CDOs were reduced by 70-85% if compared to the low traffic scenarios. This work also reported a reduction of 45 kg per CDO due to traffic congestion at ATL airport. A similar experiment at ATL airport reported a sightly different reduction of 23 to 36 kg per CDO (Cao et al.,2014). In addition to fast-time simulations, human-in-the-loop studies with ATC showed that CDOs would not be feasible at ATL airport during the busiest traffic periods, since too much efficiency would be lost to accommodate the traffic demand (Johnson,2009). In such traffic conditions, few flights would be issued the CDO clearance and even some flights would likely need to be removed from CDO procedures after being cleared in order to manage the demand. Nevertheless, results also showed that at least 15% of the total environmental benefit would still be achievable for CDOs that are terminated earlier. I.3 New concepts to enable CDOs in busy TMAs Results from previous works reveal that the predicted increase in traffic demand will negatively impact the potential benefits of CDOs due to their low predictability if compared to current descent operations. In order to face this issue, several ATM concepts have been proposed that aim to enable CDOs also in high traffic demand scenarios. I.3.1 Fixed flight path angle descents The three-degrees deceleration approach (TDDA) (De Prins et al.,2007;Sopjes et al.,2011;De Leege et al., 2009) is a concept in which the descent is performed all along a 3◦path from the TOD to the runway threshold. To perform the spacing task, two controls are given to the pilot: the first is the thrust cutback altitude, in which the throttle is set to idle and, consequently, the aircraft starts to decelerate; and the second is the flaps/slats and gear setting, which allows to manage the deceleration after the thrust cutback altitude. Figure I-4 (adapted from De Prins et al. (2007)) shows a scheme of this concept. Pradeep & Wei (2017) investigated the variability and operational feasibility aspect of CDOs with fixed path angle descents, and concluded that this kind of flight operation has higher degree of vertical predictability and lesser variability in terms of altitude and speed profiles if compared to the CDOs performed at idle thrust. This is achieved at the cost of thrust settings different from idle and with higher drag configurations, which implies more environmental impact. I.3.2 CDOs in vertical corridors Aiming to allow for more flexible vertical profiles than fixed flight path angle descents, the International Civil Aviation Organization (ICAO) published some CDO guidance material (ICAO,2010) to support air navigation service providers (ANSP) to design vertical corridors in which all descent trajectories must be contained, helping in this way to strategically separate them from other arrival and departure procedures in the vicinity. As reported by Fricke et al. (2015), however, these vertical corridors have been established without explicitly considering the aircraft type, assuming international standard atmosphere (ISA) conditions and with coarse assumptions regarding the aircraft gross mass and aircraft performance models. This leads, in the majority of cases, to too restrictive corridors that limit the CDO adherence in real operations.
I.3 New concepts to enable CDOs in busy TMAs 7 Figure I-4: Scheme of the three-degrees deceleration approach concept I.3.3 Optimised profile descents Standard arrival procedures (STARs) used in most of the airports have been historically designed to ensure that a wide variety of aircraft can fly the procedure, from those basically equipped to the most modern aircraft. In order to enable the descent to aircraft with various navigation capabilities, current STARs include altitude and speed constraints at various fixes of the route and, generally speaking, are not designed to execute the whole descent with the engines at idle. An optimised profile descent (OPD) is a procedure, normally associated with a published STAR, designed to allow maximum practical use of CDOs (ICAO,2010). Accordingly, OPD could be interpreted as one method of facilitating CDOs by using altitude constraints that restrict the vertical profile to a near-idle descent for many of the currently in use commercial aircraft. Shresta et al. (2009) assessed the potential benefits and issues of implementing OPDs at Denver international airport (DEN). Arrival flights descending along CDOs were modelled by transforming the descent trajectories of real flights, taken from radar tracks, into descent trajectories at idle thrust. Several potential conflicts were identified as a result of allowing aircraft to descend along their optimal trajectories. Implementing OPDs only in a certain altitude interval, allowing level-off above and below the delimiting altitudes, demonstrated to greatly reduce conflicts while still offering attractive benefits in terms of fuel consumption. As in previously mentioned studies (Johnson,2009), a major conclusion of the work performed by Shresta et al. (2009) was that CDOs do not need to be implemented from the cruise altitude to the stabilisation point in order to obtain reasonable reductions of fuel consumption and gaseous emissions. I.3.4 3D-Path arrival management The working principle of the solutions presented so far consists of removing level-offs to the greatest possible extent, while still supporting ATC to predict the vertical profile of aircraft. In addition, these solutions assume static routes, corresponding to the published STARs of the airport, and therefore not tailored to each individual flight. The next solutions go one step further by controlling the time at which aircraft reach one or several metering fixes along the route. For instance, using the 3D-path arrival management (3D-PAM) concept (Nagle et al.,2011;Coppenbarger et al.,2010), a ground system computes the arrival schedule of flights at a metering fix. Then, for each flight, it strategically assigns the best route to the metering fix in combination of cruise and speed advisories and (eventually) a path adjustment in the form of point-bearing-distance to meet the schedule. These instructions, which are designed to facilitate the delivery by voice using standard phraseology, are cleared prior the TOD and given to the pilot for manual entry into the FMS. This negotiation prior the TOD allows the FMS to know the remaining distance to go, compute the optimal trajectory plan complying with the cleared instructions, and uninterruptedly execute it. This is generally not feasible today due to the use of open-loop vectors. A system for the integration of this concept is the efficient descent advisor (EDA), developed by the National Aeronautics and Space Administration (NASA) (Coppenbarger et al.,2004).
8Chapter I - Introduction Figure I-5 shows a scheme of the 3D-PAM concept. Figure I-5: Scheme of the 3D-path arrival management concept It should be noted that, even if when implementing 3D-PAM some flights are actually travelling a longer distance due to the point-bearing-distance instruction (see blue segment in Fig. I-5), the descent to the metering fix is typically flown at a lower power setting, overall reducing the fuel consumption, gaseous emissions and noise nuisance. I.3.5 Tailored arrivals Tailored arrivals (TAs) are a step further of 3D-PAM to enable CDOs. TAs are similar concept to OPDs, but are not restricted to published STARs (Coppenbarger et al.,2009;Elmer,2008). In the TA concept, a customised, dynamic and conflict-free descent trajectory, tailored by ATC through use of a list of fixes with associated speed and altitude constraints, is uplinked via data-link to the FMS before starting the descent, given as a single clearance well before to TOD. These clearances provide sufficient information for the FMS to compute the trajectory plan from the current aircraft position to the runway threshold. When possible, the clearances should allow to generate a near-idle descent trajectory. It should be noted that even for ATC without data-link support, TAs could be implemented by using a pre-negotiated set of STARs stored in the database of the FMS. I.3.6 CDOs with controlled times of arrival and fixed routes A primary challenge of 3D-PAMs and TAs is the design of the route and constraints that are issued to the FMS, which shall be chosen to allow the FMS to generate a trajectory plan that can be flown at idle thrust. Ground systems in charge of computing these constraints, however, use trajectory predictors configured with simplified models to represent the dynamics and performance of the aircraft. These models will rarely reproduce the behaviour of all aircraft types without errors. Any discrepancy between the models used by the on-board and ground trajectory predictors might result in clearances selected by the 3D-PAM or TA ground automation system that, even if supposed to enable an idle descent according to the (inaccurate) ground trajectory predictor, might require the use of additional thrust or drag according to the FMS.
I.4 Motivation of this PhD 9 Another promising concept to enable CDOs in busy TMAs consists of enforcing controlled time of arrival (CTA) at one or several metering fixes while keeping the aircraft on a lateral route cleared before the TOD (Muynch et al.,2012). Note that this route could be either a published STAR or a sequence of waypoints uplinked to the FMS by the ground system before the TOD. With this type of flight operations, ATC would assign (using data-link, for instance) at least one CTA to each aircraft, ideally before they reach their TOD, to sequence and merge arrival traffic. A negotiation process of CTAs via data-link between ATC and FMSs provided with 4D capabilities was proposed by Uebbing-Rumke & Temme (2011) and Oberheid et al. (2008). The negotiation process includes the on-board computation (i.e., using accurate models bundled into the FMS) of the earliest and latest achievable CTA at the metering fix (in other words, the feasible time window) for the known lateral route that is to be flown, and the subsequent down-link of this information to the ground automation system. Based on the feasible time window and the surrounding traffic, the ground system replies with the CTA. The received CTA is entered as a required time of arrival (RTA) into the FMS, and then the on-board trajectory planner computes a new (optimal) trajectory plan starting at the current position, while satisfying the RTA and other operational constraints. Then, the resulting trajectory plan is executed by an on-board guidance system with 4D capabilities (De Prins et al.,2007). Note that, different from TAs and 3D-PAM, the constraints are not in form of altitude and speed instructions at the different waypoints, but purely CTAs. The on-board trajectory planner is who selects the optimal altitude and speed profile such that the applicable CTAs are satisfied. Accordingly, this kind of flight operations are expected to increase the likelihood to generate a trajectory plan at idle thrust that satisfies the CTA requested by ATC, supposed to safely sequence and merge arrival traffic. Summing up, vertical dispersions of the CDOs trajectories could be mitigated by appropriately selecting the location of the metering fix and conveniently assigning CTAs to each inbound aircraft, while allowing them to follow a pre-negotiated lateral route to completion without ATC interventions. This kind of flight operations, however, require a modernisation of both current ground automation and airborne systems. The focus of this PhD will be in the development of new planning and guidance strategies to accurately and efficiently satisfy CTAs during CDOs. I.4 Motivation of this PhD Sections I.4.1 and I.4.2 present the capabilities of current ground automation and on-board trajectory management systems, respectively, and the future needs for the efficient implementation of CDOs subject to CTAs. These future needs drive the motivation of this PhD thesis. I.4.1 Ground automation Arrival managers (AMANs) were developed in the late 1990s to provide automated sequencing support for ATC to achieve a smooth, predictable and optimised arrival stream with as few tactical interventions as possible. State-of-the-art AMANs use information from the flight plans, radar tracks, weather data, local airspace structure information, and approximated aircraft performance models to predict trajectories and provide the predicted time of arrival at the different fixes for every flight. Research is currently ongoing into how modern AMANs could also provide support to select appropriate CTAs (Muynch et al.,2012). ATC would be ultimately in charge to provide the CTA to the pilot, who would then use a FMS trajectory management with 4D capabilities to satisfy it. I.4.2 On-board trajectory management FMSs were developed in the early 1980s to reduce both workload of the flight crew and flight costs due to their capability to generate trajectory plans minimising fuel and time, as well as to automatically guide the aircraft along that plan by using elevator and throttle commands. In the last decades, a variety of FMSs have been designed to comply with CTAs. Most of these systems, however, have difficulties to satisfy time constraints with an accuracy of 30 s. It should be noted that even this 30 s accuracy is too large for reliable use within the TMA. In 2009, Lenz & Korn (2009) reported that only about 11% of the aircraft flying in the European airspace were equipped with FMS capable of meeting a CTA with an accuracy of 6 s.
10 Chapter I - Introduction Modern FMSs include vertical navigation (VNAV) functions that can compute the optimal trajectory plan for the descent and provide guidance to execute it. Limitations in both planning and guidance capabilities of VNAV, however, have prevented widespread adoption of time-constrained CDOs for operational use (Oseguera-Lohr et al.,2007). For instance, despite modern FMSs have the ability to meet a CTA on cruise, state-of-the-art VNAV functions do not allow pilots to easily perform speed changes during the descent while maintaining the engines at idle. I.4.2.1 Future needs for on-board trajectory management An example of concept that overcomes the main issues of current VNAV functionalities is the continuous descent approach for maximum predictability (CDA-MP), developed by Boeing (Garrido-L´ opez et al.,2009). This strategy uses a tactical controller to nullify time deviations with the elevator by following a ground speed plan, and operates with the throttle applying small amounts of thrust in order to provide a mean to reduce the energy in the face of disturbances. Since the elevator and the throttle are following, respectively, the ground speed and thrust plans, when the models used by the FMS to generate the trajectory plan do not match the reality (e.g., if the wind forecast is not representative of the actual wind conditions), the aircraft will deviate from the altitude plan. If the aircraft deviates above the plan by more than a predefined threshold, the throttle will revert to idle. Conversely, if the aircraft deviates below the plan by more than a predefined threshold, the thrust will increase to nullify the vertical deviation. This strategy is able to satisfy CTAs while being more fuel-efficient that current VNAV functionalities. If the models used by the FMS were perfect, however, the resulting fuel consumption would be higher than that of a CDO performed at idle thrust. In fact, taking into consideration the definition of CDO stated early in this chapter, the trajectory plan computed with the CDA-MP would not be considered a CDO. The CDA-MP concept assumes that the trajectory plan is computed right after receiving the CTA prior the TOD and that is frozen afterwards, meaning that it is not dynamically updated during the course of the descent. The initial trajectory plan, however, will be actually optimal and its time prediction will be correct only if the models used by the FMS to represent the aircraft dynamics, aircraft performance, and weather are accurate enough. In practise, these models will rarely reproduce the real world and are not exempted from errors. In addition, if the deviation from the initial trajectory plan is significant, it might be preferable (in terms of environmental impact) to update the trajectory plan, rather than trying to strictly follow it. Last but not least, in the future trajectory-based ATM paradigm we could envisage ATC updating the CTA once the descent has been initiated, requiring a new trajectory plan after the TOD. Based on the above discussion, future FMS shall be able to dynamically update the trajectory plan (on-board and in real-time) at any time during the course of the descent, in such a way constraints are satisfied requiring the minimum amount of additional thrust and/or speed-brakes use. For instance, Liu et al. (2018) presented a trajectory optimisation algorithm that allows to update the trajectory during the course of a CDO. Similarly, Ramasamy et al. (2013) proposed a conceptual design for the next generation of FMSs, which includes 4D trajectory planning and air-to-ground trajectory negotiation functionalities. In this novel prototype of FMS, the optimal trajectory satisfying the negotiated CTAs is computed using direct collocation methods, updating the trajectory with an execution time around 120 seconds, which might not be acceptable for real-time applications considering the typical duration of a descent. Time and energy management concepts have been also investigated in this context, showing very promising results. I.4.2.2 The time and energy management operations concept Time and energy managed operations (TEMO) is a trajectory management (planning and guidance) concept that aims to satisfy a CTA at a metering fix by means of only elevator control while keeping, as far as possible, the engines at idle during the whole descent (De Jong,2014;De Jong et al.,2015). Otherwise, optimal amounts of energy are removed or added by means of thrust or speed brakes, respectively, to satisfy the CTA and applicable operational constraints while minimising the environmental impact. A descent in which the total energy is not increased by means of thrust higher than idle nor removed by using active drag devices is known as energy-neutral. TEMO is a new concept that was developed within the management of trajectory and mission work package of the area of systems for green operations of the Clean Sky European Joint Undertaking research initiative. TEMO is in line with SESAR step 2 capabilities, since it proposes 4D trajectory management and aims to provide significant environmental benefits in the arrival phase without negatively affecting throughput, even in high density and peak-hour operations.
I.4 Motivation of this PhD 11 From an ATC point of view, the TEMO concept assumes that the AMAN will use available trajectory information to determine the preferred arrival route, landing sequence, inter-aircraft spacing, and arrival schedule based on the capabilities and constraints of the inbound aircraft, as well as the scheduled airport constraints (such as runway configuration, mixed-use runway use, dependent approaches, and weather conditions). The scheduling process will be coordinated with adjacent ATC centres and, when the schedule is frozen, a fixed area navigation (RNAV) arrival route (to the runway) with CTAs at some metering fixes will be provided. The assigned control times will be entered as RTAs by the on-board TEMO tool-set. Then, the on-board trajectory planner will compute the optimal trajectory plan such that a given cost function is minimised (typically fuel), while fulfilling the entered RTAs. The trajectory plan is computed by assuming more or less complex models for the aircraft dynamics, aircraft performance and weather conditions. When the models used by the FMS trajectory planning function do not match the reality, the initial trajectory is not longer the most optimal and some operational constraints may even be violated if errors are not nullified by the guidance system. Several guidance strategies were proposed as a part of the TEMO concept (De Jong et al.,2017,2015), which working principle will be explained in Section II.3. While all these methods have in common to manage the total energy of the aircraft (altitude and speed) such that CTAs are satisfied in an environmentally friendly way, they differ in the mechanisms implemented to correct time and energy deviations. At present, however, it is still not known which of these guidance strategies is the best in terms of environmental impact mitigation and capability to satisfy operational constraints (of utmost interest CTA). Initial batch studies of the TEMO concept were done using a conventional laptop to test its feasibility, reaching technology readiness level (TRL)-3 and proving lower fuel consumption and noise levels on ground, if compared with current step-down descents (De Jong et al.,2015). Thereafter, a human-in-the- loop study was performed in a high-fidelity flight simulator to look at human factors aspects, reaching TRL-4 and showing acceptable RTA adherence performance and operational acceptability by qualified pilots (De Jong et al.,2017). Yet, the models used by the trajectory planner were subject to several limitations. The fast optimizer for continuous descent approaches (FASTOP) project, funded by the CleanSky Joint Undertaking initiative, enhanced that version of the TEMO algorithm in order to test it in more realistic environments, aiming at the TRL-5 gate. The main improvements of the model were the consideration of realistic wind fields, non-standard atmospheres or curved routes. In addition, the TEMO software was redesigned from scratch allowing to use it in real-time on-board applications (Prats et al.,2014,2015b). In 2014, a second human in the loop study was performed in the Netherlands Aerospace Centre (NLR) full-motion flight simulator generic research aircraft cockpit environment (GRACE) to test TEMO in a more realistic simulated environment, achieving in this way TRL-5. The design of the experiment and the qualitative assessments gathered from pilots were reported by Bussink et al. (2016), while Prats et al. (2015a) presented a more detailed description of the trajectory optimisation algorithm and the quantitative results. In 2015, and aiming at TRL-6, the NLR in cooperation with the Delft university of technology and with the support of the flight operations for novel continuous descents (CONCORDE) consortium, executed some flight trials with a Cessna Citation II research aircraft. Several TEMO guidance variants were tested, including current step-down descents for benchmarking purposes. The preparation of this flight testing campaign was described by Verhoeven et al. (2016), and the results were presented by Prats et al. (2017). I.4.2.3 Wind modelling As mentioned in previous section, the accuracy of the trajectories computed by the FMS, and especially the computation of the estimated time of arrival over a fix in the flight plan, critically depends on the quality of the wind and, to a less extent, temperature data used during the trajectory planning process. Therefore, accurate weather predictions are of utmost importance to generate accurate trajectory plans and minimise the need of corrective actions during the execution of the descent, which typically restrict the accomplishment of energy-neutral descents. Historically, the flight crew have received the wind forecast prior to the flight from the flight dispatchers, and have manually entered this information into the FMS. At present, in order to simplify operations, wind forecasts can be also up-linked automatically via the aircraft communication and reporting system (ACARS). Typical FMS allow to store these wind data at several waypoints of the flight plan. Yet, at each one of these waypoints, wind data is typically available at only five altitudes. Then, the wind is interpolated between these altitudes and extrapolated beyond the limits when performing FMS trajectory predictions.
12 Chapter I - Introduction The wind data entered to the FMS are primarily based on wind charts from numerical weather prediction (NWP) models. Nowadays, the observations required to initialize NWP models are mainly gathered from radiosondes and aircraft equipped with aircraft meteorological data relay (AMDAR). The spatial distribution of the radiosondes, which are launched only two to four times a day, however, is too coarse, and wind observations gathered through AMDAR are not sufficient because not all aircraft in operation are equipped with such system (Haan,2010). Due to the relatively low spatio-temporal resolution of the data required to initialize NWP models and the computational burden of running a prediction, wind forecasts are generated over large areas of study only three to eight times per day, and are valid until +6, +12 or +24 hours beyond their issued time. Consequently, the wind forecast for the descent is not tailored for the current as-amended flight plan, and could be several hours old by the time the TOD is reached. Most flight planning suppliers currently use the world area forecast model (WAFC) that is produced to an ICAO specification by meteorological providers. Other suppliers allow to use higher resolution wind forecasts in the spatial domain, but do not provide a higher update rate due to the time required to uplink the data to the FMS (Bienert & Fricke,2013). Aiming to increase the update rate of wind forecasts, Boeing developed an advanced weather service that uplinks tailored wind data to the FMS via ACARS (Durham, 2011). This system, however, does not include additional altitudes at which wind data can be stored, but determines the best depending on the prevailing meteorological conditions of the flight plan. Some studies identified aircraft themselves as a potential network of airborne sensors emitting the sensed wind and temperature to ground receivers (Hollister et al.,1986) or directly to nearby aircraft (De Leege et al.,2013) to provide accurate, high-resolution and up-to-date weather data replacing or complementing the forecasts obtained from NWP models. I.5 Objectives of this PhD thesis The optimisation of the trajectory of an aircraft can be formulated as a multi-phase optimal control problem (Bryson & Ho,1975), which aims at finding the best control history (e.g., flight path angle and thrust) such that a given cost function is minimised (e.g., fuel consumption) while satisfying applicable operational constraints (e.g., a CTA). For real-life problems, obtaining the analytical solution is seldom. Consequently, most practical implementations transform the continuous optimal control problem into a finite-dimensional nonlinear programming (NLP) optimisation problem (Betts,2010), which can be solved using standard gradient-based optimisation techniques. The execution time required to solve the resulting NLP optimisation problem and the robustness of the algorithm, however, are critical aspects that need to be addressed when using these optimisation techniques for on-board applications. A first objective of this PhD thesis is to establish a framework for robust and real-time optimisation of aircraft trajectories, and also to provide guidelines to cope with the most common implementation issues that one’s might encounter. The paradigm of FMS considered in this PhD thesis generates an optimal trajectory plan (on-board and in real-time) whenever either the guidance system requests a trajectory plan update due to, e.g., an excessive time deviation, or the ATC up-links a new CTA. In both cases, the FMS attempts to generate an energy-neutral trajectory plan, i.e., correcting eventual deviations and/or adjusting the time of arrival at the metering fix without requiring additional engine thrust and keeping the active drag devices retracted. This might be achieved by using time and energy management concepts, using only elevator control to exchange potential energy by kinetic energy and vice-versa. For instance, if the aircraft had excess of potential energy and needed to arrive earlier than initially planned due to, e.g., a CTA update, altitude could be exchanged by speed by pitching down. The amount of time that can be absorbed or gained by using only elevator control, as a function of the state of the aircraft (i.e., altitude, speed, mass and distance to the metering fix) at the moment of updating the optimal trajectory plan, however, is still unknown. This PhD thesis aims to quantify the performance limits of energy-neutral trajectories. Any trajectory plan obtained after the optimisation process needs to be executed by the FMS in an uncertain environment, in which the models used to generate the plan are not exempt of errors. Guidance strategies build on time and energy management concepts are very attractive for executing CDOs subject to CTAs in an environmentally-friendly way. To the best of our knowledge, however, a comprehensive comparison of the existing guidance strategies based on time and energy management concepts has not been performed before. An important objective of this PhD is to assess the performance differences across various guidance strategies.
I.5 Objectives of this PhD thesis 13 In this PhD thesis, a new guidance strategy is also proposed as a part of the TEMO concept. This strategy is based on an iterative update of the optimal trajectory plan when new measurements of the current state of the aircraft and/or better knowledge of the surrounding environment are available (e.g., an improved wind prediction). This guidance strategy, commonly referred in the literature as model predictive control (MPC) - or nonlinear model predictive control (NMPC) for nonlinear systems - has been in use in the process industries in chemical plants and oil refineries since the 1980s. NMPC guidance strategies typically exploit the fact that the consecutive optimal control problems solved to update the trajectory plan are similar to each other. This allows to initialise the solution procedure efficiently by a suitably shifted guess from the previous trajectory plan, saving considerable amounts of computation time and improving convergence. The similarity of subsequent optimal control problems is even further exploited by strategies that never attempt to iterate the optimisation problem to convergence, and sensitivity-based trajectory updates that use parametric sensitivity information of the previous solution with respect to the uncertain parameters to generate acceptable approximations of the optimal solution quasi-instantaneously. To the best of our knowledge, this is the fist time that NMPC guidance and its sensitivity-based variants are implemented to satisfy time constraints during the descent of an aircraft. As explained before, one of the most uncertain parameters is the wind that the aircraft will encounter during the execution of the descent. An inaccurate wind forecast will both negatively effect CTA compliance and reduce the environmental benefits of CDO. A lot of research activities on this field are currently conducted, and future air-to-air communications to exchange accurate and up-to-date wind observations are expected to improve the quality of the wind data input to the FMS. The quantitative benefits, in terms of fuel consumption and time and energy deviations at the metering fix, of using wind observations emitted from nearby aircraft to update the trajectory plan in real-time are still unknown. Therefore, this PhD thesis proposes a method to include up-to-date weather observations in the trajectory management process under NMPC guidance, and assesses the benefits of this new concept that combines optimisation and estimation. Based on the discussions presented so far, the specific objectives of this PhD thesis can be outlined as follows: • Establish a trajectory management (planning and guidance) framework for the robust optimisation and execution of aircraft descent trajectories, on-board and in real-time. • Demonstrate the potential applicability of time and energy management concepts by quantifying the feasible time windows at a metering fix of energy-neutral descents subject to last minute CTA updates once the descent has been initiated (i.e., the aircraft has overflown the TOD). • Perform a comprehensive comparison of guidance strategies build on time and energy management concepts to satisfy CTAs during CDOs flying in an uncertain environment. The set of guidance strategies compared in this PhD thesis includes those already existing in the TEMO concept (i.e., tactical, strategic and hybrid) and the newly proposed NMPC. • Assess sensitivity-based alternatives to the generic NMPC methodology to overcome its main drawback, namely, the uncertain and in some cases excessive execution time required to update the trajectory plan. These guidance strategies make use parametric sensitivity information from the active optimal trajectory plan to frequently update it quasi-instantaneously. • Investigate the benefits of combining the NMPC stategy with wind networking concepts (i.e., sharing wind observations between aircraft), in terms of environmental impact and ability to comply with operational constraints. In addition, the sensitivity of the performance metrics to the rate of wind observations received from nearby aircraft will be also quantified.
14 Chapter I - Introduction I.6 Scope and limitations of this PhD thesis In order to accomplish the objectives of this PhD thesis, the research is subject to several assumptions and limitations that define its scope: • This PhD thesis only focuses on the on-board systems required to satisfy CTAs during CDOs. Accordingly, it is assumed that the CTAs is provided by an hypothetical ATC through a data-link communications system. Furthermore, the time constraint is computed by an ideal AMAN, which is able to accurately sequence the arrival traffic and ensure sufficient spacing. • The concept of operations investigated assumes a closed-loop path from the TOD to the runway threshold. Therefore, updates of the optimal trajectory plan due to route changes are not contemplated. Changes of the lateral route, however, could be translated to changes in the distance to go, and use the proposed algorithms with minor modifications. • All the experiments of this PhD thesis evaluate the best a single aircraft can do given a CTA and route, thus neglecting interaction with surrounding traffic. The potential effectiveness of the proposed planning and guidance algorithms are not investigated through a traffic simulation, and discussions on their benefits in the traffic stream management are not performed. • Aircraft emissions consist primarily of carbon dioxide (CO2), nitrogen oxides (NOx), and methane. Of additional concern are the generation of contrails, the sulphur oxides (SOx) and the particulate matter (PM). In the recent years, several studies have proposed trajectory optimisation strategies to mitigate the environmental impact of aviation considering some of these factors (see for instance the work performed by Hartjes & Visser (2017)). In this PhD thesis, however, fuel consumption is the only environmental impact metric considered for the optimisation of aircraft trajectories. Yet, airframe noise is indirectly considered in the optimisation algorithm by also minimising the number of speed brakes deployments. In addition, in certain experiments a noise model is used to asses, ad-hoc, the noise footprints of the trajectories. • The performance assessment of the planning and guidance strategies is done through simulation, which requires an aircraft performance model. All the results presented in this PhD thesis are only valid for the Airbus A320, which represents the vast majority of worldwide narrow-body commercial aircraft. The proposed planning and guidance algorithms, however, are generic enough to consider different aircraft types by adopting the corresponding performance model. • It it assumed that speed brakes and throttle are the only mechanisms to remove and add energy from the system to correct disturbances, respectively. Despite flaps and gear might be more efficient than deploying speed brakes or decreasing the thrust, they have not been considered in this PhD thesis as means to control the energy rate of the aircraft. In spite of that, they have been modelled in some experiments to represent realistic approaches. • The trajectory optimisation problem tackled in this PhD involves nonlinear constraints and/or cost function and is, in consequence, a non-convex optimisation problem. This means that several local minima may exist. Since global optimisation algorithms require a high computational burden, they are prohibited for real-time applications such as the on-board generation of optimal trajectory plans when using NMPC guidance. Thus, the proposed trajectory optimisation techniques involve the use of local optimisation algorithms. • Contingencies, such as downgraded guidance or engine failure, are not considered in the experiments. In other words, normal operations are assumed in all the experiments.
I.7 Outline of this PhD thesis 15 I.7 Outline of this PhD thesis The present document is organised in seven chapters, which are summarised as follows. It is worth noting that a broad state-of-the-art of the main topics addressed in this PhD thesis has been presented before. A deeper and more specific state-of-the-art for each individual topic is included at the beginning of the chapter that addresses it. •Chapter II presents a detailed framework for trajectory management, which includes optimisation techniques to generate optimal trajectory plans in real-time, and various guidance strategies to execute the resulting plans using time and energy management concepts. •Chapter III quantifies the feasible time window for energy-neutral CDO trajectories as a function of the altitude and distance to go; and analyses the sensitivity of the feasible time window to different winds, metering fix position, aircraft masses, and initial speeds. •Chapter IV compares the performance of four guidance strategies designed to satisfy time constraints during a continuous descent operation: tactical, strategic, hybrid and NMPC. •Chapter Vproposes two variants of the generic NMPC guidance strategy that rapidly update the optimal trajectory plan by using parametric sensitivities of the active trajectory plan. Then, the performance of these two variants are compared with those of an ideal NMPC, which updates the optimal trajectory plan instantaneously at each re-calculation instant, as well as those of the open-loop execution of the initial trajectory plan computed before the TOD. •Chapter VI assesses the benefits, in terms of environmental impact and CTA compliance, of combining the NMPC guidance strategy with wind networking concepts to generate accurate and up-to-date wind predictions on-board before updating the optimal trajectory plan. •Chapter VII gives the conclusions that are drawn from this work and points out some future work that could be done in the direction of the presented research.
22 Chapter II - Framework on trajectory management Note that since in the concept of operations tackled in this PhD thesis the lateral route (sequence of waypoints) is known before generating the optimal vertical profile, realistic air temperature and pressure models will be essentially function of the altitude and distance to go. A similar rationale is adopted to model the wind field. Simple models used for academic purposes assume constant wind speed and direction throughout the descent, independently of the geographical location and altitude. In practise, however, wind speed tends to increase with the altitude. This typical behaviour can be represented using various empirical and theoretical models, from linear functions to more or less complex mathematical expressions. An example of empirical model for the total wind speed (|~w|), which assumes constant longitudinal wind above a reference altitude (hr) is the well-known Hellmann exponential law: |~w|(h) = (wrif h > hr wrh hrα if h≤hr,(II.7) where wris a known wind speed at hr. The Hellmann exponent (α) depends on many factors, such as the coastal location, the shape of the terrain on the ground, and the stability of the air. Figure II-1(b) shows the Hellmann exponential law for different values of α, fixing wrto 30 kt and hrto FL300. State-of-the-art weather forecasts obtained from NWPs typically provide the wind vector at different altitudes and geographical locations. These forecasts decompose the wind vector in north component (wn), east component (we) and vertical component. Again, since the route is known beforehand in this PhD thesis, the bearing between consecutive waypoints of the route is also known, and the wind field can be projected to along-track and cross-track wind components as a function of the altitude and distance to go: ws(h, s) = wn(h, s) cos χg(s) + we(h, s) sin χg(s)(II.8a) wx(h, s) = −wn(h, s) sin χg(s) + we(h, s) cos χg(s),(II.8b) where χgis the track angle, a known function of the distance to go. In should be noted that the temporal variations of temperature, pressure and wind components are neglected in this PhD thesis because the duration of a descent is small if compared with the time scale of typical atmosphere dynamics. In addition, the trade-off between additional complexity (thus execution time) of the model and accuracy of the resulting trajectories is not appealing. Yet, readers interested on the extension of the model that considers temporal variations of the weather variables could consult the model proposed by Prats et al. (2014). For each of the experiments conducted in this PhD thesis, the expressions for δ,θ,wsand wxwill be particularised as a function of hand s, either with mathematical expressions approximating realistic weather data from NWP forecasts or using empirical or theoretical models. II.2 Optimal trajectory planning As mentioned before, the optimisation of the vertical profile of an aircraft trajectory can be formulated as a multi-phase, constrained optimal control problem (Soler et al.,2015). The formulation of a generic optimal control problem is presented in Section II.2.1. Since the trajectory optimisation problem tackled in this PhD thesis cannot be solved analytically due to the nonlinearity and complexity of the constraints and cost function, a discretisation approach is presented in Section II.2.2 to solve it numerically. II.2.1 Optimal control problem formulation An optimal control problem for a generic system during a fixed or variable continuous time horizon [tI, tF] can be formulated as follows (Bryson & Ho,1975):
II.2 Optimal trajectory planning 23 min ~u(t)J:= φ~x(tF),~ d+ZtF tI π~x(t), ~u(t),~ ddt s.t ~x(tI) = ~ X ˙ ~x =~ f~x(t), ~u(t),~ d ~ bin ~x(t), ~u(t),~ d≤0 ~ beq ~x(t), ~u(t),~ d= 0 ~ ψ~x(tF),~ d= 0, (II.9) where ~x ∈ <nxis the vector of differential states, ~u ∈ <nuis the vector of controls, and ~ d∈ <ndis the vector of fixed parameters of the model. The cost function J:<nx× <nu× <nd→ <, which is composed by a running cost (or Lagrange term) π:<nx×<nu×<nd→ <and an end cost (or Mayer term) φ:<nx×<nd→ <, is to be minimised subject to: dynamic constraints ~ f:<nx×<nu×<nd→ <nxin the form of ODEs with initial conditions ~ X∈ <nx, algebraic constraints~ beq :<nx×<nu×<nd→ <nbeq , inequality path constraints ~ bin :<nx×<nu×<nd→ <nbin , and terminal constraints ~ ψ:<nx×<nd→ <nψenforced only at tF. In some applications, free parameters that, in contrast to the controls, are constant in time might be also included as a part of the decision variables. For instance, when the final time is fixed, tFis a known and fixed parameter. For optimal control problems with variable final time, however, tFbecomes a new decision variable to be optimised, in addition to ~u(t). It should be noted that the vector of fixed parameters of the model (~ d), however, is not part of the decision variables and, consequently, must be chosen by the user and remains constant during the whole optimisation process. Yet, it might change for different problems. For instance, parameters that could be included in ~ dare the coefficients α,wr, and hrof Eq. (II.7) if the Hellmann model were considered. The optimal control problem described by Eq. (II.9) assumes that the same running cost, dynamic constraints and algebraic and path constraints apply during the whole time horizon. In addition, event constraints can be set only at the very end of the time horizon. Yet, many real-life processes can be divided into several phases (or stages), where the dynamics of the system, the running cost and the algebraic and path constraints might change. In addition, in some particular applications it is necessary to formulate interior-point constraints between two consecutive phases. Even in some applications, it could be desirable to define different states and controls vectors for each phase. In this PhD thesis, different phases could have different algebraic, path and interior-point constraints, yet they must have the same number of differential dynamic constraints, because it has been assumed that the state vector is identical in all phases. In order to extend the optimal control problem of Eq. (II.9) to multiple phases, let the time horizon [tI, tF]be divided into Pknown intervals [tj, tj+1], with j= 0, . . . , P −1. Note that t0=tIand tP=tF. Each one of these intervals corresponds to a different phase. The generic multi-phase optimal control problem defined over Pphases can be formulated as: min ~u(t)J:= P−1 X j=0 φj~x(tj+1),~ d+Ztj+1 tj πj~x(t), ~u(t),~ ddt! s.t ~x(t0) = ~ X ˙ ~x =~ fj~x(t), ~u(t),~ d;j= 0, . . . , P −1 ~ bin j~x(t), ~u(t),~ d≤0; j= 0, . . . , P −1 ~ beq j~x(t), ~u(t),~ d= 0; j= 0, . . . , P −1 ~ ϑeq j~x(tj+1),~ d= 0; j= 0, . . . , P −2 ~ ϑin j~x(tj+1),~ d≤0; j= 0, . . . , P −2 ~ ψ~x(tP),~ d= 0, (II.10)
24 Chapter II - Framework on trajectory management where ~ fj:<nx×<nu×<nd→ <nxand πj:<nx×<nu×<nd→ <are the differential dynamic constraints and running cost of the jth phase, respectively. Similarly,~ beq j:<nx×<nu×<nd→ <nbeq jand~ bin j:<nx×<nu× <nd→ <nbin jare the algebraic and path constraints, respectively, of the jth phase; and ~ ϑeq j:<nx×<nd→ <nϑeq jand ~ ϑin j:<nx×<nd→ <nϑin jrepresent applicable equality and inequality interior-point constraints, respectively, applying at the last time of the jth phase. When dealing with mutli-phase optimal control problems, a new set of constraints is added to link the state variables across two consecutive phases and guarantee continuity in the solution: ~x(t+ j+1)−~x(t− j+1) = 0; j= 0,...P −2.(II.11) Analogously to the single-phase optimal control problem, for multi-phase optimal control problems the final time of each phase could be either a known parameter or a decision variable to be optimised. For instance, if only the final time of the whole time interval were fixed (e.g., due to a CTA) but the duration of each phase were flexible, tPwould be a known parameter and tj+1,j= 0, . . . , P −2, additional decision variables of the optimisation problem, in addition to ~u(t). Basically, two different methods are available for solving optimal control problems: •Indirect methods: consist of first deriving the necessary conditions for optimality in an analytical way, and then solve them numerically by using root-finding solvers and ODE integrators. These methods are commonly referred in the literature to optimise-then-discretise. •Direct methods: consist of first discretising the time history of the controls and/or states to transcribe infinite-dimensional optimal control problem Eq. (II.9) into a finite-dimensional nonlinear programming (NLP) optimisation problem, and then attempt a direct numerical solution. These methods are commonly referred in the literature to discretise-then-optimise. In turn, indirect methods can be further classified into single-shooting and multiple-shooting, while the direct methods can be divided into single-shooting, multiple-shooting and collocation (see Fig. II-2). An excellent review of these methods was performed by Andersson (2013), illustrating the working principle of them by introducing a simple, yet industrially relevant optimal control problem. Prats (2010) also described the main features and the trade-offs of each method, in terms of complexity and computational burden, from a more practical perspective. Optimal control methods Direct methods (discretise-then-optimise) Indirect methods (optimise-then-discretise) Direct collocation Indirect Single-shooting Direct Single-shooting Direct Multiple-shooting Indirect Multiple-shooting Figure II-2: Methods to solve optimal control problems The major drawback of indirect methods is the requirement for a detailed mathematical analysis of each single optimal control problem. Even small changes in the definition of the constraints can lead to a completely different solution, often requiring a complete derivation of the necessary conditions for optimality. These methods also require a very accurate initial guess for the adjoint variables. These variables, however, do not have any physical meaning, and the optimal solution is highly sensitive to even small changes of them (Prats,2010). Last but not least, including path inequality constrains in the formulation of the optimal control problem can be prohibitely difficult when using indirect methods. For these reasons, they were not implemented in this PhD thesis.
II.2 Optimal trajectory planning 25 Direct methods are much easier to set up, if compared with the indirect ones, and allow to solve very complex problems with a minimum effort of mathematical analysis. Using direct methods, the whole time horizon is divided into several time intervals, and the control vector is parameterised with a low order polynomial on each interval (e.g., a piecewise constant control). In direct single-shooting, the initial state vector ( ~ X) is propagated forward in time from tIto tFby integrating the differential equations (~ f) with ODE/DAE integrators. As a result, the state vector at tF becomes a function of the parametrised controls on each interval and the initial state vector. The integral of the running cost (π) is also calculated to evaluate the cost incurred during the whole time horizon, which is to be minimised. Then the value of the terminal constraint (~ ψ) is evaluated and, by using a NLP optimisation algorithm, the parametrised controls and the initial state vector are adjusted in order to satisfy it. In other words, with the single-shooting approach, only the time history of the controls is parametrised, eliminating the state variables with ODE/DAE integrators. This method is effective for problems that are either simple or have an extremely good initialisation. The main drawback of direct single-shooting is that the objective and constraint are highly nonlinear dependent on ~u, thus slowing down the convergence of the NLP solution (Andersson,2013). In addition, if a guess for the state variables is available, e.g., from a previously computed trajectory, it cannot be used to initialise and speed-up the solution. Multiple-shooting methods address these issues by including the state trajectory in the NLP, essentially breaking the problem into shorter steps. The working principle of multiple-shooting consists on parametrising not only the controls, but also the state vector at the beginning of each interval. On each interval, the initial state vectors can be propagated to the beginning of the next one by integrating the differential equations with ODE/DAE integrators. As a result, the state vector at the end of each interval becomes a (nonlinear) function of the state and control vectors at its beginning. The discretised state vectors could be pieced together to form a continuous trajectory if and only if their values match at the junction points. The discontinuities at the junction points, called defects, are included in the constraints vector as a system of nonlinear equations that depend on the parametrised states and controls. Finally, the resulting finite-dimensional optimisation problem is solved by using standard NLP optimisation algorithms. A deep analysis on the proprieties of direct shooting methods was performed by Gath (2002) and Bulirsch et al. (1993). In this PhD thesis, this approach was selected to solve the trajectory optimisation problem. For a detailed explanation of the remaining direct collocation methods, the reader is referred to Betts (2010). The main difference between direct collocation and multiple-shooting methods is that the former perform an implicit integration of the dynamic system on each time interval, while multiple-shooting methods use explicit integration formulas, such as Euler or Runge-Kutta. One of the major advantages of using multiple-shooting methods, if compared to direct collocation, is the possibility of running them with ODE integrators that include step-size control techniques, resulting in an algorithm almost independent of the discretisation grid that will, if it converges, at least deliver a suboptimal solution Gath (2002). A very interesting survey on existing optimal control software packages was performed by Virtanen et al. (1999). II.2.2 Direct multiple-shooting formulation Let the continuous time horizon [tI, tF]of Eq. (II.9) be discretised into N+ 1 equidistant time samples τk, with k= 0, . . . , N. Note that τ0=tI,τN=tFand the discretisation step is ∆τ=(τN−τ0)/N. The optimal control problem discretised by using direct multiple-shooting methods and minimising Jin a time horizon of Ntime intervals can be formulated as: min ~xk,k=0,...,N ~uk,k=0,...,N−1 J:= φ~xN,~ d+ N−1 X k=0 Π~xk, ~uk,~ d, ∆τ s.t ~x0=~ X ~xk+1 =~ F~xk, ~uk,~ d, ∆τ;k= 0, . . . , N −1 ~ bin ~xk, ~uk,~ d≤0; k= 0, . . . , N −1 ~ beq ~xk, ~uk,~ d= 0; k= 0, . . . , N −1 ~ ψ~xN,~ d= 0, (II.12)
26 Chapter II - Framework on trajectory management where ~xk∈ <nxand ~uk∈ <nuare the state and control vectors discretised at τk, respectively, for k= 0, . . . , N. While the algebraic, path and terminal constraints of the continuous optimal control problem Eq. (II.9) can be directly evaluated at the discretised states and controls, the states evolution function ~ F: <nx× <nu× <nd× < → <nxmust be defined as the result of integrating ~ fduring a time interval of duration ∆τusing an appropriate integration scheme with the discretised states and controls. Similarly, Π : <nx×<nu×<nd×< → <is a quadrature function resulting from a numerical integration of πduring a time interval of duration ∆τ. For instance, using the Euler method with one integration step these functions would be: ~ F(~xk, ~uk,~ d, ∆τ) = ~xk+~ f~xk, ~uk,~ d∆τ;k= 0, . . . , N −1(II.13a) Π(~xk, ~uk,~ d, ∆τ) = π~xk, ~uk,~ d∆τ;k= 0, . . . , N −1.(II.13b) where ∆τcould be a known parameter for optimal control problems with fixed final time, or an additional decision variable to be optimised in Eq. (II.12) if the final time were free. It should be noted, however, that Eq. (II.13) represents the most simple method to integrate ~ fand πduring ∆τ. In this PhD thesis, more complex and accurate techniques have been used, which will be presented in the following chapters. Analogously to Eq. (II.10), in the discrete domain the optimal control problem can be also defined over more than one phase. First, let the continuous time horizon [tI, tF]be divided into Ptime intervals [tj, tj+1] for j= 0, . . . , P −1; each time interval corresponding to a different phase. Again, t0=tIand tP=tF. Then, each time interval (or phase) is discretised into Njequidistant time samples τk, τk+1, . . . , τk+Nj−1, where τk=tj,τk+Nj−1=tj+1 and k=Pi<j Ni, for all j= 0, . . . , P −1. The discretisation step of the jth phase is denoted by ∆τj. As a result, the whole time horizon is discretised again into N+ 1 = PP−1 j=0 Nj time samples [τ0, τ1, . . . , τN]. Let Tbe a multi-dimensional set that relates the index of each phase to the indexes of its corresponding time samples. The subset E ⊆ T only includes the index corresponding to the last time sample of each phase; and Iis defined as T\E. For instance, for a multi-phase optimal control problem with P= 3 phases and N+ 1 = 12 time samples equally distributed among them (i.e., Nj= 4 for j= 0,1,2), these subsets would be particularised as: T={(0,0),(0,1),(0,2),(0,3),(1,4),(1,5),(1,6),(1,7),(2,8),(2,9),(2,10),(2,11)} E={(0,3),(1,7),(2,11)} I={(0,0),(0,1),(0,2),(1,4),(1,5),(1,6),(2,8),(2,9),(2,10)}. Figure II-3 shows an illustrative scheme of the discretisation process for an arbitrary number of phases and time samples. Based on the definitions stated above, the discrete multi-phase optimal control problem is: min ~xk,k=0,...,N ~uk,∀(j,k)∈I J:= X (j,k)∈E φj~xk,~ d+X (j,k)∈I Πj~xk, ~uk,~ d, ∆τj s.t ~x0=~ X ~xk+1 =~ Fj~xk, ~uk,~ d, ∆τj;∀(j, k)∈ I ~ bin j~xk, ~uk,~ d≤0; ∀(j, k)∈ I ~ beq j~xk, ~uk,~ d= 0; ∀(j, k)∈ I ~ ϑeq j~xk,~ d= 0; ∀(j, k)∈ E\{(P−1, N) ~ ϑin j~xk,~ d≤0; ∀(j, k)∈ E\{(P−1, N) ~ ψ~xN,~ d= 0 ~xk−~xk+1 = 0; ∀(j, k)∈ E\{(P−1, N)}. (II.14)
II.2 Optimal trajectory planning 27 Figure II-3: Discretisation process for a generic optimal control problem In Eq. (II.14), Πj:<nx×<nu×<nd×< → <and ~ Fj:<nx×<nu×<nd×< → <nxare the quadrature and states evolution functions for the jth phase, respectively. The discretisation step of each individual phase could be considered either known parameter or variable to be optimised, depending on the problem. For instance, if the duration of the whole time horizon were fixed to a certain parameter, say a CTA, but the duration of each phase were flexible, ∆τjfor j= 0, . . . , P −1would become decision variables subject to: P−1 X j=0 (Nj−1) ∆τj−CTA = 0,(II.15) which would be appended as an equality constraints in Eq. (II.14). II.2.3 Nonlinear programming problem As mentioned above, direct methods for optimal control compute the optimal trajectory plan by formulating the generic multi-phase optimal control problem Eq. (II.14) as a parametric NLP optimisation problem. This NLP optimisation problem has the following form: min ~z f(~z, ~p) s.t ~ h(~z, ~p)=0 ~g (~z, ~p)≤0, (II.16) where ~z ∈ <nzis the vector of primal variables, ~ h:<nz×<np→ <nhis the vector of equality constraints, ~g :<nz×<np→ <ngis the vector of inequality constraints, and ~p ∈ <npis the vector of (fixed) parameters of the NLP optimisation problem. In this PhD thesis, the following definitions are considered: ~z := ~z0 ~z1 . . . ~zN ,~ h(~z, ~p) := ~ h0(~z, ~p) ~ h1(~z, ~p) . . . ~ hN(~z, ~p) , ~g (~z, ~p) := ~g0(~z, ~p) ~g1(~z, ~p) . . . ~gN(~z, ~p) ,(II.17) where:
28 Chapter II - Framework on trajectory management ~zk:= "~uk ~xk#if (j, k)∈ I ~xkif (j, k)∈ E ~ hk:= ~xk+1 −~ Fj~xk, ~uk,~ d, ∆τj ~ beq j~xk, ~uk,~ d if (j, k)∈ I "~ ϑeq j~xk,~ d ~xk−~xk+1 #if (j, k)∈ E\{(P−1, N)} ~ ψ~xk,~ dif k=N ~gk:= ~ bin j~xk, ~uk,~ dif (j, k)∈ I ~ ϑin j~xk,~ dif (j, k)∈ E\{(P−1, N)}. (II.18) According to Eq. (II.18), ~zkincludes both discretised states and controls at the time sample τk. Similarly, ~gk and ~ hkinclude the inequality and equality constraints applied at τk, respectively. In Eq. (II.16), fis the cost function of the original optimal control problem evaluated at the primal variables and NLP parameters, i.e., f(~z, ~p) = J(~z, ~p). In this PhD thesis, the vector of NLP parameters is composed of both current state of the aircraft and parameters of the model: ~p ="~ X ~ d#.(II.19) Furthermore, in order to reduce the number of NLP variables and constraints, the constraint that fix the initial conditions of the optimal control problem to the current state of the system is eliminated by substituting the variables ~x0for the vector of fixed initial conditions ~ Xin the whole NLP optimisation problem (i.e., equality constraints, inequality constraints and cost function). This allows to remove the constraint ~x0=~ Xand the variable ~x0from the optimisation problem. The Lagrangian function associated to the NLP optimization problem (II.16) is: L~z, ~p,~ λ, ~µ:= f(~z, ~p) + ~ λT~g (~z, ~p) + ~µT~ h(~z, ~p),(II.20) where ~ λ∈ <ngand ~µ ∈ <nhare the Lagrange multipliers (dual variables) vectors paired up with the constraints ~g and ~ h, respectively. An optimal primal-dual pair ~z∗,~ λ∗, ~µ∗, where (·)∗indicates optimal, satisfies the first-order necessary conditions of optimality, also known as Karush-Kuhn-Tucker conditions (KKT) (Kuhn & Tucker,1950), if: ∂L ∂~z ~z∗, ~p,~ λ∗, ~µ∗=∂f ∂~z (~z∗, ~p) + ~ λ∗T∂~g ∂~z (~z∗, ~p) + ~µ∗T∂~ h ∂~z (~z∗, ~p)=0 (II.21a) ~ h(~z∗, ~p) = 0; ~µ∗>0(II.21b) gi(~z∗, ~p) = 0; λ∗ i>0; ∀(gi, λ∗ i)∈ Gac (II.21c) gi(~z∗, ~p)<0; λ∗ i= 0; ∀(gi, λ∗ i)∈ Gin,(II.21d) where the pairs of inequality constraints and associated dual variables have been divided into two complementary sets: the active set (Gac) and the inactive set (Gin) (W¨ urth et al.,2009). Newton methods suitable to solve this kind of optimisation problems try to find a point satisfying Eq. (II.21), by using successive linearisation of the cost function and constraints. The major differences between them are on how to achieve the conditions related with the inequality constraints Eqs. (II.21c) and (II.21d). The two big families are:
II.2 Optimal trajectory planning 29 •Interior-point methods, which typically replace the original NLP optimisation problem by a series of barrier subproblems (Wachter & Biegler,2006), which are controlled by a barrier parameter. Two representative interior-point methods are the primal-dual method, and the trust-region method. •Active-set methods, which solve in turn a sequence of subproblems based on a quadratic approximation of the original NLP optimisation problem. Representative active-set methods are the sequential quadratic programming (SQP) method (Gill et al.,2005), the sequential linear-quadratic programming (SLQP) method (Byrd et al.,2003), and the gradient projection method (Lin & Mor´ e,1999). Let PNbe the NLP algorithm (either interior-point or active-set) that provides the optimal primal-dual solution (i.e., satisfying Eq. (II.21)) as a function of ~p for the next Ntime samples: ~z∗,~ λ∗, ~µ∗← PN(~p).(II.22) Remember that the vector of NLP parameters includes both current states and parameters of the model (see Eq. (II.19)). Taking into account that, by definition of Eq. (II.22), the optimal primal-dual solution is an implicit function of ~p, any change in ~p would lead to a new optimal primal-dual solution. This fact will be exploited in Section II.3.4.2 when considering parametric sensitivity methods to update the optimal trajectory plan quasi-instantaneously. II.2.4 Implementation considerations Several common issues must be considered when dealing with NLP optimisation problems. In particular, how to react in the face of an infeasible solution in which some constraints are not satisfied, how to handle discontinuities present in tabular data, and which mechanisms should be used to speed up the optimisation process when solving a sequence of similar problems. The following sections give a couple of recommendations related to these three topics, respectively. II.2.4.1 Recovering from infeasible solutions In practical applications, the solution of the NLP optimisation problem described by Eq. (II.16) may not always be feasible. Accumulated deviations caused by disturbances and modelling errors could eventually lead to the impossibility to satisfy one or several constraints to a given numerical precision. If this were the case, the NLP solver would not provide a valid solution and, consequently, the optimal trajectory plan could not be obtained. Yet, in some circumstances it may be reasonable to allow, to some extent, violations in certain constraints, so that a sub-optimal solution can still be obtained, as long as the safety requirements of the application permit. With this aim, there exist several alternatives to recover the feasibility of the solution. Constraint relaxation method An attractive method consists of not defining the admissible range of the constraints to their actual operational limits but to more conservative values. Then, these limits will be systematically relaxed whenever an infeasible solution is obtained. This approach is commonly referred as the constraint relaxation method. It is not always obvious, however, which constraints to relax, nor the amount which they should be relaxed. In order to overcome this issue, some optimisation techniques have been proposed, which determine the constraints that should be relaxed, according to some priority list, and the minimum relaxation required to recover feasibility (Vada et al.,2001). These techniques resort to optimisation problems running in parallel to the original optimal control problem. Accordingly, their algorithmic complexity may be significant. Transient constraint violation method A second alternative consists of allowing constraint violations only during a certain period of time, typically at the beginning of the new trajectory, which duration is to be minimized. In some cases, however, the time transient constraints violations could be large (Scokaert & Rawlings,1999).
30 Chapter II - Framework on trajectory management Soft-constraint method Finally, the soft-constraint method consists of defining additional slack variables in the optimisation problem, which allow to violate the path and/or terminal constraints in case of not obtaining a feasible solution to the optimal control problem (Scokaert & Rawlings,1999). Using this method, the original NLP optimisation problem, Eq. (II.16), is re-formulated as follows: min ~z,~eeq ,~ein f(~z, ~p) + ~einTWin~ein +~eeqTWeq~eeq s.t ~ h(~z, ~p) = ~eeq ~g (~z, ~p)≤~ein ~ein ≥0, (II.23) where ~ein ∈ <ngand ~eeq ∈ <nhare the vector of slack variables paired up with the inequality and equality constraints of the original NLP optimisation problem, respectively. Additional constraints ensure that the slack variables related to inequality constraints do not make the associated constraints more restrictive than their original values. The positive-definite weighting matrices Win ∈ <ng×ngand Weq ∈ <nh×nhdetermine how much the violation of inequality and equality constraints is penalised, respectively. By educatedly tuning these weights, it is expected that the optimisation algorithm will violate constraints only if necessary. Otherwise, the result shall be equivalent to that obtained from solving (II.16). The rule of thumb known as Bryson’s rule (Bryson & Ho,1975) can be used as a guideline to determine these weighting matrices. This rule states that the weights should be chosen so that the contribution of each slack variable to the cost function is approximately the same. In order to accomplish that, the inverse square of the admissible range of each constraint could be used as weight for the associated slack variable. Last but not least, it should be noted that using this technique not all the constraints need to be necessarily slacked. For these constraints which cannot be violated under any condition, the corresponding slack variable could be either removed from the set of variables, or the weight associated to the slack variable could be set to a very large value. Yet, the later solution is not recommended for numerical stability reasons. II.2.4.2 Handling tabular data in NLP algorithms On optimal control problems such as Eq. (II.9) one often encounters a cost function and/or constraints given in tabular form, as a function of the states and/or controls. For instance, the drag coefficient of the aircraft, which depends on the lift coefficient and Mach number, could be given as a look-up table with performance data obtained from the manufacturer. Similarly, wind data from NWP forecasts, which are required to compute the ground speed, are typically given in tabular form as a function of the latitude, longitude and pressure level. For use in simulations, the (linear) interpolation between discrete data points is a simple and convenient way to extract information from look-up tables. However, typical NLP solvers require the first and second derivatives of the cost functions and constraints with respect to the primal variables in order to work. The derivatives of a piecewise linear function are undefined at the data points, where two linear functions with different slopes intercept. For this reason, state-of-the-art NLP solvers have difficulties when using tabular data. In this kind of optimisation problems, the data must be approximated off-line in order to generate continuous and twice-differentiable curves by some smooth function fitting techniques. The easiest way to approximate tabular data is perhaps to use polynomial fittings. Although some of the aircraft performance data could be approximated quite well by polynomial functions, these expressions are not able to accurately approximate sharp change such as the drag rise region due to compressibility effects or complex functions such as the wind field. In these cases, high-order polynomials prone to oscillation due to the Runge’s phenomenon would be required, potentially leading to poor convergence, local minima issues and large errors between data points. In order to deal with tabular data in optimal control problems, Betts (2010) suggested to approximate the data points with tensor product splines. Splines are functions defined piecewise by polynomials of any degree, which keep sufficiently continuity and differentiability at the places where the polynomial pieces connect, called knots (de Boor,1972).
II.2 Optimal trajectory planning 31 Altitude [FL] 050 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 Tmax [kN] 20 30 40 50 60 (a) Interpolating spline Altitude [FL] 050 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 Tmax [kN] 20 30 40 50 60 (b) Smoothing spline Altitude [FL] 050 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 ∂Tmax ∂h −100 −80 −60 −40 −20 0 20 (c) Interpolating spline: first derivative with respect to altitude Altitude [FL] 050 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 ∂Tmax ∂h −50 −40 −30 −20 −10 0 10 (d) Smoothing spline: first derivative with respect to altitude Altitude [FL] 050 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 ∂2Tmax ∂h2 −7500 −5000 −2500 0 2500 5000 7500 (e) Interpolating spline: second derivative with respect to altitude Altitude [FL] 050 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 ∂2Tmax ∂h2 −300 −200 −100 0 100 200 300 400 (f) Smoothing spline: second derivative with respect to altitude Figure II-4: Comparison between interpolating and smoothing spline approximation for the maximum thrust of an A320
38 Chapter II - Framework on trajectory management accomplished by computing the error in the Lagrange sensitivity (opt) and the nonlinear constraint infeasibility (infs) at the updated pair of primal-dual variables (Wolf & Marquardt,2016): opt =k∂L∗ i: ∂~z k∞ k~ λ∗ i:k2+k~µ∗ i:k2+ 1, infs =k~g∗+ i:k∞+k~ h∗ i:k∞ k~z∗ i:k2+ 1 .(II.37) If these metrics were higher than pre-defined values, the set of first and second-order sensitivities F∗ i::= {∂2L∗ i: ∂~z2 i: ,∂2L∗ i: ∂~zi:∂~p ,∂f∗ i: ∂~zi:,∂~g∗ i: ∂~zi:,∂~g∗ i: ∂~p ,∂~ h∗ i: ∂~zi:,∂~ h∗ i: ∂~p }would be updated with the new ~z∗ i:,~ λ∗ i:, ~µ∗ i:solution and ~p, and further iterations of QN−iwould be triggered until satisfying feasibility and optimality criteria. Summing up, if ∆~p is small, it is not really necessary to solve PN−iat each time sample. In the neighbourhood of the active optimal solution, parametric sensitivities can be calculated to rapidly update the optimal solution at τiusing a first-order approximation and checking for optimality and feasibility criteria. This fast trajectory update is performed by solving QN−i, which is expected to be faster than PN−i. A virtue of this method is that F∗ i:can be evaluated in the background at the preceding optimal solution assuming the estimated parameters ~p (W¨ urth et al.,2008a;Gros et al.,2009). Then, the optimal solution is updated on-line by solving QN−ialmost instantaneously right after measuring and/or estimating ~p. This is why this method is commonly referred as sensitivity-based NMPC (SbNMPC). Algorithm (II.2) shows the main steps of SbNMPC, where τopt and τinfs are the pre-defined tolerances for the optimality and feasibility criteria, respectively. The contingency step (line 17 in Algorithm II.2) is optional. An alternative is to not update the trajectory if the optimality and feasibility criteria are not satisfied after kmax QP iterations. Algorithm II.2: Working principle of the sensitivity-based NMPC (SbNMPC) 1: ~z∗,~ λ∗, ~µ∗← PN(~p) 2: Build F∗ 1: from F∗ 3: for i= 1, . . . , N −1do 4: Measure ~ Xand estimate ~ d 5: ~p ←"~ X ~ d# 6: ∆~p ←~p −~p(τi−1) 7: ∆~z∗ i:,~ λ∗ i:, ~µ∗ i:← QN−i(∆~p) 8: ~z∗ i:←~z∗ i:+ ∆~z∗ i: 9: Evaluate F∗ i:at the new (~z∗ i:,~ λ∗ i:, ~µ∗ i:) and ~p 10: k←1 11: while (opt > τopt |infs > τinfs) & k <=kmax do 12: ∆~z∗ i:,~ λ∗ i:, ~µ∗ i:← QN−i(0) 13: ~z∗ i:←~z∗ i:+ ∆~z∗ i: 14: Evaluate F∗ i:at the new (~z∗ i:,~ λ∗ i:, ~µ∗ i:) and ~p 15: k←k+ 1 16: if k== kmax then 17: ~z∗ i:,~ λ∗ i:, ~µ∗ i:← PN−i(~p)using warm start 18: Build F∗ i+1 from F∗ i: 19: Build the guess for ~z∗ i+1:,~ λ∗ i+1: and ~µ∗ i+1: from ~z∗ i:,~ λ∗ i:and ~µ∗ i:, respectively 20: Implement ~u∗ iuntil τi+1 It should be noted that sensitivity-based trajectory updates presented in this section will work only for small perturbations of the parameters vector, since they are based on a linear approximation of the necessary conditions of optimality. For large perturbations, the linear approximation will not be accurate enough, and a large amount of QP iterations might be required to satisfy feasibility and optimality criteria.
II.3 Trajectory guidance 39 II.3.4.3 Advanced-step NMPC (AsNMPC) Consider that at any τi, i = 0, . . . , N −1, the current state vector is ~ Xand that ~u∗ iis applied to the system. If the model and associated parameters used during the trajectory optimisation process match with the reality, the system will evolve according to ~ F. Under these circumstances, at τiit is possible to predict with high accuracy the future state vector ˆ ~ X(τi+1)and, using the most recent prediction for ˆ ~ d, solve PN−i+1 in advance. If this problem can be solved in background during the discretisation interval [τi, τi+1], and under the strong assumption that the models match the reality, then ~u∗ i+1 will be already available when reaching τi+1. Accordingly, when using this guidance strategy, the discretisation interval between two consecutive time samples needs to be higher than the expected time required to solve the optimisation problem. Even if the execution time could be guessed from past experience, in some particular cases the optimisation may not finish before the next time sample. If this were the case, one of the following three contingency actions could be implemented: first, to stop the optimisation and not update the optimal trajectory at the current time sample; second, to wait until the optimisation finishes, and then shift the remaining time sample in such a way that the separation between them keeps unchanged; third, to wait until the optimisation finishes, but keeping the remaining time samples at their original values. This simple approach allows to remove the computational delay and associated issues of the INMPC. In the presence of errors in the models, however, the real state vector at τiwill not coincide with that predicted at τi−1. In this case, the optimal control cannot be computed in advance with accuracy. In order to account for errors in the models used by the FMS planning function, a fast update using parametric sensitivities can be used when measuring the actual state and parameters at τi. This strategy, which mixes an advanced full recalculation in background with an on-line sensitivity-based update when receiving measurements, is shown in Algorithm (II.3) (J¨ aschke et al.,2014;Zavala & Biegler,2009). Algorithm II.3: Working principle of the advanced-step NMPC (AsNMPC) 1: ~z∗,~ λ∗, ~µ∗← PN(~p) 2: Build ~z∗ 1:,~ λ∗ 1: and ~µ∗ 1: from ~z∗,~ λ∗and ~µ∗ 3: for i= 1, . . . , N −1do 4: Predict ˆ ~ Xand ˆ ~ d 5: ˆ ~p ← ˆ ~ X ˆ ~ d 6: ~z∗ i:,~ λ∗ i:, ~µ∗ i:← PN−iˆ ~pusing warm start 7: Evaluate F∗ i:at the new (~z∗ i:,~ λ∗ i:, ~µ∗ i:) and ˆ ~p 8: Measure ~ Xand estimate ~ d 9: ~p ←"~ X ~ d# 10: ∆~p ←~p −ˆ ~p 11: ∆~z∗ i:,~ λ∗ i:, ~µ∗ i:← QN−i(∆~p) 12: ~z∗ i:←~z∗ i:+ ∆~z∗ i: 13: Same as steps 9 to 20 of Algorithm (II.2)
You may delay, but time will not. —Benjamin Franklin I only feel angry when I see waste. When I see people throwing away things we could use. —Mother Teresa III The energy-neutral time window Arrival managers (AMANs) were designed to provide support for air traffic control (ATC) when sequencing and merging arrival traffic in the terminal maneouvering area (TMA). Existing AMANs use groundbased trajectory predictors to estimate the time of arrival of each inbound flight. These trajectory predictions are used to provide sequence information to the ATC, and also advisories of the required time to loose or gain for each flight in order to achieve the desired landing rate and spacing at the runway. The ATC is then responsible for selecting and notifying appropriate instructions (e.g., vectoring, speed changes) for the aircraft to loose or gain that time. Future AMANs will also determine the optimum landing sequence based on inbound flight plans and trajectory predictions. Then, a controlled time of arrival (CTA) will be calculated for each inbound aircraft at a specific metering fix such that the required landing rate and spacing are satisfied. When executed with sufficient timing accuracy, the CTAs will help to achieve an orderly arrival sequence of optimum descent profiles without excessive workload to the ATC. In order to efficiently allocate the CTAs, AMANs should know the earliest and latest time of arrival of each inbound flight (i.e., the feasible time window). An estimation of the feasible time window could be computed by the same ground system based on trajectory predictions. The models used by the groundbased trajectory predictor to represent the aircraft dynamics and performance, however, are typically subject to several simplifications, which limit the accuracy of the trajectory predictions. Furthermore, even if the models used by the ground-based trajectory predictor were perfect, the lack of flight intent information would make difficult to accurately predict the trajectory of the aircraft. A novel concept for synchronisation of ground-based trajectory predictors and flight management systems (FMSs) was proposed by Bronsvoort (2014), in which FMSs provide to the ground system flight intent information to improve their predictions. An alternative solution consists of using the on-board trajectory planner of the aircraft, which obviously knows the flight intent and has access accurate aircraft dynamics and performance models, to compute both the earliest and latest times of arrival at the metering and report them to the AMAN. Having this valuable information for all aircraft, the AMAN could compute the required CTAs (SESAR,2015b). 41
42 Chapter III - The energy-neutral time window The earliest and latest times of arrival depend on many factors, including aircraft performance, weather conditions and, noticeably, the current aircraft states when computing it (i.e., altitude, speed, mass and remaining distance to the metering fix). If the aircraft were in the cruise phase, the time of arrival at the metering fix could be easily modified by advancing or receding the position of the top of descent (TOD). Therefore, wide feasible time windows are expected for aircraft that have not initiated their descent yet. In some circumstances, however, the aircraft might be required to adapt the time of arrival after initiating the descent (e.g., due to a late CTA update or a large time deviation). In this situation, time and energy management concepts might still be used to efficiently re-plan the trajectory and adjust the time of arrival to the requested CTA with minimum fuel consumption and speed brakes use. The extent for which energy management can be used to adjust the time of arrival in this situation depends on the feasible time window, which, to the best of our knowledge, has not been quantified before. In this chapter, the feasible time windows for energy-neutral trajectories are quantified, assuming that the descent has been already initiated when the time of arrival at the metering fix needs to be changed by the on-board trajectory planner. In order to accomplish that, the earliest and latest times of arrival at a metering fix are computed for a wide range of aircraft states. This optimisation is restricted in such a way that neither thrust higher than idle nor speed brakes use are allowed along the whole descent. In this way, the entire trajectory is subject to optimisation to minimise or maximise the flight time, and only elevator control is left to modify the speed profile and change the time of arrival of the energy-neutral trajectory. In addition, the CTA window sensitivity to different longitudinal winds, metering fix position and masses is also investigated. III.1 State of the art Several studies have dealt with the assignment of CTA and the quantification of the feasible CTA windows when the aircraft is still in cruise, well before reaching the TOD. See for instance the works of Takeichi (2017) and Park & Clarke (2013,2016), which computed CTA windows at a metering fix by allowing the aircraft to adjust the TOD position (i.e., the duration of the cruise phase) and the descent speed profile, while keeping the thrust at idle during the whole descent. With similar purposes, Takeichi & Inami (2010) also enabled the addition (resp. omission) of waypoints to stretch (resp. shorten) the flight path length in addition to adjust the TOD under a tailored arrivals (TAs) concept. However, in the cruise phase, energy must be added to the system by means of throttle in order to maintain constant altitude and speed. This leads to considerably more fuel consumption and emissions than during an energy-neutral descent. The robustness of energy-neutral descent trajectories facing late changes to the CTA during the descent was assessed by Lindsay et al. (2009). Aiming to minimise the environmental impact of the continuous descent operation (CDO), only elevator control was permitted to adjust the time at which the metering fix was reached, assuming rather simple Mach/callibrated airspeed (CAS) profiles and allowing a single (and instantaneous) modification of the scheduled speed. Furthermore, the assessment was performed for few initial aircraft states, and the method used to generate the energy-neutral descents could not ensure the optimality (in the strict mathematical sense) of the resulting trajectories. Another important limitation of previous works is that none of them (up to our best knowledge) took into account the remaining descent between the metering fix at which the CTA was assigned and the instrumental landing system (ILS) interception or the runway threshold. Adjusting the descent speed profile to minimise or maximise flight time may result in a change to the altitude at which the metering fix is overflown. If the energy of the aircraft at this fix were too low, additional thrust might be needed after overflying it, interrupting in this way the CDO. Similarly, if the energy were too high, it might be required to use speed brakes and/or to deploy high-lifting devices or the landing gear earlier than initially planned, leading to more airframe noise. Furthermore, most previous works used the base of aircraft data (BADA) v3 performance model, which has been reported to show significant limitations for accurate trajectory planning tasks in the TMA (Senzig et al.,2009;Hoekstra,2016). More sophisticated aircraft performance models (APMs) are needed in order to compute realistic descent profiles and to obtain accurate fuel consumption and flight time, in particular time window, figures.
III.2 Setup of the experiment 43 III.2 Setup of the experiment The goal for the study was to quantify the CTA window (latest arrival time - earliest arrival time) at a certain metering fix as a function of the aircraft state when receiving the CTA after the TOD, and to study the sensitivity of the feasible CTA window to the following parameters: initial speed when receiving the CTA, aircraft mass, position of the metering fix along the route and longitudinal wind speed. In order to accomplish that, thousands of optimal trajectory plans starting at different aircraft states and configured with different values of the listed parameters were generated to minimise and maximise the time of arrival at a single metering fix. The optimisation of these aircraft trajectories was formulated as an optimal control problem as described in Section II.2.1, and was solved by using the numerical optimisation techniques described in Section II.2.3. This section describes the setup of the experiment performed in this chapter. First, Section III.2.1 presents the scenario. Then, the generic models introduced in Section II.1 are particularised for this experiment in Section III.2.2. Analogously, the generic optimal control problem described in Section II.2.1 is particularised in Section III.2.3 with specific cost function and constraints. Finally, the selected case studies are listed in Section III.2.4. III.2.1 Scenario Considering a specific aircraft type, descent procedure with associated constraints, geographical location and realistic weather conditions corresponding to a certain day and time only gives a clue of the feasible CTA window for that particular scenario, and results might not be representative of other scenarios. Since computing the CTA time window for all possible scenarios would be impracticable in combinatorial terms, this study aimed to be as generic as possible. Based on the discussion above, the following generic scenario was considered in this experiment: an Airbus A320 (one of the mostly used commercial aircraft types nowadays in Europe) attempts to perform an energy-neutral descent along an hypothetical straight-line route from the TOD to the interception of the ILS glide slope at the final approach point (FAP), and somewhere in the descent, before reaching the initial approach fix (IAF), the ATC assigns (or updates) a CTA at this metering fix. The typical speed constraint of 250 kt at FL100, enforced nowadays in most of the TMAs, was also considered in the formulation of the optimal control problem. Figure III-1 shows a schematic representation of this scenario. CTA Figure III-1: Scenario to investigate the feasible time window of energy-neutral descents III.2.2 Models The models required to formulate the optimisation of the aircraft trajectory, which were presented in Section II.1, are particularised below for the experiment conducted in this chapter.
44 Chapter III - The energy-neutral time window III.2.2.1 Aircraft dynamics model The point-mass model described by Eq. (II.1) was reduced to a γ-command model, which could be considered accurate enough for the objectives of the experiment. The γ-command model assumes continuous vertical equilibrium thorough the flight (i.e., lift balances weight): ˙γ= 0 →L=mg cos γ. (III.1) As a result, Eq. (II.1b) is removed from Eq. (II.1) and the aerodynamic flight path angle becomes an input control variable that can change instantaneously. The following assumptions were also considered: • The effect of the cross-wind on the ground speed is orders of magnitude below that of the longitudinal wind (i.e., wx= 0). • The longitudinal wind is a function of only the altitude (i.e., ws=ws(h)). • The fuel consumption during an energy-neutral descent is a very small fraction of the total mass (Clarke et al.,2004) and, consequently, the variations in mass mcan be neglected (i.e., ˙m= 0). In spite of that, the idle fuel flow qidle was computed in the model to determine the amount of fuel consumption for each trajectory, which was used as a metric during the assessment of results. Under the aforementioned assumptions, Eq. (II.1) can be simplified to: dv dt= ˙v=Tidle(v, h)−D(v, h, γ, ζ) m−gsin γ ds dt= ˙s=vcos γ+ws(h) dh dt=˙ h=vsin γ. (III.2) In this widely known γ-command model, the state vector is ~x = [v, s, h]T. Furthermore, since neither additional thrust nor speed brakes use is permitted, the control vector is composed by only the aerodynamic flight path angle (i.e., ~u = [γ]), and the thrust in Eq. (II.1a) is replaced by Tidle, which mathematical expression, function of ~x, is particularised by the APM described below. III.2.2.2 Aircraft performance model Aerodynamic and propulsive data from the Airbus performance engineering program (PEP) were taken to approximate the drag coefficient, idle thrust, and idle fuel flow tables with continuous and twicedifferentiable functions suitable for gradient-based nonlinear programming (NLP) optimisation algorithms, as suggested in Section II.2.4.2. Using these accurate data from the manufacturer allowed us to obtain realistic fuel and flight time figures. The drag coefficient was approximated by a smoothing spline of the form: CD(CL, M, h) = X i,j,k ci,j,kBi(CL)Bj(M)Bk(h).(III.3) Assuming temperature changes only with altitude, the Mach number can be expressed as a function of the altitude and speed: M(v, h) = v a(h)=v pγaRθ(h)τSSL ,(III.4) where ais the speed of sound. In Eq. (III.3), B(·)(with (·)∈ {i, j, k}) are spline basis functions; ci,j,k are control points of the spline (de Boor,1972); and CLis the lift coefficient, which can be obtained by assuming continuous vertical equilibrium (see Eq. (III.1)) and solving Eq. (II.2a) for CL:
III.2 Setup of the experiment 45 CL(v, h, γ) = 2mg cos γ SρSSLσ(h)v2.(III.5) Following the same methodology, Tidle and qidle were expressed as a function of hand M: Tidle (h, M) = X i,j ci,jBi(h)Bj(M),(III.6a) qidle (h, M) = X i,j ci,jBi(h)Bj(M).(III.6b) The control points and knots of the splines appearing in Eqs. (III.3) and (III.6) were obtained by fitting the aircraft performance data from the manufacturer using the FITPACK library (Dierckx,1993). Figure III-2 shows the smoothing splines functions approximating the aircraft performance data obtained from PEP. Altitude [FL] 050 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 qidle [kg h−1] 200 250 300 350 400 450 (a) Tidle Altitude [FL] 0 50 100 150 200 250 300 350 400 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 Tidle [kN] −15 −10 −5 0 5 (b) qidle CL[-] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 CD[-] 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 (c) CD(for h= 0) Figure III-2: Smoothing spline approximating aircraft performance data obtained from PEP
46 Chapter III - The energy-neutral time window III.2.2.3 Weather model This study aimed to be as generic as possible. Accordingly, the international standard atmosphere (ISA) model, established to provide a common reference for temperature and pressure, was considered to represent the normalised temperature θand pressure δas a function of hthrough Eqs. (II.4) and (II.5), respectively. Then, the normalised density was computed assuming the perfect gas law relationship with Eq. (II.6). With the same rationale, the empirical Hellman model defined by Eq. (II.7) was used to represent the wind speed as a function of h. The parameters of this empirical model were particularised as follows: α=1/7, which is the typical value for neutral stability conditions; hrwas set to 36,000 ft; and different values of wrwere analysed to assess the sensitivity of the feasible time window to the longitudinal wind speed (see Section III.2.4 for the detailed list of case studies). In addition, the wind vector was assumed to point in the direction of the aircraft track (i.e., ws=|~w|). III.2.3 Optimal control problem formulation The optimisation process adopted to generate earliest and latest trajectories at the metering fix is a constrained, multi-phase and nonlinear optimal control problem, as described in Section II.2.1. This section particularises the cost function and constraints of the generic multi-phase optimal control problem Eq. (II.10) for each one of the phases composing the flight profile of an energy-neutral CDO. The differential constraints (~ f) were particularised by the γ-command model given by Eq. (III.2): ~ fj= Tidle(v,h)−D(v,h,γ,ζj) m−gsin γ vcos γ+ws(h) vsin γ ;j= 0, . . . , P −1,(III.7) where ζjis the flaps/slats configuration of the jth phase. In this experiment, ζjwas set to clean (i.e., no flaps/slats deployed) for all phases except for the last one, which modelled the approach. In order to mimic a representative yet generic operation, the descent was split into P= 4 phases, where different phase-dependent path, algebraic and interior-point constraints may apply. Table III-1 wraps up the different phases and their constraints. This table is illustrated in Fig. III-3. Table III-1: Definition of the phases and constraints for an energy-neutral descent Phase (j) Description ~ bin j~ beq j~ ϑin j~ ϑeq j 0 Descent above FL100 "vCAS(v, h)−VMO GD −vCAS(v, h)#- - hh−10,000 fti 1 Descent below FL100 "GD −vCAS(v, h) vCAS(v, h)−250 kt#- - hs−sIAFi 2 Approach "GD −vCAS(v, h) vCAS(v, h)−250 kt#- - "vCAS (v, h)−GD h−hFAP # 3 Levelled deceleration "vFAP −vCAS (v, h) vCAS(v, h)−GD#h˙ hi- - All Entire descent M(v, h)−MMO γ γmin −γ - - - Several path constraints applying all along the descent were also considered in order to ensure that the Mach number remains within operational limits and that the maximum descent gradient is not exceeded. The maximum operative CAS (VMO) and the maximum operative Mach (MMO) are constant parameters that depend on the aircraft type. In this study, MMO and VMO were set to 0.80 and 340 kt, respectively, according to Airbus (1993). Furthermore, vFAP corresponds to the CAS that the aircraft should have when intercepting the ILS glide slope at the FAP, and was computed as a fraction of the stall speed, which is a function of both mand ζ(Airbus,1993). The maximum descent gradient (γmin) was set to -15◦.
III.2 Setup of the experiment 47 Phase 0 1 2 3 Figure III-3: Definition of the phases and constraints for an energy-neutral descent According to Table III-1, the first phase starts at the current state of the aircraft ( ~ X) and ends at FL100. Note that since the goal of this study was to quantify the CTA window for different initial conditions, the optimal control problem was solved for several values of ~ X. Below FL100, ATC procedures typically restrict the CAS to 250 kt. Afterwards, the aircraft heads towards the IAF, where the approach procedure begins. Few nautical miles before the FAP, the aircraft decelerates to the green dot speed (GD)1and starts configuring with flaps at constant altitude (still at idle thrust), in order to intercept the ILS glide slope at the FAP. Therefore, terminal constraints were: ~ ψ= vCAS(v, h)−vFAP h−hFAP s−sFAP ,(III.8) where hFAP is the altitude at which the FAP has to be intercepted, which value was set to 2,000 ft; and the distance of the FAP (sFAP) was fixed at 6 NM from the runway threshold (see Fig. III-1). The objective of this study was to compute trajectories minimising and maximising the time of arrival at a given metering fix, which was assumed to be the IAF of an hypothetical straight-line procedure. Accordingly, the running cost function was defined as: πj=(±1 vcos γ+ws(h)if j= 0,1 0otherwise,(III.9) where the +sign applied when computing the earliest trajectory, and −when computing the latest. The end cost for each phase was not considered, that is φj= 0 for j= 0, . . . , P −1. Furthermore, since the duration of the time horizon, tF, was not fixed but had to be optimised, tj+1 for j= 0, . . . , P −1were additional decision variables of the optimal control problem. The optimal control particularised above was transformed to a parametric NLP optimisation problem as described in Section II.2.3. Then, the resulting parametric NLP optimisation problem was formulated in the general algebraic modelling system (GAMS) software suit, and solved with CONOPT as NLP solver, which uses active-set methods. In this experiment, the trapezoidal integration scheme was adopted to obtain ~ Fjand Πjfrom ~ fjand πj, respectively. In addition, techniques to recover from infeasibility were not used, i.e., all the constraints shown in Table III-1 are hard constraints. Last but not least, the profile shown in Table III-1 and Fig. III-3 does not include speed and altitude constraints throughout the descent (excepting the typical speed limit at FL100), aiming to be as generic as possible. Consequently, it is rather optimistic and the obtained results correspond to a best-case scenario. In current operations, aircraft may encounter additional altitude and speed constraints, which negative effect on the feasible CTA window would depend on their location and severity. It might even be the case that it would not be longer possible to reach the metering fix using neither additional thrust nor speed brakes. Nevertheless, it is expected that future procedures will trade altitude and speed constraints for CTAs, giving more flexibility to FMSs for planning the optimal trajectory provided that time CTAs are satisfied. 1For the Airbus A320, the green dot speed is the minimum operating speed in managed mode and clean configuration, being approximately the best lift-to-drag ratio speed.
54 Chapter III - The energy-neutral time window of the mass), slower descents can be achieved with lighter weights because of the lower GD speeds. Consequently, there is a trade-off between maximising the feasible time window and maximising the amount of payload carried on-board the aircraft. It is worth noting that the difference of feasible time window between two different aircraft masses becomes more significant at farthest distances from the runway threshold. For instance, for an increment of the aircraft mass from 80% to 100% of MLM, at 60 NM from the runway threshold the maximum time window decreases 30 s, while at 100 NM the detriment is about 1 min. III.3.3.2 Effect of the longitudinal wind Figure III-10 shows the effect of the longitudinal wind on the CTA window in a matrix form. Each column represents a different magnitude of the wind speed. The first row shows the results for the simulations with head wind, while the second row for those with tail wind. It is important to recall that the wind speed is not constant, but changes with the altitude according to Eq. (II.7). Therefore, the values of wind speed shown in Fig. III-10 correspond to wr, not to ws. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 100 150 200 250 300 350 Altitude [FL] Infeasible region (excessive energy) Infeasible region (insufficient energy) Tail wind 0.4 0.4 0.8 1.3 1.7 2.1 Wr= 20 kt 0.4 0.8 1.2 1.6 2.0 Wr= 40 kt 0.4 0.7 1.1 1.5 1.9 Wr= 60 kt 40 60 80 100 120 Distance to go [NM] 100 150 200 250 300 350 Altitude [FL] Head wind 0.5 0.9 1.4 1.9 2.4 40 60 80 100 120 Distance to go [NM] 0.5 1.0 1.5 2.1 2.6 40 60 80 100 120 Distance to go [NM] 0.6 0.6 1.1 1.7 2.3 2.8 Maximum time window trajectory 0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 Time window [min] Figure 8 Figure III-10: Effects of longitudinal wind on the feasible time window at the metering fix According to Fig. III-10, in presence of head wind the feasible range of initial aircraft states moves towards the time-constrained metering fix and becomes narrower, but the feasible energy-neutral time window at the IAF increases, if compared with calm wind conditions (see Fig. III-6). The opposite behaviour applies in presence of tail wind. This observation is easy to understand by analysing the impact of the wson the total energy lost per unit length, which expression can be deducted by combining Eqs. (III.2) and (II.29): dEt ds=dEt dt dt ds=˙ Et ˙s=v(Tidle −D) vcos γ+ws .(III.17) Equation (III.17) leads to the following inequalities for a given weight, altitude and speed: dEt dsws>0<dEt dsws=0 <dEt dsws<0.(III.18) According to Eq. (III.18), in presence of head (resp. tail) wind, less (resp. more) distance is needed to release the same amount of energy if compared with calm winds condition. Namely, for a given initial
III.3 Results 55 energy, in presence of head wind the descent must start closer to the metering fix to satisfy all the constraints. Equation (III.18) explains why, for a given altitude, the set of initial states moves towards the IAF and becomes narrower for head wind conditions. Generally speaking, for the earliest trajectories it would be desirable to maintain the maximum allowed speed (either MMO/VMO or 250 kt, depending on the altitude) down to the IAF. Similarly, for the latest trajectories, the deceleration to GD should be performed as soon as possible. In practice, however, the capability to release energy once the IAF is overflown constrains the altitude and speed at which this fix can be reached, and the ideal conditions cannot always be met. For instance, for the earliest energy-neutral trajectories in the case studies with head wind, the optimal solution consists of reaching FL100 at VMO and decelerating to 250 kt just before the IAF, since the aircraft is capable of releasing such amount of energy in the remaining descent distance. In presence of tail wind, however, the energy at the IAF cannot be as high because of the lower |dEt/ds|, and the deceleration to 250 kt must be performed well before the metering fix. III.3.3.3 Effect of the IAF position Figure III-11 shows the effect of the position of the metering fix on the feasible CTA window. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 40 60 80 100 120 Distance to go [NM] 100 150 200 250 300 350 Altitude [FL] Infeasible region (excessive energy) Infeasible region (insufficient energy) 0.5 1.0 1.5 2.1 2.6 sIAF = 15 NM Maximum time window trajectory 60 70 80 90 100 110 120 130 Distance to go [NM] 0.3 0.7 1.0 1.3 1.7 sIAF = 50 NM 0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 Time window [min] Figure 9 Figure III-11: Effects of the position of the metering fix on the feasible time window As expected, the closer the metering fix to the runway threshold, the wider the CTA window because more distance is left to gain or lose time. It is interesting to note that, even if the time window is significantly affected by the position of the metering fix, the set of feasible initial conditions remains unchanged. This is because of the particularities of the selected scenario, in which neither speed nor altitude constraints have been considered. As a result, the set of feasible initial conditions is not affected by the IAF location, but only the value of the cost function. III.3.3.4 Effect of the initial speed Results shown in the preceding sections were computed fixing the initial speed of the aircraft to that corresponding to a kinetic energy in the middle of the feasible kinetic energy range. However, the initial kinetic energy of the aircraft plays an important role in the feasible time window, which must be also investigated. For instance, if the CTA update forces the aircraft to arrive earlier than initially planned while flying very slow, it might be not possible to accelerate enough without adding energy by means of thrust. Conversely, if the aircraft is forced to arrive later than initially planned while flying at a very high speed, it might be infeasible to decelerate (thus decrease the drag) and at the same time release the required energy without using speed brakes Figure III-12 shows the effects of the initial speed on the feasible energy-neutral time window. According to this Fig. III-12, flying at a very high speed when receiving a CTA update is typically preferred. For a given altitude and distance to go, the higher the speed, the higher the total energy. A high energy level leads to more freedom to either release or retain it by increasing (provided that the aircraft is not flying
56 Chapter III - The energy-neutral time window 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 40 60 80 100 120 Distance to go [NM] 100 150 200 250 300 350 Altitude [FL] 0.7 1.4 2.1 2.8 3.5 Free initial v 40 60 80 100 120 Distance to go [NM] Infeasible region (excessive energy) Infeasible region (insufficient energy) 0.4 0.4 0.9 1.4 1.9 2.3 Maximum initial v 40 60 80 100 120 Distance to go [NM] 0.4 0.8 1.2 1.6 2.1 Minimum initial v Maximum time window trajectory 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Time window [min] Figure 10 Figure III-12: Effects of the initial speed on the feasible time window at the metering fix at VMO/MMO/250 kt) or decreasing the speed, respectively. If the initial energy is limited, however, the only choice left consists of preserving it as long as possible to satisfy the terminal condition, not being able to fly faster and, consequently, highly penalising the earliest energy-neutral time of arrival. Results for the case with a free initial airspeed (i.e., allowing the aircraft to change its initial speed instantaneously) show feasible energy-neutral time windows up to 4 min, in contrast to the 2.8 min shown for the baseline case study. However, changing the initial speed implies adding or removing energy from the system by means of additional thrust or drag, respectively, leading to extra fuel consumption and/or noise nuisance. III.4 Discussion An efficient assignment of CTAs when sequencing and merging traffic could enable the introduction of CDOs in the TMA, where conflictive trade-offs exist between airspace capacity and environmental objectives. For this purpose, future ground automation systems will require a good knowledge of the feasible CTA windows for each inbound aircraft. This study quantified this feasible time window at the IAF for an Airbus A320 performing a CDO such that the speed profile is only adjusted by means of elevator control along the descent, and assuming the descent has been already initiated. Results for a baseline case study show that CTA windows up to 2.8 min at the IAF could be achieved for certain flight conditions. Another important remark that arises from this study is that a previously assigned CTA or even the CI could have a significant impact on the CTA window. Interestingly, it has been observed that minimum time (i.e., faster) trajectories are very similar to those of minimum fuel, provided that the optimisation takes place once in the descent (i.e., the top of descent has been overflown) and neither additional thrust nor drag devices usage are allowed. It has been also found that there is a trade-off between maximising the robustness of facing late CTA updates throughout the descent (i.e., adapting the time of arrival) and minimising the fuel consumption. Therefore, for robust and optimal CDOs in the presence of uncertainties, it would be more interesting to fly a higher extra fuel consumption if compared to that minimum achievable. Moreover, it has been observed that the lighter the aircraft, the wider the feasible time window. In addition, in presence of head wind, the set of feasible flight conditions such that an energy-neutral trajectory can be accomplished, narrows and moves towards the time-constrained metering fix, but the magnitude of the feasible time window increases. Results also show that receiving the CTA update while flying at a high speed is typically preferred. Finally, and as expected, the closer the aircraft to the metering fix when receiving the CTA update, the narrower the feasible time window The most important conclusion of this chapter is that CTAs can be updated after overflying the TOD, and that large energy-neutral time windows could still be attainable depending on the flight conditions.
He who is not satisfied with a little, is satisfied with nothing. —Epicurus The aim of science is not to open the door to infinite wisdom, but to set a limit to infinite error. —Bertolt Brecht IV Comparison of guidance strategies to meet time constraints in optimal descents In Chapter I, the assignment of time constraints to each inbound aircraft at some strategic fixes for separation, sequencing and merging tasks was identified as a potential approach to enable continuous descent operations (CDOs) in dense traffic scenarios. With this type of flight operations, air traffic control (ATC) would assign (using data-link, for instance) at least one controlled time of arrival (CTA) to each aircraft, ideally before they reach the top of descent (TOD). Chapter III quantified the feasible time windows of energy-neutral CDOs that have overflown their TOD, showing that time windows up to 2.8 minutes could still be achieved for certain flight conditions. After receiving the CTA from the ATC, either before of after the TOD, the time instruction would be entered as a required time of arrival (RTA) by the aircraft crew into the flight management system (FMS). Right after, the on-board trajectory planner would compute the optimal trajectory plan for the descent satisfying the RTA and other applicable operational constraints. In practise, however, because of modelling errors and uncertainties encountered during the execution of the trajectory plan, aircraft may deviate from the planned altitude and/or speed profile, and may also be delayed or advanced with respect to the planned time schedule. In order to safely execute the trajectory plan as well as to accurately comply with the entered RTA in such uncertain environment, corrective actions must be taken by the on-board guidance system of the aircraft. Several guidance strategies were proposed as a part of the time and energy managed operations (TEMO) concept for this purpose, as mentioned in Section I.4.2.2. The aim of this chapter is to compare the performance of the three existing TEMO guidance strategies described in Section II.3 (i.e., strategic, tactical and hybrid), in terms of noise nuisance, energy error at the metering fix, fuel consumption, speedbrakes use, as well as their ability to satisfy time constraints in presence of modelling errors. Furthermore, a new TEMO guidance strategy called nonlinear model predictive control (NMPC), which is one of the major contributions of this PhD thesis, will also take part of this comparison. 57
58 Chapter IV - Comparison of guidance strategies to meet time constraints in optimal descents IV.1 State of the art Various guidance strategies have been proposed to comply with a CTA during a CDO, which essentially differ on how the aircraft actuators are used to nullify the time deviations with respect to the trajectory plan, and whether this trajectory plan is dynamically adjusted or not. For instance, the speed profile could be continuously modified to tactically nullify a time error signal by using throttle and speed-brakes commands, while remaining on the vertical flight path with the elevator, as suggested by Rumbo et al. (2002). At present, however, the automatic control of speed brakes is seldom, implying that an additional margin must be available to control the throttle. This margin may be acquired, for instance, by generating a trajectory plan at idle thrust plus some arbitrary increment. For constant speed segments, the resulting flight path angle will be less steep, while for constant flight path angle decelerations, the resulting decelerations will be slower. The continuous action on thrust required by this guidance strategy, however, has a negative effect on noise nuisance, gaseous emissions, fuel usage and engine wear and tear. A more appealing guidance strategy, in terms of environmental impact, consists of actively controlling the speed with the elevator (speed-on-elevator) in order to tactically nullify time deviations, then using near-idle thrust variations to recover the initially planned path only when the vertical deviation exceeds a predefined threshold (Garrido-L´ opez et al.,2009). A set of additional guidance strategies were proposed as a part of the TEMO concept (De Jong et al., 2017,2015). These strategies have in common to manage the total energy of the aircraft (altitude and speed) such that time constraints are satisfied in the most environmentally friendly way. For instance, if the aircraft is late with respect to the CTA but exceeds the planned energy, the guidance system could command a steeper descent angle to accelerate, which would result in an increased energy loss rate while simultaneously reducing the time error. These guidances were described in Section II.3, and will be summarised below for the sake of completeness. One guidance option in TEMO is to use speed-on-elevator to nullify time errors, while throttle and speed brakes are continuously used to correct energy deviations (tactical guidance). Differently, the strategic guidance approach consists of strictly follow the trajectory plan using speed-on-elevator while allowing certain energy and time deviations during the execution, as long as some predefined thresholds are not exceeded. Whenever time and/or energy errors exceed these thresholds, a new trajectory plan can be computed from the current aircraft state in such a way that corrective actions are taken to satisfy all operational constraints with minimum fuel usage (De Jong et al.,2017). In previous work, this guidance technique was compared with the guidance strategy of a typical FMS that does not correct time errors, showing superior performances in terms of fuel consumption and RTA compliance in presence of wind errors (Dalmau et al.,2016). The combination of the strategic and tactical guidance variants leads to the hybrid strategy, where only the thrust plan is followed and a tactical controller is in charge to nullify time errors by using speed-on-elevator (De Jong et al.,2015). Energy deviations exceeding the allowed threshold are corrected by means of strategic updates of the trajectory plan. Finally, with NMPC (Diehl et al.,2008), the remaining descent trajectory to the metering fix is frequently updated by solving an optimal control problem, setting the current state of the aircraft as initial conditions. Different from the previous guidances strategies, the NMPC does not wait for energy and time deviations to exceed a certain threshold, nor nullifies them to strictly follow the initial trajectory plan. Thus, it could be seen as the limit of the strategic guidance, with energy and time deviation thresholds tending to zero, in which the optimal plan is tactically updated. Recent research demonstrated the feasibility of using NMPC guidance strategies to achieve precise spacing between aircraft, the main objective of interval management operations (Weitz & Bai,2018). However, due to limitations in the execution time of the trajectory optimisation algorithm, previous work implemented a receding NMPC that did not update the remaining descent trajectory at each time sample, but only computed the optimal controls over a relatively short future horizon such that the initial trajectory plan was closely followed in the most optimal way. Other relevant contributions of NMPC in the aerospace field have been performed, for instance, by Arrieta-Camacho & Biegler (2005), Gibbens & B. Medagoda (2011) and Eren et al. (2017).
IV.2 Setup of the experiment 59 IV.2 Setup of the experiment The comparison of the four TEMO strategies was performed by simulating 123 descents per guidance strategy (a total of 492) in the avionics prototyping environment for research and operations (APERO), a flight simulator located at the Netherlands Aerospace Centre (NLR). The APERO is composed by a research FMS, which allows the development testing of modules targeting applications for the optimisation of trajectories, and a high-fidelity simulator of the aircraft dynamics, which implements a six-degree-of-freedom model fed with accurate aircraft performance data. A trajectory management software developed at the Technical University of Catalonia (UPC) as a result of this PhD thesis was connected to the research FMS of the APERO, as illustrated in Fig. IV-1. This trajectory management software is composed by a trajectory planner and a variety of guidance algorithms build on time and energy management concepts. NLR High-fidelity 6DOF flight simulator Research FMS (RFMS) Trajectory planner Guidance strategy Aircraft performance model ISA Re-plan trigger Commands for the actuators Polynomial fitting Figure IV-1: Setup of the experiment to compare guidance strategies Basically, the trajectory planner was in charge of computing the initial trajectory plan complying with a CTA at a certain metering fix of a realistic procedure. The optimisation of the aircraft trajectory was formulated as an optimal control problem (see Section II.2.1), and solved by using the optimisation techniques described in Section II.2.3. Thereafter, the trajectory plan was executed with the specific guidance strategy subject of study. Using its particular logics, the guidance system provided the optimal controls to the research FMS, which was in charge of computing the commands for the actuators. The trajectory planner updated the optimal trajectory plan when requested by the guidance system, based on the state observations provided by the research FMS. As shown in Fig. IV-1, the simulations were performed by introducing intentional discrepancies between the weather forecast and aircraft performance models implemented in the trajectory planner, and the corresponding models implemented in the simulation environment. The detailed list of case studies and their respective modelling errors will be presented in IV.2.4. Fuel consumption, time and energy deviations at the metering fix, the number of speed brakes deployments, and area affected by the aircraft noise were selected as metrics for this comparison. Noise contours were calculated using the model described by Van der Wal et al. (2001). This section describes the setup of the experiment performed in this study. First, the scenario is presented in Section IV.2.1. Then, the generic models introduced in Section II.1 are particularised in Section IV.2.2. Analogously, the generic optimal control problem stated in Section II.2.1 is particularised in Section IV.2.3 with specific cost function and constraints. The case studies of the experiment are listed in Section IV.2.4, and the workflow of a single descent simulation is summarised in Section IV.2.5.
60 Chapter IV - Comparison of guidance strategies to meet time constraints in optimal descents IV.2.1 Scenario All simulations assumed the following scenario: an Airbus A320 was cruising at a certain altitude when the ATC requested a CTA at a metering fix. Right after entering the RTA in the FMS, the on-board trajectory planner computed the optimal trajectory plan and executed it by implementing any of the guidance strategies considered in the experiment. However, the trajectory plan was computed under certain weather forecast and aircraft performance modelling errors, and the guidance system reacted to correct energy and time using its specific mechanisms. All the simulations were performed to Barcelona-El Prat airport (LEBL), using four realistic standard arrival procedures (STARs) published for runway 25L. As can be observed in Fig. IV-2, each STAR ends at one of the four initial approach fixs (IAFs) of LEBL and then, all the approaches finally merge into a single intermediate fix (SOTIL). The instrumental landing system (ILS) glide slope is intercepted at the final approach point (FAP) few nautical miles later. For this experiment, SOTIL (the intermediate fix) was selected as the metering fix where CTAs were assigned by an hypothetical ATC for sequencing and merging tasks. Figure IV-2: Routes at Barcelona-El Prat chosen to compare guidance strategies (Source: Spanish AIP) IV.2.2 Models The models presented in Section II.1 are particularised below for the experiment conducted in this chapter. IV.2.2.1 Aircraft dynamics model The point-mass model described by Eq. (II.1) was implemented in the trajectory planner of this experiment. The control vector, however, was modified in the following way: the lift coefficient was exchanged for the load factor (nz), and the thrust for the throttle setting (κ); that is: ~u = [κ, nz, β]Tinstead of ~u = [T, CL, β]T. This also implies the following change in Eq. (II.1b): dγ dt= ˙γ=g v(nz−cos γ),(IV.1) where the load factor is defined as the ratio of the aerodynamic lift of an aircraft to its weigh: nz=L mg ,(IV.2) which represents a measure of the stress to which the structure of the aircraft is subject. According to Eq. (IV.2), a load factor of one corresponds to straight and level flight, where the lift is equal to the weight. Load factors different from one (or even negative) are the result of manoeuvres.
IV.2 Setup of the experiment 61 The throttle was assumed to control the revolutions of the engine fan (N1) using a linear relationship: N1(v, h, κ) = N1idle(v, h) + κ(N1max(v, h)−N1idle(v, h)) ,(IV.3) where N1idle and N1max represent, respectively, the minimum and maximum revolutions of the engine fan, which are both function of the altitude and airspeed of the aircraft. In this model, the thrust force and fuel flow were described as a function of N1, the altitude and the speed: T(v, h, κ) = T(v, h, N1(v, h, κ)) (IV.4a) q(v, h, κ) = q(v, h, N1(v, h, κ)).(IV.4b) The differential equations of Eq. (II.1) for the remaining states kept unchanged. IV.2.2.2 Aircraft performance model The aircraft performance model (APM) used in this experiment was based on the performance data provided in tabular form by NLR for the Airbus A320. Since the nonlinear programming (NLP) solver used in this study required first and second-order differentiability of the cost function and constraints, the drag coefficient was modelled as the following drag polar: CD(η, ζ, β, CL) = CD0(η, ζ) + CD2(η, ζ)C2 L+CDβ(η, ζ)β, (IV.5) which neglects the compressibility effects and consider the effects of the speed brakes deflection as an additional linear term. In Eq. (IV.5), CDβ is a coefficient that represents the increase in drag coefficient when the speed brakes are deployed, and CLcan be computed by combining Eq. (IV.2) with Eq. (II.2a): CL(nz, v, h, m) = 2nzσ(h)ρSSLSv2 mg .(IV.6) Figure IV-3 shows the drag coefficient as a function of CL, for different aerodynamic configurations. 0.0 0.5 1.0 1.5 CL[-] 0.00 0.05 0.10 0.15 CD[-] ζ= 0.0 ζ= 5.0 ζ= 10.0 ζ= 15.0 Figure IV-3: Drag polars approximating aerodynamic data for the A320 The maximum revolutions of the engine fan were set to a constant value, and the normalised revolutions of the engine at idle thrust were modelled as fifth order polynomial: N1max = 87 (IV.7a) N1idle(v, h) pθ(h)= 3 X i=0 3 X j=0 ci,jθi(h)Mj(v, h).(IV.7b)
62 Chapter IV - Comparison of guidance strategies to meet time constraints in optimal descents θ[-] 0.75 0.80 0.85 0.90 0.95 1.00 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 N1max/√θ[%] 88 90 92 94 96 98 100 (a) Normalised max. revolutions of the engine fan θ[-] 0.75 0.80 0.85 0.90 0.95 1.00 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 N1idle/√θ[%] 30 35 40 45 50 (b) Normalised min. revolutions of the engine fan N1/√θ[%] 30 40 50 60 70 80 90 100 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 T/δ[kN] 0 20 40 60 80 100 120 140 160 (c) Normalised thrust N1/√θ[%] 30 40 50 60 70 80 90 100 Mach [-] 0.3 0.4 0.5 0.6 0.7 0.8 q/δ√θ[kg h−1] 1000 2000 3000 4000 5000 6000 7000 8000 (d) Normalised fuel flow Figure IV-4: Polynomials approximating propulsive data for the A320 The same methodology was used to model the normalised fuel flow and the normalised thrust as a function of the normalised revolutions of the engine fan: q(v, h) δ(h)pθ(h)= 3 X i=0 3 X j=0 ci,j N1 pθ(h)!i Mj(v, h)(IV.8a) T δ(h)= 3 X i=0 3 X j=0 ci,j N1 pθ(h)!i Mj(v, h).(IV.8b) The coefficients of the polynomials functions appearing in Eqs. (IV.5), (IV.7) and (IV.8) where obtained as a result of a least-squares fitting process with the aerodynamic and propulsive tabular data used in the flight simulator, respectively. Figure IV-4 shows the polynomial approximation of the propulsive data. It should be noted that even if both flight simulator and trajectory planner draw upon the same aircraft type (i.e., an Airbus A320), the dynamics and performance models used by these components were sightly different. The flight simulator used as a high fidelity six-degree-of-freedom model to represent the aircraft dynamics, and very accurate look-up tables for the aircraft performance. Moreover, the trajectory planner
IV.2 Setup of the experiment 63 represented the aircraft as a three-degree-of-freedom point-mass model (similar to Eq. (II.1)), and approximated the aircraft performance data included in the look-up tables of the flight simulator with polynomial functions, aiming at achieving the continuity and second-order differentiability required by the embedded NLP solvers. IV.2.2.3 Weather model The trajectory planner used in this experiment assumed international standard atmosphere (ISA) conditions and calm winds. Accordingly, the normalised temperature θand pressure δwere modelled as a function of hthrough Eqs. (II.4) and (II.5), respectively; the normalised density was computed assuming the perfect gas law relationship with Eq. (II.6); and ws=wx= 0. IV.2.3 Optimal control problem formulation The trajectory planner used in this experiment generated the optimal descent trajectory subject to a RTA by solving a multi-phase, constrained and nonlinear optimal control problem, as described in Section II.2.1. This section particularises the cost function and the dynamic, path, algebraic and interior-point constraints for each one of the phases composing the flight profile. The differential constraints on the dynamics of the aircraft (~ f) were particularised by the following point-mass model, independently of the phase: ~ fj= T(v,h,κ)−D(v,h,m,nz,β) m−gsin γ g v(nz−cos γ) vsin γ vcos γ −q(v, h, κ) ;j= 0, . . . , P −1.(IV.9) The descent was divided into P= 4 different phases, where several phase-dependent path, algebraic and interior-point constraints may apply. Table IV-1 wraps up the different phases and their associated constraints. This table is also illustrated in Fig. IV-5. Table IV-1: Definition of the phases and constraints for the comparison of guidance strategies Phase (j) Description ~ bin j~ beq j~ ϑin j~ ϑeq j 0 Cruise phase "GD −vCAS(h, v) vCAS(h, v)−VMO# "˙ h(v, γ) ˙ M(h, v)#- - 1 Cruise deceleration "GD −vCAS (h, v) vCAS(h, v)−VMO#h˙ h(v, γ)i- - 2 Descent "GD −vCAS(h, v) vCAS(h, v)−VMO#- - hvCAS(h, v)−GDi 3 Approach "vIF −vCAS(h, v) vCAS(h, v)−GD#- - - All Entire descent M(h, v)−MMO κ κ−1.0 β β−1.0 nz−1.15 0.85 −nz - - - Phase-independent path constraints (i.e., applying all along the descent) as a function of the state and control vectors were also considered, which ensured that the Mach number remains within operational limits all along the descent, and that the bounds on the controls are not exceeded.
70 Chapter IV - Comparison of guidance strategies to meet time constraints in optimal descents F T N E S 175 178 181 184 187 0.4 0.8 1.2 1.6 2 254 256 258 260 262 15 30 45 60 75 2 4 6 8 10 MPC Tactical Hybrid Strategic NMPC (a) Aggregated results F T N E S 175 178 181 184 187 0.4 0.8 1.2 1.6 2 254 256 258 260 262 15 30 45 60 75 2 4 6 8 10 (b) Only wind forecast errors F T N E S 175 178 181 184 187 0.4 0.8 1.2 1.6 2 254 256 258 260 262 15 30 45 60 75 2 4 6 8 10 (c) Only temperature forecast errors F T N E S 175 178 181 184 187 0.4 0.8 1.2 1.6 2 254 256 258 260 262 15 30 45 60 75 2 4 6 8 10 (d) Only aircraft performance modelling errors Figure IV-9: Mean metrics per guidance strategy. E: Specific energy deviation [ft]; F: Fuel consumption [kg]; N: Area affected by aircraft noise [NM2]; S: Number of speed brakes deployments; T: Time deviation [s]. monitored, but the elevator, throttle and speed brakes continuously act together to nullify them. Yet, the corrective actions performed by the NMPC are mathematically optimal and take into account the remaining descent trajectory, which is dynamically adjusted. For this reason, the noise nuisance, fuel consumption and speed brakes use for the NMPC strategy are comparable with those of the strategic and hybrid strategies. The hybrid strategy show the best results in terms of environmental impact mitigation and speed brakes use, proving to exploit the advantages of the strategic and tactical mechanisms. However, the differences with respect to the strategic and NMPC strategies are not significant enough to conclude that the hybrid guidance could be preferred over the others. As expected, the specific energy deviation at the metering fix is higher for the strategic and hybrid strategies, which allow deviations up to 100 ft at this fix before triggering a re-plan. Similarly, the tactical strategy permits specific energy deviations up to 50 ft before applying throttle or deploying speed brakes. The NMPC shows the best results in terms of energy deviation, because every 30 seconds the energy deviation is nullified when activating the new trajectory plan. The fuel consumption shown in Fig. IV-9 does not consider the phases after the metering fix. If this fix
IV.3 Results 71 Hybrid NMPC Strategic Tactical −100 −50 0 50 100 150 Specific energy deviation at SOTIL [ft] (a) Specific energy deviation at the metering fix Hybrid NMPC Strategic Tactical −2 0 2 Time deviation at SOTIL [s] (b) Time deviation at the metering fix Hybrid NMPC Strategic Tactical 160 170 180 190 200 210 Fuel consumption [kg] (c) Fuel consumption Hybrid NMPC Strategic Tactical 255 260 265 270 275 280 Area affected by aircraft noise [NM2] (d) Area affected by aircraft noise Hybrid NMPC Strategic Tactical 0 5 10 15 20 Number of speed-brakes deployments (e) Number of speed brakes deployments Figure IV-10: Metrics distribution per guidance strategy were reached with a certain energy deviation, the missing or excessive energy would need to be added by means of additional thrust (leading to more fuel consumption and noise) or released by deploying speed brakes (leading to more airframe noise), respectively, at low altitudes and near populated areas. With NMPC guidance, the aircraft typically reaches the metering fix with negligible energy deviations. Consequently, the NMPC does not require the addition nor subtraction of energy after the metering fix. As expected, and according to Fig. IV-9(b), wind errors have the effect of increasing both absolute time and energy deviations. It is also worth to note that, at aggregated level (and considering the same number of case studies with prevalent head wind than with tail wind), the fuel consumption, noise nuisance and speed brakes use are similar to those observed in Fig. IV-9(a). In addition, Fig. IV-9(c) shows that temperature deviations also lead to significant time and energy deviations. On the one hand, all the guidance strategies considered herein execute the CAS plan, leading to true airspeed (TAS) errors in presence of temperature deviations. On the other hand, the aircraft performance strongly depends on the air temperature (see Eq. (IV.8b)). Finally, according to Fig. IV-9(d), errors in the aircraft performance models are not penalising
72 Chapter IV - Comparison of guidance strategies to meet time constraints in optimal descents 10 0 10 ¯ Ws[kt] 170 180 190 200 Fuel Consumption [kg] (a) Fuel consumption 10 0 10 ¯ Ws[kt] 2 0 2 Time Deviation [s] SHMPC hybrid strategic tactical NMPC (b) Time deviation 10 0 10 ¯ Ws[kt] 200 100 0 100 Specific Energy Deviation [ft] (c) Specific energy deviation 10 0 10 ¯ Ws[kt] 0.0 2.5 5.0 7.5 10.0 12.5 Number of Speed-brakes Deployments [-] (d) speed brakes use 10 5 0 5 10 ¯ Ws[kt] 32 34 36 38 Area Afected by the Noise [NM2] (e) Area affected by the noise 10 5 0 5 10 ¯ Ws[kt] 0 5 10 15 20 25 30 Number of Replans [-] (f) Number of re-plans Figure IV-11: Correlation between metrics and longitudinal wind speed error the capability to satisfy CTAs, but have a major impact on the environment. Figure IV-9 showed only the mean of the different metrics per guidance strategy. Figure IV-10 shows the complete distribution for each metric and guidance strategy. Essentially, the conclusions remain the same than stated before. Namely, the strategic and hybrid show the larger dispersion of energy deviations at the metering fix, due to their passive behaviour with respect to the energy deviation as long as it remains within the allowed threshold; the tactical hybrid and NMPC show excellent figures in terms of time error at the metering fix, while the time error for the strategic approach shows larger dispersion; and the tactical strategy is the one leading to more fuel consumption, noise nuisance and speed-brakes use. IV.3.3 Sensitivity to the wind error Figure IV-11 shows, for each guidance strategy, the dependence of each particular metric with respect to the average longitudinal wind speed error encountered druing the descent. Each point in the scatter plot shows the result of a single simulation. Positive (resp. negative) values of longitudinal wind speed error represent unexpected head (resp. tail) wind encountered during the execution of the descent. According to Fig. IV-11, and as expected, the higher the unexpected head wind, the larger the fuel consumption and noise nuisance. In contrast, the number of speed brakes deployments typically increases with the unexpected tail wind. The rationale behind these results was already discussed in Section III.3.3.2: the head wind leads to an increase of the energy loss rate per unit distance, while the tail wind reduces the amount of energy lost per unit distance. During the execution of the descent, the guidance system will add more energy (directly related to fuel consumption and noise) or will release it by deploying speed brakes to safely and efficiently guide the aircraft along the planned trajectory. According to Fig. IV-11(c), for the tactical and NMPC guidance strategies, the energy deviation at the metering fix is always close to zero, independently of the wind speed error. However, in strategic and hybrid guidance more dispersion is found for large wind speed errors. The same conclusion applies for the
IV.4 Discussion 73 time deviation (see Fig. IV-11(b)). In this case, however, the tactical, hybrid and NMPC guidance strategies show excellent figures of CTA adherence, while the time deviation for the strategic guidance shows a linear correlation with the longitudinal wind speed error. Figure IV-11(f) shows that the number of re-plans for the strategic and hybrid guidance strategies strongly depends on the absolute value of the longitudinal wind error. However, the number of re-plans of these guidance strategies is still much lower than that of the NMPC strategy, which generates a new descent plan starting at the current aircraft state every 30 seconds. IV.4 Discussion This chapter compared the performance of four guidance strategies based on time and energy management concepts: tactical, strategic, hybrid and NMPC. Results from a batch study demonstrated that the four guidance variants can safely accomplish a CTA during a CDO with negligible time and energy deviations, even in presence of significant weather and aircraft performance modelling errors. The tactical guidance strategy requires continuous actions on the throttle and speed brakes to remain on the planned energy level and arrive on time, leading to more fuel consumption and noise nuisance if compared with other strategies. The hybrid strategy shows the best results in terms of environmental impact mitigation, and also requires less pilot actions to deploy speed brakes. However, due to its strategic behaviour with respect to specific energy deviations, the metering fix is reached with larger specific energy deviations than for other strategies. In this context, the NMPC strategy demonstrated to be the most robust in terms of specific energy and time deviations at the metering fix, while simultaneously showing environmental impact figures comparable to those of the hybrid and strategic. The main conclusion of this study is that these four guidance strategies provide very similar performance. However, even these small differences could be significant at aggregated level. Since all the guidance strategies showed similar results in terms of fuel consumption and time and energy error at the metering fix, the additional robustness of the NMPC guidance strategy promotes its use for this application. Compared to the other guidance strategies investigated herein, NMPC provides more situation awareness for the pilot, since deviations from the (dynamic) trajectory plan are lower than for strategic and hybrid. Also for this reason, the NMPC reduces the chances of obtaining infeasible solution as a result of a trajectory update. Note that when the strategic and hybrid strategies update the trajectory plan due to an excessive time and/or energy error reaching the threshold, it might be too late to react and satisfy all the constraints. In such cases, any of the strategies discussed in Section II.2.4.1 could be used to recover feasibility as long as the operational limits of the aircraft are not overpassed, yet degrading the performance of the operation. Last but not least, another advantage of the NMPC is that, if at certain moment during execution of the descent up-to-date and accurate wind data is obtained from, e.g., the leading aircraft in the procedure, re-planning the trajectory plan may be beneficial in both environmental impact and constraints compliance terms. The combination of NMPC and wind networking concept will be further investigated in Chapter VI. A limitation of the experiment presented herein is that the remaining descent after the metering fix was not taken into account. If the energy of the aircraft at this fix were too low, additional thrust would be needed after overflying it. Similarly, if the energy were too high, it would be required to use speed brakes and/or to deploy high-lifting devices or the landing gear earlier. After the metering fix, however, the approach procedure is heavy constrained, since the aircraft has to intercept the ILS glide slope in few nautical miles and thereafter follow a constant ground flight path angle (not allowing for energy management). Last but not least, in this chapter the guidance strategies were not optimised for noise explicitly, even if noise nuisance was used as metric for comparison. In future research, descents will be optimised for noise (and emissions) explicitly, and even better results may be obtained in terms of these metrics.
A complex system that works is invariably found to have evolved from a simple system that works. —John Gaule Any sufficiently advanced technology is indistinguishable from magic. —Arthur C. Clarke V Model predictive control to meet time constraints in optimal descents In Chapter IV, the performance of various guidance strategies build on the time and energy managed operations (TEMO) concept (De Jong et al.,2015) were compared using a high-fidelity flight simulator located at Netherlands Aerospace Centre (NLR), focusing in the environmental impact mitigation of these guidance strategies as well as their ability to comply with operational constraints (of particular interest, time constraints). Hundreds of optimal descents to Barcelona-El Prat airport subject to controlled time of arrival (CTA) at a metering fix were simulated by introducing errors in the weather and aircraft performance models implemented in a research flight management system (FMS). Results showed that nonlinear model predictive control (NMPC) (Diehl et al.,2008), a guidance strategy based on a frequent update of the optimal trajectory plan during the execution of the descent (e.g., every 30 seconds), could be very robust in terms of correcting energy (speed and altitude) and time deviations at the metering fix, providing at the same time acceptable fuel consumption and noise nuisance figures, if compared to existing TEMO strategies. Another advantage of NMPC is that deviations from the (continuously adjusted) trajectory plan are always small, thus providing more situational awareness for the pilot and being less prone to infeasible solutions as a result of a trajectory plan update. Based on these conclusions, NMPC was identified as a very attractive guidance strategy to enable continuous descent operations (CDOs) in high traffic demand scenarios. A practical limitation of the basic NMPC for real-life applications, however, was that each trajectory update could take up to 10 seconds. The main reason of this excessive computation time was that the trajectory planner formulated the optimal control problem in the remaining time horizon from scratch, without taking advantage of the nonlinear programming (NLP) structure, parametric sensitivity or primaldual solution of a similar optimisation problem computed at a previous time sample (see Section II.2.4.3). In this chapter, two additional variants of the basic NMPC will be investigated, which take advantage of parametric sensitivities to update the optimal trajectory plan quasi-instantaneously. The hypothesis is that sensitivity-based strategies will provide considerable reduction in execution time (thus allowing to updated the trajectory plan more often), simplification of the optimisation algorithm (thus improving convergence and robustness), at the expense of a small reduction on the accuracy of the solution. 75
76 Chapter V - Model predictive control to meet time constraints in optimal descents In order to accomplish that, more than 4,000 descents per guidance strategy were simulated at the National Aeronautics and Space Administration (NASA), supported by their expertise in trajectory operations and, in particular, CDOs and interval management. Errors in the parameters used by the FMS to describe the wind profile were intentionally introduced in the simulations. The open-loop fuel consumption as well as time and energy errors at the metering fix, obtained by applying the initial optimal control plan (computed at the TOD) in open-loop subject to these modelling errors, were compared with those of the pure sensitivity-based NMPC (SbNMPC), which uses only parametric sensitivities to update the optimal trajectory plan at each time sample; the ideal NMPC (INMPC) (Suwartadi et al.,2017), which ideally updates the optimal trajectory plan at each time sample by solving a rigorous NLP optimisation problem without delay (i.e., zero execution time); and the advanced-step NMPC (AsNMPC), which combines an early rigorous re-calculation of the trajectory plan based on the predicted state at the next time sample, with a sensitivitybased update to correct perturbations after measuring the actual state of the aircraft and parameters. Last but not least, the sensitivity of the performance metrics to the mean longitudinal wind error experienced during the descent was also investigated for the INMPC, if compared to the open-loop execution. It should be noted that the flight simulator used in this experiment was simpler than that used in Chapter IV. In particular, the simulator consisted on a set of accurate ordinary differential equation (ODE) integrators that, given the controls of the aircraft and its current state, calculated the future position at the next time sample using CVODES, a solver for stiff and nonstiff ODE systems bundled into the SUite of Nonlinear and DIfferential/ALgebraic Equation Solvers (SUNDIALS) (Hindmarsh et al.,2005). V.1 State of the art Generic NMPC implementations, such as that used in Chapter IV, frequently update the optimal trajectory plan by solving a rigorous NLP optimisation problem. Ideally, the trajectory should be updated without delay at each time sample, right after measuring the current state of the aircraft ( ~ X) and estimating the parameters of the model (~ d). In practical real-time applications, however, solving the rigorous NLP optimisation problem may take a significant amount of time, leading to potential stability issues and/or degrading the overall performance of the descent operation (Suwartadi et al.,2017). In order to reduce the execution time, educated simplifications in the models used by the FMS have been used, for instance by Kim & Hull (1995), at the expense of compromising the accuracy of the solution. Other NMPC implementations compensate for computational delay by starting the optimisation process in advance, setting the initial conditions of the new trajectory plan to the predicted state of the aircraft at the next time sample (Zavala & Biegler,2009), or at a look-ahead time equal to the estimated execution time. In the later case, the unpredictability of the execution time remains a critical issue (Wetzel,1996). An alternative method widely used in process industries consists of computing fast trajectory updates in the neighbourhood of the last optimal solution using the theory of neighbouring extremals (Pesch,1979). Using this approach, parametric sensitivities are obtained by linearisation of the necessary conditions of optimality, also known as Karush-Kuhn-Tucker conditions (KKT) conditions, along the (optimal) active trajectory plan to rapidly update it for small perturbations in the parameters of the model. This strategy, known as SbNMPC, has shown important benefits in terms of execution time (W¨ urth et al.,2008b,2009;Wolf & Marquardt,2016). For more technical details about this strategy the reader is referred to Section II.3.4.2. V.2 Setup of the experiment In the experiment performed in this chapter, 4,143 descents per guidance strategy were simulated to asses the performance of the INMPC, SbNMPC and AsNMPC guidance strategies in the presence of errors in the wind forecast model used by the trajectory planner. All the descents were subject to CTA at a metering fix, and to several altitude and speed constraints of a realistic flight procedure. First, Section V.2.1 presents the scenario. Then, the generic models introduced in Section II.1 are particularised in Section V.2.2. Analogously, the generic optimal control problem stated in Section II.2.1 is particularised in Section V.2.3 with specific cost function and constraints. The selected case studies are listed in Section V.2.4, and the workflow of a single simulation is summarised in Section V.2.5.
V.2 Setup of the experiment 77 V.2.1 Scenario The Denver international airport (DEN) was selected as the scenario for the experiment performed in this chapter. Fig. V-1 shows the BOSSS TWO standard arrival procedure at DEN, which was the starting point to define the route and the vertical profile used in this study. The metering fix where CTA were assigned by an hypothetical air traffic control (ATC) during the simulations was the final approach point (DYMON in Fig. V-1). The aircraft type used for the assessment was an Airbus A320-214. Figure V-1: BOSSS TWO standard arrival procedure (source: US AIP) V.2.2 Models The models required to formulate the optimisation of the aircraft trajectory, which were presented in Section II.1, are particularised below for the experiment conducted in this chapter. V.2.2.1 Aircraft dynamics model A variant of the γ-command model described by Eq. (III.2) was implemented in the trajectory planner. Different from typical approaches, the independent variable of this variant is the distance to go (s) and not the time. The selection of sas the independent variable was driven by the fact that during an ideal CDO, with no intervention from the ATC except for the negotiation of the CTA, the aircraft will follow a closedloop route and the remaining distance to go will be known. In addition, this formulation replicates how constraints are defined in the current operational environment, thereby enabling more precise modelling. In the distance-based γ-command model, the state vector is ~x = [t, v, h]T, and the control vector is ~u = [γ, T, β]T. The dynamics of ~x are expressed by the following system of ODEs: dt ds=t0=1 vcos γ+ws(h) dv ds=v0=T−D(v, h, γ) m−gsin γ1 vcos γ+ws(h) dh ds=h0=vsin γ vcos γ+ws(h). (V.1) Note that wxwas neglected for the sake of simplicity.
78 Chapter V - Model predictive control to meet time constraints in optimal descents V.2.2.2 Aircraft performance model Different alternatives can be used to model Tidle,Tmax,CD,q,qidle, as mentioned in Section II.1.2. In this experiment, the base of aircraft data (BADA) v4 model was adopted (Poles et al.,2010). Unfortunately, BADA v4 still does not include a model for the effects of the speed brakes on the drag coefficient. As a workaround, the contribution of the speed brakes was modelled as an extra linear term in the generic BADA v4 drag coefficient model (CD,BADA). A similar approach was adopted in the study performed in Chapter IV: CD(CL, M, β) = CD,BADA(CL, M) + CDββ. (V.2) Figure V-2 show the propulsive functions of the BADA v4 performance model for the A320-214 used in this study. Analogously, Fig. V-3 shows the drag coefficient. h[FL] 0 50 100 150 200 250 300 350 Mach [-] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Tmax/δ[kN] 75 100 125 150 175 200 225 250 275 (a) Maximum thrust h[FL] 0 50 100 150 200 250 300 350 Mach [-] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Tidle/δ[kN] −40 −30 −20 −10 0 10 (b) Idle thrust h[FL] 0 50 100 150 200 250 300 350 Mach [-] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 qidle/δ√θ[kg h−1] −1000 −500 0 500 1000 1500 2000 2500 3000 (c) Idle fuel flow T/δ[kN] −50 0 50 100 150 200 250 Mach [-] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 q/δ√θ[kg h−1] 0 5000 10000 15000 20000 25000 (d) Nominal fuel flow Figure V-2: BADA 4 propulsive model for the A320-214 V.2.2.3 Atmosphere model In this experiment, the international standard atmosphere (ISA) model, established to provide a common reference for temperature and pressure, was considered to represent the normalised temperature θand pressure δas a function of hthrough Eqs. (II.4) and (II.5), respectively. Then, the normalised density was computed assuming the perfect gas law relationship with Eq. (II.6). The longitudinal component of the wind (ws) was modelled by a smoothing spline of the form:
V.2 Setup of the experiment 79 CL[-] 0.0 0.2 0.4 0.6 0.8 1.0 Mach [-] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 CD[-] 0.02 0.04 0.06 0.08 0.10 0.12 0.14 Figure V-3: BADA 4 aerodynamic model for the A320-214 ws(h) = nc X i=1 ciBi(h),(V.3) where Bi,i= 1, . . . , nc, are B-spline basis functions and ~c = [c1, . . . , cnc]are control points of the spline (de Boor,1972). It should be noted that the longitudinal wind was modelled as a function of the altitude only, as done in similar works (De Jong et al.,2015). Yet, it would be straightforward to consider dependence also on the distance to go, which is the independent variable of the model. V.2.3 Optimal control problem formulation The trajectory planner used in this experiment generates the optimal plan trajectory subject to a CTA by solving a multi-phase, constrained and nonlinear optimal control problem, as described in Section II.2.1. This section particularises the cost function and the dynamic, path, algebraic and/or interior-point constraints for each one of the phases composing the flight profile of the procedure illustrated in Fig. V-1. The differential constraints that enforce the dynamics of the aircraft (~ f) were particularised by the following distance-based γ-command model, independently of the phase: ~ fj= 1 T−D(v,h,γ) m−gsin γ vsin γ 1 vcos γ+ws(h);j= 0, . . . , P −1.(V.4) From all the waypoints of the route, only the altitude and speed constraints at QUAIL, BOSSS, CHAPP and DYMON were modelled. In order to accomplish that, the descent was divided into P= 4 different phases, with associated phase-dependent path, algebraic and/or interior-point constraints. It should be noted, however, that the fact of modelling few constraints of the real procedure is not a shortcoming nor a limitation of the model. The model proposed in this study could handle an unlimited number of phases and associated constraints, yet few constraints were selected aiming to represent a futuristic and less restricted procedure facilitating CDOs, as well as to ease the interpretation of the results. Table V-1 wraps up the different phases and their associated constraints. This table is also illustrated in Fig. V-4. Note that the last row of Table V-1 includes phase-independent path constraints (i.e., applying all along the descent), which ensure that the aircraft airspeed remains within operational limits, and that the maximum and minimum descent gradients, thrust and speed brakes are not exceeded. In Table V-1,β= 0 and β= 1 indicate that speed brakes are retracted and fully extended, respectively; the minimum descent gradient was set to γmin =−7◦; and the values for maximum operative CAS (VMO) and maximum operative Mach (MMO) were obtained from the BADA v4 performance model. In addition, the following terminal constraints were set at the metering fix DYMON:
86 Chapter V - Model predictive control to meet time constraints in optimal descents (a) Open-loop (b) INMPC Figure V-9: Planned and executed trajectories by guidance strategy (2018-04-21 00:00 with a look-ahead time of +1 hour) modulation process around FL260 to avoid the altitude interval with extreme and unexpected tail wind, given the best available wind profile estimation at that time sample (see Fig. V-7). However, according to the executed trajectory, a later trajectory update (not shown in Fig. V-9(b)) which had more accurate information of the actual wind profile requested to start exchanging potential for kinetic energy around 1,000 ft earlier. After QUAIL, the updated trajectory demanded a sightly lower speed than that of the initial plan, in order to compensate the unexpected head wind in the altitude interval FL170–7,000 ft. It should be noted that, using NMPC, the optimal control applied at each time sample was based on the best estimation of the wind profile at the current altitude. The actual wind conditions that the aircraft would encounter at lower altitudes, however, were still unknown. Even if the control applied at each time sample was optimal for the estimated wind, it was sub-optimal for the actual (unknown) wind conditions (provided, of course, that the estimation of the wind downstream was not perfect). This can be observed by comparing Figs. V-9(b) and V-8(b): the executed trajectory (resulting from applying the optimal control for the best estimated wind profile at each time sample), and the trajectory update computed at the TOD with accurate information of the actual wind for the entire descent, are similar yet not identical. For the wind forecast error of this illustrative example, the open-loop trajectory arrived at DYMON 27 seconds earlier than the CTA and 600 ft above the enforced specific energy level. For any of the NMPC strategies investigated in this chapter, the aircraft arrived at DYMON with negligible time and specific energy error, while simultaneously using the minimum required fuel consumption and speed brakes. V.3.3 Aggregated results This section presents the aggregated performance metrics for the 4,143 case studies per guidance strategy. Figure V-10 shows the time error at the metering fix (DYMON) with respect to the enforced CTA for the INMPC, AsNMPC and SbNMPC guidance strategies. The time error that would be achieved by applying the optimal control resulting from the initial trajectory plan in open-loop is also shown in Fig. V-10(d). Analogously, Fig. V-11 shows the specific energy (Es)4error at DYMON. This kind of plots visualise the distribution of the data. The red horizontal line represents the average value, the central black horizontal line represents the median value, and the whiskers extend out to the minimum and maximum values. According to Fig. V-10, the more the look-ahead time of the wind forecast, the more time error (in absolute value) is observed, independently of the guidance strategy. The maximum time error at the metering fix (115 s) was realised when executing a trajectory plan computed with a wind forecast of +6 h in open-loop (note that the scale of Fig. V-10(d) is different from the other subfigures). The median time error at the metering fix when considering a wind forecast of +6 h is around 15 s for the open-loop execution. For 4The specific energy is defined as the total energy of the aircraft divided by the aircraft weight. By definition, the units of the specific energy are ft.
V.3 Results 87 +1 +3 +6 Wind forecast look-ahead time [h] 0 10 20 30 Time deviation [s] (a) NMPC +1 +3 +6 Wind forecast look-ahead time [h] 0 10 20 30 Time deviation [s] (b) SbNMPC +1 +3 +6 Wind forecast look-ahead time [h] 0 10 20 30 Time deviation [s] (c) AsNMPC +1 +3 +6 Wind forecast look-ahead time [h] 0 20 40 60 80 100 120 Time deviation [s] (d) Open-loop Figure V-10: Time deviation at the metering fix for different NMPC variants a wind forecast of +1 h, this value is reduced to 8 s. For the set of descents simulated with a wind forecast of +6 h, the distribution of time errors is concentrated in the range 0-60 s. When using strategies build on NMPC concepts, the time error is drastically reduced if compared to the open-loop execution of the initial trajectory plan. In particular, for all case studies analysed herein (one year of descent operations at DEN using the BOSSS TWO procedure) the time error is always lower than 27 s. The median time error for any of the NMPC guidance strategies is lower than 2 s for the set of descents simulated with a wind forecast of +6 h. In addition, the distribution of time errors is always concentrated below 10 s. Similar results can be observed for the specific energy error. According to Fig. V-11(d), the metering fix could be achieved with specific energy errors up to 1,800 ft by implementing the control of the initial trajectory plan computed at the TOD in open-loop. Conversely, guidance strategies build on NMPC concepts achieve the metering fix with negligible specific energy error, always lower than 150 ft and with a median value around 10 ft for a wind forecast of +6 h. In addition, the distribution of specific energy errors for the open-loop execution extends up to 1,000 ft, while for any of the NMPC strategies it extends up to 100 ft. Figures V-10 and V-11 show that the performance of the SbNMPC and AsNMPC, in terms of energy and time errors, is equivalent to that of INMPC. Results agree with those shown in Section V.3.1, where the trajectory updates based on parametric sensitivities were similar to those obtained by solving the rigorous NLP optimisation problem, considering relatively small perturbations in the vector of NLP parameters. In the model proposed in this chapter, the optimal control computed at each time sample could demand to modulate energy with the elevator or to add/remove energy to/from the system by means of additional thrust or speed brakes. Ideally, the NMPC optimiser will attempt to obtain an energy-neutral trajectory. In certain conditions with significant errors in the model parameters, however, energy modu-
88 Chapter V - Model predictive control to meet time constraints in optimal descents +1 +3 +6 Wind forecast look-ahead time [h] 0 50 100 150 Specific energy deviation [ft] (a) NMPC +1 +3 +6 Wind forecast look-ahead time [h] 0 50 100 150 Specific energy deviation [ft] (b) SbNMPC +1 +3 +6 Wind forecast look-ahead time [h] 0 50 100 150 Specific energy deviation [ft] (c) AsNMPC +1 +3 +6 Wind forecast look-ahead time [h] 0 500 1000 1500 2000 Specific energy deviation [ft] (d) Open-loop Figure V-11: Specific energy deviation at the metering fix for different NMPC variants lation by means of elevator control may not be sufficient to obtain a solution satisfying all the constraints. In this case, the NMPC trajectory optimiser would calculate the optimal amount of energy to be added or removed in terms of fuel consumption and use of speed brakes such that all constraints are satisfied. Table V-2 shows the mean, maximum and standard deviation of the specific energy removed by using speed brakes. For each case study, this metric was computed as follows: ∆Esβ=1 2 S mg ZtF tI v3(t)σ(t)ρSSLCDβ β(t)dt(V.9) For each guidance strategy and forecast lookahead times, the metrics shown in Table V-2 aggregate the results of the corresponding 1,381 case studies. Note that results for the open-loop execution are not shown in Table V-2 because the initial plan was always computed without speed brakes (i.e., β(t)=0). Table V-2: Metrics of energy removed by using speed brakes for the different NMPC variants and wind forecast lookahead times Guidance strategy Mean [ft] Standard deviation [ft] Maximum [ft] +1 h +3 h +6 h +1 h +3 h +6 h +1 h +3 h +6 h INMPC 2 7 13 12 34 45 187 780 641 SbNMPC 2 8 14 16 41 52 340 977 668 AsNMPC 2 8 14 11 40 55 113 914 904 According to Table V-2, and independently of the guidance strategy, the more the look-ahead time of
V.3 Results 89 the wind forecast, the more the mean and standard deviation of energy removed by using speed brakes. Interestingly, the maximum value of ∆Esβcorresponds to a case study with a +3 h wind forecast look-ahead time. The sensitivity of the energy removed by using speed brakes with respect to the mean longitudinal wind error encountered during the execution of the descent will be presented in Section V.3.4. Table V-2 also shows very similar results for the three strategies based on NMPC, meaning that updating the trajectory by solving a rigorous NLP or by using parametric sensitivities would result in almost identical number of speed brakes deployments during the descent, at aggregated level. Not surprisingly, the worst figures correspond to the SbNMPC, followed by the AsNMPC. These two NMPC strategies update the trajectory plan by using an approximated method, which in some cases could be sub-optimal. The INMPC show the best results, yet the differences with the other strategies are very small. It is also interesting to remark the significant drop in the three metrics when comparing the results for a wind forecast look-ahead time of +3 hours with those for a wind forecast look-ahead time of +1 hour. For example, for the AsNMPC the maximum specific energy removed by using speed brakes decreases from 914 ft to 113 ft, the mean value from 14 ft to 2 ft, and the standard deviation from 40 ft to 14 ft. Analogously, Fig. V-12 shows the difference between the executed trajectory and the open-loop execution of the initial plan, in terms of specific energy added by means of thrust above idle during the descent. For each case study, this metric was computed according to the following formula: ∆EsT=1 mg ZtF tI (T(t)−Tidle(v(t), h(t))) dtNMPC −ZtF tI (T(t)−Tidle(v(t), h(t))) dtOL.(V.10) Due to the particular altitude and speed constraints of the BOSSS TWO procedure (see Table V-1), the initial trajectory plan was not computed assuming a complete engine-idle descent, but included some segments with additional thrust above idle. For instance, additional thrust was typically required in the 2nd phase, in order to keep the speed within the allowed bounds while flying at constant altitude. Note that, since the wind encountered during the execution was different from that used to generate the trajectory plan, the aircraft may deviate from the altitude and speed plan. Accordingly, the idle thrust for the actual flight conditions (which depends on the altitude and speed) may be different from that of the initial plan at the same time sample. In other words, altitude and speed deviations may lead to an actual Tidle different from that assumed in the plan. It is also important to remark that, during the execution of the descent, the thrust of the aircraft was not allowed to be lower than Tidle for the actual flight conditions: T(t)actual = T(t)plan if T(t)plan ≥Tidle(v(t), h(t))actual Tidle(v(t), h(t))actual otherwhise.(V.11) Results shown in Fig. V-12(a) demonstrate that, in most of the case studies, the specific energy added to the aircraft by means of thrust when implementing any NMPC variant is identical to that of the open-loop execution (since the median and mean values are close to zero). In addition, even if the maximum value for each combination of guidance strategy and wind forecast look-ahead time is higher than the minim value (in absolute terms), the main distribution is almost symmetric, meaning that the NMPC guidance strategies have similar chances to demand additional specific energy than to remove it by reducing the thrust, if compared to the open-loop execution. These results indicate that the NMPC strategies could take advantage of the reiterated re-planning of the trajectory plan, based on the actual state of the aircraft and a continuously improving estimation of the parameters that describe the model, to execute the descent with the same thrust than initially planned. Taking into consideration the results shown in Figs. V-10 and V-11, it can be concluded the NMPC strategies are able to satisfy operational constraints and to accurately comply with the CTA without requiring extra energy, at aggregated level. In order to further analyse this interesting conclusion, Table V-3 shows, for each NMPC variant, the percentage of case studies not requiring speed-brakes or additional thrust during the entire descent. In addition, the percentage of case studies that performed an energy-neutral trajectory while complying with the CTA are also shown (i.e., uncertainties in the execution were cancelled by energy modulation). Not surprisingly, the higher the look-ahead time of the wind forecast, the lower the percentage of case studies not requiring use of speed brakes and/or additional thrust during the descent. An interesting result of Table V-3 is that the figures for the three NMPC variants are almost identical. Another striking result is
90 Chapter V - Model predictive control to meet time constraints in optimal descents +1 +3 +6 Wind forecast look-ahead time [h] −500 0 500 1000 Specific energy added with T[ft] (a) NMPC +1 +3 +6 Wind forecast look-ahead time [h] −500 0 500 1000 Specific energy added with T[ft] (b) SbNMPC +1 +3 +6 Wind forecast look-ahead time [h] −500 0 500 1000 Specific energy added with T[ft] (c) AsNMPC Figure V-12: Additional specific energy difference with respect to the open-loop execution for the different NMPC variants Table V-3: Percentage of energy neutral trajectories for the different NMPC variants and forecast lookahead times Guidance strategy Without speed brakes Without additional thrust Energy-neutral +1 h +3 h +6 h +1 h +3 h +6 h +1 h +3 h +6 h INMPC 98.3 93.7 90.0 82.1 75.6 75.5 80.6 70.3 67.0 SbNMPC 98.3 93.2 89.2 81.9 77.4 76.6 80.5 71.3 67.4 AsNMPC 98.6 93.2 89.9 82.0 76.3 76.7 80.9 70.5 68.4 that, even for look-ahead times of +6 h, around 70% of the descents could be performed without requiring neither additional thrust nor speed-brakes use during a whole year of operations at DEN airport following the BOSSS TWO procedure, provided that the sequencing and merging tasks were performed by means of CTAs and without any kind of tactical intervention (e.g., vectoring) from ATC after the TOD. Finally, Fig. V-13 shows the fuel consumption difference with respect to the open-loop execution for the different NMPC variants. Negative values indicate fuel savings, positive values extra fuel. Results shown in Fig. V-13 agree with those of Fig. V-12: the more the additional specific energy if compared to the open-loop execution, the more the extra fuel consumption. Interestingly, Fig. V-13 shows that using NMPC guidance strategies the fuel burned during the execution of the descent is typically the same than that burned by implementing the initial trajectory plan in open-loop (the median and mean values are close to zero). In addition, the distribution of fuel consumption difference is almost symmetrical, meaning that the NMPC strategies have similar chances to burn more fuel than to save it, if compared to the open-loop execution. The maximum extra fuel for each combination of guidance strategy and look-ahead
V.3 Results 91 +1 +3 +6 Wind forecast look-ahead time [h] −10 −5 0 5 10 15 Fuel consumption difference [%] (a) NMPC +1 +3 +6 Wind forecast look-ahead time [h] −10 −5 0 5 10 15 Fuel consumption difference [%] (b) SbNMPC +1 +3 +6 Wind forecast look-ahead time [h] −10 −5 0 5 10 15 Fuel consumption difference [%] (c) AsNMPC Figure V-13: Fuel consumption difference with respect to the open-loop execution for the different NMPC variants time of the wind forecast, however, is considerable larger than the corresponding maximum fuel saving. However, the maximum value was observed for a very particular scenario (out of 4,131 case studies) with severe differences between the initial wind forecast and the actual wind. In fact, the distribution of fuel consumption differences is concentrated in the range [−5%,5%]. Taking into account that these results include simulations for one year of wind forecasts at DEN, it can be concluded that the dispersion of fuel consumption differences is relatively small, independently of the NMPC variant. The open-loop execution implements the control plan that is optimal only for the initial wind forecast assumed at the TOD. In contrast, NMPC strategies could take advantage of the frequent update of the optimal trajectory plan based on a continuously improving wind forecast to save fuel. Again, taking into consideration the results shown in Figs. V-10 and V-11, it can be concluded the NMPC strategies assessed herein are able to satisfy altitude and speed constraints as well as to accurately comply with the enforced CTA without requiring additional fuel consumption, at aggregated level. Based on the results shown in this section, it can be also concluded that the three NMPC guidance strategies assessed in this PhD thesis show very similar performance at aggregated level, even if the method used to update the optimal trajectory plan is sightly different. Recall that the INMPC solves a rigorous NLP optimisation problem at each time sample; the SbNMPC solves a quadratic programming (QP) optimisation problem to update the optimal trajectory plan; and the AsNMPC updates the trajectory in advance by solving a rigorous NLPoptimisation problem (as the INMPC) and then corrects for potential prediction errors by solving a QP optimisation problem (as the SbNMPC). Figure V-14 shows the normalised distribution of execution times for the two optimisation methods, extracted from the whole set of simulations. According to Fig. V-14, the maximum execution time for the QP algorithm is one order of magnitude
92 Chapter V - Model predictive control to meet time constraints in optimal descents 0.9 1.0 0.0 0.1 PN1QNi Algorithm to update the trajectory plan Normalised execution time [-] Median Mean Figure V-14: Normalised execution time for the NLP and QP algorithms. 1 corresponds to the maximum execution time for all trajectory updates of all simulations, which was 10 seconds below that for the NLP algorithm. In addition, the mean and median execution times for the QP algorithm are approximately half of those for the NLP algorithm. The most important conclusion arising from Fig. V- 14, however, is that the dispersion of execution times for the QP algorithm is much lower than for the NLP algorithm, thus demonstrating the additional robustness of the QP algorithm. V.3.4 Sensitivity to the wind error Previous section assessed the performance of various NMPC guidance strategies in presence of errors in the wind profile forecast, and compared the results with those obtained by executing the initial trajectory plan in open-loop. The time and energy deviations induced by errors in the longitudinal component of the wind are significantly higher than those caused by errors in the cross and vertical components. For this reason, Figs. V-15 to V-17 only show the correlation between the different metrics of the performance assessment and the mean longitudinal wind error. The mean longitudinal wind error for each case study was computed by averaging the longitudinal wind error of all time samples. Positive values indicate more tail wind than initially expected, negative values manifest a stronger head wind if compared to the forecast. Each data point in these figures corresponds to a particular case study. In order to keep this figure simple, only the results for the open-loop execution and the INMPC are shown. Yet, it should be noted that results for the remaining strategies are very similar to those of the INMPC, as shown in previous sections. According to Fig. V-15(a), when executing the initial trajectory plan in open-loop, the time error at the metering fix shows a strong and negative correlation with the mean longitudinal wind error, i.e., the stronger the unexpected head wind, the later the aircraft will arrive with respect to the CTA. Analogously, the specific energy error at the metering fix (see Fig. V-15(b)) presents a strong, positive correlation with the mean longitudinal error, meaning that the stronger the unexpected head wind, the lower the energy of the aircraft will be when arriving at the metering fix. For any of the NMPC strategies, the time error appears to be weakly correlated with the mean longitudinal wind error; and the specific energy error is not correlated at all. The differences in the correlation patterns for the OL and NMPC guidance strategies are caused by the corrective actions performed by the later, which repeatedly generates the best optimal trajectory plan to satisfy the terminal constraints at the metering fix while simultaneously minimising the fuel consumption and speed brakes usage. Regarding the extra energy removed by using speed-brakes (see Fig. V-16(b)), it can be observed that for the open-loop execution speed brakes are never used (because the initial trajectory plan is always computed without extra drag). Generally speaking, for the NMPC guidance strategies the stronger the unexpected head wind, the more energy is removed. There exist some particular case studies in which energy is removed for negative mean longitudinal wind error. Even if the mean longitudinal wind error is negative for these case studies, in some regions of the descent more head wind than expected is encountered, which requires the use of speed brakes to satisfy all the constraints.
V.3 Results 93 (a) Time error (b) Specific energy error Figure V-15: Time and specific energy error at the metering fix as a function of the mean longitudinal wind error (a) Added by means of T(b) Removed by means of β Figure V-16: Specific energy added/removed with respect to the initial plan as a function of the mean longitudinal wind error for the INMPC and OL Figure V-17: Fuel consumption difference with respect to initial plan as a function of the mean longitudinal wind error
94 Chapter V - Model predictive control to meet time constraints in optimal descents Finally, according to Figs. V-16(a) and V-17 the specific energy difference due to thrust and the fuel consumption difference show very similar patterns. Both performance metrics present a negative correlation with the mean longitudinal wind error, meaning that the stronger the unexpected tail wind, the less thrust and fuel consumption will be required if compared to the initial trajectory plan computed at the TOD. Remember that for the open-loop execution, the differences in fuel consumption and energy added by means of thrust are only caused by altitude and speed deviations from the initial plan (changing therefore the value of the idle thrust). Conversely, for the NMPC guidance strategies, these differences are also caused by trajectory updates that might modify the states and controls plans. For this reason, results for the NMPC strategies show more dispersion than for the open-loop execution. V.4 Discussion This chapter described the application of two additional NMPC guidance strategy that use parametric sensitivities to update the optimal trajectory during the descent of an aircraft subject to time constraints. Results showed that trajectory updates resulting from solving a rigorous NLP optimisation problem are very similar to those obtained by using parametric sensitivities. Yet, the later method to update the optimal trajectory is simpler, faster and more robust in terms of dispersion of execution times and convergence. Accordingly, a guidance strategy that updates the optimal trajectory at each time sample by using parametric sensitivities, called SbNMPC, showed equivalent performance (in terms of fuel consumption, extra thrust, speed brakes use and extra fuel consumption) to the INMPC, which instantaneously updates the optimal trajectory at each time sample by solving a rigorous NLP optimisation problem. A third NMPC variant, which combines an early re-calculation of the optimal trajectory with a sensitivity-based update, also showed similar performance to the INMPC. These results encourage the use of the later strategy. If compared with the generic NMPC, which solves a rigorous NLP optimisation problem right after measuring the actual state and parameters of the model, the AsNMPC has the advantage of not being affected by computational delay. This is only true, however, if the optimal trajectory plan can be computed in background before reaching the next time sample. Recall that the overall performance, capability to meet constraints and stability of the generic NMPC deteriorate as the execution time increases. Moreover, if compared to the SbNMPC, the AsNMPC is less affected by large perturbations in the vector of parameters, since the sensitivity-based update is not performed over an approximated trajectory plan but over a very accurate trajectory plan generated recently by solving a rigorous NLP optimisation problem. Any trajectory obtained by using parametric sensitivities is just an approximation to the actual optimal trajectory, and therefore includes some approximation error, which obviously increases with the magnitude of the perturbation. This error would be propagated over time if only sensitivity-based updates were performed (as in the case of the SbNMPC). The rigorous NLP optimisation problem solved by the AsNMPC at each time sample, based on the predicted state and parameters of the model at the next one, acts as a mechanisms to ”reset” these approximation errors. Moreover, since the predictions errors in the state and parameters of the model are expected to be relatively small, provided that the time interval between consecutive time samples is sufficiently small, the subsequent sensitivity-based approximation when measuring the actual state and parameters of the model is expected to be very fast and accurate. In addition, results showed that using any of the NMPC guidance strategies, energy modulation through elevator control typically suffices to compensate time and energy deviations at the metering fix caused by errors in the wind forecast. In particular, around 70% of the descents at DEN following the BOSSS TWO procedure were able to perform an energy-neutral trajectory for a wind forecast look-ahead time of +6 hours. However, in some circumstances with significant errors in the wind forecast, these guidance strategies could demand to add or remove energy by means of thrust or speed brakes, respectively, to satisfy operational constraints. In this case, the trajectory planner would compute the optimal energy to be added or removed which, at each time sample, minimises the fuel consumption and speed brakes usage while simultaneously satisfying the constraints. In this context, promising results showed that NMPC guidance strategies can accurately comply with CTAs and other operational constraints in presence of significant wind prediction errors at no extra fuel consumption (considering the aggregated results of around 16,000 simulations). Most common languages used in aviation safety critical systems are C,C++ and Ada. In the experiment of this chapter, and taking into account that the developed system was just a prototype, the formulation of
V.4 Discussion 95 the optimal control problem (i.e., definition of constraints and cost function of the trajectory optimisation problem) was performed in Python, but the resulting NLP optimisation problem was solved by using offthe-shell NLP solvers. Many of these solvers are written in Cand could be technically integrated into the avionics systems. Regarding the hardware, modern FMSs provide significant improvements to already highly reliable computers, as well as increased processing power and memory. The computational demand of the proposed algorithm is not extremely high, yet the feasibility of running it using current on-board resources deserves a more comprehensive study. What drives value in avionics, however, is compliance to mandates, reliability and safety. On-board trajectory planning and guidance algorithms must be simple for the pilot to be easily understood, as well as robust enough to have hundreds of lives depend on them. Accordingly, new features being added into conventional FMSs typically require few extra computational resources, and there is little incentive for the the latest programming paradigms or innovatory technology. Because the NMPC guidance algorithm uses a NLP solver for the trajectory optimisation process, this trajectory management system can be characterised as an adaptive system with nondeterministic feature, which conflicts with current civil aviation certification processes. In particular, NLP solvers lack from determinism and cannot guarantee to always find a solution. In addition, the execution time has a lot of variance and in some unpredictable situations could be unacceptably high. Accordingly, in order to implement any NMPC variant in a real FMS, even if having enough hardware capabilities, operational and certification issues will need to be tackled.
102 Chapter VI - Wind-networking to improve time predictability and fuel efficiency of descents 60 40 20 0 Approximated longitudinal wind speed [kt] 0 100 200 400 Altitude [FL] ⌧0 ⌧10 ⌧40 h0 h10 h40 Forecast Actual (a) Only ownship observations (µ= 0.0) 60 40 20 0 Approximated longitudinal wind speed [kt] 0 100 200 400 Altitude [FL] ⌧0 ⌧10 ⌧40 h0 h10 h40 Forecast Actual (b) µ= 0.5 60 40 20 0 Approximated longitudinal wind speed [kt] 0 100 200 400 Altitude [FL] ⌧0 ⌧10 ⌧40 h0 h10 h40 Forecast Actual (c) µ= 1.0 Figure VI-3: RAP wind forecast and analysis (17-09-2019 at 18:00 for a look-ahead time of +6 hours) VI.3.2 Aggregated results This section presents the performance metrics for all case studies of the experiment. Figure VI-4 shows the time deviation at the metering fix (DYMON) with respect to the enforced CTA for three different values of µ and for the case in which the wind forecast is not updated, i.e., the initial wind forecast obtained from RAP keeps static throughout the descent. The time deviation that would result by applying the optimal control of the initial plan in open-loop was already shown in Fig. V-10(d), which could be as high as 110 seconds, with a mean and median value around 20 and 15 seconds for the worst forecast look-ahead time (+6 hours), respectively. When using NMPC guidance, both maximum, mean and median time error drastically reduce, independently of the mechanism adopted to update the wind profile (if any). For instance, when using the static wind profile approach, the time error is lower than 30 seconds for all case studies, and the distribution shows that most of the time deviations are lower than 20 seconds, with a median value of 2.5 seconds and a mean value around 5 seconds for a forecast look-ahead time of +6 hours. Even if at first sight 30 seconds could seem a very large time error, it should be noted that the simulations performed in this chapter considered all the wind forecasts for +1, +3 and +6 h look-ahead time during one year of operations at DEN airport. For the simulations in which the wind forecast is updated with only ownship observations (i.e.,µ= 0.0), the maximum time deviation further reduces to 25 seconds, while most of the case studies show a time deviation lower than 10 seconds. Using only ownship measurements, the time error is approximately halved with respect to that using a static wind forecast. Results show that if additional wind data emitted by aircraft in the neighbourhood were used to update the wind profile, the time error would be negligible for most case studies. Figure VI-5 show the analogous results for the specific energy deviation at the metering fix. The specific energy deviation at the metering fix that would result by applying the optimal control of the initial plan in open-loop was already shown in Fig. V-11(d), with maximum specific energy errors up to 1,800
VI.3 Results 103 +1 +3 +6 Wind forecast look-ahead time [h] 0 10 20 30 Time deviation [s] (a) Static +1 +3 +6 Wind forecast look-ahead time [h] 0 5 10 15 20 25 30 Time deviation [s] (b) µ= 0.0 +1 +3 +6 Wind forecast look-ahead time [h] 0.0 2.5 5.0 7.5 10.0 12.5 15.0 Time deviation [s] (c) µ= 0.5 +1 +3 +6 Wind forecast look-ahead time [h] 0.0 2.5 5.0 7.5 10.0 12.5 15.0 Time deviation [s] (d) µ= 1.0 Figure VI-4: Time deviation at the metering fix for different values fo µ ft as well as a mean and median values around 600 ft and 500 ft, respectively, for a wind forecast lookahead time of +6 hours. Conversely, when executing the descent with NMPC, the metering fix is achieved with much smaller specific energy deviation. In particular, for the static wind profile approach, the specific energy error is lower than 300 ft for all case studies. When the wind profile is updated at each time sample yet using only ownship wind observations, the maximum specific energy error reduces to, approximately, 140 ft. Finally, if additional wind observations from nearby aircraft were used to update the wind profile at each time sample, the specific energy error would be lower than 50 ft for all case studies. As discussed in Chapter V, the optimal control computed at each time sample could claim to modulate energy with the elevator or to add/remove energy to/from the system by means of additional thrust or speed brakes. Ideally, the NMPC trajectory planner will attempt to obtain an energy-neutral trajectory. In certain conditions, however, energy modulation may not be sufficient to obtain a solution satisfying all the constraints. If this were the case, the trajectory planner would calculate the optimal amount of energy to be added or removed, in terms of fuel consumption and speed brakes use, while satisfying all the constraints. Table VI-1 shows the mean, maximum and standard deviation of the specific energy removed by using speed brakes. Similarly to Chapter V, for each case study this metric was computed according to Eq. V.9. It should be noted that, for each value of µand look-ahead time of the wind forecast, results shown in Table V-2 aggregate the results of the corresponding 1,381 case studies. According to Table V-2, and independently of µ, the more the look-ahead time of the wind forecast, the more the mean energy removed by using speed brakes. For the standard deviation, however, when using wind observations emitted by nearby aircraft to update the wind profile (i.e., µ > 0), the larger values are found for a +3 h wind forecast look-ahead time. As observed in Chapter V, the maximum value of ∆Esβ corresponds to a case study with a +3 h wind forecast look-ahead time.
104 Chapter VI - Wind-networking to improve time predictability and fuel efficiency of descents +1 +3 +6 Wind forecast look-ahead time [h] 0 100 200 300 Specific energy deviation [ft] (a) Static +1 +3 +6 Wind forecast look-ahead time [h] 0 25 50 75 100 125 150 Specific energy deviation [ft] (b) µ= 0.0 +1 +3 +6 Wind forecast look-ahead time [h] 0 10 20 30 40 50 60 Specific energy deviation [ft] (c) µ= 0.5 +1 +3 +6 Wind forecast look-ahead time [h] 0 10 20 30 40 50 60 Specific energy deviation [ft] (d) µ= 1.0 Figure VI-5: Specific energy deviation at the metering fix for different values of µ Table VI-1: Metrics of energy removed by using speed brakes for different values of µand wind forecast lookahead times Guidance strategy Mean [ft] Standard deviation [ft] Maximum [ft] +1 h +3 h +6 h +1 h +3 h +6 h +1 h +3 h +6 h Static 2 13 22 14 37 49 161 481 453 µ= 0.02 7 13 12 34 45 187 780 641 µ= 0.51 2 3 24 31 25 586 831 439 µ= 1.01 2 3 23 33 29 518 982 421 The largest mean and standard deviation, independently of the wind forecast look-ahead time, correspond to the static wind forecast. These metrics decrease as the value of µincreases. In particular, for µ= 1 the mean energy removed by using speed brakes for a +3 h wind forecast look-ahead time is only 3 ft, with a standard deviation of 29 ft. For the maximum value, however, it is more difficult to determine which value of µgives the best results. It should be noted that adding wind observations not only reduces the mean and standard deviation of the energy removed by using speed brakes, but simultaneously reduces the time and energy error at the metering fix, if compared to the static wind profile and µ= 0. Analogously, Fig. VI-6 shows the difference between the executed trajectory and the open-loop execution of the initial plan, in terms of specific energy added by means of thrust above idle during the descent. For each case study, this metric was computed according to Eq. (V.10). By comparing the distributions of additional thrust for the static wind forecast with those corresponding to µ= 0, it can be concluded that updating the wind profile using only ownship observation does not
VI.3 Results 105 have a great impact on the additional thrust used during the descent. Therefore, the time and energy error at the metering fix can be significantly reduced without requiring additional thrust by simply updating the wind profile with ownship observations. When including additional wind observations emitted by nearby aircraft, however, even if the median and mean additional thrust are still close to zero (and even decrease sightly for all wind forecast look-ahead times), the dispersion is larger, as are the maximum and minimum values. Therefore, the huge reduction of time and specific energy deviations observed in Figs. VI-4 and VI- 5for µ > 0comes at a cost: there is a trade-off between reducing the dispersion of additional thrust and reducing the time and specific energy deviations at the metering fix. +1 +3 +6 Wind forecast look-ahead time [h] −500 0 500 1000 Specific energy added with T[ft] (a) Static +1 +3 +6 Wind forecast look-ahead time [h] −500 0 500 1000 Specific energy added with T[ft] (b) µ= 0.0 +1 +3 +6 Wind forecast look-ahead time [h] −500 0 500 1000 Specific energy added with T[ft] (c) µ= 0.5 +1 +3 +6 Wind forecast look-ahead time [h] −500 0 500 1000 Specific energy added with T[ft] (d) µ= 1.0 Figure VI-6: Specific energy added by means of Twith respect to the initial plan for different values of µ Table VI-2 shows, for each strategy to update the wind profile forecast during the descent, the percentage of case studies not requiring speed-brakes or additional thrust during the entire descent. In addition, the percentage of case studies that performed at energy-neutral trajectory while complying with the CTA are also shown in this table. Table VI-2: Percentage of energy neutral trajectories for different values of µand forecast look-ahead times Wind forecast update mode Without speed brakes Without additional thrust Energy-neutral +1 h +3 h +6 h +1 h +3 h +6 h +1 h +3 h +6 h Static 97.4 87.8 80.5 76.3 71.9 73.8 74.4 61.5 55.9 µ= 0.098.3 93.7 90.0 82.1 75.6 75.5 80.6 70.3 67.0 µ= 0.599.2 98.9 97.8 79.9 75.1 74.0 79.6 74.9 73.0 µ= 1.099.5 99.0 98.6 80.8 76.2 75.2 80.5 75.9 74.5 According to Table VI-2, and as expected, the more the forecast look-ahead time, the less the per-
106 Chapter VI - Wind-networking to improve time predictability and fuel efficiency of descents centage of case studies not requiring use of speed brakes and/or additional thrust during the descent. For a given forecast look-ahead time, increasing the number of wind observations used to update the wind profile typically reduces the chances of requiring speed brakes. Regarding the use of thrust, it is difficult to determine the probability of executing a descent without additional thrust depending on the number of weather observations used to update the forecast wind profile. Finally, the probability of executing an energy-neutral descent, in which neither additional thrust nor speed brakes use are required to satisfy constraints, typically increases with the amount of wind observations used to update the wind profile forecast. This enhancement becomes even more significant for large forecast look-ahead times. In particular, when using a static wind profile forecast with a look-ahead time of +6 h, only 55.9% of the descents are able to perform an energy-neutral trajectory. Conversely, when updating the wind profile forecast with a cadence of observations corresponding to µ= 1.0, the percentage of energy-neutral trajectories increases to 74.5%. Finally, Fig. VI-7 shows the difference in fuel consumption between the executed trajectory and the open-loop execution of the initial plan, with negative values indicating fuel savings and positive values indicating that additional fuel was required. Results shown in Fig. VI-7 agree with those of Fig. VI-6: the more specific energy added by means of thrust, the more fuel consumption. As already discussed in Chapter V, when implementing NMPC strategies the fuel burned during the execution of the descent is typically the same than that achieved by implementing the initial trajectory plan in open-loop (the median and mean values are close to zero). In addition, the distribution of fuel consumption difference is almost symmetrical, meaning that the NMPC strategies have the same chances to burn more fuel than to save it, if compared to the open-loop execution. The same conclusions than for the additional specific energy added by means of thrust apply here: there is a trade-off between reducing the dispersion of additional fuel and reducing the time and specific energy deviations at the metering fix. +1 +3 +6 Wind forecast look-ahead time [h] −10 −5 0 5 10 15 Fuel consumption difference [%] (a) Static +1 +3 +6 Wind forecast look-ahead time [h] −10 −5 0 5 10 15 Fuel consumption difference [%] (b) µ= 0.0 +1 +3 +6 Wind forecast look-ahead time [h] −10 −5 0 5 10 15 Fuel consumption difference [%] (c) µ= 0.5 +1 +3 +6 Wind forecast look-ahead time [h] −10 −5 0 5 10 15 Fuel consumption difference [%] (d) µ= 1.0 Figure VI-7: Fuel consumption difference with respect to the initial plan for different values of µ
VI.4 Discussion 107 VI.4 Discussion In this chapter the NMPC guidance strategy, which repeatedly updates the optimal trajectory of the aircraft, was combined with wind networking concepts, in which aircraft collaborate by sharing their wind observations to improve the wind profile forecast on-board and in real-time. Results show that the performance of the conventional NMPC strategies assessed in Chapter V, which (at most) used ownship wind observations to predict the wind profile just before updating the optimal trajectory plan, could be significantly improved when taking advantage of up-to-date wind observations emitted by nearby aircraft, notably in terms of time and specific energy deviations at the metering fix. Another interesting conclusion is that there exists a trade-off between reducing the dispersion of additional fuel and thrust and reducing the time and specific energy deviations at the metering fix. Generally speaking, however, increasing the value of µreduces the mean and median additional fuel and thrust. It was also observed that by including wind observations emitted by nearby aircraft the number of energy-neutral descents satisfying a CTA at a metering fix could be as high as 75%, considering a wind forecast with a look-ahead time of +6 hours. Without updating the wind profile during the descent, the number of energy-neutral descents satisfying a CTA at a metering fix would be around 55%. These are very promising results considering that the case studies selected for this experiment correspond to the RAP forecasts generated for +1, +3 and +6 h look-ahead times during one year over the region of Denver. It should be noted that for those case studies in which energy modulation was not sufficient to satisfy operational constraints, the amount of energy added/removed by the NMPC guidance system was optimal in terms of fuel consumption and use of speed brakes. The integration of CDOs subject to CTAs will require changes in both on-board trajectory management and ground automation systems. In particular, the interaction between aircraft and arrival manager (AMAN) systems on-ground is envisaged to follow the next sequence: 1. the ground system sends via datalink the area navigation (RNAV) route and runway in use to the aircraft, which subsequently computes the most optimal trajectory plan (not subject to time constraints) from its current position to the runway threshold. 2. the ground system requests (also via datalink) the estimated time of arrival for a certain fix along that route. 3. the aircraft sends the estimated time on arrival and possibly the time window (energy-neutral and/or absolute) at the requested fix. 4. the ground system sequences aircraft by assigning CTAs within the time windows transmitted by the aircraft (preferably using their estimated time of arrival, which is supposed to be the most optimal from the airline point of view). A crucial aspect for the acceptance of the AMAN system in future operations is the accuracy of its trajectory predictions. One of the main factors that directly affect the accuracy of any aircraft trajectory prediction is the wind data that is used in the model. The better the wind data, the closer is to the real wind that the aircraft will encounter, and the more accurate the resulting trajectory prediction (in particular, altitude, speed and time at the different fixes). Ground systems could also take advantage of wind observations emitted by airborne aircraft to improve their trajectory predictions. Enhanced trajectory predictions would also improve the perceived accuracy of the AMAN system, which will be very important to the overall acceptance by the air traffic controllers.
We can only see a short distance ahead, but we can see plenty there that needs to be done. —Alan Turing It’s the job that’s never started as takes longest to finish. —J. R. R. Tolkien VII Concluding remarks The growth in traffic implies increased pressure on the environmental sustainability of air transport. The introduction of more fuel-efficient profiles, however, is achieved at the cost of a reduction on capacity due to the more diverse trajectories that the air traffic control (ATC) need to handle. Continuous descent operations (CDOs) with controlled time of arrivals (CTAs) were identified as a potential solution to maintain capacity at acceptable levels, or even increasing it, while simultaneously allowing fuel-efficient descent profiles. The introduction of this kind or flight operations, however, will require the modernisation of both ground and airborne systems. The main objective of this PhD thesis was the development of on-board trajectory management algorithms to comply with time constraints during the descent of an aircraft in an environmentally friendly way, even in presence of uncertainties. During this work some questions arose that were assessed and some of them are still open and could be topics of further research. A brief summary and conclusions of the achieved results, as well as hints on the possible directions for future work, are presented in what follows. VII.1 Summary of contributions The main contributions of this PhD thesis are summarised as follows: • A trajectory management (planning and guidance) framework for the robust optimisation and execution of aircraft descent trajectories was established in Chapter II. This framework includes the formulation of a generic yet very complete optimal control problem to generate optimal trajectory plans taking into account any aircraft dynamics, performance and weather model as well as any number of flight phases, constraints and cost function. Moreover, several recommendations to face the most common issues that one could encounter when implementing nonlinear optimisation algorithms were also provided in Chapter II. Finally, several guidance strategies to execute the resulting trajectory plan, which take advantage of time and energy management concepts, were are also pro- 109
110 Chapter VII - Concluding remarks posed. • The potential applicability of these time and energy management concepts was demonstrated in Chapter III by quantifying the feasible time windows at a metering fix for energy-neutral descents. Remember that during an energy-neutral descent the total energy is not increased by means of thrust nor removed by using active drag devices, if compared to the initial trajectory plan. The most important conclusion of this study was that CTAs can be updated even after overflying the top of descent (TOD), and that considerable time windows could be attainable at idle thrust and without using speed brakes depending on the flight conditions. These results should encourage the introduction of time and energy management concepts for CDOs subject to CTAs. • A comprehensive comparison of guidance strategies build on time and energy management concepts was performed in Chapter IV. The set of guidance strategies compared in this PhD thesis included those already existing in the time and energy managed operations (TEMO) concept (i.e., tactical, strategic and hybrid) and a new guidance strategy called nonlinear model predictive control (NMPC), which has been in use in the process industries in chemical plants and oil refineries since the 1980s. Interestingly, results showed that all these guidance strategies could comply with CTAs with an average time error lower than 2 seconds and very attractive fuel consumption and noise nuisance figures, even in the presence of significant aircraft performance modelling and weather forecast errors. The main conclusion of this study, however, was that the performance of the four guidance strategies assessed are very similar. In spite of that, the NMPC was selected for further research because it demonstrated to be the most robust in terms of time and energy deviations at the metering fix. • In Chapter Vsensitivity-based alternatives to the generic NMPC were proposed to overcome its main limitation: the computational delay. Sensitivity-based NMPC (SbNMPC) and advanced-step NMPC (AsNMPC), which take advantage of parametric sensitivities to update the optimal trajectory plan solving a quadratic programming (QP) optimisation problem at each time sample, were compared to the generic NMPC, which solves a rigorous nonlinear programming (NLP) optimisation problem at each time sample. Fuel consumption, energy added and removed by means of thrust and speed-brakes, respectively, as well as time and energy deviations at the metering fix were the metrics selected for the comparison. Results from 16,000 simulations showed that trajectory updates resulting from solving a rigorous NLP optimisation problem are very similar to those obtained by using parametric sensitivities. Yet, the later method to update the trajectory plan is simpler, faster and more robust. Based on these results, the AsNMPC was the most recommended. Firstly, it has the advantage of not being affected by computational delay, in contrast to the ideal NMPC (INMPC). Secondly, if compared to the SbNMPC, it is less affected by large perturbations in the vector of parameters, since the sensitivity-based update is not performed over an approximated trajectory plan but over a very accurate trajectory generated recently by solving a rigorous NLP optimisation problem. There results shall encourage the integration of AsNMPC into future flight management system (FMS) because it is more robust than INMPC and SbNMPC, thus paving the way to certification processes. • Finally, the benefits of combining NMPC with wind networking concepts (i.e., sharing wind observations between aircraft) were investigated in Chapter VI, in terms of environmental impact and ability to comply with operational constraints. In addition, the sensitivity of the performance metrics to the rate of wind observations received from nearby aircraft was also quantified. Results showed that the performance of the conventional NMPC could be significantly improved when taking advantage of up-to-date wind observations emitted by nearby aircraft, in terms of time and energy deviations at the metering fix and fuel consumption. In particular, around 75% of the simulated flights were able to execute an energy-neutral descent when taking advantage of wind networking concepts to continuously improve the wind forecast in real-time, in contrast to the 55% observed when using a static wind forecast. Unfortunately, the current and accurate automatic dependent surveillance-broadcast (ADS-B) and Mode-S data being transmitted is not currently being used to update the wind profile used by the FMS to optimise the trajectory plan. Results arising from this work should encourage the aviation community to take advantage of data provided in the wind networked concept to further optimise the trajectory plan and reduce the overall environmental impact of aviation. Ground systems could also take advantage of wind observations emitted by airborne aircraft to improve their trajectory predictions and, consequently, increase the capacity and the safety of the air traffic management (ATM) system.
VII.2 Future Research 111 VII.2 Future Research During this PhD thesis new questions and research lines arose. Taking advantage of the trajectory optimisation framework that has been developed and for the sake of completeness, the following work items are proposed for the future: • It is almost impossible to predict how difficult it is to solve a particular model with a particular algorithm, especially for NLP optimisation problems. The only reliable way to find which solver to use for a particular class of models is so far to experiment. A benchmarking of optimal control problem methods, optimisation methods, solvers, and parameters of these components for the particular problem of trajectory optimisation would be very valuable for the aviation community. • The absolute time windows shown in Chapter III could be extended by allowing the aircraft to apply thrust or speed brakes. It would be interesting to assess the sensitivity of the time window to the amount of energy added or removed by means of thrust or speed brakes, respectively, as well as to quantify the theoretical limits of this time window as a function of the state of the aircraft. • It would be also interesting to compare the energy-neutral and absolute time windows with a fixed route with that achievable by using path lengthening or stretching, even though pilots and ATC would probably prefer energy management to reduce workload and increase situation awareness. • The NMPC guidance strategy updates the optimal trajectory plan at regular time intervals; while the strategic and hybrid guidance strategies update the optimal trajectory plan only when the time and/or energy deviations exceed a pre-defined threshold. In the former case, a continuous replanning of the optimal trajectory plan could be contrapositive if the models used by the FMS are not accurate enough; in the later case, it might be too late to react if the thresholds were not selected properly. In future work, a new guidance strategy based on reinforcement learning could be used to develop a system that ”learns” the best time samples to update the optimal trajectory plan, based on the current state of the aircraft and the environment (for instance, procedure being flown and meteorological conditions). The agent (the guidance system) would ought to take actions (re-plan or not re-plan) in an environment so as to maximise some notion of cumulative reward (e.g., the total fuel consumption of the entire descent). A large number of simulations would be required to train this agent. • All the experiment performed in this PhD thesis were performed from the perspective of a single aircraft. In future work, the potential effectiveness of the CTA window could be tested by a traffic simulation, aiming to assess the effectiveness of the CTA window in the traffic stream management.
118 BIBLIOGRAPHY SENZIG, D. A. & FLEMING, G. G. . 2009. Fuel consumption modeling in support of ATM environmental decision-making. In: 8th USA/Europe Air Traffic Management Research and Development Seminar (ATM2009). Napa, CA: FAA/EUROCONTROL. 20 SENZIG, D. A. , FLEMING, G. G. & IOVINELLI, R. J. . 2009. Modeling of Terminal-Area Airplane Fuel Consumption. Journal of Aircraft, 46(4), 1089–1093.20,42 SERGEEVA, M. , DELAHAYE, D. , MANCEL, C. & VIDOSAVLJEVIC, A. . 2017. Dynamic airspace configuration by genetic algorithm. Journal of Traffic and Transportation Engineering,4(3), 300 – 314.1 SESAR. 2015a. European ATM Master Plan. Tech. rept. Single European Sky ATM Research (SESAR) Joing Undertaking, Brussels, Belgium. Edition 2015. 1,2 SESAR. 2015b. SESAR Concept of Operations Step 1 Final Edition. Tech. rept. Single European Sky ATM Research (SESAR) Joing Undertaking, Brussels, Belgium. 41 SHRESTA, S. , NESKOVIC, D. , & WILLIAMS, S. S. . 2009. Analysis of Continuous Descent Benefits and Impacts during Daytime Operations. In: 8th USA/Europe Air Traffic Management Research and Development Seminar (ATM2009). Berlin, Germany: FAA/EUROCONTROL. 7 SLATER, G. L. . 2009. Study on variations in vertical profile for CDA descents. In: 9th AIAA Aviation Technology, Integration, and Operations Conference (ATIO).Aviation Technology, Integration, and Operations (ATIO) Conferences. Hilton Head, SC: AIAA. 20 SOLER, M. , OLIVARES, A. , STAFFETTI, E. & ZAPATA, D. . 2012. Framework for Aircraft Trajectory Planning Toward an Efficient Air Traffic Management. Journal of Aircraft,49(1), 341–348.2 SOLER, M. , OLIVARES, A. & STAFFETTI, E. . 2015. Multiphase Optimal Control Framework for Commercial Aircraft Four-Dimensional Flight-Planning Problems. Journal of Aircraft,52(1), 247–286.22 SOPJES, R. , DEJONG, P. M. A. , BORST, C. , VAN PAASSEN, R. & MULDER, M. . 2011. Continuous Descent Approaches with Variable Flight-Path Angles under Time Constraints. In: AIAA Guidance, Navigation, and Control Conference.Guidance, Navigation, and Control and Co-located Conferences. Portland, OR: AIAA. 6 SUN, J. , VU, H. , ELLERBROEK, J. & HOEKSTRA, J. . 2017 (November). Ground-based Wind Field Construction from Mode-S and ADS-B Data with a Novel Gas Particle Model. In: 7th SESAR Innovation Days.98 SUWARTADI, E. , KUNGURTSEV, V. & J ¨ ASCHKE, J. . 2017. Sensitivity-Based Economic NMPC with a Path-Following Approach. Processes,5(1), 8.76 TAKEICHI, N. & INAMI, D. . 2010. Arrival-Time Controllability of Tailored Arrival Subjected to Flight-Path Constraints. Journal of Aircraft,47(6), 2021–2029.42 TAKEICHI, N. . 2017. Nominal flight time optimization for arrival time scheduling through estimation/resolution of delay accumulation. Transportation Research Part C: Emerging Technologies,77(April), 433–443.42 THOMPSON, T. , MILLER, B. , MURPHY, C. , AUGUSTINE, S. , WHITEE, T. & SOUIHI, S. . 2013. Environmental Impacts of Continuousdescent Operations in Paris and New York Regions. In: 10th USA/Europe Air Traffic Management Research and Development Seminar (ATM2013). Chicago, IL: FAA/EUROCONTROL. 3 UEBBING-RUMKE, M. & TEMME, M.-M. . 2011. Controller AIDS for integrating negotiated continuous descent approaches into conventional landing traffic. In: 9th USA/Europe Air Traffic Management Research and Development Seminar (ATM2011). Berlin, Germany: FAA/EUROCONTROL. 9 VADA, J. , SLUPPHAUG, O. , JOHANSEN, T. A. & FOSS, B. A. . 2001. Linear MPC with optimal prioritized infeasibility handling: application, computational issues and stability. Automatica,37(11), 1835 – 1843.29 VAN DER WAL,H.,VOGEL,P.&WUBBEN, F. . 2001 (July). Voorschrift voor de berekening van de lden en lnight geluidbelasting in db(a) ten gevolge van vliegverkeer van en naar de luchtaven. Tech. rept. National Aerospace Laboratorium (NLR). 59 VERHOEVEN, R. , BUSSINK, F. J. L. , PRATS, X. & DALMAU, R. . 2016. Flight Testing Time and Energy Managed Operations. In: Society of Flight Test Engineers - European Chapter Symposium.11 VIRTANEN, K. , EHTAMO, H. , RAIVIO, T. & HAMALAINEN, R. P. . 1999. VIATO - Visual Interactive Aircraft Trajectory Optimization. IEEE Transactions of systems, man, and cybernetics - Part C: Applications and reviews,29(3), 409–421.25 WACHTER, A. & BIEGLER, L. . 2006. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming,106(1), 25–57.29 WARREN, A. & TONG, K. . 2002. Development of continuous descent approach concepts for noise abatement. In: 21st IEEE/AIAA Digital Avionics Systems Conference (DASC). Irvine, CA: IEEE. 3 WAT, J. , FOLLET, J. , MEAD, R. , BROWN, J. , KOK, R. , DIJKSTRA, F. & VERMEIJ, J. . 2006. In Service Demonstration of Advanced Arrival Techniques at Schiphol Airport. In: 6th AIAA Aviation Technology, Integration and Operations Conference (ATIO).Aviation Technology, Integration, and Operations (ATIO) Conferences. Wichita, Kansas: AIAA. 3 WEITZ, L. A. & BAI, X. . 2018. Using Model Predictive Control for Trajectory Optimization and to Meet Spacing Objectives. In: 2018 AIAA Guidance, Navigation, and Control Conference.AIAA SciTech Forum. Kissimmee, FL: AIAA. 58
BIBLIOGRAPHY 119 WETZEL, T. A. . 1996. Development of a finite difference neighboring optimal control law and application to the optimal landing of a reusable launch vehicle. PhD Thesis, Iowa State University. 76 WILSON, I. & HAFNER, F. . 2005. Benefit assessment of using continuous descent approaches at Atlanta. In: 24th IEEE/AIAA Digital Avionics Systems Conference (DASC). Washington, DC: IEEE. 6 WOLF, I. J. & MARQUARDT, W. . 2016. Fast NMPC schemes for regulatory and economic NMPC - A review. Journal of Process Control, 44(August), 162–183.37,38,76 W¨ URTH, L. , HANNEMANN, R. & MARQUARDT, W. . 2008a. An Efficient Strategy for Real-Time Dynamic Optimization based on Parametric Sensitivities. IFAC Proceedings Volumes,41(2), 1928–1933.38 W¨ URTH, L. , HANNEMANN, R. & MARQUARDT, W. . 2008b. An Efficient Strategy for Real-Time Dynamic Optimization based on Parametric Sensitivities. IFAC Proceedings Volumes,41(2), 1928–1933.76 W¨ URTH, L. , HANNEMANN, R. & MARQUARDT, W. . 2009. Neighboring-extremal updates for nonlinear model-predictive control and dynamic real-time optimization. Journal of process control,19(8), 1277–1288.28,37,76 ZAVALA, V. M. & BIEGLER, L. T. . 2009. The advanced-step NMPC controller : Optimality , stability and robustness. Automatica,45(1), 86–93.39,76