Industrial Engineering Department
Simulation Basics
Week 2
Modeling&Simulation-IW-2007
Industrial Engineering Department
Types of Simulation
Static or dynamics
Stochastic or deterministic
Discrete event or continuous
Modeling&Simulation-IW-2007
Industrial Engineering Department
Static versus Dynamic Simulation
A static simulation is one that is not based on time.
It often involves drawing random samples to generate
a statistical outcome, so it is sometimes called Monte
Carlo simulation.
Dynamic simulation includes the passage of time.
It looks at state changes as they occur over time.
A clock mechanism moves forward in time and state
variables are updated as time advances.
Modeling&Simulation-IW-2007
Industrial Engineering Department
Stochastic versus Deterministic
Simulation
Random
input
(a) deterministic simulation
Random
out-put
(b) stochastic simulation
Modeling&Simulation-IW-2007
Industrial Engineering Department
Simulation in which one or more input variables are random
are referred to as stochastic or probabilistic
A stochastic simulation produces output that is itself random.
Simulation having no input components that are random are
said to be deterministic.
A deterministic simulation will always produce the exact same
outcome no matter how many times it is run.
In stochastic simulation, several randomized runs or
replication must be made to get an accurate estimate because
each run varies statistically.
Performance estimates for stochastic simulations are obtained
by calculating the average value of performance metric
accross all of the replication.
Modeling&Simulation-IW-2007
Industrial Engineering Department
Random Behavior
Natural behavior
Values vary within a given range and
according to a probability distribution.
Probability distributions are useful for
predicting the values, where these values are
random variables.
Depend on the type of distribution and their
parameters.
Modeling&Simulation-IW-2007
Industrial Engineering Department
Probability distributions may be either discrete
or continuous.
A discrete distribution represents a finite or
countable number of possible values.
The number of items in a lot
A continuous distribution represents a
continuum of values
Processing time
Modeling&Simulation-IW-2007
Industrial Engineering Department
An example of a discrete probability function
Modeling&Simulation-IW-2007
Industrial Engineering Department
An example of a continuous probability function
Modeling&Simulation-IW-2007
Industrial Engineering Department
Simulating Random Behavior
One of the most powerful features of
simulation
The ability to mimic random behavior.
Simulating random behavior requires that a
method be provided to generate random
numbers as well as for generating random
variates based on a given probability
distribution.
Modeling&Simulation-IW-2007
10
Industrial Engineering Department
Generating Random Numbers
Random behavior is imitated in simulation by
using a random number generator.
Random number generator is responsible for
producing this stream of independent and
uniformly distributed numbers.
Modeling&Simulation-IW-2007
11
Industrial Engineering Department
The U(0, 1) distribution of a random
number generator
Modeling&Simulation-IW-2007
12
Industrial Engineering Department
The numbers produced by a random number
generator are not random in the truest sense.
The generator can produce the same sequence of
numbers again and again, which is not indicative
of random behavior.
They are often referred to as pseudo-random
number generators.
Good pseudo-random number generators can
pump out long sequences of numbers that pass
statistical tests for randomness (the numbers
are independent and uniformly distributed).
Modeling&Simulation-IW-2007
13
Industrial Engineering Department
Linear Congruential Generators (LCG)
The most common method for generating
random numbers
A sequence of integers Z1, Z2, Z3, ... is defined
by the recursive formula:
Z i aZ i 1 c mod m
a : the multiplier;
m : the modulus;
c : the increment;
Z0 : the seed or starting
value;
a,
c,
m
>
0
and
integer
Modeling&Simulation-IW-2007
14
Industrial Engineering Department
The Zi values are bounded by 0 Zi m 1
and are uniformly distributed in the discrete
case.
The continuous version of the uniform
distribution with values ranging between 0 and
1 can be obtained by
Zi
U i , i 1,2,3,
m
Modeling&Simulation-IW-2007
15
Industrial Engineering Department
Example of LCG
Modeling&Simulation-IW-2007
16
Industrial Engineering Department
The maximum cycle length that an LCG can
achieve is m.
To realize maximum cycle length, select the
values of a, c, and m carefully (Pritsker 1995)
m = 2b where b is determined based on the
number of bits per word on the computer being
used (for computer with 32 bits, b = 31)
c and m such that their greatest common factor is
1.
a = 1 + 4k, where k is an integer.
Modeling&Simulation-IW-2007
17
Industrial Engineering Department
Frequently, the long sequence of random number is
subdivided into smaller segments referred to as
streams.
To subdivide the generators sequence of random
numbers into streams, it is need:
to decide how many random numbers to place in each
stream
to generate the entire sequence of random numbers (cycle
length) produced by the generator and record the Zi values
that mark the beginning of each stream.
Each stream has its own starting or seed value.
Modeling&Simulation-IW-2007
18
Industrial Engineering Department
There are two types of LCG:
Mixed LCG:
Mixed LCG
Multiplicative LCG
c>0
Multiplicative LCG
c=0
Modeling&Simulation-IW-2007
19
Industrial Engineering Department
Other kinds of random generators
(see Kelton pp.431-434)
More General Congruences
Composite Generators
Tauswarthe and Related Generators
Modeling&Simulation-IW-2007
20
Industrial Engineering Department
Testing Random Number Generators
The numbers produced by the random number
generator must satisfy two properties:
Independent
Uniformly distributed between zero and one
Generate a sequence of random numbers:
U1, U2, U3, ... .
Modeling&Simulation-IW-2007
21
Industrial Engineering Department
The hypothesis for testing the independence
property:
H0: Ui values from the generator are independent
H1: Ui values from the generator are not independent
The most common statistical method
the run test
the runs above and below the median test
the runs up and down test
Modeling&Simulation-IW-2007
22
Industrial Engineering Department
The hypothesis for testing the uniformly
property:
H0: Ui values are U(0, 1)
H1: Ui values are not U(0, 1)
The most common statistical methods:
The Kolmogorov-Smirnov test
The chi-square test
Modeling&Simulation-IW-2007
23
Industrial Engineering Department
Generating Random Variates
Methods for generating random variates (see
Kelton Chap.8)
Inverse transformation
Composition
Convolution
Acceptance-Rejection
Methods employing special properties
Modeling&Simulation-IW-2007
24
Industrial Engineering Department
Inverse Transformation Method
Continuous distribution
For a given probability density function f(x), find
the cumulative distribution function of X, i.e., F(x)
= P(X x).
Set U = F(x), where U is uniform (0, 1) and solve
for x.
Solving for x yields x = F-1(U).
The equation x = F-1(U) transforms U into a value
for x that conforms to the given distribution f(x).
Modeling&Simulation-IW-2007
25
Industrial Engineering Department
Inverse transformation for continuous distribution
Modeling&Simulation-IW-2007
26
Industrial Engineering Department
Generating variates from the exponential
distribution with mean .
The probability density function f(x):
1 x /
, x0
e
f x
0,
elsewhere
The cumulative distribution function F(x):
1 e x / , x 0
F x
elsewhere
0,
Modeling&Simulation-IW-2007
27
Industrial Engineering Department
Setting U = F(x) and solving for x yields:
U 1 e x /
e x / 1 U
ln e x / ln 1 U
x ln 1 U
x ln 1 U
The random variate x is exponentially distributed
with mean .
Modeling&Simulation-IW-2007
28
Industrial Engineering Department
U1 0.27 x1 2 ln1 0.27 0.63
U 2 0.89 x2 2 ln 1 0.89 4.41
Modeling&Simulation-IW-2007
29
Industrial Engineering Department
Discrete distribution
For a given probability mass function p(x), find
the cumulative distribution function of X, i.e.,
F(x) = P(X x).
It is assumed that X take on only the values x1,
x2, ... where x1 < x2 < ... .
Algorithm:
Generate U ~ U(0, 1)
Determine the smallest positive integer I such that U
F(xI), and return X = XI.
Modeling&Simulation-IW-2007
30
Industrial Engineering Department
Inverse transformation for discrete distribution
Modeling&Simulation-IW-2007
31
Industrial Engineering Department
Generating variates from the following
probability mass function (arbitrary discrete
distribution):
0.10, x 1
p x P X x 0.30, x 2
0.60, x 3
Modeling&Simulation-IW-2007
32
Industrial Engineering Department
The cumulative distribution function F(x)
U1 = 0.27, because 0.10 < U1 0.40 then x1 = 2
Modeling&Simulation-IW-2007
33
Industrial Engineering Department
Algorithms for generating variates from other
distributions, see Law and Kelton, Simulation
Modeling and Analysis, Chapter 8.
Modeling&Simulation-IW-2007
34
Industrial Engineering Department
Automatic teller machine (ATM) system
What type of data required ?
What performance measurement to be analyzed ?
Modeling&Simulation-IW-2007
35
Industrial Engineering Department
TIME CALCULATION
ArrivalTime (i) = ArrivalTime (i 1) + InterarrivalTime (i)
If ArrivalTime (i) < DepartureTime (i 1) then
BeginServiceTime (i) = DepartureTime (i 1)
else
BeginServiceTime (i) = Arrival (i)
DepartureTime (i) = BeginServiceTime (i ) + ServiceTime (i)
TimeInQueue (i) = BeginServiceTime (i) ArrivalTime (i)
TimeInSystem (i) = DepartureTime (i) ArrivalTime (i)
Modeling&Simulation-IW-2007
36
Industrial Engineering Department
A Spreadsheet Simulation
A single-server queuing system (i.e., ATM
system)
Queuing model : (M/M/1): (FIFO//)
Inter-arrival time ~ exponential distribution with
mean = 3.0 minutes
Service time ~ exponential distribution with mean
= 2.4 minutes.
Modeling&Simulation-IW-2007
37
Industrial Engineering Department
Inter-arrival times (X1i):
Z10 = 3
Z1i = (21 * Z1i-1 + 3) mod (128)
U1i = Z1i / 128
X1i = -3.0 ln (1 U1i)
Service times (X2i):
Z20 = Z125
Z2i = (21 * Z2i-1 + 3) mod (128)
U2i = Z2i / 128
X2i = -2.4 ln (1 U2i)
Modeling&Simulation-IW-2007
38
Industrial Engineering Department
Modeling&Simulation-IW-2007
39
Industrial Engineering Department
Simulation Replications and
Output Analysis
The spreadsheet simulation of the first 25
customers to arrive to the ATM system
The simulation output gives:
average time in queue = 1.94 minute
average time in system = 4.26 minute
These results represent only one possible value
of each performance measure.
Modeling&Simulation-IW-2007
40
Industrial Engineering Department
The inter-arrival times of customers and their
service times are random variables, then the
output of the simulation model is also random.
The 2nd replication change the seed values
Z10 = 29
Z20 = Z125 = 92
Modeling&Simulation-IW-2007
41
Industrial Engineering Department
Modeling&Simulation-IW-2007
42
Industrial Engineering Department
Summary of ATM system simulation output (for
two replications)
Modeling&Simulation-IW-2007
43
Industrial Engineering Department
Generating Random Variates
Modeling&Simulation-IW-2007
44
Industrial Engineering Department
Uniform
1.
Generate U ~ U(0,1)
Return X = a + (b a)U
Modeling&Simulation-IW-2007
45
Industrial Engineering Department
Uniform
The distribution function of a U(a, b) random
variable is easily inverted by solving u = F(x)
for x to obtain, for 0 u 1,
Algorithm:
1.
2.
Generate U ~ U(0, 1)
Return X = a + (b a)U
Modeling&Simulation-IW-2007
46
Industrial Engineering Department
Exponential
Inverse-transform algorithm:
1.
Generate U ~ U(0, 1)
2.
Return X = ln (1-U)
Modeling&Simulation-IW-2007
47
Industrial Engineering Department
m-Erlang
If X is an m-Erlang random variable with mean
, it can be written X = Y1 + Y2 + ... + Ym,
where Yis are IID exponential random
variable, each with mean /m.
Algorithm:
1.
Generate U1, U2, ..., Um as IID U(0, 1)
2.
Return
X
ln
m
U
i
i 1
Modeling&Simulation-IW-2007
48
Industrial Engineering Department
Weibull
X~
Algorithm:
1.
2.
Generate U ~ U(0, 1)
Return
X lnU
1/
Modeling&Simulation-IW-2007
49
Industrial Engineering Department
Normal
Given X ~ N(0, 1), X ~ N(, ) can be obtained
by setting X = + X.
Algorithm:
1.
2.
Generate U1 and U2 as IID U(0, 1), let Vi = 2Ui 1
for i = 1, 2; and let W = V12 + V22.
If W > 1, go back to step 1. Otherwise, let
Y 2 ln W W , X1 = V1Y and X2 = V2Y.
Then X1 and X2 are IID N(0, 1) random variables.
Modeling&Simulation-IW-2007
50
Industrial Engineering Department
Lognormal
If Y ~ N(, 2) then eY ~ LN(, 2).
Algorithm:
1.
Generate Y ~ N(, 2).
2.
Return X ~ eY.
Modeling&Simulation-IW-2007
51
Industrial Engineering Department
Generating Discrete Random Variates
Bernoulli
Discrete uniform
Binomial
Geometric
Negative binomial
Poisson
Modeling&Simulation-IW-2007
52
Industrial Engineering Department
Binomial
Algorithm:
1.
Generate Y1, Y2, ..., Yt as IID Bernoulli(p)
random variates.
2.
Return X = Y1 + Y2 + ...,+ Yt
Modeling&Simulation-IW-2007
53
Industrial Engineering Department
Bernoulli
Algorithm:
Generate U ~ U(0, 1).
If U p, return X = 1. Otherwise, return
X = 0.
Modeling&Simulation-IW-2007
54
Industrial Engineering Department
Discrete Uniform
Algorithm:
Generate U ~ U(0, 1).
Return X = i + [(j i + 1)U]
Modeling&Simulation-IW-2007
55
Industrial Engineering Department
Geometric
Algorithm:
1.
Generate U ~ U(0, 1)
2.
Return X = [ln U/ln(1 p)]
Modeling&Simulation-IW-2007
56
Industrial Engineering Department
Negative Binomial
Algorithm:
1.
Generate Y1, Y2, ..., Ys as IID geom(p) random
variates.
2.
Return X = Y1 + Y2 + ...,+ Ys
Modeling&Simulation-IW-2007
57
Industrial Engineering Department
Poisson
Algorithm:
1.
Let a = e, b = 1, and i = 0.
2.
Generate Ui+1 ~ U(0, 1) and replace b by
bUi+1. If b < a, return X = i. Otherwise, go to
step 3.
3.
Replace i by i + 1 and go back step 2.
Modeling&Simulation-IW-2007
58
Industrial Engineering Department
Group Assignment#2
Develop Simulation by using ms-Excell
2 ATM machines (identical)
1 Queuing, FIFO
50 customers
Inter-arrival time ~ exponential distribution with mean = 2.0
minutes
Service time ~ exponential distribution with mean = 2.4
minutes.
Analyze the simulation result
Report should be submitted Patrisya: Next week from
(Thursday) hardcopy & softcopy
Modeling&Simulation-IW-2007
59