全书目录

1 仿真导论

仿真（simulation）是对真实世界过程或系统运行过程的模仿。无论采用手工方式还是计算机方式，仿真都包含两个要素：生成系统的一段人造历史（artificial history），并观察这段人造历史，从而推断真实系统的运行特性。

研究系统随时间演化的行为，方法是建立仿真模型（simulation model）。该模型通常以一组关于系统如何运行的假设的形式出现。这些假设通过系统中各实体（entity，即关注的对象）之间的数学、逻辑和符号关系来表达。模型一旦建立并通过验证（validated），就可以用它来研究关于真实世界的各种"如果……会怎样"（what if）问题。可以先对系统的潜在变更进行仿真，以预测其对系统性能的影响。仿真也可用于研究尚处于设计阶段、尚未建成的系统。因此，仿真建模既可作为分析工具，用于预测对既有系统的变更所产生的效果，也可作为设计工具，用于预测新系统在各种不同情形下的性能。

在某些情况下，所建立的模型足够简单，可以用数学方法"求解"。这类解可能通过微分学、概率论、代数方法或其他数学技术获得。解通常由一个或多个数值参数构成，这些参数称为系统的性能度量（measures of performance）。然而，许多真实世界系统极其复杂，其模型几乎不可能用数学方法求解。在这些情况下，可以采用基于计算机的数值仿真来模仿系统随时间的行为。从仿真中收集数据，就如同在观察一个真实系统一样。这些由仿真生成的数据被用来估计系统的性能度量。

本章首先讨论何时应当使用仿真、仿真的优点与缺点，以及仿真的实际应用领域。随后探讨系统与模型的概念。最后给出建立并使用系统仿真模型的步骤概要。

1 仿真何时是恰当的工具

专用仿真语言的出现、运算成本不断下降所带来的巨量计算能力，以及仿真方法学的进展，使仿真成为运筹学与系统分析中被最广泛使用和认可的工具之一。何种情形下适合使用仿真，已有众多作者作过讨论，从 Naylor 等人 [1966] 到 Shannon [1998]。仿真可用于以下目的：

1. 仿真使得对复杂系统内部交互作用、或复杂系统内某个子系统的交互作用进行研究、开展实验成为可能。
2. 信息、组织和环境方面的变更可以先进行仿真，从而观察这些改变对模型行为的影响。
3. 在设计仿真模型过程中所获得的知识，对提出针对所研究系统的改进建议极有价值。
4. 改变仿真输入并观察相应的输出，可以产生有价值的洞察：哪些变量最为重要，以及变量之间如何相互作用。
5. 仿真可以用作一种教学手段，以强化解析求解方法。
6. 仿真可以在实施之前对新的设计或政策进行试验，从而为可能发生的事情做好准备。
7. 仿真可以用来验证解析解。
8. 对机器的不同能力配置进行仿真，有助于确定其需求。
9. 为培训而设计的仿真模型，使学习成为可能，且无须付出在职培训的成本与干扰。
10. 动画可以展示系统的仿真运行过程，使方案得以可视化。
11. 现代系统（工厂、晶圆厂、服务组织等）极其复杂，其内部交互作用只能通过仿真来处理。

2 仿真何时不适用

本节基于 Banks 与 Gibson [1997] 的一篇文章，该文给出了十条用于判断何时不宜使用仿真的规则。第一条规则指出：当问题可以凭常识解决时，不应使用仿真。文中举了一个汽车牌照服务设施的例子，顾客随机到达，平均到达率为 100 人/小时，而平均服务率为 12 人/小时。为了确定所需服务人员的最小数量，并不需要仿真。只须计算 100/12 = 8.33，即可得出需要 9 名或更多服务人员的结论。

第二条规则指出：若问题可以用解析方法求解，则不应使用仿真。例如，在某些条件下，上述例子中的平均等待时间可以利用 www.bcnn.net 提供的工具求得。

第三条规则指出：若直接做实验的成本更低，则不应使用仿真。文中以一家快餐店汽车穿梭通道（drive-in）为例：为了确定增设一个接单点对顾客等待时间的影响，安排一个人使用手持终端并通过语音接单，其成本低于开展仿真研究。

第四条规则指出：若成本超过收益，则不应使用仿真。如第 12 节将要讨论的，完成一次仿真包含许多步骤，且每一步都必须扎实地完成。倘若一项仿真研究花费 2 万美元，而可能带来的节省仅为 1 万美元，那么仿真就不合适。

第五和第六条规则指出：若资源或时间不具备，则不应开展仿真。如果估计仿真成本为 2 万美元，而只有 1 万美元可用，那么建议不要贸然启动仿真研究。同理，若决策需要在两周内做出，而仿真要花一个月，那么也不建议开展仿真研究。

仿真需要数据，有时候需要大量数据。若没有数据可依，连估计值都没有，则不建议使用仿真。下一条规则涉及能否验证和确认模型。若时间不够，或人手不足，仿真就不适用。

如果管理者抱有不切实际的期望，过早要求过多东西，或者高估了仿真的能力，那么仿真也可能不适用。

最后，若系统行为过于复杂或无法定义，仿真就不适用。人类行为有时极其复杂，难以建模。

3 仿真的优点与缺点

仿真对客户而言具有直觉上的吸引力，因为它模拟了真实系统中发生的事情，或模拟了处于设计阶段的系统给人的感知。仿真的输出数据应当与从真实系统中可记录到的输出直接对应。此外，还可以在不依赖可疑假设（诸如"每个随机变量服从同一统计分布"这类数学可解模型的假设）的前提下，为系统建立仿真模型。由于这些以及其他原因，仿真常常成为问题求解的首选技术。

与优化模型不同，仿真模型是被"运行"（run）而非被"求解"（solved）的。给定一组特定的输入与模型特征，模型被运行，仿真行为被观察。这一改变输入与模型特征的过程，产生一组待评估的方案（scenario）。随后针对这些方案给出推荐方案，无论是对既有系统的分析还是对新系统的设计，最终都可将推荐的解决方案付诸实施。

仿真有许多优点，但也有一些缺点。Pegden、Shannon 与 Sadowski [1995] 列举了这些优缺点。部分优点如下：

1. 可以在不干扰真实系统日常运行的前提下，探索新的政策、作业流程、决策规则、信息流、组织流程等。
2. 可以对新的硬件设计、物理布局、运输系统等进行试验，而无须为购置它们而投入资源。
3. 关于某些现象如何发生或为何发生的假设，可以就其可行性进行检验。
4. 时间可以被压缩或拉伸，从而使所研究的现象得以加速或减速。
5. 可以获得关于变量之间相互作用的洞察。
6. 可以获得关于各变量对系统性能重要性的洞察。
7. 可以进行瓶颈分析，以发现在制品、信息、材料等被过度延误的环节。
8. 仿真研究有助于理解系统实际上是如何运行的，而不仅仅停留在人们认为系统如何运行的层面上。
9. “如果……会怎样"类问题可以得到解答。这在新系统设计中尤其有用。

部分缺点如下：

1. 建模需要专门训练。这是一门需要随时间、经由经验习得的技艺。而且，若两个模型由不同的行家分别构建，它们或许有相似之处，但完全相同的可能性极低。
2. 仿真结果可能难以解释。多数仿真输出本质上是随机变量（它们通常基于随机输入），因此很难分辨某个观察结果究竟源于系统内部的相互关系，还是源于随机性。
3. 仿真建模与分析可能既耗时又昂贵。在建模与分析上削减资源，可能导致仿真模型或分析不足以胜任任务。
4. 在某些情况下，即便解析解存在甚至更可取，人们仍会使用仿真，正如第 2 节所讨论的。这在某些排队系统的仿真中可能尤为明显，因为那里已有闭式排队模型（closed-form queueing model）可用。

为仿真辩护，针对上述四点缺点分别可以作如下回应：

1. 仿真软件厂商一直在积极开发包含若干模型的软件包，这些模型只需输入数据即可运行。这类模型有"模拟器”（simulator）或"模板"（template）这样的通用称谓。
2. 许多仿真软件厂商已在其软件包中开发了输出分析功能，用以进行非常细致的分析。
3. 今天的仿真比昨天做得更快，明天还会更快，这既得益于硬件的进步（可以快速运行方案），也得益于众多仿真软件包的进步。例如，某些仿真软件包含用于物料搬运建模的构件，可使用叉车、输送带和自动导引车（AGV）这类运输设备。
4. 闭式模型无法分析实践中遇到的大多数复杂系统。在本书两位作者多年的咨询实践中，没有遇到过一个可以用闭式解求解的问题。

4 应用领域

仿真的应用极为广泛。冬季仿真会议（Winter Simulation Conference, WSC）是了解仿真应用与理论最新进展的绝佳途径。会议还提供从入门到高阶的大量教程。WSC 由六个技术学会和美国国家标准与技术研究院（NIST）赞助。这六个技术学会分别是：美国统计学会（ASA）、美国计算机协会/仿真专业兴趣组（ACM/SIGSIM）、电气电子工程师学会系统、人类与控制论学会（IEEE/SMCS）、工业工程师学会（IIE）、运筹学与管理科学学会仿真分会（INFORMS-SIM），以及国际建模与仿真学会（SCS）。关于即将召开的 WSC 的信息可从 www.wintersim.org 获取。含完整论文的 WSC 会议文集可从 www.informs-cs.org/wscpapers.html 获取。下面按领域列出最近一次 WSC 的部分报告：

制造业应用

• 选择最合适瓶颈检测方法的方法论
• 船厂制造模型的自动化开发
• 制造工程过程中的仿真（Emulation）应用
• 面向提升制造绩效的优化维护设计
• 汽车零部件行业的生产率管理
• 日本汽车制造厂的产线设计

晶圆制造

• 批次分配到机台的范式转变
• 带有重入过程的多芯片封装装配线调度
• 封装—测试设施的执行层产能分配决策
• 借助回路控制管理在制品（WIP）与周期时间

业务流程处理

• 服务请求分派问题的新策略
• 流程执行监控与调整方案
• 零售门店的店内陈列
• 小型零售店的销售预测

建筑工程与项目管理

• 多个建筑工地有限弯曲机（bar-bender）的调度
• 重复性项目的施工
• 面向公路建设项目改进规划的交通组织
• 隧道竖井施工建模模板
• 隧道施工规划决策支持工具

物流、运输与配送

• 驳船运输系统的运营策略
• 生物恐怖事件下应急医疗物资的发放方案
• 复杂邮件运输网络的分析
• 提升集装箱码头性能
• 基于实时数据的场桥调度
• 航空运营中的单元装载设备库存
• 基于预测更新策略的库存系统
• 食品配送中心的泊位分配
• 驳船运输系统的运营策略

军事应用

• 多国战区内部物流配送
• 加拿大军队行动的可持续性前景研究
• 替换 MK19 自动榴弹发射系统的可行性研究
• 为非对称作战训练联合部队
• 多目标无人机任务规划
• 由作战需求驱动的联邦的发展

医疗保健

• 缩短预约提前期与降低患者爽约率的干预措施
• 支持智慧型思维以改善医院绩效
• 急诊室流程精益改进的验证
• 缓解急诊科拥挤状况
• 易腐药品的库存建模
• 门诊手术中心的实施
• 传染病控制政策
• 手术室与麻醉后资源的平衡
• 结直肠癌筛查的成本效益分析

其他应用

• 采用多重反馈控制的劳动力资源行动管理
• 分析球洞大小对推杆入洞的影响
• 粒子滤波在野火蔓延仿真中的应用
• 封闭生境中的捕食者—被捕食者关系
• 集约化仔猪生产系统
• 呼叫中心的实时延迟估计
• 某公立大学的大流行性流感应急预案

在 ICS Newsletter（ICS 通讯）上一篇关于仿真未来的文章（Banks, 2008）中，十六位仿真研究者——其中包括该领域知名的专家——就以下问题给出了各自的回答：“从长远来看，比如说三年之后，我们将在仿真软件上看到哪些非凡的成果？“以下回答即出自该文：

• 在诸如基于智能体建模（agent-based modeling）这类全新思想与方法被软件厂商采纳之后（这将在几年内发生），仿真建模的进展将由计算能力的提升所驱动。例如，在某个时点上，我们将能够仿真超大型供应链和制造设施的详细运行过程。
• 仿真建模将更像一种"装配"活动，而非一种"从零开始建造"的活动。建模时将使用智能的、参数化的组件来拼装模型，而不是定义大量细化的逻辑。
• 在解决困难问题方面将取得进展。要为真正困难的问题提供工具，仿真软件开发者必须回到起点，认真重新思考一个根本问题：“谁提供力量，谁提供远见？“把这两者的组合摆正，对进一步进步至关重要。当软件开发者试图同时提供力量与远见时，最终用户将被迫接受开发者的范式。这对简单问题很有效，对困难问题则效果不佳。对于所有真正困难的问题，用户对问题的了解往往比软件开发者更多显而易见的解决方案是：由开发者提供经过深思熟虑的真正原语集合，以及组合和打包这些原语的方式。
• 仿真软件将与控制软件更紧密地集成。
• 建模者将拥有一个在组织内跨应用共享的单一模型。
• 仿真应用将不仅限于设计类应用，还将用于组织内部的日常运营决策。
• 运行在强大服务器上的仿真建模将通过基于网页的接口来访问。
• 对人类活动更好、更方便的建模将成为可能（例如把基于智能体的模型嵌入离散事件模型之中）。
• 可以预期会有更多协作式的仿真项目开发。
• 将采纳接口标准并融入 Web 服务，使仿真软件不仅能作为一个整合的联邦相互协作，还能标准化并简化其他应用与仿真软件之间的接口方式。
• 解析求解技术（如线性规划）将与仿真能力相结合。
• 那些"非凡"的成果将源自计算技术其他领域和软件工程方面的进展。但没有什么能超越由 SIMULA 67 引入软件主流开发之中的面向对象范式，而 SIMULA 67 恰恰源自仿真软件这一社区。

5 仿真的若干近期应用

本节介绍仿真的若干近期应用。这些案例均见于所标注的文献中，您可以查阅这些案例并了解更多细节。

案例 1

• 标题： “土耳其军队使用仿真对燃料供应系统建模与优化”
• 作者： I. Sabuncuoglu, A. Hatip
• 报道于： 2005 年 11—12 月，《Interfaces》
• 挑战： 对土耳其军队的燃料供应系统进行分析。
• 技术方法： (1) 在各种情形下测量既有系统与拟建系统的绩效；(2) 开发了一个基于遗传算法的仿真优化模型以优化系统性能；(3) 开展了大量仿真实验。
• 节约金额： 数百万美元。

案例 2

• 标题： “PLATO 帮助雅典夺得金牌：面向组织变革与资源管理的奥运会知识建模”
• 作者： D.A. Beis, P. Poucopoulos, Y. Pyrgiotis, K.G. Zografos
• 报道于： 2006 年 1—2 月，《Interfaces》
• 挑战： 开发一套用于规划和设计场馆运营的系统化流程；建立一个内容丰富的模型库，可直接迁移给未来的奥运会组委会及其他体育赛事。
• 技术方法： 基于仿真与其他决策分析方法的知识建模与资源管理技术及工具。
• 节约金额： 超过 6970 万美元。

案例 3

• 标题： “斯伦贝谢使用仿真进行陆地地震勘测的投标与实施”
• 作者： P.W. Mullarkey, G. Butler, S. Gavirneni, D.J. Morrice
• 报道于： 2007 年 3—4 月，《Interfaces》
• 挑战： 快速且准确地度量地震勘测的成本。
• 技术方法： 开发了一套仿真工具，用于评估队伍规模、勘测区域、地理区域和天气条件对勘测成本与工期的影响。
• 节约金额： 每年 150 万至 300 万美元。

案例 4

• 标题： “运筹学助力 Standard Register 公司重塑运营战略”
• 作者： S.L. Ahire, M.F. Gorman, D. Dwiggins, O. Mudry
• 报道于： 2007 年 11—12 月，《Interfaces》
• 挑战： 在竞争激烈的传统印刷市场中，最小化提供有竞争力定价所需的总成本。
• 技术方法： (1) 用回归分析估计成本与时间属性；(2) 用优化建模确定订单路由策略；(3) 对生产—配送网络进行仿真建模。
• 节约金额： 每年超过 1000 万美元。

案例 5

• 标题： “仿真为服务业实现需求驱动的劳动力调度器”
• 作者： M. Zottolo, O.M. Ulgen, E. Williams
• 报道于： 2007 年冬季仿真会议论文集，编者 S.G. Henderson, B. Biller, M.-H. Hsieh, J. Shortle, J.D. Tew, R.R. Barton。
• 实施方： PMC（www.pmcorp.com）
• 客户： 美国某大型家电公司。
• 挑战： 由于排班过程耗时、依赖手工且效率低下，客户出现了工作组排班过量与不足、所提供服务水平不一致、现场管理人员不堪重负等问题。
• 技术方法： (1) 建立仿真模型，依据一天中客户需求的变化和服务时间为不同工作组排班；(2) 确定作业标准；(3) 开发用于输入数据的接口以及排班的存储与发布功能；(4) 实施该调度工具并培训员工使用。
• 节约金额： 估计为其美国工厂每年节约 8000 万美元。

案例 6

• 标题： “仿真改进某家电生产商的产线末端分拣与物料搬运取货调度”
• 作者： N. Kale, M. Zottolo, O.M. Ulgen, E. Williams
• 实施方： PMC（www.pmcorp.com）
• 客户： 美国某大型汽车租赁公司。（原文如此印）
• 报道于： 2007 年冬季仿真会议论文集，编者 S.G. Henderson, B. Biller, M.-H. Hsieh, J. Shortle, J.D. Tew, R.R. Barton。
• 挑战： 确定把若干种成品家电（SKU）分配到 12 条分拣通道的最有效方法，以尽量减少物料搬运作业和取货所需的叉车数量。由于产品组合波动大、SKU 数量众多、管理费用高昂，这一分拣任务及其相关的物料搬运任务复杂到值得采用离散事件仿真来分析。
• 技术方法： 建立仿真模型，并就以下方面开展实验：(1) 不同 SKU 向不同通道的分配；(2) 取货策略与取货批量；(3) SKU 的批量规模。
• 节约金额： 估计每年节约 10 万美元，而一次性资本投入仅 5 万美元。

6 系统与系统环境

要对一个系统建模，必须理解系统的概念以及系统边界（system boundary）的概念。系统（system）被定义为：一组以某种有规律的交互或相互依赖关系相结合、为达成某一目的而存在的对象。汽车制造生产系统即是一例：机器、零部件和工人沿装配线协同作业，生产出高质量的整车。

系统常常受到系统之外所发生变更的影响。这类变更被称为发生在系统环境（system environment）中 [Gordon, 1978]。在为系统建模时，必须界定系统与其环境之间的边界。这一决定可能取决于研究的目的。

以工厂系统为例，控制订单到达的因素可能被视为处于工厂影响之外，因而属于环境的一部分。然而，如果要考虑供给对需求的影响，工厂产出与订单到达之间就存在一种关系，而这种关系必须视为系统的活动。同样地，以银行系统为例，可能存在可支付最高利率的上限。对于研究单家银行而言，这将被视为环境施加的约束；但在研究货币法规对银行业影响的研究中，该上限的设定则属于系统的活动 [Gordon, 1978]。

7 系统的组成要素

为了理解与分析一个系统，需要定义若干术语。实体（entity）是系统中关注的对象；属性（attribute）是实体的性质；活动（activity）代表一段规定长度的时间。如果研究的是银行，那么顾客可能是实体之一，其支票账户余额可能是属性，而存款则是一项活动。

构成某项研究之系统的实体集合，可能只是另一项研究所针对的整个系统的一个子集 [Law, 2007]。例如，如果研究上述银行是为了确定办理存取款业务所需的柜台员人数，那么系统可以定义为银行中由普通柜台员和排队等待的顾客构成的那一部分。若研究目的扩展到确定所需的专职柜台员人数（用以开具银行本票、办理商业交易等），则系统的定义必须相应扩大。

系统的状态（state）被定义为：为在任意时刻描述系统、相对于研究目标而言所必需的那一组变量。在银行研究中，可能的状态变量包括忙碌的柜台员人数、排队等待或正在接受服务的顾客人数，以及下一位顾客的到达时间。事件（event）被定义为可能改变系统状态的一次瞬时发生。内生（endogenous）一词用于描述发生在系统内部的活动和事件；外生（exogenous）一词用于描述环境中那些影响系统的活动和事件。在银行研究中，顾客到达是外生事件，顾客服务完成则是内生事件。

表 1 列出了若干系统中实体、属性、活动、事件和状态变量的示例。这里仅给出系统组成部分的部分列表。除非明确研究目的，否则无法给出完整的清单。依据研究目的不同，系统中受到关注的方面也不同，从而可以补全这一组成要素清单。

表 1 系统及其组成要素示例

8 离散系统与连续系统

系统可以被划分为离散型或连续型。“在实践中，完全离散或完全连续的系统很少，但由于对大多数系统而言总有一种变化占主导地位，因此通常可以把一个系统归类为离散型或连续型” [Law, 2007]。离散系统（discrete system）是指状态变量仅在一组离散时间点上发生变化的系统。银行就是一个离散系统的例子：状态变量"银行内的顾客人数"仅在顾客到达或顾客服务完成时才发生变化。图 1 展示了顾客人数如何仅在离散时间点上发生变化。

图 1 离散系统的状态变量。

连续系统（continuous system）是指状态变量随时间连续变化的系统。水坝后的水位即是一例。在暴雨期间及暴雨结束后的若干时间内，雨水流入坝后的湖中。水从大坝中引出，用于防洪和发电。蒸发也会使水位下降。图 2 展示了这一连续系统的状态变量——坝后水位——的变化过程。

图 2 连续系统的状态变量。

9 系统的模型

有时我们希望研究一个系统，以理解其组成部分之间的关系，或预测该系统在新政策下将如何运行。研究系统时，有时可以直接对系统本身做实验。但这并不总是可行的。一个新系统可能尚不存在，它可能只处于假想状态或设计阶段。即便系统已经存在，对其做实验也可能不切实际。就银行而言，为了研究减少柜台员人数对排队长度的影响而削减人员，可能会激怒顾客，以致他们把账户转往竞争对手。因此，对系统的研究往往要借助系统的模型来完成。

我们曾承担过一项咨询业务：为西澳大利亚州一个重新设计的港口做仿真。对于一个造价 3 亿美元的装卸泊位而言，投入如此巨资后再发现该泊位不适合承担任务，是不明智的。

模型（model）被定义为：为了研究某系统而对该系统所作的表示。对于大多数研究而言，只需考虑系统中那些影响所研究问题的方面。这些方面就体现在系统的模型中；模型按其定义就是对系统的简化。另一方面，模型又应足够详细，以便能对真实系统得出有效的结论。随着研究目的的变化，同一系统可能需要不同的模型。

正如系统的组成要素包括实体、属性和活动一样，模型也由这些要素来表示。不过，模型中只包含与研究相关的那些组成要素。

10 模型的类型

模型可以分为数学模型和物理模型。数学模型（mathematical model）使用符号记法和数学方程来表示系统。仿真模型（simulation model）是系统的一类特定数学模型。物理模型（physical model）是某个实物放大或缩小后的版本，例如放大的原子模型或缩小版的太阳系模型。

仿真模型还可以进一步分类为静态或动态、确定性或随机、离散或连续。静态仿真模型（static simulation model）有时也称为蒙特卡罗仿真（Monte Carlo simulation），表示系统在某一特定时间点的状态。动态仿真模型（dynamic simulation model）表示系统随时间的变化。例如，对银行从上午 9:00 到下午 4:00 的仿真就是一个动态仿真的例子。

不含随机变量的仿真模型被称为确定性（deterministic）的。确定性模型有一组已知的输入，从而产生一组唯一确定的输出。如果牙科诊所的所有病人都按预约时间准时到达，那么到达就是确定性的。随机仿真模型（stochastic simulation model）则含有一个或多个随机变量作为输入。随机输入导致随机输出。由于输出是随机的，因此只能被看作模型真实特性的估计值。对银行的仿真通常涉及随机的到达间隔时间和随机的服务时间。因此，在随机仿真中，输出度量——平均等待人数、顾客平均等待时间——必须被当作系统真实特性的统计估计来处理。

离散系统与连续系统已在第 7 节中定义。离散模型与连续模型的定义方式与此类似。然而，离散仿真模型并不总是用来模拟离散系统，连续仿真模型也不总是用来模拟连续系统。储罐和管道尽管我们知道流体流动是连续的，也可能被离散地建模。此外，仿真模型可以是混合型的，同时兼具离散与连续成分。究竟采用离散、连续（或者离散与连续兼有）的仿真模型，取决于系统的特征和研究目标。因此，如果每条消息的特征和移动都显得重要，通信信道就可以离散地建模；反之，若关注的是消息在信道上的总流量，那么用连续仿真来建模该系统可能更合适。本书所重点讨论的模型是动态的、随机的和离散的。

11 离散事件系统仿真

本书讨论的是离散事件系统仿真（discrete-event system simulation）。离散事件系统仿真是对状态变量仅在一组离散时间点上发生变化的系统进行建模。仿真模型采用数值方法而非解析方法进行分析。解析方法（analytical method）运用数学的演绎推理来"求解"模型。例如，可以使用微分学计算某些库存模型的最小成本策略。数值方法（numerical method）则运用计算过程来"求解"数学模型。对于采用数值方法的仿真模型而言，模型是被"运行"而非被"求解"的——也就是说，依据模型假设生成系统的一段人造历史，并收集观察值以供分析，从而估计真实系统的性能度量。真实世界的仿真模型相当庞大，所存储和操作的数据量巨大，因此这类运行通常要借助计算机来完成。

概括地说，本书讨论的是离散事件系统仿真，其中所关注的模型采用数值方法进行分析，通常借助计算机完成。

12 仿真研究的步骤

图 3 给出了一组步骤，用以指导建模者开展周密而可靠的仿真研究。其他文献中也有类似的图示与关于步骤的讨论 [Shannon, 1975; Gordon, 1978; Law, 2007]。图 3 中各符号旁的编号对应正文中更为详细的讨论。仿真研究的步骤如下：

图 3 仿真研究的步骤。

1. 问题表述（Problem formulation）

每项研究都应当以问题的陈述开始。如果问题陈述由决策者或遇到问题的人给出，分析人员必须确保对所描述的问题有清晰的理解。如果问题陈述由分析人员拟定，那么让决策者理解并认可这一表述就很重要。尽管图 3 中未予显示，但有时随着研究的推进，问题必须重新表述。在许多情形下，决策者和分析人员早在问题的性质尚未明确之前，就已意识到问题的存在。

2. 确定目标与总体项目计划（Setting of objectives and overall project plan）

目标指明仿真要回答的问题。在这一点上，应当判断：对于所表述的问题和所陈述的目标而言，仿真是否为恰当的方法学。假定已确定仿真是恰当的，那么总体项目计划应当包括对待考察备选系统的说明，以及评估这些备选系统有效性的方法。它还应当包括研究计划：涉及人员数量、研究成本、完成各阶段工作所需的天数，以及每一阶段结束时期望得到的结果。

3. 模型概念化（Model conceptualization）

构建系统模型大概是艺术与科学的结合。Pritsker [1998] 对这一步做了详尽的讨论。“尽管不可能给出一套指令，以确保在每一种情况下都能构建出成功而恰当的模型，但仍有一些可以遵循的一般性准则” [Morris, 1967]。建模艺术因以下能力而提升：抽象问题的本质特征，选择并修正刻画系统的基本假设，然后丰富和细化模型，直至得到一个有用的近似。因此，最好从一个简单的模型开始，逐步构建到更高的复杂度。不过，模型的复杂度不必超过达成模型既定目的所需的程度。违反这一原则只会增加建模和计算机方面的开销。模型与真实系统之间并不需要一一对应，只需要真实系统的精髓即可。

建议让模型使用方参与模型概念化工作。让模型使用方参与，既能提升所建模型的质量，也能增强模型使用方对模型应用的信心。（不过，只有对真实系统——而非教科书习题——的亲身经验，才能"教会"人们建模的艺术。）

4. 数据收集（Data collection）

模型的构建与所需输入数据的收集之间始终存在相互影响 [Shannon, 1975]。随着模型复杂度发生变化，所需的数据元素也会随之变化。此外，由于数据收集在完成一次仿真所需的总时间中占有很大比重，因此必须尽早开始，通常要与建模的早期阶段同步进行。

研究目标在很大程度上决定了要收集哪些数据。例如在银行研究中，若希望了解柜台员人数变化时排队长度的变化，那么所需的数据类型包括：到达间隔时间的分布（按一天中不同时段分别统计）、柜台员的服务时间分布，以及不同条件下排队长度的历史分布。最后一类数据将用于验证仿真模型。（参见 Henderson [2003] 的精彩讨论。）

5. 模型翻译（Model translation）

大多数真实世界的系统所导出的模型都需要大量的信息存储和计算，因此必须把模型录入计算机可识别的格式中。我们使用"程序"这一术语，尽管在许多情形下，或许几乎不需要真正的编码就能达成预期结果。建模者必须决定是用仿真语言（如 GPSS/H™）为模型编写程序，还是使用专用的仿真软件。在制造和物料搬运领域，有 AnyLogic®、Arena®、AutoMod™、Enterprise Dynamics®、Extend™、Flexsim、ProModel® 和 SIMUL8® 等软件。仿真语言强大而灵活。然而，若问题可用仿真软件解决，模型开发时间将大幅缩短。此外，多数仿真软件包都增加了增强其灵活性的功能，尽管灵活性的差异很大。

6. 验证了吗？（Verified?）

验证（verification）针对的是为仿真模型所编制的计算机程序。计算机程序的运行是否正确？对于复杂模型，要在不做大量调试的情况下完整成功地翻译模型是困难的，甚至几乎不可能；如果模型的输入参数和逻辑结构在计算机中得到了正确表示，验证即告完成。就大部分情况而言，完成这一步所依据的是常识。（Sargent [2007] 对此有专门讨论。）

7. 确认了吗？（Validated?）

确认（validation）通常通过模型的校准（calibration）来实现——这是一个迭代过程：把模型与实际系统行为作比较，利用两者之间的差异以及由此获得的见解来改进模型。反复进行这一过程，直到模型精度被认为可以接受。在前面提到的银行例子中，我们收集了当前条件下排队长度的数据。仿真模型能否复现这一系统度量？这是确认的一种手段。（Sargent [2007] 对此有专门讨论。）

8. 实验设计（Experimental design）

必须确定要仿真的备选方案。通常，关于仿真哪些备选方案的决策，会取决于已完成并分析过的运行结果。对每一个被仿真的系统设计，都需要就以下事项作出决策：初始化期的长度、仿真运行的长度，以及每次运行所要做的重复次数。（Sanchez [2007] 对此有非常详尽的讨论。）

9. 生产性运行与分析（Production runs and analysis）

生产性运行及其后续分析，用于估计被仿真系统设计的性能度量。（可辅助完成这一步的软件包括 AutoStat（在 AutoMod 中）、OptQuest（在若干仿真软件中）以及 SimRunner（在 ProModel 中）。）

10. 还需要更多运行吗？（More runs?）

在分析了已完成的运行之后，分析人员据此判断是否需要追加运行，以及这些追加实验应遵循何种设计。

11. 文档与报告（Documentation and reporting）

文档分为两类：程序文档与进展文档。程序文档（program documentation）的必要性体现在许多方面。如果该程序将由同一批或不同的分析人员再次使用，那么理解程序如何运行可能是必要的。这将建立对程序的信心，使模型使用方和决策者能够基于分析结果作出决策。同样，如果程序将由同一名或另一名分析人员修改，充分的文档会极大地方便这一工作。分析人员一次因文档不充分而吃了苦头，通常就足以让分析人员确信这一重要步骤的必要性。为程序编写文档还有一个理由：模型使用方可以随时更改参数，从而了解输入参数与输出性能度量之间的关系，或发现能"优化"某项输出性能度量的输入参数。

Musselman [1998] 讨论了进展报告（progress report），它提供仿真项目重要的书面历史。项目报告给出已完成工作和已作决策的时间脉络，这极有助于使项目保持在正确轨道上。Musselman 建议频繁出具报告（至少每月一次），使得即便那些未参与日常工作的人也能及时了解情况。他人的这种了解常常能促进项目的成功完成，因为能及早暴露误解，而此时问题尚易解决。Musselman 还建议维护一份项目日志（project log），以全面记录工作成果、变更请求、关键决策以及其他重要事项。

在交付物方面，Musselman 建议频繁产出阶段性成果。这些成果未必都是重大成就的产物。他的座右铭是：“与其只有一个绝对的截止期限，宁可设置许多中间里程碑。“最终报告之前的可能产出包括模型说明书、原型演示、动画、培训结果、中期分析、程序文档、进展报告和演示汇报。他建议在项目生命周期内对这些交付物作出合理的时间安排。

所有分析结果都应当在最终报告中清楚而简洁地报告。这将使模型使用方（此时即决策者）能够复核最终的问题表述、所涉及的备选系统、比较这些备选系统所依据的准则、实验结果，以及针对该问题所推荐的解决方案。此外，如果决策需要在更高层级上得到论证，最终报告还应当为模型使用方／决策者提供一种认证载体，并提升模型及建模过程的可信度。

12. 实施（Implementation）

实施阶段能否成功，取决于前十一个步骤执行得如何好，也取决于分析人员在整个仿真过程中让最终模型使用方参与的充分程度。如果模型使用方在整个建模过程中都参与其中，并且理解模型及其输出的性质，那么实施有力推进的可能性就会增大 [Pritsker, 1995]。反之，如果模型及其背后的假设没有得到恰当的沟通，那么无论仿真模型本身是否有效，实施都可能受到影响。

图 3 所示的仿真建模过程可以分解为四个阶段。第一阶段包括步骤 1（问题表述）和步骤 2（确定目标与总体设计），是一段发现或定向时期。最初的问题陈述通常相当"模糊”，最初的目标通常不得不重新设定，最初的项目计划通常也必须微调。这些重新校准与澄清可能发生在本阶段内，也可能在另一个阶段之后或进行之中发生（例如，分析人员可能不得不重新启动整个过程）。

第二阶段与模型构建和数据收集相关，包括步骤 3（模型概念化）、步骤 4（数据收集）、步骤 5（模型翻译）、步骤 6（验证）和步骤 7（确认）。各步骤之间需要持续相互配合。在此阶段将模型使用方排除在外，可能会在实施时带来严重后果。

第三阶段涉及模型的运行，包括步骤 8（实验设计）、步骤 9（生产性运行与分析）和步骤 10（是否需要更多运行）。本阶段必须针对仿真模型制定一个周密的实验计划。事实上，离散事件随机仿真就是一次统计实验。输出变量是含有随机误差的估计值，因此需要进行恰当的统计分析。这种思想与那种只做一次运行、就从这一个数据点得出推断的分析人员形成了鲜明对照。

第四阶段即实施阶段，涉及步骤 11（文档与报告）和步骤 12（实施）。能否成功实施取决于模型使用方的持续参与，以及流程中每一步的顺利完成。整个流程中也许最关键的一点是步骤 7（确认），因为一个未通过确认的模型将导致错误的结果，而这些结果一旦付诸实施，可能既危险又昂贵，甚至两者兼有。

参考文献

BANKS, J. [2008], “Some Burning Questions about Simulation.” ICS Newsletter, INFORMS Computing Society, Spring.

BANKS, J., and R. R. GIBSON [1997], “Don’t Simulate When: 10 Rules for Determining when Simulation Is Not Appropriate,” IIE Solutions, September.

GORDON, G. [1978], System Simulation, 2d ed., Prentice-Hall, Englewood Cliffs, NJ.

HENDERSON, S. G. [2003], “Input Model Uncertainty: Why Do We Care and What Should We Do About It?” in Proceedings of the Winter Simulation Conference, eds. S. Chick, P. J. Sánchez, D. Ferrin, and D. J. Morrice, New Orleans, LA, Dec. 7–10, pp. 90–100.

KLEIJNEN, J. P. C. [1998], “Experimental Design for Sensitivity Analysis, Optimization, and Validation of Simulation Models,” in Handbook of Simulation, ed. J. Banks, John Wiley, New York.

LAW, A. M. [2007], Simulation Modeling and Analysis, 4th ed., McGraw–Hill, New York.

MORRIS, W. T. [1967], “On the Art of Modeling,” Management Science, Vol. 13, No. 12.

MUSSELMAN, K. J. [1998], “Guidelines for Success,” in Handbook of Simulation, ed. J. Banks, John Wiley, New York.

NAYLOR, T. H., J. L. BALINTFY, D. S. BURDICK, and K. CHU [1966], Computer Simulation Techniques, Wiley, New York.

PEGDEN, C. D., R. E. SHANNON, and R. P. SADOWSKI [1995], Introduction to Simulation Using SIMAN, 2d ed., McGraw–Hill, New York.

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

PRITSKER, A. A. B. [1998], “Principles of Simulation Modeling,” in Handbook of Simulation, ed. J. Banks, John Wiley, New York.

SANCHEZ, S.R. [2007], “Work Smarter, Not Harder: Guidelines for Designing Simulation Experiments,” in Proceedings of the 2007 Winter Simulation Conference, eds. S. G. Henderson, B. Biller, M.-H. Hsieh, J. Shortle, J. D. Tew, and R. R. Barton, Washington, DC, Dec. 9–12, pp. 84–94.

SARGENT, R.G. [2007], “Verification and Validation of Simulation Models,” in Proceedings of the 2007 Winter Simulation Conference, eds. S. G. Henderson, B. Biller, M.-H. Hsieh, J. Shortle, J. D. Tew, and R. R. Barton, Washington, DC, Dec. 9–12, pp. 124–137.

SHANNON, R. E. [1975], Systems Simulation: The Art and Science, Prentice-Hall, Englewood Cliffs, NJ.

SHANNON, R. E. [1998], “Introduction to the Art and Science of Simulation,” in Proceedings of the Winter Simulation Conference, eds. D. J. Medeiros, E. F. Watson, J. S. Carson, and M. S. Manivannan, Washington, DC, Dec. 13–16, pp. 7–14.

习题

1. 为下列系统各列举若干实体、属性、活动、事件和状态变量：

   • (a) 一家自助餐厅
   • (b) 一家杂货店
   • (c) 一家自助洗衣店
   • (d) 一家快餐店
   • (e) 一家医院急诊室
   • (f) 一家拥有 10 辆出租车的出租车公司
   • (g) 一条汽车装配线

2. 考虑图 3 所示的仿真过程。

   • (a) 通过合并相似的活动，至少减少两个步骤。请给出理由。
   • (b) 通过拆分现有步骤或对现有步骤加以细化，至少增加两个步骤。请给出理由。

3. 拟对某大型交通路口进行仿真，目的是改善当前的交通流。请针对图 3 仿真过程中的步骤 1 和步骤 2，按复杂度递增的顺序给出三个迭代版本。

4. 拟对烹制一顿意大利面进行仿真，以确定一个人应当在什么时间开始动手，才能在晚上 7:00 前把饭菜端上餐桌。请从网上下载一份意大利面的食谱（或向朋友或亲戚索取食谱）。尽可能地追踪你认为在图 3 仿真过程的数据收集阶段所需要进行的内容，以便完成一次包含食谱中每个步骤的仿真。该系统中的事件、活动和状态变量分别是什么？

5. 使用你的支票簿所涉及的事件和活动有哪些？

6. 第 1 节中给出了大量使用仿真的理由。但可能还有其他理由。请查阅当届 WSC 论文集中的应用，看能否发现其他理由。（WSC 论文集见 www.informs-cs.org/wscpapers.html）

7. 在 www.informs-cs.org/wscpapers.html 可获取的当届 WSC 论文集中，选取一篇与你所学专业或兴趣领域相关的仿真应用文章，撰写报告，说明作者是如何完成图 3 所列各步骤的。

8. 从习题 7 的网址获取一份近期的 WSC 论文集，就你感兴趣的领域，报告其中讨论的不同应用。

9. 从习题 7 的网址获取一份近期的 WSC 论文集，报告你所发现的最不寻常的应用。

10. 找出第 12 节中的某一个案例，描述其所采用的方法。

11. 访问冬季仿真会议网站 www.wintersim.org，回答下列问题：

    • (a) 上一届 WSC 提供了哪些高阶教程，或下一届 WSC 计划提供哪些高阶教程？
    • (b) 下一届 WSC 将于何时何地举办？

12. 访问冬季仿真会议网站 www.wintersim.org，回答下列问题：

    • (a) 迄今规模最大（按参会人数计）的一届 WSC 是哪一届，共有多少人参加？
    • (b) 自 WSC 创办以来，哪一个公历年份没有召开会议？
    • (c) 自 WSC 创办以来，两次会议之间最长的时间间隔是多久？
    • (d) 从第 25 届 WSC 开始，你能否看出会议举办地点的某种规律？

13. 向 WSC 提交论文的流程是怎样的？

14. 将近期 WSC 的参展商名单按"软件供应商"和"其他"两类进行划分。

15. WSC 基金会的宗旨与历史是什么？

16. 使用你常用的搜索引擎，在网上搜索 ‘discrete event simulation output analysis’，并就搜索结果撰写报告。

17. 使用你常用的搜索引擎，在网上搜索 ‘supply chain simulation’，并就搜索结果撰写报告。

18. 使用你常用的搜索引擎，在网上搜索 ‘web based simulation’，并就搜索结果撰写报告。

2 电子表格中的仿真示例

本章借助电子表格示例来介绍仿真（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 之间的随机值，如下所示：

或者用于生成两个值之一（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 掷硬币的仿真表

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

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

每一次调用 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 服务时间的输入规格说明与仿真表

图 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 从随机数到服务时间的"随机飞镖"变换。

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

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

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

其中 $ 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()，它带两个参数：或者先给出含期望整数区间最小值的单元格、再给出含最大值的单元格，或者直接给出实际的上下限，例如：

或

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

对本例而言即为：

其中 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 到达时间的输入规格说明与仿真表

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

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

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

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

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

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

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

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

3. 一个系统状态；

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

5. 一个模型输出；

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

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

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

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

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

表 4 通用仿真表

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

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

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 掷硬币游戏的输入规格说明与仿真表

如果你查看电子表格解答中的「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 掷硬币游戏的实验工作表

表 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 领先机会的直方图。

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

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

3 电子表格中的排队仿真

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

图 3 单通道排队系统。

3.1 等待线模型

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

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

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

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

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

图 4 离开事件流程图。

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

图 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 强调时钟时间的仿真表

表 8 事件的时间顺序

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

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

图 8 系统中的顾客数。

3.2 单服务台队列的仿真

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

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

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

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

表 9 到达之间时间的分布

表 10 服务时间的分布

如何由服务时间分布与到达间隔时间分布生成服务时间和到达时间，已分别在第 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 位顾客）

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

图 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 次试验平均等待时间的频率分布。

3.3 双服务台队列的仿真

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

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

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

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

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

表 13 Able 的服务时间分布

表 14 Baker 的服务时间分布

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

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

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

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

图 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 次试验实验中平均来电者延迟的频率。

接下来我们要问的问题是：平均来电者延迟（前 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 概率型订购至水平库存系统。

为简单起见，图 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 按新闻日类型划分的每日报纸需求量分布

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

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

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

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

表 17 新闻日类型的分布

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

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

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

其余各天的利润用同样的方法即可容易地算出。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 天总利润的频率。

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

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

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

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

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

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

例 8：冰箱库存问题

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

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

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

表 19 每日需求量的分布

表 20 提前期的分布

“订购至 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 台冰箱的库存策略仿真表

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

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

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

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

在「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 轴承寿命的分布

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

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

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

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

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

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

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

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

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

轴承的总寿命为 $( 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 个炸弹落点的弹药库。

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

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

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

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

表 26 模拟的一次轰炸

下标 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()) 来生成具有所需标准差的正态变量，如下所示：

用于 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 轰炸任务命中数量的直方图。

5.3 提前期需求分布的估计

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

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

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

例 11：提前期需求

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

表 27 提前期需求的仿真表

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

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

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

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

图 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() 变换到任何指定范围，如下：

' 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 表示活动的先后顺序，每条路径代表一个人的任务，如下：

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 的活动网络。

由于连续均匀分布的均值是下限与上限的中点，可以看到图 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 次试验的项目完成时间直方图。

项目完成时间被绘制在频率图（直方图）中，如图 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. 考虑如下连续运转的作业车间。作业的到达间隔时间分布如下：

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

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

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

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

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

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

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

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

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

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

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

5. 随机变量 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 的分布如下表所示：

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

使用样本量为 100 的样本。

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

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

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

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

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

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

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

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

   (h) 让仿真运行 5 周。

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

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

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

(a) 对随机变量

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

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

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

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$。则随机变量

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

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

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

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

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

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

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

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

用仿真估计该湖的面积。（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 中，把服务时间分布改为如下：

运行 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 设计一个蒙特卡罗实验来估计 π 的值，如同圆面积的公式

提示：考虑单位正方形内一个半径为 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 老虎机，其每个转轴上的符号数量如下：

中奖赔付如下：

使用 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 文件

3 一般原则

本章通过使用离散事件仿真（discrete-event simulation），为复杂系统的建模建立一个共同的框架。本章涵盖所有离散事件仿真模型的基本构件：实体（entity）与属性（attribute）、活动（activity）以及事件（event）。在离散事件仿真中，系统是依据其在每个时点的状态来建模的；依据流经系统的实体以及代表系统资源的实体来建模；依据引起系统状态改变的活动和事件来建模。对于那些系统状态仅在离散时间点上发生变化的系统，离散事件模型是适合的。

本章的讨论聚焦于仿真语言与仿真软件（统称为仿真软件包，simulation package）的离散事件概念与方法学。

本章介绍并解释所有离散事件仿真包所蕴含的基本概念与方法学。这些概念与方法学并不依附于任何特定的软件包。许多软件包使用的术语与本章所用的不同，而且大多数软件包都配备了若干高阶构件，以便在各自的应用领域中使建模更为简单直接。例如，本章讨论实体的基本抽象概念，但不涉及某些软件包中内置的、更为具体的实体实现（如机器、输送带和车辆），这些构件是为了便利制造、物料搬运或其他领域的建模而设置的。

第 1 节讨论离散事件仿真的一般原则与概念、事件调度/时间推进算法（event-scheduling/time-advance algorithm），以及三种流行的世界观（world view）：事件调度（event scheduling）、过程交互（process interaction）和活动扫描（activity scanning）。第 2 节介绍列表处理（list processing）的一些概念，列表处理是离散事件仿真软件中更为重要的方法学之一。

1 离散事件仿真中的概念

本章专门讨论动态的、随机的（即涉及时间且包含随机成分的）并且以离散方式变化的系统。本节将展开讨论这些概念，并提出一个用于建立系统离散事件模型的框架。主要概念先给出简要定义，随后以示例加以说明：

系统（system）　一组随时间相互交互、从而实现一个或多个目标的实体（例如人和机器）的集合。

模型（model）　系统的一种抽象表示，通常包含用以状态、实体及其属性、集合、过程、事件、活动和延迟来描述系统的结构性、逻辑性或数学性关系。

系统状态（system state）　一组变量，包含了在任意时刻描述系统所需的全部信息。

实体（entity）　系统中需要在模型中显式表示的任何对象或构件（例如，一台服务台、一位顾客、一台机器）。

属性（attributes）　某一给定实体的性质（例如，等待顾客的优先级，某个工件在作业车间中的加工路线）。

列表（list）　一组（永久或临时）关联起来的实体，按某种逻辑方式排序（例如，当前所有处于排队线上的顾客，按"先到先服务"或按优先级排序）。

事件（event）　改变系统状态的一次瞬时发生（例如，一位新顾客的到达）。

事件通知（event notice）　一条关于将在当前或某个未来时刻发生的事件的记录，以及执行该事件所需的一切相关数据；该记录至少应包括事件类型和事件时间。

事件列表（event list）　未来事件的事件通知所组成的列表，按发生时间排序；也称为未来事件列表（future event list, FEL）。

活动（activity）　一段规定长度的时间（例如服务时间或到达间隔时间），其长度在开始时即已知（尽管它可能是以统计分布的形式定义的）。

延迟（delay）　一段长度不确定、无法预先确定的时间，在其结束之前无法知道其长度（例如，顾客在后进先出排队线中的延迟，在延迟开始时还取决于未来的到达）。

时钟（clock）　表示仿真时间的变量，在后续示例中称之为 CLOCK。

不同的仿真软件包对相同或相似的概念使用不同的术语——例如，列表有时被称为集合（set）、队列（queue）或链（chain）。集合或列表用于保存实体和事件通知。列表中的实体总是按某种规则排序，例如先入先出或后入先出，或者按某种实体属性（如优先级或到期日）排名。未来事件列表则总是按事件通知中记录的事件时间排名。第 2 节讨论处理列表的若干方法，并介绍一些对有序集合或列表进行高效处理的方法学。

活动通常代表服务时间、到达间隔时间，或其他任何持续时间已经由建模者刻画并定义的处理时间。活动持续时间的指定方式有若干种：

a. 确定性——例如，总是恰好 5 分钟；

b. 统计性——例如，以等概率从集合 2、5、7 中随机抽取；

c. 依赖于系统变量和/或实体属性的函数——例如，一艘铁矿石船的装载时间是该船允许载重量与每小时装载吨位装载率的函数。

无论以何种方式刻画，活动持续时间在其开始的那一瞬间都可由其规定算出。其持续时间不受其他事件发生的影响（除非如某些仿真软件包所允许的那样，模型中包含取消或推迟正在进行的活动的逻辑）。为跟踪活动及其预计完成时间，在仿真时刻某一活动持续时间开始时，生成一条事件通知，其事件时间等于该活动的完成时间。例如，若当前仿真时间为 CLOCK 100 分钟，且一项恰好 5 分钟的检验时间刚刚开始，则生成一条事件通知，规定事件类型（检验结束事件）与事件时间（100 + 5 = 105 分钟）。

与活动相反，延迟的持续时间并非由建模者预先指定，而是由系统状况决定。延迟的持续时间常常需要测量，并且是模型运行所期望的输出之一。通常，延迟在某组逻辑条件成立或一个或多个其他事件发生时结束。例如，顾客在排队线中的延迟可能取决于排在前面其他顾客的人数与服务持续时间，也可能取决于服务台与设备的可用性。

延迟有时称为条件等待，活动称为无条件等待。活动的完成是一个事件，常称为主事件（primary event），通过在 FEL 上放置一条事件通知来管理。与之相对，延迟则通过把相关实体放置到另一个列表（或许代表一条排队线）上来管理，直至系统状况允许处理该实体为止。延迟的完成有时称为条件事件或次事件（conditional or secondary event），但这类事件并不用事件通知表示，也不会出现在 FEL 上。

这里所考虑的系统都是动态的，即随时间变化的。因此，系统状态、实体属性与在活动实体数、集合的内容，以及当前正在进行中的活动和延迟，都是时间的函数，并随时间不断变化。时间本身由一个称为 CLOCK 的变量表示。

例 1：呼叫中心，再度审视

考虑《电子表格中的仿真示例》一章例 6 中的 Able–Baker 呼叫中心系统。一个离散事件模型具有以下构件：

系统状态 $L_{Q}(t)$，时刻 $t$ 正在等待接受服务的呼叫者数量；$L_{A}(t)$，取 0 或 1，指示时刻 $t$ Able 是空闲还是忙碌；$L_{B}(t)$，取 0 或 1，指示时刻 $t$ Baker 是空闲还是忙碌。

实体　除非希望得到某些呼叫者的平均值（参见例 4 和例 5），否则呼叫者和服务台都不需要显式表示，只能通过状态变量来体现。

事件 到达事件； Able 的服务完成； Baker 的服务完成。

活动 到达间隔时间，定义见《电子表格中的仿真示例》一章的表 12； Able 的服务时间，定义见《电子表格中的仿真示例》一章的表 13； Baker 的服务时间，定义见《电子表格中的仿真示例》一章的表 14。

延迟 呼叫者在队列中等待，直到 Able 或 Baker 空闲下来。

模型构件的定义给出了模型的静态描述。此外，还需要对这些构件之间的动态关系与交互作用加以描述。需要得到回答的问题包括：

1. 每个事件如何影响系统状态、实体属性和集合内容？

2. 活动是如何定义的（即，确定性的、概率性的，还是某个其他数学方程）？什么事件标志着每项活动的开始或结束？该活动能否不受系统状态限制就开始，还是它的开始以系统处于某种状态为条件？（例如，一项机加工"活动"只有在机床空闲、未发生故障且不在维修之中的情况下才能开始。）

3. 哪些事件触发每类延迟的开始（与结束）？延迟在什么条件下开始或结束？

4. 时刻 0 的系统状态是什么？应当在时刻 0 生成哪些事件来"预热"模型——也就是使仿真得以启动？

离散事件仿真，是对这样一个系统随时间进行的建模：该系统的全部状态变化都发生在离散的时间点上——即事件发生的那些时点。离散事件仿真（下文简称仿真）通过产生一系列系统快照（system snapshot，或称系统映像）来推进，这些快照表示系统随时间的演化。给定时刻的一个快照（CLOCK t）不仅包含时刻 t 的系统状态，还包含当前正在进行的所有活动及其各自结束时刻的列表（FEL）、所有实体的状态与所有集合的当前成员，以及用于在仿真结束时计算汇总统计量的累积统计量与计数器的当前值。图 1 给出了一个原型系统快照。（并非所有模型都会包含图 1 中所展示的每一个要素。本章各示例中提供了进一步的说明。）

图 1 仿真时刻 t 的原型系统快照。

1.1 事件调度/时间推进算法

推进仿真时间、并保证所有事件按正确的时序发生的机制，建立在未来事件列表（future event list, FEL）的基础之上。该列表包含所有已安排在未来某个时刻发生的事件的事件通知。安排一个未来事件，意味着在某一活动开始的那一刻计算其持续时间，或从某个统计分布中抽取一个样本；并把该活动结束事件连同其事件时间一起放到未来事件列表上。在现实世界中，大多数未来事件并非事先安排，而只是偶然发生——例如随机故障或随机到达。在模型中，这类随机事件由某项活动的结束来表示，而活动的结束又由一个统计分布来表示。

在任一给定时刻 t，FEL 包含所有先前已安排的未来事件及其相应的事件时间（在图 1 中记作 $t_{1} , t_{2} , \ldots$）。FEL 按事件时间排序，意味着事件按时间先后排列——即各事件时间满足

时刻 t 是 CLOCK 的值，即当前仿真时间。与时刻 $t_{1}$ 对应的事件称为即时事件（imminent event）；也就是说，它是接下来将要发生的事件。在仿真时刻 CLOCK t 的系统快照更新之后，CLOCK 被推进到仿真时刻 ${\mathrm{C L O C K}} = t_{1}$，即时事件的事件通知从 FEL 中被移除，然后执行该事件。执行即时事件意味着：基于时刻 t 的旧快照以及即时事件的性质，为时刻 $t_{1}$ 生成一个新的系统快照。在时刻 $t_{1} .$，可能会生成新的未来事件，也可能不会；但只要生成了，就要通过创建事件通知并把其放入 FEL 上的适当位置来加以安排。在时刻 $t_{1}$ 的新系统快照更新之后，时钟被推进到新的即时事件发生的时刻，并执行该事件。这一过程不断重复，直至仿真结束。仿真器（或仿真语言）为推进时钟并构建新系统快照而必须执行的一系列动作，称为事件调度/时间推进算法，其步骤列于图 2（并在其后加以解释）。

随着仿真的推进，FEL 的长度与内容不断变化，因此在仿真软件包中对其进行高效管理，将对模型的计算机运行时间产生重大影响。对一个列表的管理称为列表处理（list processing）。在 FEL 上执行的主要列表处理操作是移除即时事件、向列表中添加新事件，以及偶尔移除某个事件（称为事件的取消）。由于即时事件通常位于列表的顶部，移除它的效率是最高的。添加新事件（以及取消旧事件）则需要搜索该列表。这一搜索的效率取决于列表的逻辑组织方式以及搜索的实施方式。除 FEL 之外，模型中的所有集合也按某种逻辑顺序维护，向集合中添加和移除实体同样需要高效的列表处理技术。第 2 节将简要介绍仿真中的列表处理。

FEL 中事件的移除与添加在图 2 中得到说明。事件时间为 $t_{1}$ 的事件 3，比如说是服务台 3 的一个服务完成事件。由于它是时刻 t 的即时事件，因此它在事件调度/时间推进算法的第 1 步（图 2）中被从 FEL 移除。当事件 4（比如一个到达事件）在第 4 步被生成、其事件时间为 $t^{*}$ 时，确定其在 FEL 上正确位置的一种可能做法是进行自顶向下的搜索：

（在图 2 中，假定 $t^{*}$ 位于 $t_{2}$ 与 $t_{3}$ 之间。）另一种做法是进行自底向上的搜索。维护 FEL 效率最低的方式是把它留作一个无序列表（新增项任意放置在顶部或底部），这将要求在图 2 的第 1 步中，在每次推进时钟之前对整个列表进行完整搜索以找出即时事件。（即时事件是 FEL 上事件时间最低的那个事件。）

时刻 0 的系统快照由初始条件以及所谓外生事件（exogenous event）的生成所确定。所指定的初始条件确定了时刻 0 的系统状态。例如，在图 2 中，若 $t = 0$，则状态 (5, 1, 6) 可以代表系统中三个不同位置上初始的顾客数。外生事件（exogenous event）是"发生在系统之外"、并对系统产生影响的一次发生。一个重要的例子是到达排队系统的事件。在时刻 0，生成第一个到达事件并将其安排到 FEL 上（也就是说，把其事件通知放到 FEL 上）。到达间隔时间是活动的一个例子。当时钟最终被推进到这一首次到达的时刻时，生成第二个到达事件。首先，生成一个到达间隔时间，记为 $a^{*} .$

时刻 t 的旧系统快照

事件调度/时间推进算法

第 1 步. 从 FEL 中移除即时事件的事件通知（事件 3，时刻 $t_{1} )$。

第 2 步. 把 CLOCK 推进到即时事件时刻（即把 CLOCK 从 t 推进到 $t_{1}$）。

第 3 步. 执行即时事件：更新系统状态，按需改变实体属性与集合成员关系。

第 4 步. 生成未来事件（若有必要），并把其事件通知按事件时间排序后放到 FEL 上。（例如：在时刻 $t^{*}$ 发生的事件 4，其中 $t_{2} < t^{*} < t_{3} . )$。

第 5 步. 更新累积统计量与计数器。

时刻 $t_{1}$ 的新系统快照

图 2 推进仿真时间并更新系统映像。

并把它加到当前时间上，${\mathrm{C L O C K}} = t ;$ 由此得到的（未来）事件时间 $t + a^{*} = t^{*}$ 用于确定新的到达事件通知在 FEL 上的位置。这种生成外部到达流的方法称为自举（bootstrapping）；它给出了事件调度/时间推进算法第 4 步中如何生成未来事件的一个例子。自举过程在图 3 中得到说明。

图 3 通过自举生成外部到达流。

生成的前三个到达间隔时间为 3.7、0.4 和 3.3 个时间单位。到达间隔区间的结束是主事件的一个例子。

未来事件如何生成（图 2 第 4 步）的第二个例子，是排队仿真中的服务完成事件。当一位顾客在当前时刻 $\mathrm{C L O C K} = t ,$ 完成服务时，若下一位顾客已在，则要为下一位顾客生成一个新的服务时间 $s^{*}$。下一个服务完成事件将被安排在未来时刻 $t^{*} = t + s^{*}$ 发生，做法是在 FEL 上放入一条新的、类型为服务完成、事件时间为 t∗ 的事件通知。此外，若在到达事件发生的时刻，服务台组中至少有一台空闲的服务台，则也将生成一个服务完成事件并加以安排。服务时间是活动的一个例子。开始服务是条件事件（conditional event），因为其发生仅在"有顾客在场且有服务台空闲"这一条件下才被触发。服务完成是主事件的一个例子。注意，像开始服务这样的条件事件，是由某个主事件的发生以及系统中盛行的某些条件共同触发的。只有主事件才会出现在 FEL 上。

未来事件如何生成的第三个重要例子，是一台可能发生故障的机器其运行时长与停机时长的交替生成。在时刻 0，生成第一个运行时长，并安排一个运行时长结束事件。每当运行时长结束事件发生时，生成一个停机时长，并在 FEL 上安排一个停机时长结束事件。当 CLOCK 最终被推进到这一停机时长结束事件的时刻时，生成一个运行时长，并在 FEL 上安排一个运行时长结束事件。以此方式，运行时长与停机时长在整个仿真过程中不断交替。运行时长和停机时长是活动的例子，而运行时长结束和停机时长结束是主事件。

每一个仿真都必须有一个停止事件，在此称为 E，它规定了仿真将运行多长时间。停止仿真通常有两种方式：

1. 在时刻 0，把一个停止仿真事件安排在某个指定的未来时刻 $T_{E} .$。这样，在仿真之前就知道仿真将在时间区间 $[ 0 , T_{E} ]$ 上运行。例如：仿真一个作业车间，$T_{E} = 40$ 小时。

2. 运行长度 $T_{E}$ 由仿真本身决定。一般而言，$T_{E}$ 是某个指定事件 E 的发生时刻。例如：$T_{E}$ 可以是某个服务中心第 100 次服务完成的时刻；$T_{E}$ 可以是一个复杂系统发生故障的时刻；$T_{E}$ 可以是战斗仿真中脱离接触或彻底毁伤（以先发生者为准）的时刻；$T_{E}$ 可以是配送中心把某一天订单中的最后一个纸箱发运出去的时刻。

在第 2 种情形下，$T_{E}$ 事先是不知道的。事实上，它完全可能就是仿真所要产生的主要统计量之一。

1.2 世界观

在使用仿真软件包时，甚至在做手工仿真时，建模者都会采取某种世界观或视角来开发模型。最为盛行的世界观是上一节讨论过的事件调度世界观、过程交互世界观和活动扫描世界观。即便某个特定的软件包并不直接支持其中一种或多种世界观，理解这些不同的方法仍可为某一给定系统的建模提供其他思路。

概括前面的讨论，在采用事件调度方法时，仿真分析人员关注的是事件及其对系统状态的影响。第 1.3 节的手工仿真将说明这一世界观。

在使用支持过程交互方法的软件包时，仿真分析人员以过程（process）的思维来思考问题。分析人员用实体或对象及其在流经系统过程中的生命期来定义仿真模型，这些实体或对象会要求资源并排队等待资源。更确切地说，一个过程就是一个实体的生命期。这一生命期由各种事件、活动构成。有些活动可能需要使用一项或多项容量有限的资源。正是这类资源约束和其他约束导致了过程之间的交互，最简单的例子就是：某个实体因其所需要的资源正被另一实体占用，而被迫在队列中（在一个列表上）等待。过程交互方法之所以流行，是因为它具有直观上的吸引力，而且实现该方法的仿真软件包允许分析人员用高层的模块或网络构件来描述过程流，而过程之间的交互则由软件自动处理。

更确切地说，一个过程是一个按时间顺序排列的事件、活动和延迟的列表（其中包含对资源的需求），它定义了一个实体在系统内移动时的生命期。图 4 给出了一个"顾客过程"的例子。在该图中，可以看到两个顾客过程之间的交互：顾客 n + 1 被延迟，直至前一位顾客的"服务结束事件"发生。通常，一个模型中会有许多过程同时处于活动状态，过程之间的交互可能相当复杂。

在仿真软件包中实现过程交互方法的基础（但通常对建模者不可见）是：事件被安排到未来事件列表上，实体在遇到延迟时被放到各个列表上，从而使一个过程暂时挂起其执行，而其他过程继续推进。建模者对相关概念有基本理解、并对所用仿真软件包中内置但隐藏的运行规则有详细理解，是十分重要的。Schriber 与 Brunner [2003] 对这方面提供了阐释。

事件调度方法与过程交互方法都采用可变时间推进——也就是说，当某一仿真时刻的所有事件与系统状态变化都已发生之后，仿真时钟被推进到 FEL 上下一个即时事件的时刻。与之相反，活动扫描方法采用固定时间增量，并用一套基于规则的方法来决定在每个仿真时点上是否有活动可以开始。

图 4 单服务台队列中两个相互作用的顾客过程。

在活动扫描方法下，建模者关注模型中的各项活动，以及那些使某项活动得以开始的条件（无论简单还是复杂）。在每次时钟推进时，检查每项活动的条件，若条件成立，则相应的活动开始。该方法的倡导者认为活动扫描方法概念简单，并能产生模块化的模型，这类模型日后更易于被其他分析人员维护、理解和修改。然而他们也承认，为查明某项活动能否开始而进行的反复扫描，会导致计算机上的运行速度缓慢。因此，纯活动扫描方法被所谓的三阶段方法（three-phase approach）加以修改（并在概念上变得略微复杂一些）；该方法把事件调度的若干特征与活动扫描结合起来，从而既允许可变时间推进，又能在不必要的时候避免扫描，同时保留活动扫描方法的主要优点。

在三阶段方法中，事件被视为持续时间为零个时间单位的活动。按照这一定义，活动被分为两类，称为 B 类和 C 类。

B 类活动　必然要发生的活动；即所有主事件和无条件活动。

C 类活动　以某些条件成立为前提的活动或事件。

B 类活动与事件可以像在事件调度方法中那样预先安排。这就允许了可变时间推进。FEL 中只包含 B 类事件。为查明是否有 C 类活动可以开始或 C 类事件发生而进行的扫描，只在每次时间推进结束、全部 B 类事件完成之后才进行。总之，在三阶段方法下，仿真通过反复执行以下 3 个阶段而推进，直至仿真完成：

阶段 A　从 FEL 中移除即时事件，并把时钟推进到其事件时刻。把 FEL 中所有具有相同事件时间的其他事件也一并移除。

阶段 B　执行所有从 FEL 中移除的 B 类事件。（这可能释放若干资源，或以其他方式改变系统状态。）

阶段 C　检查触发各项 C 类活动的条件，并激活其中条件已满足者。反复扫描，直至没有更多的 C 类活动可以开始、也没有事件发生。

三阶段方法提高了活动扫描方法的执行效率。此外，该方法的倡导者认为，活动扫描方法与三阶段方法特别擅长处理复杂的资源问题——在这些问题中，完成不同任务需要资源的各种不同组合。这两种方法保证：给定仿真时刻被释放的所有资源，都将在任何可用资源被重新分配给新任务之前得到释放。

例 2：呼叫中心，再度登场

对于例 1 中的呼叫中心，若采用面向活动扫描世界观的三阶段方法，则阶段 C 中各项活动的开始条件如下：

若采取过程交互世界观，则应从一位呼叫者及其"生命期"的视角来看待该模型。若认为生命期始于到达，则一个顾客过程如图 4 所示。

总之，过程交互方法已被美国最为流行的仿真软件包所采用。另一方面，若干活动扫描软件包则在英国和欧洲颇受欢迎。有些软件包允许模型的一部分采用基于事件调度的形式（若该视角更方便），再与过程交互方法混合使用。最后，还有一些软件包以流程图、框图或网络结构为基础，而仔细考察之下会发现，它们其实是过程交互概念的一种具体实现。

1.3 使用事件调度进行手工仿真

在进行事件调度仿真的过程中，需要使用一张仿真表（simulation table），以记录时间推进过程中的各个相继系统快照。

例 3：单通道队列

考虑《电子表格中的仿真示例》一章例 5 中以某种特殊方法仿真的、只有一个收银台的杂货店。系统由排队线上的顾客加上正在结账的那一位（若有）构成。本例设定停止时刻为 60 分钟。模型具有以下构件：

系统状态　(LQ(t), LS(t))，其中 LQ(t) 是排队线上的顾客数，LS(t) 是时刻 t 正在接受服务的顾客数（0 或 1）。

实体　除通过状态变量体现外，服务台和顾客都不作显式建模。

事件

到达（A）； 离开（D）； 停止事件（E），安排在时刻 60 发生。

事件通知

(A, t)，表示将在未来时刻 t 发生的到达事件； (D, t)，表示将在未来时刻 t 发生的顾客离开事件； (E, 60)，表示将在未来时刻 60 发生的仿真停止事件。

活动

到达间隔时间，定义见《电子表格中的仿真示例》一章的表 9； 服务时间，定义见《电子表格中的仿真示例》一章的表 10。

延迟　顾客在排队线中花费的时间。

事件通知写作（事件类型, 事件时间）。在本模型中，FEL 将始终包含两条或三条事件通知。到达事件与离开事件的执行效果在图 5 和图 6 中有详细说明。

图 5 到达事件的执行。

图 6 离开事件的执行。

收银台的仿真表见表 1。读者应遮住除一个之外的所有系统快照（从第一个开始），并尝试根据前一个快照以及图 5 和图 6 中的事件逻辑构建出下一个快照。

初始条件为：第一位顾客在时刻 0 到达并开始接受服务。这一点在表 1 中体现为时刻零 $( \mathbf{C L O C K} = 0 )$ 的系统快照，其中 $L Q ( 0 ) = 0 , L S ( 0 ) = 1$，且 FEL 上有一个离开事件和一个到达事件。此外，仿真安排在时刻 60 停止。只收集两项统计量：服务台利用率和最大队列长度。服务台利用率定义为服务台总忙碌时间（B）除以总时间 $( T_{E} )$。总忙碌时间 B 与最大队列长度 MQ 将随着仿真的推进而累积。表中另设"备注"一栏以帮助读者（$a^{*}$ 和 $s^{*}$ 分别是生成的到达间隔时间与服务时间）。

时刻 CLOCK 0 的系统快照一完成，仿真即开始。在时刻 0，即时事件为 $( \mathsf {A} , 1 )$。CLOCK 被推进到时刻 1，$( \mathsf {A} , 1 )$ 从 FEL 中移除。由于 $L S ( t ) = 1$ 对 $0 \leq t \leq 1 ( \mathrm{i . e .}$，服务台忙碌了 1 分钟），累积忙碌时间由 $B = 0$ 增加到 $B = 1$。依据图 6 中的事件逻辑，置 $L S ( 1 ) = 1$（服务台变为忙碌）。FEL 中还剩三个未来事件，(A, 2)、(D, 4) 和 (E, 60)。仿真 CLOCK 接着被推进到时刻 2，并执行一个到达事件。表 1 剩余部分的解读留给读者完成。

表 1 中的仿真覆盖了时间区间 [0,23]。在仿真时刻 23，系统清空，但下一次到达也发生在时刻 23。在所仿真的 23 个时间单位中，服务台忙碌了 21 个，最大队列长度为 2。当然，这次仿真太短，无法得出任何可靠的结论。注意，仿真表给出的是所有时刻的系统状态，而不仅是所列出的那些时刻。例如，从时刻 11 到时刻 15，有一位顾客正在接受服务，另有一位在排队线中。

当事件调度算法在仿真软件中实现时，软件只维护一套系统状态和其他属性值，也就是说，只维护一个快照（当前的或部分更新的那个）。若打算在 Java 或其他通用语言中实现事件调度，则应遵循如下规则：新的快照只能由此前的快照、新生成的随机变量以及事件逻辑（图 5 和图 6）导出。在推进时钟时，应当忽略以往的快照。当前快照必须包含继续仿真所需的全部信息。

表 1 收银台的仿真表