# Simulating a Queue in R [closed]

A few months ago, I learned how to simulate a simple queue in R https://stackoverflow.com/questions/77734912/r-how-to-simulate-a-queue . Now, I am trying to a make a more complicated version of this queue with the following parameters:

lambda <- 5  # Arrival rate
mu <- 1      # 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 <- 100  # Initial queue size
time_step <- 1  # Time step for discretization


I simulated this code in R for 3 Servers vs 4 Servers:

      library(dplyr)
library(ggplot2)
library(tidyr)
library(purrr)
library(gridExtra)

lambda <- 5
mu <- 1
sim_time <- 200
k_minutes <- 15
num_simulations <- 100
initial_queue_size <- 100
time_step <- 1
k_values <- c(3, 4)

run_simulation <- function(seed, k) {
set.seed(seed)

events <- data.frame(
time = c(0, cumsum(rexp(ceiling(sim_time * lambda), rate = lambda))),
type = "arrival"
)
events <- events[events$time <= sim_time, ] queue <- numeric(initial_queue_size) # Initialize queue with initial_queue_size servers <- numeric(k) processed <- 0 waiting_times <- numeric() results <- data.frame( time = seq(0, sim_time, by = time_step), queue_length = initial_queue_size, processed_orders = 0, waiting_longer = 0, total_arrivals = initial_queue_size ) event_index <- 1 for (i in 1:nrow(results)) { current_time <- results$time[i]

# Process events up to current time
while (event_index <= nrow(events) && events$time[event_index] <= current_time) { event_time <- events$time[event_index]

# Process completed services
finished <- servers <= event_time
if (any(finished)) {
processed <- processed + sum(finished)
servers[finished] <- 0
}

# Process new arrival
results$total_arrivals[i] <- results$total_arrivals[i] + 1
if (any(servers == 0)) {
free_server <- which(servers == 0)[1]
servers[free_server] <- event_time + rexp(1, mu)
waiting_times <- c(waiting_times, 0)
} else {
queue <- c(queue, event_time)
}

# Update queue
while (length(queue) > 0 && any(servers == 0)) {
free_server <- which(servers == 0)[1]
wait_time <- event_time - queue[1]
waiting_times <- c(waiting_times, wait_time)
servers[free_server] <- event_time + rexp(1, mu)
queue <- queue[-1]
}

event_index <- event_index + 1
}

results$queue_length[i] <- length(queue) results$processed_orders[i] <- processed
results\$waiting_longer[i] <- sum(waiting_times > k_minutes)
}

results
}

run_simulations <- function(k_values) {
map(k_values, function(k) {
map(1:num_simulations, ~run_simulation(., k)) %>%
set_names(paste0("sim_", 1:num_simulations))
}) %>% set_names(paste0("k", k_values))
}

simulations <- run_simulations(k_values)

process_results <- function(simulations) {
map_dfr(names(simulations), function(k_name) {
k <- as.integer(gsub("k", "", k_name))
bind_rows(simulations[[k_name]], .id = "simulation") %>%
mutate(k = k, simulation = as.integer(gsub("sim_", "", simulation))) %>%
group_by(simulation, k) %>%
mutate(
cumulative_waiting_longer = cumsum(waiting_longer),
cumulative_total_arrivals = cumsum(total_arrivals),
waiting_percentage = pmin(100, pmax(0, (cumulative_waiting_longer / cumulative_total_arrivals) * 100))
) %>%
ungroup()
})
}

all_results <- process_results(simulations)

plot_waiting_percentage <- function(data, k) {
ggplot(data %>% filter(k == !!k), aes(x = time, y = waiting_percentage, group = simulation)) +
geom_line(alpha = 0.1, color = "blue") +
stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red", size = 1) +
labs(title = paste("Percentage of People Waiting >", k_minutes, "Minutes (k=", k, ")"),
subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu),
x = "Time", y = "Percentage") +
theme_minimal() +
ylim(0, 100)
}

plot_queue_length <- function(data, k) {
ggplot(data %>% filter(k == !!k), aes(x = time, y = queue_length, group = simulation)) +
geom_line(alpha = 0.1, color = "blue") +
stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red", size = 1) +
labs(title = paste("Queue Length Over Time (k=", k, ")"),
subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu, ", Initial Queue Size =", initial_queue_size),
x = "Time", y = "Queue Length") +
theme_minimal() +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA))  # Start y-axis from 0
}

plot_cumulative_orders <- function(data, k) {
ggplot(data %>% filter(k == !!k), aes(x = time, y = processed_orders, group = simulation)) +
geom_line(alpha = 0.1, color = "blue") +
stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red", size = 1) +
labs(title = paste("Cumulative Orders Processed (k=", k, ")"),
subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu, ", Initial Queue Size =", initial_queue_size),
x = "Time", y = "Cumulative Orders") +
theme_minimal() +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA))
}

plots <- map(k_values, function(k) {
list(
waiting_percentage = plot_waiting_percentage(all_results, k),
queue_length = plot_queue_length(all_results, k),
cumulative_orders = plot_cumulative_orders(all_results, k)
)
})

do.call(grid.arrange, c(unlist(plots, recursive = FALSE), ncol = 2))


Based on these results, we can see that on average, the same queue with 4 servers outperforms the 3 server queue for cumulative orders processed and queue length - but somehow the percent of customers waiting longer than 15 minutes is better (i.e. increases slower) for 3 servers than 4 servers?

• Ahoy! "Can someone please help me review and find where the mistake is" Does it seem like the mistake is in the simulation or in the code? Commented Sep 6 at 21:30
• @ sam: thank you for your reply! I think its in the code because the results of the simulation are illogical! Commented Sep 6 at 21:42
• Okay- Thank you for your honest reply. CR requires that the code be working correctly, to the best of the author's knowledge, before proceeding with a review. Until it works as expected it will be closed. Commented Sep 6 at 23:39
• thanks! I made a different version of this question codereview.stackexchange.com/questions/293622/… can you please check it out? Commented Sep 7 at 6:16