Predicting Fuel Consumption Reduction Potentials Based on 4D Trajectory Optimization with Heterogeneous Constraints

: Investigating potential ways to improve fuel efﬁciency of aircraft operations is crucial for the development of the global air trafﬁc management (ATM) performance target. The implementation of trajectory-based operations (TBOs) will play a major role in enhancing the predictability of air trafﬁc and ﬂight efﬁciency. TBO also provides new means for aircraft to save energy and reduce emissions. By comprehensively considering aircraft dynamics, available route limitations, sector capacity constraints, and air trafﬁc control restrictions on altitude and speed, a “runway-to-runway” four-dimensional trajectory multi-objective planning method under loose-to-tight heterogeneous constraints is proposed in this paper. Taking the Shanghai–Beijing city pair as an example, the upper bounds of the Pareto front describing potential fuel consumption reduction under the inﬂuence of ﬂight time were determined under different airspace rigidities, such as different ideal and realistic operating environments, as well as ﬁxed and optional routes. In the congestion-free scenario with ﬁxed route, the upper bounds on fuel consumption reduction range from 3.36% to 13.38% under different benchmarks. In the capacity-constrained scenario, the trade-off solutions of trajectory optimization are compressed due to limited available entry time slots of congested sectors. The results show that more ﬂexible route options improve fuel-saving potentials up to 8.99%. In addition, the sensitivity analysis further illustrated the pattern of how optimal solutions evolved with congested locations and severity. The outcome of this paper would provide a preliminary framework for predicting and evaluating fuel efﬁciency improvement potentials in TBOs, which is meaningful for setting performance targets of green ATM systems. This provides a preliminary but adaptable framework for predicting and evaluating city-pair fuel efﬁciency improvement potentials in the TBO context. A more comprehensive study that expands the research scope from a single city pair to the entire airport network including domestic and international ﬂights would achieve more insights for setting green development goals in future ATM. investigating


Introduction
Growth in air travel is having an increasing environmental impact. According to the 2019 International Civil Aviation Organization (ICAO) Environmental Report, one of the aviation industry's targets to improve fuel efficiency by 1.5% every year has been constantly met and exceeded, averaging a 2.4% annualized improvement since 2009. The aim is to halve net emissions by 2050 compared with 2005 [1]. Studies have shown that there is a significant proportional relationship between carbon emission and fuel consumption [2]. Therefore, with the continuous growth of demand for air transportation, increasing the fuel efficiency of aircraft has become a key goal in the aviation industry. Aircraft technology and design improvements [3,4], improvement of the socioeconomic policies related to air transportation [5,6], optimization of aircraft operation processes [7,8], and the search for alternative fuels [9,10]  Trajectory-based operation is one of the basic characteristics of the new generation of air transportation systems. It aims to realize four-dimensional (4D) trajectory planning and control before takeoff and during flight through continuous information exchange between various stakeholders. The realization not only plays a major role in enhancing the predictability of air traffic and flight efficiency but also provides a new means for saving energy and reducing emissions. The 4D trajectory optimization model is a process of designing continuous or discrete trajectories composed of 3D positions and time, which minimizes (or maximizes) some measure of performance (travel time, environmental impact, etc.) while satisfying a set of constraints (aircraft performance, capacity, etc.). At present, trajectory optimizations for different flight phases, including taxiing on the airport surface, climb and descent in the Terminal Maneuvering Area (TMA), and cruise in en-route airspace, aiming at improving fuel efficiency have been extensively studied.
In order to reduce fuel consumption and exhaust emissions during taxiing on the airport surface, push-back control strategy [11,12], taxiing trajectory and speed profile optimizations [13,14], and improvement of taxiing mode (e.g., single-engine taxiing, operational tow-outs, and electric taxiing) [15,16] have been comprehensively investigated. During climb and descent, especially in high-density terminal airspace, prior research has focused more on how to make full use of runway capacity, while minimizing noise, total fuel consumption, emissions, and delays [17][18][19], and put forward operating concepts and operating procedures such as Continuous Climbing Operation (CCO) [20], Continuous Descending Operation (CDO) [21], Point Merge (PM) [22], etc. Prats et al. [23] presented a non-linear multi-objective optimal control methodology to design noise abatement procedures aimed at reducing the global annoyance perceived by the population living around the airports. Zhu et al. [24] generated the optimal descent trajectory based on the verification that there existed a certain relationship between fuel consumption and flight altitude as well as between aircraft weight and vacuum speed during descent. Rodríguez-Díaz et al. [25] developed a multi-objective optimization model under the Constrained Position Shifting (CPS) constraint to minimize noise pollution, fuel consumption, and delays. In addition, with the continuous deepening of research and field trials, the traditional terminal standard arrival and departure procedures are shown to be relatively rigid, which are less flexible to support optimal trajectory planning and implementation in the TBO context. Some novel airspace structure and operation modes were proposed to enhance the flexibility and controllability of the user-preferred trajectory [26,27].
For trajectory optimization in the en-route airspace, considerations are mainly given to fuel consumption, emissions, condensation trails, and their impact on the climate. Lovegren et al. [28][29][30] improved the fuel efficiency of aircraft by optimizing the altitude and speed profile. Soler et al. [31] studied the influence of atmospheric uncertainty, convective indicators, and cost-index on the leveled aircraft trajectory optimization. For the horizontal profile, Patron et al. [32] considered the wind direction during cruise and used genetic algorithms to select an optimal flight path for the aircraft in the pre-tactical stage. With the objective of minimizing total fuel consumption, Hartjes et al. [33] studied the long-haul flight fuel efficiency of commercial aircraft formations and analyzed the impact of wind field on the performance. Studies have shown that the choice of flight path has a significant impact on fuel efficiency and the environment. Coordinated optimization of a flight's path and vertical profile has gradually become popular in recent years [34,35].
In addition, integrated trajectory optimization during multiple flight stage have been studied. Soler et al. proposed a multi-phase optimal control framework for commercial aircraft 4D flight trajectory planning [36]. Taking conventional cost and environmental cost Sustainability 2021, 13, 7043 3 of 33 into account, Tian et al. optimized a green 4D trajectory of single commercial aircraft [37]. Based on alternative trajectory options and traffic volume hotspot detection, Melgosa et al. presented an enhanced Demand and Capacity Balance (DCB) formulation, with the aim of minimizing operating costs [38]. Harada et al. evaluated the operational advantages of TBO in terms of fuel consumption, flight time, and flight distance of domestic flights in Japan [39].
Extensive studies on green trajectory planning have been implemented for different flight phases (takeoff, landing, and cruise) at the strategic and pre-tactical levels, but there are still some limitations. (1) Trajectory planning is more focused on individual optimization of the flight stage, while complete trajectory optimization from runway to runway, or even gate to gate, is rarely discussed. Meanwhile, trajectory optimization is more focused on a single dimension of a horizontal path or a vertical profile. (2) In the process of trajectory optimization, less consideration is given to realistic constraints such as the candidate routes and sector capacity limitations. (3) Investigating the potentials of fuel efficiency improvement trading off against travel time through 4D trajectory optimization with different constraints is rarely reported. In particular, how the fuel-saving potentials evolve with the congestion uncertainties has not been revealed, such that the emission reduction targets of air traffic management planning lack solid support. For this reason, based on the existing research results, taking Shanghai Hongqiao Airport (ZSSS)-Beijing Capital International Airport (ZBAA) as an example, and considering real meteorology data, aircraft dynamics, available route limitations, sector capacity constraints, and air traffic control restrictions on aircraft altitude and speed, a "runway-to-runway" multiobjective planning model under loose-to-tight constraints was developed in this study. The goal is to investigate potential fuel consumption and flight time reductions under different levels of operation restrictions, such as free-flow and capacity-constrained airspace sectors, as well as fixed and optional routes. The model enables the formulation of green aviation development goals.
The main work of this study includes: • A runway-to-runway 4D trajectory model and fuel consumption model were established, and a fixed-path multi-objective trajectory optimization model based on trajectory discretization was proposed. The NSGA-II algorithm was used to search the optimal speed-altitude profile between city pairs. The upper bound of fuel efficiency was investigated by examining the trade-offs between fuel consumption and flight time.

•
In order to cope with uncertainties in sector-level demand and capacity balancing that often appear in real-world operation, two 4D multi-objective optimization models based on fixed and optional routes in capacity-constrained scenarios were established, respectively, by taking into account the variation of available sector entry slots changed after takeoff. The impact of airspace congestion on fuel-saving was modeled and measured. • Further, a sensitivity analysis underpinned by a metric for evaluating the effect of multi-objective optimization was supplemented to give insight on how the fuel-saving potentials, as well as Pareto optimal front evolve with the congestion location (e.g., distance to the departure airport) and severity (e.g., level of delay).
The rest of this paper is organized as follows. The problem description is presented in Section 2. The 4D trajectory model and fuel consumption estimation model are presented in Section 3. An ideal trajectory optimization model for fixed-path scenario and two models that take into consideration capacity constraints for fixed and optional paths are presented in Section 4. The Non-dominated Sorting Genetic Algorithm with Elitist Strategy (NSGA-II) is adaptively improved for determining a multi-objective solution of fuel consumption and flight time in Section 5. The flight between ZSSS and ZBAA is taken as an example in Section 6 to verify the accuracy of the fuel consumption model and explore the effect of trajectory optimization on fuel consumption improvement under different strictness constraints. Conclusions and future research directions are indicated in Section 7.

Problem Description
In addition to being affected by its own aerodynamics and meteorological conditions, aircraft fuel consumption has significant sensitivity to flight trajectory, such as altitude, speed, and flight distance. In the TBO concept, although Shared Business Trajectory (SBT) would be coordinatively planned before takeoff, based on information such as weather conditions, scheduled flight time, scheduled route, and airspace capacity, during flight, the SBT is often difficult to implement due to disturbances of external dynamic factors such as airspace availability changes and air traffic controller intervention. The main problem to be solved in this paper is that of predicting the possible upper bound of fuel consumption reduction for a given city pair using trajectory optimization with various constraints without changing the airspace structure.
To address this problem, three multi-objective mixed-integer program models were formulated in this study while considering fuel consumption and flight time for a city pair: an ideal 4D trajectory optimization model and trajectory optimization models that considered the constraint of a sector's controlled time of arrival (CTA) for fixed and optional paths, respectively. These models were used to determine potential fuel consumption reductions in different scenarios. To simplify the problem, the vertical profile was divided into three stages: climb, cruise, and descent. The horizontal trajectory of the aircraft was discretized in 1 km increments. The time-continuous trajectory optimization problem was transformed into a Mixed-Integer Nonlinear Program (MINLP) problem, resulting in a limited set of decision variables.
For fixed-path scenarios, under various operational constraints such as aircraft properties, control rules, and flight comfort, the relationship between flight speed and altitude was optimized collaboratively to generate a speed and vertical profile that weighs the dual goals of fuel consumption and flight time. When optional routes are considered, the route network was abstracted as a network G = (V, E), where V = {j 1 , j 2 , . . . , j n } is the set of nodes, j n is the starting point for each flight segment, E = {k 1 , k 2 , . . . , k m } is the set of flight segments, and the length of each flight segment is set to 1 km. Let two sequential nodes be j n and j n+1 that are members of an adjacent matrix Θ = e[j n , j n+1 ]. If j n → j n+1 is connected, then e[j n , j n+1 ] = 1. Otherwise, e[j n , j n+1 ] = +∞. At this time, the path of aircraft α can be expressed as the set of nodes connected in series. This problem can be expressed as: under the constraints of connectivity in the network G, various operational constraints are integrated to plan the aircraft's flight altitude and speed profile simultaneously. By establishing a Pareto frontier for fuel consumption and flight time in various scenarios, the potential fuel consumption reduction can be identified.

Meteorology Model
For aircraft trajectory optimization research, the main atmospheric properties are temperature, pressure, and the density of the air. In Base of Aircraft Data (BADA) 3.11 published by EUROCONTROL [40], the relations among the atmospheric temperature, air density, and pressure are modeled as Formula (1): where T em is the atmospheric temperature at the corresponding altitude; T 0 = 288.15 (K) is the mean sea level standard atmospheric temperature; h is altitude (m); the altitude of the tropopause is h trop = 11,000 (m) in the ISA model and the temperature is con-sidered constant above the tropopause; T trop is the temperature of the tropopause (K); β = −0.0065 (K/m) is the ISA temperature gradient below the tropopause; P is the atmospheric pressure; P 0 = 101,325 (Pa) is the mean sea level standard pressure; g is the local gravity acceleration; R = 287.05287 (m 2 / K · s 2 ) is the gas constant; and ρ is the atmospheric density. The data of weather in this article are mainly derived from the GRIB2 files published by Global Forecast System (GFS) of the National Oceanic and Atmospheric Administration (NOAA), and the GFS data are further divided into GFS Analysis data and GFS Forecasts data. In this article, the weather data of the 0.5 * 0.5 latitude and longitude grid are selected in the GFS Analysis data.

Aircraft Dynamic Model
To accurately assess fuel consumption and flight time, in this paper, based on the Point Mass Model (PMM) [41] and the Total-Energy Model (TEM) [40], as well as adaptive improvement based on the problem, the aircraft dynamic model is established to simulate aircraft flight.
According to PMM, assuming that all forces act on the center of gravity of an aircraft, the motion of the aircraft is reduced to three degrees of freedom (three translations), and the equations of motion in three dimensions can be integrated directly. The forces acting on the aircraft are shown in Figure 1. In Figure 1, γ is the flight path angle; α is the flight angle of attack; ϕ T is the engine installation angle (usually very small and approximately ϕ T = 0 ); v is the true airspeed (TAS); G is the force of gravity; and L and D are the lift and drag on the aircraft, respectively, expressed as: where S is the reference wing area; C L and C D are the lift and drag coefficients, respectively, expressed as: φ is the slope of the aircraft; C D0 is the parasitic drag coefficient; and C D2 is the induced drag coefficient.
In an air reference system with wind component, the aircraft dynamics model is described as follows: where the state vector u = [s, h, v, m] is formed, respectively, by the position s along path distance (km), altitude h (m), true airspeed v (m/s), and mass m of the aircraft (kg). The control vector u = [T, γ] is composed of the thrust T on the aircraft and the flight path angle γ; v wind represents wind speed (m/s); and FF is the fuel flow of the aircraft (kg/s). In this article, after discretizing the trajectory according to the voyage, in order to facilitate calculation, we assume that, on flight segment i n → i n+1 , the T, γ, and the rate of climb/descent (Rate of Climb or Decent (ROCD)) remain unchanged. The TEM equates the rate of work done by forces acting on the aircraft to the rate of increase in potential and kinetic energy, namely:

Aircraft Performance Model
Aircraft performance is the parameters that describe the motion law for the aircraft's center of mass, including the speed, altitude, thrust, and so on [42]. The models and parameters of aircraft performance used in this paper are from BADA 3.11.
In the process of simulating aircraft movement, the acquisition of thrust is very important. This paper transforms the calculation of thrust into the solution of acceleration. That is, when ROCD i n ,i n+1 and v i are determined, the acceleration of the aircraft a i n ,i n+1 on the flight segment i n → i n+1 can be estimated according to Formula (8), and the thrust of the aircraft can be calculated in combination with Formula (7).
where d is the distance of the flight segment i n → i n+1 ; t i n ,i n+1 is the flight time of the flight segment i n → i n+1 ; and the corresponding ROCD i n ,i n+1 feasibility of different flight altitudes are determined according to the speed-altitude profile recommended by PTF files of BADA 3.11. The speed of the aircraft is an important parameter for aircraft performance as well. According to the BADA 3.11 OPF file, the calibrated air speed (CAS) is often used to determine whether the aircraft speed exceeds the operational threshold. The CAS is calculated with Formula (9) : where ρ 0 = 1.225 (kg/m 3 ) is the density values at sea level; p is the atmospheric pressure; and κ = 1.4 is the specific heat ratio of the air. According to Formula (8), the ground speed (GS) v G can be used to calculate the γ i n ,i n+1 and a i n ,i n+1 . The GS is obtained after the correction of the v wind for TAS:

Fuel Consumption Model of Aircraft
When calculating fuel consumption, this article adopts the fuel flow rate model of BADA 3.11 and makes adaptive improvements to calculate the fuel consumption of aircraft from takeoff to landing. According to the BADA 3.11, a jet aircraft fuel consumption can be calculated with the following equations: where η is the thrust specific fuel flow; C f 1 and C f 2 are, respectively, the first and second thrust-specific fuel consumption coefficients; and FF is the instantaneous fuel flow, and the total fuel consumed from 0 to t is F. Formula (11) is applicable at all flight stages except during the cruise phase and idle thrust. Idle thrust is not considered in this paper, while FF in the cruise phase is defined in Formula (12): where C f cr is the fuel flow correction coefficient in cruise phase. After discretization, the total fuel consumption can be approximated as a sum of discrete values of fuel consumption in each segment throughout the flight as following: An estimation is done for the value of the aircraft initial state, based on historical data. After the ROCD of each segment is determined, using Formulas (1)-(10), the state vector, control vector in any segment, and arrival time at each node during the flight can be estimated, thus the average fuel flow rate FF α i n ,i n+1 can be calculated using Formulas (11) and (12).

Basic Assumptions
To simplify the problem, the following basic assumptions were made: • The clean configuration of the aircraft is adopted for optimization. • Human interventions are not included in the flight trajectories.

Decision Variables
In this paper, the horizontal trajectory is discretized, and the continuous 4D trajectory optimization is transformed into a problem involving determining the flight speed and altitude at the nodes in each flight segment. Therefore, the trajectory optimization decision variables under a fixed path are the speed v α i and altitude h α i of the aircraft at node i.

Objective Function
Model I is intended to generate an ideal economic flight trajectory for the aircraft based on the aircraft performance model proposed in Section 3 and a pre-set horizontal path, while meeting the basic air traffic operation restrictions. The total fuel consumed by the aircraft in the air is the main optimization objective of Model I. In addition, the flight time between city pairs is also an important factor that airlines must consider when making flight plans and crew scheduling. Therefore, this paper establishes the dual objectives of fuel consumption and flight time to examine the trade-off relationship between the two, as well as to analyze the influence of flight time adjustments on the potential of fuel consumption improvement. The objective function is defined in Formula (14): where α is the aircraft; A is the set of aircrafts; i is a node; I is the set of nodes; and F α , respectively, represent the fuel consumption and flight time of aircraft α on the flight segment.

•
Altitude constraint The aircraft needs to maintain a specific altitude during flight. This limitation is determined by the aircraft's own properties and airspace limitation, as defined in Formula (15).
where h α max is the maximum altitude of the aircraft under the constraint of its own properties and h α i,control min and h α i,control max are the altitude limits in the controlled sector, i.e., the minimum and maximum altitude thresholds of aircraft α at node i. In particular, if the flight is a two-way route, the aircraft should comply with the "odd for eastbound and even for westbound" regulations in addition to meeting the above constraints, i.e., an aircraft flying eastward adopts odd-numbered altitude levels and an aircraft flying westward adopts even-numbered altitude levels.
• Speed constraint Similarly, the speed of aircraft α at node i must meet the limitations of its own properties and control rules of the sector.
where v α i,min and v α i,max are the minimum and maximum speed of aircraft α at node i under the constraint of its own properties and v α i,control min and v α i,control max are the minimum and maximum speed of aircraft α at node i under the air traffic control constraint in a specific sector. The thrust required to operate the aircraft is also taken as a constraint to ensure passengers' comfort: where a = 0.6096 m/s 2 is the maximum longitudinal acceleration recommended in the BADA 3.11 GPF file. • Climb and descent performance constraints The rates of climb and descent of the aircraft cannot exceed the maximum rate of climb or the maximum rate of descent. In this article, the rate of climb and the rate of descent introduced from BADA 3.11 PTF file are used as the thresholds for the climb and descent operations.
where CL α max and DE α max are the aircraft's maximum rate of climb and maximum rate of descent, respectively. State variables x α i n ,i n+1 and y α i n ,i n+1 are introduced to denote the flight state, as defined in Formula (20).
• Timeliness constraint Time efficiency is another important dimension from airline's perspective. Therefore, when planning a trajectory, the city-pair flight time should be kept within an acceptable range, as shown in Formula (21).
where SAT α is the standard airborne time (SAT) of aircraft α, which is obtained according to flight schedule and the Collaborative Decision Making (CDM) system in China; MATD is a positive number, which is the maximum acceptable time delay (MATD) of landing; and MATA is a negative number, which is the maximum acceptable time advance (MATA) of landing.

Aircraft Flight Trajectory Optimization Model for the Fixed Path Scenario under Capacity Limitation (Model II)
The CDM system is a regional air traffic flow management system that allocates flight departure slots by integrating various air and ground restrictions. However, due to the takeoff time slot freeze (generally 30 min in advance), changes in the sector capacity-flow balance may occur due to uncertainty, resulting in unavailability of an entry time slot in the original schedule. In the TBO context, the aircraft would be capable of re-planning flight trajectories based on air-ground negotiation to generate an updated available entry time slot in the restricted sector. Therefore, based on the ideal horizontal path in Model I, a 4D trajectory multi-objective optimization model for a fixed path with capacity constraints was developed to analyze the impact of airspace congestion on potential reductions in fuel efficiency for a rigid air route structure.

Objective Function
Model II is constructed based on Model I, and the meaning and expression of the objective function are the same as that in Formula (14).

Constraints
The meanings and expressions of altitude constraints, speed constraints, climb and descent constraints, and time constraints in Model II are the same as those in Model I, as shown in Formulas (15)- (21). Due to capacity limitation, this model shall guarantee that all the aircraft enter the sectors within available time slots, as shown in Formula (22).
where i is a node that enters sector w; τ α i is the time aircraft α passes node i; T w is the set of available time slots of sector w; and W is the set of sectors.

Multi-Path-Based Aircraft Flight Trajectory Optimization Model under Capacity Constraints Scenario (Model III)
Compared with Model II, Model III relaxes the rigid route structure restrictions, allowing the aircraft to independently select flight paths based on an updated sector capacity-flow relationship after takeoff. This allows the effect of airspace flexibility on aircraft fuel efficiency to be to explored.

Decision Variables
v α j and h α j have the same meaning as in Model I. To identify whether aircraft α passes through node j, the variable z α j is introduced, as shown in Formula (23).
Furthermore, to identify the connectivity between nodes that the aircraft passes, the variable γ α j n ,j n+1 is introduced to indicate whether the aircraft chooses to fly through the edge formed by connected nodes. For any two nodes j n and j n+1 , there is:

Objective Function
Model III allows the aircraft to flexibly select the flight path according to sector congestion information after takeoff. This provides integrated planning of the aircraft's horizontal trajectory, vertical profile, and flight speed, which are used to minimize flight time and fuel consumption, as shown in Formula (25).
where C(j) is the node connected with node j.

Constraints
The Model III constraint condition changes and expands the Model II constraint based on network connectivity. For ∀α ∈ A, ∀j ∈ V, ∀w ∈ W, the constraint conditions are shown in Formula (26).
Formula (26a) is the connectivity constraint, i.e., adjacent nodes in the path selected by the aircraft must be connected in the network. Formula (26b) is a unique constraint, i.e., after the aircraft passes through any node there is one and only one next node. Formulas (26c)-(26g) are the speed, altitude, climb and descent performance, timeliness, and sector available time slot constraints, respectively, which are the same as those defined in Model I.

Solution Algorithm
The complexity of the aircraft trajectory optimization problem increases geometrically as the number of decision variables increases. It is difficult to obtain and verify the solution in polynomial time. As an NP-hard problem, the solutions are independent, therefore a metaheuristic algorithm is used to find the solution [43]. A genetic algorithm is a random search algorithm that mimics natural selection, heredity, and mutation. It provides a general framework for solving complex optimization problems, and it has been successfully used in many disciplines. Research results show that the Non-Dominated Sorting Genetic Algorithm (NSGA-II) using the elite retention strategy can efficiently seek a set of different solutions and converge to the real Pareto optimal solution set. Recently, a comprehensive comparison was made to justify the overall advantages of NSGA-II in solving trajectory optimization problems in terms of solution quality and convergence speed [27,44,45]. Therefore, the NSGA-II algorithm was chosen to solve the optimization problem outlined herein, and a flowchart for this algorithm is shown in Figure 2.  To solve the problem, the NSGA-II algorithm was modified as follows: • Encoding The chromosome encoding pattern is shown in Figure 3. For the fixed path scenario, the altitude and speed of each node from initial node i s to termination node i e are defined as positive integer, as shown in Figure 3a. For the multi-path scenario, to improve the computational efficiency, a hierarchical and variable-length chromosome was designed. The first layer is encoded as a 0-1 matrix presenting the route choice based on the variable γ j n ,j n+1 defined in Section 4.3. The green gene represents the connectivity between node j n and j n+1 according to the adjacent matrix Θ = e[j n , j n+1 ] defined in Section 2. Otherwise, the links are not connected and the genes are encoded as 0, which do not change during subsequent genetic manipulation. The second layer of the chromosome is as the same as the fixed path scenario, as shown in Figure 3b. Note that the length of the second layer is determined by the route coded in the first layer.
was designed. First layer is encoded as a 0-1 matrix presenting the route choice based on the variable γ j n ,j n+1 defined in Section 4.3. The green gene represents the connectivity between node j n and j n+1 according to the adjacent matrix Θ = e[j n , j n+1 ] defined in Section 2. Otherwise, the links are not connected and the genes are encoded as 0, which do not change during subsequent genetic manipulation. Second layer of the chromosome is as the same as the fixed path scenario, as shown in Figure 3b. Noted that the length of the second layer is determined by the route coded in the first layer. • Initial population The initial population, as the iterative starting point of the heuristic algorithm, is an important factor affecting population evolution and efficiency of the algorithm. When optimizing a single path, historical trajectory data was used to form the initial population to obtain the Pareto frontier solution. When optimizing multiple paths, the Pareto frontier solutions for a single path were randomly selected and combined into multiple path optimization individuals, forming the initial population. One should note that in actual operation, the altitude and speed of an aircraft, being significantly sensitive to factors such as maneuvering actions, changes in aerodynamic configuration, and meteorological conditions, are very prone to fluctuations. Therefore, to determine the optimal altitude-velocity profile under ideal conditions, the altitude and velocity data in the historical trajectory data were smoothed [46]. • Genetic manipulation The tournament algorithm and elite retention strategy were used to select individuals from the parent population, and the elite individuals were selected for crossover and mutation to obtain the offspring population.
1. Crossover method • Initial population The initial population, as the iterative starting point of the heuristic algorithm, is an important factor affecting population evolution and efficiency of the algorithm. When optimizing a single path, historical trajectory data were used to form the initial population to obtain the Pareto frontier solution. When optimizing multiple paths, the Pareto frontier solutions for a single path were randomly selected and combined into multiple path optimization individuals, forming the initial population. One should note that, in actual operation, the altitude and speed of an aircraft, being significantly sensitive to factors such as maneuvering actions, changes in aerodynamic configuration, and meteorological conditions, are very prone to fluctuations. Therefore, to determine the optimal altitude-velocity profile under ideal conditions, the altitude and velocity data in the historical trajectory data were smoothed [46]. • Genetic manipulation The tournament algorithm and elite retention strategy were used to select individuals from the parent population, and the elite individuals were selected for crossover and mutation to obtain the offspring population.

1.
Crossover method For fixed path scenario, the one-point crossover was performed in units of decision variables (altitude and speed) corresponding to an aircraft, i.e., a random crossover point is selected and the first part from each parent is copied to the corresponding child, and the second parts are copied crosswise. For the multi-path scenario, two crossover methods with the same probability were adopted. Case 1 is the crossover of paths that lead to the crossover of altitude and speed simultaneously. For Case 2, the chromosomes with the same path are selected for crossover of the second layer, which is the same as the crossover operator in fixed-path scenario. It should be noted that, in structural airspace, there are often a limited number of bifurcation nodes in the route network. Therefore, for Case 1, in order to improve the efficiency of path search, the selection of crossover nodes is not completely random, but randomly selected from the bifurcation nodes in the route network.

2.
Mutation method In this paper, the aircraft altitude and speed profiles are optimized by mutation independently to generate more solutions along the Pareto-optimal front. Once the mutation node i m is selected, the altitude and speed are mutated with with equal probability: When the aircraft is climbing or descending, according to altitude and speed of the predecessor node and successor node of mutation node, based on Formulas (15) (27).
When the aircraft is cruising, altitude and speed were independently disturbed using the disturbance values r h and r v with equal probability, respectively, as shown in Formula (28). The mutation method of the fixed and optional path scenarios is the same.
One should note that, because the predecessor node has a greater impact on the successor node, after crossover and mutation, offspring chromosomes' feasibility must be verified. The three gene positions before, in the middle of, and after crossover and mutation are used to determine whether a chromosome meets the constraint conditions. If it does, the offspring chromosome is obtained; otherwise, correction is made using the correction operator.

• Correction operator
The altitude profile is taken as an example to demonstrate the correction operator, as shown in Figure 4. The crossover node is i c , the mutation node is i m , the blue solid line is the original trajectory, and the red solid line is the trajectory of the aircraft after crossover or mutation. Obviously, dramatic changes in altitude between adjacent nodes is not possible. The following three scenarios illustrate the correction operator. Case 1 is the crossover correction operator, while Case 2 and Case 3 are the mutation correction operators.

1.
Case 1 For the crossover operation, as shown in Figure 4a, the speed relationship between node i c and node i k should be determined firstly. On the segment and v i c > v i k means the aircraft conducted accelerated climb, climb at a constant speed and decelerated climb, respectively, until reached the same altitude as the node i k . The altitude and speed of any node from i c to i k is corrected by using Formulas (1)-(10) and the rate of climb recommended by the BADA 3.11 manual. The situation of the descent is similar, so it is not repeated here.

2.
Case 2 After i m , if there exists h i l > h i m , the correction method is illustrated in Figure 4b.
During the segment i m → i k , the altitude and speed correction method is the same as the crossover correction method (see i c → i k in Figure 4a), and is not repeated. On segment i k → i l , the trajectory is further corrected by a cruise segment at constant altitude and speed. Thus far, the correction from i m to i l is completed. In a similar way, the correction works if there exists h i l < h i m after i m .

3.
Case 3 If the aircraft cruises at a constant altitude after i m , the correction method is as shown in Figure 4c. The altitude and speed of each subsequent node after i m are calculated and corrected using Formulas (1)-(10) and the rate of climb recommended by the BADA 3.11 manual. When h i k = h i m , the climb ends, and the aircraft cruises at constant altitude and speed from i k to the top of decent (TOD). The altitude and speed of each node during the descent can be estimated using Formulas (1)-(10) and the rate of descent recommended by the BADA 3.11 manual. Thus far, the correction is completed. In a similar way, the correction works for modification of a sharp descent due to mutation.
• Non-dominant sorting and crowding distance sorting The non-dominated sorting method is used to quickly stratify all individuals in the population, forming multiple Pareto frontiers with different levels. Crowding degree is used to measure the distance from one solution to two adjacent solutions on the frontier of the same non-dominated solution. • Termination condition The algorithm is terminated by setting the maximum number of iterations and testing the convergence results in each generation of the population. When the number of Pareto optimal frontier solutions in the population exceeds 80% of the population size, the algorithm is automatically terminated.

Route Structure and Operation Restrictions
In this study, one of China's busiest routes, ZSSS-ZBAA, was selected for optimization. Its flight routes and the sectors passed by flights are shown in Figure 5. Among them, waypoint data primarily include the longitude and latitude of the waypoint and the connection between waypoints, forming a basic network structure. The runway-to-runway 4D trajectory optimization has received focus without considering details of the aircraft's airport surface operation, so departure from and arrival at airports are each regarded as one network node. To study fuel consumption under temporary changes in the capacity-flow relationship, the concept of a "hotspot" was introduced to characterize the busy area when the capacity and flow lose balance during a certain period [47]. By analyzing the historical air traffic data, Sector 5 and Sector 6 are identified as "hotspots". City-pair flight data and operating restrictions are shown in Tables A1 and A2.

Typical Aircraft Types and Performance Parameters
In 2019, the Airbus A333 was the most frequently used aircraft on the ZSSS-ZBAA city pair (accounting for 37% of total flights). Therefore, Airbus A333 was used as an example aircraft, and its basic performance parameters are shown in Table 1.

Performance Parameter Value
Reference mass (kg) 174,000 Wing area (m 2 ) 361.6 Parasitic drag coefficient C D0 0.019805 Induced drag coefficient C D2 0.031875 First thrust-specific fuel consumption coefficient C f 1 0.61503 Second thrust-specific fuel consumption coefficient C f 2 919.03 Fuel flow correction factor during cruise C f cr 0.93655

Optimization Parameters
The optimization algorithm was written in Python 3.6.3. The parameters in the NSGA-II algorithm were set as follows: the population size was set to 1000, the maximum number of iterations was set to 100, the mutation probability was set to 0.1, the crossover probability was set to 0.95, r h = 600 m, and r v = 1 m/s. The relevant parameters of an Airbus A333 are shown in Table 1. The SAT for ZSSS-ZBAA is 105 min, where the maximum allowable advance and delay are 5 and 10 min, respectively.
The flight level is one of the important factors that affect aircraft fuel efficiency. Generally, air density is lower at higher altitude, so the aircraft receives less drag and is more fuel-efficient. In actual flight, the choice of altitude must consider the restrictions of waypoints and ATC preferences during flight. The height profile of an Airbus A333 on a typical day during a week is shown in Figure 6. Clearly, before reaching TOD, the aircraft descends to a relatively low altitude and maintains a level flight for a certain distance. The cruising altitude level of an Airbus A333 for the week and level-flight altitude level before TOD are shown in Table 2. Flights that cruised at 10,400 m and higher accounted for 80.9%. Considering that the ultimate objective of optimization is to increase the aircraft's fuel efficiency, the initial cruising level was set to 10,400 m and above. At the flight level before TOD, 98.3% of aircraft used 8400 and 9200 m, which were selected as the initial flight altitudes before TOD.  13.60% --

Verification of Fuel Consumption Model
First, to verify the accuracy of the fuel calculation model, QAR data were used to estimate fuel consumption. QAR data of A333 from ZSSS to ZBAA with mass of 172,365 kg were chosen, as shown in Table A3. According to the aircraft's speed, altitude, and ROCD from takeoff to landing, the fuel consumption model established in this paper is used to estimate fuel flow, and the estimated fuel flow is compared with the actual fuel flow, as shown in Figure 7. The fuel consumption change pattern calculated with the BADA model is basically the same as the trend of the recorded QAR data, and the value is similar. In general, the Mean Relative Error (MRE) from the BADA method is 12.3%, and the R 2 is 91.2%. The fuel flow prediction accuracy during climb and level-flight phases is relatively satisfactory with the average relative errors 9.6% and 7.1%, respectively. The overall estimation accuracy is within the acceptable range [48]. Therefore, the accuracy of the BADA fuel consumption model used in this study meets the requirements for trajectory optimization. The optimal flight trajectory was planned for each of the four horizontal paths shown in Table A1. The Pareto frontier solution in each path is shown in Figure 8. The red line segment is the reference baseline of the SAT. The results show that flight time has a significant impact on fuel consumption, and, within the preset time range, the overall relationship is shown to be a linear one with similar slope, i.e., approximately 118.8 kg is required for every minute of advance. Meanwhile, without considering capacity constraints, the flight distance has an absolute impact on the Pareto frontier position, i.e., as the flight distance increases, the fuel consumption under the same flight time gradually increases, and the flight time under the same fuel consumption increases. Clearly, under the same operating constraints and optimization environment, any one of the frontier solutions with a longer horizontal path cannot be better than other frontier solutions with a shorter one in terms of fuel and flight time. In addition, as the flight distance increases, the Pareto frontier begins to show an inflection point, because the official standard flight time does not leave too much room for flights to arrive early.
To investigate the trade-off space between fuel consumption and flight time, the altitude and speed of the aircraft during flight were analyzed. Figure 9 shows a comparison of the optimal and historical trajectory profiles corresponding to the solution with minimum fuel consumption (Trajectory 1) and minimum travel time (Trajectory 2) in Route 1. The red solid line is Trajectory 1, the blue solid line is Trajectory 2, and the black dotted lines are the historical trajectories of Airbus A333 on 1 July 2019. It shows that there were many low-altitude level-flight phases in the original trajectory, and the altitude and speed fluctuated greatly during cruise, resulting in higher fuel consumption. The difference between the results for Trajectories 1 and 2 obviously comes from the change in speed. Figure 9 shows that the aircraft reaches top of climb (TOC) at 284 km and TOD at 969 km when fuel consumption is minimized. The cruising altitude is maintained at 11,600 m. To determine the minimum fuel consumption, the flight speed was reduced to 201 m/s. To meet the control preference, the cruise level drops to 8400 m at 740 km away from departure airport. When the flight time is minimized, the aircraft reaches TOC at 276 km and TOD at 958 km. The cruising altitude is maintained at 11,600 m, and the cruising speed is maintained at 230 m/s. To meet the control preference, the cruising altitude drops to 8400 m at 733 km from the departure airport. High altitude and low speed cruising can significantly reduce fuel consumption while meeting timeliness constraints.  To evaluate the optimal upper bound of fuel efficiency in ideal conditions, the total fuel consumption before and after flight trajectory optimization on a typical day were compared. The historical fuel consumption during the flight was calculated using the fuel consumption model presented in this paper. First, the minimum fuel consumption in Route 1 and standard flight time were taken as upper bounds for comparison with the historical fuel consumption in each flight to obtain the ideal upper bound on fuel consumption reduction. To further increase the practical significance of the comparison, the same or similar flight time (the MRE of flight time of selected solutions in Routes 1-4 are 0.07%, 0.07%, −0.02%, and −0.01%) corresponding to the optimal fuel consumption solution in the Pareto frontier along different paths was selected for comparison, as shown in Figure 10. Among them, a positive value represents reductions compared to the historical fuel consumption/time, and vice versa for a negative value. Specifically, as shown in Table 3, the aircraft operates under the ideal altitude-speed profile corresponding to the minimum fuel consumption in Route 1 with standard flight time. The average ideal upper bounds on fuel consumption reduction are 13.38% and 3.36%, respectively; when similar flight times were selected for Routes 1-4, the average fuel consumption reduction rates were 8.05%, 6.55%, 5.17%, and 4.06%, respectively. Under the flight time limit, choosing a flight path with a shorter horizontal distance has greater potential for fuel consumption reduction.
Furthermore, the flight time was varied from 100 to 115 min in 5 min increments. The group of flights on a typical day is shown in Table 3. Multi-factor analysis of variance was used to evaluate the significance of the influence of route and flight time on fuel consumption. First, the normality of the data was tested, and the results are shown in Table A4. The results show that the data obey a normal distribution. Second, the data were tested for homogeneity of variance. Table A5 shows a Levene test table for error variance. The results show that the p-value of the Levene homogeneity of variance test statistic is 0.886, which is >0.05, so the data can be taken to meet the homogeneity of variance. Therefore, the data satisfy the premise of performing analysis of variance. Table 4 shows the main effects test results of the analysis of variance. The results show that the p-values of the main effects of time period and path are both <0.05, indicating that flight time and path have a significant impact on fuel consumption reduction. The p-value of interaction factors between time period and path is 0.999, which is greater than 0.05, implying that the effect of time on fuel saving potential does not change significantly with the choice of path. The result seems to conform to the reality because in this case the lengths of different paths do not differ too much, such that, no matter which path is chosen, the impact of flight time on fuel saving is roughly identical.     In summary, in a traffic scenario with no capacity constraints on a fixed path, the trajectory optimization method proposed in this paper can weigh the trade-off between flight time and fuel consumption, plan an ideal altitude-speed profile, and obtain upper bounds on fuel consumption reduction under different constraints. The model provides a basis for determining the ideal city-pair flight energy savings and emission reductions determined by 4D trajectory optimization.

Scenario 2: Trajectory Optimization Analysis Based on Fixed and Available Flight Path Sets under Capacity Constraints
In addition to the highly dynamic atmospheric environment, potential reductions in aircraft fuel consumption in a real flight are affected by uncertainty in capacity and flow. The sector capacity-flow relationship from historical data was analyzed and four A333s were selected during the busy period of 10:00-14:00 on a typical day for experimental analysis. The congestion level in sectors from ZSSS to ZBAA at 10:00-14:00 is shown in Figure 11. Sector capacity was calculated using 15 min as statistical granularity and 5 min as rolling capacity in the sliding window, and sector congestion is the ratio of flow to capacity. The updated available time slots for a hotspot are shown in Table A6. The "time of enter hotspot" in the table was determined using the lowest fuel consumption trajectory of Route 1 (no capacity constraint) in Scenario 1.
analysis. The congestion level in sectors from ZSSS to ZBAA at 10:00-14:00 is shown in Figure 11. Sector capacity was calculated using 15 min as statistical granularity and 5 min as rolling capacity in the sliding window, and sector congestion is the ratio of flow to capacity. The updated available time slots for a hotspot are shown in Table A6. The "time of enter hotspot" in the table was determined using the lowest fuel consumption trajectory of Route 1 (no capacity constraint) in Scenario 1.  In the above scenario, fixed and optional path trajectory optimization for capacity constraints were carried out, as shown in Figure 12. In general, the limited capacity reduces the space for trade-off optimization of flight trajectories, and a more flexible horizontal operation space helps increase the potential reduction in fuel consumption. Compared with the fixed route (Route 1) and the historical flight trajectory, the minimum fuel consumption of the optional route trajectory optimization decreased by 1.45% and 8.25% on average. Under the same and similar flight time, fuel consumption on the optional route was decreased by 1.30% and 6.89%, respectively. The planned entry time of the four flights at the hotspot is different from the time slot for allowed entry, and the trajectory optimization potential is significantly different. Among them, the optimization results for Flights 1 and 4 are consistent with those for the fixed route, indicating that, when the planned arrival time is slightly later than the controlled arrival time, the aircraft can still obtain the lowest fuel consumption by adjusting the speed and altitude while maintaining the optimal fixed route (Route 1). In particular, when the planned arrival time is earlier than the controlled arrival time within a certain range, such as Flights 2 and 3, the Pareto frontier is very different from that of the fixed path. To further explore the causal factors driving these results, the flight trajectory was further optimized using Routes 2 and 4 as fixed paths under capacity constraints. The frontier solutions are shown in Figure 12. The Pareto optimal frontiers for different paths for Aircraft 2 are shown in Figure 12b. Under capacity constraints, the fuel efficiency from high to low is Route 2, Route 4, Route 1, and Route 3. The reason for this ranking is that the available time slots for aircraft flying over Sector 5 and Sector 6 are consistent and limited. Compared with Routes 1 and 3, Routes 2 and 4 have longer flight distances before entering the hotspot, and the aircraft can use smaller speed adjustments to remain in the time slot when entering the sector. In addition, Route 2's flight distance is shorter than that of Route 4's, so Route 4's optimal Pareto frontier and Route 2's (with capacity constraint) Pareto frontier nearly coincide. In addition, compared to the optimal Pareto frontier for Route 1 (without capacity constraint), the Pareto frontier trade-off space for all paths under capacity constraints is compressed, and flight time tends to be distributed above 110 min. This shows that the use of air speed control to delay the time required for the aircraft to enter the sector has limited space for trajectory optimization.
The Pareto optimal frontiers for different paths of Aircraft 3 are shown in Figure 12c. Under capacity constraints, the fuel efficiency from high to low is found on Route 3, Route 4, Route 2, and Route 1. The reason for this ranking is that Routes 3 and 4 have the same time slots to enter Sector 6, and Sector 6's available time slot is more relaxed than Sector 5's. In addition, Route 3 has a shorter flight distance than Route 4. Therefore, Route 4's optimal Pareto frontier and Route 3's (with capacity constraint) Pareto frontier nearly coincide. Routes 1 and 2 enter Sector 5 within the same available time slot, but Route 2 has a longer flight distance before flying over the hotspot, such that the aircraft can use smaller speed adjustments to remain in the time slot when entering the sector, resulting in better fuel efficiency compared to Route 1. In particular, the optimal frontier for Route 1 (with capacity limitation) is very tight and is quite close to the maximum flight time of 115 min, which suggests the lowest trade-off flexibility in trajectory selection. It shows that, when the aircraft does not have optional route space, and the sector is available within the allotted time slot, the method of decelerating to delay the aircraft's entry into Sector 5 and meet the timeliness constraints has almost reached its limit. This better reflects the importance of optional route space for reducing travel time.
To explore the potential of optional routes for fuel efficiency improvement in the case of airspace congestion, the minimized fuel consumption point determined with optional path trajectory optimization was used as the upper bound and compared with the minimum fuel consumption determined for Routes 1-4 held fixed under capacity constraints. To further increase the practical significance of the comparison, the optimal fuel consumption solutions corresponding to the flight times that were similar or equal to the historical data were chosen (the MRE of flight time of selected solutions in Routes 1-4 are −0.05%, −0.03%, −0.12%, and 0.02%, respectively) from the Pareto frontiers for the optional and fixed paths for comparison, as shown in Figure 13. Specifically, compared with that of Routes 1-4, the average fuel consumption reduction for the optional route at its optimal fuel consumption point is 1.45%, 1.35%, 1.17%, and 1.41%, respectively. At the same or similar flight time points, the average fuel consumption reduction is 1.31%, 1.64%, 2.39%, and 2.70%, respectively. This verifies the effectiveness of the path selection model, as shown in Table 5.   Furthermore, the fuel consumed by the four A333s during the busy period of 10:00-14:00 on a typical day for optional and fixed routes is shown in Table A7. Similarly, a multi-way analysis of variance was carried out to evaluate whether the effect of the optional route on fuel consumption is significant. First, normality of the data was tested, and the results are shown in Table A8. The results show that the data follow a normal distribution. Second, the data were tested for homogeneity of variance. Table A9 shows a Levene test table for error variance. The results show that the significance of the test statistics for Levene homogeneity of variance is 0.477 > 0.05, so variance can be considered homogeneous. Therefore, the analysis of variance can be used to analyze these data. Table 6 shows analysis of variance results, indicating that optional routes have a significant and positive effect (p = 0.019) on fuel consumption.

Fuel Consumption Reduction Potential
The upper bound on potential fuel consumption reduction for city-pair flights was determined using trajectory optimization without changing the airspace structure. Compared with historical fuel consumption data, when an aircraft operates with an ideal altitudespeed profile corresponding to minimal fuel consumption on Route 1 and standard flight time in ideal operating conditions, the average ideal upper bounds for fuel consumption optimization are 13.38% and 3.36%, respectively, in the congestion-free scenario ( Figure 14); regarding the same (similar) flight time, average fuel consumption reductions in Routes 1-4 are 8.05%, 6.55%, 5.17%, and 4.06%, respectively, as shown in Figure 14. When the aircraft operates under the capacity limit, the fuel consumption reductions for the optional route and Routes 1-4 are 8.99%, 6.92%, 7.01%, 7.18%, and 6.95%, respectively. One can see that the 4D trajectory multi-objective optimization model proposed in this paper can be used to examine trade-offs between fuel consumption and flight efficiency. By evaluating the fuel consumption reduction potential in different scenarios, different airspace rigidities, and different flight time requirements, the upper bounds on fuel consumption reduction can be determined for different comparison benchmarks. This provides a basis for determining city-pair flight energy savings and emission reduction targets based on 4D trajectory optimization, as shown in Figure 14.

Sensitivity Study
In order to further explore the impact of congestion location and uncertainty level on fuel efficiency, a sensitivity analysis was conducted using time-optimal trajectory along Route 1 in congestion-free scenario as the baseline (go through Sector 2, Sector 3, Sector 4, Sector 5, Sector 7, Sector 8 and Sector 9 in turn). Series of scenarios were designed by two independent variables: congested sector and severity of congestion (i.e., difference between the available entry slot and the planned arrival time).
Trajectory optimizations in above scenarios are shown in Figure 15. The time offset between the available time slot and the planned arrival time of the congested sector ranged from 2 to 12 min, incremented by 2 min. In general, the time offset and congested sector (interpreted as the distance to the departure airport) have a significant impact on the trajectory optimization effect. With the increase of the severity of congestion, it becomes less likely to arrive at the destination airport in advance. Meanwhile, the trade-off solution space decreases gradually. In addition, as the congested sector gets farther away from the departure airport, the aircraft has stronger ability to tolerate higher variations by adjusting speed and altitude profile. In particular, when the time offset is 12 min, solutions can be achieved if and only if the congested sector is Sector 9, which is closest to the landing airport.
In order to further quantify the impact of airspace congestion on the optimization effect, and to explore the internal relationship between aircraft fuel efficiency and temporalspatial characteristics of airspace congestion, the Hypervolume (HV) index is used to evaluate the pros and cons of the Pareto front solution [49], as shown in Formula (29). R i is a rectangle constructed by each element in the Pareto front and the reference point as the diagonal vertices, as shown in Figure 16. It should be noted that, when calculating HV, the data should be normalized.
The normalized (1,1) is selected as the reference point, and the results of HV are shown in Table 7. For the same sector, larger deviation would result in worse optimization effect. Interestingly, when the degree of congestion is small to medium (i.e., the time offset is 2-6 min), the fuel efficiency increases at first and then drops, as the congested sector gets close to the landing airport. In addition, for different time offsets, the "best" sectors varies accordingly, as shown in Table 7. It implies that, for low-to-medium variation, more balanced space on both sided of the congested sector would be more beneficial for trajectory optimization. When the congestion is getting severe (i.e., the time offset is 8-12 min), the "best" sector with best performance is always Sector 9 (see Table 7). It implies that allocating heavy delay to the airspace sector that is farthest away from the departure airport would achieve optimal solutions within predefined tolerable delay (i.e., [−5 min, 10 min]) in the model. The results are fairly fit with the nature of uncertainty, that is uncertainty is magnified over time. Therefore, more robust and refined ATFM strategies would significantly improve the performance of Pareto optimal front of trajectory optimization.

Conclusions
Studying the fuel-saving potentials and their influencing factors in the TBO context is a very first step of setting environmental performance targets in the development of the next generation of the ATM system, which is still an open topic receiving insufficient concern. This paper focuses on developing runway-to-runway 4D trajectory optimization models with heterogeneous constraints of travel time boundaries, route options, and sector capacity limitations from both airspace users' and air traffic managers' perspectives. The significance of fuel consumption reduction was examined for city-pair (ZSSS-ZBAA) trajectory optimization in ideal and congested operation scenarios with either fixed or optional routes. The congestion impact on the fuel-saving potentials was supplementarily included. In the following, our main conclusions, some implications, and future work are stated.

• Main conclusions
In the congestion-free scenario with fixed route, travel time has a significant impact on fuel consumption. Overall, linear trade-off characteristics between fuel burn and flight time within the predefined range of tolerable delay was observed (i.e., 1 min in advance requires an extra 118.8 kg of fuel). Within a predefined travel time boundary, cruising at higher altitude with lower speed would burn less fuel. In the congestionfree scenario, the flight distance simply has a negative impact on the optimal Pareto front (i.e., an extra 23.6 kg of fuel required for each kilometer), the slopes of which are similar for different route choices. The results imply that, considering a tolerable delay, choosing a flight path with shorter horizontal distance would achieve greater potential for fuel consumption optimization. In the capacity-constrained scenario, which is the common situation in daily ATFM (Air Traffic Flow Management) practice, the tradeoff space of optimal trajectory solutions is compressed due to fewer available entry time slots in congested sectors, which would significantly reduce the fuel efficiency. Nevertheless, more energy-saving potentials would benefit from flexible route options. Apparently, the trajectory optimization effect is sensitive to the congestion location and severity. Generally, within travel time constraints, the maximum acceptable sector entry time delay after takeoff (i.e., ATFM uncertainty) increases as the congested sector becomes further away from the departure airport. The best performance would be achieved when the entry time uncertainty (less than 6 min) emerges around the middle of the flight path. Therefore, in the TBO context, aircraft performance margin, and RTA constraints of individual flight trajectory should be taken into consideration in balancing capacity and air traffic flow among multiple sectors, in order to achieve minimum deviation to the planned trajectory while maintaining high operational efficiency, economy, and predictability. • Limitations and future work This paper provides a preliminary but adaptable framework for predicting and evaluating city-pair fuel efficiency improvement potentials in the TBO context. A more comprehensive study that expands the research scope from a single city pair to the entire airport network including domestic and international flights would achieve more insights for setting green development goals in future ATM. More impact factors shall be considered in addition to route choice, flight time, and coarse-grained sector congestion, in investigating fuel-saving potentials at the entire network level by incorporating airport capacity constraints, incidental capacity shortfalls, Miles (Minutes)-in-Trail regulations, ATC separation management, prevailing climate, or emerging severe weather conditions. A big data-driven fuel assumption modeling would be a promising approach to deal with the above uncertainties and even emergency and in-stable air traffic operation scenarios. In addition, a higher degree of routing flexibility would unleash fuel-saving potentials in the TBO context. Free Route Airspace (FRA) [50], compared with structural airspace with preset but limited routing options, is a more promising means of increasing trajectory-based operation performance, which would be worthwhile for advancing airspace planning and management. For the China airspace system, a more practical way to enhance the operational flexibility is to design modularized airspace that retains the original ATS (Air Traffic Service) route and adaptively configures entry and exit points, as well as free maneuvering zones based on the availability of airspace resources underpinned by Flexible Use of Airspace (FUA). Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: