How Long Might We Wait at Random?

TL;DR

This paper investigates the distribution of maximum queue length in a random customer-server model, highlighting the complexity of extending known results from simpler systems and discussing potential generalizations.

Contribution

It provides algebraic expressions for maximum queue length in a two-server system and explores the challenges of extending these to three servers, reviewing related queue wait time distributions.

Findings

Algebraic expressions derived for two-server maximum queue length.

Complexity increases significantly with three servers, making algebraic solutions difficult.

Discussion on generalizing maximum wait time distributions from M/M/1 queues.

Abstract

In discrete time, customers arrive at random. Each waits until one of three servers is available; each thereafter departs at random. We seek the distribution of maximum line length of idle customers. Algebraic expressions obtained for the two-server scenario do not appear feasible here. We also review well-known distributional results for maximum wait time associated with an M/M/1 queue and speculate about their generalization.

Equations

$$\begin{array}[c]{l}\text{{u = 0}}\\ \text{{m = 0}}\\ \text{{for t=1:n}}\\ \text{{\ \ x = rand()$<$p \ \ \ \ \ \ \ \ \# x=1 means that an arrival occurs}}\\ \text{{\ \ y1 = rand()$<$r \ \ \ \ \ \ \ \# y1+y2+y3=3 means that three departures occur}}\\ \text{{\ \ y2 = rand()$<$r \ \ \ \ \ }}\\ \text{{\ \ y3 = rand()$<$r \ \ \ \ \ }}\\ \text{{\ \ if u==0 }}\\ \text{{\ \ \ \ u = x\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \# increment is 1 or 0\ \ }}\\ \text{{\ \ else}}\\ \text{{\ \ \ \ if u==1}}\\ \text{{\ \ \ \ \ \ u = max(0,u+x-y1) \ \ \ \ \ \ \ \ \ \# increment is 1, 0 or -1}}\\ \text{{\ \ \ \ elseif u==2}}\\ \text{{\ \ \ \ \ \ u = max(0,u+x-y1-y2)\ \ \ \ \ \ \ \# increment is 1, 0, -1 or -2}}\\ \text{{\ \ \ \ else}}\\ \text{{\ \ \ \ \ \ u = max(0,u+x-y1-y2-y3)\ \ \ \ \# increment is 1, 0, -1, -2 or -3}}\\ \text{{\ \ \ \ end}}\\ \text{{\ \ end}}\\ \text{{\ \ m = max(m,u)}}\\ \text{{end}}\\ \text{{return m}}\end{array}$$

$$\begin{array}[c]{ccc}\mathbb{P}\left\{\text{time lapse between adjacent arrivals is }i\right\}=p\,q^{i-1},&&i\geq 1\end{array}$$

$$\begin{array}[c]{ccc}\mathbb{P}\left\{\text{time lapse between adjacent departures is }j\right\}=r\,s^{j-1},&&j\geq 1\end{array}$$

$$\mathbb{P}\left\{M_n<k\right\}\sim \exp \left(-\frac{{\pi }_k}{\mathbb{E}(C)}n\right)$$

$$\left({\small\begin{array}[c]{cccccccc}q&p&0&0&0&0&0&\cdots\\ qr&pr+qs&ps&0&0&0&0&\cdots\\ qr^{2}&pr^{2}+2qrs&2prs+qs^{2}&ps^{2}&0&0&0&\cdots\\ qr^{3}&pr^{3}+3qr^{2}s&3pr^{2}s+3qrs^{2}&3prs^{2}+qs^{3}&ps^{3}&0&0&\cdots\\ 0&qr^{3}&pr^{3}+3qr^{2}s&3pr^{2}s+3qrs^{2}&3prs^{2}+qs^{3}&ps^{3}&0&\cdots\\ 0&0&qr^{3}&pr^{3}+3qr^{2}s&3pr^{2}s+3qrs^{2}&3prs^{2}+qs^{3}&ps^{3}&\cdots\\ 0&0&0&qr^{3}&pr^{3}+3qr^{2}s&3pr^{2}s+3qrs^{2}&3prs^{2}+qs^{3}&\cdots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}}\right)$$

$$\begin{array}[c]{ccccc}\pi_{j}=\pi_{j}^{\ast}\pi_{3}&&\text{for}&&0\leq j\leq 2,\end{array}$$

$${\pi }_3={\left(\frac{1}{1-\omega }+{\pi }_0^{\ast }+{\pi }_1^{\ast }+{\pi }_2^{\ast }\right)}^{-1},$$

$$\begin{array}[c]{ccccc}\pi_{j}=\omega^{j-3}\pi_{3}&&\text{for}&&j\geq 4\end{array}$$

$${\pi }_2^{\ast }=\frac{1}{p{s}^2}\left[\left(1-3pr{s}^2-q{s}^3\right)-\left(3p{r}^{2}s+3qr{s}^2\right)\omega -\left(p{r}^3+3q{r}^{2}s\right){\omega }^2-\left(q{r}^3\right){\omega }^3\right],$$

$${\pi }_1^{\ast }=\frac{1}{ps}\left[\left(1-2prs-q{s}^2\right){\pi }_2^{\ast }-\left(3p{r}^{2}s+3qr{s}^2\right)-\left(p{r}^3+3q{r}^{2}s\right)\omega -\left(q{r}^3\right){\omega }^2\right],$$

$${\pi }_0^{\ast }=\frac{1}{p}\left[\left(1-pr-qs\right){\pi }_1^{\ast }-\left(p{r}^2+2qrs\right){\pi }_2^{\ast }-\left(p{r}^3+3q{r}^{2}s\right)-\left(q{r}^3\right)\omega \right]$$

$$\omega =\frac{1}{3qr}\left[\left(-3+2r+3ps\right)+(3q-r)r{\left(\frac{2}{\chi +\theta }\right)}^{1/3}-{\left(\frac{\chi +\theta }{2}\right)}^{1/3}\right]$$

$$\begin{array}[c]{ccc}\chi=-9qr^{2}+2r^{3}-27q^{2}s,&&\theta=\sqrt{4(3q-r)^{3}r^{3}+\chi^{2}}.\end{array}$$

$${\left(\frac{1}{\omega }\right)}^k=n{\left(\frac{1}{\omega }\right)}^{h+1}$$

$${\pi }_{k}n={\pi }_3{\omega }^{k-3}n={\pi }_3{\omega }^{h-2}.$$

$$\left\{\begin{array}[c]{lll}-3&&\text{with probability }qr^{3},\\ -2&&\text{with probability }pr^{3}+3qr^{2}s,\\ -1&&\text{with probability }3pr^{2}s+3qrs^{2},\\ 0&&\text{with probability }3prs^{2}+qs^{3},\\ 1&&\text{with probability }ps^{3}.\end{array}\right.$$

$${\nu }_j=p{s}^3{\nu }_{j-1}+(3pr{s}^2+q{s}^3){\nu }_j+(3p{r}^{2}s+3qr{s}^2){\nu }_{j+1}+(p{r}^3+3q{r}^{2}s){\nu }_{j+2}+q{r}^3{\nu }_{j+3},$$

$${\nu }_0=p{s}^3{\nu }_{-1}+(3pr{s}^2+q{s}^3)+(3p{r}^{2}s+3qr{s}^2){\nu }_1+(p{r}^3+3q{r}^{2}s){\nu }_2+q{r}^3{\nu }_3$$

$$F(z)$$

$$=p{s}^{3}z\sum _{j=1}^{\infty }{\nu }_{j-1}{z}^{j-1}+(3pr{s}^2+q{s}^3)\sum _{j=1}^{\infty }{\nu }_j{z}^j+\frac{3p{r}^{2}s+3qr{s}^2}{z}\sum _{j=1}^{\infty }{\nu }_{j+1}{z}^{j+1}$$

$$+\frac{p{r}^3+3q{r}^{2}s}{z^2}\sum _{j=1}^{\infty }{\nu }_{j+2}{z}^{j+2}+\frac{q{r}^3}{z^3}\sum _{j=1}^{\infty }{\nu }_{j+3}{z}^{j+3}$$

$$=p{s}^{3}z\left[F(z)+1\right]+(3pr{s}^2+q{s}^3)F(z)+\frac{3p{r}^{2}s+3qr{s}^2}{z}\left[F(z)-{\nu }_{1}z\right]$$

$$+\frac{p{r}^3+3q{r}^{2}s}{z^2}\left[F(z)-{\nu }_{1}z-{\nu }_2{z}^2\right]+\frac{q{r}^3}{z^3}\left[F(z)-{\nu }_{1}z-{\nu }_2{z}^2-{\nu }_3{z}^3\right]$$

$$\left[1-p{s}^{3}z-3pr{s}^2-q{s}^3-\frac{3p{r}^{2}s+3qr{s}^2}{z}-\frac{p{r}^3+3q{r}^{2}s}{z^2}-\frac{q{r}^3}{z^3}\right]F(z)$$

$$=p{s}^{3}z-\frac{3p{r}^{2}s+3qr{s}^2}{z}({\nu }_{1}z)-\frac{p{r}^3+3q{r}^{2}s}{z^2}({\nu }_{1}z+{\nu }_2{z}^2)-\frac{q{r}^3}{z^3}({\nu }_{1}z+{\nu }_2{z}^2+{\nu }_3{z}^3)$$

$$\left[q{r}^3+(p{r}^3+3q{r}^{2}s)z+(3p{r}^{2}s+3qr{s}^2)z^2-(1-3pr{s}^2-q{s}^3)z^3+p{s}^3{z}^4\right]F(z)$$

$$=-p{s}^3{z}^4+(3p{r}^{2}s+3qr{s}^2)z^2({\nu }_{1}z)+(p{r}^3+3q{r}^{2}s)z({\nu }_{1}z+{\nu }_2{z}^2)+q{r}^3({\nu }_{1}z+{\nu }_2{z}^2+{\nu }_3{z}^3)$$

$$(1-z)\left[q{r}^3+(3q-2r+3pr)r^{2}z+(1+s+s^2-3p{s}^2)r{z}^2-p{s}^3{z}^3\right]F(z)$$

$$=-p{s}^3{z}^4+3p{r}^{2}s{z}^3{\nu }_1+3qr{s}^2{z}^3{\nu }_1+p{r}^3{z}^2{\nu }_1+p{r}^3{z}^3{\nu }_2+3q{r}^{2}s{z}^2{\nu }_1+3q{r}^{2}s{z}^3{\nu }_2+q{r}^{3}z{v}_1$$

$$+q{r}^3{z}^2{\nu }_2+z^3\left[{\nu }_0-p{s}^3{\nu }_{-1}-(3pr{s}^2+q{s}^3)-(3p{r}^{2}s+3qr{s}^2){\nu }_1-(p{r}^3+3q{r}^{2}s){\nu }_2\right]$$

$$=z^3{\nu }_0+\left(qr+prz+3qsz\right)r^{2}z{\nu }_1+q{r}^3{z}^2{\nu }_2-p{s}^3{z}^3{\nu }_{-1}-3pr{s}^2{z}^3-q{s}^3{z}^3-p{s}^3{z}^4.$$

$${\nu }_{-j}=p{s}^3{\nu }_{-j-1}+(3pr{s}^2+q{s}^3){\nu }_{-j}+(3p{r}^{2}s+3qr{s}^2){\nu }_{-j+1}+(p{r}^3+3q{r}^{2}s){\nu }_{-j+2}+q{r}^3{\nu }_{-j+3},$$

$${\nu }_0=p{s}^3{\nu }_{-1}+(3pr{s}^2+q{s}^3)+(3p{r}^{2}s+3qr{s}^2){\nu }_1+(p{r}^3+3q{r}^{2}s){\nu }_2+q{r}^3{\nu }_3$$

$$G(z)$$

$$=\frac{p{s}^3}{z}\sum _{j=1}^{\infty }{\nu }_{-j-1}{z}^{j+1}+(3pr{s}^2+q{s}^3)\sum _{j=1}^{\infty }{\nu }_{-j}{z}^j+(3p{r}^{2}s+3qr{s}^2)z\sum _{j=1}^{\infty }{\nu }_{-j+1}{z}^{j-1}$$

Full text

00footnotetext: Copyright © 2019 & 2022 by Steven R. Finch. All rights reserved.

In queueing theory, maximum line length and maximum wait time are different sides of the same coin, one spatial (Section 1) and the other temporal (Section 2).

1 In Spatium

Let $0<r<1$ and $0<p<3r$.  Consider the Julia program:

[TABLE]

which simulates the maximum value of a Geo/Geo/$3$ queue with LAS-DA over $n$ time steps.  The Boolean expressions containing Julia’s Uniform $[0,1]$ random deviate generator ensure that $X\sim\,$Bernoulli($p$) and $Y\sim\,$Bernoulli($r$).  The word “Geometric” arises because

[TABLE]

where $q=1-p$ and

[TABLE]

where $s=1-r$.  Clearly $3s-2<q<1$.  LAS stands for “late arrival system” and DA stands for “delayed access” [1]; in particular, a customer entering an empty queue at time $t$ is not immediately eligible for service, but rather at time $t+1$. We study the asymptotic distribution of the maximum $M_{n}$.  This section is best thought of as an extended addendum to [2], written for the sake of completeness.  The only difference with our earlier work is that closed-form expressions here become unwieldy and thus our approach is more numeric and less symbolic.

The Poisson clumping heuristic [3], while not a theorem, gives results identical to exact asymptotic expressions when such exist, and evidently provides excellent predictions otherwise.  Consider an irreducible positive recurrent Markov chain with stationary distribution $\pi$.  For large enough $k$, the maximum of the chain satisfies

[TABLE]

as $n\rightarrow\infty$, where $C$ is the sojourn time in $k$ during a clump of nearby visits to $k$.

Starting with transition matrix

[TABLE]

we obtain [4]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and $0<\omega<1$ satisfies the cubic equation $\omega=(q\omega+p)(r\omega+s)^{3}$, that is,

[TABLE]

where

[TABLE]

Note that, if $k=\log_{1/\omega}(n)+h+1$, we have

[TABLE]

thus

[TABLE]

We need now to calculate $\mathbb{E}(C)$. Consider a random walk on the integers consisting of incremental steps satisfying

[TABLE]

For nonzero $j$, let $\nu_{j}$ denote the probability that, starting from $-j$, the walker eventually hits [math].  Let $\nu_{0}$ denote the probability that, starting from [math], the walker eventually returns to [math] (at some future time).  We have two values for $\nu_{0}$: when it is used in a recursion, it is equal to $1$; when it corresponds to a return probability, it retains the symbol $\nu_{0}$.  Let $j\geq 1$.  Using

[TABLE]

[TABLE]

define

[TABLE]

equivalently

[TABLE]

equivalently

[TABLE]

equivalently

[TABLE]

Examine the denominator of $F(z)$.  Only the first three of its four zeroes $z_{1}$, $z_{2}$, $1$, $z_{3}$ are of interest (the fourth is $>1$).  Substituting $z=1$, $z=z_{1}$ and $z=z_{2}$ into the numerator $N_{F}$ of $F(z)$, then setting $N_{F}=0$, gives three equations in four unknowns.

Let $j\geq 1$.  Using

[TABLE]

[TABLE]

define

[TABLE]

equivalently

[TABLE]

equivalently

[TABLE]

equivalently

[TABLE]

Examine the denominator of $G(z)$.  Only the zero $z_{4}$ of smallest modulus interests us.  Substituting $z=z_{4}$ into the numerator $N_{G}$ of $G(z)$, and setting $N_{G}=0$, gives a fourth equation (to include with the other three from earlier).  For example, if

[TABLE]

we have

[TABLE]

from earlier and

[TABLE]

after solving the simultaneous system in $\nu_{0}$, $\nu_{-1}$, $\nu_{1}$, $\nu_{2}$.  Observe that

[TABLE]

as $n\rightarrow\infty$.  The ratio within the exponential argument is

[TABLE]

and hence

[TABLE]

for sufficiently large $n$, where $\gamma$ denotes Euler’s constant [5].  Such moment formulas usually contain tiny periodic fluctuations, but we omit these from consideration.

The use of an expected maximum for performance analysis, instead of a simple average, does not appear to lead to surprising outcomes.  A corollary of the preceding numerical results is that, in a busy hospital emergency room (with $p=1/3$), one fast doctor (with $r=1/2$) outperforms three slow doctors (each with $r=1/6$).  For average queue lengths [4],

[TABLE]

corresponding to Geo/Geo/$1$ and

[TABLE]

corresponding to Geo/Geo/$3$.  This is also consistent with results in [6] governing deterministic traffic signals: we do better with an $RGRGRG...$ pattern than with $RRRGGG...$.

2 In Tempore

Let $c\geq 1$ be an integer and $0<\lambda<c\mu$.  Consider the R program:

[TABLE]

which simulates the maximum wait time associated with an M/M/$c$ queue over the time interval $[0,n]$.  More precisely, at any arrival time $t=P_{i,1}$, let $W_{sys}$ denote the $i^{\text{th}}$ patient wait time in the system (either queue or treatment) and $W_{que}$ denote the $i^{\text{th}}$ patient wait time in the queue (excluding treatment).  For $c=1$, the maximums of $W_{sys}$ and $W_{que}$ over all arrival times up to $n$ satisfy [7, 8, 9]

[TABLE]

[TABLE]

where $\rho=\lambda/\mu$.  These imply

[TABLE]

[TABLE]

for large enough $n$.  Also, if $\lambda=1/3$ and $\mu=1/2$, then

[TABLE]

[TABLE]

In particular, assuming $n=20000$, the expected maximum wait times are $43.109$ and $40.676$ respectively.

For $c\geq 2$, no analogous formulas are known.  Assuming $n=20000$, $\lambda=1/3$ and $\mu=1/(2c)$, we estimate via simulation that

[TABLE]

for $c=2$ and

[TABLE]

for $c=3$.  With regard to the expected maximum of $W_{sys}$, one fast doctor outperforms $c$ slow doctors.  With regard to the expected maximum of $W_{que}$, however, $c$ slow doctors outperform one fast doctor.  Experimental histograms of such maximums appear to be roughly Gumbel-shaped: a two-moment fit provides a fairly compelling but ultimately inconclusive match between theory and data.  Only rigorous analysis for $c\geq 2$ will clarify the situation.

The same $W_{que}$ trend occurs for simple averages.  Given $\lambda=1/3$, $\mu=1/(2c)$ as before and allowing $n\rightarrow\infty$, we have the following when $c=1$:

[TABLE]

when $c=2$:

[TABLE]

and when $c=3$:

[TABLE]

Such mean results apply not only to the FIFO discipline (first in, first out), but also to the LIFO (last in, first out) and SIRO (serve in random order) disciplines.  The underlying SIRO wait time distribution is more complicated than that for FIFO, as shown for $c=1$ in [10, 11].

Alternative treatments of M/M/$1$ maximum wait times include [12, 13], where $n$ or $k$ refer not to the length of the time interval but instead the total number of customers.  Simple averages of wait times for Geo/Geo/$c$ can be calculated by methods given in [14]; maximums remain open.

3 Acknowledgements

I am thankful to Guy Louchard for introducing me to the Poisson clumping heuristic (especially recursions for $\nu_{j}$ and $\nu_{-j}$), and to Stephan Wagner for extracting discrete Gumbel asymptotics in [15] (a contribution leading to [2, 16, 17]).  See also [18, 19] for work involving the case of Uniform$[a,b]$ service times and Deterministic$[a]$ subcase ($b\rightarrow a$).  The creators of R, Julia, Mathematica and Matlab, as well as administrators of the MIT Engaging Cluster and the MIT Supercloud Cluster, earn my gratitude every day.