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

Figure 1 depicts the basic concept of our 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 figure 1, cells A and B represent the start and end cells, respectively. 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.

Zoom In Zoom Out Reset image size

In our model, each particle has a pair of antennas in the back and front directions as shown in the right part of figure 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 particles to stop at the current location, according to the priority rules predefined at the junction. Figure 2 shows a flowchart of our method. In the following sections, (a) and (b) are described in section 2.2, and (c), (d), and (e) are explained in section 2.3.

Zoom In Zoom Out Reset image size

Figure 3 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 particles.

Zoom In Zoom Out Reset image size

Before the simulation, we establish a timetable that has five instructions for each particle: (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 particle 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 particles 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 particle at uniform random.

During the simulation, each particle 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 particle to either of the routes according to the timetable at checkpoint C. At checkpoints E or F, the particle 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 particle selects one of the parking sites on the route from the rest of the parking sites whose state is vacant. Then, the particle moves towards the parking site.

After reaching the parking site, the particle 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 (${\tau }_{{\rm{s}}},{\sigma }_{{\rm{s}}}$) of the half-normal distribution are obtained from the mean and deviation of (${\bar{\tau }}_{{\rm{s}}},{\bar{\sigma }}_{{\rm{s}}}$) of a normal distribution, as follows:

$$\begin{eqnarray}&&{\tau }_{{\rm{s}}}={\bar{\tau }}_{{\rm{s}}}+{\bar{\sigma }}_{{\rm{s}}}\sqrt{\displaystyle \frac{2}{\pi }}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\sigma }_{{\rm{s}}}={\bar{\sigma }}_{{\rm{s}}}\sqrt{\left(1-\displaystyle \frac{2}{\pi }\right)}\end{eqnarray} \tag{ 2 }$$

In the latter case, we decide the T_sp by setting ${\bar{\sigma }}_{{\rm{s}}}$ to zero.

The particles move toward checkpoint G after vacating the parking sites. The flows of the particles from both, the domestic route R_D and international route R_I, merges at checkpoint G. When some particles go through path HG on the domestic route R_D, other particles move on path IJG, and the system pauses the particles on path HG. After passing checkpoint G, the particle 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.

In this section, we describe a polyhedral geometrics algorithm for the detection of approaching antenna particles. Let the number of particles that exist on the system be N_p. We set a sequential number from 1 to N_p to the particles. We detect the collision of antennas among the i-th particle and the other j-th particle, as follows.

First, we judge whether the i-th and j-th particles are on an identical line by using the relative vector from the i-th particle to the j-th particle ${{\boldsymbol{x}}}_{\mathrm{ij}}$, the normal vector of i-th particle ${{\boldsymbol{n}}}_{{\rm{i}}}$, the angle $\phi =\left|\arccos \left\{\tfrac{({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{x}}}_{\mathrm{ij}})}{| {{\boldsymbol{n}}}_{{\rm{i}}}| | {{\boldsymbol{x}}}_{\mathrm{ij}}| }\right\}\right|$ (the angle between the normal vector ${{\boldsymbol{n}}}_{{\rm{i}}}$ and the relative vector ${{\boldsymbol{x}}}_{\mathrm{ij}}$), and the small angle . If the condition:

$$\begin{eqnarray}&&\phi \,\,\lt \,\,\epsilon \,\,\mathrm{or}\,\,\phi \,\,\gt \,\,\pi -\epsilon \end{eqnarray} \tag{ 3 }$$

is satisfied, we judge that the i-th and j-th particles 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 equation (3), we decide that the j-th particle is on an identical line as the i-th particle, which indicates that these two particles are on the same path. In this case, we can judge the collision of the antennas by simply measuring the distance between both particles. At this time, we have to identify the back and front positions of the two particles and their orientations because the system needs to issue a pause instruction to the back particle. By judging from the normal vector of the i-th particle ${{\boldsymbol{n}}}_{{\rm{i}}}$, of the j-th particle ${{\boldsymbol{n}}}_{{\rm{j}}}$, and the relative vector ${{\boldsymbol{x}}}_{\mathrm{ij}}$, four interaction types exist, as listed below.

• (A)  
   If $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{n}}}_{{\rm{j}}})\gt 0$ and $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{x}}}_{\mathrm{ij}})\gt 0$, then:The i-th particle locates the back of the j-th particle. The system issues a pause instruction to the i-th particle.
• (B)  
   If $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{n}}}_{{\rm{j}}})\gt 0$ and $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{x}}}_{\mathrm{ij}})\lt 0$, then:The j-th particle locates the back of the i-th particle. The system issues a pause instruction to the j-th particle.
• (C)  
   If $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{n}}}_{{\rm{j}}})\lt 0$ and $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{x}}}_{\mathrm{ij}})\gt 0$, then:Both particles 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 particles according to the priority rule predefined at junction G.
• (D)  
   If $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{n}}}_{{\rm{j}}})\lt 0$ and $({{\boldsymbol{n}}}_{{\rm{i}}}\cdot {{\boldsymbol{x}}}_{\mathrm{ij}})\lt 0$, then:Both particles 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 equation (3), we decide that the i-th and j-th particles are not on an identical line, which suggests that these two particles 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 particle be a(a_x, a_y) and b(b_x, b_y), and let those of the j-th particle 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

$$\begin{eqnarray}&&{{\bf{p}}}_{\mathrm{ab}}={\bf{a}}+{u}_{{\rm{i}}}({\bf{b}}-{\bf{a}})\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\bf{p}}}_{\mathrm{cd}}={\bf{c}}+{u}_{{\rm{j}}}({\bf{d}}-{\bf{c}})\end{eqnarray} \tag{ 5 }$$

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 equation (4) and that of equation (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 equation (4) and equation (5). The cross of two lines only occurs when 0 < u_i < 1 and 0 < u_j < 1:

•  If 0 < u_i < 1 and 0 < u_j < 1, then:The i-th particle and j-th particle collide with each other. The system issues pause instructions to the either of the two particles according to the priority rule defined at the joint.

Unlike the case that both particles 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 the location of the particles and pauses the particle that belongs to the rear path of the corner. Figure 4 summarizes equations and the schematics of all interaction types from (A) to (E).

Zoom In Zoom Out Reset image size

In our model, the system issues multiple instructions to the particles when they go on to the forks. (e.g., the forks around checkpoint C in figure 3). In this case, the particles pause when at least one instruction given to them suggests a pause. For better understanding, we describe the final decision-making of a particle using the following formula. Let the total number of particles approaching the i-th particle be N_i; set a serial number to the approaching particles between 1 and N_i. Donate a binary function of the m-th particle of the N_i particles as I_i^m, which returns one when the system issues pause to the i-th particle, and it returns zero for other cases. The final decision-making D_i of the i-th particle is obtained as

$$\begin{eqnarray}&&{I}_{{\rm{i}}}^{{\rm{s}}}=\sum _{m\in {N}_{{\rm{i}}}}{I}_{{\rm{i}}}^{{\rm{m}}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{D}_{{\rm{i}}}=\left\{\begin{array}{l}\mathrm{true}\,\,\,({I}_{{\rm{i}}}^{s}\ne 0)\\ \mathrm{false}\,\,(\mathrm{otherwise})\end{array}\right.\end{eqnarray} \tag{ 7 }$$

In case that D_i is $\mathrm{true}$, the i-th particle pauses. In the end, the position of the i-th particle is updated according to the velocities and routes recorded on the timetable, the priority rules, and the instruction issued by equation (7).

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

$$\begin{eqnarray}&&Q:= \displaystyle \frac{{N}_{\mathrm{exit}}}{{N}_{\mathrm{plan}}}\end{eqnarray} \tag{ 8 }$$

Here, ${N}_{\mathrm{exit}}$ indicates the number of particles that exit from outlet ${\rm{L}}$ of figure 3 during the simulation. N_plan shows the number of particles that plan to enter the system from inlet A each day.

Second, we define ${N}_{\mathrm{chg}}$ as the number of particles 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

$$\begin{eqnarray}&&S:= \displaystyle \frac{{N}_{\mathrm{chg}}}{{N}_{\mathrm{exit}}}\end{eqnarray} \tag{ 9 }$$

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

$$\begin{eqnarray}&&{U}_{{\rm{T}}}:= {U}_{{\rm{D}}}+{U}_{{\rm{I}}}\end{eqnarray} \tag{ 10 }$$

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, 1483) 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, in which arrivals and departures are repeated to ensure sufficient statistics. 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 $({\bar{\tau }}_{{\rm{s}}},{\bar{\sigma }}_{{\rm{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 $({\bar{\tau }}_{{\rm{s}}},{\bar{\sigma }}_{{\rm{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 KL, 15 substeps for path BC and path GK, and 3 substeps for other paths. Each of these values corresponds to 237.6 km h⁻¹, 118.8 km h⁻¹, and 23.7 km h⁻¹ in the physical scale, respectively. We carry out each simulation for different cases of intervals of arrival time T_in in 2 min, 3 min, and 4 min.

Figure 5 shows the dependence of throughput Q of the whole system on probability α, which distributes the particles 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}_{{\rm{D}}}+{G}_{{\rm{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.

Zoom In Zoom Out Reset image size

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 [25–27]. Besides, we confirmed from figure 5 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 6 shows the dependence of the changing rate S of the whole system on probability α that distributes the particles 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 figure 5. 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.

Zoom In Zoom Out Reset image size

An obvious reason was found attributable to these unique behaviors. Figure 7 shows the sample mean of the throughput displayed in figure 5 (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 figure 7, 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.2 because of the saturation of the parking lane, and it decreases when α > 0.2. Meanwhile, the throughput of the domestic route R_D increases because the amount of particles increases proportional to probability α. The same is almost true in the case of 3 min. These linear increases can be attributed to the fact that the traffic of particles does not surpass the capacity of the domestic lane G_D in 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.

Zoom In Zoom Out Reset image size

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 particles 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 figure 8. Although each route simply shows either monotonic decrease, monotonic increase, or saturation, their sum displays very complex behavior.

Zoom In Zoom Out Reset image size

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 [24] (For the details of the classical theoretical queueing model, refer to [28, 29]). Since our system has a finite number of service windows (the parking sites) in each route, we apply the $M/M/c$ queueing model [28] 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 each particle to the domestic route R_D with probability α, the effective arrival rates of route R_D and route R_I become αλand (1 − α) λ, respectively.

In case of c = 1, the arrival rate λ and the service rate μ respectively correspond to the inversed values of the interval of arrival time T_in and the staying time T_sp ($\lambda ={T}_{\mathrm{in}}^{-1},\mu ={T}_{{sp}}^{-1}$) as long as we ignore the effect of the congestion. In this approximation, we regard ${T}_{\mathrm{in}}^{-1}$ and ${T}_{\mathrm{sp}}^{-1}$ as the arrival rate λ and service rate μ as for the $M/M/1$ queue. Because the number of occupied sites does not suppress the maximum capacity, U_D and U_I in each route can be expressed by the classical $M/M/c$ queueing model, as follows:

$$\begin{eqnarray}&&{U}_{{\rm{D}}}=\min \left[{N}_{{\rm{D}}},\displaystyle \frac{\alpha {T}_{\mathrm{sp}}}{{T}_{\mathrm{in}}}\right]\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{U}_{{\rm{I}}}=\min \left[{N}_{{\rm{I}}},\displaystyle \frac{(1-\alpha ){T}_{\mathrm{sp}}}{{T}_{\mathrm{in}}}\right]\end{eqnarray} \tag{ 12 }$$

Here, N_Dand N_I indicate the number of parking sites in the domestic lane G_D and that in the international lane G_I, respectively. From equation (10), equation (11), and equation (12), we obtain the average number of occupied sites U_T as

$$\begin{eqnarray}&&{U}_{{\rm{T}}}=\min \left[{N}_{{\rm{D}}},\displaystyle \frac{\alpha {T}_{\mathrm{sp}}}{{T}_{\mathrm{in}}}\right]+\min \left[{N}_{{\rm{I}}},\displaystyle \frac{(1-\alpha ){T}_{\mathrm{sp}}}{{T}_{\mathrm{in}}}\right]\end{eqnarray} \tag{ 13 }$$

Figure 9 shows the comparison of the simulations with the approximations by the classical $M/M/c$ queueing models for different cases of the arrival time interval T_in between 2 min and 4 min. The solid lines with the filled symbol indicate the simulations. The dashed lines with the white-colored symbols show the approximations obtained using the $M/M/c$ queueing model. 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.

Zoom In Zoom Out Reset image size

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 10 shows the three-dimensional contour plot of the dependencies of the throughput Q on the distribution probability $\alpha$ 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 [25–27]. On the other hand, the dependence of throughput Q on arrival time interval T_in does not reach the plateau in between 4 and 2min 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.

Zoom In Zoom Out Reset image size

By devising a well-balanced model between theory and application, we have made it possible to quantitatively investigate the dependence of throughput on the usage ratio of parking sites on realistic systems. At the same time, we have found the existence of an optimal value that maximises the throughput by comparing simulations with queueing theory. Although this finding might look easy for some readers once it has been done, it surely brings many suggestions to practical problems. E.g., in the studies on air-ground transport systems, it is crucial to examine throughput on a system from every angle to enhance its turnover rate. Our simulation results provide air traffic controllers with a new perspective of maximising performance by optimising the usage ratio of parking sites which is scaled by the parameter α in our model. Furthermore, all the because of the simple mechanism on the throughput reduction, we may apply our research findings not only to junction flows consisting of multiple paths with different capacities as observed in rivers but also to the traffic flow of insects such as ants [30]. It can be said that the scope of applications ranges across a broaden area of physics.