Controlling Aircraft Inter-Arrival Time to Reduce Arrival Traffic Delay via a Queue-Based Integer Programming Approach

Abstract

Despite the importance of controlling the inter-arrival times of flights to propose strategies for efficient arrival management by the Arrival Manager (AMAN), the specific guidelines of such adjustments and their effect on reducing delays have not been explicitly considered. Accordingly, this paper proposes a novel approach, which integrates the $G_t/GI/s_t+GI$ time-varying fluid model and nonlinear integer programming to flatten the arrival rate at terminal gates. This, in turn, is achieved by minimizing the variance in inter-arrival times by penalizing any excessive change in arrival time, considering operational constraints. The results for Tokyo International Airport show potential to significantly reduce arrival traffic delays by minimizing said variance. This study may also spawn subsequent work, which builds a queuing network comprising upstream and terminal airspace and demonstrates the scope to reduce delays in the terminal airspace by controlling inter-arrival times at the upstream airspace.

1. Introduction

A significant increase in air traffic has caused severe congestion at airports and the surrounding airspaces and will remain in issue in the future. Japan Aircraft Development Corporation (2021) [1] forecasts that the number of passenger and cargo jets in service will experience 2.3% and 2.0% growth per year, respectively. This trend has also meant increasing delays. According to a joint report by Eurocontrol and FAA (2019) [2], in 2017, about Air Traffic Management (ATM)-related delays for departing flights, the total delay per flight was 1.73 min and 2.06 min in Europe and the U.S., of which the en-route-related delay per flight was 0.89 min and 0.35 min, respectively. This showed that in the U.S., 82.8% of delays were from terminals and 17.2% from en-route airspace, and severe weather was the main factor for both. Conversely, in Europe, 48.4% were from terminals and 51.6% from the en-route sector, and the Air Traffic Control (ATC) capacity and staffing constraints were no less dominant than adverse weather was. Note that the current situation with regard to air traffic in Japan is more similar to that in Europe, considering the policy of ATC and its consequences.

Relevant to the capacity constraint, the major factor behind any delay is the imbalance between supply and demand, which prevents air traffic controllers (ATCos) from matching controllable air traffic flow with existing facilities. To solve this problem and manage air traffic efficiently, Arrival Manager (AMAN), a system for supporting metering and vectoring, has been developed and implemented. Using AMAN, ATCos check the Estimated Time of Arrival (ETA) and the recommended time of arrival of flights and instruct some of them to reduce speed to prevent overloaded airspace intervals. Considering the AMAN mechanism, a guideline for controlling inter-arrival times for flights should be issued to reduce delays. This control corresponds to minimizing the variance in inter-arrival times or smoothing the arrival rate for the time horizon. However, Sekine et al. (2022) [3] indicated that previous studies focusing on speed control tend to be generally applicable to various airports and the surrounding airspaces and fail to pinpoint a viable strategy to control the speed of each flight. In response, rule-based simulation via cellular automation, multi-objective optimization and a decision tree were incorporated to determine such a strategy. Nevertheless, the issue of inter-arrival times for control was still not explicitly and quantitatively addressed, although Itoh et al. (2019) [4] suggested that a smaller variance in inter-arrival times could mitigate congestion and be realized by trajectory-based operation from a qualitative perspective.

Accordingly, this paper demonstrates potential to alleviate crowding by controlling the inter-arrival times of flights to flatten arrival rates at congested periods. In concrete terms and using Tokyo International Airport as a case study, we analyze the current terminal airspace to calculate delays in crowded intervals and verify the benefit of optimizing arrival times on delay reduction, taking operational constraints into consideration. This paper is a subsequent work of Higasa and Itoh (2022) [5], which indicates the potential scope to reduce delays by replacing a current arrival flow with a more orderly virtual flow, featuring a smaller variance in inter-arrival times at the en-route airspace of Tokyo International Airport. This study, which focuses on the actual (inter-)arrival times of flights, can directly boost AMAN development by reflecting delay mitigation strategies within the system as part of efforts to generate recommended arrival times.

The paper tackles the problem by combining two types of modeling approaches: the queuing model and mathematical programming. Concretely, the $G_t/GI/s_t+GI$ queuing model, which can consider the time-varying characteristics of the general arrival process, is applied to calculate delay times. It can check the effect of arrival intervals with less deviation and is adopted with the future work of building the queuing network in mind. Moreover, as a pioneer work to explicitly include variance in inter-arrival times for an objective function, nonlinear integer programming is employed, while penalizing any excessive changes in arrival time. The optimization model can consider operational constraints to discuss operationally feasible solutions. The problem with less than 200 decision variables is enough to be formulated as a flight-by-flight model, whose computation can become complex in the large-scale network compared with the flow-based model considering a set of flights together (see Sandamali et al. (2020) [6], in which a two-stage Air Traffic Flow Management framework is proposed with the latter as the first stage and the former as the second stage to deal with the potential increase in traffic). To the best of our knowledge, no combination of two modeling approaches has ever been determined for the purpose of improving AMAN by controlling inter-arrival times, although Jacquillat and Odoni (2015) [7] developed an integrated approach to jointly optimize flight schedules on a strategic level and airport capacity on a tactical level with the $M\left(t\right)/E_k\left(t\right)/1$ queuing model for airport queues. Note that this method can be applied to building queuing networks and for verifying that the control at the boundary of en-route airspace can mitigate congestion in terminal airspace (described in further detail in Section 5). Given the lack of prior research on jointly optimizing both airspaces, this study also contributes in that regard, although excluding the en-route airspace.

Here we review the literature on the queuing model and mathematical programming for Air Traffic Management. The queuing model has been employed to model queues in an air traffic flow to estimate operational delays and can consider two types of flow characteristics: stochasticity and nonstationarity (Shone et al., 2021 [8]). It can also test any improvement in ATC or AMAN strategies/tactics by focusing on the arrival process of the system. Among various possibilities for modeling the process, the homogeneous/nonhomogeneous Poisson arrival process has long been used (e.g., Bäuerle et al., 2007 [9], Itoh and Mitici, 2020 [10], Wang et al., 2018 [11]). However, the Poisson model has recently been questioned, since the arrival process cannot be considered appropriate for a congested arrival flow with arrivals successively adjusted by ATCos (Caccavale et al., 2014 [12]). Moreover, the model cannot confirm the positive effect of accurately anticipating arrival time and controlling it with a developed management system because the variance in inter-arrival time is the same as its mean, according to the Poisson process. In response, some researchers have suggested the Pre-Scheduled Random Arrivals (PSRA) model, which presumes the actual arrival time of each customer to constitute the sum of pre-scheduled arrival times and variance from it, following the independent identical distribution. From Guadagni et al. (2011) [13] as pioneer work, the PSRA model has been applied to a range of air traffic situations (see Caccavale et al. (2014) [12], Gwiggner and Nagaoka (2014) [14] and Lancia and Lulli (2020) [15] for instance). Note that these studies mainly focus on fitting the model to an air traffic flow and pay little attention to verifying the effect of less-deviated arrivals on delays.

Mathematical programming is another possibility with which to consider optimization problems from strategic and tactical perspectives. As Agustin et al. (2010) [16] described, most of the literature has focused on airports rather than sectors, because the problems were originally and mainly studied in the U.S., where delays mainly occur in the airport. As with Japan, where ATCos allocate metering delay to the en-route phase and the remaining delay to the ground (Gwigger and Nagaoka, 2014 [14]), the ground delay program has been conducted to let departure aircraft wait for some time to alleviate crowding in the destination airport. The ground holding problem has been formulated for single/multiple airports relevant to the ground delay program (Dell’Olmo and Lulli (2003) [17], Mukherjee and Hansen (2007) [18], Zhang et al. (2007) [19] for instance). The slot allocation problem for arriving flights in airports, needing to match the slot demand by airlines and airport capacity, has also been considered (e.g., Zografos et al., 2012 [20], Zografos and Jiang, 2019 [21], Ribeiro et al., 2018 [22] and 2019 [23]). Contrary to those studies, however, few studies have addressed en-route congestion. Barnier et al. (2001) [24] presented a sorting model for the slot allocation problem, where capacity constraints can always realize underloaded intervals; this is incorporated in our model (see Constraints (4) in Section 3.1). Flener et al. (2007) [25] considered relieving complex air traffic in the target airspace by the acceleration and deceleration of flights before entering the area. See Ball et al. (2003) [26], Lulli and Odoni (2007) [27], and Ivanov et al. (2017) [28] for another example. Then, as a recent work focusing on the concept of smoothing demand by speed control, Rosenow et al. (2022) [29] formulated the integer programming problem to minimize fuel consumption and capacity violation by controlling the arrival rate of 26 long-range flights arriving at the terminal airspace of Singapore Changi Airport. The problem is based on a work shift scheduling concept, and they found the scope to minimize it and prevent overload by adjusting speed three to four hours before entering the sector. We have reviewed the mathematical programming approaches, and found that none of them explicitly focused on minimizing the variance of inter-arrival times while considering operational constraints.

The rest of the paper is ordered as follows. Section 2 corresponds to the present analysis of terminal airspace. In Section 2.1, aspects of the $G_t/GI/s_t+GI$ queuing model, such as the assumption, parameters and functions, are described, whereupon a methodology to simulate airspace accurately via a queuing model using actual data is proposed. This approach is generally applicable to wide-ranging arrival flow around various airports. In Section 2.2, we discuss how to set up parameters in the $G_t/GI/s_t+GI$ queuing model based on real data for flights arriving at Tokyo International Airport in the terminal airspace during the evening hours when congestion peaks. In particular, how to determine capacity should be carefully considered, given that the balance between arrival rate and capacity is what dictates system behavior, and capacity setting can be relatively more arbitrary than that of other parameters. In Section 2.3, the calculation by the model indicates that an excessive arrival rate can cause a very considerable delay over four periods in the terminal airspace of Tokyo International Airport. Section 3 corresponds to the replacement of the arrival flow with an updated one, Section 3.1 presents the methods and Section 3.2 discusses applications. In Section 3.1.1, a nonlinear integer programming problem is formulated that minimizes the sum of the variance of inter-arrival times of entering aircraft and arrival time adjustment for all flights as a penalty. The constraints are for smoothing the peak of arrivals to prevent the overload of the airspace. Moreover, operationally feasible constraints are additionally considered in Section 3.1.2. In Section 3.2.1, we apply the approach to air traffic flow entering the terminal airspace to demonstrate that an early control of inter-arrival times can be effective for drastically reducing delay, leveling a peak in an originally congested period. In Section 3.2.2, we impose more realistic constraints on the extent of early arrival and adjustment, the result of which implies the limit of controlling the flow at the boundary of terminal airspace from an operational point of view. After the discussion on the distribution of inter-arrival time for each adjustment scenario in Section 4, we propose future work, which will consider more front-loading control of air traffic flow in en-route airspace using a network model. Concluding remarks are given in Section 5.

2. Estimating Arrival Queue in Terminal Airspace

2.1. Queue-Based Modeling of Aircraft Arrival Traffic Flow

The conceptual diagram of the $G_t/GI/s_t+GI$ queuing model is presented in Figure 1. It comprises the arrival process characterized by a time-varying arrival rate $\lambda \left(t\right)$($G_t$), the service process characterized by an independent and identical distribution of service-receiving time $g\left(x\right)$ ($GI$) and a time-varying service capacity $s\left(t\right)$ ($s_t$), and the independent and identical distribution of the abandon time, namely when the system is left before service is rendered, $f\left(x\right)$ ($+GI$).

This paper applies the deterministic fluid model developed by Liu and Whitt (2012) [30], which approximates the stochastic queuing model. In the fluid model, the functions and delay time, which represent the system state, are uniquely determined by the given parameters and initial conditions. Note that the service discipline is first-come–first-served. The input parameters are as follows:

• $\lambda \left(t\right)$: arrival rate at time t;
• $g\left(x\right)$: probability distribution function of service time x;
• $f\left(x\right)$: probability distribution function of abandon time x;
• $s\left(t\right)$: service capacity at time t.

The output functions, including delay time, are as follows:

• $Q\left(t\right)$: number of aircraft in a queue at time t;
• $B\left(t\right)$: number of aircraft in service at time t;
• $\sigma \left(t\right)$: number of exiting aircraft per unit time at time t;
• $v\left(t\right)$: expected delay time of an airplane entering the airspace at time t.

Figure 2 is the methodology used to analyze the current airspace with given $s\left(t\right)$, with the top left box representing the given flight data. We obtain the time of passage at each point, provided by the Trajectorized Airport Traffic Data Processing System (TAPS) (Ministry of Land, Infrastructure, Transport and Tourism [31]). These data then elicit the arrival and exit time of each aircraft, whereupon an overview of airspace in the real world can be determined by counting $c\left(t\right)$, the actual number of aircraft in the airspace. We can then roughly predict the time-varying behavior of the airspace, such as when the area has a large quantity of traffic flow and when delays are expected. Conversely, parameters in the $G_t/GI/s_t+GI$ queuing model can be set up according to the data to calculate functions and delay time. Then, to check whether we can develop the model adequately, we must validate it with the following equation:

$$c\left(t\right)\approx X\left(t\right)=B\left(t\right)+Q\left(t\right).$$

(1)

Equation (1) means that the number of aircraft in real airspace is nearly equivalent to the sum of the number of aircraft in a queue and in service in the model (see Figure 1, where $B\left(t\right)=Q\left(t\right)=2$). It ignores the amount of arrival, exit, and abandonment, which are proportionally smaller than $c\left(t\right)$, $B\left(t\right)$, and $Q\left(t\right)$. The calculation by the model should be repeated to ensure that Equation (1) can be properly satisfied while tuning parameters. Combining data analysis with queuing theory can elicit more insight into the time-varying characteristics of the airspace, and this process can also be applied to various types of arrival flow at general airports.

2.2. Stochastic Features in the Queuing Model

We obtained flight data for 25 days with 4700 airplanes, limiting the scope to flights arriving in the area after 17:00 and exiting (or leaving) it before 22:00. This period was the most congested (Itoh and Mitici, 2020 [10]). In this paper, we focus on the terminal airspace, 10–80 nm away from Tokyo International Airport (see the orange region in Figure 3), surrounded by the adjacent en-route airspace.

The parameters are determined as follows. Firstly, we set the averaged airspace arrival rate, denoted as ${\lambda }_{av}\left(t\right)$, which is the mean value of ${\lambda }_i\left(t\right)$ for each time t, and the arrival rate of each day i ($1\le i\le 25$). $av$ corresponds to the present state of the terminal area for the remainder of the paper. We assume that the arrival rate is right continuous with left limits and is constant for the time interval of arrival between preceding and subsequent flights, being the reciprocal of inter-arrival times. Accordingly, ${\lambda }_{av}\left(t\right)$ and ${\lambda }_i\left(t\right)$ are step functions. Figure 4 is the boxplot of ${\lambda }_i\left(t\right)$, showing the time/day-dependent behavior of the arrival flow. Secondly, we assume identical service and flight times, defined as the difference between the passage time of each side of the airspace. Although this only technically applies when no system queue exists, it remains broadly valid given that underloaded intervals prevail over overloaded intervals in airspace. Figure 5 is the service time distribution $g_{av}\left(x\right)$. In this paper, we presume $g_{av}\left(x\right)$ to be $N(1233.4,183.5^2)$, referencing Itoh et al. (2019) [4]. Thirdly, the abandon time distribution $f_{av}\left(x\right)$ is assumed to be $N(1233.4\times 1.4,183.5^2)$. The left-hand tail of $f_{av}\left(x\right)$ slightly overlaps the right-hand tail of $g_{av}\left(x\right)$, reflecting the rarity of emergency circumstances, such as the runway closing.

Unlike the parameters mentioned above, $s\left(t\right)$ can have a higher degree of freedom for its value in this process. Since even a few servers suffice to accommodate light traffic, but more are needed when the traffic intensifies, the ratio between $s\left(t\right)$ and $\lambda \left(t\right)$ must reflect the reality as accurately as possible. With that in mind, we have to check Equation (1) while changing $s\left(t\right)$ to accurately assess the airspace, since an excessively large/small $s\left(t\right)$ can result in the queue or delay time being under- or over-estimated.

Table 1. RMSE for three types of $s\left(t\right)$.

$\mathit{s}\left(\mathit{t}\right)$	RMSE
15.0	1.09
16.0	1.10
17.0	1.11

is the RMSE (Root Mean Squared Error) of the model for some fixed $s\left(t\right)$, defined by

$$\sqrt{\frac{1}{T+1}}\parallel \mathbf{X}\left(t\right)-\mathbf{c}\left(t\right)\parallel =\sqrt{\frac{1}{T+1}\sum _{i=0}^T{(X_i-c_i)}^2},$$

where $0\le t\le T$, $X_i=X\left(i\right)$ and $c_i=c\left(i\right)$. The possible candidates of $s\left(t\right)$ are selected based on the history of $c\left(t\right)$. From this result, we can deduce that $s\left(t\right)=15.0$ would provide the most precise analysis of airspace and that the system can accommodate all flights with constant capacity. In this paper, with the optimization discussed below in mind, we assume the constant capacity $s_{av}\left(t\right)=15.0$, which can simulate airspace with high accuracy. Note that $s\left(t\right)=14.0$ is infeasible in the model.

2.3. Arrival Delay Time in the Queue

For the calculation in the fluid model, we used Julia as the programming language and a computer with a 3.2 GHz 6 Core Intel Core i7 processor and 64 GB 2667 MHz DDR4 memory. The computational time is within 2 min for calculating the functions in the fluid model. The computational grid of the fluid model is orthogonal with 1674 × 252 points ($0\le t\le T=1673$, $0\le x\le E\left[g\right(x\left)\right]+7V\left[g\right(x\left)\right]=251$). Figure 6a shows the histories of ${\lambda }_{av}\left(t\right)$ [/10 s], $Q_{av}\left(t\right)$, $B_{av}\left(t\right)$ and ${\sigma }_{av}\left(t\right)$ [/10 s], whereas Figure 6b is on $v_{av}\left(t\right)$ [×10 s] and $X_{av}\left(t\right)=B_{av}\left(t\right)+Q_{av}\left(t\right)$, under $s_{av}\left(t\right)=15.0$, $g_{av}\left(x\right)$ and $f_{av}\left(x\right)$. $c_{av}\left(t\right)$ is also shown with them. Note that the initial time $t=0$ corresponds to 17:00:00 and the end time $t=1673$, namely the time the final aircraft arrives, corresponding to 21:38:50 in the real world. The mean value, variance, and maximum value are listed in

Table 2. Representative values of functions in the actual terminal airspace.

$\mathit{s}\left(\mathit{t}\right)=15.0$	$\mathit{E}[·]$	$\mathit{V}[·]$	max$\{·\}$
$Q\left(t\right)$	0.15	0.12	1.95
$B\left(t\right)$	13.3	5.7	15.1
$\sigma \left(t\right)$[/s]	0.010	$9.2\times {10}^{-5}$	0.014
$v\left(t\right)$[s]	12.4	77.0	154.0
$X\left(t\right)$	13.5	6.4	17.0

. The delay time occurs over four periods in the current terminal airspace: 17:45–17:58, 18:07–18:22, 20:21–21:00, and 21:24–21:36. In particular, congestion tends to peak during the period between 20:00 and 21:00, with $max\left\{v_{av}\left(t\right)\right\}$=157.3 s. These delays occur due to the excessive arrival rates at 17:30–17:40, 18:00–18:20, 20:10–20:30, and 21:20–21:30, respectively, all of which exceed 45 flights in terms of the number of arrivals per hour. This excessive arrival traffic should be dispersed by modifying the arrival time to the neighborhood time.

3. Controlling Inter-Arrival Time for Minimizing Arrival Delay

3.1. Optimization Model

3.1.1. Formulation of a Nonlinear Integer Programming Problem

As mentioned above, highly concentrated arrivals should be modified by changing their arrival times. This corresponds to the smoothing of inter-arrival time since the (piece-wise) arrival rate is the inverse of inter-arrival time. Subsequently, we develop the integer programming problem.

Firstly, the given parameters are set up:

• $t\in \mathbb{N}\cup \left\{0\right\}$: time horizon ($0\le t\le T$);
• $j\in \mathbb{N}\cup \left\{0\right\}$: j th-arriving flight ($0\le j\le J$);
• $b_j\in \{0,1,...,T\}$: arrival time of flight j before adjustment ($b_0=0,b_j=T$, see Figure 7);
• $s\in \mathbb{N}$: airspace capacity;
• ${\mu }_1\in \mathbb{N}$: flight time;
• ${\mu }_2\in \mathbb{R}$: mean inter-arrival time for all aircraft to be achieved ($=T/J$).

Flight 0 is for convenience in writing. Note that airspace capacity s, determined by the result of the queuing model, is interpreted as the maximum number of flights permissible in the airspace in the problem, although it technically corresponds to the number of flights in service in the queuing model. Moreover, we assume a constant flight time ${\mu }_1$ for all flights, as Rosenow et al. (2022) [29] adopt, while current arrival flow experiences various flight times for each aircraft. This study presumes ${\mu }_1$ as the mean flight time.

Secondly, we consider the variable corresponding to $b_j$.

• $a_j\in \{0,1,...,T\}$: arrival time of flight j after adjustment ($a_0=0,a_j=T$, see Figure 8)

We move the arrival time of each flight j so as to minimize the deviation in inter-arrival time from an ideal value ${\mu }_2$, namely, the variance in inter-arrival time. This is achieved using the decision variable $d_j$ representing the inter-arrival time.

• $d_j(=a_{j+1}-a_j)$$\in \mathbb{N}$: inter-arrival time between flight j and $j+1$ after adjustment ($0\le j\le J-1)$ (see Figure 8)

Since we assume an unchanged order of arrival, $d_j>0$. We set not $a_j$ but $d_j$ as the decision variable due to the computational tractability for theoretical reasons and interpretability for AMAN for practical reasons.

Moreover, given the penalty for excessive adjustments, such as operational inconvenience for airlines and extra fuel consumed when accelerating and decelerating, the sum of arrival adjustment should also be minimized. This penalty term can reduce computational effort, which means that combining arrival times with unrestrained changes elicits a large degree of freedom. Accordingly, we can formulate the nonlinear integer programming problem shown below:

$$Minimize\frac{1}{J-1}\sum _{j=0}^{J-1}{(d_j-{\mu }_2)}^2+\alpha \sum _{j=1}^J{Y}_j,$$

(2)

subject to

$$\begin{array}{c}\hfill \sum _{j=0}^{J-1}{d}_j=T,\end{array}$$

(3)

$$\begin{array}{c}\hfill \sum _{k=j}^{j+s-1}{d}_k>{\mu }_1(1\le j\le J-s),\end{array}$$

(4)

$$\begin{array}{c}\hfill -Y_j\le \sum _{k=0}^{j-1}{d}_k-b_j\le Y_j(1\le j\le J),\end{array}$$

(5)

$$\begin{array}{c}\hfill d_j\in \mathbb{N}(0\le j\le J-1),\end{array}$$

(6)

$$\begin{array}{c}\hfill Y_j\in \mathbb{N}\cup \left\{0\right\}(1\le j\le J).\end{array}$$

(7)

where $\alpha$ is weight.

$\alpha$ remains 1 throughout the paper for simplicity, which even elicits a significant reduction in delay time. Note that when $\alpha$ decreases, the variance in inter-arrival time carries more weight. The first term in the objective function (2) is the unbiased variance in inter-arrival time, while the second term constitutes the adjustment for all flights. The auxiliary variable $Y_j$ is used to avoid the absolute value, corresponding to $|a_j-b_j|$. Constraint (3) works as the conservation condition of the time horizon. Constraint (4) represents a situation where the number of aircraft in the airspace is less than or equal to the capacity, since we can deduce that $a_j={\sum }_{k=0}^{j-1}{d}_k(1\le j\le J)$ from the definition of $d_j$ and $a_{j+s}-a_j>{\mu }_1(1\le j\le J-s)$ represents it. Provided Constraint (4) is met, no delay is expected because it is equivalent to the number of aircraft $c\left(t\right)$ being less than or equal to capacity $s\left(t\right)$ from the perspective of the $G_t/GI/s_t+GI$ fluid model. Constraints (5) and (7) are established by introducing $Y_j$.

3.1.2. Nonlinear Integer Programming Problem with Operational Constraints

In the previous section, we assumed that any adjustments would be allowed unless the constraints were not satisfied. Consequently, congested periods saw earlier arrivals or significant changes in arrival times occur. However, early arrivals rarely arise in air traffic management in Japan today. Similarly, substantial changes in entry time are unpopular with both airlines/pilots and controllers, given the operational inconvenience and extra fuel consumption. Accordingly, we would like to add operational constraints to the problem. To restrict the adjustment, we can set constraints as follows:

$$Y_j\le n$$

(8)

which means that flights are instructed to move their arrival time by a maximum of n [×10 s] ($n\in \mathbb{N}$), whether this involves advancing or postponing the operation. Moreover, if we prohibit any early arrival of more than m [×10 s] ($m\in \mathbb{Z}$), we can impose a condition:

$$\sum _{k=0}^{j-1}{d}_k-b_j\ge -m$$

(9)

For example, if $b_j$ = 17:00 and m = 2 h, $a_j$ ≥ 15:00 and aircraft j can arrive earlier (15:00–17:00) or later (17:00–). If $b_j$ = 17:00 and m = 0 h, $a_j$ ≥ 17:00 and aircraft j cannot arrive earlier. Moreover, although it is not considered in this paper, if $b_j$ = 17:00 and m = $-2$ h, $a_j$≥ 19:00 and aircraft j can only arrive later (19:00–). Thus, if $m=0$, no flight enters the terminal airspace before the original arrival time.

3.2. Results

In this section, we demonstrate the delay-reduction effect by controlling inter-arrival times, with none of the operational constraints in Section 3.2.1 (given in Scenario 1 below) and with the constraints in Section 3.2.2 (given in Scenarios 2 to 6 below). To set scenarios, we consider severe constraints from two directions. We whittle down n and m from 10 [×10 s] independently until the model becomes infeasible and find that $n\le 8$ [×10 s] and $m\le -1$ [×10 s] are not allowed. Moreover, two kinds of constraints can only be satisfied simultaneously when $(n,m)=(9,9),(9,10),(10,9),(10,10)$. Accordingly, five types of scenarios are to be discussed (see

Table 3. Scenarios to be discussed. Scenario 0 is the original arrival flow, and - represents no limitation.

Scenario	Arrival Rate	Maximum Earlier Arrival [s]	Maximum Later Arrival [s]
0	${\lambda }_{av}\left(t\right)$	-	-
1	${\lambda }_{av,n=10}\left(t\right)$	-	-
2	${\lambda }_{av,n=9}\left(t\right)$	90	90
3	${\lambda }_{av,m=9}\left(t\right)$	90	-
4	${\lambda }_{av,m=6}\left(t\right)$	60	-
5	${\lambda }_{av,m=3}\left(t\right)$	30	-
6	${\lambda }_{av,m=0}\left(t\right)$	0	-

), compared with Scenario 1. Note that $n=10$ in Scenario 1 corresponds to $(n,m)=(10,10)$ and $n=9$ in Scenario 2 corresponds to $(n,m)=(9,9)$.

3.2.1. Delay Reduction by Control of Arrival Time in Terminal Airspace

For the actual arrival flow in the current terminal airspace, the arrival time is controlled using nonlinear integer programming to demonstrate the scope for the control to ease crowding. From the analysis mentioned in Section 2.2 and Section 2.3, we set the parameters as follows: $T=1673$ [×10 s], $J=188$, ${\mu }_1=123$ [×10 s], and ${\mu }_2=T/J=8.90$ [×10 s]. The parameters are optimized using Gurobi Optimizer with Python. We used the same computer as the one used for the fluid model, and the calculation time is within 2 s for optimizing the arrival sequence. Figure 9 represents the calculation result of Scenario 1, focusing on the adjustment for each flight j. If $a_j-b_j>0$ for flight j arriving at the airspace at time $b_j$, the flight is decelerated beforehand in more upstream airspace. Conversely, if $a_j-b_j<0$ for flight j arriving at the airspace at time $b_j$, the flight is accelerated beforehand. Figure 10 clarifies that an excessive number of arrivals are manipulated: some are advanced, while others are postponed. The updated arrival history is denoted as ${\lambda }_{av,n=10}\left(t\right)$. In the first cluster (17:30–17:40), most flights should be instructed to advance their arrival time. In the second cluster (18:00–18:20), controllers should ask pilots to slow down for later arrival. Then, for the most overcrowded arrival flow (20:10–20:30), flights are divided into two clusters, in the first and second half, respectively. The former group should move forward with their incoming schedules, while the latter should push off them.

The adjustment represents a value peaking at about 100 s. Finally, in the last cluster (21:20–21:30), most aircraft should enter the airspace behind the originally scheduled time. This spacing adjustment can reduce arrival flow deviation, conserving the total amount of arriving aircraft. This optimization reduces the standard deviation of inter-arrival time from 14.0 s (Scenario 0) to 12.5 s (Scenario 1).

Using the updated arrival history ${\lambda }_{av,n=10}\left(t\right)$, we recalculate the functions of the model. Figure 11a shows the histories of ${\lambda }_{av,n=10}\left(t\right)$ [/10 s], $Q_{av,n=10}\left(t\right)$, $B_{av,n=10}\left(t\right)$, and ${\sigma }_{av,n=10}\left(t\right)$ [/10 s], whereas Figure 11b is on $v_{av,n=10}\left(t\right)$ [×10 s] and $X_{av,n=10}\left(t\right)=B_{av,n=10}\left(t\right)+Q_{av,n=10}\left(t\right)$. This calculation is under $s_{av}\left(t\right)=15.0$, $g_{av}\left(x\right)$, and $f_{av}\left(x\right)$. Comparing Figure 6 and Figure 11, we find that no delay can be achieved by a more controlled arrival flow and that the flow leads to a more ordered exit flow, i.e., a landing flow. Smoothing arrivals helps avoid periods of heavy traffic, and it can be concluded that air traffic controllers should anticipate the arrival times beforehand as precisely as possible and instruct flights in the intervals to adjust before entering the terminal airspace.

3.2.2. Delay Reduction Considering Operational Constraints

For the relatively feasible scenarios in which arrival time is controlled by AMAN, the delay time is calculated by the $G_t/GI/s_t+GI$ fluid model.

Table 4. Calculation results for Scenarios 2 to 6 with Scenario 0 and 1.

Scenario	Standard Deviation of Inter-Arrival Time [s]	min{${\mathit{a}}_{\mathit{j}}-{\mathit{b}}_{\mathit{j}}$} [s]	max{${\mathit{a}}_{\mathit{j}}-{\mathit{b}}_{\mathit{j}}$} [s]	E[$\mathit{v}\left(\mathit{t}\right)$] [s]	max{$\mathit{v}\left(\mathit{t}\right)$} [s]
0 (${\lambda }_{av}\left(t\right)$)	14.0	-	-	12.4	154.0
1 (${\lambda }_{av,n=10}\left(t\right)$)	12.5	−100	100	0.0	0.0
2 (${\lambda }_{av,n=9}\left(t\right)$)	12.6	−90	90	0.0	0.0
3 (${\lambda }_{av,m=9}\left(t\right)$)	12.5	−90	100	0.0	0.0
4 (${\lambda }_{av,m=6}\left(t\right)$)	13.1	−60	120	0.0	0.0
5 (${\lambda }_{av,m=3}\left(t\right)$)	13.7	−30	150	0.03	6.68
6 (${\lambda }_{av,m=0}\left(t\right)$)	14.8	0	180	0.05	8.67

summarizes Scenarios 2 to 6 with Scenarios 0 and 1 already mentioned above. Firstly, the variance in inter-arrival time rises with excessive restrictions on early arrival. In particular, the arrival flow is less ordered if we prohibit early arrival (Scenario 6), although delays can be reduced under the problem structure in this paper (Constraint (4) guarantees the delay reduction). Figure 12 and Figure 13 illustrate the characteristics of arrival rate in Scenario 6, which fluctuates in the congested period by substantial adjustment (see Figure 14). There is a trade-off between fluctuation and delay reduction by maintaining an operational policy of prohibiting early arrival. Consequently, although all arrival histories, including Scenario 6, can alleviate crowding as well as Scenario 1 without any operational constraints, the possibility of postponing the arrival for a couple of minutes should be considered.

As mentioned in Section 3.1.1, Constraint (4) ensures that no delay is expected, but a little delay is observed if early arrival is allowed only within 30 s. This is probably attributable to differing assumptions between mathematical optimization and the queuing model. The actual flight time distribution $g_{av}\left(x\right)$ demonstrates a wide-ranging flight time, while we presume flight time to be constant in the integer programming. Based on these findings, we can deduce that Scenario 3, allowing for a 90 s earlier arrival at most, with an up to 100 s later arrival, constitutes the ideal control strategy, but may represent a daring attempt. From a tactical perspective and based on the concept presented in this paper, AMAN should calculate the adjustment time of each aircraft according to the ETA and instruct pilots to accelerate or decelerate. Hence, a precise ETA forecast is necessary using the latest technologies, such as the concept of four Dimension (4D) Trajectories-Based Operation (TBO) or data science using machine learning, as studied in the recent literature.

Although our approach is unique, and it is difficult to compare the results with those in conventional studies due to the difference in targeting airspace and operational conditions, here, we introduce one of the benchmark studies by Khassiba et al. (2019) [32], which proposed a two-stage (before and after entering the terminal airspace) stochastic programming model considering the uncertainty of passage time at the Initial Approach Fix. We could reduce the maximum delay time $max\left\{v\right(t\left)\right\}$[s] by 94.4 % (comparing the results in Scenario 0 and 6 in Table 4), which is the most conservative value in our study. On the other hand, the delay reduction was 73.5 % inside the terminal airspace of Paris Charles-De-Gaulle airport under the simulation-based experiments for 10 arrivals [32]. Note that 73.5 % was calculated under the condition of deviations in terminal-entry times from target times following $N(0,{\sigma }^2)$ ($\sigma =30\mathrm{s}$). The enhancement of the modeling fidelity under detailed assumptions, such as flight-dependent flight times, routes constraints, and rule-based control policies, is necessary for validating our results. As mentioned in the discussion section, fast-time simulation experiments will be conducted to confirm the delay reduction in future works.

4. Discussion

Controlling arrival intervals means a more orderly air traffic flow. This section focuses on the inter-arrival times themselves for each scenario. Figure 14 is the histogram of inter-arrival times for each scenario. Scenario 0 is the baseline, i.e., the original arrival flow without any intervention. Scenario 3 corresponds to the adjusted arrival flow with an arrival time change of up to 100 s and a maximum of 90 s earlier arrival. Note that Scenarios 1, 2, and 3 have the same shape for the histogram, although $d_j$ may differ in each case. Scenario 6 is the condition where no early arrival is allowed. For those scenarios, inter-arrival times have a maximum value of 140 s in periods with the sparsest arrivals around 19:20–19:30 (840–900 [×10 s], see Figure 13) and the behavior of flights in congested periods influences the distribution of inter-arrival time.

From this diagram, it is deduced that Scenario 3 is the best strategy for forming a narrow peak of inter-arrival times if the ATM can allow speed control for flights to shift their passage time to about 1.5 min. Before control, during congested periods, some flights would experience narrowly spaced rows, while others would be widely spaced, and their inter-arrival times would be varied. By smoothing the arrival rate by adjusting during crowded intervals, flights in the different periods are made uniform, since the flattened arrival rate equates to the flattened inter-arrival time. Conversely, the mode value of inter-arrival times is 100 s in Scenario 6 and the arrival flow becomes more varied, compared with not only Scenarios 1, 2 and 3 but also the current arrival flow in Scenario 0. This is due to the fluctuation in congested periods (see Figure 13), in which $d_j$ and $d_{j-1}$ can be totally different and there is no choice for some flights but to postpone their arrivals by more than 90 s.

The minimum inter-arrival time is 70 s in all scenarios, including Scenario 0 with the original arrival flow. This satisfies the minimum time-based separation at the Initial Approach Fix. According to Meyn and Erzberger (2005) [33], 72 s is required for the 5 NM distance-based separation. This indicates that both radar separation (3 NM) and wake vortex separation minima are satisfied considering the arrival aircraft types at the Tokyo International Airport. Note that one of the extensions when analyzing the arrival flow of each day is to incorporate the type of each aircraft into the problem as constraints, restricting the maximum value of $d_j$.

Table 5. Objective function and its breakdown for each scenario.

Scenario	Objective Function	First Term	Second Term
0	-	-	-
1	263.6	1.6	262.0
2	263.6	1.6	262.0
3	263.6	1.6	262.0
4	274.7	1.7	273.0
5	331.9	1.9	330.0
6	440.2	2.2	438.0

summarizes the values of the objective function, first term (the variance in inter-arrival time), and second term (the adjustment for all flights). For Scenarios 4, 5, and 6, all the values, particularly the second term, become worse than those in Scenarios 1, 2, and 3 due to operational constraints. These limitations on early arrival generate postponement on a very large scale for flight arrival times during congested periods, although congestion is mitigated even in more deviated arrival flows, as guaranteed by Constraint (4) in Section 3.1.1. Note that weight $\alpha$ is constant at 1 and the second term dominates in this paper, with an imbalance emerging between the first and second terms. In the current problem structure, the penalty for excess adjustment is supposed to have more influence than optimization for a more orderly flow. Although even this setting can realize alleviated crowding, there is room to improve the programming problem by changing $\alpha$ or allowing an increase in mean inter-arrival time and so on. The latter modification means the time horizon is expanded backward (e.g., from 17:00:00–21:38:50 to 17:00:00–22:00:00) and we should investigate the potential for more reasonable operation by AMAN through this modification as a topic for future work.

Moreover, it is meaningful to further analyze the characteristics of the arrival flow for each day. This paper used the mean arrival rate by averaging 25-day arrival histories. Although the given assumption leads to the solution being generally applicable, it is important to clarify how the flight-time-dependent arrival rates have effects on the queuing model and integer programming in future work.

Another approach involves verifying the delay-reduction propagation effect. Figure 15 shows a conceptual diagram of a queuing network. We learned that the arranged arrival flow helps to reduce delays in the interested area (Higasa and Itoh, 2022 [5] for the en-route area by generating a virtual flow, and in this paper for the terminal area using the real flow). Accordingly, we will verify that controlling the arrival time of certain clusters at the border of en-route airspace can mitigate congestion in the terminal airspace. According to Figure 15, this corresponds to controlling $\lambda \left(t\right)$, i.e., stabilizing the exit flow $\sigma \left(t\right)$ and pouring ${\lambda }_T\left(t\right)$ into the terminal airspace to decrease delay time. If this holds, the exit flow from the terminal airspace ${\sigma }_T\left(t\right)$ also becomes stabilized, providing a more orderly landing flow, and the smoother landing can ease crowding on the surface or at departure in an airport. Furthermore, we will verify the effectiveness of earlier interval control in the en-route airspace, assuming the increasing volume of air traffic in the future. The proposed optimization method will be applied to a virtually created aircraft arrival flow, which mimics the increasing volume of future air traffic.

Finally, we will conduct fast-time simulation experiments using AirTOP simulator, with $b_j$ and $a_j$ as input. Future work will validate whether the simple assumptions mentioned above, such as the constant flight time in the integer programming formulation, reduce aircraft arrival delays under realistic operational conditions, considering aircraft-dependent flight time, and even increasing air traffic volume.

5. Conclusions

For efficient Arrival Traffic Management (ATM), the imbalance between demand and capacity should be remedied in the surrounding airspace. Considering the mechanism of the Arrival Manager (AMAN) displaying the Estimated Time of Arrival (ETA) and the Scheduled Time of Arrival (STA), the advanced control of flights to optimize their inter-arrival times can be an effective approach to solve this problem and mitigate congestion. However, little attention has been paid to the quantitative demonstration of interval control. Thus, we developed a methodology to optimize arrival flow so that inter-arrival times are less deviated while the excess adjustment of arrival time is restricted, combining the $G_t/GI/s_t+GI$ time-varying queuing model and nonlinear integer programming. The integrated approach clarifies that crowding in the terminal airspace can be mostly alleviated by advancing and postponing flights entering the area during the congested period, even with the operational constraints for early/late arrival being at most about 90 s in the case of Tokyo International Airport.

Note that airlines and controllers may experience a much heavier workload than before because they must maintain close contact in order to accelerate or decelerate flights in cruise phase in the en-route airspace. Since the terminal airspace is far more crowded to enter than en-route airspace, more aggressive intervention is required for potentially overcrowded arrival flow, as identified by our additional analysis, which will be discussed in future work. In this future work, we will apply this process to relatively sparse en-route airspace to examine whether adjustment elicits a similar effect when controlling inter-arrival times at en-route gates and when the adjustment is modest. Furthermore, future work will include demonstrating the positive effect of advancing adjustment points by building a queuing network. It may also help to propose operational strategies for integrating AMAN (arrival management), SMAN (surface management), and DMAN (departure management).