R Simulation not matching theoretical answer

[swaheera]

I am trying to a M/M/K queue (M/M/c queue - Wikipedia) simulation in R with an arrival rate of 8, service rate of 10 and 1 server. The average queue length at steady state according to the theoretical formula should be rho/(1-rho) where rho = (lambda/mu). In my case, this should result in average queue length of 4.

I tried to do this with an R simulation.

First, I defined the queue parameters:

set.seed(123)
library(ggplot2)
library(tidyr)
library(dplyr)
library(gridExtra)

#  simulation parameters
lambda <- 8          # Arrival rate
mu <- 10               # Service rate
sim_time <- 200       # Simulation time
k_minutes <- 15       # Threshold for waiting time
num_simulations <- 100  # Number of simulations to run
initial_queue_size <- 0  # Initial queue size
time_step <- 1        # Time step for discretization

servers <- c(1)  

Next, I defined a function perform a single simulation. My approach takes the current queue length, adds random arrivals and subtracts random departures - and then repeats this process:

# single simulation
run_simulation <- function(num_servers) {
    queue <- initial_queue_size
    processed <- 0
    waiting_times <- numeric(0)
    queue_length <- numeric(sim_time)
    processed_over_time <- numeric(sim_time)
    long_wait_percent <- numeric(sim_time)
    
    for (t in 1:sim_time) {
        # Process arrivals
        arrivals <- rpois(1, lambda * time_step)
        queue <- queue + arrivals
        
        # Process departures
        departures <- min(queue, rpois(1, num_servers * mu * time_step))
        queue <- queue - departures
        processed <- processed + departures
        
        # Update waiting times
        if (length(waiting_times) > 0) {
            waiting_times <- waiting_times + time_step
        }
        if (arrivals > 0) {
            waiting_times <- c(waiting_times, rep(0, arrivals))
        }
        if (departures > 0) {
            waiting_times <- waiting_times[-(1:departures)]
        }
        
        # Record metrics
        queue_length[t] <- queue
        processed_over_time[t] <- processed
        long_wait_percent[t] <- ifelse(length(waiting_times) > 0,
                                       sum(waiting_times > k_minutes) / length(waiting_times) * 100,
                                       0)
    }
    
    return(list(queue_length = queue_length, 
                processed_over_time = processed_over_time, 
                long_wait_percent = long_wait_percent))
}

I then run this simulation:

results <- lapply(servers, function(s) {
    replicate(num_simulations, run_simulation(s), simplify = FALSE)
})

And finally, I tidy everything up into data frames:

# Function to create data frames for plotting
create_plot_data <- function(results, num_servers) {
    plot_data_queue <- data.frame(
        Time = rep(1:sim_time, num_simulations),
        QueueLength = unlist(lapply(results, function(x) x$queue_length)),
        Simulation = rep(1:num_simulations, each = sim_time),
        Servers = num_servers
    )
    
    plot_data_processed <- data.frame(
        Time = rep(1:sim_time, num_simulations),
        ProcessedOrders = unlist(lapply(results, function(x) x$processed_over_time)),
        Simulation = rep(1:num_simulations, each = sim_time),
        Servers = num_servers
    )
    
    plot_data_wait <- data.frame(
        Time = rep(1:sim_time, num_simulations),
        LongWaitPercent = unlist(lapply(results, function(x) x$long_wait_percent)),
        Simulation = rep(1:num_simulations, each = sim_time),
        Servers = num_servers
    )
    
    return(list(queue = plot_data_queue, processed = plot_data_processed, wait = plot_data_wait))
}

plot_data <- lapply(seq_along(servers), function(i) {
    create_plot_data(results[[i]], servers[i])
})

plot_data_queue <- do.call(rbind, lapply(plot_data, function(x) x$queue))
plot_data_processed <- do.call(rbind, lapply(plot_data, function(x) x$processed))
plot_data_wait <- do.call(rbind, lapply(plot_data, function(x) x$wait))

My Problem: When I calculate the average queue length, I get the following:

> mean(plot_data_queue$QueueLength)
[1] 2.46215

And this average does not match the theoretical answer.

Can someone please help me understand what is lacking in my approach and what I can to do to fix this?

Thanks!