Organisatie. Essay College 1: Introductie Monte Carlo Simulatie Discrete Event Simulatie Java. Enige DES boeken. Informatie

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 (5 en 7 januari);

I

Groepjes van twee;

I

Simulatieprogramma van een wachtrijmodel in Java;

I

Cases;

I

Verslag inleveren eind vd maand;

I

Later mondelinge toelichting;

5 Januari 2016

c Ad Ridder (VU)

Essay – January 2016

1 / 83

Informatie

c Ad Ridder (VU)


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)


Lawrence Leemis and Steve Park. Discrete-Event Simulation: A First Course, 2004.

3 / 83

c Ad Ridder (VU)


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)


5 / 83

Studie van een Stochastisch Systeem

c Ad Ridder (VU)


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


7 / 83

c Ad Ridder (VU)


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)


9 / 83

Probability Theory

c Ad Ridder (VU)


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)


11 / 83

c Ad Ridder (VU)


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


13 / 83

Random Estimation Interval

I

c Ad Ridder (VU)


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)


≈1−α

b I(`)

15 / 83

c Ad Ridder (VU)


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

(ii). Strongly consistent estimator: S2 → σ 2 ; i.e., lim

N 2 1 X S = H(Xi ) − `b . N − 1 i=1

N→∞

2

I

I

`b − ` D `b − ` D √ ∼ N(0, 1) ⇒ √ ∼ tN−1 , σ/ N S/ N

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);

b by its estimate Hence, replace in your computations SE[`] √ c = S/ N. SE

c Ad Ridder (VU)


17 / 83

Normal vs Student t

I


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

(iii).

Mathematically equivalent: S2 =

P(|S2 − σ2 | > ) = 0

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


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)


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)


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)


22 / 83

Modeling aspects

As a reference model we consider a queueing network model.

c Ad Ridder (VU)


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)


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)


25 / 83

Entities


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)


27 / 83

c Ad Ridder (VU)


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)


29 / 83

Event list or calendar


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)


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


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’, ] for the arrival event at 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, ] 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)


33 / 83

The 1-st event

c Ad Ridder (VU)


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


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:


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


Lev = {[‘A1’, 1.861], [‘A2’, 2.043], [‘D2’, 1, 2.510], [‘D1’, 1, 3.987]}

35 / 83

c Ad Ridder (VU)


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


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]}


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)


37 / 83

The k-th event (cont’d)

c Ad Ridder (VU)


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


(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)


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)


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)


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)


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)


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)


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)


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


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


t


I

47 / 83

then we wish to compute

c Ad Ridder (VU)

E[Y].


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)


49 / 83

E

Interpretation of the integral: area below graph of the function Q(t).

c Ad Ridder (VU)


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)


51 / 83

c Ad Ridder (VU)


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)


E

I

53 / 83

Limit Distribution

Question: what do we mean by random variable X(∞)?

c Ad Ridder (VU)


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)


55 / 83

c Ad Ridder (VU)


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)


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)


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


59 / 83

D

c Ad Ridder (VU)


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 + τ ].


61 / 83

Discussion

observation

c Ad Ridder (VU)


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


63 / 83

c Ad Ridder (VU)

1 t

t

Z

H X(s) ds. 0


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)


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)



N 1 X `b = Yi N i=1

67 / 83

c Ad Ridder (VU)


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 τ


batch 3 length τ

69 / 83

Confidence Interval

I

I

c Ad Ridder (VU)


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)


4. Programma schrijven en uitvoeren 5. Output statistiek 6. Verificatie

71 / 83

c Ad Ridder (VU)


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)


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


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


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)


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)


DES datastructuren

77 / 83

Dynamische vectoren


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)


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)


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)


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

c Ad Ridder (VU)


83 / 83

c Ad Ridder (VU)


82 / 83

Organisatie. Essay College 1: Introductie Monte Carlo Simulatie Discrete Event Simulatie Java. Enige DES boeken. Informatie

Recommend Documents