M/M/K队列R模拟平均队长与理论值不符，求排查修正

问题背景 #

使用R语言模拟M/M/1队列系统：到达率λ=8，服务率μ=10，服务器数量为1。根据M/M/1队列的稳态理论公式，平均队长应为ρ/(1-ρ)（其中ρ=λ/μ=0.8），计算得理论值为4。但模拟后计算所有时间点的队列长度均值为2.46215，与理论值差距较大。

用户的模拟代码如下：

1. 参数定义 #

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)  

2. 单次模拟函数 #

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

3. 执行模拟与数据整理 #

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

问题分析 #

1. 未排除暂态阶段 #

M/M队列从空队列（初始队列长度0）开始模拟，前一段时间系统处于暂态阶段，队列长度还未收敛到稳态值。直接对所有时间点（包括暂态）取均值，会拉低整体结果，因为暂态期的队列长度普遍远低于稳态值。

2. 离散时间步的近似误差 #

当前模拟采用离散时间步（time_step=1），用泊松分布近似到达和离开事件。虽然泊松分布的期望符合速率要求，但离散时间的近似会引入误差，尤其是当时间步长较大时，无法精确模拟连续时间下的指数分布间隔。

3. 离开事件模拟的逻辑瑕疵 #

对于单服务器，当队列中有顾客时，服务器持续处于忙碌状态，服务完成的事件是指数分布。当前用rpois(1, num_servers * mu * time_step)生成离开数，虽然期望正确，但在离散时间步下，可能出现单次时间步内离开多个顾客的情况（单服务器在1个时间步内完成多个服务），这不符合M/M/1队列的单服务台逻辑（一个服务器同一时间只能处理一个顾客）。

修正方案 #

方案1：排除暂态期，延长模拟时间 #

修改点：

• 延长模拟时间至1000，给系统足够时间进入稳态；
• 统计时排除前200个时间步的暂态数据，只计算稳态阶段的均值。

代码修改：

1. 更新模拟参数：

sim_time <- 1000  # 延长模拟时间
warmup_time <- 200  # 暂态期长度

2. 计算平均队长时过滤暂态数据：

# 只保留warmup_time之后的稳态数据
steady_state_data <- plot_data_queue %>% filter(Time > warmup_time)
mean(steady_state_data$QueueLength)

方案2：改用事件驱动模拟（更准确的连续时间模拟） #

事件驱动模拟直接模拟到达和离开事件的发生时间，完全符合M/M队列的连续时间假设，消除离散时间步的近似误差。

修正后的事件驱动模拟函数：

run_event_driven_simulation <- function(num_servers) {
    current_time <- 0
    queue <- initial_queue_size
    processed <- 0
    server_busy_until <- rep(-Inf, num_servers)  # 记录每个服务器的空闲时间
    queue_length_history <- data.frame(time = numeric(), length = numeric())
    
    # 初始化第一个到达事件
    next_arrival <- rexp(1, lambda)
    
    while (current_time < sim_time) {
        # 找到下一个离开事件的时间
        next_departure <- min(server_busy_until[server_busy_until > current_time])
        if (is.infinite(next_departure)) next_departure <- Inf
        
        # 确定下一个事件类型
        if (next_arrival < next_departure) {
            # 处理到达事件
            time_diff <- next_arrival - current_time
            # 记录当前队列长度在这段时间内的持续时间
            if (current_time > 0) {
                queue_length_history <- rbind(queue_length_history, 
                                             data.frame(time = time_diff, length = queue))
            }
            current_time <- next_arrival
            queue <- queue + 1
            # 生成下一个到达事件
            next_arrival <- current_time + rexp(1, lambda)
            
            # 如果有空闲服务器，立即开始服务
            free_server <- which(server_busy_until <= current_time)[1]
            if (!is.na(free_server)) {
                service_time <- rexp(1, mu)
                server_busy_until[free_server] <- current_time + service_time
            }
        } else {
            # 处理离开事件
            time_diff <- next_departure - current_time
            queue_length_history <- rbind(queue_length_history, 
                                         data.frame(time = time_diff, length = queue))
            current_time <- next_departure
            queue <- queue - 1
            processed <- processed + 1
            
            # 如果队列还有顾客，安排下一个服务
            if (queue > 0) {
                free_server <- which(server_busy_until <= current_time)[1]
                service_time <- rexp(1, mu)
                server_busy_until[free_server] <- current_time + service_time
            }
        }
    }
    
    # 计算时间加权的平均队列长度（稳态）
    # 排除前warmup_time的暂态数据
    warmup_idx <- cumsum(queue_length_history$time) >= warmup_time
    steady_state <- queue_length_history[warmup_idx, ]
    avg_queue_length <- sum(steady_state$time * steady_state$length) / sum(steady_state$time)
    
    return(list(avg_queue_length = avg_queue_length, 
                queue_length_history = queue_length_history))
}

使用事件驱动模拟的计算：

# 执行多次模拟
event_results <- replicate(num_simulations, run_event_driven_simulation(1), simplify = FALSE)
# 提取所有模拟的平均队长
avg_lengths <- sapply(event_results, function(x) x$avg_queue_length)
# 计算均值
mean(avg_lengths)

该方法模拟结果会更接近理论值4，因为它完全贴合M/M队列的连续时间假设，且排除了暂态期。

内容的提问来源于stack exchange，提问作者stats_noob