Simulating a Queue in R

General
1
  • Suppose there is a factory with 5 employees
  • Each day, boxes arrive at the factory according to some random date
  • Employees work on the boxes based on a First In First Out (FIFO) rule : Boxes that have been waiting for longer periods of time are processed before boxes waiting for shorter periods of time
  • The employees can finish a box based on some other random rate. Once they finish a box, they immediately move to the next box (employees are never idle and never take breaks)
  • In the situation that there more than two employees are sitting idle since there are no new boxes, the employee that was sitting idle the longest will be assigned in the newest box

Imagine this factory is open for 10 hours each day. I want to simulate the detailed evolution of this factory over 10 days

Here is my logic.

Step 1: Define initial parameters

library(data.table)

lambda_arrival <- 5  # Arrival rate (Poisson)
lambda_service <- 3  # Service rate (Exponential)
hours_per_day <- 10
days <- 10
workers <- 5

results <- vector("list", days * hours_per_day)
completed_boxes <- vector("list", days * hours_per_day)

box_id_counter <- 1

Step 2: Define a Queue Function - this took me the longest time to do. It was really confusing to synchronize the dynamics of the factory (i.e. ensuring that incomplete boxes from the previous hours are worked on chronologically):

simulate_hour <- function(hour, day, queue, workers_status) {
    # Generate arrivals
    arrivals <- rpois(1, lambda_arrival)
    arrival_times <- rep(hour, arrivals)
    arrival_days <- rep(day, arrivals)
    arrival_ids <- seq(box_id_counter, box_id_counter + arrivals - 1)
    box_id_counter <<- box_id_counter + arrivals
    
    # Update queue with new arrivals
    if (arrivals > 0) {
        new_boxes <- data.table(box_id = arrival_ids, arrival_hour = arrival_times, arrival_day = arrival_days, arrival_id = arrival_ids, waiting_time = 0)
        queue <- rbind(queue, new_boxes)
    }
    
    # Process boxes
    completed <- data.table(worker = character(0), box_id = numeric(0), completion_hour = numeric(0), completion_day = numeric(0), arrival_hour = numeric(0), arrival_day = numeric(0))
    for (i in 1:workers) {
        if (nrow(queue) > 0) {
            service_time <- rexp(1, lambda_service)
            queue[1, waiting_time := waiting_time + service_time]
            completed <- rbind(completed, data.table(worker = paste0("worker_", i), box_id = queue[1, box_id], completion_hour = hour, completion_day = day, arrival_hour = queue[1, arrival_hour], arrival_day = queue[1, arrival_day]))
            queue <- queue[-1]
        }
    }
    
    # Update waiting times for remaining boxes
    queue[, waiting_time := waiting_time + 1]
    
    # Store results
    results[[hour + (day - 1) * hours_per_day]] <<- data.table(current_hour = hour, current_day = day, queue)
    completed_boxes[[hour + (day - 1) * hours_per_day]] <<- completed
    
    return(queue)
}

# Simulate 10 days
queue <- data.table(box_id = numeric(0), arrival_hour = numeric(0), arrival_day = numeric(0), arrival_id = numeric(0), waiting_time = numeric(0))
for (day in 1:days) {
    for (hour in 1:hours_per_day) {
        queue <- simulate_hour(hour, day, queue, workers)
    }
}

Step 3: Run the simulation

queue <- data.table(box_id = numeric(0), arrival_hour = numeric(0), arrival_day = numeric(0), arrival_id = numeric(0), waiting_time = numeric(0))
for (day in 1:days) {
    for (hour in 1:hours_per_day) {
        queue <- simulate_hour(hour, day, queue, workers)
    }
}

print(results[[1]])
print(completed_boxes[[1]])

Problem: In this simulation, a worker only finishes one box an hour. How can I make it so that when a worker is finished a box, they start working on the next one?

Thanks!

Appendix: Here is a function to summarize the factory

summarize_factory_status <- function(results, completed_boxes, hours_per_day, days) {
    summary <- data.table(
        day = rep(1:days, each = hours_per_day),
        hour = rep(1:hours_per_day, days),
        opening_inventory = 0,
        arrivals = 0,
        completed_boxes = 0,
        closing_inventory = 0
    )
    
    for (i in 1:(days * hours_per_day)) {
        current_result <- results[[i]]
        current_completed <- completed_boxes[[i]]
        
        if (i == 1) {
            summary[i, opening_inventory := 0]
        } else {
            summary[i, opening_inventory := summary[i-1, closing_inventory]]
        }
        
        summary[i, arrivals := nrow(current_result[arrival_hour == current_result[1, current_hour] & 
                                                       arrival_day == current_result[1, current_day]])]
        summary[i, completed_boxes := nrow(current_completed)]
        summary[i, closing_inventory := nrow(current_result)]
    }
    
    return(summary)
}

factory_summary <- summarize_factory_status(results, completed_boxes, hours_per_day, days)
data.frame(factory_summary)

A possible fix?

library(data.table)

# Step 1: Define initial parameters
lambda_arrival <- 5  # Arrival rate (Poisson)
lambda_service <- 3  # Service rate (Exponential)
hours_per_day <- 10
days <- 10
workers <- 5

results <- vector("list", days * hours_per_day)
completed_boxes <- vector("list", days * hours_per_day)

box_id_counter <- 1

# Step 2: Define a Queue Function
simulate_hour <- function(hour, day, queue, workers_status) {
    # Generate arrivals
    arrivals <- rpois(1, lambda_arrival)
    arrival_times <- rep(hour, arrivals)
    arrival_days <- rep(day, arrivals)
    arrival_ids <- seq(box_id_counter, box_id_counter + arrivals - 1)
    box_id_counter <<- box_id_counter + arrivals
    
    # Update queue with new arrivals
    if (arrivals > 0) {
        new_boxes <- data.table(box_id = arrival_ids, arrival_hour = arrival_times, arrival_day = arrival_days, arrival_id = arrival_ids, waiting_time = 0)
        queue <- rbind(queue, new_boxes)
    }
    
    # Process boxes
    completed <- data.table(worker = character(0), box_id = numeric(0), completion_hour = numeric(0), completion_day = numeric(0), arrival_hour = numeric(0), arrival_day = numeric(0))
    time_remaining <- 1  # Start with 1 hour available
    
    while (time_remaining > 0 && nrow(queue) > 0) {
        for (i in 1:workers) {
            if (nrow(queue) > 0) {
                service_time <- rexp(1, lambda_service)
                if (service_time <= time_remaining) {
                    time_remaining <- time_remaining - service_time
                    completed <- rbind(completed, data.table(worker = paste0("worker_", i), box_id = queue[1, box_id], completion_hour = hour, completion_day = day, arrival_hour = queue[1, arrival_hour], arrival_day = queue[1, arrival_day]))
                    queue <- queue[-1]  # Remove the processed box from the queue
                } else {
                    # Update waiting time if the service time extends beyond the current hour
                    queue[1, waiting_time := waiting_time + time_remaining]
                    time_remaining <- 0  # No more time left in the hour
                    break
                }
            }
        }
    }
    
    # Update waiting times for remaining boxes
    queue[, waiting_time := waiting_time + 1]
    
    # Store results
    results[[hour + (day - 1) * hours_per_day]] <<- data.table(current_hour = hour, current_day = day, queue)
    completed_boxes[[hour + (day - 1) * hours_per_day]] <<- completed
    
    return(queue)
}

# Step 3: Simulate 10 days
queue <- data.table(box_id = numeric(0), arrival_hour = numeric(0), arrival_day = numeric(0), arrival_id = numeric(0), waiting_time = numeric(0))
for (day in 1:days) {
    for (hour in 1:hours_per_day) {
        queue <- simulate_hour(hour, day, queue, workers)
    }
}

# Step 4: Summarize the factory status
summarize_factory_status <- function(results, completed_boxes, hours_per_day, days) {
    summary <- data.table(
        day = rep(1:days, each = hours_per_day),
        hour = rep(1:hours_per_day, days),
        opening_inventory = 0,
        arrivals = 0,
        completed_boxes = 0,
        closing_inventory = 0
    )
    
    for (i in 1:(days * hours_per_day)) {
        current_result <- results[[i]]
        current_completed <- completed_boxes[[i]]
        
        if (i == 1) {
            summary[i, opening_inventory := 0]
        } else {
            summary[i, opening_inventory := summary[i-1, closing_inventory]]
        }
        
        summary[i, arrivals := nrow(current_result[arrival_hour == current_result[1, current_hour] & 
                                                       arrival_day == current_result[1, current_day]])]
        summary[i, completed_boxes := nrow(current_completed)]
        summary[i, closing_inventory := nrow(current_result)]
    }
    
    return(summary)
}

# Generate the summary
factory_summary <- summarize_factory_status(results, completed_boxes, hours_per_day, days)

# Print the summary data frame
print(data.frame(factory_summary))
2

I would guess you should do your loop at the required timing resolution. Presumably you cant succeeded by looping over hours, as you want events to happen in partial hours ? what is the appropriate resolution ? minutes ? seconds ? I would guess minutes. or even groups of minutes, like 5min chunks .

something to think about.

1 Like
3

@nirgrahamuk : thank you for your reply! I made change (possible fix) - what do you think about this?

4

@nirgrahamuk

Can you please take a look at this question?
Simulating Expired Orders in a Factory - General - Posit Community