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.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 [Reference Zhang, Gu, Zhao and Zhao14].

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 [Reference Hinton, Srivastava, Krizhevsky, Sutskever and Salakhutdinov15]. 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.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 [Reference Li20]. 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 [Reference Zhu21]. 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.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 [Reference Ghosh and Cabrera22].

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 [Reference Fan23]. 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.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.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.