2.1 Matheuristics in routing, scheduling and logistics

The integration of exact solvers within heuristic frameworks has proven particularly effective in logistics. Archetti and Speranza [Reference Archetti and Speranza24] classify these approaches into decomposition, improvement and column generation-based methods. Solyal and Süral [Reference Solyal and Süral25] addressed the inventory routing problem using a matheuristic that sequentially solves different MILPs, demonstrating that this hybrid approach could outperform existing algorithms and provide new best solutions. Similarly, Boschetti and Maniezzo [Reference Boschetti and Maniezzo26] applied a matheuristic to city logistics by combining a dual ascent procedure with a randomised constructive heuristic, which proved effective in real-world scenarios. Vadseth et al. [Reference Vadseth, Andersson and Stålhane27] presented an iterative matheuristic for inventory routing that emphasised the maximum level replenishment policy and outperformed state-of-the-art exact methods on large benchmark instances. Furthermore, Moura et al. [Reference Moura, Pinto, Alves and de Carvalho28] introduced a hybrid approach for the combined bin packing and vehicle routing problem, demonstrating significant improvements in operational efficiency under real-world constraints.

In the broader context of transportation networks, Chouman and Crainic [Reference Chouman and Crainic29] developed a cutting-plane matheuristic that efficiently computes tight lower bounds, surpassing current best methods in the literature. Villegas et al. [Reference Villegas, Prins, Prodhon, Saldanha-da Gama and Santos-Pereira30] proposed a hybrid method for truck and trailer routing using a set-partitioning formulation that outperformed leading methods and highlighted the robustness of matheuristics. Additionally, Angelelli et al. [Reference Angelelli, Archetti and Peirano31] focused on freight forwarder services and validated their method on real data instances, showing its capability to offer high-quality solutions for large-scale problems within reasonable computing times.

In industrial scheduling, matheuristics often bridge the gap between theory and practice. Dang et al. [Reference Dang, van Diessen, Martagan and Adan32] tackled parallel machine scheduling using a genetic algorithm combined with linear programming, showing significant improvements over current practices. Fanjul-Peyro et al. [Reference Fanjul-Peyro, Perea and Ruiz33] proposed strategies for unrelated parallel machine scheduling that significantly outperformed mathematical models on medium-sized instances. Della Croce and Salassa [Reference Della Croce and Salassa34] introduced a variable neighbourhood search-based matheuristic for nurse rostering and demonstrated that hybrid methods can handle complex resource constraints better than pure solvers. Finally, Karlsson et al. [Reference Karlsson, Rönnberg, Stenberg and Uppman35] presented a matheuristic for large-scale avionic scheduling in integrated modular systems, effectively solving large instances and providing valuable managerial insights.

2.2 Matheuristics in aviation operations

Within the aviation sector, matheuristics have been successfully adapted to various ground and airside operations. Silva et al. [Reference da Silva, Ueno, dos Santos and de Carvalho36] addressed the aircraft landing problem (ALP) using a relax-and-fix approach that proved efficient for instances with a large number of aircraft and offered significant improvements over existing mathematical formulations. In maintenance operations, Zhang et al. [Reference Zhang, Chan, Chung and Fu37] combined ant colony optimisation with set partitioning to handle flexible cruise times. Their approach significantly improved aircraft utilisation and demonstrated high accuracy when validated with Bureau of Transportation Statistics data. Karsu and Solyal [Reference Karsu and Solyal38] employed a matheuristic for gate assignment to minimise passenger walking distances, outperforming existing algorithms in computational tests. Delgado and Mora [Reference Delgado and Mora39] applied column generation-based heuristics to air-cargo recovery and significantly improved profit margins by allowing flexible cargo assignment policies. Kyriakidou [Reference Kyriakidou40] addressed aircraft landing scheduling using variable neighbourhood search combined with interior-point methods, showing substantial improvements over traditional methods through statistical analysis.

2.3 The Gap in DSP: From black-box to explainable optimisation

Despite these successes, there is currently no matheuristic that embeds an exact MILP evaluator directly within the DSP as existing work relies predominantly on pure metaheuristics. Studies by Atkin et al. [Reference Atkin, Burke, Greenwood and Reeson41], Capr $\grave{i}$ and Ignaccolo [Reference Capr and Ignaccolo42] and Benlic et al. [Reference Benlic, Brownlee and Burke43] typically enforce feasibility through forward timing rules or ad-hoc repairs rather than through an exact optimisation layer. More recently, constraint-guided search frameworks have been introduced to inject structural information [Reference Riahi, Newton, Polash, Su and Sattar44], yet these still largely lack the rigorous guarantees of a MILP for complex ground-holding couplings.

Conversely, adoption barriers persist for exact methods due to their opacity. Even when a stand-alone MILP offers a detailed solution, the intermediate reordering logic remains hidden inside the branch-and-bound tree, forcing controllers to accept a ‘black box’ output without insight into the operational decisions that led to it. Recent work highlights that such opacity limits trust [Reference Aigner, Goerigk, Hartisch, Liers and Miehlich15, Reference Hurter, Degas, Guibert, Durand, Ferreira, Cavagnetto, Islam, Barua, Ahmed, Begum, Bonelli, Cartocci, Di Flumeri, Borghini, Babiloni and Aricó45] and that controllers prioritise tools aligning with their mental models [Reference Westin, Borst and Hilburn22]. To bridge this divide, this study draws on principles of explainable metaheuristics [Reference Almeida, Soares, Lezama, Limmer, Rodemann and Vale20, Reference Fyvie, McCall, Christie, Brownlee and Singh21], aiming to combine the transparency of trajectory-based search with the precision of mathematical programming.

4.1 Evaluation of fixed sequences

The core innovation of our search strategy is the replacement of heuristic objective functions with an exact MILP evaluator. For any candidate permutation generated by the metaheuristic, we solve a simplified version of the DSP model where the binary ordering variables are fixed. This transforms the NP-hard sequencing problem into a computationally tractable linear programming (LP) or reduced MILP problem focusing solely on resource allocation (Fig. 1).

The evaluation protocol proceeds as follows:

1. 1) The metaheuristic generates a candidate neighbour order using a constraint-aware move (e.g. a CPS-compliant swap).

2. 2) This candidate order is passed to the restricted MILP, where the binary variables defining the sequence are fixed.

3. 3) The restricted MILP optimises the remaining continuous variables: reaching times to each point  $\{{t_i^p}\}$ and hold durations $\{ {s_i^p}\}$ , ensuring optimal distribution of ground holds across the taxiway network.

4. 4) The MILP returns the optimal objective value (total hold fuel) and feasibility status, which are fed back to the metaheuristic to guide the search.

5. 5) Steps 1–4 repeat until the time or iteration budget is exhausted.

Frequency of evaluation: In each iteration, the metaheuristic filters candidates using quick feasibility checks (e.g. CPS bounds). Only promising candidates trigger a restricted-MILP solve. To ensure computational efficiency compatible with real-time operations, the MILP solver is warm-started from the previous solution, and the total number of MILP calls is capped by a global budget.

Computing environment: All experiments were performed on a Windows workstation equipped with an 11th Gen Intel^® Core $^{TM}$ i7-1185G7 CPU (3.00 GHz) and 64 GB of RAM. The algorithms were implemented in Python 3.11 using Pyomo and solved with Gurobi Optimizer (v12.0.3). Default solver parameters were used, except for the TimeLimit, which was kept consistent across all runs to ensure comparability.

4.2 MILP model

We formulate the DSP on a time–space representation of the airport surface. The model minimises total ground-hold fuel while enforcing runway wake separations, multi-point ground conflict separations and CPS. Below we define all sets, indices, parameters and decision variables before presenting the constraints.

Sets and indices:

• • i $,j \in I$ : set of aircraft.

• • $P$ : set of surface/operation points (pushback, conflict points, runway); for clarity we partition $P = {P^{{\rm{conf}}}} \cup {P^{{\rm{rwy}}}} \cup {P^{{\rm{oth}}}}$ .

• • ${P_i} \subseteq P$ : ordered path of aircraft $i$ (only the points actually traversed by $i$ ). We write $q{ \prec _i}p$ if point $q$ precedes $p$ on $i$ ’s path.

• • $O$ : wake-turbulence types (Heavy/Medium/Light); ${\tau _i} \in O$ denotes the type of aircraft $i$ .

Parameters:

• • ${\alpha _i}$ [s]: system entry (ready) time of aircraft $i$ .

• • ${\gamma _i}$ [s]: prescribed gate/pushback delay for $i$ (treated as a parameter).

• • $\theta _i^p$ [s]: nominal taxi time for $i$ to reach point $p \in {P_i}$ from the previous point on $i$ ’s path.

• • $\bar s_i^p$ [s]: maximum allowable local ground hold of $i$ at point $p \in {P_i}$ .

• • ${{\rm{\Delta }}_p}$ [s]: required time separation at conflict point $p \in {P^{{\rm{conf}}}}$ (point-dependent).

• • ${W_{{\tau _i},{\tau _j}}}$ [s]: wake separation at the runway when a type- ${\tau _i}$ aircraft is followed by a type- ${\tau _j}$ aircraft.

• • ${\phi _i}$ [kg/s]: hold-fuel rate parameter used to evaluate ground holds of aircraft $i$ .

• • ${\rm{MPS}}$ : constrained position-shifting limit (max allowed deviation from scheduled order).

• • $n = \left| I \right|$ , $M \gt 0$ : a sufficiently large constant (big- $M$ ).

Decision variables:

• • $t_i^p$ [s]: time aircraft $i$ reaches point $p \in {P_i}$ .

• • $s_i^p$ [s]: local ground-hold (delay) of aircraft $i$ at point $p \in {P_i}$ .

• • $x_{ij}^p \in \left\{ {0,1} \right\}$ : ordering at point $p$ (1 if $i$ is before $j$ at $p$ ), defined only for $p \in {P_i} \cap {P_j}$ .

• • $\pi _i^p$ : position of aircraft $i$ at point $p \in {P_i}$ ; $\pi _i^{{\rm{std}}}$ scheduled position; $\pi _i^{{\rm{fin}}}$ final (runway) position.

• • ${F_i}$ [kg]: fuel consumed by aircraft $i$ due to ground holds.

Initialization and propagation on paths

(1) \begin{align} \displaystyle t_i^0 & = {\alpha _i} + {\gamma _i}\qquad\qquad &\forall i &\in I\\[-10pt]\nonumber\end{align}

(2) \begin{align}\displaystyle t_i^p & = t_i^0 + \mathop \sum \limits_{q{\rm{\;}}{ \prec _i}{\rm{\;}}p} \left( {\theta _i^q + s_i^q} \right) \qquad\qquad\qquad &\forall i &\in I,{\rm{\;\;}}\forall p \in {P_i}\\[-8pt] \nonumber\end{align}

(3) \begin{align} \displaystyle 0\quad &\le \quad s_i^p \le \bar s_i^p \qquad\qquad\qquad &\forall i &\in I,{\rm{\;\;}}\forall p \in {P_i} \\[10pt] \nonumber\end{align}

Positions, CPS, and consistency

(4) \begin{align}\displaystyle \pi_{i}^{p} & = n - \mathop\sum\limits_{\begin{matrix}{j \in I} \\ {j \ne i}\end{matrix} }x_{ij}^p \qquad &\forall i &\in I,\quad \forall p \in P_i \\[-10pt] \nonumber\end{align}

(5) \begin{align}\displaystyle \pi _i^p{\rm{\;\;}} & \le {\rm{\;\;}}\pi _i^q{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}&\forall i &\in I,{\rm{\;\;}} \forall p,q \in {P_i}:{\rm{\;\;}}p{ \prec _i}q \\[-10pt] \nonumber\end{align}

(6) \begin{align}\displaystyle \left| {\pi _i^{{\rm{fin}}} - \pi _i^{{\rm{std}}}} \right|{\rm{\;}} & \le {\rm{\;\;MPS\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}&\forall i &\in I\\[10pt] \nonumber\end{align}

Conflict and runway separations (only where aircraft meet)

(7) \begin{align}\displaystyle t_j^p - t_i^p{\rm{\;\;}} &\geqslant {\rm{\;\;}}{{\rm{\Delta }}_p} - M\left( {1 - x_{ij}^p} \right){\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}&\forall p &\in {P^{{\rm{conf}}}},{\rm{\;\;}} \forall i \ne j:{\rm{\;\;}}p \in {P_i} \cap {P_j}\\[-10pt] \nonumber\end{align}

(8) \begin{align}\displaystyle t_i^p - t_j^p{\rm{\;\;}} & \geqslant {\rm{\;\;}}{{\rm{\Delta }}_p} - M{\rm{\;}}x_{ij}^p{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}&\forall p &\in {P^{{\rm{conf}}}},{\rm{\;\;}} \forall i \ne j:{\rm{\;\;}}p \in {P_i} \cap {P_j}\\[-10pt] \nonumber\end{align}

(9) \begin{align}\displaystyle t_j^p - t_i^p{\rm{\;\;}} & \geqslant {\rm{\;\;}}{W_{{\tau _i},{\tau _j}}} - M\left( {1 - x_{ij}^p} \right){\rm{\;\;\;\;\;\;\;\;}}&\forall p &\in {P^{{\rm{rwy}}}},{\rm{\;\;}}\forall i \ne j:{\rm{\;\;}}p \in {P_i} \cap {P_j} \\[-10pt] \nonumber\end{align}

(10) \begin{align}\displaystyle t_i^p - t_j^p{\rm{\;\;}} & \geqslant {\rm{\;\;}}{W_{{\tau _j},{\tau _i}}} - M{\rm{\;}}x_{ij}^p{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}&\forall p &\in {P^{{\rm{rwy}}}},{\rm{\;\;}} \forall i \ne j:{\rm{\;\;}}p \in {P_i} \cap {P_j} \\[10pt] \nonumber\end{align}

Hold Fuel accounting and objective

(11) \begin{align}\displaystyle {F_i}{\rm{\;}} & = {\rm{\;}}{\phi _i}\mathop \sum \limits_{p \in {P_i}} s_i^p{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}} \forall i \in I\\[-10pt] \nonumber\end{align}

(12) \begin{align} &\displaystyle {\rm{min\;\;\;\;}}\mathop \sum \limits_{i \in I} {F_i}\\[10pt] \nonumber\end{align}

Model overview: Equation (1) sets each aircraft’s pushback time from ready time and a fixed gate delay (parameter). Equation (2) propagates times along the aircraft-specific path ${P_i}$ , summing nominal taxi and optimised holds only over points actually traversed. Bounds in Equation (3) cap local holds. Positions $\pi $ are derived from order binaries $x$ at the points where aircraft meet (Equation (4)), kept nondecreasing along each path (Equation (5)), and restricted by CPS (Equation (6)) to stay close to the scheduled order. Conflict-point separations use the point-dependent parameter ${{\rm{\Delta }}_p}$ (Equations (7) and (8)); runway separations use the wake matrix ${W_{{\tau _i},{\tau _j}}}$ (Equations (9) and (10)). Fuel is linear in holds (Equation (11)), and the objective minimises total hold fuel (Equation (12)). Emissions (HC/CO/NOx) are reported post-solve from ${F_i}$ and type-specific indices. Importantly, if an aircraft does not traverse a point (e.g. different apron), that point is not in ${P_i}$ , so timing and separation constraints never apply spuriously. The big- $M$ terms appear only in disjunctive separations and do not affect propagation or hold-fuel accounting.

4.3 Tabu search (TS)

The TS is a metaheuristic that explores the solution space by iteratively moving from one potential solution to another, avoiding cycles through a tabu list [Reference Glover46–Reference Cecen and Çetek48]. Our enhanced TS incorporates custom features such as dynamic, reactive tabu tenure adjustment, priority pairs based on ETD values, CPS compliance and a small elite pool with path-relinking/micro-scale large neighbourhood search (LNS) to improve its effectiveness in aircraft sequencing. Dynamic tabu tenure is updated reactively: it increases when the search stagnates and may be relaxed again after diversification steps. Priority pairs are grouped based on estimated time of departure (ETD) values, guiding the search towards feasible solutions that respect temporal constraints. CPS compliance ensures that position shifts for each aircraft do not exceed a certain limit, maintaining fairness and minimising controller workload. The TS also maintains a small elite pool of best sequences and periodically applies path-relinking and a micro-scale LNS between the current and elite solutions, which intensifies the search without violating CPS.

The system initialises the search state, best solution and parameters, generates a neighbourhood of CPS-compliant candidates using allowed swap/insert moves, and evaluates only a small top- $k$ subset of them using the MILP model based on heuristic scores, in order to respect the runtime budget. Moves are checked against the tabu list and the aspiration criterion, and the tabu list is updated to prevent cycling back to recent moves. The system adjusts tabu tenure reactively and performs random CPS-compliant shuffles, together with occasional elite-based path-relinking/micro-LNS, to escape local optima. The search stops when a global time limit, a MILP evaluation budget or the maximum number of non-improving iterations is reached. The pseudo-code of the enhanced TS approach is given in Algorithm 1.

Algorithm 1

Enhanced Tabu Search

Algorithm 1 Long description

The equation outlines an enhanced tabu search algorithm. It initializes the search state, best solution, and parameters. The algorithm generates a neighborhood of CPS-compliant candidates using allowed swap or insert moves and evaluates a small top-k subset of these candidates using a MILP model based on heuristic scores. Moves are checked against a tabu list and an aspiration criterion, and the tabu list is updated to prevent cycling back to recent moves. The system adjusts tabu tenure reactively and performs random CPS-compliant shuffles, along with occasional elite-based path-relinking or micro-LNS, to escape local optima. The search stops when a global time limit, a MILP evaluation budget, or the maximum number of non-improving iterations is reached.

Algorithm 2.

Enhanced Simulated Annealing

Algorithm 2 Long description

The equation outlines the Enhanced Simulated Annealing algorithm. It begins by setting initial parameters such as the initial state, maximum temperature, minimum temperature, alpha, maximum stagnation, maximum no improvement, time limit, and maximum evaluations. The algorithm sets the current solution, best solution, temperature, iteration, stagnation counter, no improvement counter, number of evaluations, and start time. It then enters a loop that continues while the temperature is greater than the minimum temperature, the no improvement counter is less than the maximum no improvement, the current time minus the start time is less than the time limit, and the number of evaluations is less than the maximum evaluations. Within the loop, a new solution is generated, and its cost is evaluated. The delta cost is calculated, and the current solution is updated based on the cost comparison and a random value. The best solution, stagnation counter, and no improvement counter are updated accordingly. The temperature is updated, and if the stagnation counter exceeds the maximum stagnation, a compliant shuffle is performed. The loop continues until the conditions are met, and the best solution is returned.

4.4 Simulated annealing (SA)

The SA algorithm is a simulation technique that draws inspiration from metallurgy, where materials are heated and cooled to achieve a stable crystalline structure [Reference Kirkpatrick, Gelatt and Vecchi49]. It uses a temperature parameter to accept worse solutions and escape local optima. Custom enhancements include an adaptive temperature schedule, CPS-compliant moves and priority-based moves, which balance exploration and exploitation effectively.

The enhanced SA system involves initialisation, generation of new CPS-compliant solutions, evaluation via the restricted MILP, acceptance criteria, temperature update and stagnation check. It balances exploration and exploitation through an adaptive temperature schedule, avoiding premature convergence. CPS-compliant moves ensure feasible and practical solutions, even under complex operational constraints. The system also performs a shuffle if stagnation is detected to escape local optima. In practice, the SA loop is additionally capped by a global wall-clock time limit and a maximum number of MILP evaluations to keep its runtime comparable with the other metaheuristics. The pseudo-code of the enhanced SA approach is given in Algorithm 2.

4.5 Multi-ant colony optimisation (MACO)

Ant colony optimisation (ACO) is a method inspired by ants’ behaviour, where they deposit pheromones on paths to reinforce efficient routes [Reference Dorigo, Di Caro and Gambardella50]. Multi-ACO uses one or more colonies, each searching independently but sharing information periodically, to enhance global search capability [Reference Ting and Chen51]. In our implementation, the number of colonies is kept small, which is sufficient to maintain diversity while staying within the restricted-MILP evaluation budget. Custom enhancements include multiple colonies running in parallel when needed, dynamic pheromone updating and periodic inter-colony communication to propagate good solutions across the search space. These features enhance diversity and reduce the likelihood of premature convergence.

The initialisation process involves setting up multiple colonies, each with its own pheromone matrix and initial state. Ants in a colony use pheromone information and heuristics like CPS compliance and priority pairs to construct solutions. Each constructed sequence is then evaluated either by the restricted MILP or a fast-timing surrogate, depending on the configuration. The pheromone update is a system where pheromones are updated based on the quality of a solution to guide future ants toward better solutions. Inter-colony communication involves sharing the best solutions among colonies to facilitate the spread of beneficial solutions. The global best update updates the global best solution across all colonies. The pseudo-code of the enhanced MACO approach is given in Algorithm 3.

Algorithm 3.

Enhanced Multi-Ant Colony Optimisation

Algorithm 3 Long description

The equation outlines the steps of an enhanced multi-ant colony optimization algorithm. It begins with initializing multiple colonies, each with a pheromone matrix and initial state. Ants within each colony use pheromone information and heuristics such as CPS compliance and priority pairs to construct solutions. Each constructed sequence is evaluated using either a restricted mixed-integer linear programming (MILP) model or a fast-timing surrogate, depending on the configuration. The pheromone update mechanism adjusts pheromone levels based on the quality of the solutions to guide future ants toward better solutions. Inter-colony communication involves sharing the best solutions among colonies to spread beneficial solutions. The global best update ensures that the best solution across all colonies is identified and updated. The pseudo-code for this enhanced MACO approach is detailed in Algorithm 3.

4.6 Bacterial foraging optimisation (BFOA)

BFOA is a method that uses bacteria’s foraging behaviour to find optimal solutions [Reference Passino52]. It involves tumbling and swimming bacteria to find the most nutritious regions. BFOA’s custom enhancements include CPS compliance, adaptive swimming length and dynamic reproduction and elimination-dispersion phases. These features ensure feasible solutions and maintain diversity through adaptive reproduction and elimination-dispersion phases (see Ref. [Reference Niu, Fan, Tan, Rao, Li, Huang, McGinnity, Heutte and Zhang53] for further information about BFOA).

The process involves initialising the bacteria population with feasible solutions and setting parameters for different phases. Each bacterium moves randomly through the solution space, performing a tumble, and then swims toward improved solutions. Similar to the other metaheuristics, BFOA uses the same CPS-compliant swap/insert operators and ETD-based priority pairs when generating new candidate sequences. Reproduction involves retaining and cloning the most efficient bacteria, allowing the rest to be repurposed to concentrate on the most promising areas of the solution space. The elimination-dispersion method introduces diversity by randomly eliminating and redistributing certain bacteria. The global best update updates the global best solution for all bacteria. The enhanced BFOA prevents premature convergence and ensures solutions respect operational constraints like CPS and priority pairs, making them suitable for real-world implementation. The pseudo-code of the enhanced BFOA approach is given in Algorithm 4.

Algorithm 4.

Enhanced Bacterial Foraging Optimisation

Algorithm 4 Long description

The image displays a detailed flowchart outlining the steps of the enhanced bacterial foraging optimization algorithm. It begins with initializing the bacteria population and setting parameters for different phases. Each bacterium performs random tumble and swim movements to explore the solution space. The algorithm uses CPS-compliant swap/insert operators and ETD-based priority pairs for generating new candidate sequences. Reproduction involves retaining and cloning the most efficient bacteria, while the elimination-dispersion method introduces diversity by randomly redistributing bacteria. The global best update ensures the best solution is tracked across all bacteria. This process prevents premature convergence and respects operational constraints, making it suitable for real-world applications.

4.7 Custom features

In our matheuristic approach, several custom features are shared across all algorithms, while others are tailored to specific metaheuristics based on their inherent structure and operational characteristics. Priority pair selection, CPS constraints and hold-fuel-weighted swap logic are universally relevant and applied across all four metaheuristics. In contrast, features like probabilistic acceptance are used only where they are particularly effective (e.g. in SA), while random shuffling is incorporated into algorithms like TS and SA to handle stagnation effectively. This selective application of features ensures that each algorithm operates optimally, leveraging the strengths of its design while remaining adaptable to the DSP.

Priority pair selection: First, priority pairs are determined based on the ETD closeness between aircraft. Aircraft with ETDs within a specified threshold (ETD_threshold = 90 s) are considered priority pairs. This strategy ensures that the algorithm focuses on swaps between aircraft with close departure times, which are likely to have a more significant impact on the overall solution. When no priority pairs remain, the algorithm selects two random aircraft to either swap or insert, maintaining diversity in solution space exploration. The same priority-pair generator is shared across TS, SA, MACO and BFOA so that all algorithms focus first on ETD-critical swaps.

Swap and insert operations: At each step, the algorithm selects either a swap or an insert move. A swap exchanges the positions of two aircraft in the sequence, whereas an insert removes an aircraft from its current position and reinserts it at another position.

CPS compliance: Both swap and insert operations respect CPS constraints, ensuring that the new positions of the aircraft stay within a $ \pm 3$ position range relative to their ETD-based positions. This approach helps maintain a feasible sequence while staying close to the FCFS order.

Hold-fuel weight-based swap logic: To further refine the swaps, all metaheuristics use a hold-fuel weight-based swap probability when proposing swap moves. Aircraft with higher hold-fuel weights are given priority in swaps, as these have a more substantial impact on hold-fuel consumption and overall solution quality. The swap probability is computed as the ratio of the fuel weight of the first aircraft to the combined hold-fuel weights of both aircraft:

\begin{equation*} {\rm{swap}}\_{\rm{probability}} = \frac{{{\rm{holdfuel}}\_{\rm{weigh}}{{\rm{t}}_a}}}{{{\rm{holdfuel}}\_{\rm{weigh}}{{\rm{t}}_a} + {\rm{holdfuel}}\_{\rm{weigh}}{{\rm{t}}_b}}}{\rm{\;}}.\end{equation*}

This logic ensures that more hold-fuel-efficient decisions are prioritised across all neighbourhood generators. The hold-fuel weight data for each aircraft type is obtained from the International Civil Aviation Organisation (ICAO) engine emissions databank [54], providing accurate hold-fuel consumption estimates.

Selection of aircraft for swap: Two aircraft, denoted $a$ and $b$ , are selected either randomly or based on priority pairs. Their hold-fuel weights are calculated based on their respective aircraft types, and the probability of performing a swap is then determined by the above formula. If a random number generated is less than the calculated swap probability, the swap is executed. Otherwise, the algorithm either performs a swap or an insert operation as part of its regular process.

Probabilistic acceptance of solutions for SA: The SA algorithm employs a probabilistic acceptance function based on the Boltzmann distribution to determine whether to accept a newly generated sequence. This function allows worse solutions to be accepted with a probability that decreases as the temperature drops, helping the algorithm escape local minima. The acceptance probability is computed as:

\begin{equation*} P = {e^{ - {\rm{\Delta }}E/T}}\end{equation*}

where ${\rm{\Delta }}E$ is the change in energy (hold-fuel consumption) between the current and the new solution, and $T$ is the current temperature. The probabilistic acceptance mechanism helps explore a wider solution space, increasing the chances of escaping local optima and finding the global optimum.

Stagnation and random shuffle: If the algorithm detects stagnation (i.e. no improvement over a set number of iterations), it performs a CPS-compliant random shuffle to explore new regions of the solution space. This operation swaps aircraft positions within the CPS constraints to ensure feasibility while introducing larger changes to the sequence. Random shuffling is a less-frequent operation compared to regular swaps but is crucial in preventing the algorithm from being stuck in local optima. While the regular swap makes small incremental changes, the random shuffle introduces more substantial variations when necessary.

Termination conditions: All metaheuristics incorporate three key termination criteria: a wall-clock time limit (e.g. 600 s), a cap on the total number of restricted-MILP evaluations per run and a maximum number of iterations without improvement. These criteria ensure that the algorithms run within a practical duration, maintain a fair MILP budget relative to the standalone MILP and prevent infinite loops if no improvement is observed. Once any of these conditions is met, the algorithm terminates and returns the best solution found.