Chapter 11 Queueing models

A queueing model describes a system in which customers, jobs, people, packets, or tasks arrive, wait if necessary, receive service, and then leave. The word customer is used generically. A customer could be a person standing in line, a phone call, a web request, a patient, or a machine part waiting for repair.

A typical queueing model specifies:

  1. an arrival process,
  2. a service-time distribution,
  3. the number of servers,
  4. the rule for selecting the next customer to serve.

In this chapter we focus on the simplest continuous-time queueing models. These models combine the Poisson process, exponential waiting times, and continuous-time Markov processes.

11.1 Kendall notation

Queueing models are often described with Kendall notation. The basic form is \[ A/S/k, \] where:

  • \(A\) describes the arrival process,
  • \(S\) describes the service-time distribution,
  • \(k\) is the number of servers.

The letter \(M\) stands for Markovian, which in this context means memoryless. Thus:

  • \(M/M/1\) means Poisson arrivals, exponential service times, and one server.
  • \(M/M/k\) means Poisson arrivals, exponential service times, and \(k\) identical servers.

The default service discipline in this chapter is first-come-first-served. For the basic formulas below, the exact order of service often matters less than the fact that the service times are exponential and customers are not lost.

System size. In this chapter \(X_t\) denotes the total number of customers in the system at time \(t\). This includes both customers waiting in line and customers currently receiving service.

11.2 The \(M/M/1\) queue

The \(M/M/1\) queue has:

  • arrivals according to a Poisson process with rate \(\lambda\),
  • service times that are independent exponential random variables with rate \(\mu\),
  • one server.

The rate \(\lambda\) is the average number of arrivals per unit time. The rate \(\mu\) is the average number of customers one server can serve per unit time. Thus the average service time is \(1/\mu\).

The process \(X_t\) records the number of customers in the system. If \(X_t=n\), then the possible transitions are \[ n\to n+1\quad\text{at rate }\lambda \] from an arrival, and \[ n\to n-1\quad\text{at rate }\mu,\qquad n\geq 1, \] from a completed service.

There are no service completions from state \(0\) because the system is empty.

This is a birth-death process on the state space \[ \{0,1,2,3,\ldots\}. \] The birth rate is constant, \(\lambda_n=\lambda\), and the death rate is \[ \mu_n=\mu,\qquad n\geq1. \]

\(M/M/1\) queue. The \(M/M/1\) queue is the continuous-time Markov process with transitions \[ n\to n+1\text{ at rate }\lambda, \qquad n\to n-1\text{ at rate }\mu\quad(n\geq1). \]

11.3 Stability

The traffic intensity of the \(M/M/1\) queue is \[ \rho=\frac{\lambda}{\mu}. \] This compares the arrival rate to the service capacity. If \(\lambda<\mu\), then the server can keep up with arrivals in the long run. If \(\lambda\geq\mu\), then the queue does not have a stationary distribution.

Stability condition for \(M/M/1\). The \(M/M/1\) queue has a stationary distribution exactly when \[ \rho=\frac{\lambda}{\mu}<1. \] Equivalently, the arrival rate must be smaller than the service rate.

This condition is intuitive. If customers arrive faster than they can be served, then the line should grow without bound. If the arrival rate is smaller than the service rate, then the system can settle into statistical equilibrium.

11.4 Stationary distribution for \(M/M/1\)

Assume \(\lambda<\mu\). Let \[ \pi_n=\lim_{t\to\infty}\mathbb P(X_t=n) \] be the stationary probability of having \(n\) customers in the system.

The equilibrium flow between neighboring states gives \[ \pi_n\lambda=\pi_{n+1}\mu, \] so \[ \pi_{n+1}=\frac{\lambda}{\mu}\pi_n=\rho\pi_n. \] Therefore \[ \pi_n=\rho^n\pi_0. \] Since the probabilities must sum to one, \[ 1=\sum_{n=0}^\infty \pi_n =\pi_0\sum_{n=0}^\infty \rho^n =\frac{\pi_0}{1-\rho}. \] Thus \[ \pi_0=1-\rho, \] and \[ \pi_n=(1-\rho)\rho^n, \qquad n=0,1,2,\ldots. \]

Stationary distribution for \(M/M/1\). If \(\rho=\lambda/\mu<1\), then the stationary distribution is geometric: \[ \pi_n=(1-\rho)\rho^n, \qquad n=0,1,2,\ldots. \]

The probability that the server is busy is \[ \mathbb P(X>0)=1-\pi_0=\rho. \] Thus the traffic intensity \(\rho\) is also the long-run fraction of time the server is busy.

11.5 Mean system size and waiting time

Let \(L\) denote the long-run expected number of customers in the system. Since the stationary distribution is geometric, \[ L=\mathbb E(X)=\frac{\rho}{1-\rho}. \] Let \(L_q\) denote the long-run expected number of customers waiting in line, not including the customer in service. Since the server is busy with probability \(\rho\), \[ L_q=L-\rho=\frac{\rho}{1-\rho}-\rho=\frac{\rho^2}{1-\rho}. \]

Little’s law relates average number in system to average time in system: \[ L=\lambda W, \] where \(W\) is the average amount of time a customer spends in the system. Similarly, \[ L_q=\lambda W_q, \] where \(W_q\) is the average amount of time a customer spends waiting in line before service begins.

Therefore, for the \(M/M/1\) queue, \[ W=\frac{L}{\lambda}=\frac{1}{\mu-\lambda} \] and \[ W_q=\frac{L_q}{\lambda}=\frac{\rho}{\mu-\lambda} =\frac{\lambda}{\mu(\mu-\lambda)}. \]

Main \(M/M/1\) performance formulas. For \(\rho=\lambda/\mu<1\), \[ L=\frac{\rho}{1-\rho}, \qquad L_q=\frac{\rho^2}{1-\rho}, \] \[ W=\frac{1}{\mu-\lambda}, \qquad W_q=\frac{\rho}{\mu-\lambda}. \]

These formulas show a common queueing phenomenon: as \(\lambda\) gets close to \(\mu\), the average queue length and average waiting time grow very quickly.

11.6 Example: a single-server queue

Example. Suppose customers arrive at a help desk at rate \(\lambda=4\) per hour, and the employee can serve customers at rate \(\mu=6\) per hour.

The traffic intensity is \[ \rho=\frac{4}{6}=\frac23. \] The system is stable because \(\rho<1\).

The stationary probability of having \(n\) customers in the system is \[ \pi_n=\left(1-\frac23\right)\left(\frac23\right)^n =\frac13\left(\frac23\right)^n. \] The long-run average number of customers in the system is \[ L=\frac{\rho}{1-\rho}=\frac{2/3}{1/3}=2. \] The long-run average time spent in the system is \[ W=\frac{1}{\mu-\lambda}=\frac{1}{6-4}=\frac12 \] hours, or 30 minutes.

The long-run average time spent waiting in line is \[ W_q=\frac{\rho}{\mu-\lambda}=\frac{2/3}{2}=\frac13 \] hours, or 20 minutes.

11.7 Computing \(M/M/1\) quantities in R

The following small calculation plots the stationary distribution for an \(M/M/1\) queue.

lambda <- 4
mu <- 6
rho <- lambda / mu

n <- 0:15
pi_n <- (1 - rho) * rho^n

barplot(pi_n, names.arg = n,
        xlab = "number in system",
        ylab = "stationary probability",
        main = "Stationary distribution for an M/M/1 queue")

Here is a compact function for the basic \(M/M/1\) performance measures.

mm1_summary <- function(lambda, mu){
  rho <- lambda / mu
  if (rho >= 1) stop("The M/M/1 queue is unstable when lambda >= mu.")
  L  <- rho / (1 - rho)
  Lq <- rho^2 / (1 - rho)
  W  <- L / lambda
  Wq <- Lq / lambda
  c(rho = rho, L = L, Lq = Lq, W = W, Wq = Wq)
}

mm1_summary(lambda = 4, mu = 6)
##       rho         L        Lq         W        Wq 
## 0.6666667 2.0000000 1.3333333 0.5000000 0.3333333

11.8 The \(M/M/k\) queue

The \(M/M/k\) queue has:

  • arrivals according to a Poisson process with rate \(\lambda\),
  • service times that are exponential with rate \(\mu\),
  • \(k\) identical servers.

The state \(X_t\) is again the number of customers in the system. Arrivals increase the state by one: \[ n\to n+1\quad\text{at rate }\lambda. \] Service completions decrease the state by one. If \(n\) customers are in the system, then the number of busy servers is \[ \min(n,k). \] Since each busy server works at rate \(\mu\), the total service completion rate is \[ \mu\min(n,k). \] Thus \[ n\to n-1\quad\text{at rate }\mu\min(n,k),\qquad n\geq1. \]

\(M/M/k\) queue. The \(M/M/k\) queue is the birth-death process with transitions \[ n\to n+1\text{ at rate }\lambda, \] \[ n\to n-1\text{ at rate }\mu\min(n,k),\qquad n\geq1. \]

The total service capacity of the system is \(k\mu\). The traffic intensity is \[ \rho=\frac{\lambda}{k\mu}. \] The stability condition is \[ \rho<1, \] or equivalently \[ \lambda<k\mu. \]

11.9 Stationary probabilities for \(M/M/k\)

It is useful to define \[ a=\frac{\lambda}{\mu}. \] For \(0\leq n<k\), the stationary probabilities are \[ \pi_n=\pi_0\frac{a^n}{n!}. \] For \(n\geq k\), all servers are busy, so \[ \pi_n=\pi_0\frac{a^n}{k!k^{n-k}}. \] The value of \(\pi_0\) is determined by normalization: \[ \pi_0=igg[\sum_{n=0}^{k-1}\frac{a^n}{n!} +\frac{a^k}{k!}\frac{1}{1-\rho}\bigg]^{-1}. \]

Stationary distribution for \(M/M/k\). If \(\rho=\lambda/(k\mu)<1\), then \[ \pi_n=\pi_0\frac{a^n}{n!},\qquad 0\leq n<k, \] and \[ \pi_n=\pi_0\frac{a^n}{k!k^{n-k}},\qquad n\geq k, \] where \(a=\lambda/\mu\) and \[ \pi_0=igg[\sum_{n=0}^{k-1}\frac{a^n}{n!} +\frac{a^k}{k!}\frac{1}{1-\rho}\bigg]^{-1}. \]

These formulas are more complicated than the \(M/M/1\) formulas, but the idea is the same: write down the birth and death rates, use equilibrium flow between neighboring states, and then normalize.

11.10 Waiting probability and Erlang C

In an \(M/M/k\) queue, an arriving customer must wait if all \(k\) servers are busy. The stationary probability that all servers are busy is related to the Erlang C formula.

The probability that an arriving customer has to wait is \[ C=\frac{a^k}{k!}\frac{1}{1-\rho}\pi_0. \] This is called the Erlang C waiting probability.

The expected number waiting in line is \[ L_q=C\frac{\rho}{1-\rho}. \] The expected number in service is \(a=\lambda/\mu\), so \[ L=L_q+a. \] By Little’s law, \[ W_q=\frac{L_q}{\lambda}, \qquad W=\frac{L}{\lambda}. \]

Basic \(M/M/k\) performance formulas. Let \(a=\lambda/\mu\) and \(\rho=\lambda/(k\mu)<1\). Then \[ C=\frac{a^k}{k!}\frac{1}{1-\rho}\pi_0, \] \[ L_q=C\frac{\rho}{1-\rho}, \qquad L=L_q+a, \] \[ W_q=\frac{L_q}{\lambda}, \qquad W=\frac{L}{\lambda}. \]

11.11 Example: two servers

Example. Suppose customers arrive at rate \(\lambda=8\) per hour. There are \(k=2\) servers, and each server has service rate \(\mu=5\) per hour.

The traffic intensity is \[ \rho=\frac{\lambda}{k\mu}=\frac{8}{2\cdot5}=0.8. \] The system is stable.

Here \[ a=\frac{\lambda}{\mu}=\frac85=1.6. \] The normalizing constant is \[ \pi_0=\left[1+a+\frac{a^2}{2!}\frac{1}{1-\rho}\right]^{-1}. \] The probability that an arriving customer must wait is \[ C=\frac{a^2}{2!}\frac{1}{1-\rho}\pi_0. \] The expected number waiting in line is \[ L_q=C\frac{\rho}{1-\rho}. \] These quantities are easy to compute in R.

mmk_summary <- function(lambda, mu, k){
  rho <- lambda / (k * mu)
  if (rho >= 1) stop("The M/M/k queue is unstable when lambda >= k * mu.")
  a <- lambda / mu
  n <- 0:(k - 1)
  pi0 <- 1 / (sum(a^n / factorial(n)) + (a^k / factorial(k)) * (1 / (1 - rho)))
  C <- (a^k / factorial(k)) * (1 / (1 - rho)) * pi0
  Lq <- C * rho / (1 - rho)
  L <- Lq + a
  Wq <- Lq / lambda
  W <- L / lambda
  c(rho = rho, pi0 = pi0, wait_probability = C, L = L, Lq = Lq, W = W, Wq = Wq)
}

mmk_summary(lambda = 8, mu = 5, k = 2)
##              rho              pi0 wait_probability                L 
##        0.8000000        0.1111111        0.7111111        4.4444444 
##               Lq                W               Wq 
##        2.8444444        0.5555556        0.3555556

11.12 Comparing one server and two servers

Suppose the arrival rate is \(\lambda=8\) customers per hour. Compare:

  • one fast server with \(\mu=10\),
  • two slower servers with \(\mu=5\) each.

Both systems have the same total service capacity, \(10\) customers per hour. In both cases, \[ \rho=\frac{8}{10}=0.8. \] However, the waiting behavior is not the same.

mm1_summary(lambda = 8, mu = 10)
## rho   L  Lq   W  Wq 
## 0.8 4.0 3.2 0.5 0.4
mmk_summary(lambda = 8, mu = 5, k = 2)
##              rho              pi0 wait_probability                L 
##        0.8000000        0.1111111        0.7111111        4.4444444 
##               Lq                W               Wq 
##        2.8444444        0.5555556        0.3555556

This comparison illustrates an important queueing principle: pooling servers can reduce waiting. A single shared line feeding two servers usually performs better than two separate lines, because one server can help absorb randomness when the other would otherwise be idle.

11.13 Little’s law

Little’s law says \[ L=\lambda W. \] It applies much more broadly than just the \(M/M/1\) and \(M/M/k\) queues. Here:

  • \(L\) is the long-run average number of customers in the system,
  • \(\lambda\) is the long-run arrival rate,
  • \(W\) is the long-run average time a customer spends in the system.

Similarly, \[ L_q=\lambda W_q \] for the queue alone.

Little’s law. Under broad stability assumptions, \[ L=\lambda W. \] Average number in system equals arrival rate times average time in system.

A useful way to remember Little’s law is by units: \[ \text{customers} = rac{\text{customers}}{\text{time}}\cdot \text{time}. \]

11.14 Poisson arrivals see time averages

For queues with Poisson arrivals, arriving customers see the system in its ordinary stationary distribution. This principle is often called PASTA, which stands for Poisson arrivals see time averages.

For example, in an \(M/M/1\) queue, the long-run fraction of time that the system is empty is \[ \pi_0=1-\rho. \] A Poisson arrival also sees an empty system with probability \(1-\rho\).

This principle is special to Poisson arrivals. It may fail for other arrival processes.

11.15 Markov-modulated Poisson arrivals

A Poisson process has a constant arrival rate \(\lambda\). In some applications the arrival rate changes over time. For example, a website might have high-traffic and low-traffic periods, or an insurance company might see claims arrive faster during stormy weather than during ordinary weather.

A simple way to model this is to let the arrival rate depend on a background Markov chain.

Let \(J_t\) be a continuous-time Markov process on states \[ \{1,2,\ldots,m\}. \] When \(J_t=i\), arrivals occur at rate \(\alpha_i\). Thus the arrival process has rate \[ \alpha_{J_t} \] at time \(t\).

Markov-modulated Poisson process. A Markov-modulated Poisson process, or MMPP, is a counting process whose instantaneous arrival rate depends on an underlying continuous-time Markov process. If the background state is \(i\), arrivals occur at rate \(\alpha_i\).

The MMPP is not usually a Poisson process, because the rate is random and correlated over time. Conditional on the path of the background process \(J_t\), the arrivals behave like a Poisson process with time-varying rate \(\alpha_{J_t}\).

11.16 Two-state MMPP

The simplest MMPP has two background states:

  • state \(0\): low arrival rate \(\alpha_0\),
  • state \(1\): high arrival rate \(\alpha_1\).

Suppose the background chain switches according to \[ 0\to1\quad\text{at rate }\beta, \] \[ 1\to0\quad\text{at rate }\delta. \] The stationary distribution of the background chain is \[ \pi_0^J=\frac{\delta}{\beta+\delta}, \qquad \pi_1^J=\frac{\beta}{\beta+\delta}. \] The long-run average arrival rate is therefore \[ \bar\lambda=\pi_0^J\alpha_0+\pi_1^J\alpha_1. \]

Example. Suppose a system alternates between low and high demand. In the low state, arrivals occur at rate \(\alpha_0=2\) per hour. In the high state, arrivals occur at rate \(\alpha_1=10\) per hour. The background process switches from low to high at rate \(\beta=1\) and from high to low at rate \(\delta=3\).

The stationary probabilities of the background states are \[ \pi_0^J=\frac{3}{4}, \qquad \pi_1^J=\frac{1}{4}. \] The long-run average arrival rate is \[ \bar\lambda=\frac34(2)+\frac14(10)=4. \] Even though the long-run average arrival rate is \(4\), the arrivals are more bursty than a Poisson process with constant rate \(4\).

11.17 An MMPP queue

An MMPP can be used as the arrival process for a queue. The simplest version is an MMPP/M/1 queue. The state must include both the background state and the number of customers in the system: \[ (J_t,X_t). \] If the background is in state \(i\) and there are \(n\) customers in the system, then the transitions are: \[ (i,n)\to(i,n+1)\quad\text{at rate }\alpha_i, \] from an arrival, \[ (i,n)\to(i,n-1)\quad\text{at rate }\mu,\qquad n\geq1, \] from a service completion, and \[ (i,n)\to(j,n)\quad\text{at rate }q_{ij}, \] from a change in the background Markov chain.

The stability condition is more subtle than for the ordinary \(M/M/1\) queue, but the basic intuition is the same. The average service capacity must exceed the long-run average arrival rate: \[ \bar\lambda<\mu. \] For the two-state example above, this would mean \[ \frac{\delta}{\beta+\delta}\alpha_0 +\frac{\beta}{\beta+\delta}\alpha_1 <\mu. \]

The MMPP/M/1 queue is a good example of why Markov process methods are useful. The customer count alone is not Markovian, because knowing whether the system is currently in a high-arrival or low-arrival environment matters. But the pair \((J_t,X_t)\) is Markovian.

State description for an MMPP queue. In a queue with Markov-modulated arrivals, the customer count \(X_t\) alone usually does not contain enough information. The Markov state should include both the background state and the queue size.

11.18 Simulating an \(M/M/1\) queue

The following simulation is included as optional code. It is not evaluated when compiling the book, but it can be run interactively.

sim_mm1 <- function(lambda, mu, tmax){
  t <- 0
  x <- 0
  times <- 0
  states <- 0
  
  while (t < tmax){
    arrival_rate <- lambda
    service_rate <- ifelse(x > 0, mu, 0)
    total_rate <- arrival_rate + service_rate
    
    t <- t + rexp(1, rate = total_rate)
    if (t > tmax) break
    
    if (runif(1) < arrival_rate / total_rate){
      x <- x + 1
    } else {
      x <- x - 1
    }
    
    times <- c(times, t)
    states <- c(states, x)
  }
  
  list(times = times, states = states)
}

set.seed(423)
out <- sim_mm1(lambda = 4, mu = 6, tmax = 20)
plot(out$times, out$states, type = "s",
     xlab = "time", ylab = "number in system",
     main = "Simulated M/M/1 queue")

This event-by-event simulation uses the same idea as other continuous-time Markov processes: from the current state, determine the possible transition rates, wait an exponential time with rate equal to the total rate, and then choose which transition occurs with probabilities proportional to the rates.

11.19 Summary

The \(M/M/1\) and \(M/M/k\) queues are basic birth-death processes. Their main ingredients are Poisson arrivals, exponential service times, and Markovian state dynamics. The \(M/M/1\) queue has especially simple formulas: \[ \pi_n=(1-\rho)\rho^n, \qquad L=\frac{\rho}{1-\rho}, \qquad W=\frac1{\mu-\lambda}. \] The \(M/M/k\) queue extends this idea to multiple servers, with total service rate \(\mu\min(n,k)\) when there are \(n\) customers in the system. Finally, MMPP arrivals provide a natural way to model bursty arrival streams whose rates depend on an underlying Markov environment.