r-simmer / simmer

Discrete-Event Simulation for R
https://r-simmer.org
GNU General Public License v2.0
223 stars 42 forks source link

Adding a custom arrival rate vector to simmer #311

Closed swaheera closed 3 months ago

swaheera commented 3 months ago

I made this queue simulation that shows how the percent of all customers that waited more than 15 minutes changes when more servers are added to the system:

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

# Arrival rate
lambda <- 5  
# Service rate
mu <- 1      
# Simulation time
sim_time <- 200  
# Threshold for waiting time
k_minutes <- 15  
# Number of simulations to run
num_simulations <- 100  
# Initial queue size
initial_queue_size <- 10  

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

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

simulation_results_k3 <- lapply(1:num_simulations, function(i) run_simulation(i, 3))
simulation_results_k4 <- lapply(1:num_simulations, function(i) run_simulation(i, 4))

process_results <- function(simulation_results, k_value) {
    queue_length_data <- do.call(rbind, lapply(1:num_simulations, function(i) {
        resources <- simulation_results[[i]]$resources
        resources$simulation <- i
        resources$k <- k_value
        resources
    }))

    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$k <- k_value
        arrivals
    }))

    waiting_percentage_data <- waiting_time_data %>%
        group_by(simulation, time = floor(start_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)
}

results_k3 <- process_results(simulation_results_k3, 3)
results_k4 <- process_results(simulation_results_k4, 4)

I then plotted the results in ggplot (code not shown, but I can add if someone is interested). The plots look like this - confirming that when more servers are added, the percent of customers that waited longer than 15 minutes increases at a lesser rate:

image

My Question: Suppose if instead of a constant arrival rate, I have a custom arrival rate vector:

my_arrival_vector = as.integer((rnorm(200, 500, 5)))

That is, for each value of sim_time 1,2,3,4,... we insert my_arrival_vector[1], my_arrival_vector[2], my_arrival_vector[3] etc number of new people into the queue.

Is it possible to feed this vector into the simulation directly instead?

Note: Plotting Code

combined_waiting_data <- rbind(results_k3$waiting_percentage_data, results_k4$waiting_percentage_data)

waiting_percentage_plot_k3 <- ggplot(results_k3$waiting_percentage_data, aes(x = time, y = percentage, group = interaction(simulation, k), color = factor(k))) +
    geom_line(alpha = 0.7) +
    stat_summary(fun = mean, geom = "line", aes(group = 1), color = "black", size = 1) +
    labs(title = paste("Percentage of People Waiting >", k_minutes, "Minutes (k=3)"),
         subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu),
         x = "Time", y = "Percentage", color = "Number of Servers (k)") +
    theme_minimal()

waiting_percentage_plot_k4 <- ggplot(results_k4$waiting_percentage_data, aes(x = time, y = percentage, group = interaction(simulation, k), color = factor(k))) +
    geom_line(alpha = 0.7) +
    stat_summary(fun = mean, geom = "line", aes(group = 1), color = "black", size = 1) +
    labs(title = paste("Percentage of People Waiting >", k_minutes, "Minutes (k=4)"),
         subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu),
         x = "Time", y = "Percentage", color = "Number of Servers (k)") +
    theme_minimal()

grid.arrange(waiting_percentage_plot_k3, waiting_percentage_plot_k4, ncol = 2)