Adding a custom arrival rate vector to simmer

[swaheera]

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:

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)

[Enchufa2]

Yes, it is certainly possible. You pass a function to the generator, so that function could return anything you need (as long as the numbers represent interarrival times). For example:

env %>% 
  add_generator("customer", customer, function() c(my_arrival_vector, -1))

The -1 makes the generator stop. Another option is to use absolute times instead:

env %>% 
  add_generator("customer", customer, at(cumsum(my_arrival_vector))

Note the cumsum(), which assumes that my_arrival_vector contains interarrival times, and the at() helper. There are other helpers that may be useful for other tasks.

Note: next time, please open this kind of questions in the Discussions section instead of the Issues one. I've moved this one there.

[swaheera]

@Enchufa2 : Thank you so much Dr. Ucar . In the future, I will post directly in the discussions area.

[swaheera]

As I understand, in this example I posted , the arrivals appear in single "bursts" but the services remain classical. For example, suppose 100 people appear on the first day and 83 are processed. On the second day, suppose 20 people appear - this means that we have 17 people from the first day and the 20 people from day 2 in the current queue. On the second day, we will start with the customer that has been waiting the longest (i.e. the 84th customer from the first day). Is this what the R code is doing?

[Enchufa2]

A vector of interarrivals like c(1, 1, 1, 1) will put an arrival at t = 1, 2, 3, 4, no matter how you provide them (i.e., one by one or all at once in a single call). A burst would be c(1, 0, 0, 0). This would put 4 arrivals at t = 1.

[swaheera]

@Enchufa2 Hi Dr. Ucar . Thank you so much for your reply. Could you please show me how to do this in R?

As an example, here is the R code I am using for a standard m/m/k queue that compares the queue length for 3 servers vs 4 servers ... and the percent of people that have been waiting longer than 15 minutes :

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

# [Keep the existing simulation parameters and run_simulation function]

# [Keep the existing simulation runs for k3 and k4]

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$time <- floor(resources$time / time_step) * 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)
}

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

# Combine results
combined_queue_length <- rbind(results_k3$queue_length_data, results_k4$queue_length_data)
combined_waiting_percentage <- rbind(results_k3$waiting_percentage_data, results_k4$waiting_percentage_data)

# Updated plotting functions
plot_combined_queue_length <- function(data) {
    ggplot(data, aes(x = time, y = queue, group = interaction(simulation, k), color = factor(k))) +
        geom_line(alpha = 0.1) +
        stat_summary(aes(group = k), fun = mean, geom = "line", size = 1) +
        scale_color_manual(values = c("3" = "blue", "4" = "red")) +
        labs(title = "Queue Length Over Time",
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu),
             x = "Time", y = "Queue Length", color = "Number of Servers (k)") +
        theme_minimal()
}

plot_combined_waiting_percentage <- function(data) {
    ggplot(data, aes(x = time, y = percentage, group = interaction(simulation, k), color = factor(k))) +
        geom_line(alpha = 0.1) +
        stat_summary(aes(group = k), fun = mean, geom = "line", size = 1) +
        scale_color_manual(values = c("3" = "blue", "4" = "red")) +
        labs(title = paste("Percentage of People Waiting >", k_minutes, "Minutes"),
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu),
             x = "Time", y = "Percentage", color = "Number of Servers (k)") +
        theme_minimal() +
        ylim(0, 100)
}

# Create combined plots
combined_queue_length_plot <- plot_combined_queue_length(combined_queue_length)
combined_waiting_percentage_plot <- plot_combined_waiting_percentage(combined_waiting_percentage)

# Arrange and display plots
grid.arrange(combined_queue_length_plot, combined_waiting_percentage_plot, ncol = 2)

[swaheera]

(Assuming the above code is correct), could you please show me how to modify the above code to correspond to a situation where the arrivals are predetermined? E.g..

arrivals = as.integer(rnorm(1000, 100,5))

Thank you so much Dr. Ucar.

[Enchufa2]

What do you mean? You mean to modify the simulation to introduce predetermined arrivals? That was shown in my first comment. Or you mean to modify the processing? The processing shouldn't change.

[swaheera]

Hi Dr. Ucar,

My apologies for the confusion. I tried to use your advice on how to input a custom arrivals vector, but I am still confused on whether I am doing this correctly.

Here is the code I am using:

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

# Simulation parameters
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)

# Custom arrival rate vector
set.seed(123)  # For reproducibility
my_arrival_vector <- as.integer(rnorm(sim_time, 500, 5))

# 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, at(cumsum(my_arrival_vector)))
   
    # 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 k = 3 and k = 4
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))

# 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
        resources$time <- floor(resources$time / time_step) * 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 k = 3 and k = 4
results_k3 <- process_results(simulation_results_k3, 3)
results_k4 <- process_results(simulation_results_k4, 4)

# 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))
        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 k = 3 and k = 4
cumulative_orders_k3 <- process_cumulative_orders(simulation_results_k3, 3)
cumulative_orders_k4 <- process_cumulative_orders(simulation_results_k4, 4)

# 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 = "Custom Arrival Rate, Service Rate = 1",
             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("Custom Arrival Rate, Service Rate = 1, Initial Queue Size =", initial_queue_size),
             x = "Time", y = "Queue Length") +
        theme_minimal() +
        scale_y_continuous(expand = c(0, 0), limits = c(0, NA))
}

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("Custom Arrival Rate, Service Rate = 1, Initial Queue Size =", initial_queue_size),
             x = "Time", y = "Cumulative Orders") +
        theme_minimal()
}

# Create plots
waiting_percentage_plot_k3 <- plot_waiting_percentage(results_k3$waiting_percentage_data, 3)
waiting_percentage_plot_k4 <- plot_waiting_percentage(results_k4$waiting_percentage_data, 4)
queue_length_plot_k3 <- plot_queue_length(results_k3$queue_length_data, 3)
queue_length_plot_k4 <- plot_queue_length(results_k4$queue_length_data, 4)
cumulative_orders_plot_k3 <- plot_cumulative_orders(cumulative_orders_k3, 3)
cumulative_orders_plot_k4 <- plot_cumulative_orders(cumulative_orders_k4, 4)

# Arrange and display all plots
grid.arrange(waiting_percentage_plot_k3, waiting_percentage_plot_k4,
             queue_length_plot_k3, queue_length_plot_k4,
             cumulative_orders_plot_k3, cumulative_orders_plot_k4,
             ncol = 2)

The above code runs, but the % plots are empty , which makes me think I am not using the arrivals vector correctly as you had mentioned. Is this the correct way to use the arrivals vector? Once again, my apologies for the confusion.

[Enchufa2]

You are using it correctly, but you have arrival times that are higher than your simulation time, so you don't see any arrivals:

> set.seed(123)  # For reproducibility
> my_arrival_vector <- as.integer(rnorm(sim_time, 500, 5))
> cumsum(my_arrival_vector)
  [1]   497   995  1502  2002  2502  3010  3512  4005  4501  4998  5504  6005  6507  7007  7504  8012  8514  9004  9507 10004 10498 10996
 [23] 11490 11986 12482 12973 13477 13977 14471 14977 15479 15977 16481 16985 17489 17992 18494 18993 19491 19989 20485 20983 21476 21986
 [45] 22492 22986 23483 23980 24483 24982 25483 25982 26481 26987 27485 27992 28484 28986 29486 29987 30488 30985 31483 31977 32471 32972
 [67] 33474 33974 34478 34988 35485 35973 36478 36974 37470 37975 38473 38966 39466 39965 40465 40966 41464 41967 42465 42966 43471 43973
 [89] 44471 44976 45480 45982 46483 46979 47485 47981 48491 48998 49496 49990 50486 50987 51485 51983 52478 52977 53473 53964 54462 54966
[111] 55463 55966 56457 56956 57458 57959 58459 58955 59450 59944 60444 60939 61436 61934 62443 62939 63440 63940 64435 64934 65441 65943
[133] 66443 66940 67429 67934 68426 68929 69438 69930 70433 70931 71423 71915 72406 72903 73395 73898 74408 74901 75404 75907 76408 76902
[155] 77401 77899 78401 78899 79403 79901 80406 80900 81393 81909 82406 82907 83410 83907 84409 84910 85408 85908 86407 86917 87413 87907
[177] 88407 88908 89410 89907 90401 90907 91405 91900 92398 92897 93402 93902 94405 94902 95403 95901 96401 96896 97389 97898 98401 98894
[199] 99390 99884

but sim_time is only 200.