Introduction

You need to be familiar with John’s notes and the first vignette.

Motivation

In last vignette, we modeled the initial stage of COVID-19. Now, as time passes by, people are paying more attention to the disease and will go to the hospital for treatment. They will wait in line until they are treated. In this vignette, we want to model the length of the queue and also the average waiting time, which leads us to birth-death process and birth-death queuing model. Here, we assume that every patient wears a mask now and thus the infection rate comes down a lot so that it’s enough to model the infection rate as a constant. For now, we assume the whole hospital can only treat one patient one time, and don’t consider the capacity of hospital and the death of patients. We’ll mention how to model these scenarios in next vignette.

Birth-death process

Birth-death process is a generalization of pure birth process, and has two types of state transitions: ‘birth’, which increases the current state variable by one, or ‘death’, which decreases the current state variable by one.

Postulates

In general, the process \(X(t)\) is called birth-death process if

\[ P\{X(t+h) - X(t) = 1 | X(t) = n \} = \lambda_{n}h + o(h) \\ P\{X(t+h) - X(t) = -1 | X(t) = n \} = \mu_{n}h + o(h) \\ P\{X(t+h) - X(t) = 0 | X(t) = n \} = 1- (\lambda_n+\mu_n) h + o(h) \\ P\{|X(t+h) - X(t)| > 1 | X(t) = n \} = o(h) \\ \] Then, the transition matrix is: \[ \begin{pmatrix} -\lambda_0 & \lambda_0 & 0 & 0 & \cdots \\ \mu_1 & -(\lambda_1+\mu_1) & \lambda_1 & 0 & \cdots \\ 0 & \mu_2 & -(\lambda_2+\mu_2) & \lambda_2 & \cdots \\ 0 & 0 & \mu_3 & -(\lambda_3 + \mu_3) & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{pmatrix} , \] and we can obtain the system equation as follows:

\[ P_0(t+h) = P_0(t)(1-\lambda_0h) + P_1(t)\mu_1h,\\ P_n(t+h) = P_{n-1}(t) \lambda_{n-1} h + P_{n+1}(t)\mu_{n+1} h + P_{n}(t) [1- (\lambda_n+\mu_n)h ]. \]

Differential equation and equilibrium solution

Subtract \(P_n(t)\) from each side of the system equations and further take limit as \(h\) goes to 0. Then we can get the differential equations as \[ \frac{d}{dt}P_0(t) = P_{1}(t)\mu_{1} - P_{0}(t) \lambda_0\\ \frac{d}{dt}P_n(t) = P_{n-1}(t) \lambda_{n-1} + P_{n+1}(t)\mu_{n+1} - P_{n}(t) (\lambda_n+\mu_n), \qquad n\geq 1. \]

Suppose the equilibrium probability exists, i.e. \(\underset{t \rightarrow \infty} \lim P_n(t)\) is a constant, then \(\underset{t \rightarrow \infty} \lim \frac{dP_n(t)}{dt}=0\). We can obtain the following results

\[ P_{1}(t)\mu_{1} = P_{0}(t) \lambda_0\\ P_{n-1}(t) \lambda_{n-1} + P_{n+1}(t)\mu_{n+1} = P_{n}(t) (\lambda_n+\mu_n), \qquad n \geq 1. \] which will result in the equilibrium probability as \[ P_0 = \frac{1}{1+\sum_{n=1}^{\infty}\frac{\lambda_0\lambda_1\cdots\lambda_{n-1}}{\mu_1\mu_2\cdots\mu_n}},\qquad \qquad \\ \\ P_n = \frac{\frac{\lambda_0\lambda_1\cdots\lambda_{n-1}}{\mu_1\mu_2\cdots\mu_n}}{1+\sum_{n=1}^{\infty}\frac{\lambda_0\lambda_1\cdots\lambda_{n-1}}{\mu_1\mu_2\cdots\mu_n}}, \qquad n\geq1\\ \] under condition \(\sum_{n=1}^{\infty}\frac{\lambda_0\lambda_1\cdots\lambda_{n-1}}{\mu_1\mu_2\cdots\mu_n} < \infty\).

Birth-death queuing model

A birth and death queuing model is an exponential queuing system model in which the arrival rates and departure rates only depend on the number of customers in the system. It’s the same as birth-death process if we consider the number of customers as state variable; the arrival rate as birth rate; and departure rate as death rate. Using the previous notation, when there are \(n\) customers in the queue (i.e. at state \(n\)), the arrival rate is denoted as \(\lambda_n\) and the departure rate is denoted as \(\mu_n\).

M/M/1 queuing model

M/M/1 model is the most common birth-death queuing model. In M/M/1 model, we assume there is only one server in the system. Customers arrive the system at rate \(\lambda\) according to a Poisson process, and the service times have an exponential distribution with rate \(\mu\). Then, the transition matrix can be expressed as

\[ \begin{pmatrix} -\lambda & \lambda & 0 & 0 & \cdots \\ \mu & -(\lambda+\mu) & \lambda & 0 & \cdots \\ 0 & \mu & -(\lambda+\mu) & \lambda & \cdots \\ 0 & 0 & \mu & -(\lambda + \mu) & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{pmatrix} , \]

which is equivalent to birth-death process with constant birth rate \(\lambda\) and death rate \(\mu\).

Using the previous equilibrium solution, we would gain

\[ P_0 = \frac{1}{1+\sum_{n=1}^{\infty}(\frac{\lambda}{\mu})^n} = 1-\frac{\lambda}{\mu}=1-\rho\\ P_n = \frac{ (\frac{\lambda}{\mu})^n }{1+\sum_{n=1}^{\infty}(\frac{\lambda}{\mu})^n} = (\frac{\lambda}{\mu})^n (1-\frac{\lambda}{\mu})=\rho^n(1-\rho), \]

where \(\rho= \frac{\lambda}{\mu} < 1\), otherwise the summation terms will not converge and we will not gain a stationary distribution. Then the average number of customers in the system, and in the queue will be:

\[ L = \sum_{n=0}^{\infty} nP_n = \sum_{n=0}^{\infty} n \rho^n (1-\rho) = \frac{\rho}{1-\rho} = \frac{\lambda}{\mu-\lambda}\\ L_q = \sum_{n=1}^{\infty} (n-1)P_n = \sum_{n=1}^{\infty} (n-1) \rho^n (1-\rho) = \frac{\rho^2}{1-\rho} = \frac{\lambda^2}{\mu(\mu-\lambda)}.\\ \]

Using Little’s Law, the average waiting time in the system and in queue will be: \[ W = \frac{L}{\lambda} = \frac{1}{\mu-\lambda} \\ W_q = \frac{L_q}{\lambda} = \frac{\lambda}{\mu(\mu-\lambda)}. \]

Simulation of M/M/1

Let \(\lambda_n\) denote the birth rate and \(\mu_n\) denote the death rate at state \(n\). Then the time before next event would follow an exponential distribution with parameter \(\lambda_n + \mu_n\). The probability that the next event is a birth is \(\frac{\lambda_n}{\lambda_n + \mu_n}\), and the probability that the next event is a death is \(\frac{\mu_n}{\lambda_n + \mu_n}\). Using these, we can write a function MM1_sim to simulate a M/M/1 queuing process. Notice that only births can occur if the population size is 0.

MM1_sim <- function(lambda, mu, 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
    # Either birth or death can occur if the population size is not 0
    if (rev(n)[1] > 0){
      Z[i] <- rexp(1, rate=lambda+mu)
      Z_type[i] <- sample(c(-1,1), size=1, prob=c(mu/(lambda+mu), lambda/(lambda+mu))) 
      n[i] <- n[i-1] + Z_type[i]
    }
    # only birth can occur if the population size is 0
    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))

}

Example: The relationship between \(\mu\) and \(\lambda\)

Here let’s look at a few simulated examples using several combinations of \(\mu\) and \(\lambda\). We first set (mu, lambda) =(2.5, 3), (3,3), (3.5,4), (4,4), which satisfies \(\mu \leq \lambda\). The simulated \(N(t)\) are as follows:

set.seed(123)
par(mfrow=c(2,2))
lambda0 = c(3, 3, 4, 4)
mu0 = c(2.5, 3, 3.5, 4)
for (i in 1:4){
  res <- MM1_sim(lambda=lambda0[i], mu=mu0[i], T_end=500)
  plot(cumsum(res$Z), res$n, type='l', xlab='time', ylab='population', main=paste0('lambda=', lambda0[i], ',', 'mu=', mu0[i]))
}

From the plots, we can see that the population explodes with time passing for either \(\mu = \lambda\) or \(\mu < \lambda\). Try to change the random seed and implement the above code by yourself. Sometimes you may observe the population doesn’t seem to explode for \(\lambda = \mu\). Then, you may want to extend the time length and see what happens.

Next, we set \((mu, lambda) = (1, 1.5), (2, 2.5), (3,3.5), (4,4.5)\) which satisfy \(\lambda < \mu\):

set.seed(123)
par(mfrow=c(2,2))
lambda0 = c(1, 2, 3, 4)
mu0 = lambda0+0.5
for (i in 1:4){
  res <- MM1_sim(lambda=lambda0[i], mu=mu0[i], T_end=500)
  plot(cumsum(res$Z), res$n, type='l', xlab='time', ylab='population', main=paste0('lambda', '=', lambda0[i], ',', 'mu', '=', mu0[i]))
}

From the plots, we can see the population is quiet stationary. Even though some quasi-explosion seems to exist in the first and second graph, the explosion is not consistent and the population comes down quickly as time passing by. Next, you may want to testify whether the distribution of the stationary case is consistent with the previous theoretical values.

Example: The equilibrium solutions

Here, we simulate the M/M/1 process and compute the simulated average stay time of each state. Then we compare the simulated distribution of stay time for each state with the theoretical values. Here we set (lambda, mu) = (1,4) and T_end = 500. We write a function p_mm1 to compute the theoretical distribution and sim_mm1 to compute the simulated values.

lambda=1
mu=4
rho=lambda/mu
# Write a function to compute the theoretical values
p_mm1 <- function(n,rho){
  return(rho^n*(1-rho))
}
# Using the previous simulated result to compute the simulated stay time at each state
sim_mm1 <- function(state, res){
  res$Z <- c(res$Z, 0)
  if (state==0){return(sum(res$Z[1], res$Z[which(res$n==0)+1])/sum(res$Z))}
  else{return(sum(res$Z[which(res$n==state)+1])/sum(res$Z))}
}

res <- MM1_sim(lambda=lambda, mu=mu, T_end=500)
n_min <- min(res$n)
n_max <- max(res$n)
x <- seq(n_min, n_max, 1)
sim_p <- c()
the_p <- c()
# Use sim_mm1 to compute the simulated result
for (i in 1:length(x)){
  sim_p[i] <- sim_mm1(x[i], res)
}
# Use p_mm1 to compute the theoretical result
theo_p <- p_mm1(x, rho)
result <- data.frame(simulated = round(sim_p,4), theory = round(theo_p,4))
rownames(result) = x
result
##   simulated theory
## 0    0.7403 0.7500
## 1    0.2014 0.1875
## 2    0.0486 0.0469
## 3    0.0067 0.0117
## 4    0.0030 0.0029