A New Computational Technique for the Generation of Optimised Aircraft Trajectories

1 Introduction

Trajectory optimisation has been a core research topic in the aerospace engineering domain for several decades. Most of the operational, economic and environmental performances of an aerospace vehicle’s mission are in fact directly correlated to the flown trajectory [1]. The maximisation of range or endurance performances and the minimisation of fuel-based or time-based costs are some of the objectives that have been frequently pursued in these studies. In the 1930’s, Ernst Zermelo studied the problem of wind-optimal aircraft routing [2,3]. This was one of the first and most challenging trajectory optimisation problems tackled in the aerospace domain, extensively investigated by De Jong and the Royal Netherlands Meteorological Institute in the 1970’s [4], and is still actively researched at present [5–7]. While the steady-state approximations of aircraft dynamics adopted in these and more generic aircraft performance studies have been very successful in providing the rough figures required to develop Multidisciplinary Design Optimisation (MDO) and basic flight planning algorithms, the limitations of this approach have been increasingly evident with the steady growth of air traffic, which has resulted in new constraints on route and airspeed and in an increasing congestion in low altitude terminal airspace regions. To resolve the emerging issues and simultaneously improve the predictability of traffic demand, the aviation research community has committed to better capture aircraft, airspace and airport dynamics in all trajectory planning, prediction and management algorithms that will be integrated in future avionics and Air Traffic Management (ATM) Decision Support Systems (DSS). This evolution will be enabled by shifting from legacy flight plans to 4-Dimensional Trajectory (4DT) intents which will be exchanged, negotiated and updated continuously by exploiting Air-to-ground Data-Links (ADL) communications [8–11].

When introducing aircraft dynamics in the problem statement, the formulation of the Trajectory Optimisation Problem (TOP) mathematically corresponds to the formulation of the Optimal Control Problem (OCP). Optimal control is exploited for a wide range of engineering problems, as it permits finding theoretically optimal control laws that minimise a performance criterion (cost functional) of a dynamic system through mathematical optimisation techniques. In the aerospace and avionics domain, the optimisation process involves defining aircraft dynamics, path costraints (imposed by adverse weather conditions or by ATM and Air Traffic Flow Management services) and a set of objectives (usually formulated to maximise fuel/time efficency and to reduce the environmental impacts) to compute trajectories that are optimal with respect to user-preferred performances. The aircraft dynamics model consists of a set of differential equations correlating the system states and its control inputs, which are usually accompanied by a set of constraints on the control inputs leading to a reduced set of admissible inputs or controls. The solution of the OCP consists in determining the admissible inputs that fulfill all constraints and, in doing so, minimize the cost functional [12]. The techniques available to solve aerospace OCP are reviewed in [13–16]. Traditionally, OCP were solved using indirect methods, applying the calculus of variations or Pontryagin’s maximum principle to satisfy first-order necessary conditions for optimality [13–16]. These methods are characterized by attempting an explicit solution of the optimality conditions stated in terms of the adjoint differential equations, the maximum principle, and associated boundary conditions [14]. This is practical for classical problems and some special weakly non-linear low dimensional systems. However, to obtain a solution of dynamic systems described by strongly non-linear differential equations, it is necessary to use numerical methods [13]. These methods suffer from the fact that the addition of new constraints requires deriving new necessary conditions. In many complex problems, determining the necessary conditions in a useful form can be a very difficult task [26]. As problems became more complex, indirect methods became increasingly harder to use, eventually being replaced by the more computationally intensive direct methods. Direct methods transcribe the continuous OCP into a parameter optimisation problem [8910111213141516]. Satisfaction of the system equations is accomplished by integrating them stepwise using either implicit or explicit rules; in either case, the effect is to generate non-linear constraint equations which must be satisfied by the parameters, which typically correspond to the discretised state and control histories [27]. The problem is thus converted from the original infinite dimensional OCP into a finite Non-Linear Programming (NLP) problem that can be solved using standard NLP solvers. A significant opportunity is recently identified when considering the application of very efficient direct numerical OCP solution methods for aircraft 4DT planning and management algorithms [1]. Visser, among others, analysed the application of OCP in 4D Trajectory Based Operations (TBO) [8]. Since the concept of 4D-TBO by itself required substantial research, the subsequent Programme for Harmonised ATM Research in EUROCONTROL (PHARE) research initiative focussed on the operational concept [17–19]. With the growing concerns for environmental sustainability of aviation in the last decade, the application of optimal control to compute environmentally optimal aircraft trajectories has been one of the key focus subjects of the Systems for Green Operations Integrated Technology Demonstrator (SGO-ITD) of the Clean Sky joint technology initiative for aeronautics and air transport [1,20–22]. As part of this and similar endeavours, a need for more versatile multi-objective formulations, encompassing operational, economic and environmental criteria at the same time was identified [23]. Although a large number of studies have tackled the application of optimal control techniques to aircraft flight trajectories, only a few have undertaken such investigation with the aim of implementing these aspects in avionics and ATM DSS [23]. A notable challenge consists in the transposition of the mathematically opimal trajectory computed by OCP solution methods into an operationally meaningful 4DT that can be flown by state-of-the-art Automated Flight Control Systems (AFCS) and that can be concisely described so to minimise the impacts on the data-link bandwidth and system processing [24]. Due to the aforementioned growing focus in multi-objective TOP formulations, computationally efficient bi- and multi-objective optimality techniques are an essential element of 4DT optimisation algorithms that will be integrated in future avionics and ATM DSS. The work presented in this paper investigates the adoption of custom developed techniques to solve aircraft TOP using direct methods to determine optimal trajectories. The computationally-efficient adaptive bisection ε-constraint scalarisation method is adapted to attain a computationally-efficient generation of Pareto frontiers for bi-objective optimal control problems (BOOCPs). This method is well suited for bi- and multi-objective trajectory optimisation algorithms that will be implemented in future avionics and ATM DSS. The underlying NLP solver which is used for the solution of the discretised BOOCP is the open-source large-scale non-linear solver IPOPT [28], which has been integrated with the MATLAB environment , where all the algorithms were numerically verified in representative simulation case studies. The basis of this work was initially presented in [30].

The paper is structured as follows. Section 2 formulates a generic single-objective continuous OCP that is discretised using pseudospectral techniques in Section 3. In Section 4, the discretised OCP is transcribed into a parameter optimisation problem that can be supplied to the NLP solver. In Section 5, the OCP is extended to encompass bi-objective problems, and the formulation of the adaptive bi-section ε-constrained scalarisation method is described in Section 6, supporting an efficient solution of bi-objective TOP with standard solvers. Section 7 discusses the NLP solver selection and its implementation in our algorithms. In Section 8, the techniques are validated against a classical problem of optimal control theory (i.e., the brachistochrone). Furthermore, a conventional vertical trajectory optimisation problem is evaluated. Section 9 presents the proposed implementation of these algorithms to solve realistic single- and multi-objective 4DT optimisation problems in avionics and ATM DSS for 4D-TBO. The paper terminates with some concluding remarks and directions for further research.

2 Optimal control problem formulation

In general, a Single-Objective Optimal Control Problem (SOOCP) is solved by finding the state trajectories x(t), the control trajectories u(t), and times t₀ and t_f in the interval t ∈ [t₀,t_f], that minimize a specific cost functional J. The problem is formulated as follows [31]:

where ϕ is the endpoint cost and L is the integrand cost, known as the Mayer and Lagrange cost respectively. The SOOCP is subjected to the following constraints that must be satisfied by the solution:

where ẋ represents the system dynamics in the form of differential constraints, and h and e are the path and event constraints respectively. The state, control and time variables are also bounded as follows:

3 Pseudospectral discretisation

Solving the OCP analytically is generally very difficult even for simple problems and becomes prohibitive for complex nonlinear problems. This difficulty can be eliminated adopting efficient numerical techniques on digital computers. By following this approach, the OCP is discretised such that state and control trajectories are represented by vectors of discrete values for each time interval. The solution of the resulting discretised OCP involves a number of steps, the most computationally intensive of which are those approximating the derivatives of the state trajectories at the discretisation nodes and integrating the cost functionals. Over the last few years, pseudospectral discretisation techniques have emerged as the most effective transcription methods for solving complex nonlinear OCP owing to their accuracy and speed, with an impressive convergence rate known as spectral accuracy [32]. In particular, for smooth problems, spectral accuracy implies an exponential convergence rate [33]. For these reasons, pseudospectral techniques are now very commonly adopted in aerospace TOP and have been selected for this work. The pseudospectral formulation starts with a polynomial expansion, which is well known in numerical analysis [34]. If τ₀τ₁,…,τ_N are N + 1 distinct values and f is a function whose values are given at those values, then a unique polynomial P(τ ) of degree at most N exists with:

Assume that {p_k}_{k = 0,1,…} is a system of algebraic polynomials (with degree of p_k = k) that are mutually orthogonal over the interval (–1, 1) with respect to a weight function w ≠ 0 [35]:

The classical Weierstrass theorem implies that such a system is complete in the space Lw2(−1,1). This is the space of functions v such that the following norm is finite:

The formal series of a function f∈Lw2(−1,1) in terms of the system {p_k} is:

where the coefficients of the expansion are given by:

Equation 13 represents the polynomial transform of f. For an integer N > 0, the truncated series of f of order N is the polynomial:

The polynomial P(τ)in equation 10 is known as a Lagrange polynomial and it is given by:

where

A particular class of orthogonal polynomials are the Legendre polynomials. The Legendre polynomials L_k (τ), k = 0, 1, …, N are the eigenfunctions of the singular Sturm-Liouville problem. If L_k (τ) is normalised (i.e., L_k (1) = 1) then for any k:

where [k/2] denotes the integral part of k/2. Consider the N-th degree Legendre orthogonal polynomial L_N(τ). The polynomial L_N has N – 1 extrema τ_k, i.e., LN˙(τk)=0, for k = 1,…, N – 1; that belong to the interval (–1,1). Adding the boundary points τ₀ = –1 and τ_N = 1, a total of N + 1 intervals are obtained, which are high-precision quadrature nodes termed the Legendre-Gauss-Lobatto nodes; there exists weights w_k, such that the quadrature formula

is exact for all polynomials p of degree ≤ 2N – 1. Legendre Gauss, Legendre-Gauss-Radau and Legendre-Gauss-Lobatto (LGL) series are the most common Legendre series, the choice of which determines the relative displacement of the quadrature nodes and their corresponding weights. Since equations for the quadrature nodes are not available in explicit form, the associated values have to be calculated using implicit numerical approximations. The quadrature weights can be expressed in closed form in terms of the nodes, as indicated in the following formulas for the LGL series, which is used extensively in the literature for trajectory optimisation:

If L(τ) is a general smooth function, then for a suitable N, its integral over τ ∊ [–1, 1] can be approximated as follows:

In a collocation method the fundamental representation of a smooth function f on τ ∊ [–1, 1] is in terms of its values at the discrete Gauss-type points. Derivatives of the function are approximated by analytic derivatives of the interpolating polynomial. The interpolating polynomial is denoted by I_N f and satisfies:

Since it is a polynomial of degree N, interpolated over the LGL nodes as a discrete expansion using Legendre polynomials, funtion f admits an expression of the form:

Therefore, from equation 14, the weights of the expansion are approximated as follows:

where

Since I_N f (τ) is an interpolant of f (τ) at the LGL nodes, and since the interpolating polynomial is unique, I_N f (τ) may be expressed as a Lagrange interpolating polynomial:

so that equations 24 and 27 are mathematically equivalent. It is possible to write the Lagrange basis (or characteristic) polynomials as follows:

The use of polynomial interpolation to approximate a function using the LGL points is known in the literature as the Legendre pseudospectral approximation method. Denote f ^N (τ)=I_N f (τ). Then,

It should be noted that L_k (τ_j) = 1 if k = j and L_k (τ_j) = 0, if k ≠ j, so that:

If the function f is known at the LGL quadrature points, one can compute an approximate derivative of f by differentiating the interpolant I_N f and evaluating it at the same nodes. The polynomial of degree N – 1

is called the Legendre interpolation derivative of f relative to the chosen set of quadrature nodes, since in general, it is different from the projection derivative (PNf)˙. The interpolation derivative at the nodes can be performed efficiently through matrix multiplication:

The entries (D_N) _kl can be computed by differentiating the characteristic Lagrange polynomials L_k of degree N, which are 1 at τ_l and 0 at all the other collocation points. For the commonly used Gauss-Lobatto points, the closed form for the first-derivative matrix is:

The discretisation process of the general OCP derived in Section 2 commences with the introduction of the following transformation [36]:

This results in the mapping:

The objective of the OCP solution process is now to find the state and control trajectories x(τ) and u(τ) respectively, in the interval τ ∈ [–1, 1], and times t₀ and t_f, which minimize the performance index:

subject to the following constraints and bounds:

In the Legendre pseudospectral approximation, the state and control trajectories x(τ) and u(τ) respectively, in the interval τ ∈ [–1, 1], are approximated by N^th order Lagrange polynomials x^N(τ) and u^N(τ) based on interpolation at the Legendre-Gauss-Lobatto nodes [33]:

where x^N(τ) and u^N(τ) are the Lagrange interpolating polynomials, and ϕ_k(τ) are known as Lagrange basis polynomials. It should be noted that ϕ_k (τ_j) = 1 if k = j and ϕ_k (τ_j) = 0 if k ≠ j such that:

The derivative of the state vector is approximated as follows:

where D is (N + 1) × (N + 1) differentiation matrix corresponding to the LGL nodes. By defining u^Nand x^N as the n_u × (N + 1) control and n_x × (N +1 ) state matrices storing the trajectories at the LGL nodes respectively:

Then, the trajectory of the state derivatives can be easily calculated with a simple matrix multiplication such that ẋ^N = x^ND^T to result in the ẋ^N matrix of dimensions n_x × (N + 1):

The matrix of the differential constraints, i.e. the RHS of the differential equations in the OCP formulation, of dimensions n_x × (N + 1) is calculated at each node as follows:

The matrix of differential defects at the collocation points, of dimensions n_x × (N + 1), is obtained from the subtraction of the differential constraints from the state derivatives:

The differential defects are introduced to enforce the differential constraints at the LGL nodes. A subsequent step consists in defining the matrix of path constraint function values evaluated at the LGL nodes as:

The objective function of the OCP is therefore approximated as follows:

where the weights ω_k are defined at the LGL nodes.

4 NLP formulation

In Non-linear Programming (NLP), a feasible solution of a system of equalities and/or inequalities is determined over a set of unknown variables, such that an objective function is minimized. In general, a NLP problem is defined as:

where y is the vector of decision variables of size n_y to be optimised with lower and upper bounds y_li and y_ui respectively. g_j (y) and h_k (y) are the r inequality and s equality constraints respectively. The OCP can be directly transcribed into a NLP problem by following a simple procedure. The objective function can be derived directly from equation (55), by constructing the decision vector y as a combination of the state vector x^N, the control vector u^N and the initial and final times t₀ and t_f respectively. This results in vector y having dimension n_y = n_x(N + 1) + n_u(N + 1) + 2 and this can be expressed as:

The decision vector is constrained by a lower bound vector y_li and an upper bound vector y_ui:

The equality constraint vector h_k (y) is constructed from the differential constraints such that:

with the dimension of the equality constraints vector r equal to n_x(N + 1). The equality constraints vector incorporates the path, event and time constraints with their respective bounds:

5 Generic BOOCP problem formulation

In aircraft trajectory optimisation, it is often the case that flight planners as well as aircraft designers and operators need to determine the optimal control law that minimizes the fuel consumption and the time of flight, or a compromise might be required between carbon dioxide and nitrogen oxides emissions. These are inherently BOOCPs. In the bi-objective case, the OCP is formulated similarly to the SOOCP in Section 2 with the exception for equation (1) which is now modified to cater of two objective functions:

where i = 1,2.

Following the pseudospectral discretisation and approximation processes described in Section 3, the objective functions are calculated as follows:

where i = 1,2. NLP formulation of the BOOCP will lead to the following problem:

Solving the BOOCP problem results in a Pareto set of optimal solutions. A standard NLP solver is limited to solving only single objective optimisation problems. Therefore, a scalarisation technique (also known as multi-objective articulation of preference techniques) needs to be used to reduce the BOOCP into a series of SOOCPs to populate the Pareto frontier. The three most applied scalarisation methods in the literature are the linear weighted sum (LWS), the normal boundary intersection (NBI) and the normal constraint (NC) methods. In previous works, the authors of this paper have proposed a new method, that of the adaptive bisection ε-constraint method which has been shown to perform better than the reviewed methods [30,37]. Hence, the adaptive bisection ε-constraint method was adopted to generate Pareto frontiers for the BOOCPs in this work.

6 Adaptive bisection ε-constraint method

The adaptive bisection ε-constraint method [37], is a scalarisation method conceived to efficiently generate optimal Pareto frontiers in BOOCPs. Similarly to other scalarisation methods, the Pareto frontier is generated by solving a sequence of single-objective optimisation problems in a systematic manner. In this method, one of the objective functions is selected to be optimised while the other is converted into an additional constraint, leading to a solution that can be proven to always be weakly Pareto optimal. Systematic modifications to the value of the objective function forming the additional constraint leads to the generation of an evenly distributed Pareto frontier. The method transforms the BOOCP into n_p SOOCPs, where n_p is user specified and it determines the target number of points on the Pareto frontier. The algorithm starts by obtaining the anchor points J1∗ and J2∗ of the BOOCP, corresponding to the minimum values of each of the performance indices through solving the two SOOCPs formulated in equation (67). The anchor points define the extremal points on the Pareto frontier, ensuring the ensuing scalarisation method does not fail to consider any part of the Pareto frontier in the optimisation process. The intersection of the lines J1=J1∗ and J2=J2∗ defines the utopian point J_u, which, albeit being an ideal solution, does not lie in the feasible region of the optimisation problem (Fig. 1). The remaining (n_p – 2) problems are formulated as follows:

with the additional constraint:

Fig. 1

Graphical representation of the design metric space of a BOOCP.

The value of ε_c for each of the remaining SOOCPs is calculated as follows. Once the the anchor points are determined, the utopian line, which is a straight line joining the anchor points, is bisected to obtain the first value of ε_c in equation (75). The solution of the resulting SOOCP will lead to an additional point on the Pareto frontier. Sometimes this will lead to an infeasible problem that leads to no solution. In this case, the line joining the anchor points is subdivided into four sections and the constraint ε_c is set to the value of J₂ at one-fourth the length of the line. If the problem is still infeasible, the value of J₂ at three-quarters of the line joining the anchor points is then tried. The line will continue being bisected until a solution is found or a constant K set by the user is reached. Once an additional point J3∗ is found, the euclidean distance between the point and other points on the Pareto frontier is determined. The two points with minimal euclidean distance are then used to find an additional point on the Pareto frontier by using the method of line bisection as previously described. This process is repeated until the number of Pareto points n_p requested by the user is found. For each point on the frontier, the optimal state and control trajectories can be easily deduced from the optimisation variables.

7 NLP solver implementation

A number of Nonlinear Programming solvers are publicly available, the predominant ones are SNOPT, IPOPT and fmincon. SNOPT is a software package for solving large-scale optimisation problems (linear and nonlinear programs) that was developed by Gill, Murray and Saunder [38]. It is implemented in Fortran 77 and distributed as source code against a license fee. IPOPT is an interior-point optimiser that was written by Wächter to solve large-scale nonlinear optimisation problems [39]. The algorithm is written in C++ and distributed open-source. Finally, fmincon is a function included in MATLAB’s optimisation toolbox which seeks the minimiser of a scalar function of multiple variables, with a region specified by linear constraints and bounds. In this work, both IPOPT and fmincon were adopted to demonstrate that the proposed method is non-solver-specific, as the intended application to avionics and ATM DSS will require the impementation in potentially very diverse platforms, such as ADA language on 64bit Acorn RISC Machine (ARM64) architectures. SNOPT was not used due to software licensing restrictions. The trajectory optimisation tool can be configured to use either fmincon or IPOPT. However, the latter is the preferred solver because in all test cases, IPOPT always converged faster than fmincon or converged in cases when fmincon did not find a solution. The execution time depends significantly on the length and complexity of the required trajectory.

8 Vertical trajectory optimisation case studiy

An aircraft trajectory optimisation problem was tackled to demonstrate the validity of the presented technique. The trajectory optimisation problem involves generating a Pareto frontier of optimal climb trajectories for an aircraft flying from 35 ft Above Mean Sea Level (AMSL) to a cruising altitude of 35,000 ft AMSL while covering a range of 900 kilometres. Two cost functionals were considered in this case: the minimization of flight time and the minimization of fuel consumption. Eurocontrol’s Base of Aircraft Data (BADA) coefficients of an Airbus A320 airliner were used to develop the aircraft dynamics model, following the methodology of Glover and Lygeros [41], described by the following state equations:

where X is the horizontal distance (range) covered by the aircraft on the ground, h is the altitude AMSL , V is the true airspeed (TAS) of the aircraft and m is the mass of the aircraft. X, h, V and m constitute the four state variables characterising the vertical profile of the aircraft. The control inputs are the flight path angle y, and the thrust ratio T_R, which is a fraction of the maximum thrust available from the engines at a particular altitude, T_max. The engine model integrated in the set of dynamic constraint to determine T_max as a function of altitude and airspeed is an empirical turbofan thrust model adopted by BADA. C_D is the coefficient of drag that varies with the aircraft configuration i.e. take-off, initial climb or clean. S represents the total area of the lifting surfaces on the aircraft, ρ is the air density at a particular value of h assuming standard atmospheric conditions and f is the rate of fuel consumption. Finally, g is the gravitational acceleration assumed to be constant. The initial and final conditions of the aircraft, formulated as event constraints in the BOOCP problem are defined in Table 1. Representative path constraints were applied to the trajectory optimisation problem to reproduce realistic limits on the performance of the aircraft and its propulsion system. The first path constraint is the stall speed of the aircraft which is a function of the aircraft configuration and constitutes the minimum value for the airspeed. The upper end of the speed scale is limited by V_MO and M_MO which are the maximum operating speed and the maximum operating Mach number respectively. Finally, the maximum thrust provided by the engines is also formulated as a path constraint.

The Pareto frontier of this climb case study is presented in Figure 2. The frontier is evenly distributed and illustrates the extremal solutions at the far ends of the frontier, as well as intermediate solutions that provide a compromise between fuel consumed and flight time for a particular trajectory. The altitude-range profiles of the complete Pareto set are illustrated in Figure 3 as a three-dimensional surface plot. The trajectory at the front of the surface is the minimum time trajectory, whereas the profile at the far back represents the trajectory which consumes the least fuel. The complete surface is filled with successive plotting of trajectories representing the rest of the points in the Pareto optimal set. The discussion henceforth will focus on the extremal solutions.

Fig. 2

Pareto frontier for minimum time and minimum fuel burn.

Fig. 3

Altitude profiles for Pareto set.

Figure 4 illustrates the altitude-range profile for the two extremal solutions. A total of 4.36 tonnes of fuel are consumed in the minimum time trajectory, which is flown in 1 hr 2 mins. On the other hand, 3.49 tonnes of fuel and a flight time of 1 hr 14 mins characterise the minimum fuel trajectory. This means that for this specific problem, an increased flight time of 12 mins over a 900 km leg can be traded for approximately 20% of the fuel consumed.

Fig. 4

Altitude profile for minimum time and minimum fuel burn.

The minimum fuel trajectory involves a direct climb to the cruising altitude at which the aircraft flies the rest of the leg. This behaviour is expected since modern jetliners are most efficient at high altitudes where the air density is lower, which results in less drag on the aircraft. However, it is worth noting that the upper bound on the altitude in the BOOCP formulation was that of 45,000 ft giving the aircraft the possibility to fly higher than 35,000 ft in the level flight phase followed by a descent to the requested altitude at the end of the leg. This did not materialise, however, since the cost of attaining higher altitudes was larger than the savings in fuel yielded by flying a few flight levels higher. It is interesting to note that the minimum time trajectory follows a completely different strategy. The flight starts with a slightly shallower climb up to around 24,500 ft, followed by a level flight and finally an additional shallower climb to the requested altitude of 35,000 ft. The altitude flown for most of the time by the aircraft is the cross-over altitude i.e. the altitude at which the airspeed corresponding to the maximum operating Mach number M_MO coincides with the maximum operating speed V_MO. Consequently, this happens to be the altitude at which the aircraft can fly at the maximum true airspeed, which results in minimum flight time. In both cases, the climb from 35 ft occurs at the maximum thrust provided by the engine as can be seen in Figure 5. In the minimum time trajectory it is clear that engine thrust is kept at the maximum level permissible by operational constraints to gain the maximum possible airspeed which will result in minimal flight time. When minimising fuel costs, however, this is less obvious because fuel consumption is directly related to thrust levels. Reducing the thrust during the climb while keeping a constant flight path angle will result in a decrease in the rate of climb. As a result, it will take longer to reach the target altitude. Since the total fuel consumed during the climb is the integral of the rate of fuel consumption for the duration of the trajectory, the longer the flight time, the larger the fuel consumption. The optimisation results suggest that in the absence of winds it is more convenient to climb at high thrust levels (and high rates of fuel consumption) at high rates of climb rather than by simply keeping low thrust levels for a prolonged time. This is also understandable in the context that, for a given calibrated airspeed (CAS), TAS increases disproportionately with altitude. In still air, such as the case considered, TAS is what defines the ground speed and hence the flight time. As a result, it is expected that in the absence of winds it will be always advantageous to expedite climb in order to reach higher altitudes quicker.

Fig. 5

Thrust ratio profile for minimum time and minimum fuel burn.

The true airspeed and flight path angle (FPA) profiles illustrated in Figures 6 and 7 correlate perfectly with the observations previously made. The minimum time trajectory yields a speed profile that is larger in magnitude throughout the whole flight range over the minimum fuel counterpart. The former adopts a TAS of 495 kts, equivalent to 0.86 Mach, at the crossover altitude, while the latter flies at a nearly constant TAS varying from 419 kts to 410 kts at 35,000 ft. This result is consistent with the fact that during cruise, the aircraft would be expected to climb or slow down gradually as the aircraft becomes lighter through progressive fuel burn. In this case, the latter strategy was adopted. The steep climb of the minimum fuel trajectory is reflected in Figure 7 in which a maximum FPA of 6.2° is adopted as opposed to a maximum FPA of 4.4° in the minimum time case. The minimum time trajectory exhibits a lower climb gradient than its counterpart in order to afford a quicker TAS during the climb, since in both cases the engine thrust is set to the maximum level. Furthermore, in the minimum time trajectory, the second climb at the end of the trajectory, which is shallower than the first with a maximum FPA of 1°, is explained by the fact that at higher altitudes turbine engines generate less thrust due to the thinner atmosphere. Moreover, no trade-off is occurring between kinetic and potential energy to provide the additional rate of climb until the very end of the flight.

Fig. 6

True Airspeed profile for minimum time and minimum fuel burn.

Fig. 7

Flight path angle profile for minimum time and minimum fuel burn.

In Figures 5 to 7, it can be observed that the end of the trajectory is salient. In particular, in the minimum time profile the last 5,000 ft are climbed in just 12 km. Clearly, this is not practical for both operational and safety reasons. The speed profile during this time indicates a rapid decrease in TAS from 495 kts to 363 kts. The speed change occurs due to the trading of kinetic energy (speed) with potential energy (altitude), which is operationally defined zoom-climb [43], and also due to a rapid reduction in the thrust ratio. This zoom-climb solution falls within the feasible search space of the optimiser since no final condition is set on the speed in the problem formulation. A similar but less extreme shift is adopted in the minimum fuel trajectory. In this case, the aircraft is already flying at the cruise altitude at fuel-optimal speed. Since fuel consumption is of primary importance, in the final part of the trajectory the thrust is reduced to zero to conserve fuel leaving the aircraft to decelerate slowly towards the stall speed due to the aerodynamic drag. Theoretically this is possible since no constraints are included in the problem formulation even though for practical reasons such a solution is not operationally practical.

9 Combined later and vertical trajectory optimisation case study

As mentioned in section 1, the potential benefits of implementing trajectory optimisation techniques in avionics and ATM systems has recently met a growing interest by the research community. This case study therefore addresses the online planning of optimal 4DT intents in representative ATM scenarios. In particular, the spacing of dense arrival traffic towards in the Terminal Manoeuvring Area (TMA) was evaluated as a representative case study of online tactical 4DT planning. For trajectory optimisation of medium-large transport aircraft, Three Degrees of Freedom (3-DOF) point-mass aircraft dynamics are very commonly adopted. The 3-DOF implemented in our simulation case study is

where the state vector consists of: longitudinal velocity v [m s^–1], flight path angle y [rad]; track angle 𝒳 [rad]; geographic latitude ϕ [rad]; geographic longitude λ [rad]; altitude z [m]; thrust angle of attack ∊ [rad]; aircraft mass m [kg]; and the control vector includes: thrust force T [N]; load factor N [ ]; bank angle μ [rad]. Other variables and parameters include: aircraft weight W and aerodynamic drag D [N]; wind velocity v_w in its three scalar components [m s^–1]; gravitational acceleration g [m s^–2]; Earth radius R_E [m]; fuel flow FF [kg s^–1]. For turbofan aircraft, the following empirical expressions were adopted in the development of Eurocontrol’s Base of Aircraft Data (BADA), to determine the climb thrust and the fuel flow FF, which operationally equates to the maximum thrust T_MAX in all flight phases excluding take-off [41]:

where τ is the throttle control, H_P is the geopotential pressure altitude in feet, ΔT is the deviation from the standard atmosphere temperature in kelvin, v_TAS is the true airspeed. C_T1… C_T5, C_f1… C_f4 are the empirical thrust and fuel flow coefficients, which are also supplied as part of BADA for a considerable number of currently operating aircraft [41]. While carbon dioxide (CO₂) emissions are characterised by an approximately constant emission index of EI_CO2 = 3.16 [Kg/Kg], an empirical model for carbon monoxide (CO) and unburned hydrocarbons (HC) emission indexes (EI_CO/HC) in [g/Kg] at mean sea level based on nonlinear fit of experimental data from the ICAO emissions databank is:

where the fitting parameters c_{1,2, 3} accounting for the average emissions of 165 currently operated civil turbofan engines are c = {0.556, 10.208, 4.068} for CO and c = {0.083, 13.202, 1.967} for HC [23]. The nitrogen oxides (NO_X) emission index [g/Kg] based on the curve fitting of 177 currently operated civil aircraft engines is [23]:

In the AMAN scenario, a best arrival sequence is defined by suitable scheduling algorithms. Longitudinal separation is enforced at the final approach fix to ensure sufficient separation upon landing and to prevent separation infringements in the approach phase itself. The assumed minimum longitudinal separation is 4 nmi 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 aircraft. The pseudospectral trajectory optimisation algorithm is therefore executed targeting a minimum fuel. The results of one representative simulation run are depicted in Fig. 8.

Fig. 8

Results of the pseudospectral 4D trajectory optimisation algorithm applied to the online TMA trajectory planning scenario.

Computation times consistently below 20 seconds were achieved by our simulated ATM system to plan each 4DT intent of approximately 40 nmi with multiple path constraints in a representative TMA arrival and departure traffic conditions. Validated 4DT intents for multiple conflicting arrival traffic in a complex AMAN Terminal Sequencing and Spacing (TSS) scenario were produced in an average time of 41 seconds and consistently less than 60 seconds in a Monte Carlo simulation [25].

10 Conclusion and future work

This paper has presented a new computational technique to solve multi-objective aircraft trajectory optimisation problems formulated as optimal control problems. Due to their computational efficiency, pseudospectral methods were implemented to transcribe the formulated nonlinear optimal control problem. In terms of multi-objective optimality, the adaptive bisection ε-constraint scalarisation method was adopted to yield a very efficient generation of Pareto frontiers for a variety of bi- and multi-objective aircraft trajectory optimisation problems that are being currently investigated as part ofmajor aviation modernisation programmes. This multi-objective technique demostrates the required computational efficiency for implementation into avionics and Air Traffic Management (ATM) systems to generate optimal 4-Dimensional Trajectories (4DT) addressing multiple operational, economic and environmental criteria. A bi-objective vertical flight trajectory optimisation case study was evaluated in detail, involving a climb from 35 ft to a cruise altitude of 35,000 feet covering a range of 900 km and including operationally representative path constraints in terms of airspeed and engine thrust. The strategies adopted by the optimisation process are consistent with flight performance theory, therefore validating the algorithm. In particular, to maximise fuel efficiency in the absence of winds it is preferable to attain the cruise flight level as early as possible. Conversely, for maximised time performance it is preferable to interrupt the climb and maintain the crossover altitude as long as possible, then exploit the excess in kinetic energy to quickly reach the final cruise altitude. The proposed computational technique was also applied to the optimisation of 4DT intents in realistic ATM contexts. This work forms the basis of a trajectory optimisation tool that will be used to solve more complex and challenging problems in the future. The ultimate goal of the ongoing research is to develop innovative real-time trajectory optimisation methods which have the potential to be used in the next generations of ATM 4DT intent planning/negotiation systems and avionics Flight Management Systems (FMS). Trajectory optimisation algorithms based on pseudospectral discretisation and efficient Non-Linear Programming (NLP) solvers also have the potential to meet the real-time performance requirements for emergency avoidance such as in the case of obstacle warning systems and Sense-and-Avoid (SAA) systems [44–46]. Ongoing research activities are investigating the feasibility of avionics/ATM technologies capable of contributing to the emissions reduction targets set by ICAO and by European, American and Australian governments. In this respect, significant opportunities are emerging in relation with the adoption of multiobjective optimal control-based algorithms for trajectory planning and management considering multiple operational, economic and environmental criteria.