1. Introduction
With the demand for air travel increasing after the COVID-19 pandemic, the environmental impacts caused by air traffic have become increasingly serious. Previous research [
1] has suggested that aviation accounts for 2.5% of global CO
2 emissions. Nevertheless, it has been demonstrated to contribute approximately 4% to the global warming observed to date. For the sustainable development of the aviation and commercial aircraft manufacturing industries, many new technologies that can reduce carbon emissions have been extensively explored. Among the numerous measures aimed at reducing CO
2 emissions from aviation, enhancing the environmental performance of air traffic systems is currently regarded as the most viable, significant, and promising approach.
In the National Airspace System, a Terminal Maneuvering Area (TMA) is a complex area of airspace surrounding one or more metropolitan airports, designed to aggregate arriving traffic flows from different directions. To increase the capacity of TMAs, an arrival route structure based on area navigation, called the Point Merge System (PMS), has been developed to assist controllers in sequencing and reducing open-loop radar vectoring. A PMS is a method for delaying aircraft involving less radar vectoring. Its structure includes a merge point and equidistant arc-shaped sequencing legs to the merge point. The flight path in a PMS can be lengthened or shortened as needed without radar vectoring. In the future Global Air Navigation Plan of ICAO [
2], the PMS is recognized as having the potential to reduce environmental impacts. However, a well-defined, scientific emissions benchmark for the minimum harmful input under a given air traffic flow is still lacking. Therefore, a scientifically accurate benchmark for measuring environmental performance is needed. With such a benchmark, low-carbon aircraft scheduling techniques can be applied in a step-by-step manner.
In order to better understand the environmental impacts of TMAs, both at present and in the future, one approach being adopted by Air Traffic Management (ATM) operators is quantifying their performance using measures of flight inefficiency [
3]. Flight inefficiency can be defined as anything that causes an aircraft to fly a path different from its fuel-optimal four-dimensional trajectory. Therefore, the accuracy of the optimized flight trajectory is important. A challenge posed in the evaluation of the environmental impact of air traffic, both at present and in the future, is determining the optimal trajectory and minimal fuel burn (which is the benchmark for fuel-based inefficiency analysis) for any given route, especially for arrival routes within a TMA.
This research encompasses two technical aspects. First, the arrival traffic flow with a point merge pattern is identified from the cluttered trajectories within the terminal area, and prevalent arrival routes are further extracted. Second, as continuous descent operations are recommended by the ICAO, descent profile optimization provides an important basis for CO2 emission benchmarking studies. A literature survey on these two technical aspects is provided below.
The implementation of new airspace surveillance technologies, such as Automatic Dependent Surveillance-Broadcast (ADS-B), has led to improvements in the integrity, quality, and accessibility of aircraft tracks. The analysis of track records is becoming a viable method for identifying traffic patterns and assessing the operational performance of air traffic. At present, trajectory clustering is becoming the most effective method for achieving these goals. Various distance-based clustering algorithms have been reported, but the core component of all of them is the similarity or distance measure used for the classification of data [
4].
In order to accurately represent the overall similarity between two trajectories, various spatial distance-based similarity metric models have been proposed. Almost all of the typical clustering algorithms can be used, either directly or synthetically, to perform trajectory clustering using pre-computed similarity values [
5]. Besse et al. summarized the calculation of spatial distances of trajectories for clustering, giving the distances between sets of points, between points and line segments, and between line segments. This new distance was compared to other metrics, according to the corresponding clustering results obtained using both hierarchical clustering and affinity propagation (AP) methods [
6]. Poppe and Buxbaum clustered aircraft climb trajectories using the k-means algorithm, where the resulting clusters could further be used in trajectory prediction applications [
7]. The calculated point-to-segment distance matrix—a spectral clustering algorithm—has been adopted to cluster flight trajectories [
8].
Given the limitations of distance measures for high-dimensional spatio-temporal trajectories, researchers have attempted to improve the credibility of trajectory similarity through dimensional reduction. For example, Olive and Morio used the Douglas–Pecuker algorithm to approximate the original trajectories according to a fixed step size [
9]. Xuhao et al. revealed the time-series nature of the trajectories, converted the trajectories into a time-series of headings relative to their remaining distances, and calculated the similarity between them using the dynamic time warping algorithm [
10]. In recent years, the continuous development of Deep Learning approaches has demonstrated their powerful ability to extract internal potential features from large data sets. Deep autoencoders have been developed to train deep neural networks in order to learn suitable representations from high-dimensional trajectory data. Subsequently, the output of the Deep autoencoder is input into a GMM or DNSCAN for clustering [
11]. In order to understand the spatial distribution characteristics of the arrival traffic flow in a TMA, Wang established a flight trajectory feature compression model on the basis of the t-SNE method, following which the density peaks clustering approach was used to extract the prevailing traffic flow [
12].
Trajectory clustering plays a crucial role in air traffic management, providing insights for airspace structure optimization, trajectory planning, and trajectory prediction. The application of trajectory clustering in the field of traffic pattern recognition within TMAs is evolving and expanding, including enhancement of the airspace structure around airports [
13], examination of the reliability of arrival and departure route networks [
14], and analysis of the efficiency of air traffic operations [
9]. However, due to the structural complexity of the point merging procedure and the eventual convergence of arrival traffic, which causes flight trajectories to be more mixed, aggregated, and intertwined, identification of traffic flows with PMS topology and its representative trajectories may require a new similarity-measuring model that takes into account the complex features of arrival traffic trajectories.
Trajectory optimization plays a key role in ensuring the success of trajectory-based operations. The optimization problem for 4D aircraft trajectories has been studied extensively by many scholars, who have gradually realized that the 4D trajectory optimization problem can be formulated as a multi-phase optimal control problem with process constraints, terminal points constraints, and control constraints [
15,
16]. At present, 4D trajectory optimization methods can be classified into four categories: dynamic programming, indirect, direct, and meta-heuristic algorithms.
Initially, researchers explored the potential of dynamic programming methods to address trajectory optimization problems [
17]. However, it was found that the computational complexity and issues associated with interpolation limited the applicability of these methods in practical trajectory optimization scenarios [
18]. Many scholars have attempted to address the aforementioned limitations of dynamic programming methods [
19,
20,
21]. For example, Ahmed et al. [
22] recently proposed a modified dynamic programming method which overcomes the curse of dimensionality and the extended mesh problem while enhancing the computational efficiency and feasibility of 4D trajectory optimization.
Indirect methods are based on transforming the 4D trajectory optimal control problem into a Hamiltonian boundary value problem, which is then solved using numerical methods to obtain differential or algebraic equations [
23,
24]. The primary application of this technique is to optimize relatively simple systems [
25]. Indirect methods encounter difficulties in solving the trajectory optimization problem due to the complex model dynamics and the numerous constraints involved in the 4D trajectory optimization problem [
26].
Compared with indirect methods, direct methods discretize the complex optimal control model into a finite-dimensional nonlinear planning problem, which has higher computational efficiency and numerical stability, and is suitable for optimal control problems with complex constraints [
27,
28]. Soler et al. [
29] have transformed the 4D trajectory optimization problem into a multi-stage optimal control problem and used the point collocation method to obtain a 4D trajectory with the objective of minimum fuel consumption. Based on this study, Soler et al. [
30] further analyzed the performance of aircraft, considering different take-off weights and cost indices to obtain a 4D trajectory with minimum fuel consumption. Bonami et al. [
31] transformed the 4D trajectory optimization problem into a mixed-integer nonlinear programming problem and used the branch and bound method to find a 4D trajectory with the objective of optimal fuel consumption.
Continuous Descent Operations (CDOs) have been recognized as an effective way to significantly reduce fuel burn and CO
2 emissions in aircraft operations [
32]. By keeping arriving aircraft at their cruise altitude longer, then following a continuous descent at near-idle thrust with no level-flight segments, CDO can bring substantial fuel savings and environmental benefits [
33]; therefore, it has become a globally recommended emission mitigation practice by the International Civil Aviation Organization (ICAO).
The pseudo-spectral method is a representative direct method that transforms the continuous optimal control problem into a nonlinear programming (NLP) problem, enabling more efficient parameterization of constraints. Many scholars have employed pseudo-spectral methodologies to address the descent profile optimization problem, in accordance with specified horizontal route conditions. For instance, Park and Clarke divided a CDO into a number of descent stages according to different flap settings and calibrated airspeed constraints, then modeled the trajectory optimization problem as a multi-phase optimal control problem between two previously known endpoints, and solved for the optimal continuous descent trajectory using the Gaussian pseudo-spectral method [
34]. Subsequently, Park et al. employed the Legendre–Gauss pseudo-spectral method to determine the optimal commencement of the descent point and descent profile, and proposed a CDO trajectory that simultaneously reduces flight time and minimizes fuel consumption [
35]. Similarly, Ma et al., taking delay time, fuel burn, and aircraft noise as objectives, used the Legendre–Gauss pseudo-spectral method to generate optimized trajectories [
36]. Zhu et al. developed trajectory optimization models based on the Gaussian pseudo-spectral method to obtain optimal descent trajectories for different cost index values [
37]. Most of the above studies relied on the GPOPS toolbox for solutions. However, the quality of the solution and the convergence of the computational process depends heavily on the initial estimates of the control variables, and inappropriate values can lead to solution failure [
38].
The advent of swarm intelligence optimization technology has prompted some scholars to employ heuristic algorithms in an attempt to resolve the optimal control problem of 4D trajectories. Yang et al. addressed the challenge of not being able to provide time-continuous control variables in trajectory optimization problems by transforming the optimal control problem into a nonlinear planning problem [
39]. The nonlinear planning model was constructed using the thrust and load factor as decision variables, which was solved using IPSO-SQP [
40]. They claimed that, despite its simplicity and intuitiveness, the particle swarm methodology proved to be quite effective in finding the optimal solution to orbital rendezvous optimization problems with considerable numerical accuracy. Ren et al. developed an integer planning model based on flight segment division and employed a particle swarm optimization algorithm to obtain optimal green trajectories [
41].
Although the aforementioned studies have achieved notable success in their respective fields, they have not yet addressed the CO2 emissions baseline trajectory problem of traffic flow. In order to investigate the mechanism driving the environmental inefficiency of point merge traffic flows within a busy terminal area, we combine trajectory clustering and optimal control methods to propose a 4D trajectory optimization model for benchmarking the CO2 emissions of inbound traffic flows. This is the first time that such a systematic approach has been utilized to gain insight into the emissions benchmarking of traffic flows from the perspective of a data-driven optimization model.
The differences between the research in this study and the existing achievements are specifically reflected in three aspects: First, the continuous descent operation under the special structure of the PMS presents distinct constraints from traditional arrival procedures. To satisfy these constraints, we propose a new optimal control model for the 4D arrival trajectory optimization problem. Second, a 4D trajectory nonlinear optimization method is developed on the basis of a genetic algorithm to solve the 4D trajectory optimal control problem. The method is able to model and solve the PMS vertical profile encompassing a level-flight arc leg in an ad hoc manner. Third, a spherical segment–path distance between trajectories directly using their coordinates is calculated, which overcomes the deformation defects caused by the indirect use of map projection coordinates.
To establish environmental benchmarks for arrival traffic flows utilizing PMS topology, this study introduces a 4D trajectory optimization method that integrates data-driven and optimal control models. The predominant arrival routes are identified via trajectory spectral clustering. A tailored optimal control model, specifically designed for the vertical profiles of PMS topology, is developed with the objective of minimizing fuel–time costs, using the identified arrival routes as horizontal guidance. To overcome the limitations of traditional solvers like GPOPS in addressing optimal control problems involving multiple descent processes, a genetic algorithm-based solution approach is proposed. Subsequently, the proposed method was applied to a realistic terminal airspace scenario, and its applicability was validated through two consecutive experiments: one for identifying arrival traffic flows and the other for optimizing 4D trajectories.
2. Methodology
The section outlines the process of data pre-processing, point merge operation flow identification, and continuous descent operation 4D trajectory generation, as illustrated in
Figure 1. The process consists of five steps: (1) extract arrival trajectories from the traffic flow records within the boundaries of the terminal control area; (2) calculate the distances between pairs of trajectories sequentially to obtain the spherical symmetric segment distances, and construct the similarity matrix of the arrival trajectory set; (3) use a spectral clustering algorithm to obtain the spatial pattern of arrival traffic flow, identify the core trajectories for the identified traffic flow and, finally, take the horizontal profile of the core trajectory as the prevailing arrival route; (4) establish an optimal control model for a CDO profile with PMS topology based on predominant arrival routes; and (5) solve the optimal control problem for vertical profiles using a nonlinear trajectory optimization model—called GACDO—based on a genetic algorithm.
2.1. Data Pre-Processing
Air traffic surveillance equipment, such as secondary surveillance radar and ADS-B, are capable of acquiring many kinds of aircraft navigational traffic data, such as the position, altitude, speed, and heading of aircraft in the terminal area and their changes, in real-time. The Air Traffic Control (ATC) automation system receives, processes, and displays the surveillance data and, at the same time, records all or part of the flight tracks of the aircraft.
The concepts of track and trajectory are not distinguished in ATC operations and research and, for the purpose of this study, they are defined as follows:
A basic track is the geographic position and altitude information of the kth point of an aircraft i recorded by a surveillance sensor at a certain time t, expressed as a four-dimensional vector , where , , and denote the longitude, latitude, and altitude of the aircraft in space, respectively, and t denotes the timestamp.
The trajectory is the historical trace of an aircraft’s flight over a period of time, which is actually a continuous curve. Due to the data processing cycle of monitoring sensors and the cognitive limitations of air traffic controllers, the trajectory T consists of a series of discrete tracks arranged in chronological order with certain time intervals; that is, , where k ∈ [1, K] is the chronological index of tracks and K is the total number of tracks in the trajectory T.
Suppose that the set of trajectories formed by all arriving aircraft in the terminal area in a given time period is , where each denotes an arriving flight trajectory and N denotes the total number of all trajectories.
Although the trajectory data recorded by ATC surveillance systems has been processed for correlation, noise reduction, and elimination of duplicates, it still contains a large number of redundant and irrelevant tracks in the context of analyzing the environmental performance of the arriving traffic flows.
Most trajectories also contain tracks of the en-route flight, or even of the entire flight phase from takeoff to landing. It is necessary to remove those tracks outside of the terminal airspace. The geographic boundary of the terminal airspace is usually formatted as a two-dimensional polygon. For such processing, tracks outside the boundary of this polygon can be eliminated, while those inside are retained.
Tracks within the terminal airspace may also contain overflying aircraft tracks. In addition, there are general aviation flights and training flights within the terminal airspace. When general aviation aircraft are equipped with ADS-B transmitters, the surveillance system will also record these tracks. We can eliminate these tracks through determining whether the aircraft lands at the major airport in the terminal airspace. As there are different movement patterns for arriving and departing flights, they need to be distinguished. Arrival and departure trajectories can be easily separated according to their departure and destination airports in the flight plans that correspond to the trajectories.
Finally, we obtain the arrival trajectories during a period of time. The set of trajectories formed by all arriving aircraft in the terminal area over a given time period is , where each denotes an arriving flight trajectory and N denotes the total number of all such trajectories.
2.2. Traffic Flow Identification and Predominate Arrival Route Recognition
Arrival traffic flow consists of successive arrivals with similar movement behaviors. As they have similar trajectories, they are often recognized through trajectory clustering. In general, clustering can be performed through comparing trajectory objects. Therefore, it is important that a trajectory distance measurement model reflects the trajectory shape characteristics, including the physical distance and total length of the trajectory. In addition, as aircraft fly over a large area of the Earth’s surface, the trajectory distance must take into account its spherical geometry in order to avoid the deformation caused by map projection.
2.2.1. Spherical Segment–Path Distance between Trajectories
Supposing that there are two tracks,
P1(
λ1,
ϕ1,
h1,
t1) and
P2(
λ2,
ϕ2,
h2,
t2), where
λ1 and
ϕ1 represent the longitude and latitude of
P1, and
λ2 and
ϕ2 represent the longitude and latitude of
P2, the shortest distance between them can be calculated using the Haversine formula [
42], as shown in Equation (1):
where
R0 represents the average radius of the Earth. This formula takes into account the rounding error in the proximity between the two points and calculates the great circle distance between them.
If trajectory
has a segment arc
defined by its endpoints
and
, the great circle distance from any track point
O on trajectory
to arc
can be calculated with Equation (2):
Figure 2a shows the distance
from point
O to trajectory
, which is the minimum distance between the point
O and all segments that make up
. This can be calculated using Equation (3):
where
is the number of tracks contained in
. Then, the spherical segment–path distance from
to
can be defined as Equation (4):
where
Ok is the
kth track on trajectory
and the number of track points contained in trajectory
is
. A schematic representation of the calculation of
is provided in
Figure 2b.
The inter-trajectory distance matrix
D = [
di,j]
N×N is generated by calculating the symmetric unidirectional great circle distance between trajectories
and
in the set of arrival trajectories
, as shown in Equation (5).
2.2.2. Predominant Arrival Route Recognition
Spectral clustering is an algorithm derived from graph theory, the main idea of which is to treat all data as points in a high-dimensional space connected by edges. The weight of an edge between two distant points is low, while the weight of an edge between two close points is high. Clustering can then be achieved through dividing the graph of data points into sub-graphs in such a way that the sum of edge weights between the sub-graphs is minimized and the sum of edge weights within the sub-graphs is maximized.
The matrix of trajectory similarity is constructed based on inter-trajectory distances, in order to accurately reflect the similarity relationship between data points. This results in higher similarity between similar points and lower similarity between dissimilar points. Equation (6) shows the Gaussian kernel function
used to construct the similarity between trajectories
and
:
The parameter σ in Equation (6) controls the width of the neighborhood; a larger value of σ indicates a higher similarity between the trajectories, which should be adjusted according to the specific situation. A non-negative link weight unidirectional graph G can be constructed on the basis of the similarity between the trajectories. The weight of the zero eigenvalue of the Laplace matrix L of the graph G is approximately equal to the number of connected sub-graphs of the graph G, denoted as Z. It is important to note that Z represents the estimated number of clusters.
The spectral clustering method divides the original trajectory samples into Z subsets (i.e., Z different clusters of trajectories) according to their similarities in structure and proximity to each other. However, in practice, the clusters obtained from preliminary clustering do not accurately reflect the spatial characteristics of the traffic flow. This is because the distances between trajectories within the clusters are not the same, especially when the clusters contain some abnormal trajectories. Thus, it is often necessary for a cluster to be subdivided into several smaller sub-clusters. This study proposes the use of a k-Nearest Neighbor (KNN) distance model between the trajectories to identify prevalent, anomalous, and representative trajectories obtained from the same cluster. The aim is to identify the main trajectories and their representative flight routes, which can effectively express the characteristics of the spatial distribution of traffic flow.
To obtain the average KNN distance of a trajectory
in the trajectory cluster
, the inter-trajectory distance matrix
D—given by Equation (5)—is first searched to find the distances to other trajectories in
. Then, these distances are sorted from smallest to largest, and the average KNN distance
of
is calculated using Equation (7):
The k-nearest neighbor distance matrix of also can be used to classify a given cluster of trajectories into two sub-classes using k-medoids. This process yields the prevalent traffic flow trajectories and the anomalous trajectories. Each class of data obtained using the k-medoids classification algorithm has a corresponding center that serves as a representative of the data with high similarity. The prevalent traffic flow trajectories also exhibit a core trajectory, which represents the most frequent movement patterns under specific operating conditions, such as air traffic flow, traffic mix, runway-in-use configuration, and meteorological conditions. It serves as a horizontal reference for arrival routes during the following 4D trajectory optimization.
2.3. Optimization of Continuous Descent 4D Trajectory with Point Merge Pattern
2.3.1. Continuous Descent Operation in Point Merge System
The PMS is a new type of arrival route structure based on satellite navigation, proposed by the European Organization for the Safety of Navigation (EUROCONTROL). It provides a systematic method for the sequencing and continuous descent of arriving aircraft, and has been gradually adopted worldwide in recent years. A typical point merging system consists of inner and outer sequencing arc legs and a merge point (MP), as shown in
Figure 3. An arc leg is a segment of arc with equal distance from the merge point, which is used to regulate the spacing of successive aircraft. For the purpose of establishing sufficient vertical separation, the inner arc is generally higher than the outer arc by several flight levels, and aircraft are required to maintain the same altitude and constant speed while on the arc. The merge point is an orientation point designed within the terminal area. Aircraft arriving from different directions fly into the arc leg from their arrival routes and turn directly to the merge point upon receiving turning instructions from air traffic controllers.
When an aircraft passes through the merge point, it will fly along a short segment to the final approach fix. The most important task for the pilot of the aircraft during this period is to smoothly and quickly intercept the Instrument Landing System (ILS) localizer. There is no need to optimize the trajectory at this segment and, so, the MP is used as the end of the CDO.
Figure 4 shows a typical CDO descent profile, combining the route structure of the point merge procedure and the basic requirements of the CDO. The diagram shows the relationship between flight distance to the MP, flight time to the MP, and flight altitude.
The descent can be broken down into five phases, according to standard airline operation procedures and the special structure of the point merge procedure. During phase 1, the aircraft maintains a constant Calibrated Airspeed (CAS) at the flight level while entering the terminal airspace, until it reaches the Top of Descent (TOD). The regulations for operating large commercial aircraft require that the CAS of an aircraft must not exceed 250 knots below 10,000 ft. Additionally, air traffic control mandates that the CAS of an aircraft must not be less than 250 knots above 10,000 ft. Therefore, during phase 2, the aircraft must descend to 10,000 ft and decelerate to 250 knots. In phase 3, the aircraft descends to the specified altitude of the sequencing arc and decelerates to the specified speed limit of the arc. During phase 4, the aircraft flies along the arc while maintaining a constant speed and altitude. In phase 5, the aircraft leaves the sequencing arc and gradually descends to the specified altitude and decelerates to restricted speed at the MP.
In this section, we employed spectral clustering to identify the traffic flow with the PMS pattern and to obtain the core route of this traffic flow. This core route represents the horizontal trajectory of the aircraft among this traffic flow, thereby providing a horizontal reference for 4D trajectory optimization of the arriving aircraft. Subsequently, we optimize the vertical profile of this core route based on the five phases illustrated in
Figure 4. As such, the horizontal route and its optimal vertical profile constitute the complete 4D trajectory.
2.3.2. Optimal Control Model for Point Merge Pattern
The optimization problem for the CDO vertical profile in point merge operations can be transformed into a multi-stage nonlinear optimal control problem, given the initial conditions of horizontal route, aircraft type, arrival altitude, and speed. The relevant constraints are determined according to air traffic rules depending on the given horizontal distance of the arrival route and the length of the sequencing arc. A multi-stage nonlinear optimal control model of the CDO profile is then established in order to meet the requirements of the point merge systems.
- (1)
Objective function
There is a simple proportional relationship between aircraft fuel consumption and CO
2 emissions [
43]. The CO
2 emission index of aviation kerosene is 3.15; that is, 3.15 kg of CO
2 is produced for every 1 kg of fuel burned. Therefore, the objective of minimum CO
2 emissions can be mapped to fuel consumption optimization. The objective function is given as Equation (8):
Meanwhile, air traffic control aims to expedite air traffic flow while minimizing flight time. The objective function is shown in Equation (9):
In these equations, the fuel consumption rate for phase p is represented by , while denotes the level flight true airspeed of the aircraft during the same phase. Additionally, represents the distance to be flown from the end of phase p to the MP, and indicates the end moment of the CDO.
Typically, in the actual operation of a commercial aircraft, fuel and time-related costs are considered together through setting a Cost Index (
CI) in the Flight Management System (FMS). Thus, an economic cruising speed that minimizes the total Direct Operating Cost (DOC) is obtained. The objective function is shown in Equation (10):
- (2)
Control Variables
In the optimal control model, the aircraft state
X is defined in terms of the distance to the MP
x, the altitude
h, and the true airspeed
. Therefore,
. The force analysis of the aircraft descent is used to calculate each variable of the flight state, as shown in Equation (11):
Equation (12) is used to calculate
T (engine thrust) and
D (drag), where
m is the mass of the aircraft and g denotes the gravitational acceleration:
where
Cdes represents the thrust correction coefficient,
CT1 represents the first thrust coefficient,
CT2 represents the second thrust coefficient,
CT3 represents the third thrust coefficient,
CD0 represents the zero-lift drag coefficient,
CD2 represents the induced drag coefficient,
ρ represents the atmospheric density, and
Sw represents the reference area of the wing. These coefficients are related to the aircraft type, and can all be derived from the aircraft parameters in the BADA documents [
44].
According to Equation (12), the thrust of an aircraft during descent is solely dependent on its altitude in an ideal scenario. As the descent angle γ determines the aircraft state, it may be considered as a unique control variable.
- (3)
Constraints
The values of state variables (e.g., flight distance and altitude) in each phase must be restricted to a corresponding range, as shown in Equation (13).
where
represents the initial altitude for phase
p,
represents the final altitude for phase
p,
represents the handover altitude when the aircraft enters the terminal airspace, and
represents the altitude when the aircraft is at its MP.
Equation (14) shows the corresponding speed restrictions that the aircraft must satisfy in different phases when descending within the terminal airspace:
where
represents the calibrated airspeed of the aircraft in phase
p,
is the maximum allowable CAS of the aircraft,
is the CAS at the end of phase
p, and
is the CAS at the MP.
To avoid an uncomfortable situation for passengers caused by a high descent rate, the glide path angle
γ must be limited to a certain range, as shown in Equation (15):
where
and
are the upper and lower limits of the glide path angle, respectively.
Although the optimal control model divides the complete flight into several independent phases, it is important to maintain consistency between the states before and after the critical point of the artificially divided successive segments, as shown in Equation (16):
where
and
denote the initial and final states of the aircraft in phase
p, respectively.
2.3.3. GACDO Algorithm for Point Merge Pattern
The challenge in solving the optimal control model for continuous descent operations lies in the uncertainty regarding the future flight time. To address this issue, we use the exemplar trajectory identified in
Section 2.2 as the reference route, and choose the horizontal flight distance as one of the key decision variables for modelling.
Most airlines have developed operating procedures for various aircraft types, in order to enhance their economic efficiency, safety, and passenger comfort. These procedures offer clear guidance on the CAS profile of the aircraft, as illustrated in
Figure 5. This study models the CAS profile of the arrival process using seven decision variables, denoted by
.
Figure 5 shows that
represents the extreme CAS for the second phase, while
represents the CAS when the aircraft is flying on the sequencing arc. The notation used in this context is as follows:
represents the distance of the level flight segment before the top of descent (TOD);
represents the distance flown by the aircraft as it changes from its initial speed to
;
represents the distance flown by the aircraft as it maintains
;
represents the distance flown by the aircraft as it decreases from
to 250 kn; and
represents the distance flown by the aircraft as it descends from an altitude of 10,000 ft to the sequencing arc altitude.
Thus, the optimization model for the vertical profile is established as shown in Equation (17):
The objective function is represented by , while x(3) represents the flight distance from the beginning of the sequencing arc to the MP. The minimum value of is constrained to 2 NM by x(3) − x(0) − 2. To facilitate aircraft operations, decision variables are typically constrained to integer values.
During level flight, the aircraft maintains a constant CAS, resulting in linearly varying fuel consumption with flight distance. However, during the descent phase, the fuel consumption is nonlinearly related to the flight distance due to the complex model dynamics. To solve this nonlinear objective function during aircraft descent, a numerical solution method is proposed.
The descent between the interval [
x(1),
x(3)] is used as an example to illustrate the numerical solution algorithm for fuel consumption. There are five known data points—namely, (
x(1),
(
x(a),
(
x(b),
(
x(2),250
(
x(3),
and (
x(4),
—from the aircraft velocity profile shown in
Figure 5. Similarly, dividing the interval [
x(1),
x(3)] by a sufficiently small distance
will produce a number of query points
xi,
i = [1, 2,…,
N1], where
N1 denotes the number of query points. Assuming that the CAS is uniformly decreasing, a linear interpolation method can be used to generate the calibrated airspeed corresponding to any query point
xi.
The conversion of calibrated airspeed
to true airspeed
under international standard atmospheric conditions is given by Equation (18):
where
and
hi represent the true airspeed and altitude, respectively, at the query point
.
For a sufficiently short time period
, the acceleration
can be approximated as a constant, and the displacement
is linearly related to the square of time, with the underlying kinematic shown in Equation (19):
From Equations (11) and (19), the displacement with time can be approximately described by a quadratic equation, as shown in Equation (20):
Equation (21) demonstrates how to derive the flight time differential using the horizontal flight distance increment
and the differential
of the true airspeed from Equation (20):
From Equations (11) and (20), the corresponding flight time
, glide path angle
, altitude
, and fuel flow
are calculated iteratively, as shown in Equation (22):
where
Cf1,
Cf2,
Cf3, and
Cf4 are the correction factors, which can all be obtained from the data in the BADA documents [
44].
According to the velocity–distance profile shown in
Figure 5, the state of the aircraft at each query point of the continuous descent operation can be obtained and, finally, the total fuel consumption
and the total flight time
can be calculated.
The fitness function incorporates penalty terms during evaluation, in order to ensure that the solution satisfies the constraints of the optimal control problem. The glide path angle constraint function
is defined as shown in Equation (23), in order to determine whether the glide path angle is within the range:
In the velocity profile shown in
Figure 5, the altitude
that the aircraft must reach at
(
p = 2, 3, 5) follows Equation (13). In order to ascertain whether the iteratively calculated altitude meets its constraints, the altitude constraint function
is defined as shown in Equation (24):
where
δ is an allowable height tolerance. Finally, the fitness function of the GACDO model is given by Equation (25):
where
M is a sufficiently large positive integer.
In accordance with the aforementioned models, a numerical solution for a continuous descent 4D trajectory-based genetic algorithm was developed in Matlab 2016b, which is referred to as the GACDO algorithm.