Simulation Modeling

1
Simulation Modeling

• J. M. Akinpelu

2
What is Simulation?

• Simulation is the use of a computer to evaluate a
  system model numerically, in order to estimate
  the desired true characteristics of the system.
• Simulation is useful when a real-world system is
  too complex to allow realistic models to be
  evaluated analytically.

3
Systems

• System a collection of entities, e. g., people
  or machines, that act and interact together to
  accomplish some end
• System state the collection of variables
  necessary to describe the system at a particular
  time
• Types of systems
• discrete
• continuous

4
Types of Simulations

• discrete vs. continuous
• deterministic vs. stochastic
• static vs. dynamic

5
Advantages of Simulation

• Provides a way to study complex, real-world
  systems that cannot be accurately described by a
  mathematical model that can be evaluated
  analytically.
• Allows estimation of the performance of an
  existing system under some projected set of
  operating conditions.
• Allows comparison of alternative proposed system
  designs to see which one best meets a specified
  requirement.
• Allows study of a system with a long time frame
  in compressed time, or alternatively, study of
  the detailed workings of a system in expanded
  time.

6
Disadvantages of Simulation

• Each run of a stochastic simulation model
  produces only estimates of a models true
  characteristics for a particular set of input
  parameters. Thus several independent runs are
  required for each set of input parameters to be
  studied.
• Simulations models are often expensive and
  time-consuming to develop.
• The large volume of numbers produced by a
  simulation study often creates a tendency to
  place greater confidence in a studys results
  than is justified.

7
Discrete-Event Simulation

• Discrete Event System
• a system whose state changes at discrete points
  in time due to the occurrence of asynchronous
  events
• Example M/M/1 queueing system
• State
• number of customers in system
• Events
• customer arrival
• customer departure

8
Discrete-Event Simulation

• Example M/M/1 queueing system
• Arrival times 0.5, 1, 2.5, 4.5
• Service times 1.5, 1.5, 0.5, 1

Customer 1 arrives and starts service
Customer 2 arrives and joins queue
Customer 1 finishes service customer 2 starts
service
Customer 3 arrives and joins queue
Customer 3 finishes service
Customer 4 arrives and starts service
Customer 2 finishes service customer 3 starts
service
Customer 4 finishes service
9
Discrete-Event Simulation

• Discrete-Event Simulation
• The operation of a system is represented as a
  chronological sequence of events.
• Each event occurs at an instant in time and marks
  a change of state in the system.
• Although discrete-event simulation can
  conceptually be done by hand calculations, the
  amount of data that must be stored and
  manipulated for most real-world systems dictates
  that discrete-event simulations be done on a
  digital computer.

10
Discrete-Event SimulationComponents

• System state The collection of state variables
  necessary to describe the system at a particular
  time
• Simulation clock A variable giving the current
  value of simulated time
• Event list A list containing the next time when
  each type of event will occur
• Statistical counters Variables used for storing
  statistical information about system performance

11
Discrete-Event Simulation
12
Discrete-Event Simulation
13
Discrete-Event Simulation
14
Discrete-Event Simulation
15
Discrete-Event Simulation
16
Discrete-Event Simulation
17
Discrete-Event Simulation
18
Discrete-Event Simulation
19
Discrete-Event Simulation
20
Discrete-Event Simulation

• Average number in system 6.5/5.5 1.18
• Server utilization 4.5/5.5 0.82

21
Discrete-Event SimulationComponents

• Initialization routine A subprogram to
  initialize the simulation model at time 0
• Timing routine A subprogram that determines the
  next event from the event list and then advances
  the simulation clock to the time when the next
  event is to occur
• Event routine A subprogram that updates the
  system state when a particular type of event
  occurs (there is one event routing for each type
  of event)
• Library routines A set of subprograms used to
  generate random observations from probability
  distributions that were determined as part of the
  simulation model

22
Discrete-Event SimulationComponents

• Report generator A subprogram that computes
  estimates (from the statistical counters) of the
  desired measures of performance and produces a
  report when the simulation ends
• Main program A subprogram that invokes the
  timing routing to determine the next event and
  then transfers control to the corresponding event
  routine to update the system state appropriately.
  The main program may also check for termination
  and invoke the report generator when the
  simulation is over.

23
Discrete-Event SimulationLogical Flow
24
Discrete-Event SimulationStopping Rules

• Number of events of a certain type reached a
  pre-defined value
• Example stop M/M/1 simulation after the 1000th
  departure
• Simulation time reaches a certain value this is
  usually implemented by scheduling a
  end-simulation event at the desired simulation
  stop time.

25
Monte Carlo Simulation

• A scheme employing random numbers which is used
  to solve certain stochastic or deterministic
  problems where the passage of time plays no
  substantive role.
• Common problem is estimation of
  where f is a function, x is a vector and O is
  domain of integration.
• Special case Estimate for scalar
  x and limits of integration a, b

26
Monte Carlo Simulation

• Let X be a uniform random variable on the
  interval a, b with density
• and let x1, , xn be a random sample from X. Then

27
Monte Carlo Simulation

• Example Estimate
• We approximate this by
• where x1, , xn are a sample from a uniform 0,
  b random variable.

28
Monte Carlo Simulation

• Example Estimate
• There is considerable variability in the quality
  of solution accuracy of numerical integration
  sensitive to integrand and domain of integration

29
Monte Carlo Simulation

• Example (Ross 9.27) Consider a two-out-of-three
  system of independent components, in which each
  components lifetime (in months) is uniformly
  distributed over (0, 1). The reliability function
  is given by

30
Monte Carlo Simulation

• Example (Ross 9.27) cont. Use Monte Carlo to
  estimate the expected lifetime given by
• where Fi(t) is the lifetime distribution for
  component i, i 1, 2, 3, and

31
Monte Carlo Simulation

• Example (Ross 9.27) cont.
• where pi is a sample from a uniform (0,1) R. V.

32
Sources of Randomness for Common Simulation
Applications
33
Random Number Generation

• Random variate1 generation plays a key role in
  simulation.
• The most basic step is the generation of random
  variates from a uniform distribution on (0, 1)
  (denoted U(0,1) and called random numbers) .
• Computer-based random number generators produce
  randomness from a precise algorithm initiated
  using a seed.
• Random variates having other specified
  distributions are built from the basic U(0,1)
  numbers produced with the random number generator.

1 particular outcome of a random variable
34
Random Number Generation

• Pseudo random number generators (PRNG) are
  algorithms that can automatically create long
  runs (for example, millions of numbers long) of
  random numbers with good random properties but
  eventually the sequence repeats exactly.
• One of the most common PRNG is the linear
  congruential generator, which uses the recurrence
• where m (the modulus), a (the multiplier), c (the
  increment), and Z0 (the seed) are suitably chosen
  non-negative integers. The quantity Zn/m is taken
  as an approximation to a uniform (0, 1) random
  variable.

35
Random Number Generation

• Good PRNGs have several properties
• The numbers produced appear to be distributed
  uniformly on (0, 1) and should not exhibit any
  correlation with each other.
• The generator is fast and avoids the need for a
  lot of storage.
• We are able to reproduce a given stream of random
  numbers exactly.
• There is a provision for easily producing
  separate streams of random numbers.

36
Random Number Generation

• Most computer programming languages include
  functions or library routines that purport to be
  random number generators. Such library functions
  often have poor statistical properties and some
  will repeat patterns after only tens of thousands
  of trials. These functions may provide enough
  randomness for certain tasks but are unsuitable
  where high-quality randomness is required, such
  as in cryptographic applications, statistics or
  numerical analysis.
• Much higher quality random number sources are
  available on most operating systems for example
  /dev/random on various BSD flavors, Linux, Mac OS
  X, IRIX, and Solaris, or CryptGenRandom for
  Microsoft Windows.

37
Simulating Continuous Random Variables

• The Inverse Transformation Method
• The Acceptance-Rejection Method

38
The Inverse Transformation Method

• Proposition. Let U be a uniform (0, 1) random
  variable and F be a continuous distribution
  function. Then the random variable
• X F?1(U)
• has distribution function F.
• Proof Since F(x) is a monotone function, then

39
The Inverse Transformation Method
1
F
U
0
X
40
Uniform Distribution

• If
• where u is a U(0, 1) random variate, then
• is a U(a, b) random variate.

41
Exponential Distribution

• If
• where u is a U(0, 1) random variate, then
• is an exponential random variate.

42
Exponential Distribution

• 1000 Exponential Random Numbers (? 2)

43
Erlang-m Distribution

• If X is an Erlang-m random variable with mean ß,
  then
• where the Yis are IID exponential random
  variables, each with mean ß/m. If we use the
  inverse-transform method to generate the
  exponential Yis, then

44
Erlang-m Distribution

• Hence we generate random variate x as follows
• Generate u1, u2, , um as IID U(0,1).
• Return

45
Acceptance-Rejection Method

• Let X have density function f (x). Let g(x) be
  another density function such that there exists a
  number c satisfying
• To generate a random variate with density
  function f (x),
• Generate u from U(0,1).
• Generate a random variate y from the density
  function g, independent of u.
• If u f (x)/cg(x) then set x y. Otherwise, go
  to step 1.

46
Normal Distribution

• Note that if Z has distribution N(0, 1) then its
  distribution function F satisfies
• Therefore, to generate a random variate z, we can
  first generate a nonnegative random variate x
  with density function
• and then assign x a random sign (positive or
  negative with equal probability) to get z.

47
Normal Distribution
48
Normal Distribution

• If we choose g(x) e-x, x gt 0, and
  then
• Hence, to simulate X we
• Generate u from U(0,1).
• Generate a random variate y from the density
  function g, independent of u.
• If u exp?(y ? 1)2 /2, then set x y.
  Otherwise, go to step 1.
• We now generate z by letting z be equally x or x.

49
Simulating Discrete Random Variables

• Inverse Transformation Method

mapped to x3
50
Bernoulli Distribution

• Let X be a Bernoulli random variable with success
  probability p. To generate a Bernoulli random
  variate,
• Generate u from U(0,1).
• If u p, then set x 1. Otherwise set x 0.

51
Discrete Uniform Distribution

• Let X be a discrete uniform random variable
  taking values i, i 1, , j. To generate a
  discrete uniform random variate,
• Generate u from U(0,1).
• Set x i ?(j i 1) u?.

52
Geometric Distribution

• Let X be a geometric random variable with success
  probability p (and q 1 p). Note that
• Setting this equal to u and solving for n yields
• To generate a geometric random variate,
• Generate u from U(0,1).
• Set x ?ln u / ln q?.

53
Binomial Distribution

• Let X be a binomial random variable with
  parameters n and p. Note the sum of n IID
  Bernoulli(p) random variables is a Binomial(n, p)
  random variable.
• To generate a binomial random variate,
• Generate y1, y2, , yn as IID Bernoulli(p) random
  variates.
• Set x y1 y2 yn .

54
Poisson Distribution

• Let X be a Poisson random variable with parameter
  ? and probability mass function p(n). Now
  consider an exponential random variable Y with
  rate ?. Note that
• If Y is the interarrival distribution for some
  event, then the probability that n events occur
  in 0, 1 is given by
• The number of events that occur in the interval
  0, 1 is n if and only if the sum of the first n
  interarrival times is no more than one, but the
  sum of the first n 1 interarrival times is
  greater than one.

55
Poisson Distribution

• Hence, if yi is the ith interarrival time, and ui
  is the random number used to generate it, then
• This implies that we should multiply independent
  U(0,1) random variates until the product falls
  below e??.

56
Poisson Distribution

• To generate a Poisson random variate,
• Let a e??, b 1, and n 0.
• Generate un1 U(0,1) and replace b by bun1. If
  b lt a, then let x n. Otherwise go to step 3.
• Replace n by n 1 and go to step 2 .

57
Estimation of Distribution Parameters Maximum
Likelihood Estimation

• Suppose that Y1, , Yn are continuous random
  variables with respective densities fi(y ?) that
  depend on some common parameter ? (which can be
  vector-valued). Assume that
• ? is unknown
• we observe y1, , yn.
• We want to estimate the value of ? associated
  with Y1, , Yn . Intuitively, we want to find the
  value of ? that is most likely to give rise to
  the data sample y1, , yn.

58
Estimation of Distribution Parameters Maximum
Likelihood Estimation

• Assume that the observations are independent. We
  define the likelihood function
• In maximum likelihood estimation, we choose the
  value of ? that maximizes the likelihood
  function.

59
Estimation of Distribution Parameters Maximum
Likelihood Estimation

• Furthermore, since logarithm is a monotone
  increasing function, then the value of ? that
  maximizes () also maximizes the log of the
  likelihood function
• A similar technique applies to discrete R.V.s
  employing the probability mass function.

60
Estimation of Distribution Parameters Maximum
Likelihood Estimation

• Example Suppose that we assume that some data
  y1, , yn associated with our system is
  exponentially distributed with parameter ?. To
  estimate ? we calculate the likelihood function
• and the log-likelihood function
• Solving for ? gives

61
Tests for Determining How Representative the
Fitted Distributions Are

• Histograms/Frequency Comparisons
• Probability Plots
• Quantile-quantile (Q-Q) plots
• Probability-probability (P-P) plots

62
Tests for Determining How Representative the
Fitted Distributions Are

• Histograms/Frequency Comparisons
• Probability Plots
• Quantile-quantile (Q-Q) plots graph of the
  quantiles of the fitted distribution vs. the
  quantiles of the sample distribution
• Probability-probability (P-P) plots graph of the
  fitted cumulative distribution vs. the sample
  cumulative distribution

63
Tests for Determining How Representative the
Fitted Distributions Are

• Quantile-quantile (Q-Q) plots
• A graph of the quantiles of the fitted
  distribution vs. the quantiles of the sample
  distribution
• If the two sets of quantiles come from
  populations with the same distribution, then the
  points should fall approximately along a
  45-degree reference line.
• The greater the departure from this reference
  line, the greater the evidence for the conclusion
  that distribution for the sample data is
  different from the fitted distribution.

64
Quantile-Quantile Plot
65
Tests for Determining How Representative the
Fitted Distributions Are

• Probability-probability (P-P) plots
• A graph of the fitted cumulative distribution vs.
  the sample cumulative distribution
• The P-P plot is formed by comparing, for each
  value in the ordered sample
• The proportion of points in the sample that are
  less than or equal to that sample value
• The probability of a value less than or equal to
  the sample value based on the fitted distribution
  
• Departures from a straight line indicate
  departures from the specified distribution.

66
Tests for Determining How Representative the
Fitted Distributions Are

• Probability-probability (P-P) plots
• Example
• A sample of 1000 values is fitted to an
  exponential distribution with parameter ? 2.
• The 500th value in the ordered sample is 0.3644,
  i.e., ½ of the data points do not exceed 0.137.
• For the fitted exponential distribution, the
  probability of seeing a value less than or equal
  to 0.3644 is 0.5175.
• Hence (0.5, 0.5175) is a point on the P-P plot.

67
Probability-Probability Plot
68
Homework

• Show how to use the inverse-transform method to
  generate random variates from a Weibull
  distribution.
• Use the method discussed in class to generate 500
  N(0,1) random variates. Use a P-P plot to show
  that the random variates match a normal
  distribution.
• Use the relationship between the negative
  binomial and geometric distributions to give an
  algorithm for generating negative binomial
  variates.

69
Homework

• For an M/M/1 queueing system with interarrival
  rate 2 and service rate 3
• Generate arrival times and service times for 5
  customers.
• Simulate this system until all customers depart,
  showing the system clock, system state and
  scheduled events.
• Use your simulation to calculate the average
  customer delay.
• Use Monte Carlo simulation to estimate the
  expected life of a three-component parallel
  system in which each component has a U(0, 5)
  lifetime distribution.

70
References

• Sheldon M. Ross, Introduction to Probability
  Models, Ninth Edition, Elsevier Inc., 2007.
• Averill M. Law, W. David Kelton, Simulation
  Modeling and Analysis, Third Edition, 2000.
• James C. Spall, Introduction to Stochastic Search
  and Optimization Estimation, Simulation, and
  Control, Wiley Interscience, 2003.