Machine-Learning-Assisted Optimization of Aircraft Trajectories Under Realistic Constraints

Open AccessEngineering NotesMachine-Learning-Assisted Optimization of Aircraft Trajectories Under Realistic ConstraintsYifan Xu, Sebastian Wandelt, Xiaoqian Sun, Yichen Yang, Xianfei Jin, Sureshan Karichery and Maciej DrwalYifan XuBeihang University, 100191 Beijing, People’s Republic of China*Ph.D. Candidate, School of Electronic and Information Engineering; .Search for more papers by this author, Sebastian WandeltBeihang University, 100191 Beijing, People’s Republic of China†Professor, School of Electronic and Information Engineering; .Search for more papers by this author, Xiaoqian SunBeihang University, 100191 Beijing, People’s Republic of China‡Professor, School of Electronic and Information Engineering; .Search for more papers by this author, Yichen YangSabre Corporation, Southlake, Texas 76092§Senior Data Science Engineer, 3150 Sabre Drive; .Search for more papers by this author, Xianfei JinSabre Corporation, Southlake, Texas 76092¶Principle Data Science Engineer, 3150 Sabre Drive; .Search for more papers by this author, Sureshan KaricherySabre Corporation, Southlake, Texas 76092**Senior Principle Data Science Engineer, 3150 Sabre Drive; .Search for more papers by this author and Maciej DrwalCAE Corporation, Québec City, Québec H4T 1G6, Canada††Principal Operations Research Engineer, 8585 Cote-de-Liesse Saint-Laurent; .Search for more papers by this authorPublished Online:10 Jun 2023https://doi.org/10.2514/1.G007038SectionsRead Now ToolsAdd to favoritesDownload citationTrack citations ShareShare onFacebookTwitterLinked InRedditEmail AboutI. IntroductionWith the rapid expansion of air traffic over the past decades, increasing capacity/demand imbalances in the airspace system have raised concerns about inefficient aircraft operation, particularly in light of the induced greenhouse gas emissions. For example, in 2019, the fuel consumption of the airline industry reached 95 billion gallons with 905 million tons CO2 emissions [1]. On the other hand, it is estimated that nearly 1.5 million flights (160 million passengers) cannot be accommodated by the current air transportation system by the year 2040 [2]. Although COVID-19 has caused a significant, yet temporary setback to the aviation industry, it is expected that the aviation system will recover within the next one to two years. A broad discussion of COVID-19 and its impact on aviation can be found in [3,4]. To tackle the inherent challenges of sustainable aviation growth and to gain efficiency for both individual aircraft and the entire system, the aviation industry has bundled its efforts to promote a set of operational enhancements. One representative transformation is trajectory-based operation (TBO), which enables the air traffic management system to know and modify the planned and actual trajectory before operation or during flight based on accurate shared information. Besides, to overcome the aviation sector’s efficiency, capacity, and environmental problems, free route airspace (FRA) has been implemented in Europe since 2009. Within FRA, users may plan a route freely between a defined entry and exit point, even without reference to the air traffic service (ATS) route network. Through improved strategic planning, these improvements jointly help airspace users to fly close to their preferred trajectory while obtaining measurable improvements in direct operating costs and greenhouse emission reductions. Therefore, the 4D trajectory of an aircraft, which consists of the three spatial dimensions plus time, is fundamental in air transportation operations. While improvements in aircraft manufacturing and engine performance have played important roles in reducing fuel consumption [5], the improvement of flight trajectories has proven to be effective in significantly reducing fuel consumption and noise levels in recent years, e.g., through continuous descent approach [6]. Therefore, flight trajectory optimization has been extensively studied in the literature to improve the overall performance of the air transportation system.As the foundation to model traffic flow and aircraft characteristics such as fuel consumption and speed, the aircraft performance model is widely studied and adopted in the trajectory optimization problem. Such a critical model is typically developed and maintained by aircraft manufacturers and organizations such as the Base of Aircraft Data (BADA) [7] from Eurocontrol. Emerging research efforts are devoted to developing open-source performance models based on a point-mass model like BlueSky and OpenAP [8,9], or machine-learning-based methods, including Gaussian process regression for fuel burn prediction [10] and k-nearest neighbor for various statistics, like descent rates and airspeed [11]. In general, the aircraft performance can be integrated into trajectory optimization via equations of motion [12,13] and numerical performance data approach [14,15].In the work of [16], a multiphase mixed-integer optimal control technique is proposed and transcribed into a nonlinear programming model using the collocation method. While airways are not considered in this study under the TBO assumption, a multipartite graph is considered to represent possible aircraft transitions. Similarly, the mixed-integer optimal control approach adopted in [17] models the point-mass motion of the aircraft and allows decisions on 4D routing and airspeed. To address the uncertainty in wind prediction, González-Arribas et al. [18] apply the robust optimal control approach to aircraft trajectory planning based on a randomly parameterized ordinary differential equation. The problem is solved by approximating the uncertain parameters using stochastic quadrature rules.While the equations of motion approach provides a uniform structure for aircraft performance, solving such a nonlinear differential system is not trivial and can be too time-consuming for airline real practice. Therefore, a wide range of studies use numerical performance data, which contain a set of aircraft performance tables derived from real or experimental flight data, to facilitate the computation of critical indicators such as fuel flow and true airspeed based on interpolation in aircraft status (e.g., ISA deviation and mass weight). Accordingly, the trajectory optimization problem can be formulated as a combinatorial problem and solved using graph search or enumeration methods [19]. In the work of [20], a space reduction algorithm is developed to select the optimal flight profile. Specifically, the aircraft performance database is adopted and airspeed is computed using Lagrange linear interpolation method. A beam search method is adopted in [21] to optimize vertical profile based on a numerical performance model for all flight phases.Secondly, regarding trajectory dimensions (i.e., latitude, longitude, altitude, and time), the complete trajectory optimization problem is frequently handled by separated lateral and vertical profile optimization due to tractability issues. In particular, lateral optimization aims to find a set of geographical waypoints during the cruise phase to reduce flight time. Following the grid network modeling ideas of [15,22], the search space can be constructed by discretizing the geodesic cruise trajectory into a fixed number (e.g., 15) of waypoints. Thereafter, parallel edges are built with decreasing number of waypoints for cruise phase only (given top of climb/descent). Compared to lateral profile optimization, vertical profile optimization tries to optimize the combination of speed schedules and flight levels that reduce fuel consumption and flight time cost. In the work of [14], strategies are developed to reduce the formation of contrail aside from optimizing fuel consumption. An enumeration-based algorithm developed by [23] selects the pre-optimal cruise profile and evaluates the nearby complete profile. The authors further extend their work by modeling the problem as a decision tree and considering some realistic air traffic management constraints [21].The optimization of 3D or 4D aircraft trajectories often involves solving lateral and vertical problems sequentially [19,24,25]. For the fully integrated approach that simultaneously optimizes the vertical and lateral profile, 3D networks composed of coordinates/waypoints and discrete flight levels (e.g., defined in increments of 1000 ft) are developed in multiple studies [26–28]. Similarly, a grid system is defined in [29] to facilitate the development of a dynamic programming algorithm with equations of motion. Based on a geodesic trajectory, Murrieta-Mendoza et al. [30] generate search space between the top of climb and the top of descent with a fixed separation. These approaches aim to reduce computational complexity by optimizing the problem with a partial space around an initial trajectory. Although the integration of lateral and vertical decisions can significantly reduce operational costs, it can lead to a very large-scale network, particularly for long-haul flights with refined altitude intervals. Moreover, due to the high complexity caused by the surging solution space, relevant studies typically do not thoroughly consider dynamic weather conditions and operational constraints related to FRA, compulsory routes and so other factors. Therefore, fully integrated trajectory optimization still lacks sufficient and tractable research efforts.To conclude, depending on the choice of aircraft performance modeling and problem dimensions, the existing research on flight trajectory optimization is classified and summarized in Table 1.Sequentially optimizing lateral and vertical profiles actually ignores the coupling relationship between dynamic weather conditions and aircraft performance. For instance, in the route optimization phase, the lack of vertical information and corresponding weather data makes the computation of fuel consumption inaccurate and may even result in an infeasible solution given strict flight level restrictions at part of the generated lateral trajectory. In this study, an integrated 4D trajectory optimization problem is tackled to provide cost-effective airline reference trajectories. Furthermore, the promotion of new operation concepts like FRA poses new requirements to airlines’ ability to plan trajectory through a dense airspace region that is computationally difficult to optimize given more flexible flying patterns. Therefore, the transformation from conventional flight operations to TBO certainly brings opportunities and challenges for flight trajectory optimization.To resolve the resulting computation challenges and simultaneously optimize the lateral and vertical dimensions for all three flight phases, (step) climb, cruise, and descend, along with airspeed decisions, this study proposes a graph traversal algorithm based on the A* heuristic, which is an informed search algorithm that finds the shortest path from the source node to target node based on the cost of the path and an estimated cost from every node to the target node. In our implementation, the A* algorithm is adjusted based on the multilabel shortest path algorithm [34] to handle dynamic aircraft status and operation restrictions. In this way, not only are wind, temperature, and an aircraft numerical performance model taken into account, but additional detailed operational restrictions on airway flight levels, restricted areas, Extended-range Twin-engine Operational Performance Standards (ETOPS) [35], FRA, and compulsory routes are considered as well. In particular, compulsory routes are airway segments that are available depending on the availability of other ATS route portions [36]. Besides, a backward Dijkstra algorithm is adopted to estimate accurate heuristic values from each waypoint to the destination airport (target node) based on a machine-learning-assisted regression model for aircraft fuel consumption. The search space is further reduced through a vertex pruning heuristic that utilizes the information provided by route optimization. Three major contributions of this study are summarized as follows: An integrated A* optimization algorithm is proposed for the realistic aircraft reference trajectory optimization problem. This algorithm outperforms the widely applied two-stage method, which involves optimizing lateral and vertical trajectories separately, and incorporates weather, aircraft performance, and complex restrictions. The proposed algorithm is capable of producing high-quality global feasible solutions within a reasonable computation time and has been validated on both the European and the U.S. aviation networks.Two domain knowledge-based heuristics are developed to enhance the convergence of the algorithm. To begin with, machine-learning-based regression models are utilized for each flight phase and aircraft type based on a comprehensive set of attributes such as wind, temperature, and current aircraft weight. Secondly, a vertex filtering heuristic that restricts the search space based on trajectories generated by lateral optimization is proposed. By integrating these two heuristics in a backward Dijkstra algorithm, reliable heuristic costs from each vertex to the destination airport are estimated to enhance the integrated A* algorithm.The proposed algorithm is validated with real-world scenarios, including realistic weather, aircraft performance, airway restrictions, and FRA data. Compared with the commonly used two-stage solution algorithm, the proposed algorithm is able to compute high-quality trajectory solutions within 7 minutes, which is expected to result in substantial annual savings for airlines on fuel consumption and flight time.The remainder of this paper is organized as follows. Section II introduces the underlying airspace network and a generic mixed-integer nonlinear model. Then, algorithm implementation and improvement techniques are presented in Sec. III. Subsequently, an illustrative example and the computational results are presented in Sec. IV. Finally, Sec. V concludes the work and proposes future research directions.II. Mathematical FormulationThe objective of the 4D flight trajectory optimization problem addressed in our study is to determine a series of (waypoint coordinate, altitude, speed, estimated time of arrival [ETA]) tuples for a single flight given its origin/destination airport, scheduled departure time, and the aircraft type in order to minimize the weighted cost. This cost is composed of fuel consumption (measured in kg) and the product of flight time (measured in min) and cost index coefficient CI (measured in kg/min), which trades off flight time and fuel cost.A. Underlying Airspace NetworkTo support the modeling of the feasible flight trajectories, an airspace network G(V,E) is introduced, which is a multi-edge directed graph built based on air navigation data. The vertices v∈V represent either airports or (Enroute, Navaids) waypoints with ID, coordinates, and elevation. Significant points and coordinates from FRA with a horizontal separation of 1° are also included to allow flexible flight plans and reduced route extension. An edge e∈E between two vertices v1 and v2 is constructed in four cases: 1) v1 and v2 are two endpoints of an airway segment; 2) v2 can be reached by v1 via direct flight (DCT, direct connection with great circle); 3) two points belong to a portion of a procedure; 4) v2 lies in the neighborhood of v1 (e.g., 10 nearest points) within the same FRA and meet specific requirements (e.g., not crossing boundaries). The possible procedures include standard instrumental departure procedure (SID), standard terminal arrival route (STAR), and compulsory routes. Each procedure contains one or more route portions (i.e., edges) and may require successive execution of the involved edges with specific speed and altitude restrictions.Due to heterogeneous weather conditions (e.g., wind, ISA deviation) and operational constraints (e.g., minimum sector altitude MSA, restricted areas), the connectivity of edges varies for different flight levels and regions. An intuitive example is demonstrated in Fig. 1a, which illustrates the main modeling elements of the studied problem. A simulated airspace network is presented that includes two airports, each with a line-by-sight distance of 400 km, 13 waypoints, and two flight levels (FLs). The origin and destination airports are connected with other transition vertices through SIDs or STARs, which may consist of multiple edges, respectively. Additionally, the compulsory route from WP 2 to WP 6 requires that the edge between WP 2 and WP 4 must be visited in advance before traversing the edge between WP 4 and WP 6. We have also defined an FRA, which composes of one entry point WP 6 and two exit points WP 9/WP 10. Within the polyhedron, aircraft can fly freely without considering the dependency on predefined waypoints. Furthermore, altitude restrictions on the transition waypoints WP 1 and WP 13 only allow their usage at FL 410 (41,000 ft) with ISA assumed, while WP 12 is only available from FL 390 to FL 410. Figure 1b exemplifies the FRA concept with several significant points for entry (E), exit (X), intermediate (I), arrival connecting (A), and departure connecting (D), which are represented by triangles. Different types of edges are created according to the related rules and highlighted with corresponding colors. It is worth noticing that the complete network is constructed not only on the great circle curve from the origin airport to the destination airport but also on the aircraft’s maximum cruise range (e.g., ETOPs) and related operational time restrictions on airways/waypoints during the planning horizon. Therefore, the overall network size is relatively larger than the typical practices in the literature that sample vertices along the great circle curves.Fig. 1 The airspace network structure of an illustrative example.B. Generic Mathematical FormulationBefore introducing the customized solution algorithms for the 4D flight trajectory planning problem, we depict the corresponding constraints and dependencies with a mixed integer nonlinear mathematical model. To start with, as lateral and vertical profiles are jointly considered, a node set N is defined where each node is a combination of vertex v∈V and flight level FL. In particular, node S and T denote the origin and destination airport, respectively. Likewise, an arc set A contains node pairs that can be successively flown by aircraft. For an arc between node i(vi,FLi) and node j(vj,FLj), it is feasible if (vi,vj)∈E and |FLi−FLj| must be less than the upper bound derived from the aircraft maximum climb/descend rate. A subset of arcs belonging to one procedure (i.e., SIDs/STARs/compulsory routes) is represented by Ac.Further define nonlinear functions FT and FC to compute the flying time and fuel consumption of arc (i, j) given the weather condition (ISA deviation, temperature, wind), flight level, and the performance profile. Based on the relevant parameters and decision variables defined in Table 2, where restrictions on airway segments can be derived from Aeronautical Information Publication (AIP), we model the 4D flight trajectory optimization problem as follows: min(fT−fS)+CI⋅(tT−tS)/60(1)∑j∈δi+xij−∑j∈δi−xji=Di∀i∈N(2)∑j∈δS+xSj=1(3)tj=∑i∈δj−ti+xij⋅FT(i,j,ti,vij)∀j∈N(4)fj=∑i∈δj−fi−xij⋅FC(i,j,ti,vij)∀j∈N(5)∑j∈δi+xij≤max{0,fi−FRF−ZFW}∀i∈N(6)xij≤max{0,min{ti−ot_ij,o¯tij−tj}}∀(i,j)∈E(7)xij≤xkl∀(i,j)∈Ac,(k,l)∈θij(8)fi≤MTFW∑j∈δi+xij∀i∈N(9)ti≤T¯∑j∈δi+xij∀i∈N(10)The objective function minimizes the cost of fuel and flying time given a user-specified cost index where the differences between the target node T and the source node S are computed. Flow balance is ensured at every intermediate node as specified by constraints (2). Constraint (3) ensures that the aircraft departs from the origin airport. Constraints (4) and (5) model the vertex arrival time and remaining mass weight given the selected edges. Aircraft are not allowed to fly when the remaining fuel weight is less than the FRF as specified by constraints (6) and the scheduled arrival time at each edge with opening time restrictions is ensured by constraints (7). Constraints (8) require the consecutive execution of every procedure. Finally, when a node is not selected, the corresponding mass weight and flight time are forced to take zero according to constraints (9) and (10).From the model formulation, it is evident that the model depends heavily on the possible combinations of waypoints and flight levels, which drastically increase along with the instance size. Variable xij implicitly requires the satisfaction of both operational and aircraft-performance-related constraints (such as airspace closure, and maximum climb rates). On the other hand, the value of generic function FC(⋅) can be set to an arbitrarily large number to deter the selection of infeasible variables, thereby preventing violations of the constraints. Moreover, determining the flight time and fuel consumption on a specific edge involves making decisions on airspeed, climb/descend rate (angle), and the experienced flight phases (e.g., climb then cruise) simultaneously. This, coupled with the large solution space and a high degree of nonlinearity, poses significant challenges in efficiently solving the 4D flight trajectory optimization problem.III. Optimization AlgorithmsIn this section, multiple graph search solution algorithms and acceleration techniques are presented to tackle the flight trajectory optimization problem in order to minimize fuel consumption and flight time.A. Integrated A* AlgorithmThe 4D flight trajectory optimization problem is heuristically solved using a revised A* algorithm given that a large, complex solution space, and the consideration of realistic operational restrictions makes the problem computationally intractable. The principle of optimality, which is commonly adopted in vast literature for shortest path algorithms, does not always hold due to the following reasons. Firstly, the consideration of the remaining fuel and aircraft status leads to a resource-constrained shortest path problem. Furthermore, complex dependent restrictions, e.g., compulsory routes, require the execution of backtracking for previous flight paths and states. In addition, the inclusion of vertical profile optimization where cruise altitudes typically range from FL 100 to FL 410 significantly increases the graph size and computational complexity. Therefore, the shortest subpath at a particular vertex is not necessarily a part of the complete shortest path, and applying the traditional multilabel shortest path algorithm leads to a large number of labels being generated given their cardinality.Hence, the underlying optimization problem is inherently hard, and we propose a variant of the A* algorithm to solve the problem. Different from exact graph search methods that thoroughly search the complete graph depending on the actual cost c defined in Eq. (1), A* traverses graph vertices according to both c and a heuristic cost h from this vertex to the target node to guide the search direction. With the heuristic approximation of the cost to the destination, the A* algorithm works similarly to a greedy best-first search procedure. It is important to note that the heuristic values h are critical to the efficiency and applicability of the algorithm. Given the aforementioned problem properties, the A* algorithm is combined with the multilabel shortest path algorithm as shown in Algorithm 1, where Extend is a function to generate labels and compute cost c.In Algorithm 1, a priority queue Q is maintained to pop up labels with the minimum cost that is composed of actual cost c and a heuristic value h. In addition, the label contains other attributes such as flight level, flight time, flight phase, aircraft mass, and entered airways/areas. The algorithm terminates when encountering a label at the target node for the first time or once the priority queue is empty. Consequently, the algorithm can derive an optimal solution if h(l) is a valid lower bound to c(l). Starting from the origin airport node S, the algorithm visits and extends labels along the edges until the destination airport node T. When the initial takeoff weight (TOW) is not provided, it is approximated using the Breguet range Eq. (11): Winitial=Wfinal⋅exp(Range⋅g⋅SFCV⋅(L/D))(11)where Range is calculated from the great circle distance from the origin airport to the destination airport and g represents the gravity of Earth. As the quotient of lift (L) and drag (D) is equivalent to the division of lift coefficient CL and CD, these two parameters along with specific fuel consumption (SFC) and flight velocity (V) are associated with individual aircraft. To account for the reversed fuel weight, Winitial can be estimated as 5% more than zero fuel weight (ZFW).For each label retrieved from the priority queue Q, it is extended toward every successor vertex v2 to generate a set of new labels given all possible decisions (climb, cruise, and descend). If the new label features less cost than one existing label l′ for the same flight level from Bv2, l′ is discarded and substituted by l2. The aircraft trajectory is finally generated from the best label in BT. The Extend function in line 10 of Algorithm 1 creates new labels by extending the existing labels through edges. The algorithm first checks the restrictions on airway segment (v,v2) and vertex v2 in terms of available altitudes, aircraft turning angle, requirements posed by compulsory routes, and FRA. Subsequently, cruise/climb/descend with different aircraft performance schedules are explored with different fuel costs and flight time. If the corresponding remaining fuel weight is less than the required minimum reserved fuel weight, the schedule is discarded. To capture the dynamic aircraft mass and weather condition, the edge (v,v2) is discretized during the extension according to the applied wind grid resolution (50 nm), and the fuel cost is computed iteratively from tail vertex v to the head vertex v2 with necessary linear interpolation. Since the desired flight level during climb and descent can be reached before arriving at vertex v2, the edge segment after reaching the desired altitude is flown with a cruise phase.Finally, to solve the problem efficiently, a high-quality heuristic h(l) is required. In previous studies, great-circle-distance-based estimation is applied to approximate the fuel consumption from every vertex to the final destination vertex T [28]. In our study, we first use a similar heuristic function h(l)=gcd×β, where gcd represents the great circle distance between two vertices and β is the approximated cost per nautical miles as a comparison method. Regarding the illustrative example, assuming that the distance between Airport A and Airport B is 400 km, using the performance data from the CAE Flight Plan Manager‡‡ of a B787 aircraft, the estimated consumed fuel is computed as Winitial−Wfinal=3777.54 kg. To simplify the computation process, the wind is not taken into account. During the algorithm solution process, every vertex is explored at different altitudes, and one feasible route consists of Airport A→SID2→WP 1→WP 3→WP 5→WP 7→WP 11→WP 13→STAR 2→Airport B.B. Two-Stage AlgorithmDifferent from an integrated solution algorithm that jointly optimizes lateral and vertical profiles, most existing studies apply a sequential approach that optimizes these two profiles separately, resulting in a significant reduction in computational complexity. In this context, we present a two-stage algorithm to solve the lateral and vertical optimization problems that can produce suboptimal solutions within a short computation time.To be specific, the algorithm terminates only when the priority queue is empty, thereby eliminating the need for lines 6 and 7 in Algorithm 1. Comparisons of the generated labels with existing labels are performed for dominance check. Label l1 is weakly dominated by label l2 if c(l1)≥c(l2) and all attributes like flight level, flight phase, and performance schedule remain the same. For dependent restrictions like SID/STAR and compulsory route, additional equality conditions are also needed during the extension.Since altitude decisions are not taken into account during lateral optimization, it is also essential to modify the Extend function regarding fuel consumption and different flight phases. Specifically, high flight levels are preferred in the function if available (considering altitude restrictions and aircraft performance) to reduce the fuel cost with less air drag and low air density at high altitudes since only the cruise phase is considered. Dummy edges that connect the airports with their SID/STAR transition points are created as well in order to avoid inaccuracy in selecting SIDs/STARs using cruise performance profiles. The complete solution procedures are shown in Algorithm 2. If the 2D trajectory obtained in lateral optimization can be further transformed into a 3D trajectory after vertical profile optimization satisfying altitude restrictions, it can be utilized as the output flight trajectory solution. In contrast, the last valid edge in path p is identified and deactivated to ensure its successive invalid edge during the re-optimization phase. For instance, due to the lack of consideration for altitude limit at WP 12, the resultant trajectory, i.e., Airport A,SID1,WP 1,SID1,WP 2→WP 4→WP 6→WP 10→WP 12→STAR1→Airport B, obtained through lateral optimization is not feasible for the vertical optimization phase. The minimum flight level of FL 390 is too high for the aircraft to join STAR1. To resolve this issue, it is necessary to avoid edge WP 10→WP 12 during the subsequent round of optimization, which results in a slight detour from WP 4 to WP 11.C. Acceleration TechniquesAlthough the presented integrated A* algorithm is capable of ob