Organisatie
Essay
College 1:
Introductie
Monte Carlo Simulatie
Discrete Event Simulatie
Java
A.A.N. Ridder
Department EOR
Vrije Universiteit Amsterdam
Homepage:
http://personal.vu.nl/a.a.n.ridder/essay/default.htm
I
Twee keer college (6 en 8 januari);
I
Groepjes van twee;
I
Simulatieprogramma van een wachtrijmodel in Java;
I
Cases;
I
Verslag inleveren eind vd maand;
I
Later mondelinge toelichting;
6 Januari 2015
c Ad Ridder (VU)
Essay – January 2015
1 / 83
Informatie
c Ad Ridder (VU)
Essay – January 2015
2 / 83
Enige DES boeken
Wachttijdtheorie
Hoofdstuk 5 in
Henk C. Tijms
A First Course in Stochastic Models
Wiley, 2003
Jerry Banks, John Carson, Barry Nelson and David Nicol.
Discrete-Event System Simulation, Fourth edition, Pearson 2005.
Website (Internet)
http://personal.vu.nl/a.a.n.ridder/
essay/
´
Christos Cassandras and Stephane
Lafortune. Introduction to Discrete
Event Systems, Second edition, Springer 2008.
Louis Birta and Gilbert Arbez. Modelling and Simulation. Exploring
Dynamic System behaviour, Springer 2007.
Averill M. Law and W. David Kelton. Simulation Modeling and Analysis,
Third edition, McGraw-Hill 2000.
Website (Blackboard)
https://bb.vu.nl/
c Ad Ridder (VU)
Essay – January 2015
Lawrence Leemis and Steve Park. Discrete-Event Simulation: A First
Course, 2004.
3 / 83
c Ad Ridder (VU)
Essay – January 2015
4 / 83
Onderwerpen
Voorkennis en -kunde
1. Monte Carlo simulatie
2. Discrete-event simulatie
I
Kansrekening (SLLN; CLT)
3. Voorbeeld van een discrete-event simulatie
I
Markov-ketens (DTMC; CTMC)
4. Eindige en oneindige horizon simulatie
I
Statistiek (BTI)
5. Ingredienten van een simulatiestudie
I
Simulatie uit verdelingen (inverse transformatie; acceptatie-rejectie)
6. Kenmerken van een simulatieprogramma
I
Java programmeren (Linked List)
7. Wachtrijmodellen
8. Wachttijdtheorie
c Ad Ridder (VU)
Essay – January 2015
5 / 83
Studie van een Stochastisch Systeem
c Ad Ridder (VU)
Essay – January 2015
Methoden van studie
I
Veel bedrijfseconomische en financiele systemen bevatten onzekere
(stochastische) componenten;
I
Denk aan
I
Exact (analytisch): Markov-ketens, kansgenererende functies, Laplace
transformaties, complexe analyse.
Risico; levensduur; bederf;
I
Benaderingen: numerieke algoritmes.
Volatiliteit; koersverloop;
I
Simulatie: last resort. Meestal als discrete event simulation (DES).
I
Transporttijden; operatietijden; wachttijden;
I
Congestie; ongelukken;
I
I
I
Onderzoek en studie van het systeem nodig om die onzekerheden te
quantificeren;
I
Bv, een ziekenhuis moet weten wat de gemiddelde ‘ligtijd’ is om te
beoordelen of het aantal bedden voldoende is.
c Ad Ridder (VU)
6 / 83
Essay – January 2015
7 / 83
c Ad Ridder (VU)
Essay – January 2015
8 / 83
Simulation Model
Monte Carlo Simulatie
I
The basis of simulation is the static simulation model.
I
Also called Monte Carlo simulation.
I
A probability model is determined by random variables (or process)
X = (X1 , X2 , . . .) (input variables).
I
The output of the model is the variable Y = H(X) for some function H
(output variable).
I
The goal is to compute the performance measure
Z
` = [Y] = [H(X)] =
H(x)f (x) dx,
E
E
X
where f (x) is the pdf of the random vector X;
c Ad Ridder (VU)
Essay – January 2015
9 / 83
Probability Theory
c Ad Ridder (VU)
Essay – January 2015
10 / 83
The Monte Carlo Method
SLLN
D
When X1 , X2 , . . . are iid replications of X ∼ f then
N
1 X
a.s.
H(Xi ) →
N i=1
E[H(X)]
(N → ∞).
I
Consider finitely many iid replications X1 , X2 , . . . , XN ;
I
This is called a random sample with sample size N;
I
The mean of the associated H(Xi ) variables is a random variable; i.e., a
function on the sample space; called sample mean:
I
I
Make use of the last property;
I
By executing the simulation program, you generate observations (or
data; realizations; values; ...) x1 , x2 , . . . , xN of the random sample;
I
Compute the average of the data:
`=
N
1 X
H(xi );
N i=1
I
Then (SLLN): ` ≈ `: we say that ` is an estimate of `;
N
1 X
`b =
H(Xi ).
N i=1
I
Question: how good is the estimate?
I
Actually; it is custom in the simulation literature to just consider and
analyse the estimators;
E`b = `; i.e., `b is an unbiased estimator of the performance measure `;
I
Thus; how good is the sample mean estimator?
Properties:
I
I
For large N, according to SLLN, `b is not so random; i.e., it is almost a
constant function; i.e., `b ≈ ` for almost all samples.
c Ad Ridder (VU)
Essay – January 2015
11 / 83
c Ad Ridder (VU)
Essay – January 2015
12 / 83
Statistics
Using the CLT
Let σ 2 =
Var[H(X)] be the variance of the output variable;
CLT
D
When X1 , X2 , . . . are iid replications of X ∼ f , and `b is the sample mean based
on N samples; then
√
D
N `b − ` → N(0, σ 2 )
D
D
Interpretation: `b ≈ ` + N(0, σ 2 /N) ∼ N(`, σ 2 /N);
f`b(y)
`q = is q-th quantile of
b
distribution of `:
I
√
b ≈ σ/ N;
First: the standard error of the estimator is SE[`]
I
Suppose you estimated ` by ` based on a simulated sample of size N;
I
You might consider ` as a random realisation of of the estimator `b ∼ f`b(y).
I
Thus (see graph previous slide),
√
|` − `| < 1.96 σ/ N ⇔ 0.025 <
I
P(`b ≤ `q ) = q
D
P `b ≤ `
< 0.975
In other words, when you would generate many independent estimates
(each based on sample N), about 95% of these would be within 1.96
standard errors of the true value.
Note:
`0.025
` = `0.5
c Ad Ridder (VU)
`0.975
√
`0.975 ≈ `+1.96 σ/ N
y
Essay – January 2015
13 / 83
Random Estimation Interval
I
c Ad Ridder (VU)
Essay – January 2015
14 / 83
Confidence Interval
Generally, let α ∈ (0, 1) be your uncertainty parameter such as α = 0.05
or α = 0.1.
2
Var[Y] of the output
I
Let ` be your estimate, and suppose you know σ =
variable.
I
Let zq be the q-th quantile of the standard normal distribution; i.e.,
Φ(zq ) = (N(0, 1) ≤ zq ) = q. We saw
√ |`b − `| ≤ z1−α/2 σ/ N ≈ 1 − α
| {z }
P
I
The analysis above is the basis for reporting simulation output;
I
Your estimate is ` computed as the average of N iid output data obtained
by simulating N times the output variable Y = H(X);
I
You get a realisation of the random estimation interval:
b ` + z1−α/2 SE[`]
b
I(`) = ` − z1−α/2 SE[`],
P
I
You report this interval as a 100(1 − α)% confidence interval.
I
Interpretation: when you would generate many independent estimates
(each based on sample N), and when you would calculate the
associated confidence intervals, then about 100(1 − α)% of these
confidence intervals would cover the true value.
b
SE[`]
I
Rewrite: the unknown parameter ` satisfies
b `b + z1−α/2 SE[`])
b
P ` ∈ (|`b− z1−α/2 SE[`],
{z
}
random estimation interval
c Ad Ridder (VU)
Essay – January 2015
≈1−α
b
I(`)
15 / 83
c Ad Ridder (VU)
Essay – January 2015
16 / 83
Unknown Standard Error
Properties
I
Almost always the standard error of the estimator is unknown;
I
Equivalently, the variance σ 2 =
variable is unknown;
I
Thus variance can be estimated by the sample variance
Properties of the sample variance.
Var[Y] = Var[H(X)] of the output
(i). Unbiased estimator:
E[S2 ] = σ2 ;
P
lim
N
2
1 X
S =
H(Xi ) − `b .
N − 1 i=1
2
I
N→∞
(4.6)
N
N 1 X 2
H (Xi ) − `b2 .
N − 1 N i=1
the (Student) t distribution with N − 1 degrees of freedom.
(iv). Asymptotic normality:
√
D
N `b − ` /S → N(0, 1);
Essay – January 2015
17 / 83
Normal vs Student t
I
Essay – January 2015
18 / 83
E
Summary. Suppose we wish to compute a performance measure ` = [Y]
where the output variable Y = H(X) can be calculated as a function H of a
sequence of random input variables X. Suppose σ 2 = ar[Y] also unknown.
According to property (iii), we should use the Student t distribution in
stead of the normal distribution when constructing confidence intervals:
√
√ `b − tN−1,1−α/2 S/ N, `b + tN−1,1−α/2 S/ N ,
V
Monte Carlo algorithm
1. Repeat for i = 1, . . . , N:
(i). Generate Xi independently of previous runs;
(ii). Compute output Yi = H(Xi ).
According to property (iv), there is not much difference between the
critical points for large sample size N:
α
t9,1−α/2
t49,1−α/2
t99,1−α/2
t249,1−α/2
t999,1−α/2
z1−α/2
c Ad Ridder (VU)
c Ad Ridder (VU)
The Monte-Carlo Algorithm
where tN−1,q is the q-th quantile for the t-distribution with N − 1 degrees
of freedom.
I
∀ > 0
`b − ` D
`b − ` D
√ ∼ N(0, 1) ⇒
√ ∼ tN−1 ,
σ/ N
S/ N
b by its estimate
Hence, replace in your computations SE[`]
√
c = S/ N.
SE
c Ad Ridder (VU)
P(|S2 − σ2 | > ) = 0
(iii).
Mathematically equivalent:
S2 =
I
(ii). Strongly consistent estimator: S2 → σ 2 ; i.e.,
0.15
1.574
1.462
1.451
1.443
1.441
1.440
0.1
1.833
1.677
1.660
1.651
1.646
1.645
0.05
2.262
2.010
1.984
1.970
1.962
1.960
Essay – January 2015
2. Compute sample mean estimator for `:
0.025
2.685
2.312
2.276
2.255
2.245
2.241
N
1 X
`b =
Yi ;
N i=1
3. Compute sample variance estimator for σ 2 :
S2 =
N
1 X
b 2;
(Yi − `)
N − 1 i=1
4. Report estimate, confidence interval, and/or (estimated) standard error.
19 / 83
c Ad Ridder (VU)
Essay – January 2015
20 / 83
What is DES?
Discrete-event simulation is a computer simulation model/program of a
stochastic system that evolves dynamically in time via state information which
changes at discrete time epochs.
The components of a DES are
Discrete Event Simulation
c Ad Ridder (VU)
Essay – January 2015
21 / 83
A queueing network example
I
Entities
I
Time variable (or simulation clock)
I
System state variables
I
Events and Event list (or calendar)
I
Global variables
I
Statistics collectors (or counter variables)
c Ad Ridder (VU)
Essay – January 2015
22 / 83
Modeling aspects
As a reference model we consider a queueing network model.
c Ad Ridder (VU)
Essay – January 2015
23 / 83
I
The queueing model consists of two stations (‘queues’) Q1 and Q2 .
I
At Qi jobs arrive according to a Poisson(λi ) process.
I
There are three servers at Q1 and two servers at Q2 .
I
The service times of the servers at Qi have distribution function Gi (·).
I
After service completion at Q1 the jobs enter Q2 .
I
After service completion at Q2 a job leaves the system with probability p
or re-enters (‘feeds back’ to) Q1 with probability 1 − p.
I
Both stations have infinite waiting spaces.
I
Waiting jobs are served in order of arrival (FCFS).
c Ad Ridder (VU)
Essay – January 2015
24 / 83
Performance measures
Transient vs Steady-State
Typically we are interested in
I
waiting times or system times at the two stations;
I
waiting lines or system lines at two stations;
I
time spend in the system (sojourn time);
I
throughput (per station or from the whole system);
I
utilization;
I
...
We have to specify
I
whether we wish to estimate these performance measures for a finite
period; for instance the queueing system operates daily, opening empty
at 8.00 hr in the morning, closing at 18.00 hr in the evening.
I
or whether we wish to estimate steady-state averages; then we assume
that the system operates for an infinite time.
Let’s assume that
c Ad Ridder (VU)
Essay – January 2015
25 / 83
Entities
Essay – January 2015
26 / 83
System state
I
Actually, a system is defined to be a collection of entities, e.g., people or
cars or machines, that act and interact together.
I
Without entities, nothing would happen.
I
Entities have attributes, usually given as data values.
I
In our example the entities are (attributes between brackets)
I
c Ad Ridder (VU)
I
I
I
The system state is a collection of variables necessary to describe a
system at a particular time.
The set of all states is denoted by X ; a specific state by x ∈ X ; and the
(random) state at time t by X(t) ∈ X .
In our example the state comprises
I
jobs or customers (arrival time);
I
the number of jobs (x1 , x2 ) present at the two stations;
I
servers (idle/busy).
I
two vectors a1 , a2 of the arrival times of these waiting jobs;
I
two vectors b1 , b2 specifying the status of the servers (idle/busy).
Most often in a simulation study, we do not bother too much with entities.
c Ad Ridder (VU)
Essay – January 2015
27 / 83
c Ad Ridder (VU)
Essay – January 2015
28 / 83
Events
I
I
I
Time or clock variable
An event is an instantaneous occurrence that may change the state or
trigger a state transition.
The set of all possible events is denoted by E ; a specific event by e; the
events active in state x ∈ X by E(x) ⊂ E .
I
Events occur at some point in time.
I
For this we need a variable representing the current time of the
simulation.
I
This is also called the simulation clock or system time that measures the
elapsed simulation time.
I
Clearly we need to specify its unit size (second or minute or ...); then we
set it to zero at the start of a simulation and update it every time an event
occurs.
In our example events are
I
the arrival of a new job (at Q1 or at Q2 );
I
a service completion at one of the five servers.
c Ad Ridder (VU)
Essay – January 2015
29 / 83
Event list or calendar
Essay – January 2015
30 / 83
A Set of Alarm Clocks
I
The calendar for the simulation is a list of the events that are currently
scheduled to occur.
I
There is only one event list and it consists of the scheduled event times,
sorted by the earliest scheduled time first.
I
The event list at time t is denoted by Lev (t).
I
In our example the event list comprises
I
the next arrival time of a customer at Q1 ;
I
the next arrival time of a customer at Q2 ;
I
for any busy server the next time he is ready completing the job.
c Ad Ridder (VU)
c Ad Ridder (VU)
Essay – January 2015
I
Another view of the event list is that it is a collection of alarm clocks, one
for each scheduled event. The clocks are preset at different (random)
times in the future. For instance at some t there are 3 events scheduled:
(
)
Lev (t) =
arrival Q1
31 / 83
c Ad Ridder (VU)
arrival Q2
service
completion Q2
Essay – January 2015
32 / 83
The Event-Scheduling Approach
A walk through DES
In our queueing example we let time be measured in seconds, starting at
08:00. This is system time 0.
I
The discrete-event simulation runs as follows.
I
Suppose current time is tsim and state x ∈ X .
I
Each active event ei ∈ E(x) has an associated scheduled alarm time ti .
I
We denote the type-time pair of events by [‘A1’, <time>] for the arrival
event at time <time> at Q1 (and similarly for arrival event at Q2 );
I
The simulation clock is advanced forward to the earliest time on the event list
(when the first alarm goes off).
I
The alarm belongs to a particular event that triggers several activities in the
system.
I
I
These activities make that there will be changes in the system state and in
the event list; these need to be updated.
and by [‘D1’, 1, <time>] for departure events due to service completion
at server 1 of Q1 ; similarly for the other servers, and the other queue.
I
We start tsim = 0 with an empty system and idle servers:
I
Then the simulation clock is advanced again, etc.
x = {x1 = 0, x2 = 0, a1 = [], a2 = [], b1 = (0, 0, 0), b2 = (0, 0)}.
I
There are two active events: arrivals at the queues; the associated event
times are drawn by our random number generator, resulting. Thus
Lev = {[‘A2’, 0.817], [‘A1’, 1.284]}.
c Ad Ridder (VU)
Essay – January 2015
33 / 83
The 1-st event
c Ad Ridder (VU)
Essay – January 2015
The 2-nd event
I
The simulation time is advanced to the earliest event time 0.817
belonging to the event of an arriving job at Q2 ;
I
The simulation time is advanced to the earliest event time 1.284
belonging to the event of an arriving job at Q1 .
I
Immediately, this job enters service with server 1 of Q2 ;
I
This job enters service with server 1 of Q1 .
I
We draw a sample of the service time S2 , say 1.693, to schedule a new
event ‘D2’ with departure time 0.817 + 1.693 = 2.510;
I
We draw a sample of the service time S1 , say 2.613, to schedule a new
event ‘D1’ with departure time 1.284 + 2.613 = 3.987;
I
Also we realise a new sample of the interarrival time A2 , say 1.226, and
update the scheduled time of event ‘A2’ to 0.817 + 1.226 = 2.043;
I
We realise a new sample of the interarrival time A1 , say 0.577, and
update the scheduled time of event ‘A1’ to 1.284 + 0.577 = 1.861;
Thus the new situation is:
Thus the new situation is:
tsim = 0.817
tsim = 1.284
x = {x1 = 0, x2 = 1, a1 = [], a2 = [], b1 = (0, 0, 0), b2 = (1, 0)}
x = {x1 = 1, x2 = 1, a1 = [], a2 = [], b1 = (1, 0, 0), b2 = (1, 0)}
Lev = {[‘A1’, 1.284], [‘A2’, 2.043], [‘D2’, 1, 2.510]}
c Ad Ridder (VU)
34 / 83
Essay – January 2015
Lev = {[‘A1’, 1.861], [‘A2’, 2.043], [‘D2’, 1, 2.510], [‘D1’, 1, 3.987]}
35 / 83
c Ad Ridder (VU)
Essay – January 2015
36 / 83
The 3-rd event
The k-th event
Suppose that currently the situation is
I
The simulation time is advanced to the earliest event time 1.861
belonging to the event of an arriving job at Q1 .
I
This job enters service with server 2 of Q1 .
I
We draw a sample of the service time S1 , say 0.811, to schedule a new
event ‘D1’ with departure time 1.861 + 0.811 = 2.752;
I
We realise a new sample of the interarrival time A1 , say 1.428, and
update the scheduled time of event ‘A1’ to 1.861 + 1.428 = 3.289;
tsim = 249.31
x = {x1 = 4, x2 = 8, a1 = (240.82), a2 = (238.71, 240.61, 241.01, 244.55, 246.91, 248.88),
b1 = (1, 1, 1), b2 = (1, 1)}
Lev = {[‘D2’, 2, 250.28], [‘D2’, 1, 251.43], [‘A1’, 254.38], [‘D1’, 2, 255.36],
[‘A2’, 256.42], [‘D1’, 1, 260.91], [‘D1’, 3, 263.93]}
Thus the new situation is:
I
The next event is service completion at Q2 at time 250.28.
I
There are jobs waiting at Q2 to occupy the empty seat. We draw a
sample of the service time S2 , say 0.83, and update the scheduled time
of event ‘D2’, 2 to 250.28 + 0.83 = 251.11.
I
For the job leaving Q2 we flip a coin (Bernoulli random variable) to decide
a feed back; suppose the outcome is to loop back for entering Q1 again.
I
The looped job finds all servers busy at Q1 , and thus joins the queue.
tsim = 1.861
x = {x1 = 2, x2 = 1, a1 = [], a2 = [], b1 = (1, 1, 0), b2 = (1, 0)}
Lev = {[‘A2’, 2.043], [‘D2’, 1, 2.510], [‘D1’, 2, 2.752], [‘A1’, 3.289], [‘D1’, 1, 3.987]}
c Ad Ridder (VU)
Essay – January 2015
37 / 83
The k-th event (cont’d)
c Ad Ridder (VU)
Essay – January 2015
Trace
Next situation:
I
This detailed walk through a system simulation is called a trace.
tsim = 250.28
I
In a trace a print is made of all the (important) simulation components
after each event.
I
It serves two objectives.
x = {x1 = 5, x2 = 7, a1 = (240.82, 250.28), a2 = (240.61, 241.01, 244.55, 246.91, 248.88),
b1 = (1, 1, 1), b2 = (1, 1)}
(i). In developing a simulation program it helps you to think about the events that
may happen in the system, the activities they trigger and their logic.
Lev = {[‘D2’, 2, 251.11], [‘D2’, 1, 251.43], [‘A1’, 254.38], [‘D1’, 2, 255.36],
[‘A2’, 256.42], [‘D1’, 1, 260.91], [‘D1’, 3, 263.93]}
c Ad Ridder (VU)
38 / 83
Essay – January 2015
(ii). After having written a computer program of the simulation, it provides a
check whether the program is correct.
39 / 83
c Ad Ridder (VU)
Essay – January 2015
40 / 83
Simulation run
Programming a Run
I
A simulation run of the system goes from event to event at discrete time
epochs determined by realisations of the appropriate random variables.
I
A trace covers usually a short time interval or a small number of events.
I
I
The actual execution of the computer program of the simulation is longer
but should end at some time or after some number of events: stopping
criterion.
When you write a computer program of the simulation it is convenient to
modularise your program into several subprograms or routines to clarify
the logic and the interactions.
I
Roughly a simulation run looks in pseudo-code as follows.
I
One such a simulation is called a simulation run which corresponds to
the concept of sample path of a stochastic process.
I
Later we will see that in most simulations studies you will simulate
(many) more simulation runs in order to obtain more reliable estimates.
c Ad Ridder (VU)
Essay – January 2015
initialise;
REPEAT
[e,t] = next_event_type_and_time;
switch (e)
case ’arrival’: ...
case ’departure’: ...
case ’. . .’: ...
UNTIL stopping_criterion;
41 / 83
Programming a Run (cont’d)
Essay – January 2015
42 / 83
The Goal of DES
In the subroutines for the different events you program code for updating
state, event list, and statistical counters. Similar as what you would do in the
trace. For instance:
The goal of a DES-study is to estimate performance measures; for instance :
subroutine execute_arrival_at_Q1_event(t)
begin
update_all_relevant_counter_variables();
j = find_free_server();
if (server(j)==‘idle’)
server(j) = ‘busy’;
s = generate_servicetime();
add_to_eventlist(‘D1’,j,t+s);
else
add_customer_to_queue(j,t);
end
a = generate_interarrivaltime();
add_to_eventlist(‘A1’,t+a);
end
c Ad Ridder (VU)
c Ad Ridder (VU)
Essay – January 2015
I
average waiting times (per queue or overall);
I
mean waiting line (per queue or overall);
I
average time spend in the system (sojourn time);
I
throughput (per queue or overall);
I
utilization;
I
...
We need to say whether we want these quantities measured (averaged)
during a short period (finite-horizon simulation), or in the long-run
(steady-state simulation).
43 / 83
c Ad Ridder (VU)
Essay – January 2015
44 / 83
Procedure
I
A run starts at time 0, at a specific specified state;
I
A run ends at some stopping criterion, e.g.
(i). A time horizon T;
(ii). A number of events Ne ;
Finite-Horizon
(iii). A number of arrivals Na ;
(iv). Etc.
Simulation
c Ad Ridder (VU)
Essay – January 2015
45 / 83
Illustration
I
Execute N runs, independent and (probabilistically) identical.
I
In each run we keep track of certain variables (counter variables) that at
the end of the run determine its outcomes or outputs.
I
The performance measures are estimated by averaging these N
outcomes.
I
The counter variables are initialised at the start of a simulation run, and
because nothing changes in between two consecutive events, it suffices
to update them after an event occurs.
c Ad Ridder (VU)
Essay – January 2015
46 / 83
Example
state
0
T
t
run 1: compute
output
Y1 .
I
Consider Q1 in the tandem network.
I
Suppose the goal is to estimate the performance measure “mean
average waiting time per customer from 08:00-18:00”.
I
This means:
state
0
T
t
run 2: compute
output
Y2 .
c Ad Ridder (VU)
let K be the (random) number of customers that has been served during
these ten hours at Q1 ;
I
let W1 , . . . , WK be their corresponding (random) waiting times;
P
define the output Y = (1/K) K
i=1 Wi ;
I
state
0
I
T
Essay – January 2015
t
run 3: compute
output
Y3 .
I
47 / 83
then we wish to compute
c Ad Ridder (VU)
E[Y].
Essay – January 2015
48 / 83
Example (cont’d)
I
Estimating the Average Queue Length
During a run in the simulation program you keep track of two variables:
1. nserv1 = the current total number customers who went into service at Q1 ;
I
How to compute (or estimate) the mean length at Q1 during the interval
08:00-18:00?
Whenever a new customer enters service at Q1 , these two variables are
updated.
I
At the end of the run you have a realisation y = wtotal1 /nserv1 of the output
Y.
I
Let Q(t) be the queue length at time t, 0 ≤ t ≤ T, where t = 0 represents
08:00 and T represents 18:00; measured in some time-units;
RT
Goal: [Y] for output Y = T1 0 Q(t) dt;
2. wtotal1 = the current total waiting time of these nserv1 customers.
I
I
I
I
Repeat N times and do your statistics.
I
For such finite-horizon simulation the statistics is similar as for static
stochastic simulation (see the part on the Monte Carlo simulation).
c Ad Ridder (VU)
Essay – January 2015
49 / 83
E
Interpretation of the integral: area below graph of the function Q(t).
c Ad Ridder (VU)
Essay – January 2015
50 / 83
Illustration
Q(t)
Steady-State
0 T1 T2
TR T
t
Simulation
Let 0 = T0 < T1 < T2 < · · · < TR < T ≤ TR+1 be the consecutive event times.
Y=
R
1X
1
Tr − Tr−1 Q(Tr−1 ) +
T − TR Q(TR )
T r=1
T
Note: number of events R is random.
c Ad Ridder (VU)
Essay – January 2015
51 / 83
c Ad Ridder (VU)
Essay – January 2015
52 / 83
Steady-State Definition
Looking to the Future
I
Suppose that you describe the system state at time t by X(t);
I
You observe the system at time 0 (‘now’);
I
The system is acting already a long time (from time −∞);
I
You might be interested in
(i). What is the CDF F0 (x) of X(0)?
(ii). What is the expectation
E[X(0)]?
I
You might reason as follows;
I
The system starts at time 0;
I
The system will act forever in the future;
I
The system will be in steady-state at time +∞;
I
Now we are interested in
(i). The CDF F∞ (x) of X(∞);
(ii). The expectation
H X(∞) .
(iii). More generally,
suppose
that H(x) is a function acting on the state variable,
what is
H X(0) ?
E
I
F0 is the steady-state distribution of the system;
I
The expectations are steady-state performance measures.
c Ad Ridder (VU)
Essay – January 2015
E
I
53 / 83
Limit Distribution
Question: what do we mean by random variable X(∞)?
c Ad Ridder (VU)
Essay – January 2015
54 / 83
Convergence in Mean
P(X(t) ≤ x) is a CDF for all t ≥ 0;
I
Assume that Ft (x) =
I
Equivalently, X(t) is a proper random variable;
Definition
D
Assume that X(t) → X(∞). If
E[X(∞)] = t→∞
lim E[X(t)]
Definition
Assume
we say that the sequence of random variables {X(t), t ≥ 0} converges in
mean to X(∞); and we denote
(a). The limit
F∞ (x) = lim Ft (x)
t→∞
exists for all x ∈
R (we say pointwise);
L
X(t) →1 X(∞)
(b). F∞ (·) is a CDF (recall the definition of a CDF!)
Then we say that the sequence of random variables {X(t), t ≥ 0} converges
weakly or in distribution to a (proper!) random variable X(∞); and we denote
Generally you need to check certain sufficient conditions to ensure
connvergence in distribution and/or in mean (see the counterexamples next
slides).
D
X(t) → X(∞)
c Ad Ridder (VU)
Essay – January 2015
55 / 83
c Ad Ridder (VU)
Essay – January 2015
56 / 83
Counterexample I
Counterexample II
I
I
P(X(t) = 0) = (t − 1)/t; P(X(t) = t) = 1/t.
Define random variable X(t) for t ≥ 0 by
P(X(t) = 0) = 0.5; P(X(t) = t) = 0.5.
I
I
I
Hence,
0,
Ft (x) = (t − 1)/t,
1,
Hence,
0,
Ft (x) = 0.5,
1,
I
Define random variable X(t) for t ≥ 1 by
x < 0;
0 ≤ x < t;
x ≥ t.
I
Easy to see that F∞ (x) = limt→∞ Ft (x) exists (pointwise), with
(
0,
x < 0;
F∞ (x) =
1,
x ≥ 0.
I
Conclude X(t) → X(∞), where
I
However
Easy to see that F∞ (x) = limt→∞ Ft (x) exists (pointwise), with
(
0,
x < 0;
F∞ (x) =
0.5,
x ≥ 0.
No convergence in distribution.
D
lim
t→∞
c Ad Ridder (VU)
Essay – January 2015
57 / 83
Steady-State Simulation
I
I
I
P(X(∞) = 0) = 1;
E[X(t)] = t→∞
lim 1 = 1 6= 0 = E[X(∞)].
c Ad Ridder (VU)
Essay – January 2015
58 / 83
A Steady-State Simulation Run
I
D
Suppose that in your system both (i) X(t) → X(∞) and (ii)
L
H X(t) →1 H X(∞) ;
D
NB: you may assume that also H X(t) → H X(∞)
Hence,
(i). you start simulating a run from an arbitrary initial state, for instance the
empty state;
(ii). continue long enough, say until τ0 , at which the system is about in
steady-state;
The performance measure of interest is the steady-state expectation
`=
H X(∞) ;
I
Note that for all t ≥ τ0 it holds that X(t) ≈ X(∞) and
How would you estimate ` by simulation?
I
Thus define as output of the run
Z
1 τ0 +τ
Y=
H( X(t) dt
τ τ0
I
NB: similarly for discrete processes in which you take sums.
(iii). continue simulating until τ0 + τ ;
E
I
x < 0;
0 ≤ x < t;
x ≥ t.
Notice:
I
I
D
From (i): X(∞) ≈ X(τ0 ) for large τ0 ;
From (ii): ` =
H X(∞) ≈
H X(τ0 ) for large τ0 ;
c Ad Ridder (VU)
E
E
Essay – January 2015
59 / 83
D
c Ad Ridder (VU)
Essay – January 2015
EH
X(t)
≈ `;
60 / 83
Replication-Deletion
Illustration
state
I
The output satisfies
Z
h 1 Z τ0 +τ
i
1 τ0 +τ
H( X(t) dt =
[Y] =
τ τ0
τ τ0
E
I
I
E
EH
X(t)
τ
run 2: compute
output
Y2 from state
in
t process
[τ0 , τ0 + τ ].
dt ≈ `;
τ0
0
Repeating N times independently, you get an (approximately) unbiased
estimator of ` by
N
1 X
`b =
Yi
N i=1
warm-up
observation
state
The rest (SLLN; CLT; confidence interval; accuracy; ...) is the same as in
static Monte Carlo simulations.
τ0
0
warm-up
c Ad Ridder (VU)
τ
run 1: compute
output
Y1 from state
in
t process
[τ0 , τ0 + τ ].
Essay – January 2015
61 / 83
Discussion
observation
c Ad Ridder (VU)
Essay – January 2015
More Good News
I
Suppose the performance measure is a long-run fraction or average
quantity; e.g.
I
A simple method for estimating a steady-state expectation;
I
the long-run fraction of time the server is available;
I
The initial long period of a run until steady-state is called warm-up
period, and is not used for computing output;
I
the long-run average cost per unit of time;
I
the long-run fraction of time the machine is down;
The output Y of a run is computed from the observation period;
I
the long-run average number of particles per unit of time;
I
Pro: iid observations which makes the statistics easy;
I
etc.
I
Con: inefficient because each run needs the warm-up period.
I
I
These can be modelled by
lim
t→∞
c Ad Ridder (VU)
62 / 83
Essay – January 2015
63 / 83
c Ad Ridder (VU)
1
t
t
Z
H X(s) ds.
0
Essay – January 2015
64 / 83
Almost Sure Limits
I
The Difficulty of Steady-State Simulation
If the system state process {X(t), t ≥ 0} is regenerative these long-run
average limits converge almost surely to the steady-state expectation
(this is known as an ergodic theorem):
Z
1 t
lim
H X(s) ds =
H X(∞)
a.s.
t→∞ t 0
|
{z
}
random variable
E
I
See your textbook on Markov chains;
I
Corollary: you
estimate
long-run averages/fractions similarly as
estimating
H X(∞) .
I
How large should τ0 (length of warm-up period) be before starting the
observations?
I
This is the initial transient problem;
I
There are statistical techniques, but in practice you just simulate and
experiment;
E
c Ad Ridder (VU)
Essay – January 2015
65 / 83
Steady-State Simulation Strategies
c Ad Ridder (VU)
66 / 83
C. The Batch Means Method
I
I
The replication-deletion method, described in the slides above, is just
one of many methods for steady-state simulation.
I
The most important ones are
Remedy for the inefficiency of the replication-deletion method;
(i). simulate one long run of length T;
(ii). first a warm-up period of length τ0 ;
(iii). then split the remaining time T − τ0 in N sections (called batches) of length
τ = (T − τ0 )/N.
A. Stationary process (single run)
B. Replication-deletion method (multiple runs)
I
Compute the sample mean of the i-th batch (i = 1, . . . , N)
Z
1 τ
Yi =
H X(τ0 + (i − 1)τ + t) dt
τ 0
I
Final estimator is the average of these batch means
C. Batch means method (single run)
D. Regenerative method (multiple cycles)
I
A. and D. are more involved obtaining variance estimates for the
construction of the confidence intervals; not in this course.
I
C. is an efficient alternative for B.
c Ad Ridder (VU)
Essay – January 2015
Essay – January 2015
N
1 X
`b =
Yi
N i=1
67 / 83
c Ad Ridder (VU)
Essay – January 2015
68 / 83
Illustration
Discussion
state
I
Typically, number of batches N is small.
I
The batch length (or size) τ should be large, so that
1. the correlations between the batch means Yi , i = 1, 2, . . . are neglegible;
2. each batch mean Yi is approximately Gaussian.
t
0
I
K
warm-up
c Ad Ridder (VU)
batch 1
length τ
batch 2
length τ
Essay – January 2015
batch 3
length τ
69 / 83
Confidence Interval
I
I
c Ad Ridder (VU)
Essay – January 2015
70 / 83
Ingredienten van een simulatiestudie
The sample variance to estimate v2 :
1. Random number generator
N
2
1 X
S2 =
Yi − `b .
N − 1 i=1
I
Thus the Yi ’s may be considered as independent Gaussian, with mean `
and variance v2 ;
2. Trekken uit kansverdelingen
3. Fitten van parameters
Since N is small:
√ `b − ` D
N
→ tN−1 (t → ∞),
S
where tN−1 is the Student-t variable with N − 1 df.
√ √
CI: `b − tN−1,1−α/2 S/ N, `b + tN−1,1−α/2 S/ N .
c Ad Ridder (VU)
Essay – January 2015
4. Programma schrijven en uitvoeren
5. Output statistiek
6. Verificatie
71 / 83
c Ad Ridder (VU)
Essay – January 2015
72 / 83
1. Random Number Generator (RNG)
RNG in Java
Een RNG produceert een rij getallen u1 , u2 , . . . tussen 0 en 1 die beschouwd
mogen worden als realiseringen van een rij IID uniform verdeelde
stochastische variabelen U1 , U2 , . . ..
I
Gebruik de klasse Random:
long seed = 123456; // kies je eigen seed
Random generator = new Random(seed);
Een RNG wordt als een algoritme ge¨ımplementeerd dat
I
Dan steeds u = generator.nextDouble() geeft een
(pseudo-)random getal tussen 0 (inclusief!) en 1 (exclusief).
I
Alternatief: Mersenne Twister. Een snelle en uitstekende RNG. Code
beschikbaar op Internet.
(i). snel moet zijn;
(ii). reproduceerbaar;
(iii). uitstekende statistische eigenschappen;
(iv). machine-onafhankelijk.
c Ad Ridder (VU)
Essay – January 2015
73 / 83
2. Simuleren uit verdelingen
I
Inverse-transformatie: als X verdelingsfunctie F(·) heeft, en U is de
D
F(x) = 1 − e−λx = 1 − u ⇔ x = − ln(u)/λ.
Acceptatie-rejectie: als X kansdichtheid f (·) heeft, Y kansdichtheid g(·),
en er is een eindige constante c zodat f (x) ≤ cg(x) voor alle x ∈ ; dan
R
D
α α−1 −λx
λ x
e
Γ(α)
Calculus: optimaal is µ = α/λ en c =
c Ad Ridder (VU)
I
Gegeven een geparametrizeerde kansverdelingstype, bv Gamma(α, λ),
Coxian2 (p, µ1 , µ2 ), etc van een stochastische variabele X.
I
Gegeven voldoende momenten m1 =
I
Vind de parameters zodat de momenten kloppen.
Meestal eerste twee momenten.
I
Ipv tweede moment, beter genormaliseerde variantie; die heet
(kwadratische) variatiecoefficient
c2 [X] =
en g(x) = µe−µx .
I
λ
Γ(α)
α α−1
e
Essay – January 2015
E[X], m2 = E[X2 ], . . .
I
is X = Y|(U ≤ f (Y)/(cg(Y))). Algoritme: trek y uit g(·), u uit U(0, 1); als
u ≤ f (y)/(cg(y)) neem x = y, anders herhaal. Bv (voor α > 1)
f (x) =
74 / 83
Probleem:
D
uniforme s.v. op [0, 1]; dan is X = F −1 (U) = F −1 (1 − U). Algoritme: los
F(x) = u (of F(x) = 1 − u) op naar x. Bv
I
Essay – January 2015
3. Fitten van parameters
Algemeen: inverse-transformatie methode of de acceptatie-rejectie methode.
Specifiek: convolutie, compositie, ...
I
c Ad Ridder (VU)
Var[X] = m2 − 1.
m21
m21
Zie de notitie des notes.
.
75 / 83
c Ad Ridder (VU)
Essay – January 2015
76 / 83
6. Verificatie
Verificatie is het proces dat nagaat of je simulatieprogramma correct is.
Enige tips:
I
Programmeer in modules en check deze; bijvoorbeeld een module om
trekkingen uit een Lognormale verdeling te simuleren met verwachting 5
en variantie 20. Ga experimenteel na of 100000 trekkingen inderdaad
deze gemiddeldes heeft.
I
Druk een trace af (eventlist, toestand, statistische tellers na iedere event)
en controleer ‘met de hand’.
I
Controleer bekende eigenschappen zoals Little.
I
Simuleer vereenvoudigde versie van het systeem waarvoor analytische
uitdrukkingen bekend zijn.
I
Neem extreme waarden van parameters.
c Ad Ridder (VU)
Essay – January 2015
DES datastructuren
77 / 83
Dynamische vectoren
Essay – January 2015
78 / 83
Linked List
I
In wachtrijsimulaties heb je te maken met dynamische
vectoren/matrices: geplande events; klanten in de rij (of rijen); klanten in
het systeem; etc.
I
Afhankelijk van model/probleem/simulatieprogramma heb je behalve de
lengte (ofwel ‘aantallen’) ook nodig informatie per element. Bv, de
aankomsttijden van klanten die staan te wachten.
I
Dynamisch wil zeggen dat de lengte varieert gedurende de simulatie.
Implementatie als ‘arrays’ is dan bijzonder inefficient (vrijwel onmogelijk).
c Ad Ridder (VU)
c Ad Ridder (VU)
Essay – January 2015
79 / 83
I
De datastructuur Linked List (gelinkte lijst) is de oplossing.
I
Uit Wikipedia: Een gelinkte lijst bestaat uit een reeks van
´ of twee
data-elementen, records, elk met een aantal datavelden en e´ en
verwijzingen naar een volgende en/of vorig data-element. Het is niet
mogelijk om meteen een willekeurig data-element direct te benaderen in
tegenstelling tot bijvoorbeeld een array: om een element op de N-de
positie te benaderen, moet eerst de hele lijst tot positie N doorlopen
worden.
c Ad Ridder (VU)
Essay – January 2015
80 / 83
Operaties
Zelf programmeren
Het is toegestaan om op elke plek in een gelinkte lijst data-elementen in te
voegen of te verwijderen door de verwijzing van het voorgaande element
naar respectievelijk het nieuw toe te voegen, dan wel het tweede volgend
element te laten wijzen.
Mijn advies: voor wie het kan zelf programmeren:
I class Knoop (of Event of Klant of ...):
I
I
bevat velden voor data, een wijzer naar vorige en/of een wijzer naar
volgende;
I
bevat methode(n) voor constructie
class Lijst (of EventList of Wachtrij of ...):
I
bevat velden voor lengte (vd lijst), en voor knopen current/first/last;
I
bevat methoden voor constructie, invoegen, verwijderen, zoeken, enz.
Een wachtrijsimulatie bevat meestal meer lijsten!
c Ad Ridder (VU)
Essay – January 2015
81 / 83
Gebruik van Java’s klassen
I
Java kent een klasse LInkedList van objecten. Je moet zelf klasse object
invoeren (bv Event, Knoop, etc, zie vorige slide).
I
LinkedList bevat al methoden voor invoegen en verwijderen, en een
aantal anderen. Opzoeken bij documentatie van Java
http://docs.oracle.com/javase/1.4.2/docs/api/java/
util/LinkedList.html
I
Zie ook mijn voorbeeldprogramma LLproef.java op
http://personal.vu.nl/a.a.n.ridder/essay/src.html
c Ad Ridder (VU)
Essay – January 2015
83 / 83
c Ad Ridder (VU)
Essay – January 2015
82 / 83
© Copyright 2024 ExpyDoc
