Discrete-Event Simulation

Dr. Phillip M. Feldman

1. Introduction

1.1 Continuous-time and Discrete-Event/Discrete-Time Simulation

Most dynamical system models fall into one of two general categories:
1. continuous-time: Changes of the system state occur at every moment of time. Examples of systems and phenomena that would be modeled as continuous-time include the following:

• heat flow
• orbital motion
• analog circuits

The evolution of the state of such a system is typically described by one or more differential equations.

2. discrete-event/discrete-time: Changes of system state occur at isolated instants in time. Between these event times, either the system state does not change, or it can be viewed as changing in a simple, predictable fashion. Examples of systems that would be modeled as discrete time include:

• any type of queueing system, e.g.:

• a network of computers that exchange packets
• a fast-food restaurant
• an inventory control system. A facility such as a supermarket that sells items and has a fixed warehouse capacity must decide when to reorder items and in what quantities.

• a digital circuit

• a system of hard particles that move and collide with one another. (Because the particles are hard, collisions can be modeled as instantaneous).

• warfare, e.g., an engagement between two forces, each consisting of infantry and tanks.

Combined continuous-time and discrete-event models are also possible. In such a model, some elements of the state evolve continuously while others change at discrete instants.

1.2 Analysis Versus Simulation

In application domains such as queueing, certain simple problems can be solved analytically, i.e., using mathematical methods rather than simulation. For example, for the M/G/1 queue, one can calculate the mean queue length and mean system time. But, obtaining other performance measures such as the full distribution of system time requires complex numerical calculations. Furthermore, seemingly minor changes to the model, e.g., changing the number of servers from 1 to 2, may make analytical calculation of even the most basic performance measures impractical. An analytical device such as embedded Markov chains that is applicable to one model may be inapplicable to a slightly different model.

An important advantage of discrete-event simulation (DES) is that it permits one to answer practically any question about the behavior of a queueing system using methods that are largely independent of model details (such as the number of servers) and the performance measures that we seek to estimate. A small change in the model typically requires a small change to the simulation code. (This last is not always true).

Most discrete-event simulations draw values from random distributions to model such stochastic elements as arrivals, service times, and failures. A concomitant drawback of DES is that there is a tradeoff between (a) the number of repetitions and/or the simulated time interval and (b) the accuracy of performance measure predictions. (This can be thought of as a tradeoff between CPU time and accuracy). While a mathematical analysis of a queueing model may produce exact values for performance measures, a stochastic simulation can never produce exact results.

1.3 Sources of Inaccuracy in DES Results

There are four main sources of error in results obtained from discrete-event simulations:

1. errors and approximations in the model itself. In a network simulation, for example, if packets arrive at a router in groups but are modeled as arriving individually according to a Poisson process, model predictions will be wrong.
2. random variation. A simulation of a dynamical system model with stochastic elements must use random numbers. As a consequence, performance measure estimates based on a finite sample of simulation runs of finite duration will not be precisely repeated with different random numbers, and will differ from the ensemble averages. (To match the ensemble averages, one would need to perform either an infinite number of finite-duration runs or a single run of infinite duration).

3. initialization bias. Many stochastic models have a steady-state average behavior that does not depend on the initial conditions. Simulations of such models are often used to predict the steady-state behavior. Aside from the stochastic variation mentioned above, there is also the issue that results from runs of finite duration will tend to depend to some degree on the initial conditions in the sense that estimates will differ from the true steady-state values in a consistent fashion. (More precisely, averages taken over $n$ runs will not converge to the steady-state truth as $n$ tends to infinity). This effect is known as initialization bias. In queueing simulations, for example, it is common to start the simulation with all queues empty; this tends to produce estimates of waiting time and system time that are biased low. Performing longer runs tends to reduce the initialization bias.

4. model implementation errors. If there are bugs in the software, or if the software model is inconsistent with the mathematical model, all bets are off.

It is good practice (not always followed) to use statistical methods to assess the uncertainty in simulation predictions due to random variation. This is accomplished by generating a confidence interval for each performance measure of interest. The data that goes into the confidence interval calculation may be collected in either of two ways:

1. Repeat the same simulation with different sequences of random numbers.
2. Perform a single long simulation run and divide the simulated time into non-overlapping intervals.

The second method has the advantage that the impact of the initialization bias is incurred only once, rather than in every run. On the other hand, when one uses data from a single simulation run, there is serial dependence between successive time intervals; the statistical tools that are ordinarily used for calculating confidence intervals do not account for such dependence and may produce misleading results. More specialized methods of time series can be used to generate confidence intervals from data with serial dependence.

1.4 The Concept of an Event Queue

Most discrete-event simulations store information about pending events in a data structure known as an event queue. Each item in the queue would at minimum contain the following information: (a) a timestamp and (b) an integer or string that identifies the type of event. In a queueing simulation, event types might include 'customer arrival' and 'customer departure'. Events sometimes have associated parameters, e.g., a customer departure event might specify the number of the departing customer (to facilitate the collection of statistics).

Some simulations make use of special types of events that do not affect the state of the model. Such events may be used to collect statistical information or to terminate the simulation.

The typical operations on an event queue are: (a) inserting a new event and (b) removing the next event—the one with the lowest timestamp—from the queue. It may also be necessary to (c) cancel a scheduled event. Because operations (a) and (b) must be performed efficiently, the common implementation of an event queue is a partially-ordered heap. With such a data structure, obtaining the next event requires O(1) time, while inserting an event requires O(log(n)) time, where n is the number of events currently in the queue.

In a Python-based simulation, one can implement an event queue using Python's heapq module, which is an efficient priority queue implementation.

1.5 Using Python and SimPy for Discrete-Event Simulation

Writing a discrete-event simulation in Python does not require infrastructure beyond the Python Standard Library and NumPy. (NumPy is invaluable for arrays, for random number generation, and for the generation of random deviates from a wide variety of discrete and continuous distributions). But, the use of other infrastructure can simplify the implementation of a complex DES model in code.

There are three approaches to implementing a discrete-event simulation with a general-purpose programming language such as Python:

1. There are two types of models in which no special machinery is required for managing events: (1) Events occur at regular intervals. Such models are called synchronous. (2) All events can be generated at one shot, followed by processing of the data to collect results (model performance data). In the next section, I present an example of (2).

2. For models that require an event queue, one might choose to use Python's heapq module.

3. For complex models, e.g., with multiple types of shared resources, it is often convenient to build on a simulation framework. When creating a Python-based simulation, the best available discrete-event framework is SimPy.

The following is quoted from simpy.org:
• SimPy is a process-based discrete-event simulation framework based on standard Python. Its event dispatcher is based on Python’s generators and can also be used for asynchronous networking or to implement multi-agent systems (with both, simulated and real communication).
• Processes in SimPy are defined by Python generator functions and can, for example, be used to model active components like customers, vehicles or agents. SimPy also provides various types of shared resources to model limited capacity congestion points (like servers, checkout counters and tunnels).

In the remainder of this article, I present two extended discrete-event simulation examples. One of these is implemented using only Python and NumPy (with matplotlib for graphical output); the second uses SimPy.

2. Discrete-Event Simulation Example #1: Shot Noise

2.1 Overview

This Python program simulates shot noise, such as that associated with the flow of electrons across a reverse-biased PN junction, and plots the resulting data in several different ways. Shot noise is a type of electronic noise originating from the discrete nature of electric charge.

Our model is as follows:

• Electrons arrive at and begin crossing the junction at random times; these times can be modeled as a Poisson process with specified rate lambda. The Poisson assumption implies that the number of electrons arriving during any time interval is independent of the number arriving during any non-overlapping time interval.
• The time T_transit for an electron to transit the junction is modeled as a constant. In actuality, the transit time varies, but for purposes of an initial investigation, this is a reasonable assumption. [Allowing the transit time to have a distribution with a tail causes the autocorrelation to decay smoothly to zero, and eliminates the nulls in the power spectral density.]

The current in units of Amperes is given by the following equation:

I(t)= n(t) * q_e / T_transit

In the above, n(t) is the number of electrons in transit (moving through the junction) at time t, q_e is the charge of one electron, and T_transit is the time for an electron to pass through the junction. Thus, the current is proportional to the number of electrons in transit.

2.2 Objectives

The objectives of the simulation are to verify the following:

1. The distribution of the number of electrons in transit should be approximately normal (Gaussian).
2. The autocorrelation should have a triangular shape, diminishing to zero when the lag reaches T_transit.
3. The power spectral density should be nearly white up to frequencies on the order of 0.25 / T_transit.

2.3 Implementation Approach

Because electrons do not interact with one another in our model, we can generate all of the arrival and departure times at the outset and then calculate the number of electrons in transit at each event time. The use of a DES framework such as SimPy would be overkill for this model. In detail, the following sequence of calculations is performed:
1. Generate the total number of electrons that arrive over the course of the simulation and store this value as N_arrivals. This number is a Poisson random variable with mean $\lambda * T_{sim}$.
2. Generate the arrival times. Given the number of arrivals and the total simulation time, the arrival times are i.i.d. uniform on $[0, T_{sim}]$.

3. Add the transit time T_transit to each arrival time to obtain the corresponding departure time. (This is done as a single array operation to avoid slow, Python-level looping).

4. Combine the lists of arrival and departure times, add the time zero to the combined list, and sort the result to obtain a list of all event times of interest. These are stored in the array times.

5. Remove elements of times that are greater than or equal to T_sim.

6. For each time in times, calculate the number of electrons in transit, i.e., inside the junction. Note that this is the number of arrivals minus the number of departures.

2.5 Results

Figure A—the first figure generated by the program—illustrates the time history of the number of electrons in transit. From the graph, the noise appears to be 'white', i.e., uncorrelated over time, but this is in fact not exactly the case because of the finite time required for an electron to transit the junction.

Figure B shows a histogram of the number of electrons in (i.e., crossing) the junction. It looks very nearly Gaussian, but close examination reveals that it is skewed slightly to the right. (The number of electrons in the junction cannot be less than zero, but there is no upper limit).

Figure C shows the estimated autocorrelation of the number of electrons in the junction. Note that the autocorrelation initially slopes down linearly, reaching zero at the transit time. Wiggles to the right of t=1 represent statistical variation, which would be eliminated in the limit of infinite simulation time.

Figure D shows the power spectral density of the number of electrons in transit. The power spectral density was estimated using Welch's method, which applies a taper function to sections (windows) of the data, calculates a periodogram from each section, and averages the periodograms together. Two different window sizes—2048 and 8192—were used, as well as two different taper functions. Note that the estimates based on the larger window size do a better job of capturing the null in the PSD, but are also more sensitive to noise in the data.

2.6 Conclusions

The histogram confirms that the distribution of electrons in transit is approximately normal. (It becomes closer to normal as the mean number of electrons in transit increases).

The autocorrelation graph has the expected triangular shape, except for some noise that could be reduced by increasing T_sim.

The PSD is not as flat as expected at low frequencies. Probable causes are either insufficient data or limitations of the PSD estimation algorithm.

3. Discrete-Event Simulation Example #2: Automotive Traffic

3.1 Overview

This Python script simulates automotive traffic arriving at a traffic light. The script demonstrates the following:
1. several basic principles of discrete-event simulation, including the representation of system states, the connection between events and changes of state, and the Poisson process,
2. the use of the SimPy module,
3. simulation of random deviates using NumPy, and
4. the use of Python generator functions.

3.2 The Model

Cars drive on a single-lane road and arrive at the intersection from one direction only according to a Poisson process with specified rate. Positions and velocities of the cars are not modeled, and there is thus no coordinate system.

The traffic light switches between green and red at regular intervals; there is no yellow state.

When a car arrives at a green light with no cars queued, it passes immediately through the intersection and departs the simulation.

When a car arrives at a red light or with one or more cars queued at the intersection, it stops and joins the queue.

When the light turns green, if there is a queue, cars enter the intersection one at a time. Once a car has entered the intersection, it clears the intersection, regardless of the state of the light, after a delay that is a random draw from a triangular distribution.

3.3 Input and Output Interfaces

All inputs are hardwired into the code. (This is not good programming practice, but it has the advantage that everything is contained in a single file).

The simulation generates a trace of events as it runs, and also reports estimates of the average number of cars in the queue and the average waiting time.

3.4 Possible Extensions of the Model

This model might be extended in any of several directions. Some possibilities include the following:
1. There are several sources of bias in the average waiting time: (a) The simulation starts with a green light. (b) The simulation starts with no cars queued. (c) Cars remaining in the queue at the end of the simulation do not contribute to the estimated waiting time. The effects of (a) and (b) might be mitigated by eliminating the initial transient; this would involve collecting and saving individual waiting times, choosing alternative values for the duration of the initial transient, generating estimates using only the data outside of the transient interval, and then selecting one of these results using some statistical criteria. (We want the mean of the selected result to have no statistically significant difference from results based on longer transient intervals, but otherwise we want to discard as little of the data as possible).
2. The delay of the driver at the head of the queue in responding to the change of the light from red to green could be made somewhat longer than that of drivers who are further back in the queue. This would allow for greater realism in the model without explicitly modeling positions and velocities of the cars.

3. One might have some fraction of the cars make a right-turn at the light, with the remaining cars going straight, and compare the capacity of a single-lane road, a two-lane road with no turn restrictions, and a two-lane road with a dedicated right-turn-only lane.

3.6 Output from a Run

> python traffic_sim.py

Simulation of Cars Arriving at Intersection Controlled by a Traffic Light

The light turned green at time 0.000.
Car #1 arrived to a green light with no cars waiting at time 4.706.
Car #2 arrived to a green light with no cars waiting at time 15.595.
Car #3 arrived to a green light with no cars waiting at time 18.829.

The light turned red at time 30.000.
Car #4 arrived and joined the queue at position 1 at time 30.226.
Car #5 arrived and joined the queue at position 2 at time 30.376.
Car #6 arrived and joined the queue at position 3 at time 42.595.
Car #7 arrived and joined the queue at position 4 at time 47.117.
Car #8 arrived and joined the queue at position 5 at time 62.352.
Car #9 arrived and joined the queue at position 6 at time 66.757.

The light turned green at time 70.000.
Car #4 departed at time 72.010, leaving 5 cars in the queue.
Car #5 departed at time 74.265, leaving 4 cars in the queue.
Car #6 departed at time 76.161, leaving 3 cars in the queue.
Car #7 departed at time 77.782, leaving 2 cars in the queue.
Car #10 arrived and joined the queue at position 3 at time 78.785.
Car #8 departed at time 79.802, leaving 2 cars in the queue.
Car #11 arrived and joined the queue at position 3 at time 80.976.
Car #12 arrived and joined the queue at position 4 at time 81.616.
Car #9 departed at time 81.918, leaving 3 cars in the queue.
Car #10 departed at time 83.558, leaving 2 cars in the queue.
Car #11 departed at time 85.273, leaving 1 cars in the queue.
Car #13 arrived and joined the queue at position 2 at time 85.427.
Car #12 departed at time 87.101, leaving 1 cars in the queue.
Car #13 departed at time 88.765, leaving 0 cars in the queue.
Car #14 arrived to a green light with no cars waiting at time 89.862.
Car #15 arrived to a green light with no cars waiting at time 93.045.

The light turned red at time 100.000.
Car #16 arrived and joined the queue at position 1 at time 103.176.
Car #17 arrived and joined the queue at position 2 at time 106.408.
Car #18 arrived and joined the queue at position 3 at time 115.312.
Car #19 arrived and joined the queue at position 4 at time 125.277.
Car #20 arrived and joined the queue at position 5 at time 127.870.

The light turned green at time 140.000.
Car #16 departed at time 141.863, leaving 4 cars in the queue.
Car #21 arrived and joined the queue at position 5 at time 142.561.
Car #17 departed at time 143.754, leaving 4 cars in the queue.
Car #22 arrived and joined the queue at position 5 at time 145.814.
Car #18 departed at time 145.827, leaving 4 cars in the queue.
Car #19 departed at time 148.060, leaving 3 cars in the queue.
Car #20 departed at time 149.871, leaving 2 cars in the queue.
Car #21 departed at time 151.817, leaving 1 cars in the queue.
Car #22 departed at time 153.938, leaving 0 cars in the queue.
Car #23 arrived to a green light with no cars waiting at time 156.028.
Car #24 arrived to a green light with no cars waiting at time 157.470.
Car #25 arrived to a green light with no cars waiting at time 160.066.
Car #26 arrived to a green light with no cars waiting at time 164.971.
Car #27 arrived to a green light with no cars waiting at time 165.049.

The light turned red at time 170.000.
Car #28 arrived and joined the queue at position 1 at time 172.428.
Car #29 arrived and joined the queue at position 2 at time 175.255.
Car #30 arrived and joined the queue at position 3 at time 179.434.
Car #31 arrived and joined the queue at position 4 at time 184.672.
Car #32 arrived and joined the queue at position 5 at time 184.777.

*** Statistics ***

Mean number of cars waiting: 2.345
Mean waiting time (seconds): 13.671


4. References

1. Documentation on SimPy, including tutorials, can be found at the following URL: http://simpy.readthedocs.org/en/latest/index.html
2. For general background on thermal noise and shot noise:

http://www.fis.unipr.it/~gigi/dida/strumentazione/harvard_noise.pdf

3. For background on the Poisson process:

http://en.wikipedia.org/wiki/Poisson_process

4. Bandwidth of a transistor depends on junction transit time, which in turn depends on the materials, geometry, and temperature (higher temperature reduces the junction transit time). See, for example, the following: