Throughput reduction on the air-ground transport system by the simultaneous effect of multiple traveling routes equipped with parking sites

This paper proposes a stochastic lattice model of the air-ground transport system that comprises a junction of two traveling routes: the domestic route and the international route, each of which has parking sites. The system distributes the arrived aircrafts to either of the two routes and selects one of the parking sites in the route for each aircraft, which stops at the parking site once during its travel. The polyhedral geometric algorithm is introduced to represent inter-vehicular distance and detect other approaching aircrafts. The system displays interesting behavior; remarkably, the dependence of the throughput on the distribution ratio of aircrafts to the domestic route reduces after reaching the maximum parking capacity of the domestic route. Our simulations and analysis with the queueing model describe this phenomenon and suggest the following fact: As the distribution ratio of aircrafts to the international route decreases, the throughput of the international route reduces, and simultaneously, that of the domestic route saturates. This simultaneous effect of the decrease and saturation causes a reduction in the throughput of the entire system.


INTRODUCTION
Effective control of an air-ground transport system plays an essential role in airport terminal management. Notably, the reduction of throughput becomes a crucial factor in the deterioration of efficiency in the entire airport system. However, in many cases, it is highly challenging to specify the cause of throughput reduction because of the mutually-correlated multiple factors (e.g., the delay of spotting, the congestion on aircraft taxiing, and the paths of complex geometries); numerical simulations are necessary to investigate such complex systems.
Several numerical simulations of the air-ground transport systems have been reported. Nevertheless, little studies have reported on the simulations using cellular automata [1,2], which are expected to become promising approaches since the cellular automata have solved a variety of similar systems in pedestrian dynamics [3,4,5,6,7,8,9] and automobile transport systems [10,4,11,12] for years. In comparisons with a lot of the previous analytical studies using the optimization and scheduling algorithms [13,14,15,16,17], the cellular automata models have the advantage of the easiness to examine the effect of the geometry of the targeting airports.
In this paper, we propose a stochastic lattice model of the airground traffic system that comprises a junction of two traveling routes, each of which has adjacent parking sites. The whole system distributes arrived aircrafts to either of the traveling routes and selects one of the parking sites of the route for each aircraft. The aircraft stops at the designated site of the route once during Email addresses: tsuzuki.satori@mail.u-tokyo.ac.jp (Satori Tsuzuki), tDaichi@mail.ecc.u-tokyo.ac.jp (Daichi Yanagisawa), tknishi@mail.ecc.u-tokyo.ac.jp (Katsuhiro Nishinari) its travel. When a aircraft finds that the targeted site is in use, the aircraft changes its destination to one of the other parking sites that is vacant.
Our model has the combined characteristics of both conventional cellular automata and the polyhedral methods of computational geometries [18,19]. Each aircraft hops on a checkerboard in the arbitrary one-way direction in 2D space. An advantage of using our model is that it is easy to reproduce the smooth motion of vehicles in a real-world system, as is the case with using a continuous model such as the optimal velocity model [20,21,22,23] while keeping the simplicity of the lattice-based model using fine-grained cells. Additionally, in our model, the polyhedral geometric algorithm is introduced to represent inter-vehicular distance and detect other approaching aircrafts. These characteristics suggest that our model has high flexibility in both translation and rotational motions compared to the ordinary stochastic lattice-based model.
As preliminary works, we examined the characteristics of one-dimensional cellular automata on a single traveling lane equipped with the functions of site assignments to the adjacent parking sites in [24,25]. Our current model can be said to be a combined model of our two previous models in a different scale in a broader sense; however, the existence of a single junction and the introduction of inter-vehicular distance are the distinguishing factors compared to the previous models.
The remainder of this paper is structured as follows. Section 2 describes our model with the priority rules at each intersection and an effective method to detect the approaching aircrafts. In Section 3, we investigate the characteristics of the system properties through simulations and compare them with the theoretical queueing models [26]. Section 4 summarizes our results  Figure 1 depicts the basic concept of the proposed model. A checkerboard spread over the whole system, and the selected cells from the checkerboard, called checkpoints, construct a route. In the example of the left part in Fig.1, cells A and B represent the start and end cells, respectively. A yellow-colored particle represents the center of a aircraft. The particle hops from cell A to cell B along the relative vector obtained by connecting the centers of both cells. At this time, the particle moves along the vector by a distance equal to the length of a single cell in every time step and stores itself on the closest cell from the relative vector.

Basic concept of our model
In our model, each particle has a pair of antennas in the back and front directions to represent inter-vehicular distance among aircrafts as shown in the right part of Fig.1. We detect the intersections and the contacts of two different antenna particles by using the polyhedral geometric algorithm as mentioned later in Section 2.3. When the system detects the intersection, the system instructs both or either of the aircrafts to stop at the current location, according to the priority rules predefined at the junction. Figure 2 depicts the schematic of the target system. We abstract the target system from a real-world airport, Fukuoka airport in Japan, which has one of each domestic and international terminals; the roads on the background show the coarse-grained geometries of the airport for reference. We denote the number of cells in the x-direction by L x and denote those in the y-direction by L y . The three symbols (circle, square, and triangle) indicate the checkpoints on the checkerbord. The blue line with the circular symbol shows a domestic route R D comprising the parking lane G D that has 23 parking sites, and the red line with the triangle symbol shows a international route R I comprising the parking lane G I that has 12 parking sites. The upper and lower green lines with square symbols respectively show the lanes utilized for the arrival and departure of aircrafts.

Target system
Before the simulation, we establish a timetable that has five instructions for each aircraft: (a) the arrival time t arr determined by the fixed interval of arrival time T in , (b) the type of routes (the domestic route R D or the international route R I ), (c) the index of the parking site of the route, (d) the scale and shape parameters that determine the staying time T sp as described later, and (e) the interval velocities among the checkpoints. We set the arrival time t arr of each aircraft such that all intervals of the arrival times become the same constant value T in . In additions, we assign the domestic route R D to aircrafts with probability α, and we assign the international route R I to them with the probability 1 − α. In the end, we select one of the parking sites on the route for each aircraft at uniform random.
During the simulation, each aircraft enters the system from inlet A with a fixed interval T in and moves on the path ABC towards the checkpoint C. The system distributes the aircraft to either of the routes according to the timetable at checkpoint C. At checkpoints E or F, the aircraft checks the state of the parking site designated in the timetable and changes it when the parking site is still in use. In such a case, the aircraft selects one of the parking sites on the route from the rest of the parking sites whose state is vacant. Then, the aircraft moves towards the parking site.
After reaching the parking site, the aircraft stays at the parking site during time T sp . In this paper, we investigate the target system in both cases with stochastic and deterministic parameters to determine the effect of the delay of staying time T sp . In the former case, we consider the delay of the staying time T sp by using a half-normal distribution. The half-normal distribution is selected because the delay only occurs for the positive direction in the target system. The scale and shape parameters (τ s , σ s ) of the half-normal distribution are obtained from the mean and deviation of (τ s ,σ s ) of a normal distribution, as follows: In the latter case, we decide the T sp by settingσ s to zero. The aircrafts move toward checkpoint G after vacating the parking sites. The flows of the aircrafts from both, the domestic route R D and international route R I , merges at checkpoint G. When some aircrafts go through path HG on the domestic route R D , other aircrafts move on path IJG, and the system pauses the aircrafts on path HG. After passing checkpoint G, the aircraft moves on path GKL and exits from the outlet L. Note that we describe the setting of interval velocities among the checkpoints in Section.3.2 .

Detection of approaching antenna aircrafts
In this section, we apply a polyhedral geometrics algorithm to the proposed model to detect approaching aircrafts. Each aircraft has a pair of antennas in the back and front directions as explained in Section 2.1. Let the number of aircrafts that exist on the system be N p . We set a sequential number from 1 to N p to the aircrafts. We detect the collision of antennas among the i-th aircraft and the other j-th aircraft, as follows.
First, we judge whether the i-th and j-th aircrafts are on an identical line by using the relative vector from the i-th aircraft to the j-th aircraft x ij , the normal vector of i-th aircraft n i , the angle φ = arccos (n i ·x ij ) |n i ||x ij | (the angle between the normal vector n i and the relative vector x ij ), and the small angle . If the condition: is satisfied, we judge that the i-th and j-th aircrafts are on an identical line. becomes zero in an ideal case; however, we must set it to a larger value than the scale of the "machine epsilon." In case φ meets the relationship in Eq. (3), we decide that the j-th aircraft is on an identical line as the i-th aircraft, which indicates that these two aircrafts are on the same path. In this case, we can judge the collision of the antennas by simply measuring the distance between both aircrafts. At this time, we have to identify the back and front positions of the two aircrafts and their orientations because the system needs to issue a pause instruction to the back aircraft. By judging from the normal vector of the i-th aircraft n i , of the j-th aircraft n j , and the relative vector x ij , four interaction types exist, as listed below.
(A) If (n i · n j ) > 0 and (n i · x ij ) > 0, then: The i-th aircraft locates the back of the j-th aircraft. The system issues a pause instruction to the i-th aircraft.
(B) If (n i · n j ) > 0 and (n i · x ij ) < 0, then: The j-th aircraft locates the back of the i-th aircraft. The system issues a pause instruction to the j-th aircraft.
(C) If (n i · n j ) < 0 and (n i · x ij ) > 0, then: Both aircrafts face each other. This case only occurs at junction G in the target system. The system issues a pause instruction to the either of the aircrafts according to the priority rule predefined at junction G.
(D) If (n i · n j ) < 0 and (n i · x ij ) < 0, then: Both aircrafts move in the opposite direction. Since this situation does not occur in the target system, the system ignores this type of interaction.
In case φ does not meet the relationship shown in Eq.(3), we decide that the i-th and j-th aircrafts are not on an identical line, which suggests that these two aircrafts are on different crossing paths. In such a case, we judge the collisions among these antennas by using the line-line collision detection algorithm used in polyhedral geometrics. Let the edge positions of the antennas of the i-th aircraft be a(a x , a y ) and b(b x , b y ), and let those of the j-th aircraft be c(c x , c y ) and d(d x , d y ). At this time, the arbitrary point of p ab on the line ab and that of p cd on the line cd are expressed as Here, u i and u j indicate the scalar value uniquely determined by each point. In case line ab and line cd cross, the left-hand side of Eq.(4) and that of Eq.(5) correspond to each other. At this time, we obtain the unique values of u i and u j by solving the simultaneous linear system of Eq.(4) and Eq.(5). The cross of two lines only occurs when 0 < u i < 1 and 0 < u j < 1: (E) If 0 < u i < 1 and 0 < u j < 1, then: The i-th aircraft and j-th aircraft collide with each other.
The system issues pause instructions to the either of the two aircrafts according to the priority rule defined at the joint.
Unlike the case that both aircrafts exist on the same lane, we do not divide the interaction (E) into four subtypes of interactions by their normal vectors because every corner is separated into two paths by the single checkpoint. Namely, once the interaction (E) is detected, the system identifies the exact corner by Interaction types between i-th antennas particle and j-th antennas particle (A) (B) (C) j-th particle the location of the aircrafts and pauses the aircraft that belongs to the rear path of the corner. For example, the system pauses the aircraft existing on path CD in case it detects an interaction (E) at checkpoint D in Fig.2. Figure 3 summarizes equations and the schematics of all interaction types from (A) to (E).
In our model, the system issues multiple instructions to the aircrafts when they go on to the forks. (e.g., the forks around checkpoint C in Fig.2). In this case, the aircrafts pause when at least one instruction given to them suggests a pause. For better understanding, we describe the final decision-making of a aircraft using the following formula. Let the total number of aircrafts approaching the i-th aircraft be N i ; set a serial number to the approaching aircrafts between 1 and N i . Donate a binary function of the m-th aircraft of the N i aircrafts as I m i , which returns one when the system issues pause to the i-th aircraft, and it returns zero for other cases. The final decision-making D i of the i-th aircraft is obtained as In case that D i is true, the i-th aircraft pauses.

Definitions of the physical properties
In this paper, we examine three physical properties: the throughput Q of the whole system, changing rate of schedule S , and averaged occupation of the parking sites U T of the target system. First, we define throughput Q as Here, N exit indicates the number of aircrafts that exit from outlet L of Fig.2 during the simulation. N plan shows the number of aircrafts that plan to enter the system from inlet A each day. Second, we define N chg as the number of aircrafts that changes the destination of the parking site when going through either checkpoint E or F (E for the domestic route and F for the international route). We define the changing rate of schedule S as In the end, we describe U D as the number of occupied parking sites in parking lane G D and U I as that of the occupied parking sites in parking lane G I . U D and U I are averaged by the total time step N t . We define U T as U T := U D + U I (10)

Setting of input parameters
We set the length of a square cell on the checkerboard to 2.2 m and the number of cells of (L x , L y ) to (382, 1, 483) so that our system fits the real geometry of the targetting airport. We adopt a parallel update method; the system updates all cells on the checkerboard simultaneously. We set the physical time of a time step to be 1 s and set the number of total time steps N t to 86, 400 for each case, which corresponds to 24 h of a day. We perform simulations for different cases of the distribution parameter α in between 0 and 1.0 at intervals of 0.05. We carry out 15 cases in each simulation. As mentioned in Section 2.2, we investigate the system by setting the staying time T sp in both cases, with stochastic and with deterministic parameters. In the former case, we set the pair of stochastic parameters (τ s ,σ s ) of T sp to (60 min, 10 min), aiming to reproduce the typical staying time at real airport terminals. In the latter case, we set (τ s ,σ s ) to (60 min, 0 min). We set the interval velocities by calculating the substep at each time step; 30 substeps for path AB and path  Figure 4 shows the dependence of throughput Q of the whole system on probability α, which distributes the aircrafts to the domestic route R D , for different cases of the arrival time interval T in (from 2 min to 4 min). The blue-colored circle symbol shows the results when setting T sp with deterministic parameters. The orange-colored triangle indicates results in the case of setting T sp with stochastic parameters using the half-normal distribution. The dashed line shows the ratio of the number of parking sites of the domestic route G D to that of the total system (G D + G I ). This ratio corresponds to the maximum capacity of the parking lane G D . In case of setting T in to 4 min, the throughput Q increases as probability α increases, and it reaches a plateau after α becomes more larger than a specific value. We observe similar phenomena in the case of 3 min.

Simulation results
In the case of 2 min, an important observation was made; the throughput Q reduces after α reaches the maximum capacity. The overall reduction of throughput Q in between 0.6 and 1.0 seems to be unique because the throughput plateaus after saturation, in similar systems [27,28,29]. Besides, we confirmed from Fig.4 that the existence of the delay of staying time is not the critical factor for the overall reduction because we observe it in both cases: with stochastic and with deterministic parameters. Indeed, the delay of T sp causes a substantial deterioration of throughput Q in exceptional cases; however, it has little effect on the mainstream behavior of the whole system. Figure 5 shows the dependence of the changing rate S of the whole system on probability α that distributes the aircrafts to the domestic route R D , for different cases of arrival time interval T in from 2 min to 4 min. Each symbol is the same in Fig.4. The results show further unexpected behaviors. It is not easy to find the apparent dependence on parameter α in every case. Furthermore, the curves are entirely different among the three cases of setting T in to 2 min, 3 min, and 4 min regardless of whether it is with stochastic or with deterministic parameters. Hereafter, we focus on the deterministic cases to clarify the reasons for the phenomena observed in throughput Q and changing rate S .
An obvious reason was found attributable to these unique behaviors. Figure 6 shows the sample mean of the throughput displayed in Fig.4 (the solid line with the green-square symbol) and its breakdown: the throughput of the domestic route R D (the solid line with the blue-colored circle symbol) and that of the international route R I (the solid line with the red-colored triangle symbol). As shown in Fig.6, the sum of the throughputs of the domestic route R D and the international route R I becomes the throughput of the whole system. In the case of 4 min, the throughput of the international route R I shows an almost plateau when the probability α < 0. these cases. In the case of 2 min, the throughput of the international route R I decreases similarly as that in the cases of 3 min and 4 min. On the other hand, the throughput of the domestic route R D increases until it reaches the maximum capacity of the domestic lane G D , and after that, it saturates unlike the cases of 3 min and 4 min. Simply put, the reduction of the throughput of the whole system shown in the case of 2 min is described as follows. As the distribution ratio of the aircrafts to the international route R I decreases, the throughput of the international route R I reduces while at the same time, that of the domestic route R D saturates, the simultaneous effect of which causes the reduction in the throughput of the whole system. A similar explanation can be provided for every case of the changing rate S , as shown in Fig.7. Although each route simply shows either monotonic decrease, monotonic increase, or saturation, their sum displays very complex behavior.

Analysis with the classical M/M/c queue
The simulation results in Section 3.3 indicate that the function of the queueing system of each route independently contributes to the whole system. In this section, we compare the simulation results of the average occupied time with the classical queueing model by Erlang [26] (For the details of the classical theoretical queueing model, refer to [30,31]). Since our system has a finite number of service windows (the parking sites) in each route, we apply the M/M/c queueing model [30] to each route of the target system as the first-order approximation level.
In the classical M/M/c queue, the average number of occupied sites among c sites corresponds to the ratio of the arrival rate λ to the service rate µ. Since the target system directs Here, N D and N I indicate the number of parking sites in the domestic lane G D and that in the international lane G I , respectively. From Eq.(10), Eq.(11), and Eq. (12), we obtain the average number of occupied sites U T as  It was found that the M/M/c queueing models describe the simulation results appropriately, and they strongly support the fact that the simultaneous effect of the decrease in the throughput of the international route R I and the saturation of that of the domestic route R D causes the reduction in the throughput of the whole system in the case of 2 min. For further investigations, we describe the following observations. Whereas the simulation results of route R I show good agreement with approximations in all cases from 2 to 4 min, the simulations of the domestic route R D was observed to deviate from the approximations as T in decreases because of the longer distance of route R D compared to route R I , and the lack of the congestion effect in the classical M/M/c queueing model. Figure 9 shows the three-dimensional contour plot of the dependencies of the throughput Q on the distribution probability α and the interval of arrival time T in . It was observed that throughput Q gets sluggish and reaches the plateau as T in increases, thus making a type of sigmoid curve, the tendency of which is similar to the experimental data of the same type of systems in previously reported works [27,28,29]. On the other hand, the dependence of throughput Q on arrival time interval T in does not reach the plateau in between 4 and 2 min at around where the parameter α gets close to the ratio of the number of parking sites in the domestic route to that in the whole system. It should be emphasized that, the α-dependence of the physical properties of this kind of system was first discussed in this paper. Particularly, we found that the simultaneous effect of the decrease, increase, and saturation of the physical properties in the two different routes causes unusual behaviors of the whole  system. Since the targeting system fixes the number of parking sites in this study, there might exist a possibility to find another interesting fact caused by the aforementioned simultaneous effects in different situations; these other facts should be studied in future works.

CONCLUSIONS
We introduced a two-dimensional stochastic lattice model that comprises a junction of two traveling lanes, each of which has a parking lane. In our model, each aircraft has a pair of antennas in the moving direction to detect other approaching aircrafts. As a case study, we abstracted a target system from a real-world airport that has one of each domestic and international terminals. We applied our model to the target system and investigated the physical characteristics of the proposed model thoroughly. The contribution of this paper is as follows: The proposed system was observed to display interesting behavior. The throughput shows a sudden decrease after reaching the maximum capacity of the parking lane. On the other hand, the dependence of the changing rate of schedule on the system parameter disorderly fluctuates, and it hides its dependence on the parameter. Our simulations and approximations by the M/M/c queueing model clearly explain these phenomena and support the following fact. As the distribution ratio of aircrafts to the international route decreases, the throughput of the international route reduces; simultaneously, the throughput of the domestic route R D saturates. This simultaneous effect of the decrease and saturation causes reduction in the throughput of the whole system.