15 随机服务系统模拟

15.1 概述

我们在第10章中讲到, 较为复杂的随机模型往往难以进行彻底的理论分析, 这时常常使用随机模拟方法产生模型的大量数据, 从产生的数据对模型进行统计推断。 随机服务系统就是这样的一种模型, 经常需要利用随机模拟方法进行研究。

随机服务系统在我们日常生活、工业生产、科学技术、军事领域中是经常遇到的随机模型, 比如,研究银行、理发店、商店、海关通道、高速路收费口等服务人员个数的设置和排队规则, 研究计算机网络网关、移动网络的调度规则,等等。

在概率统计理论中排队论用来研究随机服务系统的数学模型, 可以用来设计适当的服务机构数目和排队规则。 如下面的M/M/1排队系统。

例15.1 设某银行仅有一个柜员, 并简单假设银行不休息。 顾客到来间隔的时间服从独立的指数分布Exp(\(\lambda\)) (\(1/\lambda\)为间隔时间的期望值), 如果柜员正在为先前的顾客服务, 新到顾客就排队等待, 柜员为顾客服务的时间服从均值为\(1/\mu\)的指数分布, 设\(u = \lambda / \mu < 1\)。 设\(X_t\)表示\(t\)时刻在银行内的顾客数(包括正在服务的和正在排队的), 则\(X_t\)是一个连续时马氏链。

这是一个生灭过程马氏链, 有理论结果当系统处于稳定状态时 \[\begin{align*} P(X_t = i) = u^i (1-u), \ i=0,1,2,\dots. \end{align*}\] 设随机变量\(N\)服从\(X_t\)的平稳分布。 于是银行中平均顾客数为 \[\begin{align*} E N = \frac{u}{1-u}, \end{align*}\] 平均队列长度\(EQ\)等于\(E N\)减去平均正在服务人数, 正在服务人数\(Y_t\)\[\begin{align*} Y_t = \begin{cases} 1, & \text{as } X_t > 0, \\ 0, & \text{as } X_t = 0 \end{cases} \end{align*}\] 所以\(E Y_t = P(Y_t = 1) = 1 - P(N = 0) = u\), 于是平均队列长度为 \[\begin{align*} E Q = E N - E Y_t = \frac{u}{1-u} - u = \frac{u^2}{1-u}, \end{align*}\] 设顾客平均滞留时间为\(E R\), 由关系式 \[\begin{align*} E N = \lambda \cdot E R \end{align*}\] 可知平均滞留时间为 \[\begin{align*} E R = \frac{u}{\lambda(1-u)} = \frac{1}{\mu - \lambda}, \end{align*}\] 进一步分析还可以知道顾客滞留时间\(R\)服从均值为\(1/(\mu-\lambda)\)的指数分布。

※※※※※

从上面的例子可以发现, 一个随机服务系统的模型应该包括如下三个要素:

  • 输入过程: 比如,银行的顾客到来的规律。
  • 排队规则: 比如,银行有多个柜员时, 顾客是选最短一队, 还是随机选一队, 还是统一排队等候叫号, 顾客等得时间过长后是否会以一定概率放弃排队。
  • 服务机构:有多少个柜员, 服务时间的分布等。

虽然某些随机服务系统可以进行严格理论分析得到各种问题的理论解, 但是,随机服务系统中存在大量随机因素, 使得理论分析变得很困难以至于不可能。 例如,即使是上面的银行服务问题, 可能的变化因素就包括: 顾客到来用齐次泊松过程还是非齐次泊松过程, 柜员有多少个, 是否不同时间段柜员个数有变化, 柜员服务时间服从什么样的分布, 顾客排队按照什么规则, 是否VIP顾客提前服务, 顾客等候过长时会不会放弃排队, 等等。 包含了这么多复杂因素的随机服务系统的理论分析会变得异常复杂, 完全靠理论分析无法解决问题, 这时,可以用随机模拟方法给出答案。

在模拟随机服务系统时,可以按时间顺序记录发生的事件, 如顾客到来、顾客接受服务、顾客结束服务等, 这样的系统的模拟也叫做离散事件模拟(Discrete Event Simulation, DES)。

离散事件模拟算法可以分为三类:

  • 活动模拟,把连续时间离散化为很小的单位, 比如平均几秒发生一个新事件时可以把时间精确到毫秒, 然后时钟每隔一个小时间单位前进一步并查看所有活动并据此更新系统状态。 缺点是速度太慢。
  • 事件模拟。仅在事件发生时更新时钟和待发生的事件集合。 这样的方法不受编程语言功能的限制, 运行速度很快,也比较灵活,可以实现复杂逻辑, 但是需要自己管理的数据结构和逻辑结构比较复杂, 算法编制相对较难。
  • 过程模拟。把将要到来的各种事件看成是不同的过程, 让不同的过程并行地发生, 过程之间可以交换消息, 并在特殊的软件包或编程语言支持下自动更新系统状态。 在计算机实现中需要借助于线程或与线程相似的程序功能。 这类软件包有C++语言软件包C++SIM和Python语言软件包SimPy, R的simmer包,Julia的simjulia包。 优点是系统逻辑的编码很直观,程序模块化, 需要用户自己管理的数据结构和逻辑结构少。 过程模拟是现在更受欢迎的离散事件模拟方式。

事件模拟算法必须考虑的变量包括:

  • 当前时刻\(t\);
  • 随时间而变化的计数变量, 如\(t\)时刻时到来顾客人数、已离开人数;
  • 系统状态,比如是否有顾客正在接受服务、排队人数、队列中顾客序号。

这些变量仅在有事件发生时(如顾客到来、顾客离开)才需要记录并更新。 为了持续模拟后续事件, 需要维护一个将要发生的事件的列表(下一个到来时刻、下一个离开时刻), 列表仅需要在事件发生时进行更新。 其它变量可以从这三种变量中推算出来,比如, 顾客\(i\)在时间\(t_1\)时到达并排队, 在时间\(t_2\)时开始接受服务, 在时间\(t_4\)时结束服务离开, 则顾客\(i\)的滞留时间为\(t_4 - t_1\)

例15.2 用事件模拟的方法来模拟例15.1。 目的是估计平均滞留时间\(E R\)

想法是,模拟生成各种事件发生的时间, 模拟很长时间, 丢弃开始的一段时间\(T_0\)后, 用\(T_0\)后到达的顾客的总滞留时间除以\(T_0\) 后到达的顾客人数来估计平均滞留时间。 设\(T_0\)时间后每位顾客滞留时间的模拟值为 \(R_i, i=1,2,\dots,m\), 可以用\(\{ R_i \}\)作为随机变量\(R\)的样本来检验\(R\)的分布是否指数分布。

用事件模拟方法进行离散事件模拟的算法关键在于计算系统状态改变的时间, 即各个事件的发生时间, 这个例子中就是顾客到来、顾客开始接受服务、顾客离开这样三种事件, 由此还可以得到每个顾客排队的时间和服务的时间。

在没有明确算法构思时, 可以从时间0开始在纸上按照模型规定 人为地生成一些事件并人为地找到下一事件。 这样可以找到要更新的数据结构和更新的程序逻辑。

下面的算法保持了一个将要发生的事件的集合, 在每次事件发生时更新时钟, 更新时钟时从事件集合中找到最早发生的事件进行处理, 并生成下一事件到事件集合中, 如此重复直到需要模拟的时间长度。

M/M/1 事件模拟算法
{初始化当前时钟\(t \leftarrow 0\), 柜员忙标志\(B \leftarrow 0\), 当前排队人数\(L \leftarrow 0\),
最新到来顾客序号\(i \leftarrow 0\), 正在服务顾客序号\(j \leftarrow 0\), 已服务顾客数\(n \leftarrow 0\) }
从Exp(\(\lambda\))抽取\(X\), 设置下一顾客来到时间\(A \leftarrow X\)
repeat {
\(\quad\) if(\(B = 0\) or (\(B = 1\) and \(A < E\))) { # \(E\)是正在服务的顾客结束时刻
\(\qquad\) \(t \leftarrow A\)
\(\quad\) } else {
\(\qquad\) \(t \leftarrow E\)
\(\quad\) }
\(\quad\) if (\(t > T_1\)) break  # \(T_1\)是预先确定的模拟时长
\(\quad\)
\(\quad\) if(\(t == A\)) { # 待处理到达事件
\(\qquad\) \(L \leftarrow L+1\)
\(\qquad\) \(i \leftarrow i+1\), 记录第\(i\)位顾客到来时间\(a_i \leftarrow t\)
\(\qquad\) 从Exp(\(\lambda\))抽取\(X\), \(A \leftarrow t + X\)
\(\qquad\) if(\(B==0\) ) { # 不用排队,直接服务
\(\qquad\quad\) \(B \leftarrow 1\), \(L \leftarrow L-1\)
\(\qquad\quad\) \(j \leftarrow j+1\), 置第\(j\)位顾客开始服务时间\(s_j \leftarrow t\)
\(\qquad\quad\) 从Exp(\(\mu\))抽取\(Y\), 置\(E \leftarrow t + Y\)
\(\qquad\) }
\(\quad\) } else { # 待处理结束服务事件
\(\qquad\) \(B \leftarrow 0\)
\(\qquad\) \(n \leftarrow n+1\), 记录第\(n\)个顾客结束服务时间\(e_n \leftarrow t\)
\(\qquad\) if(\(L>0\)) { # 排队顾客开始服务
\(\qquad\quad\) \(L \leftarrow L-1\)
\(\qquad\quad\) \(B \leftarrow 1\)
\(\qquad\quad\) \(j \leftarrow j+1\), \(s_j \leftarrow t\)
\(\qquad\quad\) 从Exp(\(\mu\))抽取\(Y\), 置\(E \leftarrow t + Y\)
\(\qquad\) }
\(\quad\) }
}
\(I = \{ i:\ T_0 \leq s_i \leq T_1 \}\), 求\(\{ e_i - a_i, \ i \in I \}\)的平均值作为\(ER\)估计

从这个算法可以看出, 事件模拟方法需要用户自己管理待处理事件集合与时钟, 算法设计难度较大。

※※※※※

用随机服务系统进行建模和模拟研究的一般步骤如下:

  • 提出问题。比如,银行中顾客平均滞留时间与相关参数的关系。
  • 建立模型。比如,设顾客到来服从泊松过程, 设每个顾客服务时间服从独立的指数分布。
  • 数据收集与处理。 比如,估计银行顾客到来速率, 每个顾客平均服务时间,等等。
  • 建立模拟程序。一般使用专用的模拟软件包或专用模拟语言编程。
  • 模拟模型的正确性确认。
  • 模拟试验。
  • 模拟结果分析。

离散事件模拟问题一般都比较复杂, 即使借助于专用模拟软件或软件包, 也很难确保实现的模拟算法与问题的实际设定是一致的。 为此,需要遵循一些提高算法可信度的一般规则。 算法程序一定要仔细检查, 避免出现参数错误、逻辑错误。 在开始阶段, 可以利用较详尽的输出跟踪模拟运行一段时间, 人工检查系统运行符合原始设定。 尽可能利用模块化程序设计, 比如, 在M/M/1问题模拟中, 顾客到来可能遵循不同的规律, 比如时齐泊松过程、非时齐泊松过程, 把产生顾客到来时刻的程序片段模块化并单独检查验证, 就可以避免在这部分出错。 问题实际设定可能比较复杂, 在程序模块化以后, 如果有一种较简单的设定可以得到理论结果, 就可以用理论结果验证算法输出, 保证程序框架正确后, 再利用模块化设计修改各个模块适应实际复杂设定。

例15.3 15.2可以用R语言直接按照时间模拟法实现。

这是单服务员、泊松到来、指数服务时间、先进先出服务系统。 设T0为开始计入统计的时间, T1为结束模拟的时间, lambda为顾客到来速率, mu服务结束速率。

demo.mm1 <- function(T0=0, T1=10000, lambda=1, mu=1.2, DEBUG=FALSE){
  ## 初始化
  t <- 0     # 时钟
  B <- FALSE # 柜员忙标志
  L <- 0     # 排队等候人数
  i <- 0     # 当前到达顾客号
  j <- 0     # 正在接受服务顾客号
  n <- 0     # 已结束服务顾客数
  a <- numeric(0) # 顾客到来时间的集合
  s <- numeric(0) # 顾客开始服务时间的集合
  e <- numeric(0) # 顾客结束服务时间的集合
  X <- rexp(1, lambda); A <- X # 下一顾客到来时间
  E <- NA         # 下一个服务结束时间
  
  ## 时钟反复跳到下一事件发生的时间,并处理事件,更新待发生事件集合
  repeat{
    if(DEBUG){
      cat('\n')
      cat('t=', round(t,2), ' B=', as.integer(B), ' L=', L,
          ' i=', i, ' A=', round(A,2),
          ' n=', n, ' E=', round(E,2), 
          ' j=', j, 
          '\n', sep='')
      cat('到达序列', round(a,2), '\n')
      cat('开始序列', round(s,2), '\n')
      cat('结束序列', round(e,2), '\n')
    }
    if(!B || (B && A<E)){ # 不忙,或待处理下一到来事件
      t <- A
    } else {
      t <- E
    }
    if(t > T1) break
    
    if(t == A){ # 待处理到达事件
      L <- L+1
      i <- i+1; a <- c(a, t)
      X <- rexp(1, lambda); A <- t + X
      if(!B){
        B <- TRUE
        L <- L-1
        j <- j+1; s <- c(s, t)
        Y <- rexp(1, mu); E <- t + Y
      } # if !B
    } else { # 待处理结束服务事件
      B <- FALSE
      n <- n+1; e <- c(e, t)
      if(L > 0){
        L <- L-1
        B <- 1
        j <- j+1; s <- c(s, t)
        Y <- rexp(1, mu); E <- t + Y
      }
    } # end of 待处理结束服务事件
  } # end of repeat
  
  nn <- min(length(a), length(s), length(e))
  a <- a[seq(nn)]; s <- s[seq(nn)]; e <- e[seq(nn)]
  I <- a >= T0 & a <= T1
  ER <- mean(e[I] - a[I])
  ER.true <- 1/(mu-lambda)
  cat('估计的平均滞留时间ER=', ER,
      ' 期望值=', ER.true, '\n')
  c(ER.hat=ER, ER.true=ER.true)
}

测试运行:

※※※※※

15.2 用R的simmer包进行随机服务系统模拟

R扩展包simmer是仿照python语言的simpy包编写的一个用于离散事件模拟的软件包, 属于过程模拟算法, 利用这样的软件包可以使得离散事件模拟程序变得比较简单。

例15.4 用R的simmer包模拟M/M/1随机服务系统。 比较简单的版本,进行短时间的模拟, 用详细过程输出验证算法的合理性。

程序:

测试运行:

set.seed(1)
demo.mm1.simmer1()
## 0.755182: 顾客0: 进入银行
## 0.755182: 顾客0: 开始接受服务
## 0.876604: 顾客0: 结束服务
## 1.93682: 顾客1: 进入银行
## 1.93682: 顾客1: 开始接受服务
## 2.07662: 顾客2: 进入银行
## 2.30022: 顾客1: 结束服务
## 2.30022: 顾客2: 开始接受服务
## 3.32485: 顾客2: 结束服务
## 4.97159: 顾客3: 进入银行
## 4.97159: 顾客3: 开始接受服务
## 5.51127: 顾客4: 进入银行
## 5.65832: 顾客5: 进入银行
## 5.76873: 顾客3: 结束服务
## 5.76873: 顾客4: 开始接受服务
## 6.40375: 顾客4: 结束服务
## 6.40375: 顾客5: 开始接受服务
## 7.04905: 顾客6: 进入银行
## 7.43509: 顾客5: 结束服务
## 7.43509: 顾客6: 开始接受服务
## 8.31388: 顾客6: 结束服务
## 11.473: 顾客7: 进入银行
## 11.473: 顾客7: 开始接受服务
## 12.5082: 顾客8: 进入银行
## 13.0363: 顾客7: 结束服务
## 13.0363: 顾客8: 开始接受服务
## 13.163: 顾客9: 进入银行
## 13.3171: 顾客8: 结束服务
## 13.3171: 顾客9: 开始接受服务
## 13.7515: 顾客10: 进入银行
## 14.3933: 顾客11: 进入银行
## 14.6875: 顾客12: 进入银行
## 15.2533: 顾客13: 进入银行
## 15.2876: 顾客9: 结束服务
## 15.2876: 顾客10: 开始接受服务
## 15.3371: 顾客10: 结束服务
## 15.3371: 顾客11: 开始接受服务
## 15.3594: 顾客14: 进入银行
## 15.8193: 顾客11: 结束服务
## 15.8193: 顾客12: 开始接受服务
## 16.7971: 顾客12: 结束服务
## 16.7971: 顾客13: 开始接受服务
## 17.6278: 顾客13: 结束服务
## 17.6278: 顾客14: 开始接受服务
## 18.8239: 顾客14: 结束服务
## 19.3183: 顾客15: 进入银行
## 19.3183: 顾客15: 开始接受服务
## 19.3556: 顾客16: 进入银行
## 19.5883: 顾客15: 结束服务
## 19.5883: 顾客16: 开始接受服务
## 19.7579: 顾客16: 结束服务
## simmer environment: M/M/1模拟 | now: 20 | next: 20.6760773248559
## { Resource: 柜员 | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Generator: 顾客 | monitored: 1 | n_generated: 18 }

※※※※※

例15.5 上一例子的改进。 用R的simmer包模拟一个银行,多名顾客以指数间隔时间到来(速率lambda) 有一个服务柜台仅一名柜员,只排一队, 服务时间是速率为mu的指数分布随机数。 程序更容易改成其它模拟问题,还增加了一些额外的输出。

R程序:

测试:

set.seed(1)
demo.mm1.simmer2()
## Warning: Attribute retrieval through function arguments is deprecated.
## Use 'get_attribute' instead.

## Warning: Attribute retrieval through function arguments is deprecated.
## Use 'get_attribute' instead.
## 0.755182: 顾客0: 到达
## 0.755182: 顾客0: 排队时间:  0
## 0.876604: 顾客0: 离开:  0.876604105381506
## 1.93682: 顾客1: 到达
## 1.93682: 顾客1: 排队时间:  0
## 2.07662: 顾客2: 到达
## 2.30022: 顾客1: 离开:  2.3002151337181
## 2.30022: 顾客2: 排队时间:  0.223595259614148
## 3.32485: 顾客2: 离开:  3.32485017823051
## 4.97159: 顾客3: 到达
## 4.97159: 顾客3: 排队时间:  0
## 5.51127: 顾客4: 到达
## 5.65832: 顾客5: 到达
## 5.76873: 顾客3: 离开:  5.76872798965007
## 5.76873: 顾客4: 排队时间:  0.257456738084932
## 6.40375: 顾客4: 离开:  6.40375286920846
## 6.40375: 顾客5: 排队时间:  0.745435627139374
## 7.04905: 顾客6: 到达
## 7.43509: 顾客5: 离开:  7.43508916148749
## 7.43509: 顾客6: 排队时间:  0.386036790607969
## 8.31388: 顾客6: 离开:  8.31387513424638
## 11.473: 顾客7: 到达
## 11.473: 顾客7: 排队时间:  0
## 12.5082: 顾客8: 到达
## 13.0363: 顾客7: 离开:  13.0363492322359
## 13.0363: 顾客8: 排队时间:  0.528118697642798
## 13.163: 顾客9: 到达
## 13.3171: 顾客8: 离开:  13.3171271292232
## 13.3171: 顾客9: 排队时间:  0.154149957416097
## 13.7515: 顾客10: 到达
## 14.3933: 顾客11: 到达
## 14.6875: 顾客12: 到达
## 15.2533: 顾客13: 到达
## 15.2876: 顾客9: 离开:  15.2875565067429
## 15.2876: 顾客10: 排队时间:  1.53609961348276
## 15.3371: 顾客10: 离开:  15.3370891404057
## 15.3371: 顾客11: 排队时间:  0.943739658913893
## 15.3594: 顾客14: 到达
## 15.8193: 顾客11: 离开:  15.8193495265548
## 15.8193: 顾客12: 排队时间:  1.13187965742313
## 16.7971: 顾客12: 离开:  16.7971096147388
## 16.7971: 顾客13: 排队时间:  1.54377422102827
## 17.6278: 顾客13: 离开:  17.627787077442
## 17.6278: 顾客14: 排队时间:  2.2683790604488
## 18.8239: 顾客14: 离开:  18.8238581969195
## 19.3183: 顾客15: 到达
## 19.3183: 顾客15: 排队时间:  0
## 19.3556: 顾客16: 到达
## 19.5883: 顾客15: 离开:  19.5883493298518
## 19.5883: 顾客16: 排队时间:  0.232739934309695
## 19.7579: 顾客16: 离开:  19.7579412882115
##      name start_time   end_time activity_time finished replication
## 1   顾客0  0.7551818  0.8766041    0.12142227     TRUE           1
## 2   顾客1  1.9368246  2.3002151    0.36339052     TRUE           1
## 3   顾客2  2.0766199  3.3248502    1.02463504     TRUE           1
## 4   顾客3  4.9715884  5.7687280    0.79713958     TRUE           1
## 5   顾客4  5.5112713  6.4037529    0.63502488     TRUE           1
## 6   顾客5  5.6583172  7.4350892    1.03133629     TRUE           1
## 7   顾客6  7.0490524  8.3138751    0.87878597     TRUE           1
## 8   顾客7 11.4729866 13.0363492    1.56336264     TRUE           1
## 9   顾客8 12.5082305 13.3171271    0.28077790     TRUE           1
## 10  顾客9 13.1629772 15.2875565    1.97042938     TRUE           1
## 11 顾客10 13.7514569 15.3370891    0.04953263     TRUE           1
## 12 顾客11 14.3933495 15.8193495    0.48226039     TRUE           1
## 13 顾客12 14.6874699 16.7971096    0.97776009     TRUE           1
## 14 顾客13 15.2533354 17.6277871    0.83067746     TRUE           1
## 15 顾客14 15.3594080 18.8238582    1.19607112     TRUE           1
## 16 顾客15 19.3183409 19.5883493    0.27000846     TRUE           1
## 17 顾客16 19.3556094 19.7579413    0.16959196     TRUE           1
##     waiting_time
## 1   2.775558e-17
## 2  -1.665335e-16
## 3   2.235953e-01
## 4   0.000000e+00
## 5   2.574567e-01
## 6   7.454356e-01
## 7   3.860368e-01
## 8  -4.440892e-16
## 9   5.281187e-01
## 10  1.541500e-01
## 11  1.536100e+00
## 12  9.437397e-01
## 13  1.131880e+00
## 14  1.543774e+00
## 15  2.268379e+00
## 16 -1.221245e-15
## 17  2.327399e-01

※※※※※

例15.6 用R的simmer包模拟例15.1。 一个银行,多名顾客以指数间隔时间到来(速率lambda), 有一个服务柜台仅一名柜员,只排一队, 服务时间是速率为mu的指数分布随机数, 比较平均等待时间与理论值。

测试:

※※※※※

15.3 用Julia语言的SimJulia包进行离散事件模拟

待完成。

习题

习题1

选择适当的编程语言实现例15.1的算法。 考虑如下的变化情形:

  1. 银行8:00开门,17:00关门,关门后不再允许顾客进入但已进入的顾客会服务完;
  2. 顾客到来服从非齐次的泊松分布,其速率函数\(\lambda(t)\)为阶梯函数;
  3. 如果顾客到来后发现队太长,有一定概率离开; 如果顾客等待了太长时间,有一定概率离开, 这时需要求这些离开顾客的平均人数。

习题2

设计模拟如下离散事件的算法。 设某商场每天开放\(L_0\)时间, 有扒手在该商城出没, 扒手的出现服从强度\(\lambda_1\)的齐次泊松过程, 出现后作案\(X\)时间后离开, 设\(X\)服从对数正态分布, \(\ln X \sim \text{N}(\mu_2, \sigma_2^2)\)。 设有一个警察每隔\(L_3\)时间在商城巡逻\(Y\)时间, \(Y\)服从Exp\((\lambda_2)\)分布, 只要扒手和警察同时出现在商城内扒手就会被抓获。 模拟估计一天时间内有扒手被抓获的概率。

习题3

设计模拟M/M/c随机服务系统的算法。 设服务机构共有\(c\)个服务窗口, 顾客到来服从齐次泊松过程,速率为\(\lambda\), 顾客排成一队, 一旦某个窗口空闲则队头的顾客接受服务, 每个窗口的服务时间均服从期望为\(1/\mu\)的指数分布。

  1. 编程模拟上述随机服务系统;
  2. 模拟估计顾客在这个服务机构的平均滞留时间\(ER\);
  3. 在什么条件下这个服务机构的平均滞留时间能够稳定不变?