R：如何模拟队列
我正在使用 R 编程语言。
我正在尝试学习如何模拟队列。
例如：
- 假设一名顾客每 poisson(lambda1) 分钟到达一家咖啡馆
- 有一名服务员（即咖啡馆员工）
- 服务器每 poisson(lambda2) 分钟完成一个订单。
我想在一段时间内进行模拟和绘图：
- 已完成订单数
- 线长
- 刚刚得到服务的客户等待的总时间（即从输入到订单完成）
从数学中，我知道泊松到达之间的差异具有指数分布。如果我选择 lambda1 和 lambda2 的某些值 - 我可以在运行循环之前模拟所有到达时间和服务时间。这是因为到达时间和服务时间完全相互独立且相互关联：
library(ggplot2)
# Parameters
lambda1 <- 5 # customer arrival rate
lambda2 <- 7 # order completion rate
time_period <- 1000 # total time period
# Initialize
arrival_times <- cumsum(rexp(time_period, rate = 1/lambda1))
service_times <- rexp(time_period, rate = 1/lambda2)
completion_times <- numeric(time_period)
queue_length <- numeric(time_period)
total_wait_time <- numeric(time_period)
然后我可以尝试在 1000 个时间点上模拟这个队列：
# Simulation
for (i in 1:time_period) {
if (i == 1) {
completion_times[i] <- arrival_times[i] + service_times[i]
} else {
if (arrival_times[i] < completion_times[i-1]) {
completion_times[i] <- completion_times[i-1] + service_times[i]
} else {
completion_times[i] <- arrival_times[i] + service_times[i]
}
}
queue_length[i] <- sum(arrival_times <= completion_times[i]) - i
total_wait_time[i] <- completion_times[i] - arrival_times[i]
}
df <- data.frame(Time = 1:time_period,
Completed_Orders = 1:time_period,
Queue_Length = queue_length,
Total_Wait_Time = total_wait_time)
然后，我可以绘制结果：
ggplot(df, aes(Time)) +
geom_line(aes(y = Completed_Orders, color = "Completed Orders")) +
geom_line(aes(y = Queue_Length, color = "Queue Length")) +
geom_line(aes(y = Total_Wait_Time, color = "Total Wait Time")) +
scale_y_continuous(sec.axis = sec_axis(~., name = "Total Wait Time")) +
scale_color_manual(values = c("Completed Orders" = "blue", "Queue Length" = "red", "Total Wait Time" = "green")) +
labs(x = "Time", y = "Completed Orders / Queue Length",
title = "Queue Simulation", color = "Legend") +
theme_minimal()
我认为这是正确的 - 这是因为订单数量永远不会减少，只能增加/保持不变。我的图表上的蓝线始终在增加。
我的问题：现在，假设有 2 台服务器。顾客排成一队，然后前往第一个可用的服务器。如果两台服务器同时可用，则客户随机选择。
我不确定如何修改我的代码以“拆分”客户流并将其重新定向到任一服务器。
有人可以告诉我该怎么做吗？
谢谢！
我认为你的模拟不太正确。如果服务时间因客户而异，则已完成的订单不应呈直线上升。我们应该只看到这种具有固定服务时间和始终有顾客的队列的模式。
如果我们有多个服务器，处理计算的一种方法是保留一个与服务器数量相同长度的向量，以跟踪每个服务器下次可用的时间。这将以零向量开始。
当客户到达时，我们会找到下一个可用的服务器（即“下一个可用时间”最短的服务器）。如果一台或多台服务器的“下一个可用时间”小于或等于到达时间，则立即有一台服务器可用，等待时间为0。完成时间只是客户到达时间加上客户服务时间。然后，我们将所选服务器的“下一个可用时间”设置为完成时间。
如果下一个可用服务器时间是在客户到达之后，则等待时间是客户到达和下一个可用服务器时间之间的差值。然后，我们只需添加下一个可用服务器时间和服务时间即可获得完成时间。再次，所选服务器的“下一个可用时间”被设置为完成时间。
对所有客户重复此序列，之后我们就有了到达时间、完成时间和等待时间的向量。如果我们将所有这些向量放入一个数据框中，我们可以通过转为长格式、按时间排序并计算任意时间商店中的顾客数量来计算队列大小。我们需要减去服务员的数量来求出店内任意时刻未得到服务的顾客数量，从而得到队列大小。
这是封装在函数中的该算法的实现：
model_queue <- function(arrival_times, service_times, n_servers) {
if(length(arrival_times) != length(service_times)) {
stop('length of "arrivals" must equal length of "service_time"')
}
n <- length(arrival_times)
complete <- numeric(n)
waits <- numeric(n)
server_free_at <- numeric(n_servers)
for(i in seq(n)) {
server <- which.min(server_free_at)
if(arrival_times[i] > server_free_at[server]) {
waits[i] <- 0
server_free_at[server] <- arrival_times[i] + service_times[i]
complete[i] <- arrival_times[i] + service_times[i]
} else {
waits[i] <- server_free_at[server] - arrival_times[i]
complete[i] <- arrival_times[i] + waits[i] + service_times[i]
server_free_at[server] <- server_free_at[server] + service_times[i]
}
}
data.frame(arrival_times, complete, customer = seq(n), waits) |>
tidyr::pivot_longer(1:2, names_to = "event", values_to = "time") |>
dplyr::arrange(time) |>
transform(completed_orders = cumsum(event == "complete"),
customers = cumsum(event == "arrival_times")) |>
transform(queue = pmax(0, customers - completed_orders - n_servers)) |>
transform(waits = ifelse(event == "complete", NA, waits))
}
让我们看一些可重现的样本数据：
set.seed(1)
lambda1 <- 5
lambda2 <- 7
n <- 100
arrivals <- cumsum(rexp(n, 1/lambda1))
service_time <- rexp(n, 1/lambda2)
让我们使用这些数据并使用单个服务器对队列进行建模，然后检查输出：
queue_df1 <- model_queue(arrivals, service_time, n_servers = 1)
head(queue_df1)
#> customer waits event time completed_orders customers queue
#> 1 1 0.000000 arrival_times 3.775909 0 1 0
#> 2 2 7.390689 arrival_times 9.684123 0 2 1
#> 3 3 10.275599 arrival_times 10.412657 0 3 2
#> 4 4 12.719607 arrival_times 11.111633 0 4 3
#> 5 5 13.196461 arrival_times 13.291976 0 5 4
#> 6 1 NA complete 17.074813 1 5 3
并有两个服务器：
queue_df2 <- model_queue(arrivals, service_time, n_servers = 2)
head(queue_df2)
#> customer waits event time completed_orders customers queue
#> 1 1 0.000000 arrival_times 3.775909 0 1 0
#> 2 2 0.000000 arrival_times 9.684123 0 2 0
#> 3 3 2.884909 arrival_times 10.412657 0 3 1
#> 4 4 5.328917 arrival_times 11.111633 0 4 2
#> 5 5 3.782836 arrival_times 13.291976 0 5 3
#> 6 2 NA complete 13.297566 1 5 2
我们可以将结果绘制如下：
library(geomtextpath)
ggplot(queue_df1, aes(time, queue)) +
geom_step(aes(color = "Queue length")) +
geom_step(aes(y = completed_orders, color = "Completed orders")) +
geom_step(aes(y = waits, color = "Wait time"),
data = subset(queue_df1, !is.na(waits))) +
geom_textvline(aes(xintercept = max(time[event == "arrival_times"]),
label = "Last arrival"), vjust = -0.3, hjust = 0.75) +
scale_color_brewer(NULL, palette = "Set1") +
theme_minimal(base_size = 16)
ggplot(queue_df2, aes(time, queue)) +
geom_step(aes(color = "Queue length")) +
geom_step(aes(y = completed_orders, color = "Completed orders")) +
geom_step(aes(y = waits, color = "Wait time"),
data = subset(queue_df2, !is.na(waits))) +
geom_textvline(aes(xintercept = max(time[event == "arrival_times"]),
label = "Last arrival"), vjust = -0.3, hjust = 0.75) +
scale_color_brewer(NULL, palette = "Set1") +
theme_minimal(base_size = 16)
我们可以看到，使用两台服务器，队列或等待时间要少得多。