The purpose of this vignette is to demonstrate how to use queuecomputer to understand M/M/k queues. We consider one M/M/1 queue and two M/M/3 queues. The working follows the theoretical results for M/M/k queues as shown in Thomopoulos, N (2012).
library(queuecomputer)P_0_func <- function(rho, k){
sum_i <- rep(NA, k)
for(i in 0:I(k-1))
{
sum_i[i+1] <- rho^i / factorial(i)
}
p_0 <- 1/(sum(sum_i) + rho^k/(factorial(k - 1) * (k - rho)))
return(p_0)
}
P_n <- function(rho,n,k){
p_0 <- P_0_func(rho, k)
if(n <= k){
output <- rho^n / factorial(n) * p_0
} else {
output <- rho^n / (factorial(k) * k^(n-k)) * p_0
}
return(output)
}
Lq <- function(rho, k){
p_0 <- P_0_func(rho, k)
output <- p_0 * rho^{k+1} / ( factorial(k-1) * (k - rho)^2)
return(output)
}set.seed(1)
n_customers <- 1e6
lambda_a <- 1/1
lambda_s <- 1/0.8
interarrivals <- rexp(n_customers, lambda_a)
arrivals <- cumsum(interarrivals)
service <- rexp(n_customers, lambda_s)
rho <- (1/lambda_s) / (1/lambda_a)k = 1
p_0 <- P_n(rho, n=0, k)
### System lengths -----------------------
Vectorize(P_n, "n")(rho=rho, n=c(0:30), k = k)## [1] 0.2000000000 0.1600000000 0.1280000000 0.1024000000 0.0819200000
## [6] 0.0655360000 0.0524288000 0.0419430400 0.0335544320 0.0268435456
## [11] 0.0214748365 0.0171798692 0.0137438953 0.0109951163 0.0087960930
## [16] 0.0070368744 0.0056294995 0.0045035996 0.0036028797 0.0028823038
## [21] 0.0023058430 0.0018446744 0.0014757395 0.0011805916 0.0009444733
## [26] 0.0007555786 0.0006044629 0.0004835703 0.0003868563 0.0003094850
## [31] 0.0002475880### Estimated queue length -----------------
LQ <- Lq(rho, k)
LQ## [1] 3.2### Estimated units in system -----------
LQ + rho## [1] 4Ws = 1/lambda_s
Wq = LQ / lambda_a
W = Ws + Wq
Wq # Mean waiting time (time in queue)## [1] 3.2W # Mean response time (time in system)## [1] 4MM1 <- queue_step(arrivals = arrivals, service = service, servers = k)
MM1_summary <- summary(MM1)
MM1_summary$slength_sum## # A tibble: 52 x 2
## queuelength proportion
## <int> <dbl>
## 1 0 0.200
## 2 1 0.160
## 3 2 0.128
## 4 3 0.102
## 5 4 0.0817
## 6 5 0.0655
## 7 6 0.0524
## 8 7 0.0420
## 9 8 0.0337
## 10 9 0.0270
## # … with 42 more rows# Mean queue length
MM1_summary$qlength_mean## [1] 3.210141# Mean system length (number of customers in system)
MM1_summary$slength_mean## [1] 4.009881MM1_summary$mwt # Mean waiting time## [1] 3.212666MM1_summary$mrt # Mean response time## [1] 4.01304k = 3
p_0 <- P_n(rho, 0, k)
### System lengths -----------------------
Vectorize(P_n, "n")(rho=rho, n=c(0:30), k = k)## [1] 4.471545e-01 3.577236e-01 1.430894e-01 3.815718e-02 1.017525e-02
## [6] 2.713400e-03 7.235732e-04 1.929529e-04 5.145410e-05 1.372109e-05
## [11] 3.658958e-06 9.757221e-07 2.601926e-07 6.938468e-08 1.850258e-08
## [16] 4.934022e-09 1.315739e-09 3.508638e-10 9.356368e-11 2.495031e-11
## [21] 6.653417e-12 1.774245e-12 4.731319e-13 1.261685e-13 3.364493e-14
## [26] 8.971982e-15 2.392529e-15 6.380076e-16 1.701354e-16 4.536943e-17
## [31] 1.209851e-17### Estimated queue length -----------------
LQ <- Lq(rho, k)
LQ## [1] 0.01892092### Estimated units in system -----------
LQ + rho## [1] 0.8189209Ws = 1/lambda_s
Wq = LQ / lambda_a
W = Ws + Wq
Wq # Mean waiting time (time in queue)## [1] 0.01892092W # Mean response time (time in system)## [1] 0.8189209MM3 <- queue_step(arrivals = arrivals, service = service, servers = k)
MM3_summary <- summary(MM3)
MM3_summary$slength_sum## # A tibble: 13 x 2
## queuelength proportion
## <int> <dbl>
## 1 0 0.447
## 2 1 0.358
## 3 2 0.143
## 4 3 0.0380
## 5 4 0.0101
## 6 5 0.00268
## 7 6 0.000704
## 8 7 0.000193
## 9 8 0.0000399
## 10 9 0.0000120
## 11 10 0.00000201
## 12 11 0.00000276
## 13 12 0.000000106# Mean queue length
MM3_summary$qlength_mean## [1] 0.01865591# Mean system length (number of customers in system)
MM3_summary$slength_mean## [1] 0.8184002MM3_summary$mwt # Mean waiting time## [1] 0.01867057MM3_summary$mrt # Mean response time## [1] 0.8190443# Setup ----------
set.seed(2)
n_customers <- 1e6
lambda_a <- 1/1
lambda_s <- 1/2.5
interarrivals <- rexp(n_customers, lambda_a)
arrivals <- cumsum(interarrivals)
service <- rexp(n_customers, lambda_s)
rho <- (1/lambda_s) / (1/lambda_a)rho <- (1/lambda_s) / (1/lambda_a)
# MM3 queue ------------------------------
k = 3
## Theoretical -------------------
p_0 <- P_n(rho, 0, k)
### System lengths -----------------------
Vectorize(P_n, "n")(rho=rho, n=c(0:30), k = k)## [1] 0.0449438202 0.1123595506 0.1404494382 0.1170411985 0.0975343321
## [6] 0.0812786101 0.0677321751 0.0564434792 0.0470362327 0.0391968606
## [11] 0.0326640505 0.0272200421 0.0226833684 0.0189028070 0.0157523392
## [16] 0.0131269493 0.0109391244 0.0091159370 0.0075966142 0.0063305118
## [21] 0.0052754265 0.0043961888 0.0036634906 0.0030529089 0.0025440907
## [26] 0.0021200756 0.0017667297 0.0014722747 0.0012268956 0.0010224130
## [31] 0.0008520108### Estimated queue length -----------------
LQ <- Lq(rho, k)
LQ## [1] 3.511236### Estimated units in system -----------
LQ + rho## [1] 6.011236### Waiting times -----------
Ws = 1/lambda_s
Wq = LQ / lambda_a
W = Ws + Wq
Wq # Mean waiting time (time in queue)## [1] 3.511236W # Mean response time (time in system)## [1] 6.011236MM3_2 <- queue_step(arrivals = arrivals, service = service, servers = k)
MM3_2_summary <- summary(MM3_2)
MM3_2_summary$slength_sum## # A tibble: 53 x 2
## queuelength proportion
## <int> <dbl>
## 1 0 0.0450
## 2 1 0.113
## 3 2 0.140
## 4 3 0.118
## 5 4 0.0980
## 6 5 0.0820
## 7 6 0.0679
## 8 7 0.0568
## 9 8 0.0471
## 10 9 0.0391
## # … with 43 more rows# Mean queue length
MM3_2_summary$qlength_mean## [1] 3.472223# Mean system length (number of customers in system)
MM3_2_summary$slength_mean## [1] 5.971736MM3_2_summary$mwt # Mean waiting time## [1] 3.46892MM3_2_summary$mrt # Mean response time## [1] 5.966074Thomopoulos, N (2012). Fundamentals of Queuing Systems: Statistical Methods for Analyzing Queuing Models. Springer Science & Business Media