General Queueing Code in R

[swaheera]

Not exactly related to Simmer directly, but I wanted to show some generic queueing code I wrote.

This code (below) statistically simulates the progression of a MMK queue at different time points. It does NOT track individual customers at each time point directly:

# Load required libraries
library(ggplot2)
library(tidyr)
library(dplyr)
library(gridExtra)

# Set simulation parameters
lambda <- 7          # Arrival rate
mu <- 1.5               # 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

# Define the number of servers to simulate
servers <- c(2,3,4)  

# Function to run a 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))
}

# Run multiple simulations for each server count
results <- lapply(servers, function(s) {
    replicate(num_simulations, run_simulation(s), simplify = FALSE)
})

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

# Create data frames for all scenarios
plot_data <- lapply(seq_along(servers), function(i) {
    create_plot_data(results[[i]], servers[i])
})

# Combine data for plotting
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))

# Function to calculate average values
calculate_average <- function(data, value_column) {
    data %>%
        group_by(Time, Servers) %>%
        summarise(Average = mean(!!sym(value_column)), .groups = "drop")
}

# Function to create a single plot
create_plot <- function(data, value_column, title, y_label) {
    avg_data <- calculate_average(data, value_column)
    
    subtitle <- sprintf("Initial Queue Size: %d, Arrival Rate: %.1f, Service Rate: %.1f\nSimulations: %d, Simulation Time: %d",
                        initial_queue_size, lambda, mu, num_simulations, sim_time)
    
    # Generate a color palette for the number of server scenarios
    color_palette <- rainbow(length(servers))
    
    ggplot() +
        geom_line(data = data, aes(x = Time, y = !!sym(value_column), group = interaction(Simulation, Servers), color = factor(Servers)), alpha = 0.1) +
        geom_line(data = avg_data, aes(x = Time, y = Average, color = factor(Servers)), size = 1) +
        scale_color_manual(values = setNames(color_palette, servers), name = "Servers") +
        labs(title = title, subtitle = subtitle, x = "Time", y = y_label) +
        theme_minimal() +
        theme(legend.key = element_rect(fill = "white"),
              legend.background = element_rect(fill = "white", color = "black"),
              legend.position = "bottom",
              plot.subtitle = element_text(size = 8, face = "italic"))
}

# Create the 3 plots
plot1 <- create_plot(plot_data_queue, "QueueLength", "Queue Length Over Time", "Queue Length")
plot2 <- create_plot(plot_data_processed, "ProcessedOrders", "Cumulative Orders Processed Over Time", "Orders Processed")
plot3 <- create_plot(plot_data_wait, "LongWaitPercent", paste("% of People Waiting >", k_minutes, "Minutes Over Time"), "% Long Wait")

# Arrange plots side by side
grid.arrange(plot1, plot2, plot3, ncol = 3)

Maybe somebody might find this useful?

[swaheera]

I have also modified this to work in continuous time as well: in case anyone is interested

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

# Simulation parameters
lambda <- 7          # Arrival rate
mu <- 1.5            # Service rate
sim_time <- 1000     # Simulation time
k_minutes <- 15      # Threshold for waiting time
num_simulations <- 100  # Number of simulations to run
initial_queue <- 1 # Initial queue length

# Define the number of servers to simulate
servers <- c(2, 3, 4)  

# Function to run a single simulation
run_simulation <- function(num_servers) {
    queue <- initial_queue
    processed <- 0
    waiting_times <- numeric(0)
    queue_length <- numeric(sim_time)
    processed_over_time <- numeric(sim_time)
    long_wait_percent <- numeric(sim_time)
    
    # Generate inter-arrival times
    arrival_times <- cumsum(rexp(sim_time * lambda, rate = lambda))
    
    # Generate service times for each arrival (including initial queue)
    service_times <- rexp(length(arrival_times) + initial_queue, rate = mu)
    
    # Initialize servers
    server_available_times <- rep(0, num_servers)
    
    for (t in 1:sim_time) {
        # Process arrivals
        new_arrivals <- sum(arrival_times <= t & arrival_times > (t-1))
        queue <- queue + new_arrivals
        
        # Process departures
        for (s in 1:num_servers) {
            if (server_available_times[s] <= t && queue > 0) {
                queue <- queue - 1
                processed <- processed + 1
                service_time <- service_times[processed]
                server_available_times[s] <- t + service_time
                
                # Calculate waiting time (assume initial queue members arrived at t=0)
                if (processed <= initial_queue) {
                    waiting_time <- t
                } else {
                    arrival_time <- arrival_times[processed - initial_queue]
                    waiting_time <- max(0, t - arrival_time)
                }
                waiting_times <- c(waiting_times, waiting_time)
            }
        }
        
        # 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))
}

# Run multiple simulations for each server count
results <- lapply(servers, function(s) {
    replicate(num_simulations, run_simulation(s), simplify = FALSE)
})

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

# Create data frames for all scenarios
plot_data <- lapply(seq_along(servers), function(i) {
    create_plot_data(results[[i]], servers[i])
})

# Combine data for plotting
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))

# Function to calculate average values
calculate_average <- function(data, value_column) {
    data %>%
        group_by(Time, Servers) %>%
        summarise(Average = mean(!!sym(value_column)), .groups = "drop")
}

# Function to create a single plot
create_plot <- function(data, value_column, title, y_label) {
    avg_data <- calculate_average(data, value_column)
    
    subtitle <- sprintf("Arrival Rate: %.1f, Service Rate: %.1f\nSimulations: %d, Simulation Time: %d, Initial Queue: %d",
                        lambda, mu, num_simulations, sim_time, initial_queue)
    
    # Generate a color palette for the number of server scenarios
    color_palette <- rainbow(length(servers))
    
    ggplot() +
        geom_line(data = data, aes(x = Time, y = !!sym(value_column), group = interaction(Simulation, Servers), color = factor(Servers)), alpha = 0.1) +
        geom_line(data = avg_data, aes(x = Time, y = Average, color = factor(Servers)), size = 1) +
        scale_color_manual(values = setNames(color_palette, servers), name = "Servers") +
        labs(title = title, subtitle = subtitle, x = "Time", y = y_label) +
        theme_minimal() +
        theme(legend.key = element_rect(fill = "white"),
              legend.background = element_rect(fill = "white", color = "black"),
              legend.position = "bottom",
              plot.subtitle = element_text(size = 8, face = "italic"))
}

# Create the 3 plots
plot1 <- create_plot(plot_data_queue, "QueueLength", "Queue Length Over Time", "Queue Length")
plot2 <- create_plot(plot_data_processed, "ProcessedOrders", "Cumulative Orders Processed Over Time", "Orders Processed")
plot3 <- create_plot(plot_data_wait, "LongWaitPercent", paste("% of People Waiting >", k_minutes, "Minutes Over Time"), "% Long Wait")

# Arrange plots side by side
grid.arrange(plot1, plot2, plot3, ncol = 3)

Please note that neither of these codes have been formally peer reviewed and are simply draft ideas/proof of concept codes that I have been working on.

[swaheera]

In any case, I think we should all take a moment to deeply thank Dr. Ucar @Enchufa2 for all his creativity, efforts, dedication, accomplishments and kindness throughout the development of the simmer R package. It is a huge accomplishment and I am confident that many students (like myself) in the future will greatly benefit from his hard work, intelligence and dedication. Thank you Dr. Ucar!