Robust 4D climate-optimal flight planning in structured airspace using parallelized simulation on GPUs: ROOST V1.0

The climate impact of non-CO₂ emissions, which are responsible for two-thirds of aviation radiative forcing, highly depends on the atmospheric chemistry and weather conditions. Hence, by planning aircraft trajectories to reroute areas where the non-CO₂ climate impacts are strongly enhanced, called climate-sensitive regions, there is a potential to reduce aviation-induced non-CO₂ climate effects. Weather forecast is inevitably uncertain, which can lead to unreliable determination of climate-sensitive regions and aircraft dynamical behavior and, consequently, inefficient trajectories. In this study, we propose robust climate-optimal aircraft trajectory planning within the currently structured airspace considering uncertainties in standard weather forecasts. The ensemble prediction system is employed to characterize uncertainty in the weather forecast, and climate-sensitive regions are quantified using the prototype algorithmic climate change functions. As the optimization problem is constrained by the structure of airspace, it is associated with hybrid decision spaces. To account for discrete and continuous decision variables in an integrated and more efficient manner, the optimization is conducted on the space of probability distributions defined over flight plans instead of directly searching for the optimal profile. A heuristic algorithm based on the augmented random search is employed and implemented on graphics processing units to solve the proposed stochastic optimization computationally fast. An open-source Python library called ROOST (V1.0) is developed based on the aircraft trajectory optimization technique. The effectiveness of our proposed strategy to plan robust climate-optimal trajectories within the structured airspace is analyzed through two scenarios: a scenario with a large contrail climate impact and a scenario with no formation of persistent contrails. It is shown that, for a nighttime flight from Frankfurt to Kyiv, a 55 % reduction in climate impact can be achieved at the expense of a 4 % increase in the operating cost.

1 Introduction

The aviation industry has experienced strong growth in recent years (Lee et al., 2021). Air traffic is estimated to grow at a 4.3 % annual rate over the next 20 years (Scherer, 2019). Aviation's contribution to global warming through CO₂ and non-CO₂ emissions is currently responsible for 3.5 % of total anthropogenic radiative forcing (RF) (Lee et al., 2021). The non-CO₂ climate impact includes changes in ozone and methane concentrations induced by nitrogen oxides (NO_x), water vapor (H₂O), hydrocarbons (HC), carbon monoxide (CO), sulfur oxides (SO_x), and increased cloudiness due to persistent contrail formation. Accounting for the growth rate of air traffic, a critical increase in its associated climate impact is foreseen.

Mitigating the climate impact associated with aviation-induced CO₂ emissions requires progressing along with technical aspects, such as moving toward more efficient aircraft and using novel propulsion as well as sustainable aviation fuel. Technological improvements can only be moderately incorporated into the existing aircraft fleet. This is due to the aircraft's long life service and long phases of development, production, and certification. On the other hand, the climate impact caused by non-CO₂ emissions, which is about twice as large as that from CO₂ emissions (Lee et al., 2021), reveals a high spatial and temporal dependency (Grewe et al., 2014b). Such dependencies provide the possibility to mitigate their climate effects through operational strategies, particularly aircraft trajectory optimization to avoid areas sensitive to aircraft emissions, called climate-sensitive regions (e.g., Simorgh et al., 2022).

Numerous studies have been proposed to reduce the climate impacts of non-CO₂ emissions by changing aircraft maneuvers to avoid climate-sensitive regions. These studies differ mainly in (1) how the climate-sensitive areas are defined and (2) how climate-friendly trajectories are determined. The first attempts to consider climate hotspots were based on areas sensitive to the formation of persistent contrails (see Gierens et al., 2008). In order to provide information on the spatial and temporal dependency of non-CO₂ effects, climate change functions (CCFs) were developed. These CCFs provide the climate impact of aviation emissions per flown kilometer and per emitted mass of the species as five-dimensional datasets (i.e., longitude, latitude, altitude, time, type of emission) (Matthes et al., 2012; Frömming et al., 2013; Grewe et al., 2014b). Due to their computational complexity, CCFs were not suited for real-time operations. Therefore, the so-called algorithmic climate change functions (aCCFs) were developed. The aCCFs provide a very fast computation of the individual non-CO₂ climate impact, as they are based on mathematical formulas, which only need relevant local meteorological input parameters (e.g., van Manen and Grewe, 2019). The aCCFs are well-suited for trajectory optimization tools due to their computational efficiency (Matthes et al., 2017). An enhanced and consistent set (with respect to emission scenario, metrics, etc.) of aCCFs (aCCF-V1.0A) has recently been developed and introduced within the EU project FlyATM4E (see Yin et al., 2022; Dietmüller et al., 2022; Matthes et al., 2023).

As for climate-optimal trajectory planning methods, various strategies ranging from mathematical programming (e.g., Campbell et al., 2008) to meta-heuristic (e.g., Yin et al., 2018a; Yamashita et al., 2020), indirect optimal control (e.g., Sridhar et al., 2011), and direct optimal control methods (e.g., Niklaß et al., 2017; Lührs et al., 2021, 2016; Matthes et al., 2020) have been adopted. For instance, the direct optimal control approach has been employed by Hartjes et al. (2016) to minimize flight time (or distance flown) in areas sensitive to persistent contrail formation, by Lührs et al. (2021) to minimize average temperature response over the next 20 years (ATR20) associated with non-CO₂ emissions, and by Vitali et al. (2021) to minimize the global warming potential (GWP) of NO_x, H₂O, soot, SO₂, and contrails. Using aCCFs to quantify climate impacts, Yamashita et al. employed a genetic algorithm to determine climate-optimal aircraft trajectories (Yamashita et al., 2020, 2021). Regarding the optimization methodologies, the mathematical programming methods only apply to simplified aircraft trajectory optimization problems (e.g., in the study conducted by Campbell et al., 2008, the aircraft's dynamic behavior is represented with a linearized model). The meta-heuristic methods (e.g., genetic algorithm) require very fast aircraft trajectory prediction in order to find an optimal solution with a large number of iterations; thus, the flight planning problem is usually approximated with a simplified but representative enough problem (e.g., in Yamashita et al., 2020, the optimization is defined with 11 decision variables to characterize lateral path and flight altitude, and the speed profile is considered constant). Finally, with the optimal control methods, the capability to model more accurate aircraft trajectory optimization problems is provided since the problem is represented as a dynamic optimization problem. Nevertheless, there are some drawbacks associated with addressing the formulated problem. The dynamic programming method (as an optimal control approach) results in the “curse of dimensionality” for complex problems (e.g., a full 4D aircraft trajectory optimization problem). Regarding the indirect optimal control approach, deriving analytical solutions using Pontryagin's maximum principle is daunting, especially for problems with singularities (e.g., only a 2D trajectory optimization problem has been addressed in the literature; Sridhar et al., 2011). The direct optimal control approach, despite being very flexible in modeling aircraft trajectory optimization problems (e.g., considering a full 4D dynamical model with nonlinear path and boundary constraints; Lührs et al., 2021), has a high sensitivity to initial conditions, and thus local optimality is its main drawback. In addition, considering the airspace structure with indirect and direct optimal control methods is not straightforward. Interested readers are referred to Simorgh et al. (2022) for our recent comprehensive survey on climate-optimal aircraft trajectory planning, reviewing both the approaches to model climate-sensitive regions and trajectory planning methods. A classification of the most recent studies in this field is provided in Table 1.

To quantify the non-CO₂ climate effects, specific weather variables are required. In the case of aCCFs, variables such as temperature (T), potential vorticity (PV), geopotential (Φ), relative humidity over ice (r_hum), and outgoing longwave radiation (OLR) are needed. These variables are obtained from standard weather forecasts. Several factors, including incomplete understanding of the state of the atmosphere, computational complexity, and nonlinear and sometimes chaotic dynamics, affect the accuracy of weather predictions, implying that the weather forecast is inevitably uncertain (WMO, 2012). These weather-forecast-related uncertainties in the aCCFs and also in aircraft dynamical behavior (e.g., uncertainty in wind and temperature), if not accounted for within aircraft trajectory planning, can lead to inefficient trajectories. Previous research in the field of climate-optimal aircraft trajectory planning has been conducted in a deterministic manner, neglecting the inclusion of any sources of uncertainty (see Table 1) (Simorgh et al., 2022). A first step in managing and integrating meteorological uncertainties into aircraft path planning is to obtain reliable weather forecasts that can predict probable variations in meteorological conditions. To characterize weather forecast uncertainties, probabilistic weather forecasting (PWF) is typically used (AMS-Council, 2008). State-of-the-art probabilistic weather forecasting is obtained from the ensemble prediction system (EPS), which provides N_EPS possible realizations of meteorological conditions called ensemble members (Bauer et al., 2015).

When accounting for the ensemble weather forecast in solving aircraft trajectory optimization, the computational time is an important issue that arises in addition to the capability to consider such uncertainties. This is due to the fact that, instead of taking one weather forecast and solving the trajectory optimization in a deterministic manner, the optimizer should be capable of considering N_EPS (e.g., N_EPS=50) different forecasts. Several studies have been proposed in the literature to determine robust aircraft trajectories in the presence of meteorological uncertainty quantified using EPS weather forecast (though not considering climate impact; Simorgh et al., 2022). However, these studies suffer mainly from computational perspectives (i.e., the computational time of the optimizer when considering weather uncertainty) and some restrictive assumptions (e.g., inaccurate modeling of the aircraft dynamical model or ignoring important decision variables such as the flight altitude) (Simorgh et al., 2022, Sect. 5.3). For instance, in the study conducted by González-Arribas et al. (2018), a robust aircraft trajectory optimization problem is solved using the optimal control approach within fully free-routing airspace. In this approach, the effects of uncertainty are included in the optimization problem by expanding the dynamical model of the aircraft (almost) linearly to the number of ensemble members, resulting in a larger dimensional deterministic optimization problem. Thus, it requires more computational time compared to the deterministic flight planning problem. Also, the flight planning problem is limited to optimizing only the lateral path. As for the structured airspace, Franco Espín et al. (2018) proposed uncertain flight plan optimization using mixed-integer linear programming (MILP). Similarly, the optimization in this study is performed in 2D airspace. In addition, the fuel burn and its associated nonlinearities are ignored, and it has a sizable computational cost of several minutes. In this respect, developing efficient trajectory optimization solvers capable of delivering robust climate-optimal trajectories with a computational time compatible with operations has been identified as a scientific gap (see Simorgh et al., 2022).

The focus of recent studies has been restricted to planning climate-optimal trajectories considering the concept of future free-route airspace (see last column Table 1), which is thus not applicable for the structured airspace of today. The trajectory optimization problem constrained by the structure of airspace results in hybrid decision spaces (e.g., route and flight level are discrete, and speed schedule is continuous) (González-Arribas et al., 2023). The trajectory optimization problem with the combination of discrete and continuous decision variables is one of the most challenging optimization problems, typically solved using mixed-integer nonlinear programming with intensive computational cost (see, e.g., Bonami et al., 2013).

Soler et al. (2014)Grewe et al. (2014a)Hartjes et al. (2016)Lührs et al. (2016)Lim et al. (2017)Matthes et al. (2017)Niklaß et al. (2017)Yin et al. (2018b)Yin et al. (2018a)Niklaß et al. (2019)Yin et al. (2022)Yamashita et al. (2020)Matthes et al. (2020)Lührs et al. (2021)Yamashita et al. (2021)

Table 1 A classification of the recent studies in the literature proposed to reduce the climate impact of aircraft emissions with aircraft trajectory optimization.

Drawing upon the brief literature review and the presented open problems, we aim to address the problem of determining robust climate-optimal aircraft trajectories within the structured airspace in this study. Our main contributions are summarized as follows: (1) full 4D climate-optimal trajectory planning within the currently structured airspace, (2) accounting for uncertain meteorological conditions and uncertainty associated with initial flight conditions such as departure time and aircraft initial mass, and (3) determining the optimized trajectory computationally very fast. The uncertainty in weather forecast is characterized using the ensemble prediction system, and aviation's climate impacts are quantified by employing the latest version of aCCFs (V1.0A). The concept of robustness that we refer to is the determination of the aircraft trajectory considering all possible realizations of meteorological variables provided using an EPS weather forecast. In other words, instead of planning a trajectory based on one forecast in a deterministic manner, we aim to determine a trajectory that is optimal considering the overall performance obtained from using different members of an ensemble weather forecast. In this respect, from the operational point of view, the optimized trajectory is tracked as determined, and the effects of meteorological uncertainties are reflected in the flight performance variables such as flight time, fuel burn, and climate impact. Mathematically, the perturbations due to the meteorological uncertainty are considered in the dynamical model of aircraft, and the proposed trajectory optimization is generic in terms of the objective function, in which a wide range of objectives, such as flight time, fuel consumption, emissions, and climate impact, and different statistics including expected values and variance of the performance variables, can be considered. Such flexibility in defining the cost function allows for solving a multi-objective optimization problem. Moreover, by penalizing the mean and variance of the objectives, the effects of uncertainty on flight variables can be controlled. In this study, the flight planning objective is a weighted sum of the simple operating cost (as a function of flight time and fuel consumption) and climate impact. The focus is restricted to optimizing the expected performance since, as will be shown in the simulation results, minimizing the averaged performance leads to reducing the uncertainty ranges for the considered case studies.

We employ the probabilistic flight planning method firstly developed by González-Arribas et al. (2023) to determine robust climate-optimal trajectories for three phases: climb, cruise, and descent. In this approach, to account for discrete and continuous decision variables in an integrated manner, the optimization is carried out on the space of probability distributions defined over flight plans instead of directly searching for the optimal profile. Then, the probability distribution over flight plans is parameterized, allowing the generation of multiple flight plans stochastically. The augmented random search algorithm is employed and implemented on graphics processing units (GPUs) to deliver a nearly optimal solution to the resulting stochastic optimization in seconds. We have developed an open-source Python library called ROOST V1.0 (Robust Optimization of Structured Trajectories) based on the proposed robust aircraft trajectory optimization technique, which is currently available via the DOI https://doi.org/10.5281/zenodo.7121862 (Simorgh, 2022). ROOST is a tool that efficiently uses the information provided by the prototype aCCFs (implemented in our recently developed Python library CLIMaCCF V1.0 (DOI: https://doi.org/10.5281/zenodo.6977273; Dietmüller et al., 2022) for planning climate-optimal trajectories accounting for the operational constraints and uncertainty. It should be noted that CLIMaCCF V1.0 is the first release of the Python library, which includes both versions of aCCFs, i.e., V1.0 and V1.0A. For performing aircraft trajectory optimization in this study, we use V1.0A aCCFs. Users need to input the required weather variables, route graphs, and aircraft type specifications to start working with the library. In addition, users should assign values to the weighting parameters associated with different objectives in the objective functions. This paper is mainly devoted to the climate-optimal aircraft trajectory planning algorithm implemented in the library and its application to optimize different case studies. Instructions to get started with the library can be found in the repository of ROOST.

The paper is arranged as follows. The robust climate-optimal aircraft trajectory planning problem is stated and formulated in Sect. 2 and solved in Sect. 3. The potential of our flight planning algorithm to reduce aviation-induced climate impacts under uncertain meteorological conditions is explored in Sect. 4. Finally, the discussion on the obtained results and concluding remarks are presented in Sect. 5.

2 Problem statement

Aviation-induced non-CO₂ climate effects have a direct dependency on atmospheric conditions at the location and time of emissions. In this respect, generating flight plans that could potentially minimize emissions in areas with high sensitivity to aircraft emissions in terms of climate change is a step towards climate-friendly air transport. To enable such a flight planning strategy, information on climate-sensitive areas is required first. Once climate information is available, it needs to be included in flight planning tools to generate eco-efficient trajectories. It is worth mentioning that to generate reliable trajectories, different sources of uncertainty need to be identified, understood, and quantified.

The focus of this study is on determining robust flight plans taking into account meteorological uncertainty and uncertainty associated with the initial flight conditions. In this respect, we start by presenting a general formulation of the dynamic optimization problem with uncertainty in Sect. 2.1, which is then employed to model the proposed robust flight planning (in Sect. 2.2). In Sect. 2.3, the effects of meteorological uncertainty on the efficiency of climate-optimal trajectories are discussed, and the motivation to solve robust trajectory optimization is provided.

2.1 General formulation of dynamic optimization problem with uncertainty

We present a general formulation of a dynamic optimization problem with the inclusion of uncertainty effects, focusing mainly on the dynamical model and objective function.

Let us consider a class of dynamical systems with uncertainty as follows:

$$\begin{array}{}\text{(1)}& \dot{\mathit{x}}(t,\mathit{\omega })=\mathbf{f}(t,\mathit{x}(t,\mathit{\omega }),\mathit{u}(t,\mathit{\omega }),\mathit{z}(t,\mathit{\omega }),\mathit{\zeta }(\mathit{\omega }\left)\right),\end{array}$$

where $\mathit{u}\in {\mathbb{R}}^{n_u}$, $\mathit{x}\in {\mathbb{R}}^{n_x}$, and $\mathit{z}\in {\mathbb{R}}^{n_z}$ are the vectors of control inputs as well as states and algebraic variables, respectively, and f is a vector field, mapping $\mathbb{R}\times {\mathbb{R}}^{n_z}\times {\mathbb{R}}^{n_x}\times {\mathbb{R}}^{n_u}\times {\mathbb{R}}^{n_{\mathit{\zeta }}}\to {\mathbb{R}}^{n_x}$. The uncertain parameters (denoted with $\mathit{\zeta }\in {\mathbb{R}}^{n_{\mathit{\zeta }}}$) are considered continuous random variables and assumed to have known probability distribution functions $\mathit{\zeta }(\cdot ):\mathrm{\Omega }\to {\mathbb{R}}^{n_{\mathit{\zeta }}}$, where Ω is a space of possible outcomes. The random variables take different values depending on the probable outcomes (i.e., ζ(ω) for ω∈Ω). The nonlinear function f(⋅) is assumed to be a measurable function in ζ. To emphasize the effects of random variables on the system's trajectories, all the variables were denoted with dependency on possible outcomes (e.g., x(t,ω)). For the sake of improving clarity, the abbreviated notation (e.g., x(t)) will be used in the following.

A general form of the cost functional considered for dynamic optimization problems with uncertainty is

$$\begin{array}{}\text{(2)}& \begin{array}{rl}J& =\mathbb{E}\mathit{\left\{}\mathcal{M}\right(t_{\mathrm{0}},\mathit{x}\left(t_{\mathrm{0}}\right),t_{\mathrm{f}},\mathit{x}\left(t_{\mathrm{f}}\right))\\ & +\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{f}}}{\int }}\mathcal{L}(t,\mathit{x}(t),\mathit{u}(t),\mathit{z}(t),\mathit{\zeta })\text{d}t\mathit{\}},\end{array}\end{array}$$

where $\mathcal{M}:\mathbb{R}\times {\mathbb{R}}^{n_x}\times \mathbb{R}\times {\mathbb{R}}^{n_x}\to \mathbb{R}$ and $\mathcal{L}:\mathbb{R}\times {\mathbb{R}}^{n_x}\times {\mathbb{R}}^{n_u}\times {\mathbb{R}}^{n_z}\times {\mathbb{R}}^{n_{\mathit{\zeta }}}\to \mathbb{R}$ are the Mayer and Lagrangian terms, respectively. Notice that the objective function and also dynamical model are similar to what is usually considered within the context of optimal control theory. Although the objective function is written in terms of using the expectation operator, other statistics can be evaluated under this formulation. For instance, one can include the variance of a function A(ζ) as $\mathbb{V}\mathit{\left\{}A\mathit{\right\}}=\mathbb{E}\mathit{\left\{}\right(A-\mathbb{E}\mathit{\left\{}A\mathit{\right\}}{)}^{\mathrm{2}}\mathit{\}}=\mathbb{E}\mathit{\{}{A}^{\mathrm{2}}\mathit{\}}-\mathbb{E}\mathit{\{}A{\mathit{\}}}^{\mathrm{2}}$, where 𝕍{⋅} is the variance operator. Depending on the benchmark problem, several constraints such as equality and inequality path and boundary constraints are considered (see Simorgh et al., 2022).

The objective of the stated dynamic optimization problem is to find a feasible control policy (u∈𝒰) to minimize the performance index (Eq. 2) respecting a set of constraints, including dynamical constraints (Eq. 1). Depending on the benchmark problem, reformulations and approximations are normally made to address the required performance, such as computational complexity. In this study, we will slightly reformulate the optimization for the proposed path planning problem.

2.2 Robust climate-optimal flight planning problem formulation

The definition of the aircraft trajectory optimization problem mainly requires the aircraft dynamical model, flight planning objectives, and physical and operational constraints. To consider climate impact within aircraft trajectory planning, information on the climate impacts of CO₂ and non-CO₂ emissions is necessary and needs to be included in the objective function. In the following, the modeling of the mentioned elements is briefly presented.

2.2.1 Aircraft dynamical model

To determine reliable aircraft trajectories, accurate aircraft dynamical models are necessary. In this work, the point-mass model with the following equations of motion is used to represent the aircraft's dynamical behavior, as is usually considered within air traffic management studies.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& \left[\begin{array}{c}\dot{\mathit{\varphi }}\\ \dot{\mathit{\lambda }}\\ \dot{h}\\ \dot{v}\\ \dot{m}\end{array}\right]=\left[\begin{array}{c}(v\cos \mathit{\gamma }\cos \mathit{\chi }+w_y(\mathit{\omega }\left)\right)\left(R_M\right(\mathit{\varphi })+h{)}^{-\mathrm{1}}\\ (v\cos \mathit{\gamma }\sin \mathit{\chi }+w_x(\mathit{\omega }\left)\right)\left(\right(R_N\left(\mathit{\varphi }\right)+h)\cos \mathit{\varphi }{)}^{-\mathrm{1}}\\ v\sin \mathit{\gamma }\\ \left(T\right(C_T)-D(C_L\left)\right)m^{-\mathrm{1}}-g\sin \mathit{\gamma }\\ -\text{FF}\left(T\right(\mathit{\omega }),C_T)\end{array}\right],\\ & \text{State variables }\left(\mathit{x}\right):{\left[\begin{array}{ccccc}\mathit{\varphi }& \mathit{\lambda }& h& v& m\end{array}\right]}^T,\\ & \text{Control variables }\left(\mathit{u}\right):{\left[\begin{array}{ccc}{C}_{\mathrm{T}}& \mathit{\chi }& \mathit{\gamma }\end{array}\right]}^T.\end{array}\end{array}$$

Here, ϕ is the latitude, λ is the longitude, v is the true airspeed, h is the altitude, m is the mass, C_T is the thrust coefficient, γ is the climb angle, χ is the heading angle, $C_L\left(\mathit{\gamma }\right)=\left(\mathrm{2}mg\cos \mathit{\gamma }\right)/\left(\mathit{\rho }{v}^{\mathrm{2}}S\right)$, and (w_x,w_y) are wind components. The Earth's ellipsoid radii of curvature in the meridian and the prime vertical are denoted with R_M and R_N, respectively, T is the magnitude thrust force, and the drag force is denoted with D. g is the Earth's gravity, FF is the fuel flow, and S is the wetted surface of the aircraft. The BADA4.2 model (Gallo et al., 2006) is employed to provide the aerodynamic and propulsive performance of the aircraft. Notice that we used the notation presented in Sect. 2.1 (i.e., ω) to represent the uncertainty in the temperature and components of wind. The realization of the state and control variables depends on uncertainty, e.g., for mass, $m=m(t,\mathit{\omega })$; in the interest of notational clarity, we use the abbreviated notation m.

Structured airspace

As trajectory optimization in this study is performed within the structured airspace, the evolutions of aircraft states are constrained. In the following, we briefly present our proposed modeling of airspace structure and flight plan.

We consider a directed acyclic graph $G=(V,E)$ to model the airspace with V as the navigation waypoints and e∈E as the airway edges connecting waypoints. The trajectory is assumed to start at the end of the standard instrument departure procedure (SID) and end at the beginning of the standard instrument arrival procedure (STAR) to the destination airport, denoted as o∈V and d∈V, respectively. We define the flight plan ℱ with a tuple $(\mathcal{R},\stackrel{\mathrm{‾}}{\text{FL}},\stackrel{\mathrm{‾}}{\text{M}},\mathcal{C},\mathcal{D},d_{\mathcal{D}})$. In the flight plan (ℱ), the route (or lateral path) denoted as ℛ includes a sequence of waypoints, i.e., $\mathcal{R}:=(r_{\mathrm{0}},r_{\mathrm{1}},\mathrm{\cdots },r_{n_r})$. The vertical profile of the cruise, $\stackrel{\mathrm{‾}}{\text{FL}}$, is composed of an ordered sequence of tuples of the form (r_k,FL_k), indicating that, if the aircraft is in the cruise phase, the flight level will be changed to FL_k when the waypoint r_k is reached (see Fig. 1). The Mach schedule $\stackrel{\mathrm{‾}}{\text{M}}:=({\text{M}}_{\mathrm{0}},{\text{M}}_{\mathrm{1}},\mathrm{\cdots },{\text{M}}_{n_r})$ indicates the target Mach number M_k at waypoint r_k during the cruise phase. Climb and descent profiles $\mathcal{C},\mathcal{D}:\mathbb{R}\to \mathbb{R}$ are represented by continuous and piecewise-differentiable functions mapping the altitude to the target airspeed during the climb and descent phases, respectively. Finally, a scalar variable d_𝒟 shows the distance to go to the destination node at which the aircraft should end the cruise and start the descent phase.

Figure 1 Structure of airspace. The route graph is generated by processing the full airspace graph to include paths from the end of the SID to the beginning of the STAR to the destination airport.

2.2.2 Objective function

The goals of the aircraft trajectory optimization problem are interpreted mathematically and defined as an objective function to be minimized (or maximized). In addition to the climate impact, the operating cost is a crucial aspect that needs to be considered as it is one of the main interests of airliners. Generally, there is a trade-off between the operating cost and climate impacts. This is due to the fact that rerouting areas sensitive to climate increases operational costs as the aircraft tends to fly longer routes (Niklaß et al., 2021). To consider both objectives in the trajectory optimization, we define the following objective function:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\text{Objective function}& ={\mathit{\psi }}_{\mathrm{CST}}\cdot \text{Operating cost}\\ & +{\mathit{\psi }}_{\mathrm{CLM}}\cdot \text{Climate impact},\end{array}\end{array}$$

where ψ_CST and ψ_CLM are weighting parameters penalizing cost and climate impact, respectively. In the following, the proposed modeling of these two objectives is presented.

Operating cost

There are various approaches to account for operating costs within aircraft path planning. Flight time and/or fuel are common objectives. However, more realistic cost metrics exist, which include additional costs such as flight crew, cabin crew, and landing fees. Interested readers are referred to Simorgh et al. (2022) for a classification of these metrics.

In this study, we use simple operating cost (SOC) as a metric expressing cost in USD with linear relation to flight time and fuel consumption:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \text{Operating cost}={\mathit{\psi }}_t\cdot \text{ Exp. Flight time}+{\mathit{\psi }}_m\cdot \text{Exp. Fuel burnt}\\ & \text{ Exp. Flight time}:\mathbb{E}\mathit{\left\{}\text{FT}\mathit{\right\}}:=\mathbb{E}\mathit{\left\{}{t}_{\mathrm{f}}\right(\mathit{\omega })-t_{\mathrm{0}}(\mathit{\omega }\left)\mathit{\right\}}\\ & \text{ Exp. Fuel burnt}:\mathbb{E}\mathit{\left\{}\text{FB}\mathit{\right\}}:=\mathbb{E}\mathit{\left\{}{m}_{\mathrm{0}}\right(\mathit{\omega })-m_{\mathrm{f}}(\mathit{\omega }\left)\mathit{\right\}},\end{array}\end{array}$$

where t₀ and t_f are the initial flight time and final flight time weighted by ψ_t=0.75 USD s⁻¹, and m₀ and m_f are the initial mass and final mass weighted by ψ_m=0.51 [USD kg⁻¹]. The expectation operator is used here because the meteorological uncertainty (included in the aircraft dynamical model) and also uncertainty due to initial conditions directly affect the flight time and fuel consumption.

In spite of considering only flight time and fuel consumption to represent the operating cost, it was reported in Table 4 of Yamashita et al. (2020) that employing SOC and a more comprehensive metric such as cash operating cost within trajectory optimization delivered almost similar results. This is because these two metrics, in the end, consider flight time and fuel consumption as inputs to estimate the operating cost.

Climate effects induced by aviation

Numerous approaches have been proposed in the literature to consider climate impact within aircraft trajectory planning strategies (see Table 1 and Fig. 2 of Simorgh et al., 2022). In this study, the climate impact of aircraft operations is modeled using the prototype aCCFs. These aCCFs take as input specific meteorological variables (e.g., temperature and relative humidity) and estimate the climate impact in terms of average temperature change. The suitability of aCCFs for climate-optimal trajectory planning can be justified as follows:

• aCCFs account for the temporal and spatial dependency of climate impacts associated with non-CO₂ species, including ozone and methane, resulting from NO_x emissions, water vapor emissions, and persistent contrails;

• aCCFs estimate the climate impact associated with aircraft emissions computationally in real time, making it well-suited for climate-optimal trajectory planning; and

• aCCFs directly quantify climate impacts in average temperature change.

The aCCF V1.0 reported by Yin et al. (2022) is the first complete and consistent set of prototype aCCFs providing spatially and temporally dependent non-CO₂ climate effects in terms of average temperature response over the next 20 years for a pulse emission scenario (P-ATR20). In this study, we use the latest version of aCCFs, i.e., aCCF V1.0A developed by Matthes et al. (2023) within the EU project FlyATM4E, which has been calibrated to the state-of-the-art climate response model AirClim (Dahlmann et al., 2016).

The V1.0A aCCFs are multiplied by some factors in order to be more suitable for the flight planning application.

• Metric conversion. The selection of a suitable metric depends on the question to be answered (see Grewe and Dahlmann, 2015, for more details); therefore, different questions require the use of different metrics. The P-ATR20 metric was selected as a metric for the aCCFs to assess the impact of a “simple” pulse emission. However, factors are available (see Dietmüller et al., 2022) to convert P-ATR20 to other metrics: for example, assuming the future emission scenario or longer time horizons. This way, one can select the emission scenario and time horizon that are best-suited for the question. In this study, the F-ATR20 metric is used to assess the climate effect reduction obtained by steadily applying a specific routing strategy under the assumption of a future business-as-usual emission scenario. To this end, the aCCFs based on the pulse emission scenario are converted to the future scenario (F-ATR20) using values reported in Table 3 of Dietmüller et al. (2022). These conversion factors were derived by simulations with the climate response model AirClim (Dahlmann et al., 2016): one simulation with pulse emission and one with the future emission scenario. By comparing the two simulations, the factors can be derived.

• Efficacy. Efficacies were introduced to take into account the fact that the radiative forcing of some non-CO₂ forcing agents (e.g., ozone, methane, contrails) is more or less effective in changing the global mean near-surface temperature per unit forcing compared to the response of CO₂ forcing (Hansen et al., 2005; Ponater et al., 2007). Dietmüller et al. (2022) have summarized the efficacy parameters reported by Lee et al. (2021). For a detailed explanation of efficacy, the reader is referred to the state-of-the-art literature (e.g., Ponater et al., 2007; Rap et al., 2010; Bickel et al., 2020).

For the sake of compactness of notation, F-ATR20 (from aCCF V1.0A) with efficacy parameters is replaced with ATR in the following.

For the aCCFs of (daytime and nighttime) contrails, the ice supersaturation is applied using temperature and relative humidity over ice in order to predict regions where persistent contrails are expected to form, called persistent contrail formation areas (PCFAs) (Schmidt, 1941; Appleman, 1953). To represent the climate impacts using aCCFs in the average temperature change (i.e., [K]), fuel consumption rate, NO_x emissions, and distance flown through PCFAs are required.

The geographical aCCF pattern of water vapor, NO_x-induced ozone (production), and methane (destruction), as well as of contrails, is shown in Fig. 2a–c for 13 June 2018 at 00:00 UTC over the European region for FL340. Moreover, Fig. 2d shows the pattern of the merged non-CO₂ aCCFs that combines the individual aCCFs (Fig. 2a–c). Note that to generate the merged aCCFs and to compare the contribution of each species, we adopt typical transatlantic fleet mean values to unify the units of aCCFs in K kg⁻¹(fuel). The approximated conversion factors for NO_x emissions and contrails are $\mathrm{13}\times {\mathrm{10}}^{-\mathrm{3}}$ kg(NO₂) kg⁻¹(fuel) and 0.16 km kg⁻¹(fuel), respectively (Graver and Rutherford, 2018; Penner et al., 1999). It is clear from the merged aCCFs that the contrails have dominant climate effects, which is in line with related studies employing aCCFs (see, e.g., Dietmüller et al., 2022).

Figure 2 Algorithmic climate change functions of (a) water vapor, (b) NO_x, (c) contrails, and (d) the total non-CO₂ effects on 13 June 2018 at 00:00 UTC over European airspace for FL340.

To benefit from the spatial and temporal dependency of non-CO₂ climate effects identified using aCCFs in planning climate-optimal trajectories, we define the following objective expressed in the Lagrangian form for Eq. (4).

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \text{Climate impact}=\text{ Exp. ATR}:\mathbb{E}\mathit{\left\{}\text{ATR}\mathit{\right\}}\\ & :=\mathbb{E}\mathit{\{}\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{f}}}{\int }}\sum _{i=\mathrm{1}}^{\mathrm{5}}{\mathit{\psi }}_{\mathrm{ATR},i}\cdot {\text{ATR}}_i(t,\mathit{x}\left(t\right),\mathit{u}\left(t\right),\mathit{\zeta }\left)\text{d}t\mathit{\right\}}\end{array}\end{array}$$

for $i\in \mathit{\{}{\text{CH}}_{\mathrm{4}},\text{Cont.},{\text{O}}_{\mathrm{3}},{\text{H}}_{\mathrm{2}}\text{O},{\text{CO}}_{\mathrm{2}}\mathit{\}}$:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\text{ATR}}_{{\mathrm{O}}_{\mathrm{3}}}(t,\mathit{x},\mathit{u},\mathit{\zeta })={\mathrm{10}}^{-\mathrm{3}}\times {\text{aCCF}}_{{\mathrm{O}}_{\mathrm{3}}}(t,\mathit{x},\mathit{\zeta })\cdot {\dot{m}}_{nox}(t,\mathit{\zeta })\\ & {\text{ATR}}_{{\mathrm{CH}}_{\mathrm{4}}}(t,\mathit{x},\mathit{u},\mathit{\zeta })={\mathrm{10}}^{-\mathrm{3}}\times {\text{aCCF}}_{{\mathrm{CH}}_{\mathrm{4}}}(t,\mathit{x},\mathit{\zeta })\cdot {\dot{m}}_{nox}(t,\mathit{\zeta })\\ & {\text{ATR}}_{\mathrm{Cont}.}(t,\mathit{x},\mathit{u},\mathit{\zeta })={\mathrm{10}}^{-\mathrm{3}}\times {\text{aCCF}}_{\mathrm{Cont}.}(t,\mathit{x},\mathit{\zeta })\cdot v_{\mathrm{gs}}(t,\mathit{\zeta })\\ & {\text{ATR}}_{{\mathrm{H}}_{\mathrm{2}}\text{O}}(t,\mathit{x},\mathit{u},\mathit{\zeta })=-{\text{aCCF}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}(t,\mathit{x},\mathit{\zeta })\cdot \dot{m}(t,\mathit{\zeta })\\ & {\text{ATR}}_{{\mathrm{CO}}_{\mathrm{2}}}(t,\mathit{x},\mathit{u},\mathit{\zeta })=-{\text{aCCF}}_{{\mathrm{CO}}_{\mathrm{2}}}\cdot \dot{m}(t,\mathit{\zeta })\end{array}\end{array}$$

Here, $\mathit{\zeta }(\cdot ):=\left[w_x\right(\cdot \left)w_y\right(\cdot \left)T\right(\cdot \left)\text{OLR}\right(\cdot \left)\text{PV}\right(\cdot \left)\mathrm{\Phi }\right(\cdot ){\text{r}}_{\mathrm{hum}}$ ${\text{F}}_{\mathrm{in}}(\cdot )\text{q}(\cdot )t_{\mathrm{0}}(\cdot )m_{\mathrm{0}}(\cdot )]$ is a vector of uncertain variables, i.e., meteorological variables and initial flight conditions. In addition, ${\dot{m}}_{{\mathrm{no}}_x}\left(t\right)=\text{FF}(t,\mathit{u})\cdot {\text{EI}}_{{\mathrm{NO}}_x}(t,\mathit{x},\mathit{u})$, where ${\text{EI}}_{{\mathrm{NO}}_x}(\cdot )$ is the actual NO_x emission index in [g(NO₂) kg⁻¹(fuel)], and v_gs is the ground speed. The NO_x emission index varies with many factors such as aircraft type, fuel flow, flight altitude, and synoptical situation. To consider such dependencies, the Boeing fuel flow method 2 (BFFM2), calculating the actual emission index of NO_x from the reference conditions, is adopted (DuBois and Paynter, 2006; Jelinek, 2004). The reference conditions are obtained from the International Civil Aviation Organization (ICAO) data bank (https://www.easa.europa.eu/en/domains/environment/icao-aircraft-engine-emissions-databank, last access: 2 July 2023). Notice that the objective regarding climate impact is expressed in the Lagrangian form (i.e., as an integral) since the climate impact needs to be evaluated (and accumulated) along the route, unlike the operating cost, which requires only information on boundary values (e.g., initial and final flight time). Notice that we define the objectives considering the expected performance. One can use other statistics, such as variance, without loss of generality. For instance, to penalize the range of uncertainty in climate impacts, $\mathbb{V}\mathit{\left\{}\text{ATR}\mathit{\right\}}=\mathbb{E}\mathit{\left\{}{\text{ATR}}^{\mathrm{2}}\mathit{\right\}}-\mathbb{E}\mathit{\{}\text{ATR}{\mathit{\}}}^{\mathrm{2}}$ should be included in the objective function. For the considered case studies, it will be shown that minimizing the expected objective function will reduce the uncertainty in the most uncertain variable for the climate-optimal routing option, thus providing a robust solution without penalizing the deviation.

2.3 Effects of weather-induced uncertainty on climate

The dynamical model of aircraft requires weather-related variables such as wind and temperature (see, e.g., Sect. 2.2.1). In addition, the non-CO₂ climate effects of aviation included in the objective function of the optimization problem strongly rely on weather conditions. Since the required meteorological variables are obtained from standard weather forecasts, they are inevitably uncertain. It is worth mentioning that there are other sources of uncertainty affecting the efficiency of the planned climate-optimized trajectories, including uncertainty from climate science (e.g., modeling and estimating aviation-induced climate effects), emission calculation, and also inaccurate modeling of aircraft behavior, which are not within the scope of this paper but have been identified as open problems (see Matthes et al., 2023).

In this paper, the focus is on forecast-related uncertainties, which will be characterized by employing ensemble prediction system (EPS) weather forecasts, a numerical weather prediction method introduced to deal with uncertainty in weather forecast (Bauer et al., 2015). These are forecasts in which both the initial conditions and the physical parameters of a numerical weather integration model are slightly modified from one member to another and provide N_EPS (typically N_EPS = 50) different predictions known as ensemble members (Bauer et al., 2015). Each member of an ensemble represents one possible realization of meteorological situations. As the aCCFs take as inputs some meteorological variables, N_EPS different aCCFs can be calculated for an EPS weather forecast. For instance, the meteorological variables temperature and relative humidity over ice are required for the aCCF of nighttime contrails. Notice that relative humidity over ice is required for identifying ice-supersaturated areas. Feeding N_EPS probable realizations of these meteorological variables (i.e., ensemble members), N_EPS different aCCFs (i.e., ${\text{aCCF}}_{{\mathrm{Cont}}_i}$ for $i=\mathrm{1},\mathrm{\cdots },N_{\mathrm{EPS}}$) are calculated. The same applies to system dynamics (considering uncertainty in temperature and wind) and also the NO_x emission index (due to the dependency on ambient temperature and specific humidity).

To investigate the degree of uncertainty (or variability) in the meteorological variables provided by an EPS and its effects on the computed aCCFs, we take the standard deviation (SD) from 10 ensemble members of weather data obtained using the ERA5 reanalysis data products (https://cds.climate.copernicus.eu/, last access: 2 July 2023). It should be noted that the reanalysis data products are generated from post-processing with observations. Thus, the variability among the ensemble members is expected to be lower than the forecast data with more ensemble members yet still valid to illustrate. Figure 3 shows the SD of weather variables required to calculate the aCCFs and aircraft trajectory on 13 June 2018 at 00:00 UTC for FL340. The SD is taken over the normalized variables (with respect to their maximum values) for comparison purposes. The variability of geopotential, temperature, and wind is small compared to potential vorticity, outgoing longwave radiation, and relative humidity. The SDs of the calculated aCCFs based on the ensemble members are illustrated in Fig. 4. Since the aCCF of NO_x emissions (i.e., methane and ozone) depends on geopotential and temperature, its SD is small compared to the aCCFs of water vapor and nighttime contrails, which are based on potential vorticity and relative humidity (when applying the ice-supersaturated condition), respectively. Notice that the uncertainty in the climate impact of contrails is much higher than water vapor due to the variability of relative humidity in satisfying the persistency condition of contrails (see SD of PCFA in Fig. 4). Due to the considered time (i.e., 00:00 UTC), the aCCF of contrails is based on the formulation of nighttime aCCFs. In spite of negligible uncertainty in the aCCF of NO_x and also relatively low uncertainty in the aCCF of water vapor compared to aCCF of contrails, due to the dominant climate impact of contrails, the net non-CO₂ climate effect is considerably uncertain (see SD of the merged aCCFs in Fig. 4), which must be crucially taken into consideration.

Figure 3 Variability (quantified using SD) of the meteorological conditions for an ensemble weather forecast with 10 ensemble members on 13 June 2018 at 00:00 UTC over European airspace for FL340.

Figure 4 Variability (quantified using SD) of aCCFs for an ensemble weather forecast with 10 ensemble members on 13 June 2018 at 00:00 UTC over European airspace for FL340.

Figure 5 shows how uncertainty in meteorological variables can affect the performance of aircraft trajectories. As can be seen, these uncertainties are accumulated and can considerably degrade the efficiency of the optimized trajectory if not considered in the aircraft trajectory planning a priori. No recent study on the determination of climate-optimized trajectories has considered robustness in the sense of uncertainty in weather forecasts. One of the main reasons for not considering such variations is the computational time that arises within the optimization techniques. In fact, instead of considering one member, the optimization, in this case, is to consider N_EPS ensemble members (e.g., 50) and find an optimized trajectory to be optimal with respect to all probable realizations of meteorological variables. Such an increase in dimensions is daunting to cope with by employing classical dynamic optimization approaches such as direct optimal control.

In Sect. 3, we will address the problem of robust climate-optimal trajectory planning with an efficient heuristic method implemented on GPUs.

Figure 5 Propagation of uncertainties (associated with initial flight conditions and meteorological variables) within climate-optimal aircraft trajectory planning.

3 Solution approach

The aircraft trajectory optimization problem formulated in Sect. 2.2 is solved by employing the method firstly developed by González-Arribas et al. (2023), which is a heuristic optimization technique for the structured airspace and is capable of determining the optimized trajectory in four dimensions, i.e., latitude, longitude, altitude, and time.

As mentioned in Sect. 2.1, depending on the problem, some reformulations and approximations are normally made to the dynamic optimization problem to address the required objectives. Here, instead of seeking the optimal control policy (i.e., u^o), the goal is to find an optimal flight plan ℱ^o, i.e., $({\mathcal{R}}^{\mathrm{o}},{\stackrel{\mathrm{‾}}{\text{FL}}}^{\mathrm{o}},{\stackrel{\mathrm{‾}}{\text{M}}}^{\mathrm{o}},{\mathcal{C}}^{\mathrm{o}},{\mathcal{D}}^{\mathrm{o}},d_{\mathcal{D}}^{\mathrm{o}})$ (see Sect. 2.2.1) that minimizes the objective function given in Eq. (4) and satisfies the aircraft dynamical model (given in Eq. 3), airspace structure, and path and boundary constraints. Such a selection of decision space allows us to directly account for the operational restrictions, such as determining lateral routes that follow the airspace structure.

3.1 Ensemble trajectory integration

To determine the performance of a flight plan and evaluate the cost function Eq. (4), the corresponding trajectories of the aircraft are to be calculated using the aircraft dynamical model (provided in Sect. 2.2.1). As mentioned in Sect. 2.3, aircraft trajectories are affected by uncertainty in meteorological variables, including temperature and wind. Assuming a unique lateral path (constant course) to be tracked in practice with low-level controllers in real time, the uncertainty in the wind (both magnitude and direction) will affect ground speed and, consequently, the time the aircraft flies the route (see Fig. 6). In addition, uncertainty in temperature affects fuel consumption because the propulsive and aerodynamic performance of the aircraft and also airspeed have a dependency on temperature (González-Arribas et al., 2023). From Eq. (6), one can conclude that the uncertainty in flight time and flight mass can also affect the climate impacts (see also Fig. 5). In addition to the uncertainty in meteorological variables, which is characterized using ensemble weather forecasts, uncertainty associated with initial flight conditions is taken into account within the trajectory optimization problem. In this study, the initial flight time and initial mass of the aircraft are modeled as Gaussian variables, i.e., $t_{\mathrm{0}}\sim \mathcal{N}({\overline{t}}_{\mathrm{0}},{\mathit{\sigma }}_{t_{\mathrm{0}}})$ and $m_{\mathrm{0}}\sim \mathcal{N}({\overline{m}}_{\mathrm{0}},{\mathit{\sigma }}_{m_{\mathrm{0}}})$.

Figure 6 Relationship between wind, course, heading, airspeed, and ground speed.

To efficiently reflect the effects of wind uncertainty in flight performance variables, instead of time, the distance flown along the route (s) is considered to be the independent variable using $\left(\text{d}t\right)(\text{d}s{)}^{-\mathrm{1}}=v_{\mathrm{gs}}^{-\mathrm{1}}$. This is beneficial since the uniqueness of time for all possible realizations of wind means that the position of the aircraft is fixed with respect to time. In this case, the effects of uncertainty cannot be considered efficiently because the range of feasible solutions is limited. The selection of distance flown as the independent variable allows for reflecting wind uncertainty in flight time.

According to the defined objective function (Eqs. 4, 5, 6), the flight performance variables required to evaluate Eq. (4) are flight time, final mass, and climate impact. By using TI(⋅) to denote the integration of the aircraft dynamical model for a given flight plan, weather data, and initial flight conditions, we receive the expected final mass, the expected final time, and the expected ATR required to evaluate the performance of aircraft trajectory (i.e., calculate the objective function given in Eq. 4) as follows:

$$\begin{array}{}\text{(8)}& \left[\mathbb{E}\mathit{\right\{}\text{FT}\mathit{\}},\mathbb{E}\mathit{\{}\text{FB}\mathit{\}},\mathbb{E}\mathit{\{}\text{ATR}\mathit{\left\}}\right]=\mathbb{E}\mathit{\left\{}\text{TI}\right(\mathcal{F},\mathcal{W},t_{\mathrm{0}},m_{\mathrm{0}}\left)\mathit{\right\}}.\end{array}$$

In the formulated robust optimization problem, the uncertain variables are considered to be continuous random variables with known probability distribution functions. As the weather variables for an ensemble weather forecast are directly represented in a discrete fashion, we assume a discrete probability distribution for the uncertainty. Considering N_EPS probable realizations of uncertainty (i.e., $\mathit{\{}{\mathit{\zeta }}_j{\mathit{\}}}_{\mathrm{1}}^{N_{\mathrm{EPS}}}$), we have

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathit{\zeta }}^j& :=\left[{\mathcal{W}}^j{t}_{\mathrm{0}}^j{m}_{\mathrm{0}}^j\right]\\ & :=\left[w_x^j{w}_y^j{T}^j{\mathrm{\Phi }}^j{\text{r}}_{\mathrm{hum}}^j{\text{q}}^j{\text{PV}}^j{\text{OLR}}^j{\text{q}}^j{t}_{\mathrm{0}}^j{m}_{\mathrm{0}}^j\right],\end{array}\end{array}$$

where 𝒲^j(⋅) is a set of meteorological variables required for climate-optimal aircraft trajectory planning corresponding to the jth member of an EPS weather forecast, and $t_{\mathrm{0}}^j,m_{\mathrm{0}}^j\sim t_{\mathrm{0}}$, and m₀ are initial flight conditions that are sampled independently. Generally, different members of the ensemble weather forecasts are considered equally probable. This implies that a specific forecasted weather pattern that has a higher probability will be represented by a larger number of ensemble members. In this study, we use equal weights for each probable realization of weather variables, i.e., with a probability of $\mathbf{P}(\mathcal{W}={\mathcal{W}}^j)=N_{\mathrm{EPS}}^{-\mathrm{1}}$. Thus, using an unweighted average between all ensemble members, Eq. (8) can be written as

$$\begin{array}{}\text{(10)}& \left[\mathbb{E}\mathit{\right\{}\text{FT}\mathit{\}},\mathbb{E}\mathit{\{}\text{FB}\mathit{\}},\mathbb{E}\mathit{\{}\text{ATR}\mathit{\left\}}\right]\approx {\displaystyle \frac{\mathrm{1}}{N_{\mathrm{EPS}}}}\sum _{j=\mathrm{1}}^{N_{\mathrm{EPS}}}\text{TI}(\mathcal{F},{\mathcal{W}}^j,t_{\mathrm{0}}^j,m_{\mathrm{0}}^j).\end{array}$$

The expected ATR is calculated as

$$\mathbb{E}\mathit{\left\{}\text{ATR}\mathit{\right\}}\approx {\displaystyle \frac{\mathrm{1}}{N_{\mathrm{EPS}}}}\sum _{j=\mathrm{1}}^{N_{\mathrm{EPS}}}{\text{ATR}}^j;{\text{ATR}}^j=\sum _{i=\mathrm{1}}^{\mathrm{5}}{\text{ATR}}_i^j$$

for $i\in \mathit{\{}{\text{CH}}_{\mathrm{4}},\text{Cont.},{\text{O}}_{\mathrm{3}},{\text{H}}_{\mathrm{2}}\text{O},{\text{CO}}_{\mathrm{2}}\mathit{\}}$. For instance, ATR^j for ozone and contrails can be calculated as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\text{ATR}}_{{\mathrm{O}}_{\mathrm{3}}}^j& ={\mathrm{10}}^{-\mathrm{3}}\underset{t_{\mathrm{0}}^j}{\overset{t_{\mathrm{f}}^j}{\int }}{\text{aCCF}}_{{\mathrm{O}}_{\mathrm{3}}}^j\left({\mathit{x}}^j\right(t^j),t^j)\\ & \times \stackrel{{\dot{m}}_{{\mathrm{no}}_x}^j\left(t^j\right)}{\overbrace{\text{FF}\left({\mathit{x}}^j\right(t^j,t^j),{\mathit{u}}^j(t^j\left)\right)\times {\text{EI}}_{{\mathrm{NO}}_x}^j\left({\mathit{x}}^j\right(t^j),{\mathit{u}}^j(t^j,t^j)}})\text{d}{t}^j\\ {\text{ATR}}_{\mathrm{Cont}.}^j& ={\mathrm{10}}^{-\mathrm{3}}\underset{t_{\mathrm{0}}^j}{\overset{t_{\mathrm{f}}^j}{\int }}{\text{aCCF}}_{\mathrm{Cont}.}^j\left({\mathit{x}}^j\right(t^j),t^j)v_{\mathrm{gs}}^j\left(t^j\right)\text{d}{t}^j\\ & ={\mathrm{10}}^{-\mathrm{3}}\underset{\mathrm{0}}{\overset{s_{\mathrm{f}}}{\int }}{\text{aCCF}}_{\mathrm{Cont}.}^j\left({\mathit{x}}^j\right(t^j\left(s\right)),t^j(s\left)\right)\text{d}s,\end{array}\end{array}$$

where x^j(t^j) and u^j(t^j) are the state and control variables of the aircraft considering the jth realization of weather variables and the jth sampled initial conditions, $\text{d}s=v_{\mathrm{gs}}^j\cdot \text{d}{t}^j$, and

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\text{aCCF}}_{{\mathrm{O}}_{\mathrm{3}}}^j\left({\mathit{x}}^j\right(t^j),t^j)& :={\text{aCCF}}_{{\mathrm{O}}_{\mathrm{3}}}\left(T^j\right({\mathit{x}}^j\left(t^j\right),t^j),\\ & {\mathrm{\Phi }}^j\left({\mathit{x}}^j\right(t^j),t^j))\\ {\text{aCCF}}_{\mathrm{Cont}.}^j\left({\mathit{x}}^j\right(t^j),t^j)& :={\text{aCCF}}_{\mathrm{Cont}.}\left(T^j\right({\mathit{x}}^j\left(t^j\right),t^j),\\ & {\text{OLR}}^j\left({\mathit{x}}^j\right(t^j),t^j),{\text{r}}_{\mathrm{hum}}^j\left({\mathit{x}}^j\right(t^j),t^j)),\end{array}\end{array}$$

where weather variables such as T^j and Φ^j are the jth members of the ensemble weather forecast. The actual NO_x emission index, i.e., ${\text{EI}}_{{\mathrm{NO}}_x}^j\left({\mathit{x}}^j\right(t^j),{\mathit{u}}^j(t^j\left)\right)$, is calculated using BFFM2. As can be seen in Eq. (11), the climate impact due to the NO_x emissions depends on the amount of NO_x emitted in NO_x-sensitive regions (i.e., multiplied by the aCCF of NO_x emissions), while for contrails, it depends on the distance flown in persistent contrail formation areas.

Heun's method is adopted for integrating the aircraft dynamical model along discretized segments of the route through each phase, i.e., climb, descent, and cruise (González-Arribas et al., 2023). Since the calculations are similar for different members, parallelization would be beneficial in reducing computational time. Here, CUDA (Guide, 2013; Klöckner et al., 2012), a tool for general-purpose computing on the graphics processing unit, is employed to parallelize the computations.

3.2 Performance evaluation of a flight plan

The expected values obtained from Eq. (10) are for a specific flight plan. By these settings, for this flight plan, the cost function given in Eq. (4) can be evaluated with the following equation:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}J\left(\mathcal{F}\right)& ={\mathit{\psi }}_{\mathrm{CST}}[{\mathit{\psi }}_t\cdot \mathbb{E}\mathit{\{}\text{FT}\mathit{\}}+{\mathit{\psi }}_m\cdot \mathbb{E}\mathit{\{}\text{FB}\mathit{\left\}}\right]\\ & +{\mathit{\psi }}_{\mathrm{CLM}}\cdot \mathbb{E}\mathit{\left\{}\text{ATR}\mathit{\right\}}.\end{array}\end{array}$$

Figure 7 shows how the expected performance is calculated and evaluated for a given flight plan and an ensemble weather forecast.

Figure 7 Calculation and evaluation of the expected performance for a given flight plan and an ensemble weather forecast.

The objective is to find a flight plan that minimizes Eq. (13). Since the flight plan includes both discrete and continuous decision variables, the optimizer should be capable of solving the optimization within hybrid decision spaces. A classical approach to solving such optimization problems is mixed-integer nonlinear programming, which is mathematically complex and computationally intensive.

3.3 Probabilistic flight plan generation

To solve the optimization with hybrid decision variables in an efficient manner, instead of directly searching for the optimal flight plan, the optimization is conducted in the space of probability distributions defined over flight plans. In other words, instead of directly minimizing min_ℱJ(ℱ), the minimization of the following equivalent problem is considered:

$$\begin{array}{}\text{(14)}& \underset{p\left(\mathcal{F}\right)}{min}{\mathbb{E}}_{p\left(\mathcal{F}\right)}\left[J\right(\mathcal{F}\left)\right],\end{array}$$

which provides the same optimal solution; i.e., if ℱ^* is the optimal solution to the original problem, $\mathbf{P}(\mathcal{F}={\mathcal{F}}^{*})=\mathrm{1}$ provides the same results. The optimization is carried out in the space of probability distributions to move to continuous search space. In addition, it can facilitate searching by parameterizing the distribution p with a parameter vector θ∈ℝ^Θ, approximating Eq. (14) as

$$\begin{array}{}\text{(15)}& \underset{\mathit{\theta }}{min}{\mathbb{E}}_{p(\mathcal{F};\mathit{\theta })}\left[J\right(\mathcal{F}\left)\right],\end{array}$$

which is generally not identical to Eq. (14), as it relies on whether the parameterization is able to capture a distribution wherein $\mathbf{P}(\mathcal{F}={\mathcal{F}}^{*})=\mathrm{1}$. Thus, the parameterization of the distribution p using the vector θ plays an important role in approximating the original problem. In González-Arribas et al. (2023), the probabilistic-execution flight plan (𝒫ℱ) is introduced to parameterize the distribution over the space of possible flight plans. In this approach, the parameter vector θ is defined as follows:

$$\begin{array}{}\text{(16)}& \mathit{\theta }={\left[\begin{array}{cccccc}{\mathrm{{\rm Y}}}^T& {\widehat{\text{M}}}^T& {\widehat{\text{FL}}}^T& {\widehat{\mathcal{C}}}^T& {\widehat{\mathcal{D}}}^T& d_{\mathcal{D}}^T\end{array}\right]}^T\in {\mathbb{R}}^{\mathrm{\Theta }},\end{array}$$

where $\mathrm{{\rm Y}}\in {\mathbb{R}}^{n_{\mathrm{sp}}}$ is a vector assigning probability to select each airway, and $\widehat{\text{M}}$, $\widehat{\text{FL}}$, $\widehat{\mathcal{C}}$, $\widehat{\mathcal{D}}$, and d_𝒟 are the parameterized Mach schedule, flight level, climb profile, descent profile, and distance to go, respectively. For a parameter vector θ of a 𝒫ℱ, flight plans can be randomly generated. For instance, let us explain how the lateral path is sampled from a given vector of parameters Υ. First, assume that all waypoints of the graph are processed and limited to two outgoing airways. This can be done by considering virtual edges of zero length for waypoints with more than two outgoing airways. The vector composed of a set of junctions is defined as $\overline{V}=\mathit{\{}{v}_k{\mathit{\}}}_{k=\mathrm{1}}^{n_{\mathrm{sp}}}\in {\mathbb{R}}^{n_{\mathrm{sp}}}$. At kth junction waypoint, the selection of an airway is done through a random process: for kth entry of a given vector Υ (i.e., υ_k) and kth entry of a vector containing uniform random variables ${\mathit{\xi }}_k(\in \mathrm{\Xi }\sim U(\mathrm{0},\mathrm{1}\left)\right)$, the airway is selected as

$$\mathrm{\text{selected airway at }}k\text{th junction waypoint}=\left\{\begin{array}{ll}\text{airway 1}& S\left({\mathit{\upsilon }}_k\right)\le {\mathit{\xi }}_k\\ \text{airway 2}& S\left({\mathit{\upsilon }}_k\right)>{\mathit{\xi }}_k,\end{array}\right.$$

where S(⋅) is a sigmoid function: $S\left(x\right)=\mathrm{0.5}[\mathrm{1}+x(\sqrt{\mathrm{1}+x^{\mathrm{2}}}{)}^{-\mathrm{1}}]$. Therefore, for a given vector of parameters Υ and a given vector of random variables, a lateral path (restricted to the structure of airspace) can be sampled. The approach to sampling a complete flight plan from p(ℱ;θ), such as sampling Mach schedule and flight level, is presented in González-Arribas et al. (2023) (see Algorithm 1). It is worth mentioning that the operational aspects and feasibility of the aircraft trajectory are considered within the generation of flight plans. For instance, continuous Mach adjustment is avoided, and the frequency of flight level changes is limited.

3.4 Optimization: augmented random search

In the probabilistic-execution flight plan approach, to sample the flight plan from a given θ, vectors containing random variables are required. Thus, the 𝒫ℱ associated with θ is stochastic. To evaluate ${\mathbb{E}}_{\mathit{\theta }}\left[J\right(\mathcal{F}\left)\right]:={\mathbb{E}}_{p(\mathcal{F};\mathit{\theta })}\left[J\right(\mathcal{F}\left)\right]$ for a given θ, sampling of multiple flight plans is required. By generating N_FP sets of random variables, N_FP flight plans (i.e., FP_j for $j=\mathrm{1},\mathrm{\cdots },N_{\mathrm{FP}}$) can be sampled independently for a given θ. To benefit from the provided N_EPS probable realizations of meteorological variables and N_EPS samples of exogenous sources of uncertainty (i.e., initial flight time and initial flight mass), we sample N_EPS (i.e., N_FP=N_EPS) potential flight plans for a given θ. In other words, each sampled flight plan is evaluated with one realization of meteorological conditions. In this respect, one can rewrite Eqs. (10) and (13) as

$$[\widehat{\text{FT}},\widehat{\text{FB}},\widehat{\text{ATR}}]={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{EPS}}}}\sum _j\text{TI}({\text{FP}}_j,{\mathcal{W}}^j,t_{\mathrm{0}}^j,m_{\mathrm{0}}^j),$$

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathbb{E}}_{\mathit{\theta }}\left[J\right(\mathcal{F}\left)\right]& \approx \widehat{J}:={\mathit{\psi }}_{\mathrm{CST}}[{\mathit{\psi }}_t\cdot \widehat{\text{FT}}+{\mathit{\psi }}_m\cdot \widehat{\text{FB}}]\\ & +{\mathit{\psi }}_{\mathrm{CLM}}\cdot \widehat{\text{ATR}}.\end{array}\end{array}$$

The flight planning problem can now be expressed as the following stochastic optimization problem:

$$\begin{array}{}\text{(18)}& \underset{\mathit{\theta }}{min}\mathbb{E}\left[\widehat{J}\right],\end{array}$$

where the objective is to find the optimal value of θ (which parameterizes the probability distribution p(ℱ;θ)) such that a population of randomly sampled flight plans minimizes the expected cost in Eq. (18). Sampling N_EPS flight plans and using them in a pairwise manner with ensemble members for trajectory integration to evaluate Eq. (17) is computationally more efficient than sampling a different number of flight plans (N_FP≠N_EPS) and then evaluating each sampled flight plan with all the ensemble members. This is because, in this case, we only integrate the aircraft trajectory N_EPS times instead of performing N_FP×N_EPS trajectory integrations. In spite of reducing the number of computations, it provides similar results. This is due to the fact that, despite sampling multiple flight plans for the evaluation of the objective function in Eq. (17) for a given θ, as the process goes by, all flight plans converge to a unique flight plan. For instance, let us consider the process of sampling the lateral path presented in Sect. 3.3. The choice of each airway relies on two parameters: υ and ξ. For the first iterations of the optimization algorithm, all outgoing airways from waypoints are almost equally probable. However, with more iterations, the decision variable θ is improved, leading to increasing or decreasing the parameter υ. This parameter converges to a large positive or negative value for the last iterations. For instance, in the case of large υ, ${lim}_{\mathit{\upsilon }\to \mathrm{\infty }}S\left({\mathit{\upsilon }}_k\right)=\mathrm{1}$, implying S(υ_k)≥ξ_k for all the sampled flight plans. Thus, in the end, for the optimal value of θ obtained from optimization, we receive a unique flight plan. In this case, for the similar flight plans that are sampled for the last iterations, the expected flight performance variables obtained (from Eq. 17) using N_FP×N_EPS and N_EPS times trajectory integration will be almost the same.

In this work, the V1 version of the augmented random search (ARS) algorithm adopted from Mania et al. (2018) and González-Arribas et al. (2023) is employed. This is a gradient-like algorithm, which starts from an initial point θ₀ and then generates n random search directions ω∈ℝ^Θ. In the next step, it evaluates ${\widehat{J}}^{+}:=\widehat{J}(\mathit{\theta }+\mathbf{S}\mathit{\omega })$ and ${\widehat{J}}^{-}:=\widehat{J}(\mathit{\theta }-\mathbf{S}\mathit{\omega })$, where $\mathbf{S}\in {\mathbb{R}}^{\mathrm{\Theta }\times \mathrm{\Theta }}$ is a diagonal matrix, adjusting the relative variations that are allowed between decision variables. Then, the decision variables θ are improved along all search directions proportional to ${\widehat{J}}^{+}-{\widehat{J}}^{-}$. As can be concluded, at each iteration, we generate 2n different vectors of decision variables (i.e., n search directions with θ+Sω and θ−Sω), and for each decision variable, we sample N_EPS flight plans. Therefore, we need to perform trajectory evaluation (i.e., TI) 2n×N_EPS times for each iteration. As these calculations are similar (or very similar) and independent from each other, the parallelization on GPUs is beneficial, enabling very fast function evaluation. The algorithms, details on the optimization approach, and also parallelization on GPUs are provided in González-Arribas et al. (2023).

4 Results

The effectiveness of the proposed optimization algorithm to plan robust climate-optimal aircraft trajectories with respect to uncertain meteorological conditions is analyzed for a flight from Frankfurt to Kyiv for two different days and departure times.

• Scenario 1: 13 June 2018 at 00:00 UTC is a scenario in which the aircraft flies through areas favorable for the formation of persistent contrails.

• Scenario 2: 10 December 2018 at 12:00 UTC is a scenario with no formation of persistent contrails.

The dominant climate impact of contrails is the main reason for selecting these two scenarios, providing better insight into the mitigation potentials.

For the route graph, the full airspace graph of the considered days is filtered and processed to include all paths from the end of the standard instrument departures of the origin airport to the beginning of the standard instrument arrivals of the destination airport with the maximum length of 104 % of the shortest path length. The considered aircraft is an Airbus model A320-214, with the engine CFM56-5B4/P. Table 2 provides the required parameters to calculate the NO_x emission index of the considered aircraft using BFFM2. The initial flight time and initial mass are modeled as Gaussian variables: t₀∼𝒩(00:00 UTC, 660) [s] for scenario 1, t₀∼𝒩(12:00 UTC, 660) [s] for scenario 2, and m₀∼𝒩(61 600, 164) [kg]. The considered standard deviations for initial flight time and initial flight mass are adopted from the studies conducted by Dalmau et al. (2021) and Benjumea (2003), respectively.

As for meteorological input data, due to ease of availability, the ERA5 Reanalysis data products containing 10 ensemble members are adopted for this study. It is worth mentioning that forecast data with more ensemble members can be similarly employed. Simulations are launched on the NVIDIA GeForce RTX 3090 graphics card, providing 10496 CUDA cores at a clock speed of 1.4 to 1.7 GHz.

(DuBois and Paynter, 2006)

Table 2 The data obtained from the ICAO data bank to calculate the actual emission index of Airbus model A320-214, with the engine CFM56-5B4/P.

The weighting parameters of the objective function given in Eq. (4) are selected as ψ_CST=α [–] and ${\mathit{\psi }}_{\mathrm{CLM}}=(\mathrm{1}-\mathit{\alpha })K$ [USD K⁻¹], where K is a scaling (or conversion) factor determined as

$$\begin{array}{}\text{(19)}& K={\displaystyle \frac{{\text{SOC}}_{\mathrm{climate}}-{\text{SOC}}_{\mathrm{cost}}}{{\text{ATR}}_{\mathrm{cost}}-{\text{ATR}}_{\mathrm{climate}}}},\end{array}$$

where, for instance, SOC_climate is the SOC calculated when the optimization objective is only the climate impact or ATR_cost is the ATR when the objective is only SOC. $\mathit{\alpha }\in [\mathrm{0},\mathrm{1}]$ is a weighting parameter that penalizes cost versus climate impact in which α=0 is the pure cost-optimal and α=1 is the pure climate-optimal routing strategy. In the simulations, we consider five different values for α in order to explore the trade-off between operating cost and climate impact represented by SOC and ATR, respectively.

4.1 Scenario 1: formation of persistent contrails

We consider a scenario in which the aircraft flies through warming contrails for the cost-optimal routing option. Before presenting the results, the performance of the proposed optimizer in terms of convergence and computational time is analyzed. Since the optimization approach is stochastic, different results may be obtained with different executions. To explore the sensitivity of the optimization method, 50 different runs are performed with similar settings for pure cost- (i.e., α=1.0) and pure climate-optimal (i.e., α=0.0) routing options. Then, the objective gap is calculated considering the best performance obtained from different solutions (i.e., the minimum value of J) as the reference. The convergence performances with averaged values as solid lines, and 0, 10, 90, and 100 percentiles are depicted in Fig. 8. For both cases, with around 700 iterations (≈2.8 s) for the cost-optimal trajectory and 900 iterations (≈3.6 s) for the climate-optimal one, the estimated objective gaps are quickly reduced up to 1 % of the values of the objective functions J. As the climate-optimal routing option is associated with the inclusion of aCCFs calculated from meteorological variables, the optimization is much more complex, which can be validated in Fig. 8. With around 4000 iterations (16 s), the objective gap is reduced up to 0.5 % of J. Consequently, with a maximum of 4000 iterations, nearly optimal performance can be obtained with +0.5 % maximum deviation around the best-obtained value (i.e., the most optimal case).

Figure 8 Sensitivity of the convergence performance of the optimization algorithm to 50 different executions for the cost- and climate-optimal routing options (α=1 and α=0, respectively) (one iteration ≈4 ms). The objective gap is calculated considering the deviations from the best performance obtained among 50 runs (i.e., the minimum value of J) as the reference, i.e., at each iteration, the value of the objective function (used for trajectory optimization in Eq. 17) is compared with the minimum value obtained from other executions. The mean value is highlighted in as a solid black line, and 0, 10, 90, and 100 percentiles are represented with different color bands.

Now, we proceed to present the obtained results. The aircraft profiles and climate responses for different routing options are given in Fig. 9. The SOC depends on the flight time and fuel consumption. Therefore, the aircraft for routing strategies with higher values of α, such as $\mathit{\alpha }=\mathrm{1.0},\mathrm{0.8}$, tends to fly at higher altitudes within the vertical constraints because flying at higher altitudes is beneficial for reducing fuel consumption, which contributes a large part of the total operating cost (see Fig. 9a). By analyzing the lateral paths depicted in Fig. 10 with the direction and speed of wind at different flight levels, one can see that the aircraft deviates from the shortest path to benefit from stronger tail winds. For trajectories with lower climate impacts, as can be seen in Fig. 9a, the aircraft flies at relatively lower altitudes compared to the cost-optimal routing option, mainly to avoid the formation of persistent contrails (due to warming impacts during nighttime). The climate-optimal routing options reduce the warming effect of contrails. Although the warming climate effects of NO_x emissions and water vapor increase with climate-optimal trajectories, the net climate impact decreases. This is because the climate impact of contrails outweighs the impact associated with other species (as discussed in Sect. 2.2.2). The contribution of each species to total climate impact, variability of the obtained climate impacts, and SOC with ranges of uncertainty for different α values and Pareto frontiers are provided in Fig. 11. For a specific case (α=0.2), by accepting an increase of 4.3 % in cost, there is a potential to mitigate the climate impact by 53.0 % considering median performance. In Sect. 2.3, it was shown that the variability of relative humidity among ensemble members is high, leading to high uncertainty in the aCCF of contrails. As expected, the obtained climate impact of contrails is highly uncertain when the aircraft flies through areas sensitive to forming persistent contrails. In contrast, as the aircraft tends to avoid PCFAs, the ranges of uncertainty are reduced, in which, for the complete avoidance that is achieved with α=0.2, the climate impact is almost deterministic. In addition, SOC requires flight time and fuel consumption to represent operating cost in USD, and as it is affected by relatively less uncertain meteorological variables (i.e., wind and temperature compared to relative humidity for the considered case study as analyzed in Sect. 2.3), the uncertainty in its value is small.

Figure 9 Results of scenario 1 (13 June 2018, 00:00 UTC) for different routing options (i.e., α). The shaded regions show the ranges of uncertainty associated with uncertain meteorological conditions (outer lighter areas show the minimum and maximum values, while the inner darker ones represent the 95 % confidence interval).

By analyzing the contribution of each species to the net ATR for different α, one can conclude that the mitigation potential is achieved mainly by avoiding contrail-sensitive areas, which results in slight increases in the climate impact of NO_x emissions (see Fig. 11a). However, when the formation of persistent contrails is completely avoided (i.e., with α=0.2), the optimizer tends to reduce NO_x emissions mainly by reducing speed to reduce the fuel flow required to calculate the NO_x emission index and also total NO_x emissions (i.e., NO_x emissions = NO_x emission index ⋅ fuel consumption). Reducing NO_x emissions by flying at lower speeds is achieved at the expense of a considerable increase in flight time and, consequently, SOC. As can be concluded from Pareto frontiers, such a reduction in climate impact for this scenario is not cheap as only 5 % more reduction in climate impact is obtained with almost 3 % more increase in SOC (α=0.0). As the aCCF of contrails is only evaluated in areas favorable for the formation of persistent contrails, typically determined using inequality constraints, it has sharp spatial behaviors (i.e., PCFA (latitude, longitude, altitude, time) ∈ {0, 1}). In addition, contrails have dominant climate effects. Therefore, the optimizer's first choice is to avoid forming persistent contrails, which may be achieved more efficiently than reducing the impacts of other species with relatively lower climate impacts and smooth spatial behaviors. This can be validated in Fig. 11a, as the lowest priority is given to reducing the climate impact associated with NO_x emissions (for α=0).

Figure 10 Lateral paths for scenario 1 (13 June 2018, 00:00 UTC) depicted with (a) wind and (b) aCCF of contrails as color maps.

Figure 11 Overall performance of the optimized trajectories in terms of ATR and SOC for scenario 1 (13 June 2018, 00:00 UTC).

4.2 Scenario 2: no formation of persistent contrails

In the next scenario, we analyze the mitigation potential when no persistent contrails are formed with the cost-optimal routing option.

For this case, aircraft profiles and climate responses are depicted in Fig. 12a and b, respectively. As can be seen in Fig. 12a, the optimizer chooses to fly at lower altitudes for routing strategies with higher penalization on climate impact. As no persistent contrails are formed (see Fig. 13a), we depict the lateral paths with the merged aCCFs (calculated using the mean values of the obtained NO_x emission index) as the color map at different flight levels in Fig. 13b. As can be seen, flying at lower altitudes is more beneficial in reducing the climate impact of other species (mainly NO_x). In addition to lowering cruise altitude, the aircraft flies at lower speeds to reduce the fuel flow and, consequently, fuel consumption, NO_x emission index, and NO_x emissions. The variability of climate impact and SOC for different α values and Pareto frontiers are given in Fig. 14. By reducing α, the climate impact decreases at the cost of an increase in SOC. For instance, for α=0.2, by accepting a 1.1 % increase in cost, a 16 % reduction in ATR can be achieved. As in the previous case, the relative increase in SOC is considerable for α=0, since the aircraft tends to fly at a relatively lower speed for more reduction in climate impact.

In conclusion, climate impact reduction is achieved at the expense of a higher cost increase than in the previous scenario. Moreover, since no contrails are formed, the uncertainty in climate impact is almost negligible.

Figure 12 Results of scenario 2 (10 December 2018, 12:00 UTC) for different routing options (i.e., α). The shaded regions show the ranges of uncertainty associated with uncertain meteorological conditions (outer lighter areas show the minimum and maximum values, while the inner darker ones represent the 95 % confidence interval).

Figure 13 Lateral paths with (a) aCCF of contrails and (b) merged aCCFs as color maps at different flight levels for scenario 2 (10 December 2018, 12:00 UTC).

Figure 14 Overall performance of the optimized trajectories in terms of ATR and SOC for scenario 2 (10 December 2018, 12:00 UTC).

5 Discussion of results and conclusion

This paper presented a methodology to plan a robust climate-optimal aircraft trajectory under uncertain meteorological conditions. The climate-sensitive regions were identified using the prototype algorithmic climate change functions (version 1.1). The ensemble prediction system was employed to characterize uncertainty in weather forecasts. A heuristic algorithm was employed and implemented on graphics processing units to solve the proposed robust trajectory optimization in a computationally fast manner. The effectiveness of the proposed approach was explored in two scenarios. Discussion of the obtained results (mainly related to the achieved mitigation potentials, current limitations, and future lines of research) and some general remarks are presented in the following.

The obtained mitigation potentials for the considered scenarios were different due to the variability of meteorological conditions. In both cases, the climate-optimal routing options could reduce the climate impacts. The cost-optimal trajectories flew at higher altitudes compared to climate-optimal ones, as flying at higher altitudes is beneficial for reducing fuel consumption. This is also in line with related studies in the literature (e.g., Yamashita et al., 2020). Due to the dominant climate impact and non-smooth spatial behavior of contrails, the mitigation potential obtained for the scenario with contrail effects (i.e., scenario 1) was higher than the scenario with no formation of persistent contrails. The non-smooth spatial behavior of the contrail climate impact is related to the conditions of PCFA. Due to the high variability among the ensemble members of relative humidity over ice needed to determine PCFA, the climate impact of contrails was highly uncertain. However, for the cases without contrail formation, the total climate impact was almost deterministic. This is because the variability in the other weather variables was almost negligible. For both scenarios, α=0.2 seems to be a reasonable choice, since the climate impacts were reduced at the expense of acceptable increases in the operating cost, and the results were almost deterministic.

In spite of considering the ensemble members in trajectory planning, a unique (or deterministic) flight plan is determined. This reflects the operational feasibility and applicability of this method since, in the flight planning context, the requirement is to determine a unique lateral route in latitude and longitude that starts and ends at predefined points in space and follows the real structure of airspace as well as having a fixed altitude profile and a fixed airspeed schedule. In this case, the effects of uncertainty are reflected in the aircraft performance variables. For instance, let us consider the first scenario. For α=1.0, the lateral path, speed schedule, and flight level were determined in a cost-optimal manner (see Figs. 9a, 10). In our approach, the optimized trajectory is assumed to be tracked as close as possible in practice with the system's low-level controllers in real time. Aiming to optimize a unique flight plan, the proposed method provides some ranges for the aircraft performance variables, such as fuel consumption, flight time, and climate impact, due to different probable realizations of weather conditions. For instance, the climate impact is expected to lie within the determined ranges (see Fig. 9b) if the considered ensemble members could acceptably predict future weather conditions.

One of the next steps should be analyzing the feasibility of such a routing strategy for real traffic scenarios. In fact, air traffic management (ATM) is a complex multi-agent system that cannot be represented by individual elements but by their collective behavior at the network scale. It was shown in the paper that for the climate-optimal routing options, the aircraft tends to fly at relatively lower altitudes compared to the cost-optimal one. Such behavior to avoid climate-sensitive areas may result in more congested areas, raising some challenges, particularly increased workload, complexity, and conflicts. Thus, the mitigation potentials reported at the micro-level may not be achievable considering real traffic scenarios. Therefore, after generating climatically optimal flight plans, one needs to assess the resulting effects at the network scale and perform a resolution (typically modeled as an optimization problem) to re-stabilize the ATM system by compensating for the negative impacts while keeping the modified trajectories as close as possible to inputted climate-optimized ones. The assessment of manageability of climate-optimal trajectory planning in an ATM system is called macro-scale analysis and lies outside the scope of this paper (see Baneshi et al., 2023, for a study in this area).

In this study, we only considered the minimization of the expected performance, e.g., expected climate impact. However, the concept of robustness is mainly related to having less uncertain results (i.e., also minimizing the uncertainty range). In the case of robustness to meteorological uncertainty, we need to find a flight plan that avoids areas of airspace with high variability among the ensemble members. For instance, in González-Arribas et al. (2018), the dispersion of the arrival time is minimized by avoiding regions with high variability among the ensemble members of wind (characterized by SD). In this study, we observed that the most uncertain variable is the climate impact of contrails. Since, for both scenarios, only warming contrails were formed, minimization of the expected values directly led to avoiding uncertain PCFAs, and as can be seen from the results, for the climate-optimized trajectories, we obtained robust solutions (almost deterministic results). However, during the daytime, different behavior is expected for the cases with the cooling contrails. This is because, to reduce the expected climate impact, aircraft will tend to fly in uncertain PCFAs to benefit from cooling effects. In the next versions of ROOST, such scenarios will be taken into account, and controlling the dispersion of all flight performance variables will be addressed by including their variance in addition to the averaged values as objectives in the objective function.

Regarding the computational time and convergence performance, it was shown that they are scenario-dependent. For more complex problems, such as the case including climate impacts quantified by using aCCFs, the optimizer required more iterations to enhance the convergence compared to the cost-optimal routing option. It is worth mentioning that the distance between the origin and destination, available route graphs, and also parameters within the optimization algorithm can change the convergence performance and computational time. The number of iterations is a user-defined parameter that needs to be specified based on the required performance and availability of computational resources. In the performed simulations, we considered 4000 iterations. By looking at the Pareto frontiers, it is clear that the optimizer was able to find nearly optimal solutions. Thanks to the parallelization on GPUs, the computational time for achieving a nearly optimal solution is promising. There are several controlling parameters within the optimization algorithm of ROOST, including the number of search directions, the augmented random search (ARS) step size, and the Nesterov velocity factor (see González-Arribas et al., 2023, for a description of these controlling parameters) that can affect the convergence performance. In our future work, we aim to examine the effects of all these parameters and propose an optimal selection of them. Moreover, adaptive (scenario-dependent) stopping criteria will be proposed, helping to optimize aircraft trajectories more efficiently in the sense of computational time.

To explore the trade-off between climate impact and the operating cost, Pareto frontiers were generated. By changing the weighting parameter α, different Pareto-optimal solutions were obtained. However, with this approach, a specific value for α does not necessarily result in a similar cost increase and climate impact mitigation potential for different scenarios (e.g., for α=0.2, the climate impact is reduced by 53 % and 16 % by accepting a 4.3 % and 1.1 % increase in the operating cost for scenarios 1 and 2, respectively). This approach is suitable for analyzing the mitigation potential for a single flight. However, it is not the most efficient way to study the Pareto-optimal solutions of the aggregated results of a large number of flights. For such cases, having the flexibility to directly optimize flights, requesting a certain range reduction in climate impact or allowing a specific range for the increased operating costs would be beneficial. In this respect, one can aggregate climate-optimal trajectories having, for instance, a 0.5 % to 1.0 % increase in the operating cost. In future versions of ROOST, we will add this feature to the optimization tool by defining some path constraints. Such an aggregation of results is doable with α; however, one needs to generate more points in the Pareto frontier in order to classify the results based on a certain percentage increase in cost or a percentage decrease in climate impact. Scaling up this trajectory-level analysis to the network scale will increase the computational time by a factor of the considered α values.

As was explored in the paper, the mitigation of climate impact within the flight planning context is achieved only by accepting some extra costs due to avoidance of highly climate-sensitive regions, which is in line with related studies in the literature (e.g., Yamashita et al., 2020, 2021; Lührs et al., 2021; Niklaß et al., 2019). One potential approach to offsetting these additional costs is implementing fees and taxes that account for the climate impact caused by aviation in order to motivate airliners to adopt climate-optimal trajectories. Currently, there is no climate policy for the aviation-induced non-CO₂ climate effects in the planned market-based instruments such as the emission trading scheme (ETS). However, several studies have proposed frameworks to include costs of aviation-induced non-CO₂ climate effects in the operating cost, such as the concept of equivalent CO₂ emissions and the concept of climate-charged airspace (Niklaß et al., 2021). If taxes are applied to the non-CO₂ climate effects, climate-optimal trajectories can also be economically beneficial.

The aCCFs used in this study represent a prototype formulation. The aCCF algorithms were developed for meteorological summer and winter conditions with a focus on the North Atlantic flight corridor. Thus, the usage of the aCCFs for different seasons and regions needs special caution. However, further development of the aCCFs and an expansion of their geographic scope and seasonal representation represent ongoing research.