# 电子表格中的仿真示例

LLMS 索引： [llms.txt](/llms.txt)

---

本章借助电子表格示例来介绍仿真（simulation）。这些示例相对简单；其中许多既可以在电子表格中完成，也可以手工完成。我们的主要目的，是引入并图解仿真中的若干关键概念；而不是讲授如何构建或开发电子表格模型，而是通过若干简单的电子表格模型，以示例的方式说明任何仿真都具备的一些关键要素。这些关键要素涵盖从模型的组成部分直至对模型进行实验。

尽管某些现实世界的仿真模型——例如风险分析、财务分析、某些可靠性问题等——适合用电子表格来表述和求解，但电子表格对于大多数复杂的现实世界动态事件驱动仿真而言，存在着严重的局限。这类复杂模型更多地是在仿真软件包或通用编程语言中开发的，而不是在电子表格中开发的。

本章的电子表格示例既包含蒙特卡罗仿真（Monte Carlo simulation），也包含随时间跟踪事件的动态仿真。电子表格仿真的实现方式是设计一张仿真表（simulation table），它提供了一种跟踪系统状态随时间变化的方法。尽管这些仿真表包含若干共同的要素，但总体而言它们都是临时拼凑的（ad hoc），也就是说，它们是针对手头问题专门定制的。由于这种临时拼凑的性质，开展动态离散事件仿真需要一个通用的框架。

本章引入的关键要素包括：利用随机数表示不确定性；用于预测排队模型与库存模型中系统性能的描述性统计量；以及把模型理解为由输入、活动、事件、状态、输出和响应所构成所带来的好处。这些主题在此均通过示例加以介绍。

第 1 节讲述基本内容，例如生成随机样本以及仿真模型的基本组成部分；后续示例都依赖于这里学到的概念与技术。其余各节的示例涉及若干应用领域，例如掷硬币、排队或等待线、库存策略、可靠性、命中目标（轰炸或飞镖靶）以及项目网络分析。在这些示例中，我们对模型实现（implementation）只作有限关注；而把重点放在识别模型的组成部分（输入、事件、状态等）以及开展一项简单的实验上，其中的关键概念是试验（trials）或重复（replications）。

本章所有示例都配有 Excel 电子表格解答，可在本书网站 www.bcnn.net 上获取。这些电子表格包含说明与注释，为每个示例补充了技术细节，其中包括实现细节；对于修改现有模型或按照某些习题要求开发新的电子表格仿真模型的学生而言，这些细节可能颇有价值。

## 1 电子表格仿真的基本原理

本节介绍完成本章其余示例仿真所需的基本内容。我们将学习如何在电子表格仿真中生成随机数，以及手工仿真时到哪里去取得随机数。我们将学习如何针对若干简单情形生成随机样本，例如掷硬币、排队仿真中的服务时间与到达时间，以及库存模型中的随机需求。此处学到的方法，足以应付本章示例中用到的任何简单离散统计分布。

### 1.1 如何仿真随机性

本书所考虑的几乎所有仿真模型都包含一个或多个随机变量（random variable）。显然，对于掷硬币而言，需要某种机制来模拟一次随机的投掷，结果为正面或反面。在排队模型中，当到达时间与服务时间事先未知、无法预测但遵循某种统计规律时，服务时间与到达间隔时间（interarrival time）二者都可以由统计分布来定义。

我们如何取得随机量？关于随机数和随机变量生成背后的理论与一般性解释，这里暂且不论。本章采取一种务实、实用的观点：对于电子表格仿真或手工仿真，我们如何取得随机量？

生成任何随机量，都要从生成一个位于 0 与 1 之间区间内的随机数（random number）开始。"随机数"一词始终指在 0 到 1 之间的单位区间内随机生成的量。与之相对，"随机变量"一词指具有特定统计分布的任意随机生成的量。任何生成随机数序列的方法，都称为随机数发生器（Random Number Generator，RNG）。

随机数的生成方法应当产生具有两项重要统计性质的数列：

1. 这些数应当在 0 与 1 之间均匀分布。

2. 随后生成的数应当与此前生成的所有数在统计上相互独立。

均匀性除其他含义外还意味着：落入区间 (a, b)（其中 a 和 b 均介于 0 与 1 之间）的概率，等于该区间的长度 b − a。统计独立性除其他含义外还意味着：知道先前已生成数的取值，对预测随后生成数的取值毫无帮助。此外，它还意味着相邻两个数之间的相关系数为零。

由此得到的这些数通常称为伪随机数（pseudo-random number）。其原因之一在于这样一个悖论：我们希望生成的数在统计上均匀且相互独立，同时我们又希望能加以控制。我们希望在任何时候都能重现该数列，从而使得仿真实验能够得到控制并被复现；还有一个实用原因：这有助于模型的开发与调试。

一般而言，我们如何取得这些随机数？它们从哪里来？（以下文中，VBA 指 Visual Basic for Applications，即 VB 编程语言在 Excel 及其他 Microsoft Office 应用程序中的实现。）随机数有如下若干标准来源：

1. 在 Excel 中，我们可以在单元格的公式里使用内置的工作表函数（worksheet function）RAND。RAND 生成 0 与 1 之间的值。使用 RAND 的公式的典型例子包括生成 0 与 1 之间的随机值，如下所示：

$$
= \mathrm{RAND()}
$$

 或者用于生成两个值之一（0 或 1）的更复杂表达式，如下所示：

$$
= \text{IF} (\text{RAND} () <   = 0. 5, 0, 1)
$$

 注意，工作表函数 IF 在其第一个参数中检查一个条件；若条件为真，则返回第一个值，否则返回第二个值。

2. Excel 还提供工作表函数 RANDBETWEEN()，它返回指定下限与上限之间均匀分布的整数值。Excel 或其某些加载项还提供其他生成随机值的函数。这些我们一概不用。

3. 在 Excel 的工作表单元格中，我们可以调用一个用户编写的 VBA 函数，该函数进而调用 VBA 内部库函数，例如 Rnd()，或本书网站上提供的 VBA 函数（MRG32k3a）。本章的电子表格示例全都带有可从单元格中调用的 VBA 函数 Rnd01()。（Rnd01() 只是简单地调用 VBA 内部库函数 Rnd()。）

4. 仿真软件包提供多种方法，可从一长串标准统计分布以及从数据中，生成随机数和随机变量。

5. 大多数通用编程语言都提供内置的库函数来生成随机数，其中许多还提供从各种统计分布中生成样本的设施。若缺少这些，模型开发者可以自行编写一个或多个算法。

6. 对于手工（纸笔）仿真，可以使用某种物理随机数发生器，例如抛一枚或多枚硬币，或者掷骰子。由于这些做法除最简单的情形外都不切实际，另一种可能是使用随机数表或来自特定分布的随机变量表，例如附录中的表 A.1（均匀分布）与表 A.2（标准正态分布）。

工作表函数 RAND() 和 VB 函数 Rnd() 都不应当用于专业工作；它们都存在已知的缺陷。例如，VB 函数 Rnd() 返回单精度值（约 5 位小数的精度），而所有优秀的 RNG 都返回双精度（约 15 位小数）；其次，Rnd() 的周期很短——在生成相对少量的数之后便会重复。它的周期是 2^24，约 1677 万。RAND() 函数的周期更短，而且在均匀性与独立性这两项关键性质上有更严重的缺陷。（注意 Microsoft 在 Office 2003 中修改了 RAND() 函数；据说它既有错误又有缺陷，因而广受批评。）我们只在例 1 中使用 RAND()，且仅出于教学目的。既然有易于获得的更好的发生器，包括本书网站上 Excel 文件 "RandomNumberTools.xls" 中的 VBA 函数 MRG32k3a，我们深信在专业工作中没有任何理由（或借口）去使用有缺陷的发生器。

大多数 RNG 都可以用一个称为种子（seed）的用户指定数字来初始化。电子表格示例中「重置种子并运行」（'Reset Seed & Run'）按钮背后的 VBA 代码，展示了如何为 VBA 的 Rnd() 发生器设定种子。显然，Excel 的 RAND 发生器无法设定种子（这是不使用它的另一个理由）。MRG32k3a 发生器（在 "RandomNumberTools.xls" 中）附带 VBA 函数 InitializeRNSeed()，用于设定该发生器的种子。通过指定种子，用户可以让随机数发生器重复同一数列。这对调试仿真模型和重复仿真实验都很有价值。请记住，我们所仿真的是随机性；可重复性所带来的这种可控性是一项特性，而不是缺陷。

### 1.2 示例中所用的随机数发生器

每个示例电子表格都带有若干 VBA 函数，用于在若干简单情形下生成随机变量。我们建议在开发新模型时，从一个现有的电子表格模型入手；这样一来，为各个示例编写的全部 VBA 函数都可为新模型所用。

用于在工作表单元格中生成随机值的 VBA 函数与 Excel 工作表函数包括：

1. Rnd01()，用于生成 0 与 1 之间的随机数；

2. DiscreteUniform(min, max)，用于生成从最小值 min 到最大值 max 的随机整数，其中 min 与 max 可以是数字，也可以是包含该数字的单元格区域；

3. DiscreteEmp(rCumProb, rValues)，用于从由两个 Excel 区域（单元格区域）指定的离散分布中生成样本：rCumProb 为累积概率，rValues 为取值；

4. Uniform(low, high)，用于从下限 low 与上限 high 之间的连续均匀分布中生成一个实值（而不只是整数），其中 low 与 high 可以是数字，也可以是包含该数字的单元格；

5. NORMSINV(Rnd01())，用于生成均值为 0、标准差为 1 的标准正态变量。

除 NORMSINV() 之外，上述所有函数都是为本章示例编写的 VBA 函数；NORMSINV() 是 Excel 工作表函数，可以（按上述形式）用来生成标准正态变量。我们建议不要使用工作表函数 RAND() 和 RANDBETWEEN()。所有电子表格解答都可以调用其中任何一个函数，无论所提供的解答是否用到了它；因此，这些示例电子表格是完成某些习题的良好起点。

在本书的这一点上，我们并不期望读者完全理解上述所有函数的工作方式。不过，我们的目标是：让你理解每个函数所生成随机值的性质，并能在修改示例或创建新的电子表格模型时，于习题中加以运用。后续各节在示例中使用这些函数时，会对它们作更详细的说明。

### 1.3 如何使用这些电子表格

电子表格仿真模型位于各解答工作簿的「One Trial」（单次试验）工作表中。该工作表包含模型输入的定义与规格说明、包含用于生成仿真表中每一行从而推进仿真的公式的仿真表，以及若干概括统计量（通常还有一幅直方图），用以表示一次仿真运行的结果。工作表上的两个按钮控制仿真：

1. 每点击一次，「生成新试验」（'Generate New Trial'）按钮就用一列新的随机数完成一次新的仿真运行。

2. 每点击一次，「重置种子并运行」（'Reset Seed & Run'）按钮就把 RNG 种子设为用户在相邻单元格中指定的值，然后用该种子运行仿真。默认种子是任意取值 12,345。

你随时可以点击「重置种子并运行」按钮回到初始结果。

「Experiment」（实验）工作表允许用户开展一项仿真实验：反复执行「One Trial」工作表，每次使用不同的随机数，并把指定的模型响应（取自 Link 单元格）记录到响应表（Response Table）中。用户指定试验次数，例如 100，意思是「One Trial」工作表将执行 100 次，并把 100 个响应值记录到响应表中。这称为重复（replication）；重复的次数就是试验次数。

所有试验的实验结果汇总在「Experiment」工作表上的「多次试验汇总」（'Multi-Trial Summary'）表中。这些汇总统计量包括频数分布及其配套的直方图，以及样本均值、中位数、最小值、最大值等标准统计量。第 2 节的表 6 给出了一个示例。

要在「Experiment」工作表上运行实验，按如下方式使用这两个按钮：

1. 每点击一次，「运行实验」（'Run Experiment'）按钮就按指定的试验次数完成一次实验运行。随后的每次点击，都使用不同的随机数运行该实验。

2. 每点击一次，「重置种子并运行」按钮就把 RNG 种子设为相邻单元格中指定的值（默认 12,345），并完成一次实验运行。每次点击都会复现完全相同的原始结果。

本章所有示例中给出的结果，都来自点击「重置种子并运行」按钮；这保证了读者能够复现出同样的结果。每个电子表格解答都包含关于解答各个方面的注释与说明，其中包括模型逻辑的实现；有些解答还带有一个名为「Explain」（说明）的工作表，对模型的某些方面给出更详尽的解释。

### 1.4 如何仿真掷硬币

如何在电子表格模型中仿真一次掷硬币？本解答展示了如何在电子表格仿真中使用均匀随机数发生器。

### 例 1：掷硬币

我们希望仿真一系列 10 次掷硬币。该硬币是一枚"公平"硬币，也就是说出现正面的概率与出现反面的概率相同。二者各以 0.5 的概率发生。你应当多次运行该仿真，并把结果与现实生活中掷硬币的预期进行比较。这是一个蒙特卡罗仿真的例子；其中并没有跟踪任何事件或时钟时间。

解答见表 1，位于电子表格 "Example2.1CoinToss.xls" 中。事实上，该电子表格中共有四个解答，前两个使用 RAND() 随机数发生器，后两个使用所提供的 VBA 函数 Rnd01()。留作练习（见本章末尾）的是：为这四个解答在统计上的等价性给出一番论证，尽管每次点击「生成新试验」按钮时，这四组结果都不相同。

**表 1 掷硬币的仿真表**

|  | B | C | D | E | F | G | H | I |
|---|---|---|---|---|---|---|---|---|
| 11 |  | 解答 #1A | 解答 #1B |  | 解答 #2A | 解答 #2B |  |  |
| 12 | 投掷 | RAND() | 使用 C 列随机数的解答 | 使用 RAND() 的解答 | 投掷 | Rnd01() | 使用 G 列随机数的解答 | 使用 Rnd01() 的解答 |
| 13 |  |  |  |  |  |  |  |  |
| 14 |  |  |  |  |  |  |  |  |
| 15 | 1 | 0.7317089 | T | T | 1 | 0.9871114 | T | T |
| 16 | 2 | 0.8285837 | T | T | 2 | 0.0225598 | H | H |
| 17 | 3 | 0.0701168 | H | T | 3 | 0.0008356 | H | T |
| 18 | 4 | 0.2147295 | H | H | 4 | 0.2127686 | H | T |
| 19 | 5 | 0.8613179 | T | H | 5 | 0.8586159 | T | H |
| 20 | 6 | 0.4159098 | H | H | 6 | 0.5503362 | T | T |
| 21 | 7 | 0.7492317 | T | H | 7 | 0.5243730 | T | H |
| 22 | 8 | 0.9671661 | T | T | 8 | 0.8884351 | T | H |
| 23 | 9 | 0.4823089 | H | T | 9 | 0.2111118 | H | T |
| 24 | 10 | 0.6738897 | T | H | 10 | 0.7680427 | T | H |

掷硬币的电子表格逻辑很简单。第二个解答（E 列）为：

$$
= \text{IF} (\text{RAND} () <   = 0. 5, \text{"H","T"})
$$

 而第四个解答（I 列）为：

$$
= \text{IF} (\text{Rnd01} () <   = 0. 5, \text{"H","T"})
$$

 每一次调用 RAND() 或 Rnd01() 都会产生一个新的随机数。「One Trial」工作表记录 10 次投掷中正面与反面各自出现的频数，并以直方图把结果绘出。

第一个和第三个解答分作两步，这种做法在所生成的随机数要被使用多次时很有用。C 列与 G 列存放所生成的随机数，D 列与 H 列存放结果 \(( ^ {66} \mathrm{H} ^ {3} \mathrm{o r} ^ {66} \mathrm{T} ^ {3} )\)。第四个解答是例 4 中掷硬币游戏的基础。

在此后所有的电子表格仿真中，我们都使用所提供的 VBA 函数 Rnd01()，或以 Rnd01() 为基础的某个 VBA 函数，而不再使用 RAND()。不过，对于严肃的专业工作，我们建议使用文件 "RandomNumberTools.xls" 中由 VBA 函数 MRG32k3a 实现的发生器。如前所述，Excel 工作表函数 RAND() 与 VBA 函数 Rnd() 都存在严重的局限与缺陷，只应用于随意或教学目的。

### 1.5 如何仿真随机服务时间

下面的示例展示如何从任意离散分布中生成随机样本。这个示例很简单，服务时间只有 3 个可能取值。其关键原理——即从一个随机数变换为具有给定概率的服务时间——很容易推广到任何具有有限个可能取值的离散分布。

此处以及第 1.6 节中就到达时间所描述的方法，被用于本章几乎所有后续示例。

### 例 2：随机服务时间

某自动电话信息服务与每位来电者的通话时间，要么是 3 分钟，要么是 6 分钟，要么是 10 分钟。各服务时长所占的呼叫比例分别为 30%、45% 和 25%。我们希望在电子表格中仿真这些服务时间。我们真正的目的，是学习如何从离散分布中生成随机样本，为随后的排队、库存及其他示例做准备。这同样是一个蒙特卡罗仿真的例子。

表 2 给出了输入规格说明和仿真表的一部分（前 4 位来电者），取自电子表格模型 "Example2.2ServiceTimes.xls"。输入规格说明定义了各服务时间的概率，同时也给出了累积概率。累积概率总是递增到 1.0；正如我们将要看到的，生成方法用到的正是累积概率。该方法总是从一个随机数开始，以所期望的值（本例中即具有指定分布的随机服务时间）结束。电子表格中的仿真表展示出为 25 位来电者所生成服务时间的频率分布（或直方图）。

**表 2 服务时间的输入规格说明与仿真表**

|  | A | B | C | D |
|---|---|---|---|---|
| 4 |  | 服务时间 | 概率 | 累积概率 |
| 5 |  |  |  |  |
| 6 |  |  |  |  |
| 7 |  | 3 | 0.30 | 0.30 |
| 8 |  | 6 | 0.45 | 0.75 |
| 9 |  | 10 | 0.25 | 1.00 |
| 10 |  |  |  |  |
| 11 |  | 来电者人数= | 25 |  |
| 12 |  | 仿真表 |  |  |
| 13 |  | 步骤 | 活动 |  |
| 14 |  | 来电者 | 服务时间 |  |
| 15 |  |  |  |  |
| 16 |  | 1 | 6 |  |
| 17 |  | 2 | 6 |  |
| 18 |  | 3 | 10 |  |
| 19 |  | 4 | 6 |  |

图 1 说明了如何把一个随机数变换为服务时间。可以把它想成：向一块特殊的飞镖靶投掷一支特殊的飞镖；飞镖总会落在单位区间上的某处，且落在 0 到 1 之间直线上任何一点的可能性都相等。飞镖所落入的子区间，决定了所生成的值。

在图 1 中，各概率（位于上方）用单位区间上互不重叠的线段表示，且必须加总为 1。累积概率（紧贴线的下方）用线上的点表示。箭头表示这一变换。该算法的工作方式如下：首先选取一个随机数，记作 R。若 R 落在第一个区间内，也就是 R 小于或等于 0.30，则该过程生成一个 3 分钟的服务时间样本。类似地，若 R 介于 0.30 与 0.75 之间，服务时间为 6 分钟；最后，若 R 大于 0.75，则服务时间为 10 分钟。该算法很容易推广到任何离散分布。

为说明这一过程，我们首先生成 5 个随机数样本，或者在 Excel 中生成（使用 RAND 或 Rnd01 或任何其他随机数发生器），或者从诸如附录中的表 A.1 那样的随机数表中取 5 个样本。假设这 5 个样本是：

0.9871
0.0226
0.0008
0.2128
0.8586

经过图 1 所示的变换后，所得服务时间为：

10
3
3
3
10

在一个小样本（此处样本量为 5）中，我们不能指望各取值的实际出现频数接近其概率；不过，当样本量足够大时，样本频数应当相当接近这些概率。你可以试用电子表格示例 "Example2.2ServiceTimes.xls"，看看小样本与所假定的概率可以相差多远。

![图 1 从随机数到服务时间的"随机飞镖"变换](img/ch02/img_e43a07e186166847.jpg)

**图 1 从随机数到服务时间的"随机飞镖"变换。**

在 Excel 中，我们可能会以为可以实现这一变换——至少对于取值个数较少（此处为 3 个）的分布——用一个嵌套的 IF 工作表函数，例如：

$$
= \text{IF} (\text{Rnd01} () <   = 0. 30, 3, \text{IF} (\text{Rnd01} () <   = 0. 75, 6), 10)\tag{WRONG !!!}
$$

 这个求解尝试引入了一个严重错误。每一次引用 Rnd01() 函数都会产生一个新的随机数；而该变换过程依赖于只使用一个随机数。（你或许愿意把 "Example2.2ServiceTimes.xls" 中的电子表格解答改成这个错误的解答，看看由此得到的概率。）这一做法可以修改为使用两个单元格：第一个存放 Rnd01 的结果，第二个使用修改后的 IF 语句。

更好的办法是使用随所有示例一道提供的 VBA 函数 DiscreteEmp()（离散经验分布，discrete empirical，即由数据定义）。要了解该函数的使用方法，请参看电子表格解答 "Example2.2ServiceTimes.xls"；一个生成随机服务时间的典型单元格如下：

$$
= \text{DiscreteEmp} (D 7: D 9, B 7: B 9)
$$

 其中 \$ D\$7:\$ D\$9 是包含累积概率的单元格区域，\$ B\$7:\$ B\$9 是包含服务时间期望取值的单元格区域。Discrete 意味着只会生成所列出的那些取值。这两个工作表区域中的每一个都必须是列向量（即所有数据在一列中），且长度相同，也就是单元格个数相同。每个地址中行号与列号之前的"$"表示绝对地址而非相对地址；这样一来，在把公式键入「服务时间」列中某一个单元格之后，把它复制粘贴到该列的所有单元格中就更方便了。

此处所使用的、随附所提供的 VBA 函数 DiscreteEmp() 的方法，很容易推广到本章中使用简单离散分布作为模型输入的多数示例（和习题）。要使用 DiscreteEmp()，需在工作表中提供两个区域，一个用于累积概率，另一个用于期望取值。它适用于分布中任意多个取值。

在后续示例中，请记住图 1；它恰好说明了 DiscreteEmp() 的工作原理。虽然它的 VBA 实现涉及一个循环，对非编程人员而言可能显得复杂，但其基本逻辑就是图 1 的简单逻辑——向单位区间这块特殊飞镖靶投掷飞镖这一简单逻辑，以及随之而来的变换。

### 1.6 如何仿真随机到达时间

本节的目标是学习如何生成随机的到达间隔时间（相继到达之间的时间），以及如何计算顾客或其他实体到达系统的时间。到达间隔时间，与例 2 中的服务时间一样，都是由仿真生成的时间长度，称为活动时间（activity time）或活动（activity）。它们与到达时间（arrival time）相区别——后者是仿真时钟上的时刻，表示到达某一设施的时刻。

仿真时钟（simulation clock）是每一个动态离散事件仿真的关键组成部分。在所有 Excel 电子表格解答中，时钟时间都标注在仿真表相应列标题的正上方，从而便于你区分活动时间（例如服务时间或到达间隔时间）与计算得到的时钟时间。在所有情形下，一个时钟时间都表示某事件发生的时间，例如排队模型中的到达、服务开始或服务完成。

### 例 3：随机到达时间

打往该电话信息服务（其服务时间在例 2 中定义）的电话，其发生时刻是随机的，由一个离散分布定义，其中到达间隔时间的取值为 1、2、3 或 4 分钟，且概率相同。我们在此的目的，是展示如何同时生成到达间隔时间与到达时间。

由于本例只有一个事件，即到达事件，它是我们第一个动态的、事件驱动的模型示例（尽管很简单）。动态（dynamic）简单地说就是基于时间的，系统状态随时间变化；动态、事件驱动（event-based）则意味着模型跟踪事件发生随时间的推进过程。

到达间隔时间的生成方式与例 2 中服务时间的生成方式完全相同。不过，电子表格还提供了第二个 VBA 函数，以简化从离散均匀分布中随机生成数值。均匀（uniform）意味着所有取值的概率相同；本例中，所有到达间隔时间都可能是 1、2、3 或 4 分钟，其出现概率各为 0.25。对于这样的离散均匀分布，电子表格提供了 VBA 函数 DiscreteUniform()，它带两个参数：或者先给出含期望整数区间最小值的单元格、再给出含最大值的单元格，或者直接给出实际的上下限，例如：

$$
= \text{DiscreteUniform} (G 5, G 6)
$$

 或 

$$
= \text{DiscreteUniform} (1, 4)
$$

 DiscreteUniform(low, high) 函数由一行 VBA 代码实现，如下：

$$
\text{DiscreteUniform} = \text{low} + \text{Int((high - low + 1)*Rnd()}
$$

 对本例而言即为：

$$
\text{DiscreteUniform} = 1 + \text{Int} (4 * \text{Rnd()})
$$

 其中 Rnd() 生成 0 与 1 之间的随机数（在 VBA 中它总是小于 1 但可以等于 0），而 Int() 是一个 VBA 函数，它对其参数截断（向下取整）。（注意，4 \* Rnd() 生成 0.0 到 3.99999 之间的实值，因此截断后再加 1，就得到所需的 1 到 4 的范围。为什么我们期望概率相等？）

表 3 给出了到达间隔时间的输入规格说明以及仿真表的一部分，取自电子表格模型 "Example2.3ArrivalTimes.xls"。如步骤 1（第 14 行）所示，我们假定第一次到达发生在仿真时刻 0。（第一位顾客的假定到达时间取决于模型；它可以是任何其他指定的时刻，既可以是常数，也可以是随机的。）到达间隔时间用 DiscreteUniform(1, 4) 函数随机且独立地生成。到达时间则由一个简单的公式计算得出，即把当前顾客的到达间隔时间加上前一位顾客的到达时间。（你应当核对一下电子表格中计算到达时间的公式。）

**表 3 到达时间的输入规格说明与仿真表**

|  | A | B | C | D |
|---|---|---|---|---|
| 4 |  | 到达间隔时间（分钟） |  |  |
| 5 |  |  |  |  |
| 6 |  | 最小值 | 1 |  |
| 7 |  | 最大值 | 4 |  |
| 8 |  |  |  |  |
| 9 |  | 来电者人数= | 25 |  |
| 10 |  | 仿真表 |  |  |
| 11 |  | 步骤 | 活动 | 时钟 |
| 12 |  | 来电者 | 到达间隔时间 | 到达时间 |
| 13 |  |  |  |  |
| 14 |  | 1 |  | 0 |
| 15 |  | 2 | 2 | 2 |
| 16 |  | 3 | 1 | 3 |
| 17 |  | 4 | 3 | 6 |

例 3 是第一个包含事件的示例，该事件即到达事件，以及与该事件相关联的事件时间，即到达时刻的 CLOCK 时间。注意表 3 把活动时间与 CLOCK 时间区分开来。

### 1.7 电子表格仿真的一个框架

仿真表提供了结构，而下面列出的步骤则提供了开展电子表格仿真的一套通用准则。我们已经见过若干仿真表的示例，其中例 3 中的表 3 是目前最好的一例——因为它包含到达事件的时钟时间，从而使它成为最简单的动态事件驱动仿真之一。

首先，模型开发者必须定义模型输入、系统状态和模型输出。输入是那些（通常）独立于其他系统特征而定义的外生变量；例子包括掷硬币出现正面的概率、排队系统中的服务时间分布与到达间隔时间分布，以及库存模型中的需求分布。这些输入规格说明用于生成随机活动时间和其他随机变量。其他一些输入，例如策略或决策参数，可能是常数。

输出用于计算系统性能度量（也称为模型响应，或简称响应）。例如，一个输出可以是某位顾客在队列中的等待时间，或者库存系统中某笔交易的成本。相应的响应则可能是队列中的平均等待时间，或库存系统中单位时间的平均成本。

每张仿真表都是针对手头问题专门定制的。表中的每一列属于以下几种类型之一：

1. 与某个模型输入相联系的活动时间；

2. 由某个模型输入所定义的其他任何随机变量；

3. 一个系统状态；

4. 一个事件，或某个事件的时钟时间；

5. 一个模型输出；

6. 有时是一个模型响应。

对于此处的电子表格仿真，上面列出的这些项目已经足够。

一般而言，模型响应（也称为系统性能度量）是在仿真完成之后、在仿真表之外计算的。例如，在一个等待线模型中，每位顾客在队列中的延迟是一个模型输出，而所有顾客的平均延迟则是一个模型响应。模型输出在仿真表中跟踪；响应则在此之后计算。

活动时间是由常数或统计分布（例如服务时间或到达间隔时间）所定义的时间。活动是一段时间长度，与之相对的事件时间则是时钟时间。

系统状态例如包括排队模型中的服务台状态（忙碌或空闲）以及库存模型中的库存水平。一般来说，系统状态的集合包括仿真随时间推进（或从一步推进到下一步）所需的任何信息，以及任何能使仿真逻辑更简单、更清晰、因而更易理解的量或信息。在模型逻辑正确的前提下，尽可能简单有助于理解，而这种理解极有价值，尤其是在试图验证模型逻辑确实正确的时候。

**表 4 通用仿真表**

| 步骤 | 活动与系统状态 | 输出 |  |  |  |  |  |
|---|---|---|---|---|---|---|---|
| \(x_{i1}\) | \(x_{i2}\) | \(\cdots\) | \(x_{ij}\) | \(\cdots\) | \(x_{ip}\) | \(y_i\) |  |
| 1 |  |  |  |  |  |  |  |
| 2 |  |  |  |  |  |  |  |
| 3 |  |  |  |  |  |  |  |
| . |  |  |  |  |  |  |  |
| . |  |  |  |  |  |  |  |
| . |  |  |  |  |  |  |  |
| n |  |  |  |  |  |  |  |

仿真表中的每一行代表仿真中的一个步骤，通常与一个或多个事件的发生，或者某个实体（例如一位顾客）在系统中推进的过程相联系。仿真表的设计使得：给定的某一步只依赖于模型输入，和（或）当前这一步之前的一个或多个步骤，以及先前已计算出的数值。

一般而言，模型开发者按如下准则开发并运行电子表格仿真或手工仿真：

1. 确定仿真中各个输入的特征。经常遇到的输入类型包括常数以及概率分布。

2. 确定与问题相关的活动、事件和系统状态。

3. 确定模型响应或概括性的性能度量，视目标以及为实现这些目标而须回答的具体问题而定。

4. 确定计算模型响应所需的模型输出。

5. 构建一张仿真表。仿真表的一种通用表示见表 4。在该表中，有 p 个活动与系统状态，\(x_{i j} , j = 1 , 2 , \dotsc , p .\)，以及每个仿真步骤 \(i = 1 , 2 , \dots , n\) 的一个输出 y_i。通过填写步骤 1 的数据来初始化该表。在构建电子表格模型时，最初的一个或若干步骤（通常只有步骤 1）与后续步骤不同，后续步骤由公式定义。

6. 对每个步骤 i，先为各活动生成一个取值，然后计算系统状态和输出 y。随机输入值可以通过从步骤 1 中所选定的分布中抽样得到。一个输出通常依赖于输入、系统状态，也可能依赖于其他输出。

7. 仿真结束时，利用各输出计算模型响应或性能度量。

在以下各节中，我们将看到若干在专门设计的仿真表中完成的电子表格仿真示例。这些示例包括一个掷硬币游戏的示例，以及排队、库存、可靠性和网络分析等应用领域中的示例。这个掷硬币游戏提出的问题很简单，但其答案我们相信会让大多数读者感到意外。

两个排队示例分别说明单服务台和双服务台系统。第一个库存示例涉及一个有闭式解的问题；因此该仿真解答可以与数学解答作比较。第二个库存示例则涉及经典的订购至水平模型。

其余示例包括一个命中目标的模型、一个可靠性模型、一个用于估计某个复杂随机变量（即提前期需求）分布的蒙特卡罗模型，最后还有一个项目作业网络模型。命中目标模型引入了正态分布随机数的使用。

## 2 掷硬币游戏

在掷硬币的游戏中，你和对手各自领先的时间大约占一半（即大约有一半的投掷次数上双方各领先一半）的可能性，是否大于其中一方在大多数时间领先的可能性？

借助下面的示例，我们更详细地说明电子表格解答的布局、两个主要工作表（「One Trial」和「Experiment」）、用于运行仿真的按钮，以及所使用的术语：试验（或重复）与实验。模型逻辑的细节（所用的公式、调用的函数等等）在电子表格本身中都有解释，这些都可以在本书网站 www.bcnn.net 上找到。

### 例 4：掷硬币游戏

Charlize 为她的朋友 Harry 和 Tom 掷硬币，恰好掷 100 次。掷出正面时 Harry 赢 \$1.00，Tom 输 \$1.00；掷出反面时 Tom 赢 \$1，由 Harry 支付。Charlize 在他们玩这个游戏的过程中记录各自的盈亏。在任何一次投掷中，可能是 Harry 领先，也可能是 Tom 领先，还可能双方打平。

在任意一局游戏（一局游戏恰好是 100 次投掷）中，两人各自领先大约占一半投掷次数的可能性，是否大于其中一人在大多数时间领先的可能性？

注意，我们问的不是游戏结束时的输赢；也不是问正面与反面之比（正如你所知，从长远看应当接近 50:50）。这个问题关心的是在整个游戏过程中，Harry 或 Tom 领先或落后的频繁程度。还要注意，问题中的这两种可能性并不是一局游戏仅有的两种可能结果；因此，它们的概率不会加起来等于 1。

具体地说，下面哪一种更可能，以及大约相差多少？

(A) Harry 和 Tom 各自领先大约占一半的投掷次数；具体地说，Harry 在 45 到 55 次投掷中领先，因而 Tom 也在 45 到 55 次投掷中领先；或者

(B) Harry 在 95 次或更多次投掷中领先；或者

(C) Harry 在 5 次或更少次投掷中领先。

先猜一猜，再看结果。表 5 给出了输入规格说明以及仿真表中的前 14 次投掷，取自电子表格模型 "Example2.4CoinTossGame.xls"，使用默认种子（12,345）并点击「重置种子并运行」按钮。「One Trial」工作表中含有一局 100 次投掷游戏的逻辑与结果。

**表 5 掷硬币游戏的输入规格说明与仿真表**

|  | A | B | C | D |
|---|---|---|---|---|
| 4 |  | 硬币 | 概率 | 累积概率 |
| 5 |  |  |  |  |
| 6 |  |  |  |  |
| 7 |  | 正面 | 0.50 | 0.50 |
| 8 |  | 反面 | 0.50 | 1.00 |
| 9 |  |  |  |  |
| 10 |  | 投掷次数 = 100 |  |  |
| 11 |  | 仿真表 |  |  |
| 12 |  | 步骤 | 状态 | 输出 |
| 13 |  | 投掷 | 投掷结果 | Harry 的盈亏 |
| 14 |  | 1 | T | -\$1 |
| 15 |  | 2 | H | 0 |
| 16 |  | 3 | T | -1 |
| 17 |  | 4 | H | 0 |
| 18 |  | 5 | H | 1 |
| 19 |  | 6 | T | 0 |
| 20 |  | 7 | H | 1 |
| 21 |  | 8 | T | 0 |
| 22 |  | 9 | T | -1 |
| 23 |  | 10 | H | 0 |
| 24 |  | 11 | T | -1 |
| 25 |  | 12 | H | 0 |
| 26 |  | 13 | T | -1 |

如果你查看电子表格解答中的「One Trial」工作表（采用 RNG 默认种子的那次试验），你会看到在第 13 次投掷时 Harry 的盈亏为负（他在输钱），并且在这局游戏余下的过程中一直为负；事实上，在这局游戏中，Tom 在 100 次投掷中的 92 次上领先，而他最终却只赢了 \$6.00。Harry 只在 100 次投掷中的 2 次上领先，双方打平的情形有 6 次（最后一次是在第 12 次投掷时）。

现在，我们不能只靠玩一局游戏来回答所提出的问题（这相当于点击一次「生成新试验」按钮，因为每次试验就是一局 100 次投掷的游戏）。要玩很多局，我们使用「Experiment」工作表，并把试验次数（单元格 C3）设为 400。（就我们目前的目的而言，任何相当大的试验次数都可以。）然后点击标有「重置种子并运行」的按钮，它把种子设为默认值（12,345）并用该种子运行实验。运行实验意味着「One Trial」工作表被执行 400 次，每次试验之后，其结果都被记录到「Experiment」工作表中。该电子表格中包含关于公式、函数以及实现模型逻辑和收集所需输出的其他技术细节的注释与说明。

**表 6 掷硬币游戏的实验工作表**

|  | A | B | C | D | E | F |
|---|---|---|---|---|---|---|
| 3 | 试验次数：400 |  |  |  |  |  |
| 4 | 性能度量链接： |  |  |  |  |  |
| 5 | 度量名称：链接 |  |  |  |  |  |
| 6 | Harry 领先的 100 次投掷次数：31 |  |  |  |  |  |
| 7 |  |  |  |  |  |  |
| 8 |  |  |  |  |  |  |
| 9 | 响应表 |  | 多次试验汇总 |  |  |  |
| 10 | 试验 | Harry 领先的投掷次数 |  | 组（区间） | 频数 | 相对频数 |
| 11 | 1 | 2 |  | 0 到 5 | 73 | 18.3% |
| 12 | 2 | 91 |  | 45 到 55 | 28 | 7.0% |
| 13 | 3 | 1 |  | 95 到 100 | 46 | 11.5% |
| 14 | 4 | 76 |  |  |  |  |
| 15 | 5 | 17 |  |  |  |  |
| 16 | 6 | 11 |  | 平均值 | 48.2 |  |
| 17 | 7 | 94 |  | 中位数 | 50 |  |
| 18 | 8 | 95 |  |  |  |  |
| 19 | 9 | 78 |  | 最小值 | 0 |  |
| 20 | 10 | 0 |  | 最大值 | 100 |  |
| 21 | 11 | 40 |  |  |  |  |

表 6 给出了实验规格说明（400 次试验与响应）、前 11 次试验（前 11 局游戏）的响应表，以及多次试验汇总。图 2 是根据多次试验汇总中的频数结果绘制的直方图。

基于这次包含 400 次试验（或重复）的实验，结果可能会让许多人感到意外：

Harry 在大多数时间落后 在 400 局游戏中的 73 局（18.3%）里，Harry 只在 100 次投掷中的 5 次或更少次上领先。

Harry 与 Tom 差不多 在 400 局游戏中的 28 局（7.0%）里，Harry 与 Tom 大致平分领先。

Harry 在大多数时间领先 在 400 局游戏中的 46 局（11.5%）里，Harry 在 100 次投掷中的 95 次或更多次上领先。

![图 2 Harry 领先机会的直方图](img/ch02/img_80a19a57792d2bb7.jpg)

**图 2 Harry 领先机会的直方图。**

当然，这三种可能性并非全部的可能性；这正是频数之和不等于 400 局、也不等于 100% 的原因。你也许会想：是什么解释了"0 到 5"这一结果看起来比"95 到 100"更有可能，而两者又都比 Tom 与 Harry 大致平分领先更有可能？这一差异有多大？这些结果是否与硬币是"公平"的、正反面概率相等的假设相矛盾？这些结果对 Harry 或 Tom 最终的盈亏又说明了什么？

你可以多次点击「运行实验」按钮做实验，看看结果在一次实验与下一次实验之间如何变化。每次实验包含 400 局游戏，直觉上结果不会有太大变化。一次实验运行与另一次实验运行是完全一样的，只不过使用了一组新的随机数。（你随时可以点击「重置种子并运行」按钮，回到这里所给出的原始随机数和原始结果。）要验证这里基于 400 局游戏（即 400 次试验）引用的结果，你可以运行一次试验次数非常多的实验，比如 4000 次试验。随着试验次数增加，精度会提高。

## 3 电子表格中的排队仿真

我们先描述排队模型或等待线模型的组成部分。随后给出两个示例：一个单服务台队列和一个双服务台队列；两者都是动态的、事件驱动的模型。

![图 3 单通道排队系统](img/ch02/img_b391e11732e70af6.jpg)

**图 3 单通道排队系统。**

### 3.1 等待线模型

一个排队系统（queueing system）由其呼叫总体（calling population）、到达的性质、服务机制、系统容量以及排队规则（queueing discipline）加以描述。图 3 描绘了一个简单的单通道排队系统。呼叫总体的成员可以是顾客、电话呼叫，或是修理铺里的活件；这里我们称它们为单位（unit）。

在本章的简单模型中，我们假定呼叫总体在潜在意义上是无限的，且到达率是恒定的。也就是说，无论已经到达或仍留在系统中的单位有多少，下一次到达之前的到达间隔时间都服从同一个统计分布。到达的服务请求以随机的方式一个接一个地发生；一旦加入等待线，它们最终都会得到服务。此外，服务时间按某个概率分布取某个随机长度，该分布不随时间变化。系统容量没有限制，这意味着任何数量的单位都可以在排队线上等待。最后，单位按其到达顺序接受服务（通常称为 FIFO：先进先出，或 FCFS：先到先服务），服务由一个服务台或多个并行服务台完成。

服务与到达分别由服务时间的分布和到达之间时间的分布来定义，正如例 2 与例 3 所示。对于任何简单的单通道或多通道队列，总的有效到达率必须小于总服务率，否则等待线将无限增长。当队列无限增长时，称其为爆炸的或不稳定的。（在某些可重入的排队网络中，单位在最终离开系统之前会多次返回同一个服务台，此时"到达率小于服务率"这一条件未必能保证稳定性。更多解释参见 Harrison 与 Nguyen [1995]。有趣的是，这类不稳定性最初并非在理论中被发现，而是在半导体制造厂的实际制造过程中被注意到。）也可能出现更为复杂的情形——例如在短时间内到达率大于服务率，或者具有路由的队列网络。然而，本章只讨论最基本的队列。

系统状态、事件和仿真时钟的概念如何应用于排队模型？排队模型的系统状态通常由系统中的单位数以及服务台状态（忙碌或空闲）构成。事件是一组引起系统状态瞬时变化的情形。在单通道（或任何多通道）排队系统中，事件包括到达事件、服务开始事件和离开事件（也称服务完成事件）。仿真时钟用于跟踪仿真时间。

当离开事件发生时，仿真按图 4 的流程图所示的方式推进。离开事件引起系统状态的变化。如果队列中没有单位在等待，服务台状态变为空闲；如果有一个或多个单位在等待，则下一个单位开始服务。在这两种情形下，系统中的单位数都减少一个。

![图 4 离开事件流程图](img/ch02/img_4682cc5014854604.jpg)

**图 4 离开事件流程图。**

到达事件的流程图见图 5。该单位发现服务台要么空闲要么忙碌；因此，该单位要么立即开始服务，要么进入队列。该单位遵循图 6 所示的行动路线。如果服务台忙碌，该单位进入队列。如果服务台空闲且队列为空，该单位开始服务。服务台不可能在队列非空时处于空闲状态。

![图 5 到达事件流程图](img/ch02/img_5a9b7652758c5bb3.jpg)

**图 5 到达事件流程图。**

|  | 队列状态 |  |  |
|---|---|---|---|
| 非空 | 空 |  |  |
| 服务台状态 | 忙碌 | 进入队列 | 进入队列 |
| 空闲 | 不可能 | 开始服务 |  |
| 服务台结果 | 忙碌 |  | 不可能 |
| 空闲 | 不可能 |  |  |

**图 6 到达时单位可能采取的行动。**

**图 7 服务完成后的服务台结果。**

服务完成之后，服务台要么变为空闲，要么在有下一个单位的情况下保持忙碌。这两种结果与队列状态之间的关系见图 7。如果队列非空，另一个单位就占用服务台，服务台于是保持忙碌状态。如果队列为空，服务台变为空闲。这两种可能性表现为图 7 中的阴影部分。

在大多数仿真中，许多事件发生在随机时刻，这种随机性模拟着现实生活中的不确定性。例如，无法确切知道下一位顾客何时会到达杂货店收银台，也无法确切知道银行柜员完成一笔业务要花多长时间。在这些情形下，要根据收集并分析得到的数据，或者根据主观估计与假设，建立数据的统计模型。就我们的目的而言，我们假定到达间隔时间与服务时间服从简单的离散概率分布。

那么，这里所描述的事件是如何在仿真时间中发生的呢？排队系统的仿真通常需要维护一份事件列表（event list），用以判定接下来会发生什么。事件列表跟踪各类事件将要发生的未来时刻。本章通过显式地跟踪每一个单位来简化仿真。注意，这种临时拼凑的做法对于电子表格中的简单仿真而言是令人满意的，但对于更为复杂的现实世界仿真（以及更为有趣的学生仿真）而言，更为系统化的做法是必不可少的。

表 7 是专门为按先进先出（FIFO）原则为顾客服务的单通道队列设计的一张仿真表。它跟踪每个事件发生时的时钟时间。（注意表中不含输出列；这些将在下一个示例中加入。）这张仿真表依据假定的一系列到达间隔时间计算到达时间：

2 4 1 2 6

B、C 和 E 列是仿真时钟时间。B 列（到达时间）仅由所生成的到达间隔时间计算得到。E 列（离开事件或服务完成事件的时间）则可由 D 列所生成的输入以及 C 列中服务开始的时刻容易地计算出来。因此我们看到，计算 C 列（服务开始时间）是主要难点。对一个到达的单位来说，服务开始的时间是以下两个时刻中较晚的一个：其到达时间，或前一单位的离开时间。（对于双服务台以及其他多服务台队列而言，服务完成时间的计算就没那么容易了。）

如表 7 所示，第一位顾客在时钟时刻 0 到达并立即开始服务，该服务需要两分钟。服务在时钟时刻 2 完成。第二位顾客在时钟时刻 2 到达，在时钟时刻 3 服务结束。注意第四位顾客在时钟时刻 7 到达，但服务直到时钟时刻 9 才能开始。这是因为顾客 3 直到时钟时刻 9 才服务结束。

**表 7 强调时钟时间的仿真表**

| A | B | C | D | E |
|---|---|---|---|---|
| 顾客编号 | 到达时间（时钟） | 服务开始时间（时钟） | 服务时间（活动） | 服务结束时间（时钟） |
| 1 | 0 | 0 | 2 | 2 |
| 2 | 2 | 2 | 1 | 3 |
| 3 | 6 | 6 | 3 | 9 |
| 4 | 7 | 9 | 2 | 11 |
| 5 | 9 | 11 | 1 | 12 |
| 6 | 15 | 15 | 4 | 19 |

**表 8 事件的时间顺序**

| 事件类型 | 顾客编号 | 时钟时间 |
|---|---|---|
| 到达 | 1 | 0 |
| 离开 | 1 | 2 |
| 到达 | 2 | 2 |
| 离开 | 2 | 3 |
| 到达 | 3 | 6 |
| 到达 | 4 | 7 |
| 离开 | 3 | 9 |
| 到达 | 5 | 9 |
| 离开 | 4 | 11 |
| 离开 | 5 | 12 |
| 到达 | 6 | 15 |
| 离开 | 6 | 19 |

到达事件与离开事件的时间顺序见表 8 和图 8。表 8 按时钟时间排序，但事件未必按顾客编号排序。事件的时间顺序是离散事件仿真方法的基础。

图 8 描绘了系统中的顾客数随时间的变化。它是表 8 事件列表的直观图像。顾客 1 在系统中的时段是从时钟时刻 0 到时钟时刻 2。顾客 2 在时钟时刻 2 到达，在时钟时刻 3 离开。从时钟时刻 3 到时钟时刻 6，系统中没有顾客。在某些时间段内系统中有两名顾客，例如在时钟时刻 8，顾客 3 和顾客 4 都在系统中。此外，也存在事件同时发生的情形，例如在时钟时刻 9，顾客 5 到达而顾客 3 离开。

![图 8 系统中的顾客数](img/ch02/img_affc75aeda5a8636.jpg)

**图 8 系统中的顾客数。**

### 3.2 单服务台队列的仿真

第一个排队示例是单服务台队列。我们用它说明如何计算典型的输出与性能度量，例如顾客的平均等待时间以及服务台空闲的时间比例。

### 例 5：杂货店收银台——单服务台队列

一家小型杂货店只有一个收银台。顾客在随机时刻到达收银台，相邻两次到达之间的间隔在 1 到 8 分钟之间。我们假定到达间隔时间为整数值，且 8 个取值各以相等概率出现；这就是一个离散均匀分布，如表 9 所示。服务时间在 1 到 6 分钟之间变化（同样为整数值），其概率如表 10 所示。我们的目标是通过仿真 100 位顾客的到达与服务来分析该系统，并计算排队模型中若干典型的性能度量。

实际上，100 位顾客所对应的样本量可能太小，不足以得出可靠的结论。视我们的目标而定，可以通过增大样本量（顾客数）或进行多次试验（或重复）来提高结果的精度。第二个问题是初始条件。一个从空系统开始的杂货店仿真，除非意在模拟系统从启动开始的运行过程，或者意在模拟直至达到稳态运行，否则可能并不现实。这里，为使计算保持简单，起始条件设定为杂货店为空，任何相关顾虑都不予考虑。

**表 9 到达之间时间的分布**

|  | G | H |
|---|---|---|
| 4 | 到达间隔时间（分钟） |  |
| 5 |  |  |
| 6 | 最小值 | 1 |
| 7 | 最大值 | 8 |

**表 10 服务时间的分布**

|  | A | B | C | D |
|---|---|---|---|---|
| 4 |  | 服务时间（分钟） | 概率 | 累积概率 |
| 5 |  |  |  |  |
| 6 |  |  |  |  |
| 7 |  | 1 | 0.10 | 0.10 |
| 8 |  | 2 | 0.20 | 0.30 |
| 9 |  | 3 | 0.30 | 0.60 |
| 10 |  | 4 | 0.25 | 0.85 |
| 11 |  | 5 | 0.10 | 0.95 |
| 12 |  | 6 | 0.05 | 1.00 |

如何由服务时间分布与到达间隔时间分布生成服务时间和到达时间，已分别在第 1.5 节和第 1.6 节中讨论；同样的方法也适用于此处及后续各示例。对于手工仿真，请使用附录中的表 A.1 里的随机数。良好的做法是：从表 A.1 中的某个随机位置开始，沿某个方向系统地推进，并且在给定问题中绝不使用同一串数字。如果反复使用同一种模式，统计偏差或其他怪异效应可能会影响结果。

表 11 给出了单通道队列仿真表的一部分，取自 "Example2.5SingleServer.xls"；其中还给出了若干选定列的合计与平均值（位于仿真表上方）。第一步是初始化该表，即填写第一位顾客所在的各单元格——假定该顾客在时刻 0 到达并立即开始服务。这位顾客在 2 分钟时刻结束服务（第一个随机服务时间）。在第一位顾客之后，表中后续各行依据当前顾客随机生成的到达间隔时间与服务时间、前一位顾客的完成时间，以及对电子表格解答而言仿真表中各单元格公式所包含的仿真逻辑来确定。

接下来，第二位顾客在 5 分钟时刻到达并立即开始服务。因此，第二位顾客在队列中没有等待，而服务台有 3 分钟的空闲时间。

往下看，我们看到第五位顾客在 16 分钟时刻到达，发现服务台忙碌。在 18 分钟时刻，前一位顾客离开，第五位顾客开始服务，已在队列中等待了 2 分钟。这一过程对全部 100 位顾客持续进行。"Example2.5SingleServer.xls" 中的电子表格解答包含实现此处所述步骤的公式。为了最好地理解这些公式，你应当先尝试手工仿真若干步骤、弄清楚其中的逻辑，然后再研究这些公式是如何实现该逻辑的。

表中增加了 G、I 和 J 三列，用来收集三个模型输出，即每位顾客在队列中的时间、在系统中的时间，以及从前一位顾客离开到当前顾客到达之间服务台的空闲时间（如有）。这些模型输出又进一步用于计算若干性能度量，即顾客在队列中的平均时间和在系统中的平均时间，以及服务台空闲的时间比例。

为了计算系统性能度量，可以使用表 11 中已计算出的响应，该表给出了到达间隔时间、服务时间、顾客在队列中的等待时间、顾客在系统中所花的时间以及服务台空闲时间的合计与平均值。这些结果取自电子表格模型「One Trial」工作表上使用默认 RNG 种子（12,345）时的第一次试验运行。利用这些结果，我们可以说明等待线模型中若干典型的性能度量：

**表 11 食品店仿真的模型响应与仿真表（前 11 位顾客）**

|  | A | B | C | D | E | F | G | H | I | J |
|---|---|---|---|---|---|---|---|---|---|---|
| 15 |  | 合计 | 420 |  | 320 |  | 163 |  | 483 | 106 |
| 16 |  | 平均值 | 4.24 |  | 3.20 |  | 1.63 |  | 4.83 | 1.07 |
| 17 |  | 顾客人数= | 100 |  |  |  |  |  |  |  |
| 18 |  | 仿真表 |  |  |  |  |  |  |  |  |
| 19 |  | 步骤 | 活动 | 时钟 | 活动 | 时钟 | 输出 | 时钟 | 输出 | 输出 |
| 20 |  | 顾客 | 到达间隔时间（分钟） | 到达时间 | 服务时间（分钟） | 服务开始时间 | 队列中等待时间（分钟） | 服务结束时间 | 顾客在系统中所花时间（分钟） | 服务台空闲时间（分钟） |
| 21 |  |  |  |  |  |  |  |  |  |  |
| 22 |  |  |  |  |  |  |  |  |  |  |
| 23 |  | 1 | 0 | 0 | 2 | 0 | 0 | 2 | 2 |  |
| 24 |  | 2 | 5 | 5 | 2 | 5 | 0 | 7 | 2 | 3 |
| 25 |  | 3 | 5 | 10 | 4 | 10 | 0 | 14 | 4 | 3 |
| 26 |  | 4 | 4 | 14 | 4 | 14 | 0 | 18 | 4 | 0 |
| 27 |  | 5 | 2 | 16 | 3 | 18 | 2 | 21 | 5 | 0 |
| 28 |  | 6 | 8 | 24 | 2 | 24 | 0 | 26 | 2 | 3 |
| 29 |  | 7 | 7 | 31 | 3 | 31 | 0 | 34 | 3 | 5 |
| 30 |  | 8 | 8 | 39 | 5 | 39 | 0 | 44 | 5 | 5 |
| 31 |  | 9 | 5 | 44 | 1 | 44 | 0 | 45 | 1 | 0 |
| 32 |  | 10 | 2 | 46 | 6 | 46 | 0 | 52 | 6 | 1 |
| 33 |  | 11 | 1 | 47 | 4 | 52 | 5 | 56 | 9 | 0 |

1. 顾客的平均等待时间为 1.63 分钟，计算如下：

$$
\begin{array}{r l} \text{Average   waiting   time} & = \frac {\text{total   time   customers   wait   in   queue   (minutes)}}{\text{total   numbers   of   customers}} \\ & = \frac {163}{100} = 1. 63 \text{minutes} \end{array}
$$


2. 顾客需要在队列中等待的概率为 0.46，计算如下：

$$
\begin{array}{r l} \text{Probability   (wait)} & = \frac {\text{numbers   of   customers   who   wait}}{\text{total   number   of   customers}} \\ & = \frac {46}{100} = 0. 46 \end{array}
$$


3. 服务台约有 25% 的时间处于空闲，计算如下：

$$
\begin{array}{r l} \text{Probability   of   idle} & = \frac {\text{total   idle   time   of   server   (minutes)}}{\text{total   run   time   of   simulation   (minutes)}} \\ & = \frac {106}{426} = 0. 25 \end{array}
$$

 因此，服务台约有 75% 的时间处于忙碌。

4. 平均服务时间为 3.20 分钟，计算如下：

$$
\begin{array}{r l} \text{Average   service   time} & = \frac {\text{total   service   time   (minutes)}}{\text{total   number   of   customers}} \\ & = \frac {320}{100} = 3. 20 \text{minutes} \end{array}
$$

 通过求服务时间分布的均值，可以把这一结果与期望服务时间作比较，所用的公式为 

$$
E (S) = \sum_{s = 0} ^ {\infty} s p (s)
$$

 把期望值公式应用于表 10 中的分布，得到 

$$
\begin{array}{l} \text{Expected   service   time} = \\ 1 (0. 10) + 2 (0. 20) + 3 (0. 30) + 4 (0. 25) + 5 (0. 10) + 6 (0. 05) = 3. 2 \text{minutes} \end{array}
$$

 期望服务时间与本次仿真试验所得的平均服务时间完全相同。一般说来并非如此；这里纯属巧合。

5. 到达之间的平均时间为 4.24 分钟，计算如下：

$$
\begin{array}{r l} \text{Average   time   between} & = \frac {\text{sum   of   all   times}}{\text{between   arrivals   (minutes)}} \\ \text{arrivals   (minutes)} & = \frac {420}{99} = 4. 24 \text{minutes} \end{array}
$$

 由于第一次到达发生在时刻 0，我们除以 99 而不是 100。通过求端点为 \(a = 1\) 和 \(b = 8 .\) 的离散均匀分布的均值，可以把这一结果与到达之间的期望时间作比较。均值为 

$$
E (A) = \frac {a + b}{2} = \frac {1 + 8}{2} = 4. 5 \mathrm{minutes}
$$

 到达之间的期望时间略高于平均值。不过，随着仿真时间变长，到达之间时间的平均值会越来越接近理论均值 \(E ( A )\)。

6. 需要等待的那些顾客的平均等待时间为 3.54 分钟，计算如下：

$$
\begin{array}{r l} \text{Average   waiting   time   of} & \\ \text{those   who   wait} & = \frac {\text{total   time   customers   wait   in   queue   (minutes)}}{\text{total   number   of   customers   who   wait}} \\ (\text{minutes}) & = \frac {163}{46} = 3. 54 \text{minutes} \end{array}
$$

 注意，这不同于所有顾客的平均等待时间，因为许多顾客（本次仿真中有 54 位）在开始接受服务之前完全没有等待。在全部顾客中，平均等待时间为 \(\frac {163} {100} ~ = ~ 1 . 63\) 分钟。需要等待的顾客人数 46，来自仿真表中"队列中等待时间"一栏下的 G 列，方法是数出非零等待时间的个数；也可以从图 9 得到，图 9 展示了这 100 个等待时间的直方图。

7. 顾客在系统中所花的平均时间为 4.83 分钟，可用两种方式计算。第一种：

$$
\begin{array}{r l} \text{Average   time   customer} & = \frac {\text{total   time   customers   spend   in   the}}{\text{system   (minutes)}} \\ & = \frac {483}{100} = 4. 83 \text{minutes} \end{array}
$$

 第二种方式是认识到下述关系必定成立：

$$
\begin{array}{c c c c} \text{Average time} & \text{average time} & \text{average time} \\ \text{customer spends} & \text{customer spends} & \text{customer spends} \\ \text{in the system} & = & \text{waiting in the} + & \text{in service} \\ (\text{minutes}) & \text{queue (minutes)} & & (\text{minutes}) \end{array}
$$


**个体顾客队列等待时间的频率分布**

![图 9 队列中等待时间的频率](img/ch02/img_ae4cfe46c72646b8.jpg)

**图 9 队列中等待时间的频率。**

由结论 1 和结论 4 可得：

顾客在系统中所花的平均时间 = 1.63 + 3.20 = 4.83 分钟

现在我们感兴趣的是开展一项实验，以回答前 100 位顾客在队列中的平均时间如何逐日变化这一问题。把「One Trial」工作表以 100 位顾客运行一次，相当于一天。把它运行 50 次，则相当于使用「Experiment」工作表并把"试验次数"设为 50，其中每次试验代表一天。

对于例 5，图 9 给出了第一次试验中前 100 位顾客在队列中等待时间的频率。（注意：在本章余下部分的所有直方图中，组（区间）的上限都标注在 x 轴的图例上，即使图例显示在组（区间）的中央。）如直方图所示，54% 的顾客不必等待，而在需要等待的顾客中，34% 等待不到 4 分钟（但大于零分钟）。这些结果和图 9 中的直方图代表的是一次试验的结果，仅凭它们本身并不构成对我们问题的回答。

为回答我们的问题，我们使用「Experiment」工作表，把"试验次数"设为 50，然后点击标有「重置种子并运行」的按钮。这样就得到 50 个平均队列时间的估计值（每个都是对 100 位顾客取平均）。50 次试验上的总体平均等待时间为 1.32 分钟。图 10 是 50 次试验所得 50 个平均等待时间的直方图。它显示了在 50 天的样本中，平均等待时间如何逐日变化。

通过在电子表格解答上做实验并完成相关习题，你可以发现随机性以及所假定输入数据的影响。例如，如果你运行 400 次试验或 10 次试验，而不是 50 次，会怎样？图 10 中分布的形状会因试验次数而改变多少？为什么图 9 中的取值范围（0 到 12 分钟）要比图 10 中的取值范围（0.5 到 3 分钟）宽得多？

![图 10 50 次试验平均等待时间的频率分布](img/ch02/img_4c01584442458226.jpg)

**图 10 50 次试验平均等待时间的频率分布。**

### 3.3 双服务台队列的仿真

第二个排队示例是一个具有两个服务台且两者服务时间不同的队列。

### 例 6：Able–Baker 呼叫中心问题

某计算机技术支持中心由 Able 和 Baker 两个人负责，他们接听来电、尝试解答问题并解决计算机故障。两次来电之间的间隔在 1 到 4 分钟之间，其分布如表 12 所示。Able 更有经验，能够比 Baker 提供更快的服务。他们的服务时间分布分别见表 13 和表 14。

当两人都空闲时，由 Able 接听。当两人都忙碌时，来电转入等待；这基本上与任何排队模型中顾客必须在排队线上等待是一回事。一般而言，仿真中的到达时间、服务时间以及由此得出的时钟时间都应当是连续的（实值）变量；它们在此处（以及其他示例中）是离散的、整数值，纯粹是为了便于讲解以及理解前面所解释的随机变量生成技术。

**表 12 技术支持来电的到达间隔分布**

|  | B | C | D |
|---|---|---|---|
| 4 | 来电的到达间隔分布 |  |  |
| 5 | 到达间隔时间（分钟） | 概率 | 累积概率 |
| 6 |  |  |  |
| 7 |  |  |  |
| 8 | 1 | 0.25 | 0.25 |
| 9 | 2 | 0.40 | 0.65 |
| 10 | 3 | 0.20 | 0.85 |
| 11 | 4 | 0.15 | 1.00 |

**表 13 Able 的服务时间分布**

|  | F | G | H |
|---|---|---|---|
| 4 | Able 的服务时间分布 |  |  |
| 5 | 服务时间（分钟） | 概率 | 累积概率 |
| 6 |  |  |  |
| 7 |  |  |  |
| 8 | 2 | 0.30 | 0.30 |
| 9 | 3 | 0.28 | 0.58 |
| 10 | 4 | 0.25 | 0.83 |
| 11 | 5 | 0.17 | 1.00 |

**表 14 Baker 的服务时间分布**

|  | I | J | K |
|---|---|---|---|
| 4 | Baker 的服务时间分布 |  |  |
| 5 | 服务时间（分钟） | 概率 | 累积概率 |
| 6 |  |  |  |
| 7 |  |  |  |
| 8 | 3 | 0.35 | 0.35 |
| 9 | 4 | 0.25 | 0.60 |
| 10 | 5 | 0.20 | 0.80 |
| 11 | 6 | 0.20 | 1.00 |

仿真的推进方式与例 5 中的单服务台情形类似，只是其内在逻辑（仿真表中的电子表格公式）更为复杂，原因有二：有两个服务台，而且两个服务台的服务时间不同。完整的解答在工作簿 "Example2.6AbleBaker.xls" 中，表 15 展示了其中的一部分，并给出了估计来电者延迟（等待时间）所需的一些合计。该工作簿详细说明了仿真的各个步骤，包括计算每一步所依据的内在逻辑（公式）。（要理解这些步骤，你或许可以先做一次手工仿真，例如从一张包含步骤 1（第 20 行）以及到达间隔时间和服务时间各列的仿真表开始，尝试与表 15 中的结果对照。）

我们的目标是估计来电者在等待 Able 或 Baker 接听过程中所承受的延迟（等待时间）。为估计这一以及其他所期望的系统性能度量，我们对前 100 位来电者进行一次仿真。

从表 15 可见，假定来电者 1 在时钟时刻 0 到达。此时 Able 空闲，因此来电者 1 立即开始服务。服务时间 4 分钟，是按例 2 中的过程（投掷"随机飞镖"）由表 13 给出的分布生成的。因此，来电者 1 在时钟时刻 4 分钟完成服务，没有延迟。

**表 15 呼叫中心示例的仿真表（100 位来电者中的前 10 位）**

|  | B | C | D | E | F | G | H | I | J | K | L | M |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 13 | 合计 |  |  |  |  |  |  |  |  |  | 73 | 432 |
| 14 | 来电者人数= | 100 |  | 随机数种子 | 12,345 |  |  |  |  |  |  |  |
| 15 | 仿真表 |  |  |  |  |  |  |  |  |  |  |  |
| 16 | 步骤 | 活动 | 时钟 | 时钟 | 时钟 | 状态 | 活动 | 时钟 | 时钟 | 时钟 | 输出 | 输出 |
| 17 | 来电者编号 | 到达间隔时间（分钟） | 到达时间 | Able 可服务时刻 | Baker 可服务时刻 | 选中的服务台 | 服务时间（分钟） | 服务开始时间 | 服务完成时间 | 来电者延迟（分钟） | 在系统中的时间（分钟） |  |
| 18 | Able | Baker |  |  |  |  |  |  |  |  |  |  |
| 19 |  |  |  |  |  |  |  |  |  |  |  |  |
| 20 | 1 |  | 0 | 0 | 0 | Able | 4 | 0 | 4 |  | 0 | 4 |
| 21 | 2 | 1 | 1 | 4 | 0 | Baker | 4 | 1 |  | 5 | 0 | 4 |
| 22 | 3 | 1 | 2 | 4 | 5 | Able | 3 | 4 | 7 |  | 2 | 5 |
| 23 | 4 | 3 | 5 | 7 | 5 | Baker | 3 | 5 |  | 8 | 0 | 3 |
| 24 | 5 | 2 | 7 | 7 | 8 | Able | 5 | 7 | 12 |  | 0 | 5 |
| 25 | 6 | 2 | 9 | 12 | 8 | Baker | 5 | 9 |  | 14 | 0 | 5 |
| 26 | 7 | 1 | 10 | 12 | 14 | Able | 5 | 12 | 17 |  | 2 | 7 |
| 27 | 8 | 2 | 12 | 17 | 14 | Baker | 4 | 14 |  | 18 | 2 | 6 |
| 28 | 9 | 4 | 16 | 17 | 18 | Able | 2 | 17 | 19 |  | 1 | 3 |
| 29 | 10 | 1 | 17 | 19 | 18 | Baker | 3 | 18 |  | 21 | 1 | 4 |

![图 11 第一次试验中来电者延迟的频率](img/ch02/img_caeea34e8a2d2393.jpg)

**图 11 第一次试验中来电者延迟的频率。**

按例 3 中的过程由表 12 生成的 1 分钟到达间隔时间，使来电者 2 在时钟时刻 1 分钟到达。此时 Able 忙碌，因此来电者 2 由 Baker 服务。

来看来电者 3，他由 Able 从时钟时刻 4 分钟服务到时钟时刻 7 分钟；可以看到，来电者 3 在 2 分钟时刻到达时发现两个服务台都忙碌，因而不得不等待。

对于 100 位来电者的单次试验，图 11 中的频率图显示，100 位来电者中有 78 位（78%）没有延迟，16% 的延迟为 1 分钟，而少数几位不得不等待至多 3 分钟。使用不同随机数的其他试验，可能得到相似的结果，也可能得到大不相同的结果；唯有通过实验才能确定仿真 100 位顾客时所固有的变异性。

图 11 中来电者延迟的分布是右偏的。在「One Trial」工作表上反复点击「生成新试验」按钮，留意到来电者延迟分布常常有一条又长又细的尾巴伸向右侧。除了靠近零的高柱（高频数）之外，这对来电者延迟的性质意味着什么？一次试验提供的数据不足以据以形成可靠的结论，但这是使用电子表格的一大优势——变异性的影响相当明显——而且可以借助「Experiment」工作表把仿真重复任意多次。

在表 15 中，来电者延迟合计为 73 分钟，平均每位于来电者约 0.73 分钟。然而图 11 的直方图表明，单凭 0.73 分钟这一平均值，并不能很好地反映来电者要等待 Able 或 Baker 多久。虽然本次试验中有 78% 的顾客没有经历延迟，但少数顾客的延迟长达 3 分钟。这说明了根据模型的目标与目的选择恰当的系统性能度量的重要性。毕竟，不必等待的顾客不会抱怨等待，而那些等待很久的顾客则不会在乎平均值是 0.73 分钟。（实际上，在随后的一些试验中，平均延迟要高得多，分布的右偏也更加明显。你应当多次点击「生成新试验」，看看该系统可能表现出多大的随机变异。每次试验代表不同的一天。）

![图 12 400 次试验实验中平均来电者延迟的频率](img/ch02/img_df24ab5d92e25ba3.jpg)

**图 12 400 次试验实验中平均来电者延迟的频率。**

接下来我们要问的问题是：平均来电者延迟（前 100 位顾客）如何逐日变化？使用「Experiment」工作表，并把每一天视为一次试验，我们把"试验次数"设为 400（即 400 天）。所得平均来电者延迟的分布见图 12 的直方图。请注意，这幅直方图并不是像图 11 那样个别来电者延迟的直方图。

从图 12 可见，大约 22% 的平均延迟（400 个中的 87 个）超过 1 分钟，而大约 78%（400 个中的 313 个）不超过 1 分钟。只有大约 1.25%（400 个中的 5 个）超过 2 分钟。我们看到，第一次试验的结果（平均来电者延迟 0.73 分钟）落在最频繁的组（区间）中，但也存在平均延迟高得多的日子。这个例子表明，一天的（一次试验的）仿真很可能并不能代表系统的长期性能。

总而言之，Able 和 Baker 能否承接服务需求、能在多大程度上承接，取决于你的判据，以及顾客所经历延迟的成本或价值。我们可以猜测，增设一个服务台很可能把来电者延迟（等待时间）降至接近于零；不过，减少顾客延迟所带来的收益，必须与增设一个服务台的成本相权衡。

## 4 电子表格中的库存仿真

仿真问题中有一类重要的情形涉及库存系统。图 13 所示的简单库存系统称为 (M, N) 系统，它每隔 N 个时段审查一次库存，最大库存水平为 M。每经过 N 个时间单位，就检查一次库存水平，并下订单把库存补足到水平 M。

![图 13 概率型订购至水平库存系统](img/ch02/img_244436c22daf431b.jpg)

**图 13 概率型订购至水平库存系统。**

为简单起见，图 13 中把提前期（lead time，即下达订单与收到订单之间的时间长度）画为零；订单一经下达便立即到达。实际上，提前期可能变化极大，必须建模为一个已知的常数或一个随机变量。

需求通常无法确切知晓，因此订货量要由一个概率分布来建模。在图 13 的图形中，需求被画成在该时段上均匀分布。实际上，需求通常并不均匀，而是随时间波动。

在第一个审查周期结束时，下达一个订货量 \(Q_{1}\)，把库存水平补足到 M。每一次库存审查都重复这一做法。

请注意，在第二个周期中，库存量降到零以下，表明出现了缺货。在图 13 中，这些单位是延期交货（backorder）的；订单到达时，延期交货部分的需求先得到满足。为了避免缺货，就需要持有缓冲库存或安全库存。在某些库存系统中，库存耗尽时，一部分销售会直接流失，而不是延期交货。

库存系统有若干潜在的收入与成本来源：

1. 销售收入；

2. 售出商品的成本；

3. 持有成本，即持有库存所产生的成本；

4. 订货成本，即与下达订单相关的成本；

5. 销售流失或延期交货的成本；

6. 报废所造成的成本；

7. 报废品或其他损坏商品所得的残值收入。

持有成本可以归因于为购买商品而借入资金所支付的利息（或者等价地，由于这些资金无法用于其他投资用途而造成的损失）。持有成本还包括搬运和存储空间的成本、雇请守卫的成本以及诸如此类的其他成本。订货成本可能包括运输与物流费用，以及取决于订货量的折扣。

在试图改进库存系统时，会面临若干权衡。作为持有高库存（并承担高持有成本）的替代方案，可以更频繁地进行审查，从而更频繁地采购或补货。这会带来相应的成本，即订货成本。此外，货物脱销也可能带来成本，这种短缺会使顾客不满，并导致商誉的损失。更大的库存会降低出现短缺的概率。通过调整策略参数或决策参数，可以在各项成本之间进行权衡，以使库存系统的总成本最小。

库存系统的总成本（或总利润）就是性能度量。它会受到各种备选策略的影响。例如，在图 13 中，决策者可以控制最大库存水平 M 以及周期长度 N。

库存系统及其策略有若干参数，其中一些可控，一些不可控。可控变量就是决策变量或策略变量。对库存中的每一项物品而言，这些变量包括：

1. 最大库存水平 M

2. 审查周期 N

3. 订货量（用于补货），Q

4. 补货提前期 L

当提前期无法控制时，它可能是随机的；在其他情形下，提前期也可能被严格控制，实际上就是一个已知的常数（例如在准时制（JIT）制造系统中订购零件或部件时），此时它就成为一个策略变量。

在 (M, N) 库存系统的模型中，可能发生的事件包括对库存中物品的需求、对库存位置的审查以及随之作出的下达补货订单的决定，还有补货库存的到达。当提前期被假定为零时（如图 13 所示），后两个事件在同一时刻发生。

在下面的示例与习题中，我们仿真各种各样的库存系统与库存策略。

### 4.1 报刊销售商问题的仿真

在下面的示例中，一位报刊销售商必须决定购买多少份报纸。在此类问题中，只有单一长度的一个时段是相关的，并且只进行一次采购。单一时段结束时剩余的库存按废品卖出、丢弃或回收。大量现实世界问题具有类似的特征，包括某些备件的储备、易腐物品、时尚商品以及特殊季节性商品 [Hadley and Whitin, 1963]。

### 例 7：报刊销售商问题

一位报刊销售商以每份 33 美分购入报纸，以每份 50 美分售出。一天结束时未售出的报纸，按每份 5 美分作为废品出售。报纸可以按 10 份一捆购买。因此，报摊可以买 50 份、60 份或 70 份，等等。订货量 Q 是唯一的策略决策。与某些库存问题不同，订货量 Q 是固定的，因为剩余报纸作为废品处理，期末库存总是为零。

**表 16 按新闻日类型划分的每日报纸需求量分布**

|  | B | C | D | E | F | G | H |
|---|---|---|---|---|---|---|---|
| 4 |  | 报纸需求量分布 |  |  |  |  |  |
| 5 |  | 需求 | 需求概率 |  |  | 累积概率 |  |
| 6 | 好 | 一般 | 差 | 好 | 一般 | 差 |  |
| 7 | 40 | 0.03 | 0.10 | 0.44 | 0.03 | 0.10 | 0.44 |
| 8 | 50 | 0.05 | 0.18 | 0.22 | 0.08 | 0.28 | 0.66 |
| 9 | 60 | 0.15 | 0.40 | 0.16 | 0.23 | 0.68 | 0.82 |
| 10 | 70 | 0.20 | 0.20 | 0.12 | 0.43 | 0.88 | 0.94 |
| 11 | 80 | 0.35 | 0.08 | 0.06 | 0.78 | 0.96 | 1.00 |
| 12 | 90 | 0.15 | 0.04 | 0.00 | 0.93 | 1.00 | 1.00 |
| 13 | 100 | 0.07 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 |

新闻日有三种类型："好"、"一般"和"差"，但报刊销售商无法预测某一天会属于哪一种类型。表 16 给出了按日类型划分的报纸需求量分布。表 17 给出了新闻日类型的分布。

我们的目标是计算报摊应当采购的最优报纸份数（即订货量 Q）。做法是在「One Trial」工作表中仿真 20 天的需求，并对给定的订货量记录每天的销售利润。一次试验完成后，计算总利润（模型响应或性能度量）。为完成实验，我们使用「Experiment」工作表，对给定的订货量（例如 70 份报纸）重复某个次数（例如 400 次）。最后，我们在某个合理范围内改变订货量（40、50、60、70 等等），并依据 20 天期间的总利润比较不同的订货策略。

每日利润由下列关系给出：

$$
\text{Profit} = \binom{\text{revenue}}{\text{from   sales}} - \binom{\text{cost   of}}{\text{newspapers}} - \binom{\text{lost   profit   from}}{\text{excess   demand}} + \binom{\text{salvage   from   sale}}{\text{of   scrap   papers}}
$$

 由于销售收入为每份 50 美分，而销售商的成本为每份 33 美分，因此对每一份无法供应的需求，损失的利润为 17 美分。

**表 17 新闻日类型的分布**

|  | J | K | L |
|---|---|---|---|
| 4 | 新闻日类型 |  |  |
| 5 | 类型 | 概率 | 累积概率 |
| 6 |  |  |  |
| 7 | 好 | 0.35 | 0.35 |
| 8 | 一般 | 0.45 | 0.80 |
| 9 | 差 | 0.20 | 1.00 |

**表 18 采购 70 份报纸的仿真表**

|  | B | C | D | E | F | G | HI | I |
|---|---|---|---|---|---|---|---|---|
| 16 | 仿真表 |  |  |  |  |  |  |  |
| 17 | 天 | 新闻日类型 | 需求 | 销售收入 | 超额需求造成的利润损失 | 废品出售所得残值 | 每日成本 | 每日利润 |
| 18 |  |  |  |  |  |  |  |  |
| 19 |  |  |  |  |  |  |  |  |
| 20 | 1 | 一般 | 50 | \$25.00 | \$0.00 | \$1.00 | \$23.10 | \$2.90 |
| 21 | 2 | 一般 | 50 | \$25.00 | \$0.00 | \$1.00 | \$23.10 | \$2.90 |
| 22 | 3 | 一般 | 70 | \$35.00 | \$0.00 | \$0.00 | \$23.10 | \$11.90 |
| 23 | 4 | 好 | 80 | \$35.00 | \$1.70 | \$0.00 | \$23.10 | \$10.20 |
| 24 | 5 | 差 | 40 | \$20.00 | \$0.00 | \$1.50 | \$23.10 | -\$1.60 |
| 25 | 6 | 一般 | 50 | \$25.00 | \$0.00 | \$1.00 | \$23.10 | \$2.90 |
| 26 | 7 | 差 | 50 | \$25.00 | \$0.00 | \$1.00 | \$23.10 | \$2.90 |
| 27 | 8 | 一般 | 70 | \$35.00 | \$0.00 | \$0.00 | \$23.10 | \$11.90 |
| 28 | 9 | 好 | 40 | \$20.00 | \$0.00 | \$1.50 | \$23.10 | -\$1.60 |
| 29 | 10 | 好 | 100 | \$35.00 | \$5.10 | \$0.00 | \$23.10 | \$6.80 |
| 30 | 11 | 一般 | 70 | \$35.00 | \$0.00 | \$0.00 | \$23.10 | \$11.90 |
| 31 | 12 | 差 | 50 | \$25.00 | \$0.00 | \$1.00 | \$23.10 | \$2.90 |
| 32 | 13 | 一般 | 50 | \$25.00 | \$0.00 | \$1.00 | \$23.10 | \$2.90 |
| 33 | 14 | 差 | 40 | \$20.00 | \$0.00 | \$1.50 | \$23.10 | -\$1.60 |
| 34 | 15 | 好 | 80 | \$35.00 | \$1.70 | \$0.00 | \$23.10 | \$10.20 |
| 35 | 16 | 好 | 70 | \$35.00 | \$0.00 | \$0.00 | \$23.10 | \$11.90 |
| 36 | 17 | 好 | 80 | \$35.00 | \$1.70 | \$0.00 | \$23.10 | \$10.20 |
| 37 | 18 | 差 | 40 | \$20.00 | \$0.00 | \$1.50 | \$23.10 | -\$1.60 |
| 38 | 19 | 一般 | 70 | \$35.00 | \$0.00 | \$0.00 | \$23.10 | \$11.90 |
| 39 | 20 | 好 | 100 | \$35.00 | \$5.10 | \$0.00 | \$23.10 | \$6.80 |
| 40 |  |  |  |  |  | 总利润 = | \$114.70 |  |

这种缺货成本多少有些争议，但使问题有趣得多。废纸的残值为每份 5 美分。

采购 70 份报纸这一决策的仿真表（取自解答工作簿 "Example2.7NewsDealer.xls"）见表 18。在第 1 天，需求为 50 份报纸，因此有 70 份可供出售，其中 20 份作为废品处理。卖出 50 份报纸所得收入为 \$25.00，废品残值为 \$1.00，因此第一天的利润计算如下：

$$
\text{First - day   Profit} = 25. 00 - 23. 10 - 0 + 1. 00 = 2. 90
$$

 其余各天的利润用同样的方法即可容易地算出。20 天期间的利润是各日利润之和，为 \$114.70，如表 18 所示。一般而言，由于某一天的结果与之前各天无关，这类库存问题在电子表格中要比排队问题容易解决得多。

使用 "Example2.7NewsDealer.xls" 中的「Experiment」工作表，我们针对每次订购 70 份报纸的订货策略，设置了包含 400 次试验（或重复）的实验。图 14 给出了这 400 次试验的结果，每次试验都是一次 20 天的仿真运行。20 天总利润的平均值（在 400 次试验上取平均）为 \$135.49。20 天利润的最小值为 \$86.60，最大值为 \$198.00。图 14 表明，在 400 次试验中，只有 38 次的 20 天总利润超过 \$160。（请记住，组（区间）下方的数字是该组（区间）的上限。）

![图 14 报刊销售商每天订购 70 份报纸时 20 天总利润的频率](img/ch02/img_46fe946121296423.jpg)

**图 14 报刊销售商每天订购 70 份报纸时 20 天总利润的频率。**

表 18 中所示的那一次试验利润为 \$114.70，处于图 14 中分布的偏低一侧。这再一次说明，一次 20 天试验的结果可能并不代表长期性能，同时也表明有必要仿真若干次试验，并在得出结论之前进行仔细的统计分析。

在「One Trial」工作表上，注意多次点击「生成新试验」按钮时每日利润的变化。无论在标题为"每日利润频率"的直方图（展示 20 天中每一天发生了什么）中，还是在这 20 天的总利润上，结果的变化都相当大。图 15 给出了前四次试验的每日利润直方图，其形状的变化比 20 天总利润直方图的形状变化要大得多。

习题 28 要求你开展一项实验，以确定报刊销售商的最优订货报纸数量。用仿真解决该问题，需要设定每天购买某个数量报纸的策略，并针对该策略运行大量次数的实验。把策略（所购买的报纸数量）改为某个合理范围内的其他取值（比如说 40、50、60 等等），并重复仿真，直到找到最佳取值。要有把握地估计某个仿真的差异很可能是真实存在、而非随机偶然造成的，还需要额外的统计方法。

![图 15 报刊销售商问题每日利润的四个典型直方图](img/ch02/img_7749fc7916c05b27.jpg)

![图 15 报刊销售商问题每日利润的四个典型直方图](img/ch02/img_9130786e73f413d9.jpg)

![图 15 报刊销售商问题每日利润的四个典型直方图](img/ch02/img_843594f1c38bb669.jpg)

![图 15 报刊销售商问题每日利润的四个典型直方图](img/ch02/img_91432742f6a3cd49.jpg)

**图 15 报刊销售商问题每日利润的四个典型直方图。**

### 4.2 (M, N) 库存策略的仿真

下一个示例是经典的订购至水平 (M, N) 库存策略，它具有固定的审查周期（N）、最大库存水平（M）、随机需求和随机提前期。

### 例 8：冰箱库存问题

某销售冰箱的公司通过每隔固定天数（例如 N 天）进行一次审查、并就订购多少台以补充店铺库存作出决策，来维持库存。当前策略是订购至某一水平（订购至水平——设为 M），使用下列关系：

$$
\text{Order   quantity} = (\text{Order - up - to   level}) - (\text{Ending   inventory}) + (\text{Shortage   quantity})
$$

 例如，若订购至水平 M 为 11 台冰箱，审查周期 N 为 5 天，而第 5 天结束时的库存为 3 台冰箱，那么在该周期的第五天，将向供应商订购 8 台冰箱。另一方面，若第五天出现了 2 台冰箱的缺货，那么就会订购 13 台冰箱；最先收到的两台将提供给那些已经下单并愿意等待的顾客（所谓的延期交货）。 （与之相对，销售流失的情形是指在无库存可用时顾客的需求直接流失。）

每日需求量（即顾客每天购买的冰箱台数）是随机分布的，如表 19 所示。公司下达补货订单之后，提前期为 1 至 3 天，随机分布，如表 20 所示。假定订单在一天结束时下达。提前期为零时，来自供应商的订单在次日清晨到达，冰箱在第二天即可供应。提前期为一天时（表 20 中的最小可能值），来自供应商的订单在两天后的清晨到达，并在当天即可供应。

**表 19 每日需求量的分布**

|  | B | C | D |
|---|---|---|---|
| 4 | 每日需求量的分布 |  |  |
| 5 | 需求量 | 概率 | 累积概率 |
| 6 |  |  |  |
| 7 | 0 | 0.10 | 0.10 |
| 8 | 1 | 0.25 | 0.35 |
| 9 | 2 | 0.35 | 0.70 |
| 10 | 3 | 0.21 | 0.91 |
| 11 | 4 | 0.09 | 1.00 |

**表 20 提前期的分布**

|  | F | G | H |
|---|---|---|---|
| 4 | 提前期的分布 |  |  |
| 5 | 提前期（天） | 概率 | 累积概率 |
| 6 |  |  |  |
| 7 | 1 | 0.60 | 0.60 |
| 8 | 2 | 0.30 | 0.90 |
| 9 | 3 | 0.10 | 1.00 |

"订购至 11 台冰箱"这一策略的仿真表见表 21，取自电子表格解答 "Example2.8RefrigInventory.xls"。仿真以如下的假定初始条件开始：期初库存为 3 台冰箱，且此前有一笔 8 台冰箱的补货订单将在 2 天后到达。（如需要，这些值可以在电子表格解答中更改。）

跟踪仿真表中若干选取的日子，可以看出仿真的逻辑。8 台冰箱的订单在第一个周期第三天的清晨可以到货，使库存水平从零台冰箱上升到 8 台冰箱。第一个周期余下时间的需求，使第五天结束时的库存水平降至 3 台冰箱。因此下达了一个 8 台冰箱的订单；其提前期恰为 2 天。这笔订单在第 2 个周期的第 7 天清晨计入库存。

请注意，在第 10 天出现了缺货，形成了一笔延期交货，并下达了一笔 13 台冰箱的订单。由于提前期为 2 天，在订单到达之前还会发生额外的缺货，因此当这笔 13 台冰箱的订单到达时，延期交货加上当期需求使库存降至 6 台。

从 25 天（相当于 5 个库存审查周期）的仿真来看，期末库存的平均值约为 2.76（69/25）台。在 25 天中，有 5 天存在缺货状况。

**表 21 订购至 11 台冰箱的库存策略仿真表**

|  | B | C | D | E | F | G | H | I | J | K |
|---|---|---|---|---|---|---|---|---|---|---|
| 13 | 仿真表 |  |  |  |  |  |  |  |  |  |
| 14 | 步骤 | 状态 | 状态 | 状态 | 输入 | 状态 | 状态 | 状态 | 活动 | 状态 |
| 15 | 时钟 | 周期 | 周期内天数 | 期初库存 | 需求量 | 期末库存 | 缺货数量 | 在途订单（数量） | 提前期（天） | 距订单到达的天数 |
| 16 | 天 |  |  |  |  |  |  |  |  |  |
| 17 |  |  |  |  |  |  |  |  |  |  |
| 18 | 0 | 0 | 5 | - | - | 3 | 0 | 8 | 2 | 2 |
| 19 | 1 | 1 | 1 | 3 | 2 | 1 | 0 | 8 |  | 1 |
| 20 | 2 | 1 | 2 | 1 | 1 | 0 | 0 |  |  |  |
| 21 | 3 | 1 | 3 | 8 | 2 | 6 | 0 |  |  |  |
| 22 | 4 | 1 | 4 | 6 | 1 | 5 | 0 |  |  |  |
| 23 | 5 | 1 | 5 | 5 | 2 | 3 | 0 | 8 | 1 | 1 |
| 24 | 6 | 2 | 1 | 3 | 3 | 0 | 0 |  |  |  |
| 25 | 7 | 2 | 2 | 8 | 2 | 6 | 0 |  |  |  |
| 26 | 8 | 2 | 3 | 6 | 3 | 3 | 0 |  |  |  |
| 27 | 9 | 2 | 4 | 3 | 2 | 1 | 0 |  |  |  |
| 28 | 10 | 2 | 5 | 1 | 3 | 0 | 2 | 13 | 2 | 2 |
| 29 | 11 | 3 | 1 | 0 | 1 | 0 | 3 | 13 |  | 1 |
| 30 | 12 | 3 | 2 | 0 | 2 | 0 | 5 |  |  |  |
| 31 | 13 | 3 | 3 | 13 | 2 | 6 | 0 |  |  |  |
| 32 | 14 | 3 | 4 | 6 | 3 | 3 | 0 |  |  |  |
| 33 | 15 | 3 | 5 | 3 | 1 | 2 | 0 | 9 | 1 | 1 |
| 34 | 16 | 4 | 1 | 2 | 0 | 2 | 0 |  |  |  |
| 35 | 17 | 4 | 2 | 11 | 4 | 7 | 0 |  |  |  |
| 36 | 18 | 4 | 3 | 7 | 2 | 5 | 0 |  |  |  |
| 37 | 19 | 4 | 4 | 5 | 3 | 2 | 0 |  |  |  |
| 38 | 20 | 4 | 5 | 2 | 3 | 0 | 1 | 12 | 1 | 1 |
| 39 | 21 | 5 | 1 | 0 | 2 | 0 | 3 |  |  |  |
| 40 | 22 | 5 | 2 | 12 | 1 | 8 | 0 |  |  |  |
| 41 | 23 | 5 | 3 | 8 | 4 | 4 | 0 |  |  |  |
| 42 | 24 | 5 | 4 | 4 | 1 | 3 | 0 |  |  |  |
| 43 | 25 | 5 | 5 | 3 | 1 | 2 | 0 | 9 | 1 | 1 |
| 44 | 合计 |  |  |  |  | 69 | 14 |  |  |  |
| 45 | 平均值 |  |  |  | 2.04 | 2.76 | 0.56 |  |  |  |

**冰箱每日期末库存平均值的频率分布**

![图 16 25 天期末库存平均值的直方图](img/ch02/img_4c2b1a8c69009d92.jpg)

**图 16 25 天期末库存平均值的直方图。**

在例 8 中，任何时刻都只能有一笔供应商订单尚未到货；但也存在这样的情形：提前期极长，以致迄今为止所展示的关系需要作如下修改：

$$
\text{Order   quantity} = (\text{Order - up - to   level}) - (\text{Ending   inventory}) - (\text{On   order}) + (\text{Shortage   quantity})
$$

 计入在途数量，可以保证不会出现重复订货。

在「One Trial」工作表上，反复点击「生成新试验」将用不同的随机数运行下一次试验（或重复），这表明对于 25 天的仿真，每日期末库存的结果变化很大。（请注意随着新试验的运行，直方图如何变化。）相比之下，在「Experiment」工作表中把试验次数设为 100 并通过点击「运行实验」按钮反复运行实验，则在每日期末库存的平均值上产生的变异性要小得多（平均值是在一次 25 天的试验上取的）。25 天的平均值通常落在 1 到 5 台冰箱的范围内，如图 16 所示的直方图。变异性上的这一差异说明了一个众所周知的事实：平均数的波动小于单个观测值。在具体分析中该用哪一个，取决于所提出的目标与问题。

"Example2.8RefrigInventory.xls" 中的电子表格解答允许你改变订货策略（即 M 和 N 的取值）。它还允许你对每日需求量和提前期的统计分布作一些修改，但须满足某些限制：提前期必须为一天或更多，审查周期必须大于最长的提前期，可能还有其他限制；否则，实现该逻辑的公式可能不正确。

## 5 仿真的其他示例

本节包含四个杂项示例：一个可靠性问题、一次轰炸任务、提前期需求分布的估计，以及一个项目作业网络。

### 5.1 可靠性问题的仿真

下一个示例比较铣床上更换轴承的两种更换策略。该示例是一个蒙特卡罗仿真，而不是动态的事件驱动模型，因为其中不含事件和时钟时间。每个轴承的寿命都随机生成 15 次，并据此计算成本。

### 例 9：更换铣床上的轴承

一台铣机有三个不同的轴承，会在服役中失效。每个轴承寿命的分布都相同，如表 22 所示。轴承失效时，铣机停机，召唤一名机械师，由他为她装上一个新轴承（每只轴承 \$32）。机械师到达的延迟时间是随机的，其分布见表 23。估计铣机停机的成本为每分钟 \$10。机械师在现场的直接成本为每小时 \$30。机械师更换一个轴承需要 20 分钟，更换两个轴承需要 30 分钟，更换三个轴承需要 40 分钟。工程人员提出了一项新策略：每当一个轴承失效时，就更换全部三个轴承。管理层需要对这一提案作出评价，所用的性能度量为每 10,000 轴承小时的总成本。

**表 22 轴承寿命的分布**

|  | B | C | D |
|---|---|---|---|
| 3 | 轴承寿命分布 |  |  |
| 4 | 轴承寿命 | 概率 | 累积概率 |
| 5 |  |  |  |
| 6 | 1000 | 0.10 | 0.1 |
| 7 | 1100 | 0.13 | 0.23 |
| 8 | 1200 | 0.25 | 0.48 |
| 9 | 1300 | 0.13 | 0.61 |
| 10 | 1400 | 0.09 | 0.70 |
| 11 | 1500 | 0.12 | 0.82 |
| 12 | 1600 | 0.02 | 0.84 |
| 13 | 1700 | 0.06 | 0.90 |
| 14 | 1800 | 0.05 | 0.95 |
| 15 | 1900 | 0.05 | 1.00 |

**表 23 直至机械师到达的延迟时间的分布**

|  | F | G | H |
|---|---|---|---|
| 4 | 延迟时间的分布 |  |  |
| 5 | 延迟时间 | 概率 | 累积概率 |
| 6 |  |  |  |
| 7 | 5 | 0.60 | 0.60 |
| 8 | 10 | 0.30 | 0.90 |
| 9 | 15 | 0.10 | 1.00 |

**表 24 现行方法下的轴承更换**

|  | B | C | D | E | F | G | H |
|---|---|---|---|---|---|---|---|
| 17 | 仿真表 |  |  |  |  |  |  |
| 18 | 步骤 | 轴承 1 | 轴承 2 | 轴承 3 |  |  |  |
| 19 | 寿命（小时） | 延迟（分钟） | 寿命（小时） | 延迟（分钟） | 寿命（小时） | 延迟（分钟） |  |
| 20 |  |  |  |  |  |  |  |
| 21 | 1 | 1000 | 5 | 1700 | 10 | 1300 | 10 |
| 22 | 2 | 1200 | 5 | 1100 | 5 | 1100 | 5 |
| 23 | 3 | 1200 | 10 | 1000 | 10 | 1300 | 10 |
| 24 | 4 | 1500 | 5 | 1000 | 10 | 1100 | 15 |
| 25 | 5 | 1700 | 5 | 1900 | 15 | 1200 | 5 |
| 26 | 6 | 1200 | 5 | 1200 | 10 | 1500 | 10 |
| 27 | 7 | 1300 | 5 | 1500 | 5 | 1100 | 10 |
| 28 | 8 | 1700 | 5 | 1700 | 5 | 1400 | 15 |
| 29 | 9 | 1000 | 5 | 1300 | 5 | 1800 | 15 |
| 30 | 10 | 1800 | 10 | 1300 | 5 | 1200 | 5 |
| 31 | 11 | 1200 | 5 | 1100 | 5 | 1500 | 5 |
| 32 | 12 | 1100 | 5 | 1800 | 5 | 1100 | 10 |
| 33 | 13 | 1300 | 10 | 1200 | 5 | 1700 | 10 |
| 34 | 14 | 1300 | 10 | 1100 | 5 | 1300 | 10 |
| 35 | 15 | 1100 | 5 | 1300 | 5 | 1300 | 10 |
| 36 | 合计 | 19,600 | 95 | 20,200 | 105 | 19,900 | 145 |

现行策略的仿真表见表 24，取自电子表格 "Example2.9CurrentBearings.xls"，它表示在"只更换那个失效的单个轴承"的现行策略下，15 次轴承更换的仿真。请注意，存在多个轴承同时失效的情形。这在实践中不太可能发生，其原因是采用了粗略的简化假设——假定轴承寿命是 100 小时的整数倍（这样我们就能采用一个简单的离散分布，而不是更贴近现实的分布）。这里我们假定这些时刻绝不会完全相同，因此任何一次故障中至多只更换一个轴承。

由这次仿真试验（或重复）出发，现行系统的成本估计如下：

$$
\begin{array}{r l} \text{Cost of bearing} & = 45 \text{bearings} \times 32 / \text{bearing} \\ & = 1, 440 \\ \text{Cost of delay time} & = (95 + 105 + 145) \text{minutes} \times 10 / \text{minute} \\ & = 3, 450 \\ \text{Cost of downtime during repair} & = 45 \text{bearings} \times 20 \text{minutes/bearing} \times 10 / \text{minute} \\ & = 9, 000 \\ \text{Cost of mechanics} & = 45 \text{bearings} \times 20 \text{minutes/bearing} \times 30 / 60 \text{minutes} \\ & = 450 \\ \text{Total cost} & = 1, 440 + 3, 450 + 9, 000 + 450 \\ & = 14, 340 \end{array}
$$

 全部 45 个轴承的总寿命为 \(( 19 , 600 + 20 , 200 + 19 , 900 ) = 59 , 700\) 小时。因此，每 10,000 轴承小时的总成本为 (\$14,340/5.97) \$2,402。

所提方法的仿真表见表 25，取自电子表格 "Example2.9ProposedBearings.xls"。对于第一组轴承，最早的失效发生在 1,100 小时。

**表 25 所提方法下的轴承更换**

|  | B | C | D | E | F | G |
|---|---|---|---|---|---|---|
| 17 | 仿真表 |  |  |  |  |  |
| 18 | 步骤 | 轴承 1 | 轴承 2 | 轴承 3 | 首次失效（小时） |  |
| 19 | 寿命（小时） | 寿命（小时） | 寿命（小时） | 延迟（分钟） |  |  |
| 20 |  |  |  |  |  |  |
| 21 | 1 | 1300 | 1100 | 1300 | 1,100 | 10 |
| 22 | 2 | 1100 | 1200 | 1500 | 1,100 | 5 |
| 23 | 3 | 1100 | 1400 | 1800 | 1,100 | 15 |
| 24 | 4 | 1200 | 1500 | 1100 | 1,100 | 10 |
| 25 | 5 | 1700 | 1300 | 1300 | 1,300 | 5 |
| 26 | 6 | 1700 | 1100 | 1000 | 1,000 | 10 |
| 27 | 7 | 1000 | 1900 | 1200 | 1,000 | 10 |
| 28 | 8 | 1500 | 1700 | 1300 | 1,300 | 15 |
| 29 | 9 | 1300 | 1100 | 1800 | 1,100 | 5 |
| 30 | 10 | 1200 | 1100 | 1300 | 1,100 | 5 |
| 31 | 11 | 1000 | 1200 | 1200 | 1,000 | 15 |
| 32 | 12 | 1500 | 1700 | 1200 | 1,200 | 10 |
| 33 | 13 | 1300 | 1700 | 1000 | 1,000 | 10 |
| 34 | 14 | 1800 | 1200 | 1100 | 1,100 | 10 |
| 35 | 15 | 1300 | 1300 | 1100 | 1,100 | 10 |
| 36 | 合计 |  |  |  | 16600 | 145 |

在那时，三个轴承全都被更换了，尽管其余轴承还剩下更多的寿命。例如，轴承 1 本来还可以再使用 200 小时。

由该模型这次仿真试验（或重复）出发，所提策略的成本估计如下：

$$
\begin{array}{r l} \text{Cost   of   bearings} & = 45 \text{bearings} \times 32 / \text{bearing} \\ & = 1, 440 \end{array}
$$


$$
\begin{array}{r l} \text{Cost   of   delay   time} & = 145 \text{minutes} \times 10 / \text{minute} \\ & = 1, 450 \end{array}
$$


$$
\begin{array}{r l} \text{Cost   of   downtime   during   repairs} & = 15 \text{sets} \times 40 \text{minutes / set} \times 10 / \text{minute} \\ & = 6, 000 \end{array}
$$


$$
\begin{array}{r l} \text{Cost   of   mechanics} & = 15 \text{sets} \times 40 \text{minutes / set} \times 30 / 60 \text{minutes} \\ & = 300 \end{array}
$$


$$
\begin{array}{r l} \text{Total   cost} & = 1, 440 + 1, 450 + 6, 000 + 300 \\ & = 9, 190 \end{array}
$$

 轴承的总寿命为 \(( 16 , 600 \times 3 ) = 49 , 800\) 小时。因此，每 10,000 轴承小时的总成本为 (\$9,190/4.98) \$1,845。新策略每 10,000 小时的轴承寿命可节省 \$557。

例 9 的电子表格解答分别是：现行策略见 "Example2.9CurrentBearings.xls"，所提策略见 "Example2.9ProposedBearings.xls"。在每个电子表格中，用户都可以更改若干输入：轴承寿命的分布、直至机械师到达的延迟时间的分布、成本参数，以及每个轴承或每组轴承的维修时间。

### 5.2 命中目标的仿真

下一个示例仿真一次轰炸任务，其中炸弹实际落点相对于目标的偏差量由一个正态分布定义，随机取值。该电子表格解答提供了一种生成这些正态随机变量的方法；手工解答则使用附录中表 A.2 的标准正态随机数。

就我们的目的而言，请回顾一下：正态分布（钟形曲线）关于其均值对称，其可能的取值范围由标准差决定；它在下一个示例中的用途，是正态分布用来表示误差的经典用法，本例中即相对于目标的偏差。这个示例是一个纯粹的蒙特卡罗抽样实验，因为其中既不含事件也不含时钟时间；因此，它不是动态的离散事件模型。

### 例 10：一次轰炸任务

考虑一架轰炸机试图摧毁一座弹药库，如图 17 所示。（这架轰炸机携带的是常规武器而非激光制导武器。）如果炸弹落在目标内部任何位置，就算命中；否则就是未命中。（请注意，当炸弹在视觉上似乎触及某条边界线时，它可能命中也可能没有命中目标；模型利用 (X, Y) 坐标以及弹药库分段线性边界的方程，在数学上判定是否发生了命中。）

![图 17 带有 10 个炸弹落点的弹药库](img/ch02/img_0385f243f1b7ea5b.jpg)

**图 17 带有 10 个炸弹落点的弹药库。**

轰炸机沿水平方向飞行，携带 10 枚炸弹。瞄准点为 (0,0)。假定实际命中点围绕瞄准点呈正态分布，在飞行方向上标准差为 400 米，在垂直方向上标准差为 200 米。问题是仿真这一行动并估计命中目标的炸弹数量。

请回顾，标准正态变量 \(Z ,\) 的均值为 0、标准差为 1，其分布为 

$$
Z = \frac {X - \mu}{\sigma}
$$

 其中 \(X\) 是正态随机变量，\(\mu\) 是 X 的分布的均值，σ 是 X 的标准差。于是，在均值为零、且 \(\sigma_{X} = 400\) 与 \(\sigma_{Y} = 200\) 的情况下，有 

$$
\begin{array}{c} {X = 400 Z_{i}} \\ {Y = 200 Z_{j}} \end{array}
$$

 其中 (X, Y) 是模拟出的炸弹命中坐标。

**表 26 模拟的一次轰炸**

|  | F | G | H | I | J | K | L | M | N | O | P | Q |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3 |  | 各枚炸弹的命中点（1–10） |  |  |  |  |  |  |  |  |  |  |
| 4 | 炸弹 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 命中数 |
| 5 | X | 891.8 | -1257.3 | 429.6 | 24.5 | -321.0 | -77.3 | 131.4 | 304.5 | -442.7 | -403.9 |  |
| 6 | Y | -400.7 | -159.4 | 25.3 | 243.6 | 146.5 | 260.7 | -228.3 | -62.0 | 49.7 | 51.0 |  |
| 7 | 命中？ | 未命中 | 未命中 | 命中 | 命中 | 命中 | 命中 | 命中 | 命中 | 未命中 | 未命中 | 6 |

下标 i 和 j 是为了表明 x 方向和 y 方向的 Z 值应当不同。这些 Z 值是什么，又能在哪里找到？Z 值是均值为零、标准差为一的随机正态数。它们可以从均匀分布的随机数生成。就我们的目的而言，手工仿真可以使用附录中的表 A.2，电子表格模型则可以使用第 1.2 节列出的 Excel 公式。

例如，表 A.2 中最先的两个随机正态数是 0.23 和 -0.17，由此得到第一枚炸弹的 (x, y) 坐标为 (92.0, 34.0)，这是把随机正态数分别乘以 400 和 200 这两个标准差得到的，如上式所示。

表 26 给出了模拟一次轰炸运行的结果。在我们的电子表格解答 "Example2.10Target.xls" 中，我们使用第 1.2 节所述的工作表函数 NORMSINV(Rnd01()) 来生成具有所需标准差的正态变量，如下所示：

$$
= E 5 * N O R M S I N V (R n d 01 ())
$$

 用于 x 坐标（其中单元格 E5 存放 x 方向的标准差）。y 坐标可以用一个类似的公式处理。

使用 "Example2.10Target.xls" 中的电子表格解答，按假定的形状和标准差，我们运行了包含 400 次试验的实验（每次试验为 10 枚炸弹）。图 18 的直方图给出命中数量的频数，即有多少次轰炸运行取得了给定的命中数（可能的命中数从 0 到 10）。

从电子表格模型「Experiment」工作表的"多次试验汇总"中可以看到，实际结果从 2 次命中到 10 次命中，平均为 6.77 次命中。（显然，0 次命中或 1 次命中都是可能的。由于我们做了 400 次重复，实验仍表明：在具有所假定特征的轰炸任务中，1 次命中或 0 次命中极不可能发生。）如果我们只仿真一次任务（即一次试验），可能得到极具误导性的结果，而图 18 提供了有用的描述性信息。例如，在 17%（400 次中的 68 次）的轰炸运行中，命中数为 5 次或更少。在约 60% 的运行中（400 次中的 239 次），命中数为 7 次或更多。

"Example2.10Target.xls" 中的电子表格模型允许你通过编辑目标各条边两端点的 (X, Y) 坐标来改变目标的形状，只要保持形状为凸形即可。你也可以改变标准差，标准差越小则越有可能更接近目标，越大则离目标越远的可能性越大。在该电子表格解答中，模拟落点为 (X, Y) 的炸弹是否命中目标，是用初等代数导出的一组线性不等式来判定的；细节在电子表格中有说明。

**400 次轰炸任务中命中数量的频率分布**

![图 18 轰炸任务命中数量的直方图](img/ch02/img_22937124fefe9904.jpg)

**图 18 轰炸任务命中数量的直方图。**

### 5.3 提前期需求分布的估计

下一个示例的目的，是确定提前期需求（lead-time demand）的分布，即在整个提前期内所收到的总需求。提前期是从下达补货订单到收到该订单之间的时间。在这段期间内，需求继续从顾客那里到来。这一期间内顾客的总需求就是提前期需求。

这个模型是蒙特卡罗仿真的一个例子；它不是动态的事件驱动仿真。本质上，其中没有事件也没有时钟时间；它是一个纯粹的抽样实验，其目的是估计一个复杂的统计分布。提前期需求的分布常用于设计更好的库存策略，以尽可能低的成本减少延期交货和销售流失。

我们假定提前期和需求都是随机变量。因此，提前期需求是这样一个随机变量：定义为提前期内各期需求之和，即 \(\textstyle \sum_{i = 0} ^ {T} D_{i}\)，其中 T 是提前期，\(D_{i}\) 是第 i 个时间周期内的需求。提前期需求的分布是这样求得的：仿真许多个提前期周期，对提前期中的每个周期生成随机需求并求和，然后根据结果构建直方图。

### 例 11：提前期需求

某公司批发销售新闻纸卷。每日需求量由下列概率分布给出：

| 每日需求（卷） | 3 | 4 | 5 | 6 |
|---|---|---|---|---|
| 概率 | 0.20 | 0.35 | 0.30 | 0.15 |

**表 27 提前期需求的仿真表**

|  | B | C | D | E | F | G |
|---|---|---|---|---|---|---|
| 13 | 仿真表 |  |  |  |  |  |
| 14 | 周期 | 提前期 | 第 1 天需求 | 第 2 天需求 | 第 3 天需求 | 提前期需求 |
| 15 |  |  |  |  |  |  |
| 16 |  |  |  |  |  |  |
| 17 | 1 | 2 | 3 | 5 |  | 8 |
| 18 | 2 | 1 | 4 |  |  | 4 |
| 19 | 3 | 2 | 3 | 5 |  | 8 |
| 20 | 4 | 3 | 4 | 5 | 4 | 13 |
| 21 | 5 | 3 | 4 | 5 | 6 | 15 |
| 22 | 6 | 1 | 3 |  |  | 3 |
| 23 | 7 | 1 | 6 |  |  | 6 |
| 24 | 8 | 2 | 5 | 6 |  | 11 |
| 25 | 9 | 2 | 5 | 3 |  | 8 |
| 26 | 10 | 3 | 4 | 3 | 5 | 12 |
| 27 | 11 | 1 | 4 |  |  | 4 |
| 28 | 12 | 1 | 5 |  |  | 5 |
| 29 | 13 | 3 | 4 | 4 |  | 14 |
| 30 | 14 | 1 | 6 |  |  | 6 |
| 31 | 15 | 1 | 3 |  |  | 3 |
| 32 | 16 | 2 | 5 | 3 |  | 8 |
| 33 | 17 | 2 | 4 | 6 |  | 10 |
| 34 | 18 | 3 | 4 | 5 | 5 | 14 |
| 35 | 19 | 3 | 4 | 3 | 6 | 13 |
| 36 | 20 | 2 | 5 | 6 |  | 11 |

提前期（即从下达订单到公司从供应商收到订单之间的天数）是一个随机变量，由下列分布给出：

| 提前期（天） | 1 | 2 | 3 |
|---|---|---|---|
| 概率 | 0.36 | 0.42 | 0.22 |

表 27 给出了仿真表，取自 "Example2.11LeadTimeDemand.xls"。这一仿真有一处不同于先前各例的特征：一个周期中所生成的随机需求个数本身是随机的。若生成的提前期为 2 天，就要生成两个随机需求并求和；若提前期为 3 天，则必须生成三个随机需求并求和。

使用电子表格解答，你可以生成许多个周期，把每个周期内的需求求和得到提前期需求，并通过直方图看到由此得到的提前期需求分布，如图 19 所示。该结果是通过在电子表格解答 "Example2.11LeadTimeDemand.xls" 中使用「One Trial」工作表点击「重置种子并运行」来仿真 20 个提前期周期而得到的。多次点击「生成新试验」按钮，你可以看到仅仅 20 个周期之内，提前期需求会有多大的变异。

**20 个周期内提前期需求的频率分布**

![图 19 提前期需求的频率](img/ch02/img_d25cd4e240b00f6a.jpg)

**图 19 提前期需求的频率。**

把这 20 个周期的结果与电子表格模型「Experiment」工作表上的结果作比较，后者给出的是 200 次试验的结果，每次试验对应一个提前期需求。这一比较表明，由于样本量小，图 19 的直方图可能并不代表提前期需求分布的形状。这说明样本量（试验次数）必须足够大，才能对所期望的响应变量作出准确的估计。

### 5.4 作业网络的仿真

假设你有一个项目，需要完成若干项活动。每项活动需要一定的时间，或为常数，或为随机。有些活动必须依次进行，也就是说要按规定顺序；另一些活动则可以并行进行。这个项目可以用一张活动网络来表示，如图 20 所示。弧表示活动，节点表示活动的开始或结束。

在图 20 的网络中，有三条路径，每条都从标有"开始"（Start）的节点通往标有"结束"（Finish）的节点。每条路径代表一系列必须依次完成的活动。不同路径上的活动可以并行进行。每项活动的时间标注在弧上，例如 U(2, 4) 表示一个在 2 到 4 分钟之间均匀分布的随机变量。完成一条路径上所有活动所需的时间，是该路径上各活动时间之和。要完成整个项目，所有活动都必须完成；因此，项目完成时间是所有路径完成时间的最大值。

与先前所有示例都不同的是，这里我们使用的是连续均匀分布；也就是说，下限与上限之间的所有取值都是允许的，而不仅仅是整数值。随加载项文件 "BCNNvba.xla" 一起提供、并且任何电子表格模型都可以使用的 VBA 函数 Uniform()，用一行 VBA 代码实现了均匀生成，公式很简单，它把在 0 到 1 上均匀分布的 VBA 随机数发生器 Rnd() 变换到任何指定范围，如下：

```vba
' In the Uniform function, low is the lower limit,
'    high is the upper limit.
' Type ``Double'' is a real value (double precision).
Function Uniform(low As Double, high As Double) As Double
Uniform = low + (high - low) * Rnd()
End Function
```

我们希望这个简单的公式多少有些直观上的吸引力。基本上，由于 Rnd() 在单位区间 (0, 1) 上均匀分布，因此函数 Uniform() 在所指定区间 (low, high) 上也均匀分布。

Pritsker [1995] 是把仿真技术应用于项目作业网络的先驱。虽然例 12 是一个蒙特卡罗仿真，但有些作业网络模型由于复杂性更高，可能要用离散事件仿真软件包或者专门为此类模型设计的软件包来实现。

### 例 12：项目仿真

假设三位朋友想为周末来访的客人做培根、鸡蛋和烤面包作为早餐。他们把工作分开，每人准备一样东西。图 20 表示活动的先后顺序，每条路径代表一个人的任务，如下：

```txt
Top path: Start → A Crack eggs
A → B Scramble eggs
B → Finish Cook eggs
Middle path: Start → C Make toast
C → Finish Butter toast
Bottom path: Start → Finish Fry bacon
```

完成每项活动所需的时间有一定的不确定性，由图 20 中相应弧上所示的均匀分布表示。这几位朋友希望估计总准备时间的分布，以便能告知客人何时到早餐桌旁就座。有了这个分布，而不只是总准备时间的均值或平均值，他们就能估计在指定时间内准备好早餐的概率。该模型位于电子表格 "Example2.12Project.xls" 中。

![图 20 例 12 的活动网络](img/ch02/img_12468abd95c19d87.jpg)

**图 20 例 12 的活动网络。**

由于连续均匀分布的均值是下限与上限的中点，可以看到图 20 项目网络中的 3 条路径每一条的均值都是 9 分钟。由各均匀分布的下限可以看出，每条路径的最短时间是 6 分钟；同样地，最长时间是 12 分钟。这些时间对预测项目总完成时间有任何帮助吗？

就准备早餐而言，项目完成时间就是所有路径上的最大时间。使用电子表格解答 "Example2.12Project.xls"，我们可以对任务时间进行大量次数的抽样，从而了解项目完成时间可能的范围与分布。「One Trial」工作表给出了网络仿真一次试验（或重复）的解答。实现它的逻辑和公式比先前许多示例都要简单得多，我们把这些留给学生。

借助「Experiment」工作表，我们可以运行许多次试验，估计项目完成时间的可能范围以及其他统计量。我们把"试验次数"设为 400，并通过点击一次「重置种子并运行」按钮来运行实验。在 400 次试验并使用默认种子的情况下，结果如下：

样本均值：10.21 分钟

最小值：6.65 分钟

最大值：11.99 分钟

所谓的关键路径（critical path），是在某一次特定试验中最终成为最长的那条路径；也就是说，该路径的时间就是项目完成时间。在 400 次试验中，每次使用不同的随机数，每条路径都有一定的概率成为关键路径。对于这 400 次试验中的每一次，实验都确定哪条路径是关键路径，结果如下：

上路径：25.50% 的试验

中间路径：32.25%

下路径：42.25%

我们由此得出结论：培根最后才做好的概率约为 42.25%。为什么各条路径不是各占约 1/3 的时间？此外，42.25% 这个数有多准确？你可以多次点击「运行实验」按钮，以了解这一估计通常会有多大的变动；每点击一次就会变动几个百分点。这种随机变动的幅度取决于实验中的试验次数（本例中为 400）。

**项目完成时间的频率分布**

![图 21 基于 400 次试验的项目完成时间直方图](img/ch02/img_b33ddf57471d82b1.jpg)

**图 21 基于 400 次试验的项目完成时间直方图。**

项目完成时间被绘制在频率图（直方图）中，如图 21 所示。从该图我们得出以下推论：

早餐在 9 分钟或更短时间内准备好的概率约为 10.5%（400 次中的 42 次）。

耗时在 11 到 12 分钟之间的概率约为 24.25%（400 次中的 97 次）。

请注意，标注为"12.0"的组（区间）实际上是 11.0 到 12.0 分钟这一区间，因为组（区间）的标注始终是其上限。

回顾 3 条路径各自的平均时间为 9 分钟，直方图告诉我们耗时超过 9 分钟的概率约为 89.5%。换句话说，各条路径的均值几乎不能说明项目完成时间的分布——毕竟项目完成时间是由最慢的那一串任务决定的。

## 6 小结

本章的目的，是通过示例引入若干关键的仿真概念，并说明若干一般性的应用领域。

我们以仿真表——为每个问题专门定制——作为组织主线。仿真表中的各列代表活动时间或其他随机输入（例如需求）、事件发生时的时钟时间、系统状态，以及模型输出。各行代表一个仿真步骤，例如一位顾客或一天的需求。每个单元格中的公式实现仿真逻辑，或者计算模型输出。这些示例说明，需要弄清楚输入数据的特征、根据输入假设生成随机变量、定义并计算模型响应（性能度量），以及分析所得到的响应。

模型输入由策略参数和决策参数的常数值来定义，或者由服务时间、到达间隔时间、提前期和需求等变量的简单统计分布来定义。模型响应（也称系统性能度量）与仿真表分开计算，但所用的正是模型输出和其他数据，其依据是一次试验或一次仿真运行的结果。这些模型响应包括合计、平均值，以及频率分布（用于绘制直方图）。

在示例中，我们讲述了如何从四种分布类型中生成随机样本：离散经验分布、离散均匀分布、连续均匀分布和正态分布。（经验分布是指该分布基于数据和给定的概率。）本书网站 www.bcnn.net 上提供的电子表格示例，提供了生成每种分布的 VBA 函数。你应当以它们作为完成本章末尾各习题的起点。

仿真表中的事件是通过使用随机变量以及特定模型的逻辑（公式）而产生的。排队示例，特别是双通道队列的示例，说明了可能出现的某些复杂依赖关系——在所讨论的示例中，是相继来到队列的顾客之间的依赖关系。由于这些复杂性，即使是相对简单的队列网络以及其他中等复杂的模型，临时拼凑的仿真表方法也会失效或变得复杂得难以忍受。基于这一点以及其他理由，需要一种更为系统化的方法，例如事件调度方法。

我们引入了开展一项由多次试验或重复构成的实验这一想法，以及它的必要性。我们看到了由此产生的统计变异性，也看到了开展某种进一步分析的必要性。

## 参考文献

HADLEY G., AND T.M. WHITIN [1963], *Analysis of Inventory Systems*, Prentice-Hall, Englewood Cliffs, NJ.

HARRISON, J.M., AND V. NGUYEN [1995], "Some Badly Behaved Closed Queueing Networks," *Proceedings of IMA Workshop on Stochastic Networks*, eds. F. Kelly and R.J. Williams.

PRITSKER, A. A. B. [1995], *Introduction to Simulation and SLAM II*, 4th ed., John Wiley, New York.

SEILA, A., V. CERIC, AND P. TADIKAMALLA [2003], *Applied Simulation Modeling*, Duxbury, Belmont, CA.

## 习题

### 电子表格与手工习题

若想了解在 Excel 中实现的提示，请研究 www.bcnn.net 上的电子表格解答，并阅读你感兴趣的任何示例的「Explain」工作表。你可能还需要用到电子表格解答中的一个或多个 VBA 函数。因此，为了让这些 VBA 函数可用，最简便的做法是从一个现有的电子表格模型入手——比如与该习题最接近的某个示例——并把 Excel 加载项文件 "BCNNvba.xla" 放在与你的电子表格相同的文件夹中。该加载项文件可在本书网站上找到。只要某个 VBA 函数位于 Excel 的 VBA 环境中的某个 VB 模块里（无论是在电子表格自身中，还是在一个已打开或已被引用的加载项中），你就可以在电子表格单元格的公式中调用它。

对于任何需要你构建电子表格模型的习题，请从某个示例入手，在「One Trial」工作表中删除、修改或增补已有的输入与仿真表，并使用保留下来的部分。在「Experiment」工作表中，把位于"Link"一词正下方那个单元格中的响应改掉，并按你的需要修改标题和列标题中的任何文字。第 1.5 节所述全部函数的 VBA 实现都在加载项文件 "BCNNvba.xla" 中。每个函数对每一个示例电子表格都是可用的，无论该示例当前是否用到了这个函数。

1. 考虑如下连续运转的作业车间。作业的到达间隔时间分布如下：

| 到达之间的时间（小时） | 概率 |
|---|---|
| 0 | .23 |
| 1 | .37 |
| 2 | .28 |
| 3 | .12 |

作业的加工时间服从正态分布，均值 50 分钟，标准差 8 分钟。构建一张仿真表，并针对 10 位新顾客进行仿真。假定仿真开始时，有一件作业正在加工（预计 25 分钟后完成），另有一件加工时间为 50 分钟的作业在队列中。

   (a) 这 10 件新作业在队列中的平均时间是多少？

   (b) 这 10 件新作业的平均加工时间是多少？

   (c) 这 10 件新作业在系统中的最长时间是多少？

2. 一位面包师想要弄清楚每天该烤多少打贝果。贝果顾客人数的概率分布如下：

| 顾客人数/天 | 8 | 10 | 12 | 14 |
|---|---|---|---|---|
| 概率 | 0.35 | 0.30 | 0.25 | 0.10 |

顾客按下述概率分布订购 1、2、3 或 4 打贝果：

| 每位顾客订购的打数 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| 概率 | 0.4 | 0.3 | 0.2 | 0.1 |

贝果售价为每打 \$8.40，制作成本为每打 \$5.80。所有当天未售出的贝果都以半价卖给当地一家杂货店。基于 5 天的仿真，每天应当烤多少打（取整到最接近的 5 打）贝果？

3. 为一个具有 i 个通道的排队系统，设计并解读与图 4 和图 5 类似（analogous）的流程图。

4. 小城出租车（Smalltown Taxi）在上午 9:00 至下午 5:00 期间运营一辆车。目前正在考虑向车队增购第二辆车。对出租车的需求分布如下：

| 两次呼叫之间的时间（分钟） | 15 | 20 | 25 | 30 | 35 |
|---|---|---|---|---|---|
| 概率 | 0.14 | 0.22 | 0.43 | 0.17 | 0.04 |

完成一次服务所需时间的分布如下：

| 服务时间（分钟） | 5 | 15 | 25 | 35 | 45 |
|---|---|---|---|---|---|
| 概率 | 0.12 | 0.35 | 0.43 | 0.06 | 0.04 |

模拟现行系统以及增购一辆出租车后的系统各 5 个单独运行日的运作。就顾客的等待时间以及其他任何有助于说明情况的度量，对这两个系统进行比较。

5. 随机变量 X、Y 和 Z 的分布如下：

$$
\begin{array}{l} {X \sim N (\mu = 100, \sigma^{2} = 100)} \\ {Y \sim N (\mu = 300, \sigma^{2} = 225)} \\ {Z \sim N (\mu = 40, \sigma^{2} = 64)} \end{array}
$$

 对随机变量 

$$
W = \frac {X + Y}{Z}
$$

 仿真 50 个取值。用宽度为 3 的组（区间）编制所得取值的直方图。

6. 已知 A、B、C 是不相关的随机变量。变量 A 服从正态分布，\(\mu = 100\) 且 \(\sigma^{2} = 400\)。变量 B 服从离散均匀分布，其概率分布为 \(p ( b ) = 1 / 5\)，其中 b 为 0、1、2、3、4。变量 C 的分布如下表所示：

| C 的取值 | 概率 |
|---|---|
| 10 | .10 |
| 20 | .25 |
| 30 | .50 |
| 40 | .15 |

用仿真估计新变量 D 的均值，其中 D 定义为 

$$
D = (A - 25 B) / (2 C)
$$

 使用样本量为 100 的样本。

7. 用仿真估计如下运作的库存系统每周平均的销售流失量：

   (a) 每当库存水平降至或低于 10 个单位时，就下达一个订单。同一时刻只能有一笔订单尚未到货。

   (b) 每笔订单的数量等于 20 − I，其中 I 是下达订单时的库存水平。

   (c) 如果在库存水平为零的时段发生需求，则该笔销售流失。

   (d) 每日需求量服从正态分布，均值为 5 个单位，标准差为 1.5 个单位。（仿真过程中把需求量取整到最接近的整数，若结果为负值，则令其需求量为零。）

   (e) 提前期在 0 到 5 天之间均匀分布——只取整数。

   (f) 仿真开始时库存中有 18 个单位。

   (g) 为简单起见，假定订单在营业日结束时下达，并在提前期过去之后收到。因此，若提前期为一天，则该订单在下达订单之后的第二个营业日上午即可供应。

   (h) 让仿真运行 5 周。

8. 某制造厂的一台电梯恰好能载运 400 千克物料。有三种物料装在箱中运抵，乘电梯搭乘。这些物料及其到达间隔时间的分布如下：

| 物料 | 重量（千克） | 到达间隔时间（分钟） |
|---|---|---|
| A | 200 | 5 ± 2（均匀） |
| B | 100 | 6（常数） |
| C | 50 | P(2) = 0.33 |
|  |  | P(3) = 0.67 |

电梯上到二楼需要 1 分钟，卸货需要 2 分钟，返回一楼需要 1 分钟。除非装满货物，否则电梯不离开一楼。模拟该系统 1 小时的运行。物料 A 的一箱的平均运输时间是多少（从其到达直至卸下）？物料 B 的一箱的平均等待时间是多少？1 小时内共有多少箱物料 C 完成了运输？

9. 随机变量 X 和 Y 是实值（非整数）的，在如下区间上均匀分布：

$$
\begin{array}{l} X \sim 10 \pm 10 (\text{uniform}) \\ Y \sim 10 \pm 8 (\text{uniform}) \end{array}
$$


   (a) 对随机变量 

$$
Z = X Y
$$

 仿真 200 个取值。编制所得取值的直方图。Z 的取值范围是什么？其平均值是多少？

   (b) 与 (a) 相同，只是 

$$
Z = X / Y
$$


10. 考虑两块钢板的装配，每块钢板的中心都钻有一个孔。两块钢板通过一根销钉连接。装配时，两块钢板以左下角 (0,0) 为基准对齐。孔的位置在每块钢板上都居于 (3,2)。各个方向上的标准差为 0.0045。孔径服从正态分布，均值为 0.3，标准差为 0.005。销钉直径同样服从正态分布，均值为 0.29，标准差为 0.004。有多大比例的销钉能够穿过装配好的钢板？你的答案应基于 50 次观测的仿真。

$$
[ H i n t ]: \text{Clearance} = \operatorname{Min} (h_{1}, h_{2}) - [ (x_{1} - x_{2}) ^ {2} + (y_{1} - y_{2}) ^ {2} ] ^ {. 5} - p, \text{where}
$$


    h 孔径，i 表示板 1、板 2

    p 销钉直径

    x 到钢板孔中心的距离，水平方向，i 表示板 1、板 2

    y 到钢板孔中心的距离，垂直方向，i 表示板 1、板 2

11. 在上面的习题中，若销钉过松，它就会晃动。当 Min(h1, h2) − p ≥ 0.006 时发生晃动。有多大比例的装配会晃动？（条件概率——即在销钉能够穿过的前提下。）

12. 在一个圆的圆周上随机选取三个点。用蒙特卡罗抽样方法估计三点全部落在同一个半圆上的概率。执行 5 次重复。

13. 统计学中的两个定理如下：

    定理 1 设 \(Z_{1} , Z_{2} , \ldots , Z_{k}\) 是相互独立且服从正态分布的随机变量，均值 \(\mu = 0\)、方差 \(\sigma^{2} = 1\)。则随机变量 

$$
\chi^{2} = Z_{1} ^ {2} + Z_{2} ^ {2} + \dots + Z_{k} ^ {2}
$$

 称为服从自由度为 k 的卡方分布，简称 \(\chi_{k} ^ {2}\)

    定理 2 设 \(Z \sim N ( 0 , 1 )\)，V 为自由度为 k 的卡方随机变量。若 Z 与 V 相互独立，则随机变量 

$$
T = \frac {Z}{\sqrt {V / k}}
$$

 称为服从自由度为 k 的 t 分布，简称 \(t_{k}\)

    生成一个自由度为 3 的 t 分布随机变量的样本。按所示顺序，使用下列随机值中所需的尽可能多个：

| 随机数 | 随机正态数 |
|---|---|
| 0.6729 | 1.06 |
| 0.1837 | -0.72 |
| 0.2572 | 0.28 |
| 0.8134 | -0.18 |
| 0.5251 | -0.63 |

    从 t 分布中生成一个样本，需要多少个随机数和标准随机正态数？

14. 某银行有一名免下车服务柜员，另有一处可供一位顾客等待的空间。当队列已满时到达的顾客会把车停好，走进银行办理业务。到达之间的时间以及服务时间的分布如下：

| 到达之间的时间（分钟） | 概率 | 服务时间（分钟） | 概率 |
|---|---|---|---|
| 0 | 0.09 | 1 | 0.20 |
| 1 | 0.17 | 2 | 0.40 |
| 2 | 0.27 | 3 | 0.28 |
| 3 | 0.20 | 4 | 0.12 |
| 4 | 0.15 |  |  |
| 5 | 0.12 |  |  |

    模拟免下车服务柜员为 10 位新顾客服务的运作。这 10 位新顾客中的第一位在一个随机确定的时刻到达。仿真开始时，有一位顾客正在接受服务（在时刻 3 离开），另有一位在队列中。有多少顾客走进银行办理了业务？

15. 正在估算雪松沼泽湖（Cedar Bog Lake）所覆盖的acre 数。下面给出了该湖的图形，图上标出了以英尺为单位的比例尺：

![雪松沼泽湖](img/ch02/img_57b362c00eb3df3d.jpg)

    用仿真估计该湖的面积。（1 acre 等于 43,560 平方英尺。）

16. 在例 6 中，假设增设另一位服务台 Charlie 之后的平均延迟几乎为零。但另一位服务台的成本为每小时 \$20。如果平均每位来电者的延迟成本为每小时 \$200，那么在什么情况下增设另一位服务台才是可取的？

17. 在例 12 中，为什么各条路径不是各占约 1/3 的时间？

18. 在例 11 中，要得到 3 或 4 的提前期需求，唯一的可能是提前期为一天且需求量为 3 或 4。这种情况发生的概率是 0.36(0.20+0.35) = 0.198。用这一逻辑计算，提前期需求等于 11 或 12 的概率是多少？等于 18 的概率又是多少？（用计算，不要用仿真。）

19. 将例 4（掷硬币游戏）中的实验运行 4000 次试验。解释为什么"Harry 在 100 次投掷中的 5 次或更少次上领先"的概率，应当大于"Harry 在 95 次或更多次投掷中领先"的概率。为什么关于 Tom 你也能说出同样的道理——而且并不矛盾。

### 基于电子表格示例的习题

20. 修改第 1.5 节中的例 2，改用"WRONG !!"公式来生成服务时间。生成大量服务时间，比如 400 个，并编制所得频数的直方图。对于每个服务时间取值，样本频数与所指定的概率有多接近？解释为什么会不同。也就是说，解释这个公式为什么是错的。

21. 在例 5 中，令到达分布为 1 到 10 分钟之间的均匀分布（只取整数）。运行 50 次试验的实验。所显示的组（区间）频数与图 10 中的相比有差别吗？

22. 在例 5 中，把服务时间分布改为如下：

| 服务时间（分钟） | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| 概率 | 0.05 | 0.10 | 0.20 | 0.30 | 0.25 | 0.10 |

    运行 50 次试验的实验。此处的平均等待时间与修改之前的例 5 中的平均等待时间有差别吗？如果有差别，你把它归因于什么？

23. 将例 5 中的实验分别运行 25、50、100、200 和 400 次试验。（所有试验都使用不同的随机数集合。）平均等待时间的最小值与最大值有什么差异？如果存在差异，你如何解释？

24. 在例 6 中，运行 10 个实验，每个实验包含 100 次试验。在每个实验中，有多大比例的试验其来电者的平均延迟为 3 分钟或更少？对于刚刚仿真的这 1000 次试验，又有多大比例的来电者其平均延迟为 2 分钟或更少？

25. 在例 6 中，开展一次包含 400 次试验的实验。解释最小平均延迟与最大平均延迟之间为何存在如此大的差距。

26. 在例 6 中，运行 10 个各含 50 次试验的实验，以及 10 个各含 400 次试验的实验。解释你所观察到的这两类实验在输出上的差异（如果有的话）。

27. 在例 6 中，修改电子表格，使 Able 和 Baker 各自接听的电话数能够显示出来。在 400 次试验的实验中，各接听了多大比例？

28. 在例 7 中，基于 50 次试验的实验，确定订购报纸的最佳策略。

29. 在例 7 中，运行 10 个实验，每个实验含 400 次试验。你所得的每日利润最大值与最小值之间的最大差异是多少？

30. 使用例 8 的电子表格。对于 100 次试验，若 M = 11 且 N = 4、5、6，平均期末库存是多少？

31. 使用例 8 的电子表格。对于 100 次试验，若 N = 5 且 M = 10、11、12，平均期末库存是多少？

32. 以例 8 的解答为起点，修改该解答以模拟销售流失的情形；也就是说，把短缺改为销售流失。记录流失的销售量。

33. 在例 9 中，轴承成本为何值时，现行系统与所提系统的每 10,000 小时总成本实际上相同？你的估计应基于现行系统和所提系统各 10 个实验，每个实验包含 400 次试验。

34. 把例 9（现行或所提）的 Experiment 工作表上直方图的组（区间）宽度改为 \$50，起始位置低于最小值（例如，若最小值为 \$1528.46，则让第一个单元格从 \$1500 开始）。这样做有什么好处？

35. 使用例 9（所提策略）的电子表格，运行 10 个实验，每个实验含 40 次试验，并记录结果的极差（最大值减最小值）。然后计算平均极差。接着按同样的方式再做一遍，但每个实验使用 400 次试验。如果存在差异，你如何解释？

36. 在例 10 的电子表格中，把 \(\sigma_{x} = 600\) 米和 \(\sigma_{y} = 300\) 米设定好。保持目标不变。开展一次包含 200 次试验的仿真。平均命中数是多少？

37. 在例 10 的电子表格中，把 \(\sigma_{x} = 300\) 米和 \(\sigma _ {\mathrm{y}} = 600\) 米设定好。保持目标不变。开展一次包含 200 次试验的仿真。平均命中数是多少？

38. 在例 10 的电子表格中，把 \(\sigma_{x} = 300\) 米和 \(\sigma_{y} = 150\) 米设定好。保持目标不变。开展一次包含 200 次试验的仿真。平均命中数是多少？解释本习题与上一个习题之间的差别。

39. 在例 10 的电子表格中，把 \(\sigma_{x} = 2 \sigma_{y}\) 米设定好。保持目标不变。如果基于一次包含 400 次试验的实验，平均命中数约为 6.0，那么 \(\sigma_{x}\) 的取值是多少？

40. 在例 11 中，假设你希望对提前期需求的分布有更好的估计：你应当 (1) 增加「One Trial」工作表上的周期数，还是 (2) 增加「Experiment」工作表上的试验次数？解释你的答案。

41. 在例 11 中，假设你希望对提前期需求的均值有更好的估计，你应当怎么做？

42. 在例 11 中，令可能取值 3、4、5、6 的需求概率相等。运行一次包含 400 次试验的实验。比较使用原始输入数据与使用新输入数据所得的平均提前期需求。

43. 在例 11 中，令可能取值 1、2、3 的提前期概率相等。运行一次包含 400 次试验的实验。比较使用原始输入数据与使用新输入数据所得的平均提前期需求。

44. 在例 12 中，把实验运行 20 次，每次含 400 次试验。记录中间路径最长的次数。你对中间路径被选中的概率所作的最佳估计是多少？

45. 在上面的习题中，你遇到的中间路径被选中次数的最小值和最大值分别是多少？如果你只做了一次仿真，恰好得到了那个最小值（或最大值），并把它作为结果报告出来，你会如何为自己开脱？

46. 在例 12 中，假定第三条路径（即图 20 中最下面那条）改为由六个 U(1,2) 活动组成。修改电子表格以适应这一变化。现在最常出现的路径是哪一条？本习题给你什么启发？

47. 使用 Excel，用公式 =-10\*LN(RAND()) 在一列中生成 1000 个值。

    (a) 计算该列数据的描述性统计量，包括最小值、最大值、均值、中位数和标准差。

    (b) 把这些值按 10 个组（区间）制表，每个组（区间）宽度为 5，第一个组（区间）从 0 开始，第 11 个组（区间）用于溢出值（如果有）。

    (c) 该直方图像你所熟悉的某种分布吗？如果像，它叫什么名字？

    提示：在 Excel 公式中使用 FREQUENCY 或 COUNT 来统计各组（区间）的频数。

48. 使用 Excel，用公式 =RAND() 生成 12 列，每列 250 个值。

    在 M1 单元格中放入公式 =SUM(A1:L1)-6，并把它复制到 M1 列下方的 249 个单元格中。

    (a) 计算 M 列数据的描述性统计量，包括最小值、最大值、均值、中位数和标准差

    (b) 把这些值按 9 个组（区间）制表。第一个组（区间）包含所有小于或等于 3.5 的值；接下来的 6 个组（区间）宽度为 1；最后一个组（区间）包含所有大于 3.5 的值。

    (c) 该直方图像你所熟悉的某种分布吗？如果像，它叫什么名字？

    提示 1：在 Excel 公式中使用 FREQUENCY 或 COUNT 来统计各组（区间）的频数。

    提示 2：M 列中的数值可以作为表 A.2 中数值的粗略近似。

49. 使用 Excel，生成一张可用来代替表 A.1 的表，公式为 =INT(RAND()\*100000)。

50. 用 Excel 设计一个蒙特卡罗实验来估计 π 的值，如同圆面积的公式 

$$
\mathrm{A} = \pi \mathrm{r} ^ {2}
$$

 提示：考虑单位正方形内一个半径为 1 的四分之一圆。向单位正方形投掷一枚随机飞镖。

51. 一名大学生正在为来年做财务规划。支出包括下列已知成本：

    学费 \$8,400

    住宿 \$5,400

    另外还有许多不确定的支出，其范围如下：

    伙食 \$900–\$1,350

    娱乐 \$600–\$1,200

    交通 \$200–\$600

    书本 \$400–\$800

    这名学生预计来年会有收入。其中只有下列金额是确定的：

    州奖学金 \$3,000

    父母资助 \$4,000

    部分收入是可变的，其范围如下：

    餐厅侍应 \$3,000–\$5,000

    图书馆助理 \$2,000–\$3,000

    这名学生预计需要一笔贷款。开发一个电子表格模型，并用 1000 次试验来预测贷款金额的概率分布。把任何金额可变的项目都视为在所给区间上均匀分布。

52. 22 档的 Bally 1090 老虎机，其每个转轴上的符号数量如下：

| 符号 | 转轴 1 | 转轴 2 | 转轴 3 |
|---|---|---|---|
| 樱桃 | 2 | 5 | 2 |
| 橙子 | 2 | 3 | 7 |
| 李子 | 5 | 1 | 10 |
| 铃铛 | 10 | 2 | 1 |
| 条杠 | 2 | 10 | 1 |
| 7 | 1 | 1 | 1 |

    中奖赔付如下：

| 组合 | 赔付 | 组合 | 赔付 |
|---|---|---|---|
| 左樱桃 | 2 | 三个李子 | 20 |
| 右樱桃 | 2 | 李子李子条杠 | 14 |
| 任意两个樱桃 | 5 | 条杠李子李子 | 14 |
| 三个樱桃 | 20 | 三个铃铛 | 20 |
| 三个橙子 | 20 | 铃铛铃铛条杠 | 18 |
| 橙子橙子条杠 | 10 | 条杠铃铛铃铛 | 18 |
| 条杠橙子橙子 | 10 | 三个条杠 | 50 |
|  |  | 三个 7 | 100 |

    使用 Excel，模拟这台老虎机 1000 次转动，每次成本 \$1。

    盈亏是多少？

## www.bcnn.net 上可获取的资料清单

Example2.1CoinToss.xls - 掷硬币

Example2.2ServiceTimes.xls - 服务时间

Example2.3ArrivalTimes.xls - 到达时间

Example2.4CoinTossGame.xls - 掷硬币问题

Example2.5SingleServer.xls - 杂货店（单服务台队列）

Example2.6AbleBaker.xls - 呼叫中心（双通道队列）

Example2.7NewsDealer.xls - 报刊销售商的库存问题

Example2.8RefrigInventory.xls - 冰箱销售（库存订货）

Example2.9CurrentBearings.xls - 现行轴承更换策略

Example2.9ProposedBearings.xls - 所提轴承更换策略

Example2.10Target.xls - 一次轰炸任务

Example2.11LeadTimeDemand - 提前期需求分布

Example2.12Project.xls - 项目作业网络

BCNNvba.xla - 用于运行各示例的 Excel 加载项（VBA 代码）

Chap2Excel5ed.zip - 包含全部 Excel 示例的 zip 文件