Asymptotics for the Late Arrivals Problem

We study a discrete time queueing system where deterministic arrivals have i.i.d. exponential delays $\xi_{i}$. The standard deviation $\sigma$ of the delay is finite, but its value is much larger than the deterministic unit service 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 with the name of"late arrivals problem", but no solution has been found so far.


introduction
In this paper we study a one server queueing system defined in the following way. The i-th customer arrives to the system at time where ξ i are i.i.d. exponential random variables, with density f ξ (t) = λe −λt if t > 0 0 otherwise In the limit λ → 0 this arrival process weakly converges to a Poisson process whereas for λ small but fixed the arrivals are negatively autocorrelated, see [30] and references therein. In this paper we study the system for fixed λ.
In order to ensure tightness for the queueing process we assume throughout the paper an independent thinning to the arrival process: each customer can be removed/deleted before joining the queue with independent probability 1 − ρ. In this case it is not difficult to show that the arrival process converges weakly to the Poisson process as above, see again [30]. After Kendall (see below for his exact words) we will name this arrival process exponentially delayed arrivals. Service can be delivered by the unique server only at discrete times. The length of the queue n t at time t is the number of customers waiting to be served, including the customer that will be served precisely at time t, if any. Due to the thinning procedure, it is immediate to see that the traffic intensity of the system is given by ρ; see [30] for details. This model is motivated by the description of public and private transportation systems including buses, trains, aircrafts [9,30,31,33] and vessels [28,34], appointment scheduling in outpatient services [8,16,42,43] crane handling in dock operations [20,25], 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 [15]. This model 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 realization, as in [21], or on some on/off sources, as in [57]. Among queueing systems with correlated arrivals Markov modulated queueing systems (MMQS) are quite relevant. In MMQSs the parameters are driven by an independent external Markovian process, see [2,7,19,41,49,51] and references therein. As we shall see in Section 2, our model shares with MMQSs the property that one can define an external and independent Markovian process that drives the arrival rates. However, the output of this external drive also determines the evolution of the queue length. More precisely, our system can be seen as a single-server queue with deterministic service time and arrivals given by the reneging from an auxiliary queue, namely the customers that are late at time t -see (2.2) below. Due to the memoryless property of the exponential delays, each customer late at time t may arrive in the unit time slot (t, t+1] indipendently and with the same probability q = e −λ . This means that the aforesaid renenging only happens at integer times and that customers perform synchronous independent abandonments, leading to binomial transitions in the number of late customers. In Section 2 we will see that our model can be described as a bidimensional Markov chain, whose components are the queue length and the number of late customers, respectively. There exist an extensive literature about two-dimensional Markov models. Many methods for attacking the problem are available under the assumtions of either a spatial homogeneity property and that one of the two marginal chains is finite, see [10,27,29,40,44,46,48]. Unfortunately, the Markov chain defined in Section 2 does not satisfy any of the mentioned hypothesis. When both components of the Markov chain are infinite but space homogeneity is ensured, the problem is typically attacked by reduction to a Riemann-Hilbert boundary value problem. Boundary value problems represent a broadly studied subject and several techniques have been developed in the last decades to solve them: uniformization technique [37], conformal mapping [17,18,26], compensation method [3] and power series approximation (PSA) [11,12,32,38]. Usually, PSA is used to obtain the generating function in terms of a power serie in the load ρ, though different parameters may be used, as in [55]. A power series expansion in a parameter different from ρ is also the strategy we choose in this paper. The second constraint that our model fails to meet, namely the spatial homogeneity, is due to the reneging with binomial transitions mentioned above. This kind of transitions are often encountered in Mathematical Biology [6,14,22]. The functional equation for the boundary value problem (2.8) below and its counterpart in [23,24,35] are radically different due to the mismatched dimensionalities of the systems, yet it is still possible to mark some analogies. The most important is that in both equations the right hand side exhibits the generating function computed in a convex combination in the parameter q, that is the probability of each independent abandonment; this feature will play an important role in the proof of Theorem 2. Other examples of chains with binomial transitions may be found in [1,5,47,58].
Although some features of the model may be found in recent literature the problem studied in this paper was proposed and intensively studied for its remarkable importance by the founders of queueing theory in the '60s. The appearance of the stochastic point process (1.1) can be traced back to Winsten's seminal paper [56]. Winsten named such a queueing model late process and obtained results for the special case where customers can never be delayed more than two time-units, i.e. ξ i ∈ [0, 2]. At the end of [56], a discussion by Lindley, Wishart, Foster and Takács on Winsten's results is attached. According to their words, Winsten's paper has to be considered the first treatment of a queueing model with correlated arrivals [56, pg. 22-28]. It is worthy to note that in the framework depicted above Winsten surprisingly got a geometric distribution for the stationary queueing length. During the discussion of Winsten's paper, Wishart [56, pg. 24] states the following: 'If we suppose that the ξ's have the exponential distribution, then I think that we will no longer have the geometric distribution: but what we will have I do not know'. In this paper we give a contribution representing a step towards to the answer of that question. The same problem was investigated also by Kendall in [36, pg. 12], where 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 ξ i are non negative, then the process defined in (1.1) is the output of the stationary D/G/∞ queueing system; in particular, if the ξ i 's are exponentially distributed then the process defined in (1.1) is the output of the stationary D/M/∞ queueing system. Some years later, Nelsen and Williams obtained in [45] the exact characterization for the distribution of the inter-arrival time intervals and correlation coefficient between successive inter-arrival time intervals when random variables ξ i are non-negative. They also gave an explicit expression of these quantities in particular case where the ξ i 's are exponentially distributed.
In particular, in [30] we studied in a self-contained way an arrival process defined as in (1.1), assuming for ξ i a compact support distribution. We 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 a queueing system with this kind of arrivals. Thus, 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 Winsten. In this work we present an equation for the bivariate generating function of the probability distribution of the queue length of the system in the case of exponentially delayed arrivals, and we write a formula leading to the explicit expression of the generating function as of a power series in the parameter q. The paper is organized as follows. In Section 2 we derive step by step an equation for the bivariate generating function. In Section 3 we develop an iterative method to compute the explicit expression for the generating function in terms of a power series. In addition, we discuss some aspects of the solution and presents in detail the first terms of its expansion. Section 4 contains a proof of the analiticity of our generating function with respect to the parameter q. Section 5 is devoted to the discussion of some open questions and some conjectures on the system. We also compare the queue length distribution obtained using the truncated power expansion with an empirical distribution obtained via simulations of (1.1).

construction of the bivariate generating function
The process n t , describing the length of the queue at time t, is governed by the relation where m (t,t+1] is the number of arrivals in the interval (t, t + 1], and the term 1 − δ n t ,0 represents the fact that if at time t there are clients in the queue then the first of the queue is served, decreasing the length of the queue by 1.
It is evident that the quantity m (t,t+1] depends in general on the whole previous history of the system. If for some large value of T we have that m (s,s+1] = 0 for any s ∈ {t − T, t − T + 1, ..., t − 1}, then we will have a great probability that m (t,t+1] is large; conversely, if in the recent past the values of m (s,s+1] have been large, then m (t,t+1] is expected to be smaller. This is equivalent to claim that the arrival process is negatively autocorrelated, for a proof of this claim see [30]. It is then understood that the evolution (2.1) does not depend only on the present value of n t , indeed the memory of the process is infinite since T can be arbitrarily large. Nonetheless the delays are exponential, i.e. memoryless. Let us denote with l t the number of customers that have not arrived yet at time t, that is to say Then, given the value of l t , the random variable m (t,t+1] is binomially distributed. More precisely, using the notation m (t,t+1] = m t for the sake of simplicity, we have that if the customer at time t has been deleted by the thinning procedure, while Hence, let us assume that the state of the system is determined by l t and by the length of the queue n t . Relations (2.3) and (2.4) ensure the fact that the number of arrivals does not depend on the previous history of the system, because m t depends only on l t . The evolution of the pair (n t , l t ) is therefore Markovian. The resulting Markov chain is clearly ergodic, therefore the stationary measure P n,l exists and it is unique, provided ρ < 1. Consider the following bivariate generating function P(z, y) = ∑ n,l≥0 P n,l z n y l |z|, |y| ≤ 1 (2.7) We are now ready to prove the main result of this section: theorem 1 The bivariate generating function P(z, y) defined in (2.7) satisfies Proof. We can write a set of equation for the stationary probabilities P n,l computing the probabilities to have P n t+1 ,l t+1 in terms of the various P n t ,l t , and then neglecting the time dependency. We obtain P 0,l = ρ(P 1,l−1 + P 0,l−1 )b 0,l + (1 − ρ)(P 1,l + P 0,l )b 0,l (2.9) P 1,l = ρ [(P 1,l + P 0,l )b 1,l+1 + P 2,l−1 b 0,l ] + (1 − ρ) [(P 1,l+1 + P 0,l+1 )b 1,l+1 + P 2,l b 0,l ] (2.10) P 2,l = ρ [(P 1,l+1 + P 0,l+1 )b 2,l+2 + P 2,l b 1,l + P 3,l−1 b 0,l ] + (1 − ρ) [(P 1,l+2 + P 0,l+2 )b 2,l+2 + P 2,l+1 b 1,l+1 + P 3,l b 0,l ] (2.11) . . . where b j,l are defined in (2.3) and (2.4). Defining P l (z) = ∞ ∑ n=0 P n,l z n we can rewrite the system above in the compact form Eventually, defining b l (z) = ∑ ∞ n=0 b n,l z n we arrive to b l (z) = (q + zp) l . Using the latter expression into (2.12) we get (2.8).
This equation represents a non trivial boundary value problem that will be studied in the next section. We conclude the present section with the following Remark 1. The marginal distribution for l, i.e., P(1, y) can be obtained by a recursive substitution of P(1, y) in the expression (2.8). We obtain This infinite product is well known in combinatorics with the name of q-Pochhammer symbol of the pair (ρ(y − 1), q). It is the q-analog of the descending factorial (a.k.a. Pochhammer symbol), and when ρ(y − 1) is substituted by 1 it represents the well known Euler's function. It is worthy to outline that in the circle |q| < 1 the q-Pochhammer symbol is analytic for all the values of y.

expansion in powers of q
We look for the solution of the boundary value problem (2.8) in the class of the functions which are analytic in the region If a solution is found in this class then by the uniqueness of the stationary distribution P n,l it must be the unique solution of (2.8). Consider the following power expansion which is convergent for 0 ≤ q < φ. An estimate of φ is given in Section 4. The function P (k) (z, y) = [q k ]P(z, y) is the k-th coefficient of the power expansion, namely Let υ = z + q(y − z) then the equation (2.8) for the bivariate generating function can be rewritten as follows This can be combined with a Taylor's expansion of P(z, y) around y = z P(z, y) = ∑ j≥0 (y − z) j j! ∂ j ∂y j P(z, y) y=z (3.3) For brevity of notation we will denote with the symbol ∂ j y P(·, y 0 ) the j-th order derivative taken with respect to y and evaluated at y = y 0 . From (3.2) we get the following expression for the j-th derivative of P(z, y) ∂ j ∂y j P(z, y) y=z = q j ∂ j y P(0, z) + ∂ j y P(z, z) − ∂ j y P(0, z) z (1 + ρ(z − 1)) Some remarks are now due. Remark 2. If q = 0 then the number l t of customers not arrived yet at time t is identically zero; thus we can infer that the coefficient of the zero-th order, P (0) (z, y), is independent of y. Remark 3. For y = z equation (3.2) becomes P(z, z) = P(0, z) from which it follows Let us define the functions for any j, k ≥ 0, with the agreement that j = 0 means no differentiation and that a k j (z) = 0 whenever j < 0. Substituting (3.1) into (3.6) and rearranging in powers of q we get When q = 0 equation (3.10) becomes and from Remark 2 we get Equations (3.10) and (3.11) then yield to We are now ready to write the fundamental theorem of this section that will allow us to determine the coefficients P (k) (z, y) of the power expansion (3.1). theorem 2 The coefficients P (k) (z, y) satisfy P (0) (z, y) = 1 + ρ(z − 1) (3.14) Proof. Equating (3.1) and (3.3), and using (3.4) we get ∑ k≥0 q k P (k) (z, y) = ∑ j≥0 [1 + ρ(z − 1)] ∂ j y P(0, z) + ∂ j y P(z, z) − ∂ j y P(0, z) z Using the analyticity we can now substitute (3.1) into the right-hand side of (3.16) and rearrange the expression to find that (3.17) Equating the zero-th order terms and using (3.8) we find and (3.14) follows from (3.12) Consider now k ≥ 1 and operate the change of indices i + j = k in (3.17) where from (3.19) to (3.20) we have used (3.8) and the fact that a k j (z) is identically 0 whenever j < 0. Equation (3.20) easily gives (3.15).
Remark 5. Line (3.15) gives a recursive relation for P (k) (z, y) involving the derivatives of P (l) (z, y) only for l < k and P (k) (0, z). In the following we show how to find an explicit expression for P (k) (z, y) using (3.13) and (3.15).

General form of P (k) (0, z)
Let us define the following family of functions According to (3.23) we can rewrite line (3.15) as Then we have the following result lemma 3 For any k ≥ 1, Proof. From (3.13) and (3.24) The value of P (k) (0, 0) can now be used to computed P (k) (0, z). From (3.24) again, we obtain Combining equations (3.24) and (3.25) we get the general form of P (k) (z, y) Remark 8. By means of (3.22) we have that, for k ≥ 2 Therefore for any k ≥ 2 we have that in formulas (3.15), (3.24) and (3.25) the index j actually runs from 1 to k − 1.

Second order in q
From (3.24), (3.25) and Remark 8 we get which gives Note again that P (2) (1, 1) = 0 according to (3.21). Differentiating we obtain Using (3.23) we have Remark 9. Since P (1) (z, y) and P (2) (z, y) are linear in y, we have that In particular, we see that the lowest k such that P (k) (z, y) is more than quadratic in y is k = 6. Such a value could be directly deduced from the memoryless property of the delays. The most likely situation leading to l t = 3 (see (2.2) above) is that none of the customers expected to arrive at times t − 1, t − 2 and t − 3 has arrived yet. This event has probability of order q 6 . Similarly, the lowest k such that P (k) (z, y) is of order m in y is k = m (m+1) 2 .

estimate of the radius of analyticity φ
In this section we provide an estimate of the radius of analyticity of (3.1). Combining the definition (3.7) of the functions a k j (z) with formulas (3.23) and (3.28) we easily derive the following recursive relation, valid for l ≥ 1 while for l = 0 through (3.8), (3.23) and (3.25) we get lemma 4 For all k, l, m ∈ N and |z| ≤ 1 where C = max e 2 , 2 1−ρ .
Proof. By induction. For k = 0, 1 (4.3) is clearly satisfied. In addition, from Remark 8 it suffices to prove (4.3) for l ≤ k. Thus let us assume that (4.3) is fulfilled in the polyhedral lattice P k = (j, l, m) ∈ N 3 , j ≤ k, l ≤ j and show that (4.3) holds in First note that from (3.7), (3.23) and (3.28) it follows very easily that a k l (z) is a polynomial in z with degree not larger than k. Thus, using the inductive hypothesis and a Taylor expansion we have that for ∀ m < k and again m = 0 means that no differentiation is taken.
From (4.2) we now have and (4.5) infers |a k l (z)| ≤ C k provided C ≥ e 2 . The proof of (4.3) for m ≥ 1 is similar.
Lemma 4 implies then that the analyticity radius is φ = C −1 .

open problems and numerical approximations
Although equation (3.15) represents the complete solution of the problem, its use is quite involved, due to the fact that it is a recursive relation on functions. It would be nice to have an explicit expression for P (k) (z, y), or equivalently P (k) (0, y), involving only a (linear) recursion on a set of constants. Alternatively, it would be nice to have a similar tool to compute relevant quantities like the mean value of the queue and its higher momenta. Both questions remain unanswered, but it is clear that the problem has a deep combinatory structure. Here we want briefly mention that the evidence of such combinatorial structure emerges also in the easy case of the marginal distribution represented by P(1, y). According to (3.15) we can write for P (k) (1, y) the following expression And we have that where D(l, j) is the number of partition of the integer l in terms of j distinct parts, i.e. Hence (5.1) can be rewritten as which proves that (5.1) is equivalent to Last equality can be easily verified by means of Young diagrams, see for example [54]. Finally, we want to state a conjecture that seems to be very reasonable looking at the first terms of the expansion in q. We conjecture that the general form of P(z, y) should be with β k,j constants to be determined.
Another aspect of this model worthy of being studied is the convergence to the stationary measure of the system. It is very likely that the system presents the phenomenon of the cutoff ( for a discussion of the cutoff see for instance [39,53] and references therein). A proof of such cutoff can be probably achieved with the ideas in [39], using the parameter α t = n t + l t .
When the system starts from a very high value of n t , it is easy to prove (see [30] for more details) that α t decreases by one when the customer t is deleted by the thinning procedure, and it remains constant otherwise. Moreover in this case α t is very far from its stationary distribution, leading to a large total variation distance of the distribution with respect to the stationary one. Once α t arrives to the values on which its stationary distribution is supported, the convergence to the equilibrium should be very fast. This will be the subject of further studies.
The last part of this section is devoted to some numerical aspects of the model. In Section 3 we explicitly found the first four coefficients P (k) (z, y) of the expansion (3.1). A question that may naturally arise now is 'how good is approximation found this way?'. It turns out that the answer depends on the values of ρ and q. From (3.31), (3.37) and (3.46) we see that the convergence of the truncated series ∑ N k=0 P (k) (z, y) q k to P(z, y) may be slow when ρ is close to 1 due to the presence of factors (1 − ρ) −1 and (1 − ρ) −2 . On the other hand we obviously have that the smaller is q the faster the aforesaid convergence. Numerical simulations can give an idea of the quality of the approximation of the truncated expansion. The probability distribution of the queue length can be found as the marginal of the empirical distribution P n,l or by means of the following formula P n = 1 n! ∂ n ∂z n P(z, y) (z,y)=(0,1) (5.7) = 1 n! N ∑ k=0 q k ∂ n ∂z n P (k) (z, y) (z,y)=(0,1) + o(q N ) (5.8) Figure 5 shows a comparison between the distribution of the queue length empirically obtained via numerical simulations of an Exponentially Delayed Arrivals process and the distribution obtained by means of (5.8).
It is immediate to see that if the traffic index ρ is close to 1 then even small values of q are sufficient for huge discrepancies to arise between the simulated and the theoretical distributions. Conversely, if the average load of the system is kept at a moderate level then we see that the model is capable to deliver a very good accordance with the simulated data for scenarios where each customer has a great probability to arrive out of its pre-scheduled slot. Figure 1: Comparison of P n = ∑ l P n,l for interesting scenarios. For each different couple of values (ρ, q) the 'Empirical' serie represents the empirical distribution while 'Theoretical' is obtained by means of (5.8). references