4-Dimensional Trajectory Optimisation Algorithm for Air Traffic Management Systems

This is the author pre-publication version. This paper does not include the changes arising from the revision, formatting and publishing processes. The final version that should be used for referencing is: A. Gardi, M. Marino, S. Ramasamy, T. Kistan, R. Sabatini, “4-Dimensional Trajectory Optimisation Algorithm for Air Traffic Management Systems”, IEEE/AIAA 35th Digital Avionics Systems Conference (DASC), Sacramento, CA, USA, 2016. 4-Dimensional Trajectory Optimisation Algorithm for Air Traffic Management Systems Alessandro Gardi, Matthew Marino, Subramaniam Ramasamy and Roberto Sabatini RMIT University – School of Engineering Melbourne, VIC 3000, Australia Abstract—This paper presents Multi Objective Trajectory Optimization (MOTO) algorithms that were developed for integration in state-of-the-art Air Traffic Management (ATM) and Air Traffic Flow Management (ATFM) systems. The MOTO algorithms are conceived for the automation-assisted replanning of 4-Dimensional Trajectories (4DT) when unforeseen perturbations arise at strategic and tactical online operational timeframes. The MOTO algorithms take into account updated weather and neighbouring traffic data, as well as the related forecasts from selected sources. Multiple user-defined operational, economic and environmental objectives can be integrated as necessary. Two different MOTO algorithms are developed for future implementation in ATM systems: an enroute variant and a Terminal Manoeuvring Area (TMA) variant. In particular, the automated optimal 4DT replanning algorithm for en-route airspace operations is restricted to constant flight level to avoid violating the current vertical airspace structure. As such, the complexity of the generated trajectories reduces to 2 dimensions plus time (2D+T), which are optimally represented in the present 2D ATM display formats. Departing traffic operations will also significantly benefit from MOTO-4D by enabling steep/continuous climb operations with optimal throttle, reducing perceived noise and gaseous emissions. Keywords—4-Dimensional Trajectory; Air Traffic Management; Decision Support System; Trajectory Optimization I. INTRODUCTION More effective and comprehensive implementations of flight trajectory optimisation techniques are being considered as a very promising pathway to enhance the efficiency and flexibility of online air traffic operations both for short and long-haul flights. The fundamental proposition is not new, but state-of-the-art flight planning methods are still based on the optimised vertical planning techniques initially developed in the 1970s and on lateral path planning based on optimal wind routing [1-6]. The known limitations are associated to the small set of optimality criteria (currently only fuel- and time-costs) and to the fact that the initial flight plan is the static entity assumed as a reference for every subsequent amendment. As a result, the limited initial optimality is progressively compromised when strategic or tactical Air Traffic Management (ATM) and Air Traffic Flow Management (ATFM) amendments are introduced. Suitably defined models can replicate the various operational and environmental aspects The work presented in this article was supported by THALES Australia under RMIT University Contract ID 0200312837 Trevor Kistan THALES Australia – Air Traffic Management Melbourne, VIC 3000, Australia that depend on the flown aircraft trajectory, allowing for more comprehensive and real-time Multi-Objective Trajectory Optimisation (MOTO) techniques to be integrated in airborne avionics and ground-based Communication, Navigation, Surveillance and ATM (CNS/ATM) Decision Support Systems (DSS). These novel DSS integrating real-time MOTO capabilities have the potential to support a further enhanced exploitation of airspace and airport capacities, which is emerging thanks to the advances in CNS technologies. This paper describes two dedicated MOTO functionalities for integration in state-of-the-art ATM and ATFM DSS. The MOTO algorithms are conceived for the automation-assisted replanning of 4-Dimensional Trajectories (4DT) when unforeseen perturbations arise at strategic and tactical online operational timeframes. The MOTO algorithms take into account updated weather and neighbouring traffic data, as well as the related forecasts from selected sources. Multiple operational, economic and environmental objectives can be user-defined as necessary. The project specifically addresses current and short-term future ATM operational paradigms and regulations and thus two different MOTO algorithms are developed: an en-route variant and a Terminal Maneuvering Area (TMA) variant. In particular, the automated optimal 4DT replanning algorithm for en-route airspace operations is restricted to constant flight level to avoid violating the current vertical airspace structure. As such, the complexity of the generated trajectories reduces to 2 dimensions plus time (2D+T), which are optimally represented in the present 2D ATM display formats. The Air Traffic Controller (ATCo) may amend the flight level of an optimized 2D+T trajectory in the traditional manner if necessary. The constant flight level limitation will theoretically produce a sub-optimal flight trajectory, however, the computed trajectory will be more efficient than vectors as the atmospheric wind field can be exploited to maximize aircraft ground speed while reducing fuel burn and emissions. The operational, economic and environmental benefits from MOTO are maximized in the TMA variant as the full 4DT MOTO (MOTO-4D) routines are exploited to increase the operational efficiency and reduce fuel burn, emissions and noise. In the current paradigm aircraft are handled on a "first come first serve" basis where multiple objectives can be user defined and applied to the optimization problem. However, an evolution to "just in time" 4DT based operations in the TMA is supported by the MOTO functionalities. As a consequence, TMA operations are allowed to evolve and become more efficient by minimising scheduling delays, providing automated path-stretching, reducing pollutants/noise emissions and supporting continuous descent approaches. As a consequence, ATCo workload will be relieved thanks to automation-assisted deconfliction and path planning, and this is expected to increase operational safety. Departing traffic operations will also significantly benefit from MOTO-4D by enabling steep/continuous climb operations with optimal throttle, reducing perceived noise and gaseous emissions. II. OPTIMAL CONTROL FORMULATION Trajectory optimisation studies methods to determine the best possible trajectory of a dynamical system in a finitedimensional manifold, in terms of specific objectives and adhering to given constraints and boundary conditions [7]. This definition corresponds to the definition of Optimal Control Problems (OCP), and consequently the most traditional and theoretically rigorous way to pose a Trajectory Optimisation Problem (TOP) is based on the optimal control theory. Most OCP solution methods are conventionally categorised as either direct methods if based on the transcription to a finite NonLinear Programming (NLP) problem, or indirect methods if theoretical derivations based on the calculus of variation are implemented to formulate a Boundary-Value-Problem (BVP) [8-11]. An additional class of OCP solution strategies is represented by heuristic methods. The optimal control formulation of TOP is based on a scalar time and on vectors of time-dependent state variables , timedependent control variables and system parameters . Based on these, the following components are defined: dynamic constraints, path constraints, boundary conditions and cost functions. These components characterise a well-posed OCP and guide the selection of an appropriate numerical solution method and multi-objective decision logic. Dynamic constraints describe the feasible motion of the system (i.e. the aircraft, in our case) within the TOP as in: ̇ (1) All non-differential constraints insisting on the system between the initial and final conditions are classified as path constraints, as they restrict the path of states and controls of the dynamical system. A generalised expression accounting for both equality and inequality constraints is [12]: (2) Boundary conditions define the values that state and control variables of the dynamical system shall have at the initial and final times. A generalised expression also accounting for relaxed boundary conditions is: ( ) ( ) (3) Performance indexes quantify the achievement of a particular objective by means of suitably defined cost functions. The generic formulation of a performance index that takes into account both integral and terminal costs was introduced by Bolza [13-15], and is expressed as: [ ( ) ( ) ] ∫ (4) The optimisation is classified as single-objective when an individual performance index J is introduced and multiobjective when two or multiple performance indexes Ji are defined. Different objectives can be conflicting, that is the attainment of a better would lead to a worse { [ ] } . Hence, the optimisation in terms of two or more objectives generates a number of possible compromise choices, for which a trade-off decision logic must be adopted to identify an individual solution. The first step of the theoretical derivations required by indirect methods to formulate a BVP problem consists in the Lagrangian relaxation as in: ( ( ) ) ( ) ( ) ∫ { ̇ Where are the Lagrangian multipliers and the Hamiltonian function, defined as: } (5) is (6) In addition to being mathematically complex, the application of all the analytical steps involved in indirect methods lead to very different BVP depending on the initial problem statement and this aspect limits the flexibility of indirect methods. Some notable examples are presented in [15]. Finally, since considerable nonlinearities are present in most aerospace TOP, the reduction of nontrivial cases into linear quadratic BVP is typically precluded, and therefore solutions have to be attempted by either iterative methods, which conventionally exhibit non-global convergence, or by heuristic solution techniques. Direct methods, on the other hand, prescribe the immediate parameterisation of states and controls in the case of direct collocation methods, or of controls only in the case of direct shooting. This involves adopting a basis of known linearly independent functions with unknown coefficients in the general form: ∑ (7) In global collocation (pseudospectral) methods, the evaluation of state and control vectors is performed at discrete collocation points across the problem domain using suitably defined orthogonal (spectral) interpolating functions [11, 16]. Similarly to the nonlinear BVP arising in indirect solution approaches, the NLP problems encountered in the direct solution of TOP are typically solved by iterative algorithms, or through some kind of heuristics. The initial steps involve the adoption of an n-dimensional Taylor series expansion of F(x) to the third term as in: ( ) ( ( ) ( where: ) (8) ) force T [N]; load factor N [ ]; bank angle [rad]. Other variables and parameters include: aircraft weight and aerodynamic drag D [N]; wind velocity vw in its three scalar components [m s-1]; gravitational acceleration g [m s-2]; Earth radius RE [m]; fuel flow FF [kg s-1] and thrust angle of attack [rad]. The aerodynamic drag is modelled as: (9) [ ( ) ] (10) (11) [ (14) where is the local air density retrieved from weather input data grid or a weather model, S is the reference wing surface, and are the parabolic drag coefficients typically available from aircraft performance databases such as Eurocontrol’s Base of Aircraft Data (BADA). The lift coefficient can be calculated from: (15) ] An iterative NLP solution can thus be formulated so that the search direction at step k based on the n-dimensional Newton method is written: The thrust force control variable is most frequently expressed as the product of the throttle coefficient (defined as dimensionless and ranging between 0 and 1), and the maximum thrust as in: (12) Various factors have to be considered when developing computationally efficient NLP solution strategies and some further detail is given in [17, 18]. A. Optimality criteria and dynamic constraints Optimality criteria and constraints are introduced in the optimisation problem by means of suitable models. Three degrees of freedom (3-DOF) point-mass flight dynamics models are currently preferred for TOP involving mediumlarge transport aircraft and for possible avionics and ATM system implementations, as the overall size of the resulting NLP is considerably lower than in the case of six degrees of freedom (6-DOF). A fairly comprehensive 3-DOF formulation used in our MOTO algorithms assumes variable aircraft mass, constant vertical gravity and the effects of winds and is therefore written as: ̇ ̇ ̇ ̇ ̇ [( ) ] (13) ̇ { ̇ where the state vector consists of: longitudinal velocity v [m s-1], flight path angle γ [rad]; track angle χ [rad]; geographic latitude [rad]; geographic longitude [rad]; altitude z [m]; aircraft mass m [kg]; and the control vector includes: thrust (16) This allows a natural nondimensionalisation of the control variable. For turbofan aircraft, the following empirical expressions were adopted in the development of BADA, to determine the climb thrust and the fuel flow , which operationally equates to the maximum thrust in all flight phases excluding take-off [19]: ( ) (17) (18) [ ] where is the throttle control, is the geopotential pressure altitude in feet, is the deviation from the standard atmosphere temperature in kelvin, is the true airspeed. are the empirical thrust and fuel flow coefficients, which are also supplied as part of BADA for a considerable number of currently operating aircraft [19]. In order to calculate pollutant emissions as a function of the fuel flow, the emission index is introduced as per the following definition: ∫ (19) where the generic Gaseous Pollutant (GP) should be replaced by the specific one being investigated. While carbon dioxide (CO2) emissions are characterised by an approximately , an constant emission index of empirical model for carbon monoxide (CO) and unburned hydrocarbons (HC) emission indexes ( ) in [g/Kg] at mean sea level based on nonlinear fit of experimental data from the ICAO emissions databank is: (20) where the fitting parameters accounting for the emissions of 165 currently operated civil turbofan engines are { } for CO and { } for HC [20]. The nitrogen oxides (NOX) emission index [g/Kg] based on the curve fitting of 177 currently operated civil aircraft engines is [20]: (21) III. MULTI-OBJECTIVE OPTIMALITY In the aviation domain, single and bi-objective optimisation techniques have been exploited for decades but they accounted only for flight time-related costs and fuel-related costs. These techniques have also been implemented in a number of current generation FMS in terms of the Cost Index (CI), that allowed an optimal selection of Calibrated Air Speed (CAS) / Mach number based on time and fuel costs only. Conflicting objectives arise when introducing multiple environmental, economic and operational criteria in our MOTO algorithms. Furthermore, the implementation of constraints that are either unfeasible or contrasting the attainment of better optimality also has to be addressed by adopting suitable multiobjective optimality techniques. These can either be expressed a priori (i.e. beforehand) or a posteriori (i.e. afterwards). A priori methods analyse ways to articulate the preferences to identify a combined objective which is then supplied to a single-objective TOP solution algorithm, and include weighted global criterion (including the simple weighted sum), weighted min-max, weighted product, exponential weighted criterion, lexicographic and physical programming methods [7, 21]. A posteriori methods allow the user (or suitably defined decision logics) to select an individual optimal solution from a Pareto front of trade-off choices, and include normal boundary intersection, normal constraint and physical programming methods [7, 21]. In general, a specific performance objective can be defined for each route segment. This performance objective is a multiobjective generalisation of the CI. In general, the weightings can be varied dynamically among the different flight phases of the flight. Since computational times are a crucial aspect in online 4DT planning applications, an a priori articulation of preference involving the weighted sum of the various performance indexes Ji is employed to combine the multiple conflicting operational, economic and environmental objectives. For a more detailed discussion the reader is referred to [7, 9, 11]. IV. MOTO FOR TERMINAL MANOEUVRING AREA The MOTO for Terminal Sequencing and Spacing (TSS) is based on the optimal control formulation presented in section 2. The optimal control formulation shows its best potential in this context as this is a truly 4D application where multiple variables and constraints are at play at the same time. The adoption of optimal control-based techniques in the TSS context allows accurate continuous descent profiles to be natively generated, CO/HC emissions to be minimised, as well as optimised path stretching to be introduced as necessary to achieve the set constraints. The challenges associated with optimal control-based MOTO implementations lie in the characteristics of the generated trajectory. More specifically, the mathematically optimal 4DT generated as output by these algorithms is a discretised Continuous/Piecewise Smooth (CPWS) curve, which in general may not be flyable by human pilots nor by conventional Automatic Flight Control Systems (AFCS), as it includes transition manoeuvres involving multiple simultaneous variations in the control inputs. Although continuous, the variations in control input profiles can also be particularly abrupt, exceeding what can be practically and safely feasible in the modelled aircraft. Moreover, the discretised CPWS consists of a very high number of 4D waypoints, which would have unacceptable impacts on the bandwidth of the data-link where they would be exchanged. Therefore, a post-processing stage is introduced to segment the discretised CPWS trajectory in feasible flight legs, including straight and level flight, straight climbs and descents, level turns, and climbing/descending turns. The final result is a concisely described 4DT consisting of feasible flight segments. The MOTO algorithm implemented for TMA operations accounting for the combined objective is represented in Fig. 1. MOTO-4D Initial Conditions {0, 0, z0, t0 , m0, v0, 0, 0, 0} GFS Weather Field {vw, p, RH, T} (, , z, t) Combined Objective 4DT Generation Global orthogonal collocation {f, f, zf, tf ,  f , f, f} Terminal Conditions Mathematically Optimal 4DT (CPWS) OPERATIONAL SMOOTHING Straight & Level Model Manoeuvre Identification Operational Transition Criteria Likelihood Level Turn Model Feasible 4DT Optimiser Straight Straight Turning Climb Descent Climb Model Model Model Turning Descent Model Final Solution: Optimal and Feasible 4DT Fig. 1. Block diagram of the MOTO-4D algorithm for TMA operations. The currently employed solution technique is based on direct solution methods of the global orthogonal collocation family (pseudospectral), which are arguably the most effective solution methods currently available [7]. Additionally, to further enhance the algorithm stability and convergence performances, path constraints and boundary conditions are automatically strengthened on all state and control variables to restrict the search domain as much as feasible.The TMA MOTO algorithm considers the controlled time of arrival target defined by Arrival Management (AMAN) functionalities. This is used as the final time constraint by the TMA MOTO algorithm. The Estimated Time of Arrival (ETA) may be computed at multiple fixes along the flight path. A. Operational Smoothing The 4DT intent data produced as output of the operational 4DT smoothing include 4D waypoints (latitude, longitude, altitude and time), fly-by/fly-over turn parameters, as well as performance criteria and restrictions. Based on the manoeuvre identification algorithm assessment on the CPWS trajectory, the lateral path is constructed in terms of segments (straight and turns), whose geometric characteristics depend on the required course change and the predicted Ground Speed (GS) of the aircraft during turns. The ATM DSS computes turn radius R and True Air Speed based on the selected altitude and takes into account the predicted wind speed at that altitude. The bank angle is selected based on the aircraft dynamics database, on the current altitude and on local/airline operational restrictions. The radius R of a turn can be conservatively calculated based on the maximum GS during the turn as: (22) where g is the gravity acceleration module. The turn arc length is simply given by: (23) where the track change is the difference in radians between the final and the initial ground tracks. The following prescriptions on , , lead distance (from the turn initiation to the 4DT waypoint) and abeam distance (between the 4DT waypoint and the point of the circular arc abeam it) are implemented to plan fixed-radius fly-by turns that comply with RTCA DO-229D, DO-236C and DO-283B [22, 23]: 1. The fly-by transition is defined by equation (22) combined with the following equations [22]: (24) ( ) (28) (29) 3. When transitioning from one airway to another, if both require a FRT at the common waypoint, the smaller of the two radii applicable shall be selected. In the case one of the two airways does not involve a FRT, the FRT of the other shall be implemented. The criteria for implementing the definitions above in different flight phases and at various altitudes are given in Table 2-7 of RTCA DO-229D [23]. V. MOTO FOR ENROUTE OPERATIONS Even neglecting current operational restrictions, the enroute context typically involves limited or no variations in cruise altitude. Additionally, the route distance and time involved are substantially greater than the transients associated with the aircraft flight dynamics model, which can therefore be neglected since the entire phase takes place in quasi-steadystate conditions. As a result, optimal control based techniques have limited or no advantages in this context. As a result, the aviation research community is investigating alternative formulations typically based on geometric steady-state trajectory approximations. The mainstream philosophies pursued involve either the generation of a grid of discrete 2D points across the cruise space and the implementation of suitable tree search algorithms to identify the optimal path, or the implementation of evolutionary algorithms or more general NLP solvers to optimise a discrete geometric trajectory. Predefined flight level (FL) conditions are typically considered to reduce the number of optimisation parameters and also to meet current operational paradigms that prescribe semi-circular cruise levels. The MOTO algorithm implemented for enroute operations adopts a polygonal trajectory model, simplified aircraft performance equations based on the steady-state approximation of eq. 13-21 and NLP solvers. The variables considered include the latitude and longitude of each and every trajectory waypoint, the airspeed flown (which dictates the throttle setting and therefore the pollutant emissions according to the steady-state approximation of equations 13-21. VI. TMA MOTO CASE STUDY |⃗ ⃗ | { (25) { 2. (26) { (27) The geometry of the Fixed Radius Transition (FRT) is defined by the track change and the radius . The Lead Distance and the Abeam Distance are defined based on the radius and the track change as per the following equations: The sequencing of dense arrival traffic towards a single final approach path was extensively assumed as a representative case study of online tactical TMA operations. The results of one exemplary simulation run are depicted in Fig. 2. The TMA MOTO considers the selected arrival sequence. Longitudinal separation is enforced at the mergepoint to ensure safe separation upon landing based on waketurbulence category and to prevent separation infringements in the approach phase itself. The assumed minimum longitudinal separation is 4 nautical miles on the approach path for medium category aircraft approaching at 140 knots, therefore the generated time slots are characterized by a 90~160 seconds separation depending on the wake-turbulence categories of two consecutive traffics. Fig. 3 depicts the computed 4DT in the AMAN schedule display format. Waypoints and lines depicted in magenta represent the flyable and concisely-described 4DT consisting of a limited number of fly-by and overfly 4D waypoints, obtained through the smoothing algorithm. Fig. 2. Results of the TMA MOTO exploiting the operational 4DT smoothing. requirements for tactical online data-link negotiation of the 4DT in the TMA context. VII. CONCLUSIONS This paper presented two Multi-Objective Trajectory Optimisation (MOTO) algorithms being developed for state-ofthe-art Air Traffic Management (ATM) Decision Support Systems (DSS). These MOTO algorithms allow the replanning of 4-Dimensional aircraft flight Trajectories (4DT) in strategic and tactical online operations, in the presence of dense air traffic, whenever airspace or air traffic reorganization is required due to tactically changing airspace conditions as well as promoting optimal rerouting around adverse weather. Simulation case studies consistently support the viability of the two proposed implementations: a fully-4D implementation for TMA based on optimal control formulation and a 2Dimensional plus Time (2D+T) steady-state polygonal path optimisation for the enroute context. In high air traffic density conditions, the complete process of 4DT intent generation, uplink and negotiation/validation can be consistently performed within the 180 seconds limit. These results meet the 3 minutes timeframe assumed for online tactical routing/rerouting tasks and make the approach feasible for quasi-real-time applications. The future algorithm evolutions will also incorporate various Communication, Navigation, Surveillance, ATM and Avionics (CNS+A) integrity monitoring and augmentation strategies currently being researched [24, 25]. ACKNOWLEDGMENT The authors gratefully acknowledge THALES Australia for supporting the research presented in this article under RMIT University Research Contract ID 0200312837. REFERENCES Fig. 3. Traffic sequencing represented in an AMAN scheduler plot. Monte Carlo simulations resulted in an average of 11 seconds for single newly generated 4DT intents, and consistently less than 30 seconds. The 4DT post-processing allowed to reduce discretised CPWS trajectories of up to 450 points into a number of fly-by and overfly 4D waypoints consistently below 20. These results meet the set design [1] Erzberger, H., McLean, J. D., and Barman, J. F., "Fixed-range optimum trajectories for short-haul aircraft", National Aeronautics and Space Administration (NASA) NASA TN D-8115, Ames research center, Moffett Field, CA, USA, 1975. [2] Barman, J. F. and Erzberger, H., "Fixed-range optimum trajectories for short-haul aircraft", Journal of Aircraft, vol. 13, pp. 748-754, 1976. [3] Erzberger, H. and Lee, H., "Characteristics of constrained optimum trajectories with specified range", National Aeronautics and Space Administration (NASA) NASA Technical Memorandum TM-78519, Ames research center, Moffett Field, CA, USA, 1978. [4] Sorensen, J. A., Morello, S. A., and Erzberger, H., "Application of trajectory optimization principles to minimize aircraft operating costs", 18th IEEE Conference on Decision and Control, vol. 1, pp. 415-421, 1979. [5] Lee, H. Q. and Erzberger, H., "Algorithm for fixed-range optimal trajectories", National Aeronautics and Space Administration (NASA) NASA Technical Paper TP-1565, Ames Research Center, Moffett Field, CA, USA, 1980. [6] Erzberger, H. and Lee, H., "Constrained optimum trajectories with specified range", Journal of Guidance, Control, and Dynamics, vol. 3, pp. 78-85, 1980. [7] Gardi, A., Sabatini, R., and Ramasamy, S., "Multi-objective optimisation of aircraft flight trajectories in the ATM and avionics context", Progress in Aerospace Sciences, vol. 83, 2016. [8] von Stryk, O. and Bulirsch, R., "Direct and indirect methods for trajectory optimization", Annals of Operations Research, vol. 37, pp. 357-373, 1992. [9] Visser, H. G., "A 4D Trajectory Optimization and Guidance Technique for Terminal Area Traffic Management", Delft University of Technology, 1994. [10] Betts, J. T., "Survey of numerical methods for trajectory optimization", Journal of Guidance, Control, and Dynamics, vol. 21, pp. 193-207, 1998. [11] Rao, A. V., "Survey of Numerical Methods for Optimal Control", Advances in the Astronautical Sciences, vol. 135, pp. 497-528, 2010. [12] Rao, A. V., "Trajectory Optimization", in Encyclopedia of Aerospace Engineering, John Wiley & Sons Ltd, 2010. [13] Bolza, O., Lectures on the calculus of variations (orig.μ Vorlesungen über Variationsrechnung), B. G. Teubner, Leipzig, Germany, 1909. [14] Berkovitz, L. D., "Variational methods in problems of control and programming", Journal of Mathematical Analysis and Applications, vol. 3, pp. 145-169, 1961. [15] Ben-Asher, J. Z., Optimal Control Theory with Aerospace Applications, Education Series, American Institute of Aeronautics and Astronautics (AIAA), Reston, VA, USA, 2010. [16] Rao, A. V., Benson, D. A., Darby, C., Patterson, M. A., Francolin, C., Sanders, I., et al., "Algorithm 902: GPOPS, A MATLAB software for solving multiple-phase optimal control problems using the gauss pseudospectral method", ACM Transactions on Mathematical Software, vol. 37, 2010. [17] Betts, J. T., Practical methods for optimal control and estimation using nonlinear programming, 2nd ed. vol. 19, SIAM, 2010. [18] Bertsekas, D. P., Nonlinear Programming, 2nd ed., Athena Scientific, Belmont, MA, USA, 1999. [19] "User Manual for the Base of Aircraft Data (BADA) Revision 3.11", Eurocontrol Experimental Centre (EEC) Technical/Scientific Report No. 13/04/16-01, Brétigny-sur-Orge, France, 2013. [20] Gardi, A., Sabatini, R., Kistan, T., Lim, Y., and Ramasamy, S., "4Dimensional Trajectory Functionalities for Air Traffic Management Systems", Integrated Communication, Navigation and Surveillance Conference (ICNS 2015), Herndon, VA, USA, 2015. [21] Marler, R. T. and Arora, J. S., "Survey of multi-objective optimization methods for engineering", Structural and Multidisciplinary Optimization, vol. 26, pp. 369-395, 2004. [22] RTCA, "Minimum Aviation System Performance Standards: Required Navigation Performance for Area Navigation", Radio Technical Committee for Aeronautics Inc. (RTCA) RTCA DO-236C, Washington DC, USA, 2013. [23] RTCA, "Minimum Operational Performance Standards for Global Positioning System/Satellite-Based Augmentation System Airborne Equipment", Radio Technical Commission for Aeronautics Inc. (RTCA) RTCA DO-229D, Washington DC, USA, 2013. [24] Sabatini, R., Moore, T., and Hill, C., "A new avionics-based GNSS integrity augmentation system: Part 1 - Fundamentals", Journal of Navigation, vol. 66, pp. 363-384, 2013. [25] Sabatini, R., Moore, T., and Hill, C., "A new avionics-based GNSS integrity augmentation system: Part 2 - Integrity flags", Journal of Navigation, vol. 66, pp. 501-522, 2013.