11 随机服务系统模拟
11.1 笔记
这里主要学习一下R包simmer的一些用法，主要参考该包的文档.
11.1.1 基本用法
先实例化一个新的模拟环境
library(simmer)
#>
#> Attaching package: 'simmer'
#> The following object is masked from 'package:GGally':
#>
#> wrap
#> The following object is masked from 'package:dplyr':
#>
#> select
#> The following object is masked from 'package:tidyr':
#>
#> separate
set.seed(42)
env <- simmer("SuperDuperSim")
env
#> simmer environment: SuperDuperSim | now: 0 | next:
#> { Monitor: in memory }建立一个简单的轨迹：模拟患者问诊。什么是轨迹呢，实际上就是对每名对象的服务过程。
## trajectory初始化一个轨迹对象
patient <- trajectory("patients' path") %>%
## 添加一个护士事件
## 抓住一个护士资源
seize("nurse", 1) %>%
## 服务时间分布
timeout(function() rnorm(1, 15)) %>%
## 释放一个护士资源
release("nurse", 1) %>%
## 添加一个医生事件
seize("doctor", 1) %>%
timeout(function() rnorm(1, 20)) %>%
release("doctor", 1) %>%
## 添加一个管理员事件
seize("administration", 1) %>%
timeout(function() rnorm(1, 5)) %>%
release("administration", 1)往模拟系统添加资源：
env %>%
## 1个护士
add_resource("nurse", 1) %>%
## 2个医生
add_resource("doctor", 2) %>%
## 1个管理员
add_resource("administration", 1) %>%
## 患者到达的时间间隔分布
add_generator("patient", patient, function() rnorm(1, 10, 2))
#> simmer environment: SuperDuperSim | now: 0 | next: 0
#> { Monitor: in memory }
#> { Resource: nurse | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Resource: doctor | monitored: TRUE | server status: 0(2) | queue status: 0(Inf) }
#> { Resource: administration | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: patient | monitored: 1 | n_generated: 0 }测试运行我们的模拟系统：
env %>%
## 运行80个时间单位
run(80) %>%
## 验证当前的模拟时间点
now()
#> [1] 80
## 未来三个事件到来的时刻
env %>% peek(3)
#> [1] 80.69540 81.62105 81.62105也可以指定运行多少个事件：
env %>%
stepn() %>% # 1 step
print() %>%
stepn(3) # 3 steps
#> simmer environment: SuperDuperSim | now: 80.6953988949657 | next: 80.6953988949657
#> { Monitor: in memory }
#> { Resource: nurse | monitored: TRUE | server status: 1(1) | queue status: 1(Inf) }
#> { Resource: doctor | monitored: TRUE | server status: 1(2) | queue status: 0(Inf) }
#> { Resource: administration | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: patient | monitored: 1 | n_generated: 7 }
#> simmer environment: SuperDuperSim | now: 81.6210531397386 | next: 81.6210531397386
#> { Monitor: in memory }
#> { Resource: nurse | monitored: TRUE | server status: 1(1) | queue status: 2(Inf) }
#> { Resource: doctor | monitored: TRUE | server status: 1(2) | queue status: 0(Inf) }
#> { Resource: administration | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: patient | monitored: 1 | n_generated: 7 }也可以看看系统的最终状态：
env %>% peek(Inf, verbose=TRUE)
#> time process
#> 1 81.62105 patient
#> 2 86.74154 patient4
#> 3 89.36934 patient3接下来说明一下如何重置系统。一个方法是执行一个比较长的时间单位（比较一下上面的最终状态时间）：
env %>%
run(120) %>%
now()
#> [1] 120当然，也可以利用reset()函数：
env %>%
reset() %>%
run(80) %>%
now()
#> [1] 80可以利用R基本函数多次复制同一个系统，比如复制100次同一个系统：
envs <- lapply(1:100, function(i) {
simmer("SuperDuperSim") %>%
add_resource("nurse", 1) %>%
add_resource("doctor", 2) %>%
add_resource("administration", 1) %>%
add_generator("patient", patient, function() rnorm(1, 10, 2)) %>%
run(80)
})length((envs))
#> [1] 100下面这个方法可能会更好。如果单个复制品运行起来花费比较大，可以直接并行地运行他们。这里利用的是parallel包。
library(parallel)
envs <- mclapply(1:100, function(i) {
simmer("SuperDuperSim") %>%
add_resource("nurse", 1) %>%
add_resource("doctor", 2) %>%
add_resource("administration", 1) %>%
add_generator("patient", patient, function() rnorm(1, 10, 2)) %>%
run(80) %>%
wrap()
})如何获取各种属性：
## 返回到达数量
envs[[1]] %>% get_n_generated("patient")
#> [1] 9
## 获取排队人数
envs[[1]] %>% get_queue_count("doctor")
#> [1] 0
## 获取排队容量
envs[[1]] %>% get_queue_size("doctor")
#> [1] Inf
## 获取监控的资源状态
envs %>%
get_mon_resources() %>%
head()
#> resource time server queue capacity queue_size system limit replication
#> 1 nurse 8.680399 1 0 1 Inf 1 Inf 1
#> 2 nurse 18.065588 1 1 1 Inf 2 Inf 1
#> 3 nurse 22.798471 1 0 1 Inf 1 Inf 1
#> 4 doctor 22.798471 1 0 2 Inf 1 Inf 1
#> 5 nurse 27.874898 1 1 1 Inf 2 Inf 1
#> 6 nurse 36.522502 1 0 1 Inf 1 Inf 1
## 获取监控的到达状态
envs %>%
get_mon_arrivals() %>%
head()
#> name start_time end_time activity_time finished replication
#> 1 patient0 8.680399 47.82654 39.14614 TRUE 1
#> 2 patient1 18.065588 62.27708 39.47860 TRUE 1
#> 3 patient2 27.874898 77.13731 40.61481 TRUE 1
#> 4 patient0 7.639828 47.39533 39.75551 TRUE 2
#> 5 patient1 17.311110 65.52527 41.36600 TRUE 2
#> 6 patient2 30.961440 76.90152 40.23089 TRUE 211.1.2 更进一步的轨迹使用方法
前面谈到了，对于轨迹的活动，常用的有seize,timeout,release，但此外还有很多更高级的方法，全部方法如下:
methods(class="trajectory")
#> [1] [ [[ [[<-
#> [4] [<- activate batch
#> [7] branch clone deactivate
#> [10] get_n_activities handle_unfinished join
#> [13] leave length log_
#> [16] print release release_all
#> [19] release_selected release_selected_all renege_abort
#> [22] renege_if renege_in rep
#> [25] rollback seize seize_selected
#> [28] select send separate
#> [31] set_attribute set_capacity set_capacity_selected
#> [34] set_global set_prioritization set_queue_size
#> [37] set_queue_size_selected set_source set_trajectory
#> [40] stop_if synchronize timeout
#> [43] timeout_from_attribute trap untrap
#> [46] wait
#> see '?methods' for accessing help and source code许多的活动方法都可以接受动态参数，也即接受一个函数（可以很复杂，包含很多交互），特别提一下timeout，正如我们前面使用了函数作为参数timeout(function() rnorm(1, 15))，但若只使用timeout(rnorm(1, 15))，只会生成一个随机数然后是静态参数。
轨迹可以跟模拟环境进行很多焦糊，常见的就是获取模拟环境的很多属性。一个需要注意的地方是，我们最好有这样的习惯：首先实例化模拟环境，而不是创建轨迹或者其他什么的，避免难以察觉的错误。
接下来是一下对轨迹联合以及取子集的操作。
join()函数可以联合拼接轨迹，可以利用[、[[、head()、tail()等运算符或者函数取轨迹的子集。
11.2 习题
习题1
- 将时间单位看作小时。
library(simmer)
set.seed(1024)
mu <- 1.0
lambda <- 0.8
fun_rexp_mu <- function() rexp(1, mu)
fun_rexp_lambda <- function() rexp(1, lambda)
bank <- simmer("bank")
customer <- trajectory("customer") %>%
seize("counter",1) %>%
timeout(fun_rexp_mu) %>%
release("counter")
bank %>%
add_resource("counter",1) %>%
# 从8点到17点，24小时后循环
add_generator("customer",customer,from_to(8,17,fun_rexp_lambda,every = 24))
#> simmer environment: bank | now: 0 | next: 0
#> { Monitor: in memory }
#> { Resource: counter | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: customer | monitored: 1 | n_generated: 0 }
# 2天内的服务情况
bank %>%
run(48) %>%
get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 customer0 8.00000 9.217335 1.2173351 TRUE 1
#> 2 customer1 10.03479 10.534241 0.4994537 TRUE 1
#> 3 customer2 14.79583 15.765146 0.9693123 TRUE 1
#> 4 customer3 16.90828 17.084277 0.1759946 TRUE 1
#> 5 customer4 32.00000 32.243604 0.2436042 TRUE 1
#> 6 customer5 36.49778 36.861577 0.3637970 TRUE 1
#> 7 customer6 36.65448 37.902313 1.0407360 TRUE 1
#> 8 customer7 37.54700 38.997167 1.0948535 TRUE 1
#> 9 customer8 38.64765 41.565584 2.5684179 TRUE 1
#> 10 customer9 39.01735 42.006504 0.4409200 TRUE 1
#> 11 customer10 39.76817 43.491433 1.4849281 TRUE 1- 不妨假设： \[ \lambda(t)= \begin{cases} 2 & 0\le t<10\\ 1.9 & 10\le t<20\\ \cdots \\ 1.1 & t\ge 90\\ \end{cases} \]
library(simmer)
set.seed(1024)
mu <- 1.0
bank <- simmer("bank")
fun_rexp_mu <- function() rexp(1, mu)
fun_rexp_lambda <- function(){
t <- now(bank)
lambda <- ifelse(t<90,2-floor(t/10)*0.1,1.1)
rexp(1,lambda)
}
customer <- trajectory("customer") %>%
seize("counter",1) %>%
timeout(fun_rexp_mu) %>%
release("counter",1)
bank %>%
add_resource("counter",1) %>%
add_generator("customer",customer,fun_rexp_lambda)
#> simmer environment: bank | now: 0 | next: 0
#> { Monitor: in memory }
#> { Resource: counter | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: customer | monitored: 1 | n_generated: 0 }
# 测试
bank %>%
run(10)
#> simmer environment: bank | now: 10 | next: 10.3397909260835
#> { Monitor: in memory }
#> { Resource: counter | monitored: TRUE | server status: 1(1) | queue status: 9(Inf) }
#> { Source: customer | monitored: 1 | n_generated: 19 }
bank %>%
get_queue_count("counter")
#> [1] 9
bank %>%
get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 customer0 0.8139149 4.622752 3.8088375 TRUE 1
#> 2 customer1 1.4225824 4.866357 0.2436042 TRUE 1
#> 3 customer2 1.6723093 4.991716 0.1253592 TRUE 1
#> 4 customer3 2.5172887 5.705730 0.7140139 TRUE 1
#> 5 customer4 3.0019448 6.800583 1.0948535 TRUE 1
#> 6 customer5 3.0998351 7.401239 0.6006554 TRUE 1
#> 7 customer6 3.1878324 8.659349 1.2581106 TRUE 1
#> 8 customer7 4.9869445 9.805886 1.1465369 TRUE 1- 我们用上
leave、renege_in等函数：
# 输入：
# L,P.1: 队伍长度大于L，顾客有一定概率P.1离开
# T,P.2: 等待时间超过T，顾客有一定概率P.2离开
# T1: 系统运行时间
# mu,lambda: 服务系统参数
# 输出：离开的顾客人数
simmer.1.3 <- function(L=3,P.1=0.8,T=11,P.2=0.8,T1=200,mu=1.0,lambda=0.9)
{
library(simmer)
fun_rexp_mu <- function() rexp(1, mu)
fun_rexp_lambda <- function() rexp(1, lambda)
bank <- simmer("bank",log_level = 0)
tra1 <- trajectory() %>%
seize("counter",1) %>%
renege_abort() %>%
log_("I'm being attended",level = 1) %>%
timeout(fun_rexp_mu) %>%
release("counter",1) %>%
log_("finished",level = 1)
customer <- trajectory("customer") %>%
log_("here I am",level = 1) %>%
log_(function(){
paste("the length of queue is",get_queue_count(bank,"counter"))
},level = 1) %>%
# 一定概率离开
log_("I am leaving...",level = 1) %>%
leave(function(){
l <- get_queue_count(bank, "counter") # 队伍长度
ifelse(l>L,P.1,0) # 若队伍长度大于L以P.1概率离开
}) %>%
log_("It is a joke...",level = 1) %>%
# 等待时间一定概率P.2离开
renege_in(
T,
out = trajectory() %>%
log_("lost my patience. Reneging...",level = 1) %>%
leave(P.2) %>%
log_("It is a joke...",level = 1) %>%
join(tra1)
) %>%
join(tra1)
bank %>%
add_resource("counter", 1) %>%
add_generator("customer", customer, fun_rexp_lambda) %>%
run(T1) %>%
invisible %>%
get_mon_arrivals() %>%
dplyr::summarise(leave_count = sum(!finished))
}
simmer.1.3()
#> leave_count
#> 1 14以默认的参数，重复运行多次求均值：
N <- 1000
out <- vector("integer",N)
for (i in 1:N) {
s <- simmer.1.3()
out[[i]] <- s$leave_count
}
mean(out)
#> [1] 20.977习题2
没想明白，继续想想……