Software 14.5.1 (M/M/1). The M/M/1 formulas for \(\lambda = 4\), \(\mu = 5\), and a simulation of \(400\,000\) customers by Lindley’s recursion \(W_{n + 1} = \max (0, W_n + S_n - A_{n + 1})\), which gives the wait of each customer from that of the one before:

>>> import numpy as np
>>> from math import factorial
>>> lam, mu = 4.0, 5.0                     # arrivals and services an hour
>>> rho = lam / mu
>>> L, W = rho / (1 - rho), 1 / (mu - lam)
>>> Wq = rho / (mu - lam); Lq = lam * Wq
>>> rho, L, W, Wq, Lq                      # Little: L = lam W, Lq = lam Wq
(0.8, 4.000000000000001, 1.0, 0.8, 3.2)
>>> rng = np.random.default_rng(0)
>>> def fifo_waits(inter, serv):           # Lindley's recursion for the waits
...     w = np.zeros(len(serv))
...     for n in range(1, len(serv)):
...         w[n] = max(0.0, w[n - 1] + serv[n - 1] - inter[n])
...     return w
...
>>> n = 400000
>>> inter = rng.exponential(1 / lam, n); serv = rng.exponential(1 / mu, n)
>>> wq = fifo_waits(inter, serv)[n // 10:]   # discard the warm-up
>>> round(wq.mean(), 4), round((wq > 0).mean(), 4)   # Wq, P(wait) = rho
(0.7943, 0.7987)

Reading the output. \(\rho = 0.8\), \(L = 4\), \(W = 1\), \(W_q = 0.8\) hours, \(L_q = 3.2\). The simulated mean wait, after discarding the first tenth as warm-up, is \(0.794\) hours, and the fraction of customers who wait is \(0.799\), close to \(\rho \). The simulation needs many customers because successive waits are strongly correlated when \(\rho \) is high.

MATLAB/Octave: the same recursion in a loop.

Software 14.5.2 (Erlang C and pooling). The Erlang C formula, the probability of waiting and mean wait for \(3\), \(4\) and \(5\) tellers in Example 14.2.2, and the wait with three separate M/M/1 queues against one pooled queue:

>>> def erlang_c(c, lam, mu):              # P(an arrival waits) in M/M/c
...     a = lam / mu; r = a / c
...     top = a**c / factorial(c) / (1 - r)
...     return top / (sum(a**k / factorial(k) for k in range(c)) + top)
...
>>> lam, mu = 24.0, 10.0                   # a bank: 24 customers an hour
>>> for c in (3, 4, 5):
...     pw = erlang_c(c, lam, mu); wq = pw / (c * mu - lam)
...     print(c, round(pw, 4), round(60 * wq, 2), 'min')
...
3 0.6472 6.47 min
4 0.287 1.08 min
5 0.1135 0.26 min
>>> pw1 = erlang_c(1, 8.0, 10.0)           # three separate tellers, 8 each
>>> pooled = erlang_c(3, 24, 10)           # one queue for all three
>>> round(60 * pw1 / (10 - 8), 2), round(60 * pooled / (30 - 24), 2)
(24.0, 6.47)

Reading the output. With \(3\), \(4\), \(5\) tellers the probability of waiting is \(0.647\), \(0.287\), \(0.114\) and the mean wait \(6.47\), \(1.08\), \(0.26\) minutes; the fourth teller cuts the wait by a factor of six. Three separate queues give a mean wait of \(24\) minutes against \(6.47\) for one shared queue with the same tellers.

MATLAB/Octave: the same function, with factorial.

Software 14.5.3 (The Pollaczek–Khinchine formula). \(L_q\) from the formula and from simulation for one server at \(\rho = 0.8\) with deterministic, exponential and gamma service times of squared coefficient of variation \(0\), \(1\) and \(4\):

>>> lam, mu = 4.0, 5.0; rho = lam / mu
>>> def pk_lq(cv2):                  # Pollaczek-Khinchine, Var S = cv2 / mu^2
...     return (lam**2 * cv2 / mu**2 + rho**2) / (2 * (1 - rho))
...
>>> laws = (('M/D/1', 0.0, lambda k: np.full(k, 1 / mu)),
...         ('M/M/1', 1.0, lambda k: rng.exponential(1 / mu, k)),
...         ('M/G/1, cv2 = 4', 4.0, lambda k: rng.gamma(0.25, 4 / mu, k)))
>>> for name, cv2, gen in laws:
...     s = gen(n); wq = fifo_waits(rng.exponential(1 / lam, n), s)[n // 10:]
...     print(f"{name:15s} Lq: formula {pk_lq(cv2):.3f}, "
...           f"simulated {lam * wq.mean():.3f}")
...
M/D/1           Lq: formula 1.600, simulated 1.576
M/M/1           Lq: formula 3.200, simulated 3.204
M/G/1, cv2 = 4  Lq: formula 8.000, simulated 7.702

Reading the output. The formula gives \(1.6\), \(3.2\) and \(8.0\), in the ratio \(1 : 2 : 5\) of \((1 + c_s^2)/2\), and the simulations agree, the last (\(7.70\)) less closely: with very variable services, rare long jobs dominate the average, and the simulation converges slowly.

MATLAB/Octave: gamrnd(0.25, 4/mu, 1, n) for the variable service times.