Data and queueing analysis of a Japanese air-traffic flow

European Journal of Operational Research 235 (2014) 265–275

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Innovative Applications of O.R.

Data and queueing analysis of a Japanese air-traffic flow Claus Gwiggner a,⇑, Sakae Nagaoka b a b

Institute for Operations Research, University of Hamburg, Germany Electronic Navigation Research Institute, Tokyo, Japan

a r t i c l e

i n f o

Article history: Received 14 April 2012 Accepted 21 October 2013 Available online 1 November 2013 Keywords: Queueing Applied probability Air traffic management

a b s t r a c t Congestion is a major cause of inefficiency in air transportation. A question is whether delays during the arrival phase of a flight can be absorbed more fuel-efficiently in the future. In this context, we analyze Japan’s flow strategy empirically and use queueing techniques in order to gain insight into the generation of the observed delays. Based on this, we derive a rule to balance congestion delays more efficiently between ground and en-route. Whether fuel efficiency can be further improved or not will depend on the willingness to review the concept of runway pressure. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Many airports and airspaces are already congested today, but growing traffic demand will cause increasing delays in the future if no actions are taken (ICAO, 2005; JPDO, 2010; Sesar-Consortium, 2007). Congestion delays materialize either on the ground, where aircraft have to wait before accessing a runway or during the flight, where they are deviated from their intended trajectory. When the schedule buffers are tight, delays may also propagate through the whole transportation network. Congestion delays can be managed on a strategical level (by runway expansion or shorter separation standards), a pre-tactical level (by splitting flows and sectors) or a tactical level (by sequencing and re-sequencing aircraft). Fig. 1 shows the typical situation. Prior to take-off, departure slots are allocated in order to balance demand with available capacity. Before landing, flows are merged and sequenced in order to efficiently use the runway. In this context, one speaks of radar vectoring, when traffic controllers ask the pilots to stretch their flight paths in the form of specific flight angles and of holding stacks, when pilots are asked to circulate in the vicinity of the airport. Departure slot allocation and en-route sequencing operations create delays, and we call them metering delays in this paper. The Japanese flow managers currently allocate up to 10 min of metering delay to the en-route phase and the remaining delay to the ground. One can argue that en-route delays are unnecessary, once ground delays are correctly allocated. At least, speed control shall synchronize traffic flows so that radar vectoring can be reduced. But due to arrival-time uncertainties, en-route delays ⇑ Corresponding author. Tel.: +49 40 42838 6412. E-mail address: [email protected] (C. Gwiggner). 0377-2217/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ejor.2013.10.056

serve as buffers in order to keep a high pressure on the runways. On the other hand, future reduction of arrival time errors promises to reduce en-route delays, allocating more and more delay on the ground. The benefits are smoother traffic flows, less controller workload and reduced fuel consumption. In this paper we derive a decision rule to quantify this new balance. The deterministic ground delay and aircraft sequencing problems are well known in the OR and the transportation communities, see for example Erzberger (1995), Beasley et al. (2000), Bayen et al. (2004), Lulli and Odoni (2007), Artiouchine et al. (2008) and Balakrishnan and Chandran (2010). The stochastic phenomenon of congestion was pioneered by Cox (1955), defining it as (i) a flow of customers needing service, (ii) some restrictions on the availability of service, and (iii) irregularity in the flow of customers, the servicing operation or both. Moreover, when average demand is higher than available capacity (e.g. during peak hours or bad weather conditions), predictable congestion occurs. Phenomena with the above characteristics, such as telecommunication networks, passenger flows, and many more, are often modeled as stochastic or deterministic queueing systems (Newell, 1982; Wolff, 1989). The air transportation system also exhibits the classical queuing behavior: as demand approaches capacity, delays increase sharply (Ball et al., 2010). Early contributions in this regard are by Bell (1949) and Dunlay and Horonjeff (1976). Often, arguments for Poissonian arrival flows are elaborated on which steady-state control strategies are based. A homogeneous Poisson process assumes essentially that the numbers of arrivals in disjoint time intervals are independent random

266

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

Fig. 1. Metering delays for major Japanese arrival flows (Data Source: Gwiggner et al. (2009)).

variables, and that no more than one arrival takes place per time instant (see Daley and Vere-Jones (2003) for an axiomatic derivation and historic derivations from Erlang and Bateman). Also since long time, the existence of steady-states in air transportation are questioned (e.g. Barnhart et al., 2003; Janic, 2000). One of the few transient analyzes is Peterson et al. (1995), based on an algorithmic approach. We will later show that the Poisson hypothesis is not unreasonable whenever congestion is moderate, and that steady-states can be reached for pre-scheduled flows in the order of a few hours. The assumptions of homogeneous Poisson flows have been relaxed in at least three ways. Non-homogeneous queueing models allow for variations of arrival and service rates. This seems more realistic to capture high traffic peaks, particularly during bad weather. Typically, time-dependent Poisson arrivals are modeled with time-dependent service distributions (Bäuerle et al., 2007; Fan and Odoni, 2002; Hansen et al., 2009; Knessl and Yang, 2002; Pyrgiotis et al., 2013; Yang and Knessl, 1997). Models with deterministic service times serve as lower bounds and models with exponentially distributed service times as upper bounds to the real delays (Barnhart et al., 2003). Multiple arrivals per time instant (batch arrival models) are another relaxation of the Poisson assumption. Corresponding queueing models are proposed for aircraft landings (Falin, 2010; Krishnamoorthy et al., 2009). More recently, a third departure from Poisson arrivals was presented in Guadagni et al. (2011) and Guadagni et al. (2013). They clarify that the numbers of arrivals in subsequent time intervals will be negatively correlated, whenever aircraft are pre-scheduled at a critical resource but arrive with random errors. Pre-scheduled queues are also studied in the field of ‘appointment scheduling’, where customers receive initial assignments, but may arrive with a positive or negative delay. From the queueing methodology point of view, there was some activity in the 1960’s and 1970’s, but it loses its traces then (Bloomfield and Cox, 1972; Mercer, 1960, 1973; Winsten, 1959). Today, the field is widely occupied with algorithmic questions (e.g. Begen and Queyranne (2011), Green (2008) and the references therein) or numerical approaches (Nikoleris and Hansen, 2012). The remainder of the article is organized as follows: in the next section we describe a radar-data analysis of arrival traffic for a Japanese airport. We then analyze two types of queueing models as generators for metering delays: pre-scheduled random arrivals

(PSRA) and Poissonian arrivals. Based on this, we derive a rule how to balance metering delays between ground and air efficiently.

2. Data analysis We analyzed available radar data in order to get a first impression on the statistical properties of the traffic flows and the generated delays. We selected nine days of ‘normal’ traffic from the months August, October and December 2008, i.e., where no exceptional events or delays were reported. We removed outliers by hand (about 10% of missing or erroneous fields in the source data). As main results we found empirical evidence for and against a Poisson hypothesis during moderate congestion and physical evidence against a Poisson hypothesis during very high congestion. We first verified the three-dimensional structure of the problem: aircraft enter the sector on different routes and different altitudes but leave it as a single flow on the same altitude. Average ground speed at the sector entry grew with the altitude, ranging between 484 kt and 522 kt (knot). At the sector exit, the average speed was equally about 380 kt with 27 kt standard deviation. The timeseparation sm at the metering point was on average 1.4 min, with standard deviation of 4 s. This corresponds exactly to the 10 NM (nautical mile; 1 NM = 1.852 km) rule based on the ground speed. No unexpected patterns were found on this level of analysis. Fig. 2 helps to better understand why delays occur. The top panel shows the flow (aircraft (ac)/5 min) at the metering point between 7:30 and 21:00. The red1 bars are from the flight plans, the green ones from the radar data. One can see a fluctuating demand, with slightly higher periods in the morning and evening hours. Note that these flight plans are already pre-scheduled to meet the available capacity. The dotted horizontal line is the daily average arrival flow k with approximately k  0:3 ac/min (it is almost the same for the planned and the realized flights). The bold horizontal line is the capacity at the metering point l ¼ v out =sm ¼ 0:63 ac/min, where v out is the average speed of all aircraft at the sector exit. While the flight plans sometimes exceed the capacity, the radar data generally lies below it. The middle panel shows the sample sector density 1 For interpretation of color in Figs. 1, 2 and 4–9, the reader is referred to the web version of this article.

267

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

5

Demand at West gate: flight plan and radar data

3 2 0

1

# ac / 5 min

4

flight plan radar capacity hourly rate daily rate

8

9

10

11

12

13

14

15

16

17

18

19

20

21

slot (5 min) 2008/10/20

50

10

Sector density and average delays density

30 20 10

total delay (min)

6 4

0

0

2

ρ (# ac)

8

total delay

40

capacity excess

8

9

10

11

12

13

14

15

16

17

18

19

20

21

slot (5 min) 2008/10/20

0.2 −0.2

0.0

ACF

0.4

0.6

Crosscorrelation Inflow−Delay

−20

−10

0

10

20

Lag Fig. 2. Top: Planned and observed flights. Middle: Sector density and total delay. Bottom: Sample cross-correlation.

qðt þ 1Þ ¼ qðtÞ þ qin ðtÞ  qout ðtÞ (number of aircraft in sector in time slot t þ 1) (green) and total delay (black). Here, qin ðtÞ and qout ðtÞ are the corresponding sample flows (ac/5 min). One can see oscillating sector load with a period of approximately 30 min throughout the day. At the end of each oscillation, the total delays take peak values. The red bars are the demands where the capacity at the metering point is exceeded. Accordingly, the sample crosscorrelation function between inflow and total delay has two peaks of about 0:6 and 0:4 at lags 3 and 4 respectively, corresponding to the average sector traversal time of 15–20 min (bottom panel). An important question is if demand fluctuations can be predicted in advance or if they occur spontaneously. For example, peak hours should lead on average to more demand than the rest

of the day. This is one of the reasons to study non-stationary models of air traffic in the U.S. (e.g. Hansen et al., 2009, Tu et al., 2008). On the other hand, Japanese flow management distributes departure slots on a permanent basis. This is similar to Europe, but differs from the practice in the U.S., where ground delay programs are only initiated during severe weather conditions. As Odoni et al. (2011) show, the delay patterns in the European case are more regular than in the U.S. In order to identify systematic changes in the arrival rate, we plotted sample paths of the counting process. Without any structural assumptions, we are just given an increasing sequence of positive random arrival times ðT n Þn2N , defined on the events x ¼ ðx1 ; x2 ; . . .Þ 2 X of an underlying probability space ðX; F ; PÞ.

268

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

X 1T n 6t ;

C¼n

Z

1

zn ðxÞ2 wðxÞdFðxÞ;

1

2 0 0

2

4

6

8

10

theoretical quantiles (min) Fig. 4. Quantile plot for inter-arrival distributions with exponential and PSRA distributions.

where wðxÞ is a weight which allows to vary the importance of different parts of the distribution. For example, wðxÞ ¼ 1 leads to the Cramer-von-Mises test statistics (CvM), whereas wðxÞ ¼ ðFðxÞ½1  FðxÞÞ1 leads to the Anderson Darling (A-D) statistics. For comparison, the Kolmogorov–Smirnov test computes the largest gap Dn ¼ supx jzn ðxÞj (Stephens, 1993). We simulated the distributions of the test statistics under the null hypothesis by Monte Carlo, taking into account the variability of the estimated parameter. These distributions were used to determine the p-value, that is, the probability of drawing a sample that is at least as large as the observed test statistics. In case that such a probability was less than a%, we considered the null hypothesis as rejected on a level of a. Table 1 shows the results, where the entries summarize the obtained minimum and maximum values from all tests (one for each day of data). In short, in two out of nine samples, the hypothesis of exponentially distributed inter-arrival times was rejected on a 5% level (both tests rejected the same samples). These two samples are plotted in green in Fig. 4. No pattern can be detected, even after re-estimation of the sample quantiles with interpolation methods. We also tested the Poisson hypothesis of numbers of arrivals with v2 goodness-of-fit tests. The tests were performed for time intervals of sizes 5, 10 and 15 min. Our samples were always larger than 200, so inference was based on the asymptotic distribution of the test statistics. The results are shown in the lower part of Table 1. In the 5 min interval, the null hypothesis was rejected in five out of nine tests. For the 10 and 15 min intervals, the p-values were all higher than 5%. An inspection of the results in the 5 min intervals showed that high peaks caused the deviations with the Poisson distribution (see Fig. A.11 and Table A.3 in Appendix A).

100

150

Exponential Rejected PSRA High density

50

Table 1 Goodness of fit tests. Reported are minimum and maximum values among 9 tests.

0

Counts

200

250

T09: counting process

observed quantiles (min)

where 1E is the indicator function for the event E. A sample path of NðtÞ is a realization of N x for a particular outcome x. It has the shape of a step function, increasing by one each time an arrival occurs. This is shown in Fig. 3 for the nine days of data. The plotted lines are reasonably straight, showing little indication of a change in frequency of arrivals over the day. Ruling out systematic changes in the arrival rate, the question is to which extent the flow can be seen as a homogeneous Poisson process. On the one hand, aircraft enter the sector on six different routes and various altitudes. On the other hand, they are prescheduled and obey minimum separations. Fig. 4 shows quantile plots for inter-arrival distributions of the nine days of data. A quantile plot compares the distributions of two random variables graphically. The red lines are the observed inter-arrival times against an exponential distribution. They are close to a straight line and fluctuate slightly in the tail, which is not uncommon in these plots. The blue and gray lines arise from other models and will be explained below. Under the assumption that successive inter-arrival times are independent and identically distributed, we performed goodnessof-fit tests. The null hypothesis was that the observations are drawn from an exponential distribution with unknown parameter k, the alternative hypothesis was that they are not (see below for a comment). We used the family of Cramer-von-Mises tests because a parameter of our null hypothesis (the arrival rate) had to be estimated. Given a sample x1 ; . . . ; xn of a random variable X, these tests compute the difference zn ðxÞ ¼ F n ðxÞ  FðxÞ between the empirical distribution function of the sample F n ðxÞ, that is, the fraction of observations that are less than x, and the distribution function FðxÞ of X, that is, the probability that an observation is less than x. More formally, the test statistics of the Cramer-von-Mises family are

6

8

ð1Þ

nP1

4

NðtÞ ¼

10

A counting process N x ¼ ðNðtÞÞt2R is a family of (dependent) random variables

8

9

10 11 12 13 14 15 16 17 18 19 20 21

Time (h) Fig. 3. Counting processes.

Test

Statistics

p-values

#reject (5%)

A-D CvM

[0.24, 1.5] [0.03, 0.32]

[0.03, 0.92] [0.02, 0.94]

2 2

Test

Interval

Statistics

p-values

d.f.

#reject (5%)

2

5 min 10 min 15 min

[1.7–19.5] [ 3.7–10.8] [4.6–10.1]

[0.001, 0.64] [0.06, 0.6] [0.2, 0.7]

2 [4, 6] [7, 8]

5 0 0

v

269

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

By the logic of hypothesis testing, goodness-of-fit tests cannot lead to affirmations of the null distribution. Moreover, the power of these tests, i.e. the probability to chose the alternative hypothesis when it is true, is difficult to assess. Typically, Monte-Carlo simulations with a number of ‘typical’ alternative distributions are evaluated (Stephens, 1974). We thus conclude that there is empirical evidence for and against a homogeneous Poisson process as a description for our data. This makes sense for two reasons: firstly, aircraft enter our sector on six different routes and various altitudes. Moreover, there is only moderate congestion because we analyzed a gate in an en-route sector. Both suggest a high degree of randomness in the arrival instants. And secondly, it is known that pre-scheduled arrival flows converge in distribution against the Poisson process, when the uncertainties are very high (Guadagni et al., 2011). But under heavy congestion, such as at the runway threshold, traffic flows are highly regular, for which the Poisson model is no accurate description. The gray curves in Fig. 4 illustrate such a case. In any case, a Poisson process does not take into account pre-scheduled trajectories and their prediction errors, because it depends only on one parameter; the arrival rate. In the next section we analyze a model that reflects the underlying preschedule. Its inter-arrival distribution competes well with the observed data for both, moderate and heavy congestion. These are the blue lines in Fig. 4. 3. Queueing analysis The empirical analysis confirmed that demand fluctuations are the major driver for metering delays, so we modeled the metering gate as a single server queueing system. Input to the system are flows from the different routes. Output is a single flow, separated by at least sm = 10 NM. Depending on speed, wind conditions and human factors, it may vary from time to time. We expected from the analysis to further clarify which part of the observed delays is due to metering and which to operational uncertainties. We analyzed this situation by two families of queuing models; pre-scheduled random arrival (PSRA) and homogeneous Poisson arrivals with either deterministic or general service (M/D/1, M/G/1 in Kendall’s notation, meaning that arrival processes are Poissonian, service times are deterministic in the first and generally i.i.d. distributed in the second, each having one server (Kendall, 1953)). Note that in the same paper, Kendall pointed out that pre-scheduled, but randomly delayed arrivals do not exhibit independent inter-arrival times, as it is most often assumed in the GI arrival process models. As mentioned in the introduction, many other authors analyzed M/G/1 models and generalization for air traffic flow prior to us. Our contribution is that we put this model in context with the more realistic model of pre-scheduled arrivals. 3.1. Pre-scheduled random arrivals A recent model of air traffic flow is given by the so called prescheduled random arrival process (PSRA) (Guadagni et al., 2011). Aircraft are scheduled to arrive at fixed (deterministic) times, but uncertainties, such as wind conditions, lead to arrival-time errors

atai ¼ stai þ ni ;

allows an analytical treatment. The i.i.d. assumption seems strong, for example external factors such as wind or airport congestion might affect large subsets of flights in a predictable way. However, arrival time errors also depend on individual factors, such as aircraft type, weight, and conflict avoidance manoevres. This uncertainty remains at the core of the problem, and assuming independence here seems reasonable (although the distributions might not be exactly identical). We will analyze Guadagni et al.’s model in more detail and generalize it to situations where traffic is pre-scheduled on a permanent basis, as it is the case in Japan and Europe. Fig. 5 shows the PSRA process on the time axis. The points illustrate pre-scheduled arrivals. The arrows indicate the possible arrival-time errors. The number of arrivals #fg in time slot ½t; t þ 1Þ is

mt ¼ #fi : t 6 atai < t þ 1g:

ð3Þ

This arrival process is serially correlated, which can be seen as follows. We allow the three events ‘deletion’, ‘early arrival’ and ‘late arrival’ for every pre-scheduled arrival: let p be the probability of non-deletion and q ¼ Pðata < staÞ. Then the event ‘early arrival’ occurs with probability pq and ‘late arrival’ with probability pð1  qÞ. Let us further assume that these operations occur independently for every aircraft. Based on this we can obtain the conditional arrival distribution

mij ¼ Pðmtþ1 ¼ jjmt ¼ iÞ:

ð4Þ

For illustration let us further assume Pðj ni j> 1Þ ¼ 0, meaning that arrival time errors are almost surely shorter than the minimum separation (the service time), which we normalize to 1. In this case, slots ½t; t þ 1Þ and ½t þ 1; t þ 2Þ are only affected by the realization of the events of three aircraft a; b; c in Fig. 5. We can thus determine the joint probability of all 33 possible combinations. For example, deletion of all three aircraft has probability ð1  pÞ3 and is one of the 4 events contributing to Pðmtþ1 ¼ 0jmt ¼ 0Þ. In this way we obtain the exact distribution for mij . When the jump lengths are larger than 1, the situation gets more complicated. Winsten (1959) already had a similar approach but did not consider the transient properties. Also, methods from particle physics or batch processes might be helpful, as illustrated in Guadagni et al. (2011) and Lucantoni (1993). Intuitively, the more pre-scheduled arrivals miss their time slot, the more aircraft will arrive in the neighboring slots, leading to a negative autocorrelation. This allows us to describe a discrete time queuing system with PSRA input as a bi-variate homogeneous Markov chain. It depends on the queue length nt and on the state of the arrival process mt

ðn; mÞtþ1 ¼



ðnt þ mtþ1  1; mtþ1 Þ if nt > 0 ðmtþ1 ; mtþ1 Þ

else;

where nt is the length of the queue immediately before service at time t and mt is the number of arrivals in time interval ½t; t þ 1Þ, as defined previously. Our queueing system is similar to M/D/1 systems. The main difference is that mt is auto-correlated, as mentioned above. Note also that the states ðn; mÞ, with n < m are invalid, for example ð0; 2Þ can never occur. These invalid states are removed from the state space. The one-step transition matrix P ijkl ¼ Pððn; mÞtþ1 ¼ ðk; lÞjðn; mÞt ¼ ði; jÞÞ is then

ð2Þ

where stai is the pre-scheduled arrival time and the ni ’s are random variables with variance r2i . Guadagni et al. (2011) and Guadagni et al. (2013) consider the case where the initial schedule is regularly spaced, i.e. stai ¼ ki (k is the arrival rate) and generalize it by randomly deleting scheduled points. Moreover, they assume that the ni ’s are i.i.d. random variables with compact support (having probability 0 on the tails) or exponential distribution. Their model is motivated by strategic cancellations, such as in the U.S., and

ð5Þ

Fig. 5. PSRA and its counting process.

270

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

ðrÞ

lim pijkl ¼ pkl ;

r!1

8i; j; k; l:

ð6Þ

This allows us to compute the transient probability distribution of the queue length ðrÞ

k

ð0Þ r

¼k P ;

ð7Þ

ð0Þ

where k ¼ ðpij Þ is the initial probability distribution of the queue length coupled with the arrival process. Finally, taking expectations leads to the marginal queue length

E½ðn; :Þr  ¼

K X L X ðrÞ nknl ;

ð8Þ

n¼0 l¼0

where K is the total number of aircraft in one day (this is why the chain is finite) and L is the largest number of arrivals in one slot. Before we proceed, a generalization of the PSRA model is necessary. In Eq. (2) with stai ¼ ki , the pre-scheduled flow has a regular separation between two successive aircraft, which is always 1=k. Guadagni et al. (2011) add the possibility of flight cancellations by random deletions. But in Japan and Europe, pre-scheduled flows consist of flights with and without ground delay. In case of ground delay, their separation will equal the minimum separation sm . Otherwise, their flight plans will remain unchanged. So, the prescheduled separation between two consecutive flights is

 si ¼

sm

if wi > 0;

sm þ ai

if wi ¼ 0;

ð9Þ

where sm is the minimum separation, wi the initial queueing delay of aircraft i and ai þ sm is the initial separation according to the flight plans. These three flows—regular pre-schedule, random deletion and pre-schedule with ground delay—differ in their variance. A regular pre-scheduled flow has variance 0. In the case of random deletions, ai takes values sm ; 2sm ; . . . ; ism . Eq. (9) generalizes the separations for arbitrary pre-schedules. In the remainder we will speak of pre-variance in this context. Since queueing delays depend on the

16 14

The column ðn; nÞ of the first row contains m0n . The columns of rows ðn; iÞ with n > 0 contain the matrix given by Eq. (4) and shifted by Eq. (5). All other elements are 0. The matrix is finite, because at the end of the day, the airspace is closed. So the last rows are truncated. It is also irreducible which can be seen as follows. Given any two valid states ðn; iÞ; ðm; jÞ, it is always possible to move first to the state ð0; 0Þ by a sequence of zero new arrivals. From the empty queue, the states ðm; mÞ can be reached in one step. To reach all other states ðm; jÞ; m > j, a sequence of two new arrivals in every time interval leads to ð0; 0Þ ! ð2; 2Þ ! ð3; 2Þ !    ! ðm  j þ 1; 2Þ. From there, the state ðm; jÞ can be reached in one step by j new arrivals. Thus, there is a unique stationary probability pij for all states ði; jÞ. They form a row vector. Moreover, the states ðn; 1Þ; n > 0 have non-zero probability to remain in themselves because these transitions represent the event that one new arrival occurs subsequently. Thus, since the matrix is finite, all states are ðrÞ aperiodic. So the r-step transition probabilities pijkl ¼ Pððn; mÞr ¼ ðk; lÞjðn; mÞ0 ¼ ði; jÞÞ converge to their limiting stationary distribution

t (min)

18

20

Flown distance vs time

nominal h
150

160

170

180

190

200

dist (NM) W 2 , n: 894 , m: 41 Fig. 6. Distance vs. time for route W2.

variance of the processes, we expect that the PSRA model with regular separation under-estimates queueing delays in our application. In summary, although a queueing system with PSRA arrivals has similarity to M/D/1 queues, its properties are structurally different from that of a Poissonian system. 3.2. Parameter estimation ~ were determined as the difference Empirical metering delays w between the observed sector crossing times and estimations of nominal sector crossing times per route. Fig. 6 shows an example of one of the six routes. Observed crossing times are plotted against the flown distances inside sector T09. One can see a cluster of points in the range 160–170 NM, taking between 13.5 and 16 min, and a branch with correlated distance and crossing times. The red points mark the trajectories, where no radar vectors occurred. These trajectories were considered to be nominal. They took place typically in the early morning hours, where no congestion occurred. We verified the selection graphically, by plotting the corresponding lateral and vertical trajectories. The green points mark trajectories on a flight level higher than FL 350. No pattern can be detected, concerning altitude. For all six routes, the nominal crossing times lay between 14 and 15 min with standard deviations between 0.6 and 0.9 min. As an estimate, we used their sample average per route.2 This complicates the analysis because it obscures the true metering delays. As an ad hoc idea, we assume that the empirical ~ i consist of two components; the queueing delay metering delays w wi and an unknown measurement error i

~ i ¼ wi þ i ; w

ð10Þ

where i is a random variable, and i ; j are independent and identically distributed for i – j. In the remainder, we will refer to these ‘obscured’ queueing models as PSRA=D þ =1; M=D þ =1 and M=G þ =1. 2 We compared two data sources: flight plan data and radar data. It turned out that with the flight plan data, a large number of flights had negative delays suggesting that the flight plan information is inaccurate to estimate nominal sector crossing time.

271

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

Model Simulation Radar data

8 6

Correlated arrivals

2

4

ρ= 0.975 ρ= 0.95 ρ= 0.88

0

Average Queue length

Poisson arrivals

0

50

100

150

200

250

300

Time (min) Fig. 7. Transient behavior of average queue length.

We fitted the various queueing models on our instance. The major questions were (i) the transient behavior of the PSRA model, (ii) the difference between PSRA and Poissonian models, (iii) the difference between deterministic and general service times and (iii) delay patterns that cannot be explained by queueing models.

4.1. Transient analysis of the PSRA model

0.7

In order to validate the transient analysis in Section 3.1, we extracted radar data of a runway configuration change. This is an event that occurs when the wind direction changes. A number of aircraft has then to be rerouted, which causes increased queue lengths for a certain time. The question is how long it takes until the queue length returns to equilibrium. Fig. 7 shows typical results. The horizontal axis is the elapsed time, starting at the beginning of the runway configuration change. The vertical axis is the average queue length. A green dot represents the number of aircraft in the queue counted at each touchdown. The black bold lines are the average queue lengths according to our model (Eq. (8)). The black dotted lines are from a Monte Carlo simulation, where we averaged the queue lengths over 500 runs. We see that the model accurately describes the observed data. Its advantage is to predict the time to return to equilibrium for future traffic intensities. These are in the order of a few hours. For comparison, the gray line is from an M/D/1 queue. Here the equilibrium queue length is significantly larger than in our system. From theory, exponential convergence rates are known for finite Markov chains (e.g. Durrett, 2010), but more analytical work is needed to analyze the transient

properties of the PSRA process. Moreover, it turned out that nonsymmetrical jumps (q – 1=2) matched the simulations more accurately than symmetrical ones. Until now we could not find an explanation for this.

4.2. Steady state analyis of PSRA and M/D/1 The left panel of Fig. 8 shows histograms and densities of metering delays. The gray bars show the empirical data. The light-blue and yellow bars are from the obscured Poissonian and PSRA queues. The parameters were: arrival rate 0.32 (ac/min), deterministic service at 10 NM (0.63 ac/min, based on average ground speed) and obscured metering delay with   Nð0; 0:8Þ (min, i.i.d). The arrival time error of the PSRA model was n  Nð0; 5Þ (min, i.i.d.), and its pre-variance 7.5 min. Note that we also performed simulations with compact support arrival time error (uniform distribution). For an analytical solution of the queueing model, this assumption is important, but in the simulations, the Gaussian model fitted better with the empirical data. Both models fit equally well the empirical data. The ‘true’ shape of the density functions, that is with infinitesimal bins and no obscured term, is shown by the dark-blue and dotted black lines. We computed the density for the M/D/1 queue analytically from Shortle and Brill (2005) and approximated it for PSRA arrivals by a histogram with bin-size 0.1 min. Both densities include an atom at w ¼ 0 min and three segments: one segment reaches until a peak at w ¼ 1:6 min, the metering distance (the service time). Another segment is determined by another service time and lies between 1:6 < w 6 3:2 min. It is slightly concave. The third segment covers the tail. It turns out that with higher traffic intensities, the 1.0

4. Queueing results

0.4

0.6

P(W
0.4 0.3

0.2

0.2 0.1

Data M D+ε 1 PSRA D + ε 1 PSRA/D/1

0.0

0.0

Density

0.5

0.8

0.6

Data M D+ε 1 PSRA D + ε 1 M/D/1 PSRA/D/1

0

1

2

3

4

w (min) n= 3688

5

6

7

0

1

2

3

4

5

w (min) n= 2794

Fig. 8. Metering delay. Left: Relative frequencies and density functions. Right: Distribution functions.

6

7

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

1.0 0.8 0.6 0.4 0.2

90% quantile P(w
ρ = 0.5 σ = 5 ρ = 0.5 σ = 0.5 ρ = 0.95 σ = 0.5 pre−variance=5

0.0

P(W
Similar delay patterns

0

1

2

3

4

w Fig. 9. Waiting time distribution functions.

5

5.0 4.50 0 4.00 3.50

3.00

0.9 2.50

TMA 2. 00

0.8

0.7

Congestion does not increase on these lines

En−route

1.0

0

0.6 0.50

segments increase in size and get more and more concave. More details on ‘zig-zag’ densities in the M/D/1 queue can be found in Shortle and Brill (2005). The right panel of Fig. 8 shows the probability distribution functions PðW 6 wÞ. The bold gray line is the empirical data. The yellow and blue lines are again from the obscured queueing models PSRA=D þ =1 and M=D þ =1. Both fit the data quite well. The dotted black line is from the original PSRA model. We can see again the different segments. The bending of the curves differs between the first segment w  1:6 min, and the second one 1:6 < w 6 3:2 min. This time, the distribution is continuous, but at w ¼ 1:6 min (the service time), the bending changes, which corresponds to the discontinuity in the density function. In order to keep the figure interpretable, we do not show the distribution function of the M/D/1 model. As can be inferred from the densities, it resembles closely the one for the PSRA arrivals. From both figures we conclude that PSRA and Poissonian arrivals have comparable queueing behavior in the case of our instance. Fig. 9 displays the waiting time distributions under three uncertainty scenarios. The current situation is plotted in black. If the traffic intensity remains as today, but the prediction error r will reduce to 0:5 min, the metering delays will almost disappear (blue). But if the intensity increases, similar delay patterns as today will occur (red). This means that radar vectoring cannot be eliminated completely in the future. The dotted red curve has pre-variance 5 min, whereas the bold ones have pre-variance 0. We found in several simulations, that the pre-variance does not affect the delays strongly. But one question is, if two processes with same rate and pre-variance produce the same queueing delay. We generated processes with same pre-variances for Uniform and Exponential distributions. These models (pre-variance 5 and 7.5 min) produced the same delays. But more analysis of the relationship between variance and queueing delay is necessary. The queueing delay in PSRA systems depends mainly on the traffic intensity q and the arrival time errors n. A decision rule on how to balance ground and airborne delays can be obtained by looking at the relationship between these variables. Fig. 10 shows a contour-plot of the quantile function w90 ¼ f ðq; rÞ, where r is the standard deviation of n. For example, the line w90 ¼ 2:00 in the figure has the meaning that metering delay is smaller or equal than 2 min on this line in 90% of the cases. We display the 90% quantiles PðW 6 w90 Þ ¼ 0:9, but average values would produce similar results. We can see that f is an increasing function in both dimensions. In general, the contour-lines represent the combinations of

ρ

272

0.25

0.75

1

2

3

4

σ Fig. 10. Contour plot of 90% waiting-time quantiles.

traffic intensity and prediction error that keep the congestion on a constant level. As one can see, the largest quantile takes the value 5 min for all traffic and 2.5 min for en-route traffic (q < 0:75). These values were obtained for a pre-variance of 5 min. As mentioned before, the current balance between ground and airborne metering delays is 10 min airborne delay. We infer that a considerable rebalance in the order of 5 min will be enough in the future. On the other hand, the complete elimination of en-route delay is unlikely because even small uncertainties have the potential to create queues in the vicinity of a very busy airport. As long as runway pressure remains a high priority, radar vectoring remains necessary.

4.3. M/D/1–M/G/1 Until now, we used a deterministic metering distance of 10 NM between all aircraft. Although this is the operational rule, human factors and wind conditions certainly have an impact on it. In fact, closer analysis showed that aircraft pairs with large speed difference were often separated more than 10 NM. Extracted from the radar data, the distribution had an average spacing between successive aircraft of 1.64 min, which corresponds to the 10 NM, based on average ground speed. It was skewed to the right, with the majority of its density in the interval between 1.5 and 2.0 min. Its standard deviation was 4.6 s. The purpose of the analysis in this section was to assess the impact of variation of the metering distance on the generated metering delay. Standard queuing theory predicts higher delay with increase of variability of the service time (Wolff, 1989). We wanted to see if the assumption of constant service time was critical in our context. We thus sampled the service distribution from the radar data in order to understand the impact of varying service rates in our models. Table 2 summarizes the aggregated delay probabilities. For reference, the first three rows are from the radar data. For the flows from Central and South Japan, they are similar (rows 1, 2). Those from International (row 3) appear to be smaller, except for the first bin and the extreme tail. The two rows in the middle compare M/D/1 and M/G/1 models. It turns out that they behave equally. Both have problems with the core part of the real distribution and cover well the tail. The next two rows show the results from the obscured queueing models.

273

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275 Table 2 Delay frequencies per region. Origin

Central South Int’l M/D/1 M/G/1 M/D+/1 M/G+/1 PSRA/D+/1

Delay (min)

Mean

[0, 1]

(1, 2]

(2, 3]

(3, 4]

(4–5]

(5–6]

>6

0.60 0.60 0.70 0.66 0.66 0.61 0.61 0.64

0.22 0.23 0.12 0.18 0.18 0.20 0.19 0.17

0.10 0.11 0.06 0.08 0.08 0.10 0.10 0.08

0.05 0.04 0.02 0.04 0.04 0.05 0.05 0.04

0.02 0.01 0 0.02 0.02 0.02 0.02 0.02

0.01 0.01 0.02 0.01 0.01 0.01 0.01 0.01

0 0 0 0 0 0 0  0:02

Here, we added noise on the waiting times before counting frequencies (as described in subSection 3.2). Roughly, both obscured models M=D þ =1 and M=G þ =1 match better in distribution with the radar data. The average metering delay in the radar data was 1.02 min, 0.95 in the ordinary and 1.1 min in the obscured queueing models. Last but not least, the delays from the PSRA model are printed in the last row. As already seen in the figures, no major difference to the Poissonian ones exist. To conclude, the variations in the service time have no systematic impact on the metering delays. Moreover, theoretical delays under first-come-first-served are the minimal delays, because the separation of 10 NM is constant. Observing that the radar delay is almost identical to the queueing delay was thus not surprising. It also means that the controllers solve the metering task efficiently. 4.4. Other characteristics There was one interesting detail in the radar data: The flow from Central Japan has a slightly higher average delay than the one from South Japan (1.1 min vs. 1.0 min). Moreover, the flow from International destinations has only 0.8 min average delay. But if the metering strategy were first-come-first-served, the delay probabilities should be the same for every aircraft. One reason for these differences seems to be the sector geometry. Some of the six routes allow larger radar vectors than others. And aircraft from Central Japan arrive more frequently on these routes than aircraft from the South. We found a weak pattern in favor of this explanation in the radar data. One can see that the delay frequencies between 3 and 5 min for the flow from Central Japan are slightly higher than in the flow from South Japan (0.05 vs. 0.04 and 0.02 vs. 0.01, bold characters). As far as the flow from International (row 3) is concerned, large delays (5–6 min) occur more often than for the domestic flows (0.02 vs. 0.01) but very small delays occur more frequently (0.7 vs. 0.6). It looks as if controllers were penalizing International flights in some conditions and prioritizing them in other conditions. But these patterns, if statistically significant at all, are weak.

1.1 1.0 0.8 0.92 0.9 1.1 1.1 1.1

pre-scheduled random arrivals (PSRA). They are realistic pointprocess models of air-traffic flows that take into account the regular structure of flight schedules. Their main property is the negative correlation between successive numbers of arrivals. We built a bi-variate Markov Chain to capture these correlations and looked at transient properties of the average queue length. The convergence to equilibrium took place in the order of a few hours. This is encouraging for further analysis and justification for practical use. The second class were Poisson arrivals with deterministic and sampled service-time distribution. Here we reminded the discontinuities in the waiting time density. It turned out that PSRA and Poissonian models behave equivalently during moderate congestion and differ substantially during very high congestion. The advantage of the PSRA model was to assess the impact of reduction in arrival time errors. This is not possible with Poissonian flows, since they are described only by a rate parameter. Based on these results, we derived a strategic rule for delay balancing that improves the fuel efficiency. From the analysis we concluded that it will not be necessary to allocate up to 10 min enroute delay in the future anymore. A rebalance in the order of 5 min will be enough. On the other hand, the complete elimination of en-route delay is unlikely because even small uncertainties have the potential to create queues in the vicinity of a busy airport. As long as runway pressure remains a high priority, radar vectoring, like Point Merge or traditional ones, remains necessary. But before applying our results, we encourage more validation, for example by trajectory-based simulation. As far as the methodology is concerned we confirm that non-standard queueing models increase the technical difficulty considerably, although the typical queueing behavior (i.e. rapid growth of delays as soon as demand approaches capacity) is maintained. Therefore, we derived analytical results where possible and used algorithmic and simulation approaches where necessary. We believe that more descriptive and asymptotic properties of the queueing distributions need to be found in order to make work with queueing theory more practical. Also, the variance of the arrival process as driver for the queueing delays shall be further clarified. Transient analysis is more complex in nature but more results on convergence rates are important as a basis for realistic decision making.

5. Conclusions In this paper we analyzed the most congested traffic flow in Japanese airspace. The current strategy balances delays between the ground and en-route in order to take into account arrival time errors. As long-term strategy, flow managers wish to reduce the amount of en-route delay in order to smoothen the flows. Our approach was empirical, supported by standard and non-classical stochastic queuing models. A radar-data analysis revealed ‘almost’ Poissonian flow characteristics in arrival rate and distribution. The result was expected due to the Japanese pre-scheduling strategy and sector geometry. We then analyzed two classes of queueing models in order to gain more insight into the generation of delays. The first class were

Acknowledgments This work was partly supported by the Electronic Navigation Research Institute (Japan) under the study on Trajectory Modeling and Prediction. The authors are grateful to four anonymous referees for their careful reports and to Masato Fujita for comments on a draft of this paper. Appendix A. Details on goodness-of-fit test Fig. A.11 and Table A.3 show the details in the goodness-of-fit tests performed in Section 2. For the 5 min intervals, we grouped

3.0

freq 0.0

2.0

3.0

2.0

3.0

3.0

0.0

freq 0.0

2008/10/25 D: 1.7 p: 0.43

1.0

2.0

2.0

3.0

1.0

2.0

3.0

2008/10/22 D: 19.6 p: 0

0.0 0.1 0.2 0.3 0.4

freq

freq 2.0

1.0

1.0

2008/08/23 D: 15 p: 0

2008/10/20 D: 4.7 p: 0.1

0.0 0.1 0.2 0.3 0.4

1.0

0.0

freq 0.0

2008/08/24 D: 16 p: 0

0.0

3.0

0.0 0.1 0.2 0.3 0.4

freq

freq 1.0

2.0

2008/08/20 D: 5.5 p: 0.07

0.0 0.1 0.2 0.3 0.4 0.0

1.0

0.0 0.1 0.2 0.3 0.4

2.0

2008/08/19 D: 17.2 p: 0

3.0

0.0 0.1 0.2 0.3 0.4

1.0

0.0 0.1 0.2 0.3 0.4

freq 0.0

0.0 0.1 0.2 0.3 0.4

freq

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275

0.0 0.1 0.2 0.3 0.4

274

0.0

2008/12/23 D: 9.2 p: 0.01

1.0

2.0

3.0

2008/12/24 D: 5.8 p: 0.06

Fig. A.11. Distribution of number of arrivals in 5 min intervals on 9 different days.

Table A.3 Goodness of fit with re-sample (k = 4). Date

D2

p ðv2k2 Þ

p (re-sample)

^ k

n

20080819 20080820 20080823 20080824 20081020 20081022 20081025 20081223 20081224

17.17 5.46 14.99 15.95 4.67 19.55 1.68 9.16 5.76

0 0.07 0 0 0.1 0 0.43 0.01 0.06

0 0.07 0 0 0.1 0 0.45 0.01 0.06

1.31 1.61 1.52 1.38 1.57 1.56 1.58 1.53 1.58

211 259 245 222 252 251 255 247 255

the right tails of the distribution into one bin, leading to four bins in total. For 10 and 15 min intervals, we used the natural number of occurred bins. In order to verify our simulations for the distribution of the test statistics in the other tests, we computed it for the v2 as well. Here, it fitted with the v2k1l , which is the asymptotic one for l estimated parameters. It was interesting to realize that the consequence of estimated parameters is a shift in mean of the test distribution; making the same measured test statistics less likely.

References Artiouchine, K., Baptiste, P., & Dürr, C. (2008). Runway sequencing with holding patterns. European Journal of Operational Research, 189(3), 1254–1266.

Balakrishnan, H., & Chandran, B. (2010). Algorithms for scheduling runway operations under constrained position shifting. Operations Research, 58(6), 1650–1665. Ball, M., Barnhart, C., Dresner, M., Hansen, M., Neels, K., Odoni, A., et al. (2010). Total delay impact study. Technical report, The National Center of Excellence For Aviation Operations Research (Nextor). Barnhart, C., Belobaba, P., & Odoni, A. R. (2003). Applications of operations research in the air transport industry. Transportation Science, 37(4), 368–391. Bäuerle, N., Engelhardt-Funke, O., & Kolonko, M. (2007). On the waiting time of arriving aircrafts and the capacity of airports with one or two runways. European Journal of Operational Research, 177(2), 1180–1196. Bayen, A., Tomlin, C., Ye, Y., & Zhang, J. (2004). An approximation algorithm for scheduling aircraft with holding time. Proceedings of the 43rd IEEE Conference on Decision and Control, 3, 2760–2767. Beasley, J. E., Krishnamoorthy, M., & Sharaiha, Y. M. (2000). Scheduling aircraft landings – the static case. Transportation Science, 34(2), 180–197. Begen, M. A., & Queyranne, M. (2011). Appointment scheduling with discrete random durations. Mathematics of Operations Research, 36(2), 240–257. Bell, G. E. (1949). Operational research into air traffic control. Journal of the Royal Aeronautical Society, 53, 965–976. Bloomfield, P., & Cox, D. (1972). A low traffic approximation in queues. Journal of Applied Probability, 9, 832–840. Cox, D. (1955). The statistical analysis of congestion. Journal of the Royal Statistical Society. Series A, 118(3), 324–335. Daley, D., & Vere-Jones, D. (2003). An introduction to the theory of point processes (2 edition., Vol. 1). Springer. Dunlay, W., & Horonjeff, R. (1976). Stochastic properties of enroute air traffic – an empirical investigation. Journal of Aircraft, 13(5), 376–381. Durrett, R. (2010). Probability: Theory and examples (4th ed.). Cambridge University Press. Erzberger, H. (1995). Design principles and algorithms for automated air traffic management. In AGARD lecture series no. 200 on knowledge-based functions in aerospace systems. Falin, G. (2010). A single-server batch arrival queue with returning customers. European Journal of Operational Research, 201(3), 786–790. Fan, T., & Odoni, A. (2002). A practical perspective on airport demand management. Air Traffic Control Quarterly, 10(3), 285–306.

C. Gwiggner, S. Nagaoka / European Journal of Operational Research 235 (2014) 265–275 Green, L. V., & Savi, S. (2008). Reducing delays for medical appointments: A queueing approach. Operations Research, 56(6), 1526–1538. Guadagni, G., Lancia, C., Ndreca, S., & Scoppola, B. (2013). Queues with exponentially delayed arrivals. ArXiv e-prints 1302.1999. Guadagni, G., Ndreca, S., & Scoppola, B. (2011). Queueing systems with prescheduled random arrivals. Mathematical Methods of Operations Research, 73(1), 1–18. Gwiggner, C., Kimura, A., Nagaoka, S. (2009). Data and queueing analysis of a Japanese arrival flow. In: Proceedings of the Asia-Pacific International Symposium on Aerospace Technology. Gifu, Japan. Hansen, M., Lovell, D. T., Odoni, A., & Vlachou, K. (2009). Use of queuing models to estimate delay savings from 4D trajectory precision. In: Proceedings of the 8th USA/Europe ATM R&D seminar. Napa, CA, US. ICAO (2005). Global air traffic management operational concept. International civil aviation organization. Janic, M. (2000). Air transport system analysis and modeling: Capacity, quality of services and economics. Amsterdam: Gordon and Breach Science Publishers. JPDO, J.P.D.O. (2010). Concept of operations for the next generation air transportation system. Version 3.1. Kendall, D. G. (1953). Stochastic processes occurring in the theory of queues and their analysis by the method of the imbedded markov chain. The Annals of Mathematical Statistics, 24(3), 338–354. Knessl, C., & Yang, Y. P. (2002). An exact solution for an M(t)/M(t)/1 queue with time-dependent arrivals and service. Queueing Systems, 40, 233–245. Krishnamoorthy, A., Babu, S., & Narayanan, V. (2009). The MAP/(PH/PH)/1 queue with self-generation of priorities and non-preemptive service. European Journal of Operational Research, 195(1), 174–185. Lucantoni, D. M. (1993). The BMAP/G/1 queue: A tutorial. In L. Donatiello & R. Nelson (Eds.), Performance evaluation of computer and communication systems (pp. 330–358). Springer. Lecture Notes in Computer Science. Lulli, G., & Odoni, A. (2007). The European air traffic flow management problem. Transportation Science, 41(4), 431–443.

275

Mercer, A. (1960). A queueing problem in which the arrival times of the customers are scheduled. Journal of the Royal Statistical Society, Series B, 22(1), 108–113. Mercer, A. (1973). Queues with scheduled arrivals: A correction, simplification and extension. Journal of the Royal Statistical Society, Series B, 35(1), 104–116. Newell, G. F. (1982). Applications of queueing theory (2nd ed.). Chapman & Hall. Nikoleris, T., & Hansen, M. (2012). Queueing models for trajectory-based aircraft operations. Transportation Science, 46(4), 501–511. Odoni, A., Morisset, T., Drotleff, W., & Zock, A. (2011). Benchmarking airport airside performance: FRA vs. EWR. In Proceedings of the 9th Europe-USA ATM seminar. Germany: Berlin. Peterson, M., Bertsimas, D., & Odoni, A. (1995). Models and algorithms for transient queueing congestion at airports. Management Science, 41(8), 1279–1295. Pyrgiotis, N., Malone, K. M., & Odoni, A. (2013). Modelling delay propagation within an airport network. Transportation Research Part C: Emerging Technologies, 27(2), 60–75. Sesar-Consortium (2007). The ATM target concept. D3: The Sesar Consortium. Shortle, J., & Brill, P. H. (2005). Analytical distribution of waiting time in the M/iD/1 queue. Queueing Systems, 50, 185–197. Stephens, M. A. (1974). Edf statistics for goodness of fit and some comparisons. Journal of the American Statistical Association, 69(347), 730–737. Stephens, M. A. (1993). Aspects of goodness-of-fit. Technical report, Department of Statistics, Stanford University. Tu, Y., Ball, M., & Jank, W. (2008). Estimating flight departure delay distributions-a statistical approach with long-term trend and short-term pattern. Journal of the American Statistical Association, 103(481), 112–125. Winsten, C. (1959). Geometric distributions in the theory of queues. Journal of the Royal Statistical Society: Series B (Methodological), 21(1), 1–35. Wolff, R. W. (1989). Stochastic modeling and the theory of queues. New Jersey: Prentice-Hall. Yang, Y., & Knessl, C. (1997). Asymptotic analysis of the M/G/1 queue with a timedependent arrival rate. Queueing Systems, 26, 23–68.