9. Stochastic Models
Most models in the previous chapters take their data as known. This chapter makes randomness part of the model itself. It covers the probability tools OR relies on, Markov chains that describe how a random system evolves, queueing theory that predicts congestion, and stochastic programming that optimizes decisions taken before uncertainty is resolved. These models size call centers, hospital wards and server pools. They also set inventory and maintenance policies, and plan capacity when demand, yields or prices are random.
<a id="91-probability-review-for-or"></a>
9.1 Probability Review for OR
Stochastic models need only a small part of probability theory, but they use it constantly. This subchapter reviews random variables and the distributions that recur in OR, then the Poisson process and the exponential distribution, which make Markovian models tractable. It ends with conditional expectation, the main tool for computing expectations in the rest of the chapter.
Random Variables and Distributions
Why it matters. Every stochastic model starts with a choice of distribution for its random inputs: interarrival times, service times, demand, yields, lifetimes. The choice matters a great deal. Two distributions with the same mean can produce queues whose lengths differ by a factor of ten, or stock-out risks that differ by an order of magnitude. This topic reviews the vocabulary and the handful of distributions that cover most OR applications.
Theory. A random variable $X$ assigns a number to each outcome of a random experiment, such as the demand tomorrow or the duration of the next repair.
- A discrete random variable takes countably many values. It is described by its probability mass function (pmf) $p(k) = P(X = k)$.
- A continuous random variable is described by a density $f(x) \ge 0$, with $P(a \le X \le b) = \int_a^b f(x),dx$. Single points have probability zero.
- Both have a cumulative distribution function (cdf) $F(x) = P(X \le x)$. Its complement $\bar F(x) = 1 - F(x) = P(X > x)$, the survival function or tail, answers the questions OR cares about most: "how likely is a wait longer than 10 minutes?", "how likely is demand above my stock?".
The expectation (mean) is the probability-weighted average, $E[X] = \sum_k k,p(k)$ or $\int x f(x),dx$. It is linear, $E[aX + bY] = aE[X] + bE[Y]$, whether or not $X$ and $Y$ are independent. The variance measures spread around the mean:
$$ \mathrm{Var}(X) = E\big[(X - E[X])^2\big] = E[X^2] - E[X]^2, \qquad \mathrm{Var}(aX + b) = a^2,\mathrm{Var}(X). $$
For independent $X$ and $Y$, $\mathrm{Var}(X + Y) = \mathrm{Var}(X) + \mathrm{Var}(Y)$. Otherwise a covariance term $2,\mathrm{Cov}(X, Y)$ appears.
Variance has the squared unit of $X$, so it is hard to compare across problems. OR uses the dimensionless squared coefficient of variation
$$ c^2 = \frac{\mathrm{Var}(X)}{E[X]^2}. $$
Its benchmarks are $c^2 = 0$ for a deterministic quantity, $c^2 = 1$ for the exponential distribution, and $c^2 > 1$ for "highly variable" quantities. The queueing formulas of Section 9.3 depend on service and arrival variability only through $c^2$. An analogy is two commuting routes that both average 30 minutes. One always takes 30 minutes; the other takes anything from 15 to 60. The means are equal, but you leave home at a different time depending on which route you take.
The distributions below cover most OR models:
| Distribution | Type | Mean | Variance | Typical OR use |
|---|---|---|---|---|
| Bernoulli$(p)$ | discrete | $p$ | $p(1-p)$ | defective or not, show or no-show |
| Binomial$(n, p)$ | discrete | $np$ | $np(1-p)$ | defectives in a lot, passengers who show up |
| Geometric$(p)$, $k = 1, 2, \dots$ | discrete | $1/p$ | $(1-p)/p^2$ | attempts until the first success |
| Poisson$(\lambda)$ | discrete | $\lambda$ | $\lambda$ | arrivals in an interval, demand for spare parts |
| Uniform$(a, b)$ | continuous | $(a+b)/2$ | $(b-a)^2/12$ | only a range is known; random number generation |
| Exponential (rate $\lambda$) | continuous | $1/\lambda$ | $1/\lambda^2$ | interarrival times, lifetimes without aging |
| Erlang$(k, \lambda)$ | continuous | $k/\lambda$ | $k/\lambda^2$ | multi-phase service; less variable than exponential |
| Normal$(\mu, \sigma^2)$ | continuous | $\mu$ | $\sigma^2$ | aggregated demand, measurement errors |
| Lognormal, Weibull, gamma | continuous | - | - | repair times (skewed), lifetimes with aging |
Three sources guide the choice of a distribution:
- The mechanism. Counts of many independent rare events are Poisson. Sums of many small effects are approximately normal (central limit theorem). Products of many factors are lognormal. A sum of $k$ independent exponential phases is Erlang.
- Data. Fit candidate distributions and test the fit (input modeling, Section 12.2).
- Tractability. Exponential times make Markov models work (rest of this chapter). If the data say otherwise, check how sensitive the answer is to that assumption.
Example: five service-time models with the same mean. A barista's service time averages 4 minutes. The code compares five candidate models with that mean. The hyperexponential one mixes 90% quick orders with 10% long ones.
import numpy as np
rng = np.random.default_rng(seed=1)
n = 500_000
m = 4.0 # every model has a mean of 4 minutes
# name: (n sampled service times, theoretical variance)
models = {
"deterministic": (np.full(n, m), 0.0),
"uniform(2, 6)": (rng.uniform(2, 6, n), (6 - 2) ** 2 / 12),
"Erlang k=4": (rng.gamma(4, m / 4, n), m ** 2 / 4), # sum of 4 exponential phases
"exponential": (rng.exponential(m, n), m ** 2),
"hyperexp.": (np.where(rng.random(n) < 0.9, # 90% quick jobs (mean 2),
rng.exponential(2.0, n), # 10% long jobs (mean 22)
rng.exponential(22.0, n)), 88.0),
}
print(f"{'model':>14} {'mean':>6} {'var':>6} {'var(th)':>8} {'CV^2':>5} {'P(X>8)':>7} {'P(X>20)':>8}")
for name, (x, var_th) in models.items():
scv = x.var() / x.mean() ** 2 # squared coefficient of variation
print(f"{name:>14} {x.mean():6.2f} {x.var():6.2f} {var_th:8.2f} {scv:5.2f}"
f" {np.mean(x > 8):7.3f} {np.mean(x > 20):8.4f}")
# model mean var var(th) CV^2 P(X>8) P(X>20)
# deterministic 4.00 0.00 0.00 0.00 0.000 0.0000
# uniform(2, 6) 4.00 1.33 1.33 0.08 0.000 0.0000
# Erlang k=4 4.00 4.00 4.00 0.25 0.043 0.0000
# exponential 4.01 16.12 16.00 1.00 0.136 0.0069
# hyperexp. 4.00 87.57 88.00 5.48 0.086 0.0401
All five models agree on the mean, yet $c^2$ ranges from 0 to 5.5. Their tails tell different stories too. The hyperexponential model is less likely than the exponential to exceed 8 minutes, but six times more likely to exceed 20. A single summary number never captures the whole risk. Section 9.3 shows that, with Poisson arrivals, the mean queue grows in proportion to $1 + c^2$. With the same arrivals, the hyperexponential barista would have a queue more than three times longer than the exponential one.
Key takeaways
- A distribution is described by its pmf or density, its cdf, or its tail $P(X > x)$. Tails answer most service-level and risk questions.
- Expectation is always linear. Variances add only for independent (or uncorrelated) variables.
- The squared coefficient of variation $c^2 = \mathrm{Var}/\text{mean}^2$ is the standard dimensionless measure of variability: 0 deterministic, 1 exponential, above 1 highly variable.
- Choose distributions from the mechanism that generates the data, from data fitting, and with tractability in mind. Then check how sensitive the conclusions are to that choice.
Practice
- (Beginner) An airline sells 105 tickets for a 100-seat flight. Each passenger independently shows up with probability 0.92. Compute the mean and variance of the number who show up, and the probability that more than 100 show up (use
scipy.stats.binomor sum the pmf).- (Intermediate) Show that an Erlang distribution with $k$ phases and mean $m$ has $c^2 = 1/k$. A sample of repair times has mean 30 minutes and standard deviation 15 minutes: which Erlang $k$ matches the first two moments? What happens to the Erlang distribution as $k \to \infty$ with the mean fixed?
- (Interview) "Two suppliers have the same average lead time of 10 days. One is 5% cheaper. Why might the cheaper one cost more overall?" Hint: safety stock protects against demand during the lead time, so it depends on the lead time's variability, not only on its mean. Compare the two suppliers' coefficients of variation and tails.
Poisson Process
Why it matters. Arrivals of customers, calls, orders, failures and web requests drive every queueing model. When arrivals come from many independent sources, each of which contributes only rarely, the Poisson process is the natural model. The Palm-Khintchine theorem makes this precise: the superposition of many sparse independent streams is approximately Poisson. It is the "M" for arrivals in M/M/1, and its properties make many calculations of this chapter possible.
Theory. A counting process $N(t)$ counts the events in $(0, t]$. It is a Poisson process with rate $\lambda$ if
- $N(0) = 0$,
- the counts in disjoint intervals are independent (independent increments),
- the count in any interval of length $h$ is Poisson with mean $\lambda h$ (stationary increments):
$$ P\big(N(t + h) - N(t) = k\big) = e^{-\lambda h} \frac{(\lambda h)^k}{k!}, \qquad k = 0, 1, 2, \dots $$
Two equivalent definitions give useful intuition:
- Local view. In a tiny interval of length $h$, one event occurs with probability $\lambda h + o(h)$ and two or more with probability $o(h)$. Events come one at a time and the process has no memory of its past.
- Interarrival view. The times between successive events are independent exponential random variables with mean $1/\lambda$. This is how a Poisson process is simulated.
time ---|-------x-----x--x-----------x------x-------x---|--->
0 T1 T2 T3 T4 T5 T6 t
<-A1-><-A2-><A3><---A4----><-A5-><--A6->
A1, A2, ... i.i.d. Exponential(lambda); N(t) = 6; E[N(t)] = Var[N(t)] = lambda t
Properties used throughout the chapter:
| Property | Statement | Use |
|---|---|---|
| Mean equals variance | $E[N(t)] = \mathrm{Var}[N(t)] = \lambda t$ | quick data check: the dispersion index Var/mean should be near 1 |
| Superposition | independent Poisson($\lambda_1$) + Poisson($\lambda_2$) = Poisson($\lambda_1 + \lambda_2$) | merging customer streams |
| Splitting (thinning) | if each event is of type $i$ with probability $p_i$, independently, the type-$i$ events form independent Poisson($\lambda p_i$) processes | routing in queueing networks, customer classes |
| Conditional uniformity | given $N(t) = n$, the event times are distributed like $n$ sorted uniform points on $(0, t)$ | simulation, statistical tests |
| PASTA | Poisson arrivals see time averages: an arrival finds the system in its long-run distribution | queueing formulas (Section 9.3) |
The independence of split streams is surprising. Knowing that many sales calls arrived in an hour tells you nothing about the other calls in that hour.
Real arrival rates vary by time of day. A nonhomogeneous Poisson process has a rate $\lambda(t)$, and the count in $(a, b]$ is Poisson with mean $\int_a^b \lambda(t),dt$. Call centers model each half-hour as a Poisson process with its own constant rate and staff each interval separately (M/M/c topic).
A dispersion index well above 1 signals that the Poisson assumption fails. The cause may be batch arrivals (a bus unloading), a rate that changes inside the interval, or contagion (one failure triggering others).
Example: arrivals to a call center. Calls arrive at 12 per hour. The code builds the process from exponential gaps, checks the counts in 10-minute windows against the Poisson pmf, then splits off the sales calls (25%).
import numpy as np
from math import exp, factorial
rng = np.random.default_rng(seed=7)
lam = 12.0 # calls per hour
T = 20_000.0 # hours of simulated operation
# Build the process from i.i.d. exponential gaps (mean 1/lam hours)
gaps = rng.exponential(1 / lam, size=int(lam * T * 1.05))
times = np.cumsum(gaps)
times = times[times < T] # arrival epochs in [0, T)
# Counts in consecutive 10-minute windows: should be Poisson with mean lam/6 = 2
counts = np.histogram(times, bins=np.linspace(0, T, int(T * 6) + 1))[0]
print(f"10-min counts: mean {counts.mean():.3f}, variance {counts.var():.3f}")
for k in range(6):
pmf = exp(-2) * 2 ** k / factorial(k)
print(f" P(N={k}): simulated {np.mean(counts == k):.4f} Poisson(2) {pmf:.4f}")
# Splitting (thinning): each call is a sales call with probability 0.25
is_sales = rng.random(times.size) < 0.25
sales, other = times[is_sales], times[~is_sales]
g = np.diff(sales)
print(f"sales stream: rate {sales.size / T:.2f}/h, gap mean {60 * g.mean():.2f} min,"
f" gap CV {g.std() / g.mean():.3f}")
h_sales = np.histogram(sales, bins=int(T))[0] # hourly counts of each stream
h_other = np.histogram(other, bins=int(T))[0]
print(f"correlation of hourly counts (sales vs other): {np.corrcoef(h_sales, h_other)[0, 1]:+.4f}")
# 10-min counts: mean 2.002, variance 1.986
# P(N=0): simulated 0.1334 Poisson(2) 0.1353
# P(N=1): simulated 0.2707 Poisson(2) 0.2707
# P(N=2): simulated 0.2726 Poisson(2) 0.2707
# P(N=3): simulated 0.1821 Poisson(2) 0.1804
# P(N=4): simulated 0.0894 Poisson(2) 0.0902
# P(N=5): simulated 0.0352 Poisson(2) 0.0361
# sales stream: rate 3.01/h, gap mean 19.91 min, gap CV 0.998
# correlation of hourly counts (sales vs other): -0.0023
The counts have mean and variance both close to 2, and their frequencies match the Poisson pmf. The sales stream has rate $12 \times 0.25 = 3$ per hour, exponential-looking gaps (coefficient of variation 1), and hourly counts uncorrelated with the other calls. That is the splitting property at work.
Key takeaways
- A Poisson process has independent, stationary, Poisson-distributed increments. Equivalently, its interarrival times are i.i.d. exponential with mean $1/\lambda$.
- Mean and variance of the count are both $\lambda t$. A dispersion index far from 1 is evidence against the Poisson assumption.
- Merging independent Poisson streams gives a Poisson stream. Randomly splitting one gives independent Poisson streams.
- Time-varying rates are handled with a nonhomogeneous Poisson process, usually piecewise constant per interval.
- PASTA: Poisson arrivals see the system in its long-run (time-average) state.
Practice
- (Beginner) A help desk receives tickets as a Poisson process with rate 6 per hour. Compute the probability of no ticket in 20 minutes and of at least 3 tickets in 30 minutes. No ticket arrived in the last 15 minutes: what is the expected time until the next one?
- (Intermediate) In the example, condition on hours with exactly 4 sales calls and compute the mean and variance of the non-sales calls in those hours. Compare with Poisson(9) and explain the result.
- (Interview) "Our web logs show that the variance of requests per minute is five times the mean. Can we use Poisson arrivals for capacity planning?" Hint: distinguish a rate that varies over the day (fix: a nonhomogeneous model per interval) from genuine bursts or batches (fix: batch or compound arrivals). Recompute the dispersion inside short, homogeneous intervals.
Exponential Distribution and Memorylessness
Why it matters. The exponential distribution is what makes Markov models of queues and reliability tractable. If the remaining time of an activity does not depend on how long it has already lasted, then the current state alone describes the system, and the past can be forgotten. Knowing when this assumption is reasonable, and when it is badly wrong, is an essential modeling skill.
Theory. An exponential random variable with rate $\lambda$ has
$$ f(x) = \lambda e^{-\lambda x}, \quad P(X > x) = e^{-\lambda x}, \quad E[X] = \frac{1}{\lambda}, \quad \mathrm{Var}(X) = \frac{1}{\lambda^2}, \quad c^2 = 1, \qquad x \ge 0. $$
Its defining property is memorylessness:
$$ P(X > s + t \mid X > s) = \frac{e^{-\lambda (s+t)}}{e^{-\lambda s}} = e^{-\lambda t} = P(X > t). $$
An activity that has lasted $s$ minutes has the same remaining-time distribution as a fresh one. It is the only continuous distribution with this property (the geometric is its discrete counterpart). Equivalently, its hazard rate $h(x) = f(x)/P(X > x)$, the instantaneous rate of ending given survival so far, is the constant $\lambda$. A good analogy is a radioactive atom, or a light bulb that never wears out: a used one is as good as new.
The comparison with other lifetime distributions explains when the model fits:
| Hazard rate | Example distribution | Meaning | Real examples |
|---|---|---|---|
| decreasing | Weibull, shape < 1 | the longer it survived, the safer it is | infant mortality of electronics, new software |
| constant | exponential | no aging | random external shocks, time to the next call |
| increasing | Weibull shape > 1, normal, Erlang | wear-out | bearings, tires, human lifetimes |
This has direct consequences for decisions. With exponential lifetimes, preventive replacement is pointless, since a new part is no better than a working old one. Preventive maintenance pays only under wear-out (compare the maintenance MDP of Section 6.2).
Two consequences are used constantly in the rest of the chapter.
- Residual times. If service times are exponential, the remaining service of the customer currently being served is again exponential with the same rate, however long the service has lasted. So in an M/M/1 queue, an arriving customer who finds $n$ customers present waits for $n$ i.i.d. exponential services.
- Competing exponentials. If $X_1, \dots, X_k$ are independent with rates $\lambda_1, \dots, \lambda_k$, then
$$ \min_i X_i \sim \text{Exponential}\Big(\sum_i \lambda_i\Big), \qquad P(X_i \text{ is the smallest}) = \frac{\lambda_i}{\sum_j \lambda_j}, $$
and which one finishes first is independent of when it happens. Continuous-time Markov chains (Section 9.2) work exactly this way. Each possible transition has an exponential clock, and the first clock to ring decides the next state.
Finally, a sum of $k$ i.i.d. exponentials is Erlang-$k$, which is how phase-type distributions model services less variable than the exponential.
The exponential assumption is reasonable for interarrival times of unscheduled customers and for times until a random external event. It is a rough approximation for service times whose $c^2$ is near 1. It is wrong for nearly constant operations (a machine cycle, $c^2 \approx 0$) and for wear-out lifetimes.
Example: memorylessness and racing servers. Service times average 5 minutes. Part 1 compares the exponential with a uniform distribution of the same mean. Part 2 has three servers busy, with remaining times exponential at rates 0.2, 0.2 and 0.1 per minute.
import numpy as np
rng = np.random.default_rng(seed=3)
n = 1_000_000
# 1) Memorylessness: service time with mean 5 minutes
x = rng.exponential(5.0, n)
u = rng.uniform(0, 10, n) # same mean, but NOT memoryless
for name, s in [("exponential", x), ("uniform(0,10)", u)]:
fresh = np.mean(s > 3) # P(X > 3)
aged = np.mean(s[s > 5] > 8) # P(X > 8 | X > 5): 5 min already elapsed
print(f"{name:>14}: P(X>3) = {fresh:.3f} P(X>8 | X>5) = {aged:.3f}")
# 2) Competing exponentials: three busy servers, remaining times ~ Exp with these rates
rates = np.array([0.2, 0.2, 0.1]) # per minute (means 5, 5, 10)
r = rng.exponential(1 / rates, size=(n, 3)) # one column per server
first = r.min(axis=1) # time until the first service ends
who = r.argmin(axis=1)
print(f"time to first completion: mean {first.mean():.3f} (theory 1/{rates.sum():.1f} = {1 / rates.sum():.3f})")
print(f" CV {first.std() / first.mean():.3f} (exponential -> 1)")
for i in range(3):
print(f"P(server {i + 1} finishes first) = {np.mean(who == i):.3f}"
f" (theory {rates[i] / rates.sum():.3f})")
# exponential: P(X>3) = 0.550 P(X>8 | X>5) = 0.548
# uniform(0,10): P(X>3) = 0.700 P(X>8 | X>5) = 0.400
# time to first completion: mean 2.003 (theory 1/0.5 = 2.000)
# CV 0.999 (exponential -> 1)
# P(server 1 finishes first) = 0.400 (theory 0.400)
# P(server 2 finishes first) = 0.400 (theory 0.400)
# P(server 3 finishes first) = 0.200 (theory 0.200)
For the exponential, 5 minutes of elapsed service change nothing: the chance of lasting 3 more minutes is still about 0.55. For the uniform, elapsed time is informative, and the chance drops from 0.70 to 0.40 because the service is "closer to the end". With three servers, the first completion comes after an exponential time with rate $0.2 + 0.2 + 0.1 = 0.5$ (mean 2 minutes). The slow server wins the race 20% of the time, exactly its share of the total rate.
Key takeaways
- The exponential distribution has constant hazard rate $\lambda$, mean $1/\lambda$ and $c^2 = 1$. It is the only memoryless continuous distribution.
- The remaining time of an exponential activity is exponential with the same rate, whatever its age. This is why "number in system" is a sufficient state for M/M queues.
- The minimum of independent exponentials is exponential with the summed rate, and clock $i$ wins with probability $\lambda_i / \sum_j \lambda_j$.
- Use it for unscheduled arrivals and random shocks. Avoid it for near-constant operations and wear-out lifetimes, where preventive policies matter.
Practice
- (Beginner) A component's lifetime is exponential with mean 2,000 hours. What is the probability that it survives 1,000 hours? It has already worked 5,000 hours: what is the probability that it survives another 1,000? Should it be replaced preventively?
- (Intermediate) Three independent components in series have exponential lifetimes with means 1,000, 2,000 and 4,000 hours. Find the distribution and mean of the system lifetime and the probability that component 1 causes the failure. For two identical components in parallel (mean 1,000 each), use memorylessness to show that the expected system lifetime is $500 + 1000 = 1500$ hours.
- (Interview) "You arrive at a bank with two tellers. Both are busy and nobody else is waiting. Service times are exponential with the same rate. What is the probability that you are the last of the three of you to leave?" Hint: when the first teller finishes, you start service. At that moment, by memorylessness, you and the remaining customer are in a fair race.
Conditional Expectation
Why it matters. Stochastic models are analyzed by conditioning, either on what happens first or on a hidden factor: "condition on the first transition", "condition on the number of customers present", "condition on the weather". The laws of total expectation and total variance turn hard expectations into easy ones. They underlie first-step analysis of Markov chains, the Pollaczek-Khinchine formula, and safety-stock calculations with random lead times.
Theory. For a value $y$, $E[X \mid Y = y]$ is the mean of $X$ computed with the conditional distribution given $Y = y$. Viewed as a function of $Y$, the conditional expectation $E[X \mid Y]$ is itself a random variable. It is the best predictor of $X$ from $Y$ in the mean-squared-error sense.
Law of total expectation. Averaging the conditional means gives the overall mean:
$$ E[X] = E\big[E[X \mid Y]\big] = \sum_y E[X \mid Y = y], P(Y = y). $$
For example, a company's average salary is the weighted average of its department averages, with department sizes as weights.
Law of total variance. The total variability is the variability within groups plus the variability between group means:
$$ \mathrm{Var}(X) = \underbrace{E\big[\mathrm{Var}(X \mid Y)\big]}{\text{within}} + \underbrace{\mathrm{Var}\big(E[X \mid Y]\big)}{\text{between}}. $$
Random sums. Let $S = X_1 + \dots + X_N$, where the $X_i$ are i.i.d. and $N$ is a random count independent of them. Conditioning on $N$ gives $E[S \mid N] = N,E[X]$ and $\mathrm{Var}(S \mid N) = N,\mathrm{Var}(X)$, so
$$ E[S] = E[N],E[X] \quad \text{(Wald's identity)}, \qquad \mathrm{Var}(S) = E[N],\mathrm{Var}(X) + \mathrm{Var}(N),E[X]^2. $$
Random sums appear everywhere in OR. Examples are total demand during a random lead time (safety stock, Chapter 10), total claims paid by an insurer per day, and total work brought by a random number of jobs. The second variance term is the one most often forgotten.
First-step analysis. To compute an expected time or probability, condition on the first step and use the fact that the process "starts afresh". For example, let $T$ be the number of independent trials until the first success, with success probability $p$. Conditioning on the first trial gives
$$ E[T] = 1 + (1 - p),E[T] \quad \Longrightarrow \quad E[T] = \frac{1}{p}. $$
The same idea, written for all states at once, gives the absorbing-chain formulas of Section 9.2.
Example: daily insurance claims. The number of claims per day is Poisson with mean 30. Claim sizes are independent with mean 2 and standard deviation 1.5 (thousand EUR, gamma-distributed). The code checks the random-sum formulas and shows the cost of forgetting that $N$ is random.
import numpy as np
rng = np.random.default_rng(seed=11)
days = 200_000
lam = 30 # claims per day, N ~ Poisson(30)
mx, sx = 2.0, 1.5 # claim size (thousand EUR): mean 2, sd 1.5
shape, scale = (mx / sx) ** 2, sx ** 2 / mx # gamma parameters matching that mean and sd
N = rng.poisson(lam, days)
# Total paid per day: sum of N[d] independent claims. A sum of N gammas is gamma(N*shape).
S = rng.gamma(shape * np.maximum(N, 1), scale) * (N > 0)
EN, VN, EX, VX = lam, lam, mx, sx ** 2
print(f"E[S] simulated {S.mean():7.2f} formula E[N]E[X] = {EN * EX:7.2f}")
print(f"Var[S] simulated {S.var():7.2f} formula E[N]Var[X] + Var[N]E[X]^2 = {EN * VX + VN * EX ** 2:7.2f}")
print(f"Var[S] if N were fixed at 30: E[N]Var[X] = {EN * VX:7.2f}")
# Conditional expectation as a random variable: E[S | N] = 2 N
for k in (20, 30, 40):
print(f"E[S | N={k}] simulated {S[N == k].mean():6.2f} formula {k * EX:6.2f}")
# E[S] simulated 60.06 formula E[N]E[X] = 60.00
# Var[S] simulated 188.02 formula E[N]Var[X] + Var[N]E[X]^2 = 187.50
# Var[S] if N were fixed at 30: E[N]Var[X] = 67.50
# E[S | N=20] simulated 39.93 formula 40.00
# E[S | N=30] simulated 59.99 formula 60.00
# E[S | N=40] simulated 80.14 formula 80.00
The daily total averages 60 thousand EUR, as Wald's identity predicts. Its variance is 187.5, almost three times the 67.5 obtained by treating the number of claims as fixed. Reserves sized from the smaller figure would be far too thin: the standard deviation is 13.7, not 8.2. The last lines show $E[S \mid N]$ as a function of $N$: twice the number of claims.
Key takeaways
- $E[X \mid Y]$ is a random variable, the best predictor of $X$ given $Y$. Averaging it gives $E[X]$ (law of total expectation).
- Total variance = expected conditional variance (within) + variance of the conditional mean (between).
- For a random sum, $E[S] = E[N]E[X]$ and $\mathrm{Var}(S) = E[N]\mathrm{Var}(X) + \mathrm{Var}(N)E[X]^2$. Ignoring the randomness of $N$ understates risk.
- First-step analysis (condition on the first transition, then use the restart property) gives linear equations for expected times and probabilities.
Practice
- (Beginner) A store sells on average 40 units on rainy days (standard deviation 8) and 70 on sunny days (standard deviation 12). 30% of days are rainy. Compute the mean and variance of daily sales with the two laws of this topic.
- (Intermediate) The replenishment lead time is 2, 3 or 4 days with probabilities 0.2, 0.5 and 0.3. Daily demand is i.i.d. with mean 50 and standard deviation 10. Compute the mean and standard deviation of demand during the lead time, and compare with a fixed lead time of 3.1 days. Verify by simulation.
- (Interview) "How many fair coin flips do you need, on average, to see two heads in a row?" Hint: first-step analysis with two states, "last flip was not heads" and "last flip was heads". Write one linear equation per state.
<a id="92-markov-chains"></a>
9.2 Markov Chains
A Markov chain is a random process with a short memory: given the present state, the future does not depend on the past. This subchapter covers discrete-time chains (transition matrices, classification of states, steady state and absorption), then continuous-time chains and birth-death processes, which lead directly into queueing theory. The Markov decision processes of Section 6.2 are Markov chains with decisions added: under a fixed policy, an MDP is exactly a chain of this subchapter.
Discrete-Time Markov Chains
Why it matters. Many systems are described by a state that changes randomly from one period to the next. Examples are the brand a customer buys, the stock on a shelf, the condition of a machine, the credit rating of a borrower, or the page a user is viewing. A discrete-time Markov chain is the simplest model that captures the dependence between consecutive periods while remaining analytically tractable. From it one can compute long-run market shares, average costs, hitting times and the probabilities of eventual outcomes.
Theory. A discrete-time Markov chain (DTMC) is a sequence of random variables $X_0, X_1, X_2, \dots$ taking values in a finite or countable state space $\mathcal{S}$, with the Markov property
$$ P(X_{n+1} = j \mid X_n = i, X_{n-1} = i_{n-1}, \dots, X_0 = i_0) = P(X_{n+1} = j \mid X_n = i). $$
The chain is time-homogeneous if this probability does not depend on $n$. It is then written $p_{ij}$, the one-step transition probability. A homogeneous chain is completely specified by
- the initial distribution $\pi^{(0)}$, with $\pi^{(0)}_i = P(X_0 = i)$,
- the transition probabilities $p_{ij}$, collected in a matrix $P$ (next topic).
A board game like Snakes and Ladders is a Markov chain. Where you land next depends only on your current square and the die roll, not on how you reached that square.
A chain is usually drawn as a state diagram: a directed graph with one node per state and an arc $i \to j$ labeled $p_{ij}$ whenever $p_{ij} > 0$. The brand-switching example below has three brands. Each week a customer either stays with the current brand or switches:
0.80 0.75 0.80
.---. .---. .---.
| v | v | v
+-----+ 0.15 +-----+ 0.15 +-----+
| A | --------> | B | --------> | C |
| | <-------- | | <-------- | |
+-----+ 0.10 +-----+ 0.10 +-----+
^ | ^ |
| | A -> C: 0.05 | |
| +--------------------------------+ |
+------------------------------------+
C -> A: 0.10
The Markov property is a modeling decision. It holds if the state contains all the information about the past that matters for the future. If brand loyalty depends on the last two purchases, the state "brand bought this week" is not enough. Using the pair (brand last week, brand this week) as the state restores the Markov property, at the cost of a larger state space ($3^2 = 9$ states instead of 3). This state augmentation is the same trick as in dynamic programming (Section 6.1), and it carries the same risk of a combinatorial explosion.
Example: simulating brand switching and testing the Markov property.
import numpy as np
rng = np.random.default_rng(seed=5)
states = ["A", "B", "C"] # brand bought this week
P = np.array([[0.80, 0.15, 0.05], # row = this week, column = next week
[0.10, 0.75, 0.15],
[0.10, 0.10, 0.80]])
def simulate(P, x0, steps, rng):
path = [x0]
for _ in range(steps):
path.append(rng.choice(len(P), p=P[path[-1]])) # next state uses ONLY the current row
return np.array(path)
print("one customer, 15 weeks:", "".join(states[s] for s in simulate(P, 0, 15, rng)))
# A long path: estimate P from transition counts
x = simulate(P, 0, 200_000, rng)
counts = np.zeros((3, 3))
np.add.at(counts, (x[:-1], x[1:]), 1) # counts[i, j] = number of i -> j moves
print("estimated P:\n", np.round(counts / counts.sum(axis=1, keepdims=True), 3))
# Markov property check: does last week's brand matter once we know this week's?
prev, cur, nxt = x[:-2], x[1:-1], x[2:]
for p in range(3):
mask = (cur == 1) & (prev == p) # currently B, previously p
print(f"P(next=A | now=B, before={states[p]}) = {np.mean(nxt[mask] == 0):.3f}")
# one customer, 15 weeks: ABBBBAAAAACCCCCC
# estimated P:
# [[0.802 0.148 0.05 ]
# [0.099 0.75 0.15 ]
# [0.099 0.098 0.803]]
# P(next=A | now=B, before=A) = 0.101
# P(next=A | now=B, before=B) = 0.099
# P(next=A | now=B, before=C) = 0.103
The single path shows the typical behavior of a "sticky" chain: long runs on one brand, then a switch. Counting transitions in a long path recovers $P$. This is also how $P$ is estimated from real data: $\hat p_{ij}$ is the number of $i \to j$ moves divided by the number of visits to $i$. The last three lines test the Markov property. Once we know the customer bought B this week, the brand bought the week before does not change the chance of buying A next week (all about 0.10). On real panel data, a clear difference here would tell you to augment the state.
Key takeaways
- A DTMC is a sequence of states where the next state depends only on the current one: $P(X_{n+1} = j \mid \text{history}) = p_{ij}$.
- A homogeneous chain is specified by an initial distribution and the transition probabilities, and drawn as a weighted directed graph.
- The Markov property holds or fails depending on how the state is defined. Augment the state with the missing information, keeping an eye on its size.
- Transition probabilities are estimated from data by counting moves out of each state.
Practice
- (Beginner) Tomorrow is sunny with probability 0.8 if today is sunny and with probability 0.4 if today is rainy. Write $P$, draw the state diagram and simulate 30 days starting from a rainy day.
- (Intermediate) Suppose two consecutive purchases of the same brand make a customer more loyal: after AA, the probability of staying with A is 0.9. Build a chain on pairs (last brand, current brand) and fill in its rows. How many states does it have, and which transitions are impossible by construction?
- (Interview) "Can a user's navigation on a website be modeled as a Markov chain over pages? What would break the Markov property, and how would you fix it?" Hint: think of the back button, session history and time of day. Discuss state augmentation (previous page, session features) and what it does to the number of states.
Transition Matrices
Why it matters. The transition matrix turns probability questions into linear algebra. The distribution after $n$ steps is a vector-matrix product, and $n$-step transition probabilities are matrix powers. Building $P$ correctly from a verbal description of a system is also the main modeling step in any Markov analysis, and the step where most mistakes are made.
Theory. The transition matrix $P = (p_{ij})$ has one row per current state and one column per next state. Row $i$ is the conditional distribution of the next state given $X_n = i$, so $P$ is a stochastic matrix:
$$ p_{ij} \ge 0, \qquad \sum_j p_{ij} = 1 \quad \text{for every } i. $$
$n$-step transitions. The probability of going from $i$ to $j$ in exactly $n$ steps is $p^{(n)}{ij} = P(X{m+n} = j \mid X_m = i) = (P^n)_{ij}$. The reason is the Chapman-Kolmogorov equations: to go from $i$ to $j$ in $m + n$ steps, the chain must be in some state $k$ after $m$ steps, so
$$ p^{(m+n)}{ij} = \sum_k p^{(m)}{ik}, p^{(n)}_{kj}, \qquad \text{that is,} \quad P^{(m+n)} = P^{(m)} P^{(n)}. $$
Distribution after $n$ steps. With the distribution written as a row vector, $\pi^{(n)} = \pi^{(n-1)} P = \pi^{(0)} P^n$.
Building $P$ from a model takes three choices:
- When the state is observed: at the end of the week, after an order arrives, at a departure.
- What the state is: it must contain everything needed to predict the next state.
- How the state moves: a recursion $X_{n+1} = g(X_n, \text{policy}, D_{n+1})$ with an independent random input $D$. Then $p_{ij} = P\big(g(i, \cdot, D) = j\big)$.
For the inventory example below (a classic from Hillier and Lieberman), the state is the number of cameras on the shelf at the end of the week. If the shelf is empty, 3 cameras are ordered and arrive before the next week starts. Weekly demand $D$ is Poisson(1), and unmet demand is lost:
$$ X_{n+1} = \begin{cases} \max(3 - D_{n+1}, 0) & \text{if } X_n = 0, \ \max(X_n - D_{n+1}, 0) & \text{if } X_n \ge 1. \end{cases} $$
From state $i$, the stock available during the week is $s = 3$ if $i = 0$ and $s = i$ otherwise. Then $p_{ij} = P(D = s - j)$ for $1 \le j \le s$, and $p_{i0} = P(D \ge s)$. The "max with zero" collects all large demands into state 0.
Example: building and powering the inventory matrix.
import numpy as np
from math import exp, factorial
# State = cameras in stock at the end of the week (0..3).
# Policy: if stock is 0, order 3 cameras, delivered before next week. Otherwise do not order.
# Weekly demand D ~ Poisson(1); unmet demand is lost.
def pois(k, m=1.0):
return exp(-m) * m ** k / factorial(k)
S = 4
P = np.zeros((S, S))
for i in range(S):
start = 3 if i == 0 else i # stock available during the week
for j in range(1, start + 1):
P[i, j] = pois(start - j) # demand exactly start - j
P[i, 0] = 1 - sum(pois(d) for d in range(start)) # demand >= start empties the shelf
print("P =\n", np.round(P, 3))
print("row sums:", P.sum(axis=1))
# n-step transition probabilities (Chapman-Kolmogorov: P^(m+n) = P^m P^n)
for n in (2, 4, 8):
print(f"P^{n} =\n", np.round(np.linalg.matrix_power(P, n), 3))
# Evolution of the distribution, starting with 3 cameras on the shelf
pi = np.array([0, 0, 0, 1.0])
for week in range(1, 6):
pi = pi @ P # row vector times P
print(f"week {week}: P(X_n = 0..3) = {np.round(pi, 3)}")
# P =
# [[0.08 0.184 0.368 0.368]
# [0.632 0.368 0. 0. ]
# [0.264 0.368 0.368 0. ]
# [0.08 0.184 0.368 0.368]]
# row sums: [1. 1. 1. 1.]
# P^2 =
# [[0.249 0.285 0.3 0.165]
# [0.283 0.252 0.233 0.233]
# [0.351 0.319 0.233 0.097]
# [0.249 0.285 0.3 0.165]]
# P^4 =
# [[0.29 0.286 0.261 0.164]
# [0.282 0.285 0.267 0.166]
# [0.284 0.283 0.263 0.171]
# [0.29 0.286 0.261 0.164]]
# P^8 =
# [[0.286 0.285 0.263 0.166]
# [0.286 0.285 0.263 0.166]
# [0.286 0.285 0.263 0.166]
# [0.286 0.285 0.263 0.166]]
# week 1: P(X_n = 0..3) = [0.08 0.184 0.368 0.368]
# week 2: P(X_n = 0..3) = [0.249 0.285 0.3 0.165]
# week 3: P(X_n = 0..3) = [0.293 0.292 0.263 0.152]
# week 4: P(X_n = 0..3) = [0.29 0.286 0.261 0.164]
# week 5: P(X_n = 0..3) = [0.286 0.284 0.263 0.167]
Rows 0 and 3 of $P$ are identical: an empty shelf is restocked to 3, so next week looks exactly as if the week had ended with 3 cameras. Row 1 has zeros in columns 2 and 3 because stock cannot grow without an order. As $n$ grows, the rows of $P^n$ become identical. After about 8 weeks the chain has "forgotten" its starting state, and the distribution of the stock no longer depends on it. The speed of this convergence is governed by the second-largest eigenvalue of $P$ in modulus. The common limit row is the steady-state distribution, the subject of a later topic.
Key takeaways
- $P$ is a stochastic matrix: nonnegative entries, rows summing to one. Row $i$ is the distribution of the next state from $i$.
- $n$-step probabilities are entries of $P^n$ (Chapman-Kolmogorov), and distributions evolve as $\pi^{(n)} = \pi^{(0)} P^n$.
- Build $P$ by fixing when the state is observed, what it contains, and the recursion that maps (state, random input) to the next state.
- For well-behaved chains, the rows of $P^n$ converge to a common vector: the initial state is forgotten.
Practice
- (Beginner) For the brand-switching matrix of the previous topic, compute by hand (with Chapman-Kolmogorov) the probability that a customer who buys A this week buys A again two weeks from now. Check with
np.linalg.matrix_power.- (Intermediate) Change the camera policy to "order up to 4 whenever the end-of-week stock is at most 1". Build the 5-state matrix, check the row sums, and compute the distribution of the stock after 4 weeks starting from 4 cameras.
- (Interview) "How would you compute $P^n$ for $n = 10^9$, and what will it look like for a typical chain?" Hint: repeated squaring needs only $O(\log n)$ matrix products. An eigendecomposition also works. Think about the limit of the rows for an ergodic chain, and about the exceptions (next topic).
Classification of States
Why it matters. Whether a chain settles into a unique long-run behavior, keeps oscillating, or gets trapped depends on the structure of its transition graph, not on the exact probabilities. Classifying states tells you which questions make sense before you compute anything. A unique steady state is meaningful for an irreducible aperiodic chain. With several traps, the questions become which trap the chain falls into and how fast.
Theory. The definitions only look at which entries of $P$ are positive.
- State $j$ is accessible from $i$ ($i \to j$) if $p^{(n)}_{ij} > 0$ for some $n \ge 0$: there is a directed path from $i$ to $j$.
- $i$ and $j$ communicate ($i \leftrightarrow j$) if each is accessible from the other. Communication is an equivalence relation, so it partitions the states into communicating classes. In graph terms, these are the strongly connected components.
- A chain with a single class is irreducible.
- A class is closed if no arc leaves it. A state with $p_{ii} = 1$ is absorbing: it is a closed class by itself.
Recurrence. Let $f_{ii}$ be the probability that the chain, started in $i$, ever returns to $i$.
- State $i$ is recurrent if $f_{ii} = 1$: it is visited infinitely often.
- State $i$ is transient if $f_{ii} < 1$: it is visited only finitely many times, a geometric number of visits.
Recurrence is a class property. In a finite chain, a class is recurrent exactly when it is closed, and at least one closed class always exists. In infinite chains a recurrent state can be null recurrent, with an infinite mean return time. An M/M/1 queue with $\lambda = \mu$ behaves this way: it always empties eventually, but the expected time to do so is infinite.
Periodicity. The period of state $i$ is $d(i) = \gcd{n \ge 1 : p^{(n)}{ii} > 0}$. If $d(i) = 1$, the state is aperiodic. Period is also a class property. A single self-loop ($p{ii} > 0$) makes the whole class aperiodic. In a class with period $d$, the chain cycles deterministically through $d$ groups of states, so $P^n$ cannot converge.
A state that is positive recurrent and aperiodic is ergodic. A finite, irreducible, aperiodic chain is ergodic, which is the best-behaved case.
| Structure | Long-run behavior | Questions to ask |
|---|---|---|
| finite, irreducible, aperiodic | $P^n$ converges to identical rows $\pi$ | steady state, long-run averages |
| irreducible, period $d > 1$ | $P^n$ cycles; time averages still exist | long-run fractions of time |
| transient states + several closed classes | the chain eventually leaves the transient states and stays in one closed class | absorption probabilities, time to absorption |
To classify a chain on a computer, find the strongly connected components (with Tarjan's algorithm, or with the transitive closure for small chains). The closed classes are the components with no outgoing arcs. The period of a class follows from breadth-first search levels: it is the gcd of $\text{level}(i) + 1 - \text{level}(j)$ over the arcs $i \to j$ inside the class.
0.5
.---.
| v 0.5 0.4 1.0 1.0
( 0 ) ------------> ( 1 ) ----------> ( 2 ) ------> ( 3 ) ------> ( 4 )
<------------ | ^ |
0.3 | 0.3 +------------- 1.0 ----------+
v
( 5 ) --.
^ | 1.0
'-----'
classes: {0, 1} transient {2, 3, 4} recurrent, period 3 {5} absorbing
Example: classifying a six-state chain.
import numpy as np
from math import gcd
from collections import deque
P = np.array([[0.5, 0.5, 0.0, 0.0, 0.0, 0.0],
[0.3, 0.0, 0.4, 0.0, 0.0, 0.3],
[0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])
n = len(P)
# Reachability: R[i, j] = True if j can be reached from i in zero or more steps (Warshall)
R = (P > 0) | np.eye(n, dtype=bool)
for k in range(n):
R |= R[:, [k]] & R[[k], :]
def period(cls):
"""gcd of (level[i] + 1 - level[j]) over edges i -> j inside the class (BFS levels)."""
root = cls[0]
level, q = {root: 0}, deque([root])
d = 0
while q:
i = q.popleft()
for j in cls:
if P[i, j] > 0:
if j not in level:
level[j] = level[i] + 1
q.append(j)
else:
d = gcd(d, level[i] + 1 - level[j])
return d
seen = set()
for i in range(n):
if i in seen:
continue
cls = [j for j in range(n) if R[i, j] and R[j, i]] # i and j communicate
seen.update(cls)
closed = all(not (P[a, b] > 0) for a in cls for b in range(n) if b not in cls)
kind = "recurrent" if closed else "transient" # finite chain: closed <=> recurrent
extra = " (absorbing)" if len(cls) == 1 and P[cls[0], cls[0]] == 1 else ""
print(f"class {cls}: {'closed' if closed else 'open '}, {kind}{extra}, period {period(cls)}")
# Periodic class: P^n(2, 2) cycles instead of converging
print("P^n(2,2), n = 60..65:", [float(np.linalg.matrix_power(P, k)[2, 2]) for k in range(60, 66)])
# Transient states: in the long run, all mass from state 0 sits in closed classes
print("P^200 row 0:", np.round(np.linalg.matrix_power(P, 200)[0], 3))
# class [0, 1]: open , transient, period 1
# class [2, 3, 4]: closed, recurrent, period 3
# class [5]: closed, recurrent (absorbing), period 1
# P^n(2,2), n = 60..65: [1.0, 0.0, 0.0, 1.0, 0.0, 0.0]
# P^200 row 0: [0. 0. 0.286 0.124 0.162 0.429]
States 0 and 1 communicate, but the chain leaves them for good through state 2 or state 5, so they are transient. The cycle $2 \to 3 \to 4 \to 2$ is closed and has period 3: $P^n(2, 2)$ is 1 every third step and 0 otherwise, and never converges. State 5 is absorbing. Starting from 0, the chain ends up in the cycle with total probability $0.286 + 0.124 + 0.162 = 4/7$, and in state 5 with probability $3/7$. These are the ratios $0.4 : 0.3$ of the two exits from state 1. How the $4/7$ splits among states 2, 3 and 4 depends on $n \bmod 3$.
Key takeaways
- Communicating classes are the strongly connected components of the transition graph. An irreducible chain has one class.
- In a finite chain, closed classes are recurrent and all other states are transient. Transient states are eventually left forever.
- The period is the gcd of possible return times. A self-loop makes a class aperiodic. Periodic chains have no limit for $P^n$.
- Finite + irreducible + aperiodic = ergodic: a unique steady state exists and $P^n$ converges to it.
Practice
- (Beginner) Draw the graph and classify the states (classes, closed or not, recurrent or transient, period) of the chain with rows $(0.5, 0.5, 0, 0)$, $(0.5, 0.5, 0, 0)$, $(0.25, 0.25, 0.25, 0.25)$, $(0, 0, 0, 1)$. Check with the code.
- (Intermediate) The Ehrenfest model has $N$ balls in two urns. At each step a ball chosen uniformly at random moves to the other urn, and the state is the number of balls in urn 1. Show that the chain has period 2 for $N = 4$ and that $P^n$ oscillates. How does adding a "do nothing with probability 0.5" option (a lazy chain) change the classification?
- (Interview) "Why does PageRank add a small 'teleportation' probability to the random walk on the web graph?" Hint: the raw web graph has pages without outgoing links, several closed classes, and possibly periodic structure. Which property must the chain have for a unique, start-independent stationary vector?
Steady-State Probabilities
Why it matters. Many OR questions concern the long run: the average inventory cost per week, long-run market shares, the fraction of time a machine is down. For an ergodic chain, all such quantities follow from a single probability vector $\pi$, obtained by solving a linear system. This vector also underlies the average-reward criterion for MDPs (Section 6.2).
Theory. For a finite, irreducible, aperiodic chain:
$$ \lim_{n \to \infty} p^{(n)}_{ij} = \pi_j \quad \text{for every starting state } i, $$
where $\pi$ is the unique solution of the balance equations with normalization:
$$ \pi = \pi P, \qquad \sum_j \pi_j = 1, \qquad \pi_j > 0. $$
Written out, $\pi_j = \sum_i \pi_i p_{ij}$. In steady state, the probability of being in $j$ next period equals the probability flowing into $j$ from all states. Subtracting $\pi_j p_{jj}$ from both sides gives global balance: flow out of $j$ equals flow into $j$. More generally, for any set of states, the flow across the boundary is the same in both directions.
The vector $\pi$ has three interpretations:
| Interpretation | Statement | Needs aperiodicity? |
|---|---|---|
| limiting probability | $P(X_n = j) \to \pi_j$, regardless of $X_0$ | yes |
| long-run fraction of time | $\frac{1}{n}\sum_{k=1}^{n} \mathbf{1}{X_k = j} \to \pi_j$ with probability 1 | no |
| mean recurrence time | $m_{jj} = E[\text{time to return to } j \mid X_0 = j] = 1/\pi_j$ | no |
The second interpretation extends to rewards (the ergodic theorem). If being in state $j$ costs $c(j)$ per period, the long-run average cost per period is
$$ \lim_{n \to \infty} \frac{1}{n} \sum_{k=1}^{n} c(X_k) = \sum_j \pi_j, c(j) = \pi c \quad \text{with probability 1}. $$
Computation.
- Direct solve. $\pi(P - I) = 0$ has one redundant equation. Replace any one of them with $\sum_j \pi_j = 1$ and solve the resulting nonsingular system. This is the method of choice for up to a few thousand states.
- Power method. Iterate $\pi \leftarrow \pi P$ from any start. It converges geometrically for ergodic chains, and it is how PageRank is computed on billions of pages.
- Eigenvector. $\pi^\top$ is the eigenvector of $P^\top$ for eigenvalue 1, normalized to sum to 1.
- Sparse iterative methods (Gauss-Seidel, GMRES) for large structured chains.
If the chain is periodic, $\pi$ still exists and gives time fractions, but $P^n$ does not converge. If it is reducible with several closed classes, each class has its own stationary vector, and the long-run behavior depends on where the chain starts (combine with the absorption probabilities of the next topic).
Example: long-run behavior and cost of the camera inventory. The code computes $\pi$ for the inventory chain of the Transition Matrices topic and checks it against $P^{50}$ and a simulation. It then prices the policy: ordering costs 10 plus 25 per camera, each camera left at the end of the week costs 2, and each lost sale costs 50.
import numpy as np
from math import exp, factorial
P = np.array([[0.080, 0.184, 0.368, 0.368], # camera inventory chain (previous topic)
[0.632, 0.368, 0.000, 0.000],
[0.264, 0.368, 0.368, 0.000],
[0.080, 0.184, 0.368, 0.368]])
n = len(P)
# Solve pi (P - I) = 0 with sum(pi) = 1: replace one (redundant) balance equation
A = (P - np.eye(n)).T
A[-1, :] = 1.0 # last equation becomes the normalization
b = np.zeros(n); b[-1] = 1.0
pi = np.linalg.solve(A, b)
print("steady state pi :", np.round(pi, 3))
print("row of P^50 :", np.round(np.linalg.matrix_power(P, 50)[0], 3))
print("mean recurrence 1/pi :", np.round(1 / pi, 2), "weeks")
# Long-run fraction of time from one simulated path (ergodic theorem)
rng = np.random.default_rng(seed=8)
x, visits, T = 3, np.zeros(n), 100_000
for _ in range(T):
x = rng.choice(n, p=P[x])
visits[x] += 1
print("simulated fractions :", np.round(visits / T, 3))
# Long-run average weekly cost: order (10 + 25/camera) when empty, 2 per unsold camera,
# 50 per lost sale. Cost depends on the stock available at the start of the week.
pois = lambda k: exp(-1) * 1 / factorial(k)
def weekly_cost(i):
start = 3 if i == 0 else i
order = 10 + 25 * 3 if i == 0 else 0
exp_left = sum(pois(d) * max(start - d, 0) for d in range(30))
exp_lost = sum(pois(d) * max(d - start, 0) for d in range(30))
return order + 2 * exp_left + 50 * exp_lost
c = np.array([weekly_cost(i) for i in range(n)])
print("expected cost by state:", np.round(c, 2))
print(f"long-run average cost : {pi @ c:.2f} per week")
# steady state pi : [0.286 0.285 0.263 0.166]
# row of P^50 : [0.286 0.285 0.263 0.166]
# mean recurrence 1/pi : [3.5 3.51 3.8 6.01] weeks
# simulated fractions : [0.286 0.284 0.263 0.167]
# expected cost by state: [90.21 19.13 7.39 5.21]
# long-run average cost : 34.03 per week
The three routes to $\pi$ agree: the linear solve, the rows of $P^{50}$, and the time fractions along one simulated path. The shelf is empty at the end of 28.6% of weeks, so the store orders on average every $1/\pi_0 = 3.5$ weeks. Each cost is attached to an end-of-week state $i$. It covers the order placed at that moment (when $i = 0$) and the expected holding and lost-sales costs of the following week. The long-run average cost is about 34 per week. Comparing that number across alternative policies is how a Markov chain becomes a decision tool. An MDP (Section 6.2) automates the comparison.
Key takeaways
- For a finite, irreducible, aperiodic chain, $P^n$ converges to identical rows $\pi$, the unique solution of $\pi = \pi P$, $\sum_j \pi_j = 1$.
- $\pi_j$ is the limiting probability, the long-run fraction of time in $j$, and the reciprocal of the mean return time to $j$.
- Long-run average cost or reward is $\sum_j \pi_j c(j)$.
- Solve by replacing one balance equation with the normalization. Use the power method or iterative solvers for very large chains.
Practice
- (Beginner) Compute the long-run market shares of the three brands in the brand-switching chain, by hand (two balance equations and the normalization) and with code. Which brand gains share relative to an equal split?
- (Intermediate) For the camera store, evaluate the policy "order up to 3 whenever the end-of-week stock is at most 1" with the same costs. Build its chain, compute $\pi$ and the long-run average cost, and compare with 34.03.
- (Interview) "You have logs of 10 million user sessions as page sequences. How would you estimate the long-run fraction of time users spend on each page, and how would you check whether a stationary Markov model is appropriate?" Hint: estimate $P$ by counting transitions, solve for $\pi$, and compare with the empirical frequencies. Check stability of $\hat P$ across time windows and test the Markov property (does the previous page matter?).
Absorbing Chains
Why it matters. Some processes end. A loan is repaid or defaults, a patient is discharged or dies, a student graduates or drops out, a project task completes. The questions are then "how long until it ends?" and "how does it end?". Absorbing-chain theory answers both with one matrix inverse.
Theory. A chain is absorbing if it has at least one absorbing state and an absorbing state can be reached from every state. Order the states with the $t$ transient states first and the $r$ absorbing states last. The transition matrix then has the canonical form
$$ P = \begin{pmatrix} Q & R \ 0 & I \end{pmatrix}, $$
where $Q$ ($t \times t$) holds the transient-to-transient probabilities and $R$ ($t \times r$) the transient-to-absorbing ones. Because absorption is eventually certain, $Q^n \to 0$, and the fundamental matrix exists:
$$ N = (I - Q)^{-1} = I + Q + Q^2 + \cdots . $$
Everything follows from $N$:
| Quantity | Formula | Meaning |
|---|---|---|
| $n_{ij}$ | entry of $N$ | expected number of visits to transient $j$ starting from $i$, counting the start if $i = j$ |
| $\mathbf{t} = N \mathbf{1}$ | row sums of $N$ | expected number of steps before absorption |
| $B = N R$ | $t \times r$ matrix | $b_{ik}$ = probability of ending in absorbing state $k$ from $i$ |
| $(2N - I)\mathbf{t} - \mathbf{t} \circ \mathbf{t}$ | vector | variance of the number of steps ($\circ$ = elementwise product) |
The same results follow from first-step analysis (Conditional Expectation topic). From transient $i$, one step is taken, then the chain either is absorbed or continues from a transient state $k$:
$$ t_i = 1 + \sum_k q_{ik}, t_k ;\Rightarrow; (I - Q)\mathbf{t} = \mathbf{1}, \qquad b_{ik} = r_{ik} + \sum_{l} q_{il}, b_{lk} ;\Rightarrow; B = R + QB. $$
Think of a maze with several exits. $N$ counts how often you pass through each room, $\mathbf{t}$ how long you wander, and $B$ which exit you take.
Applications include accounts receivable (aging categories, with "paid" and "bad debt" absorbing), credit rating migration (with "default" absorbing), patient flow, manpower planning, gambler's ruin, and the expected duration of random projects. The technique also computes first passage times in ergodic chains: to find the expected time to reach state $j$, make $j$ absorbing and apply $\mathbf{t} = N\mathbf{1}$.
Example: progression through a three-year degree. Each year a student either repeats the year, moves on, or drops out. Graduation and dropping out are absorbing.
import numpy as np
names = ["Year 1", "Year 2", "Year 3", "Graduated", "Dropped out"]
# Canonical form: transient states first (years), then absorbing states
Q = np.array([[0.10, 0.75, 0.00], # transient -> transient
[0.00, 0.10, 0.80],
[0.00, 0.00, 0.10]])
R = np.array([[0.00, 0.15], # transient -> absorbing (graduate, drop out)
[0.00, 0.10],
[0.85, 0.05]])
N = np.linalg.inv(np.eye(3) - Q) # fundamental matrix: expected visits
t = N @ np.ones(3) # expected years until absorption
B = N @ R # absorption probabilities
print("N (expected years spent in each year, by starting year):\n", np.round(N, 3))
for i in range(3):
print(f"start in {names[i]}: {t[i]:.2f} years enrolled,"
f" P(graduate) = {B[i, 0]:.3f}, P(drop out) = {B[i, 1]:.3f}")
# Sanity check by simulation for a new student
rng = np.random.default_rng(seed=9)
P = np.block([[Q, R], [np.zeros((2, 3)), np.eye(2)]])
years, grads = [], 0
for _ in range(100_000):
s, k = 0, 0
while s < 3: # still enrolled
s = rng.choice(5, p=P[s]); k += 1
years.append(k); grads += (s == 3)
print(f"simulated: {np.mean(years):.2f} years, P(graduate) = {grads / 100_000:.3f}")
# N (expected years spent in each year, by starting year):
# [[1.111 0.926 0.823]
# [0. 1.111 0.988]
# [0. 0. 1.111]]
# start in Year 1: 2.86 years enrolled, P(graduate) = 0.700, P(drop out) = 0.300
# start in Year 2: 2.10 years enrolled, P(graduate) = 0.840, P(drop out) = 0.160
# start in Year 3: 1.11 years enrolled, P(graduate) = 0.944, P(drop out) = 0.056
# simulated: 2.86 years, P(graduate) = 0.699
A new student spends on average 1.11 years in Year 1 ($1/(1 - 0.1)$, because of repeats) and is enrolled for 2.86 years in total. That is less than 3 years because dropouts leave early. Only 70% of new students graduate, while 94.4% of those who reach Year 3 do. So retention work in the first year matters most. Raising Year 1's continuation probability from 0.75 to 0.80 (with dropouts falling to 0.10) lifts the graduation rate from 70.0% to 74.6%. Eliminating Year 3 dropouts entirely would only reach 74.1%. The simulation confirms both numbers.
Key takeaways
- Write $P$ in canonical form with transient states first. $Q$ is the transient block and $R$ leads into the absorbing states.
- The fundamental matrix $N = (I - Q)^{-1}$ gives expected visits. $N\mathbf{1}$ gives expected time to absorption, and $NR$ gives absorption probabilities.
- The same equations follow from first-step analysis: $(I - Q)\mathbf{t} = \mathbf{1}$ and $B = R + QB$.
- Making a target state absorbing turns first-passage questions in any chain into absorbing-chain questions.
Practice
- (Beginner) Gambler's ruin: you have 2 EUR and bet 1 EUR per round, winning with probability 0.4, until you have 0 or 4 EUR. Build $Q$ and $R$, then compute the probability of reaching 4 EUR and the expected number of rounds.
- (Intermediate) In the camera inventory chain, make state 0 absorbing and compute the expected number of weeks until the shelf is first empty, starting from 1, 2 and 3 cameras. Relate the answer from state 3 to the mean recurrence time $1/\pi_0 = 3.5$ weeks of the previous topic (remember that an empty shelf is restocked to 3).
- (Interview) "A bank has an annual credit-rating migration matrix with ratings from AAA to CCC plus Default. How do you compute the 5-year default probability of a BBB loan, and what assumptions does this make?" Hint: Default is absorbing, so look at the right entry of $P^5$. Then question time-homogeneity (economic cycles) and the Markov property (rating momentum: recently downgraded firms tend to be downgraded again).
Continuous-Time Markov Chains
Why it matters. In many systems, events happen at any moment rather than at regular ticks: failures, repairs, arrivals, departures. A continuous-time Markov chain (CTMC) describes such systems with rates instead of per-step probabilities. Its long-run behavior follows from "rate in = rate out" balance equations. CTMCs are the mathematical engine behind Markovian queueing models, reliability models and many performance models of computer systems.
Theory. A process ${X(t), t \ge 0}$ on a discrete state space is a CTMC if, given the present state, the future is independent of the past. Combined with time-homogeneity, this forces the following structure:
- In state $i$, the process stays for an exponential holding time with rate $\nu_i$. Exponential is the only distribution consistent with "the past does not matter" (memorylessness).
- It then jumps to $j \neq i$ with probability $p_{ij}$. These probabilities form the embedded jump chain, a DTMC.
The transition rates are $q_{ij} = \nu_i p_{ij}$. An equivalent picture is that each possible transition $i \to j$ has its own exponential clock with rate $q_{ij}$, and the first clock to ring decides the move. By the competing-exponentials property, the two pictures agree. Collected in a matrix, the rates form the generator $Q$:
$$ q_{ij} \ge 0 ;(i \ne j), \qquad q_{ii} = -\nu_i = -\sum_{j \ne i} q_{ij}, \qquad \text{so every row of } Q \text{ sums to } 0. $$
For a small interval $h$, $P(X(t+h) = j \mid X(t) = i) = q_{ij} h + o(h)$.
Transient behavior. The matrix $P(t)$ with entries $P(X(t) = j \mid X(0) = i)$ satisfies the Kolmogorov equations $P'(t) = P(t)Q = QP(t)$ with $P(0) = I$, whose solution is the matrix exponential
$$ P(t) = e^{Qt} = \sum_{n \ge 0} \frac{(Qt)^n}{n!}. $$
In practice $e^{Qt}$ is computed with scipy.linalg.expm. For large chains, uniformization is used instead: pick $\Lambda \ge \max_i \nu_i$ and write $P(t) = \sum_n e^{-\Lambda t}\frac{(\Lambda t)^n}{n!} \tilde P^n$ with the stochastic matrix $\tilde P = I + Q/\Lambda$.
Steady state. An irreducible finite CTMC always has a limit, since periodicity cannot occur in continuous time:
$$ \lim_{t \to \infty} P(X(t) = j) = \pi_j, \qquad \pi Q = 0, \qquad \sum_j \pi_j = 1. $$
Row $j$ of $\pi Q = 0$ reads
$$ \underbrace{\pi_j, \nu_j}{\text{rate out of } j} = \underbrace{\sum{i \ne j} \pi_i, q_{ij}}_{\text{rate into } j}. $$
These balance equations are usually written directly from the state diagram. The stationary distribution of the CTMC differs from that of the jump chain $\tilde\pi$: time weighs each visit by its mean length, so $\pi_j \propto \tilde\pi_j / \nu_j$.
0.20
+----+ -------> +----------+
| Up | <------- | Degraded |
+----+ 1.00 +----------+
^ | |
0.50 | | 0.05 | 0.50
| v |
+------+ <-----------+
| Down |
+------+
rates per day; each arrow is an exponential clock
Example: availability of a machine with a degraded mode. A machine is Up, Degraded (still producing) or Down. From Up it degrades at rate 0.2 per day or fails suddenly at rate 0.05. A degraded machine gets a quick fix at rate 1.0 or fails at rate 0.5. A repair takes 2 days on average.
import numpy as np
from scipy.linalg import expm
names = ["Up", "Degraded", "Down"]
# Generator Q (rates per day): off-diagonal q_ij >= 0, rows sum to 0
Q = np.array([[-0.25, 0.20, 0.05], # Up -> Degraded 0.2, Up -> Down 0.05
[ 1.00, -1.50, 0.50], # Degraded -> Up 1.0 (quick fix), -> Down 0.5
[ 0.50, 0.00, -0.50]]) # Down -> Up 0.5 (repair takes 2 days on average)
# Transient behaviour: P(t) = exp(Qt), starting from Up
for t in (1, 7, 30, 365):
print(f"t = {t:>3} days: P(state | start Up) = {np.round(expm(Q * t)[0], 4)}")
# Stationary distribution: pi Q = 0, sum(pi) = 1
A = Q.T.copy()
A[-1, :] = 1.0
pi = np.linalg.solve(A, np.array([0, 0, 1.0]))
print("stationary pi:", np.round(pi, 4), f"-> availability {pi[0] + pi[1]:.4f}")
# Simulation: exponential holding time with rate -q_ii, then jump with prob q_ij / -q_ii
rng = np.random.default_rng(seed=10)
time_in, s, clock = np.zeros(3), 0, 0.0
while clock < 200_000:
rate = -Q[s, s]
hold = rng.exponential(1 / rate)
time_in[s] += hold; clock += hold
jump = np.where(np.arange(3) == s, 0.0, Q[s]) / rate # embedded jump chain row
s = rng.choice(3, p=jump)
print("simulated time fractions:", np.round(time_in / clock, 4))
# t = 1 days: P(state | start Up) = [0.847 0.0925 0.0605]
# t = 7 days: P(state | start Up) = [0.7339 0.0981 0.168 ]
# t = 30 days: P(state | start Up) = [0.7317 0.0976 0.1707]
# t = 365 days: P(state | start Up) = [0.7317 0.0976 0.1707]
# stationary pi: [0.7317 0.0976 0.1707] -> availability 0.8293
# simulated time fractions: [0.732 0.0976 0.1704]
Within a week, the effect of starting in Up has almost vanished, and $P(t)$ approaches the stationary vector. The machine produces 82.9% of the time. Of the 17.1% downtime, the balance equations tell us where it comes from. The rate of entering Down is $0.05,\pi_{\text{Up}} + 0.5,\pi_{\text{Deg}} = 0.0366 + 0.0488$ per day, so 57% of failures come through the degraded mode. Reducing that path, for example with condition monitoring that triggers the quick fix sooner, is the most effective lever. The simulated time fractions confirm $\pi$.
Key takeaways
- A CTMC holds in each state for an exponential time with rate $\nu_i$, then jumps according to the embedded chain. Equivalently, each transition $i \to j$ has an exponential clock with rate $q_{ij}$.
- The generator $Q$ has nonnegative off-diagonal rates and rows summing to zero. $P(t) = e^{Qt}$.
- The stationary distribution solves $\pi Q = 0$, $\sum_j \pi_j = 1$: rate out of each state equals rate in.
- CTMC and jump-chain stationary vectors differ, because time weighs states by their mean holding times.
Practice
- (Beginner) A single machine fails at rate 0.1 per hour and is repaired at rate 0.5 per hour. Write $Q$, solve $\pi Q = 0$ by hand to get the availability, and compute $P(\text{up at } t = 2 \mid \text{up at } 0)$ with
expm.- (Intermediate) Two identical machines (each failing at rate 0.1 per hour) share one repairer (rate 0.5 per hour). With the number of machines down as the state, write the $3 \times 3$ generator, compute $\pi$, the expected number of machines up, and the repairer's utilization.
- (Interview) "What is the difference between the stationary distribution of a CTMC and that of its embedded jump chain? When are they equal?" Hint: the jump chain counts visits and the CTMC measures time. Relate them through the holding rates $\nu_j$, and consider what uniformization does to them.
Birth-Death Processes
Why it matters. When the state is a count (customers in a system, machines down, busy servers, a population) and changes by one unit at a time, the CTMC is a birth-death process. Its steady state has a closed-form product formula. That formula is the common source of the M/M/1, M/M/c, finite-capacity and machine-repair results in the next subchapter.
Theory. A birth-death process has states $0, 1, 2, \dots$ and only two kinds of transitions:
- births $n \to n+1$ at rate $\lambda_n$,
- deaths $n \to n-1$ at rate $\mu_n$ (with $\mu_0 = 0$).
lam_0 lam_1 lam_2 lam_{n-1} lam_n
(0) -------> (1) -------> (2) -------> ... (n-1) ----------> (n) ----------> ...
<------- <------- <------- <---------- <----------
mu_1 mu_2 mu_3 mu_n mu_{n+1}
|cut| flow up across a cut = flow down across it
Balance the flows across the cut between $n$ and $n + 1$. The only way up is a birth from $n$, and the only way down is a death from $n+1$, so
$$ \lambda_n, p_n = \mu_{n+1}, p_{n+1} \qquad \Longrightarrow \qquad p_n = p_0 \prod_{k=0}^{n-1} \frac{\lambda_k}{\mu_{k+1}}, \qquad p_0 = \Big(1 + \sum_{n \ge 1} \prod_{k=0}^{n-1} \frac{\lambda_k}{\mu_{k+1}}\Big)^{-1}. $$
These cut equations (a form of detailed balance) are much easier than the full system $\pi Q = 0$, which they imply. A steady state exists exactly when the series in $p_0$ converges. On a finite state space it always does. On an infinite one, convergence requires births to be eventually slower than deaths.
One formula covers many models, depending on the rates:
| Model | $\lambda_n$ | $\mu_n$ | Steady state |
|---|---|---|---|
| M/M/1 queue | $\lambda$ | $\mu$ | geometric, needs $\lambda < \mu$ |
| M/M/c queue | $\lambda$ | $\min(n, c),\mu$ | Erlang C, needs $\lambda < c\mu$ |
| M/M/c/K (finite room) | $\lambda$ for $n < K$, 0 at $K$ | $\min(n, c),\mu$ | truncated, always exists |
| M/M/$\infty$ (self-service) | $\lambda$ | $n\mu$ | Poisson($\lambda/\mu$) |
| Machine repair, $M$ machines, $r$ repairers | $(M - n)\lambda$ | $\min(n, r),\mu$ | finite, always exists |
| Balking customers | $\lambda, b_n$ ($b_n$ = probability of joining) | $\mu$ | depends on $b_n$ |
In steady state, the rate of births equals the rate of deaths: $\sum_n \lambda_n p_n = \sum_n \mu_n p_n$. This throughput identity is a useful check on any computation.
Example: machine repair with one or two technicians. Five machines each fail at rate 0.1 per hour while working, and each technician repairs at rate 0.5 per hour. The state is the number of machines down. The function birth_death_stationary is reused in Section 9.3.
import numpy as np
def birth_death_stationary(lam, mu):
"""lam[n] = birth rate in state n (n = 0..K-1), mu[n] = death rate in state n+1.
Returns p_0..p_K using p_{n+1} = p_n * lam[n] / mu[n]."""
w = np.concatenate(([1.0], np.cumprod(np.asarray(lam) / np.asarray(mu))))
return w / w.sum()
# Machine repair: M machines, each fails at rate f per hour while working,
# r technicians, each repairs at rate g per hour. State n = machines down.
M, f, g = 5, 0.1, 0.5
for r in (1, 2):
lam = [(M - n) * f for n in range(M)] # fewer machines up -> fewer failures
mu = [min(n, r) * g for n in range(1, M + 1)] # at most r repairs in parallel
p = birth_death_stationary(lam, mu)
down = p @ np.arange(M + 1)
busy = p @ np.minimum(np.arange(M + 1), r)
print(f"r = {r}: p = {np.round(p, 4)}")
print(f" E[down] = {down:.3f}, machines up = {M - down:.3f},"
f" technician utilization = {busy / r:.3f}")
# Rate balance check: failure rate = repair rate in the long run
print(f" failures/h = {f * (M - down):.4f}, repairs/h = {g * busy:.4f}")
# r = 1: p = [0.2849 0.2849 0.2279 0.1367 0.0547 0.0109]
# E[down] = 1.424, machines up = 3.576, technician utilization = 0.715
# failures/h = 0.3576, repairs/h = 0.3576
# r = 2: p = [0.3927 0.3927 0.1571 0.0471 0.0094 0.0009]
# E[down] = 0.891, machines up = 4.109, technician utilization = 0.411
# failures/h = 0.4109, repairs/h = 0.4109
With one technician, 1.42 machines are down on average. A second technician cuts this to 0.89 and adds 0.53 machines of average capacity, but each technician is then busy only 41% of the time. Whether the second technician pays off depends on the value of a machine-hour against the technician's wage. Note that $\lambda_n$ decreases as machines go down. In a finite-source model, a busy repair shop slows down its own arrivals, unlike a queue fed by an infinite population. In both cases, failures per hour equal repairs per hour, as the throughput identity requires.
Key takeaways
- A birth-death process moves up at rate $\lambda_n$ and down at rate $\mu_n$, one step at a time.
- Cut equations $\lambda_n p_n = \mu_{n+1} p_{n+1}$ give the product formula $p_n = p_0 \prod_{k<n} \lambda_k/\mu_{k+1}$. A steady state exists when the normalizing series converges.
- M/M/1, M/M/c, M/M/c/K, M/M/$\infty$ and machine-repair models are all special cases, differing only in the rates.
- Check results with the throughput identity: total birth rate = total death rate.
Practice
- (Beginner) Show that the M/M/$\infty$ queue ($\lambda_n = \lambda$, $\mu_n = n\mu$) has a Poisson($\lambda/\mu$) stationary distribution. Apply it to a parking lot with 30 arrivals per hour and a mean stay of 2 hours: what is the probability that more than 70 cars are parked?
- (Intermediate) At a barbershop, an arriving customer who finds $n$ customers present joins with probability $1/(n+1)$. Arrivals are Poisson at 6 per hour, there is one barber, and service is exponential at 4 per hour. Compute $p_n$ (truncate at a large $n$) and the fraction of arriving customers who leave without joining (use PASTA).
- (Interview) "How would you decide the number of technicians for a fleet of 20 machines?" Hint: compute the expected number of machines down for each $r$ with the birth-death formula. Price downtime per machine-hour and technicians per hour, and add technicians while the marginal saving exceeds the marginal cost.
<a id="93-queueing-theory"></a>
9.3 Queueing Theory
Queues form whenever demand for a service arrives faster, at least temporarily, than it can be served. Queueing theory predicts waiting times, queue lengths and utilization from the arrival and service processes, so capacity can be sized before customers start waiting. This subchapter introduces the standard notation and Little's law, and solves the classic Markovian queues as birth-death processes. It then treats general service times with the Pollaczek-Khinchine formula and ends with networks of queues. Unless stated otherwise, single-station examples use one coffee bar: 20 customers per hour (a mean gap of 3 minutes) and a mean service time of 2.5 minutes (24 customers per hour).
Kendall Notation
Why it matters. Queueing results are catalogued by model, and Kendall notation is the shared name of a model. A short code such as M/M/c states the arrival process, the service distribution and the number of servers, and points to the right formula. Before applying any formula, an analyst must also know which performance measures to compute and when a steady state exists at all.
Theory. Every queueing system has the same anatomy:
calling arrival waiting room service facility departures
population ---> process A ----> capacity K, order D ----> c parallel servers ---->
(size N) (rate lambda) [ o o o o o ] [S1] [S2] ... [Sc]
service times S (rate mu)
Kendall's notation $A/S/c/K/N/D$ lists these components. The last three are omitted when they take their defaults ($K = \infty$, $N = \infty$, FCFS):
| Position | Meaning | Common values |
|---|---|---|
| $A$ | interarrival-time distribution | M (exponential: Poisson arrivals), D (deterministic), $E_k$ (Erlang-$k$), G or GI (general, independent) |
| $S$ | service-time distribution | M, D, $E_k$, $H_k$ (hyperexponential), PH (phase-type), G |
| $c$ | number of parallel servers | 1, 2, ..., $\infty$ |
| $K$ | system capacity (in service + waiting) | $\infty$ by default |
| $N$ | size of the calling population | $\infty$ by default; finite in machine-repair models |
| $D$ | queue discipline | FCFS (default), LCFS, SIRO (random order), priority, PS (processor sharing) |
"M" stands for Markovian (memoryless). In the arrival position it means Poisson arrivals, and in the service position exponential service times. Real systems add behaviors that the basic code does not show: balking (refusing to join a long queue), reneging or abandonment (leaving after waiting too long, as in the Erlang-A model M/M/c+M), jockeying between lines, batch arrivals ($M^X$), and priorities.
Performance measures. All formulas in this subchapter refer to long-run (steady-state) averages:
| Symbol | Meaning |
|---|---|
| $\lambda$, $\mu$ | arrival rate, service rate per server ($1/\mu = E[S]$) |
| $\rho = \lambda/(c\mu)$ | utilization: fraction of time a server is busy |
| $L$, $L_q$ | mean number in the system, mean number waiting in the queue |
| $W$, $W_q$ | mean time in the system, mean waiting time before service |
| $p_n$ | probability of $n$ customers in the system |
| $P(W_q > t)$, $P(\text{wait})$ | service-level measures |
| $p_K$ (or $B$) | blocking probability in finite systems |
Stability. With unlimited waiting room, a steady state exists only if $\rho < 1$: the servers must be able to handle the arrivals on average. With $\rho \ge 1$ the queue grows without bound. Even with $\rho < 1$, the steady-state formulas describe a system after a long run under constant conditions. They describe a lunch rush that lasts 40 minutes poorly, and such transient periods call for time-dependent analysis or simulation (Chapter 12).
| Real system | Kendall code |
|---|---|
| single ATM, random customers, variable transaction times | M/G/1 |
| call center with 25 agents and unlimited hold queue | M/M/25 (or M/M/25+M with abandonment) |
| toll plaza with 3 booths, one queue, fixed 10-second processing | M/D/3 |
| hospital ward with 12 beds and no waiting room | M/G/12/12 |
| one technician serving 10 machines | M/M/1//10 |
Example: one simulator for every code. Each position of the code maps to an input of a simulator. The function below handles any FCFS G/G/c queue, and later topics reuse it to check formulas. It runs the coffee bar under different Kendall codes with the same mean arrival and service rates.
import heapq
import numpy as np
def simulate_queue(interarrival, service, c=1, n=200_000, warmup=20_000, seed=0):
"""FCFS G/G/c queue with unlimited waiting room.
interarrival(rng, n) and service(rng, n) return n samples (minutes).
Returns (mean wait in queue Wq, mean time in system W)."""
rng = np.random.default_rng(seed)
arrive = np.cumsum(interarrival(rng, n)) # arrival epochs
s = service(rng, n)
free = [0.0] * c # min-heap: time each server becomes free
wq = np.empty(n)
for k in range(n):
start = max(arrive[k], heapq.heappop(free)) # earliest free server, FCFS order
wq[k] = start - arrive[k]
heapq.heappush(free, start + s[k])
wq, s = wq[warmup:], s[warmup:] # drop the start-up (empty system) period
return wq.mean(), (wq + s).mean()
# Coffee bar: 20 customers/hour (mean gap 3 min), mean service 2.5 min
M = lambda mean: (lambda rng, n: rng.exponential(mean, n)) # Markovian (exponential)
D = lambda mean: (lambda rng, n: np.full(n, mean)) # deterministic
E = lambda k, mean: (lambda rng, n: rng.gamma(k, mean / k, n)) # Erlang-k
systems = { # Kendall code: (A, S, c)
"D/D/1": (D(3), D(2.5), 1),
"M/D/1": (M(3), D(2.5), 1),
"M/E4/1": (M(3), E(4, 2.5), 1),
"M/M/1": (M(3), M(2.5), 1),
"M/M/2": (M(3), M(5.0), 2), # two slower baristas, same capacity
}
for code, (A, S, c) in systems.items():
wq, w = simulate_queue(A, S, c)
print(f"{code:>7}: Wq = {wq:6.2f} min W = {w:6.2f} min")
# D/D/1: Wq = 0.00 min W = 2.50 min
# M/D/1: Wq = 6.38 min W = 8.88 min
# M/E4/1: Wq = 7.99 min W = 10.49 min
# M/M/1: Wq = 12.56 min W = 15.07 min
# M/M/2: Wq = 11.41 min W = 16.41 min
All five systems have utilization $\rho = 2.5/3 \approx 0.83$, yet the waits range from zero to over 12 minutes. In D/D/1, a customer arrives every 3 minutes and leaves after 2.5, so nobody ever waits: queues are caused by variability, not by a lack of capacity. Adding randomness to arrivals (M/D/1), then to service (M/E4/1, M/M/1), steadily increases the wait. The last line holds a surprise: two baristas at half speed give a shorter queue but a longer total time than one fast barista. Later topics explain both effects.
Key takeaways
- Kendall's $A/S/c/K/N/D$ names the interarrival distribution, service distribution, servers, capacity, population and discipline. M means memoryless (Poisson arrivals or exponential service).
- The standard measures are $\rho$, $L$, $L_q$, $W$, $W_q$, $p_n$, waiting-time tails and blocking probabilities, all long-run averages.
- With unlimited room, stability requires $\rho = \lambda/(c\mu) < 1$.
- Variability, not average load alone, creates waiting: a D/D/1 queue with $\rho < 1$ never waits.
Practice
- (Beginner) Write the Kendall code for: (a) a car wash with one bay, Poisson arrivals and a fixed 6-minute program; (b) a 4-bed intensive care unit with no waiting area and general lengths of stay; (c) a web server with 8 worker threads, a buffer for 100 requests in total, Poisson arrivals and exponential service; (d) one technician maintaining 10 machines with exponential failure and repair times.
- (Intermediate) Use
simulate_queueto estimate $W_q$ for D/M/1 (an arrival every 3 minutes, exponential service with mean 2.5). Compare with M/D/1 and M/M/1. Is arrival variability or service variability more harmful here?- (Interview) "Our servers run at 95% average utilization, so we have 5% spare capacity and everything is fine." How do you respond? Hint: in most queueing models delay grows like $1/(1-\rho)$. Consider also that 95% on average hides peaks above 100% during the day.
Little's Law
Why it matters. Little's law is the most general and most used result in queueing theory. It links three quantities that managers measure (how many are inside, how fast they arrive, how long they stay) and holds for almost any system in steady state. That includes any distribution, any number of servers, any discipline, and any part of a larger system. It lets you compute the hardest of the three from the two easiest.
Theory. For a system in steady state,
$$ L = \lambda, W, $$
where $L$ is the long-run average number in the system, $\lambda$ the long-run rate at which customers enter, and $W$ the average time a customer spends inside.
Why it is true. Plot the number in the system $N(t)$ over $[0, T]$, and suppose the system is empty at both ends:
3 | ####
2 | ##############
1 | ######################## ######
0 +----+-----+-----+---+---+---+------+-----+------> t
a1 a2 a3 d1 d2 d3 a4 d4
each customer contributes a horizontal strip from its arrival to its departure
The area under $N(t)$ can be counted in two ways: by time, it is $\int_0^T N(t),dt$; by customer, it is $\sum_i (d_i - a_i)$, the sum of the times in system. Dividing by $T$ and multiplying and dividing by the number of customers $n$:
$$ \underbrace{\frac{1}{T}\int_0^T N(t),dt}{\to L} = \underbrace{\frac{n}{T}}{\to \lambda} \cdot \underbrace{\frac{1}{n}\sum_i (d_i - a_i)}_{\to W}. $$
Letting $T \to \infty$ gives the law. The argument uses no distribution, no number of servers and no queue discipline.
Apply it to any "box".
| Box | Law | Consequence |
|---|---|---|
| whole system | $L = \lambda W$ | |
| waiting line only | $L_q = \lambda W_q$ | |
| servers only | busy servers $= \lambda E[S] = \lambda/\mu$ | $\rho = \lambda/(c\mu)$ is the fraction of busy servers |
| combining | $W = W_q + 1/\mu$ | $L = L_q + \lambda/\mu$ |
| one station of a network | $L_j = \lambda_j W_j$ | |
| whole network | $L = \gamma W$, with $\gamma$ the external arrival rate |
In blocked or balking systems, $\lambda$ must be the rate of customers who actually enter ($\lambda_{\text{eff}}$), not the rate of those who try.
Little's law also applies outside classic queues. In a factory, work-in-process = throughput × cycle time. In a warehouse, average inventory = demand rate × average time in stock (inventory turns). A university admitting 1,000 students per year who each stay 4 years has about 4,000 students enrolled.
Example: checking Little's law with independent measurements. In the code, $\lambda$ and $W$ come from the customers' records, and $L$ from counting customers at 20,000 random instants. The second system has uniform arrivals, lognormal service and two servers, so nothing in it is Markovian.
import heapq
import numpy as np
def run(interarrival, service, c, n, rng):
"""FCFS G/G/c: return arrival and departure times of every customer."""
a = np.cumsum(interarrival(rng, n))
s = service(rng, n)
free, d = [0.0] * c, np.empty(n)
for k in range(n):
start = max(a[k], heapq.heappop(free))
d[k] = start + s[k]
heapq.heappush(free, d[k])
return a, d
rng = np.random.default_rng(seed=13)
cases = {
"M/M/1 coffee bar": (lambda r, n: r.exponential(3, n), # minutes
lambda r, n: r.exponential(2.5, n), 1),
"U/LN/2 workshop": (lambda r, n: r.uniform(0, 4, n), # nothing Markovian
lambda r, n: r.lognormal(0.8, 0.8, n), 2),
}
for name, (A, S, c) in cases.items():
a, d = run(A, S, c, 200_000, rng)
T = d.max() # observe [0, T]: starts and ends empty
lam = len(a) / T # throughput (customers per minute)
W = np.mean(d - a) # average time in system per customer
# L measured independently: count customers present at 20,000 random instants
probes = np.sort(rng.uniform(0, T, 20_000))
L = (np.searchsorted(np.sort(a), probes) - np.searchsorted(np.sort(d), probes)).mean()
print(f"{name}: lambda = {lam:.4f}/min, W = {W:6.2f} min,"
f" lambda*W = {lam * W:6.3f}, L (random inspection) = {L:6.3f}")
# M/M/1 coffee bar: lambda = 0.3334/min, W = 14.90 min, lambda*W = 4.968, L (random inspection) = 4.992
# U/LN/2 workshop: lambda = 0.4993/min, W = 5.40 min, lambda*W = 2.694, L (random inspection) = 2.677
In both systems, the product of throughput and time in system matches the number of customers seen by an outside observer, within sampling noise. In practice the law runs in whichever direction is convenient. A manager who counts 5 customers inside on average and 20 arrivals per hour knows, without timing anyone, that the average visit lasts 15 minutes.
Key takeaways
- $L = \lambda W$ holds for long-run averages in essentially any stable system, regardless of distributions, servers or discipline.
- It follows from counting the area under $N(t)$ in two ways.
- Apply it to any box: queue ($L_q = \lambda W_q$), servers ($\lambda/\mu$ busy on average), stations, whole networks.
- Use the rate of customers who actually enter. It also covers WIP, inventory and any stock-and-flow system.
Practice
- (Beginner) An emergency department receives 12 patients per hour and holds 30 patients on average. What is the average stay? If patients spend 40% of their stay waiting for a bed after treatment, how many are boarding on average?
- (Intermediate) A warehouse holds inventory worth 2.4 million EUR on average, and its annual cost of goods sold is 18 million EUR. Compute the average time an item spends in stock and the inventory turns. Which assumption of Little's law is questionable in a fast-growing business, and how does it bias the estimate?
- (Interview) "A service handles 300 requests per second, and each request holds a database connection for 40 ms on average. How large should the connection pool be?" Hint: Little's law gives the average number of connections in use. Then size for variability with an M/M/c model (next topics) and a target probability of waiting for a connection.
M/M/1 Queue
Why it matters. The M/M/1 queue is the simplest queue with a closed-form solution. It also teaches the most important qualitative lesson of queueing theory: waiting grows nonlinearly and explodes as utilization approaches 1. It is the baseline against which multi-server, finite and non-exponential systems are compared.
Theory. Assumptions: Poisson arrivals at rate $\lambda$, exponential service at rate $\mu$, one server, unlimited waiting room, FCFS, and $\rho = \lambda/\mu < 1$. The number in the system is a birth-death process with $\lambda_n = \lambda$ and $\mu_n = \mu$ (Section 9.2). The cut equations $\lambda p_n = \mu p_{n+1}$ give a geometric distribution:
$$ p_n = (1 - \rho),\rho^n, \qquad P(N \ge n) = \rho^n, \qquad P(\text{server busy}) = \rho. $$
Summing the geometric series and applying Little's law:
$$ L = \frac{\rho}{1 - \rho}, \qquad L_q = \frac{\rho^2}{1 - \rho}, \qquad W = \frac{1}{\mu - \lambda}, \qquad W_q = \frac{\rho}{\mu - \lambda}. $$
Under FCFS the time in system is itself exponential with rate $\mu - \lambda$, so $P(W > t) = e^{-(\mu - \lambda)t}$. The waiting time has an atom at zero, $P(W_q = 0) = 1 - \rho$, and $P(W_q > t) = \rho, e^{-(\mu - \lambda)t}$.
The nonlinearity. Write the time in system as
$$ W = \frac{1}{\mu} \cdot \frac{1}{1 - \rho}, $$
a service time multiplied by a congestion factor $1/(1 - \rho)$. At 50% utilization a visit takes twice the service time, at 90% ten times, and at 99% a hundred times. The last few percent of utilization cost the most waiting. Picture a highway: traffic flows freely up to a point, then a few more cars produce a jam.
Two more facts warn against relying blindly on the formulas near $\rho = 1$. The variance of $N$, $\rho/(1-\rho)^2$, grows even faster than its mean, so the queue length fluctuates wildly. The time needed to reach steady state also grows without bound, so heavily loaded systems spend much of their life in transient states.
Example: the coffee bar as an M/M/1 queue.
import numpy as np
# simulate_queue: defined in the Kendall Notation example
def mm1(lam, mu):
"""M/M/1 steady-state measures (same time unit for lam and mu)."""
if lam >= mu:
raise ValueError("unstable: need lam < mu")
rho = lam / mu
L = rho / (1 - rho)
Lq = rho ** 2 / (1 - rho)
return dict(rho=rho, L=L, Lq=Lq, W=L / lam, Wq=Lq / lam, P0=1 - rho)
lam, mu = 20, 24 # per hour
m = mm1(lam, mu)
print({k: round(v, 3) for k, v in m.items()})
print(f"W = {60 * m['W']:.1f} min, Wq = {60 * m['Wq']:.1f} min")
print(f"P(time in system > 30 min) = {np.exp(-(mu - lam) * 0.5):.3f}")
print(f"P(at least 5 in system) = {m['rho'] ** 5:.3f}")
print("\n lam rho L Wq(min)")
for lam_i in (12, 16, 20, 22, 23, 23.5):
r = mm1(lam_i, mu)
print(f"{lam_i:5} {r['rho']:.3f} {r['L']:6.2f} {60 * r['Wq']:7.1f}")
wq, w = simulate_queue(lambda g, n: g.exponential(3, n), lambda g, n: g.exponential(2.5, n),
n=1_000_000, warmup=50_000)
print(f"\nsimulation (lam=20, mu=24): Wq = {wq:.2f} min, W = {w:.2f} min")
# {'rho': 0.833, 'L': 5.0, 'Lq': 4.167, 'W': 0.25, 'Wq': 0.208, 'P0': 0.167}
# W = 15.0 min, Wq = 12.5 min
# P(time in system > 30 min) = 0.135
# P(at least 5 in system) = 0.402
#
# lam rho L Wq(min)
# 12 0.500 1.00 2.5
# 16 0.667 2.00 5.0
# 20 0.833 5.00 12.5
# 22 0.917 11.00 27.5
# 23 0.958 23.00 57.5
# 23.5 0.979 47.00 117.5
#
# simulation (lam=20, mu=24): Wq = 12.64 min, W = 15.14 min
A customer needs 2.5 minutes of service but spends 15 minutes in the coffee bar, and 13.5% of customers stay more than half an hour. The table shows the nonlinearity. Going from 12 to 20 customers per hour multiplies the wait by 5. The next 3.5 customers per hour multiply it by almost 10 more, although the barista's utilization rises only from 83% to 98%. The simulation agrees with the formulas.
Key takeaways
- M/M/1 is a birth-death process with constant rates. The number in system is geometric: $p_n = (1-\rho)\rho^n$.
- $L = \rho/(1-\rho)$, $W = 1/(\mu - \lambda)$, $L_q = \rho^2/(1-\rho)$, $W_q = \rho/(\mu - \lambda)$, and the time in system is exponential with rate $\mu - \lambda$.
- Time in system = service time × $1/(1-\rho)$: delays explode as $\rho \to 1$.
- Near saturation, variability and transient effects dominate, so treat steady-state numbers with care.
Practice
- (Beginner) A shared printer receives 18 jobs per hour, and print times are exponential with mean 3 minutes. Compute $\rho$, $L$, $L_q$, $W$, $W_q$ and the probability that more than 3 jobs are in the system.
- (Intermediate) Waiting costs 0.50 EUR per customer-minute in the system (lost goodwill), and a server with rate $\mu$ customers per hour costs $2\mu$ EUR per hour. For $\lambda = 20$ per hour, write the total cost per hour as a function of $\mu$ and show that the optimal rate is $\mu^* = \lambda + \sqrt{\lambda, C_w / C_s}$, with $C_w$ the waiting cost per customer-hour and $C_s$ the cost per unit of rate. Compute $\mu^*$.
- (Interview) "Traffic to our service rose by 10% and average latency doubled. Is something broken?" Hint: compute the ratio of $W$ values when $\rho$ goes from 0.85 to 0.935 in an M/M/1 model.
M/M/c Queue
Why it matters. Most service systems have several parallel servers fed by one line: call centers, bank tellers, a single snake line for many checkouts, hospital beds, pools of servers in a data center. The M/M/c model, through the Erlang C formula, is the standard tool for staffing them. It also quantifies the economies of scale from pooling capacity.
Theory. Poisson arrivals at rate $\lambda$, $c$ identical exponential servers with rate $\mu$ each, one FCFS queue. The number in the system is a birth-death process with $\lambda_n = \lambda$ and $\mu_n = \min(n, c),\mu$: below $c$, only $n$ servers are working. Define
$$ a = \frac{\lambda}{\mu} \quad (\text{offered load in Erlangs, the mean number of busy servers}), \qquad \rho = \frac{a}{c} < 1. $$
The product formula gives
$$ p_0 = \Big[\sum_{n=0}^{c-1}\frac{a^n}{n!} + \frac{a^c}{c!,(1 - \rho)}\Big]^{-1}, \qquad p_n = \begin{cases} \dfrac{a^n}{n!},p_0 & n \le c, \[2mm] \dfrac{a^n}{c!,c^{,n-c}},p_0 & n \ge c. \end{cases} $$
The key quantity is the probability that an arriving customer must wait, the Erlang C formula:
$$ C(c, a) = P(\text{wait}) = \sum_{n \ge c} p_n = \frac{a^c}{c!,(1 - \rho)},p_0. $$
Everything else follows from it:
$$ L_q = C(c,a),\frac{\rho}{1 - \rho}, \qquad W_q = \frac{C(c,a)}{c\mu - \lambda}, \qquad W = W_q + \frac{1}{\mu}, \qquad L = L_q + a, $$
$$ P(W_q > t) = C(c,a), e^{-(c\mu - \lambda)t}. $$
Factorials overflow for large $c$, so software computes Erlang C from the Erlang B recursion $B_0 = 1$, $B_k = aB_{k-1}/(k + aB_{k-1})$, then $C = cB_c / (c - a(1 - B_c))$.
Pooling. Two baristas, each with a separate line, form two M/M/1 queues. One shared line feeding both is an M/M/2 queue with the same total capacity and the same utilization, and it waits much less. With separate lines, one barista can be idle while customers wait at the other. Pooling eliminates that waste. In general, at a fixed utilization, more pooled servers mean less waiting. That is why large call centers can run at 90% or more with short waits, while a small one at the same utilization cannot.
Square-root staffing. For large offered loads, the number of servers needed for a given service level is approximately
$$ c \approx a + \beta\sqrt{a}, $$
where $\beta$ depends only on the target (the Halfin-Whitt or quality-and-efficiency-driven regime). The safety capacity grows like $\sqrt{a}$, not like $a$, so utilization $a/c \to 1$ as systems grow.
Example: pooling, and staffing a call center. First, the coffee bar at 40 customers per hour three ways: two separate M/M/1 lines, one pooled line with two baristas, and one barista twice as fast. Then, a call center with 300 calls per hour and a 4-minute mean handling time, which must answer at least 80% of calls within 20 seconds.
import numpy as np
def erlang_b(c, a):
"""Blocking probability of M/M/c/c with offered load a = lam/mu (stable recursion)."""
b = 1.0
for k in range(1, c + 1):
b = a * b / (k + a * b)
return b
def mmc(lam, mu, c):
"""M/M/c steady-state measures. Requires lam < c*mu."""
a, rho = lam / mu, lam / (c * mu) # offered load (Erlangs), utilization
if rho >= 1:
raise ValueError("unstable")
B = erlang_b(c, a)
C = c * B / (c - a * (1 - B)) # Erlang C: P(arrival must wait)
Lq = C * rho / (1 - rho)
Wq = Lq / lam
return dict(rho=rho, Pwait=C, Lq=Lq, Wq=Wq, W=Wq + 1 / mu, L=Lq + a)
# Pooling: total 40 customers/hour, service rate 24/hour per barista
sep = mmc(20, 24, 1) # two separate lines, one barista each
pool = mmc(40, 24, 2) # one shared line, two baristas
fast = mmc(40, 48, 1) # one barista twice as fast
for name, r in [("2 x M/M/1", sep), ("M/M/2 pooled", pool), ("M/M/1 fast", fast)]:
print(f"{name:>13}: P(wait) = {r['Pwait']:.3f}, Wq = {60 * r['Wq']:5.2f} min,"
f" W = {60 * r['W']:5.2f} min")
# Call-center staffing: 300 calls/hour, 4-minute mean handling time (a = 20 Erlangs).
# Target: at least 80% of calls answered within 20 seconds.
lam, mu, target_t = 300 / 60, 1 / 4, 20 / 60 # per minute, minutes
print("\n c rho P(wait) Wq(s) P(Wq <= 20 s)")
for c in range(21, 27):
r = mmc(lam, mu, c)
sl = 1 - r["Pwait"] * np.exp(-(c * mu - lam) * target_t) # waiting-time tail
print(f"{c:3} {r['rho']:.3f} {r['Pwait']:.3f} {60 * r['Wq']:5.1f} {sl:.3f}")
a = lam / mu
print(f"square-root rule c = a + beta*sqrt(a), beta = 1: {a + np.sqrt(a):.1f}")
# 2 x M/M/1: P(wait) = 0.833, Wq = 12.50 min, W = 15.00 min
# M/M/2 pooled: P(wait) = 0.758, Wq = 5.68 min, W = 8.18 min
# M/M/1 fast: P(wait) = 0.833, Wq = 6.25 min, W = 7.50 min
#
# c rho P(wait) Wq(s) P(Wq <= 20 s)
# 21 0.952 0.761 182.6 0.300
# 22 0.909 0.568 68.1 0.519
# 23 0.870 0.416 33.3 0.676
# 24 0.833 0.298 17.9 0.786
# 25 0.800 0.209 10.0 0.862
# 26 0.769 0.143 5.7 0.913
# square-root rule c = a + beta*sqrt(a), beta = 1: 24.5
Pooling the two lines cuts the time in system from 15 to 8.2 minutes with no extra staff. A single barista working twice as fast is better still on total time (7.5 minutes), because nobody's service is slowed down. The pooled system, however, has a shorter queue wait (5.7 against 6.25 minutes). That explains the M/M/2 surprise in the Kendall Notation example. In the call center, 21 agents keep the system stable ($\rho = 0.95$) but answer only 30% of calls within 20 seconds. Meeting the 80/20 target needs 25 agents: the 20 Erlangs of work plus 5 agents of safety capacity, close to the square-root rule. Notice how steep the table is: near the minimum, each extra agent roughly halves the average wait (183, 68, 33, 18 seconds).
Key takeaways
- M/M/c is a birth-death process with $\mu_n = \min(n, c)\mu$. It needs $\rho = \lambda/(c\mu) < 1$.
- Erlang C gives the probability of waiting. Then $L_q = C\rho/(1-\rho)$, $W_q = C/(c\mu - \lambda)$, and $P(W_q > t) = C e^{-(c\mu - \lambda)t}$.
- Pooling servers behind one line sharply reduces waiting at the same utilization (economies of scale).
- Staffing for a service level: $c \approx a + \beta\sqrt{a}$. Larger systems run at higher utilization for the same service.
Practice
- (Beginner) A bank has 3 tellers, 40 customer arrivals per hour, and exponential service with mean 4 minutes. Compute $\rho$, the probability of waiting, $L_q$, $W_q$ and $W$.
- (Intermediate) The call center's arrival rate changes during the day: 120, 300 and 450 calls per hour in the morning, midday and evening. Find the minimum number of agents in each period for the 80/20 target, and compare with $a + \sqrt{a}$. Why does the required utilization rise with the load?
- (Interview) "Should a hospital merge two separate 10-bed units into one 20-bed unit?" Hint: Erlang formulas show the pooling gain. Also discuss what pooling can lose: specialized staff, longer service times for mixed patients, travel distances, and the risk that one patient group crowds out the other.
Finite Capacity Queues
Why it matters. Real systems have limited space: telephone lines, buffer slots, parking spaces, beds, waiting rooms. When the system is full, arrivals are blocked and lost. The question then changes from "how long do customers wait?" to "how many customers do we lose?". A finite system also has a steady state even when it is overloaded.
Theory. In an M/M/c/K queue, at most $K \ge c$ customers can be present: $c$ in service and $K - c$ waiting. Arrivals finding $K$ customers are lost. It is a birth-death process with $\lambda_n = \lambda$ for $n < K$, $\lambda_K = 0$ and $\mu_n = \min(n, c)\mu$. The state space is finite, so the product formula always gives a steady state, whatever the value of $\lambda/(c\mu)$.
For one server (M/M/1/K) with $\rho = \lambda/\mu$:
$$ p_n = \frac{(1 - \rho),\rho^n}{1 - \rho^{K+1}} \quad (\rho \ne 1), \qquad p_n = \frac{1}{K + 1} \quad (\rho = 1), \qquad n = 0, \dots, K. $$
Key measures:
- Blocking probability. By PASTA, an arrival finds the system full with probability $p_K$.
- Effective arrival rate. $\lambda_{\text{eff}} = \lambda(1 - p_K)$ customers per unit time actually enter. Little's law must use it: $W = L/\lambda_{\text{eff}}$ and $W_q = L_q/\lambda_{\text{eff}}$.
- Utilization $= \lambda_{\text{eff}}/(c\mu)$ (for one server, $1 - p_0$).
- Lost customers per unit time $= \lambda p_K$.
With no waiting room ($K = c$), we get the Erlang loss model M/M/c/c. Its blocking probability is the Erlang B formula:
$$ B(c, a) = \frac{a^c / c!}{\sum_{k=0}^{c} a^k / k!}, \qquad a = \frac{\lambda}{\mu}. $$
Erlang B is insensitive: it holds for M/G/c/c with any service-time distribution of mean $1/\mu$. This is why it is used for hospital beds, parking spaces and rental fleets, where stay lengths are far from exponential. It was originally developed to size telephone trunk lines.
Capacity creates a trade-off: a larger $K$ loses fewer customers, but those who enter wait longer.
Example: a coffee bar with limited room, and hospital beds. The code reuses birth_death_stationary from Section 9.2.
import numpy as np
def birth_death_stationary(lam, mu): # from Section 9.2
w = np.concatenate(([1.0], np.cumprod(np.asarray(lam) / np.asarray(mu))))
return w / w.sum()
def mmck(lam, mu, c, K):
"""M/M/c/K: at most K customers in the system (c in service, K - c waiting)."""
p = birth_death_stationary([lam] * K, [min(n, c) * mu for n in range(1, K + 1)])
n = np.arange(K + 1)
lam_eff = lam * (1 - p[K]) # arrivals that actually enter
L = p @ n
Lq = p @ np.maximum(n - c, 0)
return dict(block=p[K], lam_eff=lam_eff, L=L, Lq=Lq, W=L / lam_eff, Wq=Lq / lam_eff)
# Coffee bar, one barista, room for K customers inside (lam = 20/h, mu = 24/h)
print(" K P(block) lost/h L W(min)")
for K in (2, 4, 6, 10, 20):
r = mmck(20, 24, 1, K)
print(f"{K:3} {r['block']:.3f} {20 - r['lam_eff']:5.2f} {r['L']:5.2f} {60 * r['W']:6.2f}")
# Overloaded but finite: lam = 30/h > mu = 24/h is still a valid steady state
r = mmck(30, 24, 1, 6)
print(f"\noverloaded, K=6: P(block) = {r['block']:.3f}, throughput = {r['lam_eff']:.2f}/h")
# Erlang B (M/M/c/c, no waiting room): hospital ward, 2 admissions/day, 4-day stay
def erlang_b(c, a):
b = 1.0
for k in range(1, c + 1):
b = a * b / (k + a * b)
return b
a = 2 * 4 # offered load: 8 beds busy on average
for c in (8, 10, 12, 14):
print(f"beds = {c:2}: P(patient turned away) = {erlang_b(c, a):.3f},"
f" bed occupancy = {a * (1 - erlang_b(c, a)) / c:.3f}")
# K P(block) lost/h L W(min)
# 2 0.275 5.49 0.88 3.64
# 4 0.134 2.69 1.64 5.69
# 6 0.077 1.55 2.29 7.45
# 10 0.031 0.62 3.29 10.18
# 20 0.004 0.09 4.53 13.66
#
# overloaded, K=6: P(block) = 0.253, throughput = 22.41/h
# beds = 8: P(patient turned away) = 0.236, bed occupancy = 0.764
# beds = 10: P(patient turned away) = 0.122, bed occupancy = 0.703
# beds = 12: P(patient turned away) = 0.051, bed occupancy = 0.632
# beds = 14: P(patient turned away) = 0.017, bed occupancy = 0.562
With room for only 2 customers, the coffee bar turns away 27.5% of them, about 5.5 per hour, but those inside stay less than 4 minutes. With room for 20, losses almost vanish and the time in system approaches the M/M/1 value of 15 minutes. The overloaded case ($\lambda > \mu$) has a perfectly good steady state: a quarter of arrivals are blocked, and throughput (22.4 per hour) stays below the service rate because the barista is sometimes idle. In the hospital, providing exactly the average load of 8 beds turns away almost a quarter of patients. Keeping refusals below 5% needs 13 beds (3.1% blocked), at an occupancy of about 60%. A target of "95% bed occupancy" is incompatible with low refusal rates in small units.
Key takeaways
- M/M/c/K is a truncated birth-death process. Its steady state exists for any load, because the state space is finite.
- Blocking probability is $p_K$ (PASTA). Use $\lambda_{\text{eff}} = \lambda(1 - p_K)$ in Little's law.
- Erlang B gives the loss probability with no waiting room, and it is insensitive to the service-time distribution.
- More room means fewer losses but longer waits. Small loss systems need substantial idle capacity.
Practice
- (Beginner) A drive-through has space for 5 cars including the one at the window. Arrivals are Poisson at 30 per hour, and service is exponential at 36 per hour. Compute the fraction of customers lost, the throughput and the average time in system.
- (Intermediate) For the coffee bar ($\lambda = 20$, $\mu = 24$ per hour), each lost customer costs 5 EUR of profit and each place of capacity costs 3 EUR per hour (space rent). Use
mmckto find the capacity $K$ that minimizes total cost per hour.- (Interview) "Why does an overloaded finite queue ($\lambda > \mu$) have a steady state, while an infinite one does not? What does it look like?" Hint: the buffer is almost always full. Look at the throughput and at the limit of $p_K$ as $K$ grows: compare it with $1 - \mu/\lambda$.
M/G/1 and Pollaczek-Khinchine Formula
Why it matters. Service times are rarely exponential. A machine cycle is nearly constant, while a repair or a support ticket can be highly variable. With Poisson arrivals, the Pollaczek-Khinchine (P-K) formula gives the exact mean wait for any service distribution. It shows that the wait depends on the service distribution only through its mean and variance, or equivalently through $c_s^2$.
Theory. M/G/1: Poisson arrivals at rate $\lambda$, i.i.d. service times $S$ with mean $E[S] = 1/\mu$ and variance $\sigma^2$, one server, FCFS, and $\rho = \lambda E[S] < 1$. The number in system is no longer a Markov chain in continuous time, since the elapsed service matters. Its mean, however, follows from a short argument.
Mean-value derivation. By PASTA, an arriving customer sees the time-average state. It must wait for
- the $L_q$ customers already in the queue, each needing $E[S]$ on average,
- the residual service of the customer in service, if any. The server is busy with probability $\rho$, and the mean residual is $E[R] = E[S^2]/(2E[S])$.
The residual is longer than $E[S]/2$ because of the inspection paradox: a random arrival is more likely to fall inside a long service than a short one. (In the same way, if you arrive at a bus stop at a random time, you are more likely to land in a long gap between buses.) So
$$ W_q = L_q,E[S] + \rho,E[R] = \rho W_q + \rho,\frac{E[S^2]}{2E[S]} \qquad (\text{using } L_q = \lambda W_q), $$
which gives the Pollaczek-Khinchine formula:
$$ W_q = \frac{\lambda,E[S^2]}{2(1 - \rho)}, \qquad L_q = \lambda W_q = \frac{\lambda^2\sigma^2 + \rho^2}{2(1 - \rho)}. $$
With $E[S^2] = (1 + c_s^2)E[S]^2$, the formula separates three effects:
$$ W_q = \underbrace{\frac{1 + c_s^2}{2}}{\text{variability}} \times \underbrace{\frac{\rho}{1 - \rho}}{\text{utilization}} \times \underbrace{E[S]}_{\text{time}}. $$
| Service | $c_s^2$ | $W_q$ relative to M/M/1 |
|---|---|---|
| deterministic (M/D/1) | 0 | half |
| Erlang-$k$ (M/$E_k$/1) | $1/k$ | $(1 + 1/k)/2$ |
| exponential (M/M/1) | 1 | same |
| hyperexponential, heavy-tailed | > 1 | larger, proportional to $1 + c_s^2$ |
For general arrivals (G/G/1), Kingman's approximation uses the same structure with both variabilities:
$$ W_q \approx \frac{c_a^2 + c_s^2}{2} \cdot \frac{\rho}{1 - \rho} \cdot E[S]. $$
It is exact for M/G/1 and accurate in heavy traffic. Multi-server versions such as the Allen-Cunneen approximation scale the M/M/c wait by the same factor $(c_a^2 + c_s^2)/2$. This "variability × utilization × time" (VUT) decomposition is the standard rule of thumb in factory physics. To cut waiting, reduce utilization or reduce variability; the second is often cheaper.
Example: four espresso processes with the same mean. The coffee bar keeps $\lambda = 20$ per hour and $E[S] = 2.5$ minutes ($\rho = 5/6$), with four service-time distributions. The code compares the P-K formula with simulation.
import numpy as np
# simulate_queue: defined in the Kendall Notation example
lam, ES = 1 / 3, 2.5 # per minute; mean service 2.5 min (rho = 5/6)
def pk_wq(lam, ES, scv):
"""Pollaczek-Khinchine mean wait: Wq = lam E[S^2] / (2 (1 - rho))."""
rho = lam * ES
ES2 = (1 + scv) * ES ** 2 # E[S^2] = Var + mean^2 = (1 + CV^2) mean^2
return lam * ES2 / (2 * (1 - rho))
services = { # name: (sampler, squared CV)
"deterministic": (lambda g, n: np.full(n, ES), 0.0),
"Erlang-4": (lambda g, n: g.gamma(4, ES / 4, n), 0.25),
"exponential": (lambda g, n: g.exponential(ES, n), 1.0),
"hyperexp.": (lambda g, n: np.where(g.random(n) < 0.9, # 90% quick, 10% long
g.exponential(1.25, n),
g.exponential(13.75, n)), 5.5),
}
arrivals = lambda g, n: g.exponential(1 / lam, n)
print(f"{'service':>14} {'CV^2':>5} {'Wq (P-K)':>9} {'Wq (sim)':>9} {'Lq (P-K)':>9}")
for name, (S, scv) in services.items():
wq = pk_wq(lam, ES, scv)
sim, _ = simulate_queue(arrivals, S, n=2_000_000, warmup=100_000, seed=17)
print(f"{name:>14} {scv:5.2f} {wq:9.2f} {sim:9.2f} {lam * wq:9.2f}")
# service CV^2 Wq (P-K) Wq (sim) Lq (P-K)
# deterministic 0.00 6.25 6.28 2.08
# Erlang-4 0.25 7.81 7.84 2.60
# exponential 1.00 12.50 12.40 4.17
# hyperexp. 5.50 40.62 40.22 13.54
The formula matches the simulation for every distribution. A perfectly consistent barista halves the M/M/1 wait (6.25 against 12.5 minutes). A barista who is usually quick but occasionally very slow (10% of orders are complex drinks) more than triples it, to 40 minutes, with the same average speed. Here, standardizing the complex drinks or routing them to a separate station would help more than hiring a faster barista.
Key takeaways
- P-K formula: $W_q = \lambda E[S^2] / (2(1 - \rho))$. It holds for M/G/1 with any service distribution.
- Waiting depends on service variability through $1 + c_s^2$. M/D/1 waits half as long as M/M/1.
- The residual-service argument (inspection paradox) explains why variability hurts.
- Kingman's VUT approximation $W_q \approx \frac{c_a^2 + c_s^2}{2}\cdot\frac{\rho}{1-\rho}\cdot E[S]$ extends the insight to general arrivals.
Practice
- (Beginner) An automated car wash takes exactly 6 minutes per car, and cars arrive as a Poisson process at 8 per hour. Compute $L_q$, $W_q$ and $W$, then repeat with exponential wash times of the same mean.
- (Intermediate) A repair shop's jobs take 30 minutes with probability 0.8 and 2 hours with probability 0.2, and jobs arrive as a Poisson process at 0.9 per hour. Compute $E[S]$, $E[S^2]$, $\rho$ and $W_q$. Which helps more: cutting long jobs to 1.5 hours, or making them half as frequent (probability 0.1, with the short jobs taking up the rest)?
- (Interview) "Your team can either make the average service 10% faster or eliminate the slowest outliers, which halves the variance. Which reduces waiting more?" Hint: use the VUT decomposition. A faster service lowers both $\rho$ and $E[S]$ (and raises $c_s^2$ if the variance stays the same), while halving the variance halves $c_s^2$. Compare at $\rho = 0.8$ and $\rho = 0.95$.
Queueing Networks
Why it matters. Real processes have several stages. A patient visits reception, triage, a doctor and a lab. A job flows through several machines. A web request passes through a load balancer, an application server and a database. Network models find each station's load, the bottleneck, and the maximum throughput of the whole process using only flow conservation, before any waiting-time formula is needed.
Theory. A queueing network is a set of stations $j = 1, \dots, J$. Each has $c_j$ servers with rate $\mu_j$, and customers move between stations after service.
- In an open network, customers arrive from outside (rate $\gamma_j$ at station $j$) and eventually leave.
- In a closed network, a fixed population circulates forever: pallets in a manufacturing loop, a fixed number of concurrent users, a fleet of machines.
With Markovian routing, a customer leaving station $i$ goes to station $j$ with probability $r_{ij}$, or leaves with probability $1 - \sum_j r_{ij}$, independently of everything else. The traffic equations state that, in steady state, flow into a station equals flow out:
$$ \lambda_j = \gamma_j + \sum_i \lambda_i, r_{ij}, \qquad \text{in matrix form} \quad \lambda = \gamma + R^\top\lambda ;\Rightarrow; \lambda = (I - R^\top)^{-1}\gamma. $$
The inverse exists if every customer eventually leaves. The solution gives, for each station:
- the visit ratio $v_j = \lambda_j / \sum_k \gamma_k$, the mean number of visits per customer (it can exceed 1 with rework or feedback),
- the utilization $\rho_j = \lambda_j / (c_j\mu_j)$. Every station needs $\rho_j < 1$.
The station with the largest $\rho_j$ is the bottleneck. If all external arrivals are scaled up by the same factor, the bottleneck saturates first. The maximum sustainable external rate is therefore $\gamma_{\max} = \sum_k \gamma_k / \max_j \rho_j$. Adding capacity anywhere other than the bottleneck does not raise it.
Two structural facts help:
- Tandem lines. By Burke's theorem, the departure process of a stable M/M/c queue in steady state is Poisson with rate $\lambda$. A series of M/M/c stations therefore decomposes exactly into independent M/M/c queues, each fed by a Poisson stream.
- Feedback. With rework or loops, the stream entering a station is no longer Poisson. Customers who loop back create correlated arrivals. The next topic shows that the M/M/c formulas still give the right steady-state answers in Jackson networks.
Computer performance engineering uses the same ideas as operational laws: the utilization law $U_j = X,D_j$ (system throughput $X$ times service demand per customer) and the forced-flow law $X_j = v_j X$.
Example: an outpatient clinic. Patients arrive at reception at 30 per hour and then go to triage. 60% see a doctor and 40% a nurse. 20% of doctor visits require a lab test, after which the patient returns to the doctor.
0.6 +--------------+ 0.8
gamma = 30/h +-------------+ +-----------+ -----> | doctor (c=5) | -----> exit
------------> | reception | -> | triage | +--------------+
| (c=1) | | (c=2) | ^ | 0.2
+-------------+ +-----------+ 1.0 | v
| 0.4 +-------------+
v | lab (c=1) |
+-------------+ +-------------+
| nurse (c=2) | -----> exit
+-------------+
import numpy as np
names = ["reception", "triage", "doctor", "nurse", "lab"]
gamma = np.array([30.0, 0, 0, 0, 0]) # external arrivals per hour (all at reception)
R = np.array([[0, 1.0, 0, 0, 0 ], # R[i, j] = P(go to j after i)
[0, 0, 0.6, 0.4, 0 ], # triage: 60% doctor, 40% nurse
[0, 0, 0, 0, 0.2], # doctor: 20% need a lab test, 80% leave
[0, 0, 0, 0, 0 ], # nurse: leave
[0, 0, 1.0, 0, 0 ]]) # lab: back to the doctor
mu = np.array([40.0, 18, 6, 8, 9]) # service rate per server (per hour)
c = np.array([1, 2, 5, 2, 1]) # servers per station
# Traffic equations: lam = gamma + R^T lam -> (I - R^T) lam = gamma
lam = np.linalg.solve(np.eye(5) - R.T, gamma)
rho = lam / (c * mu)
print(f"{'station':>10} {'lam':>6} {'visits':>7} {'rho':>6}")
for j in range(5):
print(f"{names[j]:>10} {lam[j]:6.2f} {lam[j] / gamma.sum():7.3f} {rho[j]:6.3f}")
b = np.argmax(rho)
print(f"bottleneck: {names[b]}; max external rate before saturation = {gamma.sum() / rho[b]:.1f}/h")
print("exit probabilities 1 - row sums:", 1 - R.sum(axis=1))
# station lam visits rho
# reception 30.00 1.000 0.750
# triage 30.00 1.000 0.833
# doctor 22.50 0.750 0.750
# nurse 12.00 0.400 0.750
# lab 4.50 0.150 0.500
# bottleneck: triage; max external rate before saturation = 36.0/h
# exit probabilities 1 - row sums: [0. 0. 0.8 1. 0. ]
Only 18 patients per hour are sent to the doctors, but the lab loop raises the doctors' load to 22.5 visits per hour. Each patient who sees a doctor does so 1.25 times on average ($0.75 / 0.6$). Triage is the bottleneck at 83% utilization, so the clinic can absorb at most 36 patients per hour. Adding a sixth doctor would shorten the doctors' queue but would not raise that limit; only triage capacity does.
Key takeaways
- Traffic equations $\lambda = \gamma + R^\top\lambda$ give every station's arrival rate from external arrivals and routing alone.
- Visit ratios count mean visits per customer, and feedback loops can push them above 1.
- The bottleneck is the station with the highest utilization. It alone limits network throughput: $\gamma_{\max} = \gamma / \max_j \rho_j$.
- Burke's theorem: a stable M/M/c queue emits a Poisson stream, so tandem M/M/c lines decompose exactly.
Practice
- (Beginner) A three-station production line receives 10 jobs per hour. After station 3 (inspection), 10% of jobs go back to station 2 for rework. Write $R$, solve the traffic equations, and compute utilizations for single servers with rates 12, 13 and 14 per hour.
- (Intermediate) The clinic manager can add one server to any single station. Using the code, find the addition that maximizes the sustainable external arrival rate, and report the new bottleneck.
- (Interview) "Each request to an API gateway calls service A twice on average and service B 0.5 times. A takes 20 ms per call on 4 instances, and B takes 50 ms per call on 2 instances. What is the maximum throughput of the system?" Hint: visit ratios give the service demand per request at each service. Divide each service's capacity by its demand and take the minimum.
Jackson Networks
Why it matters. Even when every station is Markovian, a network with feedback looks hopeless to analyze. Streams merge, split and loop back, and are no longer Poisson. Jackson's theorem says that the steady-state distribution nevertheless factorizes into independent M/M/c solutions. Open networks of any size can then be evaluated with the formulas already derived in this subchapter. It is one of the most useful results in applied probability.
Theory. An open Jackson network satisfies:
- Poisson external arrivals at each station (rate $\gamma_j$),
- exponential service at each station, with $c_j$ servers of rate $\mu_j$ (more generally, any queue-length-dependent rate),
- Markovian routing $r_{ij}$, independent of the state,
- unlimited buffers, FCFS, and $\rho_j = \lambda_j / (c_j\mu_j) < 1$ at every station, with $\lambda$ from the traffic equations.
Jackson's theorem (1957, 1963). The joint steady-state distribution of the numbers of customers at the stations is a product form:
$$ P(N_1 = n_1, \dots, N_J = n_J) = \prod_{j=1}^{J} p_j(n_j), $$
where $p_j$ is the stationary distribution of an isolated M/M/$c_j$ queue with arrival rate $\lambda_j$. In steady state, at any fixed time, the stations behave as if they were independent M/M/$c_j$ queues. The consequences:
- $L_j$ and $W_j$ (per visit) come from the M/M/c formulas with $\lambda_j$.
- The total number in the network is $L = \sum_j L_j$.
- The mean time a customer spends in the network is $W = L / \sum_j \gamma_j$ (Little's law). Equivalently, $W = \sum_j v_j W_j$: visits times time per visit.
What the theorem does not say. The stream entering a station with feedback is not Poisson, and a customer's sojourn times at successive stations can be correlated. The product form concerns the joint distribution of queue lengths at a random time. Its proof checks that the product satisfies the global balance equations $\pi Q = 0$ of the network's CTMC. The deeper reason is a property called quasi-reversibility of the stations.
Extensions.
- Closed Jackson networks (Gordon-Newell) also have a product form, normalized over all states with $N$ customers. Mean Value Analysis (MVA) computes them without the normalizing constant, using the arrival theorem: a customer arriving at a station sees the network in steady state with one fewer customer. For single-server stations with service demand $D_j = v_j/\mu_j$, iterate for $n = 1, \dots, N$:
$$ R_j(n) = D_j,\big(1 + L_j(n-1)\big), \qquad X(n) = \frac{n}{\sum_j R_j(n)}, \qquad L_j(n) = X(n),R_j(n). $$
- BCMP networks extend the product form to several customer classes and to processor-sharing, LCFS-preemptive and infinite-server stations.
- Without product form (general service, finite buffers, blocking, priorities), use decomposition approximations based on Kingman-type formulas, or simulation (Chapter 12).
Example: the clinic as a Jackson network, checked by simulation. The clinic of the previous topic has a feedback loop between doctor and lab, so the doctors' arrival stream is not Poisson. The product form still predicts each station's mean queue. An event-driven simulation of the whole network checks it.
import heapq, itertools
from collections import deque
import numpy as np
# mmc: defined in the M/M/c Queue example
names = ["reception", "triage", "doctor", "nurse", "lab"]
gamma = np.array([30.0, 0, 0, 0, 0])
R = np.array([[0, 1.0, 0, 0, 0], [0, 0, 0.6, 0.4, 0], [0, 0, 0, 0, 0.2],
[0, 0, 0, 0, 0], [0, 0, 1.0, 0, 0]])
mu = np.array([40.0, 18, 6, 8, 9]); c = np.array([1, 2, 5, 2, 1])
lam = np.linalg.solve(np.eye(5) - R.T, gamma)
# Jackson's theorem: each station behaves like an independent M/M/c with arrival rate lam_j
Lj = np.array([mmc(lam[j], mu[j], c[j])["L"] for j in range(5)])
W = Lj.sum() / gamma.sum() # Little's law for the whole network
print("product form: L per station", np.round(Lj, 3), f"total L {Lj.sum():.3f},"
f" W {60 * W:.2f} min")
def simulate(T=20_000.0, seed=19):
rng = np.random.default_rng(seed)
ev, seq = [], itertools.count()
push = lambda t, kind, j, cu: heapq.heappush(ev, (t, next(seq), kind, j, cu))
queue, busy = [deque() for _ in range(5)], np.zeros(5, int)
area, last, entry, soj = np.zeros(5), 0.0, {}, []
def join(t, j, cu): # customer cu arrives at station j
if busy[j] < c[j]:
busy[j] += 1
push(t + rng.exponential(1 / mu[j]), "dep", j, cu)
else:
queue[j].append(cu)
push(rng.exponential(1 / gamma.sum()), "arr", 0, 0)
while ev:
t, _, kind, j, cu = heapq.heappop(ev)
if t > T:
break
area += (busy + [len(q) for q in queue]) * (t - last); last = t
if kind == "arr": # external arrival (all enter at reception)
entry[cu] = t
join(t, 0, cu)
push(t + rng.exponential(1 / gamma.sum()), "arr", 0, cu + 1)
else: # service completion at station j
if queue[j]:
push(t + rng.exponential(1 / mu[j]), "dep", j, queue[j].popleft())
else:
busy[j] -= 1
k = rng.choice(6, p=np.append(R[j], 1 - R[j].sum())) # route, or exit (k = 5)
if k == 5:
soj.append(t - entry.pop(cu))
else:
join(t, k, cu)
return area / last, np.mean(soj)
L_sim, W_sim = simulate()
print("simulation : L per station", np.round(L_sim, 3), f"total L {L_sim.sum():.3f},"
f" W {60 * W_sim:.2f} min")
# product form: L per station [3. 5.455 5.135 3.429 1. ] total L 18.018, W 36.04 min
# simulation : L per station [3.009 5.461 5.104 3.529 1.016] total L 18.120, W 36.20 min
The product-form numbers agree with a 20,000-hour simulation within sampling error, including at the doctors' station, which sits inside the feedback loop. A patient spends 36 minutes in the clinic on average, and triage (the bottleneck) holds the most patients. With Jackson's theorem, the full analysis takes five M/M/c evaluations instead of a simulation of the whole network. That makes it practical to compare dozens of staffing plans in seconds and simulate only the best few.
Key takeaways
- In an open Jackson network, the steady-state joint distribution is the product of independent M/M/$c_j$ distributions with the rates $\lambda_j$ from the traffic equations.
- Station measures come from M/M/c formulas. Network totals follow from $L = \sum_j L_j$ and $W = L/\gamma$.
- Internal flows need not be Poisson, yet the product form still holds. It describes queue lengths at a random time.
- Closed networks are solved with Mean Value Analysis. Non-product-form networks need approximations or simulation.
Practice
- (Beginner) Two single-server stations in tandem receive $\gamma = 4$ customers per hour, with $\mu_1 = 6$ and $\mu_2 = 5$ per hour. Use Jackson's theorem to compute $L$, $W$ and the probability that both stations are empty.
- (Intermediate) Implement MVA for a closed network of three single-server stations with service demands $D = (0.05, 0.03, 0.04)$ seconds, for populations $N = 1, \dots, 20$. Plot throughput $X(N)$ and compare it with the asymptotic bound $1/\max_j D_j$. At which $N$ does the bound become nearly tight?
- (Interview) "How can Jackson's theorem hold when the arrival stream to a station with feedback is not Poisson?" Hint: the theorem concerns the stationary joint distribution of queue lengths, not the arrival process. Think about what verifying global balance requires, and about the arrival theorem.
<a id="94-stochastic-programming"></a>
9.4 Stochastic Programming
The models of Sections 9.1 to 9.3 describe how a random system behaves under given rules. Stochastic programming optimizes decisions that must be taken before the uncertainty is resolved. This subchapter builds two-stage models with recourse and represents uncertainty with scenarios. It then handles reliability requirements with chance constraints, and measures what a stochastic model is worth with the value of the stochastic solution. It closes with robust optimization, which protects against worst cases when probabilities are not trusted. The first, second and fourth topics share one example, the farmer problem of Birge and Louveaux.
Two-Stage Stochastic Programs with Recourse
Why it matters. Many decisions come in two steps. Capacity, planting or inventory is decided now, and operations are adjusted later, once demand, yields or prices are known. Optimizing for the average scenario ignores the cost of those adjustments. It produces plans that look good on paper but perform worse on average, the "flaw of averages" of Section 1.1. The two-stage model with recourse formalizes "decide, observe, react".
Theory. Decisions are split by when they are taken relative to the random event $\xi$:
decide x now nature reveals xi react with y(xi)
(here-and-now) ----> (yields, demand, prices) ----> (recourse: buy, sell, ...)
cost c^T x cost q^T y
The two-stage stochastic linear program with recourse is
$$ \min_{x}; c^\top x + \mathcal{Q}(x) \quad \text{s.t. } Ax = b,; x \ge 0, \qquad \mathcal{Q}(x) = E_\xi\big[Q(x, \xi)\big], $$
$$ Q(x, \xi) = \min_{y}; \big{, q(\xi)^\top y ;:; W y = h(\xi) - T(\xi),x,; y \ge 0 ,\big}. $$
- $x$ is the first-stage (here-and-now) decision. It cannot depend on $\xi$, a requirement called nonanticipativity.
- $y$ is the second-stage (recourse) decision, chosen after $\xi$ is known. There is one $y$ per outcome.
- $T(\xi)$ (the technology matrix) links the stages, $W$ is the recourse matrix, and $Q(x, \xi)$ is the recourse function: the cheapest way to react to outcome $\xi$ after deciding $x$.
Useful properties and terms:
- For each $\xi$, $Q(x, \xi)$ is the optimal value of an LP as a function of its right-hand side, so it is convex and piecewise linear in $x$ (Section 2.5). The expectation $\mathcal{Q}(x)$ is then convex too, and the whole problem is convex.
- Fixed recourse: $W$ does not depend on $\xi$. Complete recourse: the second stage is feasible for every right-hand side, for example because one can always buy the missing quantity. Relatively complete recourse: it is feasible for every first-stage-feasible $x$. Without these properties, some $x$ make $Q(x, \xi) = +\infty$ for some outcomes, and the model must exclude them.
- The newsvendor is the simplest two-stage program. The first stage is the order quantity, and the second stage sells $\min(\text{order}, \text{demand})$.
An analogy is packing for a trip. You choose what to pack before knowing the weather, and you can buy anything missing at your destination, but at tourist prices. Packing for "average weather" ignores how expensive those purchases are.
The farmer problem. A farmer has 500 acres for wheat, corn and sugar beets. At least 200 tons of wheat and 240 tons of corn are needed as cattle feed; shortfalls can be bought and surpluses sold. Sugar beets sell at 36 per ton up to a quota of 6,000 tons and at 10 per ton above it, and cannot be bought.
| Wheat | Corn | Sugar beets | |
|---|---|---|---|
| Planting cost (per acre) | 150 | 230 | 260 |
| Mean yield (tons per acre) | 2.5 | 3 | 20 |
| Feed requirement (tons) | 200 | 240 | - |
| Selling price (per ton) | 170 | 150 | 36 (first 6,000 t), 10 above |
| Purchase price (per ton) | 238 | 210 | - |
Yields depend on the weather. In a good year all yields are 20% above the mean, and in a bad year 20% below; good, average and bad years are equally likely. Planting (first stage) happens before the weather is known. Buying and selling (second stage) happen after the harvest.
Example: evaluating planting plans through the recourse function. Plan A is optimal if yields are always average. Plan B is a candidate that hedges against bad years. Evaluating a plan means solving one small second-stage LP per scenario.
import numpy as np
from scipy.optimize import linprog
# Farmer problem (Birge & Louveaux). First stage: acres x = (wheat, corn, beets), <= 500 total.
plant = np.array([150.0, 230.0, 260.0]) # planting cost per acre
# Second stage: y = (buy wheat, sell wheat, buy corn, sell corn, beets at 36, beets at 10)
q = np.array([238.0, -170.0, 210.0, -150.0, -36.0, -10.0]) # cost (negative = revenue)
need = np.array([200.0, 240.0]) # tons of wheat and corn needed for cattle feed
mean_yield = np.array([2.5, 3.0, 20.0]) # tons per acre
scenarios = np.array([1.2, 1.0, 0.8])[:, None] * mean_yield # good, average, bad years
probs = np.array([1, 1, 1]) / 3
def recourse(x, t):
"""Q(x, t): minimal second-stage cost for planting x and yields t (tons/acre)."""
A_ub = np.array([[-1, 1, 0, 0, 0, 0], # wheat: t1 x1 + buy - sell >= 200
[0, 0, -1, 1, 0, 0], # corn : t2 x2 + buy - sell >= 240
[0, 0, 0, 0, 1, 1], # beets sold <= beets grown
[0, 0, 0, 0, 1, 0]]) # quota: at most 6000 t at the high price
b_ub = np.array([t[0] * x[0] - need[0], t[1] * x[1] - need[1], t[2] * x[2], 6000])
return linprog(q, A_ub=A_ub, b_ub=b_ub, bounds=(0, None)).fun
plans = {"plan A (120, 80, 300)": np.array([120.0, 80, 300]),
"plan B (170, 80, 250)": np.array([170.0, 80, 250])}
for name, x in plans.items():
Q = np.array([recourse(x, t) for t in scenarios]) # one second-stage LP per scenario
profit = -(plant @ x + Q) # profit = -(total cost)
print(f"{name}: profit good/avg/bad = {np.round(profit).astype(int)},"
f" expected = {probs @ profit:,.0f}")
# plan A (120, 80, 300): profit good/avg/bad = [148000 118600 55120], expected = 107,240
# plan B (170, 80, 250): profit good/avg/bad = [167000 109350 48820], expected = 108,390
In an average year, plan A earns 118,600, far more than plan B's 109,350. But the average-year profit is not the average profit. In bad years plan A must buy expensive feed, and in good years it produces beets beyond the quota at the low price. Averaged over the three weather scenarios, plan B earns 108,390 against plan A's 107,240. The next topic shows that plan B is in fact the optimal two-stage plan.
Key takeaways
- A two-stage program chooses $x$ before uncertainty is revealed and $y(\xi)$ after, minimizing $c^\top x + E[Q(x, \xi)]$.
- The recourse function $Q(x, \xi)$ is the optimal cost of reacting to outcome $\xi$. For LPs it is convex and piecewise linear in $x$.
- Nonanticipativity: first-stage decisions cannot depend on the outcome.
- Check recourse properties (fixed, complete, relatively complete) to make sure every plan has a finite cost in every scenario.
Practice
- (Beginner) Write the newsvendor of Section 1.1 (buy at 2, sell at 10, leftovers worthless) as a two-stage program. Identify the first-stage variable, the second-stage variables (units sold, units left over), the recourse function, and whether the recourse is complete.
- (Intermediate) Using
recourse, plot the farmer's expected profit as a function of wheat acreage $x_1 \in [0, 250]$, with corn fixed at 80 acres and beets at $420 - x_1$. Verify that the curve is concave (convex expected cost) and locate its maximum.- (Interview) "Why not just solve the deterministic problem for each scenario and average the three plans?" Hint: the average of scenario-optimal plans ignores the cost of being wrong. Each scenario's plan assumes perfect foresight. Compare with nonanticipativity, and with the wait-and-see value (Value of the Stochastic Solution topic).
Scenario-Based Modeling
Why it matters. The expectation in a two-stage program is an integral over a possibly continuous distribution, which no LP solver can handle directly. Replacing the distribution with a finite set of scenarios turns the stochastic program into one large but ordinary LP, the deterministic equivalent. Almost all practical stochastic programs are solved this way, so building good scenarios is as important as the model itself.
Theory. With scenarios $\xi_1, \dots, \xi_S$ and probabilities $p_1, \dots, p_S$, the expectation becomes a sum, and each scenario gets its own copy $y_s$ of the second-stage variables. The extensive form is
$$ \min; c^\top x + \sum_{s=1}^{S} p_s, q_s^\top y_s \quad \text{s.t.}\quad Ax = b, \quad T_s x + W y_s = h_s ;; (s = 1, \dots, S), \quad x, y_s \ge 0. $$
Its constraint matrix has a characteristic block structure:
x y_1 y_2 ... y_S
[ A ] = b first stage
[ T_1 W ] = h_1 scenario 1
[ T_2 W ] = h_2 scenario 2
[ ... ... ]
[ T_S W ] = h_S scenario S
The size grows linearly with $S$. For large $S$, decomposition exploits the structure. The L-shaped method is Benders decomposition (Section 2.6) applied to this matrix: a master problem in $x$ plus one independent LP per scenario, which returns optimality cuts. Progressive hedging instead decomposes by scenario and gradually enforces agreement on $x$.
Where scenarios come from.
| Source | When to use | Caveat |
|---|---|---|
| expert scenarios (good / average / bad) | little data, strategic decisions | probabilities are subjective |
| historical data | stable environment | the past may not repeat |
| discretizing a distribution | few random parameters | grows exponentially with the number of parameters |
| Monte Carlo sampling (SAA) | many or continuous random parameters | needs enough samples. Check stability |
| scenario reduction, moment matching | large scenario sets | keep the extremes that drive decisions |
Sample average approximation (SAA). Draw $N$ i.i.d. scenarios, give each probability $1/N$, and solve the extensive form. The solution and its value $\hat z_N$ converge to the true ones as $N \to \infty$. Two facts guide its use:
- The in-sample value $\hat z_N$ is optimistically biased: the optimizer adapts the plan to the particular sample. For a maximization problem, $E[\hat z_N] \ge z^*$.
- The honest measure of a plan is out-of-sample evaluation: fix the plan and compute its average result on a large, fresh sample. If the plan and the out-of-sample value stop changing as $N$ grows, the sample is large enough.
Multistage problems. When decisions and observations alternate over several periods, scenarios form a tree. Nonanticipativity requires that scenarios sharing a history up to time $t$ share the decision at time $t$:
t = 1 t = 2 t = 3
+-- high --> x2(H) ---+-- high --> x3(HH)
x1 -----------+ +-- low --> x3(HL)
(one decision) +-- low --> x2(L) ---+-- high --> x3(LH)
+-- low --> x3(LL)
x2(H) is shared by scenarios HH and HL: it cannot know the t = 3 outcome
The tree has $b^{T-1}$ leaves for branching factor $b$ and $T$ stages, so multistage models need aggressive scenario reduction or approximate methods such as approximate dynamic programming (Section 6.2).
Example: the farmer's extensive form, then SAA with continuous yields. First, the three-scenario extensive form. Then the yields of the three crops vary independently and uniformly within ±20% of the mean, and SAA is solved for growing sample sizes. Each plan is evaluated on 100,000 fresh scenarios. This second stage has a closed form (sell any surplus, buy any deficit), which makes the out-of-sample check fast.
import numpy as np
from scipy.optimize import linprog
# plant, q, need, mean_yield, scenarios, probs: data from the previous example
def extensive_form(T, p):
"""Deterministic equivalent of the farmer problem for scenario yields T, probabilities p."""
S = len(p)
n = 3 + 6 * S # columns: x, then y_1, ..., y_S
c = np.concatenate([plant] + [p_s * q for p_s in p])
A, b = [], []
def add(entries, rhs): # one "<=" row from {column: coefficient}
r = np.zeros(n)
for k, v in entries.items():
r[k] = v
A.append(r); b.append(rhs)
add({0: 1, 1: 1, 2: 1}, 500) # land (first stage)
for s, t in enumerate(T):
o = 3 + 6 * s # first column of scenario s
add({0: -t[0], o: -1, o + 1: 1}, -need[0]) # wheat: grown + bought - sold >= 200
add({1: -t[1], o + 2: -1, o + 3: 1}, -need[1]) # corn : grown + bought - sold >= 240
add({2: -t[2], o + 4: 1, o + 5: 1}, 0) # beets sold <= beets grown
add({o + 4: 1}, 6000) # beet quota at the high price
res = linprog(c, A_ub=np.array(A), b_ub=np.array(b), bounds=(0, None), method="highs")
return res.x[:3], -res.fun # planting plan, expected profit
x, z = extensive_form(scenarios, probs)
print(f"3 scenarios: plant {np.round(x, 1)}, expected profit {z:,.0f}")
# Sample average approximation: yields vary continuously, +-20% independently per crop
def profit(x, T): # closed-form second stage (vectorized)
net = T[:, :2] * x[:2] - need # surplus (+) or deficit (-) of wheat and corn
sell, buy = np.array([170, 150]), np.array([238, 210])
grain = np.where(net >= 0, net * sell, net * buy).sum(axis=1)
beets = T[:, 2] * x[2]
return grain + 36 * np.minimum(beets, 6000) + 10 * np.maximum(beets - 6000, 0) - plant @ x
rng = np.random.default_rng(seed=21)
test = mean_yield * rng.uniform(0.8, 1.2, size=(100_000, 3)) # out-of-sample scenarios
for N in (3, 10, 100, 1000):
T = mean_yield * rng.uniform(0.8, 1.2, size=(N, 3))
x, z = extensive_form(T, np.full(N, 1 / N))
print(f"SAA N={N:4}: plant {np.round(x, 1)}, in-sample {z:9,.0f},"
f" out-of-sample {profit(x, test).mean():9,.0f}")
# 3 scenarios: plant [170. 80. 250.], expected profit 108,390
# SAA N= 3: plant [164.1 81.1 254.8], in-sample 98,612, out-of-sample 109,494
# SAA N= 10: plant [112.1 80.6 307.4], in-sample 110,910, out-of-sample 109,226
# SAA N= 100: plant [144.4 83.6 272. ], in-sample 114,401, out-of-sample 111,118
# SAA N=1000: plant [135.5 85.2 279.3], in-sample 111,612, out-of-sample 111,254
With three scenarios, the extensive form returns plan B with expected profit 108,390, which confirms the previous topic. With continuous yields, small samples give unstable plans: wheat ranges from 112 to 164 acres for $N = 3$ and $N = 10$. Their in-sample values are unreliable, 98,612 (an unlucky sample) and 110,910 (optimistic). At $N = 1000$ the plan is about (135, 85, 279) acres, and the in-sample and out-of-sample values agree within 400. The $N = 100$ run shows the optimistic bias clearly: it promises 114,401 but delivers 111,118. The plan also differs from plan B, because the continuous model has milder deviations (uniform rather than always ±20%) and independent crops that partly offset each other, whereas the three-scenario model moves all yields together. How scenarios are generated changes the decision.
Key takeaways
- With finitely many scenarios, a two-stage program becomes one LP (the extensive form), with a copy of the second-stage variables per scenario.
- Its block structure is exploited by the L-shaped method (Benders) and by progressive hedging.
- SAA replaces the distribution with a random sample. In-sample values are optimistic, so judge plans out of sample and increase $N$ until results stabilize.
- Multistage problems need scenario trees with nonanticipativity, and their size grows exponentially with the number of stages.
Practice
- (Beginner) Count the variables and constraints of the farmer's extensive form for $S = 3$ and for $S = 1000$. Which part grows with $S$?
- (Intermediate) Run SAA with $N = 100$ for 20 different seeds. Report the mean and standard deviation of the in-sample objective and of the out-of-sample profit, and check the optimistic bias. How do the numbers change for $N = 500$?
- (Interview) "How would you choose the number of scenarios for a stochastic program used in production?" Hint: look at the stability of the solution and objective across independent samples, the gap between in-sample and out-of-sample values, the computing budget, scenario reduction, and decomposition.
Chance-Constrained Programming
Why it matters. Sometimes a violation cannot sensibly be priced as a recourse cost: a blackout, a stock-out of a critical drug, a structural load beyond its limit. Instead, the decision maker requires a constraint to hold with high probability, for example 95% of the time. Chance constraints express such reliability or service-level requirements directly.
Theory. An individual chance constraint requires one random constraint to hold with probability at least $1 - \alpha$:
$$ P\big(a(\xi)^\top x \le b(\xi)\big) \ge 1 - \alpha. $$
A joint chance constraint requires several constraints to hold simultaneously with probability at least $1 - \alpha$, which is stricter.
Deterministic equivalents. Some cases reduce to ordinary constraints:
- Random right-hand side only. $P(x \ge D) \ge 1 - \alpha$ is equivalent to $x \ge F_D^{-1}(1 - \alpha)$, the $(1-\alpha)$-quantile of $D$. For normal demand, $x \ge \mu + z_{1-\alpha}\sigma$: the safety stock formula of inventory theory (Chapter 10).
- Normal random coefficients. If $a \sim N(\bar a, \Sigma)$ and $b$ is fixed, $P(a^\top x \le b) \ge 1 - \alpha$ is equivalent to $$ \bar a^\top x + z_{1-\alpha},\big|\Sigma^{1/2} x\big|_2 \le b, $$ a second-order cone constraint, which is convex for $\alpha < 0.5$.
In general, the feasible set of a chance constraint is not convex, and joint constraints are harder still. Practical approaches:
| Approach | Idea | Property |
|---|---|---|
| Bonferroni | split the risk: each of $m$ constraints with $\alpha/m$ | conservative but simple, since $P(\text{any fails}) \le \sum \alpha_i$ |
| sample approximation with binaries | sample $N$ scenarios; binary $z_k = 1$ allows scenario $k$ to be violated; $\sum_k z_k \le \alpha N$ | exact for the sample. A MILP, so it needs tight big-M values |
| scenario approach (Calafiore-Campi) | enforce the constraint in all of $N$ sampled scenarios | convex. $N$ chosen to guarantee $1-\alpha$ with high confidence |
| CVaR approximation | replace the probability with conditional value-at-risk | convex, conservative, accounts for the size of violations |
Chance constraints have a known weakness: they ignore how badly a constraint is violated. A 5% chance of being one unit short counts the same as a 5% chance of being a thousand units short. Integrated chance constraints and CVaR-based constraints address this.
Example: stocking two critical products. A hospital pharmacy stocks two drugs with independent normal demands, $D_1 \sim N(100, 20^2)$ and $D_2 \sim N(150, 30^2)$. Units cost 3 and 2. It wants a probability of at least 90% that both demands are fully met. The code compares four ways to set the stock and checks each one on 200,000 fresh demand scenarios.
import numpy as np
from scipy.stats import norm
from scipy.optimize import milp, LinearConstraint, Bounds
cost = np.array([3.0, 2.0]) # cost per unit stocked
mu, sd = np.array([100.0, 150.0]), np.array([20.0, 30.0]) # independent normal demands
alpha = 0.10 # allowed risk: P(any stock-out) <= 10%
rng = np.random.default_rng(seed=22)
test = rng.normal(mu, sd, size=(200_000, 2)) # fresh demands to check each plan
def report(name, x):
ok = np.mean((test <= x).all(axis=1)) # joint: both products fully served
print(f"{name:>26}: x = {np.round(x, 1)}, cost = {cost @ x:6.1f}, P(no stock-out) = {ok:.3f}")
# (a) Individual chance constraints P(x_i >= D_i) >= 0.9: deterministic equivalent x_i = mu + z sd
report("individual 90% each", mu + norm.ppf(1 - alpha) * sd)
# (b) Bonferroni: 95% each guarantees >= 90% jointly (conservative)
report("Bonferroni 95% each", mu + norm.ppf(1 - alpha / 2) * sd)
# (c) Exact joint constraint (independence): Phi(z1) Phi(z2) >= 0.9, search over the split
best = min((cost @ (mu + norm.ppf([p1, (1 - alpha) / p1]) * sd), p1)
for p1 in np.linspace(0.901, 0.999, 9_801))
p1 = best[1]
report("exact joint (normal)", mu + norm.ppf([p1, (1 - alpha) / p1]) * sd)
# (d) Sample approximation with binaries: z_k = 1 lets scenario k be violated
N = 500
d = rng.normal(mu, sd, size=(N, 2))
n_viol = int(alpha * N) # at most n_viol scenarios may be violated
low = np.sort(d, axis=0)[-(n_viol + 1)] # so x_i >= (n_viol+1)-th largest demand
M = np.maximum(d - low, 0) # scenario-specific big-M (tight)
# variables: x1, x2, z_1..z_N ; constraints x_i + M_ki z_k >= d_ki and sum z <= n_viol
A = np.zeros((2 * N + 1, 2 + N))
for i in range(2):
A[i * N:(i + 1) * N, i] = 1
A[i * N:(i + 1) * N, 2:] = np.diag(M[:, i])
A[-1, 2:] = 1
lb = np.r_[d[:, 0], d[:, 1], -np.inf]
ub = np.r_[np.full(2 * N, np.inf), n_viol]
res = milp(np.r_[cost, np.zeros(N)], constraints=LinearConstraint(A, lb, ub),
integrality=np.r_[0, 0, np.ones(N)],
bounds=Bounds(np.r_[low, np.zeros(N)], np.r_[np.inf, np.inf, np.ones(N)]))
report(f"sampled joint, N={N}", res.x[:2])
# individual 90% each: x = [125.6 188.4], cost = 753.8, P(no stock-out) = 0.809
# Bonferroni 95% each: x = [132.9 199.3], cost = 797.4, P(no stock-out) = 0.902
# exact joint (normal): x = [132.6 199. ], cost = 795.9, P(no stock-out) = 0.899
# sampled joint, N=500: x = [130.1 206.6], cost = 803.5, P(no stock-out) = 0.905
Imposing 90% on each product separately is cheapest, but both demands are met only 81% of the time ($0.9^2$), so the joint requirement is violated. Bonferroni (95% each) satisfies it at a cost of 797.4. The exact joint optimum, which uses the normal distribution and independence, costs 795.9. Bonferroni is nearly optimal here, as it usually is for small $\alpha$. The sampled MILP needs no distributional formula and handles any dependence, but it optimizes for its 500 scenarios: its plan costs more and differs in shape, an effect of sampling noise. The tight, scenario-specific big-M values make it solve in seconds. With a single large M instead, the solver struggles: a run with $N = 1000$ did not finish in ten minutes.
Key takeaways
- Chance constraints require constraints to hold with probability at least $1 - \alpha$, individually or jointly.
- With randomness only in the right-hand side, they reduce to quantile constraints, such as safety stock $\mu + z_{1-\alpha}\sigma$. Normal coefficients give second-order cone constraints.
- Joint constraints are stricter than individual ones. Bonferroni is a simple conservative approximation, and sampled MILPs are exact for the sample.
- Chance constraints ignore the size of violations. Consider CVaR or integrated chance constraints when that matters.
Practice
- (Beginner) A blood bank wants enough O-negative units to cover next week's demand, normal with mean 80 and standard deviation 15, with probability 99%. How many units should it stock, and how many of them are safety stock?
- (Intermediate) Repeat the example with correlated demands (correlation +0.8, then -0.8). Simulate the joint service level of the "individual 90%" and Bonferroni plans in both cases. Under which correlation is Bonferroni more conservative, and why?
- (Interview) "When would you use a chance constraint rather than a penalty cost in a two-stage model, and what is the main weakness of chance constraints?" Hint: consider unquantifiable or regulatory consequences, nonconvexity, and the fact that only the probability of a violation counts, not its size.
Value of the Stochastic Solution
Why it matters. Stochastic models cost more to build, explain and solve than deterministic ones. Before investing, an analyst should know how much better the stochastic solution is than "plug in the mean", and how much a perfect forecast would be worth. The value of the stochastic solution (VSS) and the expected value of perfect information (EVPI) answer these two questions.
Theory. For a maximization problem with profit $f(x, \xi)$, define four quantities:
| Name | Definition | Meaning |
|---|---|---|
| RP (recourse problem) | $\max_x E_\xi[f(x, \xi)]$ | value of the stochastic program, decided here-and-now |
| EV (expected value problem) | $\max_x f(x, \bar\xi)$, with solution $\bar x$ and $\bar\xi = E[\xi]$ | what the mean-value model promises |
| EEV (expected result of EV solution) | $E_\xi[f(\bar x, \xi)]$ | what the mean-value plan actually earns on average |
| WS (wait-and-see) | $E_\xi[\max_x f(x, \xi)]$ | value with a perfect forecast before deciding |
From these:
$$ \mathrm{VSS} = \mathrm{RP} - \mathrm{EEV} \ge 0, \qquad \mathrm{EVPI} = \mathrm{WS} - \mathrm{RP} \ge 0, \qquad \mathrm{EEV} \le \mathrm{RP} \le \mathrm{WS}. $$
- VSS is the gain from modeling uncertainty instead of averaging it away. It is nonnegative because RP optimizes the true expected profit, while $\bar x$ is just one feasible plan.
- EVPI is the most one should pay for a perfect forecast. It is the same quantity as in Section 7.2, now for a decision with recourse. It is nonnegative because the clairvoyant can always copy the here-and-now plan.
- EV is not a meaningful estimate of anything. In many models (for example, minimization with randomness only in the right-hand side, by Jensen's inequality) the mean-value model is systematically over-optimistic about its own promise.
How to read the two values:
| VSS | EVPI | Interpretation |
|---|---|---|
| small | small | uncertainty hardly matters. A deterministic model is fine |
| large | any | averaging is costly: use the stochastic model |
| small | large | the stochastic model cannot hedge much, but better forecasts are valuable |
A small VSS for the expected profit does not mean the plans are equally risky. Compare their profit distributions, in particular the worst scenarios.
Example: what stochastic modeling is worth to the farmer.
import numpy as np
# scenarios, probs, mean_yield, extensive_form, profit: from the two previous examples
# RP: recourse problem (the stochastic program itself)
x_rp, RP = extensive_form(scenarios, probs)
# EV: replace yields by their mean, then EEV: evaluate that plan against the real scenarios
x_ev, EV = extensive_form(mean_yield[None, :], np.array([1.0]))
EEV = probs @ profit(x_ev, scenarios)
# WS: wait-and-see, a clairvoyant plans for each scenario separately
WS = sum(p * extensive_form(t[None, :], np.array([1.0]))[1] for t, p in zip(scenarios, probs))
print(f"EV = {EV:9,.0f} plan {np.round(x_ev)} (promise of the mean-value model)")
print(f"EEV = {EEV:9,.0f} (what the mean-value plan actually earns on average)")
print(f"RP = {RP:9,.0f} plan {np.round(x_rp)}")
print(f"WS = {WS:9,.0f}")
print(f"VSS = RP - EEV = {RP - EEV:7,.0f}")
print(f"EVPI = WS - RP = {WS - RP:7,.0f}")
# EV = 118,600 plan [120. 80. 300.] (promise of the mean-value model)
# EEV = 107,240 (what the mean-value plan actually earns on average)
# RP = 108,390 plan [170. 80. 250.]
# WS = 115,406
# VSS = RP - EEV = 1,150
# EVPI = WS - RP = 7,016
The mean-value model promises 118,600, but its plan earns only 107,240 on average: the promise overstates reality by more than 10%. The stochastic plan earns 1,150 more per year (VSS) with no new data, only a better model. A perfect weather forecast would add another 7,016 (EVPI), which caps what the farmer should pay for forecasting services. Note that the mean-value model's promise of 118,600 exceeds even what a clairvoyant could earn on average (115,406). A VSS of about 1% looks small, and expected profit is not the whole story. In the bad year the stochastic plan earns 48,820, less than plan A's 55,120, while in the good year it earns 167,000 against 148,000. A risk-averse farmer might prefer yet another plan, which a risk measure in the objective would capture.
Key takeaways
- RP: optimal expected value of the stochastic program. EEV: expected value of the mean-value plan. WS: expected value with perfect foresight.
- VSS = RP - EEV measures the benefit of modeling uncertainty. EVPI = WS - RP bounds the value of better forecasts.
- Always $\mathrm{EEV} \le \mathrm{RP} \le \mathrm{WS}$. The deterministic model's own objective (EV) is usually over-optimistic.
- Judge plans also on risk (profit distribution, worst case), not only on expected value.
Practice
- (Beginner) For the newsvendor of Section 1.1 (price 10, cost 2, demand uniform on [50, 150]), compute EV, EEV, RP and WS, then VSS and EVPI. Hint: WS $= E[8D]$.
- (Intermediate) Raise the farmer's beet quota from 6,000 to 8,000 tons and recompute VSS and EVPI. Explain the change in terms of the asymmetry between the two beet prices.
- (Interview) "A colleague reports VSS = 0.3% of the objective and proposes dropping the stochastic model. Do you agree?" Hint: compare the profit distributions, especially the worst scenarios, of the two plans. Also consider EVPI, the cost of maintaining the model, and whether the scenarios represent the uncertainty well.
Robust Optimization
Why it matters. Stochastic programming needs probability distributions, which are often unknown or unreliable. In some applications, a solution must also stay feasible for every plausible realization of the data: safety margins, contractual commitments, engineering tolerances. Robust optimization replaces the distribution with an uncertainty set and requires feasibility for every data value in it. The result is still a tractable deterministic problem.
Theory. Consider an LP whose constraint coefficients are uncertain, with row $i$ ranging over a set $\mathcal{U}_i$:
$$ \max; c^\top x \quad \text{s.t.} \quad a_i^\top x \le b_i ;;; \text{for all } a_i \in \mathcal{U}_i, ; i = 1, \dots, m. $$
The robust counterpart writes each robust constraint as $\max_{a_i \in \mathcal{U}_i} a_i^\top x \le b_i$. The inner maximization is the "adversary" choosing the worst data, the zero-sum view of nature from Chapter 8. For common uncertainty sets it has a closed form or an LP dual, so the robust counterpart is a single tractable problem:
| Uncertainty set | Description | Robust counterpart of $a_i^\top x \le b_i$ | Problem class |
|---|---|---|---|
| box (Soyster) | each $a_{ij} \in [\bar a_{ij} - \hat a_{ij}, \bar a_{ij} + \hat a_{ij}]$ | $\bar a_i^\top x + \sum_j \hat a_{ij} \lvert x_j \rvert \le b_i$ | LP |
| ellipsoid (Ben-Tal-Nemirovski) | $a_i = \bar a_i + P_i u$, $\lVert u \rVert_2 \le \Omega$ | $\bar a_i^\top x + \Omega,\lVert P_i^\top x \rVert_2 \le b_i$ | SOCP |
| budget (Bertsimas-Sim) | at most $\Gamma_i$ coefficients of row $i$ deviate | see below | LP |
The box set assumes every coefficient takes its worst value simultaneously, which is very unlikely and very conservative. The budget of uncertainty set of Bertsimas and Sim limits the adversary to $\Gamma_i$ deviating coefficients per row. Here $\Gamma_i = 0$ is the nominal problem and $\Gamma_i = n$ the box. The worst-case term $\max_{|S| \le \Gamma_i} \sum_{j \in S} \hat a_{ij}|x_j|$ is the optimal value of an LP in $S$. Replacing it with its dual gives linear constraints (for $x \ge 0$):
$$ \bar a_i^\top x + \Gamma_i, z_i + \sum_j p_{ij} \le b_i, \qquad z_i + p_{ij} \ge \hat a_{ij}, x_j, \qquad z_i, p_{ij} \ge 0. $$
If the deviations are independent and symmetric, a constraint protected with budget $\Gamma_i$ is violated with probability at most about $\exp(-\Gamma_i^2 / (2n))$. A modest budget therefore buys high reliability when rows have many uncertain coefficients.
The profit given up for protection is the price of robustness. Choosing the uncertainty set means choosing a risk attitude.
Related approaches lie between stochastic and robust optimization:
- Adjustable robust optimization adds recourse: second-stage decisions are functions of the revealed data, often restricted to affine decision rules.
- Distributionally robust optimization optimizes the expected value against the worst distribution in an ambiguity set, for example all distributions with given moments or within a Wasserstein distance of the data.
| Stochastic programming | Chance constraints | Robust optimization | |
|---|---|---|---|
| Uncertainty described by | probability distribution or scenarios | distribution | uncertainty set |
| Objective | expected value (plus risk measures) | expected or nominal value | worst case over the set |
| Feasibility | via recourse in every scenario | with probability $\ge 1 - \alpha$ | for every value in the set |
| Typical difficulty | many scenarios | nonconvexity | choosing the set (conservatism) |
Example: a production plan robust to uncertain processing times. Four products use machine and labor hours. Each processing time may deviate by up to ±20% from its nominal value, and each product's sales are capped at 60 units. The code solves the budgeted robust counterpart for several values of $\Gamma$ (the same in both rows). It then estimates each plan's probability of being infeasible when the coefficients are drawn uniformly from their ranges.
import numpy as np
from scipy.optimize import linprog
profit = np.array([12.0, 10.0, 9.0, 14.0]) # per unit of products 1..4
a_bar = np.array([[2.0, 3.0, 1.5, 4.0], # machine hours per unit (nominal)
[3.0, 2.0, 2.5, 3.0]]) # labour hours per unit (nominal)
a_hat = 0.2 * a_bar # each coefficient may deviate by +-20%
b = np.array([400.0, 450.0])
cap = 60.0 # market limit per product
m, n = a_bar.shape
def robust_plan(Gamma):
"""Bertsimas-Sim budgeted robust counterpart: at most Gamma coefficients per row deviate.
Variables: x (n), z (m), p (m*n). Row i: a_bar_i x + Gamma z_i + sum_j p_ij <= b_i,
z_i + p_ij >= a_hat_ij x_j, all variables >= 0."""
nv = n + m + m * n
P = lambda i, j: n + m + i * n + j # column of p_ij
A, rhs = [], []
for i in range(m):
r = np.zeros(nv); r[:n] = a_bar[i]; r[n + i] = Gamma
for j in range(n):
r[P(i, j)] = 1
A.append(r); rhs.append(b[i])
for j in range(n): # a_hat_ij x_j - z_i - p_ij <= 0
r = np.zeros(nv); r[j] = a_hat[i, j]; r[n + i] = -1; r[P(i, j)] = -1
A.append(r); rhs.append(0.0)
bounds = [(0, cap)] * n + [(0, None)] * (m + m * n)
res = linprog(-np.r_[profit, np.zeros(m + m * n)], A_ub=np.array(A), b_ub=rhs,
bounds=bounds, method="highs")
return res.x[:n], -res.fun
rng = np.random.default_rng(seed=24)
draws = a_bar + a_hat * rng.uniform(-1, 1, size=(100_000, m, n)) # realized coefficients
print("Gamma plan profit P(plan infeasible)")
for Gamma in (0, 0.5, 1, 2, 3, 4):
x, z = robust_plan(Gamma)
viol = np.mean(((draws @ x) > b + 1e-9).any(axis=1))
print(f"{Gamma:5} {np.round(x, 1) + 0.0} {z:7.1f} {viol:.4f}")
# Gamma plan profit P(plan infeasible)
# 0 [60. 60. 54.5 4.5] 1874.5 0.7521
# 0.5 [51.6 34.4 53.4 25.8] 1805.0 0.5253
# 1 [46.3 33.9 55.6 25.4] 1752.0 0.3046
# 2 [60. 31.1 34. 23.3] 1662.7 0.0320
# 3 [60. 30.2 33.7 20.1] 1607.3 0.0011
# 4 [60. 53.5 35.2 0. ] 1571.9 0.0000
The nominal plan ($\Gamma = 0$) earns the most on paper, 1,874.5, but it uses capacity so tightly that it becomes infeasible 75% of the time. Its stated profit is an illusion. Protecting against 2 deviating coefficients per row costs 11% of the nominal profit and cuts the infeasibility risk to 3%. With $\Gamma = 3$ it falls to 0.1%. Full box protection ($\Gamma = 4$) gains little extra safety for another 2% of profit. The product mix also changes with $\Gamma$, and not monotonically: product 4 goes from 4.5 units to over 20 and back to zero. Robustness is not simply "produce a little less of everything". The curve of profit against risk lets the decision maker choose the protection level knowingly.
Key takeaways
- Robust optimization requires feasibility for every data value in an uncertainty set. No probability distribution is needed.
- Robust counterparts are tractable: box and budget sets give LPs, ellipsoids give SOCPs.
- The Bertsimas-Sim budget $\Gamma$ tunes conservatism between the nominal problem ($\Gamma = 0$) and the box ($\Gamma = n$), trading profit (the price of robustness) for reliability.
- Adjustable and distributionally robust models bridge robust and stochastic optimization.
Practice
- (Beginner) Write the box (Soyster) robust counterpart of $2x_1 + 3x_2 \le 12$, $x \ge 0$, when each coefficient may deviate by ±10%. By what factor does the feasible region shrink along each axis?
- (Intermediate) For the example's rows with $n = 4$ uncertain coefficients, find the budget $\Gamma$ for which the bound $\exp(-\Gamma^2/(2n))$ guarantees at most 5% violation per row. Compare it with $n$ and with the simulated violation probabilities. What does this say about using the bound for short rows?
- (Interview) "You are designing a supply-chain network (which warehouses to open) under uncertain future demand. Would you use stochastic programming, chance constraints or robust optimization?" Hint: consider what data exist (history, scenarios, only ranges), whether there is natural recourse (re-routing flows after demand is known), the risk attitude of management, and model size. Hybrids such as two-stage robust and distributionally robust models are worth mentioning.
<a id="10-applied-or-models"></a>