12 Estimation of Relative Performance

This chapter discusses a few of the many statistical methods that can be used to compare the relative performance of two or more system designs. This is one of the most important uses of simulation. Because the observations of the response variables contain random variation, statistical analysis is needed to discover whether any observed differences are due to differences in design or merely to the random fluctuation inherent in the models.

The comparison of two system designs is easier than the simultaneous comparison of multiple system designs. Section 1 discusses the case of two system designs, using two possible statistica techniques: independent sampling and correlated sampling. Correlated sampling is also known as the common random numbers (CRN) technique; simply put, the same random numbers are used to simulate both alternative system designs. If implemented correctly, CRN usually reduces the variance of the estimated difference of the performance measures and thus can provide, for a given sample size, more precise estimates of the mean difference than can independent sampling. Section 2 extends the statistical techniques of Section 1 to the comparison of multiple (more than two) system designs using simultaneous confidence intervals, screening and selection of the best. These approaches are somewhat limited with respect to the number of system designs that can be considered, so Section 3 describes how a large number of complex system designs can sometimes be represented by a simpler

metamodel. Finally, for comparison and evaluation of a very large number of system designs tha are related in a less structured way, Section 4 presents optimization via simulation.

1 Comparison of Two System Designs

Suppose that a simulation analyst desires to compare two possible configurations of a system. In a queueing system, perhaps two possible queue disciplines, or two possible sets of servers, are to be compared. In a supply-chain inventory system, perhaps two possible ordering policies will be compared. A job shop could have two possible scheduling rules; a production system could have in-process inventory buffers of various capacities. Many other examples of alternative system designs can be provided.

The method of replications will be used to analyze the output data. The mean performance measure for system i will be denoted by $\theta_{i} , i = 1 , 2 .$ If it is a steady-state simulation, it will be assumed that deletion of data, or other appropriate techniques, have been used to ensure that the point estimators are approximately unbiased estimators of the mean performance measures $\theta_{i} .$. The goal of the simulation experiments is to obtain point and interval estimates of the difference in mean performance, namely $\theta_{1} - \theta_{2}$. Two methods of computing a confidence interval for $\theta_{1} - \theta_{2}$ will be discussed, but first an example and a general framework will be given.

Example 1

Software Made Personal (SMP) has a customer support call center where 7 operators handle questions from owners of their software from 8 A.M. to 4 P.M. Eastern Time. When a customer calls they use an automated system to select between the two product lines, finance and contact management. Currently each product line has its own operators (4 for financial, 3 for contact management) and hold queue. SMP wonders if they can reduce the total number of operators needed by cross-training them so that they can answer calls for any product line. This is expected to increase the time to process a call by about 10%. The current system is illustrated in Figure 1(a). and the alternative system design is shown in Figure 1(b). Before considering reducing operators, SMP wants to know what their quality of service would be if they simply cross-train the 7 operators they currently have. Incoming calls have been modeled as a nonstationary Poisson arrival process and distributions have been fit to operator time to answer questions. Comparisons will be based on the average time from call initiation until the customer’s questions are answered.

(a) Figure 1 Current SMP call center (a) and an alternative design (b).

Table 1 Simulation Output Data and Summary Measures for Comparing Two Systems

When comparing two systems, such as those in Example 1, the simulation analyst must decide on a run length $T_{E} ^ {( i )}$ for each model $i = 1 , 2 .$, a number of replications $R_{i}$ to be made of each model, or both. From replication r of system i, the simulation analyst obtains an estimate $Y_{r i}$ of the mean performance measure $\theta_{i}$. In Example 1,$Y_{r i}$ would be the average call response time observed during replication r for system i, where $r = 1 , 2 , \ldots , R_{i} ;$ and $i = 1 , 2$. The data, together with the two summary measures, the sample means ${\bar {Y}} _ {\cdot i} ,$, and the sample variances $S_{i} ^ {2}$, are exhibited in Table 1. Assuming that the estimators $Y_{r i}$ are, at least approximately, unbiased, it follows that

In Example 1, SMP is initially interested in a comparison of two system designs (same number of operators, but specialized vs. cross trained), so the simulation analyst decides to compute a confidence intervalfor $\theta_{1} - \theta_{2}$, the difference between the two mean performance measures. This will lead to one of three possible conclusions:

1. If the confidence interval (c.i.) for $\theta_{1} - \theta_{2}$ is totally to the left of zero, as shown in Figure 2(a), then there is strong evidence for the hypothesis that $\theta_{1} - \theta_{2} < 0$, or equivalently $\theta_{1} < \theta_{2}$ In Example 1 $, \theta_{1} < \theta_{2}$ implies that the mean response time for system 1, the original system, is smaller than for system 2, the alternative system.

2. If the c.i. for $\theta_{1} - \theta_{2}$ is totally to the right of zero, as shown in Figure 2(b), then there is strong evidence that $\theta_{1} - \theta_{2} > 0 .$, or equivalently,$\theta_{1} > \theta_{2}$ In Example 1,$\theta_{1} > \theta_{2}$ can be interpreted as system 2 being better than system 1, in the sense tha system 2 has a smaller mean response time.

3. If the c.i. for $\theta_{1} - \theta_{2}$ contains zero, then, with the data at hand, there is no strong statistica evidence that one system design is better than the other.

Figure 2 Three confidence intervals that can occur in the comparing of two systems.

Some statistics textbooks say that the weak conclusion $\theta_{1} = \theta_{2}$ can be drawn, but such statements can be misleading. A “weak” conclusion is often no conclusion at all. Most likely, if enough additional data were collected $( \mathrm{i} . \mathrm{e} . , R_{i}$ increased), the c.i. would shift, and definitely shrink in length, unti conclusion 1 or 2 would be drawn. Since the confidence interval provides a measure of the precision of the estimator of $\theta_{1} - \theta_{2}$, we want to shrink it until we can either conclude that there is a difference that matters, or the difference is so small that we do not need to separate them.

In this chapter, a two-sided 100(1 α)% c.i. for $\theta_{1} - \theta_{2}$ will always be of the form

where ${\bar {Y}} _ {\cdot i}$ is the sample mean performance measure for system i over all replications

ν is the degrees offreedom associated with the variance estimator,$t_{\alpha / 2 , 1}$ is the $100 ( 1 {-} \alpha / 2 )$ percentage point of a t distribution with ν degrees of freedom, and $\mathbf{s . e .} ( \cdot )$ represents the standard error of the specified point estimator. To obtain the standard error and the degrees of freedom, the analyst uses one of two statistical techniques. Both techniques assume that the basic data,$Y_{r i}$ of Table 1, are approximately normally distributed. This assumption is reasonable provided that each $Y_{r i}$ is itself a sample mean of observations from replication $r ,$ which is indeed the situation in Example 1.

By design of the simulation experiment,$Y_{r 1} ( r = 1 , \dots , R_{1} )$ are independently and identically distributed (i.i.d.) with mean $\theta_{1}$ and variance $\sigma_{1} ^ {2} ,$, say. Similarly,$Y_{r 2} ( r = 1 , \ldots , R_{2} )$ are i.i.d. with mean $\theta_{2}$ and variance $\sigma_{2} ^ {2} .$, say. The two techniques for computing the confidence interval in Equation (1), which are based on different sets of assumptions, are discussed in the following subsections.

There is an important distinction between statistically significant differences and practically significant differences in systems performance. Statistical significance answers the following question: Is the observed difference ${\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2}$ larger than the variability in ${\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2} ?$ This question can be restated as: Have we collected enough data to be confident that the difference we observed is real, or just chance? Conclusions 1 and 2 above imply a statistically significant difference, while conclusion 3 implies that the observed difference is not statistically significant, even though the systems may indeed be different. Statistical significance is a function of the simulation experiment and the output data.

Practical significance answers the following question: Is the true difference $\theta_{1} - \theta_{2}$ large enough to matter for the decision we need to make? In Example 1, we may reach the conclusion that $\theta_{1} > \theta_{2}$ and decide that system 2 is better (smaller expected response time). However, if the actual difference $\theta_{1} - \theta_{2}$ is very small—say, small enough that a customer would not notice the improvement—then it might not be worth the cost to replace system 1 with system 2. Practical significance is a function of the actual difference between the systems and is independent of the simulation experiment.

Confidence intervals do not answer the question of practical significance directly. Instead, they bound, with probability $1 - \alpha$, the true difference $\theta_{1} - \theta_{2}$ within the range

Whether a difference within these bounds is practically significant depends on the particular problem.

1.1 Independent Sampling

Independent sampling means that different and independent random number streams will be used to simulate the two systems. This implies that all the observations of simulated system 1, namely $\{Y_{r 1} , r = 1 , \ldots , R_{1} \}$, are statistically independent of all the observations of simulated system 2, namely $\{Y_{r 2} , r = 1 , \dots , R_{2} \}$. By Equation (2) and the independence of the replications, the variance of the sample mean ${\bar {Y}} _ {\cdot i}$ is given by

For independent sampling,${\bar {Y}} _ {\cdot 1}$ and ${\bar {Y}} _ {\cdot 2}$ are statistically independent; hence,

There are confidence-interval procedures that are appropriate when the two variances are equa but unknown in value; that is,$\sigma_{1} ^ {2} = \sigma_{2} ^ {2}$. However, since we would rarely know this to be true in practice, and the advantages of exploiting equal variances are slight provided the number of replications is not too small, we refer the reader to any statistics textbook for this case (for instance Hines et al. [2002]).

An approximate $100 ( 1 - \alpha ) \%$ c.i. for $\theta_{1} - \theta_{2}$ can be computed as follows. The point estimate is computed by

with ${\bar {Y}} _ {\cdot i}$ given by Equation (2), while the sample variance for system i is

Thus, the standard error of the point estimate is given by

with degrees of freedom ν approximated by the expression

rounded to an integer. The confidence interval is then given by Expression $( 1 ) ,$, using the standard error of Equation (6). A minimum number of replications $R_{1} \geq 6$ and $R_{2} \geq 6$ is recommended for this procedure.

1.2 Common Random Numbers (CRN)

CRN (also known as correlated sampling) means that, for each replication, the same random numbers are used to simulate both systems. Therefore,$R_{1}$ and $R_{2}$ must be equal, say,$R_{1} = R_{2} = R$. Thus, for each replication $r ,$ the two estimates $Y_{r 1}$ and $Y_{r 2}$ are no longer independent, but rather are correlated. However, independent random numbers are used on different replications, so the pairs $( Y_{r 1} , Y_{s 2} )$ are mutually independent when $r \neq s .$. For example, in Table 1, the observation $Y_{11}$ is correlated with $Y_{12}$, but $Y_{11}$ is independent of all other observations. The purpose of using CRN is to induce a positive correlation between $Y_{r 1}$ and $Y_{r 2}$ for each r and thus to achieve a variance reduction in the point estimator of mean difference,${\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2}$. In general, this variance is given by

where $\rho_{12}$ is the correlation between $Y_{r 1}$ and $Y_{r 2}$. [By definition,$\rho_{12} = \operatorname{c o v} ( Y_{r 1} , Y_{r 2} ) / \sigma_{1} \sigma_{2}$, which does not depend on r.]

Now, compare the variance of ${\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2}$ arising from the use of CRN [Equation $( 8 ) _ {:}$, call it $V_{C R N} ]$ to the variance arising from the use of independent sampling with equal sample sizes [Equation (3) with $R_{1} = R_{2} = R$, call it $V_{I N D} ]$. Notice that

If CRN works as intended, the correlation $\rho_{12}$ will be positive; hence, the second term on the righ side of Equation (9) will be positive, and, therefore,

That is, the variance of the point estimator will be smaller with CRN than with independent sampling. A smaller variance, for the same sample size $R ,$ implies that the estimator based on CRN is more precise, leading to a shorter confidence interval on the difference, which implies that smaller differences in performance can be detected.

To compute a $100 ( 1 - \alpha ) \%$ c.i. with correlated data, first compute the differences

which, by the definition of CRN, are i.i.d.; then compute the sample mean difference as

Thus,$\bar {D} = \bar {Y} _ {\cdot 1} - \bar {Y} _ {\cdot 2}$. The sample variance of the differences $\{D_{r} \}$ is computed as

which has degrees of freedom $\nu = R - 1$. The $100 ( 1 - \alpha ) \% \mathrm{c . i}$. for $\theta_{1} - \theta_{2}$ is given by Expression (1), with the standard error of ${\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2} = {\bar {D}}$, estimated by

Because $S_{D} / \sqrt {R}$ of Equation (13) is an estimate of $\sqrt {V_{C R N}}$ and Expression (6) is an estimate of $\sqrt {V_{I N D}}$, CRN typically will produce a c.i. that is shorter for a given sample size than the c.i. produced by independent sampling if $\rho_{12} > 0 .$. In fact, the expected length of the c.i. will be shorter with the use of CRN if $\rho_{12} > 0 . 1$, provided $R > 10$. The larger R is, the smaller $\rho_{12}$ can be and still yield a shorter expected length [Nelson 1987].

For any problem, there are many ways of implementing common random numbers. It is never enough to simply use the same seed on the random-number generator(s). Each random number used in one model for some purpose should be used for the same purpose in the second model—that is, the use of the random numbers must be synchronized. For example, if the ith random number is used to generate the call service time of an operator for the 5th caller in model 1, then the ith random number should be used for the very same purpose in model 2. For queueing systems or service facilities, synchronization of the common random numbers guarantees that the two systems face identica work loads: both systems face arrivals at the same instants of time, and these arrivals demand equa amounts of service. (The actual service times of a given arrival in the two models may not be equal; they could be proportional if the server in one model were faster than the server in the other model.) For an inventory system, in comparing different ordering policies, synchronization guarantees tha the two systems face identical demand for a given product. For production or reliability systems, synchronization guarantees that downtimes for a given machine will occur at exactly the same times and will have identical durations, in the two models. On the other hand, if some aspect of one o the systems is totally different from the other system, synchronization could be inappropriate—or even impossible to achieve. In summary, those aspects of the two system designs that are sufficiently similar should be simulated with common random numbers in such a way that the two models behave similarly; but those aspects that are totally different should be simulated with independent random numbers.

Implementation of common random numbers is model-dependent, but certain guidelines can be given that will make CRN more likely to yield a positive correlation. The purpose of the guidelines is to ensure that synchronization occurs:

1. Dedicate a random-number stream to a specific purpose, and use as many different streams as needed. (Different random-number generators, or widely spaced seeds on the same generator, can be used to get two different, nonoverlapping streams.) In addition, assign independently chosen seeds to each stream at the beginning of each replication. It is not sufficient to assign seeds at the beginning of the first replication and then let the random-number generator merely continue for the second and subsequent replications. If simulation is conducted in this manner, the first replication will be synchronized, but subsequent replications might no be.

2. For systems (or subsystems) with external arrivals: As each entity enters the system, the next interarrival time is generated, and then immediately all random variables (such as service times, order sizes, etc.) needed by the arriving entity and identical in both models are generated in a fixed order and stored as attributes of the entity, to be used later as needed. Apply guideline 1: Dedicate one random-number stream to these external arrivals and all their attributes.

3. For systems having an entity performing given activities in a cyclic or repeating fashion, assign a random-number stream to this entity. (Example: a machine that cycles between two states: up–down–up–down–. . . . Use a dedicated random-number stream to generate the uptimes and downtimes.)

4. If synchronization is not possible, or if it is inappropriate for some part of the two models, use independent streams of random numbers for this subset of random variates.

Unfortunately, there is no guarantee that CRN will always induce a positive correlation between comparable runs of the two models. It is known that if, for each input random variate X, the estimators $Y_{r 1}$ and $Y_{r 2}$ are increasing functions of the random variate X (or both are decreasing functions o X), then $\rho_{12}$ will be positive. The intuitive idea is that both models, that is, both $Y_{r 1}$ and $Y_{r 2} .$ respond in the same direction to each input random variate, and this results in positive correlation. This increasing or decreasing nature of the response variables, called monotonicity, with respec to the input random variables, is known to hold for certain queueing systems, such as the $G I / G / c$ queues, when the response variable is customer delay, so some evidence exists that common random numbers is a worthwhile technique for queueing simulations. (For simple queues, customer delay is an increasing function of service times and a decreasing function of interarrival times.) Wrigh and Ramsay [1979] reported a negative correlation for certain inventory simulations, however. In summary, the guidelines recently described should be followed, and some reasonable notion that the response variable of interest is a monotonic function of the random input variables should be evident.

Example 1: (cont.)

SMP’s call center designs in Figure 1 will be compared by simulating them with CRN. To illustrate several issues we will look at two comparisons: the current system with 7 operators vs. cross-trained system with 7 operators; and a cross-trained system with 7 operators vs. cross-trained system with 6 operators. We will also assess how much benefit was obtained by using CRN relative to independent sampling

There are four driving input processes in the SMP problem: The caller interarrival times, the caller type, the call times for financial software callers, and the call times for the contact managemen software callers. By “call time” we mean the time that an operator spends with a caller, excluding the time the caller spends on hold.

To implement CRN, we dedicate a distinct random number stream to each input process. This guarantees that precisely the same sequence of calls and caller types will be seen by each alternative on each replication. The call times could be generated in one of two ways, however: Either generate them as needed (that is, when a caller reaches an operator, generate the call time from the appropriate distribution), or generate the call time at the moment that the caller arrives, save it as an attribute of the caller, and use it when the caller reaches an operator. The latter approach (generate the call time at arrival and save it) provides better synchronization in this case; in fact, if we generate call times for the cross-trained alternative by multiplying call times from the current system by 1.1, then this approach guarantees that each caller in the cross-trained alternatives has call time that is precisely 10% greater than they would have experienced in the current design.

Table 2 shows the first 10 replications from a run of 300 replications of the three system designs (the complete data are available from www.bcnn.net). Figure 3 plots the average response times for the current system and the cross-trained system with 7 operators for these same 10 replications. The desired impact of CRN is apparent in the figure: The response averages from the two system designs tend to go up and down together, making their differences easier to observe. The sample correlation between them, based on the entire 300 replications, is 0.24, showing a modest impac of CRN. The sample correlation, again based on all 300 replications, between the proposed systems with 7 and 6 operators is a substantial 0.94, however. As shown in Figure 4, the responses of these two system designs track distinctly up and down together, revealing a strong positive correlation.

Figure 3 Average response times from the first 10 replications of the current and proposed (7 operators) SMP call center

A more quantitative way to see the impact of CRN is to look at the standard errors, because they directly impact the length of the confidence interval. For instance, the standard error for the comparison of the two proposed systems is 0.88 minutes. We can approximate what the standard error would have been with independent sampling by using the sample variances and Equation (6):

The reduction is almost 50%, making it far easier to detect differences.

Table 2 Comparison of System Designs for the SMP Call Center

Figure 4 Average response times from the first 10 replications of the proposed designs with 7 and 6 operators for the SMP call center

Understanding why the strength of the induced correlation due to CRN is different in these two cases helps in understanding when CRN will be most effective. The two proposed designs (6 and 7 operators) are structurally the same system, differing only in the number of operators. Therefore, when the random inputs are such that the system with 7 operators is congested, the system with 6 operators is certain to be at least as congested, and probably more so. The current system is structurally different from the alternatives, having a distinct queue for each type of caller. When congestion occurs in the financial software queue it may not occur in the contact management queue, so long response times for the financial callers could be averaged with shorter response times for the contact management callers. Thus, the relationship with the proposed system that has a single queue is not as strong.

Let the subscript C denote the current call center, P7 the proposal with 7 operators, and P6 the proposal with 6 operators. Then a 95% confidence interval for the difference $\theta_{C} - \theta_{P 7}$ is

$$

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

or

$$

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

Therefore, from these 10 replications we cannot say for certain which design has a lower response time, although we can say with high confidence that the difference is within 1.4 minutes.

On the other hand, and as we might expect, the 7 and 6 operator systems are distinctly different, with 95% confidence interval

$$

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

showing that having a 7th operator is at least 3 minutes faster, on average, than only having 6.

1.3 Confidence Intervals with Specified Precision

Suppose that we want the error in our estimate of $\theta_{1} - \theta_{2}$ to be less than  (the quantity  migh be a practically significant difference). Therefore, our goal is to find a number of replications R such that

We begin by making $R_{0} \geq 10$ replications ofeach system, to obtain an initial estimate ofs.e.$( {\bar {Y}} _ {\cdot 1} - {\bar {Y}} _ {\cdot 2} )$ We then approximate the total number of replications $R \ \geq \ R_{0}$ needed to achieve the half-length criterion in Equation (15). Finally, we make an additional $R {-} R_{0}$ replications (or a fresh R replications) of each system, compute the confidence interval, and check that the half-length criterion has been attained.

Example 1: (cont.)

Recall that $R_{0} = 10$ replications and CRN yielded the 95% confidence interval for the difference in expected response time of the current and 7-operator cross-trained system in Equation (14); this interval can be rewritten $\mathrm{a s} - 0 . 86 \pm 1 . 40$ minutes. Although the current call center design appears to have the smaller expected response time, the difference is not statistically significant, since the confidence interval contains 0. Suppose that a difference larger than 0.25 minutes is considered to be practically significant. We therefore want to make enough replications to obtain $H \leq \epsilon = 0 . 25$ The confidence interval used in Example 1 was of the form $\bar {D} \pm t _ {\alpha / 2 , R_{0} - 1} S_{D} / \sqrt {R_{0}}$, with the specific values $\bar {D} = - 0 . 86 , R_{0} = 10 , t_{0 . 025 , 9} = 2 . 26$, and $S_{D} ^ {2} = 3 . 87$. To obtain the desired precision, we need to find R such that

We can approximate this by finding the smallest integer R satisfying $R \geq R_{0}$ and

Substituting $z_{0 . 025} = 1 . 96$ and $S_{D} ^ {2} = 3 . 87$, we obtain

implying that 238 replications are needed, 228 more than in the initial experiment. There are 300 replications available at www.bcnn.net and the reader is encouraged to check the result.

2 Comparison of Several System Designs

Suppose that a simulation analyst desires to compare K alternative system designs. The comparison will be made on the basis of some specified performance measure $\theta_{i}$ of system i, for $i = 1 , 2 , \dots , K$ Many different statistical procedures have been developed that can be used to analyze simulation data and draw statistically sound inferences concerning the parameters $\theta_{i}$. These procedures can be classified as being either fixed-sample-size procedures or sequential-sampling (or multistage) procedures.

In the first type, a predetermined sample size (i.e., run length and number of replications) is used to draw inferences via hypothesis tests or confidence intervals. Examples of fixed-sample-size procedures include the interval estimation of a mean performance measure and the interval estimation of the difference between mean performance measures of two systems [as by Expression (1) in Section 1]. Advantages of fixed-sample-size procedures include a known or easily estimated cos in terms of computer time before running the experiments. When computer time is limited, or when a pilot study is being conducted, a fixed-sample-size procedure might be appropriate. In some cases, clearly inferior system designs may be ruled out at this early stage. A major disadvantage is that a strong conclusion could be impossible. For example, the confidence interval could be too wide for practical use, since the width is an indication of the precision of the point estimator. A hypothesis test may lead to a failure to reject the null hypothesis, a weak conclusion in general, meaning that there is no strong evidence one way or the other about the truth or falsity of the nul hypothesis.

A sequential sampling scheme is one in which more and more data are collected until an estimator with a prespecified precision is achieved or until one of several alternatives is selected, with the probability of correct selection being larger than a prespecified value. A two-stage (or multistage) procedure is one in which an initial sample is used to estimate how many additional observations are needed to draw conclusions with a specified precision. An example of a two-stage procedure for estimating the performance measure of two systems was given Section 1.3.

The proper procedure to use depends on the goal of the simulation analyst. Some possible goals are the following:

1. estimation of each parameter,$\theta_{i} ;$2. comparison of each performance measure $\theta_{i}$ to a control $\theta_{1}$(where $\theta_{1}$ could represent the mean performance of an existing system);

2. all pairwise comparisons,$\theta_{i} - \theta_{j}$, for $i \neq j ;$4. selection of the best $\theta_{i}$(largest or smallest).

The first three goals will be achieved by the construction of confidence intervals The number of such confidence intervals is $C = K , C = K - 1$, and $C = K ( K {-} 1 ) / 2$, respectively. Hochberg and Tamhane [1987] and Hsu [1996] are comprehensive references for such multiple-comparison procedures. The fourth goal requires the use ofa type ofstatistical procedure known as a multiple ranking and selection procedure. Procedures to achieve these and other goals are discussed by Kleijnen [1975, Chapters II and V], who also discusses their relative merit and disadvantages. Goldsman and Nelson [1998] and Law [2007] discuss those selection procedures most relevant to simulation. A comprehensive reference is Bechhofer, Santner, and Goldsman [1995]. The next subsection presents a fixed-samplesize procedure that can be used to meet goals 1, 2, and 3, and is applicable in a wide range of circumstances. Subsection 2.2 presents a procedure to achieve goal 4.

2.1 Bonferroni Approach to Multiple Comparisons

Suppose that C confidence intervals are computed and that the ith interval has confidence coefficient $1 - \alpha_{i}$. Let $\boldsymbol {S} _ {i}$ be the statement that the ith confidence interval contains the parameter (or difference of two parameters) being estimated. This statement might be true or false for a given set of data, but the procedure leading to the interval is designed so that statement $\boldsymbol {S_{i}}$ will be true with probability $1 - \alpha_{i}$. When it is desired to make statements about several parameters simultaneously, as in goals 1, 2 and 3, the analyst would like to have high confidence that all statements are true simultaneously. The Bonferroni inequality states that

where $\begin{array} {r} {\alpha_{E} = \sum_{j = 1} ^ {C} \alpha_{j}} \end{array}$ is called the overall error probability. Expression (16) can be restated as

or equivalently,

Thus,$\alpha_{E}$ provides an upper bound on the probability of a false conclusion. To conduct an experiment that involves making C comparisons, first select the overall error probability, say,$\alpha_{E} = 0 . 05 \mathrm{o r} 0 . 10 .$ The individual α may be chosen to be equal $( \alpha_{i} = \alpha_{E} / C )$, or unequal, as desired. The smaller the value of $\alpha_{j}$, the wider the jth confidence interval will be. For example, if two 95% c.i.s $( \alpha_{1} = \alpha_{2} =$0.05) are constructed, the overall confidence level will be 90% or greater $( \alpha_{E} = \alpha_{1} + \alpha_{2} = 0 . 10 )$. If ten 95% c.i.s are constructed $( \alpha_{i} = 0 . 05 , i = 1 , \ldots , 10 )$, the resulting overall confidence level could be as low as 50%$\textstyle ( \alpha_{E} = \sum_{i = 1} ^ {10} \alpha_{i} = 0 . 50 )$, which is far too low for practical use. To guarantee an overall confidence level of 95%, when 10 comparisons are being made, one approach is to construct ten 99.5% confidence intervals for the parameters (or differences) of interest.

The Bonferroni approach to multiple confidence intervals is based on Expression (16). A major advantage is that it holds whether the models for the alternative designs are run with independen sampling or with common random numbers.

The major disadvantage of the Bonferroni approach in making a large number of comparisons is the increased width of each individual interval. For example, for a given set of data and a large sample size, a 99.5% c.i. will be $z_{0 . 0025} / z_{0 . 025} = 2 . 807 / 1 . 96 = 1 . 43$ times longer than a 95% c.i. For smal sample sizes, say, for a sample of size $5 ,$ a 99.5% c.i. will be $t_{0 . 0025 , 4} / t_{0 . 025 , 4} = 5 . 598 / 2 . 776 = 1 . 99$ times longer than an individual 95% c.i. The width of a c.i. is a measure of the precision of the estimate. For these reasons, it is recommended that the Bonferroni approach be used only when a small number of comparisons are being made. Twenty or so comparisons appear to be the practica upper limit.

Corresponding to goals 1, 2, and 3, there are at least three possible ways of using the Bonferron Inequality (16) when comparing K alternative system designs:

1. Individual $c . i . s {\mathrm{:}}$ Construct a $100 ( 1 - \alpha_{i} ) \%$ c.i. for parameter $\theta_{i}$ by using the expression $H = t_{\alpha / 2 , R - 1} \frac {S} {\sqrt {R}} ,$ in which case the number of intervals is $C = K$. If independent sampling were used, the K c.i.s would be mutually independent, and thus the overall confidence level would be $( 1 - \alpha_{1} ) \times ( 1 - \alpha_{2} ) \times \cdot \cdot \cdot \times ( 1 - \alpha_{C} )$, which is larger, but not by much, than the right side of Expression (16). This type of procedure is most often used to estimate multiple parameters of a single system, rather than to compare systems; because multiple parameter estimates from the same system are likely to be dependent, the Bonferroni inequality typically is needed.

2. Comparison to an existing system: Compare all designs to one specific design, usually to an existing system. That is, construct a $100 ( 1 - \alpha_{i} ) \% \ \mathrm{c . i}$. for $\theta_{i} - \theta_{1} ( i = 2 , 3 , \ldots , K )$, using Expression (1). (System 1 with performance measure $\theta_{1}$ is assumed to be the existing system.) In this case, the number of intervals is $C = K - 1$. This type of procedure is most often used to compare several competitors to the present system in order to learn which are better.

3. All pairwise comparisons: Compare all designs to each other. That is, for any two system designs $i \neq j ,$ construct a $100 ( 1 - \alpha_{i i} ) \%$ c.i. for $\theta_{i} - \theta_{j}$. With K designs, the number of confidence intervals computed is $C = \dot {K} ( K - 1 ) / 2$. The overall confidence coefficient would be bounded below by $\begin{array} {r} {1 - \alpha_{E} = 1 - \sum \sum_{i \neq j} \alpha_{i j}} \end{array}$[which follows by Expression (16)]. It is generally believed that CRN will make the true overall confidence level larger than the right side of Expression (16), and usually larger than will independent sampling. The right side of Expression (16) can be thought of as giving the worst case, that is, the lowest possible overal confidence level.

Example 2

Reconsider the call center design problem of Example 1. The alternative system designs are the following:

a. C: existing system (distinct operator queues);

b. P7: cross-trained system with 7 operators (single queue);

c. P6: cross-trained system with 6 operators (single queue).

Using the data in Table 2, confidence intervals for $\theta_{C} - \theta_{P 7} , \theta_{C} - \theta_{P 6}$, and $\theta_{P 7} - \theta_{P 6}$ will be constructed, having an overall confidence level of 95%. Recall that CRN was used in all models, bu this does not affect the overall confidence level, because, as mentioned, the Bonferroni Inequality (16) holds regardless of the statistical independence or dependence of the data.

Since the overall error probability is $\alpha_{E} = 0 . 05$ and $C = 3$ confidence intervals are to be constructed, let $\alpha_{i} = 0 . 05 / 3 = 0 . 0167$. Then use Expression (1), with proper modifications, to construct $C = 3$ confidence intervals with $\alpha = \alpha_{i} = 0 . 0167$ and degrees of freedom $\nu = 10 - 1 = 9$ The value of $t _ {\alpha_{i} / 2 , R - 1} = t_{0 . 0083 , 9} = 2 . 97$ is obtained from Table A.5 in the Appendix by interpolation; the point estimates and standard errors are obtained from Table 2

The three confidence intervals, with overall confidence coefficient at least 95%, are given by

or

The simulation analyst has high confidence (at least 95%) that all three confidence statements are correct. Notice that the c.i. for $\theta_{C} - \theta_{P 7}$ again contains zero; thus, there is no statistically significan difference between the current design and the alternative with 7 cross-trained operators, a conclusion that supports the previous results in Example 1. The other confidence intervals lie completely to the left of 0, indicating that C and P7 both dominate P6.

Some of the exercises at the end of this chapter provide an opportunity to compare CRN and independent sampling and to compute simultaneous confidence intervals under the Bonferroni approach.

2.2 Selection of the Best

Suppose that there are K system designs, and the unknown expected value of the ith system’s performance is θ . We are interested in which system design is best, where “best” means having the maximum or minimum θ_i, depending on the problem. For instance, in Example 3 we will compare $K = 8$ designs for the semiconductor fabrication facility described in the Estimation of Absolute Performance chapter, in which $\theta_{i}$ is the steady-state mean cycle time of one product family for the ith design, and smaller $\theta_{i}$ is better. We want a procedure that is capable of handling large K, say, on the order of 100

Let B denote the (unknown) index of the best system design. The smaller the true differences $| \theta_{B} - \theta_{i} | , i \neq B $ are, and the more certain we want to be that we find the best system, the more replications are required to achieve our goal. Therefore, instead of demanding that we find the best design B no matter what, we compromise and ask to find B with high probability whenever the difference between system B and the others is significantly large. More precisely, we want the probability that we select the best system to be at least 1 α whenever $| \theta_{B} - \theta_{i} | \geq \epsilon $ for all $i \neq B ,$ where  depends on the problem. If there are systems that are within  of the best, then we will be satisfied to select either the best or any one of the near-best designs. Both the probability of correc selection $1 - \alpha$ and the practically significant difference  will be under our control.

The following procedure accomplishes this goal by undertaking two stages of simulation: In the first stage,$R_{0}$ replications are obtained from each of the system designs. Those systems that are

statistically significantly inferior to the others are screened out, meaning that they are eliminated from further consideration, with confidence $1 - \alpha / 2$. If more than one system survives screening, then the survivors receive enough additional replications to select the best, or a system design within  of the best, with confidence $1 - \alpha / 2$. Then, using reasoning related to the Bonferroni inequality, the selected system is the best or a near-best system with confidence at least $1 - \alpha$ for the combination of screening and selection.

The first-stage screening can be very important because, as Example 3 will illustrate, the number of replications required in the second stage to select the best can be quite large, and there is no reason to waste time and effort on system designs that are obviously not competitive. This is especially important when K is large.

The procedure below is implemented in SimulationTools.xls, which is available at www.bcnn.net. The SimulationTools.xls implementation guides the user as to what data need to be obtained for each stage.

Select-the-Best Procedure

1. Specify the desired probability ofcorrect selection $1 / k < 1 - \alpha < 1$, the practically significan difference $\epsilon > 0$, an initial number of replications $R_{0} \geq 10$, and the number of competing systems K. Set

and obtain Rinott’s constant $h = h ( R_{0} , K , 1 - \alpha / 2 )$. [Hint: Both critical values are calculated automatically by SimulationTools.xls; small tables of t and h are given in the Appendix.]

2. Make $R_{0}$ replications of each system. Calculate the first-stage sample means and sample variances

for $i = 1 , 2 , \dots , K .$3. Calculate the screening thresholds

for all $i \neq j$(a) If bigger is better, then form the survivor subset $s$ containing every system design i such that

That is, retain every system i whose sample mean is no more than max 0,$W_{i j} - \epsilon \}$ smaller than any other system’s sample mean.

(b) If smaller is better, then form the survivor subset containing every system design i such that

That is, retain every system i whose sample mean is no more than max 0,$W_{i j} - \epsilon \}$ larger than any other system’s sample mean.

4. If the set contains only one system, then stop and return that system as the best. Otherwise, for all designs i in , compute the second-stage sample sizes

where means to round up.

5. Take $R_{i} - R_{0}$ additional replications from all systems i in , or, if it is more convenient, obtain a total of $R_{i}$ replications from system i by starting over.

6. Compute the overall sample means

for all i in . If bigger is better, then select the system with the largest $\bar {\bar {Y}} _ {\cdot i}$ as best. Otherwise select the system with the smallest $\bar {\bar {Y}} _ {\cdot i}$ as best.

The critical values t and h cannot usually be obtained directly from tables so interpolation is often required. As an alternative, algorithms for computing quantiles from the t distribution can be found in most numerical analysis and spreadsheet software; unfortunately, this is not the case for Rinott’s constant h. Therefore, use of SimulationTools.xls is highly recommended for this procedure.

Nelson et al. [2001] prove that the Select-the-Best Procedure finds the system with the larges or smallest θ , or one within  of the best θ, with confidence level at least 1 α, provided the data are normally distributed and the systems are simulated independently. The procedure will also work if CRN is applied.

Example 3

Recall the FastChip, Inc. semiconductor manufacturing example from the Estimation of Absolute Performance chapter, in which FastChip, Inc. wants to use simulation to design a new “fab” by evaluating the steady-state mean cycle time (time from wafer release to completion) that they can expect under a particular loading of two distinct products, called C-chip and D-chip. The system design we evaluated in the Estimation of Absolute Performance chapter is actually a base case, and there are seven additional designs that arise by making different choices about how to implemen certain steps in the manufacturing process; these steps are called Clean, Load Quartz, Oxidize, Unload Quartz, Coat, and Develop. The question is, which design leads to the smallest steady-state mean cycle time? We will focus on the cycle time for C-chip, and we will use the warmup period and run length obtained there.

Figure 5 Input form for Select-the-Best Procedure in SimulationTools.xls

Including the base case there are $K \ = \ 8$ system designs, and the performance measures $\theta_{1} , \theta_{2} , \ldots , \theta_{8}$ are the steady-state mean cycle times for C-chip under each design. FastChip decides that even very small improvements in cycle time are important, so they set $\epsilon = 0 . 15$ hours (or 9 minutes). This is a very tight requirement when cycle times average around 45 hours, and we wil see the consequences of this choice below. FastChip would like 95% confidence they have selected the best design, or one within 9 minutes of the best. Figure 5 shows the input screen to set up this problem in SimulationTools.xls, which uses the term “indifference level” for the practically significant difference .

Table 3 contains the first-stage results obtained from $R_{0} ~ = ~ 10$ replications of each system design. The alternative that adjusts the Oxidize step seems to provide the smallest average cycle time, but all of the alternatives appear close. The complete data from this example are available at www.bcnn.net.

The t value we need for screening is

Table 3 First-Stage Data from FastChip Simulation

Table 4 Screening Threshold from FastChip Simulation

The subset will always contain the sample best system, which is system $i = 4$(Oxidize) here. To avoid being screened out, the other systems’ sample means must be small enough to satisfy $\bar {Y} _ {\cdot i} \leq \bar {Y} _ {\cdot j} + \operatorname* {m a x} \{W_{i j} - \epsilon \}$ for all $j \neq i ,$, and in particular $\bar {Y} . i \le \bar {Y} . 4 + \operatorname* {m a x} \{0 , W_{i 4} - \epsilon \}$. Table 4 contains the thresholds ${\bar {Y}} _ {\cdot 4} +$ max $\{0 , W_{i 4} - \epsilon \}$; only alternative 2 with sample mean 45.70 hours is smal enough to survive, so ${\mathcal {S}} = \{2 , 4 \}$ To determine the second-stage sample sizes for alternatives 2 and 4, we need $h = h ( R_{0} , K , 1 -$$\alpha / 2 ) = h ( 10 , 8 , 0 . 975 ) = 4 . 635$, which could be obtained from SimulationTools.xls. Then $R_{2} = \lceil ( h S_{2} / \epsilon ) ^ {2} \rceil = 349$ and $R_{4} = \lceil ( h S_{4} / \epsilon ) ^ {2} \rceil = 90$. Thus, alternative 2 requires $349 - 10 = 339$ additional replications, while alternative 4 requires $90 - 10 = 80$ additional replications. Notice that the number of replications for alternative 2 is quite large. In general, the second-stage sample size is large if the sample variance is large, making it is difficult to detect differences, or if  is small, even small differences matter, as in this example.

After obtaining the additional replications,$\bar {\bar {Y}} _ {\cdot 2} = 45 . 88$ and $\bar {\bar {Y}} _ {\cdot 4} = 45 . 00$; therefore, alternative 4, Oxidize, is selected. The procedure guarantees, with 95% confidence, that $\theta_{4}$ is the smallest steadystate mean cycle time of the 8 alternatives,$^ {\mathrm{o r ,}}$ if it is not, then it is within $\epsilon = 0 . 15$ hours of being the best.

The Select-the-Best Procedure, as presented here, does both screening and selection. This makes it particularly useful for finding the best system design from among a large number of alternatives since many of them may be screened out in the first stage. The procedure can be used for screening alone, without selection, if the objective is to eliminate designs that are not competitive. This migh be worthwhile when we want to choose from a collection of competitive alternatives based on criteria beyond mean performance alone. For screening only, the procedure stops at Step 3 and can use the critical value

which provides somewhat tighter screening.

Similarly, if the number of system designs K is small, and the goal is to select the best, then screening can be skipped and all systems receive second-stage replications. For selection only, skip Step 3 and use the somewhat smaller critical value $h = h ( R_{0} , K , 1 - \alpha )$.

3 Metamodeling

Suppose that there is a simulation output response variable Y that is related to k independent variables, say,$x_{1} , x_{2} , \ldots , x_{k}$. The dependent variable Y is a random variable, while the independent variables $x_{1} , x_{2} , \ldots , x_{k}$ are called design variables and are usually subject to control. The true relationship between the variables Y and x is represented by the often complex simulation model. Our goal is to approximate this relationship by a simpler mathematical function called a metamodel. In some cases, the analyst will know the exact form of the functional relationship between Y and $x_{1} , x_{2} , \ldots , x_{k} , \mathrm{s a y} .$$Y = f ( x_{1} , x_{2} , \dots , x_{k} )$. However, in most cases, the functional relationship is unknown, and the analys must select an appropriate f containing unknown parameters, and then estimate those parameters from a set of data Y, x . Regression analysis is one method for estimating the parameters.

Example 4

An insurance company promises to process all claims it receives each day by the end of the next day. It has developed a simulation model of its proposed claims-processing system to evaluate how hard it will be to meet this promise. The actual number and types of claims that will need to be processed each day will vary, and the number may grow over time. Therefore, the company would like to have a model that predicts the total processing time as a function of the number of claims received.

The primary value of a metamodel is to make it easy to answer “what-if” questions, such as what the processing time will be if there are x claims. Evaluating a function f, or perhaps its derivatives, at a number of values of x is typically much easier than running a simulation experiment for each value.

3.1 Simple Linear Regression

Suppose that it is desired to estimate the relationship between a single independent variable x and a dependent variable Y, and suppose that the true relationship between Y and x is suspected to be linear. Mathematically, the expected value of Y for a given value of x is assumed to be

where $\beta_{0}$ is the intercept on the Y axis, an unknown constant, and $\beta_{1}$ is the slope, or change in Y for a unit change in x, also an unknown constant. It is further assumed that each observation of Y can be described by the model

where ε is a random error with mean zero and constant variance $\sigma^{2}$. The regression model given by Equation (18) involves a single variable x and is commonly called a simple linear regression model.

Suppose that there are n pairs of observations $( Y_{1} , x_{1} ) , ( Y_{2} , x_{2} ) , \ldots , ( Y_{n} , x_{n} )$. These observations can be used to estimate $\beta_{0}$ and $\beta_{1}$ in Equation (18). The method of least squares is commonly used to form the estimates. In the method of least squares,$\beta_{0}$ and $\beta_{1}$ are estimated in such a way that the sum of the squares of the deviations between the observations and the regression line is minimized. The individual observations in Equation (18) can be written as

where $\varepsilon_{1} , \varepsilon_{2} , \ldots$ are assumed to be uncorrelated random variables.

Each $\varepsilon_{i}$ in Equation (19) is given by

and represents the difference between the observed response $Y_{i}$ and the expected response $\beta_{0} + \beta_{1} x_{i}$ predicted by the model in Equation (17). Figure 6 shows how $\varepsilon_{i}$ is related to $x_{i} , Y_{i} ,$, and $ E ( Y_{i} | x_{i} ) $ The sum of squares of the deviations given in Equation (20) is given by

and $L$ is called the least-squaresfunction. It is convenient to rewrite $Y_{i}$ as follows:

where $\beta_{0} ^ {\prime} = \beta_{0} + \beta_{1} \bar {x}$ and $\bar {x} = \textstyle \sum_{i = 1} ^ {n} x_{i} / n$. Equation (22) is often called the transformed linear regression model. Using Equation (22), Equation (21) becomes

Figure 6 Relationship of $\varepsilon_{i}$ to $x_{i} , Y_{i} ,$ and $ E ( Y_{i} | x_{i} ) $ To minimize $L ,$ find $\partial L / \partial \beta_{0} ^ {\prime}$ and $\partial L / \partial \beta_{1}$, set each to zero, and solve for $\widehat {\beta} _ {0} ^ {\prime}$ and ${\widehat {\beta}} _ {1}$. Taking the partial derivatives and setting each to zero yields

Equations (23) are often called the normal equations, which have the solutions

and

The numerator in Equation (25) is rewritten for computational purposes as

where $S_{x y}$ denotes the corrected sum of cross products of x and Y. The denominator of Equation (25) is rewritten for computational purposes as

where $S_{x x}$ denotes the corrected sum of squares of x. The value of ${\widehat {\beta}} _ {0}$ can be retrieved easily as

 The simulation model of the claims-processing system in Example 4 was executed with initia conditions $x = 100$, 150, 200, 250, and 300 claims received the previous day. Three replications were obtained at each setting. The response Y is the number of hours required to process x claims. The results are shown in Table 5. The graphical relationship between the number of claims received and total processing time is shown in Figure 7. Examination of this scatter diagram indicates that there is a strong relationship between number of claims and processing time. The tentative assumption o the linear model given by Equation (18) appears to be reasonable.

Table 5 Simulation Results for Processing Time Given x Claims

Figure 7 Relationship between number of claims and hours of processing time.

With the processing times as the $Y_{i}$ values (the dependent variables) and the number of claims as the $x_{i}$ values (the independent variables),$\widehat {\beta} _ {0}$ and ${\widehat {\beta}} _ {1}$ can be found by the following computations:

and $\bar {x} = 3000 / 15 = 200$ From Equation (26) $S_{x y}$ is calculated as

From Equation (27),$S_{x x}$ is calculated as

Then,${\widehat {\beta}} _ {1}$ is calculated from Equation (25) as

As shown in Equation (24),$\widehat {\beta_{0} ^ {\prime}}$ is just ${\bar {Y}} ,$ or

To express the model in the original terms, compute ${\widehat {\beta}} _ {0}$ from Equation (28) as

Then an estimate of the mean of Y given x,$ E ( Y | x ) $, is given by

For a given number of claims x this model can be used to predict the number of hours required to process them. The coefficient ${\widehat {\beta}} _ {1}$ has the interpretation that each additional claim received adds an expected 0.0464 hours, or 2.8 minutes, to the expected total processing time.

Regression analysis is widely used and frequently misused. Several of the common abuses are briefly mentioned here. Relationships derived in the manner of Equation (29) are valid for values of the independent variable within the range of the original data. The linear relationship that has been tentatively assumed may not be valid outside the original range. In fact, we know from queueing theory that mean processing time may increase rapidly as the number of claims approaches the capacity of the system. Therefore, Equation (29) can be considered valid only for $100 \leq x \leq 300$ Regression models are not advised for extrapolation purposes.

Care should be taken in selecting variables that have a plausible causal relationship with one another. It is quite possible to develop statistical relationships that are unrelated in a practical sense. For example, an attempt might be made to relate monthly output of a steel mill to the weight of reports appearing on a manager’s desk during the month. A straight line may appear to provide a good mode for the data, but the relationship between the two variables is tenuous. A strong observed relationship does not imply that a causal relationship exists between the variables. Causality can be inferred only when analysis uncovers some plausible reasons for its existence. In Example 4 it is reasonable tha starting with more claims implies that more time is needed to process them. Therefore, a relationship of the form of Equation (29) is at least plausible.

Testing for Significance of Regression

In Section 3.1, it was assumed that a linear relationship existed between Y and x. In Example 4, a scatter diagram, shown in Figure 7, relating number of claims and processing time was prepared, to evaluate whether a linear model was a reasonable tentative assumption prior to the calculation of $\widehat {\beta} _ {0}$ and ${\widehat {\beta}} _ {1}$. However, the adequacy of the simple linear relationship should be tested prior to using the model for predicting the response $Y_{i}$ given an independent variable,$x_{i} .$. There are several tests which may be conducted to aid in determining model adequacy. Testing whether the order of the mode tentatively assumed is correct is commonly called lack-of-fit testing.

Testing for the significance of regression provides another means for assessing the adequacy of the model. In particular, tests to establish that the estimated regression coefficients $\widehat {\beta} _ {0}$ and ${\widehat {\beta}} _ {1}$ are statistically different from 0 are essential, particularly for multiple linear regression, described below, when a large number of design variables x are included in the metamodel but it is not certain that all of them have a substantial explanatory effect.

Tests for significance and lack of fit are included in all regression software. See Box and Draper [1987], Hines et al. [2002], and Montgomery [2000] for details about these tests.

Multiple Linear Regression

If the simple linear regression model of Section 3.1 is inadequate, several other possibilities exist. There could be several independent variables, so that the relationship is of the form

Notice that this model is still linear, but has more than one independent variable. Regression models having the form shown in Equation (30) are called multiple linear regression models. Another possibility is that the model is of a quadratic form such as

Equation (31) is also a linear model which may be transformed to the form of Equation (30) by letting $x_{1} = x$ and $x_{2} = x^{2}$ Yet another possibility is a model of the form such as

which is also a linear model. The analysis of these three models with the forms just shown, and related models, can be found in Box and Draper [1987], Hines et al. [2002], Montgomery [2000], and other applied statistics texts; and also in Kleijnen [1987, 1998], which is concerned primarily with the application of these models in simulation.

3.2 Metamodeling and Computer Simulation

The material in Section 3.1 applies to regression metamodeling in general, not just computer simulation. However, there are special issues and opportunities that arise when goal is to provide a metamodel for a computer simulation; we describe these below.

Random-Number Assignment for Regression

The assignment of random-number seeds or streams is part of the design of a simulation experiment. (This section is based on Nelson [1992].) Assigning a different seed or stream to different design points (settings for $x_{1} , x_{2} , \ldots , x_{m}$ in a multiple linear regression) guarantees that the responses Y from different design points will be statistically independent. Similarly, assigning the same seed or stream to different design points induces dependence among the corresponding responses, by virtue of their all having the same source of randomness.

Many textbook experiment designs assume independent responses across design points. To conform to this assumption, we must assign different seeds or streams to each design point. However, it is often useful to assign the same random number seeds or streams to all of the design points, in other words, to use common random numbers.

The intuition behind common random numbers for metamodels is that a fairer comparison among design points is achieved if the design points are subjected to the same experimental conditions, specifically the same source of randomness. The mathematical justification is as follows: Suppose we fit the simple linear regression $Y_{i} = \beta_{0} + \beta_{1} x_{i} + \varepsilon_{i}$ and obtain least squares estimates $\widehat {\beta} _ {0}$ and ${\widehat {\beta}} _ {1}$Then an estimator of the expected difference in performance between design points i and j is

when $x_{i}$ and $x_{j}$ are fixed design points,${\widehat {\beta}} _ {1}$ determines the estimated difference between design points i and $j ,$or for that matter between any other two values of x. Therefore, common random numbers can be expected to reduce the variance of ${\widehat {\beta}} _ {1}$ and, more generally, reduce the variance of all of the slope terms in a multiple linear regression. Common random numbers typically do not reduce the variance of the intercept term $\widehat {\beta} _ {0}$The least-squares estimators ${\widehat {\beta}} _ {0}$ and ${\widehat {\beta}} _ {1}$ are appropriate regardless of whether we use common  random numbers, but the associated statistical analysis is affected by that choice. For statistica analysis of a metamodel under common random numbers, see Kleijnen [1988] and Nelson [1992].

Queueing Simulation

Many computer simulation models are of systems that involve queues or networks of queues, and metamodels of these simulations may include design variables related to capacity (numbers of servers and service rates) and customer load (arrival rates). Such models violate standard regression assumptions as described in Section 3.1. A small example illustrates the difficulties and possible solutions.

Recall the $M / M / 1$ queue: A single-server, first-come-first-served queue with a Poisson arrival process having rate λ customers/time and exponentially distributed service time having rate $\mu$ customers/time. Let $x = \lambda / \mu$, which is called the utilization $\rho .$. Then, the steady-state expected number of customers in the queue for $0 < x < 1$ is

Now, suppose that we did not know the relationship (32) and were trying to build a metamodel for L as a function of $0 < x < 1$ by running simulations of an $M / M / 1$ queue. Clearly the usual low-order polynomial models, including only the terms x and $x^{2}$ described above, will not fit very well, particularly for x close to 1, where the higher-order terms matter.

There are at least two solutions to this problem. One is to break up the design space $( 0 < x < 1$ in this example) into smaller pieces for which low-order polynomials are adequate. The second approach is to fit a nonlinear model motivated by queueing theory. Yang, Ankenman, and Nelson [2007] showed that

with $m \ \leq \ 3 .$, works well for many queueing outputs Y, including average number in the queue and average time in the queue. Notice that the unknown parameters are m,$\beta_{i} .$, and $p ,$ and therefore nonlinear least squares methods are required to fit them (Bates and Watts [1988]). When Y is the average number in queue, an estimator of L, then $m = 1 , \beta_{0} = 0 , \beta_{1} = 1$, and $p = 1$ Unfortunately, the problems go even deeper. If we estimate L by making a simulation run of length T time units and record the average number of customers in queue (call it Y), then it can be shown (Whitt [1989]) that

Thus, instead of the variance being constant across the design space, as assumed in regression, it increases, and increases explosively, as x approaches 1. As a consequence, a run of, say,$T = 100$ time units may yield a pretty good estimate of L at $x = 0 . 1$ but a very poor estimate at $x = 0 . 9$ The standard statistical inference produced by regression software will be invalid when response variance depends strongly on x, meaning we may not recognize lack of fit or poor fit.

A remedy for unequal response variance in regression is to apply a variance-stabilizing transformation to the data to make the variance more nearly constant across the design space. In simulation we can attack the problem directly by using different run lengths $T$ or different numbers of replications, depending on x.

We suggest using the $M / M / 1$ result as a rough guide for more general networks of queues. If we let x denote the network utilization, that is, the most heavily utilized queue, then Approximation (34) implies that the run lengths $T_{1}$ and $T_{2}$ at two design points $x_{1}$ and $x_{2}$ should satisfy

For instance, if we were to simulate $T_{1} = 1$ time unit at utilization $x_{1} = 0 . 1$, then we would need about $T_{2} = 10$ time units at $x_{2} = 0 . 3$, or $T_{2} = 72$ time units at $x_{2} = 0 . 5$, or $T_{2} = 876$ time units at $x_{2} = 0 . 7 .$, or $T_{2} = 101$, 944 time units at $x_{2} = 0 . 9$ to stabilize variances.

The central message is this: In queueing simulation, if the metamodel includes decision variables that affect customer load—by changing the arrival rates, service rates, or numbers of servers—then a low-order polynomial may not be a good choice for the metamodel, and equal response variance across the design space cannot be assumed without exercising run-length control

4 Optimization via Simulation

Consider the following examples. (Some of these descriptions are based on Boesel, Nelson, and Ishi [2003].

Example 6: Materials Handling System (MHS)

Engineers need to design an MHS consisting of a large automated storage and retrieval device, automated guided vehicles (AGVs),AGV stations, lifters, and conveyors.Among the design variables they can control are the number ofAGVs, the load perAGV, and the routing algorithm used to dispatch the AGVs. Alternative designs will be evaluated according to AGV utilization, transportation delay for material that needs to be moved, and overall investment and operation costs.

Example 7: Liquefied Natural Gas (LNG) Transportation

An LNG transportation system will consist of LNG tankers and of loading, unloading, and storage facilities. In order to minimize cost, designers can control tanker size, number of tankers in use, number of jetties at the loading and unloading facilities, and capacity of the storage tanks.

Example 8: Automobile Engine Assembly

In an assembly line, a large buffer (queue) between workstations could increase station utilization— because there will tend to be something waiting to be processed—but drive up space requirements and work-in-process inventory. An allocation of buffer capacity that minimizes the sum of these competing costs is desired.

Example 9: Traffic Signal Sequencing

Civil engineers want to sequence the traffic signals along a busy section of road to reduce driver delay and the congestion occurring along narrow cross streets. For each traffic signal, the length of the red, green, and green-turn-arrow cycles can be set individually.

Example 10: Online Services

A company offering online information services over the Internet is changing its computer architecture from central mainframe computers to distributed workstation computing. The numbers and types of CPUs, the network structure, and the allocation of processing tasks all need to be chosen. Response time to customer queries is the key performance measure.

What do these design problems have in common? Clearly, a simulation model could be useful in each, and all have an implied goal of finding the best design relative to some performance measures (cost, delay, etc.). In each example, there are potentially a very large number of alternative designs, ranging from tens to thousands, and certainly more than the 2 to 100 we considered in Section 2.2. Some of the examples contain a diverse collection of decision variables: discrete (number of AGVs, number of CPUs), continuous (tanker size, red-cycle length), and qualitative (routing strategy, algorithm for allocating processing tasks). This makes developing a metamodel, as described in Section 3, difficult.

All of these problems fall under the general topic of “optimization via simulation,” where the goal is to minimize or maximize some measures of system performance, and system performance can be evaluated only by running a computer simulation. Optimization via simulation is a relatively new but already vast topic, and commercial software has become widely available. In this section, we describe the key issues that should be considered in undertaking optimization via simulation, provide some pointers to the available literature, and give one example algorithm.

4.1 What Does ‘‘Optimization via Simulation’’ Mean?

Optimization is a key tool used by operations researchers and management scientists, and there are well-developed algorithms for many classes of problems, the most famous being linear programming. Much of the work on optimization deals with problems in which all aspects of the system are treated as being known with certainty; most critically, the performance of any design (cost, profit, makespan, etc.) can be evaluated exactly.

In stochastic, discrete-event simulation, the result of any simulation run is a random variable. For notation, let $x_{1} , x_{2} , \ldots , x_{m}$ be the m controllable design variables and let $Y ( x_{1} , x_{2} , \dots , x_{m} )$ be the observed simulation output performance on one run. To be concrete,$x_{1} , x_{2} , x_{3}$ might denote