Abstract
A novel study is presented aiming at characterizing and illustrating potential enhancements in flight planning predictability due to the effects of wind uncertainty. A robust optimal control methodology is employed to calculate robust flight plans. Wind uncertainty is retrieved out of Ensemble Probabilistic Forecasts. Different wind approximation functions are compared, typifying errors, and illustrating its importance for accurate solving of the robust optimal control problem. A set of key performance indicators is defined for the quantification of uncertainty in terms of flight time and fuel consumption. Two different case studies are presented and discussed. The first one is based on a representative sample of the whole 2016 year for a single origin-destination and a forecast time step of 6 hours. As for the second, we select the most uncertain day together with a multiorigin-destination set of flights with forecast time steps up to 2 days.
1. Introduction
The Air Traffic Management (ATM) System is undergoing a paradigm shift aiming at enhancing its environmental impact, capacity, safety, and efficiency. Different worldwide programmes are fostering this transformation, mainly the Single Sky ATM Research (SESAR) in Europe and Next Generation Air Transportation System (NextGen) in USA, yet also in the Asia/Pacific region, the Collaborative Actions for Renovation of Air Traffic Systems (CARATS) in Japan and the New Generation ATM System (CNAS) in China. Improvement in the ATM system predictability has been identified as paramount to achieve the above-mentioned high-level goals (see, e.g., SESAR-JU. European ATM Master Plan Edition 2015). All ATM system actors, including pilots, air traffic controllers, airlines, dispatchers, air navigation service providers, meteorological offices, and the network manager, are daily facing the effects of uncertainty. Should the capacity of the system be increased while maintaining high safety standards and improving the overall performance, uncertainty levels in ATM must be reduced and, when possible, find new strategies to deal with and eventually reduce it. For instance, it is well-known that the sector’s capacity is under-declared due to uncertain entry/occupancy counts, thus limiting the capacity of the system. Also, adding unnecessary fuel leads to a nonneglectable loss of efficiency, e.g., Hao et al. in [1] conclude that 1 min in reducing flight time dispersion could save between $120 and $452 million per year only considering US domestic airlines.
Uncertainty in ATM is nevertheless a heavily involved, multilayered, and interrelated phenomena. The analysis of uncertainty in ATM should take into the different scales, yet also the different sources that introduce uncertainty into the system. For deeper insight on this taxonomy, refer to [2] (Chapter 4). From a scale perspective, ATM uncertainty can be grouped into macroscale (the air transport network), mesoscale (traffic scale), and microscale (single flight). Previous work analyzing uncertainty effects at both macro- and mesoscales include, e.g., Goyens et al. [3], who studied delay propagation across the European network; Valenzuela et al. in [4], who studied probabilistic sector loading considering wind uncertainties; and Cook et al. [5], who applied complex network theory to the ATM system. Needless to say, the focus of the present paper is on flight uncertainty (microlevel).
Regarding flight uncertainty, [2] (Chapter 4) distinguishes four time frames, namely, strategic (flow management planning level, from months up to two/three hours before the off-block time), predeparture (flight dispatching planning level, from two/three hours up to off-block time), gate-to-gate (execution both airborne and taxi times), and postarrival (from on-block time). Source-wise, ATM uncertainty falls into the following five categories: data uncertainty, operational uncertainty, equipment uncertainty, and meteorological uncertainty. Indeed, meteorological uncertainty has been recognized as one of the most (if not the most) important sources introducing uncertainty into the ATM system (see also [2]). All five uncertainty sources obviously affect flight predictability at its different time frames, which is mainly due to [6] not perfect knowledge of the initial conditions (e.g., off-block time, departure time, and initial mass), inaccurate aircraft performance model (e.g., weight, aerodynamics, Earth models, and atmospheric variables), uncertain meteorological conditions (e.g., convection and wind), navigation and equipment errors (e.g., flight management system, avionics, and flight technical errors), and ATM-related operational uncertainty (regulations, ATC advisories). Some previous work on initial condition error estimation include, e.g., Sun et al. [7], who estimate aircraft initial mass using ADS-B data. Also, Vazquez and Rivas [8] studied the propagation of initial mass uncertainty. In [6], the authors use Polynomial Chaos Expansion (PCE) to propagate uncertainties along a given flight intent considering disturbances due to initial conditions, the aircraft performance model, the Earth and atmospheric model, the aircraft intent, and wind.
We are herein interested in both the strategic and predeparture time frames, the ones compatible with flight planning. Yet, analyzing meteorological uncertainty and its effects on flight planning, to the best of the authors’ knowledge, has received not enough attention so far. Indeed, wind uncertainty modelling in [6] considers uniform wind field distributions based on a deterministic forecast. More advanced features are demanded, such as the use of Ensemble Prediction System (EPS) forecasts, which provide probabilistic forecasts.
One of the first (if not the first) works on EPS-based wind uncertainty modelling and flight trajectories is due to Cheung et al. [9]. They used EPS information to study the impact of wind uncertainty in flight duration on a given, fixed route over the North Atlantic. A great impact in flight duration was observed, with no seasonal variation, and an important decrease when the time frame gets closer to the used forecast. Moreover, Vazquez et al. [10] analyzed the influence of along-track wind uncertainty (also based on EPS) in fuel mass distribution based on a probabilistic transformation method. Both studies can be casted as trajectory prediction problems since the flight path is known a priori, and thus, they are only capable of characterizing uncertainty. The ambition in this paper is not only to characterize but also to reduce uncertainty. For that purpose, flight planning algorithms that take into consideration predictability should be derived. Research to this end is rather recent and still needs more insight.
For flight planning purposes, a rather straightforward approach is to extend some of the probabilistic trajectory prediction approaches (e.g., [9, 10]) to consider a higher algorithmic level in which the flight path is found using the discrete optimization method. This was followed by Cheung et al. in [11], where they used a Dijkstra-based trajectory predictor and also evaluated the quality of different EPSs, and by Franco et al. [12], who addressed the optimization of a North Atlantic crossing also using a Dijkstra algorithm together with their probabilistic trajectory predictor based on the probabilistic transformation method. A different approach is due to Legrand et al. [13], who presented a robust wind (obtained from EPS) optimal trajectory design algorithm based on dynamic programming and minimum-energy trajectory clustering. It should be noted that all these three approaches for flight planning under wind uncertainty are somehow based on discrete optimization. In [14], Gonzalez-Arribas et al. presented a robust optimal control approach to flight planning under wind uncertainty (obtained from EPS) in which the flight plan can be optimized considering efficiency and predictability. The present paper is based on the latter methodology, with the aim at generalizing it to different scenarios and modelling features, and thus advance in the understanding of the effects of meteorological uncertainty.
All in all, the contribution of this paper is threefold. First, we extend the methodology in [14] to consider spline wind approximations. It is shown that the errors committed in wind approximation (using, e.g., polynomial approximation such in [14]) have a tremendous impact, leading to ill-conditioned solutions. Second, instead on analyzing a single flight (as in previous work), we generalize for different days (and thus meteorological conditions) and for different origin-destinations. We also analyze the effects of wind uncertainty at different flight planning horizons, i.e., predeparture (3 hours before departure), which could be useful for flight dispatchers, and strategic (1 day before departure), which could be useful for the network manager for probabilistic capacity-demand balancing. Third, we characterize the quality of the flight plan prediction (and its uncertainty) against a hypothetical realization of winds based on reanalysis information. We note (as will be exposed) that further insights in the latter are demanded in future research and this third contribution remains as a preliminary step.
The methodology has been sketched in Figure 1. Firstly, wind data is collected from the EPS forecast. Secondly, wind forecasts are approximated by Lipschitz-continuous functions in order to model wind is a form compatible with the robust optimal planning methodology in [14]. Thirdly, the so-termed Robust Optimal Control Problem (ROCP) is solved, obtaining a Pareto frontier with different trade-offs between fuel consumption and predictability. Then, the flight path (obtained as a solution to the ROCP and defined by a sequence of waypoints) is introduced into a standard trajectory predictor (TP) ([15], Chapter 3), which flies the precalculated route using wind information retrieved from the original source (raw, tabular), either EPS or reanalysis. Finally, flight times, fuel burns, and its associated dispersions are analyzed. Different Key Performance Indicators (KPIs) are defined to that end.
Two different case studies are presented. The first focuses on a single prediction horizon (referred to as step), choosing a single origin-destination flight (over the North Atlantic), and a systematic date analysis along year 2016 (choosing the 1st, 5th, 10th, 15th, 20th, and 25th of each month as representative days). The second focuses on a single day (5th October 2016, the one that resulted to be the most unpredictable one in the former analysis), multiple origin-destination flights (based on the schedule of an important airline), and with forecasting steps ranging from 0 to 48 hours (by 6 hours). The first study aims at characterizing the seasonal influence at the predeparture time horizon, i.e., compatible with flight dispatching time frames; the second aims at the geographical influence, yet on a more strategic vision, i.e., compatible with the network manager’s time frames.
2. Robust Flight Planning
2.1. Robust Optimal Control Methodology
Our approach is based on the robust optimal control methodology for aircraft trajectory optimization problems presented in [14]. We will here summarize the method for the sake of completeness.
2.1.1. Robust Optimal Control
Let us consider a dynamical system with a randomly parametrized differential-algebraic equation and constraints. Uncertainty is described using a standard Kolmogorov probability space , composed of a sample space of possible outcomes , a -algebra of events containing sets of outcomes, and the probability function assigning probabilities to each of these events. The uncertain parameters of the system will be modelled as a no time-dependent random variable . For each possible outcome , the random variables take a different value .
Let us denote the state vector , the control vector , the algebraic variables , and as the independent variable (typically time). For each outcome , there exist a unique trajectory path that corresponds to the realization of the random variables . The dynamics of the system are given by the functions , , and , such that valid trajectories fulfill the conditions almost surely (i.e., with probability 1) (the ≤ sign applies in an element-wise fashion in equation (3) and analogous equations.): where is the sample point on the underlying abstract probability space. Therefore, for each possible scenario or realization of the random parameters , the trajectory will follow the deterministic differential equation (1) for the corresponding fixed value of . We employ the notation and in order to emphasize the fact that the trajectory depends on the realization of the random parameters. Note that, because the random parameters are constant in time and thus the system is not described by a full stochastic differential equation (as in, for example, [16]), there is no need to consider filtered -algebras to deal with potential measurability issues.
We follow the optimal guidance scheme of [14] that relies on the notion of tracked states.
Definition 1. A state is said to be tracked if its trajectory is assumed to be independent of the realization of the random variables almost surely, i.e., with probability one.
In this concept, the controls are specific to each scenario in order to guarantee that the tracked states follow the unique computed trajectory of those states that are tracked. Naturally, the practical consideration of this framework requires a dynamical system where the controls can be computed online to ensure that the realization of the trajectory of the tracked states stays close to its precomputed value. This is the case of the aircraft, since the flight control systems of an aircraft can follow a given route and vertical profile. For an in-depth explanation of this concept and why it is adequate for a flight planning problem, we refer the reader to [14]. We now proceed to illustrate the main concepts.
Definition 2. The amount control degrees of freedom of the dynamical system is defined, in this work, as the number of controls and free algebraic variables minus the number of algebraic restrictions: .
Let be the number of tracked states; without loss of generality, we can assume that the tracked states are the first states (rearrange the state vector otherwise), i.e., where is the tracked part of the state vector and is the untracked part. Let be the identity matrix of shape and be the zero matrix (i.e., a matrix with zeroes in all its entries) of shape . We define the matrix as
This matrix transforms the state vector into the tracked states vector that contains only the states whose evolution is equal in all scenarios. In this work, will hold so that there are enough tracked states to consume all available control degrees of freedom. This is not necessary in general. With the aid of the tracking matrices, we can now define the tracking conditions (which, again, apply almost surely):
The tracking conditions enforce equality in the tracked variables between realizations: note that, is the vector of differences between the tracked states in outcome and the tracked states in outcome . The other two conditions are analogous tracking conditions for the dependent variables and the controls.
Let us define a Bolza-type functional in optimal control problems:
For each possible scenario or realization of the random parameters , we have an objective value of . In order to build a scalar objective function, we rely on both expectations and dispersions. Without generality, the latter can be characterized using the central moments, the range, or the interquartile range (IQR), just to mention a few. In our case, and in order to account for worst-case scenarios, we choose the range.
Let us then define the range of as
The robust optimal control problem with tracking is now formulated as follows (an alternative formulation could, for instance, rely on the central moment of a random variable as and then substituting the range in by .) [14]: subject to: initial conditions final conditions differential-algebraic equations and constraints (1)–(3) grouped as and the tracking conditions (6).
Where in the above, is the expectation operator associated with the probability space ; the terminal cost or the Mayer term is ; the running cost or the Lagrange term is ; dp is a weighting parameter (dispersion penalty); and function denotes the final conditions.
2.1.2. Probabilistic Discretization
In order to deal with potentially continuous random variables, we make use of discretization techniques. These methods will be referred to as stochastic quadrature rules (SQRs).
Definition 3. Given a random variable with probability measure , a stochastic quadrature rule (SQR) is defined as a function or algorithm that produces a finite set of quadrature points , and weights , that approximates the integral with the sum: with the approximation converging to the true value as .
In the present work, we employ a trivial quadrature rule (each ensemble member is a scenario with ), as uncertainty is already characterized by discrete scenarios from EPS forecasts. However, the integration of additional sources of uncertainty in the future may require the usage of a nontrivial SQR.
Let define a trajectory for the tracked states and and analogously. Suppose a SQR has been chosen, with a number of points (in our case, corresponding to the possible scenarios). For each one of these points , the tracking trajectory defines a unique trajectory given a full set of initial conditions; we will now collect each one of these trajectories in a trajectory ensemble. We define the trajectory ensemble associated with a tracking trajectory as the set of trajectories with such that the trajectory is generated by the initial conditions and the tracking trajectory with , i.e.,
We will now build a virtual dynamical system where the state vectors of all the trajectories in the trajectory ensemble are merged into a single large state vector, which represents the collective trajectory. Its state vector contains the state vector of all the trajectories in the ensemble (the control and algebraic vectors are analogous), and each individual trajectory follows the dynamics associated with each scenario (with ).
We define the differential equation, algebraic equations, and inequality constraints of this augmented dynamical system as
Let us define the approximate cost functional, divided into Mayer and Lagrange terms, using the trajectory ensemble: and discretize the boundary conditions as
For concise writing of the discretization of the tracking conditions (6), we will define the matrix as
and can be defined in analogous fashion. These matrices map the ensemble state vector to the differences in the tracked states between trajectories.
Making use of equations (14), (15), (16), (17), and (18) as well as equations (19), (20), (21), (22), and (23), we can complete now the formulation of the deterministic approximant: where
2.2. Application to Flight Planning
2.2.1. Modelling Assumptions
We consider a free-routing airspace and 3-DoF point-mass model based on BADA 3 [17]. We restrict ourselves to the analysis of the cruise phase for the sake of illustration. In addition, we assume the constant flight level and airspeed profiles. Note however that our methodology could be extended to full 4D problems. We consider an ellipsoidal Earth as in the WGS84 model (which offers slightly more precision than a spheroid model), with radii of curvature of ellipsoid meridian and prime vertical denoted by and , respectively. We take the wind and temperature fields from an EPS forecast and compute density with the ideal gas law (as the pressure is determined by the flight level). We assume the heading as a control variable instead of the bank angle, thus allowing it to change instantaneously.
All in all, aircraft dynamics can be described by the following set of differential equations expressed in the Geographic Reference System (GRS) coordinate system ( in quation (12)): where is the latitude, is the longitude, is the true airspeed, is the geodetic altitude, is the heading, and and are the zonal and meridional components of the wind. The control vector is composed by the heading .
2.2.2. System Reformulation
We reformulate the dynamical system above as a differential-algebraic system (DAE) with the addition of the ground speed as an algebraic variable and the course as a control variable, linked to the remaining variables by two new equality constraints. The reformulated system is given by the dynamical function: and the equality constraints ( in equation (12)):
We add the inequality constraint: to ensure uniqueness of and (otherwise, would produce the same left-hand side of equation (27) as ).
2.2.3. Time to Distance Transformation
We apply a coordinate transformation to make the distance flown along the route () be the new independent variable. This allows us to apply our methodology in a manner that is consistent with existing planning and flight procedures (again, see the discussion in [14]). As a consequence, the time becomes a state variable, and the new dynamical function can be obtained by dividing the time derivatives by :
All constraints remain the same as in the untransformed system of differential-algebraic equations.
2.2.4. Problem Statement
An ensemble forecast contains a set of ensemble members, each one defining a different wind forecast (and, therefore, different functions and ). Consider the ensemble contains members, then we define scenarios, each one having weight . We choose to track the course , i.e., the function is the same in every scenario (thus, we do not need to implement scenario-specific versions). As a consequence of (29), this implies that the evolution of the latitude is unique (as it only depends on the evolution of the unique variable ) and is also unique (as it only depends on and ). Therefore, the position variables also act like tracked variables, which is a desired goal (since we want to obtain a unique route).
We can define the dynamical system associated with the trajectory ensemble with the dynamical function:
The boundary conditions are where and are scalar decision variables.
We complete the definition of the discretized robust optimal control problem by adding the cost function. We apply the robust optimal control methodology to find routes that minimize a weighted sum of average fuel consumption and flight time dispersion (weighted with the dispersion penalty parameter dp). The cost functional to minimize is then
Note that is the difference between the earliest and the latest arrival time.
With this cost functional, the parameter dp regulates the solution flight plan depending on the preferences of the flight planner. High values of dp will produce more predictable trajectories by avoiding regions where the wind is more unpredictable. Low values of the parameter will produce flight plans that are more efficient in average but less predictable.
3. Trajectory Predictor
The trajectory predictor determines the aircraft trajectory for a weather scenario. The objective of the TP is to provide the aircraft variable values, i.e., flight times, aircraft mass, and ground speed, along the previously defined route (sequence of nodes/waypoints) provided by the ROC problem.
Mass at discrete node/waypoint is computed following quation (33) (note that this equation can be applied because the aircraft is flying at constant speed and altitude): where where and are first and second thrust-specific fuel consumption coefficient, respectively, is the true airspeed (TAS), is a constant to convert to knots, is the reference wing surface area, is the air density, is the parasitic drag coefficient (cruise), is the induced drag coefficient (cruise), is the gravity, and are the flight time at discrete node/waypoint and node/waypoint , respectively, and is the aircraft mass at node/waypoint .
Ground speed is computed as where
Note that is the aircraft course between node/waypoint and node/waypoint calculated as shown in equation (37) to (38).
Flight time at node/waypoint is computed as the distance () between nodes/waypoints and over . It is assumed that is constant between nodes/waypoints. Distance between nodes/waypoint is the orthodromic distance.
4. Wind Data and Modelling
4.1. Wind Data
In order to incorporate wind into the models, we make use of two meteorological products, namely, EPS, which provide ensemble-wise uncertainty information and reanalysis, considered to be a sufficiently good approximation to weather’s realization.
4.1.1. EPS
Uncertainty of wind fields will be derived from EPS. Ensemble forecasting is a prediction technique that generates a representative sample of the possible future states of the atmosphere. An ensemble forecast is a collection of typically 10 to 50 weather forecasts (referred to as members) with a common valid time, which can be obtained using different Numerical Weather Prediction (NWP) models with varying initial conditions. Each forecast is based on a model which is close, but not identical, to the best estimate of the model equations, thus representing also the influence of model uncertainties on forecast error. Thus, the spread of solutions can be used as a measure of uncertainty. In this paper, we focus on the output data of the global ensemble forecast system ECMWF EPS. Data can be accessed (among others) at the TIGGE dataset by the European Center for Medium-Range Weather Forecasts (ECMWF) [18] (http://apps.ecmwf.int/datasets/http://apps.ecmwf.int/datasets/).
(1) ECMWF EPS. The European Centre for Medium-Range Weather Forecasts EPS (ECMWF EPS) is based on 51 ensemble members with approximately 32 km resolution up to forecast day 10 and 65 km resolution thereafter, with 62 vertical levels.
4.1.2. Reanalysis
Meteorological reanalysis is based on meteorological data assimilation. The aim is at assimilating historical observational data spanning an extended period, using a given model. Reanalysis is considered to be the best estimate on many variables (such as winds and temperature).
ECMWF ERA-5 reanalysis will be used in the course of the validation activities. The data assimilation system used to produce ERA-5 includes a 4-dimensional variational analysis with a 12-hour analysis window. The spatial resolution of the dataset is approximately 31 km (in the high resolution version) on 32 vertical levels from the surface up to 0.01 hPa. ERA-5 products are normally updated once daily, with a delay of two months to allow for quality assurance and for correcting technical problems with the production, if any. ERA-5 data can be downloaded from the ECMWF Public Datasets.
4.2. Wind Approximation
Note that in order to incorporate winds in Problem (24), functions and in equation (27) must ideally be class functions, i.e., continuous and twice derivable functions (in lack of this, they should be at least Lipschitz continuous functions). This in practice translates into the need of approximating raw, tabular data from either EPS or reanalysis into smooth functions with the abovementioned properties.
We describe two possible approximations: a multiregion polynomial regression and a spline interpolation method. Approximation errors cannot be obviated; they end up distorting results as we will show in the sequel.
4.2.1. Multiregion Polynomial Regression
We seek to compute an approximation of the wind field in the area of interest. We subdivide this area into regions and consider an approximating polynomial for each region: where denotes the index of the region.
We compute the coefficients by solving the following least squares problem: where the points are generated by applying to a rectangular, equispaced grid on the affine transformation such that , where denotes the region number . The values are generated by sampling a cubic spline interpolant of the original grid values at this new grid.
We also add continuity and differentiability requirements at the boundaries between regions. Consider the boundary between region and region and parametrize it by by building the affine transformations . By composition, we can build the polynomials: that represent the value of the field approximation at the boundaries at each one of the adjacent regions. By construction, both and are polynomials of degree in ; therefore, by equating them at points, we can guarantee that they are equal. We choose equispaced points along the boundary, () and augment problem (41) with the equality constraints . Note that these are constraints on the value of the coefficients , which are the decision variables for this problem. This constraints guarantee continuity between regions.
In order to deal with differentiability constraints, we similarly define the polynomials: where are the orders of the differentiability constraints. The degree of this polynomial on is , so we only need to equate and in points in order to enforce differentiability across boundaries.
Figure 2 presents the approximation error of the global and regional regression approaches, measured in a band of 3 deg around the reference trajectory as computed in [14]. A significant improvement can be observed in both number of regions and degree of the polynomial, thus justifying the advantages of this approach. Nonetheless, a residual error always persists and that is why we resort to spline interpolation.
(a)
(b)
4.2.2. Spline Interpolation
An alternative method to do wind field spline interpolation is a well-established field of knowledge. Amodei and Benbourhim in [19] introduced a new family of vector field splines with application to wind fields. More recently, Le Guyader et al. [20] studied the use of a spline-bases approximation considering conservative vector fields as it happens with winds (winds derive from temperature potentials). We herein use the b-spline interpolation functionality implemented in CasADi [21], a symbolic framework for algorithmic differentiation and numeric optimization. The spline takes values in the region of interest and maps them to .
4.3. Wind Errors
To illustrate errors with both wind approximation models and its effects in the results to Problem (24), we bring herein an example. It corresponds to case 5th October 2016 in the single O-D flight case study, so please refer therein for more details.
Needless to say, it can be readily observed comparing Figures 3(a) and 3(b) that the wind uncertainty (computed as the standard deviation of wind norm across all members in the EPS, see equation (44)) in both multiregion polynomial regression and spline interpolation approaches substantially differ. One could qualitatively say that there is a displacement of the wind uncertainty. This is mainly due to the fact that wind approximation errors (see definition bellow) committed using the multiregion polynomial regression itself are of the same order of magnitude as the uncertainty. This fact can be observed in Figure 3(d); indeed, the higher the uncertainty, the higher the error. This is of course expected, since regression is a smoothing technique. On the contrary (see Figure 3(c)), spline interpolation errors are neglectable with orders of magnitude of . where is the wind uncertainty, and and represent the standard deviation of wind east-west and north-south components across all members in the EPS, respectively. where and are the wind east-west and north-south components of member in the EPS and and represent the average with east-west and north-south components across all members in the EPS, respectively, being is the number of members in the EPS, i.e.,
(a) Wind uncertainty (n-spline approx.)
(b) Wind uncertainty (polynomial approx.)
(c) Wind uncertainty approximation error (n-spline approx.)
(d) Wind uncertainty approximation error (polynomial approx.)
Let us define wind uncertainty error as , where is the standard deviation of wind across all members in the EPS forecast (considering raw, tabular wind data) and is the standard deviation of wind across all members in the EPS (considering the polynomial/n-spline approximation functions.)
Figure 4 provides evidence on errors’ effects in the calculated trajectories (solving Problem (24) with both wind approximations). While darker trajectories (representing the problem computed with low dp values, thus weighting less the effects of uncertainty) are very similar, brighter ones behave in rather disimilar manner. Indeed, in Figure 4(b), trajectories computed with the multiregion polynomial regression approximation and high dp values (bright colors) seek to deviate southwards to avoid a region of high uncertainty (according to the approximated wind). This region is indeed biased by approximation errors. If one looks at Figure 4(a) (approximated with splines and thus with no error), the high-uncertainty region has moved slightly south and thus trajectories do not go south anymore, rather go straight.
(a) Using n-spline wind approximation
(b) Using polynomial approximation
We bring now readers’ attention to Figure 5. Flight times and flight dispersions are represented for different dp values in boxplots, which provide qualitative information on quartiles (upper and lower box) and median (red line). The whiskers indicate variability outside the upper and lower quartiles. Outliers are plotted (should they exist) as individual points. All in all, one can readily have an idea of median (or mean) flight times and dispersions. The solution to Problem (24) should provide increasing average flight times and reduced time dispersions as dp values grow. This is as expected in Figure 5(a), where results are represented with the direct output of the optimization problem (in which wind reads as in the approximated functions). However, if ones takes the obtained discrete path () and computes flight times using a trajectory predictor (TP) over the original raw, tabular wind, results do not show this behaviour anymore. This can be noticed in Figure 5(b), e.g., comparing and . Also, looking now at Figure 5(c), which includes the same simulation with the spline interpolation (note that herein wind errors are zero), flight times are substantially different; yet, dispersions are consistent.
(a) Using wind polynomial approx
(b) Using wind polynomial approx and TP
(c) Using wind spline approx
Please refer to Table 1 for quantitative data, including also data for fuel consumptions and fuel dispersions with spline wind approximation. Remark also that this analysis has been made for the whole scenario in Case 1 and this reasoning holds for the whole set of flights. These findings motivate the use of precise wind approximation methodology, though reader should advert there is always a trade-off between accurate modelling and solvability. All in all, we resorted to spline interpolation to circumvent these problems.
5. Key Performance Indicators (KPIs)
The following variables will be analyzed: (i)Total flight time (FT)(ii)Total fuel consumption (FC)
For the sake of simplicity, henceforth, variable will be sued (denoting either FT or FC interchangeably). Note the reader that the solution to the robust flight planning provides a discrete number N (equal to the number of members in the ensemble forecast) of possible values for both flight time and fuel consumption. Thus, in order to compare the robust flight planning results with hypothetical realizations (considered to be that of the reanalysis), we should statistically characterize flight times and fuel consumptions. Figure 6 sketches the solutions provided by the robust trajectory predictor (including its statistical characterization in quartiles and mean values) and the realization (based on reanalysis).
Let us consider the following sets: (i) (set of flights)(ii) (set of dp values)(iii) (set of forecast steps in the EPS)(iv) (set of ensemble members in the EPS)
Let us consider the following variables: where super-index R-Trj refers to variables obtained from the robust trajectory predictor and computed using EPS forecasts, and super-index Rea refers to the trajectory flown under a hypothetical realization of wind (considered to be the reanalysis).
Taking mean values across all members in the EPS:
Consider also the following conditions: where is a user-specified tolerance.
Let us define the indicator functions that takes value one if the condition is true and zero elsewhere:
The flight-wise Key Validation Indicators are defined as follows: where represents the predicted flight time or fuel burnt dispersions using the robust flight planing methodology; represents the deviations of the expected value (the mean) with respect to the observation (assumed herein to be the calculated using reanalysis); denotes boolean, and thus, the associated indicator can take value true or false, depending on whether the realization fits into the defined domains (IQR; ) for both flight time and fuel consumption; last but not the least, denotes the percentage of predictions that fin within the realization a given tolerance. Notice that these indicators are calculated for all flights, all dp values and all forecast steps.
Let us also define the following aggregated KPIs: where they provide an aggregated vision over all flights in the different scenarios of the same variables as above, namely, flight time and fuel burnt dispersions, and deviations of the predicted mean flight time and fuel burnt with respect to the observation.
6. Case Studies
Two case studies have been selected in order to characterize and illustrate potential reductions in both flight time and fuel burnt uncertainty, namely, (1)A single origin-destination flight (over the North Atlantic) along a set of characteristic days along year 2016. The flights are computed considering a single step of 6 hours, i.e., a characteristic times at flight dispatching level. Please refer to Section 6.1 for additional details(2)A multiple origin-destination set of flights in a representative day, i.e., day with more time dispersion in the single origin-destination case study. The flights are computed considering multiple steps, from 0 to 48 hours, i.e., characteristic times at flow management level. Please refer to Section 6.2 for additional details
All robust flight planning have been simulated using the robust flight planning approach in Section 2.1, considering EPS forecasts as input data and approximated using the spline interpolation in Section 4. The realization of the trajectory is computed using the route provided by the robust flight planning (a sequence of (lat, long) coordinates) and considering reanalysis wind data approximated by spline interpolation. A trajectory predictor is used for that purpose. KPIs introduced in Section 5 are used to quantify results.
6.1. Single O-D Flight, Single Step, Multiple Days
We choose a flight between KJFK and LPPT, flying at constant barometric FL380, and constant Mach of 0.82. A free routing airspace is considered. Aircraft is modelled as a A330 BADA3 model with an initial mass of 200000 kg. We select a set of characteristic days along year 2016, all together 72 days (1st, 5th, 10th, 15th, 20th, and 25th of each month). The corresponding EPS forecasts produced at 00.00 of the day and with forecast of 06 step are considered. We present in Table 1 the top five days in terms of time and fuel dispersions, respectively. We proceed on discussing the results for the day with the greatest dispersion (2016-Oct-05, herein termed representative day).
6.1.1. Representative Day “2016-Oct-05”
Figure 7 includes boxplot information on flight times, fuel consumptions, time dispersions, and fuel burnt dispersions for this representative day and different values ranging 0 to 50. Furthermore, the realization computed using reanalysis is represented (as a star). It can be seen that flight time dispersions and fuel burnt consumptions systematically reduce as increases. Similarly, expected flight times and fuel burns increase. It is also interesting to remark that and decrease as dp increases; in other words, the predicted values get closer to the realization. Indeed, it is remarkable that the greater the dp value is, the more likely the realization to fit within the predicted bounds.
(a) Flight times
(b) Fuel consumptions
Quantitative information is presented in Table 1 (row 1). For the sake of illustration, it can be observed that maximum dispersions can reach up to roughly 337 sec. This dispersion could be reduced to 207 sec. (130 sec. or 38% reduction) by setting . This would come at a cost of roughly 381 sec. of extra flight time in average . Similarly, in terms of consumption, dispersion can be reduced from 487 kg to 298 kg (189 kg or 38%), all in all with a cost of roughly 550 kg of extra fuel burnt. The Pareto frontier of the problem weighting dispersion and time/consumption is presented in Figure 8. The reader is referred therein and to Table 1 for a quantitative insight on other intermediate trade-offs and additional days.
(a) Pareto front time dispersion vs. average flight time
(b) Pareto front fuel dispersion vs. average fuel consumption
Moreover, the deviations of the expected time (as predicted) with respect to the realization is decreasing as dp increases, from 132 sec. to 84 sec. (negative sign means the trajectory flown under reanalysis has greater flight time than under the EPS forecast), i.e, roughly 48.2 sec. or, in other words, predicted trajectories with are 47.6% more predictable. Similarly with fuel consumption , we reduce the deviations with respect to the expected consumption from 191 kg to 121 kg, i.e., 70 kg of reduction.
Last but not the least, it holds (as expected) that the higher dp, the greater the percentage of prediction fitting within the realization tolerances . On the contrary, it should be noted that the realization is not fitting within the predicted interquartile range (the same reasoning holds mutatis mutandis for FC), neither within the expected flight time 30 sec. range . Notice however that the realization hits closer to the IQR as dp grows (see Figure 7). This is partially due to the high uncertainty. In other words, we are able to substantially improve predictability; however, the predictability is still rather poor. The reader is referred to Section 6.3 for a discussion on accuracy of realizations (as considered herein).
6.1.2. 2016 Aggregated Results
Figure 9 shows aggregated values for dispersion and for all 72 days considered. Table 2 presents all aggregated KPIs. On average, flight time dispersion can be reduced roughly 1 min (from 151 sec. to 92 sec.) by computing flight plans with dp equal to 20. Similarly, fuel consumption dispersion could be reduced roughly 85 kg (from 219 to 133 kg). This represents 39% dispersion reduction in both times and fuel consumption. This nevertheless would come at an average cost of extra 300 sec. and 450 kg of fuel.
(a) Final time dispersion (max-min value)
(b) Fuel consumption dispersion (max-min value)
As for the indicators that take into consideration the realizations, we can observe that deviations from the expected times and fuels decrease slightly (from 29.5 to 27.3 sec. and 42.7 to 39.4 kg, respectively). It is however remarkable to look at the percentage of times the realization is within the IQR (around 40%), yet also at the percentage of time flights are within the expected value a given tolerance (30 sec. and 70 kg in this case), over 60% and almost 90% for times and consumptions, respectively. Again, readers are referred to Section 6.3 for a discussion on accuracy of realizations (as considered herein).
6.2. Multiple O-D Flights, Multiple Steps, Single Day
We keep the focus on the same day “2016-Oct-05.” We choose now a set of 164 flights flying to different destinations (according to British Airways’ schedule for that day). Figure 10 provides a snapshot of the flight set. We assume all of them fly at constant barometric altitude FL380 and constant Mach M0.82. A free routing airspace is also assumed. Different aircraft types are considered, all modeled as BADA3, and considering their corresponding reference masses as initial mass. As for wind, we take EPS forecasts produced on 2016-Oct-05 at 00.00 with forecast steps ranging from 00 to 48 hours.
For the sake of representativeness, we have clustered them into 64 short-, 46 medium-, and 54 long-haul flights. We present the top three long-haul flights in Table 3. Top three medium- and short-haul flights can be seen in Table 4 and Table 5, respectively. All in all, results show a consistent behaviour; dispersions are systematically reduced as we increase dp values; yet, they increase as we increase the forecast step. Needless to say, this corresponds to the expected behaviour. We discuss in what follows flight EGLL to SBGL (first row in Table 3).
6.2.1. Representative Flight EGLL to SBGL
First, it is interesting to look at how dispersions grow as prediction horizon (or step) does. Indeed, at one day look-ahead time, they can be more than 10 minutes; yet, at two days look-ahead time, they can go up 20 minutes. Second, it is also very important to observe how we are able to substantially reduce flight time uncertainty as we increase the dp parameter. For instance, at one day prediction step, we are able to cut it by almost two (459 to 249 sec., i.e., 46% reduction) by moving from to . All in all, this could be of course very relevant for the ATFM units at pretactical level, where they do the capacity-demand balancing the day before operation. Nonetheless, this increase of predictability would come at the cost of increasing flight time and fuel consumption (roughly 800 sec. or 1000 kg in this case). Other trade-off values (also for the other top 2 flights) can be checked both in Table 3 and Figure 11, where the Pareto frontiers of the problem are shown for this specific flight. Moreover, Figure 12 shows boxplot information, including the obtained realization. It can be seen that dispersions reduce as dp increases. On the contrary, expected flight times and fuel burns increase. It is also interesting to remark that in general, and decrease as dp increases; in other words, the predicted values get closer to the realization. Nonetheless, no conclusive results can be extracted in this sense as forecast step increases. Similar analysis holds for other flights. Of course, medium- and short-haul flights present lower dispersion values. Detailed information on them can be checked in Tables 3, 6–9.
(a) Pareto front time dispersion vs. average flight time
(b) Pareto front fuel dispersion vs. average fuel consumption
(a) Flight times (spline)
(b) Fuel consumptions (spline)
6.2.2. Aggregated Results
An aggregated analysis (for all 164 flights considered) is presented in both Table 5 and Figure 13. Also, Figures 13–16 show aggregated values considered as per haul. One can readily observe that fuel and time dispersions are systematically reduced as we increase dp values; yet, they increase as we increase the forecast step. Aggregated results in Table 5 show, e.g., for one day planning horizon, roughly 75 sec. of time dispersion reduction with an associated cost of roughly 400 sec. of flight time (580 kg of extra fuel). If one looks at forecast step 48, one would get 133 sec. of reduction in time dispersion with a flight time cost of 664 sec. (almost 1000 kg). As for the indicator associated with the realization, we can see that the realizations are typically close to the expected values, meaning that we would be quite predictable in aggregated terms.
(a) Flight time dispersion (max-min value)
(b) Fuel consumption dispersion (max-min value)
(a) Flight time dispersion (max-min value)
(b) Fuel consumption dispersion (max-min value)
(a) Flight time dispersion (max-min value)
(b) Fuel consumption dispersion (max-min value)
(a) Flight time dispersion (max-min value)
(b) Fuel consumption dispersion (max-min value)
6.3. Discussion on the realization
So far, we have observed a systematic consistency in the behaviour of the robust trajectory planning: predicted dispersions are reduced as dp grows, yet at the cost of flying longer. This is however not the case when computing the flight as with the realizations of wind based on reanalysis. This is in principle expected, since we do not know a priori what the wind would be. Nonetheless, we have observed some features that indicate further research is needed. (i)First, one should expect the realization to be within the EPS. However, this is not always the case. We have indeed found strong deviations, which are in general very dependant of the geographic coordinates. This is due to EPS usually presenting spread-error correlation, tending to be biased and underdispersed. Indeed, different techniques for the statistical postprocessing and calibration of EPS are around, e.g., based on ensemble model output statistics and Bayesian model averaging. Thus, research efforts in improving the quality of the EPS (at the barometric levels) are needed(ii)Reanalysis is not realized wind; it is just an assimilation. Aircraft-derived data should be used instead (AMDAR data). Nevertheless, these date are not uniform (they populate the network of airways), bringing in associated drawbacks (e.g., the need of estimation outside this network). Thus, using AMDAR data would lead to similar problems as using reanalysis. Initiatives such as the recent launch of Eolus by ESA, aiming at providing a worldwide gridded measurement of winds at different altitudes, will be paramount in the future
7. Conclusions
This study explored the characterization and enhancement of flight planning predictability under wind uncertainty. For these purposes two analyses (single and multiple origin-destination; single and multiple step prediction horizon) were carried out.
The single origin-destination flight (KJFK-LPPT), over a set of days along year 2016, was computed considering a forecast step used that represents the characteristic flight dispatching level. In average, the cost of increasing 1 sec. of predictability is 5 seconds of extra flight time. Also in average, the cost of reducing fuel burnt uncertainty 1 kg is 5 kg of extra fuel burnt. For some extreme days, increasing 1 sec. predictability (reducing also fuel burnt uncertainty by 1 kg) would cost 2.5 sec. (2.5 kg) of extra fuel burnt. This improvement in flight predictability has some direct relationship with airline savings via fuel reserves, yet also has important indirect effects such as less ATC tactical intervention, less delays, and thus potential miss connections. Further studies on quantifying both these direct and indirect effects should be addressed. As for the multiple origin-destination analysis, we covered a global study considering 164 flights with different origin-destinations. The forecast step used (up to two days) represents the characteristic air traffic flow management pretactical and tactical time horizons. Similarly, for 1 to 2 days forecast steps, the cost of increasing 1 sec. of predictability is between 5 and 5.5 seconds of extra flight time (also, 5 to 5.5 kg of extra fuel burnt for reducing 1 kg of fuel burnt uncertainty). Yet, we can observe that the higher the forecasting step, the higher the uncertainty: uncertainty grows between 50 and 70 secs. in flight times every 12 hours when considering minimum time trajectories. We are however able to reduce flight time uncertainty by roughly 40% in average. The direct effect of this would be a more accurate sector loading forecasting, leading to increase the capacity of the system. Indirectly, airlines would benefit from less regulations.
Data Availability
This research is using weather data from the ECMWF MARS public databases (https://www.ecmwf.int); therefore, the work can be replicated. Also, BADA data was used and it is accessible from the Eurocontrol website (https://simulations.eurocontrol.int/solutions/bada-aircraft-performance-model/).
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work is partially supported by the Spanish Government through the project entitled “Analysis and Optimization of Aircraft Trajectories under the Effects of Meteorological Uncertainty-OptMet” (TRA2014-58413-C2-2-R) (https://optmet.wordpress.com/). The very first ideas of the present paper were introduced in a talk in the Workshop on Uncertainty and Air Traffic Management organized by the OptMet project on Wednesday 25th, October 2017 in Madrid (Spain).