Stochastic Programming (two-stage recourse)

Stochastic Programming (two-stage recourse)#

A stochastic program optimizes an expected (or risk-adjusted) cost when part of the data is uncertain. In the two-stage setting a here-and-now first-stage decision \(x\) is made before the uncertainty \(\xi\) is realized, and a wait-and-see recourse \(y\) adapts per outcome [Birge and Louveaux, 2011, Kall and Wallace, 1994]:

\[\min_x\; c\cdot x + \mathbb{E}_\xi\big[\,Q(x,\xi)\,\big],\qquad Q(x,\xi)=\min_y\{ q(\xi)\cdot y : W(\xi)y \le h(\xi)-T(\xi)x \}.\]

This is the distributional sibling of discopt.ro: robust optimization protects the worst case over an uncertainty set, stochastic programming optimizes the average over a distribution [Shapiro et al., 2021]. With a finite scenario set \(\{\xi_s, p_s\}\) the expectation is \(\sum_s p_s Q(x,\xi_s)\).

The extensive form (deterministic equivalent)#

The simplest exact method builds one monolithic model — the deterministic equivalent — via a scenario-creator callback that constructs each scenario’s recourse. Here is the classic newsvendor: order \(q\) at cost \(c\), then in each demand scenario sell \(\min(q, d_s)\) at price \(r\) (revenue = negative cost).

import discopt.modeling as dm
from discopt.modeling.core import Constant
from discopt.stochastic import ScenarioSet, build_extensive_form, Expectation

scen = ScenarioSet.from_list([(0.3, {'demand': 40}),
                             (0.5, {'demand': 60}),
                             (0.2, {'demand': 80})])

def recourse(model, data, s):
    d = float(data['demand'])
    sales = model.continuous(f'sales_{s}', lb=0, ub=100)
    model.subject_to(sales <= q)            # cannot sell more than ordered
    model.subject_to(sales <= Constant(d))  # cannot sell more than demanded
    return Constant(-5.0) * sales           # recourse cost = -revenue (price 5)

m = dm.Model('newsvendor')
q = m.continuous('q', lb=0, ub=100)         # first-stage order
ef = build_extensive_form(m, scenarios=scen, recourse_builder=recourse,
                          first_stage_cost=Constant(2.0) * q,  # cost 2/unit
                          risk=Expectation())
print('scenarios:', len(scen), '| expected-value objective built')

scenarios: 3 | expected-value objective built

The objective is c·q + Σ_s p_s Q_s — the scenario probabilities are baked in. Solving it returns the classic critical-fractile order quantity. For a newsvendor with order cost \(c\) and price \(r\), the optimal order is the smallest \(q\) whose demand CDF reaches the critical ratio \((r-c)/r\). Here \((5-2)/5 = 0.6\), and the demand CDF is \(F(40)=0.3\), \(F(60)=0.8\), \(F(80)=1.0\), so the answer should be \(q^\star = 60\) — which the next cell checks rather than takes on trust.

The Birge–Louveaux farmer problem [Birge and Louveaux, 2011] is built the same way, with per-scenario planting/selling recourse.

res_nv = m.solve()

print("status    :", res_nv.status)
print("q*        :", float(res_nv.x["q"]))
print("objective :", float(res_nv.objective), "(negative = expected profit)")

# Critical fractile: smallest q with F(q) >= (r - c)/r = 0.6  ->  q = 60.
CRITICAL_FRACTILE_Q = 60.0
assert abs(float(res_nv.x["q"]) - CRITICAL_FRACTILE_Q) < 1e-6, res_nv.x["q"]

# And check the objective by hand: cost 2*q minus revenue 5*E[min(q, d)].
expected_sales = 0.3 * min(60, 40) + 0.5 * min(60, 60) + 0.2 * min(60, 80)
expected_cost = 2.0 * CRITICAL_FRACTILE_Q - 5.0 * expected_sales
print("hand-computed objective:", expected_cost)
assert abs(float(res_nv.objective) - expected_cost) < 1e-6

print(f"\nMatches the critical-fractile solution q* = {CRITICAL_FRACTILE_Q:.0f}.")

status    : optimal
q*        : 60.0
objective : -150.0 (negative = expected profit)
hand-computed objective: -150.0

Matches the critical-fractile solution q* = 60.

Risk measures: CVaR and chance constraints#

Beyond the risk-neutral expectation, discopt.stochastic provides CVaR (Conditional Value-at-Risk) via the Rockafellar–Uryasev formulation [Rockafellar and Uryasev, 2000], which penalizes the worst-tail scenarios:

from discopt.stochastic import CVaR

mc = dm.Model('newsvendor_cvar')
q = mc.continuous('q', lb=0, ub=100)
ef_cvar = build_extensive_form(mc, scenarios=scen, recourse_builder=recourse,
                               first_stage_cost=Constant(2.0) * q,
                               risk=CVaR(alpha=0.9))   # penalize worst 10%
aux = [v.name for v in mc._variables if v.name.startswith('cvar')]
print('CVaR auxiliaries (VaR + per-scenario shortfalls):', aux)

CVaR auxiliaries (VaR + per-scenario shortfalls): ['cvar_eta', 'cvar_u0', 'cvar_u1', 'cvar_u2']

A chance constraint \(P(g(x,\xi)\le 0)\ge 1-\varepsilon\) is imposed (under SAA) with a per-scenario indicator and a coverage budget:

from discopt.stochastic import chance_constraint

mch = dm.Model('chance'); x = mch.continuous('x', lb=0, ub=100)
demands = [30.0, 50.0, 70.0, 90.0]
z = chance_constraint(mch, [x - Constant(d) for d in demands],
                      probs=[0.25]*4, epsilon=0.25, big_m=1000.0)
print('indicator binaries:', len(z), '(Σ p_s z_s ≤ ε caps the violation prob.)')

indicator binaries: 4 (Σ p_s z_s ≤ ε caps the violation prob.)

Decomposition: the L-shaped method#

For many scenarios the extensive form is huge. The L-shaped method [Van Slyke and Wets, 1969] is Benders decomposition for SP: the first stage is the master, each scenario a subproblem, and cuts approximate the expected recourse. The probability weighting rides in the objective, so the decomposition engine is reused unchanged. solve_lshaped(..., solve=False) builds the model and its structure:

from discopt.stochastic import solve_lshaped

ml = dm.Model('newsvendor_lshaped')
q = ml.continuous('q', lb=0, ub=100)
res = solve_lshaped(ml, first_stage_vars=[q], scenarios=scen,
                    recourse_builder=recourse,
                    first_stage_cost=Constant(2.0) * q, solve=False)
recourse_blocks = [b for b in res.structure.blocks if 'q' not in b]
print('first-stage (complicating):', res.structure.complicating_vars)
print('separable recourse blocks:', len(recourse_blocks))

first-stage (complicating): ['q']
separable recourse blocks: 3

Progressive hedging#

Progressive hedging [Rockafellar and Wets, 1991] is a scenario-decomposition alternative: each scenario keeps its own first-stage copy, and a proximal term drives the copies to a common nonanticipative value \(\bar x = \sum_s p_s x_s\). The driver takes an injected subproblem solver; here a closed-form quadratic illustrates that PH converges to the expected-value decision:

import numpy as np
from discopt.stochastic import progressive_hedging, quadratic_subproblem_solver

targets = {0: np.array([40.0]), 1: np.array([60.0]), 2: np.array([80.0])}
ph = progressive_hedging(scen, first_stage_dim=1,
                         subproblem_solver=quadratic_subproblem_solver(targets),
                         rho=1.0)
print('x_bar =', float(ph.x_bar[0]), '(= Σ p_s a_s = 58) | converged:', ph.converged)

x_bar = 58.00000000000002 (= Σ p_s a_s = 58) | converged: True

How good is the sample? SAA statistical bounds#

When scenarios are sampled (SAA), the sample optimum is a biased-low estimate of the true optimum. Averaging independent batches gives a statistical lower bound, and evaluating a candidate on a fresh sample gives an upper bound — together an optimality-gap estimate [Mak et al., 1999]:

from discopt.stochastic import (saa_lower_bound_estimate,
                                saa_upper_bound_estimate, optimality_gap_estimate)
rng = np.random.default_rng(0)
lower = saa_lower_bound_estimate(rng.normal(100, 5, 20), confidence=0.95)
upper = saa_upper_bound_estimate(rng.normal(103, 4, 50), confidence=0.95)
gap = optimality_gap_estimate(lower, upper)
print('point gap:', round(gap['point_gap'], 2),
      '| conservative gap:', round(gap['conservative_gap'], 2))

point gap: 4.66 | conservative gap: 7.27

References#

See the References page for full citations.