随机服务系统模拟—R实现（二）
M/M/1随机服务系统的模拟
M/M/1模型是一种出生-死亡过程，此随机过程中的每一个状态代表模型中人数的数目。因为模型的队列长度无限且参与人数亦无限，故此状态数目亦为无限。例如状态0表示模型闲置、状态1表示模型有一人在接受服务、状态2表示模型有二人（一人正接受服务、一人在等候），如此类推。 此模型中，出生率（即加入队列的速率）λ在各状态中均相同，死亡率（即完成服务离开队列的速率）μ亦在各状态中相同（除了状态0，因其不可能有人离开队列）。故此，在任何状态下，只有两种事情可能发生：有人加入队列。如果模型在状态k，它会以速率λ进入状态k + 1有人离开队列。如果模型在状态k（k不等于0），它会以速率μ进入状态k − 1，由此可见，模型的隐定条件为λ < μ。如果死亡率小于出生率，则队列中的平均人数为无限大，此这种系统没有平衡点。
一、M/M/1随机服务系统的模拟
M/M/1服务系统：（1）队列长度没有限制；（2）顾客到达的时间间隔和服务时间均服从指数分布；（3）服务台数量为1。
1. 系统的理论绩效指标
模型参数符号说明
| 参数 | 平均到达率 | 平均服务率 | 系统服务强度 | 系统空闲概率 | 系统平均顾客数 | 队列平均人数 | 平均逗留时间 | 平均等待时间 |
|---|---|---|---|---|---|---|---|---|
| 符号 | \(\lambda\) | \(\mu\) | \(\rho\) | \(P_0\) | \(L_s\) | \(L_q\) | \(W_s\) | \(W_q\) |
R计算程序
Lambda <- 1
Mue <-2
Rho <-Lambda/Mue
p0=1-Rho
Lq= Rho ^ 2/(1-Rho)
Wq = Lq / Lambda
Ls <- Lq + Rho
Ws <- Ls / Lambda
R计算结果
| 参数 | 平均到达率 | 平均服务率 | 系统服务强度 | 系统空闲概率 | 系统平均顾客数 | 队列平均人数 | 平均逗留时间 | 平均等待时间 |
|---|---|---|---|---|---|---|---|---|
| 符号 | \(\lambda\) | \(\mu\) | \(\rho\) | \(P_0\) | \(L_s\) | \(L_q\) | \(W_s\) | \(W_q\) |
| 理论值 | 1 | 2 | 0.5 | 0.5 | 1 | 0.5 | 1 | 0.5 |
2. 系统的R模拟仿真
R模型构建
library(dplyr)
library(simmer)
T0=10
T1=10000
lambda=1.0
mu=2.0
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("柜台") %>%
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 <- 1/(mu-lambda)
cat('模拟的平均逗留时间Ws=', ER,
' 期望值=', ER.true, '\n')
R计算结果
cat('估计的平均逗留时间Ws=', ER,' 期望值=', ER.true, '\n')
模拟的平均逗留时间Ws=0.9900721 期望值= 1
3. R模拟可视化
mon1=get_mon_arrivals(bank)
head(mon1,6)
name start_time end_time activity_time finished replication
1 顾客0 2.501759 2.505050 0.003290978 TRUE 1
2 顾客1 2.748517 2.942109 0.193591292 TRUE 1
3 顾客2 4.491264 4.903304 0.412040757 TRUE 1
4 顾客3 4.581213 5.283519 0.380215150 TRUE 1
5 顾客4 4.783831 6.223558 0.940038339 TRUE 1
6 顾客5 5.621871 7.052889 0.829331192 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')
二、模拟仿真总结
随机服务系统在我们日常生活、工业生产、科学技术、军事领域中是经常遇到的随机模型，例如研究银行、理发店、商店、海关通道、高速路收费口等服务人员个数的设置和排队规则，研究计算机网络网关、移动网络的调度规则，等等。某些随机服务系统可以进行严格理论分析得到各种问题的理论解， 但是随机服务系统中存在大量随机因素，使得理论分析变得很困难以至于不可能。例如，即使是上面的银行服务问题可能的变化因素就包括：顾客到来用齐次泊松过程还是非齐次泊松过程，柜员有多少个，是否不同时间段柜员个数有变化，柜员服务时间服从什么样的分布，顾客排队按照什么规则，是否VIP顾客提前服务，顾客等候过长时会不会放弃排队，等等。 包含了这么多复杂因素的随机服务系统的理论分析会变得异常复杂，完全靠理论分析无法解决问题，可用随机模拟方法给出答案。
参考文献
1.(Simmer 2019带你飞 )[https://www.sohu.com/a/344940911_100040805]
2.(Simmer仿真平台高级使用技巧)[https://segmentfault.com/a/1190000019820794]
