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.