High water mark probability of an M/M/1 queue

Question

I'm familiar with the steady-state behavior of an M/M/1 queue. For example, the arrival rate of $\lambda = 2$ per minute and servicing rate of $\mu = 3$ per minute means that utilization $\rho = \lambda / \mu = \frac{2}{3}$, and therefore the steady-state probability in each state $k$ is $(1-\rho)\rho^k$, or the following:

• $\mathrm{Pr}(k=0) = 1/3$
• $\mathrm{Pr}(k=1) = 2/9$
• $\mathrm{Pr}(k=2) = 4/27$
• $\mathrm{Pr}(k=3) = 8/81$

(etc.)

But what if I want to find the probability during some large time interval $T$ (in other words, $\lambda T \gg 1$ and $\mu T \gg 1$) that the queue length never exceeds some amount $k_{\max}$, given that I know nothing about the queue length ahead of time? (assuming probability is its steady-state distribution)

This seems like it should be a well-known problem, but I can't seem to find it in the basic tutorials/textbooks on queueing theory. For a fixed value of $k_{\max}$, the probability should decrease as $T$ increases, and I would expect this to be approximately equal to $e^{-T/T_0}$ as long as $\lambda T \gg 1$ and $\mu T \gg 1$ (and likely requiring $\rho^k \ll 1$ as well), where $T_0$ is some function of $\lambda, \mu,$ and $k_{\max}$, but I don't know how to calculate $T_0$.

Answer 1

The steady-state behaviour is described by the transition rate matrix

$$ \begin{pmatrix} -\lambda & \lambda \\ \mu & -(\mu+\lambda) & \lambda \\ &\mu & -(\mu+\lambda) & \lambda \\ &&\mu & -(\mu+\lambda) & \lambda &\\ &&&\mu&\ddots \end{pmatrix}\;. $$

To find the rate at which the probability not to exceed $k_\text{max}$ decays, replace all states above $k_\text{max}$ by an absorbing state. For instance, for $k_\text{max}=2$, the modified transition rate matrix is

$$ \begin{pmatrix} -\lambda & \lambda \\ \mu & -(\mu+\lambda) & \lambda \\ &\mu & -(\mu+\lambda) & \lambda \\ &&0& 0 \end{pmatrix}\;. $$

This matrix has an eigenvalue $0$, whose eigenvector has probability $1$ in the absorbing state, and $k_\text{max}+1$ negative eigenvalues, of which you want the one closest to zero. I don’t think there’s a closed form for this in general. For low absorption rates, you can approximate the probability remaining in the non-absorbing states by a multiple of the equilibrium distribution of the unmodified system, and then the absorption rate (and thus the decay rate $\frac1{T_0}$) is approximately the rate of transitions from $k_\text{max}$ to $k_\text{max}+1$, i.e. $\lambda(1-\rho)\rho^{k_\text{max}}$.

There’s a factor in front of $\mathrm e^{-\frac T{T_0}}$, though, which is the coefficient of the eigenvector corresponding to the eigenvalue $-\frac1{T_0}$ when the truncated equilibrium distribution of the unmodified system (with all states above $k_\text{max}$ combined into the absorbing state) is decomposed into the eigenvectors of the modified transition rate matrix.

I should add that I don’t know whether this is how this problem is usually treated in queueing theory; this is just how I’d approach it.

Answer 2

Aside: I tried an approach with a truncated transition matrix as proposed by @joriki, but using a numerical simulation to evolve the probability vector $p_k$ to a new probability vector $p'_k$, where the entries $k=0$ to $k_{\max}$ reach equilibrium, and $\frac{d}{dt}p_{k_\max+1} = \lambda p_{k_\max}$, using $\lambda=2$ per minute and $\mu=3$ per minute and $k_\max=100$. (So the vector has length 102, covering $k=0$ to $k=101$.)

What's interesting is if I plot $r_k = p'_k/p_k$, the ratio of the new to the old equilibrium (where $p_k = (1-\rho)\rho^k$), it looks like this:

Graphed vs. time for different values of k, relative to $p_k(t=0)$:

The values of $r_k$ appear to settle down as follows, by around $t=3$ hours (180 minutes):

$$\begin{aligned} r_{100} &= 1/3 = 1-\rho \\ r_{99} &= 5/9 = 1-\rho^2 \\ r_{98} &= 19/27 = 1-\rho^3 \\ r_{97} &= 65/81 = 1-\rho^4 \\ r_{96} &= 211/243 = 1-\rho^5 \\ &\ldots \end{aligned}$$

Which I assume is what I'd find if I went through the algebra a little more carefully.

Anyway it looks like after all the other eigenvalues have died out, $\frac{dq}{dt} = -\lambda(1-\rho)^2\rho^{k_{\max}}q = -q/T_0$ where $q = 1 - p_{k_\max+1}$ is the sum of all the probabilities from $k=0$ to $k_\max$, and the time constant $T_0=1/(\lambda(1-\rho)^2\rho^{k_{\max}})$. For my example, that works out to $\frac{dq}{dt}=-5.465899\times 10^{-19}$ per minute when $q\approx 1$, giving a time constant $T_0$ of about 3.479 trillion years.

For $k_\max=50$, time constant $T_0 \approx 5455$ years.

For $k_\max=30$, time constant $T_0 \approx 1.64$ years (599 days), so that tells me that a queue capacity of 30 is probably enough in this case.

I'm not sure what the penultimate eigenvalue is (the one that caused a settling time of $p'_k$ after about 180 minutes for N=100) but I'm guessing I could find an approximation.

Anyway, the value of $p_{k_\max+1}=1-q \approx 1-e^{-t/T_0}$ is the probability I was asking for in the question.