27 Engs Discrete and Probabilistic Systems
Winter, 2004
LECTURE 19
Discrete Event Simulation
19.1. 19.2. 19.3. 19.4. 19.5.
Stochastic simulation, the big picture Simulating continuous random variables Simulating Poisson processes Simulating an M/M/1 queue Simulating an M/M/k queue
1 2 3 5 9
19.1. Stochastic simulation, the big picture
In many probabilistic systems we can derive a mathematical model from basic principles. The key result of such a model is a probability distribution for the quantity of interest, from which expected values, variances, and other important quantities may be calculated. When analysis fails, we turn to simulation. For example, earlier we examined the robustness of an LC circuit by the following simulation procedure: (1) Generate random values for the two components L and C. 1 (2) Calculate the resonant frequency 0 = with these values. LC (3) Repeat (1) and (2) a large number of times, collecting an ensemble of 0 values. (4) Compute a histogram of 0 to estimate the probability density function f (0 ). Calculate sample mean and sample variance of 0 to estimate the mean and variance of the resonant frequency. All simulation procedures break down into these three parts: generating random input values of appropriate distribution, processing the inputs to obtain an ensemble of output values, and analyzing the outputs statistically to draw conclusions about the process. When the process can be expressed as a simple function of the input Y = g(X), Monte Carlo simulations like the one above work well. Other kinds of systems, notably queues, have inputs which are streams of events rather than numbers. A dierent simulation approach is taken for these systems. Before introducing discrete event simulation per se, we return to a topic that we didnt get to earlier in the course.
Engs 27 Discrete and Probabilistic Systems
Winter, 2004
19.2. Simulating continuous random variables
Let X be a random variable with cumulative distribution function FX . We will show that the random variable U = FX (X) is uniformly distributed on [0, 1]. The cumulative distribution function of U is FU (u) = P (U u) = P (FX (X) u) . Now, assuming that FX is monotonically increasing (a stronger requirement than nondecreasing) 1 as well as continuous, the inverse function FX exists, and the event (FX (X) u) is equivalent 1 to the event (X FX (u)). So,
1 1 FU (u) = P (X FX (u)) = FX FX (u) = u,
0u1.
Now dierentiate FU to get the probability density function fU . fU (u) = which is uniform. We can use this to generate continuous random variables from a desired PDF fX . Use a random number generator to obtain a number U uniformly distributed on (0, 1). Then, the 1 random variable X = FX (U ) will be distributed according to fX . This is called the inverse transform method for generating random numbers. It works when a mathematical expression 1 exists for FX . For example, to generate a value from an exponential distribution, begin with the CDF, F (x) = 1 ex/ . The inverse function is obtained by solving u = 1 ex/ for x: x = log(1 u) . If U is uniformly distributed on (0, 1), then so is 1 U , so we may also write X = log U In Matlab, one line of code will generate a vector of N exponential random numbers: x = -theta * log( rand(1,N) ); Here is a sample result: 104 samples from an exponential distribution with = 2. (19.1) dFU = 1, du 0u1,
19. Discrete Event Simulation
19.2
Engs 27 Discrete and Probabilistic Systems
104 GENERATED RANDOM VARIABLES
0.5 0.45 0.4 0.35 0.3 Sample mean = 2.0049 Sample variance = 4.08
Winter, 2004
f(x)
0.25 0.2 0.15 0.1 0.05 0 0
1
2
3
4
5
6
7
8
9
10
x
19.3. Simulating Poisson processes
When simulating a Poisson process, you may want to simulate how many counts (events) occur in a given interval of length t. Or, you may want to simulate the actual event stream, by generating a sequence of interarrival times. 1 Simulating the number of events in an interval For a Poisson process with rate , we know that the interarrival times are distributed exponentially with mean 1/. To simulate the number of events in an interval T , sum exponential random variables until the sum exceeds the length of the interval, T . A particularly ecient way of doing this follows from the earlier observation that an exponential random variable is 1 proportional to the logarithm of a uniform random variable: X = log U . Let Tk be the k th interarrival time. Then,
n
T or, equivalently,
k=1
Tk =
1
n k=1
1 log Uk = log
n
Uk
k=1
n
eT
Uk .
k=1
(19.2)
Summing exponential variates until the sum exceeds T is equivalent to multiplying uniform variates until the product exceeds eT . Then the desired Poisson variate is one less than this n. Here is Matlab code that does the calculation (as usual, = T ): p = exp(-mu); n = 1;
1S. Ross, Simulation, third edition (Academic Press, 2001), Sections 4.2 and 5.4.
19. Discrete Event Simulation
19.3
Engs 27 Discrete and Probabilistic Systems u = rand; while(u >= p), u = u * rand; n = n+1; end X = n-1;
Winter, 2004
Simulating a sequence of events If instead we want to simulate a sequence of events, i.e., arrival times, we can generate a sequence of exponential random variables, stopping when their cumulative sum exceeds the length of the interval, T . The following Matlab code does the job.
t1 = []; t = 0; while (t < T), t = t - log(rand)/mu; t1 = [t1, t]; end X = length(t1) - 1;
% Set of event times. % Running sum of interarrival times % Exponential r.v. with mean 1/mu
% We also get the number of counts in T
Shown below is an example of an event stream with = 10, generated by this method. The lower plot compares the histogram of event counts, for 1000 trials, with the Poisson probability function.
TYPICAL EVENT STREAM, =10
0
0.1
0.2
0.3
0.4
0.5 TIME
0.6
0.7
0.8
0.9
1
POISSON DISTRIBUTION 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0 5 10 COUNTS 15 20
PROBABILITY FREQUENCY
19. Discrete Event Simulation
19.4
Engs 27 Discrete and Probabilistic Systems
Winter, 2004
19.4. Simulating an M/M/1 queue
Consider an M/M/1 queue with arrival rate and service rate . To simulate the system we will keep track of the time, t and the number of customers in the system, X(t). The simulation will run from t = 0 to t = T . These quantities will change when an event arrival or departure occurs. 2 With t = 0 initially, we generate an arrival time (an exponential random variable with mean 1/), update the time pointer t to this time and increment the number of customers. Then we generate a departure time, since there is a customer waiting to be served. We also generate the next arrival time. This gives us a pair {tA , tD } called the event list. If tA < tD , it means that a new customer has arrived before the rst customer has departed, and the number of customers must be incremented. On the other hand, if tD < tA , the rst customer leaves the system before the second customer enters. Whichever occurs rst, arrival or departure, a new arrival or departure time is generated, replacing the old one in the event list. If a departure leaves the queue empty, then a new departure time is not generated until a customer has arrived. In this way we proceed through time, event by event, until at last an arrival time is greater than the limit T . customer That is not added to the system, and the remaining customers empty out, one departure at a time. A sample code is shown below. It can also be easily modied for nite-capacity queues, tracking how many customers are turned away, etc.
%% DISCRETE EVENT SIMULATION of M/M/1 QUEUE %% E.W. Hansen, Engs 27 02W %% Ref: S. Ross, "Simulation", Sect. 6.2 ra = 1; rs = 1.2; T = 12*60; % A D S Event times = []; = []; = []; % Arrival and service rates (per minute) % Simulation time (12 hours)
% Series of arrival times % Series of departure (service completion) times % Series of service times
% Initial conditions n = 0; t = 0; y = -log(rand)/ra; ta = y; td = inf;
% % % % %
Current number in system Elapsed time Exponential r.v., first interarrival time First arrival time Initialize departure time
2Ross, Sections 6.1-6.2; for a dierent approach, J.J. Higgins and S. Keller-McNulty, Concepts in Probability
and Stochastic Modeling (Duxbury, 1995), pp. 318-320
19. Discrete Event Simulation
19.5
Engs 27 Discrete and Probabilistic Systems
Winter, 2004
while(min(ta, td) <= T), % Case 1: Not time for departure yet, so generate an arrival if(ta <= td & ta <= T), t = ta; % Update elapsed time if(n<C) n = n+1; % If queue isnt full, add customer if (n==1), % If only one customer, generate service time y = -log(rand)/rs; td = t + y; % Next departure time S = [S, y]; % Update service time record end A = [A, t]; % Update arrival time record end y = -log(rand)/ra; % Next interarrival time ta = t + y; % Next arrival time end % Case 2: Departure time has come if(td < ta), t = td; n = n-1; D = [D, t]; if (n==0), td = inf; else y = -log(rand)/rs; td = t + y; S = [S, y]; end end end
% % % % %
Update elapsed time Remove customer from queue Update departure time record If queue is empty set td to "infinity"
% Next service time % Next departure time % Update service time record
% Case 3: End of simulation, no new arrivals, flush queue while(n > 0) t = td; % Update elapsed time n = n-1; % Remove customer from queue if(n > 0) % If queue isnt empty y = -log(rand)/rs; % Next service time td = t + y; % Next departure time S = [S, y]; % Update service time record end D = [D, t]; % Update departure time record end Tp = max(t-T, 0); % Time past T that last customer leaves queue
The great advantage of discrete event simulation is exibility, and the power to examine systems that do not have nice analytical models. For example,
19. Discrete Event Simulation
19.6
Engs 27 Discrete and Probabilistic Systems Arrivals and departures can be made to follow any distribution, not just Poisson. Arrival and service rates can be time-varying. Servers can be connected in sequence (e.g., wait for your food, wait to pay for it). Post-simulation analysis During the simulation, we keep track of each customers arrival, departure, and service times. These are absolute times. Interarrival times are calculated by dierencing the arrival vector,
tA = diff(A);
Winter, 2004
The time each customer spends in the queue is the dierence of his departure and arrival time, and the time each customer spends waiting for service is the time in queue minus the service time,
tQ = D-A; tW = tQ-S;
Means and standard deviations can be computed from tA, tQ and tW using the mean and std functions. A running count of the queue size, X(t), can be reconstructed from the arrival and departure times. An approximate way to do this is to divide the simulation time [0, Ts ] into intervals of equal length, t. Then, for each time tm = mt, check to see how many arrivals and departures have occurred. The number of customers X(tm ) is the number of arrivals before tm minus the number of departures before tm . Equivalently, suppose that a number of customers have arrived before time tm . Some of these customers will also have departed by time tm . The length of the queue is the number of customers which havenot departed by time tm , those whose departure times are greater than tm . Here is a way to do the calculation with Matlab.
% Number of customers in system, vs t % Divide the simulation time into intervals of equal length. Number of % customers in system at time t = number of customers for which D > t > A. % You have to rebin the customers in this way to get correct statistics. dt = 0.1/rs; % 10% of the average service time t = 0:dt:(T+Tp); M = leng...