我正在尝试通过在 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
您的问题在于计算出发次数。
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 队列系统)?
–
|
simmer
,但所有的模拟都为你完成了?我想我只是假设你想创建函数。我的假设错了吗?–
–
–
–
rpois(1, lambda)
一次时,您应该1
每 60 次查询中得到 8 次。运行lapply(1:(60 * 24), \(j) rpois(1, lambda)) %>% unlist() %>% sum()/24
您的平均值是 8 吗?(提示:不是。)您需要预期的速率。您需要使用1/λ
和1/μ
。我认为您的结果可能正是您想要的…如果错了请告诉我。–
|