# Introduction

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)

## Queueing Functions

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)
}

## Setup

set.seed(1)

n_customers <- 250000

lambda_a <- 1/1
lambda_s <- 1/0.8

interarrivals <- rexp(n_customers, lambda_a)

arrivals <- cumsum(interarrivals)
arrival_df <- data.frame(ID = c(1:n_customers), times = arrivals)

service <- rexp(n_customers, lambda_s)

rho <- (1/lambda_s) / (1/lambda_a)

## MM1 queue

### Theoretical

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)

### Estimated units in system -----------
Lq(rho, k) + rho
## [1] 4
Ws = 1/lambda_s
Wq = LQ / lambda_a
W = Ws + Wq

Wq # Mean waiting time (time in queue)
## [1] 3.2
W # Mean response time (time in system)
## [1] 4

### Observed

MM1 <- queue_step(arrival_df = arrival_df, service = service, servers = k)

MM1_summary <- summary(MM1)

signif(MM1_summary$system_lengths, 4) ## 0 1 2 3 4 5 6 ## 1.986e-01 1.615e-01 1.278e-01 1.028e-01 8.173e-02 6.645e-02 5.365e-02 ## 7 8 9 10 11 12 13 ## 4.304e-02 3.386e-02 2.646e-02 2.129e-02 1.751e-02 1.377e-02 1.082e-02 ## 14 15 16 17 18 19 20 ## 8.500e-03 6.622e-03 5.275e-03 4.236e-03 3.399e-03 2.697e-03 2.361e-03 ## 21 22 23 24 25 26 27 ## 1.768e-03 1.347e-03 1.082e-03 8.410e-04 5.812e-04 3.819e-04 3.082e-04 ## 28 29 30 31 32 33 34 ## 2.494e-04 2.290e-04 1.697e-04 1.285e-04 1.197e-04 1.069e-04 7.485e-05 ## 35 36 37 38 39 40 41 ## 3.322e-05 3.442e-05 3.162e-05 4.523e-05 2.842e-05 1.681e-05 1.161e-05 ## 42 ## 3.202e-06 MM1_summary$queue_lengths %*% c(0:I(length(MM1_summary$queue_lengths) - 1)) # Mean queue length ## [,1] ## [1,] 3.146248 MM1_summary$system_lengths %*% c(0:I(length(MM1_summary$system_lengths) - 1)) # Mean system length (number of customers in system) ## [,1] ## [1,] 3.947628 MM1_summary$mwt # Mean waiting time
## [1] 3.144179
MM1_summary$mrt # Mean response time ## [1] 3.945074 ## MM3 queue ### Theoretical k = 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) ### Estimated units in system ----------- Lq(rho, k) + rho ## [1] 0.8189209 Ws = 1/lambda_s Wq = LQ / lambda_a W = Ws + Wq Wq # Mean waiting time (time in queue) ## [1] 0.01892092 W # Mean response time (time in system) ## [1] 0.8189209 ## Observed MM3 <- queue_step(arrival_df = arrival_df, service = service, servers = k) MM3_summary <- summary(MM3) signif(MM3_summary$system_lengths, 4)
##         0         1         2         3         4         5         6
## 4.462e-01 3.582e-01 1.436e-01 3.821e-02 1.014e-02 2.731e-03 7.169e-04
##         7         8         9        10        11
## 1.769e-04 3.082e-05 1.361e-05 3.202e-06 2.001e-06
MM3_summary$queue_lengths %*% c(0:I(length(MM3_summary$queue_lengths) - 1)) # Mean queue length
##            [,1]
## [1,] 0.01873655
MM3_summary$system_lengths %*% c(0:I(length(MM3_summary$system_lengths) - 1)) # Mean system length (number of customers in system)
##           [,1]
## [1,] 0.8201465
MM3_summary$mwt # Mean waiting time ## [1] 0.0187137 MM3_summary$mrt # Mean response time
## [1] 0.8196092

## MM3 queue (second)

# Setup ----------

set.seed(2)

n_customers <- 250000

lambda_a <- 1/1
lambda_s <- 1/2.5

interarrivals <- rexp(n_customers, lambda_a)

arrivals <- cumsum(interarrivals)
arrival_df <- data.frame(ID = c(1:n_customers), times = arrivals)

service <- rexp(n_customers, lambda_s)

### Theoretical

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)

### Estimated units in system -----------
Lq(rho, k) + 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.511236
W # Mean response time (time in system)
## [1] 6.011236

### Observed

MM3_2 <- queue_step(arrival_df = arrival_df, service = service, servers = k)

MM3_2_summary <- summary(MM3_2)

signif(MM3_2_summary$system_lengths, 4) ## 0 1 2 3 4 5 6 ## 4.486e-02 1.109e-01 1.377e-01 1.158e-01 9.756e-02 8.142e-02 6.784e-02 ## 7 8 9 10 11 12 13 ## 5.619e-02 4.724e-02 3.952e-02 3.270e-02 2.806e-02 2.242e-02 1.888e-02 ## 14 15 16 17 18 19 20 ## 1.519e-02 1.295e-02 1.104e-02 9.402e-03 7.739e-03 6.690e-03 5.541e-03 ## 21 22 23 24 25 26 27 ## 4.442e-03 3.859e-03 3.244e-03 2.623e-03 2.259e-03 2.232e-03 1.698e-03 ## 28 29 30 31 32 33 34 ## 1.580e-03 1.200e-03 1.099e-03 8.428e-04 7.594e-04 6.011e-04 5.342e-04 ## 35 36 37 38 39 40 41 ## 4.300e-04 3.823e-04 3.471e-04 3.062e-04 2.352e-04 2.364e-04 1.611e-04 ## 42 43 44 45 46 47 48 ## 1.964e-04 1.984e-04 1.647e-04 9.217e-05 7.935e-05 5.050e-05 4.168e-05 ## 49 50 51 52 53 54 55 ## 2.725e-05 4.008e-05 2.204e-05 2.525e-05 2.685e-05 2.244e-05 2.525e-05 ## 56 57 58 59 60 61 62 ## 1.924e-05 2.164e-05 4.248e-05 3.286e-05 4.328e-05 3.967e-05 1.202e-05 ## 63 64 65 66 67 68 69 ## 6.412e-06 8.015e-06 8.416e-06 2.405e-06 4.809e-06 9.618e-06 7.614e-06 ## 70 71 72 ## 4.408e-06 4.008e-07 1.603e-06 MM3_2_summary$queue_lengths %*% c(0:I(length(MM3_2_summary$queue_lengths) - 1)) # Mean queue length ## [,1] ## [1,] 3.646973 MM3_2_summary$system_lengths %*% c(0:I(length(MM3_2_summary$system_lengths) - 1)) # Mean system length (number of customers in system) ## [,1] ## [1,] 6.152826 MM3_2_summary$mwt # Mean waiting time
## [1] 3.640104
MM3_2_summary\$mrt # Mean response time
## [1] 6.141271

# Bibliography

Thomopoulos, N (2012). Fundamentals of Queuing Systems: Statistical Methods for Analyzing Queuing Models. Springer Science & Business Media