Analyzing Stochastic Features in Airport Surface Traffic Flow Using Cellular Automaton: Tokyo International Airport

Global air passenger transport demand is expected to increase, and there is concern that the current airport operation will not be able to cope with aircraft overcrowding. In this study, we developed a cellular automaton (CA) simulator that can model the surface traffic of the entire Tokyo International Airport in detail, including aircraft that are arriving, taxiing, parking, and departing, using actual track data. The simulator can reproduce stop-and-go aircraft taxiing based on aircraft interactions and runway rules. It can simulate the stochastic features of the surface traffic flow. To validate the developed CA simulation, the taxiing speed distribution, local delays, and taxiing times for each route were compared with the actual track data. They were in good agreement. The effects of stochastic surface traffic features, such as arrival rate, runway occupancy time, and taxiing route, on airport operations were quantitatively analyzed. This tool could lead to a better prediction of future air traffic and improve airport operations.

International Airport in Japan, a fourth runway began oper-23 ation in 2010, increasing the annual number of departures 24 The associate editor coordinating the review of this manuscript and approving it for publication was Emre Koyuncu . and arrivals by approximately 1.4 times [2]. However, runway 25 and terminal expansions take several years and are expen- 26 sive. Thus, they can only be applied to specific airports. 27 Second, research has been conducted to improve the opera- 28 tional efficiency and throughput of airports without changing 29 current facilities by optimizing the air traffic control system. 30 There are three targets for optimization of airport opera-31 tions: arrival, surface, and departure traffic. Especially, most 32 aircraft behavior in surface traffic is entrusted to the pilot. 33 Uncertainties such as aircraft separation and taxiing speed 34 make surface traffic complex and unpredictable. In previous 35 studies, genetic algorithms [3], [4], [5] and mixed-integer 36 linear programming (MILP) [6] have been used to opti-37 mize taxiing routes. Additionally, a sequential surface traffic 38 control system that advises pilots and air traffic controllers 39 on aircraft behavior and scheduling has been proposed. 40 time, and runway occupancy time (ROT), affect the waiting 66 time of departure aircraft. However, the modeling of the 67 entire airport, i.e., the connection with arrival, parking, and 68 departure aircraft, has not been done. 69 In the past few decades, modeling of traffic flows using cel-70 lular automaton (CA) simulations has been conducted [13], 71 [14], [15], [16], [17], [18]. For air traffic flow, CA simula-72 tions are widely used to optimize the alignment of arrival 73 aircraft [19], [20]. Sekine et al. [21] conducted a step-back 74 CA simulation for arrival traffic flow at Tokyo International 75 Airport to minimize the total arrival delay and total fuel con-76 sumption. In recent years, it has been applied to surface traffic 77 flow. Mori [22], [23] developed a CA model for Tokyo Inter-78 national Airport, in which the long-range interaction among 79 aircraft is considered using a variable taxiing speed based on 80 a floor field model. The analysis was performed for specific 81 dates on airport layouts used before 2010. Yang et al. [24] Table 1 summarizes the previous studies on Tokyo Inter-97 national Airport. To the best of our knowledge, this study 98 is the first to develop a cellular automaton (CA) simula-99 tor that simulates the surface traffic flow over the entire 100 area of Tokyo International Airport, that is, arrival, taxi-101 ing, parking, and departure aircraft. It can reproduce com-102 prehensive and stochastic features of the surface traffic 103 from landing to takeoff that have not been captured ear-104 lier. It could be used to simulate future air traffic and 105 improve airport operation. This study modeled the surface 106 traffic over the entire area of Tokyo International Airport 107 and compared it with actual track data to confirm the valid-108 ity of the developed simulator. In addition, initial studies 109 were conducted to quantify the effects of arrival rate, ROT, 110 and routes on airport operations such as taxiing time and 111 delays. First, data analysis of airport surface traffic flow at Tokyo 116 International Airport was performed using Collaborative 117 Actions for Renovation of Air Traffic Systems (CARATS) 118 Open Data (COD). This is a track dataset provided by the 119 Ministry of Land, Infrastructure, and Transport (MLIT). 120 From the results, the features of traffic flow required for the 121 CA simulations were obtained and used to validate the CA 122 simulations.

123
The COD contains the timestamp, flight ID, latitude, lon-124 gitude, altitude, and aircraft type for every second. At Tokyo 125 International Airport, the runway configuration is selected 126 according to the wind direction: northerly or southerly. In this 127 study, we consider the northerly wind operation, which is the 128 primary operation at this airport. Figure 1 shows an overview 129 of the runway and terminal layouts at Tokyo International 130 Airport considered here. Tokyo International Airport has four 131 runways and three terminals. In northerly wind operations, 132 the aircraft lands on runways 34L or 34R depending on the 133 departure location. It takes off from runways 05 or 34R 134 depending on the destination. Terminals 1 and 2 are for 135 domestic flights and Terminal 3 is for international flights. 136 Each terminal has parking spots. Each parking section of the 137 terminal is called an ''Area.'' For example, terminal 1 has 138 Areas 3, 4, and 5. Area 4 has 14 parking spots. In this 139 study, we randomly selected 20 days of northerly wind 140 operation from 2016 to 2018 and averaged them for the 141 data analysis. There were aircrafts parked in areas other 142 than the eight areas mentioned earlier. However, they were 143 excluded from this analysis because they are different from 144 typical passenger operations, e.g., cargo and maintenance. 145 As a result, 87.7% of all aircraft that took off and landed 146 at Tokyo International Airport on the sampling days were 147 included in this analysis. The arrival rate, that is, the num-148 ber of arrivals per hour, was 32.78. It should be noted that 149 the maximum arrival rate was 40 when all aircrafts were 150 considered.   Fig. 3, the number of arrivals is higher from 8:00 a.m. 168 to 11:00 p.m., where domestic and Asian flights are active. 169 Runway 34L is primarily used because runway 34R is also 170 used as a departure runway. Conversely, only runway 34R 171 is used in the late-night period as it is used for European 172 and U.S. flights. For the departure runways (Fig. 4), runway 173 05 is primarily used for the departure runway. The number 174 of departures increased from 6:00 a.m. to 9:00 p.m. The area 175 usage ratios are shown in Fig. 5. The ratios of domestic areas, 176 especially in Areas 1 and 4 with many parking spots, are large 177 from 6:00 a.m. to 11:00 p.m. During this time, there are many 178 take-offs and landings, as mentioned above. Contrarily, only 179 Areas 6 and 7 for international flights were used during the 180 late night period.

182
Next, taxiing speeds were calculated from the time and loca-183 tion information of each aircraft in the COD. Their distri-184 bution is shown in Fig. 6. This is the average value during 185 the peak hours from 8:00 a.m. to 9:00 p.m. Figure 6(a) 186 shows the overall view, and (b) and (c) show the enlarged 187 views around the runways. From the distributions, the taxiing 188 speeds were 10∼20 kt on the normal taxiing routes, 20∼30 kt 189 on the high-speed sections (e.g., the straight ways from run-190 way 34R to Terminal 3), less than 5 kt around the spots, 191 and more than 120 kt on the runways. There were localized 192 low-speed regions in front of departure runways 05 and 34R. 193 This is because aircrafts that cannot enter the runway due 194 to ROT form a departure queue. It is disadvantageous from 195 the viewpoint of fuel saving as they are waiting while idling. 196 Itoh et al.     impossible to directly determine the route from landing to 212 takeoff. Two methods were used to obtain the features of 213 the aircraft routes. The first is a probabilistic method. The 214 usage ratio of each route can be calculated by multiplying the 215 arrival runway usage ratio, P arr , area usage ratio from each 216 arrival runway, P arr−area , and departure runway usage ratio 217 from each area, P area−dep , using the following equation: (1) 219 VOLUME 10, 2022 For Terminal 2, 34L→Area 1→05 was the primary route.

228
However, other routes were also used relatively frequently.

229
For Terminal 1, the ratio 34L→Area 4→05 is used much  identification. In the northerly wind operation, two arrival 243 runways, eight areas, and two departure runways were con-244 sidered. Thus, 32 standard routes were identified. Figure 2(b) 245 and (c) show a particular route, for example, 34R→Area 246 7→34R is a route along the outer edge of the airport. Parking 247 time was calculated from the difference between the spot-in 248 and spot-out times of the same aircraft. The mode frequencies 249 of parking time were calculated separately for domestic and 250 international terminals. The mode frequencies were 2700 and 251 5700 s, respectively. This result is in agreement with previous 252 studies [31].

255
In this study, CA simulation on a two-dimensional lattice 256 model proposed by Tsuzuki et al. [26] was used to reproduce 257 the airport surface traffic flow at Tokyo International Airport. 258 In this section, a basic description of the CA and the newly 259 implemented velocity-setting method and runway rules are 260 presented.  (1 cell = 4.73 m). As shown in the inset in Fig. 2( The taxiing speed in this step was v tx = 2 or 3. By repeating 313 this process, the taxiing speed becomes 2.4 cells/s, which is 314 an arbitrary particle speed. Figure 9 shows the results of the 315 verification of the aforementioned procedure on runway 34L. 316 The horizontal and vertical axes represent the setting speed 317 and actual particle speed calculated from the time required for 318 the move. The actual speed reproduced the setting speed well, 319 confirming the validity of the speed determination method. 320 In the present simulations, the taxiing speed v tx was set as 321 follows according to Fig. 6: ROT is an important factor in determining airport perfor-326 mance. It is the buffer time for safe take-offs and landings. 327 Accurate prediction and optimization of the ROT can improve 328 take-off/landing throughput [12], [28], [33], [34]. In a previ-329 ous study of Tokyo International Airport, the ROT was set at 330 95 s for consecutive takeoff and 115 s for consecutive landing 331 from the entry of the preceding aircraft onto the runway 332 [29]. In this CA simulation, it took approximately 40-50 s 333 from entry to exit on the runway. Therefore, an additional 334 buffer time of 80 s was set for the total ROT of approxi-335 mately 120-130 s. This value is reasonable compared to the 336 actual data (approximately 130 s) at runway 05 [12]. For 337 simplicity, the same ROT was used for all takeoff and landing 338 combinations. Once a landing clearance has been issued, no other aircraft 341 can take off, land, or cross on the runway. To model the run-342 way constraints, this study assumes that the landing clearance 343 should be issued before reaching 1 NM from the runway 344 arrival threshold [29]. Assuming that the flight speed in the 345 final approach is the same as v tx on the runway, the time 346 required to fly 1 NM is 1 NM/140 kt 26 s. In the present 347 CA simulation, particles enter at the end of the runway and 348 remain there for 26 s. This represents the final approach of 349 an approaching aircraft that has received landing clearance. 350 During this time, no other aircraft can use or cross the runway. 351

352
Finally, the simulation procedure is described. A time table 353 containing arrival times, routes (arrival and departure run-354 ways, and areas), and parking times for each aircraft was 355 determined based on the usage ratio in Section II-E. Here, 356 VOLUME 10, 2022   The first step was to verify the validity of the CA simulation 404 by comparing it with the COD. Figure 10 presents a snapshot 405 of the CA simulation. Note that multiple snapshots were 406 combined and illustrated for the reader's clear understanding. 407 Each particle was colored according to its status. Queues of 408 several aircrafts were generated in front of departure runways 409 05 and 34R caused by the ROT. The departing aircraft could 410 not cross runway 34L and was stopped because the approach-411 ing aircraft (APP) was standing on runway 34L.
412 Figure 11 shows the taxing speed distribution obtained 413 from the CA simulation. As described in Section III-B, four 414 taxiing speeds were set in the CA simulation. However, there 415 were several areas where interactions with other aircrafts 416 reduced the speed, causing delays. The first is the blue region 417 in front of the departure runway, which corresponds to the 418 queue before takeoff. The second is the blue region (shown 419  hour and is one of the indicators of airport performance. 454 The red and blue bars represent the arrival throughput using 455 runways 34L and 34R, respectively. The green and orange 456 bars represent the departure throughput using runways 05 and 457 34R, respectively. As λ arr increases, the throughput increases 458 by the same amount, indicating that take-offs and landing at 459 the airport occur without problems as scheduled. As men-460 tioned earlier, the actual max arrival rate was 40. Thus, this 461 result is reasonable because this analysis is within the actual 462 runway allowance. Figure 14 shows the delays in arrival and 463 departure taxiing as functions of λ arr . Here, delay is defined 464 as the increment from the taxiing time, T tx,min , when the 465 aircraft travels without interactions with other aircraft, that 466 is, without delay. The delays on arrival and departure taxiing 467 are denoted as T d,arr and T d,dep , respectively. These delays 468 include stops to avoid collisions with other aircraft, stops 469 before runways owing to ROT and queues, and stops prior to 470 runway crossings. In arrival taxiing, there is approximately 471 zero delay, regardless of the arrival rate. However, in the case 472 of departure taxiing, the delay increases as the scheduled 473 arrival rate increases. Figure 15 shows the taxiing speed dis-474 tribution around runway 05 when the arrival rates are 30 and 475 35, respectively. The queue of departures is longer when the 476 arrival rate is 35. This queue causes a departure waiting time, 477 resulting in increased delays in departure taxiing. For optimal 478 airport operation, it is necessary to achieve a high throughput 479 and low delay. This result suggests that an increase in the 480 arrival rate increases the throughput and local delay. Thus, 481 the departure schedule such as off-block time needs to be 482 optimized to remove the queues. The ROT limits the entry of a departure aircraft onto the 485 runway, creating queues for departure aircrafts. Optimizing 486 the ROT within the range of safe take-offs and landings is one 487 way to improve the airport performance. Here, we evaluated 488 VOLUME 10, 2022

499
Itoh et al. [27] analyzed the traffic flow of aircraft arriving 500 at Tokyo International Airport using a queue-based model 501 to identify the optimal arrival strategy based on the distance 502 from the arrival airport. They showed that fluctuations in 503 flight time and arrival rate at each flight area are important 504 parameters in determining the delay. Therefore, in this study, 505 we evaluate the effect of fluctuations in the arrival interval 506  The terminals on the six long-distance routes were changed as 538 shown in Table 5. For example, in route 34L→Area 1→05,   are given as follows: Figure 18 shows the total taxiing times and delay times 549 as functions of terminal exchange ratio α. Here, T tx,total is 550 the average of the sum of arrival/departure taxiing times 551 among all sampling aircrafts. T * tx,total is the value from which 552 delays are excluded, T * tx,total = T tx,total − T d,arr − T d,dep . As 553 α increases, T * tx,total decreases. It indicates that the taxiing 554 time excluding delays decreases owing to route reassign-555 ment. However, as T * tx,total decreases, the departure interval 556 becomes smaller. Thus, a queue before the departure runway 557 is more likely to form. Consequently, the departure taxiing 558 delay increases with α. T tx,total changes negligibly regardless 559 of α. Similar to the results for the arrival rate mentioned 560 earlier, this problem could be solved by optimizing the depar-561 ture schedule. Spot utilization beyond the terminal-airline 562 restriction could become an important factor in optimizing 563 airport operations.

565
In this study, CARATS Open Data (COD), which is actual 566 track data, was analyzed to reveal the surface traffic features 567 at the Tokyo International Airport. Using the features-568 arrival rates, routes, taxiing speeds, and usage ratios-569 obtained from the COD, a cellular automaton (CA) simula-570 tion was developed to reproduce the complex surface traffic 571 flow over the entire area of Tokyo International Airport. 572 To validate the CA simulation, the taxiing speed distribution, 573 local delays before runways, and taxiing times for each route 574 were compared with the COD. They were in good agreement. 575 In reality, it is difficult to evaluate the effects of vary-576 ing various operational parameters at an operational airport.

577
As an initial study, the effects of surface traffic features-