Evolution mechanism and optimisation of traffic congestion

Abstract

Air route networks can no longer meet operational efficiency requirements because of the rapid growth of complex traffic flows. Machine learning is employed to investigate the evolutionary mechanism of congestion in such networks in view of their high complexity and high density, and a reasonable network optimisation scheme is presented. First, deviations between nominal and actual routes are investigated with reference to radar track data, and a network reflecting actual route operations is constructed using adversarial neural networks. Second, flight time is used to characterise congestion in route networks. Actual network operations are considered, and congestion is defined from the perspective of road traffic engineering. The effects of the operational properties of traffic flows on flight times are analysed to establish various congestion indicators. A gradient boosting model is used to select indicator characteristics and analyse patterns in the variations of indicator values for each flight segment in distinct periods. The indicator–time relationship is leveraged to explore the evolutionary mechanism of congestion in the route network. Third, on the basis of this mechanism, a multiobjective optimisation model of congestion is formulated, and a particle swarm optimisation algorithm is executed to adjust the route passage structure, thereby solving the optimisation model. Finally, calculation validation is conducted using radar track data from the control sector of the Yunnan region. The average flight time in a route segment is 10% shorter in the optimised route network than in the nonoptimised route network, which confirms that the optimisation solution is practicable.

Nomenclature

DF

Dickey and Fuller

ADF

Augmented Dickey-Fuller

ReLU

Rectified Linear Unit

GAN

Generative Adversarial Networks

GBRT

Gradient Boost Regression Tree

1.0 Introduction

To maximise the use of airspace resources, machine learning, deep learning and big data analysis are employed to analyse the factors influencing route network congestion, study the evolutionary mechanism of such congestion, and construct a route network optimisation scheme with reference to actual traffic flow and network operation characteristics.

The literature on the mechanism of congestion propagation in air route networks, which concerns the field of civil aviation, is incipient. By contrast, in-depth studies on the evolutionary mechanism of congestion in road traffic are numerous. Dong et al. [1] cellularised traffic flows and roads and used the obtained cells to simulate traffic flow changes and demonstrate the evolutionary mechanism of traffic congestion by constructing a transportation model. Polson and Sokolov [2] estimated the traffic flow density state and Lighthill–Whitman–Richards parameters by using a particle filter learning algorithm. They also explored the mechanism of traffic state change by investigating the relationship between traffic flow and density. Omer [3] developed a mixed-integer linear model based on spatial discretisation to ensure the true trajectory by considering the speed vector relative to time continuity and the restrictions of speed, acceleration and yaw rate; and confirmed that the model can solve complex situations within a few seconds without incurring more than a few kilograms of extra fuel consumption per aircraft. Li et al. [4] probed the effects of passenger flow changes on the traffic state and introduced the passenger flow load index as a congestion parameter for studying the dynamic evolution of traffic congestion. Jiang [5] examined the factors affecting the operational state of traffic flow in a terminal area; established indicators for evaluating the congestion state; and defined an evolutionary, spatiotemporal mechanism of congestion. On the basis of vehicle trajectory data recorded by the global positioning system, Yang [6] established an urban road network model to identify the number of traffic congestion zones in a city, analyse the spatiotemporal evolutionary characteristics of congestion at the micro and macro levels, and determine the spatiotemporal evolutionary pattern of traffic congestion. Chang [7] examined the causes of traffic congestion; studied the influences of road grades, the frequency of real-time information updates, and other factors on the diffusion of congestion; and established a mechanism of traffic congestion on the basis of changes in traffic information. Yu et al. [8] developed a congestion identification model based on the discriminatory criteria of free and congested states to analyse patterns in traffic state changes on the basis of traffic density trends and traffic speed upstream and downstream of a highway. Çeçen et al. [9] presented a mixed integer linear programming model (MILP), which can solve aircraft conflicts in 450 different traffic scenarios in less than a minute. Mohammad et al. [10] used vehicle speed as a basis for assessing traffic congestion and traffic density, determined the traffic congestion index at peak hours, and explored the evolutionary mechanism of congestion. Cecen et al. [11] proposed a two-step optimisation method for solving aircraft conflicts within the tactical pre time window in general free route airspace which can resolve all conflicts in less than 4min. Basturk et al. [12] proposed prediction of aircraft estimated time of arrival (ETA) using machine learning algorithms which can accurate prediction of ETA was important for management of delay and air traffic flow etc. Wu et al. [13] proposed a depth generative model based on one-dimensional convolutional neural network (Conv1D-GAN), two-dimensional convolutional neural network and short-term memory neural network (LSTM-GAN), and demonstrated that Conv1D GAN was the most suitable generative adversarial approach for long-term aircraft trajectory prediction.

In the present study, radar track data are first preprocessed. Next, correlations of track features are analysed, and an operational route network is generated using adversarial neural networks to extract nominal tracks from the radar track data. Subsequently, the definition of congestion in road traffic engineering is extended to the field of civil aviation, and flight time is employed as the basis for assessing congestion in air networks, selecting characteristics to obtain a set of critical congestion indicators, analysng the distribution relationship between these indicators and time and exploring the evolutionary mechanism of congestion as a function of time. Finally, the traffic structure of the route network is optimised on the basis of the evolutionary mechanism, thereby reducing congestion.

2.0 Generating route networks from radar track data

2.1 Radar track data preprocessing

The radar track data, which comprise discrete data points with a 4-s refresh rate, contain information such as time, aircraft position coordinates, altitude, heading, ground speed and flight number. These high-volume, typically high-dimensional data must be preprocessed prior to analysis. Raw radar track data contain errors such as missed points and data jumps; thus, data cleaning and data extraction must be performed.

2.1.1 Extraction of effective trajectories

The limitations of data collection equipment and the unpredictability of weather conditions make it difficult to guarantee the validity of each piece of data. In this study, the radar track data are first filtered by height and selected for tracks above the altitude of the terminal area (6,000m). Next, the data are preprocessed to fill in missing values and remove outliers.

The validity check of radar track data is performed on the basis of the following principles. First, a single track ${P_i}$ contains more than 500 track points, and the time interval between any two adjacent points does not exceed 40s. Second, all incorrect trajectories are screened out; no unreasonable jump in height is observed between any two adjacent points in a single trajectory ${P_i}$ .

2.1.2. Flight path resampling

The collection interval of the radar track data is short. Aircraft have a long cruise phase in actual flight, during which their characteristic information remains essentially unchanged (except the coordinate positions), which results in redundant trajectory data. A sequence of trajectories of equal length is required for generating an operational route network from actual flight trajectories, and the number of radar data points contained in different trajectories varies according to the extraction of valid trajectories from the radar data. Following this extraction process, the radar track data are processed to eliminate redundancy and ensure that data points are uniform in length while retaining flight characteristics.

Padding with zeros is performed to equalise the length of the data points, thus preserving the original trajectory characteristics and ensuring data accuracy.

2.2 Generating real-world operational route networks by using adversarial neural networks

2.2.1 Neural network construction methods

Since their introduction by Goodfellow et al. in 2014, Generative Adversarial Networks (GAN) have been extensively researched and applied in the field of deep generative models [14].

When constructing an adversarial neural network, certain methods must be implemented to improve network performance, and relevant basic theories must be analysed.

1. (a) Excitation function

In a multilayer neural network, between the output value of the network node in the previous layer and the input value of the network node in the next layer, a nonlinear excitation function is essential for enhancing the representativeness of the network.

We employ the rectified linear unit (ReLU) function and a sigmoid function. The ReLU function can solve the gradient dispersion problem. In stochastic gradient descent, the calculation speed and convergence speed can be rapidly improved by the derivative of the ReLU function (1 or 0). The sigmoid function, which converts all values into a value between 0 and 1, outputs 0 for extremely small negative numbers and outputs 1 for extremely large positive numbers.

1. (b) Dropout

Dropout, which was proposed by Hinton, is a technique for preventing the overfitting of a complex feedforward neural network (which occurs frequently when such a network is trained on a small dataset) by preventing feature detectors from acting together. This technique enhances the performance of a neural network [15]. Dropout reduces interactions between neural nodes in the hidden layer, which makes the neural network more generalisable. The standard process of dropout in neural networks involves forward-propagating the data input x and then backpropagating the errors to determine how to update the network learning parameters. Subsequently, the neural network runs as follows.

Step 1: Half of the neurons in the hidden layer of the neural network are randomly deleted. The neurons in the input and output layers are untouched.

Step 2: The data input x is propagated forward, and the outcome of the resulting loss function is propagated backward. When the propagation of the selected batch of training data is complete, the corresponding parameters $(\omega ,b)$ are updated on the remaining half of the neurons through stochastic gradient descent.

Step 3: The deleted neurons are restored to their initial state, and the other half of the neurons are updated.

Step 4: After 50% of the neurons from the hidden layer of the neural network are randomly selected for backup, they are removed.

The aforementioned steps are repeated until training is complete.

1. (c) Loss function and cross-entropy function

The loss function is used to calculate the error between the predicted and true values of the network. The smaller the value of the loss function, the more favourable is the performance of the network. A commonly used loss function is the cross-entropy function, which is expressed as follows:

(1) \begin{align}L = \frac{1}{N}\sum\limits_i {{L_i} = \frac{1}{N}} \sum\limits_i { - \left[{y_i}\log {p_i} + (1 - {y_i})\log (1 - {p_i})\right ]} \end{align}

where L is the cross-entropy value; ${y_i}$ denotes the label of sample i; positive and negative classes are designated as 1 and 0, respectively; and ${p_i}$ is the probability that sample i is predicted to be positive.

2.1.2 Extraction of nominal trajectories

On the basis of the theory of adversarial neural network construction, the ReLU function is employed as the excitation function in the intermediate layer, the sigmoid function is used as the excitation function in the discriminator output layer, the two dropout probabilities are set as 0.3, and the cross-entropy function is selected as the loss function. The trajectory data of a particular flight in a given year are selected as the sample data, and the adversarial neural network is developed to generate the nominal trajectory of the flight. The training data sample comprises the radar trajectory data of flights for 1 year; the padding method is used to obtain track data of equal length, and the number of training sessions is set as 1,000.

An adversarial neural network is constructed and trained on the 1-year radar track data for a flight. The loss function is calculated, and the iterative plot of the loss function values is displayed in Fig. 1. The loss function values drop rapidly until approximately 50 training sessions have elapsed, after which they level off and stabilise at <0.1. Therefore, the constructed neural network can appropriately generate the nominal track on the basis of the radar track data and accurately reflect the actual operating route network.

Figure 1. Iterative graph of the loss function.

3.0 Evolutionary mechanism of congestion in air route networks

3.1 Theory of congestion in air route networks

Air route network congestion has no established definition in the civil aviation sector. Traffic congestion is well defined in road traffic engineering. Shibata et al. [16] suggested that the traffic status of different road sections depends on the length of the arrival time. Zhang [17] defined congestion as a function of passing time. Specifically, congestion exists when the passing time exceeds the standard time corresponding to normal operations. Yan [18] stated that traffic congestion occurs when vehicles queues are lengthened. Du [19] defined traffic congestion according to the average travel time, taking the ratio of the average travel time to the driving time in a noncongested state as an index of traffic congestion. A larger ratio (i.e., a longer delay) indicates a greater degree of congestion. Overall, most definitions of congestion in road traffic engineering are based on trip time and trip speed.

According to the definition of traffic congestion in road traffic engineering and the actual situation of aircraft operation on airways, we consider the description of the congestion state of the air network from the perspective of flight time and analyse its evolutionary mechanism. If the flight time of a route segment exceeds the average flight time of the segment during a certain period, congestion exists in the segment.

3.2 Selection and analysis of congestion indicators

3.2.1 Selection of segment congestion indicators

On the basis of the definition of congestion in road traffic engineering and the actual operation of the air network, the air network congestion index is selected to characterise the state of congestion in the air network with regard to the flight time and its influencing factors. These factors are defined as follows.

1. (a) Number of flights in the segment

Number of flights in the segment refers to the total number of aircraft ${n_{{r_j}}}$ passing through segment ${r_j}$ in the sampling interval.

1. (b) Segment traffic flow

Segment traffic flow refers to the number of flights that pass through the route per unit of time and is reflective of the operational load on the segment. This parameter is expressed as follows:

(2) \begin{align}{Q_{{r_j}}} = \frac{{{n_{{r_j}}}}}{t}\end{align}

where ${Q_{Li}}$ is the flow of traffic in segment ${r_j}$ , ${n_{{r_j}}}$ is the number of aircraft in segment ${r_j}$ , and t is the statistical time interval.

1. (c) Number of altitude adjustments

Number of altitude adjustments refers to the number of altitude changes made by the aircraft in flight during the statistical time interval. If the first-order difference corresponding to aircraft altitude over a statistical time interval remains at a certain altitude level for 1 min, the altitude adjustment of the aircraft is complete. The relevant equations are presented as follows:

(3) \begin{align}g(\Delta h) = \left\{ {\begin{array}{*{20}{c}}{0,\Delta {h_{{t_i}}} - \Delta {h_{{t_j}}} \le 600}\\{1,\Delta {h_{{t_i}}} - \Delta {h_{{t_j}}} \gt 600}\end{array}} \right.\end{align}

(4) \begin{align}{N_C} = \sum\limits_t {g(\Delta h)} \end{align}

where $\Delta h$ is the first-order difference corresponding to aircraft altitude, ${t_i}$ is the start time of the statistical time interval, ${t_j}$ is the end time of the statistical time interval, t is the statistical time interval, and ${N_C}$ is the number of altitude adjustments the flight segment. The standard altitude layer interval for cruising is 600 m.

1. (d) Rate of approaching traffic in a segment

The proximity of all aircraft in a segment of the route network during a statistical time interval reflects the distribution density of aircraft within that segment [20]. China’s civil aviation regulations state that in a radar-controlled environment, aircraft must meet a minimum horizontal separation of >10km or a minimum vertical separation of 300 m. Because the Earth is ellipsoidal, the Euclidean distance cannot simply be used to calculate the proximity of one aircraft from another. Therefore, the relative ellipsoidal distance between two aircraft is calculated, and the rate of approaching traffic in a segment is determined. The relevant equations are presented as follows:

(5) \begin{align}{d_{ij}} = {\left\| {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over p} }_i} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over p} }_j}} \right\|_{s{i_h},s{i_v}}} = \sqrt {\frac{{{{({x_i} - {x_j})}^2} + {{({y_i} - {y_j})}^2}}}{{s{i_h}^2}} + \frac{{{{({h_i} - {h_j})}^2}}}{{s{i_v}^2}}} \end{align}

(6) \begin{align}p{r_L} = \sum\limits_{i = 1}^n {\sum\limits_{j = 1,i \ne j}^n {\exp \left( - \frac{{{d_{ij}}}}{n}\right )} } \end{align}

where ${d_{ij}}$ is the relative distance in space between aircraft i and aircraft j; ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over p} _i}$ and ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over p} _j}$ are the position vectors of aircraft i and aircraft j, respectively; $s{i_h}$ and $s{i_v}$ are the horizontal minimum safe interval and vertical minimum safe interval of the segment, respectively ( $s{i_h} = 10{\rm{km}}$ . $s{i_v} = 300{\rm{m}}$ ); $({x_i},{x_j})$ and $({y_i},{y_j})$ and $({h_i},{h_j})$ are the coordinate positions and altitude positions of aircraft i and aircraft j; $p{r_L}$ is the rate of approaching traffic in the segment; and n is the number of aircraft in the segment.

1. (e) Number of altitude layers passed by the aircraft

The number of altitude layers passed by all aircraft in leg L during the statistical time interval is expressed as follows:

(7) \begin{align}N_{fl}^L = \sum\limits_m {\frac{{h_{{t_i}}^m - h_{{t_j}}^m}}{{{h_{std}}}}} \end{align}

where $N_{fl}^L$ is the number of altitude layers above and below the aircraft in segment r, $h_{{t_i}}^m$ is the altitude of aircraft m at the beginning of the statistical time interval, $h_{{t_j}}^m$ is the altitude of aircraft m at the end of the statistical time interval, and ${h_{std}}$ is the specified altitude layer interval for cruising (i.e., 600 m).

1. (f) Segment traffic altitude adjustment rate

Segment traffic altitude adjustment rate refers to the ratio between the altitude adjustment and the initial altitude of all aircraft in segment r during the statistical time interval. This parameter reflects the magnitude of the range of altitude adjustments made by aircraft in a segment and is expressed as follows:

(8) \begin{align}N_{hr}^r = \frac{{\left| {\sum\limits_m {\frac{{h_{{t_i}}^m - h_{{t_j}}^m}}{{h_{{t_i}}^m}}} } \right|}}{n}\end{align}

where $N_{hr}^r$ is the segment traffic height adjustment ratio, $h_{{t_i}}^m$ is the altitude of aircraft m at the beginning of the statistical time interval, $h_{{t_j}}^m$ is the altitude of aircraft m at the end of the statistical time interval, and n is the total number of aircraft operating in segment r during the statistical time interval.

1. (g) Segment traffic density

Segment traffic density refers to the instantaneous aircraft count per unit length of a segment r in a route network. This parameter reflects the intensity of flight in a segment and is expressed as follows:

(9) \begin{align}{k_{{r_j}}} = \frac{{{n_{{r_j}}}}}{{{l_{{r_j}}}}}\end{align}

where ${k_{{L_i}}}$ is the density of traffic in segment ${r_j}$ , ${n_{{r_j}}}$ is the number of aircraft operating in segment ${r_j}$ at a given instant during the sampling interval, and ${l_{{r_j}}}$ is the length of segment ${r_j}$ .

3.2.2 Correlation analysis of congestion indicators

The kernel density estimation of nonparametric estimation is used for correlation analysis of congestion indicators [21]. Let ${x_1},{x_2}, \ldots ,{x_n}$ denote the n sample points of an independent identical distribution F with the probability density function f. Density estimation is conducted using a two-dimensional Gaussian kernel, where $n = 7$ . The expression $N = \{ {x_1},{x_2},{x_3},{x_4},{x_5},{x_6},{x_7}\} $ , represents the seven congestion indicators.

(10) \begin{align} {\hat{f(x)}} = \frac{{{\rm{ }}\sum\limits_{i = 1}^n {\exp \left ( - \frac{{{{({x^{\left( 1 \right)}} - x_i^{\left( 1 \right)})}^2}}}{{2{\eta _1}^2}}\right )} \exp \left ( - \frac{{{{({x^{\left( 2 \right)}} - x_j^{\left( 2 \right)})}^2}}}{{2{\eta _2}^2}}\right )}}{{\sqrt {2\pi } n{}_1{n_2}{\eta _1}{\eta _2}}} \end{align}

Where ${x^{\left( 1 \right)}}$ is index 1; ${x^{\left( 2 \right)}}$ is index 2; and $\eta $ is the bandwidth, where $\eta \gt 0$ .

By estimating the Gaussian kernel density between the indicators and analysing the distribution between the indicators, congestion in the route network can be scientifically examined. In total, 26 routes (comprising 55 flight segments) are selected, and the time-series data of 10-day traffic flow are averaged and sorted into one-day-average data (in 1-h intervals) for each flight segment. The number of aircraft, the number of upper and lower altitude levels, and the Gaussian kernel density estimation results for the two indicators and flight time are displayed in Figs 2, 3 and 4, respectively. Furthermore, linear regression is performed to determine whether the congestion indicators and flight time are linearly related (Figs 5 and 6).

Figure 2. Joint density distribution of the number of flights and the flight time.

Figure 3. Joint density distribution of the altitude adjustments and the flight time in a flight segment.

Figure 4. Joint density distribution of the altitude adjustments and the number of flights in a flight segment.

Figure 5. Linear regression between the number of flights and the flight time in a flight segment.

Figure 6. Linear regression between the altitude adjustments and the flight time.

Figure 2 presents the joint density distribution of the number of flights in a flight segment and the flight time. The concentrated dark areas correspond to 1–10 flights and flight times of 0 to 0.1 h. The distribution in Fig. 2 is convex, which indicates that a strong correlation exists between the number of flights and the flight time.

Figure 3 shows the joint density distribution of the number of altitude adjustments and flight time for a flight segment. The concentrated dark areas correspond to 1–50 altitude layers and flight times of 0 to 0.1h. The distribution in Fig. 3 is convex, which indicates that the joint density is normal distribution.

Figure 4 displays the joint density distribution of the number of altitude adjustments and the number of flights. The concentrated dark areas correspond to 1–50 altitude adjustments and 1–10 flights. The joint density distribution displayed in Fig. 4 is convex, which indicates that this distribution is regular. In practise, the higher the number of in-segment flights, the higher is the number of altitude adjustments within the segment. This phenomenon is observed possibly because the number of height adjustments made by individual aircraft remains constant but the number of aircraft increases, which increases the overall number of altitude adjustments. An alternative explanation is that to ensure flight safety when the number of aircraft within a segment increases, the number of aircraft making altitude adjustments also increases, which increases the overall number of altitude adjustments within the segment.

Figure 5 presents a linear regression plot of the number of flights against the flight time in a segment. These two variables have a linear relationship. In practise, the higher the number of flights in a segment, the higher is the congestion within the segment and thus the longer is the flight time.

Figure 6 displays a linear regression plot of the number of altitude adjustments against the flight time. These two variables are also linearly related. In practise, the more times an aircraft adjusts its altitude within a segment, the higher is the likelihood of congestion within the segment and thus the longer is the flight time.

3.3 Feature selection of congestion indicators

Feature selection is an effective dimension reduction method to solve the problem of high-dimensional data and over fitting in the process of machine learning. In this paper, the random forest algorithm and the gradient boosting tree algorithm are used to select the route congestion indicator characteristics respectively. And show the advantages and disadvantages of the two algorithms.

3.3.1 Route congestion indicator feature selection

The application of gradient boosted regression trees to feature importance assessment involves evaluating the importance of feature ${x_i}$ on each tree and then averaging its importance on all trees.

(11) \begin{align}\left\{ {\begin{array}{*{20}{c}}{I_{{x_i}}^2 = \dfrac{1}{M}\sum\limits_{m = 1}^M {I_{{x_i}}^2({T_m})} }\\{I_{{x_i}}^2({T_m}) = \sum\limits_{j = 1}^{J - 1} {{d_j}} }\end{array}} \right.\end{align}

where $I_{{x_i}}^2$ is the squared importance of feature ${x_i}$ , M is the number of decision trees in the gradient boosted tree, ${T_m}$ is a decision tree, j is the branch node, m is the number of iterations, and ${d_j}$ denotes the extent to which the independent variable ${x_i}$ is boosted as the squared error of the jth branch node. The relative importance of all independent variables adds up to 1.

Random forest is a classification integration algorithm, which generates multiple decision trees in the process of feature selection. By training the sample data, the importance of each feature is calculated and sorted, and the features with greater contribution are selected [22].

3.3.2 Assessment of the importance of congestion indicator characteristics

According to the above theory, gradient boosting trees and random forest algorithms are used for the selected congestion indicator characteristics, respectively. Because of the different dimensions of the indicator values, standard deviation standardisation (z-score) must be performed for the indicator values. Next, 180,000 randomly sampled items from the congestion indicator dataset are used to calculate the feature importance.

As presented in Tables 1 and 2, regardless of whether a route is congested or noncongested, the number of altitude adjustments, the number of flights, and the rate of approaching traffic account for 95% of the importance of all the indicators. Thus, these three indicators are integral to the examination of congestion. A regression model can be constructed to analyse the relationship between these three indicators and flight time, thereby facilitating the examination of the evolutionary mechanism of congestion.

Table 1. Importance of congestion indicators in congested conditions

Table 2. Importance of congestion indicators in noncongested conditions

The random forest algorithm is employed to assess the significance of traffic congestion indicator attributes for each segment during congested routes. Furthermore, a random forest model with a decision number of 7 is constructed to determine the importance of traffic congestion indicators for each segment under congested and smooth conditions within the route network. The contribution of each feature is presented in Tables 1 and 2.

To ensure the stability of the two models, a 10-fold cross-validation is conducted. Cross-validation is a vital technique in statistical analysis to assess the generalisability of models on independent data sets [23]. In this study, the gradient boosting regression tree model and the random forest model are separately cross-validated in ten folds for congested and noncongested routes. Figure 7 displays the outcomes of the ten-fold cross-validation.

Figure 7. Results of tenfold cross-validation for RF and GBDT.

The consistent feature selection results depicted above indicate that the gradient boosted regression tree model outperforms the random forest algorithm. Therefore, the gradient boosting regression tree algorithm is ultimately selected for feature screening.

3.4 Analysis of the evolutionary mechanism of congestion

Due to the complex relationship between important congestion indicators and flight time, in order to ensure the accuracy of the model, gradient lifting regression tree model is used for analysis to explore the distribution relationship between important congestion indicators and flight time.

3.4.1 Gradient boosted regression tree

The negative gradient value of the loss function indicates that thegradient boosted regression tree (GBRT) model fits the residual approximation of the regression problem boost tree. The optimal gradient descent step ${\rho _m}$ is added to update the model and thereby obtain the final regression tree model.

To initialise ${f_0}(x)$ , the loss function L is employed as follows:

(12) \begin{align}{f_0}(x) = \arg \mathop {\min }\limits_\rho \sum\limits_{i = 1}^N {L({y_i},\rho )} \end{align}

The aforementioned function is a squared error loss function that is commonly used in gradient boosted regression models.

(13) \begin{align}L(y,f(x)) = {(y - f(x))^2}\end{align}

The negative gradient direction ${{\rm{z}}_{im}}$ is calculated, where $m = 1,2, \ldots ,M$ and $i = 1,2, \ldots ,N$ .

(14) \begin{align}{z_{im}} = - {\left [\frac{{\partial L(y,f({x_i}))}}{{\partial f({x_i})}}\right ]_{f(x) = {f_{m - 1}}(x)}}\end{align}

By using the training set data, the regression tree model is trained on the basis of the negative gradient ${z_{im}}$ , and the size of the optimal step ${\rho _m}$ for gradient descent is calculated as follows:

(15) \begin{align}{f_m}(x) = {f_{m - 1}}(x) + \sum\limits_{j = 1}^J {{\rho _m}{b_{jm}}I\{ x \in {R_{jm}}\} } \end{align}

(16) \begin{align}{\rho _m} = \arg \mathop {\min }\limits_\rho \sum\limits_{i = 1}^N L\left ({y_i},{f_{m - 1}}({x_i})\sum\limits_{j = 1}^J {{\rho _m}{b_{jm}}I\{ x \in {R_{jm}}\} } \right )\end{align}

where $I\{ x \in {R_{jm}}\} = \left\{ {\begin{array}{l@{\quad}l}1, & {\rm{ if }}\ x \in {R_{jm}}\\0, & {\rm{otherwise}}\end{array}} \right.$ , ${R_{jm}}$ is the region into which the feature space is divided, and ${b_{jm}}$ is the variable branch point.

The gradient boosted regression tree model can be used to analyse and visualise how a dependent variable varies with an independent variable (i.e., the nonlinear effect of the independent variable on the dependent variable).

The test of the gradient boosted regression tree model must be based on the predictive effect of the model to ensure that the overall prediction error is maintained within a certain range. The common explained variance, mean absolute error, mean squared error, and R ² are selected as evaluation metrics for testing the goodness-of-fit of the model.

In the cross-validation of the present gradient boosted regression tree model, the example of a congested air network is considered. The obtained scores are presented in Fig. 8.

Figure 8. Results of tenfold cross-validation.

The first and second scores are low. The remaining scores are clustered around 0.9, which is high. The average of the ten scores is 0.89, which indicates that the constructed model is stable, has favourable generalisability, and can be well applied to new datasets.

3.4.2 Analysis of the evolutionary mechanism of congestion

A gradient boosted regression model is employed to analyse the distribution relationship between a set of congestion indicators and flight time as well as the changes in congestion with variations in flight time.

Figure 9 displays the predicted flight time and actual flight time as functions of the number of altitude adjustments, the number of flights and the rate of approaching traffic. The orange lines represent the flight time obtained from the gradient boosted regression tree model, and the blue lines denote the actual flight time. The orange lines cover the blue lines well, which indicates that the gradient boosted regression tree model accurately fits the flight time to the three congestion indicators.

Figure 9. Fitting results of the gradient boosted regression tree model.

The values of the evaluation indicators in the gradient boosted regression tree model are presented in Table 3.

Table 3. Evaluation indicator values for the gradient boosted regression tree model

4.0 Multiobjective optimisation of route network congestion

4.1 Network congestion optimisation model

According to the evolutionary mechanism of congestion, the structure and distribution of traffic flow in the route network are optimised and adjusted from a systematic perspective such that congestion can be mitigated. The optimisation model can be more widely applied through the global optimisation and control of this congestion.

4.1.1 Model assumptions

Due to the complexity of the route network structure and actual operation, in order to restore the display problem more accurately, the relevant assumptions of the established mathematical model are as follows. First, the aircraft is considered as a mass point along the centerline of the route at a fixed cruising speed, regardless of individual differences such as aircraft type. Second, the connection relationship between waypoints (i.e. the connection order of all segments and nodes of any route) remains unchanged under optimisation. For example, the connection order of waypoints along route H before and after route network optimisation is DAL-P333-YJ-JHG. Third, three zones, namely the forbidden area, restricted area, and dangerous area, cannot be crossed. Fourth, congestion in the control area is not considered.

4.1.2 Model constraints

In view of the actual operation of the air network, the constraints described in the following sections must be established for developing a congestion optimisation scheme.

1. (a) Constraint on the upper limit of the route operation

The number of flights operating within a fixed period in the route should be less than the control limit. In general, as presented in (17).

(17) \begin{align}{x_{i,j}}(k) \lt M,k \in [1,K],i \in [1,n],j \in [1,m]\end{align}

1. (b) Nonnegative constraints on variables

The number of flights in any given route must be a positive integer.

(18) \begin{align}{x_{i,j}}(k) \in {N^*},k \in [1,K],i \in [1,n],j \in [1,m]\end{align}

4.1.3 Network optimisation modeling

To ensure the optimal operation of the route network at all times, the optimisation model takes 1 h as the optimisation period and analyses the dynamic operation of the traffic flow to minimise the operating and optimisation costs, the number of altitude adjustments and the rate of approaching traffic, the number of head-to-head flights diverted is maximised. Network optimisation model is expressed as follows.

1. (a) The objective of the established optimisation model is to minimise the operating and optimisation cost of the air network. The aforementioned objective function is expressed as follows:

(19) \begin{align}\min {\rm{ }}C = {C_P} + {C_d}{\rm{ = }}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {{\omega _{i,j}}} {C_{{p_1}}}} {\rm{ + }}\sum\limits_{i = 1}^n {\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^m {{x_{i,j}}(k){L_{i,j}}} } } \end{align}

(20) \begin{align}{\omega _j} = \left\{ {\begin{array}{*{20}{l}}{1, \quad \textrm{Add a new segment next to segment } r_j} \\ {0, \quad \textrm{otherwise}} \end{array}}\right.\end{align}

1. (b) The number of flights in a segment per unit of time should be minimised. The number of flights and the number of altitude adjustments in a segment are linearly related. The number of flights is set as x, and ${N_C}$ , which represents the number of altitude adjustments obtained through regression fitting, is expressed in (21). The fitting results reveal an R ² value of 0.866. Because the number of altitude adjustments is an integer, the calculation result obtained using (21) is rounded up.

(21) \begin{align}{N_C} = 0.3{x^2} + 5.417x - 1.702 \end{align}

The correlation analysis of the rate of approaching traffic and the number of flights in a segment reveal that these two indicators are not linearly related. The rate of approaching traffic, which is obtained from the relative distance between two aircraft, is affected by changes in the number of flights and the number of altitude adjustments. The relationship between these two indicators and the rate of approaching traffic is determined by considering the number of flights and the number of altitude adjustments concurrently. Through regression fitting, the relationships amongst the rate of approaching traffic, the number of flights and the number of altitude adjustments are established.

(22) \begin{align}pr = (0.116 + 0.79x - 9.12y + 0.639xy - 5.83{x^2}) \times {10^{ - 2}}\end{align}

where is equal to ${N_C}$ , that is, the number of altitude adjustments made in the segment.

Thus, the objective function of the optimisation model can be expressed as follows:

(23) \begin{align}\sum\limits_{i = 1}^n {\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^m {({N_C}} (k){)_{i,j}}} } = \frac{1}{m}\sum\limits_{i = 1}^n {\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^m {0.3x_{i,j}^2(k)} } } + 5.417{x_{i,j}}(k) - 1.702\end{align}

(24) \begin{align}\min \frac{1}{m}\sum\limits_{i = 1}^n {{\rm{ }}\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^m {p{r_{i,j}}(k)} } } {\rm{ }}\end{align}

(25) \begin{align}\min \frac{1}{m}\left (\sum\limits_{i = 1}^n {{\rm{ }}\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^m {(p{r_{i,j}}(k)} } } {\rm{ }} + {({N_C}(k))_{i,j}})\right )\end{align}

Given that the independent variables in the objective functions of running cost, altitude adjustments and proximity rate are frame times, the objectives functions expressed in (23) and (25) can be combined into one objective function (objective function 1). Objective function 1 is expressed as follows:

(26) \begin{align}\min \frac{1}{m}\left (\sum\limits_{i = 1}^n {{\rm{ }}\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^m {(p{r_{i,j}}(k)} } } {\rm{ }} + {({N_C}(k))_{i,j}})\right ) + C\end{align}

1. (c) Moreover, the number of head-to-head flights diverted is maximised, that is, the ratio of the number of head-to-head flights diverted to the total number of head-to-head flights diverted $u$ should be as large as possible. The objective function 2 is expressed as follows:

(27) \begin{align}\max {\rm{ }}\alpha {\rm{ = }}\frac{v}{u}\end{align}

where C is the total cost; ${C_{{P_1}}}$ is the cost of a new segment; and ${C_{_d}}$ is the operating cost, which includes the operating cost of the original segment and the cost of adding a detour to the new segment. The parameters ${L_{i,j}}$ and ${x_{i,j}}(k)$ represent the length of segment $j$ of route $i$ and the number of flights in this segment at moment k, respectively. The parameter ${x_{i,j}}(k){L_{i,j}}$ is used to calculate the operating cost of segment j of route $i$ at moment k. The parameter ${N_C}$ represents the number of altitude adjustments the segment, ${({N_C}(k))_{i,j}}$ denotes the number of altitude adjustments segment $j$ of route $i$ at moment k, pr is the rate of approaching traffic in the segment, and $p{r_{i,j}}(k)$ is the rate of approaching traffic in segment $j$ of route $i$ at moment k.

4.2 Multiobjective particle swarm optimisation–based model solution for air networks

We employ the particle swarm optimisation (PSO) algorithm with global optimisation characteristics to solve the multiobjective optimisation model of the route network. The particle swarm optimisation algorithm, which was first developed by Kennedy and Eberhart [24]. Mathematical descriptions of the PSO algorithm is as follow:

In particle swarm optimisation algorithms, a particle refers to an individual with unique properties, and an assemblage of numerous particles is called a population. In an dimensional particle swarm optimisation problem, the number of populations is, and each particle position in the iterative process can be regarded as the implicit optimum in the global optimisation process, where the current position of the first particle is $i\ {x_i} = {({x_{i1}},{x_{i2}}, \ldots ,{x_{in}})^T}$ , and its velocity is ${v_i} = {({v_{i1}},{v_{i2}}, \ldots ,{v_{in}})^T}$ . The personal best position of particle i in the iterative process is ${p_i} = {({p_{i1}},{p_{i2}}, \ldots ,{p_{in}})^T}$ , which is also denoted as ${P_{best}}$ , and the best position experienced by all particle groups is ${p_g} = {({p_{g1}},{p_{g2}}, \ldots ,{p_{gn}})^T}$ , which is also designated as the global best ( ${G_{best}}$ ). After the sum of ${P_{best}}$ and ${G_{best}}$ is calculated, the iterative process by which the velocity and position of the first particle move in the first j dimensions of the search space can be determined as follows:

(28) \begin{align}v_{il}^k = wv_{il}^k + {C_1}ran{d_1}\left (({P_{best}})_{il}^k - x_{il}^k\right ) + {C_2}ran{d_2}\left (({G_{best}})_l^k - x_{il}^k\right )\end{align}

(29) \begin{align}x_{il}^{k + 1} = x_{il}^k + v_{il}^{k + 1}\end{align}

Where $v_{il}^k$ and $x_{il}^{k + 1}$ are the velocity of particle i and the position components in dimension l, respectively. The parameter w represents the inertia weight, which is a critical parameter affecting the search performance of the algorithm. The parameters ${C_1}$ and ${C_2}$ represent two acceleration constants. Specifically, ${C_1}$ is the self-learning factor of the particle, whereas ${C_2}$ is the global learning factor of the particle. Two stochastic constants $ran{d_1}$ and $ran{d_2}$ are located in the interval (0, 1), $({P_{best}})_{il}^k$ is the individual extremum and the first-dimensional component of particle $i$ in the kth iteration, and $({G_{best}})_l^k$ is the first-dimensional component of the global extremum of the population in the kth iteration.

In multiobjective optimisation, objectives are incompatible; the optimal solution under one objective may be the least favourable solution under another objective. A solution that reduces at least one other objective function while improving one objective is called a Pareto optimal solution, the set of optimal solutions of an objective function is called the Pareto set, and the boundary defined by the set of all point mapped from the Pareto optimal set is called the Pareto-optimal front. A multiobjective particle swarm optimisation algorithm uses as few computational resources as possible to obtain a noninferior set of solutions close to the true Pareto front. These solutions must cover and be uniformly distributed throughout the entire search space. The flowchart of the multiobjective particle swarm optimisation algorithm is displayed in Fig. 10.

Figure 10. Flowchart of the multiobjective particle swarm optimisation algorithm.

To enhance the credibility of the optimisation results, this paper incorporates the simulated annealing algorithm to solve the established multi-objective optimisation model of the route network [25]. The algorithm begins by setting an initial solution and generating another solution randomly within a defined range. The Metropolis sampling acceptance criterion is then utilised, enabling the objective function to fluctuate within a limited range. This process continues until a locally optimal solution is obtained. The determination of the locally optimal solution relies on the control parameter ‘t’, which resembles the temperature ‘T’ in a physical process. The specific procedure of the simulated annealing algorithm is illustrated in Fig 11.

Figure 11. Flow chart of the multiobjective simulated annealing algorithm.

5.0 Example analysis

5.1 Nominal track generation

Radar track data for 1 year from the control sector of the Yunnan region are selected for preprocessing. Next, an adversarial neural network is constructed to train the data, the nominal trajectory is extracted, and the operational route network is generated. The parameters and excitation functions set by the adversarial neural network are described in the following text.

Padding is conducted on the radar trajectory data, with each trajectory containing 1,231 trajectory points. Information on six features is available for each of these points: latitude, longitude, height, heading, velocity and climb speed. After performing padding, the embedding layer is added to the adversarial neural network, and the six-dimensional trajectory data are tiled into one-dimensional data and input into the adversarial neural network generator for training. The first hidden layer is set as a fully connected layer with 1,000 neurons, and the dropout probability is set as 0.3 to prevent the overfitting that occurs when large quantities of data are involved. The second hidden layer is set as a fully connected layer with 2,000 neurons, the excitation functions are ReLU functions and the dropout probability is set as 0.6.

The generator output is employed as the input data of the discriminator, and the ReLU function is selected as the excitation function. The neural network in the discriminator sets a fully connected layer with 2,000 neurons as a hidden layer. The ReLU function is selected as the excitation function, and the dropout probability is set as 0.6. The sigmoid function is used to judge the trained data. If the input is true, the output is 1; otherwise, the output is 0. The cross-entropy function is employed as the loss function.

The nominal trajectory generated from the 1-year radar track data is relatively complex. Mapping this trajectory on a one-day scale enables the clear visualisation of the direction of route operation. Figure 12 is the flight trajectory operation diagram of the aircraft obtained by plotting points based on the radar trajectory data. Figure 14 presents the nominal trajectory map generated using the adversarial neural network. The purple circles in the figure are crucial nodes along the route, including airports, entry or exit points and route intersection points. The red areas around the points are thermal values. A higher intensity of red indicates more traffic between two points in the segment. A comparison of Figs 12 and 13 reveals that the nominal trajectory generated using the adversarial neural network more closely reflects the actual flight direction of aircraft and the traffic on the flight path than does the nominal trajectory generated using the radar track data.

Figure 12. Flight track operation map of the Kunming control sector.

Figure 13. Nominal track map of the Kunming control sector. Figure 12. PSO process.

As shown in Fig. 14, the planned and actual route networks differ considerably. As presented in Table 4, regarding waypoints and route connections in Kunming control sector 2, A581, H121, W143 and H24 are routes in the planned route network, with W143 and H121 being cross routes. However, in the actual route network, no cross flights are observable along routes A581 and H24. Routes W143 and A581 pass through waypoint LPS, and routes H121 and H24 pass through waypoint P249. In the actual route network, routes W143 and A581 are merged, and routes H121 and H24 are merged.

Figure 14. Nominal track map of Kunming control sector 2.

Table 4. Partial plan of route connection in kunming control sector 2

5.2 Evolutionary mechanism of congestion

On the basis of the gradient boosted regression tree model, a regression model is fitted to the three integral congestion indicators of a randomly selected segment for five days by using data from the Yunnan control sector. The effects of the number of altitude adjustments, the number of flights, and the rate of approaching traffic on the congestion status of a segment are determined.

We first analyse the evolutionary pattern of the number of altitude adjustments a segment and the state of congestion state in the segment (Fig. 15). The red and blue lines in Fig. 15 denote the states of congestion and noncongestion, respectively. The previous value of the traffic corresponding to the red line exceeds the previous value of the traffic corresponding to the blue line, which indicates that the more times aircraft appear in a segment, the higher is the probability of congestion in the segment. However, in one scenario, the segment is noncongested when the number of altitude adjustments the segment is high. In another scenario, the segment is congested when the number of altitude adjustments the segment is low. In Fig. 15(a), the number of altitude adjustments at 9:00 and 10:00 is 174 and 151, respectively. At 11:00, this number is 145 but the segment is congested. The number of altitude adjustments at 20:00 is 130, and no congestion is noted. The number of altitude adjustments at 21:00 is 118, the segment is congested. The differences between the congestion index values at the aforementioned time points are not extremely large. Considering that the segment is sometimes congested when the number of altitude adjustments is relatively low but is sometimes noncongested when the number of altitude adjustments is relatively high, the combined effect of the number of altitude adjustments and the other two crucial indicators (i.e., the number of flights and the rate of approaching traffic) has an indirect influence on the state of congestion. Analysis of the altitude adjustments times corresponding to the state of congestion in the segment on different days reveals the following information. When the difference between the number of flights going up and down the segment is large, no congestion is observed. Conversely, when this difference is small, congestion is observed. In Fig. 16(a), the number of altitude adjustments at 9:00 is 174, and the segment is noncongested. However, in Fig. 13(b), the number of altitude adjustments at 11:00, 17:00 and 19:00 is 44, 46 and 46, respectively, and the segment is congested at all three times. Thus, on different days, the number of flights in a given segment cannot be used as an indicator of the state of congestion in the segment. In summary, the number of altitude adjustmentst exerts a strong influence on the evolution of congestion in the segment. The higher the congestion index value, the higher is the possibility of congestion in a segment. The congestion in a segment is also related to the number of flights and the rate of approaching traffic in the segment.

Figure 15. Evolutionary pattern of the number of altitude adjustments and the state of congestion in the segment.

Figure 16. Evolutionary pattern of the number of flights and the state of congestion in a segment.

We first analyse the evolutionary pattern of the number of altitude adjustments a segment and the state of congestion state in the segment (Fig. 15). The red and blue lines in Fig. 15 denote the states of congestion and noncongestion, respectively. The previous value of the traffic corresponding to the red line exceeds the previous value of the traffic corresponding to the blue line, which indicates that the more times aircraft appear in a segment, the higher is the probability of congestion in the segment. However, in one scenario, the segment is noncongested when the number of altitude adjustments the segment is high. In another scenario, the segment is congested when the number of altitude adjustments the segment is low. In Fig. 15(a), the number of altitude adjustments at 9:00 and 10:00 is 174 and 151, respectively. At 11:00, this number is 145 but the segment is congested. The number of altitude adjustments at 20:00 is 130, and no congestion is noted. The number of altitude adjustments at 21:00 is 118, the segment is congested. The differences between the congestion index values at the aforementioned time points are not extremely large. Considering that the segment is sometimes congested when the number of altitude adjustments is relatively low but is sometimes noncongested when the number of altitude adjustments is relatively high, the combined effect of the number of altitude adjustments and the other two crucial indicators (i.e., the number of flights and the rate of approaching traffic) has an indirect influence on the state of congestion. Analysis of the altitude adjustments times corresponding to the state of congestion in the segment on different days reveals the following information. When the difference between the number of flights going up and down the segment is large, no congestion is observed. Conversely, when this difference is small, congestion is observed. In Fig. 16(a), the number of altitude adjustments at 9:00 is 174, and the segment is noncongested. However, in Fig. 15(b), the number of altitude adjustments at 11:00, 17:00 and 19:00 is 44, 46 and 46, respectively, and the segment is congested at all three times. Thus, on different days, the number of flights in a given segment cannot be used as an indicator of the state of congestion in the segment. In summary, the number of altitude adjustmentst exerts a strong influence on the evolution of congestion in the segment. The higher the congestion index value, the higher is the possibility of congestion in a segment.

The congestion in a segment is also related to the number of flights and the rate of approaching traffic in the segment.

We analyse the evolutionary pattern of the rate of approaching traffic and the state of congestion in the flight segment. In Fig. 17, the red and blue lines indicate that the segment is in a congested and noncongested state, respectively. Overall, the rate of approaching traffic corresponding to the red line is higher than that corresponding to the blue line. In Fig. 17(b), the rate of approaching traffic at 8:00 is greater than that at 7:00 and 9:00; however, the segment is noncongested. Considering the other two indicators, the number of altitude adjustments at 8:00 is relatively small in Figs 15(b) and 16(b), and the segment is noncongested. The aforementioned results indicate that the rate of approaching traffic in a noncongested segment does not exert a substantial effect on the state of congestion in the segment. In Fig. 16(a), the number of flights at 10:00 is 19, which is the highest number of flights at a certain time during the considered day. However, the segment is noncongested. At 1:00 and 4:00, the segment is congested even though the number of flights is only ten, which is the smallest number of flights in the segment during a congested period on the considered day. In Fig. 15(a), the number of altitude adjustments in the flight segment at 1:00, 4:00 and 10:00 is 118, 119 and 151, respectively. In Fig. 17(a), the rates of approaching traffic at 1:00, 4:00 and 10:00 are $1.64 \times {10^{ - 9}}$ , $1.13 \times {10^{ - 9}}$ and $4.63 \times {10^{ - 13}}$ , respectively, and the rate of approaching traffic at 10:00 is significantly smaller than are those at 1:00 and 4:00. Thus, in a state of congestion, changes in the rate of approaching traffic strongly affect the congestion status of a segment. The aforementioned analysis indicates that of the three considered indicators, the rate of approaching traffic best reflects the evolutionary mechanism of congestion.

Figure 17. Evolution of the relationship between the rate of approaching traffic and the congestion in a segment.

The effects of individual indicators on congestion is analysed to determine its evolutionary mechanism in the route network. The evolution of the considered indicators in states of congestion and noncongestion is examined in relation to the evolution of congestion in the route network. The values of critical congestion indicators in states of congestion and noncongestion are extracted such that these indicators can be applied to the route network according to the order in which segments along the route are connected. Furthermore, the distribution relationships among the three considered indicators are determined. The evolutionary mechanism of congestion is illustrated in Figs 18 and 19, where the color gradient represents the change in the number of flights.

Figure 18. Evolutionary mechanism of congestion in a segment.

Figure 19. Distribution relationships among the congestion indicators when a segment is noncongested.

Figure 18 displays the changes in the relationships among the three considered indicators when the segment is congested. When the number of flights in the segment is 1–35 and the number of altitude adjustments is 50–400, the rate of approaching traffic is extremely small and changes negligibly. When the number of flights is 35–40 and the number of altitude adjustments is 50–350, a small increase in the rate of approaching traffic occurs. A larger increase in this rate is observed when the number of flights and altitude adjustments is 40–45 and 50–300, respectively. This phenomenon indicates that when the number of flights increases to a certain extent, the rate of approaching traffic increases rapidly. This rate does not peak when the maximum numbers of flights and altitude adjustments are modified; instead, it shifts marginally. Specifically, when the number of flights exceeds the number of flights corresponding to the highest rate of approaching traffic, the air traffic controller slows the rate of approaching traffic.

Figure 19 illustrates the distribution relationships among the three considered indicators in a state of noncongestion. Because Fig. 19 is obtained through scatter interpolation and because no original data are available for the indicators corresponding to a plane, part of the surface appears as a plane. When the number of flights and altitude adjustments is 1–10 and 100–200, respectively, the rate of approaching traffic does not increase. By contrast, this rate increases rapidly and then decreases rapidly when the number of flights and altitude adjustments is 10–20 and 50–200, respectively. When the number of flights and altitude adjustments is 20–50 and 50–200, respectively, the rate of approaching traffic essentially remains unchanged. Table 5 presents a summary of evolutionary relationships between the status and indicators of congestion in a flight segment.

Table 5. Evolutionary relationships between congestion status and congestion indicators in a segment

5.3 Route network optimisation

Kunming control sector 6 encompasses eight regional airports (i.e., Diqing Shangri-La Airport, Ningland Luguhu Airport, Lijiang Sanyi International Airport, Dali Huangcaoba Airport, Panzhihua Bao’anying Airport, Baoshan Yunrui Airport, Tengchong Tuofeng Airport, and Dehong Mangshi Airport) and contains multiple routes, including H87, J512, J514, H90, H117, W162 and X33. This sector is selected for analysis because the eight airports are close to each other and have high flight volumes. Moreover, aircraft flying on the aforementioned sector make numerous altitude adjustments, and the route network is prone to congestion. The optimisation period is set as 1 h, and the parameters of the particle swarm optimisation algorithm are set as follows: inertia weight w = 0.8, number of particles = 1000, self-learning factor c ₁ = 0.1, global learning factor c ₂ = 0.1, and number of iterations = 100.

The particle swarm optimisation algorithm is used to solve the developed multiobjective optimisation model. Figures 20 and 21 present the optimal solution sets obtained by the two fitness functions at 30 iterations, and Fig. 22 displays the Pareto front at this point. Here, input_x1 is the number of flights (an independent variable), and input_x2 is the ratio of diverted flights. Figure 22 indicates that when the particle swarm optimisation algorithm reaches 30 iterations, the blue points are approaching the optimum solution, and a Pareto front has formed.

Figure 20. Solution of objective function 1 at 30 iterations.

Figure 21. Solution of objective function 2 at 30 iterations.

Figure 22. Pareto front at 30 iterations.

Figures 23 and 24 display the optimal solution sets obtained by the two fitness functions at 100 iterations, and Fig. 25 presents the Pareto front at this point. Figure 25 indicates that at 100 iterations, the blue point almost disappears, which suggests that all the particle points are close to the optimal position and that the optimal Pareto solution has been obtained. The Pareto solution set provides various optimisation schemes for Kunming control sector 6. A suitable optimisation scheme can be selected according to actual operations and the preferences of air traffic controllers. Figure 26 presents one of these schemes. The red line represents the newly added one-way line. A one-way line is added to a part of the flight segment of P334-PAN-P259, which reduces the number of head-to-head flights. This step results in the reduction of the rate of approaching traffic in the segment and alleviates congestion in the route network. Similarly, the utilisation of the simulated annealing algorithm for obtaining a solution is presented in Figure 27.

Figure 23. Solution of objective function 1 at 100 iterations.

Figure 24. Solution of objective function 2 at 100 iterations.

Figure 25. Pareto front of objective function 1 at 30 iterations.

Figure 26. Network optimisation scheme for Yunnan control sector 6 with PSO.

Figure 27. Network optimisation scheme for Yunnan control sector 6 with SA.

A comparison of the congestion index values obtained before and after optimisation (Table 6) indicate that the number of flights, the number of altitude adjustments, and the rate of approaching traffic in a segment decrease after optimisation. Thus, the optimisation model for PSO is effective in alleviating congestion in the route network. The average preoptimisation flight time is 485 s, whereas the average postoptimisation flight time is 435 s, which corresponds to a 10% reduction.

Table 6. Congestion index values for a segment before and after optimisation

6.0 Conclusions

In this study, a flight network is generated from flight trajectory data by using deep-learning-related methods, traffic flow characteristics observed in actual flight data are examined, the evolutionary mechanism of congestion in the route network is investigated, and a multiobjective optimisation scheme for mitigating congestion in the network is established. This study is summarised in the following text.

First, a real operational route is generated. Considering the large deviation between planned and actual operational routes, radar track data are preprocessed with the preservation of track characteristics. An adversarial neural network is employed to generate the nominal track from a large volume of radar track data. To obtain the real operational route, rules concerning track distribution and operation are examined.

Second, the evolutionary mechanism of congestion in the route network is explored. By extending road traffic engineering theory to the civil aviation field, we use the flight time as a congestion indicator and select a set of congestion indicators by integrating the operational characteristics of actual aircraft with the network operation status. Under a gradient boosted regression tree model, a set of three critical congestion indicators is established. We analyse the relationships between these indicators and the flight time, investigate the patterns of change in these indicators and in the state of congestion in different flight segments at different times, and examine the evolutionary mechanism of congestion in the route network.

Third, we construct an optimisation scheme for alleviating congestion in the route network. On the basis of the evolutionary mechanism of congestion, a multiobjective optimisation model is established with the objectives of minimising operating and optimisation costs, the number of altitude adjustments, and the rate of approaching traffic as well as maximising the proportion of same-direction flights among the diverted flights. A particle swarm optimisation algorithm is executed to solve the developed model and obtain a set of Pareto optimal solutions. The model can be employed to select different optimal solutions according to the actual operation of the route network. Furthermore, different optimal solutions can be applied to different congested segments such that congestion can be mitigated to the highest extent.

Because of the high complexity of air networks and the high dimensionality and large scale of radar track data, research on methods for the intelligent analysis of historical radar track data is incipient. On the basis of the present study, further investigations can be conducted on different aspects.

First, this study is a preliminary exploration of the evolutionary mechanism of congestion in air networks. Because of limitations in time and energy, the influence of the congestion propagation law on the congestion mechanism is not considered. Future studies can investigate this influence.

Second, in the developed multiobjective optimisation model, the effects of meteorological factors, major events, and the flight performance of actual aircraft on congestion are not considered. Future studies can incorporate such data into the developed optimisation model to improve its performance.