Best Practices for Simulating Queues with Large Parameters?

[swaheera]

I find for queues with small parameters (e.g. initial queue size, service rate, arrival rate), the results of the simulations take much shorter to run and the results seem more logical (e.g. smaller queue sizes for more servers vs fewer servers). I would have thought that the size of the parameters should not affect the time required to run the simulation? However, when looking at the time/results of the simulation involving queues with larger parameters, such issues still persist.

For example, here is a MMK queue with large parameters : library(simmer)

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

# Simulation parameters
lambda <- 100  # Arrival rate
mu <- 20      # 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 <- 100000  # Initial queue size
time_step <- 1  # Time step for discretization (e.g., 1 second)
k1 <- 90  # First k value
k2 <- 100  # Second k value

# Function to run a single simulation
run_simulation <- function(seed, k) {
    set.seed(seed)
    env <- simmer("MMK_Queue")
    
    customer <- trajectory() %>%
        set_attribute("arrival_time", function() {now(env)}) %>%
        seize("server", 1) %>%
        timeout(function() {rexp(1, mu)}) %>%
        release("server", 1)
    
    env %>%
        add_resource("server", capacity = k) %>%
        add_generator("customer", customer, function() {rexp(1, lambda)})
    
    
    # Add initial customers to the queue
    env %>% add_generator("initial_customer_", customer, at(rep(0, initial_queue_size)))
    
    env %>% run(until = sim_time)
    
    resources <- get_mon_resources(env)
    arrivals <- get_mon_arrivals(env)
    
    list(resources = resources, arrivals = arrivals)
}

# Running simulations for k1 and k2
simulation_results_k1 <- lapply(1:num_simulations, function(i) run_simulation(i, k1))
simulation_results_k2 <- lapply(1:num_simulations, function(i) run_simulation(i, k2))

# Function to process simulation results
process_results <- function(simulation_results, k_value) {
    queue_length_data <- do.call(rbind, lapply(1:num_simulations, function(i) {
        resources <- simulation_results[[i]]$resources
        # Discretize time points
        resources$time <- floor(resources$time / time_step) * time_step
        # Group by discretized time and take the last queue length in each time step
        resources %>%
            group_by(time) %>%
            summarise(queue = last(queue)) %>%
            mutate(simulation = i, k = k_value)
    }))
    
    waiting_time_data <- do.call(rbind, lapply(1:num_simulations, function(i) {
        arrivals <- simulation_results[[i]]$arrivals
        arrivals$simulation <- i
        arrivals$waiting_time <- arrivals$end_time - arrivals$start_time - arrivals$activity_time
        arrivals$k <- k_value
        arrivals
    }))
    
    waiting_percentage_data <- waiting_time_data %>%
        mutate(time = floor(start_time / time_step) * time_step) %>%
        group_by(simulation, time, k) %>%
        summarise(
            total = n(),
            waiting_longer = sum(waiting_time > k_minutes),
            percentage = waiting_longer / total * 100,
            .groups = 'drop'
        )
    
    list(queue_length_data = queue_length_data, waiting_percentage_data = waiting_percentage_data)
}

# Process the results for k1 and k2
results_k1 <- process_results(simulation_results_k1, k1)
results_k2 <- process_results(simulation_results_k2, k2)

# Function to calculate cumulative orders processed
process_cumulative_orders <- function(simulation_results, k_value) {
    cumulative_orders_data <- do.call(rbind, lapply(1:num_simulations, function(i) {
        arrivals <- simulation_results[[i]]$arrivals
        arrivals$simulation <- i
        arrivals$cumulative_orders <- seq_len(nrow(arrivals))  # Cumulative count of processed orders
        arrivals$k <- k_value
        arrivals %>%
            mutate(time = floor(end_time / time_step) * time_step) %>%
            group_by(time, k, simulation) %>%
            summarise(cumulative_orders = max(cumulative_orders), .groups = 'drop')
    }))
    cumulative_orders_data
}

# Process cumulative orders for k1 and k2
cumulative_orders_k1 <- process_cumulative_orders(simulation_results_k1, k1)
cumulative_orders_k2 <- process_cumulative_orders(simulation_results_k2, k2)

# Plotting functions
plot_waiting_percentage <- function(data, k) {
    ggplot(data, aes(x = time, y = 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, aes(x = time, y = queue, 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, aes(x = time, y = cumulative_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()
}

# Create plots
waiting_percentage_plot_k1 <- plot_waiting_percentage(results_k1$waiting_percentage_data, k1)
waiting_percentage_plot_k2 <- plot_waiting_percentage(results_k2$waiting_percentage_data, k2)
queue_length_plot_k1 <- plot_queue_length(results_k1$queue_length_data, k1)
queue_length_plot_k2 <- plot_queue_length(results_k2$queue_length_data, k2)
cumulative_orders_plot_k1 <- plot_cumulative_orders(cumulative_orders_k1, k1)
cumulative_orders_plot_k2 <- plot_cumulative_orders(cumulative_orders_k2, k2)

# Arrange and display all plots
grid.arrange(waiting_percentage_plot_k1, waiting_percentage_plot_k2,
             queue_length_plot_k1, queue_length_plot_k2,
             cumulative_orders_plot_k1, cumulative_orders_plot_k2,
             ncol = 2)

Even though the length of each simulation and the total number of simulations is relatively small, this simulation is still taking a long time to run.

For such situations (i.e. large queue parameters), are there any best practices for Simmer? Could someone please comment on my code for this situation?

Thanks!

[Enchufa2]

There's not much to say/do here: the more arrivals you simulate, the more events are generated, and more time is require to complete the simulation. The simulation time grows linearly with the number of events that happen. Best practices? Keep the number of events (i.e. complexity of the trajectories) at the minimum required for your purposes.