A queue where the number of servers change

[swaheera]

This is a proof-of-concept idea I thought of: Suppose there is a MMK queue where :

• When there are less than 300 people in line, the queue has 3 servers
• When there are more than 300 people in line, the queue now has 4 servers

I was trying to see online if it is possible to do this in Simmer. I tried to modify previous Simmer code I was using to do this (i.e. simulate queue trajectories and plot different metrics):

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

# 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)
queue_threshold <- 300  # Threshold for changing number of servers

# Function to run a single simulation
run_simulation <- function(seed) {
    cat("Running simulation", seed, "of", num_simulations, "\n")
    
    set.seed(seed)
    env <- simmer("Dynamic_MMK_Queue")
    
    # Create a data frame to store capacity changes
    capacity_changes <- data.frame(time = numeric(), capacity = numeric())
    
    customer <- trajectory() %>%
        set_attribute("arrival_time", function() {now(env)}) %>%
        seize("server", 1) %>%
        timeout(function() {rexp(1, mu)}) %>%
        release("server", 1)
    
    # Function to update server capacity based on system size
    update_capacity <- trajectory() %>%
        timeout(time_step) %>%
        log_(function() {
            system_size <- get_queue_count(env, "server") + get_server_count(env, "server")
            current_capacity <- get_server_count(env, "server")
            if (system_size >= queue_threshold && current_capacity == 3) {
                env %>% add_resource("server", 1)
                capacity_changes <<- rbind(capacity_changes, data.frame(time = now(env), capacity = 4))
                return(paste("Time:", now(env), "System size:", system_size, "Capacity increased to 4"))
            } else if (system_size < queue_threshold && current_capacity == 4) {
                env %>% update_resource("server", -1)
                capacity_changes <<- rbind(capacity_changes, data.frame(time = now(env), capacity = 3))
                return(paste("Time:", now(env), "System size:", system_size, "Capacity decreased to 3"))
            }
            return(paste("Time:", now(env), "System size:", system_size, "Capacity unchanged:", current_capacity))
        }) %>%
        rollback(2)
    
    env %>%
        add_resource("server", capacity = 3) %>%
        add_generator("customer", customer, function() {rexp(1, lambda)}) %>%
        add_generator("capacity_updater", update_capacity, at(0))
    
    # Add initial customers to the queue
    for (i in 1:initial_queue_size) {
        env %>% add_generator(paste0("initial_customer_", i), customer, at(0))
    }
    
    # Record initial capacity
    capacity_changes <- rbind(capacity_changes, data.frame(time = 0, capacity = 3))
    
    env %>% run(until = sim_time)
    
    resources <- get_mon_resources(env)
    arrivals <- get_mon_arrivals(env)
    
    list(resources = resources, arrivals = arrivals, capacity_changes = capacity_changes)
}

# Running simulations
simulation_results <- lapply(1:num_simulations, run_simulation)

# Function to process simulation results
process_results <- function(simulation_results) {
    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), system = last(system), server = last(server)) %>%
            mutate(simulation = i)
    }))
    
    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
    }))
    
    waiting_percentage_data <- waiting_time_data %>%
        filter(end_time > 0) %>%  # Exclude customers who haven't finished service
        mutate(time = floor(end_time / time_step) * time_step) %>%
        group_by(simulation, time) %>%
        summarise(
            total = n(),
            waiting_longer = sum(waiting_time > k_minutes),
            percentage = waiting_longer / total * 100,
            .groups = 'drop'
        )
    
    server_capacity_data <- do.call(rbind, lapply(1:num_simulations, function(i) {
        capacity_changes <- simulation_results[[i]]$capacity_changes
        data.frame(
            time = seq(0, sim_time, by = time_step),
            capacity = approx(capacity_changes$time, capacity_changes$capacity, xout = seq(0, sim_time, by = time_step), method = "constant", rule = 2)$y,
            simulation = i
        )
    }))
    
    list(queue_length_data = queue_length_data, 
         waiting_percentage_data = waiting_percentage_data,
         server_capacity_data = server_capacity_data)
}

# Process the results
results <- process_results(simulation_results)

# Function to calculate cumulative orders processed
process_cumulative_orders <- function(simulation_results) {
    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 %>%
            mutate(time = floor(end_time / time_step) * time_step) %>%
            group_by(time, simulation) %>%
            summarise(cumulative_orders = max(cumulative_orders), .groups = 'drop')
    }))
    cumulative_orders_data
}

# Process cumulative orders
cumulative_orders <- process_cumulative_orders(simulation_results)

# Plotting functions
plot_waiting_percentage <- function(data) {
    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"),
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu),
             x = "Time", y = "Percentage") +
        theme_minimal() +
        ylim(0, 100)
}

plot_queue_length <- function(data) {
    ggplot(data, aes(x = time, y = system, group = simulation)) +
        geom_line(alpha = 0.1, color = "blue") +
        stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red", size = 1) +
        geom_hline(yintercept = queue_threshold, linetype = "dashed", color = "green") +
        labs(title = "System Size Over Time",
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu, 
                              ", Initial Queue Size =", initial_queue_size),
             x = "Time", y = "System Size") +
        theme_minimal() +
        scale_y_continuous(expand = c(0, 0), limits = c(0, NA))
}

plot_server_capacity <- function(data) {
    ggplot(data, aes(x = time, y = capacity, group = simulation)) +
        geom_line(alpha = 0.1, color = "purple") +
        stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red", size = 1) +
        labs(title = "Server Capacity Over Time",
             subtitle = paste("Queue Threshold =", queue_threshold),
             x = "Time", y = "Number of Servers") +
        theme_minimal() +
        scale_y_continuous(breaks = c(3, 4))
}

plot_cumulative_orders <- function(data) {
    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 = "Cumulative Orders Processed",
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu, 
                              ", Initial Queue Size =", initial_queue_size),
             x = "Time", y = "Cumulative Orders") +
        theme_minimal()
}

waiting_percentage_plot <- plot_waiting_percentage(results$waiting_percentage_data)
queue_length_plot <- plot_queue_length(results$queue_length_data)
server_capacity_plot <- plot_server_capacity(results$server_capacity_data)
cumulative_orders_plot <- plot_cumulative_orders(cumulative_orders)

grid.arrange(waiting_percentage_plot, queue_length_plot,
             server_capacity_plot, cumulative_orders_plot,
             ncol = 2)

The code ran and I got this plot:

Is this how we correctly add alternating servers in simmer?

Thanks!

[Enchufa2]

It's a way of doing it, but I don't think it's the best way. I would rather run a set_capacity activity right before and after the seize in the main one. I.e., something like this:

update_capacity <- function(.env, resource, limit) {
  function() {
    if (get_queue_count(.env, resource) < limit)
      return(3)
    return(4)
  }
}

customer <- trajectory() %>%
  set_capacity("server", update_capacity(env, "server", 299)) %>%
  seize("server", 1) %>%
  set_capacity("server", update_capacity(env, "server", 300)) %>%
  timeout(function() {rexp(1, mu)}) %>%
  release("server", 1)

Note also that:

• You don't need to record the arrival time for each arrival in an attribute. The simulator collects that for you in the get_mon_arrivals dataframe.
• You don't need to record every capacity change. That is also available in the get_mon_resources dataframe.

[swaheera]

Thank you Dr.Ucar! I tried to use the update_capacity()function, but for some reason the code is taking a lot longer to run?

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

# 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

# Function to update server capacity
update_capacity <- function(.env, resource, limit) {
  function() {
    if (get_queue_count(.env, resource) < limit)
      return(3)
    return(4)
  }
}

# Function to run a single simulation
run_simulation <- function(seed) {
    cat("Running simulation", seed, "of", num_simulations, "\n")
    
    set.seed(seed)
    env <- simmer("Dynamic_MMK_Queue")
    
    customer <- trajectory() %>%
        set_capacity("server", update_capacity(env, "server", 299)) %>%
        seize("server", 1) %>%
        set_capacity("server", update_capacity(env, "server", 300)) %>%
        timeout(function() {rexp(1, mu)}) %>%
        release("server", 1)
    
    env %>%
        add_resource("server", capacity = 3) %>%
        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)
    
    list(resources = get_mon_resources(env), arrivals = get_mon_arrivals(env))
}

# Running simulations
simulation_results <- lapply(1:num_simulations, run_simulation)

# Function to process simulation results
process_results <- function(simulation_results) {
    queue_length_data <- do.call(rbind, lapply(1:num_simulations, function(i) {
        resources <- simulation_results[[i]]$resources
        resources %>%
            mutate(simulation = i)
    }))
    
    waiting_time_data <- do.call(rbind, lapply(1:num_simulations, function(i) {
        arrivals <- simulation_results[[i]]$arrivals
        arrivals %>%
            mutate(
                simulation = i,
                waiting_time = end_time - start_time - activity_time
            )
    }))
    
    waiting_percentage_data <- waiting_time_data %>%
        filter(end_time > 0) %>%  # Exclude customers who haven't finished service
        group_by(simulation, time = floor(end_time)) %>%
        summarise(
            total = n(),
            waiting_longer = sum(waiting_time > k_minutes),
            percentage = waiting_longer / total * 100,
            .groups = 'drop'
        )
    
    cumulative_orders_data <- waiting_time_data %>%
        group_by(simulation, time = floor(end_time)) %>%
        summarise(cumulative_orders = n(), .groups = 'drop')
    
    list(queue_length_data = queue_length_data, 
         waiting_percentage_data = waiting_percentage_data,
         cumulative_orders_data = cumulative_orders_data)
}

# Process the results
results <- process_results(simulation_results)

# Plotting functions
plot_waiting_percentage <- function(data) {
    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"),
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu),
             x = "Time", y = "Percentage") +
        theme_minimal() +
        ylim(0, 100)
}

plot_queue_length <- function(data) {
    ggplot(data, aes(x = time, y = system, group = simulation)) +
        geom_line(alpha = 0.1, color = "blue") +
        stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red", size = 1) +
        labs(title = "System Size Over Time",
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu, 
                              ", Initial Queue Size =", initial_queue_size),
             x = "Time", y = "System Size") +
        theme_minimal() +
        scale_y_continuous(expand = c(0, 0), limits = c(0, NA))
}

plot_server_capacity <- function(data) {
    ggplot(data, aes(x = time, y = server, group = simulation)) +
        geom_line(alpha = 0.1, color = "purple") +
        stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red", size = 1) +
        labs(title = "Server Capacity Over Time",
             x = "Time", y = "Number of Servers") +
        theme_minimal() +
        scale_y_continuous(breaks = c(3, 4))
}

plot_cumulative_orders <- function(data) {
    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 = "Cumulative Orders Processed",
             subtitle = paste("Arrival Rate =", lambda, ", Service Rate =", mu, 
                              ", Initial Queue Size =", initial_queue_size),
             x = "Time", y = "Cumulative Orders") +
        theme_minimal()
}

# Generate plots
waiting_percentage_plot <- plot_waiting_percentage(results$waiting_percentage_data)
queue_length_plot <- plot_queue_length(results$queue_length_data)
server_capacity_plot <- plot_server_capacity(results$queue_length_data)
cumulative_orders_plot <- plot_cumulative_orders(results$cumulative_orders_data)

# Display plots
grid.arrange(waiting_percentage_plot, queue_length_plot,
             server_capacity_plot, cumulative_orders_plot,
             ncol = 2)

[Enchufa2]

Are you sure? I takes half the time for me.