I am trying to a M/M/K queue (M/M/c queue - Wikipedia) simulation in R with an arrival rate of 8, service rate of 10 and 1 server. The average queue length at steady state according to the theoretical formula should be rho/(1-rho) where rho = (lambda/mu). In my case, this should result in average queue length of 4.
I tried to do this with an R simulation.
First, I defined the queue parameters:
set.seed(123)
library(ggplot2)
library(tidyr)
library(dplyr)
library(gridExtra)
# simulation parameters
lambda <- 8 # Arrival rate
mu <- 10 # Service rate
sim_time <- 200 # Simulation time
k_minutes <- 15 # Threshold for waiting time
num_simulations <- 100 # Number of simulations to run
initial_queue_size <- 0 # Initial queue size
time_step <- 1 # Time step for discretization
servers <- c(1)
Next, I defined a function perform a single simulation. My approach takes the current queue length, adds random arrivals and subtracts random departures - and then repeats this process:
# single simulation
run_simulation <- function(num_servers) {
queue <- initial_queue_size
processed <- 0
waiting_times <- numeric(0)
queue_length <- numeric(sim_time)
processed_over_time <- numeric(sim_time)
long_wait_percent <- numeric(sim_time)
for (t in 1:sim_time) {
# Process arrivals
arrivals <- rpois(1, lambda * time_step)
queue <- queue + arrivals
# Process departures
departures <- min(queue, rpois(1, num_servers * mu * time_step))
queue <- queue - departures
processed <- processed + departures
# Update waiting times
if (length(waiting_times) > 0) {
waiting_times <- waiting_times + time_step
}
if (arrivals > 0) {
waiting_times <- c(waiting_times, rep(0, arrivals))
}
if (departures > 0) {
waiting_times <- waiting_times[-(1:departures)]
}
# Record metrics
queue_length[t] <- queue
processed_over_time[t] <- processed
long_wait_percent[t] <- ifelse(length(waiting_times) > 0,
sum(waiting_times > k_minutes) / length(waiting_times) * 100,
0)
}
return(list(queue_length = queue_length,
processed_over_time = processed_over_time,
long_wait_percent = long_wait_percent))
}
I then run this simulation:
results <- lapply(servers, function(s) {
replicate(num_simulations, run_simulation(s), simplify = FALSE)
})
And finally, I tidy everything up into data frames:
# Function to create data frames for plotting
create_plot_data <- function(results, num_servers) {
plot_data_queue <- data.frame(
Time = rep(1:sim_time, num_simulations),
QueueLength = unlist(lapply(results, function(x) x$queue_length)),
Simulation = rep(1:num_simulations, each = sim_time),
Servers = num_servers
)
plot_data_processed <- data.frame(
Time = rep(1:sim_time, num_simulations),
ProcessedOrders = unlist(lapply(results, function(x) x$processed_over_time)),
Simulation = rep(1:num_simulations, each = sim_time),
Servers = num_servers
)
plot_data_wait <- data.frame(
Time = rep(1:sim_time, num_simulations),
LongWaitPercent = unlist(lapply(results, function(x) x$long_wait_percent)),
Simulation = rep(1:num_simulations, each = sim_time),
Servers = num_servers
)
return(list(queue = plot_data_queue, processed = plot_data_processed, wait = plot_data_wait))
}
plot_data <- lapply(seq_along(servers), function(i) {
create_plot_data(results[[i]], servers[i])
})
plot_data_queue <- do.call(rbind, lapply(plot_data, function(x) x$queue))
plot_data_processed <- do.call(rbind, lapply(plot_data, function(x) x$processed))
plot_data_wait <- do.call(rbind, lapply(plot_data, function(x) x$wait))
My Problem: When I calculate the average queue length, I get the following:
> mean(plot_data_queue$QueueLength)
[1] 2.46215
And this average does not match the theoretical answer.
Can someone please help me understand what is lacking in my approach and what I can to do to fix this?
Thanks!