Next Article in Journal
Dynamic Burst Actuation to Enhance the Flow Control Authority of Plasma Actuators
Next Article in Special Issue
Fault Tolerant Attitude and Orbit Determination System for Small Satellite Platforms
Previous Article in Journal
Strength Evaluation and Failure Analysis of the Vortex Reducer under Overspeed Condition
Previous Article in Special Issue
Iterative Lambert’s Trajectory Optimization for Extrasolar Bodies Interception
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Objective and Multi-Phase 4D Trajectory Optimization for Climate Mitigation-Oriented Flight Planning

Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy
*
Author to whom correspondence should be addressed.
Aerospace 2021, 8(12), 395; https://doi.org/10.3390/aerospace8120395
Submission received: 16 October 2021 / Revised: 29 November 2021 / Accepted: 7 December 2021 / Published: 13 December 2021
(This article belongs to the Special Issue Aerospace Guidance, Navigation and Control)

Abstract

:
Aviation contribution to global warming and anthropogenic climate change is increasing every year. To reverse this trend, it is crucial to identify greener alternatives to current aviation technologies and paradigms. Research in aircraft operations can provide a swift response to new environmental requirements, being easier to exploit on current fleets. This paper presents the development of a multi-objective and multi-phase 4D trajectory optimization tool to be integrated within a Flight Management System of a commercial aircraft capable of performing 4D trajectory tracking in a Free Route Airspace context. The optimization algorithm is based on a Chebyshev pseudospectral method, adapted to perform a multi-objective optimization with the two objectives being the Direct Operating Cost and the climate cost of a climb-cruise-descent trajectory. The climate cost function applies the Global Warming Potential metric to derive a comprehensive cost index that includes the climate forcing produced by CO2 and non-CO2 emissions, and by the formation of aircraft-induced clouds. The output of the optimization tool is a set of Pareto-optimal 4D trajectories among which the aircraft operator can choose the best solution that satisfies both its economic and environmental goals.

1. Introduction

Air traffic provides a significant contribution to anthropogenic global warming. In 2005, aviation contributed about 5% of the overall anthropogenic radiative forcing (RF), with the largest warming contributions being due to CO2 emissions, contrail cirrus formation, and NOX net effects [1]. It is possible to reduce the environmental impact of aviation by intervening at the Air Traffic Management (ATM) level. The Single European Sky ATM Research (SESAR) framework has introduced the concept of Free Route Airspace (FRA) that allows the free planning of a route between any defined entry and exit point in the airspace, enabling shorter routes and a more efficient use of the airspace [2].
A partial implementation of the FRA concept in the European airspace from 2014 to 2018, limited to the upper airspace (above FL335 [3]) and to a portion of the European airspace, with 20% of flight time flown in free routing in 2017, produced an estimated saving of 2.6 million tons of CO2 [4]. Current flight planning in FRA sections is limited to a reduction of route track path distance, and aims to reduce route extension (RE), i.e., the difference between flight and great circle distance between two waypoints. A full implementation of FRA will further reduce economic and environmental costs by enabling the definition of 4D flight profiles during flight planning, to be optimized with up-to-date weather information to define wind-optimal lateral paths and vertical profiles that exploit cruise climb techniques in lieu of cruise at level flight [5]. In this context, airborne trajectory optimization functionalities integrated within avionic systems will be a cornerstone of the network-centric ATM that will reduce Air Traffic Control (ATC) workload during flight plan definition [6].

1.1. 4D Trajectory Optimization

The trajectory optimization problem can be formulated as an Optimal Control Problem (OCP), commonly solved using numerical methods. The main approaches used to numerically solve OCPs are direct methods, indirect methods, and dynamic programming methods [7,8]. Direct methods are extensively used for aircraft trajectory optimization since they are easier to apply than indirect methods that require the derivation of the analytical expressions of the necessary conditions [7], and more efficient than classical dynamic programming methods that suffer from the “curse of dimensionality” and require the use of expert knowledge to reduce the solution search space [9,10]. In particular, direct collocation pseudospectral methods are characterized by a fast convergence to the optimal solution [11], and can be used to solve multi-phase problems [12,13]. Multi-phase optimization allows the optimization of an entire trajectory, composed of different flight phases corresponding to different aircraft modes and configurations that lead to a variation of the dynamic model [14,15]. Direct methods are at the base of many commercial off-the-shelf and general-purpose optimization tools [16,17,18]. However, when integrating 4D trajectory optimization functionality within the Flight Management System (FMS) of an aircraft, it is preferable to have a customizable and dedicated tool over a black-box solution.
The 4D trajectory optimization tool developed in this study will be integrated in the FMS of a Boeing 747-100, implemented in the Multipurpose Aircraft Simulation Laboratory (MASLab) developed within the Clean Sky research program, and able to execute 4D trajectory tracking [19,20,21]. The mathematical model of the aircraft used in the optimization tool uses Eurocontrol’s base of aircraft data (BADA) family 3 performance model for the propulsion system [22], while the aerodynamic model is based on Boeing 747-100 data gathered by Hanke and Nordwall [23]. The optimization algorithm is based on the Chebyshev pseudospectral direct method presented by [24], modified to allow sequential multi-phase optimization through a set of continuity boundary conditions of the states between phases and an estimation of the trajectory cost of subsequent phases in the cost function. The initial guess of the optimization algorithm is generated automatically: it is based on the orthodrome trajectory, flown at maximum specific range configuration with the respect to cost function that links the initial and final waypoints of the flight arc provided in input by the user. The optimization tool performs multi-objective optimization, using a weighted squared sum method, where the two conflicting objectives are the Direct Operating Cost (DOC) and the environmental cost of the trajectory.

1.2. Climate Mitigation-Oriented Flight Planning

Currently available 4D flight planning solutions offer the possibility to optimize routes with respect to fuel consumption, operating costs or flying time [25]. Fuel-optimal trajectories guarantee minimum CO2 emissions; however, this does not imply a minimization of the environmental impact. In fact, CO2 is the main contributor to anthropogenic global warming [26], yet the largest net warming effect produced by aviation is due to persistent contrails and contrail cirri, i.e., aviation-induced clouds (AIC), followed by CO2 and NOX [27]. In particular, 66% of the aviation net equivalent radiative forcing (ERF) in 2018 is due to non-CO2 aviation climate forcings [27]. The three major aviation environmental impacts are climate impact, noise pollution, and air quality [28]. Noise pollution and air quality are generally addressed during initial and terminal phases of flight, as well as in the ground phases, where the aircraft is closer to highly populated areas [29,30]. Climate impact, defined as the sum of CO2 and non-CO2 net warming effect produced by the aircraft along the trajectory, is mainly addressed during cruise phase. Moreover, cruise is typically the only phase of flight concerned by AIC, since they form in the ice supersaturated regions (ISSR) of upper airspace [31]. For this reason, much research work focused on the reduction of AIC formation during cruise [32,33,34,35,36], facing the problem of a trade-off between contrail formation and the increased fuel consumption (with the ensuing increase in CO2 emissions) necessary to avoid ISSRs. Other studies developed climate cost functions that exploit CO2-equivalency metrics based on the ERF of aviation climate forcings to aggregate CO2 and non-CO2 net warming contributions into a single environmental cost index, in order to assess the mitigation potential of climate-optimal trajectories in the commercial airspace and their cost when compared to fuel-optimal trajectories [37].
This paper presents the development of a novel trajectory optimization tool for an effective climate mitigation-oriented 4D flight planning. Climate cost has been assumed as the environmental cost index that is provided as a choice among the other traditional cost indexes to be minimized during the definition of the flight plan (i.e., operating cost, fuel consumption, flight time). Multi-objective optimization with respect to DOC and environmental costs leads to a set of Pareto-optimal trajectories: in this way, it is possible to select, among multiple choices, the best flight plan that can satisfy both DOC and environmental goals of the aircraft operator. Multi-phase optimization allows the definition of the optimal climb-cruise-descent arc of the flight, i.e., the flight phases where it is possible to capitalize the benefits of FRA for climate mitigation. The environmental cost index is calculated in terms of CO2-equivalent emissions of greenhouse gases (GHGs) and AICs produced by the aircraft during the climb-cruise-descent arc of the trajectory. The contributions of the individual climate forcings are weighted using the Global Warming Potential (GWP) emission metric integrated over a period of 20, 50, and 100 years [38]. Aircraft emissions are calculated using the Boeing Fuel Flow Method2 (BM2) [39], adopting engine-specific emission data from the International Civil Aviation Organization (ICAO) engine exhaust emissions databank [40]. Atmosphere and weather data are based on the Global Forecast System (GFS) model and provided by the National Oceanic and Atmospheric Administration (NOAA). The AIC formation model is built upon the Schmidt-Appleman criterion (SAC) for contrail formation, combined with a ISSR detection model [41]. The output of the multi-objective optimization tool is a set of Pareto-optimal 4D trajectories, i.e., a set of trajectories where each of the two objective costs, DOC and environmental cost, cannot be improved without worsening the other [42]. The extremes of the Pareto-optimal set of solutions are the DOC-optimal trajectory, with the higher environmental cost among the Pareto solutions, and the environmental-optimal trajectory, with the higher DOC.
The rest of the paper is organized as follows: Section 2 presents the trajectory optimization tool developed within this work, describing the optimization algorithm, the mathematical model of the aircraft, the atmospheric and AIC prediction model, the environmental cost model, the definition of the multi-objective cost function, the multi-phase optimization process, the initial guess trajectory generation; Section 3 shows the results for a test case based on a transatlantic flight from Rome to New York; Section 4 draws the conclusions about the work and suggests future improvements.

2. Materials and Methods

2.1. Direct Trajectory Optimization

The optimal trajectory is obtained performing a direct trajectory optimization using a Chebyshev pseudospectral method [24]. The trajectory optimization problem can be formulated as a general optimal control problem (OCP):
min x ( t ) , u ( t ) , t f J [ x ( t ) , u ( t ) , t ] = M [ x ( t f ) , t f ] + t 0 t f L [ x ( t ) , u ( t ) , t ] d t
subject to : { (1b) f l f [ x ˙ ( t ) , x ( t ) , u ( t ) , t ] f u (1c) b l b [ x ( t 0 ) , x ( t f ) , t 0 , t f ] b u (1d) c l c [ x ( t ) , u ( t ) , t ] c u
where x ( t ) and u ( t ) are the states and control trajectories defined in the time t [ t 0 , t f ] . The solution of the OCP is given by the control function u ( · ) and the corresponding trajectory x ( · ) that minimize the Bolza cost function defined in (1a), and composed by the Mayer term M ( · ) and the Lagrange term L ( · ) . The Problem (1) is subject to the dynamic constraints given in Equation (1b), the boundary constraints b ( · ) in Equation (1c), and the algebraic path constraints c ( · ) in Equation (1d).
Setting the lower and upper bounds equal to zero for the dynamic constraints in Equation (1b) leads to:
f [ x ˙ ( t ) , x ( t ) , u ( t ) , t ] = 0
and assuming that the Jacobian f / x ˙ is nonsingular, the dynamic constraint can be represented by the controlled ordinary differential equation:
x ˙ ( t ) = f [ x ( t ) , u ( t ) , t ]
The OCP may be solved using both direct or indirect solution methods: direct methods have the advantage of not needing to solve the necessary conditions derived from the Pontryagin maximum principle [24]. To apply direct methods, the OCP problem is discretized into a finite dimension nonlinear programming (NLP) problem: the domain of the independent time variable is subdivided into a finite number of intervals and the controls and states at the discretization points are evaluated through polynomial approximation. In direct collocation methods, both control and state variables at the nodes are the unknowns of the problem.
The Chebyshev pseudospectral method uses Chebyshev–Gauss-Lobatto (CGL) nodes as discretization points for the time variable. The CGL node points lie in the interval defined by the time domain τ [ 1 , 1 ] . Using a linear transformation, the problem is transformed from the generic interval t [ t 0 , t f ] to τ [ τ 0 , τ f ] = [ 1 , 1 ] :
t ( τ ) = ( t f t 0 ) τ + ( t f + t 0 ) 2
resulting in the following reformulation of the OCP in (1):
min x ( t ) , u ( t ) , t f J [ x ( t ) , u ( t ) , t ] = M [ x ( 1 ) , t f ] + t f t 0 2 1 1 L [ x ( t ) , u ( t ) , t ( τ ) ] d τ
subject to : { (5b) f l f 2 t f t 0 x ˙ ( τ ) , x ( τ ) , τ ( t ) , t ( τ ) f u (5c) b l b [ x ( 1 ) , x ( 1 ) , t 0 , t f ] b u (5d) c l c [ x ( τ ) , u ( τ ) , t ( τ ) ] c u
where for conciseness, x ( τ ) stands for x [ τ ( t ) ] .
The N + 1 CGL interpolation points are expressed in closed form as:
τ k = cos ( π k / N ) , k = 0 , , N
Given the N-th order interpolating Chebyshev polynomials of the first kind, defined as:
T j ( τ ) = cos ( j arccos ( τ ) ) , j = 0 , , N
for τ [ 1 , 1 ] , which yields to
T j ( τ k ) = cos ( π k j / N )
the interpolation points τ k are the extrema of the interpolating polynomials, or the roots of their derivatives T ˙ j ( τ ) .
The obtained CGL nodes are sorted in an ascending order and cluster in the proximity of the endpoints of the [ 1 , 1 ] interval, eliminating the Runge’s divergence for increasing values of N that occurs with equally spaced points [43].
The state and control continuous equations are obtained through polynomial approximation with the form:
x N ( τ ) = j = 0 N x j ϕ j ( τ )
u N ( τ ) = j = 0 N u j ϕ j ( τ )
where ϕ ( τ ) are Lagrange polynomials of order N interpolating function at the CGL nodes:
ϕ j ( τ ) = ( 1 ) j + 1 N 2 c j ( 1 τ 2 ) T ˙ N ( t ) τ τ j
with
c j = 2 , j = 0 , N 1 , 0 < j < N
The interpolating polynomials show the Kronecker delta property:
ϕ j ( τ k ) = δ j k = 1 , if j = k 0 , if j k
The derivative of the states at the CGL node points can be expressed differentiating Equation (9):
d k = x ˙ N ( τ k ) = j = 0 N x j ϕ ˙ j ( τ k ) = j = 0 N D k j x j
where the differentiation matrix D is an expression in compact form of the Lagrange polynomials derivative and can be calculated as [44], adjusted for taking into account the ascending order of the interpolation points τ k :
D : = [ D k j ] : = { c k c j ( 1 ) j + k τ k τ j , j k τ k 2 ( 1 τ k 2 ) , 0 < j = k < N 2 N 2 + 1 6 , j = k = 0 2 N 2 + 1 6 , j = k = N
The definite integral term of the continuous cost function can be discretized by approximating it to a finite sum using the Clenshaw–Curtis quadrature scheme. For a given function g:
1 1 g ( τ ) d τ k = 0 N g ( τ k ) w k
where w k are the optimal quadrature weights and can be calculated as follows [45]:
w k = c k N 1 j = 1 N / 2 b j 4 j 2 1 cos 2 j k π N , k = 0 , , N
where the floor function N / 2 returns an integer for any value of N and the coefficients b j and c j are defined as:
b j = 1 , j = N / 2 2 , j < N / 2 , c k = 1 , k = 0 , N 2 , 0 < k < N
Equation (15) holds for any N > 1 .
The OCP is finally approximated as follows:
min X , U , t f J [ X , U , t f ] M [ x N , t f ] + t f t 0 2 k = 0 N L [ x k , u k , t k ] w k
subject to : { (17b) f l f 2 t f t 0 d k , x k , u k , t k f u , k = 0 , , N (17c) b l b [ x 0 , x N , t 0 , t N ] b u (17d) c l c [ x k , u k , t k ] c u , k = 0 , , N
where the coefficients
X = ( x 0 , x 1 , , x N ) , U = ( u 0 , u 1 , , u N )
and the final time t f are the unknown of the OCP that minimize the cost function, and are found using interior-point methods for NLP.

2.2. Aircraft Model

2.2.1. State Equations

The mathematical model of the aircraft is based on a conventional three degrees-of-freedom (3-DoF) rigid-body model with variable point-mass, expressed as a set of differential algebraic equations (DAE). The model is formulated assuming that the aircraft flies on the optimal trajectory in steady wings-level flight, adjusting the lateral path through steady turning flight with negligible turning rate, and performing major flight path level adjustments only during transitions to and from cruise with a negligible pull-up rate. The resulting 3-DoF flight dynamics DAE is given by the following state equations:
λ ˙ = V cos γ cos ψ + w N R E + h
ϕ ˙ = V cos γ sin ψ + w E ( R E + h ) cos λ
h ˙ = V sin γ
m ˙ = T · T S F C
V ˙ = T D m g sin γ
Equations (18) and (19) are the equations of motion on the lateral plane for latitude λ and longitude ϕ of the World Geodetic System (WGS84), and R E = f ( λ ) is the radius of Earth as a function of the latitude. The w N and w E terms represent respectively the North and East components of wind speed and are a function of latitude, longitude and altitude h. Equation (20) is the vertical equation of motion as a function of the flight path angle γ . Equation (21) is a function of the thrust T and the thrust specific fuel consumption T S F C . Equation (22) is the state equation for the true airspeed V and is a function of thrust T and aerodynamic drag D.
The control variables of the OCP are the heading angle ψ , the flight path angle γ , the throttle Π that controls the thrust: T = f ( Π ) .

2.2.2. Performance Model

The aerodynamic model of the aircraft is based on the one built by Hanke and Nordwall for the NASA Flight Simulator For Advanced Aircraft (FSAA) of the Boeing 747-100 [23]. The model interpolates gridded data to calculate the aerodynamic drag coefficient C D as a function of the lift coefficient C L and the Mach.
The lift coefficient C L is calculated as:
C L = 2 m g cos γ ρ V 2 S
The propulsive model of the aircraft is based on performance data from Eurocontrol’s base of aircraft data (BADA), family 3 [22]. The throttle Π controls the thrust T [ T i d l e , T m a x ] with the following linear law:
T ( Π ) = T i d l e + Π ( T m a x T i d l e )
where the minimum idle and maximum thrust available are a function of the flight conditions and BADA coefficients specific to the aircraft. The fuel flow F F is then calculated as a function of the thrust and the thrust specific fuel consumption T S F C :
F F = T · T S F C
The T S F C is defined in BADA as a function of engine-specific coefficients and true airspeed [22].

2.2.3. Emissions Model

Emissions are calculated to estimate the climate impact of the trajectory. The mass of a generic pollutant m g p emitted during a given flight time can be calculated as
m g p = t 0 t f F F · E I g p d t
where the E I g p is the emission index of the given pollutant and represents the mass of emission per mass of fuel burned. Emission indexes are available for carbon dioxide CO2, water vapor H2O, sulfur dioxide SO2, soot, nitrogen oxides NOX, unburned hydrocarbons HC, and carbon monoxide CO. The emissions indexes in Table 1 have been assumed to be constant throughout the trajectory, hence the mass of emission is proportional to the mass of fuel burned.
The emissions indexes of NOX, HC, and CO vary as a function of the fuel flow and the atmospheric conditions and are calculated using the Boeing Fuel Flow Method2 (BM2) [39] using emissions data from the International Civil Aviation Organization (ICAO) engine exhaust emissions databank [40]. These data, provided to ICAO by the engine manufacturers, consist of coupled values of a relative emission index (REI) and a fuel flow factor W F F at four power settings, i.e., idle/taxi (7%), approach (30%), climb (85%), and take-off (100%).
The fuel flow factor W F F is defined as:
W F F = F F ϑ a m b 3.8 δ a m b exp ( 0.2 M 2 )
where ϑ a m b and δ a m b are respectively the ratios between ambient temperature and pressure with respect to the sea level values.
Nonlinear regression with least square error can be applied to ICAO engine data to derive the following formulation of NOX REI as a function of the fuel flow factor W F F :
R E I N O X ( W F F ) = 6.5991 · W F F 2 + 3.1542 · W F F + 2.5003
The emission index for NOX can be finally obtained using the following relation:
E I N O X = R E I N O X · exp ( H ) · δ a m b 1.02 ϑ a m b 3.3 0.5
where H is a humidity factor that depends on ambient pressure and relative humidity [39].

2.3. Atmosphere Model

Atmosphere data are provided by the National Oceanic and Atmospheric Administration (NOAA) through the NOAA Operational Model Archive and Distribution System (NOMADS) based on the Global Forecast System (GFS) model. Weather data are gathered in lookup tables with a base horizontal resolution of 0.25 degrees and a mean vertical pressure altitude resolution of 25 millibars. Gridded weather data downloaded from NOMADS are the north and west wind velocity components w N and w E in m/s, ambient temperature T a m b in K, relative humidity R H defined by values between 0 and 1.
Pressure altitude QNE expressed in meters h and ambient pressure p in Pascals are directly related as in [46]:
h = 1 p 101325 0.190284 · 44,307.694
The density for moist air ρ a is derived from the CIPM-2007 formula [47]:
ρ a = p M a Z R T a m b 1 x v 1 M v M a
where M a and M v are the molar masses of dry air and water vapor, x v is the molar fraction of water vapor, Z is the compressibility factor and R is the molar gas constant in J mol−1 K−1.

Sub-CGL Grid

The spatial resolution of the GFS data grid, equal to 0.25 degrees, has an order of magnitude of ∼ 10 1 km, whereas the cruise phase of a long-haul flight spans over an arc with a length of ∼ 10 3 km. Consequently, if the nodes of the CGL grid are lower than N∼ 10 2 it means that the GFS grid resolution is higher than the CGL resolution used for pseudospectral optimization. In these cases, all the information about the atmospheric conditions between CGL nodes is lost. This problem is solved by creating a sub-CGL grid, with equidistant sub-nodes between CGL points that have a spatial resolution comparable to the one of the GFS grid, and lay on the geodesic trajectory that links consecutive CGL nodes, as shown in Figure 1. All the other state and control variables in the arcs between CGL nodes are approximated using linear interpolation.
The actual atmosphere variable ω ¯ attributed to a specific CGL node in pseudospectral optimization is calculated as follows. Defining P the set of N + 1 CGL nodes p, and S the set of the M nodes s of the sub-CGL grid, S n is defined as the set of s nodes that are closer to a given p n node, with p [ 0 , N ] than to any other CGL node. The atmosphere variable attributed to the p n CGL node is calculated by averaging the atmosphere GFS data for the coordinates of the nodes in the S n set of | S n | elements:
ω ¯ n = i = 1 | S n | ω ( s i ) | S n |
where ω is a generic atmospheric variable gathered from the GFS gridded data as a function of spatial coordinates of node s i .

2.4. Aircraft-Induced Clouds Model

With the performance data of the aircraft model and the atmospheric model it is possible to detect the formation of persistent contrails and contrail cirri, to evaluate their impact in the environmental cost.
According to Schmidt-Appleman criterion (SAC) condensation trails form if the atmospheric temperature is lower than a threshold temperature T L C [41]. Short-lived contrails dissolve in seconds to a few minutes, with little impact on Earth’s radiative forcing. Contrails that last for at least 10 min are defined as persistent contrails by the World Meteorological Organization [48]. Persistent contrails only form and evolve in ice supersaturated regions (ISSR), where the relative humidity over ice RHi is greater than 100%:
R H i = R H p l i q p i c e
where p l i q and p i c e are the saturation pressures over liquid and ice water surfaces and can be calculated using the equations in [49].
Contrail cirri form when ISSR cold and moist conditions further influence the radiative properties of the persistent contrail, along with its shape, size, and duration, making them almost indistinguishable from natural ice clouds [31].
Persistent contrails and contrail cirri altogether are defined as aircraft-induced clouds (AIC), and form if both SAC and ISSR conditions are true:
{ T a m b T L C R H i = R H p l i q p i c e 1
The threshold temperature for the SAC criterion T L C is calculated by applying the formulation in [41]:
T L C = T L M ( 1 R H ) p l i q ( T L M ) G Δ T c
where all the temperatures are expressed in °C, Δ T c is a correction factor that is equal to 0 for relative humidity R H = { 0 , 1 } and positive otherwise; T L M is the maximum threshold temperature for R H = 1 and is given by:
T L M = 46.46 + 9.43 ln ( G 0.053 ) + 0.72 [ ln ( G 0.053 ) ] 2
Finally, G is expressed in Pa K−1 and is given by:
G = c p p a m b E I H 2 O ( M v / M a ) Q f ( 1 η )
with specific heat capacity of air c p = 1004 J, Q f is fuel specific energy in J kg−1, and η is the overall propulsion efficiency of the aircraft.

2.5. Multi-Objective Cost Functional

2.5.1. Direct Operating Cost

The direct operating cost (DOC) is the monetary cost of the flight and is directly impacted by duration and fuel consumption:
J D O C = t 0 t f ( C t + F F · C f ) d t
where C t is the estimated direct operating cost per unit time excluding fuel and C f is the cost of fuel per unit mass. The cost coefficients used within the optimization tool are based on data for the year 2018, and are C t = 0.5381 $/s [50], and C f = 0.7152 $/kg [51]. The DOC is the Mayer term of the cost function of the optimization problem as defined in Equation (1a), since C t and C f are constant, hence J D O C depends only on final time and states and can be rewritten as:
J D O C = C t ( t f t 0 ) + C f ( m f m 0 )

2.5.2. Environmental Cost

The environmental cost function is composed by the climate cost and measures the contribution of a given trajectory to global warming. Aviation contributes to anthropogenic global warming through the emissions of different greenhouse gases (GHGs) and other climate forcing agents such as AIC. Emission metrics have been introduced to compare the warming contribution of components with different physical properties and express them in a common unit, as CO2-equivalent emissions [38]. The Global Warming Potential (GWP) measures how much energy the emission of 1 ton of a gas will absorb over a given time period, relative to the emission of 1 ton of carbon dioxide (CO2). Varying the reference time horizon, e.g., 20, 50, or 100 years, affects the magnitude of the metric according to the atmospheric lifetime decay of the forcing agent. By definition, CO2 has a GWP of 1 regardless of the time period used, since it is the gas being used as the reference. GWP is only one of various CO2-equivalent emission metrics proposed in the literature, and the choice of one metric over the other depends on the environmental target that must be investigated (e.g., emission reduction target, temperature target).
For a given trajectory, the climate forcing agents produced by the aircraft (i.e., GHGs and contrail cirri) can be combined by means of a weighted sum, resulting in the environmental cost J E N V :
J E N V = s = 1 N s p e c i e s E M s t 0 t f F F ( t ) · E I s ( t ) d t
where for each emission specie s, the weight E M s is given by the chosen emission metric from Table 2. The environmental cost is measured in CO2-equivalent mass, since it is the result of the mass of the emission given by the integration term, as for Equation (26), multiplied by the CO2-equivalent emission metric.
The value of the emission index of AIC, E I A I C = { 0 , 1 } , equals 1 in the presence of AIC, according to Equation (34), and is null otherwise.

2.5.3. Multi-Objective Cost Function

A nontrivial multi-objective optimization problem that must minimize two or more conflicting objectives does not have a single solution that minimize each objective at the same time, but a set of optimal solutions defined as the Pareto-optimal set of solutions that lay on the Pareto front of the feasible solutions set in the objective space. For each solution in the Pareto-optimal set, it is not possible to further optimize one of the objectives without penalizing the others.
The multi-objective trajectory optimization is solved using a quadratic weighted sum method, where the multi-objective cost function of the OCP problem defined in Equation (17a) can be rewritten using the following scalar expression [53]:
J = w D O C J D O C σ D O C 2 + w E N V J E N V σ E N V 2
where w D O C and w E N V are the weights attributed to each objective function, while σ D O C and σ E N V are scaling factors used to normalize the single objective costs with respect to a reference value.
Weights of the weighted sum method must satisfy the following conditions:
{ w D O C , w E N V 0 w D O C + w E N V = 1
and their values are assigned using the following linear form:
{ w D O C = ( 1 κ ) w E N V = κ , κ [ 0 , 1 ]
Shifting the value of the parameter κ towards 0 increases the weight of J D O C in the optimization, producing trajectories with lower operating costs, while shifting κ towards 1 increases the weight of J E N V , producing greener trajectories. The squared weighted sum method consists of solving the optimization problem for several values of κ to find a set of feasible trajectories and identify the Pareto set of optimal solutions.
The scaling factors σ D O C and σ E N V are respectively the DOC and the environmental cost of the minimum DOC trajectory, calculated by solving the OCP problem as defined in (17), with a single objective cost function consisting only in the Mayer term of Equation (17a) and corresponding to the definition of the DOC in Equation (39):
min X D O C , U D O C , t f D O C J D O C [ X , t f ]
where X D O C is the state vector of the DOC-optimal trajectory, U D O C the control vector, and t f D O C the final time. The scaling factors σ D O C and σ E N V are then calculated using Equations (39) and (40):
σ D O C = J D O C [ X D O C , U D O C , t f D O C ]
σ E N V = J E N V [ X D O C , t f D O C ]
The weighted sum method is widely used given its simplicity; however, it has some issues that must be taken into account. The solution for a given κ does not necessarily belong to the Pareto set in the objective space, and uniform spacing in κ may not correspond to uniformed solutions in the objective space. Furthermore, the method only works with convex Pareto fronts. These issues can be mitigated by proper normalization of the objectives and by solving the problem using a high number of κ values, to obtain a valid Pareto set of solutions. If these resolutions do not work, the problem may have a non-convex Pareto front and another multi-objective optimization method should be used.

2.6. Multi-Phase Trajectory Optimization

The optimization tool optimizes multi-phase trajectory composed by climb, cruise, and descent phases. The multi-phase optimization is performed sequentially, starting from the climb phase, and each phase is optimized individually. The following continuity boundary conditions are formulated to assure continuity of the aircraft states during phase changes at top-of-climb (TOC) and top-of-descent (TOD):
x 0 ( c r ) = x N c l
x 0 ( d s ) = x N c r
Discontinuity of the throttle control variable Π is permitted, while the flight path angle γ and the heading angle ψ are constrained as the states.
To obtain an optimal climb-cruise-descent arc trajectory, the cost function of each phase must include an estimation of the costs of the subsequent phases, following the Bellman’s principle of optimality [54]. In other words, the optimal climb trajectory is not the one with the minimum climb cost, but the one that produces the optimal climb-cruise-descent trajectory with the overall minimum cost. Prior to each phase optimization, the initial guess (IG) trajectory is automatically calculated by propagating the aircraft conditions at the initial waypoint WP of the given phase, assuming an aircraft configuration that guarantees the maximum specific range (SPR) along the geodesic trajectory that leads to WP2.
Inputs for the optimization problems are the start-of-climb (SOC) and end-of-descent (EOD) coordinates and altitude, defined respectively as WP1 and WP2, the initial mass of the aircraft m 0 , the time at SOC. TOC and TOD waypoints are obtained through the optimization process. The optimal trajectory is performed assuming a free route airspace (FRA) between WP1 and WP2, with continuous climb and descent operations (CCO and CDO), and idle descent. The optimization procedures for each phase are described in the following subsections.

2.6.1. Climb Phase

Climb Phase Cost Function

Climb is the first phase to be optimized, and the cost function of the optimal climb needs to also include the approximated costs of cruise and descent phases. The cost function to be minimized for the climb phase is defined as follows:
J = J c l [ X ( c l ) , U ( c l ) , t f ( c l ) , κ ] + J ˜ c r + J ˜ d s
where J c l is the actual cost of the climb trajectory as defined in Equation (17a)
J c l = J [ X ( c l ) , U ( c l ) , t f ( c l ) , κ ]
and is a function of the climb state and control variables, of the final time of climb t f ( c l ) , and of the multi-objective coefficient κ . The subscript ( c l ) indicates variables of the climb phase, J ˜ c r and J ˜ d s are the approximated costs of the cruise and descent phases.
The cost of the approximated descent phase J ˜ d s can be considered to be negligible when compared to cruise and climb cost, since with idle descent F F d s 0 .
The cost of the approximated cruise phase J ˜ c r can be estimated by making the following assumptions:
(i)
the lateral path of the approximated cruise lays on the geodesic curve that links the TOC to WP2 (i.e., minimum distance lateral path);
(ii)
the TOD of the approximated cruise corresponds to WP2 since for long-haul flights the descent track path distance is negligible when compared to cruise distance, hence J ˜ c r | T O C W P 2 J ˜ c r | T O C T O D . This assumption allows the avoidance of the estimation of the TOD for each iteration of the cost function;
(iii)
the vertical path and the true airspeed of the aircraft for each point in the cruise trajectory are such that the SPR is maximized. Hence:
max V , h S P R ( κ ) = d c r J c r ( κ )
where d c r is cruise ground distance.
The states of the approximate cruise trajectory x ˜ c r can be estimated by propagating the TOC states x N ( c l ) through integration. The integration is performed using the implicit midpoint method by discretizing the geodesic arc between TOC and WP2 with a spatial resolution of 10 kilometers, i.e., higher than the spatial resolution of the GFS grid.
The climb cost function in Equation (49) can be rewritten as:
J = J c l [ X ( c l ) , U ( c l ) , t f ( c l ) , κ ] + J ˜ c r [ x N ( c l ) , t f ( c l ) , κ ]
and it is evident that the approximation of cruise cost is a Mayer term of the OCP since it depends only on the final states and time of climb at TOC.

Climb Phase Initial Guess

Initial guess values of the climb phase [ X ( c l ) , U ( c l ) , t ( c l ) ] I G to be used as an input for the optimization algorithm are calculated by integrating the states at WP1 along a climb trajectory with lateral path on the geodesic arc between WP1 and WP2, assuming that at each step the true airspeed V is such that the rate of climb (ROC) is maximized:
max V R O C = T m a x D m V ˙ V g
The propagation of the IG climb trajectory stops when the altitude reaches a value for which the SPR of an approximated cruise towards WP2, as defined in Equation (51), is maximized.

2.6.2. Cruise Phase

Cruise Phase Cost Function

In cruise optimization, the cost function to be minimized is equivalent to the actual cost of the cruise trajectory:
J = J c r [ X ( c r ) , U ( c r ) , t f ( c r ) , κ ]
Cruise trajectory is delimited by the TOC and TOD waypoints. The states and controls at TOC are bounded to the ones obtained in climb optimization:
x N ( c l ) = x 1 ( c r )
u N ( c l ) = u 1 ( c r )
The TOD is defined by the cost function in Equation (54) in conjunction with the following constraint:
d d s ( [ x N ( c r ) ] ) d t o W P 2 ( [ x N ( c r ) ] ) 0
where d d s is the estimated range of the descent phase and d t o W P 2 is the ground distance from the TOD to WP2. The optimal descent is the one that maximizes its range, since it reduces the cruise distance. The constraint defined in Equation (57) is satisfied for every TOD whose horizontal distance from WP2 d t o W P 2 is lower or equal to d d s , but the minimization of the cost function pushes the TOD further from WP2, on the boundary of the area defined by the constraint.
The optimal descent range d d s is estimated by maximizing the lift-to-drag ratio of the aircraft in descent, taking into account the idle thrust and the acceleration. The states of the aircraft at TOD are propagated through implicit midpoint integration by discretizing the altitude between TOD and WP2. The lateral trajectory is constrained to lay on the geodesic between TOD and WP2, while the true airspeed V for each step is obtained with
min V D T i d l e m V ˙ L

Cruise Phase Initial Guess

Initial guess values for the cruise phase [ X ( c r ) , U ( c r ) , t ( c r ) ] I G are calculated by integrating the states at TOC resulting from climb optimization, along a cruise trajectory that lays on the geodesic between TOC and WP2, assuming that at each step the maximum specific range SPR as defined in Equation (51) is maximized. The initial guess for the TOD coordinates is found when d d s ( [ x N ( c r ) ] I G ) = d t o W P 2 ( [ x N ( c r ) ] ) .

2.6.3. Descent Phase

The initial and final waypoints of the descent phase are constrained to the TOD calculated in cruise optimization and to WP2. During all the idle descent, the fuel flow is at its minimum and can be considered constant, so the operating and environmental costs as defined in Equations (39) and (40) can be minimized only by reducing the descent time. TOD in cruise optimization is calculated for a maximum range idle descent, consequently the descent duration cannot be reduced without negatively impacting the descent range.
The states, controls, and final time of the descent phase are obtained by integrating the cruise states at TOD, along a cruise trajectory that lays on the geodesic directed towards WP2, assuming that at each step the true airspeed satisfies Equation (58).

3. Results

The results of the multi-objective and multi-phase trajectory optimization tool are presented for a test case involving a transatlantic commercial flight from Rome, Italy to New York, USA on a Boeing 747-100. The inputs of the test case optimization are reported in Table 3. The number of nodes of the Chebyshev-Gauss-Lobatto (CGL) grid is user-defined. In this case, 20 nodes have been used for the cruise phase, while 10 nodes have been used for the climb and descent phases.
The output of the optimization tool is the Pareto set of optimal climb-cruise-descent trajectories from WP1 to WP2, where the objectives of the multi-objective optimization are the Direct Operating Cost (DOC), J D O C , and the environmental cost, J E N V , of the trajectory. For each trajectory on the Pareto front, it is not possible to further improve one of the two objective costs without degrading the other. The Pareto set of optimal solutions includes the DOC-optimal and the environmental-optimal trajectories that represent respectively the trajectories with the higher environmental and operating costs among the Pareto-optimal solutions. Each trajectory provided in output is represented as a time series of the states and controls of the aircraft to be processed by the Automatic Flight Control System (AFCS) of the aircraft for 4D trajectory tracking. In particular, the top-of-climb (TOC) and top-of-descent (TOD) waypoints are an output of the optimization, as well as the time of arrival at WP2.
Three multi-objective optimizations are performed considering an environmental cost function based on the Global Warming Potential (GWP) emission metric integrated for a 20-, 50-, and 100-year time horizons.
The plots in Figure 2 report the results of each optimization problem in the objective space. The axes of the objective spaces are expressed in terms of percentage variation with respect to the DOC-optimal trajectory, namely Δ J D O C for the DOC and Δ J E N V for the environmental cost. It must be noted that the DOC-optimal trajectory is the same for each of the three problems, since it is obtained optimizing with a cost function where the environmental term has a null weight. Figure 2 depicts the set of Pareto-optimal trajectories (in green) as well as the other sub-optimal feasible trajectories (in blue) produced by the optimization tool. Results show that by increasing the operating costs by 6.3% with respect to the DOC-optimal trajectory it is possible to reduce environmental costs by 51.6% for GWP20, by 47.1% for GWP50, and by 38.1% for the GWP100 case. Nonetheless, in the flight planning phase the aircraft operator can choose the trajectory that better reflects its economic and environmental goals among the ones belonging to the Pareto set of optimal solutions, and submit it to the Air Traffic Controller (ATC).
Figure 3 illustrates a breakdown of the environmental cost of the Pareto-optimal set of solutions. The environmental cost is expressed in terms of CO2-equivalent mass of emissions to underline the contribution of each climate forcing produced by the aircraft in the climb-cruise-descent arc. The first aspect that can be noticed is the reduction in the absolute value of the CO2-equivalent emissions for increasing values of the GWP time horizon. This is because short-lived climate forcings such as aircraft-induced clouds (AIC) and NOX have a smaller GWP in the long term when compared to the one of long-lived species such as CO2. The other result apparent from Figure 3 is that the main strategy adopted by the optimization tool to mitigate the environmental cost of the trajectory is to reduce the formation of AIC along the path. In particular, it can be noticed that the environmental-optimal trajectory presents zero contrail formation.
Results from the GWP100 case are now presented, to analyze the behavior of the optimization tool when shifting from DOC-optimal to environmental-optimal trajectory.
Figure 4a shows the lateral path of some of the solutions belonging to the Pareto set of optimal trajectories, including the DOC-optimal (purple) and the environmental-optimal path (green). The DOC-optimal solution ignores the presence of ice supersaturated regions (ISSRs) along the track, producing AIC on its trail, pursuing a flight track with low route extension (RE) and intercepting favorable tailwind conditions, especially in the north-western Atlantic region. On the other hand, the environment optimal trajectory adopts a greater route extension, as it can also be seen in Table 4, to avoid ISSRs and avoid AIC formation. Intermediate optimal solutions show a mixed behavior between the two limiting cases.
The altitude profile of the trajectories represented in Figure 4b shows how the DOC-optimal trajectory is characterized by a climbing cruise with a smoothly increasing flight path angle, pursuing a cruise altitude that maximizes the specific range (SPR) with respect to the DOC, regardless of AIC formation along the path. Increasing the weight of the environmental cost forces the optimization tool to reduce the altitude in the northern Atlantic area to avoid the ISSR visible in Figure 4a. This ISSR region is completely avoided by the environment optimal trajectory that deviates to the north along a longer route. The environment optimal trajectory also has a lower TOC altitude. This can be explained by the fact that lowering cruise altitude reduces NOX emission, at a cost of increased CO2 emissions. As it can be noticed from Figure 3, once AIC contribution to the CO2-equivalent emissions is eliminated, a further reduction of the net CO2-equivalent emissions can be obtained by a reduction of the net NOX component, despise the increase of CO2 emissions. The true airspeed (TAS) profile depicted in Figure 4c shows the tendency of a reduction of the TAS for increasing cruise altitudes, and fluctuations due to the variation of atmospheric conditions in the Global Forecast System (GFS) weather grid.
Figure 5 and Table 4 gather other trajectory data of the Pareto-optimal solutions for the GWP100 case.
The route extension (RE) with respect to the great circle arc between WP1 and WP2 (6900 km) of the Pareto-optimal trajectories, represented in Figure 5a, shows a trend of increasing RE when moving towards environment optimal trajectories on the Pareto front. This is due to the fact that the optimization tool adopts longer routes to avoid ISSR and reduce AIC formation. The large scattering of the data is because the optimal lateral path is strictly dependent on lateral winds, so the DOC component of the cost function is not necessarily lower for small values of the RE.
Figure 5b,c show flight duration and fuel consumption of the optimal trajectories. Both tend to increase when leaning towards less cost efficient and greener trajectories, because the DOC cost function defined in Equation (39) is proportional to both flight time and fuel consumption, with scattering on fuel mass being lower due to the higher weight of fuel consumption on DOC with respect to flight operating cost per hour.
Figure 5d represents AIC length as a fraction of route length, and shows a complete elimination of consistent reduction in AIC formation when shifting from the DOC-optimal to the environmental-optimal solution. In particular, the length of along-track AIC is 1607 km for the DOC-optimal trajectory and is completely eliminated in the environmental-optimal trajectory.

4. Conclusions

This paper presents the development of a 4D trajectory optimization tool capable of generating a Pareto set of optimal trajectories in terms of Direct Operating Cost (DOC) and climate cost, based on the Global Warming Potential (GWP) emission metric, of the climb-cruise-descent arc of a commercial flight. The multi-phase and multi-objective optimization algorithm is based on a Chebyshev pseudospectral direct collocation method. The tool has been designed to be integrated within the Flight Management System (FMS) of a commercial aircraft, the Boeing 747-100, to perform 4D trajectory flight planning in a Free Route Airspace (FRA) context. The output of the trajectory optimization tool is a set of Pareto-optimal climb-cruise-descent trajectories ranging from the DOC-optimal to the climate-optimal trajectory, and it allows the aircraft operator to select the best compromise flight trajectory in the flight planning phase that satisfies both economic and environmental constraints.
The trajectory optimization tool has been employed on a test case flight from Rome to New York. Three separate multi-objective optimizations have been performed, using environmental cost functions based on the GWP integrated for a time horizon of 20, 50, and 100 years. For the GWP100 case, the environmental-optimal trajectory grants a 38.1% reduction of the climate cost at the expense of a 6.3% increase of the DOC with respect to the DOC-optimal trajectory. In absolute terms, the environmental-optimal trajectory of the GWP100 case reduces by 242.9 Ton the CO2-equivalent emissions of the aircraft climate forcings at a cost of 4.8 k$ of DOC.
The environmental cost function takes into account the climate forcings of CO2 and non-CO2 emissions, as well as aircraft-induced clouds that are the main climate pollutant produced by aviation. In particular, the optimization algorithm showed good aircraft-induced clouds (AIC) formation avoidance capabilities, by adjusting the lateral and vertical trajectory path to avoid ice supersaturated regions (ISSR). In fact, the environmental-optimal trajectory for the test case based on GWP100 managed to completely eliminate on-track AIC formation.
The accuracy of the trajectory optimization tool is affected by the reliability of weather data and forecasts, as well as by the uncertainties of the climate metrics that have been used for the estimation of the environmental cost index. However, the same tool developed in this paper can be used with more accurate weather data and climate metrics data as soon as they will be available to the scientific community. Future work will include the implementation of dynamic weather information and will take into account flight planning with adverse weather phenomena along the path, as well as the integration of the other flight phases occurring in terminal maneuvering areas in the proximity of airports that have a lower margin for climate impact mitigation but can be addressed for the reduction of noise pollution and the improvement of local air quality.

Author Contributions

Conceptualization, A.V. and M.B.; methodology, A.V. and M.B.; software, A.V.; validation, A.V.; formal analysis, A.V.; investigation, A.V.; writing—original draft preparation, A.V.; writing—review and editing, A.V., M.B. and A.L.; visualization, A.V.; supervision, M.B.; project administration, M.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

AICAircraft-induced cloudsGWPGlobal Warming Potential
ATCAir Traffic ControlISSRIce supersaturated regions
ATMAir Traffic ManagementNLPNonlinear Programming problem
BADABase of aircraft dataOCPOptimal Control Problem
BM2Boeing Fuel Flow Method 2RERoute extension
CGLChebyshev-Gauss-LobattoREIRelative Emission In
DOCDirect Operating CostRFRadiative Forcing
DOFDegrees of freedomRHRelative Humi
EIEmission IndexROCRate of Climb
ENVEnvironmental CostSACSchmidt-Appleman Criterion
ERFEquivalent Radiative ForcingSPRSpecific Range
FFFuel FlTOCTop of Climb
FMSFlight Management SystemTODTop of Descent
GFSGlobal Forecast SystemTSFCThrust Specific Fuel Consumpti
GHGGreenhouse GasWPWaypoint

References

  1. Lee, D.S.; Fahey, D.W.; Forster, P.M.; Newton, P.J.; Wit, R.C.; Lim, L.L.; Owen, B.; Sausen, R. Aviation and Global Climate Change in the 21st Century. Atmos. Environ. 2009, 43, 3520–3537. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Eurocontrol, Free Route Airspace. Available online: https://www.eurocontrol.int/concept/free-route-airspace (accessed on 12 October 2021).
  3. AIP Italia, ENAV, FRAIT—Free Route Italy. Available online: https://ec.europa.eu/transport/sites/default/files/AIC_A_2016_11.pdf (accessed on 12 October 2021).
  4. EASA. European Aviation Environmental Report. 2019. Available online: https://www.easa.europa.eu/eaer/ (accessed on 12 October 2021).
  5. Eurocontrol, Free Routes Airspace (FRA) Design Guidelines. Available online: https://www.eurocontrol.int/sites/default/files/2020-11/fra-design-guidelines-1.0.pdf (accessed on 12 October 2021).
  6. Kistan, T.; Gardi, A.; Sabatini, R.; Ramasamy, S.; Batuwangala, E. An evolutionary outlook of air traffic flow management techniques. Prog. Aerosp. Sci. 2017, 88, 15–42. [Google Scholar] [CrossRef]
  7. Betts, J.T. Survey of Numerical Methods for Trajectory Optimization. J. Guid. Control Dyn. 1998, 21, 193–207. [Google Scholar] [CrossRef]
  8. Rao, A.V. Survey of Numerical Methods for Optimal Control. Adv. Astronaut. Sci. 2010, 135, 497–528. [Google Scholar]
  9. Hagelauer, P.; Mora-Camino, F. Flight Management Systems and Aircraft 4D Trajectory Optimization. IFAC Proc. Vol. 1997, 30, 351–356. [Google Scholar] [CrossRef]
  10. Bryson, A.; Ho, Y.C. Numerical solution of optimal programming and control problems. In Applied Optimal Control: Optimization, Estimation, and Control, 1st ed.; Taylor & Francis: Abingdon, UK, 1975; pp. 212–243. [Google Scholar]
  11. Fornberg, B.; Sloan, D.M. A review of pseudospectral methods for solving partial differential equations. Acta Numer. 1994, 3, 203–267. [Google Scholar] [CrossRef]
  12. The Optimal Control Problem. Practical Methods for Optimal Control and Estimation Using Nonlinear Programming, 2nd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2010; pp. 123–217. [Google Scholar]
  13. Ross, I.M.; Fahroo, F. Pseudospectral Knotting Methods for Solving Nonsmooth Optimal Control Problems. J. Guid. Control. Dyn. 2004, 27, 397–405. [Google Scholar] [CrossRef]
  14. Betts, J.T.; Cramer, E.J. Application of Direct Transcription to Commercial Aircraft Trajectory Optimization. J. Guid. Control Dyn. 1995, 18, 151–159. [Google Scholar] [CrossRef]
  15. Soler, M.; Olivares, A.; Staffetti, E. Multiphase Optimal Control Framework for Commercial Aircraft Four-Dimensional Flight-Planning Problems. J. Aircr. 2015, 52, 274–286. [Google Scholar] [CrossRef]
  16. Rao, A.V.; Benson, D.A.; Darby, C.; Patterson, M.A.; Francolin, C.; Sanders, I.; Huntington, G.T. Algorithm 902: GPOPS, A MATLAB Software for Solving Multiple-Phase Optimal Control Problems Using the Gauss Pseudospectral Method. ACM Trans. Math. Softw. 2010, 37, 1–39. [Google Scholar] [CrossRef]
  17. Ross, I.M. Enhancements to the DIDO Optimal Control Toolbox. arXiv 2020, arXiv:math.OC/2004.13112. Available online: https://arxiv.org/pdf/2004.13112.pdf (accessed on 12 October 2021).
  18. Rutquist, E.; Edvall, M.M. PROPT-Matlab Optimal Control Software. Available online: https://www.tomopt.com/docs/TOMLAB_PROPT.pdf (accessed on 12 October 2021).
  19. Cassaro, M.; Gunetti, P.; Battipede, M.; Gili, P. Overview of the Multipurpose Aircraft Simulation Laboratory experience. In Proceedings of the 2013 Aviation Technology, Integration, and Operations Conference, Los Angeles, CA, USA, 12–14 August 2013. [Google Scholar]
  20. Sirigu, G.; Battipede, M.; Gili, P.; Cassaro, M. FMS and AFCS Interface for 4D Trajectory Operations; SAE Technical Paper-2015; SAE International: Warrendale, PA, USA, 2015. [Google Scholar]
  21. Mazzotta, D.G.; Sirigu, G.; Cassaro, M.; Battipede, M.; Gili, P. 4D Trajectory Optimization Satisfying Waypoint and No-Fly Zone Constraints. WSEAS Trans. Syst. Control 2017, 12, 221–231. [Google Scholar]
  22. Nuic, A.; Poles, D.; Moulliet, V. BADA: An advanced aircraft performance model for present and future ATM systems. Int. J. Adapt. Control Signal Process. 2010, 24, 850–866. [Google Scholar] [CrossRef]
  23. Hanke, C.; Nordwall, D. The simulation of a jumbo jet Transport aircraft Volume II: Modeling data. In D6-30643; The Boeing Company: Wichita, KS, USA, 1970. [Google Scholar]
  24. Fahroo, F.; Ross, I.M. Direct Trajectory Optimization by a Chebyshev Pseudospectral Method. J. Guid. Control Dyn. 2002, 25, 160–166. [Google Scholar] [CrossRef] [Green Version]
  25. Lufthansa Systems, Lido Flight 4D. Available online: https://www.lhsystems.com/solutions-services/operations-solutions/lidoflightplanning/lidoflight4D (accessed on 12 October 2021).
  26. Montzka, S.A.; Dlugokencky, E.J.; Butler, J.H. Non-CO2 greenhouse gases and climate change. Nature 2011, 476, 43–50. [Google Scholar] [CrossRef]
  27. Lee, D.; Fahey, D.; Skowron, A.; Allen, M.; Burkhardt, U.; Chen, Q.; Doherty, S.; Freeman, S.; Forster, P.; Fuglestvedt, J.; et al. The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018. Atmos. Environ. 2021, 244, 117834. [Google Scholar] [CrossRef] [PubMed]
  28. ICAO. Assessing Current Scientific Knowledge, Uncertainties and Gaps in Quantifying Climate Change, Noise and Air Quality Aviation Impacts. Available online: https://www.icao.int/environmental-protection/Documents/CaepImpactReport.pdf (accessed on 12 October 2021).
  29. Park, S.G.; Clarke, J.P. Vertical Trajectory Optimization to Minimize Environmental Impact in the Presence of Wind. J. Aircr. 2016, 53, 725–737. [Google Scholar] [CrossRef] [Green Version]
  30. Sirigu, G.; Clarke, J.P.; Battipede, M.; Gili, P. Hybrid Particle Swarm Optimization with Parameter Fixing: Application to Automatic Taxi Management. J. Air Transp. 2020, 28, 36–48. [Google Scholar] [CrossRef]
  31. Kärcher, B. Formation and radiative forcing of contrail cirrus. Nat. Commun. 2018, 9, 117834. [Google Scholar] [CrossRef] [PubMed]
  32. Sridhar, B.; Ng, H.K.; Chen, N.Y. Aircraft Trajectory Optimization and Contrails Avoidance in the Presence of Winds. J. Guid. Control Dyn. 2011, 34, 1577–1584. [Google Scholar] [CrossRef]
  33. Soler, M.; Zou, B.; Hansen, M. Contrail Sensitive 4D Trajectory Planning with Flight Level Allocation Using Multiphase Mixed-Integer Optimal Control. In Proceedings of the AIAA Guidance, Navigation, and Control (GNC) Conference, Boston, MA, USA, 19–22 August 2013. [Google Scholar]
  34. Irvine, E.A.; Hoskins, B.J.; Shine, K.P. A simple framework for assessing the trade-off between the climate impact of aviation carbon dioxide emissions and contrails for a single flight. Environ. Res. Lett. 2014, 9, 064021. [Google Scholar] [CrossRef]
  35. Rosenow, J.; Fricke, H. Individual Condensation Trails in Aircraft Trajectory Optimization. Sustainability 2019, 11, 6082. [Google Scholar] [CrossRef] [Green Version]
  36. Teoh, R.; Schumann, U.; Majumdar, A.; Stettler, M.E.J. Mitigating the Climate Forcing of Aircraft Contrails by Small-Scale Diversions and Technology Adoption. Environ. Sci. Technol. 2020, 54, 2941–2950. [Google Scholar] [CrossRef] [PubMed]
  37. Matthes, S.; Lührs, B.; Dahlmann, K.; Grewe, V.; Linke, F.; Yin, F.; Klingaman, E.; Shine, K.P. Climate-Optimized Trajectories and Robust Mitigation Potential: Flying ATM4E. Aerospace 2020, 7, 156. [Google Scholar] [CrossRef]
  38. Myhre, G.; Shindell, D.; Pongratz, J. Intergovernmental Panel on Climate Change, Anthropogenic and Natural Radiative Forcing. In Climate Change 2013–The Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2014; pp. 659–740. [Google Scholar]
  39. DuBois, D.; Paynter, G.C. Fuel Flow Method2 for Estimating Aircraft Emissions. In Non-Conference Specific Technical Papers-2006; SAE International: Warrendale, PA, USA, 2006. [Google Scholar]
  40. European Union Aviation Safety Agency. ICAO Aircraft Engine Emissions Databank. Available online: https://www.easa.europa.eu/domains/environment/icao-aircraft-engine-emissions-databank (accessed on 12 October 2021).
  41. Schumann, U. A contrail cirrus prediction model. Geosci. Model Dev. 2012, 5, 543–580. [Google Scholar] [CrossRef] [Green Version]
  42. Kim, I.; de Weck, O. Adaptive weighted-sum method for bi-objective optimization: Pareto front generation. Struct. Multidiscipl. Optim. 2005, 29, 149–158. [Google Scholar] [CrossRef]
  43. Trefethen, L.N. Equispaced Points, Runge Phenomenon. In Approximation Theory and Approximation Practice, Extended Edition; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2019; pp. 95–102. [Google Scholar]
  44. Trefethen, L.N. Chebyshev Differentiation Matrices. In Spectral Methods in MATLAB; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2000; pp. 51–59. [Google Scholar]
  45. Waldvogel, J. Fast Construction of the Fejér and Clenshaw–Curtis Quadrature Rules. BIT Numer. Math. 2006, 46, 195–202. [Google Scholar] [CrossRef] [Green Version]
  46. National Weather Service: Pressure Altitude. Available online: https://www.weather.gov/media/epz/wxcalc/pressureAltitude.pdf (accessed on 12 October 2021).
  47. Picard, A.; Davis, R.S.; Gläser, M.; Fujii, K. Revised formula for the density of moist air (CIPM-2007). Metrologia 2008, 45, 149–155. [Google Scholar] [CrossRef]
  48. World Meteorological Organization (WMO). Cloud Atlas. Available online: https://cloudatlas.wmo.int/aircraft-condensation-trails.html (accessed on 12 October 2021).
  49. Sonntag, D. Advancements in the field of hygrometry. Meteorol. Z. 1994, 3, 51–66. [Google Scholar] [CrossRef]
  50. Federal Aviation Administration, Benefit-Cost Analysis: Aircraft Operating Costs. Available online: https://www.faa.gov/regulations_policies/policy_guidance/benefit_cost/media/econ-value-section-4-op-costs.pdf (accessed on 12 October 2021).
  51. Federal Aviation Administration, Benefit-Cost Analysis: Economic Values Related to Aircraft Performance. Available online: https://www.faa.gov/regulations_policies/policy_guidance/benefit_cost/media/econ-value-section-6-perf-factors.pdf (accessed on 12 October 2021).
  52. Schumann, U.; Penner, J.E.; Chen, Y.; Zhou, C.; Graf, K. Dehydration effects from contrails in a coupled contrail—Climate model. Atmos. Chem. Phys. 2015, 15, 11179–11199. [Google Scholar] [CrossRef] [Green Version]
  53. Yang, X.S. Chapter 14—Multi-Objective Optimization. In Nature-Inspired Optimization Algorithms; Yang, X.S., Ed.; Elsevier: Oxford, UK, 2014; pp. 197–211. [Google Scholar]
  54. Bellman, R.E. The Structure of Dynamic Programming Process. In Dynamic Programming; Princeton University Press: Princeton, NJ, USA, 1957; pp. 81–115. [Google Scholar]
Figure 1. Scheme of the averaging method used to pass along-trajectory weather information to the CGL nodes p of the optimization algorithms. The number of sub-CGL nodes s is chosen to have a spatial resolution similar to the GFS grid with weather data. The gradient color map represents different values of a generic GFS data variable.
Figure 1. Scheme of the averaging method used to pass along-trajectory weather information to the CGL nodes p of the optimization algorithms. The number of sub-CGL nodes s is chosen to have a spatial resolution similar to the GFS grid with weather data. The gradient color map represents different values of a generic GFS data variable.
Aerospace 08 00395 g001
Figure 2. Pareto-optimal (green) and dominated sub-optimal solutions (blue) of the three multi-objective optimization problems, with environmental cost function based on: (a) GWP20, (b) GWP50, (c) GWP100.
Figure 2. Pareto-optimal (green) and dominated sub-optimal solutions (blue) of the three multi-objective optimization problems, with environmental cost function based on: (a) GWP20, (b) GWP50, (c) GWP100.
Aerospace 08 00395 g002
Figure 3. Individual climate forcing components of the environmental cost J E N V , expressed in CO2-equivalent Ton, for environmental cost function based on: (a) GWP20, (b) GWP50, (c) GWP100. Negative values for SO2 indicate a cooling forcing.
Figure 3. Individual climate forcing components of the environmental cost J E N V , expressed in CO2-equivalent Ton, for environmental cost function based on: (a) GWP20, (b) GWP50, (c) GWP100. Negative values for SO2 indicate a cooling forcing.
Aerospace 08 00395 g003
Figure 4. Some of the Pareto-optimal trajectories for the GWP100 case, ranging from the optimal DOC trajectory (purple) to the climate-optimal trajectory (light green): (a) lateral path plotted on a map, with cyan segments indicating along track AIC formation, blue areas representing ISSR, and wind speeds expressed by a contour and vector map; (b) altitude profile plotted against time with along track AIC formation; (c) TAS profile. The atmospheric conditions represented in (a) are evaluated at an indicative cruise level, given by the average altitude of each plotted trajectory as a function of the longitude.
Figure 4. Some of the Pareto-optimal trajectories for the GWP100 case, ranging from the optimal DOC trajectory (purple) to the climate-optimal trajectory (light green): (a) lateral path plotted on a map, with cyan segments indicating along track AIC formation, blue areas representing ISSR, and wind speeds expressed by a contour and vector map; (b) altitude profile plotted against time with along track AIC formation; (c) TAS profile. The atmospheric conditions represented in (a) are evaluated at an indicative cruise level, given by the average altitude of each plotted trajectory as a function of the longitude.
Aerospace 08 00395 g004
Figure 5. The plots show trajectory data of the Pareto-optimal solutions for the environmental cost function based on GWP100: (a) route extension with respect to the great circle distance between WP1 and WP2; (b) flight duration; (c) fuel consumption; (d) fraction of AIC length with respect to route length, plotted against Δ J D O C . Green dots represent trajectories plotted in Figure 4.
Figure 5. The plots show trajectory data of the Pareto-optimal solutions for the environmental cost function based on GWP100: (a) route extension with respect to the great circle distance between WP1 and WP2; (b) flight duration; (c) fuel consumption; (d) fraction of AIC length with respect to route length, plotted against Δ J D O C . Green dots represent trajectories plotted in Figure 4.
Aerospace 08 00395 g005
Table 1. Values of the constant emission indexes [27].
Table 1. Values of the constant emission indexes [27].
EmissionEmission Index (EI)
CO23.159 kg/kg fuel
H2O1.231 kg/kg fuel
SO21.2 g/kg fuel
Soot0.03 g/kg fuel
Table 2. Global Warming Potentials of the climate forcings considered in the environmental cost function for time horizons of 20, 50, 100 years [27]. AIC values have been adjusted considering a global mean flight fraction with contrail formation of 15.6% [52], while values reported in [27] consider AIC based on global flight distance of flights scheduled in 2018.
Table 2. Global Warming Potentials of the climate forcings considered in the environmental cost function for time horizons of 20, 50, 100 years [27]. AIC values have been adjusted considering a global mean flight fraction with contrail formation of 15.6% [52], while values reported in [27] consider AIC based on global flight distance of flights scheduled in 2018.
EmissionGWP20GWP50GWP100
CO2111
AIC (Tg CO2 basis)14.876.994.04
AIC (km basis)25612271
Net NOX619205114
Soot428820181166
SO2−832−392−226
Water vapor0.220.100.06
Table 3. Inputs of the optimization tool for the test case flight from Rome to New York. WP1 and WP2 are respectively the start-of-climb and end-of-descent waypoints. Aircraft mass and time at WP2 are an output of the optimization.
Table 3. Inputs of the optimization tool for the test case flight from Rome to New York. WP1 and WP2 are respectively the start-of-climb and end-of-descent waypoints. Aircraft mass and time at WP2 are an output of the optimization.
LatLonAltitudeAircraft MassDateTime (UTC)
WP141.902812.49641000 m340 Ton25 July 202100:00
WP240.7306−73.93521000 m---
Table 4. The table reports trajectory data for the Pareto-optimal solutions for the environmental cost function based on GWP100. The trajectories are numbered for increasing direct operating cost J D O C and decreasing environmental cost J E N V , i.e., moving from left to right along the Pareto front in Figure 2c. The RE is calculated with respect to the great circle arc between WP1 and WP2 (6900 km), while AIC percentage length is expressed as a fraction of the track path distance. Boldface rows indicate the trajectories plotted in Figure 4.
Table 4. The table reports trajectory data for the Pareto-optimal solutions for the environmental cost function based on GWP100. The trajectories are numbered for increasing direct operating cost J D O C and decreasing environmental cost J E N V , i.e., moving from left to right along the Pareto front in Figure 2c. The RE is calculated with respect to the great circle arc between WP1 and WP2 (6900 km), while AIC percentage length is expressed as a fraction of the track path distance. Boldface rows indicate the trajectories plotted in Figure 4.
Op. CostsEnv. Costs (GWP100)Route FuelAIC
JDOC Δ JDOCJENV Δ JENVExtensionTimeMassLength
ID[k$][%][CO2-eq Ton][%][km][%][min][Ton][km][%]
177.11-637.4-218+3.16%54985.6160722.6%
277.18+0.10%589.9−7.4%91+1.32%52485.6122917.6%
377.19+0.10%556.6−12.7%150+2.18%47685.7100514.3%
477.20+0.11%555.4−12.9%85+1.23%51785.796613.8%
577.20+0.12%551.8−13.4%86+1.24%53285.696913.9%
677.25+0.18%535.2−16.0%97+1.41%51585.885012.1%
778.26+1.50%526.6−17.4%91+1.32%53486.775310.8%
878.41+1.68%486.3−23.7%187+2.71%50186.55758.1%
978.89+2.31%471.1−26.1%106+1.53%51886.84075.8%
1079.20+2.71%469.1−26.4%117+1.70%49587.33675.2%
1179.29+2.82%467.5−26.7%255+3.69%50487.34115.8%
1279.90+3.62%461.9−27.5%77+1.11%51887.93635.2%
1380.04+3.80%438.0−31.3%75+1.09%52987.82113.0%
1480.12+3.90%437.1−31.4%372+5.39%50787.92203.0%
1580.15+3.94%426.8−33.0%197+2.86%52688.01552.2%
1680.22+4.03%426.7−33.1%170+2.47%50988.11231.7%
1780.63+4.57%398.9−37.4%222+3.22%51788.300.0%
1881.25+5.37%398.7−37.4%427+6.19%53988.900.0%
1981.93+6.25%394.5−38.1%357+5.17%54389.300.0%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Vitali, A.; Battipede, M.; Lerro, A. Multi-Objective and Multi-Phase 4D Trajectory Optimization for Climate Mitigation-Oriented Flight Planning. Aerospace 2021, 8, 395. https://doi.org/10.3390/aerospace8120395

AMA Style

Vitali A, Battipede M, Lerro A. Multi-Objective and Multi-Phase 4D Trajectory Optimization for Climate Mitigation-Oriented Flight Planning. Aerospace. 2021; 8(12):395. https://doi.org/10.3390/aerospace8120395

Chicago/Turabian Style

Vitali, Alessio, Manuela Battipede, and Angelo Lerro. 2021. "Multi-Objective and Multi-Phase 4D Trajectory Optimization for Climate Mitigation-Oriented Flight Planning" Aerospace 8, no. 12: 395. https://doi.org/10.3390/aerospace8120395

APA Style

Vitali, A., Battipede, M., & Lerro, A. (2021). Multi-Objective and Multi-Phase 4D Trajectory Optimization for Climate Mitigation-Oriented Flight Planning. Aerospace, 8(12), 395. https://doi.org/10.3390/aerospace8120395

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop