Data 604: Homework 2
Judd Anderman
September 6, 2017
Kelton et al. (2014). Simio and Simulation: Modeling, Analysis, Applications (3rd Ed.).
Section 2.6, Problems 2, 3, 5 & 12
2.6.1) For an \(M/M/1\) queue with mean interarrival time 1.25 minutes and mean service time 1 minute, find all five of \(W_q\), \(W\), \(L_q\), \(L\), and \(\rho\). For each, interpret in words. Be sure to state all your units (always!), and the relevant time frame of operation.
\(\lambda = 1 / 1.25 = 0.8 \textrm{ entities/minute}\)
\(\mu = 1/E(S) = 1 \textrm{ entity/minute}\)
\(\rho = \lambda/\mu = 0.8 / 1 = 0.8\)
\(L = \rho/(1 - \rho) = 0.8/(1 - 0.8) = 4 \textrm{ entities}\)
\(L = \lambda W\)
\(W = L / \lambda = 4 / 0.8 = 5 \textrm{ minutes}\)
\(W = W_q + E(S)\)
\(5 = W_q + 1\)
\(W_q = 4 \textrm{ minutes}\)
\(L_q = \lambda W_q = 0.8 \times 4 = 3.2 \textrm{ entities}\)
\(W_q\), the steady-state or long-run average time in queue is 4 minutes. \(W\), the steady-state average time in system, or the sum of the expected time in queue and the expected service time, is 5 minutes. \(L_q\), the steady-state average number of entities in queue, or expected number of entities waiting to be served, is 3.2. The steady-state average number of entities in the system, \(L\), is 5. \(\rho\), the steady-state utilization of the single server in the system, or the average rate of utilization of the server over an infinte time horizon, is 0.8.
2.6.2) Repeat problem 1, except assume that the service times are not exponentially distributed, but rather (continuously) uniformly distributed between \(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. Compare all five of your numerical results to those from Problem 1 and explain intuitively with respect to this change in the service-time distribution (but with its expected value remaining at 1). Hint: In case, you’ve forgotten or your calculus has rusted completely shut, or you haven’t already found it with a web search, the standard deviation of the continuous uniform distribution between \(a\) and \(b\) is \(\sqrt{(b-a)^2/12}\) (that’s right, you always divide by 12 regardless of what \(a\) and \(b\) are … the calculus just works out that way).
\(\lambda = 1 / 1.25 = 0.8 \textrm{ entities/minute}\)
\(\mu = 1/E(S) = 1 \textrm{ entity/minute}\)
\(\rho = \lambda/\mu = 0.8 / 1 = 0.8\)
\(\theta = \sqrt{(b-a)^2/12} = \sqrt{(1.9 - 0.1)^2/12} = \sqrt{1.8^2 / 12} =\) \(\sqrt{0.27} \textrm{ minutes}\)
In an \(M/G/1\) queueing system, \(W_q = \frac{\lambda(\theta^2 + 1/\mu^2)}{2(1 - \lambda/\mu)}=\) \(\frac{0.8(0.27 + 1)}{2(1 - 0.8)} = \frac{1.016}{0.4} = 2.54 \textrm{ minutes}\)
\(L_q = \lambda W_q = 0.8 \times 2.54 = 2.032 \textrm{ entities}\)
\(W = W_q + E(S) = 2.54 + 1 = 3.54 \textrm{ minutes}\)
\(L = \lambda W = 0.8 \times 3.54 = 2.832 \textrm{ entities}\)
See below for comparison of steady-state average conditions in this system with those in single-server queue systems with exponentially distributed and triangularly distributed service times with equal expected values for service time and the inter-arrival time distribution in all three systems.
2.6.3) Repeat problem 1, except assume that the service times are triangularly distributed between \(a = 0.1\) and \(b = 1.9\), and with the mode at \(m = 1.0\). Compare all five of your results to those from Problems 1 and 2. Hint: The expected value of a triangular distribution between \(a\) and \(b\), and with mode \(m\) \((a < m < b)\), is \((a + m + b)/3\), and the standard deviation is \(\sqrt{(a^2 + m^2 + b^2 - am - ab - bm)/18}\) … do you think maybe it’s time to dust off that calculus book (or, at least hone your web-search skills)?
\(\lambda = 1 / 1.25 = 0.8 \textrm{ entities/minute}\)
\(E(S) = (a + m + b)/3 = (0.1 + 1 + 1.9)/3 = 1 \textrm{ minute}\)
\(\mu = 1/E(S) = 1 \textrm{ entity/minute}\)
\(\rho = \lambda/\mu = 0.8 / 1 = 0.8\)
\(\theta = \sqrt{(a^2 + m^2 + b^2 - am - ab - bm)/18} =\) \(\sqrt{(0.1^2 + 1^2 + 1.9^2 - 0.1 - 0.1(1.9) - 1.9)/18} =\) \(\sqrt{2.43/18} = \sqrt{0.135} \textrm{ minutes}\)
\(W_q = \frac{\lambda(\theta^2 + 1/\mu^2)}{2(1 - \lambda/\mu)}=\) \(\frac{0.8(0.135 + 1)}{2(1 - 0.8)} = \frac{0.908}{0.4} = 2.27 \textrm{ minutes}\)
\(L_q = \lambda W_q = 0.8 \times 2.27 = 1.816 \textrm{ entities}\)
\(W = W_q + E(S) = 2.27 + 1 = 3.27 \textrm{ minutes}\)
\(L = \lambda W = 0.8 \times 3.27 = 2.616 \textrm{ entities}\)
MG1 <- data.frame(Service_Time_Distribution = c("exponential", "uniform", "triangular"),
W_q = numeric(3),
W = numeric(3),
L_q = numeric(3),
L = numeric(3),
p = numeric(3)
)
lambda <- 1 / 1.25
mu <- 1
MG1$p <- lambda / mu
MG1[1, ]$L <- MG1[1, ]$p / (1 - MG1[1, ]$p)
MG1[1, ]$W <- MG1[1, ]$L / lambda
MG1[1, ]$W_q <- MG1[1, ]$W - (1 / mu)
MG1[1, ]$L_q <- MG1[1, ]$W_q * lambda
for (i in 2:3) {
a <- 0.1
b <- 1.9
m <- 1
if (MG1[i, ]$Service_Time_Distribution == "uniform") {
theta <- sqrt((b - a)^2 / 12)
} else {
theta <- sqrt((a^2 + m^2 + b^2 - (a * m) - (a * b) - (b * m)) / 18)
}
MG1[i, ]$W_q = (lambda * (theta^2 + 1 / mu^2)) / (2 * (1 - lambda/ mu))
MG1[i, ]$L_q = lambda * MG1[i, ]$W_q
MG1[i, ]$W = MG1[i, ]$W_q + (1 / mu)
MG1[i, ]$L = lambda * MG1[i, ]$W
}
knitr::kable(MG1, caption = "M/G/1 queue systems comparison")| Service_Time_Distribution | W_q | W | L_q | L | p |
|---|---|---|---|---|---|
| exponential | 4.00 | 5.00 | 3.200 | 4.000 | 0.8 |
| uniform | 2.54 | 3.54 | 2.032 | 2.832 | 0.8 |
| triangular | 2.27 | 3.27 | 1.816 | 2.616 | 0.8 |
The queue system with the triangular service time distribution has the shortest steady-state average time in queue, \(W_q\), and shortest steady-state average time in the system, \(W\), while the system with exponentially distributed service time has the largest values for those two metrics of the three systems under comparison. Similarly, the system with triangularly distributed service time has the smallest steady-state average number of entities in queue, \(L_q\), and the smallest steady-state average number of entities in the system, \(L\), while the system analyzed in Problem 1 with exponentially distributed service time has the largest values for \(L_q\) and \(L\). The three systems have equal values of \(\rho\), the average utilization of the single server.
The proability density functions and cumulative distribution functions for each of the service time distributions are plotted below. It makes intuitive sense that the queueing system with exponentially distributed service time would have the longest average wait times and average number of entities in queue given that the likelihood of a service time greater than the average interarrival time is much greater in that system than in the two others, and a longer queue and longer wait times are more likely to develop in those conditions. It also seems reasonable that the system with triangularly distributed service times would have the best metrics for steady-state average wait time and queued entities since service times would be more tightly clustered about the expected value of 1 minute than in the system with continuous uniformly distributed service times. This relative reduction in the variability of service times should improve peformance with regard to the timely flow of entities through the system.
library(ggplot2)
x <- seq(0, 5, 0.01)
dtri <- function(x, a, m, b) {
out <- numeric(length(x))
for (i in 1:length(x)) {
if (x[i] >= a & x[i] <= m) {
out[i] <- 2 * (x[i] - a) / ((b - a) * (m - a))
} else if (x[i] > m & x[i] <= b) {
out[i] <- 2 * (b - x[i]) / ((b - a) * (b - m))
} else {
out[i] <- 0
}
}
return(out)
}
ptri <- function(x, a, m, b) {
out <- numeric(length(x))
for (i in 1:length(x)) {
if (x[i] >= a & x[i] <= m) {
out[i] <- (x[i] - a)^2 / ((b - a) * (m - a))
} else if (x[i] > m & x[i] <= b) {
out[i] <- 1 - (b - x[i])^2 / ((b - a) * (b - m))
} else if (x[i] < a) {
out[i] <- 0
} else {
out[i] <- 1
}
}
return(out)
}
ggplot() +
geom_line(aes(x, dexp(x, 1))) +
geom_line(aes(x, dunif(x, 0.1, 1.9)), color = "red") +
geom_line(aes(x, dtri(x, 0.1, 1, 1.9)), color = "blue") +
labs(x = "t", y = "p(t)")
ggplot() +
geom_line(aes(x, pexp(x, 1))) +
geom_line(aes(x, punif(x, 0.1, 1.9)), color = "red") +
geom_line(aes(x, ptri(x, 0.1, 1, 1.9)), color = "blue") +
labs(x = "t", y = "p(X <= t)")
The simulation results shown below match closely with the steady-state average conditions in each system calculated above.
Simio Simulation Results (1 56-day long iteration for each service time distribution):
Problem 6.2.1: exponentially distributed service times
Problem 6.2.2: uniformly distributed service times
Problem 6.2.3: triangurlarly distributed service times
2.6.5) Repeat Problem 1, except for an \(M/M/3\) queue with mean interarrival time 1.25 minutes and mean service time 3 minutes at each of the three servers. Hint: You might want to consider creating a computer program or spreadsheet, or use mmc.exe.
lambda <- 1 / 1.25
mu <- 1 / 3
mmc <- function(lambda, mu, c) {
p <- lambda / (c * mu)
f <- function (c, p) {
x <- 0
for (n in 0:(c - 1)) {
x <- x + ((c * p)^n / factorial(n))
}
return(x)
}
p_0 <- 1 / ((c * p)^c / (factorial(c) * (1 - p)) + f(c, p))
L_q <- (p * (c * p)^c * p_0) / (factorial(c) * (1 - p)^2)
W_q <- L_q / lambda
L <- L_q + (lambda / mu)
W <- L / lambda
return(data.frame(W_q, W, L_q, L, p))
}
MM3 <- mmc(lambda, mu, 3)
knitr::kable(MM3[, 1:5], caption = "M/M/3 queue system")| W_q | W | L_q | L | p |
|---|---|---|---|---|
| 3.235955 | 6.235955 | 2.588764 | 4.988764 | 0.8 |
2.6.12) In the urgent-care clnic of Figure 2.1, suppose that the patients arrive from outside into the clinic (coming from 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 as shown in Figure 2.1. 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. Will this clinic “work”, i.e., be able to handle the external patient load? Why or why not? If you could add a single server to the system, and add it any of the five stations, where would you add it? Why? Hint: Unless you like using your calculator, a spreadsheet or computer program might be good, or perhaps use mmc.exe.
As shown below, the “local” traffic intensity, or steady-state utilization \(rho\), at each clinic station is less than 1 and as such, the clinic will be able to handle the specified external patient load. On the basis of that metric alone, one might be inclined to add an additional server to either the Exam or Treatment stations since those have the highest average server utilization rates in the clinic. Moreover, since 90% of patient arrivals will be seen in an Exam room, that seems to be a logical station to add a server to in order to improve patient flow in the clinic. However, the addition of a single Trauma room would dramatically reduce average waiting time at that station and would mean that the most acute patients arriving at the clinic would be triaged and treated much more rapidly on average. As such, one would expect improved health outcomes in that group of patients with the highest acuity, especially in terms of avoiding significant injury or worse in circumstances where time to treat is a critical factor.
sign_in_lambda <- 1 / 6
clinic <- cbind(node = c("sign_in", "registration", "trauma", "exam", "treatment"),
rbind(mmc(sign_in_lambda, (1 / 3), 2),
mmc(0.9 * sign_in_lambda, (1 / 5), 1),
mmc(0.1 * sign_in_lambda, (1 / 90), 2),
mmc(0.9 * sign_in_lambda, (1 / 16), 3),
mmc(0.64 * sign_in_lambda, (1 / 15), 2)
)
)
knitr::kable(clinic[, c("node", "p")], caption = 'Urgent-Care Clinic "Local" Traffic Intensities')| node | p |
|---|---|
| sign_in | 0.25 |
| registration | 0.75 |
| trauma | 0.75 |
| exam | 0.80 |
| treatment | 0.80 |
knitr::kable(clinic, caption = "Steady-State Conditions at Urgent-Care Clinic Stations")| node | W_q | W | L_q | L | p |
|---|---|---|---|---|---|
| sign_in | 0.20000 | 3.20000 | 0.0333333 | 0.5333333 | 0.25 |
| registration | 15.00000 | 20.00000 | 2.2500000 | 3.0000000 | 0.75 |
| trauma | 115.71429 | 205.71429 | 1.9285714 | 3.4285714 | 0.75 |
| exam | 17.25843 | 33.25843 | 2.5887640 | 4.9887640 | 0.80 |
| treatment | 26.66667 | 41.66667 | 2.8444444 | 4.4444444 | 0.80 |
exam_clinic <- cbind(node = c("sign_in", "registration", "trauma", "exam", "treatment"),
rbind(mmc(sign_in_lambda, (1 /3), 2),
mmc(0.9 * sign_in_lambda, (1 / 5), 1),
mmc(0.1 * sign_in_lambda, (1 / 90), 2),
mmc(0.9 * sign_in_lambda, (1 / 16), 4),
mmc(0.64 * sign_in_lambda, (1 / 15), 2)
)
)
knitr::kable(exam_clinic, caption =
paste("Steady-State Conditions at Urgent-Care Clinic Stations",
"with Fourth Exam Room Added",
sep = "\n"
)
)| node | W_q | W | L_q | L | p |
|---|---|---|---|---|---|
| sign_in | 0.200000 | 3.20000 | 0.0333333 | 0.5333333 | 0.25 |
| registration | 15.000000 | 20.00000 | 2.2500000 | 3.0000000 | 0.75 |
| trauma | 115.714286 | 205.71429 | 1.9285714 | 3.4285714 | 0.75 |
| exam | 2.870432 | 18.87043 | 0.4305648 | 2.8305648 | 0.60 |
| treatment | 26.666667 | 41.66667 | 2.8444444 | 4.4444444 | 0.80 |
treatment_clinic <- cbind(node = c("sign_in", "registration", "trauma", "exam", "treatment"),
rbind(mmc(sign_in_lambda, (1 /3), 2),
mmc(0.9 * sign_in_lambda, (1 / 5), 1),
mmc(0.1 * sign_in_lambda, (1 / 90), 2),
mmc(0.9 * sign_in_lambda, (1 / 16), 3),
mmc(0.64 * sign_in_lambda, (1 / 15), 3)
)
)
knitr::kable(treatment_clinic, caption =
paste("Steady-State Conditions at Urgent-Care Clinic Stations",
"with Third Treatment Room Added",
sep = "\n"
)
)| node | W_q | W | L_q | L | p |
|---|---|---|---|---|---|
| sign_in | 0.200000 | 3.20000 | 0.0333333 | 0.5333333 | 0.2500000 |
| registration | 15.000000 | 20.00000 | 2.2500000 | 3.0000000 | 0.7500000 |
| trauma | 115.714286 | 205.71429 | 1.9285714 | 3.4285714 | 0.7500000 |
| exam | 17.258427 | 33.25843 | 2.5887640 | 4.9887640 | 0.8000000 |
| treatment | 2.933537 | 17.93354 | 0.3129106 | 1.9129106 | 0.5333333 |
trauma_clinic <- clinic <- cbind(node = c("sign_in", "registration", "trauma", "exam", "treatment"),
rbind(mmc(sign_in_lambda, (1 /3), 2),
mmc(0.9 * sign_in_lambda, (1 / 5), 1),
mmc(0.1 * sign_in_lambda, (1 / 90), 3),
mmc(0.9 * sign_in_lambda, (1 / 16), 3),
mmc(0.64 * sign_in_lambda, (1 / 15), 2)
)
)
knitr::kable(trauma_clinic, caption =
paste("Steady-State Conditions at Urgent-Care Clinic Stations",
"with Third Trauma Room Added",
sep = "\n"
)
)| node | W_q | W | L_q | L | p |
|---|---|---|---|---|---|
| sign_in | 0.20000 | 3.20000 | 0.0333333 | 0.5333333 | 0.25 |
| registration | 15.00000 | 20.00000 | 2.2500000 | 3.0000000 | 0.75 |
| trauma | 14.21053 | 104.21053 | 0.2368421 | 1.7368421 | 0.50 |
| exam | 17.25843 | 33.25843 | 2.5887640 | 4.9887640 | 0.80 |
| treatment | 26.66667 | 41.66667 | 2.8444444 | 4.4444444 | 0.80 |