Biostat 213

Midterm

Author

Jiachen Ai; UID: 206182615

Load necessary libraries

library(tidyverse)
library(dplyr)
library(ggplot2)
library(gridExtra)
library(grid)
library(reshape2)
library(gtsummary)
library(kableExtra)
set.seed(2023)

Problem Statement

Conduct a simulation on your Urgent Care Center business. Here are some simplifying assumptions. Feel free to change if you wish:

Assume that your business opens at 8am and closes at 5pm. Customers arrive at a rate of 4/hr.
Between the hours of noon (12) and 1 that rate increases to 6/hr.
Suppose you have 2 tellers. There are, in general, 3 types of queuing systems:
- 1. 1 line with tandem servers.
- 1. oneline service (In-N-out example)
- 1. where each teller has their own line (grocery store).

Questions:

What is the expected waiting time for each type of system. Simulate 1 week (do 7 times and take the average).
Which system do you think is the best?
Does it matter where you place the “slow” teller?
Does you “optimal” depend on your choices?

Solution

Generate customer arrival times based on varying rates

First, we will simulate the customer arrivals times for the 1 line with two tellers.

# Parameters
lambda1 <- 4 / 60       # Arrival rate (4 per hour) before 12 pm and after 1 pm
lambda2 <- 6 / 60       # Arrival rate (6 per hour) between 12 pm and 1 pm
T <- 8 * 60             # Total simulation time per day in minutes (8 hours)
midday_start <- 4 * 60  # Start time of increased rate (12 pm in minutes from 8 am)
midday_end <- 5 * 60    # End time of increased rate (1 pm in minutes from 8 am)

# Function to generate customer arrivals with two different rates
generate_arrivals <- function(lambda1, lambda2, T, midday_start, midday_end) {
  current_time <- 0
  arrival_times <- c()
  
  while (current_time < T) {
    # Choose rate based on time of day
    if (current_time >= midday_start & current_time < midday_end) {
      interarrival_time <- rexp(1, rate = lambda2)
    } else {
      interarrival_time <- rexp(1, rate = lambda1)
    }
    
    # Update current time
    current_time <- current_time + interarrival_time
    if (current_time < T) {
      arrival_times <- c(arrival_times, current_time)
    }
  }
  return(arrival_times)
}

Plot customer arrival times on a single line.

# Generate 7 sets of arrival times and store them in a data frame
arrival_data <- do.call(rbind, lapply(1:7, function(day) {
  arrival_times <- generate_arrivals(lambda1, lambda2, T, midday_start, midday_end)
  data.frame(
    Day = paste("Day", day),
    time = arrival_times,
    customer = 1:length(arrival_times),
    label = ifelse(arrival_times >= midday_start & arrival_times < midday_end, "12-1 pm", "other times")
  )
}))


# Split arrival_data by day to get arrival times for each day
arrival_times_by_day <- split(arrival_data$time, arrival_data$Day)

# Plot using ggplot2 with facet_wrap for shared axes
ggplot(arrival_data, aes(x = time, y = customer, color = label)) +
  geom_point() +
  facet_wrap(~ Day, ncol = 3) +
  theme_bw() +
  labs(title = "Customer Arrival Times Over 7 Days", 
       x = "Time (minutes)", 
       y = "Customer", 
       color = "Time Period") +
  theme(legend.position = "top", 
        plot.title = element_text(size = 10),
        axis.text.x = element_text(angle = 45, hjust = 1))

Case 1: All servers have same rate and do type 1, 2, 3.

Assumptions:

Customers arrive following a Poisson process with two arrival rates (4/hr generally, increasing to 6/hr between 12-1 pm).
Each teller (representing a stage) has an exponential service time with rate \(\mu\) hours.
Customers first see Teller 1, then move to Teller 2 once Teller 1’s service is complete.

Type 1: One line with two servers (Tandem)

In this system, customers first go to Teller 1, then move to Teller 2 once Teller 1’s service is complete.

Code

# Set up the data for the two tellers
teller_data <- data.frame(
  x = c(1, 2),
  y = c(0, 0),
  label = c("Teller 1", "Teller 2")
)

# Create the plot
ggplot() +
  # Add the line connecting the tellers
  geom_segment(aes(x = 0.5, y = 0, xend = 2.5, yend = 0), size = 1, color = "black") +
  # Add points for each teller location
  geom_point(data = teller_data, aes(x = x, y = y), shape = 1, size = 8, color = "black") +
  # Add labels for each teller
  geom_text(data = teller_data, aes(x = x, y = y + 0.1, label = label), vjust = -1, hjust = 0.5, size = 5) +
  # Add an arrow to the end of the line
  geom_segment(aes(x = 2.5, y = 0, xend = 3, yend = 0),
               arrow = arrow(length = unit(0.3, "cm")), size = 1, color = "black") +
  # Set plot limits and remove axis lines
  coord_cartesian(xlim = c(0, 3), ylim = c(-0.5, 0.5)) +
  theme_void() %>% suppressMessages()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Here, I will the service rate is \(\mu = 5 \text{ customers/hr}\).

mu <- 5 / 60            # Service rate for each teller in minutes

# Function to simulate 1 line with tandem servers for one day
simulate_tandem_servers_day <- function(arrival_times, mu) {

  wait_times_stage1 <- c()
  wait_times_stage2 <- c()
  service_times_stage1 <- c()
  service_times_stage2 <- c()
  
  # Initialize tellers' availability times
  teller1_free_at <- 0
  teller2_free_at <- 0
  
  # Process each arrival in the line with tandem servers
  for (arrival_time in arrival_times) {
    
    # For each customer, first go to Teller 1, then Teller 2
    # Calculate the wait time and service time for each teller
    
    # Stage 1: Teller 1
    wait_time1 <- max(0, teller1_free_at - arrival_time)
    service_time1 <- rexp(1, rate = mu)
    teller1_free_at <- arrival_time + wait_time1 + service_time1
    wait_times_stage1 <- c(wait_times_stage1, wait_time1)
    service_times_stage1 <- c(service_times_stage1, service_time1)
    
    # Stage 2: Teller 2
    arrival_time_stage2 <- teller1_free_at
    wait_time2 <- max(0, teller2_free_at - arrival_time_stage2)
    service_time2 <- rexp(1, rate = mu)
    teller2_free_at <- arrival_time_stage2 + wait_time2 + service_time2
    wait_times_stage2 <- c(wait_times_stage2, wait_time2)
    service_times_stage2 <- c(service_times_stage2, service_time2)
  }
  
  # Calculate total wait times
  total_wait_times <- wait_times_stage1 + wait_times_stage2
  avg_wait_time <- mean(total_wait_times)
  
  return(list(avg_wait_time = avg_wait_time, 
              arrival_times = arrival_times, 
              wait_times_stage1 = wait_times_stage1,
              wait_times_stage2 = wait_times_stage2, 
              total_wait_times = total_wait_times))
}

Now, we will run the simulation for 7 days and calculate the average wait time over the week.

# Apply the simulate_tandem_servers_day function to each day's arrival times
tandem_weekly_results <- lapply(arrival_times_by_day, function(arrival_times) {
  simulate_tandem_servers_day(arrival_times, mu)
})

# Summary of results
tandem_avg_weekly_wait_time <- mean(tandem_weekly_results %>% 
                               sapply(`[[`, "avg_wait_time"))
cat("Average wait time over 7 days for tandem servers:", 
    tandem_avg_weekly_wait_time / 60, "\n")

Average wait time over 7 days for tandem servers: 1.405915

Finally, we will plot the wait times for each teller and the total wait time for each customer over the week.

Code

# Combine data for all days into a single data frame
all_days_data <- do.call(rbind, lapply(1:7, function(i) {
  data.frame(
    Day = paste("Day", i),
    ArrivalTime = tandem_weekly_results[[i]]$arrival_times,
    WaitTimeStage1 = tandem_weekly_results[[i]]$wait_times_stage1,
    WaitTimeStage2 = tandem_weekly_results[[i]]$wait_times_stage2,
    TotalWaitTime = tandem_weekly_results[[i]]$total_wait_times
  )
}))

# Melt data into a long format for plotting
all_days_data_melted <- melt(all_days_data, id.vars = c("Day", "ArrivalTime"), 
                             variable.name = "Teller", value.name = "WaitTime")

# Custom color palette in the style of Nature or Science journals
custom_colors <- c("WaitTimeStage1" = "#1f78b4",  # Blue
                   "WaitTimeStage2" = "#33a02c",  # Green
                   "TotalWaitTime" = "#e31a1c")   # Red

# Plot using ggplot2 with facet_wrap for shared axes
ggplot(all_days_data_melted, aes(x = ArrivalTime, y = WaitTime, color = Teller)) +
  geom_step(size = 0.8) +
  facet_wrap(~ Day, ncol = 3) +
  labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
       color = "Wait Time Type") +
  scale_color_manual(values = custom_colors, 
                     labels = c("Teller 1", "Teller 2", "Total")) +
  theme_bw() +
  theme(legend.position = "top", 
        plot.title = element_text(size = 10),
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank()) +
  ggtitle("Customer Wait Times Over 7 Days for tandem Servers") %>% suppressMessages()

Since the simulation results can vary due to randomness, we can run the simulation multiple times (1000 times) to get a more accurate estimate of the average wait time over 7 days.

# Perform 1000 simulations over 7 days, for each day's arrival times
tandem_sim_results <- replicate(1000, {
  # Apply simulate_tandem_servers_day for each day's arrival times
  lapply(arrival_times_by_day, function(arrival_times) {
    simulate_tandem_servers_day(arrival_times, mu)
  })
}, simplify = FALSE)

# Extract and calculate average wait times for each simulation over 7 days
tandem_avg_total_wait_times <- sapply(tandem_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$total_wait_times)))
})
tandem_avg_wait_times_stage1 <- sapply(tandem_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$wait_times_stage1)))
})
tandem_avg_wait_times_stage2 <- sapply(tandem_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$wait_times_stage2)))
})

# Create a data frame for the summary
tandem_results_df <- data.frame(
  Avg_Total_Wait_Time = tandem_avg_total_wait_times / 60,      # Convert to hours
  Avg_Wait_Time_Server1 = tandem_avg_wait_times_stage1 / 60,   # Convert to hours
  Avg_Wait_Time_Server2 = tandem_avg_wait_times_stage2 / 60    # Convert to hours
)

# Create a summary table with gtsummary
tandem_summary_table <- tandem_results_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    label = list(
      Avg_Total_Wait_Time ~ "Average Total Wait Time (hours)",
      Avg_Wait_Time_Server1 ~ "Average Wait Time for Server 1 (hours)",
      Avg_Wait_Time_Server2 ~ "Average Wait Time for Server 2 (hours)"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places
  ) %>%
  modify_header(label = "**Statistic**") %>%
  modify_caption("**Tandem Service Wait Times Over 7 Days**") 

# Print the summary table
tandem_summary_table

**Tandem Service Wait Times Over 7 Days**
Statistic	N = 1,000¹
Average Total Wait Time (hours)	1.18 ± 0.17
Average Wait Time for Server 1 (hours)	0.71 ± 0.15
Average Wait Time for Server 2 (hours)	0.47 ± 0.13
¹ Mean ± SD

Type 2: One Line, Multiple Servers (Customers Choose)

In this case, there is also only one queue. Customer will go to either server 1 or server 2. And customers will get their food while ordering. For example, in an In-N-Out restaurant, and customers go to the first server by default to place their order if both servers are available. If one server is busy, the customer will go to the other server. If both servers are busy, the customer will wait in the queue and go to the first server when available.

Code

# Data for tellers and customers
teller_data <- data.frame(
  x = c(1, 3),
  y = c(4, 4),
  label = c("Teller 1", "Teller 2")
)

customer_data <- data.frame(
  x = rep(2, 4),
  y = c(2.5, 2, 1.5, 1),
  label = rep("Customer", 4)
)

# Create the plot
ggplot() +
  # Add circles for tellers
  geom_point(data = teller_data, aes(x = x, y = y), shape = 1, size = 8, color = "black") +
  # Add labels for tellers
  geom_text(data = teller_data, aes(x = x, y = y + 0.2, label = label), vjust = -1, hjust = 0.5, size = 5) +
  
  # Add squares for customers
  geom_point(data = customer_data, aes(x = x, y = y), shape = 15, size = 5, color = "black") +
  
  # Add arrows from tellers to the first customer position
  geom_segment(aes(x = 2, y = 2.5, xend = 1.2, yend = 3.8), linetype=2,
               arrow = arrow(length = unit(0.3, "cm")), color = "black") +
  geom_segment(aes(x = 2, y = 2.5, xend = 2.8, yend = 3.8), linetype=2,
               arrow = arrow(length = unit(0.3, "cm")), color = "black") +
  
  # Add "or" text between the arrows
  annotate("text", x = 2, y = 3, label = "or", size = 5, fontface = "italic") +
  
  # Set plot limits and remove axis lines
  coord_cartesian(xlim = c(0.5, 3.5), ylim = c(0.5, 5)) +
  theme_void()

Since, in this case, the customer will get their food while ordering, the service speed will be slower than the previous case.

mu_slower <- 2.5 / 60            # Service rate for each teller in minutes

# Function to simulate oneline service with 1 line and 2 tellers for one day
simulate_oneline_service_day <- function(arrival_times, mu) {
  
  wait_times <- c()
  service_times <- c()
  
  # Initialize tellers' availability times
  teller1_free_at <- 0
  teller2_free_at <- 0
  
  # Process each customer arrival
  for (arrival_time in arrival_times) {
    
    # Determine which teller will serve the customer
    
    # Both tellers are available, go to Teller 1 by default
    if (arrival_time >= teller1_free_at && arrival_time >= teller2_free_at) {
      wait_time <- 0
      service_time <- rexp(1, rate = mu)
      teller1_free_at <- arrival_time + service_time
      
      # Teller 1 is available
    } else if (arrival_time >= teller1_free_at) {
      wait_time <- 0
      service_time <- rexp(1, rate = mu)
      teller1_free_at <- arrival_time + service_time
      
      # Teller 2 is available
    } else if (arrival_time >= teller2_free_at) {
      wait_time <- 0
      service_time <- rexp(1, rate = mu)
      teller2_free_at <- arrival_time + service_time
      
      # Both tellers are busy, customer waits in queue and for the next available teller
    } else {
      next_free_time <- min(teller1_free_at, teller2_free_at)
      wait_time <- next_free_time - arrival_time
      
      # Assign to the teller who becomes available first
      if (teller1_free_at <= teller2_free_at) {
        service_time <- rexp(1, rate = mu)
        teller1_free_at <- next_free_time + service_time
      } else {
        service_time <- rexp(1, rate = mu)
        teller2_free_at <- next_free_time + service_time
      }
    }
    
    # Record wait and service times
    wait_times <- c(wait_times, wait_time)
    service_times <- c(service_times, service_time)
  }
  
  # Calculate total wait times
  total_wait_times <- wait_times
  avg_wait_time <- mean(total_wait_times)
  
  return(list(avg_wait_time = avg_wait_time, 
              arrival_times = arrival_times, 
              wait_times = wait_times, 
              service_times = service_times,
              total_wait_times = total_wait_times))
}

Now, we will run the simulation for 7 days and calculate the average wait time over the week. Here, in order to keep the consistency, we will use the same arrival times for each day.

# Apply simulate_oneline_service_day function to each day's arrival times
oneline_weekly_results <- lapply(arrival_times_by_day, function(arrival_times) {
  simulate_oneline_service_day(arrival_times, mu_slower)
})

# Calculate the average wait time over 7 days
oneline_avg_weekly_wait_time <- mean(sapply(oneline_weekly_results, `[[`, "avg_wait_time"))

# Print the average weekly wait time in hours
cat("Average wait time over 7 days for oneline service (hours):", oneline_avg_weekly_wait_time / 60, "\n")

Average wait time over 7 days for oneline service (hours): 0.8161425

I think the the serving speed \(\mu\) is important for the wait time.

Code

# Combine data for all days into a single data frame
all_days_data <- do.call(rbind, lapply(1:7, function(i) {
  data.frame(
    Day = paste("Day", i),
    ArrivalTime = oneline_weekly_results[[i]]$arrival_times,
    WaitTime = oneline_weekly_results[[i]]$total_wait_times
  )
}))

# Plot for all_days_data using ggplot2 with facet_wrap
ggplot(all_days_data, aes(x = ArrivalTime, y = WaitTime)) +
  geom_step(color = "#e31a1c", size = 0.8) +
  facet_wrap(~ Day, ncol = 3) +
  labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
       title = "Wait Times for oneline Service Over 7 Days") +
  theme_bw(base_size = 12) +
  theme(
    legend.position = "none",  # No legend needed for single color
    panel.grid.major = element_blank(),  # Remove major grid lines
    panel.grid.minor = element_blank(),  # Remove minor grid lines
    strip.text = element_text(size = 10),  # Bold facet labels
    plot.title = element_text(size = 14, hjust = 0.5),  # Centered and bold title
    axis.text = element_text(size = 10),  # Adjust axis text size
    axis.title = element_text(size = 12),  # Bold axis titles
    axis.text.x = element_text(angle = 45, hjust = 1)  # Rotate x-axis labels
  ) %>%
  suppressMessages()

Since the simulation results can vary due to randomness, we can run the simulation multiple times (1000 times) to get a more accurate estimate of the average wait time over 7 days.

# Perform 1000 simulations over 7 days, using each day's arrival times
oneline_sim_results <- replicate(200, {
  # Apply simulate_oneline_service_day for each day's arrival times
  lapply(arrival_times_by_day, function(arrival_times) {
    simulate_oneline_service_day(arrival_times, mu_slower)
  })
}, simplify = FALSE)

# Extract and calculate average total wait times for each simulation over 7 days
oneline_avg_total_wait_times <- sapply(oneline_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$total_wait_times)))
})

# Convert average total wait times from minutes to hours
oneline_avg_total_wait_times_hours <- oneline_avg_total_wait_times / 60

# Create a data frame for the summary
oneline_results_df <- data.frame(
  Avg_Total_Wait_Time = oneline_avg_total_wait_times_hours  # Average wait time in hours
)

# Create a summary table with gtsummary
oneline_summary_table <- oneline_results_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    label = list(
      Avg_Total_Wait_Time ~ "Average Total Wait Time (hours)"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places
  ) %>%
  modify_header(label = "**Statistic**") %>%
  modify_caption("**Oneline Service Wait Times Over 7 Days**") 

# Print the summary table
oneline_summary_table

**Oneline Service Wait Times Over 7 Days**
Statistic	N = 200¹
Average Total Wait Time (hours)	0.61 ± 0.15
¹ Mean ± SD

Type 3: Each teller has their own line (grocery store)

In this scenario, each teller has their own line, and customers choose the shortest line to join. This is similar to the setup in a grocery store, where each checkout lane has its own queue.

Code

# Data for tellers
teller_data <- data.frame(
  x = c(1, 3),
  y = c(5, 5),
  label = c("Teller 1", "Teller 2")
)

# Data for customers in Teller 1's and Teller 2's queues
customer_data <- data.frame(
  x = rep(c(1, 3), each = 3),
  y = c(4, 3, 2, 4, 3, 2),
  label = "Customer"
)

# Create the plot
ggplot() +
  # Add circles for tellers
  geom_point(data = teller_data, aes(x = x, y = y), shape = 1, size = 8, color = "black") +
  # Add labels for tellers
  geom_text(data = teller_data, aes(x = x, y = y + 0.2, label = label), vjust = -1, hjust = 0.5, size = 5) +
  
  # Add squares for customers in each teller's queue
  geom_point(data = customer_data, aes(x = x, y = y), shape = 15, size = 5, color = "black") +
  
  # Set plot limits and remove axis lines
  coord_cartesian(xlim = c(0.5, 3.5), ylim = c(1, 6)) +
  theme_void()

mu_slower <- 2.5 / 60            # Service rate for each teller in minutes

# Function to simulate service with each teller having their own line for one day
simulate_separate_lines_day <- function(arrival_times, mu) {
  
  # Initialize variables for each teller
  wait_times_teller1 <- c()
  wait_times_teller2 <- c()
  service_times_teller1 <- c()
  service_times_teller2 <- c()
  
  # Initialize tellers' availability times and queue lengths
  teller1_free_at <- 0
  teller2_free_at <- 0
  queue_length_teller1 <- 0
  queue_length_teller2 <- 0
  
  # Initialize an empty data frame to store results
  results <- data.frame(
    ArrivalTime = numeric(0), 
    WaitTime = numeric(0), 
    ServiceTime = numeric(0), 
    Teller = character(0),
    QueueLength_Teller1 = numeric(0),
    QueueLength_Teller2 = numeric(0)
  )
  
  # Process each customer arrival
  for (arrival_time in arrival_times) {
    
    # Record the current queue lengths
    current_queue_length1 <- queue_length_teller1
    current_queue_length2 <- queue_length_teller2
    
    # Determine which teller the customer will go to based on the shorter queue
    if (queue_length_teller1 < queue_length_teller2 || 
        (queue_length_teller1 == queue_length_teller2 && runif(1) < 0.5)) {
      
      # Assign to Teller 1
      if (arrival_time >= teller1_free_at) {
        # Teller 1 is free, no wait time
        wait_time <- 0
        service_time <- rexp(1, rate = mu)
        teller1_free_at <- arrival_time + service_time
      } else {
        # Teller 1 is busy, customer waits
        wait_time <- teller1_free_at - arrival_time
        service_time <- rexp(1, rate = mu)
        teller1_free_at <- teller1_free_at + service_time
        queue_length_teller1 <- queue_length_teller1 + 1
      }
      
      # Record the result for Teller 1
      results <- rbind(results, data.frame(
        ArrivalTime = arrival_time, 
        WaitTime = wait_time, 
        ServiceTime = service_time, 
        Teller = "Teller 1",
        QueueLength_Teller1 = current_queue_length1,
        QueueLength_Teller2 = current_queue_length2
      ))
      
    } else {
      
      # Assign to Teller 2
      if (arrival_time >= teller2_free_at) {
        # Teller 2 is free, no wait time
        wait_time <- 0
        service_time <- rexp(1, rate = mu)
        teller2_free_at <- arrival_time + service_time
      } else {
        # Teller 2 is busy, customer waits
        wait_time <- teller2_free_at - arrival_time
        service_time <- rexp(1, rate = mu)
        teller2_free_at <- teller2_free_at + service_time
        queue_length_teller2 <- queue_length_teller2 + 1
      }
      
      # Record the result for Teller 2
      results <- rbind(results, data.frame(
        ArrivalTime = arrival_time, 
        WaitTime = wait_time, 
        ServiceTime = service_time, 
        Teller = "Teller 2",
        QueueLength_Teller1 = current_queue_length1,
        QueueLength_Teller2 = current_queue_length2
      ))
      
    }
    
    # Update queue lengths after the patient is served (if applicable)
    if (arrival_time >= teller1_free_at) {
      queue_length_teller1 <- max(0, queue_length_teller1 - 1)
    }
    if (arrival_time >= teller2_free_at) {
      queue_length_teller2 <- max(0, queue_length_teller2 - 1)
    }
  }
  
  # Calculate average wait times for each teller
  avg_wait_time_teller1 <- mean(results$WaitTime[results$Teller == "Teller 1"])
  avg_wait_time_teller2 <- mean(results$WaitTime[results$Teller == "Teller 2"])
  
  return(list(
    results = results,
    avg_wait_time_teller1 = avg_wait_time_teller1,
    avg_wait_time_teller2 = avg_wait_time_teller2
  ))
}

Now, we will run the simulation for 7 days and calculate the average wait time over the week. Here, in order to keep the consistency, we will use the same arrival times for each day.

# Apply simulate_separate_lines_day function to each day's arrival times
separate_weekly_results <- lapply(arrival_times_by_day, function(arrival_times) {
  simulate_separate_lines_day(arrival_times, mu_slower)
})

# Calculate the average weekly wait time for each teller
avg_weekly_wait_time_teller1 <- mean(sapply(separate_weekly_results, `[[`, "avg_wait_time_teller1"))
avg_weekly_wait_time_teller2 <- mean(sapply(separate_weekly_results, `[[`, "avg_wait_time_teller2"))

# Convert from minutes to hours for reporting
avg_weekly_wait_time_teller1_hours <- avg_weekly_wait_time_teller1 / 60
avg_weekly_wait_time_teller2_hours <- avg_weekly_wait_time_teller2 / 60

# Print the average weekly wait times for each teller
cat("Average weekly wait time for Teller 1 (hours):", round(avg_weekly_wait_time_teller1_hours, 2), "\n")

Average weekly wait time for Teller 1 (hours): 0.5

cat("Average weekly wait time for Teller 2 (hours):", round(avg_weekly_wait_time_teller2_hours, 2), "\n")

Average weekly wait time for Teller 2 (hours): 0.55

Plotting the distribution of wait times for each teller over the week.

Code

# Combine data for all days into a single data frame for both tellers
all_days_data <- do.call(rbind, lapply(1:7, function(i) {
  day_data <- separate_weekly_results[[i]]$results
  day_data$Day <- paste("Day", i) 
  return(day_data)
}))

# Convert 'Day' column to a factor for ordering in facets
all_days_data$Day <- factor(all_days_data$Day, levels = paste("Day", 1:7))

# Plot using ggplot2 with facet_wrap
ggplot(all_days_data, aes(x = ArrivalTime, y = WaitTime, color = Teller)) +
  geom_step(size = 0.8) +
  facet_wrap(~ Day, ncol = 3) +
  labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
       title = "Wait Times for Separate Lines Over 7 Days") +
  scale_color_manual(values = c("Teller 1" = "#1f78b4", "Teller 2" = "#33a02c")) +  # Colors for tellers
  theme_bw(base_size = 12) +
  theme(
    legend.position = "top", 
    panel.grid.major = element_blank(),  # Remove major grid lines
    panel.grid.minor = element_blank(),  # Remove minor grid lines
    strip.text = element_text(size = 10),  # Bold facet labels
    plot.title = element_text(size = 14, hjust = 0.5),  # Centered and bold title
    axis.text = element_text(size = 10),  # Adjust axis text size
    axis.title = element_text(size = 12),  # Bold axis titles
    axis.text.x = element_text(angle = 45, hjust = 1)  # Rotate x-axis labels
  ) %>%
  suppressMessages()

Finally, since the simulation is stochastic, we can run it multiple times to get a sense of the variability in the results.

# Perform 100 simulations over 7 days, using each day's arrival times
separate_sim_results <- replicate(1000, {
  # Apply simulate_separate_lines_day for each day's arrival times
  lapply(arrival_times_by_day, function(arrival_times) {
    simulate_separate_lines_day(arrival_times, mu_slower)
  })
}, simplify = FALSE)

# Extract and calculate average wait times and queue lengths for each simulation over 7 days
separate_avg_wait_time_teller1 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$WaitTime[day_result$results$Teller == "Teller 1"])))
})

separate_avg_wait_time_teller2 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$WaitTime[day_result$results$Teller == "Teller 2"])))
})

separate_avg_queue_length_teller1 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$QueueLength_Teller1[day_result$results$Teller == "Teller 1"])))
})

separate_avg_queue_length_teller2 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$QueueLength_Teller2[day_result$results$Teller == "Teller 2"])))
})

# Convert average wait times from minutes to hours
separate_avg_wait_time_teller1_hours <- separate_avg_wait_time_teller1 / 60
separate_avg_wait_time_teller2_hours <- separate_avg_wait_time_teller2 / 60

# Create a data frame for the summary
separate_results_df <- data.frame(
  Avg_Wait_Time_Teller1 = separate_avg_wait_time_teller1_hours,        # Average wait time for Teller 1 in hours
  Avg_Wait_Time_Teller2 = separate_avg_wait_time_teller1_hours,        # Average wait time for Teller 2 in hours
  Avg_Queue_Length_Teller1 = separate_avg_queue_length_teller1,        # Average queue length for Teller 1
  Avg_Queue_Length_Teller2 = separate_avg_queue_length_teller2         # Average queue length for Teller 2
)

# Create a summary table with gtsummary
separate_summary_table <- separate_results_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    label = list(
      Avg_Wait_Time_Teller1 ~ "Average Wait Time for Teller 1 (hours)",
      Avg_Wait_Time_Teller2 ~ "Average Wait Time for Teller 2 (hours)",
      Avg_Queue_Length_Teller1 ~ "Average Queue Length for Teller 1",
      Avg_Queue_Length_Teller2 ~ "Average Queue Length for Teller 2"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places
  ) %>%
  modify_header(label = "**Statistic**") %>%
  modify_caption("**Separate Line Wait Times and Queue Lengths Over 7 Days**")

# Print the summary table
separate_summary_table

**Separate Line Wait Times and Queue Lengths Over 7 Days**
Statistic	N = 1,000¹
Average Wait Time for Teller 1 (hours)	0.75 ± 0.20
Average Wait Time for Teller 2 (hours)	0.75 ± 0.20
Average Queue Length for Teller 1	4.38 ± 0.48
Average Queue Length for Teller 2	4.38 ± 0.48
¹ Mean ± SD

Conclusion

Assumptions:

The simulation assumes that customers arrive according to a Poisson process, and service times are exponentially distributed. The simulation also assumes that customers are served on a first-come, first-served basis.

In this simulation, I assume that the service rate for the oneline System system, each teller has a service rate of \(\mu = 5 \text{ customers per hour}\), while for the Separate Lines and One Line systems, each teller has a service rate of \(\mu = 2.5 \text{ customers per hour}\). Because, in the oneline System system, both tellers serve customers simultaneously and they divide the total service rate of \(\mu = 5 \text{ customers per hour}\), while in the Separate Lines and One Line systems, each teller need to finish the service before the next teller can start serving the next customer.

System Type	Service Rate per Teller (μ)	Description
Tandem System	5 customers per hour	Each teller has a service rate of 5 customers per hour, as they serve customers simultaneously.
Separate Lines	2.5 customers per hour	Each teller has a service rate of 2.5 customers per hour, as they must finish serving one customer before starting the next.
One Line	2.5 customers per hour	Similar to Separate Lines, each teller has a service rate of 2.5 customers per hour, with one teller needing to finish before the other starts.

Based on the simulation results, the One Line System has the lowest average total wait time, followed by the Separate Lines system and the Tandem system. The One Line System also has the lowest average queue length for both tellers, indicating efficient utilization of resources.

Code

# Combine the data frames into a single data frame for a cohesive summary
summary_data <- data.frame(
  oneline_Systems = oneline_results_df$Avg_Total_Wait_Time,
  Separate_Systems = separate_results_df$Avg_Wait_Time_Teller1,
  tandem_Systems = tandem_results_df$Avg_Total_Wait_Time
)

# Create a summary table using gtsummary with Nature or Science style formatting
summary_table <- summary_data %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",   # Display mean and standard deviation
    label = list(
      Separate_Systems ~ "Separate Lines",
      oneline_Systems ~ "One Line",
      tandem_Systems ~ "Tandem"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places for clarity
  ) %>%
  modify_header(label = "**System Type**", stat_0 = "**Average Wait Time (hours)**") %>%
  modify_caption("**Summary of Average Wait Times Across Systems**") %>%
  as_gt()

# Display the summary table
summary_table

**Summary of Average Wait Times Across Systems**
System Type	Average Wait Time (hours)¹
One Line	0.61 ± 0.15
Separate Lines	0.75 ± 0.20
Tandem	1.18 ± 0.17
¹ Mean ± SD

The service rate \(\mu\) is an important parameter in the simulation, as it determines the speed at which tellers can serve customers; thereby affecting the wait times and queue lengths in each system.

Case 2: Servers have different rates

Type 1: One line with two servers (Tandem)

In this system, customers first go to Teller 1, then move to Teller 2 once Teller 1’s service is complete. In this system, Teller 1 has a service rate of \(\mu_1 = 5 \text{ customers per hour}\), and Teller 2 has a service rate of \(\mu_2 = 3 \text{ customers per hour}\). The simulation will calculate the average total wait time and average queue length for each teller over a 7-day period.

mu1 <- 5 / 60            # Service rate for Teller 1 (customers per minute)
mu2 <- 3 / 60            # Service rate for Teller 2 (customers per minute)

# Function to simulate 1 line with tandem servers for one day
simulate_tandem_servers_day <- function(arrival_times, mu1, mu2) {

  wait_times_stage1 <- c()
  wait_times_stage2 <- c()
  service_times_stage1 <- c()
  service_times_stage2 <- c()
  
  # Initialize tellers' availability times
  teller1_free_at <- 0
  teller2_free_at <- 0
  
  # Process each arrival in the line with tandem servers
  for (arrival_time in arrival_times) {
    
    # For each customer, first go to Teller 1, then Teller 2
    # Calculate the wait time and service time for each teller
    
    # Stage 1: Teller 1
    wait_time1 <- max(0, teller1_free_at - arrival_time)
    service_time1 <- rexp(1, rate = mu1)
    teller1_free_at <- arrival_time + wait_time1 + service_time1
    wait_times_stage1 <- c(wait_times_stage1, wait_time1)
    service_times_stage1 <- c(service_times_stage1, service_time1)
    
    # Stage 2: Teller 2
    arrival_time_stage2 <- teller1_free_at
    wait_time2 <- max(0, teller2_free_at - arrival_time_stage2)
    service_time2 <- rexp(1, rate = mu2)
    teller2_free_at <- arrival_time_stage2 + wait_time2 + service_time2
    wait_times_stage2 <- c(wait_times_stage2, wait_time2)
    service_times_stage2 <- c(service_times_stage2, service_time2)
  }
  
  # Calculate total wait times
  total_wait_times <- wait_times_stage1 + wait_times_stage2
  avg_wait_time <- mean(total_wait_times)
  
  return(list(avg_wait_time = avg_wait_time, 
              arrival_times = arrival_times, 
              wait_times_stage1 = wait_times_stage1,
              wait_times_stage2 = wait_times_stage2, 
              total_wait_times = total_wait_times))
}

Now, we will run the simulation for 7 days and calculate the average wait time over the week.

# Apply the simulate_tandem_servers_day function to each day's arrival times
tandem_weekly_results <- lapply(arrival_times_by_day, function(arrival_times) {
  simulate_tandem_servers_day(arrival_times, mu1, mu2)
})

# Summary of results
tandem_avg_weekly_wait_time <- mean(tandem_weekly_results %>% 
                               sapply(`[[`, "avg_wait_time"))
cat("Average wait time over 7 days for tandem servers:", 
    tandem_avg_weekly_wait_time / 60, "\n")

Average wait time over 7 days for tandem servers: 2.799473

Finally, we will plot the wait times for each teller and the total wait time for each customer over the week.

Code

# Combine data for all days into a single data frame
all_days_data <- do.call(rbind, lapply(1:7, function(i) {
  data.frame(
    Day = paste("Day", i),
    ArrivalTime = tandem_weekly_results[[i]]$arrival_times,
    WaitTimeStage1 = tandem_weekly_results[[i]]$wait_times_stage1,
    WaitTimeStage2 = tandem_weekly_results[[i]]$wait_times_stage2,
    TotalWaitTime = tandem_weekly_results[[i]]$total_wait_times
  )
}))

# Melt data into a long format for plotting
all_days_data_melted <- melt(all_days_data, id.vars = c("Day", "ArrivalTime"), 
                             variable.name = "Teller", value.name = "WaitTime")

# Custom color palette in the style of Nature or Science journals
custom_colors <- c("WaitTimeStage1" = "#1f78b4",  # Blue
                   "WaitTimeStage2" = "#33a02c",  # Green
                   "TotalWaitTime" = "#e31a1c")   # Red

# Plot using ggplot2 with facet_wrap for shared axes
ggplot(all_days_data_melted, aes(x = ArrivalTime, y = WaitTime, color = Teller)) +
  geom_step(size = 0.8) +
  facet_wrap(~ Day, ncol = 3) +
  labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
       color = "Wait Time Type") +
  scale_color_manual(values = custom_colors, 
                     labels = c("Teller 1", "Teller 2", "Total")) +
  theme_bw() +
  theme(legend.position = "top", 
        plot.title = element_text(size = 10),
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank()) +
  ggtitle("Customer Wait Times Over 7 Days for Tandem Servers") %>%
  suppressMessages()

Since the simulation results can vary due to randomness, we can run the simulation multiple times (1000 times) to get a more accurate estimate of the average wait time over 7 days.

# Perform 1000 simulations over 7 days, for each day's arrival times
tandem_sim_results <- replicate(1000, {
  # Apply simulate_tandem_servers_day for each day's arrival times
  lapply(arrival_times_by_day, function(arrival_times) {
    simulate_tandem_servers_day(arrival_times, mu1, mu2)
  })
}, simplify = FALSE)

# Extract and calculate average wait times for each simulation over 7 days
tandem_avg_total_wait_times <- sapply(tandem_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$total_wait_times)))
})
tandem_avg_wait_times_stage1 <- sapply(tandem_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$wait_times_stage1)))
})
tandem_avg_wait_times_stage2 <- sapply(tandem_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$wait_times_stage2)))
})

# Create a data frame for the summary
tandem_results_df <- data.frame(
  Avg_Total_Wait_Time = tandem_avg_total_wait_times / 60,      # Convert to hours
  Avg_Wait_Time_Server1 = tandem_avg_wait_times_stage1 / 60,   # Convert to hours
  Avg_Wait_Time_Server2 = tandem_avg_wait_times_stage2 / 60    # Convert to hours
)

# Create a summary table with gtsummary
tandem_summary_table <- tandem_results_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    label = list(
      Avg_Total_Wait_Time ~ "Average Total Wait Time (hours)",
      Avg_Wait_Time_Server1 ~ "Average Wait Time for Server 1 (hours)",
      Avg_Wait_Time_Server2 ~ "Average Wait Time for Server 2 (hours)"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places
  ) %>%
  modify_header(label = "**Statistic**") %>%
  modify_caption("**Tandem Service Wait Times Over 7 Days**") 

# Print the summary table
tandem_summary_table

**Tandem Service Wait Times Over 7 Days**
Statistic	N = 1,000¹
Average Total Wait Time (hours)	2.69 ± 0.37
Average Wait Time for Server 1 (hours)	0.70 ± 0.16
Average Wait Time for Server 2 (hours)	1.98 ± 0.38
¹ Mean ± SD

Type 2: One Line, Multiple Servers (Customers Choose)

Since, in this case, the customer will get their food while ordering, the service speed will be slower than the previous case. So, the Teller 1 will have \(\mu_{1_{slower}} = 2.5 \text{ customers per hour}\) and the Teller 2 will have \(\mu_{2_{slower}} = 1.5 \text{ customers per hour}\).

mu1_slower <- 2.5 / 60            # Service rate for Teller 1 (customers per minute)
mu2_slower <- 1.5 / 60            # Service rate for Teller 2 (customers per minute)

# Function to simulate oneline service with 1 line and 2 tellers for one day
simulate_oneline_service_day <- function(arrival_times, mu1, mu2) {
  
  wait_times <- c()
  service_times <- c()
  
  # Initialize tellers' availability times
  teller1_free_at <- 0
  teller2_free_at <- 0
  
  # Process each customer arrival
  for (arrival_time in arrival_times) {
    
    # Determine which teller will serve the customer
    
    # Both tellers are available, go to Teller 1 by default
    if (arrival_time >= teller1_free_at && arrival_time >= teller2_free_at) {
      wait_time <- 0
      service_time <- rexp(1, rate = mu)
      teller1_free_at <- arrival_time + service_time
      
      # Teller 1 is available
    } else if (arrival_time >= teller1_free_at) {
      wait_time <- 0
      service_time <- rexp(1, rate = mu1)
      teller1_free_at <- arrival_time + service_time
      
      # Teller 2 is available
    } else if (arrival_time >= teller2_free_at) {
      wait_time <- 0
      service_time <- rexp(1, rate = mu2)
      teller2_free_at <- arrival_time + service_time
      
      # Both tellers are busy, customer waits in queue and for the next available teller
    } else {
      next_free_time <- min(teller1_free_at, teller2_free_at)
      wait_time <- next_free_time - arrival_time
      
      # Assign to the teller who becomes available first
      if (teller1_free_at <= teller2_free_at) {
        service_time <- rexp(1, rate = mu1)
        teller1_free_at <- next_free_time + service_time
      } else {
        service_time <- rexp(1, rate = mu2)
        teller2_free_at <- next_free_time + service_time
      }
    }
    
    # Record wait and service times
    wait_times <- c(wait_times, wait_time)
    service_times <- c(service_times, service_time)
  }
  
  # Calculate total wait times
  total_wait_times <- wait_times
  avg_wait_time <- mean(total_wait_times)
  
  return(list(avg_wait_time = avg_wait_time, 
              arrival_times = arrival_times, 
              wait_times = wait_times, 
              service_times = service_times,
              total_wait_times = total_wait_times))
}

Now, we will run the simulation for 7 days and calculate the average wait time over the week. Here, in order to keep the consistency, we will use the same arrival times for each day.

# Apply simulate_oneline_service_day function to each day's arrival times
oneline_weekly_results <- lapply(arrival_times_by_day, function(arrival_times) {
  simulate_oneline_service_day(arrival_times, mu1_slower, mu2_slower)
})

# Calculate the average wait time over 7 days
oneline_avg_weekly_wait_time <- mean(sapply(oneline_weekly_results, `[[`, "avg_wait_time"))

# Print the average weekly wait time in hours
cat("Average wait time over 7 days for one line service (hours):", oneline_avg_weekly_wait_time / 60, "\n")

Average wait time over 7 days for one line service (hours): 0.757806

Code

# Combine data for all days into a single data frame
all_days_data <- do.call(rbind, lapply(1:7, function(i) {
  data.frame(
    Day = paste("Day", i),
    ArrivalTime = oneline_weekly_results[[i]]$arrival_times,
    WaitTime = oneline_weekly_results[[i]]$total_wait_times
  )
}))

# Plot for all_days_data using ggplot2 with facet_wrap
ggplot(all_days_data, aes(x = ArrivalTime, y = WaitTime)) +
  geom_step(color = "#e31a1c", size = 0.8) +
  facet_wrap(~ Day, ncol = 3) +
  labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
       title = "Wait Times for One Line Service Over 7 Days") +
  theme_bw(base_size = 12) +
  theme(
    legend.position = "none",  # No legend needed for single color
    panel.grid.major = element_blank(),  # Remove major grid lines
    panel.grid.minor = element_blank(),  # Remove minor grid lines
    strip.text = element_text(size = 10),  # Bold facet labels
    plot.title = element_text(size = 14, hjust = 0.5),  # Centered and bold title
    axis.text = element_text(size = 10),  # Adjust axis text size
    axis.title = element_text(size = 12),  # Bold axis titles
    axis.text.x = element_text(angle = 45, hjust = 1)  # Rotate x-axis labels
  ) %>%
  suppressMessages()

Since the simulation results can vary due to randomness, we can run the simulation multiple times (1000 times) to get a more accurate estimate of the average wait time over 7 days.

# Perform 1000 simulations over 7 days, using each day's arrival times
oneline_sim_results <- replicate(1000, {
  # Apply simulate_oneline_service_day for each day's arrival times
  lapply(arrival_times_by_day, function(arrival_times) {
    simulate_oneline_service_day(arrival_times, mu1_slower, mu2_slower)
  })
}, simplify = FALSE)

# Extract and calculate average total wait times for each simulation over 7 days
oneline_avg_total_wait_times <- sapply(oneline_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$total_wait_times)))
})

# Convert average total wait times from minutes to hours
oneline_avg_total_wait_times_hours <- oneline_avg_total_wait_times / 60

# Create a data frame for the summary
oneline_results_df <- data.frame(
  Avg_Total_Wait_Time = oneline_avg_total_wait_times_hours  # Average wait time in hours
)

# Create a summary table with gtsummary
oneline_summary_table <- oneline_results_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    label = list(
      Avg_Total_Wait_Time ~ "Average Total Wait Time (hours)"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places
  ) %>%
  modify_header(label = "**Statistic**") %>%
  modify_caption("**One Line Service Wait Times Over 7 Days**") 

# Print the summary table
oneline_summary_table

**One Line Service Wait Times Over 7 Days**
Statistic	N = 1,000¹
Average Total Wait Time (hours)	1.00 ± 0.21
¹ Mean ± SD

Type 3: Each teller has their own line (grocery store)

In this scenario, each teller has their own line, and customers choose the shortest line to join. This is similar to the setup in a grocery store, where each checkout lane has its own queue.

The Teller 1 will have \(\mu_{1{\text{slower}}} = 2.5 \text{ customers per hour}\) and Teller 2 will have \(\mu_{2{\text{slower}}} = 1.5 \text{ customers per hour}\).

mu1_slower <- 2.5 / 60  # Service rate for Teller 1 (slower)
mu2_slower <- 1.5 / 60  # Service rate for Teller 2 (slower)

# Function to simulate service with each teller having their own line for one day
simulate_separate_lines_day <- function(arrival_times, mu1, mu2) {
  
  # Initialize variables for each teller
  wait_times_teller1 <- c()
  wait_times_teller2 <- c()
  service_times_teller1 <- c()
  service_times_teller2 <- c()
  
  # Initialize tellers' availability times and queue lengths
  teller1_free_at <- 0
  teller2_free_at <- 0
  queue_length_teller1 <- 0
  queue_length_teller2 <- 0
  
  # Initialize an empty data frame to store results
  results <- data.frame(
    ArrivalTime = numeric(0), 
    WaitTime = numeric(0), 
    ServiceTime = numeric(0), 
    Teller = character(0),
    QueueLength_Teller1 = numeric(0),
    QueueLength_Teller2 = numeric(0)
  )
  
  # Process each customer arrival
  for (arrival_time in arrival_times) {
    
    # Record the current queue lengths
    current_queue_length1 <- queue_length_teller1
    current_queue_length2 <- queue_length_teller2
    
    # Determine which teller the customer will go to based on the shorter queue
    if (queue_length_teller1 < queue_length_teller2 || 
        (queue_length_teller1 == queue_length_teller2 && runif(1) < 0.5)) {
      
      # Assign to Teller 1
      if (arrival_time >= teller1_free_at) {
        # Teller 1 is free, no wait time
        wait_time <- 0
        service_time <- rexp(1, rate = mu1)
        teller1_free_at <- arrival_time + service_time
      } else {
        # Teller 1 is busy, customer waits
        wait_time <- teller1_free_at - arrival_time
        service_time <- rexp(1, rate = mu2)
        teller1_free_at <- teller1_free_at + service_time
        queue_length_teller1 <- queue_length_teller1 + 1
      }
      
      # Record the result for Teller 1
      results <- rbind(results, data.frame(
        ArrivalTime = arrival_time, 
        WaitTime = wait_time, 
        ServiceTime = service_time, 
        Teller = "Teller 1",
        QueueLength_Teller1 = current_queue_length1,
        QueueLength_Teller2 = current_queue_length2
      ))
      
    } else {
      
      # Assign to Teller 2
      if (arrival_time >= teller2_free_at) {
        # Teller 2 is free, no wait time
        wait_time <- 0
        service_time <- rexp(1, rate = mu1)
        teller2_free_at <- arrival_time + service_time
      } else {
        # Teller 2 is busy, customer waits
        wait_time <- teller2_free_at - arrival_time
        service_time <- rexp(1, rate = mu2)
        teller2_free_at <- teller2_free_at + service_time
        queue_length_teller2 <- queue_length_teller2 + 1
      }
      
      # Record the result for Teller 2
      results <- rbind(results, data.frame(
        ArrivalTime = arrival_time, 
        WaitTime = wait_time, 
        ServiceTime = service_time, 
        Teller = "Teller 2",
        QueueLength_Teller1 = current_queue_length1,
        QueueLength_Teller2 = current_queue_length2
      ))
      
    }
    
    # Update queue lengths after the patient is served (if applicable)
    if (arrival_time >= teller1_free_at) {
      queue_length_teller1 <- max(0, queue_length_teller1 - 1)
    }
    if (arrival_time >= teller2_free_at) {
      queue_length_teller2 <- max(0, queue_length_teller2 - 1)
    }
  }
  
  # Calculate average wait times for each teller
  avg_wait_time_teller1 <- mean(results$WaitTime[results$Teller == "Teller 1"])
  avg_wait_time_teller2 <- mean(results$WaitTime[results$Teller == "Teller 2"])
  
  return(list(
    results = results,
    avg_wait_time_teller1 = avg_wait_time_teller1,
    avg_wait_time_teller2 = avg_wait_time_teller2
  ))
}

Now, we will run the simulation for 7 days and calculate the average wait time over the week. Here, in order to keep the consistency, we will use the same arrival times for each day.

# Apply simulate_separate_lines_day function to each day's arrival times
separate_weekly_results <- lapply(arrival_times_by_day, function(arrival_times) {
  simulate_separate_lines_day(arrival_times, mu1_slower, mu2_slower)
})

# Calculate the average weekly wait time for each teller
avg_weekly_wait_time_teller1 <- mean(sapply(separate_weekly_results, `[[`, "avg_wait_time_teller1"))
avg_weekly_wait_time_teller2 <- mean(sapply(separate_weekly_results, `[[`, "avg_wait_time_teller2"))

# Convert from minutes to hours for reporting
avg_weekly_wait_time_teller1_hours <- avg_weekly_wait_time_teller1 / 60
avg_weekly_wait_time_teller2_hours <- avg_weekly_wait_time_teller2 / 60

# Print the average weekly wait times for each teller
cat("Average weekly wait time for Teller 1 (hours):", round(avg_weekly_wait_time_teller1_hours, 2), "\n")

Average weekly wait time for Teller 1 (hours): 1.43

cat("Average weekly wait time for Teller 2 (hours):", round(avg_weekly_wait_time_teller2_hours, 2), "\n")

Average weekly wait time for Teller 2 (hours): 2.08

Plotting the distribution of wait times for each teller over the week.

Code

# Combine data for all days into a single data frame for both tellers
all_days_data <- do.call(rbind, lapply(1:7, function(i) {
  day_data <- separate_weekly_results[[i]]$results
  day_data$Day <- paste("Day", i) 
  return(day_data)
}))

# Convert 'Day' column to a factor for ordering in facets
all_days_data$Day <- factor(all_days_data$Day, levels = paste("Day", 1:7))

# Plot using ggplot2 with facet_wrap
ggplot(all_days_data, aes(x = ArrivalTime, y = WaitTime, color = Teller)) +
  geom_step(size = 0.8) +
  facet_wrap(~ Day, ncol = 3) +
  labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
       title = "Wait Times for Separate Lines Over 7 Days") +
  scale_color_manual(values = c("Teller 1" = "#1f78b4", "Teller 2" = "#33a02c")) +  # Colors for tellers
  theme_bw(base_size = 12) +
  theme(
    legend.position = "top", 
    panel.grid.major = element_blank(),  # Remove major grid lines
    panel.grid.minor = element_blank(),  # Remove minor grid lines
    strip.text = element_text(size = 10),  # Bold facet labels
    plot.title = element_text(size = 14, hjust = 0.5),  # Centered and bold title
    axis.text = element_text(size = 10),  # Adjust axis text size
    axis.title = element_text(size = 12),  # Bold axis titles
    axis.text.x = element_text(angle = 45, hjust = 1)  # Rotate x-axis labels
  ) %>%
  suppressMessages()

Finally, since the simulation is stochastic, we can run it multiple times to get a sense of the variability in the results.

# Perform 100 simulations over 7 days, using each day's arrival times
separate_sim_results <- replicate(200, {
  # Apply simulate_separate_lines_day for each day's arrival times
  lapply(arrival_times_by_day, function(arrival_times) {
    simulate_separate_lines_day(arrival_times, mu1_slower, mu2_slower)
  })
}, simplify = FALSE)

# Extract and calculate average wait times and queue lengths for each simulation over 7 days
separate_avg_wait_time_teller1 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$WaitTime[day_result$results$Teller == "Teller 1"])))
})

separate_avg_wait_time_teller2 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$WaitTime[day_result$results$Teller == "Teller 2"])))
})

separate_avg_queue_length_teller1 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$QueueLength_Teller1[day_result$results$Teller == "Teller 1"])))
})

separate_avg_queue_length_teller2 <- sapply(separate_sim_results, function(simulation) {
  mean(sapply(simulation, function(day_result) mean(day_result$results$QueueLength_Teller2[day_result$results$Teller == "Teller 2"])))
})

# Convert average wait times from minutes to hours
separate_avg_wait_time_teller1_hours <- separate_avg_wait_time_teller1 / 60
separate_avg_wait_time_teller2_hours <- separate_avg_wait_time_teller2 / 60

# Create a data frame for the summary
separate_results_df <- data.frame(
  Avg_Wait_Time_Teller1 = separate_avg_wait_time_teller1_hours,        # Average wait time for Teller 1 in hours
  Avg_Wait_Time_Teller2 = separate_avg_wait_time_teller1_hours,        # Average wait time for Teller 2 in hours
  Avg_Queue_Length_Teller1 = separate_avg_queue_length_teller1,        # Average queue length for Teller 1
  Avg_Queue_Length_Teller2 = separate_avg_queue_length_teller2         # Average queue length for Teller 2
)

# Create a summary table with gtsummary
separate_summary_table <- separate_results_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    label = list(
      Avg_Wait_Time_Teller1 ~ "Average Wait Time for Teller 1 (hours)",
      Avg_Wait_Time_Teller2 ~ "Average Wait Time for Teller 2 (hours)",
      Avg_Queue_Length_Teller1 ~ "Average Queue Length for Teller 1",
      Avg_Queue_Length_Teller2 ~ "Average Queue Length for Teller 2"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places
  ) %>%
  modify_header(label = "**Statistic**") %>%
  modify_caption("**Separate Line Wait Times and Queue Lengths Over 7 Days**")

# Print the summary table
separate_summary_table

**Separate Line Wait Times and Queue Lengths Over 7 Days**
Statistic	N = 200¹
Average Wait Time for Teller 1 (hours)	2.13 ± 0.45
Average Wait Time for Teller 2 (hours)	2.13 ± 0.45
Average Queue Length for Teller 1	5.42 ± 0.46
Average Queue Length for Teller 2	5.41 ± 0.45
¹ Mean ± SD

Conclusion

Assumptions:

In this simulation, I assume that the service rate for the Tandem system, teller 1 has a service rate of \(\mu = 5 \text{ customers per hour}\), teller 2 has a service rate of \(\mu = 3 \text{ customers per hour}\), and they serve customers simultaneously. While for the Separate Lines and One Line System systems, teller 1 has a service rate of \(\mu = 2.5 \text{ customers per hour}\), and teller 2 has a service rate of \(\mu = 2.5 \text{ customers per hour}\), and they must finish serving one customer before starting the next.

System Type	Service Rate (μ) per Teller	Description
Tandem Line	Teller 1: 5 customers per hour Teller 2: 3 customers per hour	Teller 1 and Teller 2 serve customers simultaneously, each with their respective service rates.
Separate Lines	Teller 1: 2.5 customers per hour Teller 2: 2.5 customers per hour	Each teller has a service rate of 2.5 customers per hour and must finish serving one customer before starting the next.
One Line System	Teller 1: 2.5 customers per hour Teller 2: 2.5 customers per hour	Similar to Separate Lines, with each teller having a service rate of 2.5 customers per hour, but one teller must finish before the next teller can start.

# Combine the data frames into a single data frame for a cohesive summary
summary_data <- data.frame(
  oneline_Systems = oneline_results_df$Avg_Total_Wait_Time,
  Separate_Systems = separate_results_df$Avg_Wait_Time_Teller1,
  tandem_Systems = tandem_results_df$Avg_Total_Wait_Time
)

# Create a summary table using gtsummary with Nature or Science style formatting
summary_table <- summary_data %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",   # Display mean and standard deviation
    label = list(
      Separate_Systems ~ "Separate Lines",
      oneline_Systems ~ "One Line",
      tandem_Systems ~ "Tandem"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places for clarity
  ) %>%
  modify_header(label = "**System Type**", stat_0 = "**Average Wait Time (hours)**") %>%
  modify_caption("**Summary of Average Wait Times Across Systems**") %>%
  as_gt()

# Display the summary table
summary_table

**Summary of Average Wait Times Across Systems**
System Type	Average Wait Time (hours)¹
One Line	1.00 ± 0.21
Separate Lines	2.13 ± 0.45
Tandem	2.69 ± 0.37
¹ Mean ± SD

Based on the simulation results, the One Line System has the lowest average total wait time, followed by the Separate Lines system and the tandem Line system. The One Line System also has the lowest average queue length for both tellers, indicating efficient utilization of resources.

Code

# Combine the data frames into a single data frame for a cohesive summary
summary_data <- data.frame(
  oneline_Systems = oneline_results_df$Avg_Total_Wait_Time,
  Separate_Systems = separate_results_df$Avg_Wait_Time_Teller1,
  tandem_Systems = tandem_results_df$Avg_Total_Wait_Time
)

# Create a summary table using gtsummary with Nature or Science style formatting
summary_table <- summary_data %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",   # Display mean and standard deviation
    label = list(
      Separate_Systems ~ "Separate Lines",
      oneline_Systems ~ "One Line",
      tandem_Systems ~ "Tandem"
    ),
    digits = all_continuous() ~ 2  # Round to 2 decimal places for clarity
  ) %>%
  modify_header(label = "**System Type**", stat_0 = "**Average Wait Time (hours)**") %>%
  modify_caption("**Summary of Average Wait Times Across Systems**") %>%
  as_gt()

# Display the summary table
summary_table

**Summary of Average Wait Times Across Systems**
System Type	Average Wait Time (hours)¹
One Line	1.00 ± 0.21
Separate Lines	2.13 ± 0.45
Tandem	2.69 ± 0.37
¹ Mean ± SD

Case 3: Three Servers with Different Service Rates

Permutations of Service Rates

There are \(3! = 6\) possible permutations of service rates for the three servers: Fast, Medium, and Slow. The service rates are as follows:

# Define service rates for Fast, Medium, and Slow servers
fast_rate <- 7 / 60    # Fast server: 7 customers per hour
medium_rate <- 5 / 60  # Medium server: 5 customers per hour
slow_rate <- 3 / 60    # Slow server: 3 customers per hour
fast_rate_half <- fast_rate / 2  # Half the rate of the fast server
medium_rate_half <- medium_rate / 2  # Half the rate of the medium server
slow_rate_half <- slow_rate / 2  # Half the rate of the slow server

# Generate all permutations of service rates
service_rate_permutations <- list(
  c(fast_rate, medium_rate, slow_rate),
  c(fast_rate, slow_rate, medium_rate),
  c(medium_rate, fast_rate, slow_rate),
  c(medium_rate, slow_rate, fast_rate),
  c(slow_rate, fast_rate, medium_rate),
  c(slow_rate, medium_rate, fast_rate)
)

service_rate_permutations_half <- list(
  c(fast_rate_half, medium_rate_half, slow_rate_half),
  c(fast_rate_half, slow_rate_half, medium_rate_half),
  c(medium_rate_half, fast_rate_half, slow_rate_half),
  c(medium_rate_half, slow_rate_half, fast_rate_half),
  c(slow_rate_half, fast_rate_half, medium_rate_half),
  c(slow_rate_half, medium_rate_half, fast_rate_half)
)

# Define descriptive names for each permutation based on rate order
rate_labels <- c("fast-medium-slow", "fast-slow-medium", "medium-fast-slow", 
                 "medium-slow-fast", "slow-fast-medium", "slow-medium-fast")

Type 1: One line with two servers (Tandem)

In this system, customers first go to Teller 1, then move to Teller 2 once Teller 1’s service is complete. In this system, Teller 1 has a service rate of \(\mu_1 = 7 \text{ customers per hour}\), Teller 2 has a service rate of \(\mu_2 = 5 \text{ customers per hour}\), and Teller 3 has a service rate of \(\mu_3 = 3 \text{ customers per hour}\). The simulation will calculate the average total wait time and average queue length for each teller over a 7-day period.

# Function to simulate 1 line with tandem servers for one day with three tellers
simulate_tandem_servers_day <- function(arrival_times, mu1, mu2, mu3) {

  wait_times_stage1 <- c()
  wait_times_stage2 <- c()
  wait_times_stage3 <- c()
  service_times_stage1 <- c()
  service_times_stage2 <- c()
  service_times_stage3 <- c()
  
  # Initialize tellers' availability times
  teller1_free_at <- 0
  teller2_free_at <- 0
  teller3_free_at <- 0
  
  # Process each arrival in the line with tandem servers
  for (arrival_time in arrival_times) {
    
    # Stage 1: Teller 1
    wait_time1 <- max(0, teller1_free_at - arrival_time)
    service_time1 <- rexp(1, rate = mu1)
    teller1_free_at <- arrival_time + wait_time1 + service_time1
    wait_times_stage1 <- c(wait_times_stage1, wait_time1)
    service_times_stage1 <- c(service_times_stage1, service_time1)
    
    # Stage 2: Teller 2
    arrival_time_stage2 <- teller1_free_at
    wait_time2 <- max(0, teller2_free_at - arrival_time_stage2)
    service_time2 <- rexp(1, rate = mu2)
    teller2_free_at <- arrival_time_stage2 + wait_time2 + service_time2
    wait_times_stage2 <- c(wait_times_stage2, wait_time2)
    service_times_stage2 <- c(service_times_stage2, service_time2)
    
    # Stage 3: Teller 3
    arrival_time_stage3 <- teller2_free_at
    wait_time3 <- max(0, teller3_free_at - arrival_time_stage3)
    service_time3 <- rexp(1, rate = mu3)
    teller3_free_at <- arrival_time_stage3 + wait_time3 + service_time3
    wait_times_stage3 <- c(wait_times_stage3, wait_time3)
    service_times_stage3 <- c(service_times_stage3, service_time3)
  }
  
  # Calculate total wait times across all stages
  total_wait_times <- wait_times_stage1 + wait_times_stage2 + wait_times_stage3
  avg_wait_time <- mean(total_wait_times)
  
  return(list(
    avg_wait_time = avg_wait_time, 
    arrival_times = arrival_times, 
    wait_times_stage1 = wait_times_stage1,
    wait_times_stage2 = wait_times_stage2,
    wait_times_stage3 = wait_times_stage3,
    service_times_stage1 = service_times_stage1,
    service_times_stage2 = service_times_stage2,
    service_times_stage3 = service_times_stage3,
    total_wait_times = total_wait_times
  ))
}

Function to perform the simulation, create step plots, and generate results data frames

run_permutation_simulation <- function(arrival_times_by_day, permutations, 
                                       rate_labels, num_simulations = 1000) {

  results_list <- list()
  
  for (i in seq_along(permutations)) {
    rates <- permutations[[i]]
    mu1 <- rates[1]
    mu2 <- rates[2]
    mu3 <- rates[3]
    
    # Run 1000 simulations over 7 days for the current permutation
    tandem_sim_results <- replicate(num_simulations, {
      lapply(arrival_times_by_day, function(arrival_times) {
        simulate_tandem_servers_day(arrival_times, mu1, mu2, mu3)
      })
    }, simplify = FALSE)
    
    # Extract average wait times for each simulation
    avg_total_wait_times <- sapply(tandem_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$total_wait_times)))
    })
    avg_wait_times_stage1 <- sapply(tandem_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$wait_times_stage1)))
    })
    avg_wait_times_stage2 <- sapply(tandem_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$wait_times_stage2)))
    })
    avg_wait_times_stage3 <- sapply(tandem_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$wait_times_stage3)))
    })
    
    # Create a data frame for the summary
    results_df <- data.frame(
      Avg_Total_Wait_Time = avg_total_wait_times / 60,      # Convert to hours
      Avg_Wait_Time_Server1 = avg_wait_times_stage1 / 60,   # Convert to hours
      Avg_Wait_Time_Server2 = avg_wait_times_stage2 / 60,   # Convert to hours
      Avg_Wait_Time_Server3 = avg_wait_times_stage3 / 60    # Convert to hours
    )
    
    # Collect data for step plot for the first simulation run of the 7 days
    first_run_data <- do.call(rbind, lapply(1:7, function(day) {
      data.frame(
        Day = paste("Day", day),
        ArrivalTime = tandem_sim_results[[1]][[day]]$arrival_times,
        WaitTimeStage1 = tandem_sim_results[[1]][[day]]$wait_times_stage1,
        WaitTimeStage2 = tandem_sim_results[[1]][[day]]$wait_times_stage2,
        WaitTimeStage3 = tandem_sim_results[[1]][[day]]$wait_times_stage3,
        TotalWaitTime = tandem_sim_results[[1]][[day]]$total_wait_times
      )
    }))
    
    # Melt data into long format for plotting
    first_run_data_melted <- melt(first_run_data, id.vars = c("Day", "ArrivalTime"), 
                                  variable.name = "Teller", value.name = "WaitTime")
    
    # Custom color palette
    custom_colors <- c("WaitTimeStage1" = "#1f78b4", 
                       "WaitTimeStage2" = "#33a02c", 
                       "WaitTimeStage3" = "#ff7f00", 
                       "TotalWaitTime" = "#e31a1c")
    
    # Generate step plot
    step_plot <- ggplot(first_run_data_melted, aes(x = ArrivalTime, y = WaitTime, color = Teller)) +
      geom_step(size = 0.8) +
      facet_wrap(~ Day, ncol = 3) +
      labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
           color = "Wait Time Type") +
      scale_color_manual(values = custom_colors, 
                         labels = c("Teller 1", "Teller 2", "Teller 3", "Total Wait Time")) +
      theme_bw(base_size = 12) +
      theme(
        legend.position = "top", 
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 9),
        strip.text = element_text(size = 10),  # Bold facet labels
        plot.title = element_text(size = 14, hjust = 0.5),  # Centered and bold title
        axis.text = element_text(size = 10),  # Adjust axis text size
        axis.title = element_text(size = 12),  # Bold axis titles
        axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis labels
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank()
      ) +
      ggtitle(paste("Customer Wait Times for", rate_labels[i], "Over 7 Days"))
    
    # Store the results
    results_list[[rate_labels[i]]] <- 
      list(
      results_df = results_df,
      step_plot = step_plot
    )
  }
  
  return(results_list)
}

# Run the function for all permutations
tandem_permutation_results <- run_permutation_simulation(arrival_times_by_day, 
                                                           service_rate_permutations,
                                                           rate_labels, num_simulations = 1000)

# Loop to access and display results
for (perm_name in names(tandem_permutation_results)) {
  cat("\n", perm_name, "Summary Table:\n")
  
  # Apply gtsummary to the results_df for a detailed summary table
  summary_table <- tandem_permutation_results[[perm_name]]$results_df %>%
    tbl_summary(
      statistic = all_continuous() ~ "{mean} ± {sd}",
      label = list(
        Avg_Total_Wait_Time ~ "Average Total Wait Time (hours)",
        Avg_Wait_Time_Server1 ~ "Average Wait Time for Server 1 (hours)",
        Avg_Wait_Time_Server2 ~ "Average Wait Time for Server 2 (hours)",
        Avg_Wait_Time_Server3 ~ "Average Wait Time for Server 3 (hours)"
      ),
      digits = all_continuous() ~ 2
    ) %>%
    modify_header(label = "**Statistic**") %>%
    modify_caption(paste("**Tandem Service Wait Times for", perm_name, "**"))
  
  # Extract the gtsummary table as a data frame
  summary_df <- as_tibble(summary_table)
  
  # Print the data frame in a clean format
  print(summary_df)
  
  print(tandem_permutation_results[[perm_name]]$step_plot)
}


 fast-medium-slow Summary Table:
# A tibble: 4 × 2
  `**Statistic**`                        `**N = 1,000**`
  <chr>                                  <chr>          
1 Average Total Wait Time (hours)        2.79 ± 0.34    
2 Average Wait Time for Server 1 (hours) 0.26 ± 0.07    
3 Average Wait Time for Server 2 (hours) 0.63 ± 0.15    
4 Average Wait Time for Server 3 (hours) 1.90 ± 0.35


 fast-slow-medium Summary Table:
# A tibble: 4 × 2
  `**Statistic**`                        `**N = 1,000**`
  <chr>                                  <chr>          
1 Average Total Wait Time (hours)        2.78 ± 0.35    
2 Average Wait Time for Server 1 (hours) 0.26 ± 0.07    
3 Average Wait Time for Server 2 (hours) 2.30 ± 0.37    
4 Average Wait Time for Server 3 (hours) 0.22 ± 0.07


 medium-fast-slow Summary Table:
# A tibble: 4 × 2
  `**Statistic**`                        `**N = 1,000**`
  <chr>                                  <chr>          
1 Average Total Wait Time (hours)        2.76 ± 0.35    
2 Average Wait Time for Server 1 (hours) 0.70 ± 0.16    
3 Average Wait Time for Server 2 (hours) 0.17 ± 0.05    
4 Average Wait Time for Server 3 (hours) 1.88 ± 0.36


 medium-slow-fast Summary Table:
# A tibble: 4 × 2
  `**Statistic**`                        `**N = 1,000**`
  <chr>                                  <chr>          
1 Average Total Wait Time (hours)        2.79 ± 0.35    
2 Average Wait Time for Server 1 (hours) 0.71 ± 0.16    
3 Average Wait Time for Server 2 (hours) 2.00 ± 0.35    
4 Average Wait Time for Server 3 (hours) 0.09 ± 0.03


 slow-fast-medium Summary Table:
# A tibble: 4 × 2
  `**Statistic**`                        `**N = 1,000**`
  <chr>                                  <chr>          
1 Average Total Wait Time (hours)        2.78 ± 0.34    
2 Average Wait Time for Server 1 (hours) 2.47 ± 0.37    
3 Average Wait Time for Server 2 (hours) 0.09 ± 0.03    
4 Average Wait Time for Server 3 (hours) 0.22 ± 0.07


 slow-medium-fast Summary Table:
# A tibble: 4 × 2
  `**Statistic**`                        `**N = 1,000**`
  <chr>                                  <chr>          
1 Average Total Wait Time (hours)        2.77 ± 0.34    
2 Average Wait Time for Server 1 (hours) 2.46 ± 0.36    
3 Average Wait Time for Server 2 (hours) 0.23 ± 0.07    
4 Average Wait Time for Server 3 (hours) 0.09 ± 0.03

Results

# Extract the average wait times for each permutation and combine them into a single data frame
tandem_summary_df <- tibble(
  `fast-medium-slow` = tandem_permutation_results$`fast-medium-slow`$results_df$Avg_Total_Wait_Time,
  `fast-slow-medium` = tandem_permutation_results$`fast-slow-medium`$results_df$Avg_Total_Wait_Time,
  `medium-fast-slow` = tandem_permutation_results$`medium-fast-slow`$results_df$Avg_Total_Wait_Time,
  `medium-slow-fast` = tandem_permutation_results$`medium-slow-fast`$results_df$Avg_Total_Wait_Time,
  `slow-fast-medium` = tandem_permutation_results$`slow-fast-medium`$results_df$Avg_Total_Wait_Time,
  `slow-medium-fast` = tandem_permutation_results$`slow-medium-fast`$results_df$Avg_Total_Wait_Time
)

# Create a summary table using gtsummary
tandem_summary_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",   # Display mean and standard deviation
    digits = all_continuous() ~ 2,  # Round to 2 decimal places for clarity
  ) %>%
  modify_header(label = "**Permutation**", stat_0 = "**Average Wait Time (hours)**") %>%
  modify_caption("**Summary of Average Wait Times Across Permutations (Tandem)**")

**Summary of Average Wait Times Across Permutations (Tandem)**
Permutation	Average Wait Time (hours)¹
fast-medium-slow	2.79 ± 0.34
fast-slow-medium	2.78 ± 0.35
medium-fast-slow	2.76 ± 0.35
medium-slow-fast	2.79 ± 0.35
slow-fast-medium	2.78 ± 0.34
slow-medium-fast	2.77 ± 0.34
¹ Mean ± SD

Type 2: One Line, Multiple Servers (Customers Choose)

In this case, there is also only one queue. Customer will go to either server 1, server 2, or server 3. And customers will get their food while ordering. For example, in an In-N-Out restaurant, and customers go to the first server by default to place their order if both servers are available. If one server is busy, the customer will go to the other server. If both servers are busy, the customer will wait in the queue and go to the first server when available.

Since, in this case, the customer will get their food while ordering, the service speed will be slower than the previous case. So, the Teller 1 will have \(\mu_{1_{slower}} = 3.5 \text{ customers per hour}\) and the Teller 2 will have \(\mu_{2_{slower}} = 2.5 \text{ customers per hour}\), and the Teller 3 will have \(\mu_{3_{slower}} = 1 \text{ customers per hour}\).

# Function to simulate oneline service with three servers for one day
simulate_oneline_service_day <- function(arrival_times, mu1, mu2, mu3) {
  
  wait_times <- c()
  service_times <- c()
  
  # Initialize tellers' availability times
  teller1_free_at <- 0
  teller2_free_at <- 0
  teller3_free_at <- 0
  
  # Process each customer arrival
  for (arrival_time in arrival_times) {
    
    # Both tellers are available, go to Teller 1 by default
    if (arrival_time >= teller1_free_at && arrival_time >= teller2_free_at && arrival_time >= teller3_free_at) {
      wait_time <- 0
      service_time <- rexp(1, rate = mu1)
      teller1_free_at <- arrival_time + service_time
      
    } else if (arrival_time >= teller1_free_at) {
      # Teller 1 is available
      wait_time <- 0
      service_time <- rexp(1, rate = mu1)
      teller1_free_at <- arrival_time + service_time
      
    } else if (arrival_time >= teller2_free_at) {
      # Teller 2 is available
      wait_time <- 0
      service_time <- rexp(1, rate = mu2)
      teller2_free_at <- arrival_time + service_time
      
    } else if (arrival_time >= teller3_free_at) {
      # Teller 3 is available
      wait_time <- 0
      service_time <- rexp(1, rate = mu3)
      teller3_free_at <- arrival_time + service_time
      
    } else {
      # All tellers are busy, customer waits for the next available teller
      next_free_time <- min(teller1_free_at, teller2_free_at, teller3_free_at)
      wait_time <- next_free_time - arrival_time
      
      # Assign to the teller who becomes available first
      if (teller1_free_at <= teller2_free_at && teller1_free_at <= teller3_free_at) {
        service_time <- rexp(1, rate = mu1)
        teller1_free_at <- next_free_time + service_time
      } else if (teller2_free_at <= teller3_free_at) {
        service_time <- rexp(1, rate = mu2)
        teller2_free_at <- next_free_time + service_time
      } else {
        service_time <- rexp(1, rate = mu3)
        teller3_free_at <- next_free_time + service_time
      }
    }
    
    # Record wait and service times
    wait_times <- c(wait_times, wait_time)
    service_times <- c(service_times, service_time)
  }
  
  # Calculate total wait times
  total_wait_times <- wait_times
  avg_wait_time <- mean(total_wait_times)
  
  return(list(avg_wait_time = avg_wait_time, 
              arrival_times = arrival_times, 
              wait_times = wait_times, 
              service_times = service_times,
              total_wait_times = total_wait_times))
}

# Now, run the simulation for 7 days with half rates for each server
oneline_weekly_results <- lapply(arrival_times_by_day, function(arrival_times) {
  simulate_oneline_service_day(arrival_times, fast_rate_half, medium_rate_half, slow_rate_half)
})

# Calculate the average wait time over 7 days
oneline_avg_weekly_wait_time <- mean(sapply(oneline_weekly_results, `[[`, "avg_wait_time"))

# Print the average weekly wait time in hours
cat("Average wait time over 7 days for one line service (hours):", oneline_avg_weekly_wait_time / 60, "\n")

Average wait time over 7 days for one line service (hours): 0.07955179

# Function to run permutation simulation for oneline service
run_oneline_permutation_simulation <- function(arrival_times_by_day, permutations, 
                                              rate_labels, num_simulations = 1000) {
  results_list <- list()
  
  for (i in seq_along(permutations)) {
    rates <- permutations[[i]]
    mu1 <- rates[1]
    mu2 <- rates[2]
    mu3 <- rates[3]
    
    # Run simulations for the current permutation
    oneline_sim_results <- replicate(num_simulations, {
      lapply(arrival_times_by_day, function(arrival_times) {
        simulate_oneline_service_day(arrival_times, mu1, mu2, mu3)
      })
    }, simplify = FALSE)
    
    # Extract average wait times for each simulation
    avg_total_wait_times <- sapply(oneline_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$total_wait_times)))
    })
    
    # Create a data frame for the summary
    results_df <- data.frame(
      Avg_Total_Wait_Time = avg_total_wait_times / 60  # Convert to hours
    )
    
    # Collect data for step plot for the first simulation run of the 7 days
    first_run_data <- do.call(rbind, lapply(1:7, function(day) {
      data.frame(
        Day = paste("Day", day),
        ArrivalTime = oneline_sim_results[[1]][[day]]$arrival_times,
        WaitTime = oneline_sim_results[[1]][[day]]$total_wait_times
      )
    }))
    
    # Generate step plot for the first simulation run
    step_plot <- ggplot(first_run_data, aes(x = ArrivalTime, y = WaitTime)) +
      geom_step(color = "#e31a1c", size = 0.8) +
      facet_wrap(~ Day, ncol = 3) +
      labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
           title = paste("Wait Times for", rate_labels[i], "Over 7 Days")) +
      theme_bw(base_size = 12) +
      theme(
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text = element_text(size = 10),
        plot.title = element_text(size = 14, hjust = 0.5),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )
    
    # Store the results
    results_list[[rate_labels[i]]] <- list(
      results_df = results_df,
      step_plot = step_plot
    )
  }
  
  return(results_list)
}

# Run the function for all permutations
oneline_permutation_results <- run_oneline_permutation_simulation(arrival_times_by_day, 
                                                                service_rate_permutations_half,
                                                                rate_labels, num_simulations = 1000)

# Loop to access and display results
for (perm_name in names(oneline_permutation_results)) {
  cat("\n", perm_name, "Summary Table:\n")
  
  # Apply gtsummary to the results_df for a detailed summary table
  summary_table <- oneline_permutation_results[[perm_name]]$results_df %>%
    tbl_summary(
      statistic = all_continuous() ~ "{mean} ± {sd}",
      label = list(
        Avg_Total_Wait_Time ~ "Average Total Wait Time (hours)"
      ),
      digits = all_continuous() ~ 2
    ) %>%
    modify_header(label = "**Statistic**") %>%
    modify_caption(paste("**One Line Service Wait Times for", perm_name, "**"))
  
  # Extract the gtsummary table as a data frame
  summary_df <- as_tibble(summary_table)
  
  # Print the data frame in a clean format
  print(summary_df)
  print(oneline_permutation_results[[perm_name]]$step_plot)
}


 fast-medium-slow Summary Table:
# A tibble: 1 × 2
  `**Statistic**`                 `**N = 1,000**`
  <chr>                           <chr>          
1 Average Total Wait Time (hours) 0.12 ± 0.05


 fast-slow-medium Summary Table:
# A tibble: 1 × 2
  `**Statistic**`                 `**N = 1,000**`
  <chr>                           <chr>          
1 Average Total Wait Time (hours) 0.13 ± 0.05


 medium-fast-slow Summary Table:
# A tibble: 1 × 2
  `**Statistic**`                 `**N = 1,000**`
  <chr>                           <chr>          
1 Average Total Wait Time (hours) 0.12 ± 0.05


 medium-slow-fast Summary Table:
# A tibble: 1 × 2
  `**Statistic**`                 `**N = 1,000**`
  <chr>                           <chr>          
1 Average Total Wait Time (hours) 0.13 ± 0.05


 slow-fast-medium Summary Table:
# A tibble: 1 × 2
  `**Statistic**`                 `**N = 1,000**`
  <chr>                           <chr>          
1 Average Total Wait Time (hours) 0.13 ± 0.05


 slow-medium-fast Summary Table:
# A tibble: 1 × 2
  `**Statistic**`                 `**N = 1,000**`
  <chr>                           <chr>          
1 Average Total Wait Time (hours) 0.14 ± 0.05

Results

# Extract the average wait times for each permutation and combine them into a single data frame
oneline_summary_df <- tibble(
  `fast-medium-slow` = oneline_permutation_results[["fast-medium-slow"]]$results_df$Avg_Total_Wait_Time,
  `fast-slow-medium` = oneline_permutation_results[["fast-slow-medium"]]$results_df$Avg_Total_Wait_Time,
  `medium-fast-slow` = oneline_permutation_results[["medium-fast-slow"]]$results_df$Avg_Total_Wait_Time,
  `medium-slow-fast` = oneline_permutation_results[["medium-slow-fast"]]$results_df$Avg_Total_Wait_Time,
  `slow-fast-medium` = oneline_permutation_results[["slow-fast-medium"]]$results_df$Avg_Total_Wait_Time,
  `slow-medium-fast` = oneline_permutation_results[["slow-medium-fast"]]$results_df$Avg_Total_Wait_Time
)

# Create a summary table using gtsummary
oneline_summary_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",   # Display mean and standard deviation
    digits = all_continuous() ~ 2,  # Round to 2 decimal places for clarity
  ) %>%
  modify_header(label = "**Permutation**") %>%
  modify_caption("**Summary of Average Wait Times Across Permutations (One Line)**")

**Summary of Average Wait Times Across Permutations (One Line)**
Permutation	N = 1,000¹
fast-medium-slow	0.12 ± 0.05
fast-slow-medium	0.13 ± 0.05
medium-fast-slow	0.12 ± 0.05
medium-slow-fast	0.13 ± 0.05
slow-fast-medium	0.13 ± 0.05
slow-medium-fast	0.14 ± 0.05
¹ Mean ± SD

Type 3: Each teller has their own line (grocery store)

In this scenario, each teller has their own line, and customers choose the shortest line to join. This is similar to the setup in a grocery store, where each checkout lane has its own queue.

The Teller 1 will have \(\mu_{1{\text{slower}}} = 3.5 \text{ customers per hour}\) and Teller 2 will have \(\mu_{2{\text{slower}}} = 2.5 \text{ customers per hour}\) and Teller 3 will have \(\mu_{3{\text{slower}}} = 1.5 \text{ customers per hour}\).

# Function to simulate service with each teller having their own line for one day (3 tellers)
simulate_separate_lines_day <- function(arrival_times, mu1, mu2, mu3) {
  
  # Initialize variables for each teller
  wait_times_teller1 <- c()
  wait_times_teller2 <- c()
  wait_times_teller3 <- c()
  service_times_teller1 <- c()
  service_times_teller2 <- c()
  service_times_teller3 <- c()
  
  # Initialize tellers' availability times and queue lengths
  teller1_free_at <- 0
  teller2_free_at <- 0
  teller3_free_at <- 0
  queue_length_teller1 <- 0
  queue_length_teller2 <- 0
  queue_length_teller3 <- 0
  
  # Initialize an empty data frame to store results
  results <- data.frame(
    ArrivalTime = numeric(0), 
    WaitTime = numeric(0), 
    ServiceTime = numeric(0), 
    Teller = character(0),
    QueueLength_Teller1 = numeric(0),
    QueueLength_Teller2 = numeric(0),
    QueueLength_Teller3 = numeric(0)
  )
  
  # Process each customer arrival
  for (arrival_time in arrival_times) {
    
    # Record the current queue lengths
    current_queue_length1 <- queue_length_teller1
    current_queue_length2 <- queue_length_teller2
    current_queue_length3 <- queue_length_teller3
    
    # Determine which teller the customer will go to based on the shortest queue
    if ((queue_length_teller1 <= queue_length_teller2) && (queue_length_teller1 <= queue_length_teller3)) {
      
      # Assign to Teller 1
      if (arrival_time >= teller1_free_at) {
        # Teller 1 is free, no wait time
        wait_time <- 0
        service_time <- rexp(1, rate = mu1)
        teller1_free_at <- arrival_time + service_time
      } else {
        # Teller 1 is busy, customer waits
        wait_time <- teller1_free_at - arrival_time
        service_time <- rexp(1, rate = mu1)
        teller1_free_at <- teller1_free_at + service_time
        queue_length_teller1 <- queue_length_teller1 + 1
      }
      
      # Record the result for Teller 1
      results <- rbind(results, data.frame(
        ArrivalTime = arrival_time, 
        WaitTime = wait_time, 
        ServiceTime = service_time, 
        Teller = "Teller 1",
        QueueLength_Teller1 = current_queue_length1,
        QueueLength_Teller2 = current_queue_length2,
        QueueLength_Teller3 = current_queue_length3
      ))
      
    } else if ((queue_length_teller2 < queue_length_teller1) && (queue_length_teller2 <= queue_length_teller3)) {
      
      # Assign to Teller 2
      if (arrival_time >= teller2_free_at) {
        # Teller 2 is free, no wait time
        wait_time <- 0
        service_time <- rexp(1, rate = mu2)
        teller2_free_at <- arrival_time + service_time
      } else {
        # Teller 2 is busy, customer waits
        wait_time <- teller2_free_at - arrival_time
        service_time <- rexp(1, rate = mu2)
        teller2_free_at <- teller2_free_at + service_time
        queue_length_teller2 <- queue_length_teller2 + 1
      }
      
      # Record the result for Teller 2
      results <- rbind(results, data.frame(
        ArrivalTime = arrival_time, 
        WaitTime = wait_time, 
        ServiceTime = service_time, 
        Teller = "Teller 2",
        QueueLength_Teller1 = current_queue_length1,
        QueueLength_Teller2 = current_queue_length2,
        QueueLength_Teller3 = current_queue_length3
      ))
      
    } else {
      
      # Assign to Teller 3
      if (arrival_time >= teller3_free_at) {
        # Teller 3 is free, no wait time
        wait_time <- 0
        service_time <- rexp(1, rate = mu3)
        teller3_free_at <- arrival_time + service_time
      } else {
        # Teller 3 is busy, customer waits
        wait_time <- teller3_free_at - arrival_time
        service_time <- rexp(1, rate = mu3)
        teller3_free_at <- teller3_free_at + service_time
        queue_length_teller3 <- queue_length_teller3 + 1
      }
      
      # Record the result for Teller 3
      results <- rbind(results, data.frame(
        ArrivalTime = arrival_time, 
        WaitTime = wait_time, 
        ServiceTime = service_time, 
        Teller = "Teller 3",
        QueueLength_Teller1 = current_queue_length1,
        QueueLength_Teller2 = current_queue_length2,
        QueueLength_Teller3 = current_queue_length3
      ))
    }
    
    # Update queue lengths after the customer is served (if applicable)
    if (arrival_time >= teller1_free_at) queue_length_teller1 <- max(0, queue_length_teller1 - 1)
    if (arrival_time >= teller2_free_at) queue_length_teller2 <- max(0, queue_length_teller2 - 1)
    if (arrival_time >= teller3_free_at) queue_length_teller3 <- max(0, queue_length_teller3 - 1)
  }
  
  # Calculate average wait times for each teller
  avg_wait_time_teller1 <- mean(results$WaitTime[results$Teller == "Teller 1"])
  avg_wait_time_teller2 <- mean(results$WaitTime[results$Teller == "Teller 2"])
  avg_wait_time_teller3 <- mean(results$WaitTime[results$Teller == "Teller 3"])
  
  return(list(
    results = results,
    avg_wait_time_teller1 = avg_wait_time_teller1,
    avg_wait_time_teller2 = avg_wait_time_teller2,
    avg_wait_time_teller3 = avg_wait_time_teller3,
    queue_length_teller1 = queue_length_teller1,
    queue_length_teller2 = queue_length_teller2,
    queue_length_teller3 = queue_length_teller3
  ))
}

run_permutation_simulation_separate_lines <- function(arrival_times_by_day, permutations, 
                                                      rate_labels, num_simulations = 100) {
  
  results_list <- list()
  
  for (i in seq_along(permutations)) {
    rates <- permutations[[i]]
    mu1 <- rates[1]
    mu2 <- rates[2]
    mu3 <- rates[3]
    
    # Run simulations for each day with the current permutation
    separate_sim_results <- replicate(num_simulations, {
      lapply(arrival_times_by_day, function(arrival_times) {
        simulate_separate_lines_day(arrival_times, mu1, mu2, mu3)
      })
    }, simplify = FALSE)
    
    # Extract average wait times and queue lengths for each teller across simulations
    avg_wait_times_teller1 <- sapply(separate_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$results$WaitTime[day_result$results$Teller == "Teller 1"], na.rm = TRUE)))
    })
    avg_wait_times_teller2 <- sapply(separate_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$results$WaitTime[day_result$results$Teller == "Teller 2"], na.rm = TRUE)))
    })
    avg_wait_times_teller3 <- sapply(separate_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$results$WaitTime[day_result$results$Teller == "Teller 3"], na.rm = TRUE)))
    })
    
    avg_total_wait_times <- avg_wait_times_teller1 + avg_wait_times_teller2 + avg_wait_times_teller3
    
    avg_queue_lengths_teller1 <- sapply(separate_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$results$QueueLength_Teller1, na.rm = TRUE)))
    })
    avg_queue_lengths_teller2 <- sapply(separate_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$results$QueueLength_Teller2, na.rm = TRUE)))
    })
    avg_queue_lengths_teller3 <- sapply(separate_sim_results, function(simulation) {
      mean(sapply(simulation, function(day_result) mean(day_result$results$QueueLength_Teller3, na.rm = TRUE)))
    })
    
    # Create a data frame for the summary
    results_df <- data.frame(
      Avg_Total_Wait_Time = avg_total_wait_times / 60,      # Convert to hours
      Avg_Wait_Time_Teller1 = avg_wait_times_teller1 / 60,  # Convert to hours
      Avg_Wait_Time_Teller2 = avg_wait_times_teller2 / 60,  # Convert to hours
      Avg_Wait_Time_Teller3 = avg_wait_times_teller3 / 60,  # Convert to hours
      Avg_Queue_Length_Teller1 = avg_queue_lengths_teller1,
      Avg_Queue_Length_Teller2 = avg_queue_lengths_teller2,
      Avg_Queue_Length_Teller3 = avg_queue_lengths_teller3
    )
    
    results_df <- na.omit(results_df)
    
    # Collect data for step plot for the first simulation run of the 7 days
    first_run_data <- do.call(rbind, lapply(1:7, function(day) {
      day_data <- separate_sim_results[[1]][[day]]$results
      day_data$Day <- paste("Day", day)
      day_data
    }))
    
    # Melt data into long format for plotting
    first_run_data_melted <- melt(first_run_data, id.vars = c("Day", "ArrivalTime", "Teller"), 
                                  measure.vars = "WaitTime", variable.name = "WaitTime_Type", value.name = "WaitTime")
    
    # Custom color palette
    custom_colors <- c("Teller 1" = "#1f78b4", 
                       "Teller 2" = "#33a02c", 
                       "Teller 3" = "#ff7f00")
    
    # Generate step plot
    step_plot <- ggplot(first_run_data_melted, aes(x = ArrivalTime, y = WaitTime, color = Teller)) +
      geom_step(size = 0.8) +
      facet_wrap(~ Day, ncol = 3) +
      labs(x = "Arrival Time (minutes)", y = "Wait Time (minutes)", 
           color = "Teller") +
      scale_color_manual(values = custom_colors) +
      theme_bw(base_size = 12) +
      theme(
        legend.position = "top", 
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 9),
        strip.text = element_text(size = 10),  # Bold facet labels
        plot.title = element_text(size = 14, hjust = 0.5),  # Centered and bold title
        axis.text = element_text(size = 10),  # Adjust axis text size
        axis.title = element_text(size = 12),  # Bold axis titles
        axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis labels
        panel.grid.major = element_blank(),  # Remove major grid lines
        panel.grid.minor = element_blank()
      ) +
      ggtitle(paste("Customer Wait Times for", rate_labels[i], "Over 7 Days"))
    
    # Store the results
    results_list[[rate_labels[i]]] <- 
      list(
        results_df = results_df,
        step_plot = step_plot
      )
  }
  
  return(results_list)
}

# Run the function for all permutations
separate_line_permutation_results <- run_permutation_simulation_separate_lines(arrival_times_by_day, 
                                                                               service_rate_permutations_half,
                                                                               rate_labels, num_simulations = 10)

# Loop to access and display results
for (perm_name in names(separate_line_permutation_results)) {
  cat("\n", perm_name, "Summary Table:\n")
  
  # Apply gtsummary to the results_df for a detailed summary table
  summary_table <- separate_line_permutation_results[[perm_name]]$results_df %>%
    tbl_summary(
      statistic = all_continuous() ~ "{mean} ± {sd}",
      label = list(
        Avg_Total_Wait_Time ~ "Average Total Wait Time (hours)",
        Avg_Wait_Time_Teller1 ~ "Average Wait Time for Teller 1 (hours)",
        Avg_Wait_Time_Teller2 ~ "Average Wait Time for Teller 2 (hours)",
        Avg_Wait_Time_Teller3 ~ "Average Wait Time for Teller 3 (hours)",
        Avg_Queue_Length_Teller1 ~ "Average Queue Length for Teller 1",
        Avg_Queue_Length_Teller2 ~ "Average Queue Length for Teller 2",
        Avg_Queue_Length_Teller3 ~ "Average Queue Length for Teller 3"
      ),
      digits = all_continuous() ~ 2
    ) %>%
    modify_header(label = "**Statistic**") %>%
    modify_caption(paste("**Separate Lines Service Wait Times for", perm_name, "**"))
  
  # Extract the gtsummary table as a data frame
  summary_df <- as_tibble(summary_table)
  
  # Print the data frame in a clean format
  print(summary_df)
  print(separate_line_permutation_results[[perm_name]]$step_plot)
}


 fast-medium-slow Summary Table:
# A tibble: 7 × 2
  `**Statistic**`                        `**N = 10**`
  <chr>                                  <chr>       
1 Average Total Wait Time (hours)        0.94 ± 0.16 
2 Average Wait Time for Teller 1 (hours) 0.14 ± 0.03 
3 Average Wait Time for Teller 2 (hours) 0.26 ± 0.09 
4 Average Wait Time for Teller 3 (hours) 0.54 ± 0.15 
5 Average Queue Length for Teller 1      1.49 ± 0.27 
6 Average Queue Length for Teller 2      1.26 ± 0.26 
7 Average Queue Length for Teller 3      1.10 ± 0.28


 fast-slow-medium Summary Table:
# A tibble: 7 × 2
  `**Statistic**`                        `**N = 10**`
  <chr>                                  <chr>       
1 Average Total Wait Time (hours)        0.93 ± 0.35 
2 Average Wait Time for Teller 1 (hours) 0.14 ± 0.05 
3 Average Wait Time for Teller 2 (hours) 0.55 ± 0.25 
4 Average Wait Time for Teller 3 (hours) 0.24 ± 0.09 
5 Average Queue Length for Teller 1      1.40 ± 0.26 
6 Average Queue Length for Teller 2      1.34 ± 0.30 
7 Average Queue Length for Teller 3      0.94 ± 0.27


 medium-fast-slow Summary Table:
# A tibble: 7 × 2
  `**Statistic**`                        `**N = 10**`
  <chr>                                  <chr>       
1 Average Total Wait Time (hours)        0.88 ± 0.20 
2 Average Wait Time for Teller 1 (hours) 0.26 ± 0.06 
3 Average Wait Time for Teller 2 (hours) 0.14 ± 0.05 
4 Average Wait Time for Teller 3 (hours) 0.48 ± 0.16 
5 Average Queue Length for Teller 1      1.64 ± 0.32 
6 Average Queue Length for Teller 2      1.24 ± 0.32 
7 Average Queue Length for Teller 3      1.13 ± 0.31


 medium-slow-fast Summary Table:
# A tibble: 7 × 2
  `**Statistic**`                        `**N = 10**`
  <chr>                                  <chr>       
1 Average Total Wait Time (hours)        1.15 ± 0.22 
2 Average Wait Time for Teller 1 (hours) 0.30 ± 0.09 
3 Average Wait Time for Teller 2 (hours) 0.71 ± 0.17 
4 Average Wait Time for Teller 3 (hours) 0.15 ± 0.06 
5 Average Queue Length for Teller 1      1.75 ± 0.23 
6 Average Queue Length for Teller 2      1.58 ± 0.21 
7 Average Queue Length for Teller 3      1.12 ± 0.22


 slow-fast-medium Summary Table:
# A tibble: 7 × 2
  `**Statistic**`                        `**N = 10**`
  <chr>                                  <chr>       
1 Average Total Wait Time (hours)        0.89 ± 0.18 
2 Average Wait Time for Teller 1 (hours) 0.54 ± 0.15 
3 Average Wait Time for Teller 2 (hours) 0.13 ± 0.04 
4 Average Wait Time for Teller 3 (hours) 0.22 ± 0.06 
5 Average Queue Length for Teller 1      1.74 ± 0.28 
6 Average Queue Length for Teller 2      1.28 ± 0.28 
7 Average Queue Length for Teller 3      1.02 ± 0.28


 slow-medium-fast Summary Table:
# A tibble: 7 × 2
  `**Statistic**`                        `**N = 10**`
  <chr>                                  <chr>       
1 Average Total Wait Time (hours)        1.09 ± 0.26 
2 Average Wait Time for Teller 1 (hours) 0.64 ± 0.19 
3 Average Wait Time for Teller 2 (hours) 0.29 ± 0.09 
4 Average Wait Time for Teller 3 (hours) 0.16 ± 0.08 
5 Average Queue Length for Teller 1      1.97 ± 0.27 
6 Average Queue Length for Teller 2      1.58 ± 0.28 
7 Average Queue Length for Teller 3      1.19 ± 0.29

Results

# Extract the average wait times for each permutation and combine them into a single data frame
separate_summary_df <- tibble(
  `fast-medium-slow` = separate_line_permutation_results[["fast-medium-slow"]]$results_df$Avg_Total_Wait_Time,
  `fast-slow-medium` = separate_line_permutation_results[["fast-slow-medium"]]$results_df$Avg_Total_Wait_Time,
  `medium-fast-slow` = separate_line_permutation_results[["medium-fast-slow"]]$results_df$Avg_Total_Wait_Time,
  `medium-slow-fast` = separate_line_permutation_results[["medium-slow-fast"]]$results_df$Avg_Total_Wait_Time,
  `slow-fast-medium` = separate_line_permutation_results[["slow-fast-medium"]]$results_df$Avg_Total_Wait_Time,
  `slow-medium-fast` = separate_line_permutation_results[["slow-medium-fast"]]$results_df$Avg_Total_Wait_Time
)

# Create a summary table using gtsummary
separate_summary_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",   # Display mean and standard deviation
    digits = all_continuous() ~ 2,  # Round to 2 decimal places for clarity
  ) %>%
  modify_header(label = "**Permutation**") %>%
  modify_caption("**Summary of Average Wait Times Across Permutations (Separate Lines)**")

**Summary of Average Wait Times Across Permutations (Separate Lines)**
Permutation	N = 10¹
fast-medium-slow	0.94 ± 0.16
fast-slow-medium	0.93 ± 0.35
medium-fast-slow	0.88 ± 0.20
medium-slow-fast	1.15 ± 0.22
slow-fast-medium	0.89 ± 0.18
slow-medium-fast	1.09 ± 0.26
¹ Mean ± SD

Conclusion

Code

# Extract the average wait times for each permutation and combine them into a single data frame for tandem
tandem_summary_df <- tibble(
  `fast-medium-slow` = tandem_permutation_results$`fast-medium-slow`$results_df$Avg_Total_Wait_Time,
  `fast-slow-medium` = tandem_permutation_results$`fast-slow-medium`$results_df$Avg_Total_Wait_Time,
  `medium-fast-slow` = tandem_permutation_results$`medium-fast-slow`$results_df$Avg_Total_Wait_Time,
  `medium-slow-fast` = tandem_permutation_results$`medium-slow-fast`$results_df$Avg_Total_Wait_Time,
  `slow-fast-medium` = tandem_permutation_results$`slow-fast-medium`$results_df$Avg_Total_Wait_Time,
  `slow-medium-fast` = tandem_permutation_results$`slow-medium-fast`$results_df$Avg_Total_Wait_Time
)

# Create summary table for tandem system
tandem_summary <- tandem_summary_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    digits = all_continuous() ~ 2
  ) %>%
  modify_caption("**Summary of Average Wait Times Across Permutations**")

# Repeat for oneline
oneline_summary_df <- tibble(
  `fast-medium-slow` = oneline_permutation_results[["fast-medium-slow"]]$results_df$Avg_Total_Wait_Time,
  `fast-slow-medium` = oneline_permutation_results[["fast-slow-medium"]]$results_df$Avg_Total_Wait_Time,
  `medium-fast-slow` = oneline_permutation_results[["medium-fast-slow"]]$results_df$Avg_Total_Wait_Time,
  `medium-slow-fast` = oneline_permutation_results[["medium-slow-fast"]]$results_df$Avg_Total_Wait_Time,
  `slow-fast-medium` = oneline_permutation_results[["slow-fast-medium"]]$results_df$Avg_Total_Wait_Time,
  `slow-medium-fast` = oneline_permutation_results[["slow-medium-fast"]]$results_df$Avg_Total_Wait_Time
)

# Create summary table for oneline system
oneline_summary <- oneline_summary_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    digits = all_continuous() ~ 2
  ) %>%
  modify_caption("**Summary of Average Wait Times Across Permutations**")

# Repeat for separate lines
separate_summary_df <- tibble(
  `fast-medium-slow` = separate_line_permutation_results[["fast-medium-slow"]]$results_df$Avg_Total_Wait_Time,
  `fast-slow-medium` = separate_line_permutation_results[["fast-slow-medium"]]$results_df$Avg_Total_Wait_Time,
  `medium-fast-slow` = separate_line_permutation_results[["medium-fast-slow"]]$results_df$Avg_Total_Wait_Time,
  `medium-slow-fast` = separate_line_permutation_results[["medium-slow-fast"]]$results_df$Avg_Total_Wait_Time,
  `slow-fast-medium` = separate_line_permutation_results[["slow-fast-medium"]]$results_df$Avg_Total_Wait_Time,
  `slow-medium-fast` = separate_line_permutation_results[["slow-medium-fast"]]$results_df$Avg_Total_Wait_Time
)

# Create summary table for separate lines system
separate_summary <- separate_summary_df %>%
  tbl_summary(
    statistic = all_continuous() ~ "{mean} ± {sd}",
    digits = all_continuous() ~ 2
  ) %>%
  modify_caption("**Summary of Average Wait Times Across Permutations**")

# Combine all three tables side-by-side
combined_summary <- tbl_merge(
  tbls = list(tandem_summary, oneline_summary, separate_summary),
  tab_spanner = c("**Tandem**", "**One Line**", "**Separate Lines**")
)

# Print the combined summary table
combined_summary

**Summary of Average Wait Times Across Permutations**
Characteristic	Tandem	One Line	Separate Lines
Characteristic	N = 1,000¹	N = 1,000¹	N = 10¹
fast-medium-slow	2.79 ± 0.34	0.12 ± 0.05	0.94 ± 0.16
fast-slow-medium	2.78 ± 0.35	0.13 ± 0.05	0.93 ± 0.35
medium-fast-slow	2.76 ± 0.35	0.12 ± 0.05	0.88 ± 0.20
medium-slow-fast	2.79 ± 0.35	0.13 ± 0.05	1.15 ± 0.22
slow-fast-medium	2.78 ± 0.34	0.13 ± 0.05	0.89 ± 0.18
slow-medium-fast	2.77 ± 0.34	0.14 ± 0.05	1.09 ± 0.26
¹ Mean ± SD

The simulation results show that the average total wait time for customers varies depending on the order of service rates at the tellers. The One Line system generally outperforms the Separate Lines and Tandem system in terms of average total wait time.
The results suggest that the order of service rates at the tellers can have a significant impact on customer wait times. In the fast-medium-slow permutation is the most efficient, while the slow-fast-medium permutation is the least efficient.
In the One Line and Tandem systems, the average total wait time is more consistent across different permutations, while in the Separate Lines system, the order of service rates can lead to differences in average total wait time.