Problem 1, 2, 3, and 5 will be done together since it calls for comparing the results between each. Problem 12 will be done separately.

Note: we will repeat problem 1 here so the results are available for comparison.

Problem 1, 2, 3, & 5

For this serie of problems, we are called to determine \({ W }_{ q }\), W, \({ L }_{ q }\), L, and \(\rho\) for various configurations of queue system where;

\({ W }_{ q }\) denotes the steady-state average time in the queue,
W denotes the steady-state average time in the system,
\({ L }_{ q }\) denotes the steady-state average number of entities in the queue,
L denotes the steady-state average number of entities in the system, and
\(\rho\) denotes the steady-state utilization of the server

Finally, \(\lambda\) will denote the arrival rate and \(\mu\) the service rate in all the scenarios.
\(\mu=\frac { 1 }{ E(S) }\), where S is an RV representing an entity’s service time in a individual server.

For the subsequent calculations, we will make use of the little law and other relationships, mainly:

\(L=\lambda\cdot W\)
\({ L }_{ q }=\lambda\cdot { W }_{ q }\)
\(W={ W }_{ q }+E(S)\)

\(\rho =\frac { \lambda }{ c\cdot \mu }\) where c denotes the number of servers

M/M/1 Queue (problem# 1)

The first case we will consider is for an M/M/1 queue with mean interarrival time 1.25 minutes and mean service time 1 minute.
In this case we have;
\(\lambda=\frac { 1 }{ 1.25 }=0.8\)
\(\mu=\frac { 1 }{ 1.0 } =1.0\)

and L is given by:
\({ L }_{ q }=\frac { { \rho }^{ 2 } }{ 1-\rho }\)

M/G/1 Queue (problem# 2)

In this case, we also have
\(\lambda=\frac { 1 }{ 1.25 }=0.8\)
\(\mu=\frac { 1 }{ 1.0 } =1.0\)

However, the service time are uniformely distributed between a=01 and b=1.9. The expected value of this distribution is (a+b)/2=1 and the standard deviation is given by:
\(\sigma =\sqrt { \frac { { (b-a) }^{ 2 } }{ 12 } }\)

and \({ L }_{ q }\) is given by:
\({ L }_{ q }=\frac { { \rho }^{ 2 }(1\quad +\quad { \sigma }^{ 2 }{ \mu }^{ 2 }) }{ 2(1-\rho ) }\)

M/G/1 Queue (problem# 3)

In this case, we also have
\(\lambda=\frac { 1 }{ 1.25 }=0.8\)
\(\mu=\frac { 1 }{ 1.0 } =1.0\)

However, the service time are triangularly distributed between a=01 and b=1.9 and with a node at m=1.0. The expected value of this distribution is (a+m+b)/2=3 and the standard deviation is given by:
\(\sigma =\sqrt { \frac { { a }^{ 2 }+m^{ 2 }+{ b }^{ 2 }-am-ab-bm }{ 18 } }\)

and \({ L }_{ q }\) is given by:
\({ L }_{ q }=\frac { { \rho }^{ 2 }(1\quad +\quad { \sigma }^{ 2 }{ \mu }^{ 2 }) }{ 2(1-\rho ) }\)

M/M/3 Queue (problem# 5)

In this case we have a M/M/3 queue with mean interarrival time 1.25 minutes and mean service time 3 minute for each server.
In this case we have;
\(\lambda=\frac { 1 }{ 1.25 }=0.8\)
\(\mu=\frac { 1 }{ 3.0 } =0.3333333\)

and { L }_{ q } is given by:
\({ L }_{ q }=\frac { \rho { (c\rho ) }^{ c }\cdot p(0) }{ c!\cdot { (1-\rho ) }^{ 2 } }\),
where p(0) denotes the following:
\(p(0)={ \left[ \frac { { (c\rho ) }^{ c } }{ c!(1-\rho ) } +\sum _{ n=0 }^{ c-1 }{ \frac { { (c\rho ) }^{ n } }{ n! } } \right] }^{ -1 }\)

Calculations for all scenarios

We will use R to calculate \({ W }_{ q }\), W, \({ L }_{ q }\), L, and \(\rho\). The results of each scenarios will be stored in a table to faciliate comparison of the results.
(Where ever possible, code will be re-used).

## Validation of inputs to any functions will not be done.  We expect input to be valide for calculations

#Function to calculate server utilization 
f_rho <- function(c, lambda, mu) {
  rho <- lambda/(c*mu)
  return(rho)
}

#Function to calculate P(0)
f_p_0 <- function(c, rho) {
  #We will calculate first term separately from sum
  p_0_1<- ((c*rho)^c)/(factorial(c)*(1-rho))
  
  # Calculate sum with for loop
  p_0_2 <- 0
  for (i in 1:c){
    p_0_i <- ((c*rho)^(i-1))/factorial(i-1)
    p_0_2 <- p_0_2 + p_0_i
  }
  
  # Calculate p_0 
  p_0 <- (p_0_1 + p_0_2)^(-1)
  
  return(p_0)
}

#Functions to Calculate Lq for M/M/c queue system
f_Lq_M <- function(c, rho, p_0){
  
  Lq <- ((rho*(c*rho)^c*p_0))/(factorial(c)*(1-rho)^2)
  
  return(Lq)
}

#Function to Calcualte Exptected Value and Standard Deviation for Uniform and triangular distribution

f_uniform_esd <- function(a, b){
  u_exp_v <- (a+b)/2
  u_sd <- sqrt((b-a)^2/12)
  u_result <- list("exp_v"=u_exp_v, "sd"=u_sd)
  return(u_result)
}

f_triangular_esd <- function(a, m, b){
  u_exp_v <- (a+m+b)/3
  u_sd <- sqrt((a^2+m^2+b^2 -a*m -a*b -m*b)/18)
  u_result <- list("exp_v"=u_exp_v, "sd"=u_sd)
  return(u_result)
}

# Function to Calculate Lq for M/G/1 queue system
f_Lq_G <- function(rho, exp_v, sd_){
  
  Lq <- ((rho^2)*(1+(sd_^2*exp_v^2))/(2*(1-rho)))
  
  return(Lq)
}

#Function to calculate other values; W, Wq, L from Lq and Lambda

f_little_law <- function (lambda, exp_v, lq){
  # Apply Little Law
  wq <- lq/lambda
  w <- wq+exp_v
  l <- w*lambda
  
  return(list("Wq"=wq, "W"=w, "L"=l))
}

f_mmc <- function(c, avg_arvl, avg_srv){
  
  #Function to determine Rho, Lq, Wq, L, W for MMC queue
  #Function will take as input parameters
  # c=number of server
  # avg_arvl = Average arrival time
  # avg_srv = Average service time at each server
  
  # This function will return Rho, Wq, W, Lq, L as 1 row data frame
  
  # Calculation lambda and mu
  r_lambda <- 1/avg_arvl
  r_mu <- 1/avg_srv
  
  # Using previously defined functions, calculating rho, p(0), and Lq
  r_rho <- f_rho(c,r_lambda, r_mu)
  r_p_0 <- f_p_0(c, r_rho)
  r_lq <- f_Lq_M(c, r_rho, r_p_0)
  
  # Using Little Law
  r_list <- f_little_law(r_lambda, avg_srv, r_lq)
  
  r_result <- data.frame(r_rho, r_list[[1]], r_list[[2]], r_lq, r_list[[3]])
  
  return(r_result)
}

f_mg1 <- function (c, avg_arvl, avg_srv, sd){
  
  #Function to determine Rho, Lq, Wq, L, W for MG1 queue
  #Function will take as input parameters
  # c=number of server (for consistency)
  # avg_arvl = Average arrival time
  # avg_srv = Average service time at each server
  
  # This function will return Rho, Wq, W, Lq, L as 1 row data frame
  
  # Calculation lambda and mu
  r_lambda <- 1/avg_arvl
  r_mu <- 1/avg_srv
  
  # Using previously defined functions, calculating rho, p(0), and Lq
  r_rho <- f_rho(c,r_lambda, r_mu)
  
  # Calculation of Lq
  r_lq <- f_Lq_G(r_rho, avg_srv, sd)
  
  # Using Little Law
  r_list <- f_little_law(r_lambda, avg_srv, r_lq)
  
  r_result <- data.frame(r_rho, r_list[[1]], r_list[[2]], r_lq, r_list[[3]])
  
  return(r_result)
}

#-----------------------------------------------#
# Building data frames will all results

# M/M/1 queue
d_mm1 <- f_mmc(1, 1.25, 1)

# Uniform distribution
sd_u <- f_uniform_esd(0.1, 1.9)[[2]]
d_mg1_u <- f_mg1(1, 1.25, 1, sd_u)

# Triangular distribution
t_list <- f_triangular_esd(0.1, 1, 1.9)
exp_v_t <- t_list[[1]]
sd_t <- t_list[[2]]

d_mg1_t <- f_mg1(1, 1.25, exp_v_t, sd_t)

# M/M/3 queue
d_mm3 <- f_mmc(3, 1.25, 3)

# Building a result data frame

d_r <- do.call("rbind", list(d_mm1, d_mg1_u, d_mg1_t, d_mm3))


d_queue <- c("mm1", "mg1_u", "mg1_t", "mm3")

d_result <- cbind(d_queue,d_r)

colnames(d_result) <- c("Model", "rho", "Wq", "W", "Lq", "L")

knitr::kable(d_result)
Model rho Wq W Lq L
mm1 0.8 4.000000 5.000000 3.200000 4.000000
mg1_u 0.8 2.540000 3.540000 2.032000 2.832000
mg1_t 0.8 2.270000 3.270000 1.816000 2.616000
mm3 0.8 3.235955 6.235955 2.588764 4.988764

Simulations for all scenarios

m/m/1

m/m/1

m/g(uniform)/1

m/g(uniform)/1

m/g(triangular)/1

m/g(triangular)/1

m/m/3

m/m/3

Problem 12

Given Model below, we are looking to compute the traffic intensity at each station, denoted \({ \rho }_{ Station }\).

Model of Clinic

Model of Clinic

\({ \rho }_{ Station }=\frac { { \lambda }_{ station } }{ { c }_{ Station }\cdot { \mu }_{ Station } }\)

The average arrival inter-arrival time into the system is 6 minutes, this will conrespond to a arrival rate \(\lambda =\frac { 1 }{ 6 }=0.1666667\)

For each station, we will determine the arrival rate.
For simplicity, arrival rate for Registration = \({ \lambda }_{ R }\), for Trauma Rooms = \({ \lambda }_{ T }\), for Exam Room = \({ \lambda }_{ E }\), for Treatment Rooms = \({ \lambda }_{ Tr }\).

Based on the branching probabilities, we can compute each arrival rate as follows:
\({ \lambda }_{ R }=0.9\times \lambda\)
\({ \lambda }_{ T }=0.1\times \lambda\)
\({ \lambda }_{ E }=1\times {\lambda }_{R}\) \({ \lambda }_{ Tr }={ 0.6\times \lambda }_{ E }\quad +\quad 1\times { \lambda }_{ T }\)

stations <- c("Sign-in", "Registration", "Trauma Rooms", "Exam Rooms", "Treatment Rooms")

s_avg_arrival <- 6
s_lambda <- 1/s_avg_arrival

all_lambda <- c(s_lambda, 0.9*s_lambda, 0.1*s_lambda, 0.9*s_lambda, (0.1*s_lambda)+(0.6*0.9*s_lambda))

all_Server_capacity <- c(2, 1, 2, 3, 2)

all_server_svc_avg <- c(3, 5, 90, 16, 15)

all_server_svc_rte <- 1/all_server_svc_avg

all_rho <- all_lambda/(all_Server_capacity*all_server_svc_rte)

d<-data.frame(all_rho)
dt <-t(d)
colnames(dt)<- stations

knitr::kable(dt)
Sign-in Registration Trauma Rooms Exam Rooms Treatment Rooms
all_rho 0.25 0.75 0.75 0.8 0.8

This clinic will “work” based on each server utilization percentage, we do not have a situation where the server would not be able to handle the load. This would be indicated by a \({ \rho }_{ Station }\) > 1.

If we could add a server, we would add it to the Exam Rooms, this is one of the station where the utilization is the highest and where arrival rate is one of the highest.