我正在尝试通过在 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

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

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

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

8

  • 1
    也许我做了另一个不该做的假设——你可以使用包simmer,但所有的模拟都为你完成了?我想我只是假设你想创建函数。我的假设错了吗?


    – 

  • @Kat:谢谢你的回复!我听说过 simmer 包,也试过,效果不错!作为学习示例,我正在尝试自己编写一个模拟函数流程!:)


    – 

  • @ Kat:早上好!对此有什么想法吗?我花了一整晚的时间尝试不同的方法来调试这个问题…… 🙁


    – 


  • @ Kat:我想我可能已经找到问题了……也许问题是步长太大?自然界中的模拟是一个连续过程,但我使用离散过程来近似它……也许我们需要将 time_step <- 1 更改为 time_step <- 0.01?


    – 

  • 2
    我 99% 确定问题出在您的到达和离开的分布上。为简单起见 – 一个时间单位 = 1 分钟,速率为每小时 /,这意味着到达的次数约为每小时 8 次。当您调用rpois(1, lambda)一次时,您应该1每 60 次查询中得到 8 次。运行lapply(1:(60 * 24), \(j) rpois(1, lambda)) %>% unlist() %>% sum()/24 您的平均值是 8 吗?(提示:不是。)您需要预期的速率。您需要使用1/λ1/μ我认为您的结果可能正是您想要的…如果错了请告诉我。


    – 


最佳答案
1

您的问题在于计算出发次数

departures <- min(queue, rpois(1, num_servers * mu * time_step))

根据理论,服务时间指数分布,速率为1/mu(而不是您假设的泊松分布)。

您可以通过假设每次出发都具有指数分布(即与服务时间相同)来解决这个问题,但我们可以简化您的模拟,因为我们知道系统中的客户数量是几何分布的,参数为1 − rho

由于您正在查看 M/M/1 队列,因此我将简化您的代码。

# single simulation
run_simulation <- function(sim_time) {
  queue <- initial_queue_size
  queue_length <- numeric(sim_time)

  for (t in 1:sim_time) {
    rho <- lambda/mu
    queue_length[t] <- rgeom(1, 1 - rho)
  }
  
  return(queue_length)
}

set.seed(123)
nsims <- 100
results <- run_simulation(nsims)
mean(results)
# [1] 3.78

这更接近理论值 4(rho /(1 – rho))。

服务器被占用时间的平均比例为:

sum(results>0)/nsims
# [1] 0.83

这与 rho (0.8) 非常接近。

2

  • 非常感谢!是否可以将其扩展到多台服务器?


    – 

  • 应该可以。但我的回答的哪方面不能满足你的问题(模拟 M/M/1 队列系统)?


    –