我正在尝试通过在 R 中模拟系统来对其进行建模，其中到达率为 8，服务率为 10，服务器为 1。根据理论公式，稳定状态下的平均队列长度应为 rho/(1-rho)，其中 rho = (lambda/mu)。在我的例子中，这应该导致平均队列长度为 4。

首先，我定义了队列参数：

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)  

接下来，我定义了一个函数来执行一次模拟。我的方法是采用当前队列长度，添加随机到达并减去随机离开 – 然后重复此过程：

# 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))
}

然后我执行多次模拟：

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

最后，我将所有内容整理到数据框中：

# 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))

我的问题：当我计算平均队列长度时，我得到以下结果：

> mean(plot_data_queue$QueueLength)
[1] 2.46215

而这个平均值和理论答案并不一致。

有人可以帮助我了解我的方法缺少什么以及我该怎么做才能解决这个问题吗？

• 注意：这是我的完整代码：