“What-if” analysis in computer simulation models: A comparative survey with some extensions

Math1 Cornput. Moaklling, Printed in Great Britain

Vol. 14, pp. 101-106,

1990

08957177/90 $3.00 + 0.00 Pergamon Press plc

“WHAT-IF” ANALYSIS IN COMPUTER A COMPARATIVE

H. Arsham and A B. Kahn School of Business, University

Abstract.

SIMULATION MODELS: WITH SOME EXTENSIONS

SURVEY

of Baltimore,

MD 21201 USA

Baltimore,

The simulation models are often subject to errors caused by the estimated parameter(s)

underlying input distribution function. to small changes in the parameters

of

“What-if” analysis is needed to establish confidence with respect

of the input distributions.

However traditional

“what3

analysis

requires a separate simulation run for each input value. Four

different

methods

extrapolation/interpolation

for

estimating

are presented.

performance By simulating

function

for

several

a simple reliability

scenarios

using

model a comparative

experimental study on the efficiency of the these methods is provided.

Kevwords. Perturbation

analysis; Model validity; Score function, Derivative estimation; What-if analysis;

Output analysis; Non-linear variance reduction technique; Extrapolation/interpolation

by simulation.

INTRODUCI-ION The simulationist Consider a stochastic system with a parameter v. Let

must meet managerial demands to consider

Model validation

alternatives.

requires

us to cope with

uncertainty in the estimation of v. Adaptation

I

J(v) = E PQI =

LOI)W;Vy

be an expected performance vector with known probability

new environments

requires

motivation

from

comes

an

adjustment

development

of a model to in v. Further

of

techniques

for

optimization, response surface mapping, problem inversion and other applications. measure where Y is a random density function

An obvious solution to the “what if’ problem is the Crude

(pdf), f(y;v)

Monte Carlo (CMC) method, which estimates J(v+ sv) for each

depends on v, and L is the performance measure.

6v separately by rerunning the system for each vf 6~. The costs

We resort to Monte Carlo simulation when L is either unknown

in CPU

or too complicated to calculate analytically.

simulations and increase when high accuracy requires increases

Estimation of J(v)

time

can

be prohibitive

even

for deterministic

is obtained by averaging over n independent replications L(yi),

the number of replications of a dynamic stochastic simulation.

i=1,2,...,n;

Recent

information

^ J(v) =

developments obtained

neighborhood.

; I&). i=l

(2)

in from

“What-iP’ analysis a single

run

Furthermore, we also consider the use of results

We refer to the former as perturbation

motivations

and the latter as

interpolation.

law of large numbers. are strong

the

from runs at two or more points over the intervening in&rval.

This is an unbiased estimator for J(v) and converges to J(v) by

There

extend

to the closed

This paper summarizes and compares four general and three for estimating

the cxpectcd

specialized recent methods, some of which have been published

performance measure J(v) for a small change in v to V+ 6v, that

only as working papers.

is to solve the so-called “what if’ problem.

The General Methods are:

Likelihood Ratio, Linear Taylor

Approximation,

and Interpolation

runs). 101

Exponential,

(IWO or more

Proc. 7th Int. Conf. on Mathematical and Computer Modelling

102

The Soeeial Methods are: Simterpulatio” Perturbation

Analysis, and All bpe

.

for Markovian Queues, Parameter

Method for

where J(v+m)

= I%@

Semi Markov Model. Since 0 is unknown , we can estimate 2 by

LIKELIHOOD

^ _s S* = S,, - 2J(v+sv). St* + J(v+sv)S,

RATIO (LR)

(6)

where Arsham et al. (1987)

proposed

a model based on Radon-

Nikodym theorem to estimate J(v+ bv) for stochastic systems in a single run.

S,,= & Wi-Li&“-l), i==l and

sz=

~,-&“-I), i=l

S12= ““(I+ WI -Lw)(W, - IF) /(n-l) i=l J(v+ bv) =

L (y).f(y;v+sv)dy I

=

which could be used to construct the desired confidence intervals

f(y;v+6v) f(y;v)dy L(y) ___ f&v)

for J(v+sv).

I

= E { L](y;v)].W )

(3) EXPONENTIAL

where the likelihood

ratio W = f(y;v+sv)/f(y;v)

adjusts the

sample path, provided fw,v) + 0. Notice that by this change of

In the statistical literature the efftcient score function is defined

probability space, we are using the common realization as J(v).

to be the gradient

The generated random vector y is roughly representative with f(v,).

Each of these random

hypothetically came from f(v+6v).

observations,

of Y, d S(y) = - 1” ffy;v) dv

could also

W weights the observations

(7)

according to this phenomena. We consider the following exponential The “What-if” estimate based on (3) is

^

(approximation)

for J(v+6v) in a first derivative neighborhood

n

J(v+ 6v) = E { L(y) e(@s(vH&)}

J(v+ 6v) = C I+. WJn i=l

(4

model

of v

(8)

where ek6.j z E { e@+%))

(9)

6v = v-v@

(10)

which is based on only one sample path of the system with parameter

v and the simulation for the system with v+ 6v is

and

not required. Unfortunately CMC.

LR produces a larger variance compared with

However, the following variance reduction

Equation (8) can be written as:

techniques

(VRT) may improve the estimate. n

n

J(v+ sv) =

(11)

C Li Wi / C Wi i=l i=l

By virtue of the central limit theorem one can show that:

to estimate J using (11) based on n independent replications we (5)

have

J(v+ 6~) = ;A i=l

/ ;B; f=l

(12)

ttlh (J(v+sv) - j(v+,v)) =r

N(O,l) as “a

where A= L(y&exp[@v)S(yi)]

(13)

and B= exp[(6v)S(y,)].

PJ)

Proc. 7th Int. Co&

Denote

A=

zA&,

on Mathematical

g = EBJn,

where S=fO;v)/f@,v)

then by a similar argument given in the previous section, we

103

and Computer Modelling

estimated

derivative

is the Efficient Score Function. could be incorporated

This

in Optimization

[Arsham (1988a), L’Ecuyer et al. (1989), Rubinstein (1989)] and have

Inverse Problem [Arsham, (1988b). (1989a)].

nH (J(v+ 6v) - J(v+ 6v)) => N&l)

as na

(15)

(21)

J’(v0) =1/n lLb].Sh] i=l where

where J = G

evaluated

(16)

at v=vh

and

Y, are

the same

n independent

replications used in (2) to estimate the nominal performance J(Q).

where

The estimated derivative based on (21) converges to the true

n S,, = C { A, - A l2 /(n-l) i=l

Sn =

(22)

SfkJ = ~“cv,,V) I fc_vJ)

and

derivative in the sense of mean square error (Arsham et al. [1989]).

c”{B, - El2 /(n-l>

TWO-POINT BASED INTERPOLATION

(‘7)

i=l SIMUTERP, proposed by Polster (1988), combines interpolation with simulation.

Given two points, v1 and v2 (scalars only)

sufficiently close, Simuterp simulates at the two points then interpolates It can be shown that for the exponential

distribution

including Markovian processes this approach

family

is equivalent to

Simuterp

Criterion”

to maintain

Although a valid

confidence interval, it is difficult to apply for the general case. Arsham

the Likelihood Ratio, Arsham (198913).

for any desired points in between. has a “Spacing

(19894) approach

assumes the given v1 and vz are

sufficiently close and looks for the “best” linear interpolation

in

the sense of minimum error on the interval. Clearly, LINEAR TAYLOR APPROXIMATION

I

J(h) = The technique proposed in A&am feasibility

of

using

Taylor’s

L (y).f(y;v&V=

et al (1989) is based on the expansion

of

J(v)

in

the

neighborhood of v,, J(v) = J(vo) + sv J’(vo) +.... where prime (‘) denotes derivative wrt v and dv=v-v,.

(18) We will

restrict our model to the first order expansion

I

e LCy).f(y;vo)dy+

^

I I

JL(v) = J(v,,) + bv J’(vO)

(19)

(l-~)WWWoY’y=

f(Y%)

4 L(Y)

I

I

f(wddy +

~

fW1)

f(Y;vo)

(l-4 L(Y)___

with variance

^

^

fWi)dy.

f(Y?zJ

One obvious choice is )= f(j’;vl)/[f(y;vl)+f(y;vz)]. This method

Var[JL(v)J= Var[J(vo)] + (bv)’ Var[J’(v,)] f ^ ^ m ~~[J(~o)~J’(~o)l

can easily extended to k-point interpolation.

where

For 2-point interpolation, then the linear interpolated

J’(v)=

L(y).r(y;v)dy=E[L(y).S] I

if we let .$ to be constant 0 i 0 d 1, “what-if” estimated value is:

(20) J(w9

= #J&o) + (1-e) J&)

(24)

104

where

Proc. 7th Int. Conf on Mathematical and Computer Model&?

the hvo estimates

on

the RHS

of (24) are two

indeuendent Likelihood Ratio extrapolations

using the two end

EXPERIMENT

Consider a coherent reliability shown in Fig. 1. Assume all

points. We define e* as the e in this convex combination minimum error in the estimate. 2

^

with the

-L--r>

That is, to minimize 2

^

varp(v,~)l=ovar[J,(v~)]+(l-0)Var[J,(v~)]

(25)

Fig. 1. A reliability network

By the first order necessary and sufficient conditions, the optimal components work according to the following Gamma density:

e is:

^ ~~~~J&)l (26) We will apply the independent replication method to evaluate the life performance J(v). Application of our approach to other Thus the “best linear” interpolation

for any point in interval

methods e.g. batch means, is quite similar. Clearly, the density function for the system is

[vr, vz] is: ^

1

^

f(wx2so) =

J(Q=)=i J&o) + (1-d) J&s)

(27)

which is the optimal interpolation variance.

in the sense of minimum

Our numerical results are very encouraging therefore,

~

vs

x,~~,exp[-(x,+~,+x,)/v]

The Likelihood Ratio Function W and Score Function S at the nominal v=vO are

we believe this approach deserves more work.

W= [ vd(v~+~)16exp[(xl+x,+x,)6v~~(v,+6v)] OTHER RELATED MODELS

and

A literature search reveals three related models which deal with special cases or require special conditions.

We shall comment

very briefly on each.

Since the analytic solution for this system is not available, we have used a Crude Monte Carlo simulation with n=60,000 runs for each point v,,=O.75(0.05)1.2.5. Our experiments are performed on a VAX 11/780 computer using streams 1 through 3 (and 4

Simternolation

proposed

by Reiman

et al. (19S7) provides

“What-if” analysis for a steady state performance

measure, but

limited to a class of open Markovian queueing systems.

through 6 for interpolation)

SIMSCRIPT II.5 to generate XI, X1,

and X,. Results are presented in Table 0. Based on these ‘true’ values we considered the Relative Absolute Error (RAE) as a measure of efficiency

Finite Perturbation

Analvsis Technicme developed by Ho, et al.

RAE=

ITrue-Estimatelnrue

(1983) uses the same random numbers at a different parameter

in comparing different methods.

setting. This approach

As a second measure for our comparison

is limited to cases where the inverse

distribution functions are available in closed form.

Coefficient of variation i.e. C.V.= lOO(Standard Error)/Estimate.

The General (Parameter1 Method presented by Cassandras and Strickland (1989) permits perturbation

of any type of system

parameters (m just probability distribution parameters as in all other proposed

techniques)

in estimating

probabilities in Semi-Markov Processes.

of stationary

state

we also considered

Proc. 7th Int. Conf: on Mathematical and Computer Modelling

Table 0: Average system function of v, n=60,000 V

J(V)

Sam

Var

v

life

time

J(v)

Sam

as

a

Table 2: Exponential nominal v,=l, n=1000

2.2969 2.1821 2.0672 1.9524 1.8375 1.7227

Table 1: Extrapolation

V

1.00 0.95 0.90 0.85 0.80 0.75

1.737 1.570 1.410 1.255 1.112 0.977

1.00 1.05 1.10 1.15 1.20 1.25

2.2969 2.4118 2.5266 2.6415 2.7563 2.8712

Likelihood Ratio with nominal v*=l,

J& 2.3463 2.2215 2.0946 1.9669 1.8392 1.7122

C.V.%

RAE%

1.76 1.45 1.24 1.20 1.34 1.63

2.15 1.80 1.33 0.74 0.09 0.61

1.737 1.915 2.101 2.297 2.501 2.714

J(v)

C.V.%

RAE%

1.00 0.95 0.90 0.85 0.80 0.75

2.3463 2.2406 2.1436 2.0545 1.9726 1.8971

1.76 1.70 1.67 1.67 1.70 1.75

2.15 2.68 3.70 5.23 7.35 10.10

1.00 1.05 1.10 1.15 1.20 1.25

2.3463 2.4613 2.5861 2.7209 2.8655 3.0193

1.76 1.87 2.03 2.23 2.49 2.79

2.15 2.05 2.35 3.01 3.96 5.16

(LHD) n=lOOO

Table 3: Linear n=lOOO

2.3463 2.4679 2.5848 2.6959 2.7997 2.8952

1.76 2.14 2.56 3.02 3.49 3.98

V

J&

C.V.%

RAE%

1.00 0.95 0.90 0.85 0.80 0.75

2.3463 2.2353 2.1231 2.0103 1.8971 1.7839

1.76 1.69 1.67 1.68 1.75 1.87

2.15 2.44 2.71 2.97 3.24 3.55

1.00 1.05 1.10 1.15 1.20 1.25

2.3463 2.4556 2.5627 2.6670 2.7680 2.8655

1.76 1.86 1.99 2.15 2.31 2.49

2.15 1.82 1.43 0.96 0.43 0.20

With

with

vo=l,

h

2.15 2.33 2.31 2.06 1.56 0.84

Table l-l: Likelihood extrapolation with VRT (LHD/VR),

Extrapolation

J(v)

C.V.%

RAE%

1.00 0.95 0.90 0.85 0.80 0.75

2.3463 2.2361 2.1258 2.0156 1.9053 1.7951

1.76 1.88 2.07 2.34 2.68 3.11

2.15 2.47 2.83 3.23 3.68 4.20

1.00 1.05 1.10 1.15 1.20 1.25

2.3463 2.4566 2.5668 2.6770 2.7873 2.8975

1.76 1.71 1.71 1.75 1.83 1.91

2.15 1.85 1.59 1.34 1.15 0.91

V

1.00 1.05 1.10 1.15 1.20 1.25

Extrapolation

Var V

1.00 0.95 0.90 0.85 0.80 0.75

105

ratio n=lOOO

Table 4: Two-point v,=O.75, v2=1.25, n=lOOO

V

Jh

interpolation (each)

C.V.%

RAE%

1.00 0.95 0.90 0.85 0.80

2.2480 2.1373 2.0319 1.9327 1.8364

1.30 1.45 1.62 1.73 1.70

2.13 2.06 1.71 1.01 0.09

1.00 1.05 1.10 1.15 1.20

2.2480 2.3622 2.4780 2.5937 2.7073

1.30 1.21 1.22 1.33 1.51

2.13 1.92 1.92 1.81 1.78

Proc. 7th Int. Conf. on Mathematical and Computer Modelling

106

CONCLUSIONS

AND FUTURE RESEARCH

As always, since these experiments numerical

example,

generalizations,

one should

were done on a specific

be careful

in making

any

Comparison of the single run simulation results

suggest that even up to a 20% change in v the coefficient of variations are reasonably small, Although interpolation has limitation encouraging.

that

v must

be a scalar,

In most practical

method

the results

cases, the estimates

Arsham, H. (1989c). Sensitivity and optimization of computer simulation models. Proceedines of Modeline and Simulation, 19, 11351142.

are

having

coefficient of variation less than 5% are generally acceptable. We noted that the efficiency in terms of C.V. and/or RAE not a single method is always superior. However, for our numerical example we notice that:

Arsham, H. (1989d). “What-if” analysis in system performance modeling: Simulation-based extrapolation and interpolation, Manuscript. Universitv of Baltimore. Arsham, H. (1988a). The what-if problem in simulation of stochastic systems. Proceedines of Southeast TIMS, 306-308. Arsham, H. (1988b). The inverse problem in dynamic system simulation. Manuscript. Universitv of Baltimore. Arsham, H., Feuerverger A. McLeish D.L. Kreimer, J. and R. Y. Rubinstein (1989). Sensitivity Analysis and the ‘What It’ problem in Simulation Analysis. Mathematical and Comuuter Modelling, 12, 193-219.

while increasing Y these methods can be ordered as: Interpolation,

Linear, LHD/VR,

Exponential,

and LHD,

in

terms of C.V. LHD/VR, Linear, Interpolation,LHD,

Exponential, in terms of

RAE. While decreasing v we have LHD=LHD/VR=Exponential=Interpolation,

Linear, in terms of

C.V. and LHD,Interpolation,Linear,LHD/VR,

Exponential,

in terms of

RAE..

Items for further research include:

0

to incorporate the result of this study in a random search optimization

ii)

technique;

to introduce efficient variance reduction and bias reduction techniques with a view to improving the accuracy of the existing and the proposed methods;

iii)

to extend our methodology to higher order approximation

iv)

to adopt these methods for larger and more complex

Arsham, H., Xreimer J. and R. Y. Rubinstein (1987). Application of Radon-Nikodym theorem for simulation of queueing model. Proceedines of the Euronean Simulation Multiconference, 95-99. Cassandras, C. G., and S. G. Strickland,(l989). A general approach for sensitivity analysis of discrete event systems with Markov properties, IEEE Transactions on Automatic m. (To Appear). L’Ecuyer, P., Giroux, N. and P. W. Glynn (1989). Stochastic optimization by simulation: Some experiments with a simple steady-state queue. Technical Report No. 13, Stanford University. Glynn, P. (1987). Likelihood ratio gradient estimation: An overview. Proceedines of the 1987 Winter Simulation 366-375. Ho, Y. C., Cao, X. and C. Cassandras, (1983). Infinitesimal and finite perturbation analysis for queueing networks. Automatica. 19, 439-445. Reiman, MI., Simon, B., and J. S. Willie, (19S7) Simterpolations: Estimating an entire queueing function from a single sample path. Proceedines of the 1987 Winter Simulation, 358363.

for linear and exponential models;

system.

Polster, R. S. (1988). A technique for interpolating system response to parametric perturbations in simulation. D. SC. Dissertation. The George Washincton Universirv. Rubinstein, R. Y., (1986). Monte Carlo Ontimization. Simulation and Sensitivltv of Oueueine Networks, John Wiley, New York.

Acknowledements-

The authors appreciate the helps from Mr.

Jae Bong Yu. This work received support from the University of Baltimore Educational

Rubinstein, R. Y. (1988). How to optimize discrete event system from a single sample path by the score function method. working paper, Technion, Israel.

Foundation. Rubinstein, R. Y. (1989). Sensitivity Analysis and Performance Extrapolation for Computer Simulation Models. Onerations Research 31, 72-81. -I REFERENCES

Arsham, H. (1989a). On the inverse problem in Monte-Carlo experiments. Inverse Problems, In Press. Arsham, H. (1989b). Exponential Tangential Extrapolation of computer performance based on the efficient score function, Manuscript. Universitv of Baltimore.