\spacedallcapsAsymptotics for the Late Arrivals Problem

\spacedallcapsAsymptotics for the Late Arrivals Problem

Carlo Lancia Leiden University Gianluca Guadagni University of Virginia Sokol Ndreca Universidade Federal de Minas Gerais Benedetto Scoppola Università di Roma ‘Tor Vergata’

We study a discrete time queueing system where deterministic arrivals have i.i.d. exponential delays . The standard deviation of the delay is finite, but much larger than the deterministic unit interarrival time. We describe the model as a bivariate Markov chain, we prove that it is ergodic and then we focus on the unique joint equilibrium distribution. We write a functional equation for the bivariate generating function, finding the solution of such equation on a subset of its set of definition. This solution allows us to prove that the equilibrium distribution of the Markov chain decays super-exponentially fast in the quarter plane. Finally, exploiting the latter result, we discuss the numerical computation of the stationary distribution, showing the effectiveness of a simple approximation scheme in a wide region of the parameters. The model, motivated by air and railway traffic, was proposed many decades ago by Kendall [46] with the name of “late arrivals problem”, but no complete solution has been found so far.

1 Introduction

In this paper we consider a single-server queue with deterministic service time, which is assumed of unitary length for the sake of simplicity. The th customer arrives to the system at time


where are i.i.d. exponential random variables with parameter .

In the limit the point process (1.1) weakly converges to a Poisson process of parameter , whereas for fixed the arrivals are negatively autocorrelated, see [21, 38] and references therein.

Remark 1.

Although the results in [38] are stated under the hypothesis that the probability density function of the delays has compact support, this assumption does not play any role in establishing the convergence to a Poisson process. As a consequence, the very same result applies here too.

We study the system described above for fixed and we assume that arriving customers might balk with independent probability . In other words, each customer can be deleted independently with probability before joining the queue. Besides being a mathematical expedient that ensures the existence of a stationary state111See Lemma 2 below., the balking is mainly a way to model empty intervals in a constant stream of customers. Again, the point process (1.1) with balking weakly converges to a Poisson process, but with parameter  [38]. In what follows, we refer to the balked version of (1.1) as Exponentially Delayed Arrivals (EDA).

Service can be delivered by the unique server only at discrete times. The length of the queue at time is ; it represents the number of customers waiting to be served, including the customer that will be served precisely at time , if any. Due to the balking, it is immediate to see that the traffic intensity of the system is given by ; see [38] for details. Using Kendall’s notation we hereafter refer to the queue model described so far as .

The model is motivated by the description of public and private transportation systems, including buses, trains, aircraft [11, 38, 39, 42] and vessels [35, 43], appointment scheduling in outpatient services [10, 18, 53, 54] crane handling in dock operations [23, 28], and in general any system where scheduled arrivals are intrinsically subject to random variations. Preliminary results show that the model described above fits very well with actual data of inbound air traffic over a large hub, see [17, 50].

The appearance of the stochastic point process (1.1) can be traced back to Winsten’s seminal paper [66]. Winsten named such a queueing model late process and obtained results for the special case and service time exponentially distributed. At the end of [66] there is a discussion on Winsten’s results by Lindley, Wishart, Foster, and Takács, where they state that Winsten’s paper can be considered as the first treatment of a queueing model with correlated arrivals [66, pages 22-28].

The same problem was investigated also by Kendall. In [46, page 11] he remarked the great importance of systems with arrivals like (1.1): “[…] perhaps too much attention has been paid to rather uninteresting variations on the fundamental Poisson stream. As soon as one considers variations dictated by the exigencies of the real world, rather than by the pursuit of mathematical elegance, severe difficulties are encountered; this is particularly well illustrated by the notoriously difficult problem of late arrivals.” Kendall also provided the following elegant interpretation: if the random variables are non-negative then the process defined by (1.1) is the output of the stationary queueing system [46]. In particular, if the random variables are exponentially distributed then can be viewed as a 2-stage tandem queueing network

However, this is not the approach followed in this work.

Some years later, under the hypothesis that , Nelsen and Williams exactly characterised in [56] the distribution of the inter-arrival time intervals and their correlations. They also gave an explicit expression of these quantities in the particular case of ’s exponentially distributed.

After the ’70s only approximations of the arrival process (1.1) [15, 64] or numerical studies of its output [5, 11, 60] seem to have appeared in the literature. In particular, in [38] the authors presented a self-contained study of an arrival process like (1.1), assuming for a compact-support distribution. They also proposed an approximation scheme that keeps the correlation of the arrivals and is able to compute in a quite accurate way the quantitative features of the queue. To the best of our knowledge, a queueing system with arrivals described by (1.1) still remains an open problem and the best results obtained so far are due to Winstein in 1959 [66].

is an example of a queueing system with correlated arrivals, a subject broadly studied in past years. There are many ways to impose a correlation to the arrival process. For instance, the parameters of the process may depend on their past realisation, as in [24], or on some on/off sources, as in [67]. Another relevant example of a queue model with correlated arrivals is the so-called Markov Modulated Queueing System. In Markov Modulated Queueing Systems the parameters are driven by an independent external Markovian process, see [2, 9, 22, 52, 61] and references therein. Our model shares with Markov Modulated Queueing Systems the property that one can define an external and independent Markovian process that drives the arrival rates. However, we see in Section 2 that the output of this external drive has also a deterministic effect on the queue length. More precisely, can be interpreted as an independent drawing from an external pool of customers late at time , see (2.2) below. Due to the memoryless property of the exponential delays, each customer late at time will be still in the pool at time independently with probability . This leads to binomial transitions in the number of late customers.

In Section 2 we show that can be described as a bivariate Markov chain representing the queue length and the number of late customers. We prove that such a bivariate chain is ergodic and write the balance equations of the stationary distribution, finding a functional equation for the bivariate generating function.

There exists an extensive literature about two-dimensional Markov models. Many methods for attacking the problem are available under two assumptions, namely, spatial homogeneity and finiteness of at least one marginal chain, see [12, 31, 36, 51, 55, 57, 59]. Unfortunately, the Markov chain defined in Section 2 does not satisfy any of the mentioned requirements.

When both components of the Markov chain are infinite but space homogeneity is still ensured, the problem is typically attacked by reduction to a Riemann-Hilbert Boundary Value Problem. These may be solved, for example, by the uniformisation technique [47], conformal mappings [19, 20, 29], the compensation method [3], or the Power Series Approximation [13, 14, 48, 41]. The aforesaid binomial transitions are responsible for the lack of spatial homogeneity and are often encountered in Mathematical Biology [8, 16, 25].

To the best of our knowledge, the functional equation (2.15) below has never previously appeared in the literature. Yet it is possible to mark some analogies with the functional equation in [26, 27, 45], the most important being that in both equations the right hand side exhibits the generating function computed in a convex combination in the parameter . Other examples of chains with binomial transitions may be found in [1, 6, 58, 69].

In Section 3 we study the marginal distribution of the number of late customers and we obtain its exact analytical expression, which reveals the rich combinatorial structure of the problem. This intermediate result allows us to show that the stationary distribution of the queue has a super-exponential decay. Finally, in Section 4 we show that such a super-exponential behaviour enables a simple, yet very effective, numerical approximation scheme of the system balance equations. For a wide range of the system parameters, including typical values for real traffic applications of the model, we give a very good a priori estimate of the total-variation distance between the true and the approximate solution.

2 Stationary distribution: generating function and balance equations

Let us consider the process , which describes the length of the queue at time . This process is governed by the stochastic recursion


where is the number of arrivals in the interval according to the arrival process (1.1), and is the usual Kronecker’s delta. The term represents the action of the server in decreasing the queue length by one if at time the queue is non-empty. Since the service time is deterministic, we focus on the so-called embedded process by observing the system at , i. e. at departure instants.

The quantity depends in general on the whole previous history of the system. Indeed, if for some large value of , for any then is large with great probability. Conversely, if in the recent past the values of have been large then is expected to be small. This suggests that the arrival process is negatively autocorrelated, as proven in [38]. Hence, the recursion (2.1) does not depend only on the present value of , and the memory of the process is infinite since can be arbitrarily large.

Let us now denote by the number of customers that are late at time , that is to say,


Let us next define and . Given the value of , the random variable is binomial with parameters and . According to the memoryless property of the exponential delays , the number of unit time intervals in which a customer is late is a geometric random variable with parameter . In other words, a customer will be late for consecutive time intervals with probability . Hence, the embedded process , , is a discrete-time Markov chain. If the customer scheduled in the interval has balked then



All in all,


We describe the queue by the embedded process , , a discrete-time Markov chain with the following transition probabilities:


Figure 1 displays the transitions of the embedded chain in the quarter plane.

Figure 1: Transitions of the queueing system in the quarter plane. They happen along the lines of Cartesian equation and if , and along the lines of Cartesian equation and if . Green transitions happen with probability (no balking), red transitions happen with probability (balking).

We now address the existence of the stationary state of the embedded chains and .


Lemma 1 The chain is ergodic if and only if .


If , the chain is clearly irreducible. In order to prove its positive recurrence, we use Foster’s criterion [63, Cor. 8.7] setting as the Lyapunov function. Thus, we need to show that there exist suitable positive constants such that

  1. for ;

  2. for ;

  3. the set is finite.

First, by (2.4),


Therefore, point 1 is satisfied, for instance, by the choice and . Next, Figure 1 shows that at each iteration increases at most by one unit and point 2 is satisfied:


By definition of point 3 is also fulfilled. On the other hand, for the chain is not ergodic because it is no longer irreducible222For , the chain satisfies in fact .. ∎


Lemma 2 The bivariate chain is ergodic if and only if and .


If and , the bivariate chain is irreducible, see Figure 1. Let us consider the process , which represent the diagonal in the quarter plane where the point lies on, see Figure 1. The bivariate chain has the property that . Equations (2.5)-(2.6) yield


In order to prove the positive recurrence of , we use again Foster’s criterion setting with . From (2.7) and (2.9)-(2.12),

Then, using a little algebra, it can be shown that the first point of Foster’s criterion is satisfied by and (for example, ). Point 2 of Foster’s criterion holds because

where we have used the property that has only nearest-neighbour transitions and equation (2.8). Point 3 is fulfilled by simply considering the definition of . For , the bivariate chain is not ergodic because it is no longer irreducible333For , the process satisfies in fact .. ∎

For , Lemma 2 guarantees both the existence and the uniqueness of the stationary distribution


Let us consider the following bivariate generating function:


We are now ready to prove the main result of this section.


Theorem 3 The bivariate generating function (2.14) satisfies


where .

Remark 2.

The functional equation (2.15) does not admit simple or immediate solutions. It is radically different from the functional equations typically studied in the literature [20, 30] and it is rather special in this respect. A simple solution can be found only in the particular case , see Section 3 below.

Proof of Theorem 2.

For each , the balance equations of are the following:


where are given by (2.3) and we agree that whenever . The special cases and respectively lead to


To show that (2.16)–(2.18) hold, it suffices to write in terms of and then neglect the time dependency. Take for example (2.18): the system is found at time in state , i. e. with empty queue and no late customers, only if at time it was either in state or in state , and the th scheduled customer444Cf. formula (1.1). is deleted by thinning. Indeed, if at time the system was in state then nothing happens and the state remains unchanged, whereas if it was in state then the customer in queue is served and at time the system is in state . Similarly, there are four cases such that the system is found at time in state , i. e. with an empty queue and customers late. In the first two cases the system is in state or in state at time , the th customer is deleted, and no one of the late customers arrives in the interval (this event has in fact probability ). In the remaining cases the system is in state or in state at time , the th customer is not deleted555Therefore, the th customer is added to the set of the customers that are already late., and no one of the late customers arrive in the interval . The latter argument gives (2.17) while an easy generalisation to the case leads to (2.16).

Let us take (2.16), multiply both sides by , and then sum over and . The summation of all terms multiplied by yields

or equivalently,

The change of indices and yields

to which we still have to sum the contribution

All in all, the sum of all terms multiplied by is

In a completely analogous way we can compute the sum of the terms multiplied by , which turns out to be

Summing up the two contributions, we get (2.15). ∎

Remark 3.

We end this section with a discussion of the special case . In this regime, the right-hand side of equation (2.15) does not depend on anymore, and . The number of late customers is in fact as the th customer can not have a delay . Then, equation (2.15) yields directly


where is the stationary probability of a void queue. Therefore, equation (2.19) is equivalent to

which is the classical result of a queue with balking.

3 The marginal distribution of late customers

In this section we focus on , the marginal distribution of late customers. First, we iterate the functional equation (2.15) to obtain the generating function of in the form of an infinite product. Then, we invert the generating function and find the exact analytical expression of . Finally, we derive the asymptotic behaviour of and use it to infer asymptotics for .

The marginal distribution of late customers is

and its generating function is

Setting into equation (2.15) yields


Evaluating the last equation in yields

Iterating (3.1) times,

The limit of for exists for each and , see [4]. Moreover, Therefore, we have proven the following


Corollary 4 For and ,

Remark 4.

The infinite product (3.2) has an interesting combinatorial interpretation:


where is the infinite -Pochhammer symbol, also known as infinite -ascending factorial in . For , , where is the well-known Euler function.

Remark 5.

For , the right-hand side of (3.3) is analytic for each . Therefore, the power series , convergent for each , can be analytically continued in the whole complex plane. As such, we expect the marginal distribution to decrease super-exponentially fast in .

Following the insight given by Remark 5, we shift now the focus to the asymptotic behaviour of and . Expanding the product and rearranging it in powers of yields


where is the number of partitions of in distinct parts.

The following theorem holds:


Theorem 5 Let be the equilibrium marginal distribution of the number of late customers and its generating function. Then,


Theorem 3 is a direct consequence of the following two results from number theory.


Lemma 6 [68] If then the number of partitions of in distinct parts equals the number of partitions of into at most parts (not necessarily distinct).


Lemma 7 [7, 40] Let be the number of partitions of in parts that do not exceed . Then equals the number of partitions of into at most parts and

Proof of Theorem 3.

Using Lemma 3 and 3 we can recast (3.4) as

where, as usual, and when .

The second part of the Theorem is proved from (3.5) as follows:

Now we use (3.6) to obtain an upper bound on : let be the -Pochammer symbol of the pair . The following inequalities hold:



where from (3.7) to (3.8) we have used the properties of -ascending factorials and -binomial coefficients [32]. If is sufficiently large then , and (3.8) yields

Remark 6.

Theorem 3 and (3.9) show that, asymptotically in , the leading order of is . This fact can be directly implied from arrival process (1.1). In fact, the most likely way to have late customers ( large) is that each of the customers originally scheduled in the interval do not balk and are late at time , an event of probability .

Since , we have just obtained the following asymptotic result:


Theorem 8 Uniformly in , the equilibrium distribution decays super-exponentially fast in . More precisely,


As a matter of fact the super-exponential decay of may be proved asymptotically in . Let us consider the auxiliary process , which we have already encountered in the proof of Lemma 2. There we have interpreted as the diagonal in the quarter plane where the point lies on. Under equilibrium conditions, the probability of finding the system on the th diagonal is just

It is straightforward to prove that the generating function of is .

Substituting into (2.15) yields

Remark 7.

Equation (3.11) gives an interesting connection between the equilibrium distribution of the quantity and the stationary probability of having late customers given that the queue is void. Figure 1 shows that the latter event drives the dynamic of through an independent Bernoulli random variable with parameter , which explains the factor .

From (3.11) we can compute as follows in terms of :


For , formulas (3.10) and (3.12) then yield


Therefore, the following asymptotic result holds:


Theorem 9 The equilibrium distribution decays super-exponentially fast as either or . More precisely,


4 Numerical approximation of the joint stationary measure

In this final section we examine the possibility to approximately compute the joint stationary distribution . Due to the very broad range of applications of the queueing model , an efficient approximate computation of the solution may prove itself crucial in contexts where practical solutions are needed.

In Section 3 we have shown that the joint stationary probability decreases super-exponentially fast in the limit of either . The natural question arising is then whether a bare truncation of the infinite linear system (2.16)-(2.18) is sufficient to obtain a satisfactory numerical expression of . As we will see, in this case the answer is positive due to (3.14).

We truncate the infinite system of balance equations (2.16)-(2.18) by fixing an integer and imposing for . For the purpose of simplifying the notation, we map the quarter plane onto the set of non-negative integers. This way, we can relabel the unknowns as and recast (2.16)-(2.18) as . For the details of both mapping and relabeling, see Appendix A. Next, we map onto the set of non-negative integers , where . We want to consider the truncated system


The idea we present is not new and has been already discussed, for instance in [65] for stationary distributions with geometric tail. As shown in Appendix B, there exists a sequence such that

The following a priori estimate of the error introduced by the truncation can be obtained from perturbation theory [34, §2.6.2]:


where is the matrix


is the usual Kronecker’s delta, and

is the norm-1 condition number of the matrix .

From Appendix B,


Figure 2 shows the behaviour of as a function of for . Therefore, unless the condition number of the matrix is very large, we expect that a truncation at the level will give a very good approximation of the stationary probabilities of .

Figure 2: Behaviour of as a function of for .

Estimating the condition number of a matrix is a notably difficult problem and a vast literature exists on this topic, see e. g. [33, 62, 44]. Since the aim of the present section is showing that a bare truncation of the balance equations (2.5)–(2.6) may be sufficient for an approximate computation of , we fall back on numerical computations. Figure 3 displays the value of in the -plane when . We see that the condition number is not larger than for .

Figure 3: for when and vary between and in steps. The condition number is of order at most in this region of the parameters.

Table 1 gives the numerical value of the right-hand side of (4.4) for between and , and . Comparison of Table