Kelton et al. (2014). Simio and Simulation: Modeling, Analysis, Applications (3rd Ed.)

Problems 2.9 & 4.15

library(ggplot2)
library(dplyr)
library(tidyr)

2.9) Find all five of the steady-state queueing metrics for an \(M/D/1\) queue, where \(D\) denotes a deterministic “distribution”, i.e., the associated RV (in this case representing service times) is a constant with no variation at all (also called a degenerate distribution). State parameter conditions for your results to be valid; use the same meaning for \(\lambda\), \(\mu\), and \(\rho\) as we did in this chapter. Compare your results to those if \(D\) were replaced by a distribution with mean equal to the constant value from the original \(D\) distribution, except having at least some variation.

For \(D\), \(\theta^2 = 0\)

\(\rho = \lambda / \mu\)

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

\(L_q = \lambda W_q = \frac{\lambda^2 / \mu^2}{2(1 - \lambda / \mu)}\)

\(W = W_q + E(S) = \frac{\lambda / \mu^2}{2(1 - \lambda / \mu)} + 1 / \mu\)

\(L = \lambda W = \frac{\lambda^2 / \mu^2}{2(1 - \lambda / \mu)} + \lambda / \mu\)

For \(D\), \(\theta^2 \neq 0\)

\(\rho = \lambda / \mu\)

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

\(L_q = \lambda W_q = \frac{\lambda^2(\theta^2 + 1 / \mu^2)}{2(1 - \lambda / \mu)}\)

\(W = W_q + E(S) = \frac{\lambda(\theta^2 + 1 / \mu^2)}{2(1 - \lambda / \mu)} + 1 / \mu\)

\(L = \lambda W = \frac{\lambda^2(\theta^2 + 1 / \mu^2)}{2(1 - \lambda / \mu)} + \lambda / \mu\)

In the queue system with some variability in the distribution of service times, the steady-state average time in queue and average number of entities in queue will be greater than those for a system with a degenerate distribution of service times in direct proportion to the variance \(\theta^2\) of the service time distribution.

4.15) Build a Simio model to confirm and cross-check the steady-state queueing-theoretic results from your solutions to the \(M/D/1\) queue of Problem 9 in Chapter 2. Remember that you Simio model is initialized empty and idle, and that it produces results that are subject to statistical variation, so design and run a Simio Experiment to deal with both of these issues; make your own decisions about things like run length, number of replications, and Warm-Up Period, possibly after some trial and error. For each of the five steady-state queueing metrics, first compute numerical values for the queueing-theoretic steady-state output performance metrics \(W_q\), \(W\), \(L_q\), \(L\), and \(\rho\) from your solutions to Problem 9 in Chapter 2, and then compare these with your simulation estimates and confidence intervals. All time units are in minutes, and use minutes as well throughout your Simio model. Take the arrival rate to be \(\lambda\) = 1 per minute, and the service rate to be \(\mu\) = 1/0.9 per minute.

lambda <- 1
mu <- 1 / 0.9

rho <- lambda / mu
W_q <- (lambda / mu^2) / (2 * (1 - (lambda / mu)))
L_q <- lambda * W_q
W = W_q + (1 / mu)
L = lambda * W

sim <- read.csv("MD1_Model_091717_Model_Experiment1_ResultsSummary.csv")
avg <- sim %>% 
  filter(Statistic.Type %in% c("Average","Average (Minutes)")) %>%
  filter(!(Data.Source %in% c("Processing", "[DestroyedEntities]")) &
           Data.Item %in% c("NumberInSystem", 
                          "TimeInSystem",
                          "UnitsUtilized",
                          "NumberInStation",
                          "TimeInStation")) %>%
  select(Data.Item, Average) %>%
  spread(Data.Item, Average) %>%
  select(TimeInStation, NumberInStation, TimeInSystem, NumberInSystem, UnitsUtilized)

avg <- cbind("simulated", avg)
colnames(avg) <- c("Results", "W_q", "L_q", "W", "L", "rho")

theory <- data.frame("theoretical", W_q, L_q, W, L, rho)
colnames(theory) <- c("Results", "W_q", "L_q", "W", "L", "rho")
theory[, 2:6] <- as.numeric(theory[, 2:6])

tab <- rbind(theory, avg)

knitr::kable(tab,
             caption = paste("M/D/1 queue system steady-state conditions with", 
                                        "$\\theta$", "= 1 per minute and ",
                                        "$\\mu$", "= 1/0.9 per minute"))
M/D/1 queue system steady-state conditions with \(\theta\) = 1 per minute and \(\mu\) = 1/0.9 per minute
Results W_q L_q W L rho
theoretical 4.050000 4.050000 4.950000 4.950000 0.9000000
simulated 4.077048 4.081291 4.977057 4.982218 0.9009273 
sim <- sim %>% 
  filter(Statistic.Type %in% c("Average","Average (Minutes)")) %>%
  filter(!(Data.Source %in% c("Processing", "[DestroyedEntities]")) &
           Data.Item %in% c("NumberInSystem", 
                          "TimeInSystem",
                          "UnitsUtilized",
                          "NumberInStation",
                          "TimeInStation")) %>%
  select(Data.Item, Average, Half.Width)

sim <- sim %>% mutate(CI.Lower = Average - Half.Width,
               CI.Upper = Average + Half.Width) %>%
  arrange(Data.Item)

knitr::kable(sim, 
             caption = paste("Estimates and confidence intervals from Simio",
                             "experiment with 10 iterations of 56-day long", "
                             simulation with 60-hour warm-up period"))
Estimates and confidence intervals from Simio experiment with 10 iterations of 56-day long simulation with 60-hour warm-up period
Data.Item Average Half.Width CI.Lower CI.Upper
NumberInStation 4.0812911 0.0711307 4.0101604 4.1524218
NumberInSystem 4.9822184 0.0723526 4.9098658 5.0545710
TimeInStation 4.0770484 0.0660214 4.0110270 4.1430699
TimeInSystem 4.9770574 0.0660271 4.9110304 5.0430845
UnitsUtilized 0.9009273 0.0021112 0.8988161 0.9030385
Results of experiment with 10 iterations of 56-day long simulation with 60-hour warm-up period in Simio

Results of experiment with 10 iterations of 56-day long simulation with 60-hour warm-up period in Simio