A novel stochastic matheuristic approach for the dynamic utilisation of arrival route structures: modelling of Istanbul Airport point merge system

Abstract

Congested airports benefit from parallel-point merge systems (P-PMSs) for efficient arrival route control. However, the decline in air traffic due to COVID-19 has curtailed its optimal utilisation, especially with the reduced need for long sequencing legs. As air traffic is poised to rebound, the evident volatility seen during and post COVID-19, as well as the daily fluctuations between peak and off-peak hours, underscore the importance of the dynamic utilisation of sequencing legs in P-PMSs. EUROCONTROL proposes various leg configurations to manage fluctuating traffics, ensuring both efficiency and safety. First, we proposed two additional leg configurations for the Istanbul Airport, offering continuous descent with the engines operating at idle thrust during leg flights; partially overlapped and fully dissociated. While they offer an alternative for controllers during low to medium traffic scenarios, current long and fully overlapped parallel legs may be used in high traffic due to the volatility of traffic density throughout a day. Therefore, we suggest an approach that provides dynamic utilisation of these configurations. We first modeled and analysed the configurations for various traffic numbers and scenarios. Then, we introduced a new stochastic matheuristic model that considers the configuration changes throughout the day and provides feasible and robust sequences applicable to all configurations by combining the benefits of mathematical models with the adaptability and speed of heuristic methods. Several test problems were evaluated using the terminal manoeuvering area structure of Istanbul Airport as a case study. The results indicate that by changing configurations, an average of 35 kg in fuel savings per aircraft can be achieved. The results also show that the proposed approach outperforms traditional stochastic mathematical models and the first-come first-serve (FCFS) strategy, ensuring efficient air traffic management in terms of fuel and delay with robust sequencing by eliminating the need for re-sequencing during configuration changes.

Nomenclature

ATC

Air Traffic Control

ATCo

air traffic controller

AST

average solution time

ETA

expected time of arrival

FCFS

first come first serve

GAMS

General Algebraic Modeling System

MM

mathematical model

MH

matheuristic model

PMS

point merge system

P-PMS

parallel point merge system

RNAV

area navigation

SA

simulated annealing

SL

sequencing leg

MP

merge point

TMA

terminal manoeuvering area

LTFM

İstanbul Airport

VES

value of expected solution

1.0 Introduction

The point merge system (PMS) was developed by EUROCONTROL in 2006 as an alternative to the traditional vector technique for terminal manoeuvering areas (TMAs) to allow air traffic controllers (ATCos) to implement systematic sequences even during peak hours [1]. The PMS involves directing an aircraft to sequencing legs (SLs) and then to the merge point (MP) with a single instruction after providing the necessary separation. Range rings are used in the system to help ATCos understand that the required separation has been achieved between successive aircraft pairs. PMS utilises the area navigation (RNAV) route structure, allowing for continuous descent operations during the flight between legs and merge point and significantly reducing controller workload and fuel savings [2]. A representation of a traditional single PMS is given in Fig. 1.

Figure 1. Traditional single PMS [1].

Several studies have indicated that the PMS offers several advantages compared to the traditional vector technique. These advantages include reducing the number of instructions given to the pilot, fuel consumption, controller workload, controller-pilot bi-directional frequency occupancy time, risk of communication errors, adverse environmental effects, flight time and flight distance. The PMS also enables predictions regarding landing sequences, efficient use of runway capacity, and collecting arrival routes into a narrower area. Through continuous descent approaches, the PMS system enhances situational awareness using range rings and enables optimised noise level control over densely populated areas. However, it is important to note that while PMS facilitates continuous descent operations, its fixed intercept of final approach may result in concentrated air traffic over-populated areas. Additionally, the PMS increases safety and separation accuracy, while initial training for controllers and pilots is straightforward [3–6].

However, the COVID-19 pandemic and the resulting significant reduction in air traffic have hindered the system’s full potential. It is observed that the long sequencing legs have become unnecessary for the limited traffic [7]. Nevertheless, with the expected increase in traffic volumes, using these legs will likely be reinstated. Even then, their use may differ during peak and off-peak hours. Long parallel legs may be necessary during peak hours, while separated short legs could be more useful during less busy hours.

The positioning of sequencing legs in air traffic management can significantly impact vertical performance. EUROCONTROL proposed three options for positioning sequencing legs shown in Fig. 2, fully overlapping, fully dissociated or partially overlapping. Overlapping legs provide larger delay absorption capacity but require vertical separation during level-off segments, while dissociated legs offer improved vertical profiles but use more airspace horizontally [1].

Figure 2. Different sequencing legs configurations of PMS [1].

Shorter, dissociated legs may be used during off-peak hours to facilitate continuous descents, while fully overlapping legs can be used during peak hours to increase capacity. Although the capacity of system increases when using the overlapping longer legs, specific safety requirements must be maintained including published vertical restrictions, adherence monitoring and specific controller training. The inner leg must also be higher than the outer leg to avoid separation losses. Therefore, the system’s efficiency can be improved by changing the configuration throughout the day [1]. However, this change may alter the sequencing making considering the previous configuration. Furthermore, with the increased usage of parallel point merge systems (P-PMS), which includes the combination of two PMS in busy airports, the dynamic use of sequencing legs may be more complex for the air traffic controllers at these airports.

To address the issue of sequencing that may become invalid in the face of configuration changes, especially in P-PMSs, a stochastic matheuristic approach is proposed as a more efficient alternative to classical stochastic model that produce sequencing for a single configuration especially for high traffic demands.

The TMA of LTFM including P-PMS is considered as a case study. Optimisation of the arrival sequences in the TMA including multiple runway assignments is handled which is defined as the runway scheduling problem (RSP) in the literature. Harikiopoulo D. proved in that the multiple runway aircraft scheduling problem is not NP-complete [8]. However, the run time can be computationally prohibitive as the number of classes or number of runways increases [9]. Therefore, a matheuristic approach is integrated to solve the problem.

Matheuristics have raised quite new techniques in terms of operational research and provide an advantage compared to exact methods or heuristics. It uses the benefit of both approaches by finding optimal results in reasonable times, even for complex problems. Matheuristics are optimisation algorithms that utilise mathematical programming techniques to generate heuristic solutions, independent of specific problem domains. Unlike traditional methods, matheuristics integrate problem-dependent elements mainly within lower-level mathematical programming, local search or constructive components. Matheuristics provide adaptable and efficient problem-solving strategies by integrating mathematical insights into heuristic frameworks, making them suitable for various optimisation challenges [10, 11].

The stochastic programming approach is extended to the stochastic-matheuristic method to provide robust solutions for dynamically changing routes even with high traffic. Thanks to this approach, it will be possible to make numerous configuration changes in the legs of P-PMSs throughout the day, even during dense operations, without requiring any sequencing changes. This approach may enable the full utilisation of the benefits of configuration changes while relieving the controllers from a significant workload during such transitions.

1.1 State of art

PMS is discussed in many studies testing its potential in simulation environments using simple or complex scenarios. After the simulations, its applications to several airports, including Oslo [12], Dublin [13] and other airports in four different continents [1]. Then the system has become more common in Europe and the world. Increased traffic revealed the necessity of optimal usage of the system. Therefore, several studies handled the PMS optimisation, including various objective functions such as minimising fuel [14], delay [15], flight time [16], makespan [17], etc. Some of the optimisation studies used exact algorithms, including mixed integer programming [18], multi-objective programming [2, 19] and stochastic programming [20, 21]. Although these approaches provided optimal results and sequences regarding objectives, the number of aircraft handled in these algorithms was relatively small. To overcome this issue, heuristics approaches are provided by testing its meta-heuristic algorithms, including SA [19, 22–24], genetic algorithm [25] and particle swarm optimisation [26]. However, these approaches only provide near-optimal solutions for the problem. On the other hand, most of the studies handled a single PMS. However, increasing traffic numbers and the increasing usage of parallel runways required the P-PMSs. A few studies considered P-PMS optimisation in the literature [2]. They provided a deterministic multi-objective mathematical model for P-PMS considering the objectives of delay, fuel and flight time minimisations. They have further increased the previously proposed model to a stochastic multi-objective model by considering the system’s wind speed and direction uncertainty [21]. However, these studies did not consider the dynamic change of the leg use in the system. Today’s demand volatility requires the active use of the PMS.

One of the important study presented by Ref. (27) examines the opportunities for introducing dynamic arrival route structures in terminal areas, utilising recorded data from 50,000 flights across four busy European airports. It employs statistical clustering techniques to identify and characterise different operational regimes, highlighting the potential for implementing varied route structures in terminal areas with sufficient traffic level variability, considering operational feasibility constraints. This provides a foundational basis for developing arrival management systems that dynamically respond to traffic density, enhancing operational flexibility.

In this current study, a stochastic matheuristic model is developed for P-PMS with various leg configurations. The matheuristic approach provides feasible solutions within a reasonable time horizon, even in dense traffic scenarios. The originality of the current study can be summarised as follow.

• • Handles the P-PMS by considering increased air traffic demand

• • Suggests new designs for PMS legs in airports

• • Considers the dynamic usage of sequencing legs of the system

• • Develops a new stochastic-matheuristic approach which is a brand-new approach in terms of operational research

• • Provides feasible solutions in reasonable times, even for dense traffic numbers

• • Provides robust sequences which are feasible for all the scenarios

• • Reduces the controller workload by eliminating the necessity of re-sequencing even in configuration changes during operations

• • Reduces the delays compared to the FCFS and deterministic approach

2.0 Methodology

2.1 Case study

Istanbul airport which has four PMS used as a case study in this work. Two PMSs are positioned in the north while others are south of TMA. Only north-side operations were handled in this study, meaning 16R, 17L and 18 runways served the arrival traffic. The P-PMS is modeled based on AIP Türkiye [28]. The entry points, TMA routes, legs, merge points and common points are represented in Fig. 3.

Figure 3. The current structure of P-PMS for Istanbul Airport (fully overlapping configuration) [1].

The P-PMS of Istanbul Airport is shown in Fig. 3 and has seven entry points (DRAMO, RILEX, RIXEN, ATPIX, ERSEN, SISPI, and INBET), four sequencing legs, two merge points (INSTA and DIVDI), three common points and three parallel runways for arrivals. In Fig. 3, the current operation layout is presented, which indicates a fully overlapping configuration. This configuration allows for greater path stretching capacity but causes level flight from the beginning of the sequencing legs to the end of the legs. Additionally, a 1000-ft vertical separation is required to prevent conflicts between aircraft. There are 2 nm distance between inner and outer legs of both PMS. In this study, several leg configurations are considered as we conceptualise alternative possibilities for the existing situation given in Figs 4 and 5.

Figure 4. Partially overlapping leg configuration design for Istanbul Airport.

Figure 5. Fully dissociated leg configuration design for Istanbul Airport.

Figure 4 represents the partially overlapping configuration, which provides an idle descent from the beginning of the leg to the overlapping waypoints. This offers vertical profile efficiency compared to the previous configuration; however, it presents less path stretching capacity and still requires vertically separated operations after the overlapping waypoints. Figure 5 indicates a full dissociation configuration, which provides continuous descent and does not require level flight, offering significant efficiency. However, this last configuration provides the least path stretching capacity for aircraft. This last configuration also requires a minimum horizontal separation between the last way points of the legs. In our case we assume that there were approximately 7 nm horizontal separation.

In the first scenario (fully overlap), the leg lengths were 30.8 and 35.5 nm for PMS-1 and 38.5 and 28.4 nm for PMS-2. In the second scenario (partially overlap), these lengths were 23.1 and 21.6 nm for both PMS-1 and PMS-2. And for the last scenario (fully dissociated), the leg length was shorter, at 15.4 and 14.2 nm for both PMSs.

To compare these configurations in terms of delay, fuel consumption and flight distance, we first modeled and analysed the configurations for various traffic numbers and scenarios. Then, we introduced a stochastic model that considers the configuration changes throughout the day and provides feasible and robust sequences applicable to all configurations. This stochastic model was then enhanced to a stochastic matheuristic model to improve solution time performance while still providing robust solutions.

2.2 Mathematical model

2.2.1 Stochastic programming

Uncertainties are addressed in stochastic programming (SP) through two primary avenues: expert opinion and historical data [29, 30]. To represent these uncertain parameters, two distinct methods are employed: the scenario-based stochastic approach and the continuous probability distribution method. The comprehensive depiction of SP is provided in the following, where ${\rm{f}}\left( {\rm{x}} \right){\rm{\;}}$ denotes the scenario-independent cost and ${\rm{Q}}\left( {{\rm{x}},{\rm{\;S}}} \right)$ represents the cost associated with scenario S. The objective function (1) encompasses the value of ${\rm{f}}\left( {\rm{x}} \right){\rm{\;}}$ alongside the expected value of ${\rm{Q}}\left( {{\rm{x}},{\rm{\;}}{{\rm{S}}^l}} \right)$ for scenario l. Multiple constraints within the model may either be scenario-independent, as expressed in Equation (2), or scenario-dependent, as delineated in Equation (3).

(1) \begin{equation}{\rm{Mi}}{{\rm{n}}_{\rm{x}}}{\rm{\;}}\bigg({\rm{f}}\left( {\rm{x}} \right) + \frac{1}{L}\mathop \sum \limits_{l = 1}^L \left[ {{\rm{Q}}\left( {{\rm{x}},{\rm{\;}}{{\rm{S}}^l}} \right)} \right]\end{equation}

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

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

We employed the scenario-based approach to incorporate leg usage configurations as scenarios into the model. The flight time limits for the legs are listed in Table 1 and span three different scenarios to represent the three different sequencing leg (SLs) setups for each PMS. We converted the lengths of legs from distances to times to align with our time-based approach used in the mathematical models because aircraft maintain a constant speed during flying the sequencing legs. A more accurate and efficient analysis is made possible by this transition into a temporal dimension. Although there are differences between the inner and outer legs’ length in each scenario, we chose to neglect these differences when converting these lengths to time based approach to simplify the modeling. This significantly reduces the number of scenarios, which notably affects the solution-finding capability of the stochastic model. Scenarios 1, 2 and 3 given in Table 1, illustrate the flight time limits at legs under the current configuration, where the legs are fully overlapped in both PMSs, partially overlapped, and fully dissociated, respectively.

Table 1. SLs configuration scenarios in seconds

A stochastic mathematical model based on Ref. (2) is initially modelled to address the RSP. The problem is solved using the exact solution method utilising the CPLEX solver in the general algebraic modeling system (GAMS). GAMS is a program for mathematical simulating and resolving mixed-integer, nonlinear and linear optimisation problems. The sets, parameters, variables and constraints are as follows.

Sets

i and p

Describe the set of aircraft

$n$

Describes the set of merge point of PMSes

s

Describes the set of scenarios

$t$

Describes the set of common points

${0_1}\;and\;{0_2}$

Describe the set of performance category

Parameters

${r_i}$

The assigned sequencing leg of i^th aircraft

$p{m_i}$

The assigned PMS of i^th aircraft

$typ{e_i}$

Type of i^th aircraft {small, medium, heavy}

$sece{p_i}$

Sector entry point of i^th aircraft

$s{e_i}$

Sector entry time of i^th aircraft

$c{t_i}$

Flight time of i^th aircraft from the sector entrance to the merge point

$tcd{a_i}$

CDA time from sequencing leg to merge point of the i^th aircraft

$le{g_{s,n}}$

The leg length of n^th PMS in scenario s

$m{f_{n,t}}$

The flight time from n^th merge point to t^th common point

$se{p_{{O_2},{O_1}}}$

Wake turbulence separation parameter between consecutive aircraft types

$M$

Large enough number

$\varepsilon $

Small enough number

$rada{r^1}$

The radar separation at PMS entry points

$rada{r^2}$

The radar separation at merge points

${q_s}$

Probability of s^th scenario

$pos_i^{eta}$

Represent the position of ${i^{th}}$ aircraft considering the expected arrival times – FCFS order

$MPS$

Maximum allowed position shifting

Variables

$z$

Objective function

${d_{i,s}}$

The PMS entry time of i^th aircraft in the scenario s

$penalty\_dela{y_{i,s}}$

Penalty delay at sector entry of i^th aircraft in scenario s

$arrdela{y_{i,s}}$

Total delay of i^th aircraft in scenario s

$m{p_{i,\;s}}$

Being time at merge point of i^th aircraft in s^th scenario

$faf{t_{i,s}}$

Being time at the common point of i^th aircraft in s^th scenario

$totdela{y_s}$

Total delay of the scenario s

${a_{i,s}}$

Flight time of i^th aircraft on sequencing leg in s^th scenario

${e_{i,p}}$

If i^th aircraft is assigned to leg before p^th aircraft takes 1, otherwise 0

${b_{i,p}}$

If i^th aircraft is assigned to merge point before p^th aircraft takes 1, otherwise 0

${x_{it}}$

If i^th aircraft is assigned to t^th common point takes 1, otherwise 0

${v_{i,p}}$

If i^th aircraft is assigned to common point before p^th aircraft 1, otherwise 0

$pos_i^{final}$

represent the position of ${i^{th}}$ aircraft at the common points

Constraints

(4) \begin{equation}{d_{i,s}} = s{e_i} + c{t_i} + dela{y_{i,s}}\qquad\qquad\qquad\qquad\forall \;i,\;s\end{equation}

(5) \begin{equation}penalty\_dela{y_{i,s}} \le 300\qquad\qquad\qquad\forall \;i,\;s\end{equation}

(6) \begin{equation}{d_{i,s}} - {d_{p,s}}\; \geq rada{r^1} - M\cdot{e_{i,p}}\qquad\qquad\forall \;i,p,\;s,\;i \ne p,\;{r_i} = {r_p}\end{equation}

(7) \begin{equation}{d_{p,s}} - {d_{i,s}} \geq rada{r^1} - M\cdot\left( {1 - {e_{i,p}}} \right)\qquad\qquad\forall \;i,p,\;s,\;i \ne p,\;{r_i} = {r_p}\end{equation}

(8) \begin{equation}m{p_{i,s}} = {d_{i,s}} + {a_{i,s}} + tcd{a_i}\qquad\qquad\forall \;i,\;s\end{equation}

(9) \begin{equation}{a_{i,s}} \le le{g_{s,\;n}}\qquad\qquad\forall \;i,\;s,\;n,\;\;n = P{M_i}\end{equation}

(10) \begin{equation}m{p_{i,s}} - m{p_{p,s}} \geq rada{r^2} - M\cdot{b_{i,p}}\qquad\qquad\forall \;i,p,\;s,\;i \ne p,\;p{m_i} = p{m_p}\end{equation}

(11) \begin{equation}m{p_{p,s}} - m{p_{i,s}} \geq rada{r^2} - M\cdot\left( {1 - {b_{i,p}}} \right)\qquad\qquad\forall \;i,p,\;s,\;i \ne p,\;p{m_i} = p{m_p}\end{equation}

(12) \begin{equation}\mathop \sum \limits_t {x_{it}} = 1\qquad\qquad\forall \;i\end{equation}

(13) \begin{equation}faf{t_{i,s}} = \mathop \sum \limits_t m{f_{n,t}}\cdot{x_{it}}\, + m{p_{i,\;s}}\qquad\qquad\forall \;i,n,\;s,\,n = p{m_i},\;\end{equation}

(14) \begin{align}faf{t_{i,s}} - faf{t_{p,s}} \geq se{p_{O,\;{0_2}}} - M\cdot{V_{i,p}} - M\cdot\left( {2 - {x_{i,t}} - {x_{p,t}}} \right) &\qquad\qquad \forall i,p,t,{0_1},{0_2},s,\,\,i \ne p,\nonumber\\& o1 = type_i,\;o2 = typ{e_p}\end{align}

(15) \begin{align}faf{t_{p,s}} - faf{t_{i,s}} \geq se{p_{O,\;{0_2}}} - M\cdot(1 - {v_{i,p}}) - M\cdot\left( {2 - {x_{i,t}} - {x_{p,t}}} \right) & \qquad\qquad\forall i,p,t,o1,o2,s,\;\;i \ne p,\nonumber\\& o1 = typ{e_i},\;\;o2 = typ{e_p}\end{align}

(16) \begin{equation}pos_i^{final} = n - \;\mathop \sum \limits_p^I {v_{i,p}}\qquad\qquad \forall i,p\;\;\end{equation}

(17) \begin{equation}pos_i^{final} - pos_i^{eta} \le MPS\;\;\qquad\qquad\forall i\;\end{equation}

(18) \begin{equation}pos_i^{eta} - pos_i^{final} \le MPS\;\;\qquad\qquad\forall i\;\end{equation}

(19) \begin{equation}arrdela{y_{i,s}} = {a_{i,s}} + penalty\_dela{y_{i,s}}\qquad\qquad\forall \;i,s\end{equation}

(20) \begin{equation}totdel{a_s} = \mathop \sum \limits_i arrdela{y_{i,s}}\qquad\qquad\forall s\end{equation}

(21) \begin{equation}z = \mathop \sum \limits_s totdela{y_s}\cdot{q_s}\end{equation}

Equation (4) is used to integrate the sector entry times, sector entry delays, and flight times up to the PMS entrances to determine the entry time for the i^th aircraft to the PMS for each scenario. The maximum of penalty delay at sector entries allowed by Equation (5) is 300 seconds. The penalty delay is used only if it resolves conflicts at intersection points or if the PMS capacity is insufficient for a feasible solution to the problem. Equations (6) and (7) are crucial in preserving the necessary separation between sequential aircraft at the PMS entries. Equation (8) allows to estimate how long it will take each aircraft to get to its specific merge point. Equation (9) guarantees that the flight times within the legs remain within the limits prescribed forth for each scenario. Equations (10) and (11) work together to prevent potential air traffic conflicts at the merge points by maintaining a safe separation between sequential aircraft that use the same PMS. Radar separations, on the other hand, are applied at PMS entrance and PMS merge points, considering the true airspeed of the aircraft.

Equation (12) ensures that each aircraft is assigned to one of the common points, which also serve as initial approach fixes for each parallel runway. The time needed for the i^th aircraft to arrive at the t^th common point is precisely calculated by Equation (13). Equations (14) and (15) guarantee the required spacing between successive aircraft at the common points, ensuring the preservation of the wake turbulence safety margins. The time-based wake turbulence separations are set according to the ICAO minimum wake vortex separation (DWVS) standard, as defined by air traffic management - ICAO DOC 4444, at common points. Table 2 provides the separation minimums at common points [31].

Table 2. Time based wake turbulence separations (sec.)

Equation (16) calculates the final position of aircraft. Equations (17) and (18) are designed to limit the maximum number of position shifts to three. For instance, if an aircraft starts at the ${10^{th}}$ position, according to these constraints, its final position could be ${7^{th}},\;{\rm{\;}}{8^{th}},{\rm{\;\;}}{9^{th}},{\rm{\;}}{10^{th}},{\rm{\;\;}}{11^{th}},{\rm{\;\;}}{12^{th}}{\rm{\;or\;}}{13^{th}}$ . The introduction of constrained position shifting (CPS) constraints to the model aims to mitigate the occurrence of excessive position changes, which could lead to unfair outcomes [32]. Equation (19), which accounts for both the sector entry delays and the time spent in the legs, calculates the total delay time for each aircraft in each scenario. The total delay for each scenario is calculated with the aid of Equation (20). The objective function, as expressed in Equation (21), aims to reduce the overall delay while taking the probability of each possible scenario into account.

2.2.2 Fuel calculations

It is imperative to determine the thrust-specific fuel consumption $\eta $ ( $\left( {kg/min \cdot kN} \right)$ ) to accurately estimate aircraft fuel consumption. This can be achieved through the following formula in Equation (22):

(22) \begin{equation}\eta = {C_{_1}} \times \left( {1 + \frac{{{V_{TAS}}}}{{{C_{{f_2}}}}}} \right)\end{equation}

Here, ${C_{{f_1}}}$ , ${C_{{f_2}}}$ , ${C_{{f_3}}}$ , and ${C_{{f_4}}}$ represent aircraft-specific coefficients, while ${V_{TAS}}$ denotes the true airspeed at pressure altitude ( ${H_p}$ ). Subsequently, the nominal fuel flow, ${f_{nom}}$ [kg/min] is computed as follows in Equation (23):

(23) \begin{equation}{f_{nom}} = \eta \times Thr\end{equation}

Note that $Thr$ is the value of thrust. Moreover, the cruise fuel flow, ${f_{cr}}$ [kg/min], is determined by incorporating, $Thr$ , and ${C_{fcr}}$ in Equation (24):

(24) \begin{equation}{f_{cr}} = {C_{fcr}} \times \eta \times Thr\end{equation}

In the optimisation model, which involves constant altitude flight operations, the required thrust is assumed to equal the drag D during level flight. This dependency is governed by the drag coefficient, ${C_D}$ , air density $\rho $ , wing reference area S, and true airspeed ${V_{TAS}}$ , according to Equation (25):

(25) \begin{equation}D = \frac{{{C_D} \cdot \rho \cdot {V_{TAS}}^2 \cdot S}}{2}\end{equation}

During descent operations, aircraft engines operate in idle configurations above the transition altitude. The fuel flow $\left( {{f_{min}}} \right){\rm{\;}}$ at the transition altitude is estimated using the formula in Equation (26):

(26) \begin{equation}{f_{min}} = {C_{f3}}\left( {1 - \frac{{{H_p}}}{{{C_{f4}}}}} \right)\end{equation}

Utilising the subsequent equation, the fuel consumption for the lower transition altitude during descent operations is computed in Equation (27):

(27) \begin{equation}{\rm{\;}}{f_{app/ld}} = MAX\left( {{f_{nom}},{f_{min}}} \right)\end{equation}

This methodology encompasses two performance categories: narrow-body and wide-body. Specifically, Airbus A321 and A333 aircraft performance values are selected to represent these categories.

2.3 Matheuristic approach

The exact solution methods may be unable to deliver a solution in a reasonable amount of time as the dataset size grows. On the other hand, using heuristics results with near optimal solutions. As a result, a matheuristic approach is presented in this study to solve the problem efficiently. Matheuristics are hybrid methods that combine mathematical models and heuristic approaches to solve complex problems. In matheuristic methods, the robustness of mathematical models is harmonised with the adaptive efficacy of heuristic methods. This combination enables the attainment of high-quality results in shorter and fixed times compared to exact methods. Scheduling and routing problems, being NP-Hard, typically require extended solution times to find optimal outcomes. Successful applications of matheuristic methods to solve complex problems are demonstrated in works by Refs (33–35). The matheuristic method hybridised with SA is employed to address the problem at hand. SA, a stochastic search method, takes inspiration from the annealing process in solids, where substances are heated to a high temperature and then gradually cooled. This concept was introduced by Ref. (36) and has since been efficaciously applied to a wide range of problem types. The flowchart of the proposed matheuristic method is shown in Fig. 6.

Figure 6. Flow chart of the mat-heuristic method.

SA is an optimisation approach which starts the search using an initial solution. The initial solution is generally determined using FCFS method depending on ETAs. The initial solution is selected as a current solution (X _now). Then, a neighbourhood solution (X _new) is obtained from current solution by swapping the sequence of aircraft. After that, the objective function values of these solutions are calculated in mathematical model. If the X _new has a lower value than X _now, X _new becomes current solution since our aim is to minimise the total delay. Otherwise, if X _new displays a weaker performance than the X _now, the new current solution is obtained utilising the probability function (P). In addition, a random value (a) is generated between 0 and 1. If a is lower than P value, X _new is chosen as a current solution. Otherwise, no change occurs in the current solution. Furthermore, the new current solution’s objective function value is compared to the best solution found so far. If X _now has a lower value, the algorithm updates the best solution and check whether the algorithm reach the stopping criteria or not. The proposed model continues the optimisation process until reaching the stopping criteria. The new temperature ( ${T_{new}}$ ) is calculated in Equation (28) using cooling rate ( $t$ ) and the initial temperature $\left( T \right)$ .

(28) \begin{equation}{T_{new}}\; = t \cdot T\end{equation}

This study structures the problem into two stages. In the first stage, the final position sequence of the aircraft is determined via SA. This schedule informs the mathematical model in the second stage, which calculates other variables and delays. Different sequencing is produced for each iteration, and based on these sequences, delay values and other variables are altered, leading to differing objective function values. The final position variable is defined as ${v_{i,p}}$ . This variable, obtained as per the result of SA, are used as parameters in the mathematical model, which then solves for the outcome. Aside from the parameterisation of variable ${v_{i,p}}$ , the mathematical model remains unchanged within the mat-heuristic algorithm and is used to solve for the optimal solution.

Figure 7. Neighbourhood structure.

In our study, the algorithm initiates with a FCFS approach and iteratively explores for better solutions using the SA heuristic technique. The algorithm typically yields a solution in approximately 900 seconds when the parameters we set, particularly the cooling rate, are applied. It’s crucial to note that we didn’t directly impose a 900-second time limit; instead, we set the parameters to generally lead to solutions within this timeframe. This observation guides us in controlling the algorithm’s solution process and determining an appropriate planning timeframe. The selection of approximately 900-second execution time as the stopping criterion is based on the considerations of air traffic management planning horizon, where a minimum preparation time of 15 minutes is deemed sufficient initially. Given the dynamic nature of air traffic and the need for timely decision-making, it’s essential to balance computational efficiency with solution quality. A 15-minute planning horizon provides controllers with adequate time to assess and implement sequencing decisions while ensuring operational flexibility.

In our approach, we initially attempted to find suitable solutions using mathematical models within GAMS. However, due to the complexity of the problem, no feasible solutions could be obtained. Consequently, the stochastic matheuristic approach was employed as an alternative, leveraging the SA heuristic to explore solution space iteratively. While matheuristic methods do not ensure optimality, they offer practical and effective solutions, aligning well with the operational demands of air traffic management. The Python Pyomo framework is utilised for modeling the matheuristic algorithm, complemented by the Gurobi solver for solving the samples.

We used a basic neighbourhood structure considering the CPS constraints of the model. Firstly, one aircraft is randomly selected in each iteration. Then, one of the six neighbourhood structures is selected depending on the generated random number. As shown in Fig. 7, a new ranking is obtained by changing the location of the selected aircraft between +3 and –3 according to the selected neighbourhood structure. For example, if the fifth aircraft is selected randomly and the second neighbourhood structure is randomly chosen, the selected aircraft moves to the front sequence.

3.0 Results

In this section, we analyse the samples representing small, medium and high traffic demand cases across three distinct leg configuration scenarios. Initially, we investigate demand samples involving 70 aircraft by resolving the mathematical model. Subsequently, we assess demand samples involving 90 aircraft, wherein we maintain a fixed sequence based on the FCFS order, as the model fails to yield optimal solutions for this level of traffic. Finally, we evaluate the efficacy of our proposed stochastic matheuristic algorithm in handling a demand of 100 aircraft within a one-hour timeframe.

3.1 Low traffic demand results (70 aircraft in one hour)

In this section, our primary focus revolved around examining each leg configuration scenarios across diverse traffic samples. We conducted a comparative analysis of these scenarios based on fuel consumption, delay and flight time metrics. Initially, we showcased the outcomes derived from the samples involving 70 aircraft, which corresponds to 35 aircraft demand per PMS. To ensure robustness, we solved 10 samples for 70 aircraft demand. Figure 1 illustrates the average total fuel consumption, average delay, and average flight time per aircraft in this samples across the three leg configuration scenarios. (It is noteworthy that all solutions were proven optimal up to the 70 aircraft scenarios.)

As depicted in Fig. 8, scenarios 2 and 3 (partially overlap and fully dissociated) notably reduced the average total fuel consumption compared to scenario 1 (fully overlap), showcasing approximately a 6.5% reduction in average fuel usage. On average, this translates to a fuel saving of 13.6 kg per aircraft, with some samples achieving savings of up to 21.4 kg. The rationale behind this reduction lies in the transition of leg flights from level-off to idle in scenarios 2 and 3, contrasting with scenario 1.

Figure 8. Seventy aircraft demand samples’ results for each sample in the three scenario.

In scenario 1, there were no idle flights, while in scenarios 2 and 3, there were no level-off flights in these samples. This deviation in flight behaviour stems from structural constraints in scenarios 1 and 3, whereas in scenario 2, the model exhibited a preference for avoiding level-off after idle segments, with all aircraft opting to turn before reaching this stage. Notably, the outcomes indicate similar results in terms of average delay and flight time across the scenarios.

For traffic scenarios of this level, choosing to use the fully dissociated leg configuration in the scenario 3 might prove most effective, as it eliminates the need for vertical separation between aircraft, potentially reducing the workload for air traffic controllers while yielding fuel savings. However, it’s crucial to acknowledge that the route structure in scenario 3 might not accommodate higher traffic volumes effectively. To investigate this further, we conducted simulations with 90 aircraft samples, detailed as follows.

3.2 Medium traffic demand results (90 aircraft in one hour)

To observe the disparities between leg scenarios, we conducted simulations for samples involving a demand of 90 aircraft. However, due to the model’s inability to find optimal solutions for these samples, we resolved the model by maintaining fixed sequences based on FCFS order. Consequently, we compared the leg scenarios under the FCFS sequence. The results from 10 samples are presented in Fig. 9.

As evident from the figures, the savings achieved in scenarios 2 and 3 surpassed those of scenario 1, particularly when compared to previous instances involving lower traffic samples. On average, the fuel saving per aircraft amounted to approximately 35 kg in scenarios 2 and 3 compared to scenario 1 with some samples achieving savings of up to 55 kg. However, in this medium traffic range, it’s worth noting that some aircraft experienced penalty delays at sector entry (SE Delay) in scenario 2 and 3, indicative of overcapacity conditions (Fig. 10).

Figure 9. Ninty aircraft demand samples’ results for each sample in the three scenarios.

Figure 10. Delay distributions in 90 AC demand samples.

As depicted in Fig. 10, penalty delays at sector entry were observed in some samples within scenarios 2 and 3, suggesting instances of overcapacity for these scenarios. Specifically, one sample exhibited penalty delays in scenario 2, while six samples encountered penalty delays in scenario 3. This underscores the potential exclusion of scenario 2 and 3 in these samples as a viable alternative for this level of traffic, owing to its limited capacity. However, note that scenario 2 and 3 provides significant fuel reduction in rest of the samples without penalty delay.

Particularly noteworthy is the high penalty delay observed in sample 4, as illustrated in Fig. 11. The subsequent figure provides an overview of the leg flights in sample 4.

Figure 11. Sample-4 delay distribution between aircraft.

This observation suggests that opting for scenario 2 may confer advantages in terms of both capacity and efficiency in scenarios characterised by medium traffic volumes except sample 4. There was no penalty delay for the scenario 2, except in sample 4. The capacity of scenario 2 appears capable of accommodating up to 45 aircraft per hour in each PMS. Nonetheless, it’s imperative to acknowledge that as traffic volume escalates, the capacity of scenario 2 may prove insufficient. It’s crucial to emphasise that these capacity assessments solely pertain to arrival operations. Once departures are factored into the system, capacity saturation may occur at lower traffic levels.

3.3 High traffic demand results (100 aircraft in one hour) with stochastic application

In this section, we present the results of high traffic demand samples involving 100 aircraft within a one-hour timeframe, employing a stochastic approach. Previous analyses indicated that the classical mathematical model struggled to find optimal solutions even in deterministic cases, particularly evident in scenarios with 90 aircraft. However, in stochastic approach, the model is solved for all three scenarios simultaneously, aiming to identify a single sequence applicable across all scenarios. In this case, a solution cannot be found using classical stochastic modeling. Therefore, we were able to compare our developed stochastic matheuristic approach with the FCFS scheduling.

This approach facilitates controller application of the sequence, even in situations where last-minute changes to leg configurations occur, thereby ensuring operational efficiency is maintained. Despite the challenges encountered with the classical stochastic mathematical model, our stochastic matheuristic model endeavors to discover a single sequence feasible for each scenario, ultimately offering improved efficiency compared to FCFS methodologies.

In Fig. 12, we present the results from 10 sample runs for the 100 aircraft demand scenarios.

Figure 12. Comparison of matheuristic and FCFS.

While the average of all samples appears to be close to each other, when examined on a sample-by-sample basis, there were improvements of up to 3.4% for delay and up to 0.5% per aircraft for fuel consumption. Figure 13 illustrates the sample-based improvements of stochastic matheuristic model compared to FCFS.

Figure 13. Sample-based improvements.

Note that the presented neighbourhood structure is a rather simple approach, especially considering that the allotted time for the matheuristic model is only 15 minutes. Given these constraints, the observed improvements can be deemed reasonable. We conducted a sensitivity analysis to explore what the model could achieve over longer durations. However, it will be addressed in detail in another study to enhance the neighbourhood structure. Randomly, Sample-9 was chosen for the sensitivity analysis. The results of the solution time sensitivity analysis are provided in Fig. 14.

Figure 14. Sensitivity analysis results (per AC).

This sensitivity analysis indicates the potential for improvements in average delay per aircraft and average fuel consumption per aircraft by extending the solution time constraint of the matheuristic (MH). Considering the cumulative savings, these findings provide optimism for future studies. Note that further enhancements to current neighbourhood structure may yield even greater improvements.

The delay distributions and the sequence obtained by the matheuristic model in 15 minutes are depicted in Figs 15 and 16 for Sample-9.

Figure 15. Delay distributions in Sample-9.

Figure 16. Position change in Sample-9 compared to FCFS.

Figure 16 illustrates the position shifts compared to the FCFS in the stochastic matheuristic model. As depicted in Fig. 12, the positions of 15 aircraft are altered, resulting in a 1% and 0.23% reduction in delay and fuel consumption, respectively. Given the significant demand for aircraft, the classical stochastic mathematical model fails to find any feasible solution. In contrast, our stochastic matheuristic approach offers feasible and efficient solutions compared to FCFS, providing hope for further developments.

4.0 Discussion and conclusion

This study introduces a novel stochastic matheuristic approach tailored for the dynamic utilisation of arrival route structures, specifically focusing on the modeling of Istanbul Airport’s P-PMS. The utilisation of P-PMSs in congested airports like Istanbul Airport has proven pivotal in enhancing air traffic control efficiency. However, the advent of the COVID-19 pandemic has disrupted air traffic patterns, rendering traditional long sequencing legs inefficient. As air traffic rebounds, the dynamic utilisation of sequencing legs becomes imperative, necessitating adaptable solutions to manage varying traffic densities throughout the day.

To tackle the intricate nature of dynamic sequencing legs, a stochastic matheuristic model is proposed, combining the strengths of mathematical models with heuristic adaptability. Traditional stochastic models may falter in providing feasible sequencing solutions, especially under high traffic demands and configuration changes. In contrast, the matheuristic approach offers robust and efficient sequencing, even amidst configuration alterations, ensuring optimal air traffic management.

The study conducted a comprehensive analysis across varying traffic demand scenarios, comparing outcomes derived from classical stochastic mathematical modeling and stochastic matheuristic model. Notably, under low traffic demand scenarios, scenarios incorporating partial overlap and full dissociation of sequencing legs showcased significant fuel savings per aircraft, indicating the efficacy of dynamic utilisation. As traffic demands escalated, the matheuristic approach outperformed traditional FCFS strategies, demonstrating improvements in delay and fuel consumption.

The examination of various leg configuration scenarios under low traffic demand scenarios revealed notable reductions in fuel consumption, particularly in scenarios 2 and 3 compared to scenario 1. These reductions stemmed from distinct flight behaviours influenced by structural constraints inherent in each scenario. Scenario 3, with fully dissociated legs, demonstrated potential advantages in reducing fuel usage and air traffic controller workload. However, it’s essential to acknowledge potential limitations in accommodating higher traffic volumes effectively, necessitating further investigation.

The analysis of medium traffic demand scenarios highlighted significant fuel savings in scenarios 2 and 3 compared to scenario 1. However, high number of instances of penalty delays in scenario 3 indicated capacity limitations, particularly in accommodating overcapacity conditions. The observed advantages of scenario 2, both in terms of capacity and efficiency, suggest its suitability for medium traffic volumes, albeit with considerations for potential capacity saturation as traffic levels increase.

The introduction of a stochastic matheuristic approach for high traffic demand scenarios demonstrated promising improvements in operational efficiency compared to traditional FCFS scheduling. Despite challenges encountered with classical stochastic mathematical models, the developed stochastic matheuristic approach aimed to identify feasible sequences applicable across scenarios, leading to reductions in delay and fuel consumption. Despite the promising results, challenges remain in refining the stochastic matheuristic approach, particularly in enhancing the neighbourhood structure. Sensitivity analyses considering the solution time limits indicate potential for further improvements, warranting future research endeavors to optimise the model’s efficiency and efficacy.

The implications of the proposed stochastic matheuristic approach extend beyond theoretical frameworks, offering tangible benefits in operational air traffic management. By providing feasible sequencing solutions adaptable to configuration changes, air traffic controllers can effectively manage fluctuating traffic densities without compromising efficiency or safety, thus optimising airport operations and mitigating potential congestion issues.

In conclusion, the study underscores the importance of dynamic utilisation in arrival route structures, particularly in the context of congested airports like Istanbul Airport. The proposed stochastic matheuristic approach emerges as a promising solution to address the complexities associated with dynamic sequencing legs, offering efficient, robust and adaptable sequencing solutions. As air traffic management continues to evolve, leveraging innovative approaches like the one presented herein will be essential in ensuring safe, efficient and sustainable airport operations.

5.0 Limitations and future works

Despite the promising results obtained in this study, there are certain limitations that should be acknowledged. These limitations provide opportunities for future research to further enhance and expand upon the findings presented here.

The current study exclusively focused on arrival sequencing and did not consider departures to isolate the effect of leg usages on arrivals. Including departures, on the other hand, in the sequencing process is crucial, as departure aircraft interact with arrival flows and impact the overall system dynamics. Future studies should aim to incorporate departure operations and explore the interdependencies between arrivals and departures to develop more realistic and efficient sequencing algorithms.

The probability of each leg usage scenario is assumed as equal in this study. However, the actual probability distributions of different configurations may vary across different airports and operational contexts. Future research could collect empirical data on configuration implementation frequencies to accurately model the probability distributions, allowing for more precise optimisation of the objective function.

The proposed stochastic matheuristic approach demonstrated improved performance compared to traditional stochastic mathematical model, but there is still room for further optimisation. Future studies could explore alternative algorithms such as genetic algorithm or particle swarm optimisation, fine-tune the parameter settings, or investigate hybrid approaches that combine the strengths of different optimisation techniques to achieve even better solution performance. Furthermore, our current neighbourhood strategy represents a basic approach. However, the development of several strategies has the potential to enhance the efficiency of the matheuristic approach.

The findings of this study were based on simulations and generic scenarios due to the lack of real data. However, validating the proposed approaches using real-world data from operational environments would provide more robust evidence of their effectiveness and applicability. Collaborations with air traffic control organisations and access to actual operational data would be invaluable in verifying the performance of the models and algorithms.

This study focused on leg usage uncertainties. However, other uncertainties are inherent in air traffic management, such as weather disruptions or unexpected events [37]. Future research should investigate the incorporation of these uncertainties into the modeling and optimisation processes to develop more robust and adaptive sequencing approaches that can handle unexpected disruptions effectively.

In summary, while this study has made significant contributions to the understanding and improvement of arrival sequencing in air traffic management, there are several avenues for future research. Addressing the limitations discussed earlier, including more realistic scenarios, considering departures, refining the optimisation algorithms, validating with real-world data and incorporating other uncertainties, will further enhance the applicability and effectiveness of the proposed approaches in real-world operational settings.