Improvement of Airport Surface Operation at Tokyo International Airport Using Optimization Approach

: Congestion and delays occur on airport surfaces as a result of a rapid increase in the demand for air transport. The aim of this study is to determine the differences between optimized and observed operations to improve airport surface operation at Tokyo International Airport by using mixed-integer linear programming to minimize the total ground movement distance and time based on real-time ﬂight information. Receding horizon schemes are considered to adapt to dynamic environments. The model obtains results that reduce the taxi distance by 18.54% and taxi time by 29.77% compared with the observed data. A comparison of taxiway usage patterns between the optimization results and observed data provides insight into the optimization process, for example, changes in runway cross strategies and taxiway direction rules. Factors such as the objective function weights and airline–terminal relationship were found to signiﬁcantly affect the optimization result. This study suggests improvements that can be made at airports to achieve a more efﬁcient surface operation.


Introduction
An increasing number of hub airports have reached their capacity owing to a rapid increase in the demand for air transport, resulting in congestion and delays on the airport surface, which can cause a series of major problems [1]. Although the number of air transport passengers has reduced during the COVID-19 pandemic, it is expected to recover by 2023 and continue increasing [2]. The Tokyo International Airport, commonly known as Haneda Airport, is one of the two primary airports that serve the Greater Tokyo Metropolitan Area. It handled 85.51 million passengers in 2019 and ranked fifth among all airports worldwide [3].
Two methods can be implemented to improve airport capacity and efficiency. One method is the expansion of the current infrastructure of airports. For example, Tokyo International Airport started operating its fourth runway in 2010 [4]. The newly built runway, taxiway, and apron area can improve airport capacity. The other method is air traffic control optimization [5], which involves a service controlling aircraft on the ground and in the airspace. This method can improve the surface operation efficiency of the airport without the addition of new airport facilities, which is more scalable and economical compared with the first method. The main operations on the airport surface include ground movement, runway management, and gate assignment [6]. These three operations are correlated and can affect each other. Indicators such as the total ground movement distance and time are affected by these operations, and improving these indicators can address the aforementioned problems [7].
Several papers have been published about airport surface operations, where ground movement is the primary area of interest. To optimize aircraft taxi times and distances on the ground, the mixed integer linear programming (MILP) formulation and solver approach is widely used [8]. Gotteland and Durand [5] and Pesic et al. [9] have implemented genetic algorithms. Ravizza [10] proposed a sequential graph-based algorithm to address the ground movement problem more realistically. The estimation of taxi time is also important in this problem. Brownlee [11] proposed an adaptive Mamdani fuzzy rule-based system to estimate taxi times. Other optimization objectives such as fuel consumption and economic implications are also considered in the multi-objective approach [12]. Some studies have considered integrating other operations with ground movement. Clare and Richards [13] considered the integrated optimization of taxiway routing and departure sequencing. Behrends and Usher [14] integrated the aircraft ground movement problem with the gate assignment problem, and Yu et al. [15] solved an integrated gate re-assignment and taxi scheduling problem.
Past studies focused on proposing an optimization method and validating the output results or execution time with the existing approach. They lacked observed data as input and validated the optimization results with observed data to evaluate and understand the optimization process. The impact of runway and terminal selection and management on taxi distance and time was not considered as a whole. Some studies ignored and simplified details on the airport surface, such as the apron area and the runway cross area. In this study, we apply an optimization approach for airport surface operations to minimize the total ground movement distance and time at Tokyo International Airport.
We apply MILP optimization with receding horizon on a real airport using secondslevel observed tracking data (CARATS Open Data [16]) as input. The advantage is that it has rich aircraft taxiing details, such as push-back and taxi route, pause time, and runway selection, with which we can better apply dynamic optimization and, more importantly, all aircraft movements in the optimization results are compared with the observed operations in CARATS Open Data to determine the differences between them, to understand the optimization process and improve airport surface operations.
In addition, the effects of runway selection and gate/terminal assignment on aircraft taxiing are also considered together, in the application of the optimization model. The model emulates the actual operation rules of the airport more closely than in previous studies; it includes details such as aircraft push-back area, multi-runway operations, and runway crossing rules at the Tokyo International Airport. Factors such as objective function weights and airline-terminal relationships that affect the optimization results are also discussed.
In Section 2, the current surface operation at the Tokyo International Airport is illustrated for finding the rule and patterns from observed data. Model development such as formulations of the optimization model are explained in Section 3. Differences between optimized and observed operations and factors that affect the optimization results are explained in Section 4. Finally, the conclusions are presented in Section 5. Figure 1 shows the runway and terminal map of the Tokyo International Airport. It comprises four runways-16L-34R, 16R-34L, 04-22, and 05-23-which are used for takeoff or landing depending on the wind direction. Basically, landing aircraft use runways 22 and 23, while take-off aircraft use runways 16L and 16R during south wind operation. However, under certain conditions, landing aircraft may use runway 16L. During north wind operation, landing aircraft use runways 34L and 34R, while take-off aircraft use runways 05 and 34R. Under certain conditions, runways 34L and 04 can be used by take-off aircraft during north wind operation. The Tokyo International Airport has three passenger terminals: Terminal 1 (used by Japan Airlines (JAL), Skymark Airlines and Star Flyer) and Terminal 2 (used by All Nippon Airways (ANA), Air Do and Solaseed Air), which serve domestic flights; and the International Terminal (renamed as Terminal 3 on 14 March 2020), which serves international flights.

Overview of the Tokyo International Airport
The rules of runway and terminal operation in the Tokyo International Airport were in effect in the period before 2020, which is considered in this study. To simplify the model development and implementation, south wind operation is selected in the analysis because runway 04-22 was not used during north wind operation. Late night and early morning hours are not considered since it has a different rule of runway operation under the non-congested circumstance.

CARATS Open Data
The Collaborative Actions for Renovation of Air Traffic Systems [16] is provided by the Ministry of Land, Infrastructure, Transport and Tourism in Japan as CARATS Open Data, which is in CSV format and includes the GPS ground tracking data (in seconds) of aircraft moving on the surface of Tokyo International Airport. Every second data contain the timestamp, flight ID, latitude, longitude, altitude, and aircraft type. The following information can be extracted from the CARATS Open Data for each aircraft: the flight ID, start time, end time, take-off or landing, runway ID, gate ID, and taxiway route. This information plays two roles in this study. First, it is used to identify the current surface operation pattern in this section. It is also used to evaluate and analyze the results of our optimization model in Section 4.
The flight ID, start time, and end time can be directly obtained from the CARATS Open Data. Regarding the taxiway route, every taxiway used by an aircraft needs to be identified. Four points are used to determine one taxiway or runway. These four points should form a convex polygon, and if the GPS tracking point of an aircraft enters this convex polygon, it is considered that the aircraft has passed through this taxiway or runway. One point is used to determine each gate. Using this point as the center, a circle is created for each gate with a radius of 15 m. If an aircraft appears in the circle for more than 20 s, the gate corresponding to this circle can be considered as the gate of the aircraft. The judgment of the aircraft take-off or landing is determined based on the speed. If the speed in the first 5 s of the tracking data of an aircraft exceeds 70 knots, the aircraft is considered as a landing aircraft. Otherwise, it is considered as a take-off aircraft. The Tokyo International Airport has three passenger terminals: Terminal 1 (used by Japan Airlines (JAL), Skymark Airlines and Star Flyer) and Terminal 2 (used by All Nippon Airways (ANA), Air Do and Solaseed Air), which serve domestic flights; and the International Terminal (renamed as Terminal 3 on 14 March 2020), which serves international flights.
The rules of runway and terminal operation in the Tokyo International Airport were in effect in the period before 2020, which is considered in this study. To simplify the model development and implementation, south wind operation is selected in the analysis because runway 04-22 was not used during north wind operation. Late night and early morning hours are not considered since it has a different rule of runway operation under the non-congested circumstance.

CARATS Open Data
The Collaborative Actions for Renovation of Air Traffic Systems [16] is provided by the Ministry of Land, Infrastructure, Transport and Tourism in Japan as CARATS Open Data, which is in CSV format and includes the GPS ground tracking data (in seconds) of aircraft moving on the surface of Tokyo International Airport. Every second data contain the timestamp, flight ID, latitude, longitude, altitude, and aircraft type. The following information can be extracted from the CARATS Open Data for each aircraft: the flight ID, start time, end time, take-off or landing, runway ID, gate ID, and taxiway route. This information plays two roles in this study. First, it is used to identify the current surface operation pattern in this section. It is also used to evaluate and analyze the results of our optimization model in Section 4.
The flight ID, start time, and end time can be directly obtained from the CARATS Open Data. Regarding the taxiway route, every taxiway used by an aircraft needs to be identified. Four points are used to determine one taxiway or runway. These four points should form a convex polygon, and if the GPS tracking point of an aircraft enters this convex polygon, it is considered that the aircraft has passed through this taxiway or runway. One point is used to determine each gate. Using this point as the center, a circle is created for each gate with a radius of 15 m. If an aircraft appears in the circle for more than 20 s, the gate corresponding to this circle can be considered as the gate of the aircraft. The judgment of the aircraft take-off or landing is determined based on the speed. If the speed in the first 5 s of the tracking data of an aircraft exceeds 70 knots, the aircraft is considered as a landing aircraft. Otherwise, it is considered as a take-off aircraft.  Figure 2 shows the regular taxiway selection patterns in the observed data obtained from Tokyo International Airport. Patterns of taxiway selection indicate the most used taxiway routes, from the gate to the runway and vice versa, selected by air traffic control. The taxiway route of each aircraft is extracted from the CARATS Open Data, from which the taxiway selection pattern is obtained by finding the common sub-taxiway route of all aircraft taxiway routes. It is observed that there is usually a pattern when aircraft moves from different terminals to different runways during operation. When there are two parallel taxiways in some areas, they are used differently based on the taxi direction or destination. Considering the take-off flight in Figure 2a as an example, when an aircraft from Terminal 1 is required to take off from runway 16R, it is likely to use the red taxiway route in the figure.

Taxiway Selection Patterns
When it is required to take off from runway 16L, it always uses the blue taxiway route.  Figure 2 shows the regular taxiway selection patterns in the observed data obtained from Tokyo International Airport. Patterns of taxiway selection indicate the most used taxiway routes, from the gate to the runway and vice versa, selected by air traffic control. The taxiway route of each aircraft is extracted from the CARATS Open Data, from which the taxiway selection pattern is obtained by finding the common sub-taxiway route of all aircraft taxiway routes. It is observed that there is usually a pattern when aircraft moves from different terminals to different runways during operation. When there are two parallel taxiways in some areas, they are used differently based on the taxi direction or destination. Considering the take-off flight in Figure 2a as an example, when an aircraft from Terminal 1 is required to take off from runway 16R, it is likely to use the red taxiway route in the figure. When it is required to take off from runway 16L, it always uses the blue taxiway route.

Model Development
This section describes a dynamic optimization model that is used for airport surface operations. The model dynamically outputs the ground movement route and schedule for each aircraft based on real-time flight information. Additionally, it can provide a runway choice using the runway decider module. The model aims to reduce the total ground movement distance and time of all aircraft. All mathematical notations are explained in Table 1.

Model Development
This section describes a dynamic optimization model that is used for airport surface operations. The model dynamically outputs the ground movement route and schedule for each aircraft based on real-time flight information. Additionally, it can provide a runway choice using the runway decider module. The model aims to reduce the total ground movement distance and time of all aircraft. All mathematical notations are explained in Table 1. Table 1. Definition of mathematical notation in Section 3.

Symbol Description
N Set of nodes in graph model.

N t
Set of taxiway intersection nodes in graph model.

N entrance
Set of runway entrance nodes in graph model.

N exit
Set of runway exit nodes in graph model.  Identified nodes in graph model.
Length of shortest path between n 1 and n 2 . g(n 1 , n 2 ) Graph matrix of graph model. v_max(n 1 , n 2 ) Maximum taxi speed of taxiway between n 1 and n 2 . t separation_node (n) Time separation requirement at node n. t separation_cross Time separation requirement of runway cross node.
t Destination node of aircraft a. t start (a) Start time of aircraft a in current planning horizon. n last (a) Last node passed by aircraft a before current planning horizon. n start (a) Node ahead of aircraft a in current planning horizon. L rest (a) Rest distance of aircraft a to n start (a) in current planning horizon. K Number of steps that can be planned in one horizon. T Duration of the planning horizon. k, k 1 , k 2 Step number. x(a, k, n 1 , n 2 ) Decision Variable. x = 1 means aircraft a taxis from n 1 to n 2 at k th step. t(a, k) Decision Variable. The time when aircraft a passing k th step's node n 1 .
The estimated taxi time from n 1 to n 2 . t last (n 1 ) Time when aircraft last passed node n 1 in the previous planning horizon. gap(a 1 , a 2 , k 1 , k 2 ) Time gap between aircraft a 1 at step k 1 and a 2 at step k 2 . con f lict(a 1 , a 2 , k 1 , k 2 , n) Equals 1 when step k 1 of aircraft a 1 and step k 2 of aircraft a 2 use the same node, n. f irst(a 1 , a 2 , k 1 , k 2 ) Equals 1 when the k 1 step of aircraft a 1 is performed earlier than the k 2 step of aircraft a 2 .
cross(a 1 , a 2 , k 1 , k 2 , c) Equals 1 if aircraft a 1 crosses the runway cross point c at step k 1 when aircraft a 2 takes off or lands on the same runway at step k 2 . f irst takeo f f (a 1 , a 2 ) Equals 1 when a 1 takes off before a 2 .

Mixed-Integer Linear Programming and Receding Horizon
This study adopts the mixed-integer linear programming (MILP) optimization, which is applicable to hybrid problems. However, the typical MILP method faces two main drawbacks in the optimization of the current research problem. First, the optimization of the airport surface operation is a type of dynamic traffic assignment problem, in which multiple coupled non-convex decisions of routing and timing are simultaneously optimized; this can be considered as an NP-hard problem. Applying a global MILP solver to handle all aircraft in one day or one hour can considerably increase the time complexity, which decelerates the solving process and inhibits the application of the approach to practical operation. Second, the optimization of airport surface operations is a dynamic problem. The prediction of when the aircraft start and arrive may be inaccurate, the predicted positions for the current aircraft may be wrong, or the predictions of when new aircraft are ready to push back from the gates or land may be wrong. Furthermore, there is a possibility that an aircraft cannot arrive, push back, or taxi as planned if the plan is determined too early, and unexpected events may occur, which can result in unsuitability of the previous plan and replanning.
A receding horizon (RH) approach, which has been extensively used in trajectory planning, is implemented to solve the above problems by dividing the MILP into several planning horizons that are sequentially optimized in the MILP solver [18,19]. A basic RH framework was first presented by Schouwenaars et al. [19]; it has been applied to the Aerospace 2022, 9, 145 6 of 21 airport surface operation problem in some studies [8,13] and is adopted as the strategy of MILP optimization herein.
MILP with RH performs better for the aforementioned problem compared with the traditional MILP approach. In this approach, the global result is not obtained through a single MILP optimization. Contrarily, the planning process is divided into several horizons, similar to an iteration loop. Applying RH in MILP ensures that the computation is spread over time because of the nonlinear increase in the solving time with the problem size. On the other hand, the optimization approach could consider dynamic characteristics of airport operations, thus ensuring that the computation is not wasted on detailed plans for the distant future, which may be revised before the aircraft arrives at its destination.
The model plans for up to K future steps for each aircraft at each planning horizon. After a planning horizon, the aircraft movement is simulated based on the previously established plan. The final position and status of each aircraft are considered as the input of the next planning horizon. The duration of the planning horizon follows an already set horizon time length, T (s). The MILP process is as follows: plan K steps, simulate and calculate for T seconds, then plan the next K steps based on the last results, and simulate T seconds, until every aircraft arrives at its destination.

Airport Graph Model
A graph, G(N, E), of the Tokyo International Airport was constructed to help in the development of the MILP optimization model. N represents the node set, and E represents the directed edge set of the graph. N = n n ∈ N t or n ∈ N entrance or n ∈ N exit or n ∈ N g (1) ∀n 1 , n 2 ∈ N, if there is a taxiway can be used from n 1 to n 2 : All the taxiway intersections N t , runway entrances N entrance , runway exits N exit , and gates N g are considered as nodes. The pairs (n 1 , n 2 ) in the set, N, are considered as directed edges in E if there is a taxiway connection from n 1 to n 2 . An overview of the graph model (excluding the gate nodes) developed after these steps is shown in Figure 3. , is built to record the shortest distance between nodes, which is an | | × | | matrix calculated using the Floyd algorithm [20] shown in Algorithm 1: An adjacency matrix, g, is built to construct the model. g is an |N| × |N| 0-1 matrix in which g(n 1 , n 2 ) = 1 if and only if (n 1 , n 2 ) ∈ E. A shortest-path matrix, L SP , is built to record the shortest distance between nodes, which is an |N| × |N| matrix calculated using the Floyd algorithm [20] shown in Algorithm 1: Additionally, for each edge (n 1 , n 2 ), a parameter v_max(n 1 , n 2 ) is added to represent the maximum taxi speed on that taxiway, which is decided according to the comments of operational experts from Japan Airlines. A speed of 10.29 m/s (20 knots) is set for a regular taxiway, 5.14 m/s (10 knots) is set for the taxiway near the apron area, and 25.72 m/s (50 knots) is set for the rapid-exit taxiway. For each node, n, t separation_node (n) is set to 20 s to represent the time separation requirement at each node (see Section 3.5.2), and t separation_cross (n) is set to 20 s to represent the time separation requirement for the second type of runway conflict (see Section 3.5.3).

Virtual Pushback Node
A virtual pushback node system is developed for the aircraft to find its start or end nodes in the optimization model. In a few previous studies, the departure aircraft of a specific terminal usually starts at the same node in the graph model. The scale of the model can be reduced by using these merged nodes. However, the taxiway usage in the area near the gate and the specific situation when multiple aircraft are pushed back cannot be simulated. In order to overcome these drawbacks, this study proposes a pushback node system as the terminal and apron areas of the Tokyo International Airport are closely related to the taxi area. A new virtual pushback node is inserted in the middle of each taxiway that is near the terminal area and is connected to the gates. It divides the original taxiway into two sub-taxiways. This node is assigned as the pushback point to the gates connected to the specific taxiway. The departure and arrival aircraft using these gates consider these pushback nodes as their origins and destinations, respectively. Using Figure 4 as an example, a new node X is added between node 60 and node 77, which divides edge (60, 77) into two edges: (60, X) and (X, 77). Node X is set as the pushback node for nearby gates 56 and 57. All the virtual pushback nodes, N p , are also considered as nodes in N. Concurrently, the edge set, E, is adjusted accordingly depending on the insertion of these nodes and the division of the edges. to the taxi area. A new virtual pushback node is inserted in the middle of each taxiway that is near the terminal area and is connected to the gates. It divides the original taxiway into two sub-taxiways. This node is assigned as the pushback point to the gates connected to the specific taxiway. The departure and arrival aircraft using these gates consider these pushback nodes as their origins and destinations, respectively. Using Figure 4 as an example, a new node X is added between node 60 and node 77, which divides edge (60, 77) into two edges: (60, X) and (X, 77). Node X is set as the pushback node for nearby gates 56 and 57. All the virtual pushback nodes, , are also considered as nodes in . Concurrently, the edge set, , is adjusted accordingly depending on the insertion of these nodes and the division of the edges.

Runway Virtual Node and Cross Node
As the take-off, landing, and sequencing processes are considered in the model, a corresponding node system is designed in the graph model. Two virtual nodes, and , are added for each runway, and they are connected by a virtual directed edge ( , ).

Runway Virtual Node and Cross Node
As the take-off, landing, and sequencing processes are considered in the model, a corresponding node system is designed in the graph model. Two virtual nodes, r 1 and r 2 , are added for each runway, and they are connected by a virtual directed edge ( r 1 , r 2 ).
Assuming Figure 5 represents a take-off runway, virtual edges (1, r 1 ) and (2, r 1 ) are added to connect the runway entrance node to the runway virtual node. Aircraft enter the runway through the runway entrance node and line up at node r 1 ; they take-off at the edge ( r 1 , r 2 ) and arrive at their destinations at node r 2 in the optimization model. Assuming Figure 5 represents a landing runway, virtual edge ( r 2 , 3) is added to connect the runway virtual node to the runway exit node. Landing aircraft start at node r 1 in the optimization model, land on edge ( r 1 , r 2 ), exit the runway at node r 2 , and taxi to runway exit node 3.

Runway Virtual Node and Cross Node
As the take-off, landing, and sequencing processes are considered in the model, a corresponding node system is designed in the graph model. Two virtual nodes, and , are added for each runway, and they are connected by a virtual directed edge ( , ).
Assuming Figure 5 represents a take-off runway, virtual edges (1, ) and (2, ) are added to connect the runway entrance node to the runway virtual node. Aircraft enter the runway through the runway entrance node and line up at node ; they take-off at the edge ( , ) and arrive at their destinations at node in the optimization model. Assuming Figure 5 represents a landing runway, virtual edge ( , 3) is added to connect the runway virtual node to the runway exit node. Landing aircraft start at node in the optimization model, land on edge ( , ), exit the runway at node , and taxi to runway exit node 3. In the Tokyo International Airport, there are runway cross points on runway 16R-34L and runway 04-22, where aircraft can taxi across the runway. To design the constraints to avoid cross conflict in Section 3.5.3, a runway cross node system as that shown in Figure  5 is designed. For each cross node (node 4 in Figure 5) in runway cross node set , , (node and node in Figure 5) represent the two runway virtual nodes used for take-off or landing on the runway, while , (node 5 and node 6 in Figure 5) represent the two nodes of the taxiway that cross the runway. In the Tokyo International Airport, there are runway cross points on runway 16R-34L and runway 04-22, where aircraft can taxi across the runway. To design the constraints to avoid cross conflict in Section 3.5.3, a runway cross node system as that shown in Figure 5 is designed. For each cross node c (node 4 in Figure 5) in runway cross node set C, c r1 , c r2 (node r 1 and node r 2 in Figure 5) represent the two runway virtual nodes used for take-off or landing on the runway, while c t1 , c t2 (node 5 and node 6 in Figure 5) represent the two nodes of the taxiway that cross the runway.
All the runway virtual nodes and cross nodes are considered as nodes in N. Concurrently, the edge set, E, is adjusted accordingly depending on the insertion of these nodes and edges.

Aircraft Data
For the entire optimization model, set F includes all the aircraft considered in the optimization. Each aircraft a in F has some static parameters given as follows: gate(a) represents the gate node of aircraft a; start_time(a) represents the pushback ready time of a departure aircraft and the time of arrival at the runway of an arrival aircraft; rwy(a) represents the runway of aircraft a; type(a) represents the aircraft type of aircraft a; and destination(a) represents the destination node of aircraft a in the optimization model. For a take-off aircraft, destination(a) is the runway virtual node r 2 of rwy(a) (see Section 3.2.3), and for a landing aircraft, destination(a) is the corresponding virtual pushback node of gate(a) (see Section 3.2.2).
In a specific planning horizon, the active aircraft set A includes the aircraft that has a start time later than the time of the current horizon and has not reached its destination. For each aircraft a in A, t start (a), n last (a), n start (a), and L rest (a) are used to define the start time and position in the current planning horizon. At time t start (a), aircraft a taxis from n last (a) to n start (a), at a distance, L rest (a) away from n start (a). For a newly added aircraft to the planning horizon, n last (a) is the gate node gate(a), n start (a) is the corresponding pushback node, L rest (a) is equal to L SP (n last (a), n start (a)) , and t start (a) equals to start_time(a). If an aircraft has already been added in the previous planning horizons, n last (a), n start (a), and L rest (a) are calculated from the result of the last horizon, and t start (a) is equal to the start time of the current planning horizon.

Decision Variables
In an MILP planning horizon, each aircraft performs K step movements, which means taxiing through K edges in the graph, G(N, E). As the model determines the taxi route and schedule of the aircraft, the information of each movement includes the target node and time of arrival at the node of each step. Clare, Richards and Sharma [21] proposed to include the above information into MILP by implementing two decision variables. Variable x represents the taxi route, which is applied as a binary variable x(a, k, n 1 , n 2 ) to MILP, in which a ∈ A, 1 ≤ k ≤ K, n 1 , n 2 ∈ N. x(a, k, n 1 , n 2 ) equals 1, which means that aircraft a taxis from start node n 1 to target node n 2 at the k th step in the plan. Another non-negative integer variable, t(a, k), determines the taxi schedule, where a ∈ A, 1 ≤ k ≤ K + 1. t(a, k) represents the time when a passes the start node, n 1 , of its k th step. Particularly, t(a, 1) is the start time of a when it is at the start location of the current planning horizon, and t(a, K + 1) is the time when a passes the target node n 2 of its final step.

Objective Function
Two parts are considered in the final objective function to define an objective for each MILP planning horizon. First, the total taxi distance of all the aircraft in the current horizon are minimized. In conventional global optimization, the total taxi distance of all the aircraft is directly set in one day or in several hours. However, in the model, the aircraft cannot arrive at their destinations in one horizon, and the regular objective function, which is given in Equation (3), cannot include the remaining taxi distances in the future; thus, the MILP model loses its direction.
This is a common problem when applying RH into MILP optimization. The design of the objective function in this study is based on the cost-to-go function computed by using the shortest-path algorithm [22]. In addition, Clare, Richards and Sharma [21] applied this approach to solve the airport ground movement optimization problem. Our model uses the real taxi distance that is already planned in the current horizon plus an estimated taxi length from the current position to the destination using the shortest-path length.
In Equation (4), obj distance is an alternative total taxi distance. The first term of Equation (4) represents the total taxi distance in the current planning horizon. The second term of Equation (4) represents the estimated taxi distance of a from the current planning horizon until its destination using the shortest-path length.
T SP (n 1 , n 2 ) =    L SP (n 1 ,n 2 ) v_max(n 1 ,n 2 ) + 120, n 2 is runway virtual node but n 1 is not L SP (n 1 ,n 2 ) v_max(n 1 ,n 2 ) , otherwise x(a, K, n 1 , n 2 ) × T SP (n 2 , destination(a)) (Minimize) obj = obj distance + w * obj time In addition to the total taxi distance, the taxi time must also be considered, thus avoiding wasting time during taxiing. In Equation (6), obj time is an alternative total taxi time using the same estimated method as for obj distance . The first term of Equation (6) represents the total taxi time length of aircraft a in the current planning horizon. The second term of Equation (6) indicates the estimated taxi time length of a from the current planning horizon until its destination. The estimated taxi time length is calculated from Equation (5). Basically, the estimated taxi time is calculated as the shortest-path length divided by the maximum speed. However, to improve the performance when a take-off aircraft does not reach the runway virtual node and ready-for-take-off point in the current planning horizon (see Section 4.5), an additional virtual long waiting time (120 s) is added to obj time . This is considered as part of the estimated taxi time length from the current planning horizon to the destination. The addition of this part can help the MILP model face congestion on the runway directly, without making meaningless detours.
The final objective function in Equation (7) is a weighted combination of obj distance and obj time , which need to be minimized. A non-negative weight parameter, w is used in obj time to control the weights of two sub-objective functions. From the equation, it can be inferred that small w values can produce better optimization of the global taxi distance, whereas large w values can reduce waiting times for the aircraft. The effect of this weight, owing to the application of RH, on the results is discussed comprehensively in Section 4.4. Herein, w = 1.5.

Constraints for Taxi Rules on the Graph Model
The taxiway and airport rules are added to the constraints. It should be noted that Equations (8)-(10), (12), and (14) refer to the model proposed in the same study [21] as the one used to determine the decision variables in Section 3.3. Equations (8) and (9) ensure that each aircraft only taxis on one existing edge at each step. Equation (10) ensures that the aircraft starts from the node where it arrived at the previous step. Equation (11) ensures that the aircraft does not use the edge of the previous step in the opposite direction. x(a, k, n 1 , n 2 ) = 1 a ∈ A, 1 ≤ k ≤ K g(n 1 , n 2 ) − x(a, k, n 1 , n 2 ) ≥ 0 n 1 , n 2 ∈ N, a ∈ A, 1 ≤ k ≤ K ∑ n 2 ∈N x(a, k − 1, n 2 , n 1 ) = ∑ n 2 ∈N x(a, k, n 1 , n 2 ) n 1 ∈ N, a ∈ A, 2 ≤ k ≤ K (10) x(a, k − 1, n 1 , n 2 ) + x(a, k, n 2 , n 1 ) ≤ 1 n 1 , n 2 ∈ N, a ∈ A, 2 ≤ k ≤ K Here, n last (a), n start (a), and L rest (a) are used to define the start position of aircraft a in each planning horizon, as discussed in Section 3.2.4. At the start time of the current planning horizon, aircraft a taxis from n last (a) to n start (a), and is located at a distance L rest (a) from n start (a). t start (a) defines the start time of aircraft a in current planning horizon. These parameters are constrained by Equations (12) and (13). An aircraft taxiing along the taxiway must follow the maximum taxi speed of the current taxiway, which is constrained by Equation (14).

Constraints to Avoid Conflict on Taxiways
To define the constraints for aircraft conflict on taxiways, some temporary variables are added to the model. The temporary variables, gap(a 1 , a 2 , k 1 , k 2 ), are used in Equation (15) and correspond to the time gap between aircraft a 1 at step k 1 and a 2 at step k 2 . In Equation (16), the binary variable, con f lict(a 1 , a 2 , k 1 , k 2 , n), equals 1 when step k 1 of aircraft a 1 and step k 2 of aircraft a 2 use the same node, n. In Equation (17), the binary variable f irst(a 1 , a 2 , k 1 , k 2 ) equals 1 when the k 1 step of aircraft a 1 is performed earlier than the k 2 step of aircraft a 2 . Figure 6 shows two types of aircraft conflict on a taxiway. The first type of conflict is a general conflict at each node. It is used as a constraint in Equation (18) as the time separation is required between each aircraft using the same node. Additionally, a second type of constraint is considered in Equations (19) and (20) to avoid head-to-head conflict between two aircraft using the same taxiway from different directions, where ' → ' represents the indicator constraints that can be solved in an MILP solver such as Gurobi. An indicator constraint z = f → a T x ≤ b states that if the binary indicator variable z is equal to f , then the linear constraint should hold. Equation (21) is used to prevent an aircraft from conflicting with the last planning horizon aircraft, where t last (n 1 ) represents the time when the aircraft last crossed node n 1 in the previous planning horizon. ∀a 1 , a 2 ∈ A, 2 ≤ k 1 , k 2 ≤ K, n ∈ N : con f lict(a 1 , a 2 , k 1 , k 2 , n) = 1 → gap abs (a 1 , a 2 , k 1 , k 2 ) ≥ t separation_node (n). (18) ∀a 1 , a 2 ∈ A, a 1 < a 2 , 1 ≤ k 1 , k 2 ≤ K, n ∈ N : f irst(a 1 , a 2 , k 1 , k 2 ) = 0 → gap abs (a 1 , a 2 , k 1 , k 2 + 1) ≥ t separation_node (n) − BIGM * (2 − x(a 1 , k 1 , n 1 , n 2 ) − x(a 2 , k 2 , n 2 , n 1 )) (19) f irst(a 1 , a 2 , k 1 , k 2 ) = 1 → gap abs (a 1 , a 2 , k 2 , k 1 + 1) ≥ t separation_node (n) − BIGM * (2 − x(a 1 , k 1 , n 1 , n 2 ) − x(a 2 , k 2 , n 2 , n 1 )) (20) ∀a ∈ A, 2 ≤ k ≤ K, n 1 ∈ N : t(a, k) ≥ ∑ n 2 ∈N t last (n 1 ) + t separation node (n 1 ) * x(a, k, n 1 , n 2 ) Aerospace 2022, 9, x FOR PEER REVIEW 12 of 23 ( , , , ) = 0 1 ∀ , ∈ , , 1 ≤ , ≤ .

Constraints to Avoid Conflict on Taxiways
To define the constraints for aircraft conflict on taxiways, some temporary variables are added to the model. The temporary variables, ( , , , ), are used in Equation (15) and correspond to the time gap between aircraft at step and at step . In Equation (16), the binary variable, ( , , , , ), equals 1 when step of aircraft and step of aircraft use the same node, . In Equation (17), the binary variable ( , , , ) equals 1 when the step of aircraft is performed earlier than the step of aircraft . Figure 6 shows two types of aircraft conflict on a taxiway. The first type of conflict is a general conflict at each node. It is used as a constraint in Equation (18) as the time separation is required between each aircraft using the same node. Additionally, a second type of constraint is considered in Equations (19) and (20) to avoid head-to-head conflict between two aircraft using the same taxiway from different directions, where '→' represents the indicator constraints that can be solved in an MILP solver such as Gurobi. An indicator constraint = → ≤ states that if the binary indicator variable is equal to , then the linear constraint should hold. Equation (21)

Constraints to Avoid Conflict in the Runway
Two types of conflict on the runway are considered in this study. First, during taxiing, aircraft can cross runway 16R-34L or runway 04-22. It is necessary to prevent these aircraft from conflicting with aircraft that are currently taking off or landing on the corresponding runway.
A temporary binary variable, cross(a 1 , a 2 , k 1 , k 2 , c), is shown in Equation (22). It is equal to 1 if aircraft a 1 crosses the runway cross point c at step k 1 when aircraft a 2 takes off or lands on the same runway at step k 2 . In Equation (23), a temporary binary variable, f irst takeo f f (a 1 , a 2 ), is used to determine the take-off order between aircraft a 1 and a 2 . It is equal to 1 when a 1 takes off before a 2 .
Equation (25) is used to prevent another type of conflict in which two take-off aircraft use the same runway. These aircraft must follow a time separation requirement set by parameter t separation_takeo f f , which is obtained from the official document "Control service processing procedure" adopted by the Tokyo International Airport. The document divides aircraft into groups A to G; different groups of aircraft have different time separation requirements when using the same runway successively, based on safety rules and on rear turbulence phenomena.

Runway Decider
Runway schedule and sequencing are determined in the MILP model based on the aforementioned constraints. However, runway management at the Tokyo International Airport also includes runway selection. According to Section 2, the airport provides two runway options for aircraft during take-off or landing. The difference in the selection can affect the taxi distance and the time to the target terminal gate. Thus, this study proposes the application of a runway decider for runway reselection. As the optimization model is developed only for airport surface operation, the landing aircraft is considered after it arrives and leaves the runway. Therefore, it is difficult to change the runway selection for the landing aircraft. For the take-off aircraft, the selection of the take-off runway is mainly based on the flight destination and its air route during actual operation. However, runway reselection is still attempted based on the runway-terminal relationship.
The runway decider is a two-stage workflow, as shown in Figure 7. The first stage is performed before the entire optimization model when the aircraft data are input. During the first stage, the best choice of runway for a take-off aircraft is selected based on its gate location. If the selection cannot meet the load capacity of the runway, the second-best choice is selected. The first stage can provide a runway selection result beforehand. During the second stage, the time when each aircraft arrives at the runway entrance is estimated, and its take-off time is calculated at both the available runways based on the queuing situation and the take-off time separation. If the other runway choice is better and its capacity is not exceeded, it is reselected in the second stage. This is a detachable module for the above optimization model, which is convenient for subsequent analysis.
choice is selected. The first stage can provide a runway selection result beforehand. During the second stage, the time when each aircraft arrives at the runway entrance is estimated, and its take-off time is calculated at both the available runways based on the queuing situation and the take-off time separation. If the other runway choice is better and its capacity is not exceeded, it is reselected in the second stage. This is a detachable module for the above optimization model, which is convenient for subsequent analysis.  Figure 8 shows the model implementation workflow, of which the central part is the model developed in Section 3. The CARATS Open Data are used as the model input. It contains the GPS ground tracking data of the aircraft in seconds, as obtained from the Tokyo International Airport. After the data processing detailed in Section 2, the following information is extracted for each flight: flight ID, start time, take-off or landing, runway ID, gate ID, and taxiway route. A one-hour range from 17:00 to 18:00 was selected on May 9, 2016, for this implementation. In this time range, the details of 34 take-off aircraft and 40 landing aircraft were recorded to be used as input data. For the MILP model, we assumed = 5 and = 50 s.

Implementation
The model was built in Python, and the MILP optimization part was built with the help of the Gurobipy Python library. Gurobi 9.11 was used as the MILP solver. The model was run on an Intel Core i5-10400 4.3 GHz PC with 16 GB RAM. An optimized result that includes the taxi route and schedule was obtained from the developed optimization model.  Figure 8 shows the model implementation workflow, of which the central part is the model developed in Section 3. The CARATS Open Data are used as the model input. It contains the GPS ground tracking data of the aircraft in seconds, as obtained from the Tokyo International Airport. After the data processing detailed in Section 2, the following information is extracted for each flight: flight ID, start time, take-off or landing, runway ID, gate ID, and taxiway route. A one-hour range from 17:00 to 18:00 was selected on 9 May 2016, for this implementation. In this time range, the details of 34 take-off aircraft and 40 landing aircraft were recorded to be used as input data. For the MILP model, we assumed K = 5 and T = 50 s.

Result Comparison
The optimized taxiway route and schedule for all the aircraft was obtained from the model results. The observed taxiway route and schedule in the history operation was obtained from the input CARATS Open Data. As the length of each taxiway is fixed, the taxi distance of each aircraft in the observed data and the optimization results can be calculated. The following figures and tables show the comparison results between our optimization results and observed data with or without runway reselection by the runway decider module.
In Figures 9 and 10, each point represents an aircraft in the one-hour time range. The x-axis represents the observed taxi distance in the historical data, while the y-axis represents the taxi distance of the optimized results. An auxiliary line, = , was used to analyze the comparison results. It was observed that very few aircraft taxi distances become longer; contrarily, most of them decrease. Table 2 shows that the rate of decrease in taxi time of the observed data without runway reselection was 9.82%. It can be observed that the space for optimization is limited, which may be because a smaller distance must be selected as far as possible in the real operation of the airport. The optimization was found to be improved after runway reselection for the take-off aircraft. The total taxi distance was found to decrease by 18.54% compared with the observed results. The model was built in Python, and the MILP optimization part was built with the help of the Gurobipy Python library. Gurobi 9.11 was used as the MILP solver. The model was run on an Intel Core i5-10400 4.3 GHz PC with 16 GB RAM. An optimized result that includes the taxi route and schedule was obtained from the developed optimization model.

Result Comparison
The optimized taxiway route and schedule for all the aircraft was obtained from the model results. The observed taxiway route and schedule in the history operation was obtained from the input CARATS Open Data. As the length of each taxiway is fixed, the taxi distance of each aircraft in the observed data and the optimization results can be calculated. The following figures and tables show the comparison results between our optimization results and observed data with or without runway reselection by the runway decider module.
In Figures 9 and 10, each point represents an aircraft in the one-hour time range. The x-axis represents the observed taxi distance in the historical data, while the y-axis represents the taxi distance of the optimized results. An auxiliary line, y = x, was used to analyze the comparison results. It was observed that very few aircraft taxi distances become longer; contrarily, most of them decrease. Table 2 shows that the rate of decrease in taxi time of the observed data without runway reselection was 9.82%. It can be observed that the space for optimization is limited, which may be because a smaller distance must be selected as far as possible in the real operation of the airport. The optimization was found to be improved after runway reselection for the take-off aircraft. The total taxi distance was found to decrease by 18.54% compared with the observed results.   The total taxi times of both the observed data and optimized results were calculated as shown in Figures 11 and 12 and Table 3, similar to the total taxi distance. However, ensuring the accuracy of the calculation of the total taxi time based on the taxi route and schedule was difficult compared with that of the total taxi distance. This is because it is difficult to simulate the change in the taxi speed of the aircraft, especially when the aircraft turns, waits, and pushes back. Only an attempt can be made to ensure that the calculation method of the taxi time is accurate and as close to the observed data as possible. The results obtained reduced the total taxi time by 22.00% without runway reselection and   The total taxi times of both the observed data and optimized results were calculated as shown in Figures 11 and 12 and Table 3, similar to the total taxi distance. However, ensuring the accuracy of the calculation of the total taxi time based on the taxi route and schedule was difficult compared with that of the total taxi distance. This is because it is difficult to simulate the change in the taxi speed of the aircraft, especially when the aircraft turns, waits, and pushes back. Only an attempt can be made to ensure that the calculation method of the taxi time is accurate and as close to the observed data as possible. The re-  The total taxi times of both the observed data and optimized results were calculated as shown in Figures 11 and 12 and Table 3, similar to the total taxi distance. However, ensuring the accuracy of the calculation of the total taxi time based on the taxi route and schedule was difficult compared with that of the total taxi distance. This is because it is difficult to simulate the change in the taxi speed of the aircraft, especially when the aircraft turns, waits, and pushes back. Only an attempt can be made to ensure that the calculation method of the taxi time is accurate and as close to the observed data as possible. The results obtained reduced the total taxi time by 22.00% without runway reselection and 29.77% with runway reselection. Particularly, the results indicated that optimization of the take-off aircraft was better than that of the landing aircraft. This is compatible with the actual operation because the airport usually assigns higher priority to landing flights. 29.77% with runway reselection. Particularly, the results indicated that optimization of the take-off aircraft was better than that of the landing aircraft. This is compatible with the actual operation because the airport usually assigns higher priority to landing flights.    29.77% with runway reselection. Particularly, the results indicated that optimization of the take-off aircraft was better than that of the landing aircraft. This is compatible with the actual operation because the airport usually assigns higher priority to landing flights.

Identification of Operation Differences
The results indicated that the developed optimization model improves the airport surface operation. The differences in the operation between the optimized and observed data were identified to obtain a better understanding of the improvement process. A taxiway usage frequency map was used to identify the operation differences. It presents the usage pattern of the optimized taxiway and the differences from the real-world scenario to clearly depict the improvement in the taxi distance and time. Figures 13 and 14 show that there are several differences between the observed and optimized operations. First, the runway cross strategy has changed. Based on the observed operation using the CARATS Open Data, landing aircraft from runway 22 to Terminal 1 or Terminal 2 primarily use point 4 to cross runway 16R during taxiing, followed by point 1; points 2 and 3 are used less frequently. However, in the optimized operation, it is observed that aircraft primarily use point 1, followed by points 2 and 3, while point 4 is used less frequently. This change may reduce the total taxi distance and waiting time.

Identification of Operation Differences
The results indicated that the developed optimization model improves the airport surface operation. The differences in the operation between the optimized and observed data were identified to obtain a better understanding of the improvement process. A taxiway usage frequency map was used to identify the operation differences. It presents the usage pattern of the optimized taxiway and the differences from the real-world scenario to clearly depict the improvement in the taxi distance and time. Figures 13 and 14 show that there are several differences between the observed and optimized operations. First, the runway cross strategy has changed. Based on the observed operation using the CARATS Open Data, landing aircraft from runway 22 to Terminal 1 or Terminal 2 primarily use point 4 to cross runway 16R during taxiing, followed by point 1; points 2 and 3 are used less frequently. However, in the optimized operation, it is observed that aircraft primarily use point 1, followed by points 2 and 3, while point 4 is used less frequently. This change may reduce the total taxi distance and waiting time.

Identification of Operation Differences
The results indicated that the developed optimization model improves the airport surface operation. The differences in the operation between the optimized and observed data were identified to obtain a better understanding of the improvement process. A taxiway usage frequency map was used to identify the operation differences. It presents the usage pattern of the optimized taxiway and the differences from the real-world scenario to clearly depict the improvement in the taxi distance and time. Figures 13 and 14 show that there are several differences between the observed and optimized operations. First, the runway cross strategy has changed. Based on the observed operation using the CARATS Open Data, landing aircraft from runway 22 to Terminal 1 or Terminal 2 primarily use point 4 to cross runway 16R during taxiing, followed by point 1; points 2 and 3 are used less frequently. However, in the optimized operation, it is observed that aircraft primarily use point 1, followed by points 2 and 3, while point 4 is used less frequently. This change may reduce the total taxi distance and waiting time.  The second difference observed was in the congestion area strategy. Based on the analysis of the observed data, it can be observed that point 6 is a congestion area. In the optimized operation, the landing aircraft uses different nodes after exiting a rapid-exit taxiway from the runway, while in the observed operation, the aircraft always uses the same node after exiting the runway. Meanwhile, lesser taxiways are used by the landing aircraft, and it is easier to provide more space for the departure aircraft.
Last, the results indicate that there is no strict rule for direction. Based on the pattern obtained earlier, Tokyo International Airport always applies direction rules on some taxiways to ensure security and efficiency. Considering point 5 as an example, the north side taxiway, K, is always used from left to right, and the south side taxiway, J, is always used from right to left. However, the approach for the optimized result is different. It always uses a closer taxiway to reach an optimized result. Therefore, it does not set taxiway direction or purpose-related rules. When there is a conflict, the optimization system uses a strategy to avoid conflict and ensure safety.

MILP Weight Parameter Analysis
The objective function discussed in Section 3.4 is a weighted objective function, in which parameter w controls the weight of obj time . For a typical global MILP optimization model, the change in w can be directly reflected in the total taxi distance and total taxi time of the final result. However, the MILP model is an optimization model with many planning horizons. The weight parameter can only affect the results of each planning horizon. Therefore, determining the actual impact of the weight parameter on the final results is crucial.
Similar to the implementation settings in Section 4.1, a one-hour time range was set, and the flight information from CARATS Open Data were used as input. Conversely, the optimization model was executed several times for different values of the weight parameter w from 0.5 to 20 in steps of 0.5. This experiment was conducted at different hours and on different dates to ensure the reliability of the results. Figure 15 shows the results of the weight parameter analysis experiments. Based on the model design, it was assumed that the taxi distance increases with the increase in w, whereas the taxi time decreases with the increase in w. The results indicate the occurrence of a sudden change when 5 ≤ w ≤ 7. Additionally, the total taxi time has a minimum value when w is approximately 7.0, after which it begins to increase. The decrease in the total taxi time was only observed for values of w between 0 and 8; therefore, w must not exceed 8. The second difference observed was in the congestion area strategy. Based on the analysis of the observed data, it can be observed that point 6 is a congestion area. In the optimized operation, the landing aircraft uses different nodes after exiting a rapid-exit taxiway from the runway, while in the observed operation, the aircraft always uses the same node after exiting the runway. Meanwhile, lesser taxiways are used by the landing aircraft, and it is easier to provide more space for the departure aircraft.
Last, the results indicate that there is no strict rule for direction. Based on the pattern obtained earlier, Tokyo International Airport always applies direction rules on some taxiways to ensure security and efficiency. Considering point 5 as an example, the north side taxiway, K, is always used from left to right, and the south side taxiway, J, is always used from right to left. However, the approach for the optimized result is different. It always uses a closer taxiway to reach an optimized result. Therefore, it does not set taxiway direction or purpose-related rules. When there is a conflict, the optimization system uses a strategy to avoid conflict and ensure safety.

MILP Weight Parameter Analysis
The objective function discussed in Section 3.4 is a weighted objective function, in which parameter controls the weight of . For a typical global MILP optimization model, the change in can be directly reflected in the total taxi distance and total taxi time of the final result. However, the MILP model is an optimization model with many planning horizons. The weight parameter can only affect the results of each planning horizon. Therefore, determining the actual impact of the weight parameter on the final results is crucial.
Similar to the implementation settings in Section 4.1, a one-hour time range was set, and the flight information from CARATS Open Data were used as input. Conversely, the optimization model was executed several times for different values of the weight parameter from 0.5 to 20 in steps of 0.5. This experiment was conducted at different hours and on different dates to ensure the reliability of the results. Figure 15 shows the results of the weight parameter analysis experiments. Based on the model design, it was assumed that the taxi distance increases with the increase in , whereas the taxi time decreases with the increase in . The results indicate the occurrence of a sudden change when 5 ≤ ≤ 7. Additionally, the total taxi time has a minimum value when is approximately 7.0, after which it begins to increase. The decrease in the total taxi time was only observed for values of between 0 and 8; therefore, must not exceed 8.

Phenomenon due to Receding Horizon
The weight analysis in Section 4.4 reveals the occurrence of a phenomenon because the approach with RH is an approximation of the global MILP optimal solution and each

Phenomenon due to Receding Horizon
The weight analysis in Section 4.4 reveals the occurrence of a phenomenon because the approach with RH is an approximation of the global MILP optimal solution and each planning horizon is only 50 s. From the results presented in Section 4.4, we find that if w is larger than 8, the total taxi time begins to increase, although it should decrease in the global optimization model. This occurs because the MILP model with RH cannot predict future congestion and only considers the immediate benefit in the current horizon. We name this phenomenon "Pretend to ignore congestion". Figure 16 explains this phenomenon. The blue aircraft is located at its start position of the current planning horizon. When K = 5, if the MILP model selects the five-step blue route, it will cause a 90 s wait in the runway entrance in the current planning horizon because of runway congestion. However, if it selects the five-step red route, although it has a longer way, the 90 s wait will be avoided, which leads to a better result in terms of the objective function. Thus, the MILP model certainly selects the second choice if the weight of obj time is large. However, the second choice does not mean that congestion will be avoided; the model just pretends to ignore congestion in this horizon.
Aerospace 2022, 9, x FOR PEER REVIEW 20 of 23 planning horizon is only 50 s. From the results presented in Section 4.4, we find that if is larger than 8, the total taxi time begins to increase, although it should decrease in the global optimization model. This occurs because the MILP model with RH cannot predict future congestion and only considers the immediate benefit in the current horizon. We name this phenomenon "Pretend to ignore congestion". Figure 16 explains this phenomenon. The blue aircraft is located at its start position of the current planning horizon. When = 5, if the MILP model selects the five-step blue route, it will cause a 90 s wait in the runway entrance in the current planning horizon because of runway congestion. However, if it selects the five-step red route, although it has a longer way, the 90 s wait will be avoided, which leads to a better result in terms of the objective function. Thus, the MILP model certainly selects the second choice if the weight of is large. However, the second choice does not mean that congestion will be avoided; the model just pretends to ignore congestion in this horizon. Two solutions can be considered to reduce the impact of this phenomenon on the results. The first method is directly related to the results presented in Section 4.4, indicating that must not exceed 8. The second method is to improve in the final objective function, which is already discussed in Section 3.4. However, the latter method only works on runway congestion, not on taxiway congestion. Because taxiway congestion does not occur at the end of the aircraft taxi process, it is impossible to add a long virtual waiting time.

Gate Re-Assignment and Analysis of Airline-Terminal Relationships
The selection of the aircraft take-off and landing runways is largely affected by the air routes. The flight route in the air space was changed by a runway decider in this study. As this study was only based on surface optimization, the gate and terminal of the flight can be re-assigned to achieve a similar optimization result without changing the selection of the runway. Table 4 shows the distribution of the airlines in each terminal before and after the gate re-assignment. Based on the operation mode of Tokyo International Airport before 2020, for domestic flights, Terminal 1 is used by JAL, Skymark Airlines, and others, while Terminal 2 is used by ANA and others. Based on the taxi distance between the different runways and terminals, the terminals can be re-assigned for unsuitable flights. Observed data Take-off  17  0  0  11  3  Landing  13  0  0  20  5 After re-assignment Take-off  4  4  13  7  3  Landing  8  13  5  7  5 The domestic flights ANA and JAL are assigned to both Terminals 1 and 2 after gate re-assignment. Table 5 shows that the optimized results with gate re-assignment are better Two solutions can be considered to reduce the impact of this phenomenon on the results. The first method is directly related to the results presented in Section 4.4, indicating that w must not exceed 8. The second method is to improve obj time in the final objective function, which is already discussed in Section 3.4. However, the latter method only works on runway congestion, not on taxiway congestion. Because taxiway congestion does not occur at the end of the aircraft taxi process, it is impossible to add a long virtual waiting time.

Gate Re-Assignment and Analysis of Airline-Terminal Relationships
The selection of the aircraft take-off and landing runways is largely affected by the air routes. The flight route in the air space was changed by a runway decider in this study. As this study was only based on surface optimization, the gate and terminal of the flight can be re-assigned to achieve a similar optimization result without changing the selection of the runway. Table 4 shows the distribution of the airlines in each terminal before and after the gate re-assignment. Based on the operation mode of Tokyo International Airport before 2020, for domestic flights, Terminal 1 is used by JAL, Skymark Airlines, and others, while Terminal 2 is used by ANA and others. Based on the taxi distance between the different runways and terminals, the terminals can be re-assigned for unsuitable flights. The domestic flights ANA and JAL are assigned to both Terminals 1 and 2 after gate re-assignment. Table 5 shows that the optimized results with gate re-assignment are better than the optimized results without gate re-assignment. In fact, gate re-assignment results in a 29.71% reduction in the total taxi distance and 37.34% reduction in the total taxi time compared with the observed data. The results of the take-off aircraft taxi time were slightly worse than those of the runway reselection optimization, but the overall results were better than those of the runway reselection optimization. In Section 4.2, Figures 11 and 12 are used to show the optimization results of each aircraft. It was found that the optimization for take-off aircraft is better than that for landing aircraft, and there is only limited space for optimization of the total taxi time of landing aircraft. Table 5 shows that gate re-assignment leads to a better optimization for landing aircraft compared with optimization with runway reselection.
Although gate re-assignment does not affect the airspace side, this analysis disregards the following facts. The gate assignment of the aircraft involves the connection of the preceding and following aircraft, but the take-off and landing aircraft in this study are independent and lose their connection data. Because flights with only a 1 h time range were considered, we do not consider the gate occupation by previous and subsequent aircraft and cannot track the same numbers of arrival and departure aircraft of the same airline in the same terminal. Further, the actual gate assignment is completed very early and cannot be based on real-time flight operation situations, whereas this analysis was based on real-time data. Therefore, the analysis presented in this section may lead to overestimation.

Conclusions
This study developed a dynamic optimization model of surface operation to apply for the Tokyo International Airport. It solves the problems of runway crossing operations, multi-runways, and multi-terminal operations. The obtained results indicate that the developed model reduces the taxi distance by 18.54% and taxi time by 29.77% compared with the observed data. This model can also identify the operation difference between optimization and reality by comparing the taxiway usage patterns. A better understanding of the improvement process is obtained from this comparison; for example, the runway cross strategy and taxiway direction rule are changed in the optimized operation. Additionally, it is observed that factors such as the MILP weights and airline-terminal relationships can affect the optimization result.
The results of this study can bring improvements to real airport surface operations. The comparison between the optimization results and the observed operations can offer suggestions to the air traffic control, such as the strategy of runway cross points selection and the taxiway direction rules. The re-selection of runway and the cooperation of terminal usage between airlines can improve the ground operation efficiency of the airport. Simultaneously, the dynamic model with a time window makes feasible the idea of applying the optimization model to practical operation.
The developed model has certain limitations. It cannot accurately simulate the speed change of the aircraft, and the comparison of the computed total taxi time with that of the observed data is not very accurate. Gate re-assignment in Section 4.6 shows a good result and presents an adequate outline of the airline-terminal relationships, but it lacks the real data generated by gate occupation in the terminal; owing to this, the results may be overestimated. In the future, the developed model can be enhanced by considering airspace side operations as the choice of runway on the ground and the landing time are closely related to the air route. Furthermore, this study focuses on optimizing the total taxi distance and time, which are direct indicators. Further benefits and indicators