Nomenclature

A

arrival

ASSP

aircraft sequencing and scheduling problem

ATCOs

air traffic control officers

CPS

constrained position shifting

D

departure

EV

expected value

FAP

final approach point

FCFS

first-come first-served

LB

lower bound

LTFJ

Istanbul Sabiha Gökçen International Airport

MILP

mixed-integer linear programming

NM

nautical miles

PMS

point merge system

PSO

particle swarm optimisation

RP

here-and-now solution

SA

simulated annealing

SAA

sample average approximation

SP

stochastic programming

TMAs

terminal manoeuvring areas

TS

tabu search

UB

upper bound

VSS

value of stochastic solution

2.0 Literature review

This paper aims to obtain resilient aircraft sequencing while considering the uncertainties of flight durations in the TMAs. The ASSP consists of both aircraft landing and take-off problems. Comprehensive literature reviews were presented in [Reference Bianco, Dell’olmo and Giordani3–Reference Zhang, Zhao, Zhang, Dai and Sui5]. The literature can be divided into two main groups, using either exact or metaheuristic approaches. Commonly used exact solution approaches are mixed-integer programming [Reference Dönmez, Çetek and Kaya6–Reference Cecen and Çetek14], the branch and bound algorithm [Reference Sölveling and Clarke15] and dynamic programming [Reference Montoya, Rathinam and Wood16–Reference Atkin, Burke, Greenwood and Reeson18]. The most common metaheuristic algorithms are the tabu search [Reference D’ariano, Pistelli and Pacciarelli19–Reference Hancerliogullari, Rabadi, Al-salem and Kharbeche22] (TS), simulated annealing [Reference Liang, Delahaye and Marechal23–Reference Cecen, Cetek and Kaya26] (SA), variable neighborhood search [Reference Salehipour, Moslemi and Kazemipoor27], ant colony optimisation [Reference Jiang, Xu, Xu, Liao and Luo28, Reference Zhan, Zhang, Li, Liu, Kwok, Ip and Kaynak29], the genetic algorithm [Reference Beasley, Sonander and Havelock30–Reference Liu, Liang, Zheng, Chu and Chu34], scatter search [Reference Pinol and Beasley35] and particle swarm optimisation [Reference Hong, Choi and Kim36, Reference Hong, Cho, Kim and Choi37] (PSO). Moreover, the heuristic algorithm [Reference Sadovsky and Windhorst38] is also implemented for this problem. In addition to the different solution approaches, the objectives in related studies are not the same, these include the following: minimising total cost [Reference Beasley, Krishnamoorthy, Sharaiha and Abramson7, Reference Briskorn and Stolletz9, Reference Faye11, Reference Murça and Müller12, Reference Lieder, Briskorn and Stolletz17, Reference Pinol and Beasley35, Reference Hong, Cho, Kim and Choi37, Reference Balakrishnan and Chandran39–Reference Solveling, Solak, Clarke and Johnson46], minimising total delay [Reference Bianco, Dell’olmo and Giordani3, Reference Dönmez, Çetek and Kaya6, Reference Furini, Persiani and Toth8, Reference Cecen and Çetek14, Reference Montoya, Rathinam and Wood16, Reference Atkin, Burke, Greenwood and Reeson18–Reference Furini, Kidd, Persiani and Toth20, Reference Hancerliogullari, Rabadi, Al-salem and Kharbeche22, Reference Liang, Delahaye and Marechal23, Reference Rodríguez-Díaz, Adenso-Díaz and González-Torre25–Reference Zhan, Zhang, Li, Liu, Kwok, Ip and Kaynak29, Reference Hu and Chen31, Reference Liu, Liang, Zheng, Chu and Chu34, Reference Bennell, Mesgarpour and Potts42, Reference Gupta, Malik and Jung47–Reference Eun, Hwang and Bang50], minimising the deviations from the target landing time [Reference Samà, D’ariano, D’ariano and Pacciarelli51], minimising workload [Reference Furini, Kidd, Persiani and Toth20], minimising fuel consumption [Reference Dönmez, Çetek and Kaya6, Reference Jones, Lovell and Ball13, Reference Sama, D’ariano, D’ariano and Pacciarelli21] and minimising completion time [Reference Zhang, Zhao, Zhang, Dai and Sui5, Reference Dönmez, Çetek and Kaya6, Reference Ghoniem, Sherali and Baik10, Reference Sölveling and Clarke15, Reference Montoya, Rathinam and Wood16, Reference Sama, D’ariano, D’ariano and Pacciarelli21, Reference Sadovsky and Windhorst38, Reference Solveling, Solak, Clarke and Johnson46, Reference Gupta, Malik and Jung47, Reference Samà, D’ariano and Pacciarelli49, Reference Salehipour52–Reference Kaplan and Çetek55]. As this paper focuses on the uncertainties in TMAs and aircraft operations, the following studies are worth mentioning in detail as they consider stochastic ASSP. Sölveling et al. [Reference Solveling, Solak, Clarke and Johnson46] proposed two-stage stochastic programming (SP) for runway scheduling, including the uncertainties caused by pushback and taxiing movements. The results showed that the stochastic method improved the results of both deterministic and first-come first-served (FCFS) approaches. Sölveling and Clarke [Reference Cecen and Çetek14] developed a stochastic branch-and-bound algorithm to reach feasible solutions quickly. The outcomes indicated that the algorithm could decrease the makespan time by 7% compared to the deterministic approach. Choi et al. [Reference Choi, Robinson, Mulfinger and Capozzi56] developed a heuristic approach based on mixed-integer linear programming (MILP) for the aircraft scheduling problem using an efficient route structure in the extended TMAs compared with FCFS strategy solutions. Uncertainties were added to their model using extra separation buffers. Heidt et al. [Reference Heidt, Helmke, Kapolke, Liers and Martin57] proposed a stochastic mixed-integer programming approach for the runway scheduling model. They first presented a deterministic model then added uncertainties to the model. The results revealed that stochastic aircraft sequencing could generate a more stable traffic flow. Ng et al. [Reference Ng, Lee, Chan and Qin58] proposed a robust aircraft sequencing and scheduling model for parallel runway operations using the artificial bee colony algorithm. The min-max regret approach was taken into consideration for the generated scenarios. Their algorithm generated a feasible solution in 5 minutes and showed a better performance than other metaheuristic algorithms. Murça [Reference Murça59] presented a MILP model for metering aircraft departures considering uncertainty in the taxi-out process. The model reduced runway delay and increased the predictability of take-off time. Solak et al. [Reference Solak, Solveling, Clarke and Johnson60] presented two different stochastic programming models and used two different solution methodologies based on Lagrangian-Based scenarios and the sample average approximation (SAA). Their results indicated that both models could generate feasible and high-quality solutions. Liu et al. [Reference Liu, Liang, Zheng, Chu and Chu34] presented a two-stage SP model for ASSP using the uncertainties of arrival time with a single runway. The first stage in their approach searches for the aircraft weight class assignment, followed by the second stage, which aims to assign each flight a position in the sequence. Similarly, Khassiba et al. [Reference Khassiba, Bastin, Cafieri, Gendron and Mongeau61] presented a two-stage SP model for the ASSP. It used simulation-based validation, which showed that the proposed approach could reduce the potential conflicts by 70% while increasing the arrival rate by 7.7%. Their findings revealed that their model significantly reduced the need for a holding stack by a few hours through integrating upstream linear holding. Hong et al. [Reference Sadovsky and Windhorst38] developed a two-stage stochastic MILP model for the aircraft landing problem (ALP). They implemented the PSO algorithm because of the complexity of the model and utilised the point merge system (PMS) to manage arrival traffic. The model used the uncertain descent times from sequencing legs to the final approach point for each aircraft. Cecen et al. [Reference Cecen, Cetek and Kaya26] presented a stochastic MILP model for a single-runway airport. They used the SA algorithm to obtain an optimum or near optimum solution quickly. The model used the PMS to regulate arrival traffic and integrated the uncertainties of two different wind directions for the sequencing legs’ upper and lower flight levels. The results revealed that the stochastic approach reduced the total aircraft delay compared to the deterministic approach. According to the results, as the number of aircraft increased, the value of stochastic solution (VSS) value also increased. Ng et al. [Reference Ng, Chen and Lee54] presented a robust optimisation model for the terminal traffic flow problem. They considered the effects of wind direction and hazardous weather on the airspeed value as an uncertainty. They enhanced the Benders-Decomposition method and compared the results with the conventional Benders-Decomposition method. The results showed that the proposed method had increased the qualities of the solution. Increases in flight duration reduce airspace capacity, create longer flight times, greater fuel consumption, and changes in the arrival sequence. The main problems here are that the re-sequencing process generates extra delays for the aircraft set and increased fuel consumption. In addition, uncertain flight durations and delays affect the airport’s gate assignment plans (i.e. gate planning, taxi route assignment). Furthermore, the delays can also lead to some changes in cabin crew schedules, which considerably raises the economic costs of re-scheduling.

Optimising fuel consumption per aircraft by considering flight duration uncertainties can help airlines and governments decrease financial costs and environmental effects. Less fuel consumption reduces the amount of petroleum needed and gas emissions from the aircraft (i.e. carbon dioxide, carbon monoxide, nitrogen oxide, and hydrocarbon). As environmental pollution negatively affects human health and nature, it directly causes governments to spend more money, which can be reduced as fuel consumption is lowered. Ensuring resilient sequencing allows for lower fuel consumption which reduces the airline’s operating costs Unlike Hong et al. [Reference Sadovsky and Windhorst38] who used uncertain aircraft arrival times from sequencing legs to the merge point and Liu et al. [Reference Liu, Liang, Zheng, Chu and Chu34] who considered uncertain aircraft arrival times, we concentrated on the uncertain flight times for each aircraft using historical data.

Our contributions are presented as follows. The model is based on individual flight duration uncertainties for arrival operations, and we propose a MILP model using the SAA method to solve this problem. To the best of our knowledge, individual flight duration uncertainties and the relationship with fuel consumption in the TMA have not been investigated in the ASSP literature.

This study presents a stochastic MILP model for a single and mixed operations runway to minimise average fuel consumption per aircraft. As aircraft sequencing depends on many variables, we focused on the uncertain flight duration; therefore, these uncertainties affected the operation times directly. Fuel consumption was estimated using three different operations. The first one is fuel consumption in descent operations. The second is fuel consumption caused by vector manoeuvres to maintain safe separation among the aircraft set. Suitable vector manoeuvres absorb airborne delays yet leads to more economical costs. The last operation is fuel consumption caused by ground holding duration.

This paper is organised as follows: Problem description is explained in Section 3, and Section 4 presents the stochastic programming approach. Then, the mathematical model is defined in Section 5, and the TS algorithm is presented in Section 6. Next, Section 7 discusses the computational results, and lastly, the conclusions and further research are described in Section 8.

3.0 Problem description

3.1. Air traffic control

This paper studies the stochastic version of the aircraft sequencing problem for a single-runway airport. Aircraft flight durations are uncertain, and this situation directly affects the expected time of arrival/departure and the delay for each aircraft. This uncertainty causes a lack of punctuality, weakens efficiency, and increases the cost of flight operations. For these reasons, obtaining a resilient aircraft sequence can help to reduce the fuel consumption per aircraft in the TMAs. ATCOs aim to effectively regulate airborne and ground aircraft traffic to maintain conflict-free aircraft operations by giving appropriate instructions to each aircraft. These instructions include altitude changes, speed changes, and heading angle changes. Safe separations minima are assumed to be 5 nautical miles (NM) on the horizontal plane and 1000 feet on the vertical plane. Each aircraft has an expected arrival/departure time. This time determines the initial runway scheduling, which is defined as the FCFS sequence, the commonly preferred method for ATCOs.

The FCFS method can lead to undesirable airborne and ground delays because of the aircraft separation manoeuvres. ATCOs can re-sequence the aircraft using vectoring and speed-changing techniques for arrivals and ground holding for departures. The vectoring method allows aircraft to extend or reduce their flight path, and speed changing also controls flight time. ATCOs can modify the initial aircraft sequencing using these techniques to minimise aircraft delays as much as possible. Dear [Reference Dear62] presented the constrained position shifting (CPS) approach to ease the workload of ATCOs by allowing limited position changing ( $\rho $ ) in the initial aircraft sequencing. As well as this approach, Bennell et al. [Reference Bennell, Mesgarpour and Potts4] presented the concept of maximum time-shifting that limits the aircraft delay. In this way, an increase in flight time was avoided as the CPS approach may cause long delays. In this study, we utilised the maximum time-shifting constraint for departure aircraft. While the separation times between two consecutive aircraft are either departure (D) or arrival (A) operations, aircraft performance categories (Heavy, Medium) are used to estimate vortex separation minima that are given in Tables 1 and 2 according to the ICAO documentation.

Table 1.

The minimum safety separation standards defined in terms of distance [Reference Balakrishnan and Chandran39] (NM).

Table 2.

Time separations [Reference Balakrishnan and Chandran39], in seconds.

By considering the aircraft speed and leading-trailing sequence, the time separation is calculated for each combination and presented in Table 2.

3.2. Traffic analysis and assumptions

This study selected Istanbul Sabiha Gökçen International Airport (LTFJ) and its TMA for the mathematical model. The arrival route structure was employed at LTFJ. The airport has a single mixed operations runway, designated 06/24, and is one of the busiest single-runway airports globally, with 35.5 million passengers in 2019 [63]. August 2018, which is one of the most active months for operations so far, was evaluated to achieve a realistic traffic distribution. The analysis found that August 17, 2018, represented the peak traffic demand [Reference Kaplan and Cetek64]. The hourly traffic data for the day of peak demand is presented in Fig. 1.

Figure 1.

Hourly distribution of landing, take-off, and total operations on August 17, 2018 [Reference Kaplan and Cetek64].

According to the traffic data, more than 95% of the flight belonged to the medium aircraft category. As a result, in our study, we assumed that the heavy and medium aircraft distribution is 10% and 90%, respectively, to evaluate the effects of the performance category mix. At LTFJ, aircraft use seven different entry points into the TMA: IBADU, DRAMO, IZMAL, ETAMP, TOKER, GUMRU and GINLI (Fig. 2). In addition, in this study, we assumed that both landing and take-off operations take place only on runway 06.

Each arrival fight follows a pre-defined arrival path with the optimal fuel velocity. According to Fig. 2, seven different arrival paths exist, and it is assumed that arrival aircraft should perform a suitable vector mauver before TMA entrance to avoid congestion in the TMA. Extra fuel consumption caused by vector manoeuvres is calculated for each aircraft flight and added to our model. Generally, ATCOs prefer to implement vector manoeuvre during peak hours within the TMA. This situation pauses descending operation for a while and changes descending trajectories. However, aircraft neither change their arrival path nor their pre-defined descent angle during their descent operation in this study. In this way, we aim to analyse how effective the vector manoeuvre is for maintaining the required delay durations. Extended TMA operations allow ATCOs to regulate air traffic before entering conventional TMA borders. Therefore, aircraft can perform a suitable vector manoeuvre before entering the TMA, and traffic congestion can be avoided. Saez et al. [Reference Sáez, Prats, Polishchuk and Polishchuk48] stated that applying some manoeuvres during the en-route phase to alter the TMA entrance time can prevent many aircraft conflicts.

According to their entry point, aircraft fly to one merge point into the TMA. The distance between merge points and the final approach point (FAP) is 6 NM, and the length from FAP to the runway is 10 NM. Westbound and eastbound traffic are separated to prevent high demand at the merge points. In addition, radar separation minima was set as 3 NM within the TMA.

Each departure has an operation time window that limits the latest operation time. With this approach, significant delays are avoided for departures because fuel consumption is lower during ground holding. On-time performance is a commonly accepted indicator for understanding punctuality for transportation. In the airline industry, an airline departure or arrival is considered to be on time if it occurs within 15 minutes of the scheduled time [65]. As this study aims to minimise fuel consumption, the model reduces airborne delays rather than departure delays. We added the maximum departure delay to the model to find the balance between arrival and departure operations and calculated fuel consumption using the ground holding duration for each departure aircraft.

Flight duration is significant for the estimation time of arrival because this value may increase or decrease the total airborne delay according to the traffic situation. Therefore, aircraft arrival times are assumed to be stochastic due to these uncertain flight durations. In this model, flight durations are independent uncertainties. The expected departure time is also considered to be stochastic and determined using historical data for the departure operations. The deviation of the ground operations time for departure aircraft is calculated as a difference between the scheduled departure time and the actual departure time. Therefore, the uncertainty of ground operations was considered in the model. To obtain the distribution of each flight, we selected a specific period, between January 2017 and December 2019, for the chosen airport and analysed all flight durations. The aircraft set was determined using historical flight data, including several arrival and departure destinations. All aircraft flight duration data were taken from the Flight Radar 24 database. Since LTFJ serves as a regional hub airport, it presents many flight options. The arrival and departure aircraft set used in our analysis contained 29 different destinations, shown in Table 3. The number of arrival operations to LTJF was 27, and the number of departure operations from LTFJ airport was 19. The arrival set had 9 different international and 18 different domestic flight destinations. Similarly, the departure set included 7 international and 12 domestic flight destinations. The mean and standard deviations of each operation duration were estimated.

Table 3.

The LTFJ airport destinations used in this study.

Figure 2.

The route structure of LTFJ within the Istanbul TMA used in this study.

To illustrate the flight duration distribution, two different routes are presented in Figs. 3 and 4.

Figure 3.

The histogram of flight duration for Samsun-Istanbul flights.

Figure 4.

The histogram of flight duration for Ordu-Istanbul flights.

The mean and standard deviation values 66.9 minutes and 5.1 for Fig. 3 and 75.4 minutes and 5.48 for Fig. 4. In this study, the expected time of departure/arrival for all aircraft is obtained using the exponential distribution. For the arrival operation, we calculated the expected TMA entrance time using the distance from the entry point to the runway and average airspeed according to the entry point for each arrival aircraft. Furthermore, the expected TMA entry times were added to the model as a parameter. Descent operations were performed according to the BADA 3.11 [66]. The fuel consumption of descending with idle configurations, approach, and landing phases were also calculated and added to the model as a parameter. All calculations were implemented considering the entry point location because it affects the descent operation. Because aircraft neither change their nominal velocity and extend nor shorten the pre-defined route path during descend operation. When calculating the fuel consumption of the departure aircraft during ground holding, it was assumed that their engines operate in idle mode, and the unit fuel consumption value was taken from the ICAO database [67].

4.0 Stochastic programming

Stochastic programming (SP) considers the uncertainties [Reference Birge and Louveaux68, Reference Rockafellar and Wets69]. There are two types of information data for uncertain parameters: historical data and expert opinion. These uncertain parameters are represented using two different approaches, continuous probability distribution and the scenario-based stochastic approach. Each scenario has the same probability in this model. The general formulation of SP is given as follows, where ${\rm{f}}\!\left( {\rm{x}} \right)$ is the scenario-independent cost and ${\rm{Q}}\!\left( {{\rm{x}},{\rm{\;S}}} \right)$ is the cost when scenario $S$ is realised. The objective function (1) of the problem includes the value of f(x) and the expected value of ${\rm{Q}}\!\left( {{\rm{x}}, {{\rm{S}}^n}} \right)$ for scenario n. Some constraints of the model may be scenario independent, as given in Equation (2), or scenario dependent, as given in Equation (3).

(1) \begin{align}{\rm{Mi}}{{\rm{n}}_{\rm{x}}} ({\rm{f}}\!\left( {\rm{x}} \right) + \frac{1}{N}\mathop \sum \limits_{n = 1}^N \left[ {{\rm{Q}}\!\left( {{\rm{x}}, {{\rm{S}}^n}} \right)} \right]\end{align}

(2) \begin{align}{\rm{subject\;to\;g}}\!\left( {\rm{x}} \right) = {\rm{a}}\end{align}

(3) \begin{align}{\rm{h}}\!\left( {{\rm{x}},{\rm{S}}} \right) = {\rm{b}}\!\left( {\rm{S}} \right)\forall S\end{align}

This study applied the SAA technique for stochastic programming to solve many scenarios based on Monte Carlo simulations [Reference Lima, Conejo, Giraldi, le Maitre, Hoteit and Knio70]. The idea is to generate an M independent SAA problem, and the size of N. A large scenario size also increases the efficiency of the model. SAA also allows us to calculate the optimal gap. Using this value enables us to evaluate the quality of the solution. The detailed steps of the SAA technique are given as follows:

1. 1) For each independent SAA problem: 1, 2, 3, … … , M

   1. a) Generate N scenario

   2. b) Solve the corresponding problem

      1. $v_n^j = {\rm{Mi}}{{\rm{n}}_{\rm{x}}} ({\rm{f}}\!\left( {\rm{x}} \right) + \frac{1}{N}\mathop \sum \nolimits_{n = 1}^N \left[ {{\rm{Q}}\!\left( {{\rm{x}}, {{\rm{S}}^n}} \right)} \right]$ also $v_n^j$ can be defined as the optimum value for each independent SAA problem where j=1, 2, 3, …, M.

2. 2) Calculate the average of each independent SAA problem and the variance of this estimator

(4) \begin{align}\overline {{v_{N,M}}} = \frac{1}{M}\mathop \sum \limits_{j = 1}^M v_n^j\qquad\qquad\qquad\qquad\end{align}

(5) \begin{align}\sigma _{N,M}^2 = \frac{1}{{M\!\left( {M - 1} \right)}}\mathop \sum \limits_{j = 1}^M {\left( {v_n^j - \overline {{v_{N,M}}} } \right)^2}\end{align}

1. 3) Calculate the objective function for the original problem

   1. a) Generate N ^′ different scenario (N ^′ > N)

   2. b) Select a feasible solution $\left( {\hat x} \right)$ for the original problem and calculate the objective function for this solution using N^′ scenario

(6) \begin{align}{v_{N^{\prime}}} = {\rm{f}}\!\left( {\rm{x}} \right) + \frac{1}{{N^{\prime}}}\mathop \sum \limits_{n = 1}^{N^{\prime}} {\rm{Q}}\!\left( {\hat x, {{\rm{S}}^n}} \right)\end{align}

1. 4) Calculate the variance of this estimator

(7) \begin{align}\sigma _{{v_{N^{\prime}}}}^2 = \frac{1}{{N\!\left( {N^{\prime} - 1} \right)}}\mathop \sum \limits_{n = 1}^{N^{\prime}} {\left( {{\rm{f}}\!\left( {\rm{x}} \right) + {\rm{Q}}\!\left( {\hat x,{\rm{\;\;}}{{\rm{S}}^n}} \right) - {v_{N^{\prime}}}} \right)^2}\end{align}

1. 5) Calculate the lower bound (LB), upper bound (UB), and optimality gap of the solution and the variance of the gap

(8) \begin{align}UB = {v_{N^{\prime}}}\end{align}

(9) \begin{align}LB = \overline {{v_{N,M}}} \end{align}

(10) \begin{align}gap = \frac{{{v_{N^{\prime}}} - \overline {{v_{N,M}}} }}{{{v_{N^{\prime}}}}}\end{align}

(11) \begin{align}\sigma _{gap}^2 = \sigma _{{v_{N^{\prime}}}}^2 + \sigma _{N,M}^2\end{align}

To evaluate the proposed model, the expected value (EV), here-and-now solution (RP), and the value of the stochastic solution (VSS) were calculated. Firstly, we solved the deterministic ASSP assuming that the flight duration is known in advance. Therefore, we included the expected value of the flight durations as a parameter in the model. Secondly, we considered flight duration uncertainties while solving the model. Finally, various scenarios were generated for the stochastic ASSP model to obtain the optimum aircraft sequencing and calculate the objective function (RP). The difference between EV and RP indicates the advantage of the SP, which is referred to as VSS.

(12) \begin{align}VSS = EV - RP\end{align}

5.0 Mathematical model

The aim is to minimise the average fuel consumption per aircraft, and the entire mathematical model is presented as follows. It was also assumed that aircraft perform a suitable vector manoeuvre when the air traffic controllers give the related instructions for a specific altitude to maintain safe separation among arrival aircraft.

Indices:

1. I set of aircraft $\;\;i,\;{i_1}$ , ${i_2} \in I$

2. L set of scenarios $l, \in L$

Parameters:

${r_i}:$

The entry point into the TMA for aircraft i

${p_i}:$

The operation type of aircraft i( ${p_i} = $ 1 departure, ${p_i} = $ 2 arrival)

${c_i}:$

The performance category of aircraft i

${s_i}:$

The approach side of aircraft i for the merge points

${T_{{p_{{i_1}}}{p_{{i_2}}}{c_{{i_1}}}{c_{{i_2}}}}}:$

The separation time on the runway between operation types of ${p_{{i_1}}}$ and ${p_{{i_2}}}$ with aircraft performance categories of ${c_{{i_1}}}$ and ${c_{{i_2}}}$

${g_{i,l}}:$

The expected TMA entry time or departure time of aircraft i in scenario l

$f{h_i}$ :

The fuel consumption per second of aircraft i for vector manoeuvres

$f{l_i}$ :

The fuel consumption per second of aircraft i for ground holdings

${k_{{r_i}}}:\;$

The fuel consumption value for the descent operations of aircraft i, according to its entry point for arrival operations, takes 0 value for departure operations.

$B:\;$

The required time separation at the merge points

$G:\;$

The maximum ground delay in scenario l for each departure aircraft i

H:

The travel time between merge points and runway for arrival aircraft

M1:

A positive number that is sufficiently high enough

O:

The number of aircraft

N:

The scenario size

Decision variables:

$w{h_{i,l}}:$

The delay duration for aircraft i in scenario l for arrival operations

$w{l_{i,l}}:$

The delay duration for aircraft i in scenario l for departure operations

${e_{i,l}}:$

The actual time of crossing the merge point of aircraft i in scenario l

${q_{i,l}}:$

The arrival/departure time of aircraft i in scenario l

${x_{{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 0

${d_{i,l}}$ :

The total delay for aircraft i in scenario l

(13) \begin{align}\min \frac{1}{O} \cdot \frac{1}{N}\;\mathop \sum \limits_i \mathop \sum \limits_l w{h_{i,l}} \cdot \;f{h_i} + w{l_{i,l}} \cdot \;f{l_i}\; + {k_{{r_i}}}\end{align}

Subject to

(14) \begin{align}{e_{i,l}} = {g_{i,l}} + w{h_{i,l}}\;\;\forall \;i,l\;\;|\;\;{p_i} = 2\end{align}

(15) \begin{align}{e_{{i_{1,l}}}} - {e_{{i_2},l}} \ge B - \left( {{x_{{i_1}{i_2}}}} \right)M1\;\;\;\;\;\;\;\;\forall \;{i_1},{i_2},l\;\;|\;{i_1} \ne {i_2},\;\;{s_{{i_1}}} = {s_{{i_2}}},\;{p_{{i_1}}} = 2,{p_{{i_2}}} = 2\end{align}

(16) \begin{align}{e_{{i_2},l}} - {e_{{i_1},l}} \ge B - \left( {1 - {x_{{i_1}{i_2}}}} \right)M1\;\;\;\;\;\;\;\;\;\;\;\forall \;{i_1},{i_2},l\;\;|\;{i_1} \ne {i_2},\;{s_{{i_1}}} = {s_{{i_2}}},\;{p_{{i_1}}} = 2,{p_{{i_2}}} = 2\end{align}

(17) \begin{align}{q_{i,l}} = {e_{i,l}} + H\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall \;i,l\;\;|\;\;{p_i} = 2\;\end{align}

(18) \begin{align}{q_{i,l}} = {g_{i,l}} + w{l_{i,l}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall \;i,l\;\;|\;\;{p_i} = 1\end{align}

(19) \begin{align}{q_{{i_1},l}} - {q_{{i_2},l}} \ge {T_{{p_{{i_2}}}{p_{{i_1}}}{c_{{i_2}}}{c_{{i_1}}}}} - \left( {{x_{{i_1}{i_2}}}} \right)M1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall \;{i_1},{i_2},l\;\;|\;{i_1} \ne {i_2}\;\end{align}

(20) \begin{align}{q_{{i_2},l}} - {q_{{i_1},l}} \ge {T_{{p_{{i_1}}}{p_{{i_2}}}{c_{{i_1}}}{c_{{i_2}}}}} - \left( {1 - {x_{{i_1}{i_2}}}} \right)M1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall \;\;{i_1},{i_2},l\;\;|\;{i_1} \ne {i_2}\end{align}

(21) \begin{align}{d_{i,l}} = \;{q_{i,l}} - {g_{i,l}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall \;\;i,l\end{align}

(22) \begin{align}w{l_{i,l}} \le G\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall \;\;i,l\;|\;\;{p_i} = 1\end{align}

(23) \begin{align}{e_{i,l}},{q_{i,l}},{d_{i,l}},w{h_{i,l}},w{l_{i,l}}\; \ge 0 \qquad \qquad \forall \;i \in I\;\forall \;l \in L\end{align}

(24) \begin{align}{x_{{i_1},{i_2}}} \in \left\{ {0,1} \right\}\qquad\qquad\qquad \forall \;{i_1},{i_2}\; \in I\;\;\;\end{align}

The objective function given in Equation (13) is to minimise the average fuel consumption per aircraft. Constraint (14) estimates the crossing time for arrival aircraft at merge points for each scenario. Moreover, Constraints (15) and (16) maintain collision avoidance at the merge points if aircraft ${i_1}$ and aircraft ${i_2}$ use the same merge point for the same scenario. Constraints (17) and (18) calculate the runway use time for arrival and departure operations, respectively, for each scenario. Constraints (19) and (20) check the vortex separations minima on the runway for each aircraft pair for the same scenario. Constraint (21) calculates the delay of each aircraft for each scenario. Constraint (22) limits the maximum ground delays. Constraints (23–24) are sign constraints.

6.0 Tabu search algorithm

Glover presented the Tabu Search (TS) algorithm in 1990 [Reference Glover71]. It can be adapted for various problems (i.e. routing problem [Reference Alinaghian, Tirkolaee, Dezaki, Hejazi and Ding72], assignment problem [Reference Zhang, Liu, Zhou and Zhang73], job scheduling problem [Reference Vela, Afsar, Palacios, González-Rodríguez and Puente74], hub location problem [Reference Sangaiah and Khanduzi75]) because the algorithm can be implemented easily. The TS algorithm uses the initial solution, shown in Fig. 5, and it was applied in our study to obtain optimum or near optimum aircraft sequencing in a short computational time. Figure 5 shows that aircraft1 is assigned to the first position in the sequence and aircraft 2 is assigned to the third position.

Figure 5.

Initial aircraft sequencing representation.

Additionally, the TS algorithm uses an interchange operation to produce neighbourhood solutions which are shown in Fig. 6.

Figure 6.

Generation of neighbourhood aircraft sequencing.

Figure 6 shows that aircraft 2 and 3 are in a different position to generate other solutions for the first neighbourhood solution. Similarly, aircraft 1 and aircraft 3 were swapped to create another solution. Any number of neighbourhood solutions can be generated using this technique. The algorithm also uses the memory structure to avoid repetitions and discover another part of the solution space for a pre-defined number of iterations. The algorithm calculates each solution’s objective functions for each iteration, chooses the best solution, and puts the selected changes in the tabu list. After a certain number of iterations, aircraft are removed from the tabu list, and tabu list size can be determined according to the problem. The main objective of the tabu list is to ensure that the algorithm searches for the optimum value by avoiding repetitions. Suppose any solution in the tabu list generates a better solution than the best-found solution so far. In that case, the algorithm allows the selection of this solution to improve the quality of the solutions. This operation is called the aspiration criterion.

7.0 Computational results

This section presented the test results that showed the effects of flight duration uncertainties. The stochastic MILP model was tested in both the GAMS/CPLEX solver and the TS algorithm. It is significant to analyse the performance of the TS algorithm; therefore, six independent traffic sets were generated, and the results of GAMS/CPLEX solver and TS algorithms were compared. This comparison was made to show whether TS algorithm can reach an optimum solution or not within a reasonable duration. The results are presented in Table 4. The first column shows the test problem numbers, the second and third columns show the objective function (z) and CPU time (sec) for the GAMS/CPLEX solver. Similarly, the fourth and fifth columns give the objective function (z) (kg) and CPU time (sec) for the TS algorithm.

Table 4.

The results of test problems.

Table 5.

The parameters of the TS algorithm.

As presented in Table 4, the TS algorithm reached optimum aircraft sequencing in test problems more quickly than the CPLEX solver. The parameters used in the TS algorithm are presented in Table 5.

A computer with a 2.3 GHz Intel Core i7 processor and 16GB of RAM was used in all the computations that were processed in MATLAB software. Exponential distribution was used to determine expected runway use time for each air traffic set because it can accommodate the time between consecutive events as time flows continuously. The mean of exponential distribution shows the operation rate per unit time; therefore, the hourly traffic was used to obtain the mean value. After finding the performance of the TS, 20 (M) sets of different air traffic were randomly generated with sizes of 50 (N) and 100 (N) scenarios by only changing the expected runway use time using the normal distribution for each air traffic set. Furthermore, uniform distribution was used to generate the operation type for each aircraft. In addition, the TMA entry point assignments of aircraft were determined using historical data, and each one has a specific probability. The test problem results are statistically presented for air traffic set (M) and with size (N) in Table 6. The first column demonstrates the scenario size. Columns two and three indicate the statistical upper bound (UB) (E) (kg.) and its standard deviation (Std). Columns four and five display lower bound (LB) (E) (kg.) and its standard deviation (Std). The last column shows the optimality gap (R) (in percentage).

Table 6.

UB and LB of the SAA method found with M = 20, N=50 and N=100, and N^′ = 400.

According to the results, the standard deviations decreased for the upper and lower bound when we increased the scenario number. Similarly, the optimality gap also declined as we increased the scenario number. As Ertem et al. [Reference Ertem, Ozcelik and Saraç76] indicated, increasing the scenarios maintains tighter confidence intervals on the LB value. In the same way, a higher number of scenarios provides tighter confidence intervals on the UB value. The EV (kg), RP (kg), and VSS (kg) of the test problems are given in Table 7, where the first column shows the air traffic set, and the second to fifth columns show the EV (kg), RP (kg), VSS (kg), and CPU (sec) times for the 36 aircraft for N=50. Similarly, the results of a N=100 are presented in the following columns.

Table 7.

Results of the TS algorithm for average fuel consumption (kg) per aircraft for 36 aircraft per hour.

The results show that the proposed model could find VSS for each air traffic set with the size of 50 and 100. The enhancement rates were 9.11% and 8.78% for the N=50 and N=100, respectively. Also, the minimum VSS results are 7.2 kg and 4 kg for the N=50 and N=100. The maximum VSS value was found to be 67.4 kg in the N=50, whereas this was 57.7 kg for the N=100. These results indicate that as the scenario size increases, the difference between EV and RP decreases. In addition to this, the average CPU times were 33.2 sec and 107.8 sec for the N=50 and N=100, and the maximum CPU times values were 73.4 and 176.1 sec. This shows that our proposed model could obtain a feasible solution for the N=50 and N=100 in an acceptable amount of time. VSS demonstrated the benefits of SP because determining aircraft sequencing with stochastic situations provides more resilient traffic sequencing than the deterministic approach. In addition, the outputs also showed that the stochastic model reduced the average delay per aircraft significantly. This decreased from 285.6 sec to 230.7 sec for the N=50 and 289.8 sec to 235 sec for the N=100. The minimum delay per aircraft was 19.9 sec and 16.9 sec for the N=50 and N=100. Additionally, the maximum delay per aircraft value was 121.3 sec in the N=50 and 105.6 sec in the N=100. Furthermore, the average number of position-shifted aircraft, when compared to the initial sequence, was found to be 19.8 and 20, respectively. The model found the same aircraft sequencing in the six air traffic sets for both N=50 and N=100.