a linear program assumes you know the data. stochastic programming admits that you do not — some coefficients are random — but insists you know their distribution, and asks for the decision that is best on average. 𐃏 the structural insight that makes this a paradigm rather than a hack: split the decision in two. commit to \(x\) now, before the coin is flipped; after uncertainty resolves, take a corrective recourse action \(y\) that adapts to whatever happened. the objective charges you for both, weighting the second stage by expectation.¹

the “distribution is known” clause is a genuine modelling commitment — those probabilities have to come from somewhere (data, physics, or a prior you are prepared to defend (Jaynes, Edwin T., 2003)). when you can only bound the uncertainty set and refuse to weight it, you are doing robust optimisation instead: worst case rather than average case. the two paradigms bracket most practical uncertainty modelling.

two-stage programs with recourse#

formulation#

the canonical two-stage stochastic linear program:

\begin{align*} \text{minimise} \quad & c^{\top} x + \mathbb{E}_{\xi}\!\left[ Q(x, \xi) \right] \\ \text{subject to} \quad & A x = b, \quad x \ge 0, \end{align*}

where the second-stage value function \(Q(x, \xi)\) is itself the optimal value of a linear program, solved after observing the random element \(\xi = (q, T, W, h)\):

\begin{align*} Q(x, \xi) = \min_{y} \quad & q^{\top} y \\ \text{subject to} \quad & W y = h - T x, \quad y \ge 0. \end{align*}

vocabulary worth pinning down:

• here-and-now decision \(x\): chosen under ignorance, shared by every future.
• recourse decision \(y\): chosen with hindsight, one per realisation of \(\xi\). \(T\) (the technology matrix) says how first-stage choices constrain the second stage; \(W\) (the recourse matrix) says what corrections are available.
• fixed recourse: \(W\) deterministic. complete recourse: \(\{Wy : y \ge 0\} = \mathbb{R}^m\), so the second stage is feasible whatever you did — every mistake is patchable at a price. relatively complete recourse asks this only for first-stage-feasible \(x\).
• \(\mathcal{Q}(x) = \mathbb{E}_{\xi}[Q(x,\xi)]\), the expected recourse cost, is convex in \(x\) (a minimum of an LP is convex in its right-hand side, and expectation preserves convexity); with finitely many scenarios and fixed recourse it is piecewise linear. so the whole problem is a convex program — large, but convex.

scenario trees#

continuous \(\xi\) makes \(\mathcal{Q}(x)\) an integral of LP values — usually hopeless to evaluate exactly. the standard move is to discretise: approximate the distribution by finitely many scenarios \(\xi_1, \dots, \xi_S\) with probabilities \(p_1, \dots, p_S\). the resulting structure is a tree: one root where \(x\) lives, one branch per scenario, one recourse decision per leaf. 𐃏

the deterministic equivalent#

with finitely many scenarios the expectation is a weighted sum, and the whole thing collapses into one ordinary — if enormous — LP, the deterministic equivalent:

\begin{align*} \text{minimise} \quad & c^{\top} x + \sum_{s=1}^{S} p_s \, q_s^{\top} y_s \\ \text{subject to} \quad & A x = b \\ & T_s x + W y_s = h_s \qquad s = 1, \dots, S \\ & x \ge 0, \quad y_s \ge 0 \qquad s = 1, \dots, S. \end{align*}

its constraint matrix has the signature dual block-angular shape — a spine of first-stage columns coupling \(S\) otherwise-independent diagonal blocks:

\begin{equation} \begin{bmatrix} A & & & & \\ T_1 & W & & & \\ T_2 & & W & & \\ \vdots & & & \ddots & \\ T_S & & & & W \end{bmatrix} \begin{bmatrix} x \\ y_1 \\ y_2 \\ \vdots \\ y_S \end{bmatrix}#

\begin{bmatrix} b \\ h_1 \\ h_2 \\ \vdots \\ h_S \end{bmatrix} \end{equation}

the structure is the algorithm. fix \(x\) and the problem shatters into \(S\) small independent second-stage LPs — one per scenario, embarrassingly parallel. the L-shaped method (benders decomposition specialised to this matrix, named after the shape of the nonzero pattern) exploits exactly this: a master problem proposes \(x\), the scenario subproblems return dual prices, and those duals become optimality cuts — supporting hyperplanes of the convex function \(\mathcal{Q}(x)\) — plus feasibility cuts when a proposal leaves some scenario unpatchable. the master re-solves against its growing outer approximation of \(\mathcal{Q}\) until the bounds close.² with 10 first-stage variables and a million scenarios, the deterministic equivalent is a million-row LP, but each subproblem stays tiny.

sample average approximation#

when scenarios are not given but the distribution can be sampled, replace the expectation with an empirical average over \(N\) iid draws \(\xi^{1}, \dots, \xi^{N}\):

\begin{equation} \min_{x \in X} \; \hat{f}_N(x) = c^{\top} x + \frac{1}{N} \sum_{i=1}^{N} Q(x, \xi^{i}), \end{equation}

and solve that — it is just a deterministic equivalent with equal weights \(1/N\). this is sample average approximation (SAA): monte carlo on the objective, then honest deterministic optimisation.

why it works: for each fixed \(x\) the law of large numbers gives \(\hat f_N(x) \to f(x)\) almost surely, and under mild conditions (compact \(X\), integrable and well-behaved \(Q\)) the convergence is uniform in \(x\) — the same machinery that pushes the empirical cdf uniformly onto the true cdf (Wasserman, Larry, 2010). uniform convergence of objectives drags the minimisers along: the SAA optimal value \(\hat v_N\) converges to the true \(v^{*}\) and SAA solutions accumulate on the true solution set.³

the honest caveats:

• the in-sample value flatters you. \(\mathbb{E}[\hat v_N] \le v^{*}\): minimising over the sample overfits the sample, so the reported optimum is biased optimistic (the bias shrinks monotonically in \(N\)). always re-evaluate the candidate \(\hat x_N\) on a fresh, much larger sample to estimate its true cost.
• rare events break it. if the expensive scenarios have probability \(10^{-4}\), a thousand samples will likely contain none of them and SAA will confidently recommend a decision that has never seen a bad day.
• no error bars for free. one SAA run is a point estimate; confidence intervals need replications (solve \(M\) independent SAA instances) or the optimality-gap estimators built from the bias inequality above.
• sample size is problem-dependent. accuracy in the solution (not just the value) can demand far larger \(N\) than the central-limit heuristic suggests, especially with flat objectives near the optimum.

chance constraints#

sometimes recourse is the wrong metaphor — some constraints must simply hold with high probability, no cleanup available. a chance constraint (charnes and cooper, 1959):

\begin{equation} \mathbb{P}\left( a(\xi)^{\top} x \le b \right) \ge 1 - \alpha, \end{equation}

for a small tolerance \(\alpha\): a bridge that carries the load on 99.9 percent of days, a portfolio that breaches its risk limit at most one day in twenty.

when they convexify#

the gaussian case is the classic success. suppose \(a \sim \mathcal{N}(\bar a, \Sigma)\). then \(a^{\top} x \sim \mathcal{N}(\bar a^{\top} x,\; x^{\top} \Sigma x)\), so

\begin{equation} \mathbb{P}\left( a^{\top} x \le b \right) = \Phi\!\left( \frac{b - \bar a^{\top} x}{\sqrt{x^{\top} \Sigma x}} \right) \ge 1 - \alpha \quad\Longleftrightarrow\quad \bar a^{\top} x + \Phi^{-1}(1-\alpha)\, \lVert \Sigma^{1/2} x \rVert_2 \le b, \end{equation}

where \(\Phi\) is the standard normal cdf. for \(\alpha \le 1/2\) the quantile \(\Phi^{-1}(1-\alpha) \ge 0\) and this is a second-order cone constraint — convex, and solvable at scale by the interior-point machinery of conic programming. the probabilistic statement has been traded for a deterministic one: back off from the mean constraint by a safety margin of \(\Phi^{-1}(1-\alpha)\) standard deviations. squint and it is robust optimisation with an ellipsoidal uncertainty set — the two paradigms meeting in the middle.

when they do not#

• the wrong tolerance. \(\alpha > 1/2\) flips the sign of the quantile and the same set turns nonconvex (a reverse cone). rarely wanted, but it shows convexity is a property of the confidence level, not just the distribution.
• general distributions. for arbitrary \(\xi\) the feasible set \(\{x : \mathbb{P}(g(x,\xi) \le 0) \ge 1-\alpha\}\) can be nonconvex and even disconnected — probability mass can sit in awkward lumps, and unioning “mostly-feasible” regions does not preserve convexity.
• finite scenarios go combinatorial. with \(S\) scenarios, “hold in at least \((1-\alpha)S\) of them” needs a binary indicator per scenario and big-M constraints: \(a_s^{\top} x \le b + M z_s\), \(\sum_s z_s \le \alpha S\). that is an integer program — you are literally choosing which scenarios to sacrifice.
• joint chance constraints. asking several inequalities to hold simultaneously with probability \(1-\alpha\) is harder still: even in the gaussian case the multivariate normal cdf has no usable closed form. the workhorse fixes are conservative convex surrogates — replace the probability by a cvar (expected-shortfall) constraint, which is convex and implies the chance constraint, at the price of some feasibility headroom.³

the newsvendor#

the fruit-fly of stochastic programming: every result in the field can be demonstrated on it, and it is the only interesting stochastic program with a closed-form solution. 𐃏

a vendor orders \(q\) papers before knowing demand \(D \sim F\). each unit of unmet demand costs \(c_u\) (foregone margin); each unsold paper costs \(c_o\) (writedown). this is a two-stage program: first stage “order \(q\)”, second stage “observe \(D\), absorb the mismatch” — with the recourse LP so trivial it evaluates in closed form:

\begin{equation} C(q) = \mathbb{E}\left[ c_o (q - D)^{+} + c_u (D - q)^{+} \right] = c_o \int_{0}^{q} (q - t)\, f(t)\, dt + c_u \int_{q}^{\infty} (t - q)\, f(t)\, dt. \end{equation}

deriving the critical fractile#

differentiate under the integral sign (leibniz rule — the boundary terms vanish because the integrands are zero at \(t = q\)):

\begin{align*} C’(q) &= c_o \left[ (q-t) f(t) \Big|_{t=q} + \int_{0}^{q} f(t)\, dt \right] + c_u \left[ -(t-q) f(t) \Big|_{t=q} - \int_{q}^{\infty} f(t)\, dt \right] \\ &= c_o F(q) - c_u \big( 1 - F(q) \big). \end{align*}

so the marginal paper is worth ordering exactly while \(c_u \mathbb{P}(D > q)\), the expected gain from covering one more unit of demand, exceeds \(c_o \mathbb{P}(D \le q)\), the expected loss from one more leftover. moreover \(C’’(q) = (c_o + c_u) f(q) \ge 0\), so \(C\) is convex and the first-order condition is the global minimum:

\begin{equation} \boxed{\; F(q^{*}) = \frac{c_u}{c_u + c_o} \quad\Longleftrightarrow\quad q^{*} = F^{-1}\!\left( \frac{c_u}{c_u + c_o} \right) \;} \end{equation}

the critical fractile. read it as a service level: stock to cover demand with probability \(c_u / (c_u + c_o)\). expensive shortages push \(q^{*}\) into the right tail; expensive leftovers pull it back. for discrete or non-invertible \(F\), take the generalised inverse (smallest \(q\) with \(F(q) \ge\) fractile) — which is exactly why the SAA solution below is a sample quantile.

numerical verification#

exponential demand with mean \(\mu\) keeps everything analytic: \(F(q) = 1 - e^{-q/\mu}\), so \(q^{*} = \mu \ln\!\big( (c_u + c_o)/c_o \big)\). with \(c_u = 3\), \(c_o = 1.5\) the fractile is \(2/3\) and \(q^{*} = \mu \ln 3 \approx 109.86\). four independent routes to the same number — closed form, SAA (which for the newsvendor is the empirical \(2/3\)-quantile), a brute-force monte-carlo sweep over \(q\), and an honest two-stage deterministic-equivalent LP with 2000 scenarios and 4001 variables:

import numpy as np
from scipy.optimize import linprog

rng = np.random.default_rng(7)

# newsvendor: order q before demand D is seen.
# underage cost cu per unit short, overage cost co per unit left over.
cu, co = 3.0, 1.5                       # lost margin vs writedown
mu = 100.0                              # exponential demand, mean mu
frac = cu / (cu + co)                   # critical fractile = 2/3

# 1. closed form: F(q) = 1 - exp(-q/mu)  =>  q* = -mu ln(1 - frac)
q_closed = -mu * np.log(1.0 - frac)

# 2. SAA: empirical quantile of a demand sample
d = rng.exponential(mu, size=200_000)
q_saa = np.quantile(d, frac)

# 3. brute force: monte-carlo expected cost on a grid of q
qs = np.linspace(60, 160, 401)
cost = np.array([(co * np.maximum(q - d, 0) + cu * np.maximum(d - q, 0)).mean()
                 for q in qs])
q_brute = qs[cost.argmin()]

# 4. two-stage deterministic equivalent LP on S scenarios:
#    min  sum_s p_s (co*u_s + cu*v_s)  s.t.  u_s >= q - d_s,  v_s >= d_s - q
S = 2_000
ds = rng.exponential(mu, size=S)
p = np.full(S, 1.0 / S)
# variables: [q, u_1..u_S, v_1..v_S]
c_lp = np.concatenate(([0.0], co * p, cu * p))
# -q + u_s >= -d_s   ->  q - u_s <= d_s
A1 = np.hstack([np.ones((S, 1)), -np.eye(S), np.zeros((S, S))])
#  q + v_s >= d_s    -> -q - v_s <= -d_s
A2 = np.hstack([-np.ones((S, 1)), np.zeros((S, S)), -np.eye(S)])
res = linprog(c_lp, A_ub=np.vstack([A1, A2]),
              b_ub=np.concatenate([ds, -ds]),
              bounds=[(0, None)] * (1 + 2 * S), method="highs")
q_lp = res.x[0]

print(f"critical fractile   cu/(cu+co) = {frac:.4f}")
print(f"closed form         q* = mu ln 3      = {q_closed:.3f}")
print(f"SAA quantile        (n = 200000)      = {q_saa:.3f}")
print(f"brute-force sweep   (grid, same MC)   = {q_brute:.3f}")
print(f"det-equivalent LP   (S = 2000)        = {q_lp:.3f}")
print(f"expected cost at q*  (MC)             = {np.interp(q_closed, qs, cost):.3f}")

critical fractile   cu/(cu+co) = 0.6667
closed form         q* = mu ln 3      = 109.861
SAA quantile        (n = 200000)      = 109.682
brute-force sweep   (grid, same MC)   = 109.750
det-equivalent LP   (S = 2000)        = 111.678
expected cost at q*  (MC)             = 164.447

everything lines up. the 200k-sample estimates land within 0.2 of the truth; the LP, fed only 2000 scenarios, lands within 2 — and its error is pure SAA sampling error (its exact optimum is the \(2/3\) sample quantile of its own 2000 draws), a live demonstration of the caveats above. note also how flat the objective is near \(q^{*}\): convex expected costs forgive small ordering errors, which is why the crude 0.25-wide grid still gets three digits.

progressive hedging#

the L-shaped method decomposes by stage; progressive hedging decomposes by scenario. give every scenario its own full copy of the first-stage decision, \(x_s\), and solve each scenario’s deterministic problem independently — trivially parallel, but now scenario 3’s plan assumes it knows scenario 3 is coming, which is cheating. the fix is to state the no-cheating rule explicitly as a constraint, nonanticipativity: \(x_s = \bar x\) for all \(s\), where \(\bar x = \sum_s p_s x_s\) is the consensus decision. progressive hedging enforces it softly via an augmented lagrangian: at each iteration, solve the per-scenario problems with an added penalty \(w_s^{\top} x_s + \tfrac{\rho}{2} \lVert x_s - \bar x \rVert^2\), recompute the consensus \(\bar x\), and update the multipliers \(w_s \leftarrow w_s + \rho\, (x_s - \bar x)\). rockafellar and wets proved convergence for convex problems (the scheme is a proximal point algorithm in disguise);⁴ in practice it is used far beyond that warranty, as a remarkably effective heuristic when the scenario subproblems are mixed-integer and the L-shaped method’s convexity assumptions have left the building.

see also#

• linear programming — the deterministic equivalent is one of these, wearing a very wide coat
• robust optimisation — the same uncertainty without the probabilities: worst case, not average case
• integer programming — where scenario-counting chance constraints and integer recourse end up
• non-linear programming — the general machinery behind the conic and cvar reformulations

References#

Jaynes, Edwin T. (2003). Probability Theory: The Logic of Science, Cambridge University Press.

Wasserman, Larry (2010). All of Statistics: A Concise Course in Statistical Inference, Springer.