Fundamental Framework to Plan 4D Robust Descent Trajectories for Uncertainties in Weather Prediction

Abstract

Aircraft trajectory planning is affected by various uncertainties. Among them, those in weather prediction have a large impact on the aircraft dynamics. Trajectory planning that assumes a deterministic weather scenario can cause significant performance degradation and constraint violation if the actual weather conditions are significantly different from the assumed ones. The present study proposes a fundamental framework to plan four-dimensional optimal descent trajectories that are robust against uncertainties in weather-prediction data. To model the nature of the uncertainties, we utilize the Global Ensemble Forecast System, which provides a set of weather scenarios, also referred to as members. A robust trajectory planning problem is constructed based on the robust optimal control theory, which simultaneously considers a set of trajectories for each of the weather scenarios while minimizing the expected value of the overall operational costs. We validate the proposed planning algorithm with a numerical simulation, assuming an arrival route to Leipzig/Halle Airport in Germany. Comparison between the robust and the inappropriately-controlled trajectories shows the proposed robust planning strategy can prevent deteriorated costs and infeasible trajectories that violate operational constraints. The simulation results also confirm that the planning can deal with a wide range of cost-index and required-time-of-arrival settings, which help the operators to determine the best values for these parameters. The framework we propose is in a generic form, and therefore it can be applied to a wide range of scenario settings.

1. Introduction

Aircraft trajectory planning, theoretically leading to a trajectory optimization problem, requires a variety of models and input data to predict a trajectory. They include the aircraft performance models, e.g., the equations of motion (EoM), meteorological models and predicted weather data. Reliable flight planning is based on accurate predictions and calculations of aircraft trajectories. However, those models and input data generally contain errors, and they usually bring uncertainties to trajectory planning. One of the most serious uncertainties comes from weather prediction. Conventional deterministic planning algorithms only assume a fixed future weather situation based on an available weather forecast and optimize the trajectory for the assumed scenario. However, the weather situation can sometimes change rapidly within a very short time. If the assumed future weather scenario is greatly different from the true weather situation that the aircraft actually experiences, the objectives for the optimization can be deteriorated and the planned trajectory is no longer the optimal trajectory. In addition to that, such a trajectory plan can also be infeasible in the sense that it violates operational limitations, such as speed limits, which is critical for flight safety. To this end, the present study proposes a framework for planning four-dimensional (4D) optimal descent trajectories that are robust against uncertainties in weather prediction.

The descent phase of flight, together with continuous descent operations (CDO), holds potential to make air traffic more efficient and environmentally friendly [1]. Nowadays, societies are concerned about the impacts of aircraft operations on the climate. To deal with this concern, worldwide collaborative projects are being carried out to renovate the current air traffic management (ATM) systems, aiming at reducing negative environmental impacts. In particular, reductions in fuel use, noise and gaseous emissions are the key performance indicators in Europe’s vision for aviation, Flightpath 2050 [2]; and the related research programs by Single European Sky ATM Research (SESAR) (in Europe) [3], Collaborative Actions for Renovation of Air Traffic Systems (CARATS) in Japan [4] and Next Generation Air Transportation System (NextGen) in the United States [5]. In Europe, the European Union (EU) formulated the Single European Sky (SES) performance scheme to comply with the above expectations [6]. In these long-term plans, continuous climb/descent operations (CCO/CDO) are considered to be the crucial components contributing to these goals: one goal is enabling low to idle thrust settings by having continuous descents, and avoiding level flight segments especially at low altitudes [1,7,8]. In 2018, the European Organisation for the Safety of Air Navigation (Eurocontrol) assessed the benefit pool of CDO for Europe and reported it in [9]. According to the report, traffic data of 2017 showed that the average time in level flight from top of descent (ToD) performed by non-CDO flights was 217 s, and the per-flight fuel savings by potential CDO implementation was estimated to be 46 kg, equivalent to 145 kg of CO${}_2$ emissions or 20 € of fuel costs. The report also estimated that CDO can reduce noise emissions by 1–5 dB per flight. Considering these benefit estimations, CDO and the optimized profile descents (OPDs), the corresponding concept of NextGen, are listed in the road maps as one of the key measures to allow the ATM systems to become more efficient and eco-friendly [3,4,5]. These more recent types of descent have been implemented in several regions throughout the world following these plans [10,11,12].

Trajectory optimization contributes to facilitating the CDO and maximizing its potential benefits. Studies on simulating and optimizing descent trajectories, which eventually lead to CDO trajectories, have been carried out intensively in the last decade [11,13,14]. The necessity of re-optimizing trajectories has also been raised in [15,16,17,18,19] by examining the deviation of an aircraft’s energy state from the planned reference trajectory as a useful trigger. However, to the best of the authors’ knowledge, these prior works assumed deterministic scenarios and ignored the impacts of uncertainties. Robust aircraft trajectory planning has recently become a field for intensive research, and several applications for ATM have been proposed. Franco et al. proposed a probabilistic trajectory predictor based on the probabilistic transformation method to transform uncertainties from ground speed to estimated flight time and fuel consumption [20]. The method was extended by Hernández-Romero, who developed methods to solve the probabilistic aircraft conflict detection and resolution (CD&R) problem [21]. Franco et al. also applied mixed-integer linear programming to analyze arrival time uncertainty [22]. Legrand et al. solved a robust planning problem with dynamic programming [23].

Another major way to formalize robust trajectory planning is to utilize robust optimal control theory. In prior publications, terminology was not fixed. The relevant methodologies were referred to as either “robust”, “stochastic” or “probabilistic” optimal control, depending on the authors. This class of methods transforms a stochastic optimal control problem to an equivalent deterministic problem. The transformation is typically carried out by either applying types of the polynomial chaos (PC) expansion or utilizing a set of scenarios representing the considered uncertainties (e.g., a dataset from ensemble weather forecasing). The theoretical aspects of this approach have been intensively studied in recent years from the control-theory perspective. Fisher et al. applied PC to solutions for constrained mechanical systems, e.g., trajectory generation, and examined its performance [24]. Cottrill proposed a solution algorithm for stochastic optimal control problems, combining the Gauss pseudospectral method (GPM) and the generalized PC (gPC) [25]. Applications of these methods to aircraft trajectory optimization were also studied by Li et al. [26] and Piprek [27]. However, there is a gap between these control-theory-oriented research works and their applications to practical operational problems of ATM, as these works mainly focused on theoretical properties of the method, and did not consider detailed models and constraints that characterize realistic ATM scenarios. Several works have been reported to address ATM applications. González-Arribas et al. applied robust optimal control for cruise flight with an uncertain wind condition by employing ensemble prediction systems (EPS) [28] and confirmed the method could reduce the impact of uncertain wind on the planning, taking as reference one year of trajectory data and multiple origin–destination pairs [29]. This robust planning method was extended in [30] to the consideration of exposure to convection and cost-index based profiles. Recently, the same authors proposed a heuristic method based on parallel graphics processing unit (GPU) computation [31], finding robust trajectories in computational times that are compatible with real operations (∼1 s), and considering not only the level flight, but step-climbs and step-descents. In addition to these works, Matsuno et al. applied stochastic optimal control to CD&R problems and solved them with gPC and GPM [32]. These prior works have revealed the importance of robust trajectory planning; however, they are limited to the en-route phase of flight. Robust trajectory planning has the potential to contribute to the descent phase as well. For example, it can improve the predictability of CDO, as robust planning is able to manage uncertainties in the aircraft’s position at certain times. On the contrary, the descent phase typically requires the optimizer to deal with more operational constraints than the cruise phase (e.g., speed/altitude constraints and waypoint location), which may limit the room for finding optimal trajectories. Therefore, research should also be carried out for descent trajectory to comprehend the benefits and limitations of robust planning. In the authors’ previous publication, we applied robust optimal control to CDO trajectory planning and evaluated its performance with various parameter settings [33], considering trajectories in a vertical plane.

The present study aims to allow robust trajectory planning to be utilized for more operation-oriented analysis. To this end, it extends the methodologies published in [33] so that they can be applied to 4D trajectory planning. The ability to optimizing 4D trajectories makes it possible to consider existing horizontal arrival route structures in planning. Unlike the prior studies mentioned above, we formalize the planning problem as a multiple-phase optimal control problem to include waypoints and speed/altitude limitations prescribed for the arrival route. We assume the algorithm is used for in-flight trajectory planning on a flight deck and consider a pre-tactical planning horizon of at least 1–2 h prior to reaching ToD. It is worth noting that unconsidered uncertainties can make the aircraft deviate from the planned trajectory at execution. It is out of the scope of this study to investigate the impact of adhering to the planned trajectory or tactically re-optimizing the trajectory during execution. Studies on these topics for deterministic planning can be found in [15,16].

2. Models for Robust Trajectory Planning

This section discusses the models for trajectory planning. The International System of Units’ units (SI units) are used, unless explicitly stated differently.

2.1. Uncertainty Models of Weather Prediction

The present study deals with uncertainties in weather prediction, which cause errors between the forecast and the true weather that the aircraft actually experiences during its flight execution. In order to model those uncertainties, we utilize the weather forecast data from the so-called Global Ensemble Forecast System (GEFS). It is a set of weather forecast scenarios provided by the National Oceanic and Atmospheric Administration (NOAA) [34] and is available online. We use version 11.0 for the present study (the newer version, 12.0, was released recently with 32 ensemble members). GEFS consists of 21 separate weather forecasts called weather ensemble members (Figure 1). Each member is obtained by perturbing the initial conditions of the weather dynamics to represent uncertainties in the meteorological models used for the forecast [35]. Small perturbations in the initial conditions evolve and lead to different weather states. The forecasts are updated every 6 h. A detailed investigation of GEFS and its application to cruise-trajectory planning can be found in [17]. The data are provided in a GRIdded Binary (GRIB2) format, which involves three dimensional grids specified by latitude, longitude and pressure level. The grid resolution is 1${}^{\circ }$ for the latitude and the longitude, and 50 hPa for the pressure level.

We utilize the GEFS data to obtain the ambient temperature $T_{mp}$, the pressure p and the wind speed vector, which holds the horizontal components u and v. The wind-speed components are positive for the west-to-east wind flow and the south-to-north flow, respectively. GEFS does not contain the vertical wind component. When creating the weather models, we make the following assumptions:

Assumption 1

(Horizontal wind field). We assume the wind field only has horizontal components and ignore the vertical one. This assumption is justified as the vertical component is usually one order of magnitude lesser than the horizontal components [36].

Assumption 2

(Stationary atmosphere). We assume the atmosphere is stationary; i.e., its time evolution is not considered. Therefore, the atmospheric quantities, i.e., the wind components, the temperature, the pressure and the density, are treated as functions of the aircraft 3D position.

Under these assumptions, the GRIB2-style data of GEFS are symbolically expressed as:

(1)$u_G=u_G\left({\varphi }_G,{\lambda }_G,p_G\right),v_G=v_G\left({\varphi }_G,{\lambda }_G,p_G\right),T_{m{p}_G}=T_{m{p}_G}\left({\varphi }_G,{\lambda }_G,p_G\right),$

where $\varphi$ and $\lambda$ are the latitude and longitude in degrees, respectively; and the subscript ${}_G$ denotes the quantities are for a grid point. We interpolate these discrete data so that the weather model will be available at any aircraft position, or any combination of ($\varphi$, $\lambda$, h). We use the three dimensional b-spline function denoted with Equation (A1) for interpolation (see Appendix A).

One remark for integrating these weather models into trajectory planning is that the aircraft dynamics normally use the geometric altitude h instead of the pressure level p to describe the height (see Section 2.2). Therefore, we need a conversion law between the two quantities. To this end, we use the conversion formula that NOAA officially provides:

(2)$p\left(h\right)=1013.25\times {\left(1-\frac{0.0065}{288.15}h\right)}^{5.2561}.$

It is worth noting that this conversion is equivalent to assuming the International Standard Atmosphere (ISA) for the air pressure [37]. If we apply Equation (2) to substitute p in Equation (A1), the b-spline weather model for the w-th member of the weather ensemble can be rewritten as functions of the aircraft 3D position: $f^{bs\left(w\right)}(\varphi ,\lambda ,p\left(h\right))$. We apply this b-spline interpolation to each member of GEFS to obtain member-specific weather models (Figure 1). The b-spline-based weather models for the weather scenario w are symbolically expressed as:

(3)$\begin{array}{ccc}\hfill u^{\left(w\right)}\left(\varphi ,\lambda ,h\right)& =& u^{bs\left(w\right)}\left(\varphi ,\lambda ,p\left(h\right)\right)\hfill \\ \hfill v^{\left(w\right)}\left(\varphi ,\lambda ,h\right)& =& v^{bs\left(w\right)}\left(\varphi ,\lambda ,p\left(h\right)\right)\hfill \\ \hfill T_{mp}^{\left(w\right)}\left(\varphi ,\lambda ,h\right)& =& T_{mp}^{bs\left(w\right)}\left(\varphi ,\lambda ,p\left(h\right)\right)\mathrm{for}w=1,\cdots ,N_w.\hfill \end{array}$

The air density $\rho$, also required for trajectory calculation, can be obtained by the gas equation

(4)${\rho }^{\left(w\right)}\left(\varphi ,\lambda ,h\right)=\frac{p\left(h\right)}{R{T}_{mp}^{\left(w\right)}\left(\varphi ,\lambda ,h\right)},$

where $R=287.058$ J kg${}^{-1}$ K${}^{-1}$ is the gas constant of dry air.

2.2. Flight Performance Models

We consider an aircraft as a dynamic system whose states are governed by the point-mass EoMs and are controlled through control inputs. When deriving the EoMs, we make the following assumptions.

Assumption 3

(Geodetic coordinate system). We express the aircraft position in the geodetic coordinates, where the aircraft position is specified with the angular coordinates, i.e., ϕ and λ, and the vertical position parameter, i.e., h. (Three different coordinate systems are usually used as reference frames fixed to the elliptical earth. They include the geodetic coordinates, the earth-centered earth-fixed Cartesian coordinates and the earth-surface Cartesian coordinates [38]. The third system is often referred to as the north–east–down (NED) coordinates, since normally its x-axis points to the north, its y-axis to the east and its z-axis downward.) This coordinate system is selected to align the aircraft dynamics with the weather model (Section 2.1). The components of wind u and v are positive for the west-to-east wind flow and the south-to-north flow, respectively.

Assumption 4

(WGS-84 Earth model). We utilize the WGS-84 Earth model defined by National Imagery and Mapping Agency (NIMA) [39]. It is an ellipsoidal earth model and is widely used for Flight Management System (FMS) [40]. The derivation of the aircraft EoM for this earth model is beyond the theme of this publication (see [38] for details).

Assumption 5

(Engine thrust is controlled by the thrust coefficient). We employ the Base of Aircraft Data family 4 (BADA4) engine thrust model [41], where the thrust T is parametrized with the thrust coefficient $C_T$ and p:

(5)$T=T\left(C_T,p\right).$

We treat $C_T$ as a control variable representing thrust control and limit it to the following range:

(6)$C_{T_{idle}}(M,p)\le C_T\le C_{T_{MCRZ}}(M,p),$

where $C_{T_{idle}}$ and $C_{T_{MCRZ}}$ are the thrust coefficients for the idle and the maximum-cruise engine ratings for turbofan engines. They are functions of Mach number M and p.

Under these hypotheses, the aircraft’s 4D dynamics are expressed as:

(7)$\frac{d}{dt}{\left\{V_{TAS},\gamma ,\psi ,\varphi ,\lambda ,h,m\right\}}^T={\left\{f_V,f_{\gamma },f_{\psi },f_{\varphi },f_{\lambda },f_h,f_m\right\}}^T,$

where

$f_V=\frac{T-D}{m}-gsin\gamma ,f_{\gamma }=\frac{Lcos{\varphi }_b-mgcos\gamma }{m{V}_{TAS}},f_{\psi }=\frac{Lsin{\varphi }_b}{m{V}_{TAS}cos\gamma },$

$f_{\varphi }=\frac{V_{TAS}cos\gamma cos\psi +v}{R_{\varphi }+h},f_{\lambda }=\frac{V_{TAS}cos\gamma sin\psi +u}{(D_{\varphi }+h)cos\varphi },f_h=V_{TAS}sin\gamma ,$

(8)$f_m=-FC.$

The aircraft state variables are $\mathit{x}={\left\{V_{TAS},\gamma ,\psi ,\varphi ,\lambda ,h,m\right\}}^T$, where $V_{TAS}$ is the true air speed (TAS), $\gamma$ is the flight path angle, $\psi$ is the heading and m is the gross mass. The control inputs, needed to close the set of differential equations, are $\mathit{u}={\left\{C_T,{\varphi }_b,C_L,{\delta }_{SB}\right\}}^T$, where ${\varphi }_b$ is the bank angle, $C_L$ is the lift coefficient and ${\delta }_{SB}$ is the speed brake control. $R_{\varphi }$ stands for the radius of curvature of the meridian’s ellipse [38] and $D_{\varphi }$ for the radius of curvature in the prime vertical [39] of the WGS-84 Earth model given as a function of $\varphi$:

(9)$R_{\varphi }\left(\varphi \right)=\frac{\mathit{a}(1-{\mathit{e}}^2)}{{(1-{\mathit{e}}^2{sin}^2\varphi )}^{3/2}},D_{\varphi }\left(\varphi \right)=\frac{\mathit{a}}{{(1-{\mathit{e}}^2{sin}^2\varphi )}^{1/2}},$

where $\mathit{a}=6.378\times {10}^6$ m is the equatorial radius and $\mathit{e}=8.18\times {10}^{-2}$ is the eccentricity. D is the drag force and is written as:

(10)$D=\frac{1}{2}\rho V_{TAS}^{2}S\left(C_{D_{\mathrm{clean}}}+C_{D_{SB}}{\delta }_{SB}\right),$

where S is the area of the wing, $C_{D_{\mathrm{clean}}}$ is the drag coefficient in the clean configuration and $C_{D_{SB}}$ is the coefficient of the speed brake effect. The two coefficients are taken from BADA4 [41]. The second term on the right-hand side of Equation (10) represents the effect of speed brake control on drag. We neglect the effect of flap and slat controls. $FC$ represents the fuel consumption and we used the BADA4 model: $FC=FC(V_{TAS},C_T,p,T_{mp})$.

To simplify the aircraft dynamics and make the problem easier and faster to solve, we make additional assumptions.

Assumption 6

(Heading as a control input). We assume ψ can be changed instantaneously instead of following the rotational dynamics (Equations (7) and (8)). We neglect the corresponding dynamics and treat ψ as a control variable, since it still appears in other dynamics (Equation 8). This also allows us to ignore ${\varphi }_b$ and remove it from the control inputs.

Assumption 7

(Vertical equilibrium). We assume the vertical forces acting on the aircraft are balanced:

(11)$L=mgcos\gamma .$

This assumption leads to $d\gamma /dt=0$ in Equation (7) and allows us to neglect the γ-dynamics. γ is excluded from the state variables and we include it in the control inputs instead. This type of model is called the γ-command model and is widely used for trajectory calculation for ATM purposes [13,15,42,43,44]. It also allows us to remove $C_L$ from the control inputs. This is because $C_L$ can be calculated automatically if γ is specified in Equation (11) as:

(12)$C_L=\frac{2L}{\rho V_{TAS}^{2}S}=\frac{2mgcos\gamma }{\rho V_{TAS}^{2}S}.$

Under these assumptions, the EoMs are simplified as

(13)$\frac{d}{dt}{\left\{V_{TAS},\varphi ,\lambda ,h,m\right\}}^T={\left\{f_V,f_{\varphi },f_{\lambda },f_h,f_m\right\}}^T,$

where the aircraft state variables are $\mathit{x}={\{V_{TAS},\varphi ,\lambda ,h,m\}}^T$ and the control variables are $\mathit{u}={\{C_T,\gamma ,\psi ,{\delta }_{SB}\}}^T$.

Currently, we have the time t as the independent variable for our formalization. From the optimal control viewpoint, this means t is required to be in the same domain $[t_0,t_f]$ for any of the state variables. In other words, the initial time $t_0$ and the final time $t_f$ must be the same for the trajectories of any state and control variables. As we discuss in Section 3.2, we eventually aggregate the EoM to cover the dynamics for all the weather scenarios to create a robust optimal control problem. t being in the same domain means the arrival time at the final position must be common for any member of the trajectory ensemble. However, this is normally not true. A different weather situation, especially a different wind speed, causes different ground speed, and consequently the arrival time should be varied. Therefore, we need a different independent variable to properly model the variations in trajectory with respect to time. We make the following assumption:

Assumption 8

(Along-track distance as the independent variable). Instead of t, we treat the along-track distance s, another strictly increasing variable, as the independent variable. As seen in Figure 2, s is defined as the distance along the aircraft path projected onto the earth surface and expressed using the ground speed $V_{GS}$:

(14)$s={\int }_{t_0}^{t_f}{V}_{GS}dt,$

where $V_{GS}$ is expressed as:

(15) $V_{GS}=\sqrt{{(V_{TAS}cos\gamma cos\psi +v)}^2+{(V_{TAS}cos\gamma sin\psi +u)}^2}.$

The feature of s being in the same domain $[s_0,s_f]$ for any member of trajectory ensemble is justified by the commonality constraints we will introduce in Section 3.2, which require the aircraft to fly the common 3D route regardless of the weather scenario. The selection of s as the alternative independent variable can also be seen in prior studies [28,43,44]. We perform this transformation, applying the chain rule $d/dt=d/ds·ds/dt$. Here, $ds/dt=V_{GS}$ due to Equation (14), and therefore the time derivative can be written as $d/dt=V_{GS}·d/ds$. Applying this formula to Equation (13), we can change the independent variable from t to s. This conversion removes t from the EoMs. Instead, we treat it as an additional state variable and include its evolution in the set of EoMs. The resulting EoMs are:

(16)$\frac{d}{ds}{\left\{V_{TAS},\varphi ,\lambda ,h,m,t\right\}}^T=\frac{1}{V_{GS}}{\left\{f_V,f_{\varphi },f_{\lambda },f_h,f_m,1\right\}}^T,$

where the independent variable is s, the state variables are $\mathit{x}={\{V_{TAS},\varphi ,\lambda ,h,m,t\}}^T$ and the control variables are $\mathit{u}={\{C_T,\gamma ,\psi ,{\delta }_{SB}\}}^T$. We use these EoMs for the robust trajectory planning. For facilitating the following discussion, we define the aircraft dynamics vector as $\mathit{f}={\{f_V,f_{\varphi },f_{\lambda },f_h,f_m,1\}}^T/V_{GS}$.

3. Robust Descent Trajectory Planning

3.1. Deterministic Formalization as a Basis

Before formalizing a robust trajectory planning problem, we first create the corresponding deterministic problem to be a basis of the robust formalization. We split the trajectory into phases so that we can consider phase-specific constraints, and we express the total number of phases with $N_p$. An overview of the phase-specific constraints and symbols is summarized in Figure 3. To deal with these phases, we apply the multiple-phase optimal control theory to the planning problem [33,45]. The variables in the optimal control problem are the state vector $\mathit{x}$ and the control vector $\mathit{u}$, and the independent variable s. We add the superscript ${}^{<p>}$ to $\mathit{x}$ or $\mathit{u}$ to denote the variable is for the phase p ($=1,\cdots ,N_p$). $s_0$ and $s_f$ are used to represent the s of the initial and final positions, respectively. We also use $s_p$ to describe the s of the edge between the p-th and the $p+1$-th phases. According to the definition, $s_f$ and $s_{N_p}$ are equal.

The deterministic multiple-phase optimal control problem consists of the following constraints. We discuss them using generic expressions in this section and show the details in Section 4.1.

Constraint D1: Aircraft EoM

The aircraft trajectory for the phase p is expressed as a set of differential equations:

(17)$\frac{d{\mathit{x}}^{<p>}}{ds}=\mathit{f}\left({\mathit{x}}^{<p>},{\mathit{u}}^{<p>}\right)\mathrm{for}p=1,\cdots ,N_p.$

We assume $\mathit{f}$ to be unchanged across the considered phases.

Constraint D2: Initial/final conditions

The initial and final conditions can be specified as:

(18)${\mathit{g}}_{init}\left({\mathit{x}}^{<1>}\left(s_0\right)\right)=\mathbf{0},{\mathit{g}}_{fin}\left({\mathit{x}}^{<N_p>}\left(s_f\right)\right)\le \mathbf{0},$

where ${\mathit{g}}_{init}$ and ${\mathit{g}}_{fin}$ are sets of functions defining the initial and final conditions, respectively. The final conditions are allowed to be inequality constraints.

Constraint D3: Path constraints

The path constraints define limitations inside each phase, not at the phase edges:

(19)${\mathit{g}}_{low}^{<p>}\le {\mathit{g}}_{path}^{<p>}\left({\mathit{x}}^{<p>},{\mathit{u}}^{<p>}\right)\le {\mathit{g}}_{up}^{<p>}\mathrm{for}p=1,\cdots ,N_p,$

where ${\mathit{g}}_{path}^{<p>}$ is a set of functions defining the constraints, and ${\mathit{g}}_{low}^{<p>}$ and ${\mathit{g}}_{up}^{<p>}$ are the corresponding lower and upper bounds. The path constraints include operational limitations, e.g., CAS and altitude limitations, representing the flight envelope and the arrival route structure.

Constraint D4: Phase-edge event conditions

The event conditions represent requirements the trajectory has to meet at phase edges:

(20)${\mathit{g}}_{event}^{<p>}\left(\mathit{x}\left(s_p\right)\right)\le \mathbf{0}\mathrm{for}p=1,\cdots ,N_p-1.$

These conditions typically represent speed and altitude limits at waypoints.

Constraint D5: Phase-linkage conditions

The phase-linkage conditions are specific to multiple-phase formalization and connect the trajectories of the state variables at each phase edge, ensuring continuity of trajectory:

(21)${\mathit{x}}^{<p>}\left(s_p\right)-{\mathit{x}}^{<p+1>}\left(s_p\right)=\mathbf{0}\mathrm{for}p=1,\cdots ,N_p-1.$

If we define ${\mathit{g}}_{linkage}^{<p>}={\mathit{x}}^{<p>}\left(s_p\right)-{\mathit{x}}^{<p+1>}\left(s_p\right)$, the linkage conditions are expressed as:

(22)${\mathit{g}}_{linkage}^{<p>}\left({\mathit{x}}^{<p>},{\mathit{x}}^{<p+1>},s_p\right)=\mathbf{0}\mathrm{for}p=1,\cdots ,N_p-1.$

This study does not impose these conditions on the control inputs. This means we allow instant change of control at phase edges.

With the above constraints, we aim to minimize the following objective functional.

Objective: Minimum operational costs

We set operational costs as the objective for the deterministic planning to minimize: $J_{det}$ (€). The operational costs are defined as a combination of the time-related costs $J_{time}$ (€) and the fuel costs $J_{fuel}$ (€). The time-related costs are defined as the costs charged for the flight time for the considered flight phases represented by the final time $t^{<N_p>}\left(s_f\right)$. A linear relationship is assumed between the time costs and the final time through the coefficient $C_{time}$ (€/s) as: $J_{time}=C_{time}{t}^{<N_p>}\left(s_f\right)$. We model the coefficient utilizing the cost index $C_I=60·C_{time}/C_{fuel}$ (kg/min) (the unit is aligned Airbus’s $C_I$ definition [46]), where $C_{fuel}$ (€/kg) denotes the jet fuel price. If we substitute $C_{time}$ with $C_I$, we obtain

(23)$J_{time}=\frac{C_I{C}_{fuel}}{60}{t}^{<N_p>}\left(s_f\right).$

On the other hand, we define the fuel costs as the costs charged for the fuel burn $FB$ (kg) for the considered flight phases: $J_{fuel}=C_{fuel}FB$. The fuel burn can be expressed as the difference between the initial and the final aircraft mass: $FB=m^{<1>}\left(s_0\right)-m^{<N_p>}\left(s_f\right)$, as shown in Appendix B. Thus, the fuel costs are written as:

(24)$J_{fuel}=C_{fuel}\left\{m^{<1>}\left(s_0\right)-m^{<N_p>}\left(s_f\right)\right\},$

leading to:

(25)$J_{det}=J_{time}+J_{fuel}=C_{fuel}\left\{\frac{C_I}{60}{t}^{<N_p>}\left(s_f\right)+m^{<1>}\left(s_0\right)-m^{<N_p>}\left(s_f\right)\right\}.$

To summarize, the deterministic multiple-phase optimal control problem is expressed as:

$\begin{array}{ccccccc}\hfill & \mathrm{Find}\hfill & \hfill & {\mathit{u}}^{<1>},\cdots ,{\mathit{u}}^{<N_p>},s_1,\cdots ,s_f\hfill & \hfill & & \\ \hfill & \min \hfill & \hfill & \mathrm{operational}\cos \mathrm{ts}\hfill & \hfill \mathrm{Equation}\left(25\right)\\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \mathrm{aircraft}\mathrm{EoM}\hfill & \hfill \mathrm{Equation}\left(17\right)\\ \hfill & \hfill & & \mathrm{initial}/\mathrm{final}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(18\right)\\ \hfill & \hfill & & \mathrm{path}\mathrm{constraints}\hfill & \hfill \mathrm{Equation}\left(19\right)\\ \hfill & \hfill & & \mathrm{phase}-\mathrm{edge}\mathrm{event}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(20\right)\\ \hfill & \hfill & & \mathrm{phase}-\mathrm{linkage}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(22\right)\end{array}$

To facilitate the following discussion, we rewrite the optimal control problem further. We gather the variables for all phases into a single vector as $\tilde{\mathit{x}}={\{{\mathit{x}}^{<1>},\cdots ,{\mathit{x}}^{<N_p>}\}}^T$, $\tilde{\mathit{u}}={\{{\mathit{u}}^{<1>},\cdots ,{\mathit{u}}^{<N_p>}\}}^T$ to cover the entire trajectory. We use the tilde symbol to denote the variables covering all the phases. Using this all-phase expression, the EoM can be expressed as:

(26)$\frac{d\tilde{\mathit{x}}}{ds}=\mathit{f}\left(\tilde{\mathit{x}},\tilde{\mathit{u}}\right).$

If we also apply the all-phase principle with the tilde symbol to the path constraints, the phase-edge event condition and the phase-linkage conditions, they are rewritten as:

(27)${\tilde{\mathit{g}}}_{low}\le {\tilde{\mathit{g}}}_{path}\left(\tilde{\mathit{x}},\tilde{\mathit{u}}\right)\le {\tilde{\mathit{g}}}_{up},$

(28)${\tilde{\mathit{g}}}_{event}\left(\tilde{\mathit{x}},s_1,\cdots ,s_{N_p-1}\right)\le \mathbf{0},$

(29)${\tilde{\mathit{g}}}_{linkage}\left(\tilde{\mathit{x}},s_1,\cdots ,s_{N_p-1}\right)=\mathbf{0}.$

As a result, the deterministic optimal control is finally expressed as:

$\begin{array}{ccccc}\hfill & \mathrm{Find}\hfill & \hfill & \tilde{\mathit{u}},s_1,\cdots ,s_f\hfill & \\ \hfill & \min \hfill & \hfill & \mathrm{operational}\cos \mathrm{ts}\hfill & \hfill \mathrm{Equation}\left(25\right)& \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \mathrm{aircraft}\mathrm{EoM}\hfill & \hfill \mathrm{Equation}\left(26\right)\\ \hfill & \hfill & & \mathrm{initial}/\mathrm{final}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(18\right)\\ \hfill & \hfill & & \mathrm{path}\mathrm{constraints}\hfill & \hfill \mathrm{Equation}\left(27\right)\\ \hfill & \hfill & & \mathrm{phase}-\mathrm{edge}\mathrm{event}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(28\right)\\ \hfill & \hfill & & \mathrm{phase}-\mathrm{linkage}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(29\right)\end{array}$

3.2. Robust Trajectory Planning

In the development of robust optimization, several types of measure of robustness have been proposed [47]. The most common measure is achieved by optimizing the problem for the worst case among the anticipated scenarios. This strategy takes a scenario that produces the highest costs and optimizes the problem for it. Such costs are expressed with the max function, and thus this type of method is also called the min–max method. However, a solution optimized for the worst case, which usually happens rarely, is often too conservative for the other cases, which are more probable. To deal with this issue, another definition of robustness has also been proposed. It regards a robust solution as the one that is optimized for the moments of the overall costs covering any realization of the uncertainty sets: such as the expected value (the first moment) or the variance (the second moment), or even higher moments [28]. This study takes the latter definition and optimizes the trajectory in terms of the expected value of the operational costs.

We formulate the robust planning problem based on the deterministic formalization discussed in Section 3.1. We take a set of $N_w$ GEFS members as possible realizations of the uncertain weather during the flight. We define robust trajectories as (1) trajectories that are feasible in any GEFS weather scenarios in the sense that they do not violate the operational limitations, and (2) trajectories that are optimized for the expected value of the costs. Feasibility of the trajectories is achieved by making the EoM and the operational limitations hard constraints regardless of weather conditions within the uncertainty range provided by GEFS data. We apply robust optimal control similarly to [28,33]. Unlike prior works, such as [26], assuming common control inputs to different weather scenarios, we basically consider that both the aircraft states and control inputs are different for different weather scenarios. We explicitly express the variables as ${\tilde{\mathit{x}}}^{\left(w\right)},{\tilde{\mathit{u}}}^{\left(w\right)}$ with the superscript to show the variables are for the weather scenario w. On the contrary, we assume the independent variable s are in a common domain $s\in [s_0,s_f]$ for any w. We also assume $s_p$ for the intermediate phase edges $p=1.\cdots ,N_p-1$ are common to the different weather scenarios. Therefore, s, $s_0$, $s_f$ and $s_p$s do not need the superscript.

The robust optimal control problem holds the following constraints based on the corresponding deterministic constraints.

Robust constraint R1: Aggregated aircraft EoM

We apply the deterministic EoM (Equation (26)) to each weather scenario $w=1,\cdots ,N_w$ as $d{\tilde{\mathit{x}}}^{\left(w\right)}/ds=\mathit{f}({\tilde{\mathit{x}}}^{\left(w\right)},{\tilde{\mathit{u}}}^{\left(w\right)})$. We assume the aircraft dynamics do not change across scenarios, and therefore $\mathit{f}$ does not need the superscript. We call this set of $N_w$ trajectories the trajectory ensemble. The trajectory for the weather scenario w is the w-th member of the trajectory ensemble. Then we aggregate the set of EoMs as:

(30)$\frac{d}{ds}\left[\begin{array}{c}{\tilde{\mathit{x}}}^{\left(1\right)}\\ ⋮\\ {\tilde{\mathit{x}}}^{\left(N_w\right)}\end{array}\right]=\left[\begin{array}{c}\mathit{f}\left({\tilde{\mathit{x}}}^{\left(1\right)},{\tilde{\mathit{u}}}^{\left(1\right)}\right)\\ ⋮\\ \mathit{f}\left({\tilde{\mathit{x}}}^{\left(N_w\right)},{\tilde{\mathit{u}}}^{\left(N_w\right)}\right)\end{array}\right].$

We call this set of EoMs aggregated aircraft EoM. With this aggregated dynamic constraint, the trajectory ensemble is computed and optimized simultaneously.

Robust constraint R2: Commonality constraints

If we solve the aggregated EoM, a different trajectory is calculated for each weather scenario. From the operational point of view, however, such a set of trajectories will make pilots unsure about which trajectory to follow, as they do not know which weather scenario to expect. To give them single reference, this study proposes the commonality constraints limiting the aircraft’s 3D positional trajectories $\tilde{\varphi }\left(s\right)$, $\tilde{\lambda }\left(s\right)$ and $\tilde{h}\left(s\right)$ common to all weather scenarios:

(31)${\tilde{\varphi }}^{\left(w\right)}\left(s\right)={\tilde{\varphi }}^{\left(1\right)}\left(s\right),{\tilde{\lambda }}^{\left(w\right)}\left(s\right)={\tilde{\lambda }}^{\left(1\right)}\left(s\right),{\tilde{h}}^{\left(w\right)}\left(s\right)={\tilde{h}}^{\left(1\right)}\left(s\right)\mathrm{for}w=2,\cdots ,N_w.$

These commonality constraints let the problem be different from just repeating the deterministic optimization $N_w$ times, as each of the trajectory ensemble members needs to share the common positional profile and thus cannot be optimized individually.

Robust constraint R3: Initial/final conditions

We assume the trajectory planning starts with the common initial aircraft state regardless of the weather scenario. Considering Equation (18), the initial aircraft state ${\tilde{\mathit{x}}}^{\left(w\right)}\left(s_0\right)$ and the final state ${\tilde{\mathit{x}}}^{\left(w\right)}\left(s_f\right)$ must satisfy:

(32)${\mathit{g}}_{init}\left({\tilde{\mathit{x}}}^{\left(w\right)}\left(s_0\right)\right)=\mathbf{0},{\mathit{g}}_{fin}\left({\tilde{\mathit{x}}}^{\left(w\right)}\left(s_f\right)\right)\le \mathbf{0}\mathrm{for}w=1,\cdots ,N_w.$

${\mathit{g}}_{init}$ and ${\mathit{g}}_{fin}$ without the superscript $\left(w\right)$ mean these functions are independent of w.

Robust constraint R4: Path constraints

We apply Equation (27) to each weather scenario w

(33)${\tilde{\mathit{g}}}_{low}\le {\tilde{\mathit{g}}}_{path}\left({\tilde{\mathit{x}}}^{\left(w\right)},{\tilde{\mathit{u}}}^{\left(w\right)}\right)\le {\tilde{\mathit{g}}}_{up}\mathrm{for}w=1,\cdots ,N_w.$

The path constraints, including the flight envelope, operational limitations, etc., are hard and must be met regardless of the weather scenario. Therefore, the lower and upper bounds are independent of the weather scenario w.

Robust constraint R5: Phase-edge event conditions

The events at phase edges, such as speed and altitude limitations, must be hard constraints as well. Based on Equation (28), the constraints are expressed as:

(34)${\tilde{\mathit{g}}}_{event}\left({\tilde{\mathit{x}}}^{\left(w\right)},s_1,\cdots ,s_{N_p-1}\right)\le \mathbf{0}\mathrm{for}w=1,\cdots ,N_w.$

Robust constraint R6: Phase-linkage conditions

The continuity of the trajectory at any phase edge is naturally required for any weather scenario. Considering Equation (29), the robust linkage conditions are:

(35)${\tilde{\mathit{g}}}_{linkage}\left({\tilde{\mathit{x}}}^{\left(w\right)},s_1,\cdots ,s_{N_p-1}\right)=\mathbf{0}\mathrm{for}w=1,\cdots ,N_w.$

In addition to the above constraints, we introduce an additional constraint related to the arrival time at the final state.

Robust constraint R7: Required time of arrival

We consider the required time of arrival (RTA) for the final position. The trajectories of t are uncommon to trajectory ensemble members and therefore the final time $t^{<N_p>\left(w\right)}\left(s_f\right)$ is different for each w. We define the robust RTA as a constraint for the expected value of the arrival time: $E\left[t^{<N_p>}\left(s_f\right)\right]=t_{RTA}$, where the time constraint $t_{RTA}$ is a given parameter. This study assumes each GEFS member is equally probable. Under this assumption, this constraint is:

(36)$\frac{1}{N_w}\sum _{w=1}^{N_w}{t}^{<N_p>\left(w\right)}\left(s_f\right)=t_{RTA}.$

We also derive the robust objective functional based on the deterministic one.

Robust objective: Expected value of the operational costs

As discussed above, we minimize the expected value of the operational costs:

(37)$J_{robust}=E\left[J_{det}\right]=\frac{C_{fuel}}{N_w}\sum _{w=1}^{N_w}\left\{\frac{C_I}{60}{t}^{<N_p>\left(w\right)}\left(s_f\right)+m^{<1>\left(w\right)}\left(s_0\right)-m^{<N_p>\left(w\right)}\left(s_f\right)\right\}.$

To summarize the discussion, the robust optimal control problem for robust trajectory planning can be described as:

$\begin{array}{cccc}\hfill & \mathrm{Find}\hfill & \hfill & {\tilde{\mathit{u}}}^{\left(1\right)},\cdots ,{\tilde{\mathit{u}}}^{\left(N_w\right)},s_1,\cdots ,s_f\hfill \\ \hfill & \min \hfill & \hfill & \mathrm{expected}\mathrm{operational}\cos \mathrm{ts}\hfill & \hfill \mathrm{Equation}\left(37\right)\\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill & \mathrm{aggregated}\mathrm{aircraft}\mathrm{EoM}\hfill & \hfill \mathrm{Equation}\left(30\right)\\ \hfill & \hfill & & \mathrm{commonality}\mathrm{constraints}\hfill & \hfill \mathrm{Equation}\left(31\right)\\ \hfill & \hfill & & \mathrm{initial}/\mathrm{final}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(32\right)\\ \hfill & \hfill & & \mathrm{path}\mathrm{constraints}\hfill & \hfill \mathrm{Equation}\left(33\right)\\ \hfill & \hfill & & \mathrm{phase}-\mathrm{edge}\mathrm{event}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(34\right)\\ \hfill & \hfill & & \mathrm{phase}-\mathrm{linkage}\mathrm{conditions}\hfill & \hfill \mathrm{Equation}\left(35\right)\\ \hfill & \hfill & & \mathrm{required}\mathrm{time}\mathrm{of}\mathrm{arrival}\hfill & \hfill \mathrm{Equation}\left(36\right)\end{array}$

Figure 4 shows the overview of the process of the proposed robust planning strategy. The trajectory optimization starts with the common initial state ${\mathit{x}}_0$ and each weather scenario is applied to the corresponding part of the aggregated EoM to make the members of the trajectory ensemble. The ensemble can provide the expected range of trajectory variation. Each ensemble member leads to its own value of the objective, and we minimize the mean value. This means of optimization simultaneously considers $N_w$ members of the trajectory ensemble.

It is well known that optimal control problems generally do not have closed-form analytic solutions, except for a few known, simplified problems, such as the linear quadratic optimal control problem [48]. We utilize the co-called pseudospectral method, a numerical computation method where the continuous optimal control problem is discretized into so-called pseudospectral nodes. Among several kinds, we take the Legendre–Gauss pseudospectral method (LGPM). The resulting discrete optimization problem is categorized as a nonlinear programming (NLP) problem, for which a variety of solvers are widely available. We utilize IPOPT (interior point optimizer) [49], an open-source software package for large-scale nonlinear optimization. We use an open-source optimization tool CasADi [50] to interface IPOPT with Python. An equal number of pseudospectral nodes is assumed in this study for any phase p.

4. Case Studies

4.1. Scenario Settings

We apply the proposed planning algorithm to an aircraft descending to Leipzig/Halle Airport in Germany (ICAO code: EDDP) through the transition procedure for runway 08R (Figure 5). This transition route is designed for CDO operation and used only at night time (CDO NIGHT). The relevant waypoints are listed in Table 1.

The optimal trajectory is calculated from the en-route waypoint LUXAR, and the initial conditions for the simulation are as listed in Table 2. In the initial state, the aircraft is assumed to pass LUXAR with level flight at a cruising altitude of 35,000 ft with 63,700 kg of total mass. We split the descent trajectory into several phases. We take the relevant waypoints as phase edges that the aircraft flies over, though they are actually defined as fly-by waypoints. The phase-specific constraints for the optimization are listed in Table 3. We require the CAS to stay slower than 250 kt below 10,000 ft. This constraint comes from the FAA regulation (14 CFR § 91.117). $M_{Min}$ and $M_{MO}$ are the minimum Mach number and the maximum operating Mach number, respectively.

In addition to the phase-specific constraints, we also have inter-phase limitations for the control inputs as described in Table 4. In the table, $C_{T_{idle}}$ and $C_{T_{MCRZ}}$ are the thrust coefficients for the idle rating and maximum-cruise rating. To calculate the constraints in Table 3 and Table 4, we utilize BADA4’s CAS model for $V_{CAS}$; the thrust coefficient models for $C_{T_{idle}}$ and $C_{T_{MCRZ}}$; and the values of minimum/maximum Mach number for $M_{min}$ and $M_{MO}$ [41]. We take the GEFS data published on the September 22nd, 2020 at 0:00. We use 746.73 $/ton for the jet fuel price published by IATA [52] for October 8th, 2021. With relevant unit changes, the unit fuel costs are calculated as $C_f=$ 0.64 €/kg.

4.2. Robust vs. Scenario-Optimal Trajectories

In order to validate the proposed planning algorithm, we calculate deterministic trajectories optimized for each weather ensemble member. We call these trajectories scenario-optimal trajectories. They are obtained by applying the deterministic planning algorithm (Section 3.1) to each weather scenario.

Figure 6 shows the lateral paths of the robust (in red) and the scenario-optimal (in blue) trajectories for the case with $C_I=30$ kg/min without RTA, together with the waypoints and the ToD locations, and Figure 7 shows the wind speed components along the robust trajectories. Both robust and scenario-optimal trajectories reach 10,000 feet altitude at MAXEB, so the end of $d_1$ phase and the beginning of $d_3$ phase are found to be equal. With these variations in u and v trends, the scenario-optimal trajectories are varied between LUXAR and MAXEB, whereas the robust trajectories share the common lateral path due to the commonality constraints (Section 3.2). This single weather-scenario-independent path facilitates the pilot following the trajectory plan no matter the weather, which is one of the benefits of the proposed robust planning strategy. $N_w$ repeats of the scenario-optimal deterministic planning is not able to calculate such trajectories. There was no significant difference in lateral path found after MAXEB.

The corresponding vertical trajectories are shown in Figure 8. Gaps between the robust and the scenario-optimal trajectories are mainly in the $d_1$ phase except for a case where the aircraft flies lower in the following descent phases. A temporary horizontal segment is shown in the $d_3$ phase, where the aircraft flies from MAXEB at 10,000 ft height to DP808. This is mainly due to the corresponding altitude restriction. The altitude is limited between 5500 ft and 10,000 ft in this phase, and the optimal trajectories maintain 10,000 ft, the highest altitude allowed, as long as it is possible to minimize the drag.

Figure 9 depicts the arrival time at GAMKO, fuel burn and operational costs that the robust trajectory results in for each weather scenario. The average arrival time is 7566 s and its standard deviation is 25.0 s, whereas the average fuel burn is 4176 kg and its standard deviation is 17.3 kg. These arrival times and fuel burns lead to 4933 € for the average operational costs and 18.4 € for its standard deviation.

4.3. Robust vs. Inappropriately-Controlled Trajectories

We introduce another type of trajectory named the inappropriately-controlled trajectory. This is the type if the pilot follows one of the scenario-optimal trajectories in a weather scenario different from the assumed one. The inappropriately-controlled trajectories represent situations where a deterministic situation is assumed for trajectory planning and the pilot has to follow a trajectory optimized for a certain weather scenario but the true weather is different. We calculate the inappropriately-controlled trajectory for the weather scenario w by applying the scenario-optimal control inputs for the other scenarios ${\tilde{\mathit{u}}}^{\left(i\right)}(i=1,\cdots ,N_w\mathrm{and}i\ne w)$. This leads to $N_w-1$ inappropriately-controlled trajectories for each weather scenario. An inappropriately-controlled trajectory can lead to better costs than those with the robust trajectory when the assumed weather scenario i for the applied scenario-optimal control is close to the true scenario w.

Figure 10 shows the costs of the robust trajectories (red dots) and the inappropriately-controlled trajectories (box plots) for each weather ensemble member (the case with $C_I=30$ kg/min and RTA not imposed). The robust trajectories achieve better operational costs than the average costs of the corresponding inappropriately-controlled trajectories (green triangles) in all weather scenarios. Except for two scenarios (8 and 14) among 21, the robust trajectory even results in lower costs than any of the corresponding inappropriately-controlled trajectories. This means that the robust trajectories are expected to achieve better performances in all objectives of the optimization. Trajectories planned in a deterministic way can result in lower costs only when the assumed weather situation is close to the actual weather.

Violation of operational constraints is another issue when the pilot only has a deterministic trajectory plan. We can see some limitations and specifications for $V_{CAS}$ are violated in the case with $C_I=20$ kg/min and without RTA. Figure 11 shows the $V_{CAS}$ profiles of the robust and inappropriately-controlled trajectories for the $d_3,\cdots ,d_6$ phases. We can see the inappropriately-controlled trajectories contain higher $V_{CAS}$ than the robust ones, and some of them violate the 250-kt limitation in the $d_3$ phase. Other inappropriately-controlled trajectories fail to meet the 180-kt final condition at the end of $d_6$ phase. The robust trajectory meets all the constraints.

4.4. Effects of Cost-Index and RTA Variations

We also test the proposed planning framework with varied parameter settings—$C_I$ and $t_{RTA}$—to investigate how robust trajectories change with their variations. We first vary $C_I$ from 10 to 80 kg/min, while not imposing RTA. Figure 12 shows the relationship between the average arrival time at the final state (GAMKO: FAF) and the average fuel burn, depending on the values of $C_I$. When a larger $C_I$ is applied to the trajectory planning, the aircraft flies faster and arrives at GAMKO earlier to reduce the time costs. This method of flight leads to more drag and requires more thrust. Consequently, more fuel is necessary. Compared to the case with $C_I=10$ kg/min, the flight time is reduced by 564 s (7.2%) but with 310 kg (7.6%) more fuel being burnt at $C_I=80$ kg/min.

As for the RTA, we use the symbol $t_f^{*}$ to represent the average arrival time at GAMKO of the robust trajectories without RTA imposed. We vary $t_{RTA}$ from $t_f^{*}-60$ (s) to $t_f^{*}+60$ (s), and we keep $C_I=30$ kg/min when $t_{RTA}$ is varied. Figure 13 depicts the trend of the average fuel burn vs. the expected operational costs (Equation (37)) with various RTAs. Earlier arrival times than $t_f^{*}$ require more fuel burn, whereas later arrival times require less, as the aircraft needs to increase/decrease the overall speed to meet the arrival time constraint, leading to more/less thrust, respectively. When the aircraft is required to arrive at GAMKO 60 s earlier than $t_f^{*}$ in the trajectory plan, 33.4 kg more fuel is burnt. When the arrival time is 60 s later, the fuel burn is 26.5 kg less. In terms of the expected value of the total operational costs, the case without RTA achieves the minimum costs. When the average aircraft arrival time is required to be 60 s earlier than the optimal average arrival time $t_f^{*}$, the costs become 2.0 € higher than for the case without RTA. If $t_f^{*}+60$ s of RTA is required, 2.1 € more costs are necessary. The change in the operational costs is much smaller than that in the arrival time and fuel burn with the formalization and the weather data we utilize for the present study. A different definition of the cost function (Equation (37)) and a different set of GEFS data could provide larger differences. Yet, as the expected trend is shown in Figure 13 (right): the case without RTA leads to the minimum costs.

Figure 12 and Figure 13 also confirm the proposed robust planning framework is able to deal with these parameter changes. The ability to carry out this type of analysis is beneficial for operators. If we interpolate the dots on Figure 12, we can obtain a curve of arrival-time and fuel-burn trend. Such a figure provides airlines with information on how much additional fuel or additional flight time they need to expect depending on their $C_I$ setting, and is helpful for determining the value of $C_I$ to take. The same kind of interpolation can also be applied to the figures in Figure 13, providing ATCs with insights on how much additional costs and fuel burn they have to impose on the aircraft with a certain value of RTA. This information allows the ATC and the pilot to negotiate for their best RTA, considering benefits for both sides.

The simulation results we have discussed so far prove that the planning algorithm can work for our scenario’s settings and weather data. In order to validate the algorithm more thoroughly, we need to test it with a variety of conditions, including weather data, arrival routes, initial conditions and $C_I$s and RTAs. Further validation is one of the authors’ main future intentions. Study in practical applications will also be authors’ future work.

5. Conclusions

This study proposes a fundamental framework and computational algorithms for planning optimal descent trajectories that are robust against uncertainties in weather prediction. We utilized a set of weather scenarios provided by GEFS to model the uncertainties related to weather prediction. A deterministic multiple-phase optimal control problem was built to form the basis of the robust optimal control problem, considering 4D aircraft dynamics. We formalized a robust optimal control problem by aggregating the state and control variables and the aircraft dynamics for each of the weather scenarios to create an aggregated EoM. We introduced constraints characterizing the robust planning, i.e., the commonality constraints, so that the robust trajectories can have the common 3D positional profiles over the considered weather scenarios. Unlike conventional worst-case optimization, the proposed robust optimal control minimizes the expected value of the costs for the considered set of weather scenarios. We applied the proposed algorithm to descent trajectory planning for Leipzig/Halle Airport. Comparison of robust and inappropriately-controlled trajectories allowed us to confirm the proposed planning algorithm can prevent high costs and infeasible trajectories that violate operational constraints in any weather scenarios. Optimization results with different $C_I$ and RTA showed the proposed robust optimal control can cope with these varying operational settings. The analyses with different parameter settings support operators, i.e., pilots, airlines and ATCs, in determining the best values for those parameters. This study discussed a robust planning problem in a generic form and proposed a fundamental framework which can be applied to a wide range of scenario settings.

One of our main future intentions is to test the proposed planning algorithm with a variety of weather scenarios, arrival routes, initial conditions and $C_I$ and RTA settings to thoroughly validate its performance. We will also utilize the robust planning strategy to solve operation-oriented problems, such as a RTA-assignment problem for a metering FIX, where the time to reach the FIX is uncertain. In addition to these, we will apply the planning algorithm to more advanced air/space structures such as trombones and point merging to enhance its applicability.

Author Contributions

Conceptualization, S.K. and J.R.; methodology, S.K.; software, S.K.; validation, S.K. and J.R.; writing—original draft preparation, S.K.; writing—review and editing, S.K., J.R., H.F. and M.S.; visualization, S.K.; supervision, J.R., H.F. and M.S.; project administration, J.R. and H.F.; funding acquisition, H.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research was carried out as a part of the project "Optimized CDO under Uncertain Environmental and Mission Conditions" funded by DFG, Deutsche Forschungsgemeinschaft (German Research Foundation)—project number: 327114631.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The GEFS data are publicly available: https://www.ncei.noaa.gov/products/weather-climate-models/global-ensemble-forecast, accessed on 25 November 2021.

Conflicts of Interest

The authors declare no conflict of interest. This product/document has been created with or contains elements of Base of Aircraft Data (BADA) Family 4 Release 4.2 which has been made available by EUROCONTROL to TU DRESDEN. EUROCONTROL has all relevant rights to BADA. ©2021 The European Organisation for the Safety of Air Navigation (EUROCONTROL). All rights reserved. EUROCONTROL shall not be liable for any direct, indirect, incidental or consequential damages arising out of or in connection with this product or document, including with respect to the use of BADA.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. GEFS consists of [Forumla omitted. See PDF.] different ensemble members ([Forumla omitted. See PDF.] for this study). A different b-spline approximation surface (blue surface) is obtained for a different GEFS member, based on the corresponding GRIB2 data (red dots).

Figure 2. Relationship between the true air speed [Forumla omitted. See PDF.] and the ground speed [Forumla omitted. See PDF.].

Figure 3. Overview of the constraints and symbols specific to phases. The constraints are categorized as phase-interior constraints and phase-edge constraints. Continuity of the trajectory is assumed for [Forumla omitted. See PDF.] at phase edges.

Figure 4. Overview of the process of the proposed robust planning strategy.

Figure 5. The horizontal route for the simulation. The arrival route for the CDO operation (CDO NIGHT) of RWY08R at Leipzig/Halle Airport, Germany (original source in black: AIP Germany [51], modified for research purpose in purple). The dotted lines and the symbols [Forumla omitted. See PDF.] show the phase separation. The lengths of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] phases do not necessarily reflect reality in this figure.

Figure 6. Lateral paths of the robust (red) and the scenario-optimal (blue) trajectories (with [Forumla omitted. See PDF.] kg/min, RTA not imposed). The red dots correspond to the waypoints and the black dots show the ToD locations of the trajectories.

Figure 7. Trends of u and v along the robust trajectory for each weather scenario. The components are positive for west-to-east wind and south-to-north wind, respectively.

Figure 8. Vertical paths of the robust (red) and the scenario-optimal (blue) trajectories (with [Forumla omitted. See PDF.] kg/min, RTA not imposed). The trajectories cover the last part of c-phase and the following descent phases. The vertical lines corresponds to phase edges of the robust trajectory.

Figure 9. Arrival time at GAMKO, fuel burn and operational costs of the robust trajectory for each weather scenario ([Forumla omitted. See PDF.] kg/min, RTA not imposed).

Figure 10. Total operational costs for each weather scenario (with [Forumla omitted. See PDF.] kg/min, RTA not imposed). The box plots show the ranges of costs of the inappropriately-controlled trajectories for the corresponding weather scenarios, and their mean values are depicted with the green triangles. Each box plot is based on the costs of [Forumla omitted. See PDF.] inappropriately-controlled trajectories. The red dots show the costs for the robust trajectory.

Figure 11. CAS profiles in the last 4 phases (with [Forumla omitted. See PDF.] kg/min, RTA not imposed). The black profiles are inappropriately-controlled trajectories and the red are the robust trajectories. The vertical lines are the phase edges and the horizontal dashed red lines show the limitations for the corresponding phases. The red diamond at the end of [Forumla omitted. See PDF.] phase depicts the final condition.

Figure 12. Trend of the average arrival time at GAMKO vs. the average fuel burn of the robust trajectories with [Forumla omitted. See PDF.] kg/min. RTA is not imposed.

Figure 13. (Left) Trend of the average arrival time at GAMKO vs. fuel burn: ±60 s from the average arrival time without RTA [Forumla omitted. See PDF.], i.e., [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] kg/min. (Right) Trend of the average arrival time at GAMKO vs. the robust operational costs [Forumla omitted. See PDF.] (Equation (37)) with the same RTA range.

Table 1

Waypoints used for the simulation.

Name	Coordinates	Altitude	Speed
MAXEB	N51 12.4 E012 13.9	FL80 or above	IAS 250 kt
DP808	N51 19.3 E011 48.9	5500 ft or above	IAS 230 kt
DP807	N51 23.5 E011 48.4	5000 ft or above	IAS 210 kt
DP442	N51 23.8 E011 53.5	-	-
GAMKO (FAF)	N51 24.1 E011 59.7	3000 ft or above	IAS 180 kt

Table 2

Initial states for the waypoint LUXAR.

States	Symbols	Conditions
Along-track distance	$s_0$	0 m
Coordinates	$({\varphi }_0,{\lambda }_0)$	(N49${}^{°}$55${}^{\prime }$48.00${}^{″}$, E021${}^{°}$10${}^{\prime }$31.00${}^{″}$)
Altitude	$h_0$	35,000 ft
Mass	$m_0$	63,700 kg
Time	$t_0$	0 s

Table 3

Phases and phase-specific constraints.

Phase	Description	Path Constraints	Final Conditions
c	Cruise	$M_{min}\le M\le M_{MO}$
		$h=h_0$, $\gamma =0$ ${}^{°}$, ${\delta }_{SB}=0$
$d_1$	ToD to 10,000 ft	$M_{min}\le M\le M_{MO}$	h = 10,000 ft
		$V_{CAS}\le 350\mathrm{kt}$	230 kt $\le V_{CAS}\le$ 250 kt
$d_2$	10,000 ft to MAXEB	230 kt $\le V_{CAS}\le$ 250 kt	$(\varphi ,\lambda )={(\varphi ,\lambda )}_{MAXEB}$
			$h\ge$FL80, $V_{CAS}=$250 kt
$d_3$	MAXEB to DP808	230 kt $\le V_{CAS}\le$ 250 kt	$(\varphi ,\lambda )={(\varphi ,\lambda )}_{DP808}$
			$h\ge$5500 ft, $V_{CAS}=$230 kt
$d_4$	DP808 to DP807	210 kt $\le V_{CAS}\le$ 230 kt	$(\varphi ,\lambda )={(\varphi ,\lambda )}_{DP807}$
			$h\ge$ 5000 ft, $V_{CAS}=$ 210 kt
$d_5$	DP807 to DP442	180 kt $\le V_{CAS}\le$ 210 kt	$(\varphi ,\lambda )={(\varphi ,\lambda )}_{DP442}$
			$h\ge$ 3000 ft
$d_6$	DP442 to GAMKO	180 kt $\le V_{CAS}\le$ 210 kt	$(\varphi ,\lambda )={(\varphi ,\lambda )}_{GAMKO}$
			$h=$ 3000 ft, $V_{CAS}=$ 180 kt

Table 4

Limitations for the control inputs across the phases.

Controls	Symbols	Lower Bounds	Upper Bounds
Thrust coefficient	$C_T$	$C_{T_{idle}}$	$C_{T_{MCRZ}}$
Flight path angle	$\gamma$	−4${}^{°}$	0${}^{°}$
Heading	$\psi$	−180 deg	180 deg
Speed brake	${\delta }_{SB}$	0	1

Appendix A. Three Dimensional B-Spline Function

We interpolate the GRIB2-style GEFS data (Equation (1)) using the three-dimensional tensor-product b-spline function:(A1)$$f^{bs\left(w\right)}(\varphi ,\lambda ,p)=\sum _i\sum _j\sum _k{c}_{i,j,k}^{\left(w\right)}{B}_i^{\left(w\right)}\left(\varphi \right)B_j^{\left(w\right)}\left(\lambda \right)B_k^{\left(w\right)}\left(p\right)\mathrm{for}w=1,\cdots ,N_w,$$ where f is either $u,v$ or $T_{mp}$, $c_{i,j,k}^{\left(w\right)}$ is the b-spline coefficient; and $B_i^{\left(w\right)}$, $B_j^{\left(w\right)}$ and $B_k^{\left(w\right)}$ are the b-spline basis functions for the w-th member of the weather ensemble. The b-spline coefficients and functions are calculated from GEFS data and are ensemble-member specific.

Appendix B. Derivation of Fuel Burn

The fuel burn $FB$ is calculated from the fuel consumption $FC$:(A2)$$FB={\int }_{t^{<1>}\left(s_0\right)}^{t^{<N_p>}\left(s_f\right)}FCdt=\sum _{p=1}^{N_p}{\int }_{t^{<p>}\left(s_{p-1}\right)}^{t^{<p>}\left(s_p\right)}FCdt,$$ where the fuel burns for each phase p are summed up. Considering the transformation of the independent variable (Section 2.2), we change the variable in the integral:(A3)$$FB=\sum _{p=1}^{N_p}{\int }_{s_{p-1}}^{s_p}FC\frac{dt}{ds}ds=\sum _{p=1}^{N_p}{\int }_{s_{p-1}}^{s_p}\frac{FC}{V_{GS}}ds.$$

$FC/V_{GS}$ can be substituted by $-dm/ds$ due to Equation (16). As a result, we can obtain:(A4)$$\begin{array}{ccc}\hfill FB& =& \sum _{p=1}^{N_p}{\int }_{s_{p-1}}^{s_p}\left(-\frac{dm}{ds}\right)ds=-\sum _{p=1}^{N_p}{\int }_{m^{<p>}\left(s_{p-1}\right)}^{m^{<p>}\left(s_p\right)}dm\hfill \\ & =& \sum _{p=1}^{N_p}\left\{m^{<p>}\left(s_{p-1}\right)-m^{<p>}\left(s_p\right)\right\}\hfill \end{array}$$

The phase-linkage condition (Equation (22)) can simplify the summation as:(A5)$$FB=m^{<1>}\left(s_0\right)-m^{<N_p>}\left(s_f\right).$$

References

1. Fricke, H.; Seiß, C.; Herrmann, R. Fuel and energy benchmark analysis of continuous descent operations. Proceedings of the 11th USA/Europe Air Traffic Management Research and Development Seminar (ATM Seminar); Lisbon, Portugal, 23–26 June 2015.

2. European Commission. Flightpath 2050: Europe’s Vision for Aviation: Maintaining Global Leadership and Serving Society’s Needs; Publications Office: Brussels, Belgium, 2011; [DOI: https://dx.doi.org/10.2777/50266]

3. SESAR Joint Undertaking. European ATM Master Plan: Digitalising Europe’s Aviation Infrastructure: Executive View: 2020 Edition; Publications Office: Brussels, Belgium, 2020; [DOI: https://dx.doi.org/10.2829/650097]

4. Study Group for the Future Air Traffic Systems. Long-Term Vision for the Future Air Traffic Systems—Changes to Intelligent Air Traffic Systems; The Ministry of Land, Infrastructure, Transport and Tourism (MLIT): Tokyo, Japan, 2010.

5. Federal Aviation Administration (FAA). NextGen Implementation Plan 2018–2019; Office of NextGen: Washington, DC, USA, 2018.

6. European Union (EU). Single European Sky Performance Scheme; EU regulation No. 290/2013 European Union: Brussels, Belgium, 2013.

7. International Civil Aviation Organization (ICAO). Continuous Descent Operations (CDO) Manual; Doc. 9931 International Civil Aviation Organization: Montreal, QC, Canada, 2010.

8. Continuous Climb and Descent Operations. European Organisation for the Safety of Air Navigation (Eurocontrol). Available online: https://www.eurocontrol.int/concept/continuous-climb-and-descent-operations (accessed on 25 November 2021).

9. European Organisation for the Safety of Air Navigation (Eurocontrol). European CCO/CDO Action Plan; Publications Office: Brussels, Belgium, 2020.

10. Vertical Flight Efficiency at Airports. European Organisation for the Safety of Air Navigation (Eurocontrol). Available online: https://ansperformance.eu/efficiency/vfe/ (accessed on 25 November 2021).

11. Toratani, D.; Wickramasinghe, N.K.; Hirabayashi, H. Simulation techniques for arrival procedure design in continuous descent operation. Proceedings of the 2018 Winter Simulation Conference; Gothenburg, Sweden, 9–12 December 2018; [DOI: https://dx.doi.org/10.1109/WSC.2018.8632350]

12. Clarke, J.P.; Brooks, J.; Nagle, G.; Scacchioli, A.; White, W.; Liu, S.R. Optimized profile descent arrivals at Los Angeles international airport. J. Aircr.; 2013; 50, pp. 360-369. [DOI: https://dx.doi.org/10.2514/1.C031529]

13. Park, S.G.; Clarke, J.P. Vertical trajectory optimization for continuous descent arrival procedure. Proceedings of the AIAA Guidance, Navigation, and Control (GNC) Conference; Minneapolis, Minnesota, 13–16 August 2012; [DOI: https://dx.doi.org/10.2514/6.2012-4757]

14. Dalmau, R.; Prats, X. Fuel and time savings by flying continuous cruise climbs estimating the benefit pools for maximum range operations. Transp. Res. Part D Transp. Environ.; 2015; 35, pp. 62-71. [DOI: https://dx.doi.org/10.1016/j.trd.2014.11.019]

15. De Jong, P.M.A. Continuous Descent Operations Using Energy Principles. Ph.D. Thesis; Delft University of Technology: Delft, The Netherlands, 2014; [DOI: https://dx.doi.org/10.4233/uuid:69cd815f-a8bd-4602-bb37-57109a880495]

16. Dalmau, R.; Prats, X. Controlled time of arrival windows for already initiated energy-neutral continuous descent operations. Transp. Res. Part C Emerg. Technol.; 2017; 85, pp. 334-347. [DOI: https://dx.doi.org/10.1016/j.trc.2017.09.024]

17. Lindner, M.; Rosenow, J.; Zeh, T.; Fricke, H. In-flight aircraft trajectory optimization within corridors defined by ensemble weather forecasts. Aerospace; 2020; 7, 144. [DOI: https://dx.doi.org/10.3390/aerospace7100144]

18. Lindner, M.; Rosenow, J.; Fricke, H. Aircraft trajectory optimization with dynamic input variables. CEAS Aeronaut. J.; 2020; 11, pp. 321-331. [DOI: https://dx.doi.org/10.1007/s13272-019-00430-0]

19. Rosenow, J.; Lindner, M.; Scheiderer, J. Advanced flight planning and the benefit of in-flight aircraft trajectory optimization. Sustainability; 2021; 13, 1383. [DOI: https://dx.doi.org/10.3390/su13031383]

20. Franco, A.; Rivas, D.; Valenzuela, A. Probabilistic aircraft trajectory prediction in cruise flight considering ensemble wind forecasts. Aerosp. Sci. Technol.; 2018; 82–83, pp. 350-362. [DOI: https://dx.doi.org/10.1016/j.ast.2018.09.020]

21. Hernàndez-Romero, E. Probabilistic Aircraft Conflict Detection and Resolution under the Effects of Weather Uncertainty. Ph.D. Thesis; Universidad de Sevilla: Sevilla, Spain, 2020.

22. Franco, A.; Rivas, D.; Valenzuela, A. Optimal aircraft path planning in a structured airspace using ensemble weather forecast. Proceedings of the 8th SESAR Innovation Days; Salzburg, Austria, 3–7 December 2018.

23. Legrand, K.; Puechmorel, S.; Delahaye, D.; Zhu, Y. Robust aircraft optimal trajectory in the presence of wind. IEEE Aerosp. Electron. Syst. Mag.; 2018; 33, pp. 30-38. [DOI: https://dx.doi.org/10.1109/MAES.2018.170050]

24. Fisher, J.; Bhattacharya, R. Optimal trajectory generation with probabilistic system uncertainty using polynomial chaos. J. Dyn. Syst. Meas. Control; 2011; 133, 014501. [DOI: https://dx.doi.org/10.1115/1.4002705]

25. Cottrill, G.C. Hybrid Solution of Stochastic Optimal Control Problems Using Gauss Pseudospectral Method and Generalized Polynomial Chaos Algorithms. Ph.D. Thesis; Air Force Institute of Technology: Kaduna, Nigeria, 2012.

26. Li, X.; Nair, P.B.; Zhang, Z.; Gao, L.; Gao, C. Aircraft robust trajectory optimization using nonintrusive polynomial chaos. J. Aircr.; 2014; 51, pp. 1592-1603. [DOI: https://dx.doi.org/10.2514/1.C032474]

27. Piprek, P. Robust Trajectory Optimization Applying Chance Constraints and Generalized Polynomial Chaos. Ph.D. Thesis; Technische Universität München: München, Germany, 2020.

28. González-Arribas, D.; Soler, M.; Sanjurjo-Rivo, M. Robust aircraft trajectory planning under wind uncertainty using optimal control. J. Guid. Control. Dyn.; 2018; 41, pp. 673-688. [DOI: https://dx.doi.org/10.2514/1.G002928]

29. García-Heras, J.; Soler, M.; González-Arribas, D. Characterization and enhancement of flight planning predictability under wind uncertainty. Int. J. Aerosp. Eng.; 2019; 2019, 6141452. [DOI: https://dx.doi.org/10.1155/2019/6141452]

30. Soler, M.; González-Arribas, D.; Sanjurjo-Rivo, M.; García-Heras, J.; Sacher, D.; Gelhardt, U.; Lang, J.; Hauf, T.; Simarro, J. Influence of atmospheric uncertainty, convective indicators, and cost-index on the leveled aircraft trajectory optimization problem. Transp. Res. Part C Emerg. Technol.; 2020; 120, 102784. [DOI: https://dx.doi.org/10.1016/j.trc.2020.102784]

31. González-Arribas, D.; Andrés-Enderiz, E.; Soler, M.; Jardines, A.; Garcıa-Heras, J. Probabilistic 4D flight planning in structured airspaces through parallelized simulation on GPUs. Proceedings of the 9th International Conference for Research in Air Transportation (ICRAT 2020); Tampa, FL, USA, 23–26 June 2020.

32. Matsuno, Y.; Tsuchiya, T.; Wei, J.; Hwang, I.; Matayoshi, N. Stochastic optimal control for aircraft conflict resolution under wind uncertainty. Aerosp. Sci. Technol.; 2015; 43, pp. 77-88. [DOI: https://dx.doi.org/10.1016/j.ast.2015.02.018]

33. Kamo, S.; Rosenow, J.; Fricke, H.; Soler, M. Robust CDO trajectory planning under uncertainties in weather prediction. Proceedings of the 14th USA/Europe Air Traffic Management Research and Development Seminar (ATM Seminar); Virtual Event, 20–23 September 2021.

34. Global Ensemble Forecast System (GEFS). National Oceanic and Atmospheric Administration (NOAA). Available online: https://www.ncei.noaa.gov/products/weather-climate-models/global-ensemble-forecast (accessed on 25 November 2021).

35. Zhou, X.; Zhu, Y.; Hou, D.; Kleist, D. A comparison of perturbations from an ensemble transform and an ensemble kalman filter for the NCEP Global Ensemble Forecast System. Weather. Forecast.; 2016; 31, pp. 2057-2074. [DOI: https://dx.doi.org/10.1175/WAF-D-16-0109.1]

36. Rödel, W. Physik unserer Umwelt, die Atmosphäre; 2nd ed. Springer: Berlin/Heidelberg, Germany, 2000.

37. International Civil Aviation Organization (ICAO). Manual of the ICAO Standard Atmosphere—Extended to 80 kilometres/262,500 feet; Doc. 7488 International Civil Aviation Organization: Montreal, QC, Canada, 1993.

38. Peters, M.; Konyak, M. The Engineering Analysis and Design of the Aircraft Dynamics Model for the FAA Target Generation Facility; Technical Report; Air Traffic Engineering Co., LLC: Lakewood, NJ, USA, 2012.

39. National Imagery and Mapping Agency (NIMA). Department of Defense World Geodetic System 1984, Its Definition and Relationships With Local Geodetic Systems; National Imagery and Mapping Agency: Bethesda, MD, USA, 1997.

40. Walter, R. Flight management systems. The Avionics Handbook; CRC Press LLC: Boca Raton, FL, USA, 2001; Chapter 15

41. Nuic, A.; Mouillet, V. User Manual for the Base of Aircraft Data (BADA) Family 4; 12/11/22-58, Version 1.3 European Organisation for the Safety of Air Navigation (Eurocontrol): Brussels, Belgium, 2016.

42. Bronsvoort, J. Contributions to Trajectory Prediction Theory and its Application to Arrival Management for Air Traffic Control. Ph.D. Thesis; Universidad Politècnica de Madrid: Madrid, Spain, 2014.

43. Dalmau, R.; Prats, X.; Baxley, B. Using wind observations from nearby aircraft to update the optimal descent trajectory in real-time. Proceedings of the 13th USA/Europe Air Traffic Management Research and Development Seminar (ATM Seminar); Vienna, Austria, 17–21 June 2019.

44. Sáez, R. Traffic Synchronization with Controlled Time of Arrival for Cost-Efficient Trajectories in High-Density Terminal Airspace. Ph.D. Thesis; Universitat Politécnica de Catalunya: Barcelona, Spain, 2021; [DOI: https://dx.doi.org/10.13140/RG.2.2.22429.64485]

45. Kamo, S.; Rosenow, J.; Fricke, H. CDO sensitivity analysis for robust trajectory planning under uncertain weather prediction. Proceedings of the 2020 AIAA/IEEE 39th Digital Avionics Systems Conference (DASC); San Antonio, TX, USA, 11–15 October 2020; [DOI: https://dx.doi.org/10.1109/dasc50938.2020.9256699]

46. Airbus. Getting to Grips with the Cost Index; Flight Operations Support and Line Assistance, Issue II, STL 945.2369/98 Airbus Customer Services: Blagnac, France, 1998.

47. Beyer, H.G.; Sendhoff, B. Robust optimization—A comprehensive survey. Comput. Methods Appl. Mech. Eng.; 2007; 196, pp. 3190-3218. [DOI: https://dx.doi.org/10.1016/j.cma.2007.03.003]

48. Todorov, E. Optimal Control Theory. Bayesian Brain: Probabilistic Approaches to Neural Coding; The MIT Press: Cambridge, MA, USA, 2006; [DOI: https://dx.doi.org/10.7551/mitpress/9780262042383.003.0012]

49. Wächter, A.; Biegler, L.T. On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Math. Program.; 2006; 106, pp. 25-27. [DOI: https://dx.doi.org/10.1007/s10107-004-0559-y]

50. Andersson, J.A.E.; Gillis, J.; Horn, G.; Rawlings, J.B.; Diehl, M. CasADi—A software framework for nonlinear optimization and optimal control. Math. Program. Comput.; 2019; 11, pp. 1-36. [DOI: https://dx.doi.org/10.1007/s12532-018-0139-4]

51. AIP Germany. Available online: https://www.eisenschmidt.aero/ (accessed on 1 December 2020).

52. Jet Fuel Price Monitor. International Air Transport Association (IATA). Available online: https://www.iata.org/en/publications/economics/fuel-monitor/ (accessed on 25 November 2021).