Effectiveness of Aircraft Inter-Arrival Control in Upstream Traffic Flow via a Combined Tandem Fluid Queue Model and Integer Programming Approach

Conventionally, several studies indicated that controlling aircraft arrival time in the en-route airspace mitigates arrival aircraft congestion in the terminal airspace. Further research is required to clarify how to leverage this idea to design an air traffic management system, a so-called Extended Arrival MANager (E-AMAN), to reduce the arrival traffic flow while assisting air traffic controllers and boosting their effectiveness quantitatively. Under these circumstances, this research proposed aircraft inter-arrival time control within the en-route airspace and clarified its effectiveness in reducing arrival delay based on mathematical modeling and simulation evaluations. In this paper, we developed the $G_{t}/GI/s_{t}+GI$ tandem fluid model to analyze the time-varying delay time of flights in both en-route and terminal airspace and demonstrated the effect of inter-arrival time control in the upstream arrival traffic flow in the en-route airspace, combining the model with the nonlinear integer programming problem. The calculation results for 3,074 aircraft over 21 days, arriving at Tokyo International Airport between 17:00 and 22:00, show the possibility for the control to reduce the mean and maximum delay time for flights by 18.8% (5.0 s) and 16.5% (37.6 s) on average within the en-route airspace. Moreover, fast-time simulation by AirTOP is conducted to validate the control, revealing the scope to reduce mean and maximum delay times in the terminal airspace by 11.5% (36.5 s) and 19.2% (148.8 s) on average.


I. INTRODUCTION A. LIMITATION OF CONVENTIONAL ARRIVAL MANAGER
To match the existing facilities in airspaces and airports and the demand of arrival traffic flow, the Arrival Manager (AMAN) has been used as a decision support tool that aids the arrival controller's sequencing, scheduling, and metering in the Terminal Maneuvering Area (TMA, or terminal airspace) [1]. It has been developed in Europe over sev- The associate editor coordinating the review of this manuscript and approving it for publication was Md Asaduzzaman . eral decades [2] and has also been implemented since fiscal 2019 in Japan as one of the policies of Collaborative Actions for Renovation of Air Traffic Systems, CARATS [3], [4].
However, air traffic demand has been soaring worldwide and the traditional AMAN cannot handle the flow efficiently via the procedures mentioned above. This results in heavy congestion, not only within the terminal airspace but also in the surrounding en-route airspace. The Japan Aircraft Development Corporation [5] forecasts that over 2022-2041, the growth rate per year will be 2.7% and 1.8% for the number of passengers and cargo jets in service, respectively. Moreover, capacity may be violated within the en-route airspace to the same extent and as frequently as the terminal airspace in Europe [6] and Japan, which operates a similar management system to that in Europe. In Europe in 2017, [7] reported that regarding Air Traffic Management-related delays on departing flights, 48.4% and 51.6% of delays were from the terminal and the en-route airspace respectively, and Air Traffic Control (ATC) capacity and staffing constraints were as the main reasons as adverse weather. The need to control air traffic flow to match the capacity and demand has emerged, in both terminal and en-route airspace in Europe and Japan.

B. POTENTIAL BENEFITS OF EXTENDED-AMAN
Accordingly, the Extended-AMAN (E-AMAN) has been studied. In the E-AMAN, to ease crowding in both areas, the horizons are extended further upstream to arrange arrival traffic flow much earlier than before [8]. According to previous studies, controlling flights ahead of time by E-AMAN may help reduce delays. The work in the context of arrival traffic in the U.S. [9] indicated the scope for most delays to be transferred to the relatively sparse en-route airspace by assigning the Controlled Time of Arrival at 500 NM away from the airport. The work implied the potential to mitigate terminal crowding by earlier control. Moreover, [10] deduced that starting the scheduling and sequencing 2-3 hours before landing could efficiently carry flights to the landing phases. They proposed a two-stage stochastic programming problem (before/after entering the TMA) for 14 flights entering the terminal area of Paris-Charles-De-Gaulle airport between 6:10 and 6:30 AM, whereupon they found that the total deviation in landing time from ideal and unconstrained values could be reduced by 71%. Furthermore, [11] extracted optimal speed-control strategies at 150-200 NM from Tokyo International Airport by combining the Cellular Automatonbased simulator, multi-objective optimization, and decision trees. This work contributed to the future implementation of the E-AMAN, starting arrival spacing and sequencing more than 350 NM away from the airport.

C. CONTROLLING AIRCRAFT INTER-ARRIVAL TIME
Towards implementing E-AMAN, conventional studies focused on optimizing the Scheduled Time of Arrival (STA) and asked Air Traffic Controllers (ATCos) to minimize errors between the aircraft's Estimated Time of Arrival (ETA) and STA. However, these operations often increase ATCos's workloads under congested air traffic operations. Additionally, it causes difficulty for ATCos to minimize the error due to the ETA accuracy, especially under uncertainties in the actual operation. In response, [12], [13] focused on improving ATCo's tactical control strategies in the en-route airspace to reduce arrival delays at a single runway. They indicated that minimizing the variance in aircraft inter-arrival times in the en-route airspace following the First-Come First-Served arrival sequences could reduce delay time in the entire arrival aircraft at a runway via theoretical analysis by steady-state queuing models. However, the model had to be improved for further analysis because the arrival rate in the arrival air traffic is time-varying.
Our previous study [14] theoretically demonstrated the scope to reduce the delay by replacing an arrival flow with a more orderly flow with fewer deviating inter-arrival times at the en-route airspace of Tokyo International Airport applying time-varying fluid queuing models. However, in their work, they generated the virtual arrival flow by random sampling of inter-arrival times and did not verify the effect for the actual arrival flow. Accordingly, in subsequent work, [15] covered the actual arrival flow at the terminal airspace of Tokyo International Airport and demonstrated the effect. Instead of minimizing errors between STA and ETA, controlling the aircraft inter-arrival time is expected to achieve ATCofriendly operation; they adjust relative time-spacing between consecutive arrivals rather than attempting to control the arrival times of all aircraft.

D. MITIGATING DELAY PROPAGATION FROM THE UPSTREAM TO DOWNSTREAM
However, even the previous work did not demonstrate that controlling arrival intervals could ease crowding in the en-route airspace as well as terminal airspace. As mentioned above, congestion of the en-route airspace should be tackled with that of terminal airspace due to increasing demand. Moreover, as well as proving this effect, the concept of mitigating congestion in both upstream and downstream areas by controlling the upstream input rate, namely, input intervals, has not been discussed.
Note that it has been considered in the context of surface management on airports. The previous study [16] described a field test at Boston Logan International Airport in 2010, where controlling the pushback rate from gates on the part of controllers was found to reduce fuel by 12,250-14,500 kg with an increase of 4.4 minutes holding at gates. Another work [17] modeled the taxi-out process, comprising gates, ramps, and runways by a tandem fluid model to demonstrate the effect of controlling the pushback rate (two components: ramp queue and spot server and runway queue and runway server). See [18] and [19] for the subsequent work and [20] and [21] for the work based on the same philosophy.

E. RESEARCH OBJECTIVE AND STRUCTURE
This paper verifies the delay reduction propagation effect: namely how controlling inter-arrival times at the boundary of the en-route airspace can ease crowding in the upstream en-route airspace and accordingly in the downstream TMA. This work follows [15], which applied a combination of the G t /GI /s t + GI fluid model and nonlinear integer programming. The former technique is applied for analyzing current and updated airspaces with a more orderly flow, whereas the latter is for optimizing arrival spacing. For Tokyo International Airport as a case study, two consecutive airspaces are set out and the tandem-queuing network model is applied, VOLUME 11, 2023 referring to [17]. Moreover, nonlinear integer programming can be used to explicitly consider how to minimize any variance in arrival intervals as the form of the objective function. To the best of our knowledge, this is the first study to demonstrate the scope for earlier control of flights by extending the AMAN boundary to reduce delays in the en-route and terminal airspace of congested airports by building a queuing network with the general output/input process, explicitly considering the optimization of inter-arrival times of the aircraft by nonlinear integer programming. This approach is supposed to directly offer insights into how E-AMAN is developed.
The rest of the paper is as follows. In Section II, the single and tandem-queuing model is explained. Section II-A summarizes the G t /GI /s t + GI queuing model and the general procedure for analysis by the model, while Section II-B focuses on two consecutive queuing models and the process of demonstrating delay reduction in both airspaces. In Section III, delay times are calculated for two current airspaces. After describing the data used and the aircraft targeted in Section III-A, time-varying characteristics of arrival air traffic flow are considered for the en-route and the terminal airspace in Section III-B, mainly focusing on delay times. In Section IV, inter-arrival times are optimized to verify the delay reduction in both areas. Parameters, decision variables, objective function, and constraints of nonlinear integer programming are introduced in Section IV-A, followed by updated arrival traffic flows into the tandem-queuing network model to calculate the delay times in Section IV-B. Additionally, in Section V, the fast-track simulation software AirTOP is introduced in Section V-B and used to calculate flight times for original and updated flow in Section V-C. This is to complement the operational constraints in real-world operations such as separation minima (detail in Section V-A). Section VI-A compares the result of our combined model with that of previous studies, whereas Section VI-B compares and discusses the result of our macroscopic model and relatively microscopic simulation. After describing some possible future extensions in Section VI-C, concluding remarks are presented in Section VII.

II. DESCRIPTION AND FORMULATION OF MODEL A. SINGLE QUEUING MODEL
In this section, we review the description of the G t /GI /s t + GI queuing model and the methodology used to analyze en-route and terminal airspace based on [15]. Fig. 1 illustrates the concept of the G t /GI /s t + GI queuing model. This model comprises the time-dependent general arrival process (G t ), the general independent service process (GI ) with the time-varying service capacity (s t ), and the general independent abandonment process (+GI ). The arrival process is characterized by a time-varying arrival rate, which influences the number of aircraft in service when there is no queue and that in a queue when there is a delay. The service time of aircraft is independent and identically dis- tributed and a time-dependent service capacity defines the maximum number of aircraft that can receive the service simultaneously at each time. Finally, abandonment refers to leaving the system before the service when there is a queue facing congestion and its time is independent and identically distributed. Our study applies a deterministic fluid model with the first-come-first-served discipline advanced by [22], approximating the original stochastic queuing model. In the fluid model, by setting parameters and initial conditions (in this case, the system is empty at the starting point), the history of the state is uniquely calculated in the form of functions including the delay time. Parameters in this model include: • w(t): actual delay time of the airplane waiting at the front at t Note that we focus not on w(t) but on v(t) for the discussion below, since knowing the potential delay time at the stage of entry is a higher priority for aircraft and controllers in an operational context. Moreover, in the rest of this paper, we adopt a flight-based delay time v j rather than a time-based delay time v(t) to compare performance between original and optimized flows. The main reason is to maintain integrity in comparison and validation phases with the trajectory-based simulation (details in Section V). In fact, for each flight j arriving at airspace at t j .
The method of analyzing the present airspace is summarized in Fig. 2. Below are the procedures when the flight data is given (the details of the flight data described in Section III-A). From a real-world perspective, we can perceive the outline of the current airspace from the data with the arrival and exit times of each flight by counting the existing number of airplanes at every time (c(t), the top right box in Fig. 2). Although we can briefly grasp the behavior of airspace, including the expected period of any delay, the time-varying crowding the flights face is not quantitatively and qualitatively understood, such as when the delay peaks during an expected period of heavy congestion and to what extent it lags behind the heavy input flow. Then the queuing model can be used to bridge the gap between simplicity and a low degree of accuracy. Parameters in the G t /GI /s t + GI queuing model can be tuned according to the flight data (the bottom left box in Fig. 2), whereupon the functions and delay time are calculated for each time. Note that not all parameters can be valid. For example, if the calculation result for an excessive number of capacities shows no delay over the time horizon but the actual airspace is supposed to experience a period of overcapacity with delay times expected, we should not discuss fact-based operational improvement of the airspace based on the result of this model. The parametersetting process can eventuate in the problem of finding the most appropriate solution to the following equation: Equation (2) can bridge the actual flow of aircraft in the real world and that simulated by the model. It argues that the sum of the number of aircraft in a queue and service in the model approximates the actual number of flights as a baseline. It should technically consider the amount of arrival λ(t), exit σ (t), and abandonment, but they are sufficiently smaller than c(t), B(t), and Q(t). Parameters should be tuned repeatedly for (2) to be satisfied. With the appropriate configuration of airspaces, this method can be applied to various airports and the surrounding airspaces. Fig. 3 is a schematic diagram of delay propagation, which can be demonstrated by a tandem-queuing model. By the arrivals at the en-route airspace (λ E (t) is taken as the arrival rate, see Section III-B1 for a detailed definition), delays occur in the area (v E (t) as a delay time, see Section II-A), and flights exit the airspace (σ E (t) as an exit rate, see Section II-A). Note that the features on the en-route airspace are denoted as E for the rest of the paper. Subsequently, the output flow σ E (t) proceeds into the terminal sector as the input flow, i.e., σ E (t) = λ T (t), and the delay time v T (t) and exit rate σ T (t) are calculated in the same way. The features on the  terminal airspace are denoted as T for the rest of the paper.

B. TANDEM QUEUING NETWORK MODEL
In the tandem fluid model, by instructing congested flights to advance or postpone their arrival times for less deviated inter-arrival times between flights, i.e. by flattening the arrival rate shown in the solid line in Fig. 4, we investigate whether v E (t) is reduced and likewise v T (t) is also reduced as a result of the smoothed output/input flow σ E (t) = λ T (t) (both in Section IV-B).

A. AIRPORT AND DATA FOR ANALYSIS
In this section, we describe the current operation of Tokyo International Airport and the flight data used for this paper. Tokyo International Airport (RJTT) had been the fifth busiest airport in passenger number terms until 2019, the year before the COVID-19 pandemic [12]. Although arrivals come from both the southwest and north, the former outnumbers the latter about threefold [12]. Flights from the southwest comprise domestic routes in the western part of Japan and East and Southeast Asia [23] and have remained on the up, even after previous studies were published. In this paper, the en-route airspace is defined as the southwest area between 100-30 NM from a benchmark waypoint XAC (the blue region in Fig. 5, same as [14]). Red lines represent the trajectories of flights arriving at Tokyo International Airport from the southwest. Moreover, we define the terminal airspace as the area between 80-10 NM from Tokyo International Airport (the orange region in Fig. 6, same as [15]). Since XAC is located 50 NM away from the airport, these two airspaces can be deemed contiguous, and the boundary of the en-route airspace is 150 NM away from the airport, referring to [11], where the cruise speed is controlled.
The data contain the passage time of each flight at each distance from the reference point (XAC or Tokyo International Airport). For the en-route airspace, the Trajectorized En-route Traffic Data Processing System (TEPS) is used, whereas the Trajectorized Airport Traffic Data Processing System (TAPS) VOLUME 11, 2023  is for the terminal airspace in this case study, both of which are provided by the Japan Civil Aviation Bureau [24], [25]. Note that TEPS is used for generating arrival flow for both airspace and parameters for the en-route airspace in the fluid model, whereas TAPS is only for parameters for the terminal airspace in the model, rather than the arrival flow for the terminal airspace. This is because the arrival flow directly enters the terminal airspace after the en-route airspace in the tandem queuing model. We obtained actual flight data for 21 days from September 2019 to February 2020, focusing on the dominant flow from the southwest. Moreover, we concentrate on the flights that entered the en-route airspace after 17:00 and exited it before 22:00, which was the most congested window for arriving air traffic in Tokyo International Airport [12]. The total number of aircraft is 3,074 for 21 days.

B. AIRCRAFT DELAY TIME IN EACH AIRSPACE 1) ESTIMATION OF PARAMETERS FOR CALCULATION
In this section, for the current en-route airspace, parameters are set up to calculate functions in the G t /GI /s t +GI queuing model using TEPS data. The process below is also applied to the terminal airspace using the TAPS data, which is not described in this paper. We denote arrival flow of each day as . Note that the features on the original airspace are denoted as O and U is for those on the updated airspace by the optimization. λ O,E,i (t) (and λ U ,E,i (t)) are defined as the piecewise constant right continuous with left limits: if the preceding flight enters the airspace at t = t 1 and the subsequent one at t = t 2 for instance, Note that the arrival rate is aggregated for every ten seconds from flight data (unit: one second). This is to reduce the computational cost of the model and improve the readability of results. Figs  Secondly, the service time distribution g O,E (x) is assumed to be N (500.7, 49.7 2 ) for all 21 days with the outlier of flight time removed, the normal distribution whose mean and variance are by the maximum likelihood estimate, referring [15] and references therein (see Fig. 9). We presume that service time equates to flight time, defined as the passage-time difference between the entry and exit times of the airspace. Technically, it is untrue if there is an overload capacity with a queue, but it can work as an accurate approximation of the current airspace operated under capacity for most of the rush hour and guaranteed by (2). Moreover, we assume the abandonment time distribution f O,E (x) to be N (500.7×1.4, 49.7 2 ) for all 21 days, considering the scarcity of abandonment. Abandonment can occur in an emergency such as diversion, but it rarely does.
Finally, service capacity s(t) is selected to describe the behavior of airspace most accurately considering c(t) ≈ X (t) = B(t)+Q(t) from RMSE (Root Mean Squared Error)'s point of view (see [15]). The most appropriate s(t) is determined by comparing the sum of the error between c(t) and X (t) for all time t.

2) DELAY IN CURRENT EN-ROUTE AIRSPACE
After setting parameters as mentioned above, by using the G t /GI /s t + GI queuing model, delays in the current en-route airspace are calculated under g O,E (x) = N (500.7, 49.7 2 ) and f O,E (x) = N (500.7 × 1.4, 49.7 2 ). 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 2,667 MHz DDR4 memory. The computational time is within two minutes for calculating the functions in the model. The computational grid of the model is orthogonal with 1,800 × 85 points in the en-route airspace ( differs from the exit rate σ (t)) and the terminal airspace as well (the maximum of x differs from that of the terminal airspace due to difference in the service time distribution). Table 1 , it emerges that the fluid model can accurately simulate the airspace by s O,E,7 (t) = 6.0. As with the case of terminal airspace (see [15]), the en-route area has a symptom: the air traffic flow exceeds the capacity of air traffic controllers with a lag behind the arrival of a huge number of flights. In specific terms, rather than a momentary large entry, ongoing and chronic entries subsequently affect the behavior of the airspace. For instance, by comparing the last period of delay around 1400×10 s with the last-minute arrival rate, it is deduced that continuous entry of flights leads to an increase in delay with λ(t) exceeding 0.10/10 s, corresponding to 36 flights per hour. Note that the value of this threshold differs day by day, but the tendency itself holds for all days.
data. Note that we estimated from σ O,E,i (t) the time when flights from the en-route airspace traverse the boundary of the terminal airspace since the fluid model cannot output the actual flight time of each flight. The order of exit remains unchanged due to First-Come First-Served (FCFS) discipline. Table 2   instructed to follow ATCos take considerable time, part of which is the delay.
From the queuing model's perspective, as a result of the excessive demand, the number of aircraft in service B(t) becomes equal to the number of service capacity s(t), resulting in the increase of aircraft in a queue Q(t). In this overloaded regime, it is extracted that s(t) flights can smoothly fly in the airspace without metering or vectoring, whereas the rest of the flights are still instructed to follow such operations.
Comparing the results with those in Section III-B2 (see Fig. 10 and 11) reveals that the delays are propagated from the en-route airspace to the terminal airspace with a lag of about the mean flight time in the en-route airspace. In the example of Day 7, both delays are amplified in five consecutive periods with an overloaded regime in the first half of the five-hour time horizon and in the one period around 21:00 (1,440 × 10 s). This is the consequence of long flight time comprising the unrestricted flight time with deceleration for landing and the time to accept metering and vectoring, the latter of which explains why the number of aircraft in airspace c(t) exceeds the capacity s(t). The former is inevitable but the latter can be reduced by shifting those adjustments from the terminal airspace to the en-route airspace, or preventing the occurrence of delays at an early upstream stage. To prevent delays, not only at the en-route area but also the terminal area, the arrival time should be controlled to flatten demand by advancing or postponing the entry time, using mathematic programming. In the next section, as a concrete method, the nonlinear integer programming model will be formulated.

IV. CONTROLLING INTER-ARRIVAL TIMES FOR DELAY REDUCTION A. OPTIMIZING INTER-ARRIVAL TIMES AT BOUNDARY OF EN-ROUTE AIRSPACE
To control the air traffic flow, air traffic controllers should order the acceleration or deceleration of each flight to decrease the variance in inter-arrival times, defined as the difference between the entry time of preceding and subsequent aircraft, namely, the reciprocal of the piecewise arrival rate. This type of operation can be modeled as an integer programming problem. Although it is already formulated in [15], we review it here. Firstly, the given parameters are below: • s ∈ N : airspace capacity • µ 1 ∈ N : flight time • µ 2 ∈ R : mean inter-arrival time for all aircraft to be achieved (= T /J ) 0 th-arriving aircraft are set virtually but used to conserve the time horizon with the J th-arriving aircraft. Airspace capacity s is set according to the result of the queuing model in which the capacity is the maximum number of flights in service, but in the problem, it is interpreted as the maximum number of flights allowed to be in the airspace. We set the common flight time µ 1 , the mean value of flight-time distribution, for all flights (see [26]). µ 2 is a mean inter-arrival time for J flights, maintaining a constant before and after the adjustment, since the end time T remains unchanged.
Within the problem structure, some originally scheduled (or expected) arrival times b j are rearranged. Below is the revised arrival time of each flight: • a j ∈ {0, 1, . . . , T }: arrival time of flight j after adjustment (a 0 = 0 and a J = T , see Fig. 15) To ease the highly nonlinear function and explicitly handle inter-arrival times, we re-declare the decision variable to transform the problem into a simpler form.
• d j (= a j+1 − a j ) ∈ N: inter-arrival time between flight j and j + 1 after adjustment (0 ≤ j ≤ J − 1, see Fig. 15) Moreover, the auxiliary variable Y j is used to avoid the absolute value. Since a j = j−1 k=0 d k (1 ≤ j ≤ J ), the nonlinear integer programming problem is formulated as follows: Minimize subject to where α is weight. The first term in the objective function (4) is the variance of inter-arrival time, and the second term is the sum of the amount of adjustment for all flights. Without the second term, not only are computational costs high but also the result is not realistic such as the excessive change of arrival time. Provided that excess adjustments lead to operational inconvenience for airlines or extra fuel consumption VOLUME 11, 2023  due to acceleration and deceleration, the sum of arrival adjustments should also be considered in the form of an objective function as a penalty term. (5) and (6) are based on a 0 = 0, a J = T , and a j+s − a j > µ 1 (1 ≤ j ≤ J − s), which limits the number of aircraft in the airspace so that the demand is below the capacity for all times. This refers to the concept of the sorting model by [27]. (7) and (9) are about Y j , and (7) are from the processing of the amount of adjustment |a j − b j |.
Finally, (8) are from a j < a j+1 (0 ≤ j ≤ J − 1), which are for the order of arrival to remain unchanged.
Since we mainly focus on the demonstration the delay-reducing effect by incorporating the variance in the objective function rather than the detailed structure of the Pareto-optimal points, we set α = 1 and the notation is abbreviated for the whole paper (see [28] as a reference). We also consider the fact that the structure may be one of the extensions of our work given that the multi-objective optimization technique can provide more insights. Note that α characterizes the influence of the variance in inter-arrival times.

B. DELAY REDUCTION IN EACH AIRSPACE
In this section, inter-arrival times of flights are adjusted so that arrival rates are smoothed to reduce delays in the en-route and terminal airspace by solving the optimization problem.

1) FLATTENED ARRIVAL RATE BY OPTIMIZATION
Here we introduce the effect of optimizing arrival intervals using Day 7 as an example. The optimization is implemented by Gurobi Optimizer with Python in less than two seconds of computation, using the same computer as that for the fluid model. In Fig. 16, the horizontal axis is the originally scheduled (or ETA from air traffic control's perspective) arrival time of flight j on day 7, whereas the vertical axis is the amount of adjustment (or the difference between ETA and STA) for each flight j. a j − b j > 0 is the late-arrival, a j − b j < 0 is the early-arrival, and a j − b j = 0 means that the change is not required. For example, the flight originally supposed to arrive at about 18:56:47 (b 7,j = 7, 007 s) is instructed to postpone its entry by 238 s (19:00:45, a 7,j = 7, 245 s), whereas the flight due to arrive at 19:09:45 (b 7,j = 7, 785 s) is made to revise the time so that it can enter the airspace at  19:04:47 (a j = 7, 487 s). Fig. 17 compares the original and revised arrival rates, in a unit time of ten minutes. When the flow is expected to become congested, the first half of the cluster is instructed to cross the border in advance, whereas the postponement of the second is ordered. Consequently, the peak arrival rate declines over a five-hour rushed period and the fluctuation of the arrival rate is moderated. In other words, the arrival flow is flattened. Note that early arrivals are allowed in this paper, which is an ideal condition from an operational point of view.

2) DELAY REDUCTION IN EN-ROUTE AIRSPACE
In the phase of demonstrating the reduction effect of delay time, statistical hypothesis testing is conducted. We focus on 21 pairs of representative values of arrivals for 21 days (before and after the optimization for 21 days each). We applied a Wilcoxon signed-rank test for all sets, one of the non-parametric hypothesis testing methods for matched samples with an unknown distribution of the population (significance level: 5%). This is because the assumption of normally-distributed samples is questioned by the Shapiro-Wilk test. The null hypothesis is that the median of the population is constant before and after optimization, and a two-sided test is conducted.
The updated arrival flow can mitigate congestion in the enroute airspace. Figs. 18-21 summarize the reduction effect in the en-route airspace. In Figs. 18 and 19, mean/maximum delay times in the original/updated en-route airspace are compared, whereas Figs. 20 and 21 are these ratios with the increase and decrease in orange and green bars respectively. Table 3 summarizes the mean and maximum delay time values for the day i in improved en-route airspace, v U ,E,i,j , and their ratio to those in original en-route airspace with pvalue. We assume that the service time distribution g O,E (x) and abandon time distribution f O,E (x) remain unchanged, even after the arrival flow is optimized. Moreover, we set s U ,E,i (t) = s O,E,i (t) as well. Note that those three assumptions do not technically hold since the arrival adjustment affects the behavior of each flight inside the airspace, but it cannot be estimated in advance of the calculation in only the combined model.
From the findings, among 20 days other than day one with no delay in original and updated flows, the mean delay time over whole flights in a day is reduced in 19 days and the maximum delay time in the flights is improved in 18 days. Consequently, the average times over 21 days are decreased by 18.8% (5.0 s, mean) and by 16.5% (37.6 s, maximum), respectively and the amount of reduction is up to 18.0 s (mean in Day 12) and 181.6 s (maximum in Day 21). Among 21 days, three days show an unusual phenomenon. Day 13 can completely prevent en-route jam-up but there is almost no delay, even in the original situation. Day 6 is the exception: both the mean and maximum delay times of the flights deteriorate. Day 10 is another exception: mean delay time is reduced whereas the maximum delay time deteriorates. These adverse effects can be caused by the simple assumption of constant flight time or the setting of airspace capacity in mathematical programming. In other words, input flow as the same flight time for all flights in the programming problem would rather thicken congested flights of the stochastically-distributed flight time in the fluid model. Finally, the Wilcoxon signed-rank test can demonstrate that the optimization operation significantly reduces the median of mean and maximum time distribution of population, concluding that it can be effective to mitigate arrival flow congestion.    is on w U ,E,7 (t) [×10 s], v U ,E,7 (t) [×10 s], and X U ,E,7 (t) = B U ,E,7 (t) + Q U ,E,7 (t). Comparing it with Fig. 10 and Fig. 11, it is clear that delays in the first half of the 5 hours are completely reduced.

3) DELAY REDUCTION IN TERMINAL AIRSPACE
Figs. 24-27 is the summary of the reduction effect in the terminal airspace. Table 4    The findings indicate that among the 21 days, both the mean delay time over whole flights in a day and the maximum delay time in the flights are improved over 15 days. Consequently, the average times over 21 days decline by 6.8% (mean) and by 10.0% (maximum), respectively. Concerning the amount of reduction, Day 6 reduces both the mean delay time of the targeted flights and the maximum delay time of a flight by 29.6 s and 180.6 s, respectively. However, we cannot conclude that the median of the population of mean values declines significantly due to interval control, whereas it is       (see Figs. 12 and 13), the absolute value of delay is reduced. With the smoothing adjustment in the en-route airspace, the arrival rate in the terminal airspace can also be flattened,  concluding that controlling λ E (t) to realize less deviated flow can ease crowding in the terminal airspace can be proven.

V. VALIDATION OF FLIGHT TIME REDUCTION BY FAST-TIME SIMULATION A. BACKGROUND OF SIMULATION
In previous sections, our combined model with the G t /GI /s t + GI tandem fluid model and the nonlinear integer programming model demonstrates that early control of arrival intervals can reduce delays in both terminal and en-route airspaces: 18.8% (5.0 s) of mean delay time and 16.5% (37.6 s) of maximum delay time for an average for 21 days and mean delay time of up to 18.0 s and maximum delay time of 181.6 s in en-route airspace. With regard to that in the terminal airspace, 10.0% (32.9 s) of maximum delay time is reduced for an average of 21 days and up to 29.6 s of mean delay time and 180.6 s of maximum delay time are cut respectively. However, the result should be validated in terms of the quantitative amount of delay reduction or its qualitative tendency because the model is so macroscopic that flight-byflight characteristics are not taken into consideration: actual operation by air traffic controllers is not explicitly regarded. Moreover, the integer programming formulation is based on a simple assumption of constant flight time for all flights, despite the day/time/path-dependent flight time in the real world. Thus, in this section, we directly calculate the flight time of arrivals by a trajectory-based simulation with rules on ATC to demonstrate that it is improved by controlling the inter-arrival times at the borderline of the en-route airspace as well. Flight times are assumed to be the sum of the fastest passage time possible and the reducible delay time by the control, whereupon the gap between the flight time before and after optimization can be linked to the reduced delay time in our model.

B. SIMULATION ENVIRONMENT
A multiagent-based simulator (AirTOP software) is applied to validate the results obtained from the queue-based model [29]. AirTOp can seamlessly model and simulate air traffic flow in the airport and airspace, which enables users to assess and improve airport and airspace capacity and complexity. Compared to the queue-based model, AirTOP explicitly takes into account wind speed and the characteristics of airspeed according to the aircraft type and state (e.g. climb, cruise and descent), incorporating the Base of Aircraft DAta (BADA) model [30]. Although rule-setting for high-fidelity traffic flow is challenging, AirTOP has the benefit of outcomes close to the actual operation, including 4D-trajectories, en-route and terminal delay time and fuel consumption.

1) TRAFFIC SCENARIO
The traffic scenario is created from the flight plan (FP) and radar track (RD). An FP is a plan that notifies the air traffic control agency when an aircraft flies. An FP includes the call sign, aircraft type, the name of the captain, flight rules, departure fix, departure time, cruise altitude, route, cruise speed, radio equipment, alternative airport, fuel load, the total number of passengers, and other information. From the FP, the call sign, cruise altitude and routing structure are extracted. RD is the time series trajectory data created from the actual data and the flight plan data from the radar information processing system for air traffic control. RD includes the time recorded by approximately ten seconds, the call sign, the latitude, longitude, altitude and aircraft type. From RD, the aircraft type, geographical point of departure including latitude and longitude and departure time are extracted. Fig. 30 demonstrates the airways to entry fixes of the Tokyo Approach Control Area (TACA), which is the terminal airspace of RJTT, in a day. The purple line is the boundary between the en-route airspace and Tokyo Approach Control Area (TACA, or terminal airspace in this paper). The blue dot indicates XAC whereas the blue line represents a boundary 100 NM away from XAC. Red lines are airways to entry fixes (black dots) of TACA. A geographical point of departure is defined as the passing point 100 NM away from XAC corresponding to the analysis area in the queue-based model. The time when each aircraft passed the departure geographical point corresponds to the departure time. After the departure, the aircraft fly to RJTT according to the routing structure defined in FP.
To validate the result obtained by optimizing inter-arrival times at the boundary of the en-route airspace, simulations are executed by assigning different numerical values to departure time. As described in Section IV, b j ∈ {0, 1, . . . , T } is the original departure time whereas a j ∈ {0, 1, . . . , T } is the departure time after the optimization. Accordingly, b j and a j are directly assigned to the departure time for this simulation.

2) AIRSPACE AND ROUTES
Information on sectors, waypoints, and runways is obtained from the Aeronautical Information Publication [31] issued by the Japan Civil Aviation Bureau. The en-route airspace in Fukuoka FIR and Tokyo Approach Control Area (TACA) is modeled in the simulation. The boundary of TACA is shown as the purple line in Fig. 31. In TACA, Standard Terminal Arrival Route (STAR) to 34L at RJTT is configured following AIP, with speed and altitude constraints at the designated waypoints. Although TACA has six entry fixes, STAR from SPENS and SELNO are modeled, since the analytical area in the queue-based model is airspace in the southwest direction to RJTT.

3) SEPARATION RULES
The separation rules due to the delay time inside and outside Tokyo Approach Control Area (TACA) are defined in the simulation.
Inside TACA, speed adjustments on STAR are performed to keep 5 NM in-trail separation. Just before the final approach, vectoring is instructed in the red polygon in Fig. 31. The separation value is based on the wake-turbulence category of current ICAO standards. In case the delay cannot be completely absorbed in the vectoring area, the holding operation is conducted at the blue dots in Fig. 31 just before vectoring area. At the blue dot, the aircraft can turn right to make a circle, with a maximum length defined as 4 minutes.
Outside TACA, the in-trail separation rule due to implicit flow management is adopted. Before entering TACA, the aircraft is instructed to engage in path-stretching to ensure 10 NM and 12 NM inter-aircraft spacing at SPENS and SELNO, respectively. The path-stretching direction is right following actual operation (as shown in Figs. 32 and 33) The maximum deviation angle is 80 degrees from the designated airways in Fig. 30, which allows aircraft to conduct radar vectoring as much as possible.

4) WEATHER CONDITION
This study assumes weather conditions that do not affect air traffic operations such as runway-change operations. Under the above conditions, the wind direction and speed affect the ground speed of the arriving aircraft. For meteorological conditions, the Mesoscopic Numerical Prediction Model GPV (MSMGPV) provided by the Japan Meteorological Agency is utilized in the simulation. This data includes atmospheric characteristics such as wind and temperature on a 3D grid covering all of Japan and is published every three hours in the database of the Research Institute for Sustainable Humanosphere, Kyoto University [32]. GPVs are located at  the time can be significantly reduced by the operation (the maximum flight time is reduced by 6.7% at most (Day 8) for ten days and constant for nine days and 1.0% (6.5 s) as an average for 21 samples and the maximum reduction effect is 51.0 s (Day 8)).

A. COMPARISON WITH PREVIOUS STUDIES
In this paper, we applied the time-varying queuing model and showed the delay reduction effect in the terminal area and surrounding upstream area. The qualitative result is consistent with the previous study which adopted the time-constant queuing model in a similar setting [13]. It introduced the G/G/c queuing model with a general arrival/service process and a constant number of servers in the terminal airspace of Tokyo International Airport. By data-driven analysis, they suggested that reducing the variance in inter-arrival times at 50-60 NM away from the airport could limit the arrival delay in the area. Our combined model with the AirTOP simulation not only validates their findings qualitatively but also enables the time-varying behavior of each flight to be analyzed. The G/G/c queuing model cannot track the time/flight-dependent behavior, making it impossible to determine which flight experiences heavy congestion and when it happens. Note that the G/G/c queuing model analytically develops the expression for the expected delay time, which is determined by the number of servers, the mean and variance in inter-arrival times, and those of service times (see [12]).
Moreover, as cited in [15] and [10] can be one of the benchmark studies, which introduced a two-stage (before and after entering the terminal airspace) stochastic programming approach, considering the uncertainty of entry time at the Initial Approach Fix (IAF) of Paris Charles-De-Gaulle airport. They estimated that the maximum time to lose within the terminal airspace was cut down by 73.5% under simulation-based experiments for ten arrivals (seven flights of wake-turbulence category H and three flights with M) from three directions. Decision variables include the continuous ones for the target time at the IAF, the binary ones for the aircraft sequence over the IAF and the continuous ones for the target landing time and it is assumed that the condition of deviations in terminal-entry times from target times follows N (0, σ 2 ) (σ =30 s). Their approach involves minimizing the sum of the length of the landing sequence, flight time before entering the area and the difference between unconstrained and target landing times respectively. Although they outperform us, the computation time for ten flights arriving within a 20-minute time window is 59.1 s (100 scenarios conducted in the second stage), much larger than 2 s of our integer programming model for about 135-151 flights. Note that our fluid model can calculate functions for the same input flow arriving for five hours within two minutes. The results of our combined model and the previous model [10] are summarized in Table 7.  [10] about reduction rate of maximum delay time, number of arrivals and its span, and computation time. Both approaches have drawbacks and advantages: our model realizes fast calculation (two minutes at most) and large-scale modeling (135-151 flights for five hours) but needs to improve the strategy for controlling intervals such as performing the operation even at the terminal airspace boundary as well. Conversely, the model in [10] achieves a higher flight-time reduction effect on a moderate calculation but has flaws in the time horizon and number of aircraft targeted. An important question is to verify whether the model retains these advantages in the same scale of arrival traffic flow as ours, or to apply our combined model to the data set [10] to estimate its performance for more direct comparison as one of the extensions of this paper. Table 8 is the summary in the Sections III-B, IV-B and V-C. The figures show the reduction effect of each time for an average of 21 samples and those in bold are to be shown as significantly reduced from the result of the Wilcoxon signed-rank test for the samples. Although the tendency is qualitatively demonstrated to ease crowding by the interval control at the entry of the en-route airspace, not only in the area but also in the terminal area, here we concentrate on the quantitative aspect and discuss the difference between the value in the terminal airspace, where the value differs quite significantly.

B. COMPARISON BETWEEN COMBINED MODEL AND SIMULATION
Since the flight time in the AirTOP Simulation includes the vectoring delay time to ensure the separation at the Initial Approach Fix (SPENS or SELNO) and the point-merge vectoring and holding delay time to adjust intervals for the final approach, the time can be broken down into two components: Flight Time (Simulation) = Time without Vectoring/Holding +Time with Vectoring/Holding.
According to the definition of the fluid model (see Fig. 1), on the other hand, the time in the model can be resolved into  We presume that the service time corresponds to the first term in (10), and the delay time to the second term. Accordingly, the amount of delay reduction is the difference between the second term of (10) (simulation) or (11) (model) before and after optimization if the first term (unconstrained flight time) remains unchanged: the first term corresponds to the unconstrained flight time with no operation. From this discussion, we would like to calculate delay time in the simulation, namely, the flight time with vectoring or holding in (10).
Only the terminal airspace is targeted in this paper and the procedure for calculating the delay time is the same as that for the flight time. Figs. 34-37 summarize the mean and maximum delay times in the AirTOP simulation. The black bar in Day 2 in Fig. 37 represents no change between the operation, namely, a ratio of 100%. Table 9 shows the mean and maximum delay times for 21 days in the terminal airspace. The median of the mean and maximum delay time distribution of population is significantly decreased: average mean delay time by 11.5% (36.5 s) and average maximum delay time by 19.2% (148.8 s), respectively. Note that the assumption of fixed unconstrained time is reasonable, since it is nearly unchanged from the result VOLUME 11, 2023  It is concluded again that crowding can be alleviated by controlling the inter-arrival times at area upstream of the enroute airspace, but the analysis of delay time by AirTOP simulation reveals the limitation of our current combined model. Table 10 shows a direct comparison of estimated and calculated delay times in the terminal airspace. It is derived that our combined model cannot be accurate quantitatively although it can determine the qualitative tendency and the AirTOP itself cannot be a exact simulation of arrival air traffic flow. The main reason why the amount takes different values is that we assumed that service time distribution equates to flight-time distribution in the model, which technically contradicts (11) above. Concretely, service time distribution is assumed to g O,T (x) = N (1187.0, 176.6 2 ), but the actual and original mean service times are estimated as 915.6 s (subtracting 316.2 s in Table 9    and calculating delay time in the fluid model, which will be described in the next section.

C. POTENTIAL FUTURE WORKS
Several tasks remain outstanding, as reflected in the discussion of the previous section. Firstly, the delay time in the en-route airspace should be estimated in the AirTOP simulation as well as that in the terminal airspace in Section VI-B.
Secondly, service time distribution should be refined. Excluding the possible modification mentioned above, for the arrival flow of each day i, the service time distribution g(x) can be formulated from the arrivals only on that day, rather than all arrivals for 21 days. Although the accuracy in the current model is guaranteed by (2), the update can improve the accuracy of the model more because date-dependent arrival histories λ E,i (t) match the behavior inside the airspace each day, defined not by g(x) but by g i (x).
Thirdly, the iteration in the combined model referring to the result of the AirTOP simulation can be effective. As mentioned above, identifying service time distribution with flight-time distribution in the fluid model is one of the reasons explaining the deviation between the reduction effect of the model and the simulation, even though (2) ensures accuracy to some extent. One possible procedure is as follows: after estimating service and delay times under conditions where both distributions are identical and that of flight time in the simulation thereafter, the mean (or standard deviation) of service time distribution in the model is adjusted by the result of flight/delay time estimation in the AirTOP simulation and then re-calculating the functions in the model. It can lead to a more appropriate estimation of time-varying functions including delay times.
Regarding the extension toward a microscopic approach, the queuing network can be elaborated by trajectory-based segmentation. Currently, the queuing network is distancebased segmentation: two tandem queues are set according to the distance from Tokyo International Airport and the benchmark waypoint XAC in the terminal and en-route airspace, respectively. This segmentation is based on previous studies which indicated that inter-arrival and service times follow mathematically tractable distributions [12] by distance-based segmentation. However, given that the actual arrival air traffic flow has several branches (see Figs. 5 and 6) and the amount of flow or the flight status can vary with routes, it may be more appropriate to divide the area by routes to build a parallel network.
Additionally, we can incorporate into our study insights from the analysis and control of bursty traffic in a computer network; the validity of presuming not the Poisson arrival process but the general arrival process can be supported by establishing whether the arrival traffic flow has self-similarity. To judge whether the arrival process of the flow has a long-range dependence regardless of timescales [33], we can apply the autocorrelation function-based approach [34] or the phase transition-based approach [35]. Not only is the assumption of the general arrival process verified, but also the technique for controlling the self-similar traffic in the data network can be applied to the arrival traffic flow.
Finally, to enhance the utility of the combined model for the implementation for various arrival traffic, evaluating computational performances should be required. Concretely, a computation load in a larger-scale arrival flow with a larger number of decision variables should be determined to discuss the scalability or applicability of the optimization problem. As [36], where a proposed model was tested by a larger-scale benchmark instance set to solve in a reasonable time to demonstrate superiority over conventional approaches, the combined model should be employed for day-long arrivals at Tokyo International Airport or more congested flow in other airports. Comparing performances with those in this paper can provide insights into the issue.

VII. CONCLUSION
This paper demonstrated the scope to prevent delay propagation from the en-route phase to the terminal phase by controlling arrival intervals to minimize the variance in inter-arrival times, using the tandem-queuing network model. For 3,074 arrivals of 21 days flying the en-route and terminal airspace around Tokyo International Airport during the most congested period between 17:00 and 22:00, the G t /GI /s t + GI queuing model and the nonlinear integer programming model were applied: the former for analyzing the time-varying char-acteristics of the current and updated airspace, including delays and the latter for the optimization problem, which explicitly considers the variance in inter-arrival times under the capacity constraints. Analysis of current airspace showed that delays are mainly attributable to relatively congested arrivals over several periods rather than a momentary rush and that delays in the en-route airspace are propagated or amplified in the terminal airspace. Moreover, we found that flight control at upstream airspace to sequence flights at the boundary of en-route airspace can mitigate congestion, not only in the en-route airspace but also the downstream terminal area. From the Wilcoxon signed-rank test for 21 pairs of representative values, delay times are significantly reduced by 18.8% (5.0 s) of mean delay times and 16.5% (37.6 s) of maximum ones in the en-route airspace and by 10.0% (32.9 s) of maximum ones in the terminal airspace, which should be improved in terms of accuracy. Furthermore, to complement the operational constraints around the airport, a fasttime simulation by AirTOP is conducted to qualitatively and quantitatively validate the combined model, concluding that arrival control from an early stage significantly reduces the mean flight times in the en-route airspace by 1.4% (7.2 s) and the mean and maximum delay times in the terminal airspace by 11.5% (36.5 s) and 19.2% (148.8 s).