Tilburg University Generalization of simulation results van den Burg, AJ

16 Pages·2017·2.59 MB·English

Tilburg University
Generalization of simulation results
van den Burg, A.J.; van der Ham, R.Th.; Kleijnen, J.P.C.
Published in:
European Journal of Operational Research
Publication date:
1979
Link to publication in Tilburg University Research Portal
Citation for published version (APA):
van den Burg, A. J., van der Ham, R. T., & Kleijnen, J. P. C. (1979). Generalization of simulation results:
Practicality of statistical methods. European Journal of Operational Research, 3(1), 50-64.
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners
and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.
• Users may download and print one copy of any publication from the public portal for the purpose of private study or research.
• You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal
Take down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately
and investigate your claim.
Download date: 14. mrt.. 2023
Generalization of answer a variety of questions. Besides answers to a
few ct'fic:~ps questions, we shall try to gain a general
simulation results- unders,anding of the system's behavior. eW also wish
to be able to answer "inverse" questions, i.e., in our
case given a specific value for the original response
variable, inventory capacity, which yearly through-
P racticality of
put is possible? eW could try to fmd this value via
statistical methods trial and error, but having available an explicit formu-
la (or nomogram) gives an immediate answer. As we
have already indicated, we also wish to know how
sensitive our solutions are to the model's specifica-
J.P.C. KLEIJNEN tions. Sensitivity studies per factor are traditional in
Kctholieke ,loohcsegoH ,grubliT Netherlands practice, but we may further examine interactions
among model parameters. The larger the area over
A.J. van den BURG which factors vary, the more desirable it seems to
Municipality of Zoetermeer {formerly with ~CT) include interactions in our investigation. (Simple
linear functions provide good approximations locally;
R.Th. van der HAM see Section 5.)
Europe Container Terminus ,)TCE( Rotterdam As Fig. I demonstrates real-life systems can be
modeled by means of a simulation program. The rela-
devieceR 4 January 1978 tionship between the inputs and the outputs of this
desiveR 26 April 1978
simu!ation can in turn be modeled, namely by a
line: regression model. The regression model acts as
The major purpose of this paper is to evaluate the prac- a metamodel or auxiliary model, so that simulation
ileal use of statistical techniques in both the generaliTation
results can be generalized. eW shall derive the regres.
or analysis of simulation results, and the design of simuiation
sion estimates fl with their covariance matrix, when
experiments. This problem is investigated with the help of a
real-life system, namely the container termhus of ECT in using either ordinary or generalized least squares.
Rotterdam. This system is modeled by a simulation program. Generalized least squares may be relevant because
The relationship between the simulation response and its
the simulation responses have different variances. The
input variables is modeled by a linear regresdon model:
validity of the fitted regression model will be tested
metamodel or auxiliary model. The paper summarizes
by an F-statistic, which compare: the residuals), -)~
regression analysis including generalized lea~t squares which
might be used for simulation responses with non-constant to the "natural" deviations specified by var 0'). Non-
v~ances. The validity of the postulated regression metamodel significant effects/) are set zero. The model is further
is tested statistically: F- and t-statistics. The selection of the
verified by comparing some new simulation response
situations to be simulated, is done through experimental
with the response predicted by the metamodel, using
design methodology, permitting both quantitative and quali-
a Student t-test. Next we discuss the selection of the
tative factors. The statistical techniques apply not only to
simulation but also to real-life experiments. design matrix. eW conclude with the actual results of
a practical study performed at Europe Container
Torminus at Rotterdam. More details on our study
I. Introduction can be found in the complete report 15.
The following notation is used. Stochastic (proba-
If our simulation model were developed to answer
a few ad hoc questions, running the shnulation pro-
gram a few times might give the required answers. eW
would still have to perform a number of simulation
runs to ees how sensitive those answers are to the
model specifications. Actu~y, a rather expensive
NOITALUMIS LEDOM
L ,
research tool like simulation is usually meant to
L, E D O M A T E M ~-
© North-Holland Publishing Company
European Journal of Operational Research 3 (1979) 50-64 Fig. .1 Stepwise modeling.
50
J.P.C. Kleinen et .Ia / Generalization of simulation results 51
bilistic) variables are in boldface and matrices in-
eluding vectors are denoted by an arrow. Note that
the techniques of this paper are applicable not only
to simulation experiments but also to real-life exp'~ri-
ments .1
2. A problem at Europe Container Terminus ~, Day t
{a} Required stacking capc~c=ty t~c day t. or tY
Europe Container Terminus (ECT) is a company
in the harbor of Rotterdam that provides facilities
Frequency fi
for loading and discharging containers from deep.sea
vessels into smaller ships called feeders, and into
trucks and railway cars, and vice versa. Also ECT pro-
vides facilities for the storage of containers. It is the
worlds' largest container terminus. Its management
is interested in the handling and storage capacities of
the ECT terminal, i.e., what is the capacity of the
present terminal, and which measures have to be
(b) Frequency of required stack,~gcopoc,t v y,. or f
taken when more cargo is to be handed? More speci-
Fig. 2. Required stacking capacity.
fically, our research team wants to know the relation.
ships between handling capacity and quay length,
arm between storage capacity and yearly "through.
put" (production). Since the handling capacity is
large relative to the storage capacity, it was assumed
I
that the storage or "stackhlg" system could be
bi = h) (i = o, l, ..., 0 (2.1)
studied separately. In preliminary studies it was 0
further found that the storage problem is easier to
where/i is the number of days at whichyt = .i The
study than the transhipment issue, which involves
corresponding cumulative distribution is
more factors and non-linearities. So in this report
we discuss only the storage oeontainers, not their i
handling. (2.2)
/=0
The relationship between stacking capacity and
yearly throughput depends on how this throughput
The 90% quantile, denoted by Yo.9o, is the value not
is realized: A specific yearly production can be gener- exceed by Yt in 90% of the cases, i.e., during 90% of
ated by many small ships or by a few big ships. the time. Hence
Other factors are relevant too. For instance, if the
(2.3)
containers remain longer at ECT, then more stacking P(Yt ~< 3'0.90) = 0.90.
capacity is needed. How do we measure the required Besides the 90% quantile ofy we estimate its 95%
capacity? Let the required storage capacity (capacity quantile since 90% is a rather arbitrary value, its
utilization) on day t be denoted by yt (t = 1, 2, ..., T) maximum value y 1.oo, and its average value 6 (y),
where stochastic variables are boldface. These Yt which equals Yo.5o for symmetrically distributed y.
form a time series; see Fig. 2. Theyt can be organized Consequently, our original question can be formu-
into a frequency diagram, yielding estimates of lated as: Consider the required stacking capacity and
P(Yt = :O determine its average, its 90% and 95% quantfles, and
its maximum, given a specific yearly throughput
1 The main difference is that in real-life experiments no realized by specific ship sizes and ship arrivals, and
subruns exist so that the experimental error variance 0 2 qualified for specific factors such as container °°dwell
must be estimated by replication. Moreover, independence
time" (number of days spent by cont~uers at ECT).
of the observations does not follow from the use of inde-
Reversely, the question can be asked: Given a certain
pendent sequences of random numbers but must be
capacity, how large a yearly productiori can e~
assured by good experimental practice.
52 £t'.C. nen~OelK et el. / Generalization of amulation remits
h~dled? This "inverse" problem t~,med out to cre- the users so that they would be more willing to ac-
am extra complicatiom. cept and implement the results of our study. Unfor-
tunately, ~ulation requires much development
effort and computer nm.ning time. Until recently,
3, Analytical and simulation models simulation models were used in an ad hoc manner at
ECF (as in most other companies), so that it was dif-
Next we are confronted with the choice between ficult to fred relationships among the variables of
an ~ a l and a simulaticn model, We decided interest. This report veal show how we attacked this
that a simulation model was prefered since an analyt- problem!
ical approach requ~ed too drastic simplifications The simulation model was de_nrmargorp in a pro-
such as Poisson service times. Moreover, we felt that cess-oriented simulation language called PROSIM.
a simulation model could be better understood by This language was developed at the University of
u
inn
I,-
|Ge~ATS A PiHS
-
'
.~-.1
DLOH'I , I EUEUO - EROFEB -LAVIRRA
.__._.,._ I
OLOH{ 0~ o,vs ' ., f
~U.~|PK3E~TOR--
I
~.~veioue,us ~,BEi:one AnnlVAL ,. ~
t-
KR60, nOe S, le ='.OALaAVA ? -- >~i
SE(
lie Nor i, tT. JOiN oueu`: ,,WA!Ti~GnOOU- ii
eTAV~S~LAe _ p
~ _ HTnEe r ~ J
r
~OLD ~:~ "OUR i _. _
eTAvlSm~L .... "
-
L . . . . . . . L
~F ! ---|nEACTtVA'i:~ ENARC ~RENNALPr ' ~,
ASS '~N CRANES a r
~ ¢ ~mm~AS S i VAT `: } . . . . . .
ETALUCLACI ROF LLA ~FIHS 0
|t,O~nA'rlON ON:. Oe OPeR "~--~-I"eACT,V'ATS ¢.~Ne.~A.Ne.'; ' t
/
ETAVITCAER SPIHS TA ESEHT LI
STN.~MOM| SEMIT:`MOS( | |
GNILURREVO NA ENITTNEVE 1~ i"D" LEAVe HTREB
, ri i i ,'
ETAVITCAER/ SHIPS NI ,MOOR.GNITIAW |F ~VNA"
'--'---'-"--- C R A N S P L A N N ~ J
i
|JOIN RET-FA----~UEUo SAILING," J
, _ . ' , .
:~T~.LUCLA¢ SUTATS FO ~ A V S
EHT LANIMRET GNISU __ - , , m
SEUEUOI AB _DNA SA ILeAVe 6oEUe -AFTEn SAILING.
t
r ~, i
IF ASKED,' .,TNIRP SUTATS I
iLEAVE SYSTe, T~."MIN A'r,:) ....
- i r i ,SM!
I
l,oL~ s seuoe.e 7
~ ~ - | N S P S G T D R - - - -
Fig. 3. The simulation model.
J.P.C Kleinen et .la / noitaziiare~,eG of simulation results 35
Technology in Delft (Netherlands); see 14. Process factors which would make the analytic~ restdts ap-
simulation makes it easier for the user to read and plicable over a wider range. Unfortunately tiffs ap-
underst~d the simulation program. proach turned out to be too rigid to give satisfactory
As far as arrivals and departures f~c containers are remits for all situations. At that point of time the
concerned, we should distinguish between two types staff of ECT came hat. contact with a new working
of events, namely ship versus truck arrivals and de- group of the Dutch society on operations research,
paxtures. Ships carry relatively large quantifies of studying the application azi simulation of the staffs-
containers, ~hereas trucks transport only one or two tical theory on experimental design. Experimental
containers. Hence, the arrival of a truck has a much design theory has been developed since the 1920"s
smaller effect on the state of the terminal than a and has beer, widely applied in experiments in agri-
ship has. For that reason we modeled the arrivals and culture, chemistry, etc. Its application in management
departures of containers by truck as a continuous and social sciences, however, is still in its infancy. In
stream, coupled with the arrivals and departures of s.el.technical systems the scientific design of experi-
ships. Arrivals of ships, however, are modeled as dis- ments is difficult and expensive (disruption of the
crete events. For more details on the simulation pro- organization). However, in a simulation model of
gram we refer to I 6. Fig. 3 gives a simplified flow such a social system, the experimental factors are
chart of the resulting simulation model. A further completely under the scientist's control. Bonini 1
discussion of the various processes in Fig. 3 can be performed a pioneering study using experimental de-
found in 15. sign techniques in his simulation model of a firm.
Note that collecting data on the real system More recently the use of statistical design and analy.
formed a quite sizable part of our study. These data sis techniques in simulation has been propagated by
are needed for Naylor 11, Mihram I0 and Kleijnen 6. In the
(i) the validation of the simulation model relative foliowing sections we shall hwestigate in detail, how
to the real-world system, more general conclusions can be derived from simula.
(ii) the determination of factor values in subse- ,ion experiments, by the combination of simulation
quent simulation experiments, aimed at investigating models with analytical models and statistical tech-
future developments. .seuqA,t
The simulation output consists of a snapshot every
eight simulated hours of the system state, measured
by a number of variables including utilized storage. 5. Regression metamodels
The transient, initializing phase of the simulation run
was first removed. As Fig. 3 illustrates, from the time Let x denote a factor influencing the outputs of
sefiesYt we compute the average and quantile~ y, the teal-word system. This f~ctor may be qualitative
• ~0.90',~'0.9S arid y 1.o0. or quantitative, continuous or discrete, in Kleijnen
5, pp. 299-313 it is shown h~w ew can represent
a qualitative factor by dummy variables assuming
4. Utilizing simulation results only the values zero or one; if a qualitative factor
assumes only two values then this factor is repre-
As we mentioned in Section ,1 we are not satisfied sented by a single dummy variable with the values
with a few ad hoe simulation runs. eW wish to study minus and plus one. The factors may further be
the sensitivity of our results to the model's specifica- partitioned into decision and environmental variables,
tions, to gain a general understanding of the real-life i.e., variables either under management control or not.
system behavior, to answer "inverse" questions, and The response or output of the real-world system is a
to estimate possible interactions. time series. eW( shall concentrate on a single response
Originally, we tried to combine an analytical solu- variable; in the ca~ of multiple outputs we apply
tion with simulation results as follows. A simple our procedure to each variable separately.) In order
analytical model was developed based on queuing to compare different system configurations, ew
and inventory theory, besides common sense. This character~.e a whole time path by a smgte or a few
an~ytical model yielded results which were valid mea~ares such sa its average, standard deviation, c,r.
only under quite restrictive assumptions. Next, simu- relation coefficients, peak, slope of a fitted linear
!ation experiments were used to determine correction trend, etc. Let y denote a measure which character-
54 .C.PL Kleiinen et ai. / noitazilareneG of simulation remits
izcs a time path of the real-world system. Hence the Ely}
response variable y is a fuaetion of the factors x: Xl = +1
(5.0
y = f (x ..., . Xl = -1
This system is approximated by a simulation model,
wherey is a function elk factors x/q = 1, ..., k),
plus a vector of random numbers 7(remember that
vectors and matrices are denoted by -+). Hence
X 2
(a) 2131 = O: no interaction ¢>
(5.2)
r-3
E|yt
-1
X t =
ereh:,,,, k is much sm~er than the unknown m in eq.
(5.1) and g symbolizes the joint effect of all factors
x ,'n eq. (5.1) not explicitly represented in eq. (5.2). X I = -1
The simulation model is specified by a computer
program which is denoted by.f~ in ¢q. (5.2). The
simulation model may be approximated in turn by a
metamodel, within a specific experimental area E.
I~x 2
io} ~1~ > :O complin~ntary
eW decided to use a metamodd that is linear in
its parameters ~, but not necessarily in its indepea-
~:{y}
dent variables since terms like x~ and x can be
Xl =
.1
utilized. A very simple metzanodel to express the
effects of the k factors would be: X 1 = -1
f5.3)
Yt = o3f + ~txi~ + "'" + f3~z,~ + ~e
(i = I, ..., n)
where in simulation run i (observation 0, factor has X 2
(el ~2 t :0 suDstitvhon _ _ _ _ _ _
the value xi/. These values xii determmeyi linearly
Fig. 4. Interactions.
except for ei, the noise (disturbance, error) term
which has expectation zero.
Such a simple metamodel implies ~hat a change ah
when more ofxl si available which can be substituted
x has a constant effect on the expected response,
for x2.
6(y), namely gj. A more general metamodel assu~es
As we shall see later on, we start by assuming a
that the effect of factor j also depends on the values
metarnodel like eq. (5.4), but next we test statisti-
of the other factors ' (j° ~ ). This can be formalized
cally whether this assumption was realistic! If the
as in eq. (5.4) where for illustration purposes we
assumption turns out to be unreasonable, we have
take k = 3:
several alternatives:
(1) Make the metamodel more complicated by
+ (f312xilxi2 + f3~3xnxl3 + f323xi2xi3) + ei, adding terms such as three-factor interactions. Since
rite intuitive interpretation of such a model ~i diffi-
(5.4)
(i = 1, ..., n).
cult, we prefer the next alternative.
Here the parameters or coefficients ~t2,/313 and 323 (2) Look for transformations of x. For instance,
denote interactions between the factors 1 and 2, t in our study we found that the response variable y
and 3, and 2 and 3. A graphical illustration of inter- reacted non.linearly to the variable x2 (inter-arrival
action ah the case of two factors, is shown in Fig. 4: time). However, the simple transformation I/x2
In case (a) :he curves are paraUe, i.e. the effect of x2 (interarfival rate) resulted in a linear response fu~c-
on 60') does not depend on the level of x -a In case tion (keeping all other variables constant). eW
(b) the two factors are complementary, i.e., as x2 strongly recommended to look for transformations,
increases the increa~ ah E(.y) is stimulated when the from the very beginning of the study. If all rise fails,
increase in x2 si accompanied by an hlcrease in x -1 we may select option (3).
In ca~ (c) the marginal output of x2 si much smaller (3) Reduce the experimental area .E In that ease
£P.C Kleijnen etal. /Generalization o rimulation results 55
the conclusions of the simulation experL, nents have can be made statistically independent indeed, hy
to be le~ general. Note that, if the only purpose of using different sequences of random numbers per
the metamodel is to f'md the optimal values of the simulation run. Unfortunately, the observat~onsyi
x's, then a small area E is used, a metamedd fitted, may show completely different (heterogeneous)
and the direction of better x-values is determined; variances .2/o Hence we introduce a diagonal matrix
see 5 for details on this so-called Response Surface D with elements ~o on the main diagonal. It is straight-
Methodology (RSM). forward to derive that ,vi this case:
After we have used the metamodd to meet the
d~ands of sensitivity analysis, optimization, and so ~2~=(X' • X-~ )-'X'" ~ D" ~- X" (X'- ,~0 -1 . (6.4)
on, we can return to the original simulation model
This rather complicated looking expression can be
to study the system behavior in detail. More on vari-
easily calculated using a computer program that reads
ous types of metamodds can be found in 7. A ..4, ..+
in the values of X and D. Note that if we select an
practical application besides the prbsent paper, is 8.
orthogonal matrix X, i.e.
The parameters ~t in the above equations can be
--~ -9,
x' x = hie" (6.5)
est~ated and tested, using the techniques of regres-
sion analysis. The main ideas of such an analysis will
then the estimated effects ~ remain correlated when
be summarized in the next section.
the responses y have non-constant variances;see 15,
p.33. An example of an orthogonal matrix X will be
given in Table 1 of Section 8.
6. Ordinary and generalized least squares
Note ~hat numerical accuracy aspects are discussed
in 9. These aspects become relevant if X is ill-con-
In metamodels such as equation (5.4) we have one
ditioned (almost singular). In that case the variances
overall or grand mean ,oO k main effects/~, and
of the estimators are large. See also 4, p.166 and
k(k- 1)/2 two-factor interactions ,#~( q </'), all
13, pp.4!-53.
together (say) q parameters. It is convenient to de-
An alternative to ordinary least squares (OLS) si
note these parameters by #t through #q. The values
generalized least squares (GLS). If the "classical"
of the independent variables in the n simulation runs
assumptions hold, then OLS is known to yield the
are denoted by the n × q matrix of independent
best linear unbiased estimators (BLUE), where "best"
variables X:
I means minimum variance, if these assumptions do
not hold, then BLUE result if GLS is u~d:
I xt, ... x,p,, (xl.k-lx,,kY
(37'
1 x:l ... xv,, (x tx22)... p
X = (6.1)
Their covariance matrix is
1 x.l ... ny ~-)-t (6.7)
~t
Observe ~at the last k(k-1)/2 columns of X follow
from the specification of the columns 2 through If ,je'f is a diagonal matrix D, then GLS can be sim-
(k + 1). The latter k columns form the so-called plified° The resulting GLS procedure is known as
design matrix; see Section 8. Using standard regres- weighted least squares, the weights for the observa-
sion notation, the familiar least-squares estimators tionyi being inversely proportional to o~; see 4.
follow: Unfortunately, GLS assumes that the covariance ma-
..¢,
trix ~y is known. Actually this covariance matrix is
~_- (.,~,. ~)-t. ,~, .~. (6.2)
unknown o~ that we have two options:
(1) Estimate the covariance mat,'L~ and substitute
In classical regression analysis it is assumed that the
this estimator ~ in eq. (6.6). As Schmidt ! 3, pp.
observationsyi are statistically independent with
71-72 ~hows this new estimator has the same
constant (homogeneous) variances 2.,.o = o .2 Under
asymptotic distribution as the regular GLS estimator
these "classical" assumptions the covariance matrix
and remains uv.biased (ander certain mild technical
of~ is
conditions). However, its small-sample proporties
"W Y)-'. )s.6(
are unknown .v
(2) Continue using OLS, even when the classical
In a simulation experiment the observationsyi
:~¢.~J Klalnen et La / Generalization of simulation retult~
~mpt~ are not met. The OLS ~t~ator of eq.
(6,2) ren'~,s unbiased.
. y = f~x)
Whether ~ ~panying tom remain valid, w~
be discussed in the next sections. eW ~ give results
for both options. In general, one might apply both
OLS and GI~, and ~ that the two approach~
yield ~ same qualitative conclusions, More specif-
ically, after we finished this study, a Monte Carlo
experiment was performed for this p~lar case.
study. The experiment showed
(i) GLS with estimated covariance matrix gives
point estimators fbr fl which have smaller variances
than OLS.
(ii) This GLS yields estimated variances which
D,X
deviate strongly from the asymptotic formula, eq.
Fig. 5. Lack-of-fit F-tests.
Hence a rule-of-thumb is: Use GLS for point esti-
mators, and OLS for tests and confidence intervals .2
Z .1 The F-test for model lack-of-fit
Observe that other criteria than Least Squares
can be used. For instance, we may wish to mimimize
Strictly speaking, this test applies only under the
the sum of the absolute deviations. Intuitively ap-
classical assumptions of regression analysis, i.e., inde.
:nilaep g is the minimization of the sum of relative
pendent observationsy~ with constant variances
errors, wh~h leads to a linear programming problem.
= o ~. Moreover, the observations should be nor-
Unforturtately, the properties of the resulting estima-
really distributed. The derivation of this test is pre-
tors are unknown, whereas for OLS and GLS we
sented in 15. The underlying idea is illustrated in
know that the estimators are BLUE, and we have a
Fig. 5. The dots in that figure show how for a specific
battery of statistica~ tests available. The criterion
x-value the y-values spread around their true mean
also affects the sensitivity of the resulting estimates
value) ~(y) -- f@). The cur!y brackets denote the
to oufliers, i.e., wild observations ony oz x. Instead
residuals, i.e., file deviations of the average response
of selecting an esthnafion algoriflu~, we can select a
at a s pec!fic xffalue, from its fitted regression model,
matrix of independent variables X that minimizes
say>y =~o +~1" x.
the sensitivity of our estimates towards outliers. For
Some authors argae that this F-test is not "very"
references on these topics we refer to vaa den Burg
sensitive to non-normality and heterogeneous varian-
et aI. 15, p.341.
ces. It still remains that the test may have small
"power", i.e., it may have a low probability of detect.
gah deviations from the nul-hypothesis, namely, the
7. Tests for model valkla~n
fitted regression model does provide good fit. In view
of these results, we decided to compute this F-statis-
In this section we shall present two statistical
tic, but we shall not base our acceptance of the fitted
tests for checking the v~dity of the metamodel, i.e.,
model exclusively on this test. A different approach
for comparing the predictions of the regression model
ll.a~ he discussed in a moment. Let us first examine
to its "true" values provided by the shnulation model.
the role of q, the number of parameters in the model.
The f'u~t test is .~tandard in experimental design, but
In general, the philosophy of science is to try and
not m regression sisylmx~ as applied by management
explain a phenomenon as simply as possible. More
scientists.
specifically, we often hypothesize that the param.
eters fl of the regression metamodel are zero. The
2 The estimates of the present stlzdy, namely #artd ~/2 sign6ffcanee of an estimated regression parameter i i
served as "true" values in the Monte Carlo experiment. As can be tested by means of the Student t-test:
X-matrices we used a 16 x I3 lct,.a 26 X 13 matrix based on
Tables ! ~d 2 of the p~ese~ srady. The number of Monte
oir_aC repliealions was 250. For the exact input data see
!5, pp.59-60, 63 ~
J*
.C.P.d Kleinen et La / noitazilareneG of simulation stluser 57
where ~ denotes the hypothesized value, usually (ii) Extra observations may also be provided
zero. Since there are a great many degrees of freedom by previous debugging runs.
(o = 128), we replace the t-statistic by the standard (iii) Moreover, we should simulate some condi-
normal variable. We can reestimate the metamodel, tions corresponding to the values that the factors
after one or more effects ~i are set zero. usually assume in practice. These values will not
eW may further test whether interactions should be deduJo~-,i in X since X represents (reasonable)
be included in the regression metamodel, tinder the extreme conditions (x i = -l, xj = + 1) as we shall
classical assumptions an F-statistic can be used to test see in the next section.
all two-factor interactionsointly. Note that the (iv) It is further wise to generate new observa-
above tests can also be applied in GLS. For more tions at the "center" of the design (x/= 0 for all
details see 15, pp.34-38. factors), in order to cheek whether pure quadratic
effects of quantitative factors are zero; see Fig. 6.
Z2. The t-test.for model validation (v') A trick to obtain "new" observations, is to
delete one observation i from the n "old" obser-
For statistical reasons, and for its appeal to model vations; compute the estimators p 0~ from the
builders and users, the following procedure is recom- remairfing (n- )1 observations~ ,~ (i' :/:/), and pre-
mended. diet Yi using these esttrr,,,~,.. " -*---~ ~ ~I i';see the discus-
(1) Generate one or more new observations.ve sion oa jackkn~mg in 5.
using the simulation model. Each response Ye is then Incorporating the "new" observafionsyg, we
compared to the predicted value-fie which is based on might estimate some new parameters, provided the
the regression metamodel estimated from the old ob- new matrix X remains non-singular. Hopefully, the
servationsYi. estimates of the new parameters turn out to be insig-
(2) If the new responsesyg agree with the model, nificant. Note that the t-test can also be used to com-
then the estimates of the model parameters 1 era pare the simulation output (~) to the real-life output
made more accurate by adding the new observations 0'). For a further discussion of statistical validation
,~,,~ era recomputing the estimators .I The procedure procedures we refer to 15.
of step (1) is traditional in the validation of models.
Actually, it may not be necessary to generate new
observationsye: 8. Selecting the experknental design matrix
(i) In specifying and realizing the matrix of
independent variables X we are liable to make In our case, as in most other situations, an exten-
errors. Indeed we made such errors in our expert- sive pilot phase preceded the phase of formal experi-
merit. Some of these "wrong" simulation runs mental design and analysis: approximate analytical
could be used in the correct specification of X; models and solutions, debugging and verification
the other runs provided the new observationsyg runs of the simulation model, ad hoe simulation runs,
for validation. general theory on inventories and queuing, common
sense. This research showed that over a relatively
large area the response variable (stacking capacity
Ely)
utilization measured by its average and quantfles) is a
linear function of the independent variables x2 and
x3 respectively, where x2 denotes the mean inter-
arrival rate of ships and x 3 denotes the mean num-
ber of days spent by containers ah the stacking area.
si~,~" linearity was found, keeping all other variables
constant. In the selection of a formal experimental
design we start by assuming that no quadratic effects
(and higher effects) are important, but leaving open
the possibility of interactions among factors. Having
specified the tentative form of the metamodel used
t --'~'X to interpret our simulation runs, we have to specify
-1 0 1¢
the values of the independent variables that occur in
.giF .6 gnitseT quadratic .stceffe
~8 .C.~J Kleilnen et .la / no~taz~iareneG o) sunulatton remits
g~ m~el In ~ present section we shall present Table 1
eh~ general theoo, of experh~ental design focussing Expe~mentsl de~ for six factors
on our s ~ situation, and in the next section we
Run 2 4 5 6 1(=56) 3(=45)
~all discass how the result~ design was applied in
the ac~val siraation.
I + ÷ + + + -!-
As we shall see in the next section (Table 3), we 2 - + + + + +
w~ to ~dy six factors, i.e. k = 6. Based on analyti- 3 + - + + + -
cal results for a simp|ified model, preceding explora- 4 - - + + + -
+ + - + - -
5
tory ~,nulafioa rare, and conunon sense, we believed
6 - + - + - -
1hat the following six mteracaons might be impor-
7 + - - ¢ - +
tant: ~t2, ~13, ~23,/~, ~2s, ~2s- Itence together with
8 - - - + - +
the six ~ effects ~! md the overall mean ,o~ we 9 + + + - - +
need o~ esttmate thh-teen parameters/3. At this point 10 - + + - - +
II + - + - -
we might ch,~eage the reader to yfice~l~ his design
12 - - + - - -
matrix, before reading the rest of this section! This
13 + + - - + -
is the traditional domain of experimental design
14 - + - - + -
theory. The results of this theory have been reported 15 + - - - + +
in numerous articles, textbooks, and tables. An exam- 16 . . . . + +
ple of a recent textbook is 3. Results of irranedhate
relevance to simulation experiments are summarized
in 5,6. The selection of experiments can be based meet a number of requirements that may be formu.
on the following approaches. hated for experimental designs. First we shall present
(1) Common sense, wha~ever that may mean: A the design ew actu~y used in our simulation experi.
popular belief is that a scientific experiment requires ment.
that all facto:s except one constant
should be kept Based on 5, pp.320-372 we propose the design
when proceeding frcm one observation to the next specified in Table .1 In ~ table the two levels of a
one: cewris paribus or one-factor-at-a-tkme method. factor are denoted by + and -, shorthand for +1 and
"factorial" experiments, howeveL ~ levels cf a -1. eW start by writing down all 24 combinations of
factor are combined with all levels of each other fac- the factors 2, 4, 5, ,ra d 6. The levels for factor 1 are
tor. Actu~y, factorial experiments eza more efficient specified by pairwise multiplication of the elements
since they yield more reiable estimator3 of the maku in the columns of the factors 5 and 6. For factor 3 the
effects; moreover, such experiments can ~ve esti- columns of the factors 4 and 5 are used. The specifi-
mates of the factor interactions; ees 5, pp.289-290. cation of the elements ofxl through x6, i e. the
Unfortunately, the ~umber of combinations grows design matrix, determines the remaining independent
dramatically as the number of factors increases. In variables x tx2~ ..., x2x6. All thirteen cdunms are
our case, even ff we examine only two levels per fac- orthogonal. The reader can fmther check that, for
:or, we have 26 = 64 combimfions. Sinr.e we con- instmce, the column for x4 is identical to that for
jeeture t~tat only thirteen parameters are important, ~e interaction x3x s but this interaction was assumed
we are satisfied with much less then sixty-four obser- to be zero from the start .3
vations. eW might select the factor combinations Let us consider the advantages and disadvantages
using common sense, or foUowing one of the next
approaches.
(2) Ad-hoe optimization ofdesCgn: eW may use ~: For readers familiar with basic experimental design metho-
the computer to fred an "optimal" design, say, a dology we mention tha~ the generators 1 = 56 and 3 = 45
lead to the defining relation I = 156 = 345 = 1346. Hence
design ~ g the Mean Squared Error ~ISE) of
~he aline pattern is 1 = 56, 3 = 45, 4 = 35, 5 = 16 = 34,
the regression estimators where ESM equals vadzuce
6 = 15, ignoring interactions among more than two factors.
plus squared bias; see 12.
Selecting three-factor interactions as generators gives less
(3) Standard experimenta: :s.ngised Over die past, desirable confounding in our particular case. For instance,
say, fifty years standard designs have been derived selecting a 26~72 design with generators 5 = 123 and
hcua~ fractional factorial defigns. At e~ff end 6 = ~24 results hal= i235 = 1246 = 3456 so that, for
as 2 ~-p
instance, 13 = 25, i.e., two possibly important interactions
of this section we shall ees how well these desi~
ate confounded.