M/M/1 Queue Simulation

What the simulation shows

Customers arrive one at a time, wait in a single first-in, first-out line and are served by one server. Arrivals form a Poisson process with rate λ; each service takes an exponential time with mean 1/μ. The ratio ρ = λ/μ, the traffic intensity, is the long-run fraction of time the server is busy.

• The line forms even though the server is fast enough on average. With the defaults the server could handle 25% more work than arrives, yet a queue is usually present. Randomness alone creates it: bursts of arrivals and long services pile up, and idle time lost during quiet spells is never recovered.
• The path is spiky and strongly correlated. The step chart of N(t) shows long excursions. The running average (orange) settles toward the dashed theoretical L only slowly.
• The number in system is geometric. After warm-up the share of time with n customers present follows (1 − ρ)ρⁿ: the server is idle a fraction 1 − ρ of the time and each extra customer is a factor ρ less likely.
• At ρ ≥ 1 there is no steady state. Set λ above μ: the chart climbs along the drift line (λ − μ)t, the theory readouts go blank and the averages just reflect how long you ran.

Press Play to watch the line from an empty start (the pace speeds up as time passes), or Step to process one arrival or departure at a time.

Kendall notation and the steady-state formulas

Kendall’s A/S/c notation names a queue by its arrival process, its service-time distribution and its number of servers. M stands for Markovian (memoryless): Poisson arrivals, exponential service times. So M/M/1 means Poisson arrivals, exponential services, one server; by convention the waiting room is unlimited and service is first come, first served. M/D/1 would have deterministic service, M/G/1 a general one, and M/M/c several servers sharing one line, as in the M/M/c multi-server queue.

Because both distributions are memoryless, N(t) is a birth–death process that moves up at rate λ and down at rate μ. For ρ < 1 its stationary distribution and means are

The four means are tied together by Little’s law, L = λW and L_q = λW_q, which holds for far more general queues. The factor 1/(1 − ρ) is the whole story of congestion: going from 80% to 90% utilization doubles the mean time in system, and 90% to 95% doubles it again.

When ρ > 1 the queue grows on average by λ − μ customers per unit time, forever. At exactly ρ = 1 it has no drift but still no steady state: it wanders like a random walk reflected at zero, reaching levels of order √t. The lab reports no theoretical values in either case rather than plugging into formulas that do not apply.

How the next-event simulation works

Nothing changes between arrivals and departures, so the simulation never ticks through time in fixed steps. It keeps a clock, the number in system, the arrival times of the customers in line, the time of the next arrival and the time the current service ends, and jumps straight from one event to the next:

1. Pick the earlier of the next arrival and the next departure; advance the clock to it.
2. Add N × (time elapsed) to the area under N(t), and the matching amounts to queue area and busy time.
3. Arrival: schedule the next arrival an exponential time later; start service at once if the server is idle, otherwise join the line.
4. Departure: record the customer’s time in system and time waiting; move the first customer in line into service and draw its service time.
5. Stop when the requested number of customers has left; divide the areas by the measured time to get L, L_q and utilization.

The discrete-event simulation walkthrough shows the same method with an explicit future event list and an event log, and SimPy examples writes queues like this one as Python processes.

Warm-up, run length and why ρ → 1 is hard

The theory describes a queue in steady state, but every run starts empty. Early on, N(t) is systematically below its long-run level, so statistics that include the start are biased low. This initialization bias is why the lab discards the first customers (shaded on the chart). Warm-up deletion trades a little data for less bias; it cannot fix a run that is simply too short.

The harder problem is autocorrelation. A customer who waited a long time is usually followed by others who wait a long time, because they are stuck behind the same backlog. Successive observations carry much less information than independent ones would, and the effect grows quickly as the server fills up. For the M/M/1 queue, measuring time in units of the mean service time, the variance of the time-average number in system over a run of length T is approximately

Pinning L down to ±1% (one standard error) therefore takes about 240,000 time units at ρ = 0.5 but about 4.2 million at ρ = 0.9. The lab estimates the uncertainty of W by batch means: it splits the measured customers into 20 consecutive batches and uses the spread of the batch averages. Compare that standard error with the gap between simulation and theory at ρ = 0.8 and at ρ = 0.97.

Example result

With λ = 0.8 and μ = 1 (ρ = 0.800), seed 42 and 20,000 customers of which the first 1,000 are discarded, the server was busy 79.8% of the time. The simulated mean number in system is L = 4.141 against the theoretical 4.000, and the mean time in system is W = 5.171 (batch-means standard error 0.409) against 5.000. The run covered 24,960.7 time units, and at no point were more than 34 customers present.

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.

// Next-event time advance: jump the clock to the earlier of the next arrival and the
// next departure, add the time-weighted statistics for the interval just skipped,
// then update the state. Arrivals and services are exponential (the two Ms).
export function step(s) {
  const departure = s.service ? s.service.end : Infinity;
  const isArrival = s.nextArrival < departure;
  const te = isArrival ? s.nextArrival : departure;
  if (s.measuring) {
    const dt = te - s.t;
    s.area += s.n * dt; // ∫ N(t) dt  → L
    if (s.n > 1) s.areaQ += (s.n - 1) * dt; // ∫ Nq(t) dt → Lq
    if (s.n > 0) s.busy += dt; // server busy time → utilization
    s.timeAt[s.n] += dt; // time spent with n in system
  }
  s.t = te;
  if (isArrival) {
    s.arrivals++;
    s.n++;
    s.nextArrival = te + s.rng.exponential(s.lambda);
    if (s.service === null) s.service = { id: s.arrivals, arrival: te, start: te, end: te + s.rng.exponential(s.mu) };
    else s.waiting.push(te);
  } else {
    const c = s.service;
    s.departures++;
    s.n--;
    if (s.departures > s.warmup) record(s, te - c.arrival, c.start - c.arrival);
    if (s.departures === s.warmup) { s.measuring = true; s.t0 = te; }
    if (s.head < s.waiting.length) {
      const arrival = dequeue(s);
      s.service = { id: c.id + 1, arrival, start: te, end: te + s.rng.exponential(s.mu) };
    } else s.service = null;
  }
  if (s.n > s.maxN) { s.maxN = s.n; s.timeAt.push(0); }
  return isArrival ? 'arrival' : 'departure';
}

def step(self):
    """Jump to the next event (arrival or departure), accumulate statistics, update the state."""
    departure = self.service[2] if self.service else math.inf
    is_arrival = self.next_arrival < departure
    te = self.next_arrival if is_arrival else departure
    if self.measuring:
        dt = te - self.t
        self.area += self.n * dt
        if self.n > 1:
            self.area_q += (self.n - 1) * dt
        if self.n > 0:
            self.busy += dt
        self.time_at[self.n] += dt
    self.t = te
    if is_arrival:
        self.arrivals += 1
        self.n += 1
        self.next_arrival = te + self.rng.exponential(self.lam)
        if self.service is None:
            self.service = (te, te, te + self.rng.exponential(self.mu))
        else:
            self.waiting.append(te)
    else:
        arrival, start, _ = self.service
        self.departures += 1
        self.n -= 1
        if self.departures > self.warmup:
            self.record(te - arrival, start - arrival)
        if self.departures == self.warmup:
            self.measuring, self.t0 = True, te
        if self.head < len(self.waiting):
            a = self.waiting[self.head]
            self.head += 1
            self.service = (a, te, te + self.rng.exponential(self.mu))
        else:
            self.service = None
    if self.n > self.max_n:
        self.max_n = self.n
        self.time_at.append(0.0)

Run the full model from a terminal:

node simulation.js --lambda 0.8 --mu 1 --seed 42
python simulation.py --lambda 0.8 --mu 1 --seed 42

Download the model

Assumptions & limitations

• Memoryless arrivals and services. Inter-arrival and service times are independent exponentials. Real service times are usually less variable than exponential, which lowers waits; Little’s law compares several service distributions through the Pollaczek–Khinchine formula.
• One server, one line, infinite patience. First come, first served; unlimited waiting room; nobody balks at a long line or leaves before service.
• Steady-state theory, finite-run simulation. Each run starts empty and stops after a fixed number of departures. Warm-up deletion reduces initialization bias but does not remove it.
• Approximate error bars. The batch-means standard error assumes the 20 batch averages are nearly independent. Near ρ = 1, short runs make them correlated and the standard error too small.
• Abstract time units. Rates are per unit time; read the unit as minutes, seconds or hours as long as λ and μ use the same one.

Related simulations

• 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-021Discrete-event simulationClock, future event list, event routines: the method behind queue simulators, one event at a time.
• 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 · How Accurate Is a Monte Carlo Estimate?

References

1. Kendall, D. G. (1953). Stochastic processes occurring in the theory of queues and their analysis by the method of the imbedded Markov chain. Annals of Mathematical Statistics, 24(3), 338–354. doi.org/10.1214/aoms/1177728975
2. Kleinrock, L. (1975). Queueing Systems, Volume 1: Theory. Wiley.
3. Law, A. M. (2015). Simulation Modeling and Analysis, 5th ed. McGraw-Hill.
4. Little, J. D. C. (1961). A proof for the queuing formula: L = λW. Operations Research, 9(3), 383–387. doi.org/10.1287/opre.9.3.383