Assignment 2 - Queueing Theory

Output Metric Definitions

\(W_q\) = the steady-state average time in queue (excluding service times) of entities.

\(W\) = the steady-state average time in system (including service times) of entities.

\(L_q\) = the steady-state average number of entities in queue. This is a time average.

\(L\) = the steady-state average number of entities in the system. This is a time average.

\(\rho\) = the steady-state utilization of a server or group of parallel identical servers.

Problem 2.6.2

Repeat Problem 1, except assume that the service times are not exponentially distributed, but rather (continuously) uniformly distributed betweeen a = 0.1 and b = 1.9. Note that the expected value of this uniform distribution is \((a + b)/2 = 1\), the same as the expected service time in Problem 1.

The standard deviation of the continuous uniform distribution between \(a\) and \(b\) is:

\[ \sqrt{(b - a)^2/12} \]

For this problem, we’re going to use the M/G/1 model which supports other types of service-time distributions.

\[ W_q = \frac{\lambda(\sigma^2 + 1/\mu^2)}{2(1 - \lambda/mu} \]

MG1 <- function(mean_interval_time, mean_service_time, sd) {
  
  (lambda <- 1/mean_interval_time)
  (mu <- 1/mean_service_time)
  (rho <- lambda/mu)
  
  L <- rho/(1-rho)

  Wq <-  lambda*(sd^2+1/(mu)^2)/(2*(1-lambda/mu))
  W <- Wq + mean_service_time
  
  Lq <- lambda * Wq
  
  return (data.frame(lambda = lambda, 
                     mu     = mu, 
                     W      = W, 
                     Wq     = Wq,  
                     L      = L, 
                     Lq     = Lq, 
                     rho    = rho))
}   
  
a <- 0.1
b <- 1.9
mean_service_time <- (a + b)/2
sd <- sqrt((b-a)^2/12)

metrics <- MG1(mean_interval_time = 1.25, mean_service_time = mean_service_time, sd = sd)
lambda mu W Wq L Lq rho
0.8 1 3.54 2.54 4 2.032 0.8

Compared to the metrics from problem 1:

lambda mu W Wq L Lq rho
0.8 1 5 4 4 3.2 0.8

We see that the M/G/1 model has lower average wait times in queue and in the system (\(W_q\) and \(W\)). This can be attributed to the difference in the service time distribution in the M/G/1 model which is continuously uniformly distributed.0.

Problem 2.6.3

Repeat Problem 1, except assume that the service times are triangularly distributed betweeen a = 0.1 and b = 1.9, and with mode at \(m = 1.0\). Note that the expected value of this uniform distribution is \((a + m + b)/3 = 1\).

The standard deviation of the triangular distribution between \(a\) and \(b\) is:

\[ \sqrt{(a^2 + m^2 + b^2 - am - ab -bm)/18} \]

For this problem, we’re going to again use the M/G/1 model which supports other types of service-time distributions.

\[ W_q = \frac{\lambda(\sigma^2 + 1/\mu^2)}{2(1 - \lambda/mu} \]

a <- 0.1
b <- 1.9
m <- 1

mean_service_time <- (a + m + b)/3 
sd <- sqrt((a^2 + m^2 + b^2 - a*m - a*b - b*m)/18)

metrics <- MG1(mean_interval_time = 1.25, mean_service_time = mean_service_time, sd = sd)
lambda mu W Wq L Lq rho
0.8 1 3.27 2.27 4 1.816 0.8

Problem 2.6.5

Repeat Problem 1, except for an M/M/3 queue with mean interarrival time of 1.25 minutes and mean service time 3 minutes at each of the three servers.

MMC <- function(lambda, mu, c, label = NULL) {

  #lambda <- 1/mean_interval_time
  #mu <- 1/mean_service_time
  
  #c <- num_servers
  ind <- 0:(c-1)

  rho <- lambda/( mu)

  pi_0 <-1/(sum(rho**ind/factorial(ind)) + (rho**c/factorial(c))*(1/(1-rho/c)))  # p(0)
  L <- rho +  ((rho**c/factorial(c))*pi_0)*((rho/c)/((1-rho/c)**2))              # NumInSystem
  Lq <- L - rho                                                                  # NuminQueue

  Wq <- Lq/lambda                                                                # TimeInQueue

  W <- Wq + 1 / mu                                                               # TimeInSystem
  W <- L/lambda
  
  rho <- lambda/(c* mu)
  
  return (data.frame(label = as.character(label),
                     lambda = lambda, 
                     mu     = mu, 
                     W      = W, 
                     Wq     = Wq,  
                     L      = L, 
                     Lq     = Lq, 
                     rho    = rho,
                     p0     = pi_0, stringsAsFactors = F))
}

  

metrics <- MMC(lambda = 1/1.25, mu=1/3, 3, "M/M/3")
label lambda mu W Wq L Lq rho p0
M/M/3 0.8 0.3333333 6.235955 3.235955 4.988764 2.588764 0.8 0.0561798

Problem 2.6.12

In the urgent-care clinic below, suppose that the patients arrive from outside into the clinic (coming form the upper right corner of the figure and always into the Sign In station) with interarrival times that are exponentially distributed with mean 6 minutes. The number of individual servers at each station and the branching probabilities are all shown in the figure.

The service times at each node are exponentially distributed with means (all in minutes) of 3 for Sign In, 5 for Registration, 90 for Trauma Rooms, 16 for Exam Rooms, and 15 for Treatment Rooms. For each of the five stations, compute the “local” traffic intensity \(\rho_{Station}\) there.

station <- data.frame()

## ---------------------------------------------
## Sign In : M/M/2
## ---------------------------------------------

lambda <- 1/6       # 100% probability of arrival to sign in
mu <- 1/3           # service rate = mean 3 minutes
c <- 1              # 1 server 

station <- MMC(1/6, 1/3, 2, "Sign In")

## ---------------------------------------------
## Registration : M/M/1
## ---------------------------------------------

lambda <- 1/6 * 0.9 # 90% probability of going a patient going to registration
mu <- 1/5           # service rate = mean 5 minutes
c <- 1              # 1 server 

station <- rbind(station,  MMC(lambda, mu, c, "Registration"))

## ---------------------------------------------
## Trauma Room : M/M/2
## ---------------------------------------------

lambda <- 1/6 * 0.1 # 10% probability of going a patient going to registration
mu <- 1/90

station <- rbind(station,  MMC(lambda, mu, 2, "Trauma Room"))

## ---------------------------------------------
## Exam Room: M/M/3
## ---------------------------------------------

lambda <- 1/6 * 0.9 # 90% probability of going a patient going to an exam room since 10% go to a trauma room
mu <- 1/16          # service rate = mean 16 minutes
c <- 3              # 3 servers

station <- rbind(station,  MMC(lambda, mu, c, "Exam Room"))

## ---------------------------------------------
## Treatment Rooms: M/M/2
## ---------------------------------------------

# Receives a % of patients going to the exam rooms + a % of patients going to trauma rooms

# Exam Room (1):  Sign In --> Registration --> Exam Room --> Treatment = 0.9 * 1.0 * 0.6 = 0.54
# Exam Room (2):   Sign In --> Trauma Room --> Treatment = 0.1 * 1.0 = 0.1

# Probability of going to a Treatment Room = 0.54 + 0.1 = 0.64

lambda <- 1/6 * 0.64 # 64% probability of going a patient going to a treatment room
mu <- 1/15           # service rate = mean 15 minutes
c <- 2               # 2 servers

station <- rbind(station,  MMC(lambda, mu, c, "Treatment Room"))

Local Traffic Intensity \(\rho_{Station}\):

  • Sign In: \(\lambda_{SignIn}\) = 0.1666666/ (2 * 0.3333333) = 0.25
  • Registration \(\rho_{Registration}\) = 0.9 * \(\lambda_{SignIn}\)/(\(\mu_{Registration}\)) = (0.9 * 0.1666666)/0.2 = 0.75
  • Trauma Rooms \(\rho_{Trauma}\) = 0.1 * \(\lambda_{SignIn}\)/(2 * \(\mu_{Trauma}\)) = (0.1 * 0.1666666)/(2 * 0.01111111) = 0.75
  • Exam Rooms \(\rho_{Exam}\) = 0.9 * \(\lambda_{SignIn}\)/(2 * \(\mu_{Exam}\)) = (0.9 * 0.1666666)/(3 * 0.06250000) = 0.80
  • Treatment Rooms \(\rho_{Treatment}\) = 0.64 * \(\lambda_{SignIn}\)/(2 * \(\mu_{Treatment}\)) = (0.64 * 0.1666666)/(2 * 0.06666667)=0.80

These values for \(\rho_{Station}\) are consistent with the calculated \(\rho\) values below.

Results

label lambda mu W Wq L Lq rho p0
Sign In 0.1666667 0.3333333 3.20000 0.20000 0.5333333 0.0333333 0.25 0.6000000
Registration 0.1500000 0.2000000 20.00000 15.00000 3.0000000 2.2500000 0.75 0.2500000
Trauma Room 0.0166667 0.0111111 205.71429 115.71429 3.4285714 1.9285714 0.75 0.1428571
Exam Room 0.1500000 0.0625000 33.25843 17.25843 4.9887640 2.5887640 0.80 0.0561798
Treatment Room 0.1066667 0.0666667 41.66667 26.66667 4.4444444 2.8444444 0.80 0.1111111

This clinic will work in the sense that it will be able to handle the patient load. The value of \(\rho\) is less than 1 for all stations indicating a stable steady-state utilization. Of the stations, the exam and treatment rooms have the highest utilization at 80%. If we were to add an additional server to the system, consider 1.) the trauma room with the highest average time in queue \(W_q\) of 115 minutes and a utilization of 75%; 2.) the exam room with utilization of 80% and low probability of having zero entities in queue of 5%; and 3.) the treatment room with utilization of 80% and an average time in queue of 26 minutes.

Let’s re-run the models adding a server to each of the three stations:

station2 <- data.frame()

lambda <- 1/6 * 0.1 # 10% probability of going a patient going to registration
mu <- 1/90          # 90 minute mean service time
c <- 3              # 3 servers

station2 <- rbind(station2,  MMC(lambda, mu, 3, "Trauma Room (2)"))

## ---------------------------------------------
## Exam Room: M/M/3
## ---------------------------------------------

lambda <- 1/6 * 0.9 # 90% probability of going a patient going to an exam room since 10% go to a trauma room
mu <- 1/16          # service rate = mean 16 minutes
c <- 4              # 4 servers

station2 <- rbind(station2,  MMC(lambda, mu, c, "Exam Room (2)"))

## ---------------------------------------------
## Treatment Rooms: M/M/3
## ---------------------------------------------

# Receives a % of patients going to the exam rooms + a % of patients going to trauma rooms

# Exam Room (1):  Sign In --> Registration --> Exam Room --> Treatment = 0.9 * 1.0 * 0.6 = 0.54
# Exam Room (2):   Sign In --> Trauma Room --> Treatment = 0.1 * 1.0 = 0.1

# Probability of going to a Treatment Room = 0.54 + 0.1 = 0.64

lambda <- 1/6 * 0.64 # 64% probability of going a patient going to a treatment room
mu <- 1/15           # service rate = mean 15 minutes
c <- 3               # 3 servers

station2 <- rbind(station2,  MMC(lambda, mu, c, "Treatment Room (2)"))

Looking at the results, we’ll focus on the trauma room and the treatment rooms. From a pure efficiency perspective, one could argue that adding the additional server to the treatment rooms would result in a utilization drop from 80% down to 53% and the highest station average time in queue dropping from 26 minutes to close to 3 minutes.

label lambda mu W Wq L Lq rho p0
Trauma Room 0.0166667 0.0111111 205.71429 115.714286 3.428571 1.9285714 0.7500000 0.1428571
Trauma Room (2) 0.0166667 0.0111111 104.21053 14.210526 1.736842 0.2368421 0.5000000 0.2105263
Treatment Room 0.1066667 0.0666667 41.66667 26.666667 4.444444 2.8444444 0.8000000 0.1111111
Treatment Room (2) 0.1066667 0.0666667 17.93354 2.933537 1.912911 0.3129106 0.5333333 0.1871658

However, the \(W_q\) value (average time in queue) for the trauma room seems high at almost two hours. Considering these are critical minutes for severally injured people, I would allocate the additional server to this station. The original value for \(\rho(0)\) shows that this station only has an 86% probability of having someone in it, indicating a consistent usage of the trauma rooms.

References:

http://homepages.ecs.vuw.ac.nz/~schukova/SCIE201/Lectures/Lecture9_10_final.html