A chance-constrained programming model for airport ground movement optimisation with taxi time uncertainties

Airport ground movement remains a major bottleneck for air traffic management. Existing approaches have developed several routing allocation methods to address this problem, in which the taxi time traversing each segment of the taxiways is fixed. However, taxi time is typically difficult to estimate in advance, since its uncertainties are inherent in the airport ground movement optimisation due to various unmodelled and unpredictable factors. To address the optimisation of taxi time under uncertainty, we introduce a chance-constrained programming model with sample approximation, in which a set of scenarios is generated in accordance with taxi time distributions. A modified sequential quickest path searching algorithm with local heuristic is then designed to minimise the entire taxi time. Working with real-world data at an international airport, we compare our proposed method with the state-of-the-art algorithms. Extensive simulations indicate that our proposed method efficiently allocates routes with smaller taxiing time, as well as fewer aircraft stops during the taxiing process.


Introduction
The International Air Transport Association forecasts current passenger numbers could double to 8.2 billion no later than 2040 (International Air Transport Association, 2020a). Despite the global pandemic of 2020, it is believed that the longterm growth of aviation is likely to resume in time (International Air Transport Association, 2020c). However, this increasing trend for air transport, and the economic benefits driven by aviation, could be curtailed due to limited airport infrastructure capacities (International Air Transport Association, 2020b). This bottleneck has been recognised as one of the greatest challenges to future aviation industry (European Commission, 2020), and thus there is an urgent need to develop a more efficient decision support system for airport operations.
As a major component of airport operations, the airport ground movement (AGM)  which allocates aircraft routes from the gate to the runway or vice versa, is closely related to other airport operation subproblems, e.g., runway scheduling (Atkin et al., 2007) and gate allocation (Benlic et al., 2017;Bagamanova and Mota, 2020a,b), and thus has a huge impact on flights throughput. Although only a small part of air traffic management procedure, AGM contributes a large portion of fuel consumption and emissions for flight operations. This is reasonable since the aircraft engine is designed for cruising speed in the air, while it operates inefficiently on the airport ground at low speed. In accordance with the industry report in Ganev (2017), the fuel consumption due to aircraft taxiing at the congested airport makes up around 6% of the whole fuel consumption, leading to 5 million tonnes of fuel cost throughout the world. The AGM on the taxiways is therefore identified as a critical problem and has received increasing attention in recent years (Atkin et al., 2010;Weiszer et al., 2015;Brownlee et al., 2018a). A comprehensive literature review in AGM optimisation will be provided in the next section.
Although various methods have been applied in AGM optimisation, one critical issue remaining largely unaddressed by existing research (notable exceptions in Koeners and Rademaker, 2011;Lee and Balakrishnan, 2012;Brownlee et al., 2018a) is the taxi time uncertainties inherent in AGM optimisation. The aircraft taxi time is typically difficult to predict, since various factors, e.g., outside environment and airport layout, may heavily impact the actual taxi time (Wang et al., 2021). The majority of the existing AGM optimisation studies neglected the taxi time uncertainties, and simply modelled the taxi time for each taxiway as a fixed value. However, the simplification of taxi time modelling can lead to conflicts when an aircraft arrives at certain taxiway point before or after the predicted time. This issue has been partially addressed by a series of research for AGM online scheduling (Roling and Visser, 2008;Roling, 2009;Evertse and Visser, 2017;Brownlee et al., 2018b;Zou et al., 2018), in which the AGM has been modelled as a dynamic routing problem, and the routing results are updated through iteratively solving the routing problem with the latest information. These reactive scheduling models are promising solutions to address the AGM with taxi time uncertainties, however, the online scheduling requirements are difficult to satisfy, especially when incorporating real-time practical speed profiles with respect to operational constraints . Also, online scheduling methods need to be computationally fast, which may not be the case if more complex routing methods are to be employed in future (e.g., multi-objective routing).
Instead, developing proactive methods, which address uncertainties with distribution information in advance, could be an alternative for AGM optimisation with taxi time uncertainties. Previously, simulation based models were established to analyse the influence of the uncertainties of AGM, and it was established that they lead to delays (Koeners and Rademaker, 2011;Lee and Balakrishnan, 2012). Scala et al. (2021) integrated the near-airspace landing aircraft sequencing and AGM optimisation together, and designed a closed-loop feedback simulation framework under uncertainty to address the integrated airport capacity management problem. Brownlee et al. (2018a) first presented the taxi time uncertainties with fuzzy membership functions, and extended the quickest path problem with time windows (QPPTW) algorithm (Ravizza et al., 2014a) to generate multiple routes for different uncertainty levels. However, the modified Fuzzy-QPPTW algorithm typically estimates taxi time with worst cases, and therefore allocates longer and over-conservative routes to avoid aircraft conflicts. This leads to wasted capacity at the airport and unnecessary increased delays and emissions.
To better incorporate the taxi time uncertainties into the AGM model and generate more effective aircraft routes, this paper introduces the chance-constrained programming (CCP) model with sample approximate method (Charnes and Cooper, 1959;Luedtke and Ahmed, 2008), and proposes a modified CCP-QPPTW algorithm with a local heuristic method to allocate more efficient aircraft routes that are still robust. Note the CCP-QPPTW and Fuzzy-QPPTW are both developed based on the deterministic QPPTW routing method. However, they address the taxi time uncertainties in different ways. CCP-QPPTW adopts the CCP model to approximate the taxi time uncertainties and then applies QPPTW under each generated scenario, while Fuzzy-QPPTW relates to an adaptive Mamdani fuzzy rule based system to handle the fuzzy taxi time estimates; the fuzzy time estimates are propagated along the length of the route and used to add extra reserved times to the edge reservations. More specifically, the scheduled aircraft routes are obtained considering different conservative degrees, i.e., the confidence level in the CCP-QPPTW algorithm, thus avoiding over-conservative routes and reducing total taxi times. Besides, the proposed method is computationally efficient as compared to the state-of-the-art algorithm, which we will show later. To the best of our knowledge, we are the first to introduce CCP in the domain of AGM optimisation, while CCP models have been widely adopted in general scheduling problems (Wang and Ning, 2018;Zandkarimkhani et al., 2020) and engineering applications (Pozo and Contreras, 2012;Nadeem et al., 2018;Han et al., 2020).
Therefore the overarching original contributions of this work can be summarised as follows: (1) For the first time we introduce the CCP to address the AGM problem with taxi time uncertainties. Building upon the CCP model, a set of scenarios is generated to approximately represent the uncertainties, and a modified sequential quickest path routing method based on the multi-scenarios is proposed to tackle the uncertain AGM problem. (2) The proposed CCP model and the routing method are validated using real-world data at an international airport. (3) By comparing with the state-of-the-art algorithms, we show that the proposed method achieves superior performance both in terms of entire taxi times, number of aircraft stops, and computational efficiency. This paper is organised as follows: Section 2 provides a comprehensive literature review for the AGM optimisation; Section 3 describes the AGM optimisation problem in detail and formulates a CCP model to address the AGM with taxi time uncertainties; a modified sequential shortest path searching method based on CCP is developed in Section 4, followed by the computational results and discussions presented in Section 5; finally, conclusions are drawn in Section 6, in which the important contributions of our work and potential future directions are highlighted.

Literature review
To clearly outline the state of the art of AGM studies, the reviewed literature has been classified into three different modelling techniques: mixed-integer linear programming (MILP), graph-based and simulation-based models. An overview of the AGM research is listed in Table 1.
The MILP models, in which the scheduling & routing results can be efficiently obtained by commercial solvers, have been contributed to the AGM over the last decade. Roling and Visser (2008), Roling (2009) formulated the AGM as a MILP model, X. Wang et al.

Table 1
Overview of the studies in airport ground movement.

Madrid International
Roling&Visser (Roling and Visser, 2008;Roling, 2009)  in which the uncertain aircraft landing/pushback times are dynamically adapted to produce robust solutions. The capabilities of the proposed model were illustrated in a simple hypothetical airport. Evertse and Visser (2017) considered an online version of AGM, aiming to minimise the emissions resulting from the taxiing movements. The problem was described as MILP formulations with weighted multiple objectives, e.g., aircraft departure/arrival time deviations, taxi time, emissions, and fuel consumption. Similar to the above studies, Zou et al. (2018) established an integer programming model, in which the iterative two-stage scheduling strategy was proposed to minimise the taxi time of AGM. All aircraft were assigned the initial 4 dimensional routes by an integer programming model in the first stage, and the landing aircraft which have unavailable initial routes due to uncertainties were rescheduled using a shortest path based algorithm. Simulation results of the Beijing Capital International Airport demonstrated superior performance of the proposed method.
Another stream of approaches to address AGM optimisation is graph-based modelling techniques. In accordance with a decision support system based on graph theory, García et al. (2005) developed an integrated approach composed of a genetic algorithm (GA) and a time-space dynamic flow management algorithm for real-world airport ground operations.
To address the airport environmental issues, Ravizza et al. (2013b) modelled the AGM with the trade-off between taxi time and fuel consumption. A combination of a graph-based routing algorithm (Climaco and Martins, 1982) and a population adaptive immune algorithm (Chen and Mahfouf, 2006) to generate Pareto speed profiles was designed for the multi-objective model. Ravizza et al. (2014a) further addressed the AGM problem with a sequential graph based algorithm, which can absorb as much waiting time for delay as possible at the stand. A lower bound was obtained as well. Weiszer et al. (2014) considered a multi-component optimisation model for the AGM, consisting of routing & scheduling problem and the speed profile optimisation problem, which aims to minimise the taxi time and fuel consumption. A heuristic with efficient computation procedure was designed to approximate the trade-off curve between the two objectives. Furthermore, Weiszer et al. (2015) combined several different airport ground problems, including aircraft routing & scheduling (i.e., AGM), 4-dimensional trajectory optimisation, runway scheduling and airport bus scheduling, and proposed a holistic economic optimisation framework. A multi-objective GA was adopted and verified on real data from an international hub airport. Weiszer et al. (2018Weiszer et al. ( , 2020 also proposed a preference-based multi-objective evolutionary optimisation framework for AGM. The evolutionary search algorithm used modified crowding distance to take into account cost of delay and fuel price. As a result, the search algorithm achieved faster convergence and potentially better solutions. Based on an iterated local search method in the modelled airport networks, a fast heuristic solution respecting all safety regulations was proposed in Lobato et al. (2015). The computational results demonstrated its superior performance compared with a MILP approach. Adacher et al. (2018) addressed the AGM optimisation, aiming to minimise the routed taxi delay time and the pollution emissions, which are presented with engine idle times. The alternative graph model was introduced to solve the problem, in which the taxi delay time was optimised in the first phase, followed by the minimisation of the emissions. Brownlee et al. (2018b) applied the A * algorithm to find the shortest taxi-time paths on AGM optimisation, where the aircraft allocation orderings are optimised through a rolling window approach incorporating GA. Simulations of three international airports demonstrated the GA reduced overall taxi time with respect to other approaches. Similarly, Dabachine et al. (2018) first introduced three routing algorithms, i.e., Dijkstra, bidirectional Dijkstra, and A * with different strengths and weaknesses. A new hybrid A * algorithm was therefore further designed to integrate the advantages of the three algorithms.
AGM optimisation has also been addressed with simulation-based methods, which are typically problem specific and can provide practical solutions with aircraft trajectories. Mori (2013) focused on the airport movement modelling accuracy and developed a new simulation method building upon the Nagel-Schreckenberg model, which is originally a car congestion model. The proposed model was validated at Tokyo International Airport. Zhang et al. (2016Zhang et al. ( , 2018) established a refined trajectory-based AGM optimisation model, and developed a dynamic routing and timing algorithm to minimise the shortest taxi time. The simulation results showed the advantages of the proposed method over existing approaches in Ravizza et al. (2014a). Yang et al. (2017) analysed the flow characteristics of airport ground network on both mesoscopic and macroscopic level, and proposed an efficient modelling approach for simulating congestion on taxiway and apron networks. The proposed model was shown to be an efficient and accurate method for supporting the AGM. Zhang et al. (2019) considered a multi-objective optimisation scheme for AGM under different operating environments, in which the taxi time, fuel consumption, and emissions may vary due to environments. The proposed optimisation problem was solved by a GA, which has been validated and compared based on three different China's airports.
Throughout the above literature review, the research gap for current AGM studies can be identified as follows: although various approaches have been proposed for AGM optimisation, limited research has proactively considered the uncertainties from the aircraft taxi time, which can heavily influence the practical AGM performance. Gotteland and Durand (2003) proposed an airport ground traffic simulation tool, in which the uncertainties were defined as a fixed percentage of the initial taxiing speed. Clare and Richards (2011) considered a coupled problem of airport taxiway routing (i.e., the ground movement) and runway scheduling, and applied a routing horizon during which the routes were determined and then updated iteratively. Koeners and Rademaker (2011) considered time margins to overcome the possible taxi time uncertainties, which is similar to the buffer taxi time in Ravizza et al. (2014a). Khadilkar and Balakrishnan (2014) regarded the AGM as a congestion control problem, and further developed a stochastic model to address the random taxi times. Notice these studies focused more on the analysis of the taxi times, and then reactively updated aircraft routes with real-time requirements. Besides, adding buffer time to absorb taxi time uncertainties is sensitive to parameter setting. If not set properly, it will result in worse performance with increased delays and emissions.
Recently, Brownlee et al. (2018a) applied the fuzzy rule-based systems to estimate aircraft taxi time and their uncertainties. Furthermore, a modified Fuzzy-QPPTW was developed for AGM optimisation. However, the results are over-conservative to some extent, as Fuzzy-QPPTW considers the worst scenarios for aircraft route allocation (We will identify the over-conservatism in the simulation section). How to properly incorporate the taxi time uncertainties into the AGM modelling, and proactively generate more practical routing results remains unexplored. We aim to fill this gap in the present work.

Chance-constrained programming for AGM
In this section, the problem of AGM with fixed taxi times is initially described and modelled in an undirected graph, followed by an extension building upon the CCP sample approximation model which addresses the taxi time uncertainties.

Problem description with fixed taxi times
The AGM allocates routes from the runway to the gate for arrival aircraft, or vice versa for departure flight. To clearly state the AGM problem, several assumptions and simplifications are provided as follows.
(i) The objective of the AGM optimisation is to minimise the sum of taxi time for all flights over a routing horizon.  (ii) The runway and gate are assumed to be predetermined by gate allocation and runway scheduling methods. Changes to runway sequence or gate allocations can be accommodated by iteratively using route allocation algorithm, which is outside the scope of this work. (iii) In the process of AGM optimisation, each edge on the taxiway graph can only have one aircraft at a time, and the aircraft must respect safe distances with other aircraft moving on the airport. (iv) In the premise that the aircraft can follow the runway schedule, the aircraft should be held at the gate as long as possible, so that the fuel consumption and emissions can be efficiently reduced. (v) The taxi time for aircraft traversing the taxiways is deterministic and can be calculated ahead of time according to airport operational and environmental conditions.
The AGM optimisation problem can now be established in line with Ravizza et al. (2014a), Brownlee et al. (2018a). As illustrated in Fig. 1, the airport layout is described as an undirected graph ( , ), where edges denote the taxiways and vertices the gates, junctions and intermediate taxiing points. For each edge ∈ , an associated set of taxi times is provided as inputs. As indicated in Assumption (v), the taxi time for each aircraft can be different due to environmental and operational conditions, and these different values have been included in . To maintain the safety constraints in Assumption (iii), each edge is restricted with no more than one aircraft at any time. The threshold of safety distance is set as 60 metres in accordance with previous research (Ravizza et al., 2014a;Brownlee et al., 2018a). Consequently, longer edges are divided into several shorter edges no more than 60 metres through adding intermediate taxiing points in advance. Notice that several aircraft can still taxi simultaneously on a long taxiway, as it has been divided into several shorter edges.
To realise the safety distance constraints in the graph, a set of time windows is further introduced to represent the edge availability. Here a time window is an interval when the edge is free to be traversed by an aircraft. Denote the set of aircraft as and the allocated conflict-free route as for each aircraft ∈ . When a conflict-free route is allocated to aircraft with specific taxi times for each edge along the route, the corresponding time windows list would be updated through deducting the already allocated taxiing interval. Moreover, the minimum safety distance needs to be satisfied through checking other edges nearby. In doing so, the graph is preprocessed to identify conflict edges set within minimum safety distance for each edge ∈ . When the edge is occupied by an aircraft within a certain time window, the corresponding time windows set + for each + ∈ is updated to block the same taxiing duration as well. Therefore the AGM optimisation problem can be loosely formulated as where ( ) is the entire taxi time for aircraft , ( )∕ ( ) is the route start/end time, ∕ is the predetermined aircraft landing/taking-off time, / is the set of arrival/departure flights, ( , ) indicates the number of aircraft on edge at time , is the scheduling horizon. The objective function (1) aims to minimise the total taxi time. The predetermined runway scheduling results are complied with constraints (2) and (3). The safety distance constraints are satisfied with (4) and (5), through iteratively updating and checking the time windows sets. The decision variables are the routes to be allocated for each aircraft. Note the AGM model above is not a strictly MILP formulation. Hence, the decision variable is not binary; instead, it abstractly represents the route to be scheduled for each aircraft , including the route edges to be selected with required route start/end time ( )∕ ( ) for arrive/departure flight (see constraints (2) and (3)). Each route is constructed iteratively as part of Algorithm 1, rather than selected from predefined routes. Given deterministic taxi time for each route edge, constrained route start/end time ( )∕ ( ) for arrive/departure flight and scheduled route edges, we can then obtain the route taxi time ( ) for each aircraft.
It is worth noting that we have mathematically established the AGM optimisation problem as a graph-based sequential searching model. Nevertheless, the CCP extension as shown in the next section is general and in theory could be applied to not only sequential approaches. Other routing algorithms, such as the ones based on MILP formulations could also be used. However, such MILP based approaches which could consider more aircraft at the same time with uncertain taxi times and remain computationally inexpensive requires more investigation for future work.

CCP extension
In the above subsection, the taxi time on the taxiways (i.e., the edges in the graph) is assumed to be fixed, while the uncertainties of taxi time are inevitable in practice (Chen et al., 2011;Ravizza et al., 2013a). To address this issue, various methods have been developed to improve the taxi time prediction accuracy (Ravizza et al., 2014b;Lordan et al., 2016). Nevertheless, due to the uncertain nature of the taxi time, one can hardly estimate the exact time. Instead, it is more practical to present the uncertainties with a bounded range, which can be obtained through statistical data processing (Ravizza et al., 2014b) and has been applied in Brownlee et al. (2018a) as inputs for AGM optimisation. Therefore the fixed taxi time values for each edge in taxi time set are now replaced by a set of bounded ranges, i.e., = {[ 1 , 1 ], [ 2 , 2 ], … }. To properly incorporate the taxi time uncertainties into the AGM optimisation, the CCP, which ensures that the probability of meeting constraints is above a certain level, is introduced in this section. Thus, we extend the previous deterministic AGM model to a CCP model, where 1 − represents the predefined confidence level. The objective function (1) is modified and the chance constraint is added as where {⋅} denotes the event probability, and the taxi time ( ) for scheduled route could vary due to uncertain traversing time at each taxiway. is defined as the entire taxi time under certain confidence level, as explained over the remainder of this section. The CCP objective function (6) is to minimise the entire aircraft taxi time under predefined confidence level. The chance constraint (7) restricts the objective function value. The CCP model for AGM optimisation with taxi time uncertainties is then formulated as: minimise , subject to constraints (2) to (5) and (7).
The CCP model is advantageous for its flexible uncertainty level settings. When the confidence level 1 − is set as 1, the CCP model aims to find the routes with minimal taxi time upper bound over all possible probabilities. The idea of this setting is similar to Brownlee et al. (2018a), where the Fuzzy-QPPTW algorithm was developed to search minimal taxi time considering the worst scenarios. Although robust solutions can be provided in this setting, we will demonstrate that the robust solutions are over-conservative to some extent in Section 5. For the CCP model, the over-conservative property can be easily accommodated via lowering its confidence level. For instance, when 1 − equals 0.5, the CCP model would seek for the routes solutions that have a minimal upper bound of the taxi time in half of the entire possible scenarios. Note the scheduled aircraft routes can always ensure operational constraints (2) to (5) through holding on certain taxiways; different setting of only influences the calculation of in (7), which is not a hard operational constraint. Therefore, decision makers can balance the robustness and effectiveness of the aircraft routes through adjusting the predefined confidence level. More explanations and clarifications of the CCP uncertainty level can be found in Charnes and Cooper (1959), Luedtke and Ahmed (2008).
However, the proposed CCP model for AGM cannot be directly solved since we have difficulty in calculating the event probability in the chance constraint (7). A sample approximation approach, which generates a set of scenarios to present the taxi time uncertainties, is thus adopted to approximately calculate the chance constraint (Luedtke and Ahmed, 2008). Denote the scenarios as = { 0 , 1 , 2 , … }, where every ∈ corresponds to an undirected graph . For each edge in , its associated taxi time ranges in set are randomly sampled, resulting in a sampled taxi time set as = {̂1,̂2, … }, in whicĥ∈ [ , ]( = 1, 2, … ). Consequently, the bounded ranges of taxi time due to uncertainties have been substituted by the sampled values in different scenarios.
Given the predefined confidence level 1 − , the inequality ∑ ∈ ( ) ≤ in (7) needs to be checked in each scenario and is allowed to be infeasible at most | | ⋅ scenarios, where | | is the number of scenarios. The binary variable is introduced for each scenario ∈ : = 1 if the entire taxi time of allocated aircraft routes ∑ ∈ ( ) is larger than and = 0 otherwise. In doing so, the chance constraint (7) now can be replaced by the sample approximation constraints as X. Wang et al.
where ∑ ∈ ( ) denotes the entire taxi time under scenario and is a large enough number. When = 0, constraints (8) indicate that the entire taxi time must be lower than or equal to ; otherwise the entire taxi time is allowed to be larger than since the inequality always holds with the big . Constraints (9) ensure that the number of scenarios, whose entire taxi time is larger than , should be at most equal to | | ⋅ . Therefore, the CCP sample approximation model for AGM with taxi time uncertainties is to minimise , subject to constraints (2) to (5), (8) and (9).
The sample approximation constraints can also be analysed as follows: constraint (9) shows that we only need to consider part of the scenarios with the lower simulated taxi time ∑ ∈ ( ). The number of scenarios that not considered is determined by the entire sampled scenarios as well as the confidence level. Besides their corresponding variables are set as 1, so that constraints (8) for these scenarios are always satisfied, as is a large number. Therefore, when we calculate the objective function, only the remaining scenarios with = 0 are considered. Then in line with constraints (8), the objective function value equals the largest simulated taxi time ∑ ∈ ( ) among the remaining scenarios. Here we define the objective function of the CCP sample approximation model as the CCP taxi time, which depends on the predefined confidence level over the sampled scenarios. To clearly state the meaning of the CCP taxi time , a toy instance is listed in Table 2, where and ( ) denote the scenario index and aircraft taxi time under scenario respectively. As shown in the table, the taxi time varies over the ten sampled scenarios. When the confidence level 1 − is 1, the CCP taxi time is identified with the largest taxi time value, i.e., 20. While is 0.2, the worst two cases (10 * (1 − 0.8) = 2) with taxi time 20 and 17 (the numbers with underline) are not considered and the remaining largest taxi time 16 is the CCP taxi time. Similarly, the CCP taxi time equals 15 when = 0.5 (neglecting the scenarios with italic numbers).

Solutions
Building on the sample approximation CCP model, a CCP-QPPTW based framework is proposed to address the AGM optimisation problem with taxi time uncertainties. The structure of the proposed framework is introduced at first, and a modified sequential QPPTW algorithm (CCP-QPPTW) is developed in accordance with sampled multi-scenarios. To further improve the algorithm performance, a local heuristic is designed to identify the peak hours and randomly change the aircraft routing orders within the identified duration.

Framework
Algorithm 1 The overall CCP-QPPTW based framework. 1: Sort aircraft to natural ordering (arrival flights before departures, then by runway time) 2: for all aircraft ∈ do 3: Building on the sample approximation CCP model, allocate a route for aircraft using CCP-QPPTW (Algorithm 2) respecting already allocated other aircraft routes 4: end for 5: Identify peak hours set over the scheduling horizon 6: for all peak duration ∈ do 7: Remove allocated aircraft routes within duration 8: Reallocate aircraft routes within duration using heuristic (Algorithm 3) 9: end for The overall framework of our proposed method is illustrated in Algorithm 1. As the core of our proposed framework is CCP-QPPTW, which is a modified sequential QPPTW algorithm, the aircraft need to be sorted in line 1. Following Ravizza et al. (2014b), Brownlee et al. (2018a), a sequential loop in lines 2 to 4 for all aircraft is considered to allocate the its corresponding route, while respecting the routes that have been allocated to other aircraft. The time windows list for each edge needs to be updated at each iteration.
Once the routes have been allocated to each aircraft, a local heuristic method is invoked, in which the peak hours are identified initially. The detailed identification will be described in Section 4.3. Then within each identified peak duration, the already allocated aircraft routes are removed, and their orderings are iteratively resorted. With each aircraft orderings, the CCP-QPPTW is again applied for routes allocation. It is guaranteed that the taxi time of the finally reallocated routes is not longer than that of the original routes.
Moreover, a flowchart of the CCP-QPPTW based framework is illustrated in Fig. 2 consisting of three major components: aircraft sorting, iterative CCP-QPPTW method application (Algorithm 2), and a local heuristic to address peak hours (Algorithm 3).

Algorithm 2
The CCP-QPPTW algorithm to find the route for current aircraft.

CCP-QPPTW
To address the established CCP sample approximation model for AGM optimisation, one should carefully design a routing method over multiple sampled scenarios. The QPPTW method is originally developed for deterministic AGM optimisation problem with fixed taxi times. The general idea of QPPTW is to allocate the aircraft routes in sequences with given aircraft orderings. For each aircraft to be allocated, the QPPTW applies a labelling algorithm to gradually extend the path from the starting node, until the ending node is found. The generated route is ensured to be conflict-free through checking the available time windows set associated with each node, and the time windows are iteratively updated when current route is determined. For more details of QPPTW, the readers are referred to Ravizza et al. (2014a).
Since QPPTW can search optimal solutions with fixed taxi time in predefined aircraft orderings, it is promising to further extend QPPTW for our developed uncertain model. Moving QPPTW from single deterministic scenario to multiple sampled scenarios is the key motivation for the adaptation. A modified CCP-QPPTW algorithm is thus designed and the structure of CCP-QPPTW is provided in Algorithm 2.
The CCP-QPPTW algorithm consists of three parts: candidate routes generation (lines 2 to 5), routes selection (lines 6 to 12), and time windows update (lines 13 to 15). The selected route and associated taxi time are initially set as empty and a big number M, respectively. To generate possible routes for current aircraft, a conflict-free route is obtained under every scenario . Considering that the taxi time has been sampled as fixed value for each scenario, the QPPTW algorithm is then repeatedly used for the route allocation. Note that the QPPTW algorithm is replaceable by any other deterministic routing methods in the proposed framework. Given the same routes could be generated under different sampled scenarios using QPPTW, the generated routes under different scenarios can only be added into the routes set when the route is different from any existing routes in . In the next part, the recorded route with the minimal CCP taxi time would be selected for current aircraft. Through comparing the CCP taxi time, the aircraft route can be eventually determined. The time windows need to be updated in the third part of CCP-QPPTW. While QPPTW only updates the time windows in one deterministic graph, CCP-QPPTW executes update process over all sampled multi-scenarios, respectively.
X. Wang et al. Algorithm 3 The heuristic to reduce taxi time via adjusting aircraft orderings.

Local heuristic
The routing solutions have been obtained using the modified CCP-QPPTW algorithm, in which the aircraft are sorted with the natural orderings, i.e., arrival flights before departures and then by runway time. It is clear that different aircraft orderings could allocate different routes, resulting in different taxi times. Thus, a local heuristic is introduced in this section to further reduce the entire taxi time via adjusting the aircraft allocation orderings.
As stated in Algorithm 1, the peak duration is identified at first for the local heuristic, and the aircraft routes will then be reallocated within each peak duration. This identification process is motivated by the trade-off between the solution space (i.e., the number of aircraft to be routed at a time) and the algorithm efficiency; only within suitable solution subspace, the heuristic designed later will have better performance. Moreover, the peak duration could include more aircraft and its corresponding potential taxi time reduction could be larger than that within peak-off duration with less aircraft.
To initiate the identification process, a time slot unit, e.g., 15 mins, is set to divide the entire routing horizon. Next, the average number of aircraft taxiing within the time slot is calculated over all sampled scenarios, since the aircraft taxi time varies over multiscenarios due to uncertainties. The peak duration can now be defined as the divided duration whose average number of taxiing aircraft is greater than a predefined threshold.
Once the peak duration has been identified, a heuristic is proposed in Algorithm 3 to adjust the aircraft orderings within each duration, aiming to further reduce the aircraft taxi time. Denote the set of aircraft within current duration as . The taxi time of the entire aircraft set with natural orderings is first calculated to provide upper bound for the heuristic. For every attempt to change the aircraft orderings, the new orderings are determined by indicator , which depends on the average waiting time for current allocated route under all scenarios, and a random indicator , . The former indicator represents that the aircraft with longer waiting time is of high priority for the next round route allocation, and the latter value ensures the proposed heuristic is capable of jumping out of the local optimum. The random property of the heuristic is controlled by the equation is the waiting time of route , (Here the waiting time is referred to the time spent on the taxiways when the aircraft has to stop and wait to avoid conflicts with other aircraft; this includes the 2-minute waits for a landing or departing aircraft to clear the runway), and is an input parameter. The larger is, the high randomness of the heuristic has. After the aircraft within certain peak duration are resorted, the existing routes as well as occupied time windows are removed for the next step, i.e., aircraft routes reallocation.
If the resorted orderings have not been found before, the CCP-QPPTW algorithm is invoked with the new orderings, generating a set of new routes for the all aircraft in . If the reallocated routes have less entire taxi time than the upper bound taxi time , , the new routes are selected and the upper bound , is updated as well. Eventually, the heuristic returns the updated routes with the smallest taxi time. Note the number of attempts to resort the aircraft and generate new routes is limited by a fixed threshold number or a maximal heuristic running time.

Experimental setup
Manchester Airport (MAN), the third busiest airport in the UK, is selected to test the performance of our proposed model and algorithm. As illustrated in Fig. 3, MAN has two runways, three terminals and 148 aircraft gates. Typically MAN operates in two X. Wang et al.  working modes: for the single runway mode, both of the arrival and departure flights use runway 05L/23R; for twin runway mode, arrival aircraft land on 05L/23R and departures take off using 05R/23L. Note the twin mode only operates in daytime peak hours. To establish the undirected graph for AGM modelling, the layout of MAN is preprocessed with the GM Tools developed in our previous work (Brownlee et al., 2018a) (available at http://hdl.handle.net/1893/30962), and is constructed using 698 nodes and 739 edges. Moreover, a set of realistic aircraft movement data for one day operations is provided by MAN. We have also generated extra sets of movement data through replicating original data and modifying the scheduled runway time. In doing so, a set of scenarios with different numbers of aircraft movements have been generated, in which the traffic levels vary from low density (0.8 times of the original movements), the realistic density, to the high density (1.3 times). The number of aircraft movements for each density level is listed in Table 3. In order to address the taxi time uncertainties at MAN, an adaptive fuzzy rule-based system is introduced to describe the taxi time based on realistic aircraft movement data (Brownlee et al., 2018a). The upper/lower bounds of the taxi times = {[ 1 , 1 ], [ 2 , 2 ], … } for each edge ∈ can be determined using the historical data. Moreover, a nominal value within each taxi time range [ , ]( = 1, 2, … ) is obtained as well. Specifically, the collected taxi times range from 4.7 to 130.9 s, and the average nominal taxi time for each edge is 40.9 seconds. We will use these values as inputs for the simulations later.
To conduct the AGM simulations at MAN and test the performance of the proposed model and method, a simulation platform developed in Brownlee et al. (2018a) is adopted in this work. The motivation for introducing the simulation platform is to reflect airport operational conditions as realistically as possible. The platform is implemented in Java and follows several rules and simplifications. First, the simulated aircraft on the platform always respect the minimum safety distance, and is required to X. Wang et al.

Table 4
Comparisons between CCP-0.5-QPPTW and the exact method. CCP taxi time and run time are the average values over 100 simulations, gap denotes the CCP taxi time relative difference with the optimal solution (i.e., the CCP taxi time obtained by the exact method), and opt is the number that CCP-0.5-QPPTW achieves global optimal solutions over 100 random simulations. stop to avoid conflict with other aircraft. This ensures that scheduled routes always satisfy the operational constraints during the simulation. Second, the position of aircraft is updated iteratively with 0.1 s increment. Third, all edges crossing the runway are considered to be occupied for 2 mins prior to scheduled runway time for arrivals and 1 min after runway time for departures. Fourth, to present the taxi time uncertainties on edges, the taxi time for certain edge is randomly generated within the taxi time range, which is provided as inputs. Moreover, to test the algorithm performance on different situations, parameter ∈ [0, 0.5] is introduced to control the taxi time uncertainty levels. Specifically, the taxi time for each edge is uniformly sampled within the range [ − 2 ( − ), + 2 ( − )]( = 1, 2, … ). When equals 0.5, the taxi time has maximal uncertainty and the time is selected in [ , ]( = 1, 2, … ). Five scenarios with different uncertainty levels are thus generated via adjusting the value of ( = 0.1, 0.2, 0.3, 0.4, 0.5). For each simulation scenario, is assigned to all aircraft with the same value. All the simulations are conducted on Queen Mary's Apocrita high-performance computing facility (http://doi.org/10.5281/zenodo.438045). A more detailed description of the simulation tool, and validation of its underlying models, can be found in Brownlee et al. (2018a).

Results and discussions
The CCP based routing algorithms developed in this work are denoted as CCP--QPPTW ( = 0.5, 0.8, 1.0), where indicates the confidence level. Here the confidence level acts as a parameter for CCP-QPPTW, and different setting of only influences which routes would be scheduled, while all the aircraft operational constraints are satisfied in the simulations. We generated 30 scenarios to employ the CCP-QPPTW routing algorithm. Note these scenarios are different from those to evaluate the allocated routes in the comparison simulations. For each aircraft, we use QPPTW to generate a path with the quickest taxi time in the multi-scenarios respectively, and for each generated path, we further calculate its required taxi time for the route in every scenario. Then the candidate route with the least CCP taxi time is selected. The selection of 30 scenarios, i.e., samples, are further discussed and justified later.
To verify the effectiveness of the proposed CCP-QPPTW, we first compare the entire CCP taxi time obtained via CCP-QPPTW with an exact routing method. Note although the uncertain AGM optimisation problem has been transferred to a deterministic optimisation problem over multiple scenarios, the problem is still NP-hard (Ravizza et al., 2013b). Hence, the comparisons are conducted in a small instance with five flights. The exact method is designed based on CCP-QPPTW, which sequentially searches local optimal solutions for each aircraft. The CCP-QPPTW is iteratively applied with all 5! = 120 aircraft orderings, and the routes with the least CCP taxi times are selected as the optimal solution. Note in general, traversing all possible aircraft orderings with CCP-QPPTW still cannot ensure optimal solutions, since two aircraft may comprise with each other and detour to avoid waiting times, leading to a better solution. While in this five-flight small instance at Manchester Airport, we manually checked that aircraft detour cannot further generate better solutions. Therefore, an optimal solution has been obtained through iteratively using CCP-QPPTW. Specifically, the aircraft traffic level and taxi time uncertainty level are set as 1.0 and 0.3, respectively. The confidence level of CCP-QPPTW is set as 0.5, and the local heuristic in CCP-QPPTW is not activated for fair comparisons. Five flights are randomly selected within a 15-min peak duration.
The comparison experiments between CCP-0.5-QPPTW and the exact method are repeated for 100 times, and the results are reported in Table 4. Clearly CCP-0.5-QPPTW realises near-optimal solutions in terms of the CCP taxi time, which is the objective function value of the established AGM optimisation model. Moreover, it obtains optimal solutions in 26 out of 100 simulations. CCP-0.5-QPPTW also has a relatively low run time, which is 100 times faster than that of the exact method. Even though in a five-flight instance, the run time of the exact method is already over 200 s. It indicates that the exact method is not applicable to the large-scale AGM optimisation problem due to limitations of its computational efficiency.
To validate the performance of proposed CCP-QPPTW algorithms with three different confidence levels, two other routing methods, the deterministic QPPTW (Ravizza et al., 2014a) and Fuzzy-QPPTW (Brownlee et al., 2018a), are compared in different airport traffic densities and taxi time uncertainties. Note the local heuristic, which is part of the proposed CCP-QPPTW based framework, is not activated by default for fair comparisons between the CCP-QPPTW and other routing algorithms. For each combination of the traffic density and taxi time uncertainty, the allocated routes for each algorithm is repetitively simulated for 30 times (the same 30 seeds are generated for each routing algorithm). Therefore we entirely have 5 × 5 × 6 × 30 = 4500 runs (5 algorithms, 5 taxi time uncertainty levels, 6 traffic density levels, and 30 random runs) for the comparison simulations.
Four performance indicators are introduced to comprehensively evaluate the routing method performances. The first indicator is the total variations of simulated times on the minimum times. Here minimum times correspond to the time to complete the shortest route (allocated by Dijkstra's algorithm (Dijkstra, 1959)) without considering aircraft conflicts, which can be regarded as a lower bound for the taxi time. The second indicator is the total of variations of simulated times on routed times, in which routed times refer X. Wang et al. to the expected times following the allocated routes by each respective routing method. The number of aircraft delayed by more than 30 s over the routed times and the number of aircraft stops due to conflicts are also recorded as the performance indicators.
The comparison results are reported in Figs. 4 to 7. Overall when the confidence level for CCP model is set as 0.5, the CCP-0.5-QPPTW has evident advantages in terms of the total variations on minimum times in Fig. 4, number of aircraft delayed in Fig. 6 and number of aircraft stops in Fig. 7. When the confidence level of the CCP-QPPTW increases, it means that we expect to allocate more conservative routes under taxi time uncertainties, so that the aircraft would detour to avoid potential conflicts even though there is a relatively lower possibility for conflicts. Consequently, the detour could reduce the simulated taxi time under the scenario where the conflict would indeed happen if the detour is not executed. However, for the majority of the simulated scenarios, the simulated taxi time would increase as the conflict does not occur due to a low possibility. Therefore the CCP-QPPTW algorithm degrades with the increase of the CCP confidence level. Note the results of CCP-0.5-QPPTW have worse exceptions under a lower uncertainty level (see Figs. 4(f) and 6(f) with 0.3 uncertainty level), where the taxi time uncertainty ranges are relatively small. This is because the results of CCP-0.8/1.0-QPPTW are not that conservative and may provide slightly better results with lower uncertainty levels. While the uncertainty level is greater than 0.3, the results of CCP-0.5-QPPTW are always better than that of CCP-0.8/1.0-QPPTW in terms of the total taxi times and number of stops.
Since QPPTW is originally designed for the deterministic routing problem, it is unsurprising to see QPPTW has comparative and even best performance in terms of all indicators with lower densities and taxi time uncertainty levels. However, when these two parameters increase, other methods that can address the uncertainty have clear advantages over QPPTW. In contrast, Fuzzy-QPPTW has a relatively poor performance with lower densities and taxi time uncertainty levels, especially in terms of the total variations on minimum times and number of aircraft stops. This is logical; Fuzzy-QPPTW considers the worst possibilities when selecting the routing path, therefore leading to over-conservative results with less complicated environments. However, even if the airport considers more aircraft and has larger taxi time uncertainties, CCP-0.5-QPPTW has overall better performance than Fuzzy-QPPTW, especially for the total variations on minimum times and number of aircraft stops. Note that Fuzzy-QPPTW indeed has slightly better performance when considering the number of aircraft delayed. This is because Fuzzy-QPPTW partially aims to reduce actual taxi time X. Wang et al.  variations on predicted routed times; it has considered possible delays to the routed time due to conflicts, and consequently a larger taxi time including stops is predicted to ensure the route has fewer number of delayed aircraft. Therefore it is understandable that allocated routes of Fuzzy-QPPTW has the least number of delayed aircraft compared to other routing methods. This observation is also consistent with the results for the total variation times in Figs. 4 and 5. As shown in Fig. 4(f), the average simulated routing taxi times of CCP-0.5-QPPTW is less than that of Fuzzy-QPPTW; while in Fig. 5(f), the Fuzzy-QPPTW has evident less variations on routed times, i.e., the predicted routing times, indicating the routed times of Fuzzy-QPPTW is clearly larger than that of CCP-0.5-QPPTW. These longer routes for Fuzzy-QPPTW are, ironically, also more prone to conflicts that only add a small delay, so the aircraft stops more often, albeit the stops are very short. This is clearly undesirable from a fuel efficiency point of view.
To further compare the routing methods on AGM optimisation with taxi time uncertainties, the average simulated taxi times and the total taxi distances for allocated routes in different traffic densities are listed in Tables 5 and 6. As shown in Table 5, CCP-0.5-QPPTW achieves 0.92% taxi time decrease on average compared to the deterministic QPPTW method, and the performance X. Wang et al. Fig. 6. Number of aircraft delayed (per the simulations) by more than 30 s over the routed times as uncertainty increases at each traffic density. Each group of five bars has the same level of uncertainty. of CCP-QPPTW gradually becomes worse when the confidence level increases. As for Fuzzy-QPPTW, although it shortens the taxi time compared to QPPTW under 1.3 traffic density, its average taxi time over all traffic densities has a 0.40% increase. The overconservative strategy applied in Fuzzy-QPPTW can be also illustrated in terms of the total taxi distances in Table 6. Compared to QPPTW, both of the CCP-QPPTW and Fuzzy-QPPTW methods allocate longer routes to address the taxi time uncertainties. This is because QPPTW aims to minimise the deterministic total taxi times, corresponding to the shortest taxi distances. However, these routes may not necessarily lead to the minimum simulated total taxi times when uncertainties are introduced (it indicates that simulated actual taxi time over one edge could change). Under such situations, although the allocated taxi distance remains, the allocated routes could experience more conflicts and have to stop and wait, leading to longer simulated total taxi times. On the other hand, the Fuzzy-QPPTW and CCP-QPPTW methods have considered the uncertainties in advance, and could detour to avoid potential conflicts in advance, thus leading to longer taxi distances and reduced taxi times. However, the Fuzzy-QPPTW increases the taxi distances over 1%, while CCP-0.5-QPPTW chooses less conservative routes with an average 0.36% distance increase. Meanwhile, X. Wang et al.  the conservative degree of the CCP-QPPTW method can be adjusted by improving the confidence level, and the routes distance difference with QPPTW would increase from 0.36% to 0.60%. To compare the computational efficiency of various routing methods, we report the average run time for each traffic density over five allocation methods in Table 7. Clearly QPPTW has the shortest run time since it does not address the taxi time uncertainties. On the other hand, Fuzzy-QPPTW generates taxi time fuzzy functions for each segment, and the functions are further combined for fuzzy operations, leading to a very large number of object creation and destruction operations and their associated overhead. This means that Fuzzy-QPPTW has the worst computational efficiency for the routes allocation. While CCP-QPPTW has better performance than Fuzzy-QPPTW, its average run time is much less than the Fuzzy-QPPTW computational time, and the variations of the confidence level do not heavily influence the run time. In conclusion, the proposed QPPTW-CCP method can better address the taxi time uncertainties with higher computational efficiency in comparison with Fuzzy-QPPTW.
The number of samples (i.e., scenarios) may also impact the accuracy and efficiency of the proposed framework. Ideally a larger sampling number better represents the approximation accuracy, while we have to strike a balance between the approximation and X. Wang et al. computational efficiency. Here, we set the number of samples as 30 in line with a rule of thumb (Hogg et al., 1977). To further justify the selection of the sampling number and demonstrate the approximation accuracy, we fix the traffic density and uncertainty level as 1.0 and 0.3 respectively, and simulated 1000 scenarios using 1000 different random seeds. The density distributions of simulated total taxi times and four performance indicators are illustrated in Fig. 8(a) to (e). As shown in the figure, the simulated results between 30 and 1000 samples (i.e., scenarios) have similar density distributions in terms of the taxi times and four performance indicators. We also compared the mean values of these five indicators over 30 and 1000 samples, and their differences are 0.07%, 2.59%, 2.80%, 0.59% and 1.93%, respectively. Moreover, the selection of 30 samples is justified from the statistic view. We introduce the two-sample Kolmogorov-Smirnov test (Massey, 1951;Xiao, 2017), which is one of the most useful and general non-parametric statistic methods, to check whether the distribution functions between two sample groups differ. For each indicator, the null hypothesis is that the two samples between 30 and 1000 simulated scenarios are from the same distributions. Using the two-sample Kolmogorov-Smirnov test, the -value for each indicator is 0.8106, 0.2102, 0.4220, 0.9219, 0.7932, respectively. All of them are X. Wang et al.  greater than the 0.05 significance level, thus the null hypothesis cannot be rejected. Therefore 30 samples are generated to represent the scenario uncertainties.
Finally we consider the effect of the local heuristic in the proposed routing framework. The parameters of the local heuristic are given as follows: time slot unit 15 mins, = 2, orderings change attempt 10, maximal heuristic run time 500 s. Preliminary simulations indicate that the local heuristic performance does not significantly vary from different parameter settings, therefore we only report the results in default parameters. The comparison results with and without heuristics are illustrated in Figs. 9 and 10, where the confidence level is set as 0.5 and CCP-QPPTW-H/CCP-QPPTW denotes for the CCP-QPPTW method with/without the local heuristic. Note the local heuristic reaches the maximal run time for each case, and the allocated routes are repetitively simulated for 30 times.
As seen in Fig. 9 (density 1.0) and Fig. 10 (density 1.3), CCP-QPPTW with local heuristic does not effectively increase the performance in terms of all indicators. Especially for the simulations with lower uncertainties, the local heuristic even leads to worse results compared to the CCP-QPPTW without local heuristic. It is true the CCP-QPPTW-H can realise better results with the additional heuristic over the same scenarios applied to select aircraft routes by the CCP-QPPTW method. However, different scenarios are simulated to evaluate the allocated aircraft routes, which could lead to increased total taxi times. Although applying more scenarios to the CCP-QPPTW/CCP-QPPTW-H could better approximate the uncertainties and improve the robustness of the proposed algorithm, the allocated aircraft routes could still have slightly worse performances under simulated scenarios. We recommend to employ the local search heuristic within a short scheduling horizon. Typically, the taxi times can be estimated more accurately within a shorter horizon, so that the updated routes using the local heuristic have a higher possibility to reduce the total taxi times. Overall, the proposed local heuristic, which resorts the aircraft orderings, has limited contributions to address the uncertainties and reduce the taxi time. This finding is also consistent with the results in Ravizza et al. (2014a), Brownlee et al. (2018b) that the heuristics to reorder the aircraft do not have huge impact on the deterministic AGM optimisation.

Conclusions
To address the aircraft taxi time uncertainties, a chance-constrained programming (CCP) model for airport ground movement (AGM) optimisation has been developed in this work. Knowing that the taxi time is typically difficult to be exactly predicted, we generate a set of scenarios to approximately sample the taxi time distributions, in which the uncertainties can be well represented. Building upon the sampled multi-scenarios, a modified sequential quickest routing method with a local heuristic (CCP-QPPTW) is then designed to solve the proposed model. In doing so, the scheduled aircraft routes can be generated considering different conservative degrees, i.e., the confidence level in the CCP-QPPTW algorithm, thus avoiding over-conservative routes and reducing total taxi times. Besides, the proposed method is computationally efficient compared to the state-of-the-art algorithm, e.g., the Fuzzy-QPPTW. Using the benchmark at Manchester Airport, extensive simulations are conducted to compare the proposed CCP-QPPTW routing method with the-state-of-art routing methods. The results indicate that CCP-QPPTW with confidence level 0.5 outperforms other algorithms in terms of the average simulated taxi times. Besides, significantly fewer number of aircraft stops has been achieved by the CCP-QPPTW, implying that the allocated routes would probably save much fuel consumption at the same time as stop-and-go has been identified as one of the main sources towards increased fuel burn (Khadilkar and Balakrishnan, 2012). However, the local heuristic which can resort the aircraft orderings has limited contributions to reduce the taxi time.
Future research on AGM optimisation with taxi time uncertainties may be oriented in three directions. First, the objective function in this study minimises the total taxi time, while other objectives, e.g., fuel consumption and emissions, are also important. Therefore, multi-objective optimisation for uncertain AGM may be taken into consideration for the future research. Second, to further address the uncertain AGM optimisation, it is promising to combine the proposed CCP-QPPTW and reactive real-time routing methods (Evertse and Visser, 2017). Third, all aircraft have successfully scheduled routes respecting operational constraints in our simulations, while at certain extreme peak hours in practice, some of the aircraft may not be able to allocate routes in time. Thus how to schedule aircraft taxiing routes to minimise the delay at peak hours deserves to be further investigated.