Example : Basic Simmer Simulation

[swaheera]

Hello Everyone,

I just wanted to share this Simmer simulation I made as a reference in case anyone wanted to use this as an example.

Thanks

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

# Simulation 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 <- 50  # Initial queue size
time_step <- 1  # Time step for discretization (e.g., 1 second)
k1 <- 3  # First k value
k2 <- 4  # 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
    for (i in 1:initial_queue_size) {
        env %>% add_generator(paste0("initial_customer_", i), customer, at(0))
    }
   
    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)

[Enchufa2]

Thanks. One detail: instead of

    # Add initial customers to the queue
    for (i in 1:initial_queue_size) {
        env %>% add_generator(paste0("initial_customer_", i), customer, at(0))
    }

which works fine, it's easier to do this:

    # Add initial customers to the queue
    env %>% add_generator("initial_customer_", customer, at(rep(0, initial_queue_size)))

[swaheera]

Hi Dr. Ucar, good idea! I made the recommended change:

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

# Simulation 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 <- 50  # Initial queue size
time_step <- 1  # Time step for discretization (e.g., 1 second)
k1 <- 3  # First k value
k2 <- 4  # 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)