M/M/c Multi-Server Queue Simulation

What the simulation shows

Arrivals form a Poisson process with rate λ. An arriving customer takes an idle server if there is one and otherwise joins a single first-come, first-served line. Each server works at rate μ, so the pool can serve up to cμ customers per unit time. Two numbers describe the load: the offered load a = λ/μ, in erlangs, is the average number of busy servers; the per-server load ρ = a/c is the fraction of time each one is busy.

• Waiting is a yes-or-no event first. A customer waits only if all c servers are busy, which on the chart means N(t) is above the green line. The lab counts how often that happens and compares it with the Erlang C probability.
• Bigger pools wait less at the same utilization. The right chart holds ρ fixed and varies c. Separate lines (orange) wait as long as a single M/M/1 queue no matter how many there are; the shared line (blue) waits less with every server added.
• Stability depends on total capacity. The queue settles only when λ < cμ. At or above it, the theory readouts go blank and the line grows along the drift line.

Try c = 1 with λ = 0.8, then c = 10 with λ = 8: every server is equally busy in both, but the waits are very different.

The Erlang C formula

With exponential services the number in system is a birth–death process: it moves up at rate λ and down at rate min(n, c)·μ. Solving for its stationary distribution and adding up the states where all servers are busy gives the probability that an arriving customer has to wait (by the PASTA property, Poisson arrivals see time averages):

The powers and factorials overflow for large c, so the code uses the Erlang B recursion B(0) = 1, B(k) = aB(k−1)/(k + aB(k−1)) and converts with C = cB/(c − a(1 − B)). Every intermediate value stays between 0 and 1. With c = 1 the formula reduces to the M/M/1 queue: C = ρ.

Why one shared line beats separate lines

Split the same arrivals evenly at random over c separate servers, each with its own line. Each becomes an M/M/1 queue with arrival rate λ/c and the same load ρ, so its mean wait is ρ/(μ(1 − ρ)) whatever c is. With one shared line the mean wait is C(c, a)/(cμ − λ); the ratio of the two works out to a/C(c, a).

At ρ = 0.8 that ratio is 2.25 for 2 servers, 7.22 for 5, 19.55 for 10 and 62.48 for 20. The mechanism is simple: with separate lines a customer can be stuck behind a long service while a server in the next line stands idle. A shared line never wastes capacity that way, and the more servers share it, the more likely one of them is about to finish.

The same effect explains why large call centers can run their agents at higher utilization than small ones for the same service level. Customers who join the shortest line do better than random splitting; a single line simply guarantees that nobody waits while a server is free.

Example result

With λ = 4, μ = 1 and c = 5 servers (ρ = 0.800, offered load 4.00 erlangs), seed 42 and 20,000 customers of which the first 1,000 are discarded, 54.7% of customers had to wait, against the Erlang C value of 55.4%. The mean wait was 0.534 (theory 0.554) and the mean time in system 1.531 (theory 1.554). Serving the same arrivals with five separate M/M/1 lines would raise the mean wait to 4.000, 7.22 times longer.

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.

// Erlang C: probability that an arriving customer must wait in an M/M/c queue with
// offered load a = λ/μ < c. Computed through the Erlang B recursion, which never forms
// the large factorials and powers of the textbook formula:
//   B(0) = 1,  B(k) = a·B(k−1) / (k + a·B(k−1)),  C = c·B / (c − a·(1 − B)).
export function erlangC(c, a) {
  let b = 1;
  for (let k = 1; k <= c; k++) b = (a * b) / (k + a * b);
  return (c * b) / (c - a * (1 - b));
}

// Steady-state M/M/c results (null when λ ≥ cμ), plus c separate M/M/1 lines at λ/c each.
export function theory(lambda, mu, c) {
  const a = lambda / mu, rho = a / c;
  if (rho >= 1) return { a, rho, stable: false, pWait: null, Wq: null, W: null, L: null, WqSep: null, gain: null };
  const pWait = erlangC(c, a);
  const Wq = pWait / (c * mu - lambda);
  const W = Wq + 1 / mu;
  const WqSep = rho / (mu * (1 - rho));
  return { a, rho, stable: true, pWait, Wq, W, L: lambda * W, WqSep, gain: WqSep / Wq };
}

def erlang_c(c, a):
    """P(wait) in M/M/c with offered load a = lambda/mu < c, via the Erlang B recursion."""
    b = 1.0
    for k in range(1, c + 1):
        b = a * b / (k + a * b)
    return c * b / (c - a * (1 - b))


def theory(lam, mu, c):
    """Steady-state M/M/c results (None when lambda >= c*mu) and c separate M/M/1 lines."""
    a = lam / mu
    rho = a / c
    if rho >= 1:
        return {"a": a, "rho": rho, "stable": False, "p_wait": None, "Wq": None, "W": None, "L": None,
                "Wq_sep": None, "gain": None}
    p_wait = erlang_c(c, a)
    wq = p_wait / (c * mu - lam)
    w = wq + 1 / mu
    wq_sep = rho / (mu * (1 - rho))
    return {"a": a, "rho": rho, "stable": True, "p_wait": p_wait, "Wq": wq, "W": w, "L": lam * w,
            "Wq_sep": wq_sep, "gain": wq_sep / wq}

Run the full model from a terminal:

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

Download the model

Assumptions & limitations

• Poisson arrivals, exponential services. Independent and memoryless, with identical servers. Service-time variability matters: less variable services shorten waits.
• One line, first come, first served, infinite patience. No abandonment, no retrials, no priority classes. Call-center models often add abandonment (Erlang A), which caps the damage of overload.
• Server choice. An arriving customer takes the lowest-numbered idle server. That loads server 1 most (see the example CSV) but does not change any queue statistic.
• Separate-lines comparison is theoretical. It assumes arrivals are split at random, so each line is an independent M/M/1 queue; it is not simulated.
• Finite runs from empty. The first warm-up customers are discarded, but some initialization bias and sampling error remain, especially near ρ = 1.

Related simulations

• OS-018M/M/1 queueOne server, random arrivals, random service: why waits explode as utilization ρ = λ/μ approaches 1.
• 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

References

1. Kleinrock, L. (1975). Queueing Systems, Volume 1: Theory. Wiley.
2. 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
3. Law, A. M. (2015). Simulation Modeling and Analysis, 5th ed. McGraw-Hill.