Introduction

You need to be familiar with previous vignettes about birth-death process.

Basic concepts

Queuing model is mathematical models of waiting lines, or queues, so that the queue lengths and waiting time could be predicted.

Three components

A basic queuing system is a service system where ‘customers’ arrive to one or several servers and require some service from one of them. In studies, the system is usually broken into three parts: arrivals, queue and service, and different assumptions about these basic components will generate different queuing models.

  • Arrival: Arrival models the process in which customers arrive at the queue. Factors to consider include but not limited to: whether the customer source is finite or infinite, whether the customer comes one by one or in groups, whether the arrival of the customers are independent.

  • Queue: Queue models how the customers wait in line after arriving. Different scenarios may involve the capacity of the line: whether a new customer will wait when there are too many customers in the queue and queuing discipline like FCFS (first come first served), LCFS (last come first served), for instance.

  • Service: Service models the process in which the customers is served. Some common assumptions may include how many servers are in the system, and how long it takes to serve one customer for each server.

Kendall’s Notation

Kendall’s notation is typically used to describe and classify a queuing model. A complete notation may include six components A/S/c/K/N/D, where

  • \(A\): the arrival process
  • \(S\): The service time distribution
  • \(c\): The number of servers
  • \(K\): The number of places in the queue
  • \(N\): The calling population
  • \(D\): The queue’s discipline

In this vignette, we only consider the first four components A/S/c/K. Besides, a code is often used to describe the arrival and service process, and some common codes are

  • \(M\): an exponential inter-arrival time or service time.
  • \(D\): a deterministic or fixed inter-arrival time or service time.
  • \(E_k\): an Erlang distribution with k as the shape parameter.
  • \(G\): a general distribution.

For example, in previous vignette, M/M/1 model refers that the inter-arrival time and service time follow an exponential distribution and there is only one server in the system.

For consistency with the previous vignette, we still use the following notations:

  • \(\lambda\): the arrival rate
  • \(\mu\): the service rate
  • \(L\): the number of customers in the system, either waiting in line or in service
  • \(L_q\): the number of customers in the queue
  • \(W\): total waiting time in the system, including the service time
  • \(W_q\): waiting time in the line.

M/M/s model

Motivation

In last vignette, we used M/M/1 system to model the queuing in the hospital, where we assume there’s only one server in the hospital. In reality, doctors in the hospital can treat patients independently, i.e. there are more than one server in the system. This scenario leads us to M/M/s model.

Definition and key formulas

M/M/s queuing model has the same assumption about the inter-arrivals time and service time of each server, which follow exponential distributions with parameter \(\lambda\) and \(\mu\), respectively. But now, there are \(s\) servers in the system. As each customer can only be served by one server, only \(n\) servers will be used if \(n \leq s\), where \(n\) is the number of customers in the system. When \(n > s\), all \(s\) servers will work at an overall service rate \(s\mu\). The arrival and service rate at state \(n\) are

\[ \lambda_n = \lambda \] and \[ \mu_n = \begin{cases} n\mu, \qquad n=1,2,\cdots,s \\ s\mu, \qquad n=s, s+1, \cdots\\ \end{cases} \]

Denote \(\rho = \frac{\lambda}{\mu}\) and \(\rho_s=\frac{\rho}{s} = \frac{\lambda}{s\mu}\). Using the equilibrium solution in previous vignette, we could derive the equilibrium distribution for M/M/s model as

\[ P_n = \begin{cases} \frac{\rho^n}{n!}P_0, \qquad n=1,2,\cdots,s \\ \frac{\rho^n}{s!s^{n-s}}P_0, \qquad n\geq s \end{cases} \] where \[ P_0 = \frac{1}{\sum_{n=0}^{s-1}\frac{\rho^n}{n!}+\frac{\rho^s}{s!(1-\rho_s)}}\\ \]

Then the average length in the queue, in service and in the system are: \[ L_q = \sum_{n=s+1}^{\infty} (n-s) P_n = \frac{\rho^s \rho_s}{s!(1-\rho_s)^2}P_0\\ L_s = \sum_{n=0}^{s-1}nP_n + s\sum_{n=s}^{\infty}P_n = \rho\\ L = L_s + L_q = \frac{\rho^s \rho_s}{s!(1-\rho_s)^2}P_0 + \rho \] Interestingly, we find that the number of customer in service (\(L_s\)) is independent from the number of servers \(s\).

For M/M/s model, we can still use Little’s Law and obtain the waiting time in the system and in the queue as \[ W = \frac{L}{\lambda}, \qquad W_q = \frac{L_q}{\lambda} \]

Simulating M/M/s model

Here, we’ll use a similar logic as previous vignette to simulate a M/M/s system. But notice that now the departure rate is related to the number of customers in the system.

MMs_sim <- function(lambda, mu, s, T_end){
  Z <- c()
  n <- c()
  Z_type <- c()
  Z[1] <- rexp(1, rate=lambda)
  Z_type[1] <- 1
  n[1] <- 1
  i=1
  while (sum(Z) < T_end){
    i=i+1
    if (rev(n)[1] > 0){
      # The departure rate is related to the number of customers in the system
      mu_n <- ifelse(rev(n)[1]<s, rev(n)[1], s) * mu 
      Z[i] <- rexp(1, rate=lambda+mu)
      Z_type[i] <- sample(c(-1,1), size=1, prob=c(mu_n/(lambda+mu_n), lambda/(lambda+mu_n))) 
      n[i] <- n[i-1] + Z_type[i]
    }
    else if (rev(n)[1] == 0){
      Z[i] <- rexp(1, rate=lambda)
      Z_type[i] <- 1
      n[i] <- n[i-1] + Z_type[i]
    }
  }
  return(list(Z=Z, n=n, Z_type=Z_type))
}

Here we set lambda=4, mu=2 and s = c(1,2,3,4), and plot the simulated results as follows.

set.seed(123)
par(mfrow=c(2,2))
lambda=4
mu=2
T_end=500
s <- c(1,2,3,4)
for (i in 1:4){
  res <- MMs_sim(lambda, mu, s[i], T_end)
  plot(cumsum(res$Z), res$n, type='l', xlab='time', ylab='# of customers', main=paste0('lambda=', lambda, ',', 's*mu=', s[i]*mu))
}

From the plots, it’s not hard to conclude that the queue is explosive if \(\lambda \geq s\mu\), and is stationary if \(\lambda < s\mu\)

M/M/s/K model

Motivation

Previously, we assume the capacity of a hospital is finite. However, the space of a hospital is usually limited and can only hold a certain number of patients. Next, we’ll model this scenario using M/M/s/K model.

Definition and key formulas

Denote \(K\) as the capacity of the system. When \(n>K\), customers will leave directly instead of waiting in line. Here, the arrival and service rate at state \(n\) are \[ \lambda_n = \begin{cases} \lambda, \qquad n < K\\ 0, \qquad n \geq K \end{cases} \] \[ \mu_n = \begin{cases} n\mu, \qquad 0 \leq n < s \\ s\mu, \qquad s \leq n < K\\ \end{cases} \] Similarly, we could derive the equilibrium solution as

\[ P_n = \begin{cases} \frac{\rho^n}{n!}P_0, \qquad 0 \leq n < s\\ \frac{\rho^n}{s!s^{n-s}}P_0, \qquad s\leq n \leq K \end{cases} \] where \[ P_0 = \begin{cases} (\sum_{n=1}^{n-1} \frac{\rho^n}{n!} + \frac{\rho^s(1-\rho_s^{K-s+1})}{s!(1-\rho_s)} )^{-1}, \qquad \rho_s \neq 1\\ (\sum_{n=1}^{n-1} \frac{\rho^n}{n!} + \frac{\rho^s}{s!}(K-s+1) )^{-1}, \qquad \rho_s = 1 \end{cases} \\ \]

Using basic linear algebra, the average lengths and waiting times in the system and of the queue can be computed accordingly. For brevity, we ommit the derivation and just give the results.

\[ L_q = \begin{cases} \frac{P_0\rho^s \rho_s }{s! (1-\rho_s)^2} [1-\rho_s^{K-s+1} - (1-\rho_s)(K-s+1)\rho_s^{K-s}], \qquad \rho_s \neq 1 \\ \frac{P_0 \rho^s (K-s) (K-s+1) }{2s!}, \qquad \rho_s = 1 \end{cases} \\ L = L_q + s + P_0 \sum_{n=0}^{n-1} \frac{(n-s)\rho^n}{n!},\\ W = \frac{L}{\lambda(1-P_k)}, \qquad W_q = W - \frac{1}{\mu}\\ \]

Simulating M/M/s/K model

Here, we can simulate M/M/s/K model using the corresponding arrival rate and departure rate. Notice that now both the departure rate and the arrival rate is related to current state. Only departure occurs if the system if full

MMsK_sim <- function(lambda, mu, s, K, T_end){
  Z <- c()
  n <- c()
  Z_type <- c()
  Z[1] <- rexp(1, rate=lambda)
  Z_type[1] <- 1
  n[1] <- 1
  i=1
  while (sum(Z) < T_end){
    i=i+1
    if (rev(n)[1] > 0 & rev(n)[1] < K){
      mu_n <- ifelse(rev(n)[1]<s, rev(n)[1], s) * mu # adjust the current mu
      Z[i] <- rexp(1, rate=lambda+mu)
      # -1 for departure 1 for arrival
      Z_type[i] <- sample(c(-1,1), size=1, prob=c(mu_n/(lambda+mu_n), lambda/(lambda+mu_n))) 
      n[i] <- n[i-1] + Z_type[i]
    }
    # The system is empty
    else if (rev(n)[1] == 0){
      Z[i] <- rexp(1, rate=lambda)
      Z_type[i] <- 1
      n[i] <- n[i-1] + Z_type[i]
    }
    # The system is full
    else if (rev(n)[1] >= K){
      mu_n <- s*mu
      Z[i] <- rexp(1, rate=mu_n)
      Z_type[i] <- -1
      n[i] <- n[i-1] + Z_type[i]
    }
  }
  return(list(Z=Z, n=n, Z_type=Z_type))
}

We use the same setting of parameter as the previous section, and set K=10. The simulated results are as follows:

set.seed(123)
par(mfrow=c(2,2))
lambda=4
mu=2
K = 10
T_end=500
s <- c(1,2,3,4)
for (i in 1:4){
  res <- MMsK_sim(lambda, mu, s[i], K, T_end)
  plot(cumsum(res$Z), res$n, type='l', xlab='time', ylab='# of patients', main=paste0('lambda=', lambda, ',', 's*mu=', s[i]*mu))
}

It’s not surprising that for the previous explosive case ( \(\lambda > s \mu\)), the system is more likely to reach its maximum capacity \(K\), and many customers have to leave instead of entering the line.

Exercise: M/M/1 model with impatient customer (Adapted from Ross Example 8.7)

Unluckily, in daily life, some patients may die in the process of queuing. Here, we consider a simple case: Suppose there is only one server in the system, where the inter-infection times follow an exponential distribution with parameter \(\lambda\) and the service times follow an exponential distribution with parameter \(\mu\). Now, suppose that each patient will only survive an exponential time with parameter \(\alpha\) and leave the line because of death. Then the system can be modeled as

\[ \lambda_n = \lambda, \qquad n \geq 0 \\ \mu_n = \mu + (n-1)\alpha, \qquad n \geq 1 \] Using the previous framework, we can similarly derive the equilibrium solution, which will be left as an exercise for the readers.