Birth-Death Processes and the Queue Length Equilibrium

A single-server queue is the smallest model that exhibits the behaviour a discrete-event simulator exists to measure: the state is a count, arrivals and departures are the only events, and the stationary distribution has a closed form that the simulation can be checked against. The birth-death structure is what makes both the analysis and the verification possible, and it is the pattern that generalizes to every queueing network a simulator might build.

The state space and its two transition types

The state is the number of customers in the system, a count that increases by one on an arrival and decreases by one on a service completion. When the state is $n$, the arrival rate is $\lambda_n$ and the service rate is $\mu_n$, and the process moves from $n$ to $n+1$ at rate $\lambda_n$ and from $n$ to $n-1$ at rate $\mu_n$ (with $\mu_0 = 0$). Only neighbours are connected, which is the defining property of a birth-death process, and it is why the stationary distribution is solvable by a one-step recursion rather than by solving the full linear system.

The transition-rate diagram is a path. Each edge carries the rate of the corresponding event, and the flow across any cut in the path must balance at stationarity: the rate of upward crossings from below $n$ to above $n$ must equal the rate of downward crossings. This gives the local balance equation

$$\lambda_{n-1} \pi_{n-1} = \mu_n \pi_n,$$

which is the detailed balance condition specialized to a one-dimensional state space, and it holds because the only way to leave the set ${0, 1, \dots, n-1}$ is upward across the edge $(n-1, n)$ and the only way to enter it is downward across the same edge.

The geometric tail of the M/M/1 queue

When arrivals are Poisson with rate $\lambda$ and service times are exponential with rate $\mu$, independent of the state, the rates are $\lambda_n = \lambda$ and $\mu_n = \mu$ for $n \ge 1$. The recursion is $\pi_n = \rho \pi_{n-1}$ with $\rho = \lambda/\mu$, so

$$\pi_n = (1 - \rho)\rho^n,$$

a geometric distribution, provided $\rho < 1$. The load $\rho$ is the fraction of time the server is busy, and the condition $\rho < 1$ is the stability requirement: at $\rho \ge 1$ the arrival rate matches or exceeds the service rate and the queue grows without bound, so no stationary distribution exists.

The geometric tail is the source of every striking property of the queue. The mean number in the system is $\rho/(1-\rho)$, which diverges as $\rho \to 1$. At $\rho = 0.9$ the mean is 9; at $\rho = 0.99$ it is 99. The queue is not a linear function of the load but a ratio that blows up as the system approaches saturation, and a simulator that runs at $\rho = 0.99$ is measuring a system whose typical state is two orders of magnitude above the $\rho = 0.9$ case.

The variance matters more than the mean for simulation. The number in the system has variance $\rho/(1-\rho)^2$, which grows as the square of the mean, so the distribution is highly skewed near saturation. A steady-state estimator that averages the queue length over a run samples a distribution whose standard deviation is comparable to its mean, and the time to reach a stable estimate grows with the variance of the sampled quantity.

Service-time distribution changes the tail, not the existence

The exponential service assumption is what makes the Markov property hold on the count alone. With general independent service times the queue is the M/G/1 model, and the mean queue length is given by the Pollaczek-Khinchine formula

$$L = \rho + \frac{\rho^2(1 + C_s^2)}{2(1-\rho)},$$

where $C_s$ is the coefficient of variation of the service-time distribution. The mean depends on the service distribution only through its first two moments, which is a rare instance of a queueing quantity being moment-determined, and it shows that variability is a first-class input: doubling $C_s^2$ raises the mean queue length at the same load.

The practical reading is that a simulator must not replace a random service time with its mean. A deterministic service time has $C_s = 0$ and the smallest possible queue for a given mean; an exponential has $C_s = 1$; a heavy-tailed service distribution can have an enormous or even undefined variance, in which case the Pollaczek-Khinchine formula itself fails and the queue length is no longer characterized by its mean. A simulation that assumes exponential services measures a different system from one that models the true service distribution, and the difference is invisible in the mean service time.

Transience versus stationarity in a discrete-event run

A discrete-event simulator produces a trajectory, not a distribution, and the trajectory starts empty. The time until the queue length reaches something like its stationary level is the warm-up period, and observations taken during it are biased downward because the system has not yet filled. The bias is not small: at $\rho = 0.9$ the mean queue length is 9, and a run that is terminated too early averages a smaller number and reports a lower load than the input parameters imply.

The distinction is between the two kinds of average a run can report. A time average of the queue length over the simulated horizon converges to the stationary mean as the horizon grows, but a single long run starting from empty is one path, and its time average is a random quantity with a variance that depends on the autocorrelation of the process, not just on the length. Independent replications, each with its own warm-up discarded, produce an ensemble average whose standard error is the familiar $\sigma/\sqrt{n}$, at the cost of throwing away the warm-up from every replication.

import math, random

def mm1_queue(arrival, service, horizon):
    """Simulate an M/M/1 queue over a fixed time horizon by next-event
    scheduling; returns (mean queue length, utilization)."""
    t = 0.0
    n = 0
    integral = 0.0
    busy = 0.0
    next_arrival = random.expovariate(arrival)
    next_departure = float('inf') if n == 0 else t + random.expovariate(service)
    while t < horizon:
        if next_arrival < next_departure:
            dt = next_arrival - t
            integral += n * dt
            busy += dt if n > 0 else 0.0
            t = next_arrival
            n += 1
            next_arrival = t + random.expovariate(arrival)
            if n == 1:
                next_departure = t + random.expovariate(service)
        else:
            dt = next_departure - t
            integral += n * dt
            busy += dt
            t = next_departure
            n -= 1
            next_departure = t + random.expovariate(service) if n > 0 else float('inf')
    return integral / t, busy / t

The scheduler holds at most two pending events, the next arrival and the next departure, and advances time to whichever comes first. The queue-length integral accumulated over the horizon is what the time average divides by, and the busy fraction is the estimate of $\rho$. This is the next-event mechanism in miniature: the only state change occurs at event times, and nothing is computed in between.

State-dependent rates keep the structure and break the closed form

The birth-death recursion survives any dependence of the rates on the state; only the geometric solution depends on constant rates. When arrivals are discouraged by a long queue, as with balking, the arrival rate $\lambda_n$ falls with $n$, and when impatient customers leave, the service rate $\mu_n$ rises with $n$ as the effective departure rate includes reneging. The stationary probabilities are still

$$\pi_n = \pi_0 \prod_{k=1}^{n} \frac{\lambda_{k-1}}{\mu_k},$$

with $\pi_0$ fixed by normalization, but the product has no closed form for an arbitrary rate function, and the simulator is the tool that evaluates it.

This is the point where simulation stops competing with analysis and becomes the only option. A finite-capacity queue with balking at a threshold, reneging after an exponentially distributed patience, and a server that changes speed under load produces a stationary distribution that a birth-death recursion can describe but not summarize, and the quantity of interest, the abandonment rate or the fraction of customers who wait beyond a deadline, is a function of the transient behaviour as much as the equilibrium. The simulator measures it directly, and the closed-form M/M/1 cases are the special limits that verify the code.

What the stationary distribution is for in a simulator

The closed form is a verification target before it is an answer. A correctly implemented simulator of an M/M/1 queue should produce a time average near $\rho/(1-\rho)$ at low load, and a utilization near $\rho$, and a departure process that is Poisson with rate $\lambda$ in stationarity. A simulator that fails any of these checks has a bug in the event scheduling, the random variate generation, or the accounting, and the failure is visible precisely because the theory gives the expected value in closed form.

The same closed form is the basis of the inverse problem. Given a measured queue length, the implied load can be recovered by inverting $L = \rho/(1-\rho)$, and discrepancies between the implied load and the configured arrival and service rates point to a mismatch between the model and the system being modeled. A simulator is calibrated against a known queue, then trusted on the queue that has no closed form.

The limits of the birth-death model are the limits of the one-step structure. A state space with feedback, where a served customer can immediately re-enter the queue, still has a birth-death structure if the re-entry is folded into the arrival rate; a network of two queues has a state space that is a product of counts, and the flow balance becomes a set of coupled equations rather than a single recursion. The detailed balance condition is the general statement of the balance that the single queue satisfies in one line, and the synchronization of common random numbers is the technique that lets two queue configurations be compared with fewer replications than either could be measured alone.