A Decision Support Framework for Aircraft Arrival Scheduling and Trajectory Optimization in Terminal Maneuvering Areas

Abstract

This study introduces a decision support framework that integrates aircraft trajectory optimization and arrival scheduling to facilitate efficient management of descent operations for arriving aircraft within terminal maneuvering areas. The framework comprises three modules designed to tackle specific challenges in the descent process. The first module formulates and solves a trajectory optimization problem, generating a range of candidate descent trajectories for each arriving aircraft. The options for descent operations include step-down descent operation, Continuous Descent Operation (CDO), and CDO with a lateral path stretching strategy. The second module addresses the assignment of conflict-free trajectories to aircraft, determining precise arrival times at each waypoint. This is achieved by solving an aircraft arrival scheduling problem. To overcome computational complexities, a novel variable neighborhood search algorithm is proposed as the solution approach. This algorithm utilizes three neighborhood structures within an extended relaxing and solving framework, and incorporates a tabu search algorithm to enhance the efficiency of the search process in the solution space. The third module focuses on comparing the total cost incurred from flight delays and fuel consumption across the three descent operations, enabling the selection of the most suitable operation for the descent process. The decision support framework is evaluated using real air traffic data from Guangzhou Baiyun International Airport. Experimental results demonstrate that the framework effectively supports air traffic controllers by scheduling more cost-efficient descent operations for arrival aircraft.

1. Introduction

Air transport demand has been growing rapidly in recent decades. Despite the significant impact of the COVID-19 outbreak on air transportation, high levels of air traffic are still projected in the near future [1,2]. It is estimated that the demand will surge by 50% within the next 17 years, accompanied by a corresponding increase in delay frequencies estimated to be 15 times greater than the current levels [3]. Short-term infrastructure development for capacity expansion is often limited by financial, geopolitical, and ecological factors. Therefore, optimizing the efficient use of existing infrastructure and managing air traffic within airport Terminal Maneuvering Areas (TMAs) play a significant role in mitigating airport congestion. Additionally, fuel efficiency is a top priority for airlines as fuel consumption costs constitute a major expense for them [4]. Furthermore, the environmental and health impacts of air traffic have prompted increased scrutiny from both the government and the public, primarily due to concerns related to noise and emissions. Consequently, aviation authorities continuously set ambitious targets to address various aspects of Air Traffic Management (ATM), including increasing airport capacity, minimizing aircraft delays, and reducing the environmental impact of air traffic while adhering to ATM safety regulations [5,6].

However, it is worth noting that some of these targets may conflict with each other when put into practice. One such example is the objective of minimizing aircraft delays, which may require increasing aircraft speed and cruising altitude to avoid flight conflicts in congested air segments. Unfortunately, this can result in higher fuel consumption and increased emissions. Furthermore, Air Traffic Controllers (ATCs) currently have limited support from automated systems, which hinders their ability to make well-informed decisions [7]. At present, ATCs primarily rely on their experience and tactical instructions, such as heading vectors, speed adjustments, and path stretching, to ensure safe separation between aircraft and maximize the throughput within the TMA. This situation is primarily attributed to the lack of granularity and efficiency in aircraft trajectory management and scheduling techniques.

Trajectory Based Operations (TBO) is an advanced concept in ATM that has been extensively explored by initiatives such as SESAR [8] and NextGen [9]. The primary objective of TBO is to transition the current air traffic control system into a collaborative one that enables precise coordination of aircraft operations and achieves higher accuracy in all dimensions. By efficiently negotiating trajectories within TBO, both airborne and ground equipment can be synchronized, allowing for the establishment of a descent plan with prescribed landing times before an aircraft enters the TMA. This approach ensures equitable and advantageous trajectories for aircraft. The management of aircraft landings within the TMA using TBO involves several stages. Firstly, the Flight Management System (FMS) on board the aircraft calculates a set of trajectories based on its descent operation, including the travel time for each air segment on the assigned landing route. This information is then transmitted to the ground system. Subsequently, the ground system sequences the aircraft and examines for any potential conflicts between trajectories. Finally, ATCs plan the descent operation, taking into account the appropriate decision criteria to ensure a satisfactory outcome.

Within the framework of TBO, the Continuous Descent Operation (CDO) has gained recognition as a vital operational technique. This procedure involves guiding arrival aircraft through an optimized profile descent using low engine thrust settings and a low drag configuration. Unlike the conventional Step-down Descent Operation (SDO), CDO aircraft descend continuously from the Top of Descent (TOD), aiming to avoid level segments. This approach leads to reduced fuel consumption and pollutant emissions [5]. Nevertheless, managing the separation between consecutive aircraft during the CDO process can pose challenges for ATCs as they need to closely monitor aircraft trajectories and increase separation buffers to avoid flight conflicts. As a result, the capacity of both airspace and runways in the TMA can be negatively affected, particularly during peak hours [10]. In our previous work [11], we developed a small proof of concept to demonstrate the successful implementation of a tactical lateral path stretching strategy in conjunction with the CDO for arrival aircraft in high-density traffic. This approach, termed ps-CDO (CDO with a lateral path stretching strategy), provides alternative lateral paths to each arrival aircraft, offering greater flexibility in avoiding conflicts within the TMA. By providing multiple options for the lateral path of each aircraft during the CDO phase, considering different arrival times, an arrival aircraft scheduling problem can be formulated and systematically resolved. This enables effective sequencing of the aircraft and prevents potential conflicts.

In this paper, we present a decision support framework that builds upon our previous work on collaborative aircraft trajectory optimization and arrival scheduling. The primary objective of this framework is to assist ATCs in effectively managing descent operations for arrival traffic within the TMA, while concurrently optimizing the utilization of the TMA’s capacity. Our proposed framework consists of three synthesis modules: descent trajectory generation, aircraft arrival scheduling, and optimal trajectory selection. These modules specifically focus on three distinct descent operations: the SDO, the CDO, and the ps-CDO. To elaborate, the first module formulates and solves an Aircraft Trajectory Optimization Problem (ATOP), which generates descent trajectories for each arriving aircraft considering all three descent operations. These descent trajectories define the time windows for each air segment across all available descent routes. Based on the generated descent trajectories, the second module assigns each aircraft a conflict-free trajectory along with the corresponding descent route, resulting in the determination of precise arrival times at each waypoint. This process involves solving a trajectory-based Aircraft Arrival Scheduling Problem (AASP) formulated as a Mixed Integer Programming (MIP). To find solutions, we employ a Variable Neighborhood Search (VNS) algorithm. Lastly, the third module determines the optimal descent operation selection among the three options by considering both total delay and total fuel consumption costs. This selection process aims to strike a balance between minimizing overall delays and reducing fuel consumption.

1.1. Literature Review

The aviation industry has developed several decision support tools to assist ATCs in managing arriving flows in the TMA, such as the Center-TRACON Automation Systems [12] and the Arrival Manager systems [13]. These systems primarily offer landing aircraft with arrival time decisions and runway assignments while aiding the ATCs in guiding aircraft to fly along specific waypoints along the established air route. These aforementioned decision support tools, however, have difficulty providing more complex trajectory management solutions, which may result in an increase in the ATCs’ workload. Developments in communication, navigation, and surveillance technologies have contributed to increasing trajectory accuracy, leading to expectations for a more advanced decision support tool based on trajectories. The tool attempts to sequence and merge arrival flow by performing time and energy management operations on descending trajectories.

The efficient management of aircraft descent trajectories is crucial for reducing fuel consumption, gas emissions, and noise during aircraft operations. Although effective trajectory management during aircraft descent is successful in airports with orderly and steady traffic, it is challenging in heavy congestion traffic situations. Recent studies have focused mainly on trajectory optimization technologies to address the challenges of managing aircraft descent trajectories in congested traffic situations. Itoh et al. [14] and Callantine et al. [15] investigate the application of flight-deck interval management for aircraft descent trajectories in a TMA. Turgut et al. [16] introduce a fixed flight path angle descent (FFPAD) technology to help management aircraft descent trajectories in a congested traffic environment. The feasibility of managing the arrival aircraft trajectories in a busy TMA with FFPAD is fully investigated in [17]. Sun et al. [18] use the four-dimensional waypoints technology to achieve a conflict-free trajectory in high-density airspace, and evaluate the feasibility of operations with a leader-follower example. Nonetheless, the aforementioned methods are only practical for particular operational equipment and aircraft and limited in terms of broader practical applications.

On the other hand, an effective trajectory management solution based on the TBO concept involves sequencing and merging arrival traffic by assigning required times of arrival (RTAs) at one or several fixes along a predetermined route. Pawelek et al. [19] propose a two-step approach to implement the full CDOs with RTA in a TMA. First, they formulated and solved a trajectory optimization problem to generate a set of candidate CDO trajectories with times of arrival at a specific metering fix for arrival aircraft. Next, they built an aircraft scheduling model to optimally assign RTAs to the arriving traffic, ensuring safe time separation. As a result, aircrew can fly optimal descent trajectories using the FMS while meeting the required times of arrival. However, the feasibility of this approach was only evaluated under a low air traffic environment. Sáez et al. [10] develop a four-dimensional trajectory negotiation and synchronization process to complete the CDO procedure with a path stretching strategy known as the tromboning paradigm in dense TMAs with the same approach. They take into account assigning optimal routes of the trombone and RTAs to the arriving traffic simultaneously. Nevertheless, the aircraft scheduling model used by [10,19] was established with only a certain metering fix, which did not include all descent waypoints in the TMA. In addition, Samà et al. [20] propose an integrated approach for aircraft trajectory optimization and scheduling for both arrival and departure aircraft, and their method involves having arrival aircraft descend using the SDO and determining optimal trajectories and sequencing for all aircraft simultaneously to demonstrate its benefits. Yet they fail to account for the scenario where an aircraft could have multiple descent routes, nor do they apply the CDO procedure.

The classical aircraft arrival scheduling problem aims to allocate a runway to each aircraft for landing and schedule a corresponding landing time for each. The classical aircraft arrival scheduling problem is commonly referred to as the Runway Scheduling Problem (RSP). The RSP is typically formulated and solved as a job-shop scheduling model [21]. Bennell et al. [22] conducted an extensive literature review of the RSP, while a more recent comprehensive survey of relevant studies was made by [23]. There are a number of algorithms that have been developed to solve the RSP since it has been proven to be a typical NP-hard problem [24]. The particular nature of the RSP prompts researchers to compute high-quality solutions with low computational time. Most of the research thus favors heuristic and meta-heuristic algorithms. The commonly used meta-heuristic algorithms are the genetic algorithm [25,26,27], the simulated annealing [28,29], and the ant colony optimization [30,31]. Additionally, the Variable Neighborhood Descent (VND), which belongs to the family of VNS, has received much attention as well. Salehipour et al. (2009) [32] develop VND to address the RSP, wherein four neighborhood structures were designed to improve the incumbent solutions. Salehipour et al. (2013) [33] consider the same problem as [32] and use VND in a hybrid algorithm by applying the simulated annealing framework. Furthermore, VND was also used in [34,35] to address the RSP. Recently, Salehipour and Ahmadian [36] propose a novel hybrid optimization approach to tackle the RSP. This approach can be decomposed into two steps. Firstly, an initial sequence is generated using the target landing time first rule. Then, the sequence is improved by a local search algorithm. With the sequence fixed, the runway scheduling problem is formed as a linear program, and the landing time is obtained by solving the linear program with the solver CPLEX. Ahmadian and Salehipour [37] extend this research by introducing a novel framework called “Relax-and-Solve” (R&S) as the sequence local search. They called it “Relax-and-Solve” because it improves the aircraft sequence at the runway by iteratively destructing (relaxing) a sub-sequence of the current sequence and reconstructing (solving) a feasible sequence using optimization techniques. In this paper, we incorporate the R&S framework into our VNS algorithm to address the AASP.

In a more realistic context, Bianco et al. [38] take into account the restrictions of various TMA resources for aircraft arrival scheduling, such as holding circles and air segments in addition to runways. Specifically, this model is formulated as a no-wait job shop scheduling problem, where the runways and air segments are interpreted as the machines, and the aircraft to be scheduled are the jobs. The separation time between each pair of aircraft corresponds to the sequence-dependent set-up time, and the earliest and latest landing time are interpreted as the release date and due date, respectively. Solving the resulting model is expected to yield the start time of each aircraft at every resource along its descent route, such that the aircraft have conflict-free descent trajectories satisfying all restrictions. Based on this effort, D’Ariano et al. [39] formulate an aircraft scheduling problem with fixed routes, which is represented by an alternative graph. Refs. [40,41,42,43] extend the work by including more assumptions and objectives to improve the level of accuracy in modeling the aircraft scheduling problem with fixed routes. Our previous work [11] covers similar ground. Samà et al. (2014) [44] propose a new formulation that considers alternative routes for each aircraft. In order to solve the aircraft scheduling problem with alternative routes, Samà et al. (2017) [45] present a VNS algorithm. They employ five distinct neighborhood structures by capturing the characteristic of the objective function to improve the search for executing the routes for each aircraft. The AGLIBRARY solver, a state-of-the-art optimization solver designed to deal with complex routing and scheduling problems, is then used to solve the problem with the searched routes. Our approach differs from [45] in both modeling and solution approaches. Specifically, our AASP formulation aims to minimize a linear combination of the total delay and the total difference between the scheduled landing time and the minimum fuel-burning landing time. This contrasts with their study, which primarily focuses on minimizing the maximum delay, with the VNS implementation adhering to this objective as well. In addition, our approach encompasses three descent operations, modeled with different equations of motion for the trajectory optimization, whereas [45] considers only the step-down descent procedure. The algorithm proposed by [45] is not applicable in our framework due to the different models used, and we thus develop a new VNS algorithm that employs an extended R&S framework to establish neighborhood structures during the search for improved routes for each aircraft.

1.2. Our Contributions

This paper introduces a novel decision support framework for managing aircraft descent operations in the TMA and makes several contributions to the existing literature:

(1) We develop a comprehensive decision support framework that can assist ATCs in determining the optimal descent operations and arrival schedules for approaching aircraft in TMAs. To the best of our knowledge, this is the first approach that integrates the aircraft trajectory optimization and arrival scheduling and incorporates various options for selecting descent operations, including SDO, CDO, and ps-CDO. By extensively comparing three different descent operations and their impacts on flight delays and fuel consumption in various traffic flows, our proposed tool can provide its users with an important reference for selecting the most suitable operations for arrival aircraft.

(2) We propose a novel VNS algorithm that employs three neighborhood structures to enhance the efficiency and quality of solutions for solving the AASP. This algorithm leverages an extended R&S framework to establish the neighborhood structures, and incorporates a Tabu Search (TS) algorithm to efficiently search the solution space.

(3) We implement extensive experiments using real traffic data from Guangzhou Baiyun International Airport (GBIA), considering both arrival and departure traffic flows. The computational results demonstrate the effectiveness and efficiency of our proposed decision support framework and our proposed solution approach.

1.3. Organization of This Paper

The remainder of this paper is organized as follows: Section 2 provides an overview of the decision support framework’s architecture and describes the mathematical models used in the framework. Section 3 discusses the solution methodologies employed for addressing the problem. The numerical experiments conducted and the corresponding results are presented in Section 4. Finally, Section 5 offers a conclusion summarizing the findings of this study and suggests potential directions for future research.

2. Decision Support Framework

In this section, we present our proposed decision support framework, starting with the introduction of the concept of operations proposed in this work (Section 2.1). We then provide a detailed description of the mathematical models proposed for the three modules within the decision support framework (Section 2.2, Section 2.3 and Section 2.4).

2.1. Concept of Operations

In Figure 1, we provide a simplified scenario that serves as an illustration of the problem addressed in this paper. The scenario considers a time-limited decision horizon where multiple aircraft seek to enter an extended-TMA (E-TMA) from different directions through various entry points. The ATCs are responsible for scheduling the descent process for each aircraft.

For each approach, there are predefined descent routes available, and the scheduling process aims to assign a suitable descent route to each aircraft along with a scheduled arrival time at each waypoint. This ensures that every aircraft follows a conflict-free descent trajectory through the different air segments and eventually lands on a runway. The traffic flows may converge at various metering waypoints, where the aircraft trajectories merge until they ultimately reach the runways.

Well before the TOD point, while still in the cruise phase, the ATC requires the aircraft to compute its travel time window for each air segment. The travel time window is determined by calculating the difference between the earliest and latest possible times for flying through the air segments, corresponding to the earliest and latest descent trajectories. The FMS is responsible for calculating these trajectories, taking into account various factors such as aircraft performance, flight envelope, and weather conditions. The width of the travel time window can be influenced by the choice of descent operation employed by the aircraft. In this paper, we consider three descent operations: the SDO, the CDO, and the ps-CDO.

Figure 2a illustrates the vertical profiles of the earliest and latest trajectories for the three descent operations. It also presents the travel time window for each descent operation on an air segment. It can be observed that the SDO offers a higher degree of operational flexibility, resulting in a wider time window. However, this flexibility comes at the cost of increased fuel consumption. On the other hand, the ps-CDO extends the width of the CDO time window through a lateral path stretching strategy, which may also lead to increased fuel consumption. The CDO is the most restrictive descent operation, with the narrowest time window. However, it is also the most fuel-efficient descent operation.

After entering the E-TMA during the cruising phase, the aircraft engages in a trajectory and scheduled time negotiation process with the ATCs. In order to facilitate this synchronization process, we propose a decision support framework, which is illustrated in Figure 3. The framework consists of three modules, outlined as follows:

• Module I: Descent trajectories generation.
  This module is responsible for generating descent trajectories for arriving aircraft, considering predefined descent operations such as the SDO, the CDO, and the ps-CDO. These trajectories encompass the earliest and latest descent trajectories, as well as the minimum fuel-burning descent trajectory. The travel time window is determined based on the earliest and latest descent trajectories, while the minimum fuel-burning landing time is associated with the minimum fuel-burning descent trajectory. To achieve this, an ATOP is formulated and solved, taking into account factors such as aircraft performance, flight envelope, and weather conditions (refer to Section 2.2). It is important to note that multiple descent routes are considered for each aircraft, and the descent trajectories differ for each route. In practice, these trajectories would be generated using advanced functionality in the FMS for each aircraft.
• Module II: Aircraft arrival scheduling.
  In the second module, based on the computed descent trajectories from the previous module, a descent route with a conflict-free trajectory is assigned to each aircraft operating under one of the three descent operations. This assignment process ensures a safe separation between aircraft during their descent procedure and determines the exact arrival time at each waypoint. The objective of this model is to minimize a linear combination of the total delay and the total difference between the minimum fuel-burning landing time and the scheduled landing time. The preferences for the total delay and the total difference are typically provided by the decision makers, in this case, the ATCs in our studied system. The AASP is formulated as a MIP (see Section 2.3) and a VNS algorithm is developed (see Section 3.2) to solve the problem.
• Module III: Optimal trajectory selection.
  This module is initiated by the decision makers (i.e., the ATCs) who choose the optimal descent operation for every arrival aircraft within the decision time horizon. We define a cost function as a linear combination of the total delay cost and total fuel consumption cost to determine the priority among the three descent operations (see Section 2.4). In this model, the total delay and the required landing time for all arriving aircraft can be calculated by solving the AASP in Module II. Additionally, the ATOP in Module I can be utilized to compute the minimum fuel-burning descent trajectory that satisfies the required landing time, from which the resulting fuel consumption can be determined. Consequently, the approaching aircraft will follow the optimal descent operation with the minimum fuel-burning descent trajectories that meet the required landing time. Notably, in the AASP, if decision makers increase the weight of the total delay indicator in the objective function to achieve a smaller total delay, it will inevitably lead to an increase in the total difference between the minimum fuel-burning landing time and the scheduled landing time. This increased difference signifies a greater deviation from the minimum fuel-burning descent trajectory, thereby leading to higher fuel consumption for each aircraft. Conversely, decreasing the weight of the total delay indicator will have the opposite effect. This approach allows for the selection of the most suitable descent operation based on the cost function and the specific requirements of each aircraft.

Without loss of generosity, we make the following assumptions to expedite the discussion of our proposed framework: (1) All waypoints are predefined in the descent procedure for each arriving aircraft. (2) There are no operational failures during descent, such as missed approaches, emergency landings, or engine failures, and operating time is not taken into account. (3) The aircraft maintains its original cruise speed until it reaches the TOD. (4) Only one aircraft is allowed to land or take-off on each runway at a time.

2.2. Descent Trajectories Generation

This section presents the methodology employed to generate aircraft descent trajectories as they approach their destination. The descent trajectories consist of two segments: the cruise segment from an entry waypoint to the TOD, and the descent segment from the TOD to the runway. The descent segment is further divided into multiple phases, which correspond to different changes in aircraft dynamics. To optimize the aircraft’s vertical profile, including altitude and speed, while maintaining a fixed lateral route and distance, we formulate the problem as a multi-phase optimal control problem. This formulation takes into consideration altitude and speed constraints, as well as aircraft performance bounds. In this study, we utilize a point mass model [46] to describe the aircraft’s equations of motion. The trajectory of each arriving aircraft along a pre-designed descent route is planned such that the aircraft follows the point mass equations in each air segment $p\in \mathcal{P}=\{1,2,\cdots ,|\mathcal{P}\left|\right\}$ and satisfies boundary conditions between any two consecutive air segments. The ATOP mathematical formulation is presented below.

$$\begin{array}{}\mathrm{(1a)}& \left(\mathbf{ATOP}\right):\hfill & min\hfill & J\hfill \mathrm{(1b)}& & \mathrm{s}.\mathrm{t}.& {\dot{V}}_T^{\left(p\right)}=\frac{T^{\left(p\right)}-D^{\left(p\right)}(V_T,h,\beta )}{m}-g{\gamma }^{\left(p\right)}\forall p\in \mathcal{P}\hfill \mathrm{(1c)}& & & {\dot{s}}^{\left(p\right)}=V_T^{\left(p\right)}cos{\gamma }^{\left(p\right)}+W_s^{\left(p\right)}\left(h\right)\forall p\in \mathcal{P}\hfill \mathrm{(1d)}& & & {\dot{h}}^{\left(p\right)}=V_T^{\left(p\right)}sin{\gamma }^{\left(p\right)}\forall p\in \mathcal{P}\hfill \mathrm{(1e)}& & & {\dot{h}}_{min}^{\left(p\right)}\le \frac{dh}{dt}\le {\dot{h}}_{max}^{\left(p\right)}\forall p\in \mathcal{P}\hfill \mathrm{(1f)}& & & T_{min}^{\left(p\right)}\le T^{\left(p\right)}\le T_{max}^{\left(p\right)}\forall p\in \mathcal{P}\hfill \mathrm{(1g)}& & & 0\le {\beta }^{\left(p\right)}\le {\beta }_{max}^{\left(p\right)}\forall p\in \mathcal{P}\hfill \mathrm{(1h)}& & & {\gamma }_{min}^{\left(p\right)}\le {\gamma }^{\left(p\right)}\le 0\forall p\in \mathcal{P}\hfill \mathrm{(1i)}& & & V_{CAS,min}^{\left(p\right)}\le V_{CAS}^{\left(p\right)}\le V_{CAS,max}^{\left(p\right)}\forall p\in \mathcal{P}\hfill \mathrm{(1j)}& & & M_{min}^{\left(p\right)}\le M^{\left(p\right)}\le M_{max}^{\left(p\right)}\forall p\in \mathcal{P}\hfill \mathrm{(1k)}& & & {[V_T,h]}^{⊺}\left(t_0^{\left(1\right)}\right)={[{V_T}_0,h_0]}^{⊺}\hfill \mathrm{(1l)}& & & {[V_T,s,h]}^{⊺}\left(t_f^{\left(\right|\mathcal{P}\left|\right)}\right)={[{V_T}_f,s_f,h_f]}^{⊺}\hfill \mathrm{(1m)}& & & {[V_T,s,h]}^{⊺}\left(t_f^{(p-1)}\right)={[V_T,s,h]}^{⊺}\left(t_0^{\left(p\right)}\right)\forall p\in \mathcal{P}\setminus \left\{1\right\}\hfill \end{array}$$

For each air segment p, the variables have the following definitions: $V_T^{\left(p\right)}$ is the true airspeed; $s^{\left(p\right)}$ is the along-track distance from the runway; $h^{\left(p\right)}$ is the altitude; ${\gamma }^{\left(p\right)}$ is the aerodynamic flight path angle; $T^{\left(p\right)}$ and $D^{\left(p\right)}$ represent the thrust and the aerodynamic drag, respectively; ${\beta }^{\left(p\right)}$ is the speed-brakes deflection; $W_s^{\left(p\right)}\left(h\right)$ is the longitudinal wind; $V_{CAS}^{\left(p\right)}$ and $M^{\left(p\right)}$ are the calibrated airspeed and Mach number, respectively; $t_0^{\left(p\right)}$ and $t_f^{\left(p\right)}$ denote the initial and final time, respectively. The parameters g and m are the aircraft gravity acceleration and mass, respectively.

Constraints (1b)–(1d) represent the point-mass dynamics of an aircraft in the vertical plane, where ${[V_T,s,h]}^{⊺}$ is the state vector and ${[T,\beta ,\gamma ]}^{⊺}$ is the control vector. Constraints (1e)–(1j) specify the value range requirements of $\frac{dh}{dt}$, T, $\beta$, $\gamma$, $V_{CAS}$ and M to ensure flight envelope protection and passenger comfort, where ${\dot{h}}_{min}$ and ${\dot{h}}_{max}$ are the minimum and maximum descent rate; $T_{\min }$ and $T_{max}$ are the idle and maximum thrust; ${\beta }_{max}$ is the maximum speed-brakes deflection; ${\gamma }_{min}$ is the minimum descent gradient; $V_{CAS,min}$ and $V_{CAS,max}$ are the minimum and maximum calibrated airspeed; $M_{min}$ and $M_{max}$ are the minimum and maximum Mach number. Constraints (1k) and (1l) represent the boundary condition at the TOD and the runway. Constraint (1m) is the link constraint.

Module I calculates the earliest and latest descent trajectories, as well as the minimum fuel-burning descent trajectory, for each arrival aircraft. Each trajectory serving its specific objectives:

Earliest descent trajectory, which is to minimize the flight time,

$$\begin{array}{}\mathrm{(2a)}& \left({\mathbf{ATOP}}_{\mathbf{E}}\right):& min& J_{time}=\underset{\mathrm{cruise}\mathrm{segment}\mathrm{flight}\mathrm{time}}{\underbrace{\frac{s\left(t_0^{\left(1\right)}\right)-s_{max}}{V_{cr}}}}+\underset{\mathrm{descent}\mathrm{segment}\mathrm{flight}\mathrm{time}}{\underbrace{t_f^{\left(\right|\mathcal{P}\left|\right)}}}\hfill \mathrm{(2b)}& & \mathrm{s}.\mathrm{t}.& \left(\mathrm{1b}\right)–\left(\mathrm{1m}\right)\hfill & \end{array}$$

where $s(t_0^{\left(1\right)})$ is an along-track distance at TOD; $s_{max}$ is the along-track distance of entry waypoint; $V_{cr}$ is the cruise speed.

Latest descent trajectory, which is to maximize the flight time:

$$\begin{array}{}\mathrm{(3a)}& \left({\mathbf{ATOP}}_{\mathbf{L}}\right):& max& J_{time}\mathrm{(3b)}& & \mathrm{s}.\mathrm{t}.& \left(\mathrm{1b}\right)–\left(\mathrm{1m}\right)\end{array}$$

Minimum fuel-burning descent trajectory, which is to minimize the fuel consumption:

$$\begin{array}{}\mathrm{(4a)}& \hfill \left({\mathbf{ATOP}}_{\mathbf{F}}\right):& min& J_{fuel}=\underset{\mathrm{cruise}\mathrm{segment}\mathrm{fuel}\mathrm{consumption}}{\underbrace{\frac{f_{cr}(s\left(t_0^{\left(1\right)}\right)-s_{max})}{V_{cr}}}}+\underset{\mathrm{descent}\mathrm{segment}\mathrm{fuel}\mathrm{consumption}}{\underbrace{\sum _{p\in \mathcal{P}}{\int }_{t_0^{\left(p\right)}}^{t_f^{\left(p\right)}}{f}_{idle}(V_T,h)dt}}\hfill \mathrm{(4b)}& & \mathrm{s}.\mathrm{t}.& \left(\mathrm{1b}\right)–\left(\mathrm{1m}\right)\hfill \end{array}$$

where $f_{cr}$ and $f_{idle}$ are the cruise fuel flow rate and idle fuel flow rate in the descent, respectively.

The ATOP may encounter different operational constraints for each aircraft in each air segment depending on the chosen descent operation. The parameters related to aircraft performance, which are utilized in both the point mass model and the fuel burn function, are derived from the Base of Aircraft Data (BADA) [47]. It is important to note that the additional thrust and speed-brakes are only permitted for the SDO. For the CDO and the ps-CDO, the point-mass model is reduced to a so-called $\gamma -command$ model, which assumes continuous vertical equilibrium through the descent process. Specifically, in order to obtain an environmentally friendly CDO trajectory, the idle thrust is imposed and the use of speed brakes is not allowed throughout the descent. As a result, the flight path angle ($\gamma$) is the only control variable in this dynamical system to manage the energy of the aircraft and achieve different travel time through the descent air segments.

We solve the $\left({\mathbf{ATOP}}_{\mathbf{E}}\right)$ formulation to obtain the earliest travel time for each aircraft in each air segment p, which is $t_f^{\left(p\right)}-t_0^{\left(p\right)}$. Similarly, solving the $\left({\mathbf{ATOP}}_{\mathbf{L}}\right)$ formulation provides us with the aircraft’s latest travel time in each air segment p. The minimum fuel-burning landing time (which is $t_f^{\left(\right|\mathcal{P}\left|\right)}$) for each aircraft is determined by solving the $\left({\mathbf{ATOP}}_{\mathbf{F}}\right)$ formulation. These times are then utilized in Module II.

2.3. Aircraft Arrival Scheduling

In this section, we formulate and solve the AASP to assign a conflict-free trajectory with a descent route to each aircraft. This scheduling process aims to determine the precise arrival time of each aircraft at every waypoint, ensuring compliance with safety separation requirements. The AASP takes into account various inputs, including the travel time window for each air segment and the minimum fuel-burning landing time (both obtained from the ATOP). Additionally, estimated earliest and latest arrival times for aircraft entering the TMA and the estimated target landing time are considered. It is worth mentioning that the estimated target landing time remains consistent regardless of the chosen descent operation.

The AASP is formulated using a directed graph $\mathcal{G}=(\mathcal{V},\mathcal{A})$ with a set of waypoints $\mathcal{V}$ and a set of arcs $\mathcal{A}$. For each aircraft $i\in \mathcal{F}$, a descent air route can be chosen in a set of alternative descent air routes ${\mathcal{R}}_i$, where $\mathcal{F}$ is the set of all arrival aircraft. Specifically, a descent air route $r_i$ for aircraft i contains a set of waypoints $(o_i,\cdots ,k,l,\cdots ,d_i)\in {\mathcal{V}}_i^{r_i}$. The origin $o_i$ and the destination $d_i$ represent the entry waypoint and the runway, respectively. Each arc $(k,l)\in {\mathcal{A}}_i^{r_i}$ indicates the air segment between two adjacent waypoints in the directed graph. Therefore, the set of waypoints that can be used by aircraft i is ${\mathcal{V}}_i={\bigcup }_{r_i\in {\mathcal{R}}_i}{\mathcal{V}}_i^{r_i}$. Similarly, the set of arcs that can be used by aircraft i is ${\mathcal{A}}_i={\bigcup }_{r_i\in {\mathcal{R}}_i}{\mathcal{A}}_i^{r_i}$. In this connection, ${\mathcal{V}}_i\subset \mathcal{V}$, ${\mathcal{A}}_i\subset \mathcal{A}$ in the digraph $\mathcal{G}$.

We introduce two sets of binary decision variables to respectively model the descent air route selection and the sequential relationship between aircraft at shared waypoints. The variable $y_i^{r_i}$ represents the selection of descent air route $r_i$ for aircraft i. The variable $z_{iju}$ indicates the sequential relationship between aircraft i and j as they traverse the same waypoint $u\in {\mathcal{V}}_i^{r_i}\cap {\mathcal{V}}_j^{r_j}$. This is equal to 1 if aircraft i arrives at waypoint u before aircraft j and 0 otherwise. We use the continuous variable $t_{ik}^{r_i}$ to represent the arrival time at waypoint k for aircraft i following descent air route $r_i$.

The objective of the AASP is to minimize a linear combination of the total delay and the total difference between the scheduled landing time and the minimum fuel-burning landing time. The sets, parameters, and decision variables used in the formulation of the AASP are summarized in

Table 1. Sets, parameters, and variables in the AASP.

Sets with indices	Explanation
$\mathcal{F}$	A set of aircraft (index $i,j$).
${\mathcal{R}}_i$	A set of alternative descent air routes for aircraft i (index $r_i$).
$\mathcal{V}$	A vertex set of waypoints in the TMA (index $o_i,\cdots ,k,l,\cdots ,d_i$), where ${\mathcal{V}}_i\in \mathcal{V}$.
$\mathcal{A}$	An air segment set of descent air route in the TMA, where ${\mathcal{A}}_i\in \mathcal{A}$.
$\mathcal{G}$	A directed graph $\mathcal{G}=(\mathcal{V},\mathcal{A})$.
Parameters	Explanation
$i,j$	Aircraft ID.
$r_i$	Descent air route for aircraft i.
$k,l,u,v$	Transit waypoint.
$o_i$	The entry waypoint for aircraft i.
$d_i$	The runway for aircraft i.
${\Gamma }_{i{o}_i}^{\mathrm{early}},{\Gamma }_{i{o}_i}^{\mathrm{late}}$	Estimated earliest (latest) arrival time for aircraft i enter the TMA.
${\Gamma }_{i{d}_i}^{\mathrm{early}},{\Gamma }_{i{d}_i}^{\mathrm{target}},{\Gamma }_{i{d}_i}^{\mathrm{late}}$	Estimated earliest (target, latest) landing time at the runway $d_i$ for aircraft i.
${\Gamma }_{i{d}_i}$	The reference landing time at the runway $d_i$ for aircraft i, where ${\Gamma }_{i{d}_i}=max\{{\Gamma }_{i{d}_i}^{\mathrm{early}},{\Gamma }_{i{d}_i}^{\mathrm{target}}\}$.
${\Gamma }_{i{d}_i}^{\mathrm{fuel}}$	The landing time at the runway $d_i$ with the minimum fuel consumption for aircraft i.
${\underline{\Theta }}_{ik\_il}^{r_i},{\overline{\Theta }}_{ik\_il}^{r_i}$	The minimum (maximum) travel time in air segment $(k,l)$ for aircraft i.
$S_{iju}$	The minimum time-based separation for a preceding aircraft i and another trailing aircraft j in a same waypoint u.
B	Large artificial variable.
Decision variables	Explanation
$y_i^{r_i}$	1, if aircraft i uses descent air route $r_i\in R_i$; 0, otherwise.
$z_{iju}$	1, if aircraft i flies through the same waypoint u before aircraft j; 0, otherwise.
$t_{ik}^{r_i}$	The arrival time of waypoint k by descent air route $r_i$ for aircraft i, where $t_{ik}^{r_i}\ge 0$.
${\alpha }_i$	The delay of aircraft i to its reference landing time, where ${\delta }_i=max\{0,t_{i{d}_i}^{r_i}-{\Gamma }_{i{d}_i}\}$.
${\mu }_i$	The absolute value of the difference $|t_{i{d}_i}^{r_i}-{\Gamma }_{i{d}_i}^{\mathrm{fuel}}|$.

. The AASP is formulated as follows:

$$\begin{array}{}\mathrm{(5a)}& \left(\mathbf{AASP}\right):& min& \lambda \sum _{i\in \mathcal{F}}{\alpha }_i+(1-\lambda )\sum _{i\in \mathcal{F}}{\mu }_i\hfill \mathrm{(5b)}& & \hfill \mathrm{s}.\mathrm{t}.\hfill & \sum _{r_i\in {\mathcal{R}}_i}{y}_i^{r_i}=1\forall i\in \mathcal{F}\hfill \mathrm{(5c)}& & & t_{ik}^{r_i}\le B{y}_i^{r_i}\forall i\in \mathcal{F},\forall k\in {\mathcal{V}}_i^{r_i},\forall r_i\in {\mathcal{R}}_i\hfill \mathrm{(5d)}& & & \begin{array}{c}{\Gamma }_{i{o}_i}^{\mathrm{early}}{y}_i^{r_i}\le t_{i{o}_i}^{r_i}\le {\Gamma }_{i{o}_i}^{\mathrm{late}}{y}_i^{r_i}\hspace{1em}\forall i\in \mathcal{F},\forall r_i\in {\mathcal{R}}_i\hfill \end{array}\mathrm{(5e)}& & & \begin{array}{c}{\underline{\Theta }}_{ik\_il}^{r_i}{y}_i^{r_i}\le t_{il}^{r_i}-t_{ik}^{r_i}\le {\overline{\Theta }}_{ik\_il}^{r_i}{y}_i^{r_i}\hfill \\ \forall i\in \mathcal{F},\forall (k,l)\in {\mathcal{A}}_i^{r_i},\forall r_i\in {\mathcal{R}}_i\hfill \end{array}\mathrm{(5f)}& & & \begin{array}{c}{t}_{ju}^{r_j}\ge t_{iu}^{r_i}+S_{iju}-B(3-z_{iju}-y_i^{r_i}-y_j^{r_j})\hfill \\ \forall i,j\in \mathcal{F},i<j,\forall u\in {\mathcal{V}}_i^{r_i}\cap {\mathcal{V}}_j^{r_j},\forall r_i\in {\mathcal{R}}_i,\forall r_j\in {\mathcal{R}}_j\hfill \end{array}\mathrm{(5g)}& & & \begin{array}{c}{t}_{iu}^{r_i}\ge t_{ju}^{r_j}+S_{jiu}-B(2+z_{iju}-y_i^{r_i}-y_j^{r_j})\hfill \\ \forall i,j\in \mathcal{F},i<j,\forall u\in {\mathcal{V}}_i^{r_i}\cap {\mathcal{V}}_j^{r_j},\forall r_i\in {\mathcal{R}}_i,\forall r_j\in {\mathcal{R}}_j\hfill \end{array}\mathrm{(5h)}& & & \begin{array}{c}{z}_{iju}-z_{ijv}\ge \sum _{\begin{array}{c}{r}_i\in {\mathcal{R}}_i\\ (u,v)\in {\mathcal{A}}_i^{r_i}\end{array}}{y}_i^{r_i}+\sum _{\begin{array}{c}{r}_j\in {\mathcal{R}}_j\\ (u,v)\in {\mathcal{A}}_j^{r_j}\end{array}}{y}_j^{r_j}-2\hfill \\ \forall i,j\in \mathcal{F},i<j,\forall (u,v)\in {\mathcal{A}}_i^{r_i}\cap {\mathcal{A}}_j^{r_j}\hfill \end{array}\mathrm{(5i)}& & & \begin{array}{c}{z}_{ijv}-z_{iju}\ge \sum _{\begin{array}{c}{r}_i\in {\mathcal{R}}_i\\ (u,v)\in {\mathcal{A}}_i^{r_i}\end{array}}{y}_i^{r_i}+\sum _{\begin{array}{c}{r}_j\in {\mathcal{R}}_j\\ (u,v)\in {\mathcal{A}}_j^{r_j}\end{array}}{y}_j^{r_j}-2\hfill \\ \forall i,j\in \mathcal{F},i<j,\forall (u,v)\in {\mathcal{A}}_i^{r_i}\cap {\mathcal{A}}_j^{r_j}\hfill \end{array}\mathrm{(5j)}& & & {\alpha }_i\ge t_{i{d}_i}^{r_i}-{\Gamma }_{i{d}_i}\forall i\in \mathcal{F},\forall r_i\in {\mathcal{R}}_i\hfill \mathrm{(5k)}& & & {\alpha }_i\ge 0\forall i\in \mathcal{F}\hfill \mathrm{(5l)}& & & {\mu }_i\ge t_{i{d}_i}^{r_i}-{\Gamma }_{i{d}_i}^{\mathrm{fuel}}\forall i\in \mathcal{F},\forall r_i\in {\mathcal{R}}_i\hfill \mathrm{(5m)}& & & {\mu }_i\ge {\Gamma }_{i{d}_i}^{\mathrm{fuel}}-t_{i{d}_i}^{r_i}\forall i\in \mathcal{F},\forall r_i\in {\mathcal{R}}_i\hfill \mathrm{(5n)}& & & t_{ik}^{r_i}\ge 0\forall i\in \mathcal{F},\forall k\in {\mathcal{V}}_i^{r_i},\forall r_i\in {\mathcal{R}}_i\hfill \mathrm{(5o)}& & & z_{iju}\in \{0,1\}\forall i,j\in \mathcal{F},i<j,\forall u\in {\mathcal{V}}_i^{r_i}\cap {\mathcal{V}}_j^{r_j}\hfill \mathrm{(5p)}& & & y_i^{r_i}\in \{0,1\}\forall i\in \mathcal{F},\forall r_i\in {\mathcal{R}}_i\hfill \end{array}$$

where $\lambda \in [0,1]$ is a weight parameter that allows for the adjustment of the relative importance of total delay and total difference of landing times. The variable ${\alpha }_i$ denotes the delay of aircraft i with respect to its reference landing time ${\Gamma }_{i{d}_i}$. The variable ${\mu }_i$ denotes the absolute value of the difference $|t_{i{d}_i}^{r_i}-{\Gamma }_{i{d}_i}^{\mathrm{fuel}}|$.

Constraint (5b) ensures that exactly one descent air route can be chosen for each aircraft. Constraint (5c) enforces that $t_{ik}^{r_i}$ is a non-zero value when the descent air route $r_i$ is selected, where B is a sufficiently large number such as the sum of the maximum travel time of all aircraft flying through all arcs. Constraint (5d) restricts that each aircraft i to enter the TMA between the estimated earliest arrival time ${\Gamma }_{i{o}_i}^{\mathrm{early}}$ and the latest arrival time ${\Gamma }_{i{o}_i}^{\mathrm{late}}$. Constraint (5e) allows the travel time in each air segment $(k,l)$ to take any value within the range $\left[{\underline{\Theta }}_{ik\_il}^{r_i},{\overline{\Theta }}_{ik\_il}^{r_i}\right]$. These lower and upper bounds are obtained by solving the aforementioned $\left({\mathbf{ATOP}}_{\mathbf{E}}\right)$ and $\left({\mathbf{ATOP}}_{\mathbf{L}}\right)$ formulations, respectively. Constraints (5f) and (5g) represent separation constraints, ensuring that if aircraft j follows aircraft i at the same waypoint $u\in {\mathcal{V}}_i^{r_i}\cap {\mathcal{V}}_j^{r_j}$, the arrival time at waypoint u of aircraft j must be greater than or equal to the arrival time of aircraft i plus the minimum time-based separation $S_{iju}$. Constraints (5h) and (5i) ensure that there is no overtake between aircraft i and j traveling through the common air segment $(u,v)\in {\mathcal{A}}_i\cap {\mathcal{A}}_j$, where ${\sum }_{r_i\in {\mathcal{R}}_i,(u,v)\in {\mathcal{A}}_i^{r_i}}{y}_i^{r_i}=1$ (${\sum }_{r_j\in {\mathcal{R}}_j,(u,v)\in {\mathcal{A}}_j^{r_j}}{y}_j^{r_j}=1$) if aircraft i (j) uses arc $(u,v)$, 0, otherwise. Constraints (5j) and (5k) define that the delay of aircraft i equals to $max\{0,t_{i{d}_i}^{r_i}-{\Gamma }_{i{d}_i}\}$, where ${\Gamma }_{i{d}_i}=max\{{\Gamma }_{i{d}_i}^{\mathrm{early}},{\Gamma }_{i{d}_i}^{\mathrm{target}}\}$. ${\Gamma }_{i{d}_i}^{\mathrm{early}}$ and ${\Gamma }_{i{d}_i}^{\mathrm{target}}$ are the earliest landing time and target landing time at the runway $d_i$ of aircraft i, respectively, where ${\Gamma }_{i{d}_i}^{\mathrm{early}}={\Gamma }_{i{o}_i}^{\mathrm{early}}+{\sum }_{(k,l)\in {\mathcal{A}}_i^{r_i}}{\underline{\Theta }}_{ik\_il}^{r_i}$. Constraints (5l) and (5m) define the absolute value of difference between $t_{i{d}_i}^{r_i}$ and ${\Gamma }_{i{d}_i}^{\mathrm{fuel}}$. The minimum fuel-burning landing time ${\Gamma }_{i{d}_i}^{\mathrm{fuel}}$ is obtained by solving the $\left({\mathbf{ATOP}}_{\mathbf{F}}\right)$ formulation. Constraints (5n)–(5p) specify the domains of the decision variables.

Note that the $\left(\mathbf{AASP}\right)$ can also be employed to allocate a feasible ascent route and the start time at each waypoint for every departing aircraft within a specified decision horizon in a TMA. In this context, each take-off aircraft follows an ascent trajectory, flying from the runway towards its assigned exit point. In the problem formulation, the origin $o_i$ refers to the runway, while the destination $d_i$ represents the exit waypoint.

2.4. Optimal Trajectory Selection

This module helps the ATCs select the optimal descent operation for each approaching aircraft. A cost function is established to prioritize the descent operations, which is formulated as the sum of the cost of total delay and the cost of total fuel consumption, as shown below:

$$\begin{array}{c}\Omega =c_{delay}{\Omega }_{delay}+c_{fuel}{\Omega }_{fuel}\end{array}$$

(6)

where $c_{delay}$ and $c_{fuel}$ are the unit delay cost (€/s) and unit fuel cost (€/kg), respectively. ${\Omega }_{delay}$ denotes the total delay experienced by all aircraft during a given decision time horizon, which is computed by solving the $\left(\mathbf{AASP}\right)$. ${\Omega }_{fuel}$ denotes the total fuel consumption of all aircraft within the decision time horizon. To calculate the fuel consumption of each aircraft, an ATOP is formulated with a fixed required time of arrival at the runway. The objective of the ATOP is to minimize fuel consumption, and the formulation is as follows:

$$\begin{array}{}\mathrm{(7a)}& \hfill \left({\mathbf{ATOP}}_{{\mathbf{F}}^{\prime }}\right):& min& J_{fuel}\hfill \mathrm{(7b)}& & \mathrm{s}.\mathrm{t}.& \left(\mathrm{1b}\right)–\left(\mathrm{1m}\right)\hfill \mathrm{(7c)}& & & t_f^{\left(\right|\mathcal{P}\left|\right)}=t_{RTA}\hfill \end{array}$$

where $t_{RTA}$ is the required time of arrival at the runway, which is obtained by solving the $\left(\mathbf{AASP}\right)$.

3. Solution Methods

The central component of our decision support framework consists of the proposed ATOP and AASP models. The ATOP is formulated as a Non-Linear Programming (NLP) optimization problem with interior point constraints, which can be solved numerically using two primary methods: direct and indirect [20,48]. In this study, a pseudospectral method [49] is employed as a direct method for solving the ATOP. On the other hand, the AASP is formulated as a MIP problem. While it can be solved using commercial solvers such as Gurobi, the high number and complexity of decision variables and constraints for a busy TMA can result in extremely long computational times. To address this issue, we have developed a VNS algorithm specifically tailored for solving the AASP. The VNS algorithm provides efficient solutions within a reasonable time frame, allowing for real-time decision-making. In practice, air traffic flow management in a TMA operates in a dynamic and evolving environment. To accommodate this dynamism, we adopt a rolling horizon approach to solve the dynamic cases of our decision support framework. This approach allows for continuous adaptation and optimization of trajectories and schedules based on updated information and changing conditions, ensuring effective and timely decision-making.

3.1. Pseudospectral Method

The pseudospectral method discretizes the time interval (time grid) and approximates the state and control variables in an optimal control problem. By converting the original infinite-dimensional problem into a finite-dimensional nonlinear optimization problem, it can be solved using standard NLP solvers like SNOPT [50]. In the pseudospectral method, the state and control inputs are approximated as piecewise polynomials, and collocation points are chosen using a quadrature rule over each time step. This approach offers high accuracy and efficiency, and for a large number of collocation points, it can yield an exact optimal solution [49].

For solving the ATOP in our study, we employ the GPOPS 5.0 package [51], which is a Matlab-based (R2021a) software specifically designed for solving multi-phase optimal control problems using the pseudospectral method. The SNOPT solver is integrated in the GPOPS as the underlying NLP solver, which is known for its robustness and efficiency in handling NLP problems. By utilizing GPOPS with the SNOPT solver, we can effectively solve the ATOP and obtain optimal trajectory solutions for aircraft descent operations in our decision support framework.

3.2. VNS Algorithm

The AASP consists of three main operations: routing, sequencing, and scheduling. Routing involves selecting the descent route for each aircraft, sequencing determines the arrival order of aircraft at each waypoint, and scheduling involves determining the arrival time for each aircraft at each waypoint. When a feasible sequence is given, which includes the selected descent route for each aircraft and the execution order of the aircraft at each waypoint, the optimal schedule for that sequence can be obtained in a polynomial time. This paper proposes a VNS algorithm that starts with an initial sequence and iteratively explores different neighborhood structures to search for improved sequences. Note that given a feasible sequence, obtaining an optimal schedule for each aircraft can be achieved using available algorithms or commercial solvers such as Gurobi.

The VNS algorithm follows a three-phase process: shaking, local search, and move. It starts with an initial sequence $\Pi$ and its corresponding objective value $z(\Pi )$. The algorithm alternates between these phases to iteratively improve the sequence. In the shaking phase, the algorithm generates a neighbor sequence ${\Pi }^{\prime }$ by applying one of the available neighborhood structures ${\mathcal{N}}_q$, for $q=1,2,\cdots ,q_{max}$, or a fixed structure ${\mathcal{N}}_1$. This neighbor sequence ${\Pi }^{\prime }$ is evaluated by calculating its objective value $z\left({\Pi }^{\prime }\right)$. Next, in the local search phase, the algorithm manipulates the neighbor sequence ${\Pi }^{\prime }$ using the available neighborhood structures to further improve it. If an improved sequence ${\Pi }^{″}$ with a lower objective value ($z\left({\Pi }^{″}\right)<z(\Pi )$) is obtained, it replaces the current best sequence, and the local search continues. If no improvement is achieved, the algorithm returns to the shaking phase. This process continues until the stopping criterion is met. The algorithm terminates based on two criteria: the maximum number of iterations ($ite{r}_{max}$) and the maximum number of iterations without improvement ($noimprovin{g}_{max}$). Once the algorithm terminates, the improved sequence $\Pi$ with its objective value $z(\Pi )$ is reported as the final solution. Algorithm 1 provides an overview of the VNS algorithm, while Algorithm 2 summarizes the steps involved in the local search phases. The generation of initial solutions and the definition of neighborhood structures are discussed in the following sections.

Algorithm 1 The variable neighborhood search (VNS) algorithm for AASP

Input:  Integers $ite{r}_{max}$, $noimprovin{g}_{max}$ and $q_{max}$. A set of neighbourhood structures ${\mathcal{N}}_q$ for $q=1,2,\cdots ,q_{max}$. An initial sequence $\Pi$ and its corresponding objective value $z(\Pi )$.   1: $iter\leftarrow 1$   2: $\xi \leftarrow 0$   3: while $iter\le ite{r}_{max}$ and $\xi <noimprovin{g}_{max}$ do   4:       $/\ast$ Shaking: $\ast /$   5:       ${\Pi }^{\prime },z\left({\Pi }^{\prime }\right)\leftarrow$${\mathcal{N}}_1(\Pi )$   6:       $/\ast$ Local search: $\ast /$   7:       ${\Pi }^{″},z\left({\Pi }^{″}\right)\leftarrow$$LS({\Pi }^{\prime },q)$    $/\ast$ See Algorithm 2 $\ast /$   8:       $/\ast$ Move or not: $\ast /$   9:       if $z\left({\Pi }^{″}\right)<z(\Pi )$ then 10:     $\Pi \leftarrow {\Pi }^{″}$ 11:     $z(\Pi )\leftarrow z\left({\Pi }^{″}\right)$ 12:     $\xi \leftarrow 0$ 13:        else 14:     $\xi \leftarrow \xi +1$ 15:     $iter\leftarrow iter+1$ 16:        end if 17: end while Output:  An improved sequence $\Pi$ and its objective value $z(\Pi )$.

Algorithm 2 The local search procedure $LS({\Pi }^{\prime },q)$ in VNS

1: $TL\leftarrow \varnothing$   2: $q\leftarrow 1$   3: while  $q\le q_{amx}$do   4:     $c{c}_q\leftarrow$ Randomly select one from all positions of ${\Pi }^{\prime }$.   5:     if $c{c}_q$ not in $TL$ then   6:      $temp,z\left(temp\right)\leftarrow {\mathcal{N}}_q({\Pi }^{\prime },c{c}_q,c{r}_q)$   7:      if $z\left(temp\right)<z\left({\Pi }^{\prime }\right)$ then   8:       ${\Pi }^{\prime }\leftarrow temp$   9:       $z\left({\Pi }^{\prime }\right)\leftarrow z\left(temp\right)$ 10:       $q\leftarrow 1$ 11:       TabuMemoryFilter($TL$) 12:      else 13:       $q\leftarrow q+1$ 14:       $TL\leftarrow$TabuMemoryUpdate($TL$, $T{L}_{max}$, $c{c}_q$, ${\Pi }^{\prime }$) 15:      end if 16:     end if 17: end while

3.2.1. Initial Solution Generation Method

The AASP exhibits a pattern where the order of aircraft at each waypoint corresponds to the order of aircraft landing on the runway based on their scheduled landing times. This alignment is a result of designing descent routes to converge aircraft arrivals towards the runway (see in Figure 1) and the prohibition of overtaking in any air segment (see constraints (5i) and (5j)). Therefore, we can represent the sequence $\Pi =\left({\pi }_1^{r_1},{\pi }_2^{r_2},\cdots ,{\pi }_{\left|\mathcal{F}\right|}^{r_{\left|\mathcal{F}\right|}}\right)$ for the AASP, where the order list $\{1,2,\cdots ,|\mathcal{F}\left|\right\}$ denotes the execution order of the aircraft at each waypoint, and the set $\{r_1,r_2,\cdots ,r_{\left|\mathcal{F}\right|}\}$ represents the descent route allocations.

To generate an initial solution for the AASP, we propose an algorithm outlined in Algorithm 3. The algorithm starts by sorting all aircraft at each waypoint in non-decreasing order based on their estimated earliest landing time ${\Gamma }_{i{d}_i}^{\mathrm{early}}$. This sorting ensures that the aircraft are initially arranged in a feasible order at each waypoint. We then formulate the problem $\left({\mathbf{AASP}}_{\mathbf{I}}\right)$, which fixes the aircraft order at each waypoint, and use the Gurobi 9.5.0 solver to solve it. The default parameters of Gurobi are used except for the time limit, which may be adjusted as needed. The initial sequence $\Pi$ is obtained by sorting all aircraft in non-decreasing order of their scheduled time at the runway. This sequence, along with its objective value $z(\Pi )$, is used to initialize the local search.

Algorithm 3 The algorithm for generating the initial solution

1: Sort the aircraft at each waypoint in non-decreasing order of $\Gamma \left(i\right)={\Gamma }_{i{d}_i}^{\mathrm{early}}$. 2: Solve the problem $\left({\mathbf{AASP}}_{\mathbf{I}}\right)$ utilizing a commercial solver like Gurobi. 3: Obtain the initial sequence $\Pi$ and its objective value $z(\Pi )$. $//$The problem$\left({\mathbf{AASP}}_{\mathbf{I}}\right)$ is formulated as follows: $$\begin{array}{}\mathrm{(8a)}& \begin{array}{ccc}& \left({\mathbf{AASP}}_{\mathbf{I}}\right):\hfill & min\hspace{1em}\left(\mathrm{5a}\right)\hfill \end{array}\mathrm{(8b)}& \begin{array}{cc}& \mathrm{s}.\mathrm{t}.\hspace{1em}{t}_{ju}^{r_j}\ge t_{iu}^{r_i}+S_{iju}-B(2-y_i^{r_i}-y_j^{r_j})\\ & \forall i,j\in \mathcal{F},\Gamma \left(i\right)<\Gamma \left(j\right),\forall u\in {\mathcal{V}}_i^{r_i}\cap {\mathcal{V}}_j^{r_j}\end{array}\mathrm{(8c)}& \hfill \left(\mathrm{5b}\right)–\left(5e\right),\left(\mathrm{5j}\right)–\left(\mathrm{5n}\right),\left(\mathrm{5p}\right)\hfill \end{array}$$

It is worth noting that we conducted preliminary experiments to determine the optimal rule for sorting aircraft at each waypoint. We tested three sorting rules based on estimated target landing time (${\Gamma }_{i{d}_i}^{\mathrm{target}}$), earliest estimated landing time (${\Gamma }_{i{d}_i}^{\mathrm{early}}$), and latest estimated landing time (${\Gamma }_{i{d}_i}^{\mathrm{late}}$). Our findings suggest that sorting the aircraft by the earliest estimated landing time yields the best solution.

3.2.2. Improvement Algorithm

Once a feasible initial sequence $\Pi$ has been obtained for the AASP, the VNS algorithm attempts to improve it by iteratively exploring its neighborhood structures for a designated number of iterations. Our VNS algorithm employs three unique neighborhood structures denoted as ${\mathcal{N}}_1$, ${\mathcal{N}}_2$, and ${\mathcal{N}}_3$.

The neighborhood structure ${\mathcal{N}}_1$ is derived from the R&S framework that allows for the exploration of improved solutions by relaxing routing and sequencing constraints for a selected sub-sequence of aircraft. This neighborhood structure differs significantly from traditional approaches. The process of the neighborhood structure ${\mathcal{N}}_1$ is illustrated in Figure 4. Given an initial sequence $\Pi$, the relaxing neighborhood selects a sub-sequence ${\Pi }^R$ from $\Pi$ based on two parameters: a randomly selected center $c{c}_1$ and a radius $c{r}_1$. Every sub-sequence ${\Pi }^R$ is relaxed through relaxing the routing and sequencing constraints of the aircraft in the sub-sequence, therefore letting the aircraft change their descent route and order in the sequence. The remaining aircraft in the sequence, denoted as ${\Pi }^{NR}=\Pi \setminus {\Pi }^R$, are not subject to changes in routing or sequencing. We then reduce the $\left(\mathbf{AASP}\right)$ with a given ${\Pi }^R$ and ${\Pi }^{NR}$ to problem $\left({\mathbf{AASP}}_{\mathbf{R}}\right)$ (see Appendix A) and solve it using a commercial solver (e.g., the Gurobi). In the $\left({\mathbf{AASP}}_{\mathbf{R}}\right)$ formulation, the descent routes and the execution orders at each waypoint for all aircraft in the sequence ${\Pi }^{NR}$ are kept unchanged (see constraints (A1c)–(A1e)). However, for the aircraft in the sequence ${\Pi }^R$, the routing and sequencing constraints are relaxed, allowing for potential adjustments (see constraints (A1f)–(A1i)). As a result, solving the $\left({\mathbf{AASP}}_{\mathbf{R}}\right)$ would involve re-routing, re-sequencing, and re-scheduling for the aircraft in ${\Pi }^R$, while only re-scheduling would occur for the aircraft in ${\Pi }^{NR}$. The sequence is updated if an improved solution is found, and the search continues.

Note that, while optimization solvers often require significant time to find the optimal route, sequence, and schedule for large instances, they display high efficiency in small-sized instances. Therefore, by relaxing the routing and sequencing constraints for a select few aircraft, one can leverage the solver’s capabilities and achieve an optimal or nearly optimal solution for the relaxed aircraft.

Figure 5 shows an example of the ${\mathcal{N}}_1$ for 10 aircraft (labeled a to j) with 3 descent routes (1 to 3). The disjunctive arcs, represented by dashed orange lines, indicate that these aircraft are subject to possible re-routing, re-sequencing, and re-scheduling operations. On the other hand, the remaining aircraft in the sequence follow conjunctive arcs represented by solid black lines. These aircraft are only subject to re-scheduling operations, meaning their routing and sequencing remain unchanged.

The neighborhood structure ${\mathcal{N}}_2$ involves sequentially swapping the location of a randomly selected center $c{c}_2$ with the positions within a radius $c{r}_2$ for a given sequence $\Pi$. By applying the swap operation to the current sequence $\Pi$, a new sequence ${\Pi }^{\prime }$ is generated. To evaluate its quality, the AASP problem is reformulated using ${\Pi }^{\prime }$ and converted into a problem denoted as $\left({\mathbf{AASP}}_{\mathbf{S}}\right)$ (see Appendix B). This reformulation captures the changes in the sequence and enables efficient solving using a commercial solver such as Gurobi, with the time limit parameter set appropriately. During the exploration process, only swaps that result in an improved solution are accepted and retained. An example of the swap operation in neighborhood structure ${\mathcal{N}}_2$ is illustrated in Figure 6.

The neighborhood structure ${\mathcal{N}}_3$ involves changing the current descent route for each aircraft within a sub-sequence ${\Pi }^C$ to a randomly selected alternative descent route, if one exists. The sub-sequence ${\Pi }^C$ is randomly selected from a given sequence $\Pi$ using a center $c{c}_3$ and a radius $c{r}_3$. A new sequence ${\Pi }^{\prime }$ is generated by applying the change operation to the current sequence $\Pi$ within this neighborhood. Then, we obtain a solution from solving the $\left({\mathbf{AASP}}_{\mathbf{S}}\right)$ with the fixed sequence ${\Pi }^{\prime }$. Only the change operations that improve the current solution are accepted. Figure 7 shows an example of a change operation within neighborhood structure ${\mathcal{N}}_3$.

Within the local search phase, a TS [52] approach is used to avoid being trapped in local optima and revisiting the same solution, as described in Algorithm 2. The TS efficiently employs a memory structure to move away from previous solutions. To facilitate this process, we introduce two operations: the tabu memory filter and the tabu memory update. The former operation eliminates all of the tabu moves from the current neighborhood (no aspiration criteria is used), while the latter operation updates the list of tabu moves during the search. The tabu list ($TL$) utilized in the local search phase has a maximum of $T{L}_{max}$ moves, with each implemented move storing two pieces of key information: the current sequence and the value of the central parameter. It is important to note that in some cases, the application of neighborhood structures ${\mathcal{N}}_2$ and ${\mathcal{N}}_3$ may result in infeasible solutions for problem $\left({\mathbf{AASP}}_{\mathbf{S}}\right)$. When this occurs, these operations are skipped, and the VNS algorithm proceeds to the next process.

We also note that the VNS algorithm only uses the neighborhood structure ${\mathcal{N}}_1$ in its shaking phase. The parameters $c{c}_1$ and $c{r}_1$ within ${\mathcal{N}}_1$ demonstrate variable values between the shaking and local search phases. In the shaking phase, we employ ${\mathcal{N}}_1$ to infuse diversification into our search, thereby facilitating an exploration of novel sectors within the solution space. In the local search, however, we deploy ${\mathcal{N}}_1$ to intensify the search process, thereby capitalizing on the incumbent solution to realize superior outcomes.

3.3. Rolling Horizon Approach

The rolling horizon approach is a dynamic strategy implemented in our decision support framework to handle real-time changes in the air traffic flow management. As illustrated in Figure 8, it involves breaking down the overall decision-making process into smaller decision time windows, each with a duration of $W_{decision}$. At the start of each time window, the decision support framework is executed to determine the optimal descent operation for each aircraft within the given time window based on the available information at that time. After completing the decision-making process for the current decision time window, the rolling horizon approach moves to the next time window with a time difference of $R_{roll}$. The rolling horizon approach iterates through the decision-making procedure for each subsequent decision time window until the end of the overall traffic decision. This iterative process ensures that the decision support framework remains responsive to the evolving traffic conditions and provides up-to-date recommendations for the optimal descent operations of the arriving aircraft.

In the rolling horizon approach, each decision time window consists of two groups of aircraft: ongoing aircraft from the previous window and new aircraft arriving in the current window. The ongoing aircraft’s descent processes continue from the previous window, while their descent routes and arrival times serve as constraints for the new aircraft. The decision support framework considers these constraints to optimize the descent operations for both groups, ensuring safe separation between them. This approach allows for efficient and safe management of the entire set of aircraft within each decision time window.

4. Experimental Results

This section presents the numerical experiments conducted to evaluate the performance of the proposed decision support framework and the VNS algorithm. We implemented the pseudospectral technique to solve the ATOP problem using the Matlab R2021a programming language, and implemented the VNS algorithm to solve the AASP problem using the Python 3.9 programming language. The problems $\left(\mathbf{AASP}\right)$, $\left({\mathbf{AASP}}_{\mathbf{I}}\right)$, $\left({\mathbf{AASP}}_{\mathbf{R}}\right)$, and $\left({\mathbf{AASP}}_{\mathbf{S}}\right)$ are solved using the commercial solver Gurobi 9.5.0 [53] with default parameters, except for the time limit. All experiments are performed on a PC with an Intel(R) Core(TM) i7-10710U [email protected] ARM 16G under the Windows 10 operating system.

4.1. Test Instances and Parameters Setting

4.1.1. Traffic Instances in Gbia

GBIA is the largest transport hub in China and serves as the primary hub for China Southern Airlines. According to the published Aeronautical Information Publication (AIP), Figure 9a displays the designed north orientation arrival routs utilized to approach the hub by incoming aircraft. Figure 9b highlights the path stretching routes for the ps-CDO procedure. We consider an E-TMA with a radius of 220 NM around the GBIA. The GBIA features three parallel runways expressly designed for arriving and departing procedures: runway 01/19, runway 02L/20R, and runway 02R/20L, all of which are shown in Figure 9c. Among them, Runway 01/19 is operational for both arrivals and departures, runway 02R/20L is exclusively reserved for arrivals, and runway 02L/20R is solely used for departures. This study assumes that the independent instrument approaches between runway 01/19 and runway 02R/20L, and the segregated parallel approaches and/or departures between runway 02R/20L and runway 02L/20R.

In our study, we conducted simulations using actual aircraft routing and scheduling data from GBIA TMA for a full day in November 2019. The test data set consists of 1290 aircraft, with 653 arrivals and 637 departures, including three types of aircraft: 267 Airbus A332 (heavy), 452 Airbus A320 (medium), and 571 Boeing B738 (medium). For all arrival aircraft, we assume a north orientation approach. To analyze different air traffic scenarios, we divided the day into three distinct time periods: low-density, normal-density, and high-density. Figure 10 illustrates the hourly movements during the simulated day, with the three highlighted periods representing the identified air traffic situations. We created six separate test instances based on these time periods. The instances labeled A1, A2, and A3 considered only arrival traffic, while the instances labeled AD1, AD2, and AD3 included both arrival and departure traffic.

Our study focuses on analyzing the effects of arrival traffic exclusively within the E-TMA phase while disregarding enroute fluctuations. We set the arrival time at the entry waypoint identical to the estimated target arrival time for every arriving aircraft. For departure traffic, a designated runway is assigned for takeoff within a given time window. We set the departure window at 900 s, and assume the minimum fuel-burning take-off time is equal to the estimated target take-off time. The take-off air segments from runway exit for the departure aircraft is not considered in this paper. Therefore, the runways are the only potential points of conflict between the arrival and departure traffic.

4.1.2. Parameters Setting

In our numerical experiments, we computed the descent trajectories for every arriving aircraft, starting from the last part of the cruise prior to the TOD and ending at the runway. The aircraft maintain a cruise altitude of 35,000 feet until the TOD, where they enter the descending phase. The runway altitude is set to 50 feet. We refer to [11,48] for more details on altitude and speed constraints at each air segment in the ATOP for our three analyzed descent operations, including the SDO, the CDO, and the ps-CDO. The altitude-dependent wind profile $W_s$ (h) used in our ATOP is taken from the meteorological terminal aviation routine weather report of the GBIA. The remaining performance parameters of the aircraft, consisting of the aerodynamic data, the fuel flow data, and the geometric data, are available in BADA.

In the AASP, the minimum time-based separation between two consecutive aircraft in the airspace is calculated by incorporating their angles and speeds. The calculation algorithm is based on the work of [2]. The minimum time-based separation matrices for operations on the same runway or close parallel runways are provided in

Table 2. The minimum time-based separation matrix for operations on the same runway (unit: s).

Preceding	AH	AM	AL	DH	DM	DL
AH	96	157	196	75	75	75
AM	60	69	131	75	75	75
AL	60	69	82	75	75	75
DH	60	60	60	90	120	120
DM	60	60	60	60	60	60
DL	60	60	60	60	60	60

[54] and

Table 3. The minimum time-based separation matrix for operations on the close parallel runways (unit: s).

Preceding	Trailing
	AH	AM	AL	DH	DM	DL
AH	-	-	-	68	68	80
AM	-	-	-	62	62	80
AL	-	-	-	48	55	80
DH	54	58	80	-	-	-
DM	54	58	80	-	-	-
DL	54	58	80	-	-	-

Note: Symbol “-” indicates that no valid operation combination is possible.

[55], respectively. Aircraft involved in arrivals and departures are symbolized by “A” and “D”, whereas aircraft weight categories are represented by “H”, “M”, and “L” for heavy, medium, and light aircraft, respectively. For example, “AH” denotes the arrival of a heavy aircraft.

Table 4. Values of the parameters for the VNS algorithm.

Parameters	Number of Aircraft
	$\mathbf{\left|}\mathcal{F}\mathbf{\right|}\mathbf{\le }\mathbf{30}$	$\mathbf{30}\mathbf{<}\mathbf{\left|}\mathcal{F}\mathbf{\right|}\mathbf{\le }\mathbf{70}$	$\mathbf{\left|}\mathcal{F}\mathbf{\right|}\mathbf{>}\mathbf{70}$
${iter}_{max}$	3	5	$max\{5,[\frac{\left|\mathcal{F}\right|}{40}]\}$
${noimproving}_{max}$	2	2	$max\{3,[\frac{\left|\mathcal{F}\right|}{25}]\}$
$c{r}_1$ (shaking)	2	4	$min\{[\frac{\left|\mathcal{F}\right|}{20}],6\}$
$c{r}_1$ (local search)	4	6	$min\{[\frac{\left|\mathcal{F}\right|}{15}],10\}$
$c{r}_2$	2	2	$max\{3,[\frac{\left|\mathcal{F}\right|}{35}]\}$
$c{r}_3$	2	2	$max\{3,[\frac{\left|\mathcal{F}\right|}{35}]\}$
${timelimit}_{initial}$ (s)	1	1	10
${timelimit}_{LS}$ (s)	1	1	3
$T{L}_{max}$	10	15	$max\{15,[\frac{\left|\mathcal{F}\right|}{10}]\}$

Note: ${timelimt}_{initial}$ and ${timelimt}_{LS}$ refer to the time limits set for the Gurobi solver during the initial solution generation phase and the local search phase, respectively.

lists a summary of the parameter values used in the VNS algorithm. We classify all test instances in Table 4 into three classes based on the number of aircraft: up to 30 aircraft, 30 to 70 aircraft, and more than 70 aircraft. We solve all test instances within each group using these parameter values.

Table 5. Values of the unit delay cost (€/s).

Delay (s)	Arrivals			Departures
	A332	A320	B738	A332	A320	B738
(0, 300]	1.25	0.67	0.67	0.91	0.4	0.4
(300, 900]	1.64	0.89	0.91	1.29	0.63	0.64

displays the delay cost per unit $c_{delay}$, which is calculated based on [56] full tactical costs but excludes the fuel cost incurred by the delay. The fuel cost per unit $c_{fuel}$ is set to 0.8 €/kg, in line with the base scenario in [56].

4.2. Framework Decision Solutions

In this section, we investigate the performances of our proposed framework with respect to three descent operations (i.e., SDO, CDO, and ps-CDO). We arbitrarily select the results with respect to $\lambda =0.5$.

Table 6. The computational results on instances A1–AD3.

Indicators		A1	A2	A3	AD1	AD2	AD3
SDO	Total cost (€)	12,842.31	58,119.74	75,851.6	12,979.91 (12,882.92)	63,028.02 (61,399.59)	80,465.05 (77,989.02)
	Total delay cost (€)	547.89	11,018.2	12,392.15	678.67 (600.07)	15,839.66 (14,522.25)	17,510.85 (15,508.58)
	Total fuel consumption cost (€)	12,294.42	47,101.54	63,459.45	12,301.24 (12,282.85)	47,188.36 (46,877.34)	62,954.2 (62,480.44)
	Total delay (s)	655	11,984	13,480	807 (687)	18,200 (15,915)	19,991 (16945)
	Total difference (s)	1222	13,655	20,651	1356 (1236)	14,875 (12,590)	21,044 (17,998)
	Total fuel consumption (kg)	15,368.03	58,876.92	79,324.31	15,376.55 (15,353.57)	58,985.45 (58,596.68)	78,692.75 (78,100.55)
CDO	Total cost (€)	11,483.37	-	-	11,589.12 (11,485.14)	-	-
	Total delay cost (€)	1554.65	-	-	1618.84 (1534.64)	-	-
	Total fuel consumption cost (€)	9928.72	-	-	9970.28 (9950.5)	-	-
	Total delay (s)	1833	-	-	1996 (1862)	-	-
	Total difference (s)	969	-	-	1174 (1040)	-	-
	Total fuel consumption (kg)	12,410.89	-	-	12,462.85 (12,438.13)	-	-
ps-CDO	Total cost (€)	11,675.87	56,729	74,834.28	11,902.97 (11,799)	63,362.42 (61,272.22)	81256.84 (78,649.21)
	Total delay cost (€)	1444.2	15,652.89	20047.09	1657.49 (1573.29)	22,109.76 (20,401.98)	26,463.63 (24,355.8)
	Total fuel consumption cost (€)	10,231.67	41,076.11	54,787.19	10,245.48 (10,225.71)	41,252.66 (40,870.24)	54,793.21 (54,293.41)
	Total delay (s)	1798	16,566	21,115	2013 (1879)	23,138 (20,242)	28545 (25485)
	Total difference (s)	1109	10,078	14,264	1248 (1114)	12,284 (9388)	15,866 (12,806)
	Total fuel consumption (kg)	12,789.59	51,345.14	68,483.99	12,806.86 (12,782.14)	51,565.82 (51,087.79)	68,491.51 (67,866.76)

Note: Symbol “-” indicates that the aircraft scheduling problem obtains no feasible solution. A number within each parenthesis indicates the value for the arrival aircraft.

displays the computational results for instances A1–AD3. In Table 6, the rows are organized in blocks of six rows for the three descent operations, including the values of the total cost, the total delay cost, the total fuel consumption cost, the total delay, the total difference, and the total fuel consumption.

Based on the analysis of total cost, the results show that different descent optimization methods yield optimal solutions depending on the specific scenarios and traffic densities. For instances A1 and AD1, which represent low traffic density scenarios, the CDO approach is found to be the optimal choice. For instances A2 and A3, which represent normal and high traffic density scenarios considering only arrival traffic, the ps-CDO approach proves to be the optimal choice. In scenarios that involve both arrivals and departures, such as instances AD2 and AD3, the SDO approach is identified as the best choice.

Note that the implementation of CDO in instances A2, A3, AD2, and AD3, which represent scenarios with dense traffic, is found to be infeasible. In such cases, ensuring safety while implementing CDO becomes challenging due to the need for extra separation buffers. As a result, arrival aircraft are unable to follow the CDO trajectories in these scenarios. With the path stretching approach, on the other hand, each aircraft approaching based on the ps-CDO is provided with various options for its lateral path with respect to different arrival times. An AASP can thus be solved to properly sequence the aircraft in order to avoid aircraft conflicts. However, the inclusion of path stretching within the ps-CDO leads to increased travel distances for aircraft during descent, resulting in higher fuel consumption and reduced travel time efficiency.

While the SDO approach performs well in terms of reducing total delay compared to the CDO and ps-CDO approaches, it leads to increased fuel consumption for arrival aircraft. On the other hand, the CDO-based approach demonstrates a fuel consumption reduction of approximately 15% compared to the SDO approach. As traffic flows increase and TMA capacity remains constant, total delay and its associated costs significantly increase. However, the increase in total fuel consumption and associated costs is negligible. The SDO approach, which reduces total delay, is thus more advantageous in terms of total cost and recommended for instances AD2 and AD3, which involve both arrival and departure traffic and higher traffic flows.

These findings highlight the trade-offs and considerations associated with different descent optimization approaches in terms of safety, fuel consumption, travel time efficiency, and delay reduction, particularly in scenarios with dense traffic. The choice of the optimal approach depends on the specific operational requirements and priorities in each scenario.

4.3. Sensitivity Analysis on Weight Parameter

In this section, we compare the performance of the SDO against the ps-CDO with varying values of the weight parameter $\lambda$ within our framework. We increase the weight parameter $\lambda$ from 0 to 1 with a step size of 0.1. Figure 11 illustrates the effect of altering the weight parameter $\lambda$ on the total delay and fuel consumption of the SDO (solid lines) and the ps-CDO (dashed lines) for instances A1-AD3. The corresponding cost results of the SDO and the ps-CDO are presented in Figure 12, including the total cost, the total delay cost, and the total fuel consumption cost.

As the weight parameter $\lambda$ increases, the importance of total delay in the AASP model increases, resulting in solutions with smaller total delay for both the SDO and the ps-CDO approaches. However, total fuel consumption and its associated costs gradually increase. The ps-CDO approach outperforms the SDO in terms of total fuel consumption and its associated cost, while the SDO performs better in terms of total delay and its associated costs, regardless of the value of $\lambda$.

In terms of total cost, during low or normal traffic scenarios (instances A1, A2, A3, and AD1), the ps-CDO consistently yields better total cost results compared to the SDO. However, during high traffic scenarios (instances AD2 and AD3), the SDO generates smaller total cost compared to the ps-CDO when $\lambda =0.5$ or $0.6$ in AD2 and $\lambda \ge 0.5$ in AD3, respectively. This detailed trade-off information can assist ATCs in determining the optimal compromise between different performance indices using our framework.

4.4. Effectiveness of the VNS Algorithm

In this section, we evaluate the effectiveness of our proposed VNS algorithm for the AASP model. Here, we arbitrarily set $\lambda =0.5$ and provide the findings in

Table 7. The computational results for the AASP on instances A1–AD3.

Index			A1	A2	A3	AD1	AD2	AD3
SDO	Gurobi	Obj. (s)	938.5	12,819.5	17,065.5	1081.5	16,706.5	20,590.5
		CPU time (s)	0.16	1800.1	1800.18	0.14	1800.21	1800.23
		Gap (%)	0	0.11	0.08	0	12.61	5.44
	Our VNS	Obj. (best) (s)	938.5 (10)	12,819.5 (10)	17,065.5 (10)	1081.5 (10)	16,537.5 (2)	20,517.5 (5)
		Obj. (avg.) (s)	938.5	12,819.5	17,065.5	1081.5	16,667.1	20,554
		CPU time (avg.) (s)	0.34	24.83	29.42	0.56	106.4	121.55
ps-CDO	Gurobi	Obj. (s)	1453.5	13,322	17,689.5	1630.5	17,713	22,205.5
		CPU time (s)	0.11	1800.07	1800.12	0.12	1800.11	1800.23
		Gap (%)	0	1.56	1.11	0	15.74	10.16
	Our VNS	Obj. (best) (s)	1453.5 (10)	13,322 (10)	17,689.5 (10)	1630.5 (10)	17,711 (8)	22,205.5 (10)
		Obj. (avg.) (s)	1453.5	13,322	17,689.5	1630.5	17,711.4	22,205.5
		CPU time (avg.) (s)	0.32	26.6	33.54	0.41	93.01	124.22

Note: A number within each parenthesis indicates the frequency of obtaining the best solution among the 10 runs of the VNS algorithm.

, in which we compare the VNS algorithm against the Gurobi solver with a time limit of 1800 s. The first three rows of the table list the outcomes of the Gurobi for the SDO, including the best objective function value, CPU time, and optimality gap at the time limit. Rows 4–6 present the results of our VNS algorithm for the same instance with the SDO run ten times, including the best objective function value, the average objective function value, and the average CPU time. Rows 7–12 provide the statistics for the ps-CDO. The number within each parenthesis indicates the frequency of obtaining the best solution among the 10 runs of the VNS algorithm.

Our VNS algorithm performed significantly better than the Gurobi solver in terms of solution quality. In all instances, our algorithm produced solutions that were either equal to or better than the best solutions obtained by Gurobi. Notably, the average solution quality achieved by our VNS algorithm was comparable to the best solutions generated by Gurobi.

The VNS algorithm can produce near-optimal solutions for all instances within a few seconds, whereas the Gurobi algorithm fails to provide better solutions within the specified time limit for all instances except A1 and AD1. Furthermore, the reliability of our VNS algorithm is evident as it consistently delivers the best solutions after multiple runs, while maintaining reasonably short CPU times. As the number of aircraft increases, resulting in more variables and constraints, the iteration times naturally increase. However, the CPU time of our VNS algorithm does not increase proportionately, showcasing its efficiency in solving the AASP problem.

4.5. Decision Solutions for Daily Operations

In this section, we employ the rolling horizon approach to explore the decision solutions of our decision support framework over a full day. The complete day instances can be found in Section 4.1.1 (see Figure 10). We consider two cases: (i) arrival traffic only, and (ii) both arrival and departure traffic. The rolling horizon method is employed with parameter settings of $W_{decision}=2$ h and $R_{roll}=1$ h. The results are analyzed based on a weight parameter of $\lambda =0.5$. Figure 13a shows the optimal descent operation solutions for arriving aircraft hourly throughout the day in case (i), while the outcomes for (ii) are presented in Figure 13b.

The effectiveness of the CDO procedure for arriving aircraft is most prominent during periods of low traffic density, specifically between 02:00 to 09:00 in both case (i) and case (ii). However, its practicality diminishes during the remaining periods of the day, as it fails to satisfy the minimum time-based separation requirement in the TMA. Implementing the CDO during these times would render the AASP problem infeasible. Regarding the ps-CDO, in case (i), it is the optimal choice for arriving aircraft between 09:00–14:00 and 15:00–24:00. For case (ii), it remains optimal during periods of normal traffic density, namely 09:00–12:00, 13:00–14:00, 15:00–18:00, and 21:00–22:00. However, during high traffic density periods, following the ps-CDO would lead to significantly increased delay costs, making it less favorable in terms of total cost. Thus, it is advisable to utilize the SDO procedure during high traffic density periods, as it minimizes the total cost. It should be noted that the ps-CDO cannot be utilized from 00:00–01:00, as it cannot provide conflict-free trajectories for arriving aircraft during this time.

5. Conclusions

In this paper, we present an innovative collaborative decision support framework for future trajectory-based operations that integrates aircraft trajectory optimization, arrival scheduling, optimal trajectory selection, and their interactive processes within a single framework. First, the framework is capable of computing the aircraft descent trajectories by performing three distinct descent operations (SDO, CDO, and ps-CDO). Second, based on the descent trajectories, the framework allocates a conflict-free trajectory to each aircraft while minimizing the total delay and total fuel consumption. Finally, the framework selects the best descent operation out of the three considering minimizing the total cost of total delay and total fuel consumption. We evaluate the performance of our framework for traffic handling in the GBIA TMA, considering both arrival and departure traffic. The computational results demonstrate that our framework can support the ATCs in efficiently managing the descent operations in the TMA. Additionally, our VNS algorithm featuring three novel neighborhood structures proves to be efficient for solving the AASP. Therefore, we propose that this framework will serve as a technical enabler for future air traffic flow management.

One avenue for future research is to extend our decision support framework to include the point merge system [57] and vectoring maneuvers [58]. Our proposed framework heavily relies on accurate trajectory planning, yet uncertainties such as weather [59], navigation accuracy, and pilot operations may undermine its performance and pose potential safety risks. Addressing these uncertainties in future research would be valuable for enabling us to develop a more robust decision support framework for realistic scenarios.