Simulating Expired Orders in a Factory

[swaheera]

I have this Operations Research problem :

• Suppose there is a factory that on day 1 has 100 boxes
• New boxes arrive according to some random rate (e.g. Poisson Distribution)
• There are 5 employees at the factory - each employee finishes boxes at the same random rate (e.g. Exponential Distribution)
• Boxes are processed according to the FIFO (First In, First Out)
• As soon as an employee finishes a box, he starts the next box right away
• In the event where there no boxes remaining and 2 or more employees are idle, the employee who has been idle the longest starts the next box

I am trying to answer the following question: After 1000 days - what will be the average queue length and the average number of boxes that have been waiting for more than 5 days (i.e. expired) ? I am not sure if there is a mathematical solution to this problem so I tried to simulate it.

I tried to set the code up like this:

Step 1 : Define parameters

library(dplyr)
library(ggplot2)
library(RColorBrewer)

set.seed(123)

initial_boxes <- 100  
num_workers <- 5  
# rate for exponential distribution (service rate)
rate <- 1  
arrival_rate <- 5  

Step 2: Define queue function

simulate_process <- function(iterations) {
    time <- 0
    queue <- data.frame(box = paste0("box_", 1:initial_boxes), time_entered = time)
    worker_histories <- vector("list", num_workers)
    processed_boxes <- data.frame(time_started = numeric(), time_finished = numeric(), box = character(), worker = integer(), time_taken = numeric())
    queue_length <- data.frame(time = numeric(), queue_length = integer())
    waiting_longer_than_5 <- data.frame(time = numeric(), percent_waiting_longer_than_5 = numeric())
    cumulative_boxes <- data.frame(time = numeric(), cumulative_boxes = integer())
    cumulative_processed_count <- 0
    
    for (i in 1:iterations) {
        num_new_boxes <- rpois(1, arrival_rate)
        new_boxes <- paste0("box_", (nrow(queue) + nrow(processed_boxes) + 1):(nrow(queue) + nrow(processed_boxes) + num_new_boxes))
        queue <- rbind(queue, data.frame(box = new_boxes, time_entered = time))
        
        queue_length <- rbind(queue_length, data.frame(time = time, queue_length = nrow(queue)))
        
        for (worker in 1:num_workers) {
            if (nrow(queue) == 0) break
            box <- queue$box[1]
            time_entered <- queue$time_entered[1]
            queue <- queue[-1,]
            time_started <- ifelse(length(worker_histories[[worker]]) == 0, time, max(sapply(worker_histories[[worker]], function(x) x$time_finished)))
            time_taken <- rexp(1, rate)
            time_finished <- time_started + time_taken
            worker_histories[[worker]] <- c(worker_histories[[worker]], list(list(time_started = time_started, time_finished = time_finished, box = box, time_taken = time_taken)))
            processed_boxes <- rbind(processed_boxes, data.frame(time_started = time_started, time_finished = time_finished, box = box, worker = worker, time_taken = time_taken))
            
            queue_length <- rbind(queue_length, data.frame(time = time, queue_length = nrow(queue)))
        }
        
        if (nrow(queue) > 0) {
            waiting_times <- time - queue$time_entered
            percent_waiting_longer_than_5 <- sum(waiting_times > 5) / nrow(queue) * 100
        } else {
            percent_waiting_longer_than_5 <- 0  
        }
        waiting_longer_than_5 <- rbind(waiting_longer_than_5, data.frame(time = time, percent_waiting_longer_than_5 = percent_waiting_longer_than_5))
        
        cumulative_processed_count <- cumulative_processed_count + sum(sapply(worker_histories, function(w) {
            if (length(w) > 0) {
                return(sum(sapply(w, function(x) x$time_finished <= time)))
            } else {
                return(0)
            }
        }))
        
        cumulative_boxes <- rbind(cumulative_boxes, data.frame(time = time, cumulative_boxes = cumulative_processed_count))
        
        time <- time + 1
    }
    
    return(list(queue_length = queue_length, waiting_longer_than_5 = waiting_longer_than_5, cumulative_boxes = cumulative_boxes))
}

Step 3: Simulate Queue 100 times and calculate summary stats

num_simulations <- 100

all_simulations <- vector("list", num_simulations)
for (sim in 1:num_simulations) {
    all_simulations[[sim]] <- simulate_process(50)
}

queue_length_data <- bind_rows(lapply(all_simulations, `[[`, "queue_length"), .id = "simulation_id")
waiting_longer_than_5_data <- bind_rows(lapply(all_simulations, `[[`, "waiting_longer_than_5"), .id = "simulation_id")
cumulative_boxes_data <- bind_rows(lapply(all_simulations, `[[`, "cumulative_boxes"), .id = "simulation_id")

average_queue_length <- queue_length_data %>%
    group_by(time) %>%
    summarise(avg_queue_length = mean(queue_length))

average_waiting_longer_than_5 <- waiting_longer_than_5_data %>%
    group_by(time) %>%
    summarise(avg_percent_waiting_longer_than_5 = mean(percent_waiting_longer_than_5))

average_cumulative_boxes <- cumulative_boxes_data %>%
    group_by(time) %>%
    summarise(avg_cumulative_boxes = mean(cumulative_boxes))

Step 4: Visualize results

ggplot() +
    geom_line(data = queue_length_data, aes(x = time, y = queue_length, group = simulation_id), 
              color = "blue", alpha = 0.2) +  # Faded individual trajectories
    geom_line(data = average_queue_length, aes(x = time, y = avg_queue_length), 
              color = "blue", size = 1.5, alpha = 0.9) +  
    theme_minimal() +
    labs(title = "Queue Length Over Time Across Simulations", 
         x = "Time", y = "Queue Length") +
    theme(
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12)
    )

ggplot() +
    geom_line(data = waiting_longer_than_5_data, aes(x = time, y = percent_waiting_longer_than_5, group = simulation_id), 
              color = "purple", alpha = 0.2) +  
    geom_line(data = average_waiting_longer_than_5, aes(x = time, y = avg_percent_waiting_longer_than_5), 
              color = "purple", size = 1.5, alpha = 0.9) +  
    theme_minimal() +
    labs(title = "Percentage of Boxes Waiting Longer Than 5 Time Units Over Time", 
         x = "Time", y = "Percent Waiting Longer Than 5 Units") +
    theme(
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12)
    )

ggplot() +
    geom_line(data = cumulative_boxes_data, aes(x = time, y = cumulative_boxes, group = simulation_id), 
              color = "green", alpha = 0.2) + 
    geom_line(data = average_cumulative_boxes, aes(x = time, y = avg_cumulative_boxes), 
              color = "green", size = 1.5, alpha = 0.9) + 
    theme_minimal() +
    labs(title = "Cumulative Boxes Processed Over Time Across Simulations", 
         x = "Time", y = "Cumulative Boxes Processed") +
    theme(
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12)
    )

Can someone please help me understand if I have done this correctly? I can provide more clarifications explaining my code is needed

Thanks!