随机服务系统模拟—R实现（三）
M/M/c随机服务系统的模拟
M/M/1服务系统：（1）队列长度没有限制；（2）顾客到达的时间间隔和服务时间均服从指数分布；（3）服务台数量为c。
一、M/M/c随机服务系统的模拟
在M/M/c排队系统中，服务台为c个。设系统的到达率为λ，每个用户的服务率为μ。当系统的用户数n>c时，用户离开的速率为cμ,（因为只有c个服务员），当n≤c时，用户离开速率为nμ（因为顾客数小于服务员数）。此时的系统状态（既系统中的用户数）转移图如下图所示。
1. 系统的理论绩效指标
模型参数符号说明
| 参数 | 平均到达率 | 平均服务率 | 系统服务强度 | 系统空闲概率 | 系统平均顾客数 | 队列平均人数 | 平均逗留时间 | 平均等待时间 |
|---|---|---|---|---|---|---|---|---|
| 符号 | \(\lambda\) | \(\mu\) | \(\frac {\lambda}{c\times\mu}\) | \(P_0\) | \(L_s\) | \(L_q\) | \(W_s\) | \(W_q\) |
R计算程序
Lambda <-9
Mue <-4
c=3
Rho <- Lambda / Mue
P0inv <- Rho^c /(factorial(c)*(1-(Rho/c)))
for (i in 0:(c-1)) {
P0inv = P0inv +(Rho^i)/ factorial(i)
}
P0 =1/P0inv
Lq=P0*(Rho^c*Rho/c)/(factorial(c)*(1-(Rho/c))^2)
Wq =Lq/Lambda
Ls <- Lq+Rho
Ws <- Ls/Lambda
R计算结果
| 参数 | 平均到达率 | 平均服务率 | 系统服务强度 | 系统空闲概率 | 系统平均顾客数 | 队列平均人数 | 平均逗留时间 | 平均等待时间 |
|---|---|---|---|---|---|---|---|---|
| 符号 | \(\lambda\) | \(\mu\) | \(\frac {\lambda}{c\times\mu}\) | \(P_0\) | \(L_s\) | \(L_q\) | \(W_s\) | \(W_q\) |
| 理论值 | 9 | 4 | 0.75 | 0.0748 | 3.9533 | 1.7033 | 0.4393 | 0.1893 |
2. 系统的R模拟仿真
R模型构建
library(dplyr)
library(simmer)
library(simmer.plot)
T0=10
T1=10000
lambda=9
mu=4
c=3
set.seed(1234)
## 建立模拟环境
bank <- simmer("bank")
## 用trajectory()建立顾客，并指定顾客的一系列活动
## seize()获取柜台服务资源，如果正在忙，就进入排队
## 服务时间用timeout指定，为了生成多个随机服务时间，
## timeout的参数是返回随机服务时间的而函数而不是时间值本身
customer <-
trajectory("顾客") %>%
seize("柜台") %>%
timeout( function() rexp(1, mu)) %>%
release("柜台")
## 用add_resource生成柜台资源
## 用add_generator()生成顾客到来列
bank %>%
add_resource("柜台",capacity=c) %>%
add_generator("顾客", customer, function() {rexp(1, lambda)} )
## 用run()执行模拟到指定结束时刻
bank %>%
run(until=T1)
R计算程序
## 用get_mon_arrivals()获取各个顾客到来的时间、离开时间、活动时间等，结果是数据框
## 用dplyr::mutate()对数据框增加新变量
resd <- bank %>%
get_mon_arrivals() %>%
dplyr::mutate(waiting_time = end_time - start_time - activity_time,
stay_time = end_time - start_time)
stay_times <- resd %>%
dplyr::filter(start_time >= T0, end_time < T1) %>%
dplyr::select(stay_time)
ER <- mean(stay_times[[1]])
ER.true <- 0.4393
cat('模拟的平均逗留时间ER=', ER,
' 期望值=', ER.true, '\n')
R计算结果
cat('模拟的平均逗留时间ER=', ER,' 期望值=', ER.true, '\n')
模拟的平均逗留时间ER= 0.4537 期望值= 0.4393
3. R模拟可视化
mon1=get_mon_arrivals(bank)
head(mon1,6)
name start_time end_time activity_time finished replication
1 顾客0 0.2779732 0.2796187 0.001645489 TRUE 1
2 顾客1 0.3053908 0.4021865 0.096795646 TRUE 1
3 顾客2 0.4990293 0.7050497 0.206020379 TRUE 1
4 顾客3 0.5090237 0.7185338 0.209510080 TRUE 1
5 顾客4 0.5315368 1.0015560 0.470019170 TRUE 1
6 顾客5 0.6160291 1.1197153 0.414665596 TRUE 1
mon = get_mon_resources(bank)
aggregate(cbind(server, queue) ~ resource, mon, mean)
library(ggplot2)
ggplot(mon, aes(x=server, fill=resource)) +
geom_histogram(binwidth = 0.5) +
facet_grid(.~resource, scales = 'free')
二、模拟仿真总结
随机服务系统在我们日常生活、工业生产、科学技术、军事领域中是经常遇到的随机模型，涉及到泊松过程与指数分布在排队模型中的应用，文中详细地从理论上计算了马尔可夫过程的M/M/c模型绩效指标，又通过随机模拟给出了这个系统绩效指标且作了对比，读者可以大致掌握排队论的应用范围与方法。对于如非指数分布，混线，排队网络，队列数量限制等，不再赘述，感兴趣地读者可以翻阅相关参考资料。
参考文献
1.(Simmer 2019带你飞 )[https://www.sohu.com/a/344940911_100040805]
2.(Simmer仿真平台高级使用技巧)[https://segmentfault.com/a/1190000019820794]
