Discussion Problem

Pick a, b, c, or d from Problem 14 of 4.10 in Kelton and solve. Discuss your Simio Model. Include a snapshot of your model as well as an interpretation of results. Discuss any issues you had building the model.

Problem 4.10.14

Build Simio models to confirm and cross-check steady-state queuing theoretic results for the four specific queuing models whose exact steady-state output performance metrics are given in Section 2.3. Remember that your Simio models are initialized empty and idle, and that they produce results that are subject to statistical variation, so design and run Simio Experiments to deal with both of these issues; make your own decision about things like run length, number of replications, possibly after some trial an error. In each case, first compute numerical values for the queuing-theoretic steady-state output performance metrics \(W_q\), \(W\), \(L_q\), \(L\), and \(\rho\) from the results in Section 2.3, 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.

\[\lambda = \frac{N}{T_a}, \quad\quad \mu = \frac{N}{T_s}, \quad\quad \rho = \frac{\lambda}{c\mu}\]

The realizations of the Simio outputs closely approximate the theoretic results. Model Entity > Number in System > Average approximates the mean number of customers in the queue \(L_q\); Model Entity > Time in System > Average approximates the mean wait time for customers in the queue excluding service times \(W_q\); Server > UnitsUtilized > Average approximates the utilization rate of the server \(\rho\); Server > InputBuffer > Number in System > Average approximates the mean number of customers in the queues \(L\); and Server > InputBuffer > Time in System > Average approximates the mean wait time for customers in the queues including service times \(W\). The upper and lower bounds were calculated using the estimates \(\pm\) halfwidth.

simio_eval <- function(averages, halfwidths, theoretical) {
  avg <- averages; half <- halfwidths
  L_sim <- c(avg[1], avg[1] - half[1], avg[1] + half[1])
  W_sim <- c(avg[2], avg[2] - half[2], avg[2] + half[2])
  rho_sim <-c(avg[3], avg[3] - half[3], avg[3] + half[3]) / theoretical[6]
  Lq_sim <- c(avg[4], avg[4] - half[4], avg[4] + half[4])
  Wq_sim <-c(avg[5], avg[5] - half[5], avg[5] + half[5])
  bools <- c(theoretical[1] >= L_sim[2] & theoretical[1] <= L_sim[3],
             theoretical[2] >= W_sim[2] & theoretical[2] <= W_sim[3],
             theoretical[3] >= rho_sim[2] & theoretical[3] <= rho_sim[3],
             theoretical[4] >= Lq_sim[2] & theoretical[4] <= Lq_sim[3],
             theoretical[5] >= Wq_sim[2] & theoretical[5] <= Wq_sim[3])
  return(cbind(rbind(Lq_sim, Wq_sim, rho_sim, L_sim, W_sim), bools))
}

Problem 4.10.14.a

\(M/M/1\) queue with arrival rate \(\lambda=1\) per minute, and the service rate \(\mu=1/0.9\) per minute.

\[L_q = \frac{\rho^2}{1 - \rho}, \quad\quad L = \frac{\lambda}{\mu - \lambda}, \quad\quad W_q = \frac{\rho}{\mu\left(1 - \rho\right)}, \quad\quad W = \frac{1}{\mu\left(1 - \rho\right)}\]

T_a <- 1
T_s <- 9 / 10
lambda <- 1 / T_a
mu <- 1 / T_s
c <- 1
rho = lambda / (c * mu)
L_q <- rho^2 / (1 - rho)
L <- lambda / (mu - lambda)
W <- 1 / (mu * (1 - rho))
W_q <- rho / (mu * (1 - rho))

One Source with Interarrival Time of Random.Exponential( 1 ). One Server with Processing Time of Random.Exponential( 9/10 ). One unmodified Sink. Experiment conducted with 200 Replications of 50 hour runs using Warm-up Period of 25 Hours.

Simulation Model Mitigate Empty/Idle Time Pivot Grid 
avg <- c(8.9423, 8.8946, 0.9003, 8.0420, 7.9946)
half <- c(0.5134, 0.4935, 0.0043, 0.5102, 0.4918)
simio <- simio_eval(avg, half, c(L, W, rho, L_q, W_q, c))
Metric Simio Estimate Lower Bound Upper Bound Theoretic Captured
\(L\) 8.9423 8.4289 9.4557 9 TRUE
\(W\) 8.8946 8.4011 9.3881 9 TRUE
\(\rho\) 0.9003 0.896 0.9046 0.9 TRUE
\(L_q\) 8.042 7.5318 8.5522 8.1 TRUE
\(W_q\) 7.9946 7.5028 8.4864 8.1 TRUE

Problem 4.10.14.b

\(M/M/4\) queue with arrival rate \(\lambda=2.4\) per minute, and the service rate \(\mu=0.7\) per minute for each of the four individual servers (the same parameters used in the mmc.exe command line program shown in Figure 2.2).

\[\rho \left( 0 \right) ={ \left[ \frac { { \left( c\rho \right) }^{ c } }{ c!\left( 1-\rho \right) } +\sum _{ n=0 }^{ c-1 }{ \frac { { \left( c\rho \right) }^{ n } }{ n! } } \right] }^{ -1 },\quad \quad { L }_{ q }=\frac { { \rho \left( c\rho \right) }^{ c }\rho \left( 0 \right) }{ c!{ \left( 1-\rho \right) }^{ 2 } }, \quad\quad L=c\rho+L_q, \quad\quad W=\frac{L}{\lambda}, \quad\quad W_q=W-\mu^{-1}\]

rho_ini <- function(c, rho) {
  sum1 <- (c * rho)^c / (factorial(c) * (1 - rho))
  iteration <- numeric(c)
  for (n in 0:(c - 1)) {iteration[n + 1] <- (c * rho)^n / factorial(n)}
  sum2 <- sum(iteration)
  return((sum1 + sum2)^(-1))
}
T_a <- 5 / 12
T_s <- 10 / 7
lambda <- 1 / T_a
mu <- 1 / T_s
c <- 4
rho <- lambda / (c * mu)
rho0 <- rho_ini(c, rho)
L_q <- (rho * (c * rho)^c * rho0) / (factorial(c) * (1 - rho)^2)
L <- c * rho + L_q
W <- L / lambda
W_q <- W - mu^(-1)
data.frame(rho0)
##         rho0
## 1 0.01744393

One Source with Interarrival Time of Random.Exponential( 5/12 ). One Server with Processing Time of Random.Exponential( 10/7 ) and Initial Capacity 4. One unmodified Sink. Experiment conducted with 200 Replications of 50 hour runs using Warm-up Period of 25 Hours.

Simulation Model Mitigate Empty/Idle Time Pivot Grid 
avg <- c(7.7971, 3.2486, 3.4295, 4.3676, 1.8180)
half <- c(0.2000, 0.0794, 0.0111, 0.1922, 0.0777)
simio <- simio_eval(avg, half, c(L, W, rho, L_q, W_q, c))
Metric Theoretic Simio Estimate Lower Bound Upper Bound Captured
\(L\) 7.6468494 7.7971 7.5971 7.9971 TRUE
\(W\) 3.1861873 3.2486 3.1692 3.328 TRUE
\(\rho\) 0.8571429 0.857375 0.8546 0.86015 TRUE
\(L_q\) 4.218278 4.3676 4.1754 4.5598 TRUE
\(W_q\) 1.7576158 1.818 1.7403 1.8957 TRUE

Problem 4.10.14.c

\(M/G/1\) queue with arrival rate \(\lambda=1\) per minute, and the service-time distribution’s being \(\Gamma\left(2.00, 0.45\right)\) (shape and scale parameters respectively). You may need to do some investigation about properties of the gamma distribution, perhaps via the web links1,2 in Section 6.1.3.

\[\mu = k\theta = \frac{\alpha}{\beta}, \quad \sigma^2 = k\theta^2 = \frac{\alpha}{\beta^{2}}, \quad L_q = \frac{\rho^2\left(1 + \sigma^2\mu^2\right)}{2\left(1 - \rho\right)}, \quad L=\rho+L_q, \quad { W }_{ q }=\frac { \lambda \left( { \sigma }^{ 2 }+{ \mu }^{ -2 } \right) }{ 2\left( 1+\rho \right) }, \quad W=W_q+\mu^{-1}\]

T_a <- 1
lambda <- 1 / T_a
# Arrival rate must be less than service rate.
# Must have T_a<mu for steady-state to exist.           
# |rho|<=1 => lambda<=mu => k,theta form
k <- 2
theta <- 9 / 20
T_s <- k * theta
mu <- 1 / T_s
sigma <- sqrt(k * theta^2)
c <- 1
rho = lambda / (c * mu)
L_q <- rho^2 * (1 + sigma^2 * mu^2) / (2 * (1 - rho))
L <- rho + L_q
W_q <- lambda * (sigma^2 + mu^(-2)) / (2 * (1 - rho))
W <- W_q + mu^(-1)
data.frame(mu, sigma)
##         mu     sigma
## 1 1.111111 0.6363961

One Source with Interarrival Time of Random.Exponential( 1 ). One Server with Processing Time of Random.Gamma( 2, 9/20 ). One unmodified Sink. Experiment conducted with 10 Replications of 90 day runs using Warm-up Period of 14 days.

Simulation Model Mitigate Empty/Idle Time Pivot Grid 
avg <- c(7.0028, 6.9980, 0.9008, 6.1020, 6.0978)
half <- c(0.1512, 0.1459, 0.0018, 0.1505, 0.1459)
simio <- simio_eval(avg, half, c(L, W, rho, L_q, W_q, c))
Metric Theoretic Simio Estimate Lower Bound Upper Bound Captured
\(L\) 6.975 7.0028 6.8516 7.154 TRUE
\(W\) 6.975 6.998 6.8521 7.1439 TRUE
\(\rho\) 0.9 0.9008 0.899 0.9026 TRUE
\(L_q\) 6.075 6.102 5.9515 6.2525 TRUE
\(W_q\) 6.075 6.0978 5.9519 6.2437 TRUE

Problem 4.10.14.d

\(G/M/1\) queue with interarrival time distribution’s being continuous uniform between 1 and 5, and service rate \(\mu=0.4\) per minute (the same situation shown in Figure 2.3).

\[z=\int _{ a }^{ b }{ { e }^{ -\mu t\left( 1-z \right) }\left( \frac { 1 }{ b-a } \right) dt } =-\frac { { e }^{ -\mu b\left( 1-z \right) }-{ e }^{ -\mu a\left( 1-z \right) } }{ \mu \left( b-a \right) \left( 1-z \right) } ,\quad W=\frac { 1 }{ \mu \left( 1-z \right) } ,\quad { W }_{ q }=W-{ \mu }^{ -1 },\quad { L }_{ q }=\frac { { W }_{ q } }{ T_{ a } } ,\quad L=\frac { W }{ T_{ a } } \]

f <- function(z, mu, a, b) {
  Fb <- exp(-mu * b * (1 - z))
  Fa <- exp(-mu * a * (1 - z))
  k <- mu * (b - a) * (1 - z)
  z + ((Fb - Fa) / k)
}
a <- 1
b <- 5
T_a <- (a + b) / 2
lambda <- 1 / T_a
T_s <- 5 / 2
mu <- 1 / T_s
c <- 1
rho = lambda / (c * mu)
epsilon <- 1e-15
z <- uniroot(f, interval=c(0+epsilon, 1-epsilon), mu=mu, a=a, b=b)$root
f_z <- f(z, mu, a, b)
W <- (mu * (1 - z))^(-1)
W_q <- W - 1 / mu
L_q <- W_q / T_a
L <- W / T_a
data.frame(z, 'f_z'=format(f_z, scientific=FALSE))
##           z            f_z
## 1 0.7234921 0.000001249264

One Source with Interarrival Time of Random.Uniform( 1, 5 ). One Server with Processing Time of Random.Exponential( 5/2 ). One unmodified Sink. Experiment conducted with 200 Replications of 50 hour runs using Warm-up Period of 25 Hours.

Simulation Model Mitigate Empty/Idle Time Pivot Grid 
avg <- c(2.9917, 8.9683, 0.8301, 2.1616, 6.4766)
half <- c(0.1530, 0.4526, 0.0056, 0.1486, 0.4412)
simio <- simio_eval(avg, half, c(L, W, rho, L_q, W_q, c))
Metric Theoretic Simio Estimate Lower Bound Upper Bound Captured
\(L\) 3.0137781 2.9917 2.8387 3.1447 TRUE
\(W\) 9.0413344 8.9683 8.5157 9.4209 TRUE
\(\rho\) 0.8333333 0.8301 0.8245 0.8357 TRUE
\(L_q\) 2.1804448 2.1616 2.013 2.3102 TRUE
\(W_q\) 6.5413344 6.4766 6.0354 6.9178 TRUE

Assignment Problems

Complete Chapter 2 Problem 9 from Kelton followed by Chapter 4 Problem 15.

Problem 2.6.9

Find all five of the steady-state queuing 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.

\[\lambda = \frac{N}{T_a}, \quad\quad \mu = \frac{N}{T_s}, \quad\quad \rho = \frac{\lambda}{c\mu}\]

\[\sigma^2 = 0 \Rightarrow L_q = \frac{\rho^2\left(1 + \sigma^2\mu^2\right)}{2\left(1 - \rho\right)} = \frac{\rho^2}{2\left(1 - \rho\right)}, \quad L=\rho+L_q, \quad { W }_{ q } = \frac { \lambda \left( { \sigma }^{ 2 }+{ \mu }^{ -2 } \right) }{ 2\left( 1+\rho \right) } = \frac { -\lambda { \mu }^{ -2 }}{ 2\left( 1+\rho \right) }, \quad W=W_q+\mu^{-1}\] No variation implies that \(\sigma^2 = 0\). Above we can see how equations for an \(M/D/1\) queue with some variance, \(\sigma^2 \neq 0\), are reduced (degenerated) when variance is equivalent to zero.

Problem 4.10.15

Build a Simio model to confirm and cross-check the steady-state queuing theoretic results to the \(M/D/1\) queue of Problem 9 in Chapter 2. Remember that your 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 decision about things like run length, number of replications, and Warm-up Period, possibly after some trial an error. For each of the five steady-state queuing metrics, first compute numerical values for the queuing-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.

T_a <- 1
T_s <- 9 / 10
lambda <- 1 / T_a
mu <- 1 / T_s
sigma <- 0
c <- 1
rho = lambda / (c * mu)
L_q <- rho^2 * (1 + sigma^2 * mu^2) / (2 * (1 - rho))
L <- rho + L_q
W_q <- lambda * (sigma^2 + mu^(-2)) / (2 * (1 - rho))
W <- W_q + mu^(-1)

One Source with Interarrival Time of Random.Exponential( 1 ). One Server with Processing Time of ( 9/10 ). One unmodified Sink. Experiment conducted with 200 Replications of 50 hour runs using Warm-up Period of 25 Hours.

Simulation Model Mitigate Empty/Idle Time Pivot Grid 
avg <- c(5.1113, 5.0698, 0.9025, 4.2088, 4.1699)
half <- c(0.2275, 0.2098, 0.0032, 0.2251, 0.2099)
simio <- simio_eval(avg, half, c(L, W, rho, L_q, W_q, c))
Metric Theoretic Simio Estimate Lower Bound Upper Bound Captured
\(L\) 4.95 5.1113 4.8838 5.3388 TRUE
\(W\) 4.95 5.0698 4.86 5.2796 TRUE
\(\rho\) 0.9 0.9025 0.8993 0.9057 TRUE
\(L_q\) 4.05 4.2088 3.9837 4.4339 TRUE
\(W_q\) 4.05 4.1699 3.96 4.3798 TRUE

References

http://rpubs.com/janderman/309871

https://rpubs.com/josezuniga/305494

http://s1.daumcdn.net/editor/fp/service_nc/pencil/Pencil_chromestore.html

https://www.youtube.com/watch?v=Rs50G2CWT70&list=PLki7GH7osYZkN9lippFxp8E3kMOwgnJLG&index=3

https://www.youtube.com/watch?v=1Q2nIIM-L2Y&list=PLki7GH7osYZkN9lippFxp8E3kMOwgnJLG&index=4

Simio and Simulation: Modeling, Analysis, Applications 3d Ed. by W. David Kelton, Jeffrey S. Smith and David T. Sturrock with Simio software.

Discrete-Event Systems Simulation, 5th Edition (2010), by Jerry Banks, John S. Carlson, Barry L. Nelson,and David M. Nicol.