NOMENCLATURE
- A
arrival
- AD
airborne delay
- ADP
aircraft departure problem
- ALP
aircraft landing problem
- ASSP
aircraft sequencing and scheduling problem
- ATCos
air traffic controllers
- D
departure
- ETAs
estimated times of arrival
- EV
expected value
- FCFS
first-come first-served
- GD
ground delay
- GH
ground holding
- H
heavy
- LTFJ
Istanbul Sabiha Gökçen International Airport
- M
medium
- MILP
mixed-integer linear programming
- PMS
point merge system
- PSA
position-shifted aircraft
- RP
here and now
- SA
simulated annealing
- SP
stochastic programming
- TAS
true airspeed
- TMA
terminal manoeuvring area
- VSS
value of stochastic solution
1.0 INTRODUCTION
The volume of air transportation has grown considerably due to recent economic developments and the increased number of low-cost airlines. The aviation industry is resistant to external effects, and thus the amount of air traffic has more than doubled in the past 20 years(1). This increasing traffic demand also raises the importance of issues regarding airport and airspace management. If the development and enhancement of airspaces and airports are not prioritised, there will be a significant increase in both airborne and ground delays. The aircraft sequencing and scheduling problem (ASSP) is connected with the reduction of these delays using different air traffic flow management and separation techniques, such as ground holding (GH), the point merge system (PMS), vectoring and speed reduction. These techniques allow air traffic controllers (ATCos) to re-sequence aircraft to maintain safe and orderly air traffic flow. Although re-sequencing results in an additional workload for ATCos, its effective use contributes to significant reductions in delay and congestion, not only at a specific airspace or airport but also across the entire air traffic system.
In this study, we aim to obtain resilient aircraft sequencing that takes wind direction changes into account in the terminal manoeuvring area (TMA). The ASSP includes two sub-problems: the aircraft landing problem (ALP) and the aircraft departure problem (ADP). The majority of previous studies, on the other hand, have presented solution approaches for either ALP or ADP, however, airports provide services for both arrival and departure operations. Bennell etal.(Reference Bennell, Mesgarpour and Potts2) and Bianco etal.(Reference Bianco, Dell’Olmo and Giordani3) provided comprehensive reviews regarding the ALP and ADP. Different deterministic and stochastic meta-heuristic algorithms have been applied to these problems to achieve good solutions in a short computational time. The frequently used algorithms are tabu search (Atkin etal.(Reference Atkin, Burke, Greenwood and Reeson4), D’Ariano etal.(Reference D’ariano, Pistelli and Pacciarelli5), Furini etal.(Reference Furini, Kidd, Persiani and Toth6), Samà etal.(Reference Samà, D’ariano, D’ariano and Pacciarelli7)), simulated annealing (Hancerliogullari etal.(Reference Hancerliogullari, Rabadi, Al-Salem and Kharbeche8), Salehipour etal.(Reference Salehipour, Moslemi and Kazemipoor9), Liang etal.(Reference Liang, Delahaye and Marechal10)), variable neighbourhood search (Salehipour etal.(Reference Salehipour, Modarres and Naeni11)), ant colony optimisation (Jiang etal.(Reference Jiang, Xu, Xu, Liao and Luo12), Zhan etal.(Reference Zhan, Zhang, Li, Liu, Kwok, Ip and Kaynak13)), the genetic algorithm (Beasley etal.(Reference Beasley, Sonander and Havelock14), Hu and Chen(Reference Hu and Chen15), Hu and Di Paolo(Reference Hu and Di Paolo16), Hu and Di Paolo(Reference Hu and Di Paolo17)), scatter search (Pinol and Beasley(Reference Pinol and Beasley18)), and particle swarm optimisation (Hong etal.(Reference Hong, Cho, Kim and Choi19)). The following reviewed studies focus on stochastic approaches to the ASSP: Solveling etal.(Reference Solveling, Solak, Clarke and Johnson20) developed two-stage stochastic programming for stochastic runway scheduling including uncertainties in pushback delay, taxiing and deviations from estimated times of arrival (ETAs). Their findings indicated that a stochastic approach performed better than deterministic or first-come first-served (FCFS) approaches. Solveling and Clarke(Reference Sölveling and Clarke21) presented a stochastic branch-and-bound algorithm to obtain an optimum or near-optimum solution that aimed to minimise the total makespan. The findings revealed that such models can reduce the production time by between 5% and 7% compared with the deterministic approach. Choi etal.(Reference Choi, Mulfinger, Robinson and Capozzi22) developed a mixed-integer linear programming model for the design of the route structure in extended TMA operations. They also proposed a genetic algorithm to minimise CPU time. Bosson and Sun(Reference Bosson and Sun23) proposed a multi-stage stochastic programming model and the sample average approximation method to obtain a solution. According to their results, stochastic optimisation methods provide better sequencing than the FCFS policy. Heidt etal.(Reference Heidt, Helmke, Kapolke, Liers and Martin24) presented a stochastic runway scheduling model. The results showed that their robust model made fewer changes to aircraft sequences and produced more stable traffic plans. Ng etal.(Reference Ng, Lee, Chan and Qin25) presented a stochastic programming model using the min–max regret approach for mixed operations runway. The artificial bee colony algorithm was applied to obtain a feasible solution. Solak etal.(Reference Solak, Sölveling, Clarke and Johnson26) presented two different stochastic runway scheduling models and tested them using realistic data. Their conclusion was that stochastic models are very useful and can be implemented easily to solve the ASSP. Liu etal.(Reference Liu, Liang, Zheng, Chu and Chu27) proposed a two-stage stochastic programming model. The first stage obtains the aircraft weight class sequence, then the second stage assigns the individual flight to the aircraft sequence. Similarly, Khassiba etal.(Reference Khassiba, Bastin, Gendron, Cafieri and Mongeau28) proposed a two-stage stochastic programming model for a single-runway airport. Their findings indicated that their stochastic model reduced the need for a holding stack by a few hours by offering more upstream linear holding. Hong etal.(Reference Hong, Choi and Kim29) proposed a two-stage stochastic algorithm using mixed-integer linear optimisation for the ALP. Due to the complexity of the problem, particle swarm optimisation was implemented. The problems of the first and second stages are sequencing and scheduling. Also, they used the PMS to regulate arrival traffic where uncertainties emerge in descent times from sequencing legs to the final approach point.
In this study, a stochastic mixed-integer linear programming model based on the simulated annealing algorithm is presented to minimise total aircraft delay for a single-runway airport. The ASSP can easily be affected by external disturbances, so we propose an algorithm that includes wind direction uncertainties for both upper and lower altitudes of sequencing legs. While Hong etal.(Reference Hong, Choi and Kim29) only considered uncertain aircraft arrival times from sequencing legs to the merge point for a single runway, we focus on the entire TMA, including all flight levels, and utilise a mixed operation runway.
Our contributions are as follows: (1) the model takes wind direction uncertainties based on real historical data into consideration and integrates them into the model. To the best of our knowledge, the effects of wind direction have not been investigated in the ASSP literature. (2) The model aims to obtain resilient aircraft sequencing, thus avoiding the need for re-sequencing. (3) Two independent wind directions for above and below the sequencing legs are adapted to the model, and an appropriate number of scenarios are generated to test the stochastic programming.
The remainder of this manuscript is structured as follows: Section 2 defines the problem and mathematical model; Sections 3 and 4 provide the simulated annealing algorithm and computational results, then Section 5 presents the conclusions.
2.0 PROBLEM DEFINITION
The aircraft sequencing and scheduling problem is a real challenge for ATCos, especially when they aim to determine optimum sequencing that maximises or minimises their objectives. Obtaining proper sequencing can reduce their workload because the need to re-schedule aircraft disappears during TMA operations. ATCos can regulate the arrival traffic using vectoring and speed reduction techniques, holding stacks or the point merge system within TMAs. The PMS is a systematic method that allows ATCos to sequence arrival traffic flows without using vectoring techniques(Reference Favennec, Rognin, Trzmiel, Vergne and Zeghal30). In this system, ATCos separate arriving aircraft systematically in sequencing legs and then direct them to a merge point. The benefits of the PMS can be listed as: reduced ATCo workload, increased efficiency and predictability for aircraft trajectories, and therefore improved traffic management within TMAs(Reference Boursier, Favennec, Hoffman, Trzmiel, Vergne and Zeghal31). Numerous studies have presented the implementation of the PMS technique to various ASSPs, by Boursier etal.(Reference Boursier, Favennec, Hoffman, Trzmiel, Vergne and Zeghal31), Ivanescu etal.(Reference Ivanescu, Shaw, Tamvaclis and Kettunen32), Hong etal.(Reference Hong, Choi and Kim29) Liang etal.(Reference Liang, Delahaye and Marechal10) and Sahin etal.(Reference Sahin, Usanmaz and Turgut33). In this study, Istanbul Sabiha Gökçen International Airport (LTFJ) and its surrounding TMA is selected for modelling, based on the airspace and runway layout of 2019(34). PMS has already been implemented at LTFJ. The airport has a single runway 06/24, which denotes the runway directions, as the angles on the horizontal plane from magnetic north are 60° and 240°. This runway serves for mixed operations of arrivals and departures. It is one of the busiest single-runway airports in the world, with passenger use of 35.5 million in 2019(35). Arrivals enter the TMA from six different entry points: ATVEP, GELBU, TURCO, PAZAR, EVNOT and ELVON (Fig. 1). In this study, we assume that both arrivals and departures use only runway orientation 06. The PMS consists of two sequencing legs separated vertically by 1000feet and a merge point. All arriving aircraft fly with a constant airspeed along the sequencing legs located at FL100 and FL110 according to the location of their entry points. Before entering the sequencing legs, safe separation among the arrival aircraft on the same leg can be maintained using vectoring or airspeed reduction. In addition to these methods, sequencing legs can also provide the required separations between arrival (A) and departure (D) operations. Two different aircraft performance categories are also considered in this study according to their wake turbulence categories: heavy (H) and medium (M). Aircraft are directed to the merge point after the time separation between an aircraft pair is achieved. The PMS includes three important times: the entry time to the sequencing leg, 
${q_i}$; the travelling time along the sequencing legs 
${x_i}$; and the touchdown time on the runway 
${t_i}$.
PMS route structure of LTFJ within the Istanbul TMA.
Airspeed values, the distance between the entry point and sequencing legs and the distance between the sequencing legs and merge point are very important parameters to detect aircraft conflict situations. In the method proposed herein, such conflicts are detected by checking the crossing times of aircraft at each critical point along the sequencing legs and on the runway. The wind direction is important in the calculation of these crossing times, because the airspeed of aircraft may increase or decrease according to the wind direction. Airspeeds are expected to be stochastic because of the wind direction.
It is assumed that a set of I = (1, 2… n) of i aircraft are in the TMA and that the airport is used for arrivals and departures. Aircraft can use the runway for either departure (D) or arrival (A) operations without violating the vortex separation times given in Table 1.
Time separations, in seconds(Reference Sölveling and Clarke21)
2.1 Stochastic programming
The air traffic system is easily affected by uncertainties and external disturbances. These situations cause extra airborne and ground delays as well as re-sequencing in TMA operations. Nevertheless, punctuality is very important for an airport to maintain safe, efficient and economic flight operations. Therefore, steady scheduling is of great importance for both airport and airspace operations. Stochastic programming (SP) takes such uncertainties into considerations and produces more resilient solutions to handle these disturbances. Uncertain parameters can be obtained using historical data or expert opinion. In general, the stochastic programming approach makes decisions considering the uncertainties in the future. Uncertainties can be represented in two different ways. The first one is the continuous probability distribution where numerical integration is employed over the random continuous probability space. This approach limits the model size but the problem presents nonlinearities and computational difficulties. The second approach is the scenario-based stochastic approach, where the random space is considered to comprise discrete events. The main difficulty of this approach is the considerable increase in computational requirements with increasing number of uncertain parameters(Reference Al-Qahtanii and Elkamel36).
As indicated by Kaya etal.(Reference Kaya, Bagci and Turkay37), the scenario-based stochastic optimisation approach generates different sets of variables for each particular scenario S. In this approach, which was originally proposed by Dantzig(Reference Dantzig38), the probabilities for each scenario are assigned and a solution for each scenario is obtained (see Birge and Louveaux(Reference Birge and Louveaux39), Kall and Wallace(Reference Kall and Wallace40)). Nevertheless, the optimum solution for a particular scenario S might not be the optimum for the general problem. We consider a scenario-based stochastic programming approach, in which the uncertainties are modelled using several scenarios that can occur, with their respective probabilities. The general formulation of a stochastic programming is given below, where f(x) is the scenario-independent cost and Q(x,S) is the cost when the scenario S is realised. The objective function in Equation (1) of the problem consists of the value of f(x) and the expected value of Q(x,S) with respect to the possible scenarios. Some of the constraints of the model might be scenario independent, as given in Equation (2), or scenario dependent, as given in Equation (3).
\begin{equation}{\rm{Mi}}{{\rm{n}}_{\rm{x}}}{\rm{}}( {{\rm{f}}( {\rm{x}} ) + {{\rm{E}}_{\rm{s}}}[ {{\rm{Q}}( {{\rm{x}},{\rm{S}}} )} ]} )\end{equation}
\begin{equation}{\rm{subject\ to\ g}}( {\rm{x}} ) = {\rm{a}}\end{equation}
\begin{equation}{\rm{h}}( {{\rm{x}},{\rm{S}}} ) = {\rm{b}}( {\rm{S}} )\,\forall S\end{equation}
In our problem, two different altitude layers are considered, in which the wind directions show independent uncertainties. First, we determined the wind directions based on weather sounding observation data at 12:00 pm for Istanbul between 2015 and 2019(41). Then we obtained the histogram for two altitude layers, i.e., layer 1 (25,000 to 11,000ft) and layer 2 (11,000 to 3,000 t), as shown in Fig. 2.
The wind direction histograms: (a) layer 1 and (b) layer 2.
In this study, eight different wind directions were determined to represent all wind directions (20°, 70°, 110°, 160°, 200°, 250°, 290° and 330°) with reasonable discrete wind directions. The probabilities of the wind direction are given in Table 2.
The wind direction probabilities according to wind direction
As a result, the model has 64 different scenarios for the stochastic programming model. The wind speed is calculated using Equation (4) presented by Turgut and Usanmaz(Reference Turgut and Usanmaz42), where the altitude is in metres:
\begin{equation}Wind\ speed\!\left( {m/s} \right) = - 9.2 \times {10^{ - 8}}\!\left[ {Altitud{e^2}} \right]+3.7 \times {10^{ - 3}}\!\left[ {Altitude} \right]+6.5\end{equation}
The ground speed of the aircraft can be expressed as
\begin{equation}{V_{ground}} = {V_{TAS}} + {V_{wind}}\cos \!\left( {Aircraft_{{heading}} - Win{d_{direction}}} \right)\end{equation}
In Equation (5), 
${V_{TAS}}$ is the true airspeed, 
${V_{wind}}$ is the wind speed, and 
$Aircraft_{{heading}}$ and 
$Win{d_{direction}}$ are the aircraft heading angle and wind direction, respectively.
To test our model, we use the expected value (EV), the here-and-now solution (RP) and the value of the stochastic solution (VSS) to analyse the benefits of the stochastic programming model. Firstly, without taking wind direction uncertainties into considerations, we solve the ASSP using deterministic programming, assuming that the wind direction is known and fixed. We use the expected value of the possible wind directions as the deterministic parameter for the wind direction in this model and obtain an aircraft sequencing as a result. The values in column EV of Table 5 show the actual expected value of the objective function with the obtained sequence when wind direction uncertainties exist. These values are calculated by finding the objective function value of the stochastic ASSP model for the fixed aircraft sequence decisions obtained from the deterministic model. Secondly, we take wind direction uncertainties into account when developing and solving the model, instead of assuming a fixed wind direction. Considering various scenarios for possible wind directions, the stochastic ASSP model determines the best aircraft sequencing and calculates the expected total delay considering different scenarios according to the probability of the wind directions. We solve the stochastic ASSP model directly by considering the randomness of the wind direction, and the objective function value of this model is shown in column RP. The difference between EV and RP is shown in column VSS, and this value reflects the benefits of considering the randomness in wind directions via stochastic programming, instead of using the expected value of the wind direction as a deterministic parameter, and using decisions obtained by deterministic programming.
\begin{equation}VSS = EV - RP\end{equation}
2.2 Mathematical model
The objective is to minimise the total aircraft delay. The entire mathematical model is given below for the ASSP.
Indices
- I:
set of aircraft in the sector
$i,{i_1},{i_2} \in l$- M:
set of entry points m M
- K:
set of aircraft performance categories
${k_1}$ ${k_2} \in$K- O:
set of operation types
${o_1}$, ${o_2}\in U$- L:
set of wind directions
$l,{l_1},{l_2}\in L$
Parameters
- M:
very large positive number
- ${g_i}$:
entry time into the TMA for aircraft i
- $p{x_i}$:
side the aircraft approaches from i
- ${r_i}$:
entry point into the sector for aircraft i
- ${p_i}$:
operation type of aircraft i
- $ap{c_i}$:
performance category of aircraft i
- ${v_{m,l}}$:
airspeed value between sector entry point and sequencing legs for an aircraft using entry point m in scenario l
- $v{d_l}$:
airspeed value between the sequencing leg and the merge point for an aircraft in scenario l
- $vmf$:
airspeed value between the merge point and the runway
- $d{h_m}$:
distance from entry point m to sequencing leg
- dl:
distance from the sequencing leg to the merge point
- dmf:
distance from the merge point to the runway
- $p{h_l}$:
probability of the wind direction of l for layer 1
- $p{l_l}$:
probability of the wind direction of l for layer 2
- $se{p_{{o_1},{o_2},{k_1},{k_2}}}$:
separation time between operation type of
${o_1}$ and ${o_2}$ with aircraft performance category ${k_1}$ and ${k_2}$
Decision variables
- ${q_{i,{l_1},{l_2}}}$:
entry time onto the sequencing leg for aircraft i with wind direction at layer 1 being
${l_1}$ and ${l_2}$ at layer 2- ${t_{i,{l_1},{l_2}}}$:
runway use time of aircraft i with wind direction at layer 1 being
${l_1}$ and ${l_2}$ at layer 2- ${h_{i,{l_1},{l_2}}}$:
delay time before entering the sequencing leg for aircraft i with wind direction
${l_1}$ for layer 1 wind direction and wind direction ${l_2}$ for layer 2- ${x_{i,{l_1},{l_2}}}$:
travel time on the sequencing leg for aircraft i with wind direction at layer 1 of
${l_1}$ and ${l_2}$ at layer 2- ${w_{i,{l_1},{l_2}}}$:
ground delay time for aircraft i with wind direction
${l_1}$ at layer 1 and ${l_2}$ at layer 2- $dela{y_{i,{l_1},{l_2}}}$:
total delay time of aircraft i with wind direction
${l_1}$ at layer 1 and ${l_2}$ at layer 2- ${c_{{i_1},{i_2}}}$:
0–1 variable that takes a value of 1 if aircraft
${i_1}$ uses the runway before aircraft ${i_2}$; otherwise, it is zero
\begin{equation}\min {\mathop \sum \limits_i \mathop \sum \limits_l \mathop \sum \limits_{{l_2}} p{h_l} \cdot p{l_{{l_2}}} \cdot dela{y_{i,l,{l_2}}}}\end{equation}
\begin{equation}{q_{i,{l_1},{l_2}}} = {g_i} + {h_{i,{l_1},{l_2}}} + {{d{h_m}} \over {{v_{m,{l_1}}}}}\forall i,{l_1},{l_2},m|{p_i}=2,m = {r_i}\end{equation}
\begin{equation}{q_{{i_2},{l_1},{l_2}}} - {q_{{i_1},{l_1},{l_2}}} \ge sreq - {c_{{i_1},{i_2}}}M\end{equation}
\begin{align}&\forall {i_1},{i_2},{l_1},{l_2}|{i_1} \ne {i_2},p{x_{{i_1}}} = p{x_{{i_2}}},{p_{{i_1}}}=2,{p_{{i_2}}}=2 \nonumber\\[3pt] & {q_{{i_1},{l_1},{l_2}}} - {q_{{i_2},{l_1},{l_2}}} \ge sreq - (1 - {c_{{i_1},{i_2}}})M \end{align}
\begin{align}& \forall {i_1},{i_2},{l_1},{l_2}|{i_1} \ne {i_2},p{x_{{i_1}}} = p{x_{{i_2}}},{p_{{i_1}}}=2,{p_{{i_2}}}=2 \nonumber\\[3pt] & {t_{i,{l_1},{l_2}}} = {q_{i,{l_1},{l_2}}} + {x_{i,{l_1},{l_2}}} + {{dl} \over {v{d_{{l_2}}}}} + {{dmf} \over {vmf}}\forall i,{l_1},{l_2},m|{p_i}=2,m = {r_i}\end{align}
\begin{equation}{t_{i,{l_1},{l_2}}} = {g_i} + {w_{i,{l_1},{l_2}}}\forall i,{l_1},{l_2},m|{p_i}=1\end{equation}
\begin{equation}{t_{{i_2},{l_1},{l_2}}} - {t_{{i_1},{l_1},{l_2}}} \ge se{p_{{o_1},{o_2},{k_1},{k_2}}} - {c_{{i_1},{i_2}}}M\end{equation}
\begin{align} &\forall {i_1},{i_2},{l_1},{l_2},{o_1},{o_2},{k_1},{k_2}|{i_1} \ne {i_2},{p_{{i_1}}} = {o_1},{p_{{i_2}}} = {o_2},ap{c_{{i_1}}} = {k_1},ap{c_{{i_2}}} = {k_2} &\quad {t_{{i_1},{l_1},{l_2}}} - {t_{{i_2},{l_1},{l_2}}} \ge se{p_{{o_1},{o_2},{k_1},{k_2}}} - (1 - {c_{{i_1},{i_2}}})M\end{align}
\begin{align} &{i_1},{i_2},{l_1},{l_2},{o_1},{o_2},{k_1},{k_2}|{i_1} \ne {i_2},{p_{{i_1}}} = {o_1},{p_{{i_2}}} = {o_2},ap{c_{{i_1}}} = {k_1},ap{c_{{i_2}}} = {k_2} \nonumber\\[3pt] &\quad dela{y_{i,{l_1},{l_2}}} = {x_{i,{l_1},{l_2}}} + {h_{i,{l_1},{l_2}}} + {w_{i,{l_1},{l_2}}}\quad\forall i,{l_1},{l_2} \end{align}
\begin{equation}{q_{i,{l_1},{l_2}}} \ge 0\quad \forall i,{l_1},{l_2}\end{equation}
\begin{equation}{w_{i,{l_1},{l_2}}} \ge 0\quad \forall i,{l_1},{l_2}\end{equation}
\begin{equation}dela{y_{i,{l_1},{l_2}}} \ge 0\quad \forall i,{l_1},{l_2}\end{equation}
\begin{equation}{t_{i,{l_1},{l_2}}} \ge 0\quad \forall i,{l_1},{l_2}\end{equation}
\begin{equation}{x_{i,{l_1},{l_2}}} \ge 0\quad \forall i,{l_1},{l_2}\end{equation}
\begin{equation}{c_{{i_1},{i_2}}}\, \in \left\{ {0,1} \right\}\quad \forall {i_1},{i_2}\ \ \end{equation}
The objective of Equation (7) is to minimise the total delay. The constraint set in Equation (8) calculates the entry time onto the sequencing legs and the airborne delays before entering the sequencing leg for each aircraft with a wind direction of 
$l_1$ in layer 1 and 
${l_2}$ in layer 2. The constraint sets in Equations (9) and (10) check the conflicts for arriving aircraft on the sequencing legs according to their routes. The constraint set in Equation (11) calculates the runway touchdown time of the arrival and the travel time on the sequencing leg for each arrival. Similarly, the constraint set in Equation (12) calculates the ground delay of each departure for a wind direction of 
${l_1}$ in layer 1 and 
${l_2}$ in layer 2. The constraint sets in Equations (13) and (14) control the conflicts for both arrival and departure operations on the runway according to the vortex separation times given in Table 1. The constraint set in Equation (15) calculates the delay time of each aircraft for a wind direction of 
${l_1}$ in layer 1 and 
${l_2}$ in layer 2. The delay is calculated as the sum of three different delay components. The first two components are the two airborne delays of arrivals, while the third component reflects the ground delay (
${w_{i,{l_1},{l_2}}}$) for departures. Airborne delays can occur before an aircraft enters the sequencing leg (
${h_{i,{l_1},{l_2}}}$) and/or during its flight along the sequencing leg 
$\left( {{x_{i,{l_1},{l_2}}}} \right)$. An aircraft suffers airborne delays if it is an arrival, or a ground delay if it is a departure. The constraint sets in Equations (16)–(21) are sign constraints.
3.0 SIMULATED ANNEALING
It is hard to reach an optimum solution using exact algorithms such as the GAMS/CPLEX solver for large-scale problems in a short time. Therefore, the simulated annealing (SA) algorithm is implemented to solve the ASSP. SA is a single solution-based metaheuristic algorithm and uses a local search approach to obtain optimum or near-optimum solutions. It begins with an initial solution, then it creates neighbourhood solutions from the current solution S. The best neighbourhood solution Sʹ is then selected and compared with S. If Sʹ has a better fitness value than S, Sʹ it is selected as the current solution. However, if Sʹ shows weaker performance than the current solution, the solution is determined according to the probability:
\begin{equation}P = {e^{{\raise0.7ex\hbox{$_{ - \Delta }$} \!\mathord{\left/ {\vphantom {\int{ - \Delta } T}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$_T$}}}}\end{equation}
A random value between 0 and 1 is generated to select either S or Sʹ for the next iteration. If the probability is greater than this random value, Sʹ is selected as the new current solution. In Equation (22), Δ is the difference between the fitness function values of Sʹ and S such that, 
$\Delta\: f(S) - f( {S'} )$ and T is the temperature. The temperature steadily reduces from an initial value to zero. At each step, the algorithm randomly choses a solution close to the current one, evaluates its quality and proceeds to it according to the temperature-dependent probability of choosing better or worse solutions. The temperature in each iteration is estimated using the formula
\begin{equation}{T_{i+1}} = \beta {T_i}\end{equation}
The cooling schedule rate and iteration number are represented using 
$\beta$ and i in Equation(23). The temperature reduces slowly, which enables poorer solutions to be preferred in the early iterations. When the temperature reaches zero, the probability of preferring poorer solutions decreases considerably. Two different neighbourhood generation processes are applied in this problem. The first one switches the sequence of two aircraft randomly, while the second one reverses the entire sequence of aircraft between a randomly selected aircraft pair. These neighbourhood generation processes are presented in Figs 3 and 4.
Examples of generating a neighbourhood solution using only two aircraft.
Examples of generating a neighbourhood solution by reversing the entire sequence between a randomly selected aircraft pair.
4.0 COMPUTATIONAL RESULTS
The aircraft sequencing and scheduling model is formulated using stochastic mixed-integer linear programming (MILP) and tested using the GAMS/CPLEX solver for three different air traffic flow rates: low, medium and high traffic demand (i.e. 30, 37 and 44 aircraft per hour, respectively). Ten different small-scale test problems are generated, including route, entry time, entry point, aircraft performance category and aircraft operation type information. While the entry times of each aircraft are obtained using an exponential distribution, the aircraft performance category and route assignments are generated using historical data. The GAMS/CPLEX solver reaches optimum results for all test problems, then the SA metaheuristic is applied to the test problems to evaluate the performance of the proposed algorithm. The parameters used in the SA are presented in Table 3. A computer with a 2.3 GHz Intel Core i7 processor and 16 GB of RAM is used in all the computations. MATLAB software is used for the metaheuristic algorithm.
Simulated annealing parameters
The results for the test problems are presented in Table 4, where the first column shows the small-scale test problem number, and the second and third columns show the GAMS/CPLEX solver total aircraft delay (Z) and CPU times in seconds. Similarly, the final columns show the SA total aircraft delay (Z) and CPU times in seconds for the test problems. For small-scale problems (ten aircraft), GAMS/CPLEX can find a global optimum solution in a short time, as shown in Table 4. However, as the size of the problem increases (i.e. 30 aircraft), as also seen in Table 4, the GAMS/CPLEX solver returns a solution with an objective function of more than 106 in 300s. On the other hand, the SA algorithm can provide much better solutions in much shorter CPU times. Therefore, we use the SA algorithm for the larger-sized problems to obtain better solutions in reasonable time limits.
Test problem results for 10 aircraft per 20min and 30 aircraft per hour
The results of the large-scale test problems are given in Table 5, where the first column shows the large-scale test problem number, and the second to fifth columns show the EV, RP, VSS and CPU times for 30 aircraft. Similarly, the results for 37 and 44 aircraft are presented in the following columns.
Results of the SA algorithm for three different air traffic flow rates per hour
The proposed stochastic programming metaheuristic algorithm found a VSS value of 8, 15 and 15 for 30, 37 and 44 aircraft. The model a VSS value of zero for some scenarios. This indicates that the stochastic programming model finds the same aircraft sequencing as the deterministic programming model. Also, this result indicates that there is no real benefit in considering wind direction uncertainties for these scenarios. The average VSS results for 30, 37 and 44 aircraft per hour are 80.2, 399.8 and 727.3 in the current version. These results are shown in the last row of Table 5. This indicates that, when the traffic increases, the average value of VSS also rises. The average percentage gains in the total aircraft delay are 6.7%, 16.7% and 15.2% when comparing the average results of EV and RS. The benefits of stochastic programming are more obvious with high traffic demand, but the use of stochastic programming reduced the total aircraft delay for all three traffic demand levels. Table 6 presents how the ground and airborne delays are distributed among the given aircraft sets in the stochastic programming solutions. The values of PSA, on the other hand, indicate the number of aircraft for which the sequence was changed when compared with the deterministic programming aircraft sequence.
GD, AD and PSA results of stochastic programming
The proposed algorithm results in an average airborne aircraft delay that is higher than the average ground aircraft delay for all traffic rates. The average airborne delays are 55.4%, 22.9% and 11.8% higher than the average total ground delays. Also, it is apparent that the use of stochastic programming results in some changes in aircraft sequencing and that the amount of position shifting rises as the number of flights is increased.
5.0 CONCLUSION
The ASSP is a real challenge for ATCos, and providing resilient aircraft sequencing reduces their workload. A stochastic programming model using the SA algorithm is presented herein, considering wind direction uncertainties. The performance of the stochastic algorithm is evaluated using an appropriate number of scenarios representing all wind directions. The results of the stochastic and deterministic approaches are then compared, revealing that the stochastic version of the model can noticeably reduce airborne delays. Additionally, stochastic programming is flexible to tolerate sudden wind direction changes and eliminates the need for re-sequencing. Furthermore, the proposed algorithm can produce a practical aircraft sequencing solution and, as a result, be used by ATCos to obtain fast and feasible sequencing. Because of the robust aircraft sequencing obtained from stochastic programming, it is easy to calculate the landing or departure times of each aircraft using the proposed model by taking into account the specific wind direction for the specific time window. ATCos can direct the aircraft set in the TMAs to complete their operations without the need for re-sequencing if any wind direction change occurs. As the aircraft sequencing is resilient to changing wind direction, the model can easily find a new schedule for each aircraft using the current wind direction value. However, the algorithm needs to be improved to be implemented in real-time operations. In future studies, a dynamic version of this stochastic model can be studied.
ACKNOWLEDGEMENTS
This research received no specific grant from any funding agency in public, commercial or non-profit sectors.
CONFLICTS OF INTEREST
The authors declare that they have no conflicts of interest.
