Discrete-Event Simulation, Step by Step

How a discrete-event simulation works

Many systems change only at separate instants: a customer arrives, a service ends, a machine breaks. Between those instants nothing that matters happens, so there is no reason to advance time in small fixed steps. A discrete-event simulation (DES) has four parts:

• State — here, who is waiting at each station and who is being served.
• Clock — the current simulated time. It only ever jumps forward, to the time of the next event.
• Future event list — every event already scheduled, sorted by time. The lab draws it on a time axis and lists it in a table.
• Event routines — one per event type. Each changes the state and schedules the future events it causes.

The main loop is short:

1. Remove the first event from the future event list.
2. Add the time since the last event, weighted by the current state, to the statistics (busy time, area under the number in the shop).
3. Set the clock to the event’s time.
4. Run the event routine: update the state, record per-customer results, schedule new events.
5. Stop when the end condition holds; otherwise repeat.

Press Reset and then Step: the first event is always the first arrival, and the status line lists the events each step adds to the list.

The coffee-shop model

Customers arrive as a Poisson process with rate λ per minute, queue first-come, first-served at an order counter (exponential service, rate μ₁), then queue for a single barista (rate μ₂) and leave. Three event types are enough:

With these three event types the list never holds more than three events: one arrival and at most one completion per station. That is typical: the list stays small while the model’s history can be millions of events long. Larger models keep hundreds of pending events, which is why production simulators store the list in a heap or calendar queue rather than the sorted array used here.

The same model written in the process-interaction style, where each customer is a process that requests resources and the library keeps the event list, is in SimPy examples.

Checking the output against Jackson-network theory

A simulation is only trustworthy once it reproduces a case with a known answer. This one has one. The departures from a stable M/M/1 queue form a Poisson process with the same rate λ (Burke’s theorem), so the barista also sees Poisson arrivals, and Jackson showed in 1957 that in steady state an open network of such stations behaves as if each were an independent M/M/1 queue:

The theory needs every station stable, λ < μ₁ and λ < μ₂. The station with the highest utilization is the bottleneck: with the defaults the barista is busier than the counter and accounts for most of the waiting. Push λ above μ₂ and the barista’s line grows without bound while the counter still behaves; the lab then reports no theoretical values for the barista or the total.

Differences between simulation and theory come from finite run length, correlated waits and the empty start (see warm-up and run length). The Little’s law page shows a second check that holds on any single run.

Example result

With λ = 0.8 per minute, μ₁ = 1.2 and μ₂ = 1 (theoretical utilizations 66.7% and 80.0%), seed 42 and 20,000 customers with the first 1,000 discarded, the simulation processed 60,021 events over 24,996.6 simulated minutes. Customers waited 1.683 minutes on average at the counter (Jackson: 1.667) and 3.875 for the barista (4.000), and spent 7.384 minutes in the shop against the theoretical 7.500. The example CSV is the event log of the first 200 events.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// The event loop: remove the most imminent event, advance the clock to it, update the
// time-weighted statistics, change the state and schedule the events this one causes.
export function step(s) {
  const ev = s.fel.shift();
  const dt = ev.time - s.clock;
  if (s.measuring) {
    if (s.busy1) s.b1 += dt;
    if (s.busy2) s.b2 += dt;
    s.area += s.inShop * dt;
  }
  s.clock = ev.time;
  s.events++;
  const c = ev.cust, p = s.p, rng = s.rng;
  if (ev.type === 'arrival') {
    s.inShop++;
    c.arrive = s.clock;
    schedule(s, s.clock + rng.exponential(p.lambda), 'arrival', { id: s.nextId++ });
    if (!s.busy1) startService(s, 1, c); else s.q1.push(c);
  } else if (ev.type === 'end1') {
    s.busy1 = null;
    if (s.h1 < s.q1.length) startService(s, 1, s.q1[s.h1++]);
    c.arrive2 = s.clock; // walks straight to the barista's queue
    if (!s.busy2) startService(s, 2, c); else s.q2.push(c);
  } else {
    s.busy2 = null;
    if (s.h2 < s.q2.length) startService(s, 2, s.q2[s.h2++]);
    s.inShop--;
    s.departed++;
    if (s.departed > s.warmup) {
      s.counted++;
      s.sumW += s.clock - c.arrive;
      s.sumWq1 += c.start1 - c.arrive;
      s.sumWq2 += c.start2 - c.arrive2;
    }
    if (s.departed === s.warmup) { s.measuring = true; s.t0 = s.clock; }
  }
  compact(s);
  return ev;
}

function startService(s, station, c) {
  if (station === 1) {
    s.busy1 = c; c.start1 = s.clock;
    schedule(s, s.clock + s.rng.exponential(s.p.mu1), 'end1', c);
  } else {
    s.busy2 = c; c.start2 = s.clock;
    schedule(s, s.clock + s.rng.exponential(s.p.mu2), 'end2', c);
  }
}

def step(self):
    """Remove the most imminent event, advance the clock, update statistics and state."""
    time, kind, c = self.fel.pop(0)
    dt = time - self.clock
    if self.measuring:
        if self.busy1:
            self.b1 += dt
        if self.busy2:
            self.b2 += dt
        self.area += self.in_shop * dt
    self.clock = time
    self.events += 1
    p, rng = self.p, self.rng
    if kind == "arrival":
        self.in_shop += 1
        c.arrive = self.clock
        self.schedule(self.clock + rng.exponential(p["lambda"]), "arrival", self.new_customer())
        if not self.busy1:
            self.start(1, c)
        else:
            self.q1.append(c)
    elif kind == "end1":
        self.busy1 = None
        if self.h1 < len(self.q1):
            self.h1 += 1
            self.start(1, self.q1[self.h1 - 1])
        c.arrive2 = self.clock
        if not self.busy2:
            self.start(2, c)
        else:
            self.q2.append(c)
    else:
        self.busy2 = None
        if self.h2 < len(self.q2):
            self.h2 += 1
            self.start(2, self.q2[self.h2 - 1])
        self.in_shop -= 1
        self.departed += 1
        if self.departed > self.warmup:
            self.counted += 1
            self.sum_w += self.clock - c.arrive
            self.sum_wq1 += c.start1 - c.arrive
            self.sum_wq2 += c.start2 - c.arrive2
        if self.departed == self.warmup:
            self.measuring, self.t0 = True, self.clock
    self.compact()
    return time, kind, c

def start(self, station, c):
    if station == 1:
        self.busy1, c.start1 = c, self.clock
        self.schedule(self.clock + self.rng.exponential(self.p["mu1"]), "end1", c)
    else:
        self.busy2, c.start2 = c, self.clock
        self.schedule(self.clock + self.rng.exponential(self.p["mu2"]), "end2", c)

Run the full model from a terminal:

node simulation.js --lambda 0.8 --mu1 1.2 --seed 42
python simulation.py --lambda 0.8 --mu1 1.2 --seed 42

Download the model

Assumptions & limitations

• Memoryless arrivals and services. Poisson arrivals, exponential service at both stations, all independent. This is what makes the Jackson check exact; the simulation method itself works with any distributions.
• Single servers, unlimited lines, no abandonment. Everyone who orders waits for a drink; nobody leaves a long line.
• Deterministic tie-breaking. Events at exactly the same time run in the order they were scheduled. With continuous random times, ties have probability zero.
• Finite runs from an empty shop. The first warm-up customers are discarded; the remaining initialization bias and sampling error show up as differences from the theory.

Related simulations

• OS-018M/M/1 queueOne server, random arrivals, random service: why waits explode as utilization ρ = λ/μ approaches 1.
• OS-019M/M/c queueSeveral servers, one shared line: Erlang C waiting probability and the pooling effect.
• OS-020Little's lawAverage number in system = arrival rate × average time in system, read off one picture.
• OS-015Poisson processEvents at random times with a constant rate: Poisson counts, exponential gaps and the waiting-time paradox.

Background: SimPy Examples: Queues in Python

References

1. Banks, J., Carson, J. S., Nelson, B. L. & Nicol, D. M. (2010). Discrete-Event System Simulation, 5th ed. Pearson.
2. Law, A. M. (2015). Simulation Modeling and Analysis, 5th ed. McGraw-Hill.
3. Jackson, J. R. (1957). Networks of waiting lines. Operations Research, 5(4), 518–521. doi.org/10.1287/opre.5.4.518
4. SimPy documentation. simpy.readthedocs.io