# 相对性能估计

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

---

本章讨论可用于比较两个或多个系统设计的相对性能（relative performance）的众多统计方法中的若干种。这是仿真最重要的用途之一。由于响应变量的观测值含有随机变异，因此需要借助统计分析来判别：所观察到的差异究竟源于设计上的差异，还是仅仅源于模型中固有的随机波动。

两个系统设计的比较比多个系统设计的同时比较更为容易。第 1 节讨论两个系统设计的情形，使用两种可能的统计技术：独立抽样（independent sampling）与相关抽样（correlated sampling）。相关抽样也称为公共随机数（common random numbers, CRN）技术；简言之，就是用同一组随机数去仿真两个备选的系统设计。如果实施得当，CRN 通常会降低性能度量之差的估计量的方差，从而在给定样本量下，能够比独立抽样给出更精确的均值之差估计。第 2 节把第 1 节的统计技术推广到多个（多于两个）系统设计的比较，使用同时置信区间（simultaneous confidence intervals）、筛选（screening）与选择最佳方案（selection of the best）。这些方法在所能考虑的系统设计数量方面有一定局限，因此第 3 节说明如何用更简单的元模型（metamodel）来表示数量庞大的复杂系统设计。最后，对于以较无结构的方式彼此关联的极大量系统设计的比较与评价，第 4 节给出基于仿真的优化（optimization via simulation）。

## 1 两个系统设计的比较

假设仿真分析人员希望比较一个系统的两种可能配置。在排队系统中，或许要比较两种可能的排队规则，或两组可能的服务台。在供应链库存系统中，或许要比较两种可能的订货策略。一个作业车间可以有两种可能的调度规则；一个生产系统可以有容量各异的在制品缓冲区。还可以举出许多其他备选系统设计的例子。

将采用重复法（method of replications）来分析输出数据。系统 \(i\) 的均值性能度量记为 \(\theta_{i} , i = 1 , 2 .\) 如果是稳态仿真，则假定已采用数据删除或其他恰当的技术，以确保点估计量是均值性能度量 \(\theta_{i}\) 的近似无偏估计量。仿真实验的目标是得到均值性能之差，即 \(\theta_{1} - \theta_{2}\) 的点估计与区间估计。下面将讨论计算 \(\theta_{1} - \theta_{2}\) 的置信区间的两种方法，但首先给出一个例子和总体框架。

### 例 1

Software Made Personal（SMP）公司设有一个客户支持呼叫中心，7 名话务员在东部时间上午 8:00 至下午 4:00 处理其软件用户的问题。客户来电时，会通过一套自动系统在两条产品线——财务软件与客户关系管理软件——之间进行选择。目前每条产品线都有自己的话务员（财务软件 4 名，客户关系管理 3 名）和各自的等候队列。SMP 想知道：能否通过交叉培训让话务员可以接听任一产品线的来电，从而减少所需的话务员总数。预计这样做会使处理一次通话的时间增加约 10%。现有系统如图 1(a) 所示，备选系统设计如图 1(b) 所示。在考虑减少话务员之前，SMP 想知道：如果只是对现有 7 名话务员进行交叉培训，其服务质量会怎样。来电被建模为非平稳泊松到达过程，话务员回答问题时间的分布也已拟合完成。比较将基于从发起呼叫到客户问题得到解答的平均时间。

![](img/ch12/img_3e90f609b6a8597e.jpg)

**图 1 (a) 现有 SMP 呼叫中心与 (b) 一种备选设计。**

**表 1 用于比较两系统的仿真输出数据与汇总度量**

| 系统 | 重复 | 样本均值 | 样本方差 |  |  |  |
|---|---|---|---|---|---|---|
| 1 | 2 | \(\cdots\) | \(R_{i}\) |  |  |  |
| 1 | \(Y_{11}\) | \(Y_{21}\) | \(\cdots\) | \(Y_{R_{1}1}\) | \(\bar{Y}_{.1}\) | \(S_{1}^{2}\) |
| 2 | \(Y_{12}\) | \(Y_{22}\) | \(\cdots\) | \(Y_{R_{2}2}\) | \(\bar{Y}_{.2}\) | \(S_{2}^{2}\) |

在比较两个系统（如例 1 中的两个系统）时，仿真分析人员必须为每个模型 \(i = 1 , 2 .\) 确定运行长度 \(T_{E} ^ {( i )}\)、每个模型要做的重复次数 \(R_{i}\)，或二者都确定。从系统 \(i\) 的第 \(r\) 次重复中，仿真分析人员得到均值性能度量 \(\theta_{i}\) 的一个估计 \(Y_{r i}\)。在例 1 中，\(Y_{r i}\) 是在系统 \(i\) 的第 \(r\) 次重复中所观察到的平均呼叫响应时间，其中 \(r = 1 , 2 , \ldots , R_{i} ;\) 且 \(i = 1 , 2\)。这些数据连同两个汇总度量——样本均值 \({\bar {Y}} _ {\cdot i}\) 与样本方差 \(S_{i} ^ {2}\)——在表 1 中给出。假定估计量 \(Y_{r i}\) 至少近似无偏，则有 

$$
\theta_{1} = E (Y_{r 1}), \quad r = 1, \ldots , R_{1}; \quad \theta_{2} = E (Y_{r 2}), r = 1, \ldots , R_{2}
$$

 在例 1 中，SMP 最初关心的是两个系统设计的比较（话务员人数相同，但专业化分工与交叉培训之别），因此仿真分析人员决定计算 \(\theta_{1} - \theta_{2}\)（即两个均值性能度量之差）的置信区间。这将导致三种可能的结论之一：

1. 若 \(\theta_{1} - \theta_{2}\) 的置信区间（c.i.）完全位于零的左侧，如图 2(a) 所示，则有强有力的证据支持 \(\theta_{1} - \theta_{2} < 0\) 的假设，或者说 \(\theta_{1} < \theta_{2}\)。在例 1 中，\(\theta_{1} < \theta_{2}\) 意味着系统 1（现有系统）的平均响应时间小于系统 2（备选系统）。

2. 若 \(\theta_{1} - \theta_{2}\) 的置信区间完全位于零的右侧，如图 2(b) 所示，则有强有力的证据支持 \(\theta_{1} - \theta_{2} > 0 .\)，或者说 \(\theta_{1} > \theta_{2}\)。在例 1 中，\(\theta_{1} > \theta_{2}\) 可解释为系统 2 优于系统 1，即系统 2 的平均响应时间更小。

3. 若 \(\theta_{1} - \theta_{2}\) 的置信区间包含零，则依据现有数据，没有强有力的统计证据表明某一个系统设计优于另一个。

![](img/ch12/img_88bf2aa5bad19a8d.jpg)

**图 2 两系统比较中可能出现的三种置信区间。**

有些统计学教科书说可以得出 \(\theta_{1} = \theta_{2}\) 这一弱结论，但这类说法可能具有误导性。"弱"结论往往根本不是结论。最可能的情况是：如果收集到足够多的补充数据（即增大 \(R_{i}\)），置信区间就会移动，并且其长度肯定会缩短，直至得出第 1 条或第 2 条结论。由于置信区间给出了 \(\theta_{1} - \theta_{2}\) 之估计量精度的度量，因此我们希望把它压缩到这样一种程度：要么能断定存在一个值得关注的差异，要么差异小到我们无须把二者区分开来。

在本章中，\(\theta_{1} - \theta_{2}\) 的双侧 \(100 ( 1 - \alpha ) \%\) 置信区间始终具有以下形式 

$$
\bar {Y} _ {. 1} - \bar {Y} _ {. 2} \pm t_{\alpha / 2, \nu} \mathrm{s.e.} (\bar {Y} _ {. 1} - \bar {Y} _ {. 2})\tag{1}
$$

 其中 \({\bar {Y}} _ {\cdot i}\) 是系统 \(i\) 在所有重复上的样本均值性能度量 

$$
\bar {Y} _ {\cdot i} = \frac {1}{R_{i}} \sum_{r = 1} ^ {R_{i}} Y_{r i}\tag{2}
$$

 ν 是与方差估计量相联系的自由度，\(t_{\alpha / 2 , 1}\) 是自由度为 ν 的 t 分布的 \(100 ( 1 {-} \alpha / 2 )\) 百分位点，而 \(\mathbf{s . e .} ( \cdot )\) 表示所给点估计量的标准误差。为得到标准误差与自由度，分析人员使用两种统计技术之一。两种技术都假定表 1 中的基础数据 \(Y_{r i}\) 近似服从正态分布。只要每个 \(Y_{r i}\) 本身是第 \(r\) 次重复中观测值的样本均值，这一假定就是合理的，而例 1 的情形正是如此。

由仿真实验的设计所决定，\(Y_{r 1} ( r = 1 , \dots , R_{1} )\) 独立同分布（i.i.d.），其均值为 \(\theta_{1}\)、方差设为 \(\sigma_{1} ^ {2}\)。同理，\(Y_{r 2} ( r = 1 , \ldots , R_{2} )\) 独立同分布，其均值为 \(\theta_{2}\)、方差设为 \(\sigma_{2} ^ {2}\)。计算式 (1) 中置信区间的两种技术基于不同的假设集，将在以下各小节中讨论。

系统性能上的统计显著性差异与实际显著性差异之间存在着重要区别。统计显著性回答这样一个问题：所观察到的差异 \({\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2}\) 是否大于 \({\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2}\) 中的变异性？这个问题也可以表述为：我们是否已经收集到足够的数据，可以有把握地断定所观察到的差异是真实的，而不仅仅是偶然？上述第 1 条与第 2 条结论意味着存在统计上显著的差异，而第 3 条结论则意味着所观察到的差异在统计上不显著，即便这些系统确实可能存在差异。统计显著性是仿真实验与输出数据的函数。

实际显著性回答这样一个问题：真实的差异 \(\theta_{1} - \theta_{2}\) 是否大到足以影响我们需要作出的决策？在例 1 中，我们可能得出 \(\theta_{1} > \theta_{2}\) 的结论，并判定系统 2 更好（期望响应时间更小）。然而，如果实际的差异 \(\theta_{1} - \theta_{2}\) 非常小——比如说小到客户根本察觉不到改善的程度——那么用系统 2 替换系统 1 可能就不值得付出相应的代价。实际显著性是系统之间实际差异的函数，且与仿真实验无关。

置信区间并不直接回答实际显著性的问题。相反，它们以 \(1 - \alpha\) 的概率把真实差异 \(\theta_{1} - \theta_{2}\) 界定在如下范围 

$$
\bar {Y} _ {. 1} - \bar {Y} _ {. 2} - t_{\alpha / 2, \nu} \text{s.e.} (\bar {Y} _ {. 1} - \bar {Y} _ {. 2}) \leq \theta_{1} - \theta_{2} \leq \bar {Y} _ {. 1} - \bar {Y} _ {. 2} + t_{\alpha / 2, \nu} \text{s.e.} (\bar {Y} _ {. 1} - \bar {Y} _ {. 2})
$$

 处于这些界限之内的差异是否具有实际显著性，取决于具体的问题。

### 1.1 独立抽样

独立抽样（independent sampling）是指用不同且相互独立的随机数流去仿真两个系统。这意味着仿真系统 1 的全部观测值，即 \(\{Y_{r 1} , r = 1 , \ldots , R_{1} \}\)，与仿真系统 2 的全部观测值，即 \(\{Y_{r 2} , r = 1 , \dots , R_{2} \}\)，在统计上相互独立。由式 (2) 以及各重复之间的独立性，样本均值 \({\bar {Y}} _ {\cdot i}\) 的方差为 

$$
\mathrm{V} (\bar {Y} _ {\cdot i}) = \frac {\mathrm{V} (Y_{r i})}{R_{i}} = \frac {\sigma_{i} ^ {2}}{R_{i}}, \qquad i = 1, 2
$$

 对于独立抽样，\({\bar {Y}} _ {\cdot 1}\) 与 \({\bar {Y}} _ {\cdot 2}\) 在统计上独立；因此，

$$
\begin{array}{r l} \mathrm{V} (\bar {Y}. _ {1} - \bar {Y}. _ {2}) & = \mathrm{V} (\bar {Y}. _ {1}) + \mathrm{V} (\bar {Y}. _ {2}) \\ & = \frac {\sigma_{1} ^ {2}}{R_{1}} + \frac {\sigma_{2} ^ {2}}{R_{2}} \end{array}\tag{3}
$$

 当两个方差相等但数值未知时，即 \(\sigma_{1} ^ {2} = \sigma_{2} ^ {2}\)，存在一些适用的置信区间程序。然而，由于在实践中我们很少知道这一点成立，而且只要重复次数不太少，利用方差相等所带来的好处也很有限，因此这一情形请读者参阅任一统计学教科书（例如 Hines 等人 [2002]）。

\(\theta_{1} - \theta_{2}\) 的近似 \(100 ( 1 - \alpha ) \%\) 置信区间可如下计算。点估计由下式计算 

$$
\bar {Y} _ {. 1} - \bar {Y} _ {. 2}\tag{4}
$$

 其中 \({\bar {Y}} _ {\cdot i}\) 由式 (2) 给出，而系统 \(i\) 的样本方差为 

$$
S_{i} ^ {2} = \frac {1}{R_{i} - 1} \sum_{r = 1} ^ {R_{i}} (Y_{r i} - \bar {Y} _ {. i}) ^ {2}
$$


$$
= \frac {1}{R_{i} - 1} \left(\sum_{r = 1} ^ {R_{i}} Y_{r i} ^ {2} - R_{i} \bar {Y} _ {\cdot i} ^ {2}\right)\tag{5}
$$

 因此，点估计的标准误差由下式给出 

$$
\text{s.e.} (\bar {Y} _ {. 1} - \bar {Y} _ {. 2}) = \sqrt {\frac {S_{1} ^ {2}}{R_{1}} + \frac {S_{2} ^ {2}}{R_{2}}}\tag{6}
$$

 自由度 ν 由下式近似 

$$
\nu = \frac {(S_{1} ^ {2} / R_{1} + S_{2} ^ {2} / R_{2}) ^ {2}}{[ (S_{1} ^ {2} / R_{1}) ^ {2} / (R_{1} - 1) ] + [ (S_{2} ^ {2} / R_{2}) ^ {2} / (R_{2} - 1) ]}\tag{7}
$$

 并四舍五入为整数。于是，置信区间由表达式 \(( 1 )\) 给出，并使用式 (6) 的标准误差。建议该程序至少使用 \(R_{1} \geq 6\) 和 \(R_{2} \geq 6\) 的重复次数。

### 1.2 公共随机数（CRN）

CRN（也称为相关抽样，correlated sampling）是指：在每一次重复中，用同一组随机数去仿真两个系统。因此，\(R_{1}\) 与 \(R_{2}\) 必须相等，设为 \(R_{1} = R_{2} = R\)。于是，对每一次重复 \(r ,\) 两个估计量 \(Y_{r 1}\) 与 \(Y_{r 2}\) 不再独立，而是彼此相关。不过，不同的重复使用独立的随机数，因此当 \(r \neq s\) 时，各对 \(( Y_{r 1} , Y_{s 2} )\) 相互独立。例如，在表 1 中，观测值 \(Y_{11}\) 与 \(Y_{12}\) 相关，但 \(Y_{11}\) 与所有其他观测值独立。使用 CRN 的目的在于对每一个 \(r\) 在 \(Y_{r 1}\) 与 \(Y_{r 2}\) 之间诱发正相关，从而在均值之差的点估计量 \({\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2}\) 上实现方差缩减。一般而言，该方差由下式给出 

$$
\begin{array}{c} \mathrm{V} (\bar {Y}. _ {1} - \bar {Y}. _ {2}) = \mathrm{V} (\bar {Y}. _ {1}) + \mathrm{V} (\bar {Y}. _ {2}) - 2 \mathrm{cov} (\bar {Y}. _ {1}, \bar {Y}. _ {2}) \\ = \frac {\sigma_{1} ^ {2}}{R} + \frac {\sigma_{2} ^ {2}}{R} - \frac {2 \rho_{12} \sigma_{1} \sigma_{2}}{R} \end{array}\tag{8}
$$

 其中 \(\rho_{12}\) 是 \(Y_{r 1}\) 与 \(Y_{r 2}\) 之间的相关系数。（按定义，\(\rho_{12} = \operatorname{c o v} ( Y_{r 1} , Y_{r 2} ) / \sigma_{1} \sigma_{2}\)，它不依赖于 \(r\)。）

现在，把使用 CRN 所产生的 \({\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2}\) 的方差［式 \(( 8 )\)，记为 \(V_{C R N}\)］与使用等样本量独立抽样所产生的方差［式 (3) 中 \(R_{1} = R_{2} = R\)，记为 \(V_{I N D}\)］作比较。注意 

$$
V_{C R N} = V_{I N D} - \frac {2 \rho_{12} \sigma_{1} \sigma_{2}}{R}\tag{9}
$$

 如果 CRN 如预期般发挥作用，相关系数 \(\rho_{12}\) 将为正；因此，式 (9) 右侧的第二项为正，从而 

$$
V_{C R N} <   V_{I N D}
$$

 也就是说，在点估计量上使用 CRN 的方差将小于使用独立抽样的方差。在相同样本量 \(R\) 下，方差更小意味着基于 CRN 的估计量更精确，从而得到更短的差异置信区间，这又意味着能够探测到更小的性能差异。

要用相关数据计算 \(100 ( 1 - \alpha ) \%\) 置信区间，首先计算差值 

$$
D_{r} = Y_{r 1} - Y_{r 2}\tag{10}
$$

 由 CRN 的定义，这些差值是独立同分布的；然后计算样本均值之差为 

$$
\bar {D} = \frac {1}{R} \sum_{r = 1} ^ {R} D_{r}\tag{11}
$$

 于是，\(\bar {D} = \bar {Y} _ {\cdot 1} - \bar {Y} _ {\cdot 2}\)。差值 \(\{D_{r} \}\) 的样本方差计算如下 

$$
S_{D} ^ {2} = \frac {1}{R - 1} \sum_{r = 1} ^ {R} (D_{r} - \bar {D}) ^ {2}
$$


$$
= \frac {1}{R - 1} \left(\sum_{r = 1} ^ {R} D_{r} ^ {2} - R \bar {D} ^ {2}\right)\tag{12}
$$

 其自由度为 \(\nu = R - 1\)。\(\theta_{1} - \theta_{2}\) 的 \(100 ( 1 - \alpha ) \% \mathrm{c . i}\). 由表达式 (1) 给出，其中 \({\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2} = {\bar {D}}\) 的标准误差估计为 

$$
\mathrm{s.e.} (\bar {D}) = \mathrm{s.e.} (\bar {Y}. _ {1} - \bar {Y}. _ {2}) = \frac {S_{D}}{\sqrt {R}}\tag{13}
$$

 由于式 (13) 中的 \(S_{D} / \sqrt {R}\) 是 \(\sqrt {V_{C R N}}\) 的估计，而表达式 (6) 是 \(\sqrt {V_{I N D}}\) 的估计，因此若 \(\rho_{12} > 0\)，CRN 在给定样本量下所产生的置信区间通常会比独立抽样所产生的置信区间更短。事实上，只要 \(R > 10\)，当 \(\rho_{12} > 0 . 1\) 时，使用 CRN 的置信区间期望长度就会更短。\(R\) 越大，能产生更短期望长度所要求的 \(\rho_{12}\) 就可以越小 [Nelson 1987]。

对任何一个问题，实施公共随机数的方法都有很多种。仅仅在随机数发生器上使用相同的种子是远远不够的。在一个模型中用于某种目的的每一个随机数，都应当在第二个模型中用于同样的目的——也就是说，随机数的使用必须是同步的（synchronized）。例如，若第 \(i\) 个随机数在模型 1 中用于生成第 5 位来电者的话务员服务时间，那么第 \(i\) 个随机数在模型 2 中就应当用于完全相同的目的。对于排队系统或服务设施，公共随机数的同步保证两个系统面对完全相同的工作负荷：两个系统在相同的时刻面临到达，且这些到达要求等量的服务。（在给定到达在两个模型中的实际服务时间未必相等；如果一个模型中的服务台比另一个模型中的服务台更快，它们可以成比例。）对于库存系统，在比较不同的订货策略时，同步保证两个系统对某种给定产品面对完全相同的需求。对于生产系统或可靠性系统，同步保证给定机器的停机在两个模型中恰好在相同时刻发生，且持续时间相同。另一方面，如果其中一个系统的某个方面与另一个系统完全不同，那么同步可能就是不恰当的——甚至根本无法实现。总之，两个系统设计中那些足够相似的方面，应当用公共随机数来仿真，使两个模型表现相似；而那些完全不同的方面，则应当用独立的随机数来仿真。

公共随机数的实施依赖于具体模型，但可以给出若干准则，使 CRN 更有可能产生正相关。这些准则的目的在于确保同步得以实现：

1. 把一条随机数流专门用于某个特定目的，并按需要使用尽可能多的不同流。（可以用不同的随机数发生器，或同一发生器上相隔很远的种子，来获得两条不同的、不重叠的流。）此外，在每一次重复开始时为每条流指定独立选取的种子。仅在第一次重复开始时指定种子，然后让随机数发生器在第二次及以后的重复中继续往下走，这是不够的。如果以这种方式进行仿真，第一次重复将是同步的，但后续的重复可能就不是了。

2. 对于带有外部到达的系统（或子系统）：每当一个实体进入系统，就生成下一个到达间隔时间，然后立即按固定顺序生成该到达实体所需的、且在两个模型中相同的所有随机变量（如服务时间、订货量等），并把它们作为该实体的属性存储起来，以便之后按需使用。应用准则 1：把一条随机数流专门用于这些外部到达及其全部属性。

3. 对于有实体以循环或重复方式执行给定活动的系统，为这个实体指定一条随机数流。（例如：一台机器在两个状态之间循环：运行—停机—运行—停机—……。用一条专用随机数流来生成运行时间与停机时间。）

4. 如果同步不可能实现，或者对两个模型的某一部分而言同步并不恰当，那么就对这一子集的变量使用独立的随机数流。

遗憾的是，并不能保证 CRN 总能在两个模型的可比运行之间诱发正相关。已知的是：若对于每个输入随机变量 \(X\)，估计量 \(Y_{r 1}\) 与 \(Y_{r 2}\) 都是 \(X\) 的增函数（或都是 \(X\) 的减函数），则 \(\rho_{12}\) 为正。直观的想法是：两个模型，即 \(Y_{r 1}\) 与 \(Y_{r 2} ,\) 对每个输入随机变量都沿相同方向作出响应，这就导致了正相关。响应变量相对于输入随机变量的这种递增或递减性质称为单调性（monotonicity）；已知它对某些排队系统成立，例如当响应变量为顾客延误时的 \(G I / G / c\) 排队，因此有证据表明公共随机数对排队仿真是一种值得采用的技术。（对于简单排队，顾客延误是服务时间的增函数、到达间隔时间的减函数。）然而，Wright 与 Ramsay [1979] 报告了某些库存仿真出现了负相关。总之，应当遵循上述准则，并且应当能够看出所关心的响应变量是随机输入变量的单调函数，这一点要有合理的把握。

### 例 1（续）

图 1 中 SMP 呼叫中心的两种设计将用 CRN 进行仿真比较。为说明若干问题，我们考察两个比较：现有 7 名话务员的系统与交叉培训的 7 名话务员系统；以及交叉培训的 7 名话务员系统与交叉培训的 6 名话务员系统。我们还将评估使用 CRN 相对于独立抽样所获得的收益有多大。

SMP 问题中有四个驱动性的输入过程：来电者的到达间隔时间、来电者类型、财务软件来电者的通话时间，以及客户关系管理软件来电者的通话时间。这里所谓"通话时间"，是指话务员与来电者相处的时间，不包括来电者在线等待的时间。

为了实现 CRN，我们为每个输入过程专门指定一条不同的随机数流。这保证在每一次重复中，每个备选方案都会看到完全相同的来电序列与来电者类型。不过，通话时间可以用两种方式之一生成：要么按需生成（即当来电者接通话务员时，从相应的分布生成通话时间），要么在来电者到达的瞬间生成通话时间，把它作为来电者的一个属性保存起来，并在来电者接通话务员时使用它。后一种做法（在到达时生成通话时间并保存）在本例中提供了更好的同步；事实上，如果交叉培训方案下的通话时间是用现有系统的通话时间乘以 1.1 得到的，那么这种做法就保证交叉培训方案中每一位来电者的通话时间恰好比其在现有设计下所经历的通话时间多 10%。

表 2 给出了三种系统设计在 300 次重复运行中前 10 次重复的结果（完整数据可从 www.bcnn.net 获取）。图 3 绘出了同样这 10 次重复中现有系统与交叉培训 7 人系统的平均响应时间。CRN 的预期效果在图中显而易见：两种系统设计的响应均值倾向于同升同落，使其差异更易于观察。基于全部 300 次重复，它们之间的样本相关系数为 0.24，表明 CRN 的影响并不大。然而，同样基于全部 300 次重复，7 人与 6 人两个备选系统之间的样本相关系数高达 0.94。如图 4 所示，这两种系统设计的响应明显地同升同落，呈现出强正相关。

![](img/ch12/img_0cb0fd894074c58c.jpg)

**图 3 SMP 呼叫中心现有方案与备选方案（7 名话务员）前 10 次重复的平均响应时间。**

观察 CRN 影响的更定量化的方式是考察标准误差，因为它们直接影响置信区间的长度。例如，两个备选系统比较的标准误差为 0.88 分钟。我们可以利用样本方差和式 (6) 来近似独立抽样情形下本应有的标准误差：

$$
\mathrm{s.e.} = \sqrt {\frac {4 . 86}{10} + \frac {23 . 99}{10}} = 1. 70
$$

 降低幅度接近 50%，使得差异更容易被探测出来。

**表 2 SMP 呼叫中心各系统设计的比较**

|  | 平均响应时间 | 差值 |  |  |  |
|---|---|---|---|---|---|
| 重复 | 现有方案 | 备选 7 人 | 备选 6 人 | 现有 − 备选7 | 备选7 − 备选6 |
| 1 | 6.24 | 6.19 | 8.35 | 0.05 | -2.16 |
| 2 | 9.06 | 10.64 | 18.03 | -1.59 | -7.39 |
| 3 | 8.02 | 9.53 | 16.17 | -1.51 | -6.64 |
| 4 | 5.93 | 6.15 | 7.40 | -0.22 | -1.25 |
| 5 | 8.31 | 7.83 | 12.70 | 0.48 | -4.87 |
| 6 | 5.91 | 6.09 | 8.26 | -0.17 | -2.17 |
| 7 | 8.74 | 7.62 | 12.32 | 1.12 | -4.70 |
| 8 | 7.78 | 7.03 | 11.40 | 0.75 | -4.37 |
| 9 | 7.15 | 12.79 | 23.04 | -5.64 | -10.24 |
| 10 | 5.72 | 7.57 | 14.30 | -1.85 | -6.73 |
| 样本均值 | 7.29 | 8.14 | 13.20 | -0.86 | -5.05 |
| 样本方差 | 1.60 | 4.86 | 23.99 | 3.87 | 7.72 |
| 标准误差 |  |  |  | 0.62 | 0.88 |

![](img/ch12/img_cdc547ec58bd5a56.jpg)

**图 4 SMP 呼叫中心备选设计中 7 人与 6 名话务员前 10 次重复的平均响应时间。**

理解这两种情形下由 CRN 诱发的相关强度为何不同，有助于理解 CRN 何时最为有效。两个备选设计（6 人与 7 人）在结构上是同一系统，仅在话务员人数上不同。因此，当随机输入使得 7 人系统发生拥堵时，6 人系统必定至少同样拥堵，而且很可能更加拥堵。现有系统与各备选方案在结构上不同，它为每种来电者设置了独立的队列。当财务软件队列发生拥堵时，客户关系管理队列未必发生拥堵，于是财务来电者的长响应时间可能会与客户关系管理来电者的短响应时间相平均。因此，它与采用单一队列的备选方案之间的关系就不那么强。

令下标 C 表示现有呼叫中心，P7 表示 7 名话务员的备选方案，P6 表示 6 名话务员的备选方案。那么 \(\theta_{C} - \theta_{P 7}\) 的 95% 置信区间为 

$$
- 0. 86 \pm 2. 26 (0. 62)
$$

 即 

$$
- 2. 26 \leq \theta_{C} - \theta_{P 7} \leq 0. 54\tag{14}
$$

 因此，根据这 10 次重复，我们不能断定哪一种设计的响应时间更低，尽管我们可以很有把握地说，差异在 1.4 分钟以内。

另一方面，正如我们所预料的，7 人与 6 人系统明显不同，其 95% 置信区间为 

$$
- 7. 02 \leq \theta_{P 7} - \theta_{P 6} \leq - 3. 06
$$

 这表明，拥有第 7 名话务员平均至少比只有 6 名要快 3 分钟。

### 1.3 具有规定精度的置信区间

假设我们希望 \(\theta_{1} - \theta_{2}\) 的估计误差小于 \(\epsilon\)（\(\epsilon\) 可以是具有实际显著性的差异）。因此，我们的目标是找到重复次数 \(R\)，使得 

$$
H = t_{\alpha / 2, \nu} \text{s.e.} (\bar {Y}. _ {1} - \bar {Y}. _ {2}) \leq \epsilon\tag{15}
$$

 我们首先对每个系统做 \(R_{0} \geq 10\) 次重复，以得到 \(\mathrm{s . e} . ( {\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2} )\) 的初始估计。然后近似出达到式 (15) 半长准则所需的总重复次数 \(R \ \geq \ R_{0}\)。最后，再对每个系统做 \(R {-} R_{0}\) 次追加重复（或重新做 \(R\) 次重复），计算置信区间，并检验半长准则是否已经达到。

### 例 1（续）

回顾一下，\(R_{0} = 10\) 次重复且采用 CRN 得到了式 (14) 中现有系统与 7 名交叉培训话务员系统期望响应时间之差的 95% 置信区间；该区间可改写为 \(- 0 . 86 \pm 1 . 40\) 分钟。虽然现有呼叫中心设计似乎具有更小的期望响应时间，但由于置信区间包含 0，该差异在统计上并不显著。假设大于 0.25 分钟的差异被认为具有实际显著性，因此我们希望做足够多的重复以获得 \(H \leq \epsilon = 0 . 25\)。例 1 中所用的置信区间形式为 \(\bar {D} \pm t _ {\alpha / 2 , R_{0} - 1} S_{D} / \sqrt {R_{0}}\)，具体取值为 \(\bar {D} = - 0 . 86 , R_{0} = 10 , t_{0 . 025 , 9} = 2 . 26\)，以及 \(S_{D} ^ {2} = 3 . 87\)。为获得所需的精度，我们需要找到满足下式的 \(R\) 

$$
\frac {t_{\alpha / 2 , R - 1} S_{D}}{\sqrt {R}} \leq \epsilon
$$

 我们可以通过寻找满足 \(R \geq R_{0}\) 且满足下式的最小整数 \(R\) 来近似求解 

$$
R \geq \left(\frac {z_{\alpha / 2} S_{D}}{\epsilon}\right) ^ {2}
$$

 代入 \(z_{0 . 025} = 1 . 96\) 与 \(S_{D} ^ {2} = 3 . 87\)，得到 

$$
R \geq \frac {(1 . 96) ^ {2} (3 . 87)}{(0 . 25) ^ {2}} = 238
$$

 这意味着需要 238 次重复，比初始实验多 228 次。www.bcnn.net 上提供了 300 次重复的数据，鼓励读者去检验这一结果。

## 2 多个系统设计的比较

假设仿真分析人员希望比较 K 个备选的系统设计。比较将基于系统 \(i\) 的某个指定性能度量 \(\theta_{i}\) 来进行，\(i = 1 , 2 , \dots , K\)。已经发展出许多不同的统计程序，可用于分析仿真数据并作出关于参数 \(\theta_{i}\) 的统计上可靠的推断。这些程序可分为固定样本量（fixed-sample-size）程序与序贯抽样（sequential-sampling，或称多阶段）程序。

在第一类程序中，使用预先确定的样本量（即运行长度与重复次数），通过假设检验或置信区间作出推断。固定样本量程序的例子包括：均值性能度量的区间估计，以及两个系统均值性能度量之差的区间估计［如第 1 节中的表达式 (1)］。固定样本量程序的优点包括：在开展实验之前，其以计算机时间计的成本是已知的或易于估计的。当计算机时间有限，或正在进行一项预研（pilot study）时，固定样本量程序可能是恰当的。在某些情况下，明显处于劣势的系统设计可以在这一早期阶段就被排除。一个主要缺点是：可能得不出强有力的结论。例如，置信区间可能过宽而无实用价值，因为区间宽度是点估计量精度的指示。假设检验可能导致无法拒绝原假设，这一般是一个弱结论，意味着无论原假设为真还是为假，都没有强有力的证据。

序贯抽样方案是指：不断收集越来越多的数据，直至得到一个具有预先规定精度的估计量，或者直至从若干备选方案中选出一个，且正确选择的概率大于预先规定的值。两阶段（或多阶段）程序是指：用初始样本来估计还需要多少补充观测，才能以规定的精度得出结论。第 1.3 节给出了一个用于估计两个系统性能度量的两阶段程序的例子。

应当使用哪一种合适的程序，取决于仿真分析人员的目标。可能的目标如下：

1. 估计每一个参数 \(\theta_{i}\)；

2. 把每一个性能度量 \(\theta_{i}\) 与一个对照 \(\theta_{1}\) 作比较（其中 \(\theta_{1}\) 可以表示现有系统的均值性能）；

3. 所有两两比较 \(\theta_{i} - \theta_{j}\)，\(i \neq j\)；

4. 选择最佳的 \(\theta_{i}\)（最大或最小）。

前三个目标将通过构造置信区间来实现。这类置信区间的个数分别为 \(C = K\)、\(C = K - 1\) 和 \(C = K ( K {-} 1 ) / 2\)。Hochberg 与 Tamhane [1987] 以及 Hsu [1996] 是关于这类多重比较（multiple comparison）程序的详尽参考文献。第四个目标需要使用一类称为多重排序与选择（multiple ranking and selection）程序的统计程序。Kleijnen [1975，第 II 章与第 V 章] 讨论了实现这些目标及其他目标的程序，并讨论了它们的相对优缺点。Goldsman 与 Nelson [1998] 以及 Law [2007] 讨论了与仿真最相关的那些选择程序。Bechhofer、Santner 与 Goldsman [1995] 是一份详尽的参考文献。下一小节给出一个可用于实现目标 1、2、3 且在很广的情形下都适用的固定样本量程序。第 2.2 小节给出一个实现目标 4 的程序。

### 2.1 多重比较的 Bonferroni 方法

假设计算了 C 个置信区间，第 \(i\) 个区间的置信系数为 \(1 - \alpha_{i}\)。令 \(\boldsymbol {S} _ {i}\) 表示"第 \(i\) 个置信区间包含所估计的参数（或两个参数之差）"这一陈述。对给定的一组数据，这一陈述可能为真也可能为假，但导出该区间的程序设计使得陈述 \(\boldsymbol {S_{i}}\) 将以概率 \(1 - \alpha_{i}\) 为真。当希望对若干参数同时作出陈述时（如目标 1、2、3），分析人员希望对所有陈述同时为真有较高的把握。Bonferroni 不等式指出 

$$
P (\text{all   statements} \mathcal {S} _ {i} \text{are   true}, i = 1, \dots , C) \geq 1 - \sum_{j = 1} ^ {C} \alpha_{j} = 1 - \alpha_{E}\tag{16}
$$

 其中 \(\begin{array} {r} {\alpha_{E} = \sum_{j = 1} ^ {C} \alpha_{j}} \end{array}\) 称为总体错误概率（overall error probability）。表达式 (16) 可以改写为 

$$
P (\text{one   or   more   statements} \mathcal {S} _ {i} \text{is   false}, i = 1, \dots , C) \leq \alpha_{E}
$$

 或等价地，

$$
\begin{array}{l} P (\text{one   or   more   of   the} C \text{confidence   intervals   does   not} \\ \text{contain   the   parameter   being   estimated}) \leq \alpha_{E} \end{array}
$$

 于是，\(\alpha_{E}\) 提供了作出错误结论之概率的上界。要开展一项涉及 C 次比较的实验，首先选定总体错误概率，比如 \(\alpha_{E} = 0 . 05\) 或 \(0 . 10\)。各个 \(\alpha_{j}\) 可以选为相等的 \(( \alpha_{i} = \alpha_{E} / C )\)，也可以按需取为不相等。\(\alpha_{j}\) 的值越小，第 \(j\) 个置信区间就越宽。例如，若构造两个 95% 置信区间 \(( \alpha_{1} = \alpha_{2} =\)0.05)，总体置信水平将为 90% 或更高 \(( \alpha_{E} = \alpha_{1} + \alpha_{2} = 0 . 10 )\)。若构造十个 95% 置信区间 \(( \alpha_{i} = 0 . 05 , i = 1 , \ldots , 10 )\)，所得到的总体置信水平可能低至 50%\(\textstyle ( \alpha_{E} = \sum_{i = 1} ^ {10} \alpha_{i} = 0 . 50 )\)，这远低于实用的要求。为在进行 10 次比较时保证 95% 的总体置信水平，一种做法是为所关心的参数（或差值）构造十个 99.5% 的置信区间。

多重置信区间的 Bonferroni 方法以表达式 (16) 为基础。它的一个主要优点是：无论各备选设计的模型是用独立抽样运行还是用公共随机数运行，它都成立。

Bonferroni 方法在进行大量比较时的主要缺点是每一个单独区间的宽度都会增加。例如，对于给定的一组数据和较大的样本量，99.5% 置信区间将是 95% 置信区间长度的 \(z_{0 . 0025} / z_{0 . 025} = 2 . 807 / 1 . 96 = 1 . 43\) 倍。对于小样本量，比如说样本量为 \(5 ,\) 99.5% 置信区间将是单个 95% 置信区间的 \(t_{0 . 0025 , 4} / t_{0 . 025 , 4} = 5 . 598 / 2 . 776 = 1 . 99\) 倍。置信区间的宽度是估计精度的度量。由于这些原因，建议只在比较次数较少时使用 Bonferroni 方法。20 次左右似乎是其实际的上限。

对应于目标 1、2、3，在比较 K 个备选系统设计时，使用 Bonferroni 不等式 (16) 至少有三种可能的方式：

1. 单独的置信区间：利用表达式 \(H = t_{\alpha / 2 , R - 1} \frac {S} {\sqrt {R}}\) 为参数 \(\theta_{i}\) 构造 \(100 ( 1 - \alpha_{i} ) \%\) 置信区间，此时区间个数为 \(C = K\)。如果使用的是独立抽样，这 K 个置信区间将相互独立，因此总体置信水平将是 \(( 1 - \alpha_{1} ) \times ( 1 - \alpha_{2} ) \times \cdot \cdot \cdot \times ( 1 - \alpha_{C} )\)，它比表达式 (16) 的右侧更大，但大不了多少。这类程序最常用于估计单一系统的多个参数，而不是用于比较系统；由于来自同一系统的多个参数估计很可能是相关的，因此通常需要 Bonferroni 不等式。

2. 与现有系统比较：把所有设计与某一个特定设计作比较，通常是与现有系统作比较。即，利用表达式 (1) 为 \(\theta_{i} - \theta_{1} ( i = 2 , 3 , \ldots , K )\) 构造 \(100 ( 1 - \alpha_{i} ) \% \ \mathrm{c . i}\)。（假定系统 1 的性能度量 \(\theta_{1}\) 为现有系统。）在这种情形下，区间个数为 \(C = K - 1\)。这类程序最常用于把若干个竞争方案与现行系统作比较，以了解哪些更好。

3. 所有两两比较：把所有设计彼此作比较。即，对任意两个系统设计 \(i \neq j ,\) 为 \(\theta_{i} - \theta_{j}\) 构造 \(100 ( 1 - \alpha_{i i} ) \%\) 置信区间。有 K 个设计时，所计算的置信区间个数为 \(C = \dot {K} ( K - 1 ) / 2\)。总体置信系数的下界为 \(\begin{array} {r} {1 - \alpha_{E} = 1 - \sum \sum_{i \neq j} \alpha_{i j}} \end{array}\)（这由表达式 (16) 推出）。通常认为，CRN 会使真实的总体置信水平大于表达式 (16) 的右侧，且通常大于独立抽样时的水平。可以把表达式 (16) 的右侧看作给出了最坏的情形，即可能达到的最低总体置信水平。

### 例 2

重新考虑例 1 的呼叫中心设计问题。备选的系统设计如下：

a. C：现有系统（各自独立的话务员队列）；

b. P7：7 名交叉培训话务员的系统（单一队列）；

c. P6：6 名交叉培训话务员的系统（单一队列）。

利用表 2 中的数据，将为 \(\theta_{C} - \theta_{P 7} , \theta_{C} - \theta_{P 6}\) 和 \(\theta_{P 7} - \theta_{P 6}\) 构造置信区间，使总体置信水平为 95%。回顾一下，所有模型都使用了 CRN，但这并不影响总体置信水平，因为如前所述，无论数据在统计上独立还是相关，Bonferroni 不等式 (16) 都成立。

由于总体错误概率为 \(\alpha_{E} = 0 . 05\)，且要构造 \(C = 3\) 个置信区间，取 \(\alpha_{i} = 0 . 05 / 3 = 0 . 0167\)。然后用经过适当修改的表达式 (1) 构造 \(C = 3\) 个置信区间，其中 \(\alpha = \alpha_{i} = 0 . 0167\)，自由度 \(\nu = 10 - 1 = 9\)。值 \(t _ {\alpha_{i} / 2 , R - 1} = t_{0 . 0083 , 9} = 2 . 97\) 由附录表 A.5 插值得到；点估计与标准误差取自表 2。

三个置信区间（总体置信系数至少为 95%）由下式给出 

$$
\theta_{C} - \theta_{P 7}: - 0. 86 \pm 2. 97 (0. 62)
$$


$$
\theta_{C} - \theta_{P 6}: - 5. 91 \pm 2. 97 (1. 42)
$$


$$
\theta_{P 7} - \theta_{P 6}: - 5. 05 \pm 2. 97 (0. 88)
$$

 或 

$$
\theta_{C} - \theta_{P 7}: - 0. 86 \pm 1. 84
$$


$$
\theta_{C} - \theta_{P 6}: - 5. 91 \pm 4. 22
$$


$$
\theta_{P 7} - \theta_{P 6}: - 5. 05 \pm 2. 61
$$

 仿真分析人员有较高的把握（至少 95%）相信三条置信陈述都正确。注意 \(\theta_{C} - \theta_{P 7}\) 的置信区间仍然包含零；因此，现有设计与拥有 7 名交叉培训话务员的备选方案之间不存在统计上显著的差异，这一结论支持了例 1 中先前的结果。其余置信区间完全位于 0 的左侧，表明 C 与 P7 都优于 P6。

本章末的部分习题提供了比较 CRN 与独立抽样、以及在 Bonferroni 方法下计算同时置信区间的机会。

### 2.2 选择最佳方案

假设有 K 个系统设计，第 \(i\) 个系统性能的未知期望值为 \(\theta_{i}\)。我们关心哪一个系统设计最佳，其中"最佳"是指 \(\theta_{i}\) 最大或最小，取决于具体的问题。例如，在例 3 中我们将比较《绝对性能估计》一章中所描述的半导体制造设施的 \(K = 8\) 个设计，其中 \(\theta_{i}\) 是第 \(i\) 个设计下某一个产品族的稳态平均周期时间，而 \(\theta_{i}\) 越小越好。我们希望有一种能够处理很大 K（比如说 100 量级）的程序。

令 B 表示最佳系统设计的（未知）下标。真实差异 $| \theta_{B} - \theta_{i} | , i \neq B $ 越小，并且我们越想确保找到最佳系统，那么实现目标所需的重复次数就越多。因此，我们不去要求无论如何都必须找到最佳设计 B，而是作出折中：当系统 B 与其他系统之间的差异足够大时，以高概率找到 B。更准确地说，我们希望：只要 $| \theta_{B} - \theta_{i} | \geq \epsilon $ 对所有 \(i \neq B\) 成立（其中 \(\epsilon\) 取决于具体问题），选到最佳系统的概率就至少为 \(1 - \alpha\)。如果存在与最佳方案相差在 \(\epsilon\) 以内的系统，那么选到最佳方案或其中任一接近最佳的设计，我们都会感到满意。正确选择的概率 \(1 - \alpha\) 与具有实际显著性的差异 \(\epsilon\) 都在我们的掌控之下。

下面的程序通过两个阶段的仿真来实现这一目标：在第一阶段，从每一个系统设计各获得 \(R_{0}\) 次重复。那些在统计上显著劣于其他系统的系统被筛选掉，即它们以 \(1 - \alpha / 2\) 的置信度被排除，不再作进一步考虑。如果筛选后仍有不止一个系统存活，那么存活者获得足够的追加重复，以便以 \(1 - \alpha / 2\) 的置信度选出最佳者，或与最佳者相差在 \(\epsilon\) 以内的某个系统设计。于是，利用与 Bonferroni 不等式有关的推理，所选出的系统在筛选与选择两者的组合下，以至少 \(1 - \alpha\) 的置信度是最佳系统或接近最佳的系统。

第一阶段的筛选可能非常重要，因为正如例 3 将要说明的，第二阶段为选出最佳方案所需的重复次数可能相当大，而没有理由把时间和精力浪费在明显没有竞争力的系统设计上。当 K 很大时，这一点尤为重要。

下面的程序在 SimulationTools.xls 中实现，该文件可从 www.bcnn.net 获取。SimulationTools.xls 的实现会引导用户了解每一阶段需要获得哪些数据。

### 选择最佳方案程序（Select-the-Best Procedure）

1. 指定所需的正确选择概率 \(1 / k < 1 - \alpha < 1\)、具有实际显著性的差异 \(\epsilon > 0\)、初始重复次数 \(R_{0} \geq 10\)，以及竞争系统的个数 K。令 

$$
t = t _ {1 - (1 - \alpha / 2) ^ {\frac {1}{k - 1}}, R_{0} - 1}
$$

 并取得 Rinott 常数 \(h = h ( R_{0} , K , 1 - \alpha / 2 )\)。［提示：这两个临界值都由 SimulationTools.xls 自动计算；附录中给出了 t 和 h 的小规模数值表。］

2. 对每个系统做 \(R_{0}\) 次重复。计算第一阶段的样本均值与样本方差 

$$
\bar {Y} _ {\cdot i} = \frac {1}{R_{0}} \sum_{r = 1} ^ {R_{0}} Y_{r i}
$$


$$
S_{i} ^ {2} = \frac {1}{n_{0} - 1} \sum_{r = 1} ^ {n_{0}} \left(Y_{r i} - \bar {Y} _ {. i}\right) ^ {2}
$$

 对 \(i = 1 , 2 , \dots , K\)。

3. 计算筛选阈值 

$$
W_{i j} = t \left(\frac {S_{i} ^ {2} + S_{j} ^ {2}}{R_{0}}\right) ^ {1 / 2}
$$

 对所有 \(i \neq j\)。

(a) 若越大越好，则构造存活子集 \(s\)，使其包含满足下式的每一个系统设计 \(i\) 

$$
\bar {Y} _ {\cdot i} \geq \bar {Y} _ {\cdot j} - \max \{0, W_{i j} - \epsilon \} \quad \text{for   all} j \neq i
$$

 也就是说，保留每一个样本均值比其他任何系统的样本均值小不超过 \(\max \{0 , W_{i j} - \epsilon \}\) 的系统 \(i\)。

(b) 若越小越好，则构造存活子集，使其包含满足下式的每一个系统设计 \(i\) 

$$
\bar {Y} _ {. i} \leq \bar {Y} _ {. j} + \max \{0, W_{i j} - \epsilon \} \quad \text{for   all} j \neq i
$$

 也就是说，保留每一个样本均值比其他任何系统的样本均值大不超过 \(\max \{0 , W_{i j} - \epsilon \}\) 的系统 \(i\)。

4. 若该集合只包含一个系统，则停止并返回该系统作为最佳系统。否则，对 \(\mathcal{S}\) 中的所有设计 \(i\)，计算第二阶段样本量 

$$
R_{i} = \max \left\{R_{0}, \lceil (h S_{i} / \epsilon) ^ {2} \rceil \right\}
$$

 其中 \(\lceil \cdot \rceil\) 表示向上取整。

5. 从 \(\mathcal{S}\) 中的所有系统 \(i\) 再取 \(R_{i} - R_{0}\) 次追加重复；或者，如果更方便的话，从头开始从系统 \(i\) 获得总共 \(R_{i}\) 次重复。

6. 计算总样本均值 

$$
\bar {\bar {Y}} _ {\cdot i} = \frac {1}{R_{i}} \sum_{r = 1} ^ {R_{i}} Y_{r i}
$$

 对 \(\mathcal{S}\) 中的所有 \(i\)。若越大越好，则选择 \(\bar {\bar {Y}} _ {\cdot i}\) 最大的系统为最佳。否则选择 \(\bar {\bar {Y}} _ {\cdot i}\) 最小的系统为最佳。

临界值 t 与 h 通常无法直接从表中查得，因此常常需要插值。作为一种替代，多数的数值分析软件和电子表格软件都提供了计算 t 分布分位数的算法；遗憾的是，Rinott 常数 h 却不是这样。因此，强烈建议使用 SimulationTools.xls 来执行这一程序。

Nelson 等人 [2001] 证明：只要数据服从正态分布且各系统独立仿真，选择最佳方案程序就能以至少 \(1 - \alpha\) 的置信水平找到具有最大或最小 \(\theta_{i}\) 的系统，或找到与最佳 \(\theta_{i}\) 相差在 \(\epsilon\) 以内的系统。如果应用了 CRN，该程序同样有效。

### 例 3

回顾《绝对性能估计》一章中 FastChip 公司的半导体制造例子。FastChip 公司希望用仿真来设计一个新的"晶圆厂"（fab），评估在两种不同产品（称为 C 芯片与 D 芯片）的某种特定装载下所能期望得到的稳态平均周期时间（从晶圆投料到完工的时间）。我们在《绝对性能估计》一章中所评价的那个系统设计实际上是一个基准情形，此外还有七个备选设计，它们来自对制造过程中某些步骤实施方式的不同选择；这些步骤称为清洗（Clean）、装石英（Load Quartz）、氧化（Oxidize）、卸石英（Unload Quartz）、涂胶（Coat）和显影（Develop）。问题是：哪一个设计带来最小的稳态平均周期时间？我们将关注 C 芯片的周期时间，并使用在那里得到的预热期和运行长度。

![](img/ch12/img_1dff73211e3b246d.jpg)

**图 5 SimulationTools.xls 中选择最佳方案程序的输入界面。**

包括基准情形在内共有 \(K \ = \ 8\) 个系统设计，性能度量 \(\theta_{1} , \theta_{2} , \ldots , \theta_{8}\) 是各设计下 C 芯片的稳态平均周期时间。FastChip 认定即使是很小的周期时间改善也很重要，因此他们设定 \(\epsilon = 0 . 15\) 小时（即 9 分钟）。当周期时间平均在 45 小时左右时，这是一个非常严格的要求，我们将在下面看到这一选择的后果。FastChip 希望有 95% 的把握选到最佳设计，或选到与最佳设计相差 9 分钟以内的设计。图 5 显示了在 SimulationTools.xls 中设置该问题的输入界面，其中用"无差异水平"（indifference level）一词来表示具有实际显著性的差异 \(\epsilon\)。

表 3 包含从每个系统设计的 \(R_{0} ~ = ~ 10\) 次重复中得到的第一阶段结果。调整氧化步骤的那个备选方案似乎给出了最小的平均周期时间，但所有备选方案看起来都很接近。本例的完整数据可从 www.bcnn.net 获取。

筛选所需的 t 值为 

$$
t = t _ {1 - (1 - \alpha / 2) ^ {\frac {1}{k - 1}}, R_{0} - 1} = t _ {1 - (0. 975) ^ {\frac {1}{7}}, 9} = t_{0. 0036, 9} = 3. 455
$$


**表 3 FastChip 仿真的第一阶段数据**

| 方案 \(i\) | 基准 1 | 清洗 2 | 装石英 3 | 氧化 4 | 卸石英 5 | 涂胶 6 | 步骤 7 | 显影 8 |
|---|---|---|---|---|---|---|---|---|
| \(\bar{Y}_{.i}\) | 46.86 | 45.70 | 47.23 | 45.13 | 46.80 | 47.81 | 47.41 | 46.94 |
| \(S_{i}^{2}\) | 0.21 | 0.37 | 0.29 | 0.09 | 0.28 | 0.98 | 0.12 | 0.28 |

**表 4 FastChip 仿真的筛选阈值**

| 方案 \(i\) | 基准 1 | 清洗 2 | 装石英 3 | 氧化 4 | 卸石英 5 | 涂胶 6 | 步骤 7 | 显影 8 |
|---|---|---|---|---|---|---|---|---|
| \(\bar{Y}_{.i}\) | 46.86 | 45.70 | 47.23 | 45.13 | 46.80 | 47.81 | 47.41 | 46.94 |
| \(W_{i4}\) | 0.60 | 0.74 | 0.68 |  | 0.67 | 1.13 | 0.50 | 0.67 |
| \(\leq Y_{.4} + \max\{W_{i4} - 0.15\}\) | 45.58 | 45.72 | 45.66 |  | 45.65 | 46.11 | 45.48 | 45.65 |

该子集总会包含样本最优系统，这里是系统 \(i = 4\)（氧化）。为避免被筛选掉，其他系统的样本均值必须足够小，以满足 \(\bar {Y} _ {\cdot i} \leq \bar {Y} _ {\cdot j} + \operatorname* {m a x} \{W_{i j} - \epsilon \}\) 对所有 \(j \neq i ,\)，特别是满足 \(\bar {Y} . i \le \bar {Y} . 4 + \operatorname* {m a x} \{0 , W_{i 4} - \epsilon \}\)。表 4 给出了阈值 \({\bar {Y}} _ {\cdot 4} +\) max \(\{0 , W_{i 4} - \epsilon \}\)；只有样本均值为 45.70 小时的备选方案 2 小到足以存活，因此 \({\mathcal {S}} = \{2 , 4 \}\)。为了确定备选方案 2 和 4 的第二阶段样本量，我们需要 \(h = h ( R_{0} , K , 1 -\)\(\alpha / 2 ) = h ( 10 , 8 , 0 . 975 ) = 4 . 635\)，它可以由 SimulationTools.xls 获得。于是 \(R_{2} = \lceil ( h S_{2} / \epsilon ) ^ {2} \rceil = 349\)，\(R_{4} = \lceil ( h S_{4} / \epsilon ) ^ {2} \rceil = 90\)。因此，备选方案 2 需要追加 \(349 - 10 = 339\) 次重复，而备选方案 4 需要追加 \(90 - 10 = 80\) 次重复。注意备选方案 2 的重复次数相当大。一般而言，如果样本方差大（使得差异难以探测），或者 \(\epsilon\) 很小（连很小的差异也要计较，如本例所示），第二阶段样本量就会很大。

在获得追加重复之后，\(\bar {\bar {Y}} _ {\cdot 2} = 45 . 88\)，\(\bar {\bar {Y}} _ {\cdot 4} = 45 . 00\)；因此选择备选方案 4，即氧化。该程序保证：以 95% 的置信度，\(\theta_{4}\) 是 8 个备选方案中最小的稳态平均周期时间；或者，如果它不是，那么它与最佳值相差在 \(\epsilon = 0 . 15\) 小时以内。

这里所介绍的选择最佳方案程序同时执行筛选与选择。这使它特别适用于从大量备选方案中寻找最佳系统设计，因为其中许多方案可能在第一阶段就被筛选掉。如果目标是淘汰那些没有竞争力的设计，该程序也可以只用于筛选而不作选择。当我们希望基于均值性能之外的准则，从一组有竞争力的备选方案中作出选择时，这样做可能是值得的。若只作筛选，程序在第 3 步停止，并可使用临界值 

$$
t = t _ {1 - (1 - \alpha) ^ {\frac {1}{k - 1}}, R_{0} - 1}
$$

 它给出的筛选稍紧一些。

类似地，如果系统设计个数 K 较小，且目标是选出最佳者，那么可以跳过筛选，让所有系统都获得第二阶段重复。若只作选择，则跳过第 3 步，并使用稍小一些的临界值 \(h = h ( R_{0} , K , 1 - \alpha )\)。

## 3 元建模

假设存在一个仿真输出响应变量 Y，它与 k 个自变量，比如说 \(x_{1} , x_{2} , \ldots , x_{k}\) 相关。因变量 Y 是一个随机变量，而自变量 \(x_{1} , x_{2} , \ldots , x_{k}\) 称为设计变量（design variable），通常受控。变量 Y 与 x 之间的真实关系由那个往往相当复杂的仿真模型所表示。我们的目标是用一个更简单的数学函数——称为元模型（metamodel）——来近似这一关系。在某些情况下，分析人员知道 Y 与 \(x_{1} , x_{2} , \ldots , x_{k}\) 之间函数关系的确切形式，比如说 

$$
Y = f ( x_{1} , x_{2} , \dots , x_{k} )
$$

。然而在大多数情况下，函数关系是未知的，分析人员必须选取一个含有未知参数的合适的 f，然后从一组数据 \(( Y , \mathbf{x} )\) 中估计这些参数。回归分析（regression analysis）就是估计这些参数的一种方法。

### 例 4

一家保险公司承诺：每天收到的所有索赔都在次日结束之前处理完毕。该公司为其拟建的索赔处理系统建立了一个仿真模型，用以评估兑现这一承诺有多困难。每天需要处理的索赔的实际数量和类型各不相同，而且数量可能随时间增长。因此，该公司希望有一个模型，能把总处理时间预测为所收到索赔数量的函数。

元模型的主要价值在于使"如果……会怎样"（what-if）类问题易于回答，例如在索赔数量为 x 时处理时间会是多少。在若干个 x 值上计算函数 f（或许还包括它的导数），通常要比对每一个取值都做一次仿真实验容易得多。

### 3.1 简单线性回归

假设希望估计单个自变量 x 与因变量 Y 之间的关系，并假设 Y 与 x 之间的真实关系被推测为线性的。在数学上，对于给定的 x 值，Y 的期望值被假定为 

$$
E (Y | x) = \beta_{0} + \beta_{1} x\tag{17}
$$

 其中 \(\beta_{0}\) 是 Y 轴上的截距，一个未知常数；\(\beta_{1}\) 是斜率，即 x 变化一个单位时 Y 的变化量，也是一个未知常数。进一步假定 Y 的每一次观测都可以由如下模型描述 

$$
Y = \beta_{0} + \beta_{1} x + \varepsilon\tag{18}
$$

 其中 ε 是一个均值为零、方差为常数 \(\sigma^{2}\) 的随机误差。由式 (18) 给出的回归模型只含单个变量 x，通常称为简单线性回归模型（simple linear regression model）。

假设有 n 对观测值 \(( Y_{1} , x_{1} ) , ( Y_{2} , x_{2} ) , \ldots , ( Y_{n} , x_{n} )\)。可以用这些观测值来估计式 (18) 中的 \(\beta_{0}\) 和 \(\beta_{1}\)。通常采用最小二乘法（method of least squares）来构造估计。在最小二乘法中，对 \(\beta_{0}\) 和 \(\beta_{1}\) 的估计要使观测值与回归线之间偏差的平方和最小。式 (18) 中的各次观测可写为 

$$
Y_{i} = \beta_{0} + \beta_{1} x_{i} + \varepsilon_{i}, \qquad i = 1, 2, \dots , n\tag{19}
$$

 其中 \(\varepsilon_{1} , \varepsilon_{2} , \ldots\) 被假定为不相关的随机变量。

式 (19) 中的每一个 \(\varepsilon_{i}\) 由下式给出 

$$
\varepsilon_{i} = Y_{i} - \beta_{0} - \beta_{1} x_{i}\tag{20}
$$

 它表示观测到的响应 \(Y_{i}\) 与由式 (17) 中模型所预测的期望响应 \(\beta_{0} + \beta_{1} x_{i}\) 之间的差。图 6 展示了 \(\varepsilon_{i}\) 与 \(x_{i} , Y_{i} ,\) 和 $ E ( Y_{i} | x_{i} ) $ 之间的关系。

![](img/ch12/img_adbdefaa7ad105bb.jpg)

**图 6 \(\varepsilon_{i}\) 与 \(x_{i} , Y_{i} ,\) 及 $ E ( Y_{i} | x_{i} ) $ 之间的关系。**

式 (20) 所给出的偏差平方和为 

$$
L = \sum_{i = 1} ^ {n} \varepsilon_{i} ^ {2} = \sum_{i = 1} ^ {n} (Y_{i} - \beta_{0} - \beta_{1} x_{i}) ^ {2}\tag{21}
$$

 而 \(L\) 称为最小二乘函数（least-squares function）。把 \(Y_{i}\) 改写如下会比较方便：

$$
Y_{i} = \beta_{0} ^ {\prime} + \beta_{1} (x_{i} - \bar {x}) + \varepsilon_{i}\tag{22}
$$

 其中 \(\beta_{0} ^ {\prime} = \beta_{0} + \beta_{1} \bar {x}\)，\(\bar {x} = \textstyle \sum_{i = 1} ^ {n} x_{i} / n\)。式 (22) 常称为变换后的线性回归模型（transformed linear regression model）。利用式 (22)，式 (21) 变为 

$$
L = \sum_{i = 1} ^ {n} [ Y_{i} - \beta_{0} ^ {\prime} - \beta_{1} (x_{i} - \bar {x}) ] ^ {2}
$$

 为使 \(L\) 最小，求 \(\partial L / \partial \beta_{0} ^ {\prime}\) 和 \(\partial L / \partial \beta_{1}\)，令二者为零，解出 \(\widehat {\beta} _ {0} ^ {\prime}\) 和 \({\widehat {\beta}} _ {1}\)。取偏导数并令其为零，得到 

$$
n \widehat {\beta_{0} ^ {\prime}} = \sum_{i = 1} ^ {n} Y_{i}
$$


$$
\widehat {\beta} _ {1} \sum_{i = 1} ^ {n} (x_{i} - \bar {x}) ^ {2} = \sum_{i = 1} ^ {n} Y_{i} (x_{i} - \bar {x})\tag{23}
$$

 方程组 (23) 常称为正规方程组（normal equations），其解为 

$$
\widehat {\beta} _ {0} ^ {\prime} = \bar {Y} = \sum_{i = 1} ^ {n} \frac {Y_{i}}{n}\tag{24}
$$

 以及 

$$
\widehat {\beta} _ {1} = \frac {\sum_{i = 1} ^ {n} Y_{i} (x_{i} - \bar {x})}{\sum_{i = 1} ^ {n} (x_{i} - \bar {x}) ^ {2}}\tag{25}
$$

 为便于计算，式 (25) 的分子改写为 

$$
S_{x y} = \sum_{i = 1} ^ {n} Y_{i} (x_{i} - \bar {x}) = \sum_{i = 1} ^ {n} x_{i} Y_{i} - \frac {\left(\sum_{i = 1} ^ {n} x_{i}\right) \left(\sum_{i = 1} ^ {n} Y_{i}\right)}{n}\tag{26}
$$

 其中 \(S_{x y}\) 表示 x 与 Y 的修正交叉乘积和。为便于计算，式 (25) 的分母改写为 

$$
S_{x x} = \sum_{i = 1} ^ {n} (x_{i} - \bar {x}) ^ {2} = \sum_{i = 1} ^ {n} x_{i} ^ {2} - \frac {\left(\sum_{i = 1} ^ {n} x_{i}\right) ^ {2}}{n}\tag{27}
$$

 其中 \(S_{x x}\) 表示 x 的修正平方和。\({\widehat {\beta}} _ {0}\) 的值可容易地由下式还原 

$$
\widehat {\beta_{0}} = \widehat {\beta_{0} ^ {\prime}} - \widehat {\beta_{1}} \bar {x}\tag{28}
$$


$$
\widehat {\boldsymbol {\beta}} _ {0}
$$


$$
\widehat {\beta} _ {1}
$$

 例 4 中索赔处理系统的仿真模型在前一天收到 x = 100、150、200、250 和 300 件索赔的初始条件下运行。每种设定下取得 3 次重复。响应 Y 是处理 x 件索赔所需的小时数。结果如表 5 所示。所收到索赔数量与总处理时间之间的图形关系如图 7 所示。考察这一散点图可知，索赔数量与处理时间之间存在很强的关系。对由式 (18) 给出的线性模型所作的初步假定看来是合理的。

**表 5 给定 x 件索赔时处理时间的仿真结果**

| 索赔数量 x | 处理时间（小时）Y |
|---|---|
| 100 | 8.1 |
| 100 | 7.8 |
| 100 | 7.0 |
| 150 | 9.6 |
| 150 | 8.5 |
| 150 | 9.0 |
| 200 | 10.9 |
| 200 | 13.3 |
| 200 | 11.6 |
| 250 | 12.7 |
| 250 | 14.5 |
| 250 | 14.7 |
| 300 | 16.5 |
| 300 | 17.5 |
| 300 | 16.3 |

![](img/ch12/img_2ca68684a51ac09e.jpg)

**图 7 索赔数量与处理时间（小时）之间的关系。**

以处理时间作为 \(Y_{i}\) 值（因变量），以索赔数量作为 \(x_{i}\) 值（自变量），\(\widehat {\beta} _ {0}\) 和 \({\widehat {\beta}} _ {1}\) 可由下列计算求得：

$$
n = 15
$$


$$
\sum_{i = 1} ^ {15} x_{i} = 3000
$$


$$
\sum_{i = 1} ^ {15} Y_{i} = 178
$$


$$
\sum_{i = 1} ^ {15} x_{i} ^ {2} = 675, 000
$$


$$
\sum_{i = 1} ^ {15} x_{i} Y_{i} = 39080
$$

 以及 \(\bar {x} = 3000 / 15 = 200\)。由式 (26)，\(S_{x y}\) 计算为 

$$
S_{x y} = 39, 080 - \frac {(3000) (178)}{15} = 3480
$$

 由式 (27)，\(S_{x x}\) 计算为 

$$
S_{x x} = 675, 000 - \frac {(3000) ^ {2}}{15} = 75, 000
$$

 然后，由式 (25) 计算 \({\widehat {\beta}} _ {1}\) 为 

$$
\widehat {\beta} _ {1} = \frac {S_{x y}}{S_{x x}} = \frac {3480}{75 , 000} = 0. 0464
$$

 如式 (24) 所示，\(\widehat {\beta_{0} ^ {\prime}}\) 就是 \({\bar {Y}} ,\) 即 

$$
\widehat {\beta_{0} ^ {\prime}} = \frac {178}{15} \approx 11. 8667
$$

 为了用原始的项来表达模型，由式 (28) 计算 \({\widehat {\beta}} _ {0}\) 为 

$$
\widehat {\beta_{0}} = 11. 8667 - 0. 0464 (200) = 2. 5867
$$

 于是，给定 x 时 Y 的均值 $ E ( Y | x ) $ 的估计由下式给出 

$$
\widehat {y} = \widehat {\beta_{0}} + \widehat {\beta_{1}} x = 2. 5867 + 0. 0464 x\tag{29}
$$

 对于给定的索赔数量 x，该模型可用于预测处理它们所需的小时数。系数 \({\widehat {\beta}} _ {1}\) 的解释是：每多收到一件索赔，期望总处理时间就增加约 0.0464 小时，即 2.8 分钟。

回归分析被广泛使用，也经常被误用。这里简要提及若干常见的误用之处。以式 (29) 的方式导出的关系，只对处于原始数据范围之内的自变量取值有效。初步假定的线性关系在原始范围之外可能并不成立。事实上，我们从排队论知道：当索赔数量接近系统容量时，平均处理时间可能迅速上升。因此，式 (29) 只能认为对 \(100 \leq x \leq 300\) 有效。不建议把回归模型用于外推。

在选择变量时应当谨慎，要让它们之间具有貌似合理的因果关系。完全有可能建立起在实际意义上毫不相干的统计关系。例如，也许有人试图把某钢铁厂的月产量与该月出现在某经理办公桌上的报告重量联系起来。一条直线对这些数据而言可能看似提供了很好的模型，但这两个变量之间的关系是脆弱的。观察到的强关系并不意味着变量之间存在因果关系。只有当分析揭示了其存在的某些合理原因时，才能推断出因果关系。在例 4 中，"从更多的索赔开始意味着需要更多时间来处理它们"是合理的。因此，式 (29) 这种形式的关系至少是貌似合理的。

### 回归的显著性检验

在第 3.1 节中，曾假定 Y 与 x 之间存在线性关系。在例 4 中，为在计算 \(\widehat {\beta} _ {0}\) 和 \({\widehat {\beta}} _ {1}\) 之前评估线性模型是否是一个合理的初步假定，绘制了如图 7 所示的、联系索赔数量与处理时间的散点图。然而，在利用该模型由自变量 \(x_{i}\) 预测响应 \(Y_{i}\) 之前，应当检验简单线性关系的充分性。可以进行若干检验来帮助确定模型的充分性。检验所初步假定的模型阶次是否正确，通常称为失拟检验（lack-of-fit testing）。

回归的显著性检验提供了评估模型充分性的另一种手段。特别地，检验所估计的回归系数 \(\widehat {\beta} _ {0}\) 和 \({\widehat {\beta}} _ {1}\) 是否在统计上不同于 0 是必不可少的；对于下面将要介绍的多元线性回归而言尤其如此——那时元模型中包含大量的设计变量 x，但并不确定它们是否都具有实质性的解释作用。

显著性与失拟的检验包含在所有回归软件之中。关于这些检验的细节，参见 Box 与 Draper [1987]、Hines 等人 [2002] 以及 Montgomery [2000]。

### 多元线性回归

如果第 3.1 节中的简单线性回归模型不充分，还存在若干其他可能。可能存在多个自变量，使关系具有如下形式 

$$
Y = \beta_{0} + \beta_{1} x_{1} + \beta_{2} x_{2} + \dots + \beta_{m} x_{m} + \varepsilon\tag{30}
$$

 注意，这个模型仍然是线性的，但含有多于一个自变量。具有式 (30) 所示形式的回归模型称为多元线性回归模型（multiple linear regression model）。另一种可能是模型具有二次形式，例如 

$$
Y = \beta_{0} + \beta_{1} x + \beta_{2} x^{2} + \varepsilon\tag{31}
$$

 式 (31) 也是一个线性模型，令 \(x_{1} = x\) 和 \(x_{2} = x^{2}\) 即可把它变换为式 (30) 的形式。还有一种可能，模型具有如下形式 

$$
Y = \beta_{0} + \beta_{1} x_{1} + \beta_{2} x_{2} + \beta_{3} x_{1} x_{2} + \varepsilon
$$

 它同样是一个线性模型。对这三种形式及相关模型的分析，可见 Box 与 Draper [1987]、Hines 等人 [2002]、Montgomery [2000] 以及其他应用统计学教材；也可参见 Kleijnen [1987, 1998]，后者主要关注这些模型在仿真中的应用。

### 3.2 元建模与计算机仿真

第 3.1 节的内容适用于一般的回归元建模，而不仅仅适用于计算机仿真。然而，当目标是为一次计算机仿真提供元模型时，会出现一些特殊的问题与机遇；下面对此加以说明。

### 回归中的随机数分配

随机数种子或流的分配是仿真实验设计的一部分。（本节基于 Nelson [1992]。）给不同的设计点（多元线性回归中 \(x_{1} , x_{2} , \ldots , x_{m}\) 的设定）分配不同的种子或流，可以保证来自不同设计点的响应 Y 在统计上独立。同理，给不同的设计点分配相同的种子或流，会由于它们都具有相同的随机性来源而在相应响应之间诱发相关性。

许多教科书中的实验设计都假定各设计点之间的响应独立。为符合这一假定，我们必须给每个设计点分配不同的种子或流。然而，把相同的随机数种子或流分配给所有设计点——换句话说，使用公共随机数——往往是有用的。

在元模型中使用公共随机数的直觉是：如果各设计点受到相同的实验条件，特别是相同的随机性来源的作用，那么就能在各设计点之间实现更公平的比较。其数学依据如下：假设我们拟合简单线性回归 \(Y_{i} = \beta_{0} + \beta_{1} x_{i} + \varepsilon_{i}\) 并得到最小二乘估计 \(\widehat {\beta} _ {0}\) 和 \({\widehat {\beta}} _ {1}\)。那么，设计点 i 与 j 之间性能期望之差的一个估计量为 

$$
\widehat {\beta} _ {0} + \widehat {\beta} _ {1} x_{i} - (\widehat {\beta} _ {0} + \widehat {\beta} _ {1} x_{j}) = \widehat {\beta} _ {1} (x_{i} - x_{j})
$$

 当 \(x_{i}\) 和 \(x_{j}\) 是固定的设计点时，\({\widehat {\beta}} _ {1}\) 决定了设计点 i 与 \(j\) 之间（或者就此而言，x 的任意其他两个取值之间）的估计差异。因此，可以预期公共随机数会降低 \({\widehat {\beta}} _ {1}\) 的方差，更一般地说，会降低多元线性回归中所有斜率项的方差。公共随机数通常不会降低截距项 \(\widehat {\beta} _ {0}\) 的方差。无论是否使用公共随机数，最小二乘估计量 \({\widehat {\beta}} _ {0}\) 和 \({\widehat {\beta}} _ {1}\) 都是恰当的，但与之相联系的统计分析会受这一选择的影响。关于公共随机数下元模型的统计分析，参见 Kleijnen [1988] 和 Nelson [1992]。

### 排队仿真

许多计算机仿真模型所针对的系统包含排队或排队网络，这些仿真的元模型可能包含与容量（服务台数量和服务率）以及顾客负荷（到达率）有关的设计变量。这类模型违背了第 3.1 节所描述的标准回归假定。一个小例子可以说明这些困难和可能的解决办法。

回顾 \(M / M / 1\) 排队：一个单服务台、先来先服务的排队，到达过程为速率 λ 顾客/时间的泊松过程，服务时间服从速率 \(\mu\) 顾客/时间的指数分布。令 \(x = \lambda / \mu\)，即所谓利用率 \(\rho\)。那么，当 \(0 < x < 1\) 时，队列中顾客的稳态期望数为 

$$
L = \frac {x}{1 - x} = x + x^{2} + x^{3} + \dots\tag{32}
$$

 现在假设我们并不知道关系式 (32)，而试图通过运行 \(M / M / 1\) 排队的仿真来为 \(0 < x < 1\) 上的 L 建立一个元模型。显然，通常的低阶多项式模型（只含上面提到的 x 和 \(x^{2}\) 项）拟合得不会很好，特别是当 x 接近 1 时，此时高阶项变得重要。

这一问题至少有两种解决办法。一种是把设计空间（本例中为 \(0 < x < 1\)）分解成若干较小的片段，使得低阶多项式在每一片段上是充分的。第二种办法是拟合一个由排队论所启发的非线性模型。Yang、Ankenman 与 Nelson [2007] 指出 

$$
Y = \frac {\sum_{i = 0} ^ {m} \beta_{i} x^{i}}{(1 - x) ^ {p}} + \varepsilon\tag{33}
$$

 其中 \(m \ \leq \ 3\)，对许多排队输出 Y（包括队列平均人数和队列平均时间）都效果良好。注意，未知参数是 m、\(\beta_{i}\) 和 \(p ,\)，因此需要非线性最小二乘法来拟合它们（Bates 与 Watts [1988]）。当 Y 是队列中平均人数（即 L 的一个估计量）时，\(m = 1 , \beta_{0} = 0 , \beta_{1} = 1\)，且 \(p = 1\)。

遗憾的是，问题还要更深一层。如果我们通过做一次长度为 \(T\) 时间单位的仿真运行来估计 L，并记录队列中的平均顾客数（记为 Y），那么可以证明（Whitt [1989]） 

$$
\operatorname{var} (Y) \approx \frac {2 x (1 + x)}{(1 - x) ^ {4} T}\tag{34}
$$

 因此，方差并非如回归中所假定的那样在整个设计空间上保持恒定，而是随着 x 趋近于 1 而增大，且是爆炸式地增大。其结果是，一次长度为 \(T = 100\) 时间单位的运行在 \(x = 0 . 1\) 处可能给出 L 的相当好的估计，但在 \(x = 0 . 9\) 处给出的估计则非常差。当响应方差强烈依赖于 x 时，回归软件所产生的标准统计推断将是无效的，这意味着我们可能察觉不到失拟或拟合不良。

针对回归中响应方差不相等的一种补救办法，是对数据施加一个方差稳定化变换，使方差在整个设计空间上更接近恒定。在仿真中，我们可以直接处理这一问题：根据 x 的不同而采用不同的运行长度 \(T\) 或不同的重复次数。

我们建议把 \(M / M / 1\) 的结果用作更一般的排队网络的粗略指南。若令 x 表示网络利用率，即负荷最重的那个队列的利用率，则近似式 (34) 意味着：在两个设计点 \(x_{1}\) 和 \(x_{2}\) 处的运行长度 \(T_{1}\) 和 \(T_{2}\) 应满足 

$$
\frac {T_{2}}{T_{1}} = \frac {x_{2}}{x_{1}} \left(\frac {1 + x_{2}}{1 + x_{1}}\right) \left(\frac {1 - x_{1}}{1 - x_{2}}\right) ^ {4}
$$

 例如，如果我们要在利用率 \(x_{1} = 0 . 1\) 处仿真 \(T_{1} = 1\) 个时间单位，那么为了稳定方差，在 \(x_{2} = 0 . 3\) 处就大约需要 \(T_{2} = 10\) 个时间单位，在 \(x_{2} = 0 . 5\) 处需要 \(T_{2} = 72\) 个时间单位，在 \(x_{2} = 0 . 7\) 处需要 \(T_{2} = 876\) 个时间单位，在 \(x_{2} = 0 . 9\) 处需要 \(T_{2} = 101\), 944 个时间单位。

其核心信息是：在排队仿真中，如果元模型包含影响顾客负荷的决策变量——通过改变到达率、服务率或服务台数量——那么低阶多项式可能不是元模型的好选择；而且若不施行运行长度控制，就不能假定响应方差在整个设计空间上相等。

## 4 基于仿真的优化

考虑下面的例子。（其中部分描述基于 Boesel、Nelson 与 Ishii [2003]。）

### 例 6：物料搬运系统（MHS）

工程师需要设计一个 MHS，它由一台大型自动化存储与检索设备、自动导引车（AGV）、AGV 站点、升降机和输送带组成。他们可以控制的设计变量包括 AGV 的数量、每辆 AGV 的载重，以及用于调度 AGV 的路由算法。备选设计将根据 AGV 利用率、需要搬运的物料的运输延迟，以及总体投资与运营成本来评价。

### 例 7：液化天然气（LNG）运输

一个 LNG 运输系统由 LNG 运输船以及装载、卸载和储存设施组成。为了使成本最小，设计者可以控制运输船的尺寸、投入使用的运输船数量、装卸设施的泊位数量以及储罐的容量。

### 例 8：汽车发动机装配

在一条装配线上，工作站之间的大缓冲区（队列）可以提高工作站的利用率——因为往往会有东西等着被加工——但会推高对空间的需求和在制品库存。我们希望找到一种缓冲区容量的分配方案，使这些相互竞争的成本之和最小。

### 例 9：交通信号配时

土木工程师希望对某繁忙路段沿线的交通信号进行排序，以减少驾驶员延误以及狭窄交叉街道上发生的拥堵。对每一个交通信号，红灯、绿灯和绿色转向箭头信号的周期长度都可以分别设定。

### 例 10：在线服务

一家通过互联网提供在线信息服务的公司正在把其计算机体系结构从集中式大型机改为分布式工作站计算。CPU 的数量与类型、网络结构以及处理任务的分配都需要作出选择。对客户查询的响应时间是关键的性能度量。

这些设计问题有什么共同点？显然，仿真模型在每一个问题中都可能有用，而且它们都有一个隐含的目标：相对于某些性能度量（成本、延误等）找到最佳设计。在每一个例子中，备选设计的数量可能都非常庞大，从几十个到几千个不等，肯定多于我们在第 2.2 节中所考虑的 2 到 100 个。其中有些例子包含种类繁杂的决策变量：离散的（AGV 的数量、CPU 的数量）、连续的（运输船的尺寸、红灯周期的长度）和定性的（路由策略、分配处理任务的算法）。这使得如第 3 节所述的建立元模型变得困难。

所有这些问题都属于"基于仿真的优化"（optimization via simulation）这一总主题之下，其目标是最小化或最大化某些系统性能度量，而系统性能只能通过运行计算机仿真来评价。基于仿真的优化是一个相对较新但已经十分庞大的主题，商业软件也已得到广泛使用。在本节中，我们描述开展基于仿真的优化时应当考虑的关键问题，给出一些可查阅文献的指针，并给出一个示例算法。

### 4.1 "基于仿真的优化"意味着什么？

优化是运筹学与管理科学工作者使用的一个关键工具，对许多类别的问题都有成熟的算法，其中最著名的是线性规划。优化方面的多数工作处理的是这样一些问题：系统的所有方面都被视为确定已知的；最关键的是，任何设计的性能（成本、利润、制造周期等）都可以被精确评价。

在随机的离散事件仿真中，任何一次仿真运行的结果都是一个随机变量。为便于记号表达，令 \(x_{1} , x_{2} , \ldots , x_{m}\) 为 m 个可控的设计变量，令 \(Y ( x_{1} , x_{2} , \dots , x_{m} )\) 为一次运行中所观察到的仿真输出性能。具体地说，在例 6 中，\(x_{1} , x_{2} , x_{3}\) 可以分别表示 AGV 的数量、每辆 AGV 的载重以及用于调度 AGV 的路由算法，而 \(Y ( x_{1} , x_{2} , x_{3} )\) 可以是 MHS 的购置与运营总成本。

相对于 \(x_{1} , x_{2} , \dots , x_{m}\) 去"优化 \(Y ( x_{1} , x_{2} , \ldots , x_{m} )\)"意味着什么？由于 Y 是一个随机变量，我们无法优化 Y 的实际取值。最常见的优化定义是 

$$
\text{maximize   or   minimize} E \left(Y (x_{1}, x_{2}, \dots , x_{m})\right)\tag{35}
$$

 换句话说，最大化或最小化性能的数学期望，即长期平均值。这是我们所知的所有商业软件包使用的默认优化定义。在我们的例子中，\(E \left( Y ( x_{1} , x_{2} , x_{3} ) \right)\) 是在有 \(x_{1} \mathrm{A G V s} , x_{2}\) 每辆 \(\operatorname{A G V}\) 的载重以及路由算法 \(x_{3}\) 的条件下，运营该 MHS 的期望（即长期平均）成本。

应当指出，式 (35) 并不是唯一可能的定义。例如，我们可能希望选择那个最有把握使购置与运营成本低于 \$ D 的 MHS 设计，从而把目标改为 

$$
\text{maximize} \operatorname * {P r} (Y (x_{1}, x_{2}, x_{3}) \leq D)
$$

 我们可以通过定义一个新的性能度量 

$$
Y^{\prime} (x_{1}, x_{2}, x_{3}) = \left\{\begin{array}{l l} 1, & \text{if} Y (x_{1}, x_{2}, x_{3}) \leq D \\ 0, & \text{otherwise} \end{array} \right.
$$

 并转而最大化 \(E \left( Y^{\prime} ( x_{1} , x_{2} , x_{3} ) \right)\)，把这个目标纳入式 (35) 的形式。

当我们希望选择那个最可能成为最佳方案的系统设计时，就出现了一个更复杂的优化问题。当所关心的是一次性（one-shot）而非长期平均性能时，这样的目标是恰当的。例子包括一次航天飞机的发射，或某笔独特的、大批量的订单的交付。Bechhofer、Santner 与 Goldsman [1995] 在"多项选择"（multinomial selection）这一主题下讨论了这一问题。

我们一直假定一个系统设计 \(x_{1} , x_{2} , \ldots , x_{m}\) 可以用单一的性能度量 Y（例如成本）来评价。显然，情况未必总是如此。在 MHS 的例子中，我们可能也对某种系统生产率的度量（例如吞吐量或周期时间）感兴趣。目前，基于仿真的多目标优化还不成熟。因此，通常采用下列三种策略之一：

1. 把所有性能度量合并成一个度量，最常见的是成本。例如，MHS 中每件完工产品所产生的收入可以代表生产率，并作为一项负成本计入。

2. 针对一个关键性能度量进行优化，然后对排名靠前的解按次要性能度量进行评价。例如，可以针对期望成本优化 MHS，然后对前 5 个设计比较其周期时间。这种方法要求保留的信息不仅仅是最佳解。

3. 针对一个关键性能度量进行优化，但只考虑那些在其他性能度量上满足某些约束的备选方案。例如，可以针对期望成本优化 MHS，但只限于那些期望周期时间低于给定阈值的备选方案。

最后，系统设计变量 \(x_{1} , x_{2} , \ldots , x_{m}\) 通常还存在约束。在例 10 的在线服务中，会有一个预算约束限制所能购买的工作站数量。复杂的、尤其是非线性的约束会使可行空间更难以搜索。在某些情况下，约束和目标函数一样无法被确定地评价，而必须由仿真来估计。例如，这家在线服务公司可能希望最小化总体客户响应时间，但他们可能必须在"某个特定（比如说政府）客户的平均响应时间低于某一水平"这一约束之下来做到这一点。

### 4.2 为什么基于仿真的优化很困难？

即使不存在不确定性，如果设计变量很多，如果问题包含种类繁杂的设计变量类型，或者如果对性能函数的结构知之甚少，优化也可能非常困难。基于仿真的优化又增加了一重复杂性：某一个特定设计的性能无法被精确评价，而只能被估计。因为我们只有估计值，所以无法有把握地断定一个设计优于另一个设计，而这种不确定性会挫败那些试图朝改进方向推进的优化算法。原则上，可以通过在每个设计点做极多次的重复或极长的运行，使性能估计实质上没有方差，从而消除这一复杂性。但在实践中，这可能意味着由于仿真每一个设计都要花费时间，因而只能探索极少的备选设计。

抽样变异性的存在迫使基于仿真的优化作出折中。以下是几种标准的折中：

1. 保证预先规定的正确选择概率。第 2.2 节中的选择最佳方案程序就是这类方法的一个例子，它允许分析人员规定希望有多大机会做对。这类算法通常要么要求仿真每一个可能的设计，要么要求各设计之间存在很强的函数关系（如元模型）。其他算法可见 Goldsman 与 Nelson [1998] 以及 Kim 与 Nelson [2006]。

2. 保证渐近收敛。有许多算法可以保证：当仿真工作量（重复次数、重复长度）趋于无穷时，算法收敛到全局最优解。这些保证是有用的，因为它们表明该算法往往会到达分析人员希望它去的地方。然而，收敛可能很慢，而且当算法在有限时间内终止时（在实践中必定如此），对所报告的解究竟有多好，往往没有任何保证。关于适用于离散或连续变量问题的具体算法，参见 Andradóttir [1998]。

3. 对确定性对应问题最优。这里的思路是：使用一种算法，如果每个设计的性能都能够被确定地评价，那么该算法就能找到最优解。一个例子可能是把标准的非线性规划算法应用于仿真优化问题。通常需要由分析人员来确保投入了足够的仿真工作量（重复次数或运行长度），以保证这类算法不被抽样变异性所误导。不建议把一种假定确定性评价的算法直接应用于随机仿真。

4. 鲁棒启发式。许多启发式方法是针对确定性优化问题发展起来的，它们并不保证找到最优解，但已被证明在困难的实际问题上非常有效。其中有些启发式把随机性作为其搜索策略的一部分，因此可以说它们对抽样变异性不像其他类型的算法那么敏感。尽管如此，仍然重要的是要确保投入足够的仿真工作量（重复次数或运行长度），以保证这类算法不被抽样变异性所误导。

鲁棒启发式是商业化的基于仿真的优化软件中最常见的算法。我们在下一节中就其使用给出一些指导。关于优化理论与实践的全面讨论见 Fu [2002]；关于鲁棒启发式的更多内容见 Olafsson [2006]。

### 4.3 使用鲁棒启发式

所谓"鲁棒启发式"（robust heuristic），我们指的是这样一种程序：它不依赖强的问题结构（例如 \(E \left( Y ( x_{1} , \ldots , x_{m} ) \right)\) 的连续性或凸性）就能有效，可以应用于具有混合类型决策变量的问题，并且在理想情况下能够容忍一定的抽样变异性。遗传算法（genetic algorithm, GA）和禁忌搜索（tabu search, TS）是两个突出的例子，但还有许多其他算法及其众多变体。这类启发式构成了大多数商业实现的核心。为了让读者对这些启发式有所感受，下面介绍 GA 与 TS。我们要提醒读者：这里只给出每种程序最简单版本的高层次描述，商业实现要复杂得多。

假设基于仿真的优化问题有 \(k\) 个可能的解。令 \(\mathbf{X} = \{\mathbf{x} _ {1} , \mathbf{x} _ {2} , \ldots , \mathbf{x} _ {k} \}\) 表示这些解，其中第 \(i\) 个解 \(\mathbf{x} _ {i} = ( x_{i 1} , x_{i 2} , \dots , x_{i m} )\) 给出 m 个决策变量的具体设定。解 \(\mathbf{x} _ {i}\) 处的仿真输出记为 \(Y ( \mathbf{x} _ {i} )\)，它可以是一次重复的输出，也可以是若干次重复的平均。我们的目标是找到使 \(E \left( Y ( \mathbf{x} ) \right)\) 最小的解 \(\mathbf{x} ^ {*}\)。

在每一次迭代（称为一代，generation）中，GA 对一个由 \(p\) 个解组成的"种群"（population）进行操作。把第 \(j\) 次迭代上的解的种群记为 \(\mathbf{P} ( j ) = \{\mathbf{x} _ {1} ( j ) , \mathbf{x} _ {2} ( j ) , \ldots , \mathbf{x} _ {p} ( j ) \}\)。\(\mathbf{P} ( j )\) 中可能有同一个解的多份副本，也可能包含在先前迭代中所发现的解。从一次迭代到下一次迭代，这个种群以如下方式演化：好的解倾向于存活并产生新的、有望更好的解，而劣解倾向于被从种群中移除。下面给出基本的 GA。

### 基本遗传算法（GA）

步骤 1. 置迭代计数器 \(j = 0\)，并（也许是随机地）选择一个由 \(p\) 个解组成的初始种群 \(\mathbf{P} ( 0 ) = \{\mathbf{x} _ {1} ( 0 ) , \ldots , \mathbf{x} _ {p} ( 0 ) \}\)。

步骤 2. 运行仿真实验，得到 \(\mathbf{P} ( j )\) 中所有 \(p\) 个解 \(\mathbf{x} ( j )\) 的性能估计 \(Y ( \mathbf{x} )\)。

步骤 3. 从 \(\mathbf{P} ( j )\) 中的解中选择一个由 \(p\) 个解组成的种群，使得 \(Y ( \mathbf{x} )\) 值较小的解更可能（但并非必然）被选中。把这个解的种群记为 \(\mathbf{P} ( j + 1 )\)。

步骤 4. 通过交叉（crossover，它把两个解 \(\mathbf x_{i} ( j + 1 )\) 和 \(\mathbf{x} _ {\ell} ( j + 1 )\) 的各一部分连接起来，形成一个新的解）与变异（mutation，它随机地改变一个解 \(\mathbf{x} _ {i} ( j + 1 ) ,\) 的一部分），对 \(\mathbf{P} ( j + 1 )\) 中的解进行重组。

步骤 5. 令 \(j = j + 1\)，转到步骤 2。

GA 可以在经过规定次数的迭代之后终止，可以在种群中几乎看不到或完全看不到改进时终止，也可以在种群中包含同一个解的 \(p\) 份副本时终止。在终止时，把最后一个种群中 \(Y ( \mathbf{x} )\) 值最小的解 \(\mathbf{x} ^ {*}\) 选为最佳；或者，也可以选择所有迭代中 \(Y ( \mathbf{x} )\) 值最小的那个解。

GA 几乎适用于任何优化问题，因为选择、交叉和变异这三种操作都可以用非常通用的方式定义，而不依赖于问题的具体细节。然而，当这些操作没有针对具体问题进行调优时，一个 GA 的进展可能非常缓慢。商业版本往往是自调优的，即它们在搜索过程中更新选择、交叉和变异的参数。有一些证据表明，GA 对 \(Y ( \mathbf{x} )\) 中的抽样变异性是容忍的，因为它们维持着一个解的种群，而不是专注于改进某一个当前最优解。换句话说，GA 是否把种群中各解的排序排得完美并不关键，因为下一次迭代取决于整个种群，而不是取决于某一个解。

另一方面，禁忌搜索（TS）在每一次迭代上识别出一个当前最优解，然后试图改进它。改进通过"移动"（move）改变解来实现。例如，解 \(( x_{1} , x_{2} , x_{3} )\) 可以通过"给第一个决策变量加 1"这一移动变成解 \(( x_{1} + 1 , x_{2} , x_{3} )\)（也许 \(x_{1}\) 表示例 \(^ {6 }\) 中 AGV 的数量，那么这一移动就是再增加一辆 AGV）。解 x 的"邻居"是所有可以通过合法移动到达的解。TS 找出最优的邻居解并移动到它。然而，为了避免那些使搜索回到先前访问过的解的移动，某些移动可能在若干次迭代中成为"禁忌"（tabu，即不可使用）。从概念上讲，想想你会如何走出迷宫：如果你走过一条通向死胡同的路，你就会避免再走那条路——它将是禁忌的。

下面给出基本的 TS 算法。该描述基于 Glover [1989]。

### 基本禁忌搜索（TS）

步骤 1. 置迭代计数器 \(j = 0\)，禁忌移动列表为空。在 X 中选择一个初始解 \(\mathbf{x} ^ {*}\)，也许是随机地选择。

步骤 2. 在 \(\mathbf{X} ^ {*}\) 的所有不能通过禁忌移动到达的邻居中，找出使 \(Y ( \mathbf{x} )\) 最小的解 \(\mathbf{x} ^ {\prime}\)，运行完成这一优化所需的任何仿真。

步骤 3. 若 \(Y ( \mathbf{x} ^ {\prime} ) < Y ( \mathbf{x} ^ {*} )\)，则 \(x^{*} = x^{\prime} ;\)；把当前最优解移动到 \(x^{\prime}\)。

步骤 4. 更新禁忌移动列表，转到步骤 2。

TS 可以在完成了规定次数的迭代时终止，可以在经过了若干次迭代而 \(\mathbf{X} ^ {*}\) 未发生变化时终止，也可以在没有更多可行移动时终止。在终止时，把解 \(\mathbf{x} ^ {*}\) 选为最佳。

TS 本质上是一个针对离散决策变量的优化器，但连续决策变量可以离散化，如第 4.4 节所述。TS 积极地追寻改进的解，因此往往进展迅速。然而，它对 \(Y(\mathbf{x})\) 中的随机变异性更为敏感，因为 \(\mathbf{x}^*\) 被当作迄今为止的真正最优解，并且要在此基础上加以改进。不过，也存在一些概率版本的 TS，它们应当不那么敏感。商业 TS 实现（而基本 TS 不具备）的一个重要特征是：当这样做有利时，有一种机制可以越过禁忌表（tabu list）。

接下来，我们就使用采用 GA、TS 或其他鲁棒启发式的商业产品提出三条建议：控制抽样变异性、重启，以及清理。

### 控制抽样变异性

在许多情况下，将由用户来决定在每个候选解上进行多少抽样（重复次数或运行长度）。一般而言这是一个困难的问题。理想情况下，抽样量应当随着启发式逼近更好的解而增加，原因很简单：区分期望性能相近的解，要比区分相差很大的解困难得多。在搜索的早期，启发式可能很容易识别出好的解和搜索方向，因为明显处于劣势的解正在与好得多的解作比较；然而到了搜索的后期，情况可能就不是这样了。

如果分析人员必须指定在整个搜索过程中使用的每个解的固定重复次数，那么应当先做一次预实验。仿真若干个设计，有些位于解空间的极端处，有些靠近其中心。用第 1 节中的方法比较这些设计中表面上的最优者与表面上的最劣者。利用第 1.3 节中描述的技术，求出要宣称这些设计在统计上显著不同所需的最小重复次数。这就是所应当使用的最小重复次数。

### 重启

由于鲁棒启发式并不保证收敛到基于仿真的优化的最优解，因此有理由把优化运行两次或更多次，看哪一次运行给出最好的解。每一次优化运行都应当使用不同的随机数种子或流，并且在理想情况下应当从不同的初始解出发。尝试从解空间极端处的解、空间中心处的解以及随机生成的解开始优化。如果熟悉该系统的人怀疑某些设计会很好，一定要把它们作为启发式可能的起始解包括进来。

### 清理

在一次或多次优化运行结束之后，关键是要对启发式所识别出的排名靠前的设计再做一组实验，以避免两类错误：未能识别出搜索过程中所仿真的最佳系统设计；以及对最终所选设计性能的估计很差。之所以会出现这些错误，是因为任何优化算法都不可能在每一步都维持统计误差控制的同时还取得进展。因此，我们建议在搜索之后进行"清理"（clean up），也就是运用第 2 节中的比较技术，特别是第 2.2 节中的选择最佳方案程序，做一次严格的统计分析，以评价在搜索过程中所发现的那些设计里哪些是最佳的或接近最佳的。搜索过程中所仿真的系统设计中排名前 5% 到 10% 的那些，都应当接受这一受控实验。"清理"这一概念由 Boesel、Nelson 与 Kim [2003] 提出。

### 4.4 一个示例：随机搜索

在本节中，我们给出一种称为随机搜索（random search）的基于仿真的优化算法。具体的实现基于 Andradóttir [1998] 中的算法 2，该算法在某些条件下保证渐近收敛。因此，只要允许它运行足够长的时间，它就会找到真正的最优解。然而在实践中，收敛可能很慢，而且这一特定版本的随机搜索对内存的需求可能相当大。尽管随机搜索并不是一种鲁棒启发式，我们还将用它来演示我们在配合这类启发式时会采用的若干策略，并演示为什么即使对于看似并不复杂的算法，基于仿真的优化也很棘手。

我们所给出的随机搜索算法要求可能存在的系统设计数目有限，尽管这个数目可能很大。这似乎排除了具有连续决策变量的问题，例如输送带速度。然而在实践中，表面上的连续决策变量往往可以用合理的方式离散化。例如，如果输送带速度可以是从每分钟 60 到 120 英尺的任意值，那么把可能的输送带速度取为 \(60 , 61 , 62 , \cdots , 120\) 英尺每分钟（共 61 个可能取值），可能并不会损失什么。但要注意，也有专门为连续变量问题设计的算法（Andradóttir [1998]）。

再次地，令基于仿真的优化问题的 \(k\) 个可能解记为 \(\{\mathbf{x} _ {1} , \mathbf{x} _ {2} , \ldots , \mathbf{x} _ {k} \}\)，其中第 \(i\) 个解 \(\mathbf{x} _ {i} = ( x_{i 1} , x_{i 2} , \ldots , x_{i m} )\) 给出 m 个决策变量的具体设定。解 \(\mathbf{x} _ {i}\) 处的仿真输出记为 \(Y ( \mathbf{x} _ {i} )\)；它可以是一次重复的输出，也可以是若干次重复的平均。我们的目标是找到使 \(E \left( Y ( \mathbf{x} ) \right)\) 最小的解 \(\mathbf{x} ^ {*}\)。

在随机搜索算法的每一次迭代中，我们把当前的一个好的解与一个随机选出的竞争者作比较。如果竞争者更好，那么它就成为当前的好解。当我们终止搜索时，我们选择的是那个被访问次数最多的解——这意味着我们预期会多次回访某些解。

### 随机搜索算法

步骤 1. 初始化计数器变量 \(C ( i ) = 0 \mathrm{f o r} i = 1 , 2 , \ldots , k\)。选择一个初始解 \(i^{0} .\)，并令 \(C ( i^{0} ) = 1 .\)［提示：\(C ( i )\) 记录我们访问解 \(i\) 的次数。］

步骤 2. 从除 \(i^{0}\) 之外的所有解组成的集合中选取另一个解 \(i^{\prime}\)，使每一个解都有相等的被选中机会。

步骤 3. 在两个解 \(i^{0}\) 和 \(i^{\prime}\) 处运行仿真实验，得到输出 \(Y ( i^{0} )\) 和 \(Y ( i^{\prime} )\)。若 \(Y ( i^{\prime} ) < Y ( i^{0} )\)，则令 \(i^{0} = i^{\prime}\)。（见步骤 4 之后的注。）

步骤 4. 令 \(C ( i^{0} ) = C ( i^{0} ) + 1\)。若尚未完成，则转到步骤 2。若已完成，则选择使 \(C ( i^{*} )\) 为最大计数的 \(x _ {i^{*}}\) 作为估计的最优解。

注意，如果问题是最大化问题，则用下式替换步骤 3：

步骤 3′. 在两个解 \(i^{0}\) 和 \(i^{\prime}\) 处运行仿真实验，得到输出 \(Y ( i^{0} )\) 和 \(Y ( i^{\prime} )\)。若 \(Y ( i^{\prime} ) > Y ( i^{0} )\)，则令 \(i^{0} = i^{\prime}\)。

许多基于仿真的优化算法的一个难题是知道何时停止。（那些保证正确选择概率的算法除外。）典型的规则可能是：经过一定次数的迭代后停止；当最优解在若干次迭代中没有多大变化时停止；或者当可用于求解该问题的全部时间已经耗尽时停止。无论使用哪一种规则，我们都建议对表面最优解中排名前 5% 到 10% 的那些应用一个统计选择程序，例如第 2.2 节中的选择最佳方案程序。这样做是为了以有保证的置信度评价它们之中哪一个是真正的最佳者，并精确估计它们的性能。如果搜索过程中的原始数据已被保存，那么这些数据可以用作两阶段选择程序的第一阶段样本。

### 例 11：实施随机搜索

假设一个制造系统由 4 个串联的工作站组成。第 0 个工作站总有原材料可用。当第 0 个工作站完成对一个零件的加工，它就把该零件传给第 1 个工作站，然后第 1 个工作站把零件传给第 2 个工作站，如此下去。工作站 0 与 1、1 与 2、2 与 3 之间的缓冲空间总共限制为 50 个零件。例如，如果工作站 2 加工完一个零件，而工作站 3 前面没有可用的缓冲空间，那么工作站 2 就被阻塞（blocked），即它无法再做任何工作。问题是如何分配这 50 个位置，以使一个班次内每件的期望周期时间最小。

令 \(x_{i}\) 为工作站 \(i\) 前面的缓冲区位置数。于是决策变量为 \(x_{1} , x_{2} , x_{3}\)，其约束为 \(x_{1} + x_{2} + x_{3} = 50\)（分配的缓冲区位置数少于我们可用的数量是没有意义的）。这意味着总共有 1326 个可能的设计（你能算出这个数字是怎么来的吗？）。

为了简化随机搜索算法的描述，把解 \(( x_{1} , x_{2} , x_{3} )\) 的计数器记为 \(C ( x_{1} , x_{2} , x_{3} )\)。

### 随机搜索算法

步骤 1. 初始化 1326 个计数器变量 \(C ( x_{1} , x_{2} , x_{3} ) ~ = ~ 0\)，每一个对应一个可能的解 \(( x_{1} , x_{2} , x_{3} )\)。选择一个初始解，比如说 \(( x_{1} = 20 , x_{2} = 15 , x_{3} = 15 )\)，并令 \(C(20, 15, 15) = 1\)。

步骤 2. 从除 (20, 15, 15) 之外的所有解组成的集合中选取另一个解，使每一个解都有相等的被选中机会。假设选中了 (11, 35, 4)。

步骤 3. 在两个解处运行仿真实验，得到期望周期时间 Y(20, 15, 15) 与 Y(11, 35, 4) 的估计。假设 \(Y ( 20 , 15 , 15 ) ~ < ~ Y ( 11 , 35 , 4 )\)。那么 (20, 15, 15) 仍作为当前的好解。

步骤 4. 令 \(C ( 20 , 15 , 15 ) = C ( 20 , 15 , 15 ) + 1\)。

步骤 2. 从除 (20, 15, 15) 之外的所有解组成的集合中选取另一个解，使每一个解都有相等的被选中机会。假设选中了 (28, 12, 10)。

步骤 3. 在两个解处运行仿真实验，得到期望周期时间 Y(20, 15, 15) 与 Y(28, 12, 10) 的估计。假设 \(Y ( 28 , 12 , 10 ) \ : < \ : Y ( 20 , 15 , 15 )\)。那么 (28, 12, 10) 成为当前的好解。

步骤 4. 令 \(C ( 28 , 12 , 10 ) = C ( 28 , 12 , 10 ) + 1 ,\) 步骤 2. 从除 (28, 12, 10) 之外的所有解组成的集合中选取另一个解，使每一个解都有相等的被选中机会。假设选中了 (0, 14, 36)。

步骤 3. 继续……

当搜索终止时，我们选择使 \(C ( x_{1} , x_{2} , x_{3} )\) 计数最大的那个解 \(( x_{1} , x_{2} , x_{3} )\)。正如我们前面所讨论的，随后应当对前 5 到 10 个解做一次单独的统计分析，以高置信度确定它们之中哪一个是真正的最佳者。在本例中，计数最大的那些解将接受第二次分析。

尽管随机搜索算法看似简单，我们却略过了在具有可证明性能的算法中经常出现的一个微妙问题。在步骤 2 中，算法必须随机地选取一个解，使所有解（当前解除外）都有相等的可能性被选中。在例 11 中如何做到这一点？约束 \(x_{1} + x_{2} + x_{3} = 50\) 意味着 \(x_{1} , x_{2}\) 和 \(x_{3}\) 不能被独立地抽样。有人可能会想把 \(x_{1}\) 作为 0 到 50 上的离散均匀随机变量来抽样，然后把 \(x_{2}\) 作为 \(0 \mathrm{t o} 50 - x_{1}\) 上的离散均匀变量来抽样，最后令 \(x_{3} = 50 - x_{1} - x_{2}\)。但这种方法并不能使所有解等可能，下面的例子可以说明这一点：假设 \(x_{1}\) 被随机抽到 50。那么试验解必定是 \(( 50 , 0 , 0 )\)；只有一种选择。但如果 \(x_{1} = 49\)，那么 (49, 1, 0) 和 \(( 49 , 0 , 1 )\) 都是可能的。因此，如果要使所有满足 \(x_{1} + x_{2} + x_{3} = 50\) 的解等可能，那么 \(x_{1} = 49\) 就应当比 \(x_{1} = 50\) 更可能出现。

## 5 小结

本章对基于仿真运行所收集数据的备选系统设计比较评价作了基本介绍。这里假定已经选定一组固定的备选系统设计以供考虑。重点放在基于置信区间的比较以及公共随机数的使用上。还简要介绍了元模型——其目的是描述设计变量与输出响应之间的关系——以及基于仿真的优化——其目的是从数量庞大且种类繁杂的系统设计中选出最佳者。在仿真相关的统计分析技术这一领域，还有许多超出本书范围的、可能令人感兴趣的其他主题。其中一些主题是：

a. 实验设计模型（experimental design model），其目的是发现哪些因子对备选系统的性能具有显著影响；

b. 除重复法与批均值法之外的输出分析方法；

c. 方差缩减技术，即提高仿真实验统计效率的方法（公共随机数就是一个重要的例子）。

关于这些主题以及与仿真相关的其他主题，请读者参阅 Banks [1998] 和 Law [2007]。

最重要的思想是：仿真输出数据需要经过统计分析才能被正确解释。特别是，统计分析可以提供仿真所产生结果之精度的度量，并提供用以达到规定精度的技术。

## 参考文献

ANDRADOTTIR, S. [1998], “Simulation Optimization,” Chapter 9 in Handbook ofSimulation, J. Banks, ed., Wiley, New York.

ANDRADOTTIR, S. [1999], “Accelerating the Convergence ofRandom Search Methods for Discrete Stochastic Optimization,” ACM TOMACS, Vol. 9, pp. 349–380.

BANKS, J., ed. [1998], Handbook ofSimulation, Wiley, New York.

BATES, D. M., AND D. G. WATTS [1988], Nonlinear Regression Analysis and its Applications, Wiley, New York.

BECHHOFER, R. E., T. J. SANTNER, AND D. GOLDSMAN [1995], Design and Analysis for Statistical Selection, Screening and Multiple Comparisons, Wiley, New York.

BOESEL, J., B. L. NELSON, AND N. ISHII [2003], “A Framework for Simulation-Optimization Software,” IIE Transactions, Vol. 35, pp. 221–229.

BOESEL, J., B. L. NELSON,AND S. KIM [2003], “Using Ranking and Selection to ‘Clean Up’After Simulation Optimization,” Operations Research, Vol. 51, pp. 814–825.

BOX, G. E. P., AND N. R. DRAPER [1987], Empirical Model-Building and Response Surfaces, Wiley, New York.

FU, M. C. [2002], “Optimization for Simulation: Theory vs. Practice,” INFORMS Journal on Computing, Vol. 14, pp. 192–215.

GLOVER, F. [1989], “Tabu Search—Part I,” ORSA Journal on Computing, Vol. 1, pp. 190–206.

GOLDSMAN, D., AND B. L. NELSON [1998], “Comparing Systems via Simulation,” Chapter 8 in Handbook ofSimulation, J. Banks, ed., Wiley, New York.

HINES, W. W., D. C. MONTGOMERY, D. M. GOLDSMAN, AND C. M. BORROR [2002], Probability and Statistics in Engineering, 4th ed., Wiley, New York.

HOCHBERG, Y., AND A. C. TAMHANE [1987], Multiple Comparison Procedures, Wiley, New York.

HSU, J. C. [1996], Multiple Comparisons: Theory and Methods, Chapman & Hall, New York.

KIM, S., AND B. L. NELSON [2006], “Selecting the Best System,” Chapter 17 in Handbooks in Operations Research and Management Science: Simulation, S. G. Henderson and B. L. Nelson, eds., North-Holland, New York.

KLEIJNEN, J. P. C. [1975], Statistical Techniques in Simulation, Parts I and II, Dekker, New York

KLEIJNEN, J. P. C. [1987], Statistical Toolsfor Simulation Practitioners, Dekker, New York.

KLEIJNEN, J. P. C. [1988], “Analyzing Simulation Experiments with Common Random Numbers,” Management Science, Vol. 34, pp. 65–74.

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

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

MONTGOMERY, D. C. [2000], Design and Analysis ofExperiments, 5th ed., Wiley, New York.

NELSON, B. L., J. SWANN, D. GOLDSMAN, AND W.-M. T. SONG [2001], “Simple Procedures for Selecting the Best System when the Number of Alternatives is Large,” Operations Research, Vol. 49, pp. 950–963.

NELSON, B. L. [1987], “Some Properties of Simulation Interval Estimators Under Dependence Induction,” Operations Research Letters, Vol. 6, pp. 169–176.

NELSON, B. L. [1992], “Designing Efficient Simulation Experiments,” Proceedings ofthe 1992 Winter Simu lation Conference, J. J. Swain, D. Goldsman, R. C. Crain, and J. R. Wilson, eds., Arlington, VA, Dec. 13–16, pp. 126–132.

OLAFSSON, S. [2006], “Metaheuristics,” Chapter 21 in Handbooks in Operations Research and Management Science: Simulation, S. G. Henderson and B. L. Nelson, eds., North-Holland, New York.

WHITT, W. [1989], “Planning Queueing Simulations,” Management Science, Vol. 35, pp. 1341–1366.

WRIGHT, R. D., AND T. E. RAMSAY, JR. [1979], “On the Effectiveness of Common Random Numbers,” Management Science, Vol. 25, pp. 649–656.

YANG, F., B. ANKENMAN, AND B. L. NELSON [2007], “Efficient Generation of Cycle Time-Throughput Curves through Simulation and Metamodeling,” Naval Research Logistics, Vol. 54, pp. 78–93.

## 习题

1. 某公司拥有自卸卡车，这些卡车反复经历三项活动：装载、称重和运输。假设共有 8 辆卡车，且在时刻 0 时 8 辆卡车都在装载机处。每台卡车在唯一一台地磅上的称重时间在 1 到 9 分钟之间均匀分布，每台卡车的运输时间服从均值为 85 分钟的指数分布。装载机前和地磅前都允许无限排队。所有卡车可以同时处于运输状态。管理层希望把一台快速装载机与目前正在使用的两台较慢装载机作比较。每台慢速装载机装满一辆卡车的时间在 1 到 27 分钟之间均匀分布。新的快速装载机装满一辆卡车的时间在 1 到 19 分钟之间均匀分布。比较的依据是 40 小时时间跨度上的平均系统响应时间，其中响应时间定义为从一辆卡车到达装载机队列到该卡车驶离地磅的持续时间。请用公共随机数对这两个选项作一次统计上有效的比较。

2. 在《绝对性能估计》一章的习题 6 中，考虑如下备选的 (M, L) 策略：

|  | L |  |  |  |
|---|---|---|---|---|
| 低 30 | 高 40 |  |  |  |
| M | 低 | 50 | (50,30) | (50, 40) |
| 高 | 100 | (100,30) | (100, 40) |  |

请利用《绝对性能估计》一章习题 6 中所建立的仿真模型作适当修改，考察这些策略的相对成本。以长期平均月成本为依据比较这四个系统设计。首先，对每个 (M, L) 策略做 4 次重复，尽可能充分地采用公共随机数。每次重复应当包含 12 个月的初始化阶段，随后是 100 个月的数据收集阶段。计算每个策略平均月成本的置信区间，使其总体置信水平为 90%。然后估计为了使各置信区间不再重叠所需要的追加重复次数。就哪一个策略最佳得出结论。

3. 重新考虑《绝对性能估计》一章的习题 7。比较习题 2 中所研究的四种库存策略，并在计算月成本时把紧急订货的成本考虑在内。

4. 重新考虑《绝对性能估计》一章的习题 8。比较习题 2 中所研究的四种月度库存策略，并把易腐物品的销售价格考虑在内。

5. 在《绝对性能估计》一章的习题 9 中，考察订货量对长期平均日成本的影响。每笔订货用送货卡车上的一个托盘运到，因此允许的订货量 Q 是 10 的倍数；即 Q 可以等于 10、20、30……。在习题 9 中考察的是 Q = 20 这一策略。

(a) 首先，考察 Q = 10 与 Q = 50 这两个策略。使用《绝对性能估计》一章习题 9 中所建议的运行长度等设定。在这些运行的基础上，判断最优的 Q（记为 Q∗）是位于 10 与 50 之间，还是大于 50。（成本曲线作为 Q 的函数应当具有什么样的形状？）

(b) 利用 (a) 中的结果，再提出两个 Q 的取值并对这两个策略作仿真。得出结论，并分析你所得结论的强度。

6. 在《绝对性能估计》一章的习题 10 中，求出卡片店店主应当采购的卡片数量 Q，以使利润最大化，且（总利润的）误差至多约为 5.00 美元。首先，对订购 Q = 250 张和 Q = 350 张的策略各做一次运行。利用这些结果，再加上习题 10 中订购 \(Q = 300\) 张的结果，确定一个值得进一步考虑的 Q 的范围：\(200 \leq Q \leq 250\)，或 \(250 \leq Q \leq 300 .\)，等等。然后在这个受限范围内再考察两个策略：Q = 200 与 \(Q = 225 ,\)，或 Q = 265 与 Q = 285，等等。

7. 在《绝对性能估计》一章的习题 11 中，考察目标水平 M 与检查周期 N 对平均月成本的影响。考虑两个目标水平 M，由习题 11 中所用的目标水平加减 10 确定，并考虑 1 个月和 3 个月的检查周期 N。根据这些仿真，哪一个 (N, M) 组合最佳？

8. 重新考虑《绝对性能估计》一章的习题 13 和习题 14，它们涉及一个作业车间中的调度规则（即排队规则）：先进先出（FIFO）与按类型优先（PR）。除了这两条规则，再考虑一种最短即时工序（SIO）调度规则：对某一个给定的工作站，具有最小平均加工时间的那一类型的所有作业都被赋予最高优先级。例如，当在工作站 1 使用 SIO 规则时，作业按如下顺序加工：先类型 II，再类型 I，最后类型 III。同一类型的两个作业按 FIFO 处理。请设计一个仿真实验，以所有作业的平均总响应时间为依据比较 FIFO、PR 和 SIO 三条规则。

9. 在《绝对性能估计》一章的习题 13（采用 FIFO 规则的作业车间）中，求出每个工作站为避免瓶颈所需要的最少工人数。当某个工作站在一段时间上的平均队列长度持续增加时，就出现了瓶颈。［提示：不要把由于服务台数量不足而引起的平均队列长度增加，与由于初始化偏差而引起的平均队列长度增加相混淆。在前一种情形下，平均队列长度会无限持续增加，且服务台利用率为 1.0。在后一种情形下，平均队列长度最终会趋于平稳，且服务台利用率小于 1。］请按作业类型以及全部类型报告工人的利用率和一个作业通过整个作业车间所需的总时间。

［提示：如果某个工作站的服务台利用率为 1.0，且平均队列长度随仿真运行长度的增加而趋于线性增长，那么该工作站很可能是不稳定的，因而是一个瓶颈。在这种情况下，该工作站至少需要增加一名工人。利用排队论，即 \(\lambda / c_{i} \mu < 1\)，给出工作站 1 所需最少工人数的建议。回顾一下，λ 是到达率，\(1 / \mu\) 是一名工人加工一个作业的总体平均服务时间，\(c_{i}\) 是工作站 \(i\) 的工人数。尝试用同一个基本条件 \(\lambda / c_{i} \mu < 1\) 给出工作站 \(i\)（\(i = 2 , 3 , 4\)）初始服务台数量的建议。］

10. (a) 对 PR 调度规则重做习题 9（参见《绝对性能估计》一章的习题 14）。

(b) 对 SIO 调度规则重做习题 9（参见习题 8）。

(c) 比较每种调度规则所需的最少工人数：FIFO 对 PR 对 SIO。

11. 在上面习题 9 与习题 10 中为《绝对性能估计》一章习题 13 的作业车间所求得的最少工人数的基础上，考虑给整个车间增加一名工人。这名工人只能被培训去处理某一个工作站的加工。这名工人应当安排在哪个工作站？这名额外工人如何影响所有作业的平均总响应时间？如何影响类型 I 作业的平均总响应时间？对每一种调度规则（FIFO、PR、SIO），分别考察有这名额外工人与没有这名额外工人的作业车间。

12. 在《绝对性能估计》一章的习题 17 中，假设在每个工人前面建造一个容量为一件物品的缓冲区。设计一个实验，考察系统设计上的这一改变是否对各工人的利用率 \(\rho_{1} , \rho_{2} , \rho_{3}\) 和 \(\rho_{4}\) 具有显著影响。至少要计算 \(\rho_{1} ^ {0} - \rho_{1} ^ {1}\) 和 \(\rho_{4} ^ {0} - \rho_{4} ^ {1}\) 的置信区间，其中 \(\rho_{i} ^ {s}\) 是缓冲区容量为 \(s\) 时工人 \(i\) 的利用率。

13. 某小型州立大学招生办公室的一名职员处理索取招生材料的请求。处理请求的时间取决于所感兴趣的专业（例如工业工程、管理科学、计算机科学等）以及该专业的层次（学士、硕士、博士）。假设处理时间可以很好地用正态分布建模，均值为 7 分钟、标准差为 2 分钟。每天开始工作时，这名职员需要一些时间才能准备就绪、开始处理请求；假设这一时间可以很好地用指数分布建模，均值为 20 分钟。招生办公室每天通常收到 40 到 60 份请求。

令 x 为一天中收到的申请数量，令 Y 为处理它们所需的时间（包括准备时间）。通过在设计点 \(x = 40 , 50\), 60 上各做 \(n\) 次重复，为 $\operatorname{E} ( Y | x ) $ 拟合一个元模型。注意，在这种情形下，我们知道正确的模型是 

$$
\operatorname{E} (Y | x) = \beta_{0} + \beta_{1} x = 20 + 7 x
$$

（为什么？）先从每个设计点 \(n = 2\) 次重复开始，估计 \(\beta_{0}\) 和 \(\beta_{1}\)。逐步增加重复次数，观察需要多少次重复才能使估计接近真实值。

14. 用 CRN 重做上一题。结果有何变化？

15. 考虑简单线性回归模型 \(Y = \beta_{0} + \beta_{1} x + \varepsilon\)。找一本涵盖回归建模的统计学教材。在那里你会找到关于假设 \(\beta_{1} = 0\) 的一个检验（拒绝这一假设就是 x 与 Y 线性相关的证据）。如果我们使用 CRN，则该检验背后的通常假定并不成立。它们在哪里失效了？

16. Riches and Associates 公司主要以定期存单（CD）的形式持有现金储备，年利率为 8%。然而，公司必须定期从这些存单中提取资金，以支付供应商款项等。这些现金流出通过一个不计息的支票账户进行。对现金的需求无法被确定地预测。从存单到支票账户的转账可以即时完成，但从存单提前支取会受到"重罚"。因此，R&A 公司使用其支票账户的透支保护也许是合理的，该透支按每美元每天 0.00033 美元（即年利率 12%）计收透支利息。

R&A 公司喜欢简单的策略：每年转账固定次数、每次固定金额。目前，它每年转账 6 次，每次 18,250 美元。你的任务是找到一种能够降低其长期日均成本的策略。

从历史规律判断，现金需求的到达率约为每天 1 次，其到达可以很好地用泊松过程建模。满足每一次需求所需的现金量可以用一个均值为 300 美元、标准差为 150 美元的对数正态随机变量合理地表示。

不同存单的提前支取罚金不同。每次支取平均为 150 美元，与金额无关，但实际罚金可以用一个范围为 100 美元到 200 美元的均匀分布随机变量来建模。

用支票账户中的现金水平来确定初始化阶段的长度。做足够多的重复，使你关于长期日均成本之差的置信区间不包含零。务必在实验设计中使用 CRN。

17. 如果你能使用商业化的基于仿真的优化软件，请测试随着仿真输出变异性的增大，它的效果如何。用一个简单的模型，例如 \(Y = x^{2} + \varepsilon ,\) 其中 ε 是服从 \(\mathrm{N} ( 0 , \sigma^{2} )\) 分布的随机变量，且最优解已知（在本例中，最小化问题的解为 \(x = 0\)）。看看随着 \(\sigma^{2}\) 的增大，该软件能否（以及多快能够）找到真正的最优解。接着尝试具有多于一个设计变量的更复杂的模型。

18. 对于例 11，说明为什么共有 1326 个解。然后导出一种对 \(x_{1} , x_{2}\) 和 \(x_{3}\) 抽样的方法，使得 \(x_{1} + x_{2} + x_{3} = 50\) 且所有结果等可能。

19. 某种关键电子元件，其平均失效时间为 x 年，售价为 2x 千美元（因此元件越可靠就越贵）。x 的取值限制在 1 到 10 年之间，实际失效时间用指数分布建模。使用该元件执行的任务持续一年；如果该元件在不到一年的时间内失效，则会产生 20,000 美元的早期失效成本。应当选择什么样的 x 值，以使期望总成本（购置成本加早期失效成本）最小？

为了求解这一问题，建立一个仿真，对平均失效时间为 x 年的元件生成总成本。这需要抽样一个均值为 x 的指数分布随机变量，然后把总成本计算为 2000x 加上（若失效时间小于 1 年）20,000。拟合一个关于 x 的二次元模型，并用它求出使所拟合模型最小的 x 值。［提示：在 1 到 10 之间选取若干个 x 值作为设计点。在每个 x 值处，令响应变量 Y(x) 为至少 30 次总成本观测值的平均。］

20. 用基于仿真的优化软件求解习题 19。如果你无法使用这类软件，就使用随机搜索算法，并令 x 的可能取值为 1.00, 1.25, 1.50, . . . , 10.00。

21. 假设某种产品的需求服从均值为 10 件的泊松分布。用基于仿真的优化软件（如果无法使用这类软件，就用随机搜索算法）求出使需求等于 x 件的概率最大的订货量 x，假定 \(0 \leq x \leq 100\) 且 x 为整数。（本题基于 Andradóttir [1999]。）［提示：对某个试验值 x 的仿真，每一次重复生成一个均值为 10 的泊松随机需求。如果该随机需求等于 x，则响应为 Y(x) = 1；否则为 Y(x) = 0。］

22. 在习题 16 中，考虑以下 5 种策略：6 次转账、每次 18,250 美元；5 次转账、每次 21,900 美元；4 次转账、每次 27,375 美元；3 次转账、每次 36,500 美元；2 次转账、每次 54,750 美元。对每种策略做 10 次重复，并为所有策略对的期望成本之差构造同时置信区间。

23. 在习题 16 中，考虑以下 5 种策略：6 次转账、每次 18,250 美元；5 次转账、每次 21,900 美元；4 次转账、每次 27,375 美元；3 次转账、每次 36,500 美元；2 次转账、每次 54,750 美元。对每种策略做 10 次重复，然后应用选择最佳方案程序选出期望日均成本最低的策略。取 \(\epsilon = \$2\) 每天。

## www.bcnn.net 上可用材料清单

SimulationTools.xls