HW4_604
jbrnbrg
September 19, 2017
library(tidyverse)
library(stringr)
library(RCurl) Complete Chapter 2 Problem 9 from Kelton followed by Chapter 4 Problem 15.
2.6.9:
pg36.
\(W_{q}=\frac{\lambda(\sigma^2+1/\mu^2)}{2(1-\lambda/\mu)}\) where \(\sigma\) is the standard deviation of the service times. Since this problem sets the service time to a constant, \(D\), there is no standard deviation in the service time. This leads to a \(W_{q}\) defined by:
\(W_{q}=\frac{\lambda}{2\mu(\mu-\lambda)}\) for \(M/D/1\) systems where \(D\) is a constant. Further, since I know that \(L_{q}=\lambda W_{q}\) for an \(M/M/1\) system, I know that for an \(M/D/1\) system, it’ll look like: \(L_{q}=\lambda W_{q}=\frac{\lambda^2}{2\mu(\mu-\lambda)}\). Continuing on, I know that \(W=W_{q}+E(S)\) and \(E(S)=1/\mu\), which would lead to: \(W=W_{q}+1/\mu=\frac{\lambda}{2\mu(\mu-\lambda)}+\frac{1}{\mu}=\frac{2\mu-\lambda}{2\mu(\mu-\lambda)}\). Now that I’ve got the calculations for \(W\),\(W_{q}\), and \(L_{q}\) and we know that \(\rho=\lambda/\mu\), I can find the remaining value for \(L\) via \(L=\lambda W=\frac{\lambda (2\mu-\lambda)}{2\mu(\mu-\lambda)}\).
Answer:
\(\therefore\) for an \(M/D/1\) system, the five steady-state queuing metrics are defined by:
\(W_{q}=\frac{\lambda}{2\mu(\mu-\lambda)}\)
\(W=\frac{2\mu-\lambda}{2\mu(\mu-\lambda)}\)
\(L_{q}=\frac{\lambda^2}{2\mu(\mu-\lambda)}\)
\(L=\frac{\lambda (2\mu-\lambda)}{2\mu(\mu-\lambda)}\)
\(\rho=\lambda/\mu\)
Now, looking ahead a bit, I select \(\lambda=1\) \(D=\frac{1}{0.9}\) as they’re coming up in the next problem and calculate the respective results using the above calcs for an \(M/D/1\) and an \(M/M/1\):
MD1 calc:
my_lambda <- 1
my_mu <- 1/0.9
Wq <- my_lambda/(2*my_mu*(my_mu - my_lambda))
W <- (2*my_mu - my_lambda)/(2*my_mu*(my_mu - my_lambda))
Lq <- my_lambda^2/(2*my_mu*(my_mu - my_lambda))
L <- my_lambda*(2*my_mu - my_lambda)/(2*my_mu*(my_mu - my_lambda))
rho <- my_lambda/my_mu
x <- (t(c("M/D/1 calc", my_lambda, round(my_mu,4), Wq,
W, Lq, L, rho, 1-rho))) %>%
as.data.frame(stringsAsFactors=F) %>%
rename(name = V1, lambda=V2, mu=V3, Wq = V4,
W = V5, Lq = V6, L = V7, rho = V8, p0 = V9)
row.names(x) <- NULLMM1 calc:
my_MMC_func <- function(t_inter, t_svc, c, prob = 1, nm){
lambda <- (1/t_inter)*prob # inter-arrival times
mu <- 1/t_svc # svc rate
rho <- lambda/(c*mu) # utilization per enti
i <- (0:(c-1)) # the index to sum over
# prob system is empty i.e. p(0)
p_empty <- 1/(((c*rho)^c)/(factorial(c)*(1-rho)) + sum(((c*rho)^i)/factorial(i)))
Lq <- (rho*((c*rho)^c)*p_empty)/(factorial(c)*(1-rho)^2)#steady-st avg enti in q
Wq <- Lq/lambda # avg time in queue
L <- Lq + rho*c # steady-state avg of enti
W <- Wq + t_svc # steady-st avg time in sy
Lq <- lambda*Wq # steady-st avg enti in q
results <- t(c(nm, round(lambda,3), round(mu,3), round(Wq,3), round(W,3),
round(Lq,3), round(L,3), round(rho,3), round(p_empty,3))) %>%
as.data.frame(stringsAsFactors=FALSE) %>%
rename(name = V1, lambda=V2, mu=V3, Wq = V4,
W = V5, Lq = V6, L = V7, rho = V8, p0 = V9) %>%
mutate_each(funs(if(is.numeric(.)) as.numeric(.) else .))
return(results)
}
my_md1<-my_MMC_func(1, .9, 1, 1, "M/M/1")
combo_mm1_md1 <- rbind(x, my_md1)
row.names(combo_mm1_md1) <- NULL| name | lambda | mu | Wq | W | Lq | L | rho | p0 |
|---|---|---|---|---|---|---|---|---|
| M/D/1 calc | 1 | 1.1111 | 4.05 | 4.95 | 4.05 | 4.95 | 0.9 | 0.1 |
| M/M/1 | 1 | 1.111 | 8.1 | 9 | 8.1 | 9 | 0.9 | 0.1 |
As shown in the table above, the \(M/D/1\) configuration is twice as fast as the \(M/M/1\) configuration.
4.10.15:
pg127.
Answer:
Below, please find a screen capture of my JaamSim sim for this M/D/1 system:
An M/D/1 in JaamSim
Above, I used a uniform distribution with limits set to [1/.9, 1/.9] to force it to a constant.
Before I could even begin: To review the report file out of JaamSim, I needed to create a function to read the data from multiple runs, and then combine them to get the averages - this wasn’t super easy but I learned a lot. Also in the code below is the creation of a helper function for calculating the confidence intervals.
If you’re interested in following along, my the data can be accessed here.
multi_rep_reader <- function(file, lambda){
# A function that combines the information from multiple runs
# from the standard JaamSim .rep file
# adds columns for t-confidence interval.
xx <- read.delim(file,
header = F,
stringsAsFactors = F) %>%
mutate(run_n = ifelse(grepl("#", V1),V1, NA)) %>%
fill(run_n) %>%
mutate(run_n = str_extract(run_n, "[:digit:]")) %>%
filter(!grepl("\\*|##|Simulation", V1),
!grepl("QueueLengthTimes", V2)) %>%
mutate(V3 = as.numeric(V3)) %>%
group_by(V1, V2, V4) %>%
summarise(r_n = n(),
r_min = min(V3),
r_max = max(V3),
r_avg = mean(V3, na.rm = T),
r_sd = sd(V3))
return(xx)
}
z <- multi_rep_reader("newMD1_24hr_50reps.rep", 1)
my_ct <- function(s, n, ci){
# helper function to get the student t half-tails.
z_t <- qt( 1-(1-ci)/2, df = n-1)
h <- z_t * s/sqrt(n)
return(h)
}
#View(z)
#unique(z$V1)Next, I create a function to calculate all of the Little’s Law values - one set with confidence intervals is created and one set is for presentation:
rep_metrics <- function(sim_rpt, inter_t, t_svc, conf_i){
g <- sim_rpt %>%
filter(grepl("Server1|Queue1", V1),
grepl("AverageQueueTime|QueueLengthAverage|NumberProcessed|Utilisation", V2))
lam <- 1/inter_t
mu <- 1/t_svc
n <- as.numeric(g[2,4])
ci <- .95
Wq_sd <- as.numeric(g[1, 8])
Lq_sd <- as.numeric(g[3, 8])
rho_sd <- as.numeric(g[5, 8])
rho <- as.numeric(g[5, 7])
Wq <- as.numeric(g[1, 7])
L <- lam * Wq
W <- Wq + t_svc
Lq <- as.numeric(g[3, 7])
p0 <- 1 - rho
rs <- data.frame(metric = c("rho", "L", "Lq", "W", "Wq"),
sim_avgs = c(rho, L, Lq, W, Wq),
sim_h = c(my_ct(rho_sd, n, ci),
my_ct(Lq_sd, n, ci),
my_ct(Lq_sd, n, ci),
my_ct(Wq_sd, n, ci),
my_ct(Wq_sd, n, ci)),
stringsAsFactors = F)
rt <- t(c("M/D/1 JaamSim", lam, round(mu,4), round(Wq,4), round(W,4),
round(Lq,4), round(L,4), round(rho,8), round(p0,8))) %>%
as.data.frame(stringsAsFactors=FALSE) %>%
rename(name = V1, lambda=V2, mu=V3, Wq = V4,
W = V5, Lq = V6, L = V7, rho = V8, p0 = V9)
return(list("rs"= rs, "rt" = rt))
}
b <- rep_metrics(z, 1, .9, .95)
combo_all <- rbind(b$rt, combo_mm1_md1)First, we’ll review the 95% confidence intervals of all of the performance metrics from the simulation that I ran (24hrs, 10min initialization, 50 reps)
| metric | sim_avgs | sim_h |
|---|---|---|
| rho | 0.9974439 | 0.0010766 |
| L | 79.0795597 | 6.4645415 |
| Lq | 79.9464861 | 6.4645415 |
| W | 79.9795597 | 5.8059479 |
| Wq | 79.0795597 | 5.8059479 |
Finally, we’ll compare the performance metrics of my M/D/1 simulation with calculated M/D/1 and an M/M/1 system:
| name | lambda | mu | Wq | W | Lq | L | rho | p0 |
|---|---|---|---|---|---|---|---|---|
| M/D/1 JaamSim | 1 | 1.1111 | 79.0796 | 79.9796 | 79.9465 | 79.0796 | 0.99744387 | 0.00255613 |
| M/D/1 calc | 1 | 1.1111 | 4.05 | 4.95 | 4.05 | 4.95 | 0.9 | 0.1 |
| M/M/1 | 1 | 1.111 | 8.1 | 9 | 8.1 | 9 | 0.9 | 0.1 |
my_MMC_func <- function(t_inter, t_svc, c, prob = 1, nm){
lambda <- (1/t_inter)*prob # inter-arrival times
mu <- 1/t_svc # svc rate
rho <- lambda/(c*mu) # utilization per enti
i <- (0:(c-1)) # the index to sum over
# prob system is empty i.e. p(0)
p_empty <- 1/(((c*rho)^c)/(factorial(c)*(1-rho)) + sum(((c*rho)^i)/factorial(i)))
Lq <- (rho*((c*rho)^c)*p_empty)/(factorial(c)*(1-rho)^2)#steady-st avg enti in q
Wq <- Lq/lambda # avg time in queue
L <- Lq + rho*c # steady-state avg of enti
W <- Wq + t_svc # steady-st avg time in sy
Lq <- lambda*Wq # steady-st avg enti in q
results <- t(c(nm, round(lambda,3), round(mu,3), round(Wq,3), round(W,3),
round(Lq,3), round(L,3), round(rho,3), round(p_empty,3))) %>%
as.data.frame(stringsAsFactors=FALSE) %>%
rename(name = V1, lambda=V2, mu=V3, Wq = V4,
W = V5, Lq = V6, L = V7, rho = V8, p0 = V9) %>%
mutate_each(funs(if(is.numeric(.)) as.numeric(.) else .))
return(results)
}
my_md1<-my_MMC_func(1, .9, 1, 1, "M/M/1")
#rbind(x, my_md1, jaamsim_results)
#names(x)