Farmer: Further Modeling Examples - Optimization for Decision Science

Optional extensions to the core farmer example: uniform-yield sampling, dependence, profit floors, other risk measures, and uncertain prices. For the main Lecture 10 reading, use SAA, CVaR, and zero VSS.

import sys
from pathlib import Path
if "google.colab" in sys.modules:
    from urllib.request import urlretrieve
    base = "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/"
    for name in ["helper.py", "farmer.py"]:
        urlretrieve(base + name, name)
    import helper
    helper.easy_install()
else:
    sys.path.insert(0, str(Path.cwd().parent))
    import helper
import farmer as f
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
helper.set_plotting_style()

from farmer import *
mean_model=solve_farmer(build_farmer([correlated_scenario(1)]))
mean_plan=plan_of(mean_model)
RP=solve_farmer(build_farmer(THREE))
rp_plan=plan_of(RP)
rp_expected=pyo.value(RP.expected_profit)
rp_profits=np.array([pyo.value(RP.scenarios[i].profit) for i in RP.SCENARIOS])

Sample average approximation: train and evaluate¶

For the continuous-uniform model (§1.1d, pp. 11–15), replace the expectation by an average over independent sampled yields. Reuse build_farmer with one Block per sample. Choose acreage using a training sample, then fix acreage and optimize trading on a new evaluation sample.

The table uses a common evaluation sample, independent of all training samples, to compare plans. Each reported standard error estimates uncertainty in that fixed plan’s mean profit. It is not a confidence bound on the optimality gap. Analyze whether bigger training samples improve these particular decisions; monotonic improvement in any finite sequence is not guaranteed.

def sample_farmer_yields(n, seed):
    """Independent uniform yields, T/acre; Birge & Louveaux, Section 1.1d."""
    rng = np.random.default_rng(seed)
    return [
        {"yield": {c: rng.uniform(0.8, 1.2) * MEAN_YIELD[c] for c in CROPS}}
        for _ in range(n)
    ]


def build_sampled_farmer(n, seed):
    """Build an n-scenario SAA model with reproducible, equally weighted draws."""
    # Equal weights 1/n are the builder's default.
    return build_farmer(sample_farmer_yields(n, seed))

evaluation_sample = sample_farmer_yields(2000, seed=20260922)
trained = solve_farmer(build_sampled_farmer(200, seed=12))
evaluated, profits = grade_plan(plan_of(trained), evaluation_sample)
se = profits.std(ddof=1) / np.sqrt(len(profits))

print("N train   wheat    corn   beets   train mean   eval mean   eval SE ($)")
for n_train, seed in [(50, 11), (200, 12), (1000, 13)]:
    trained = solve_farmer(build_sampled_farmer(n_train, seed))
    plan = plan_of(trained)
    evaluated, profits = grade_plan(plan, evaluation_sample)
    se = profits.std(ddof=1) / np.sqrt(len(profits))
    print(
        f"{n_train:7d} {plan[0]:7.1f} {plan[1]:7.1f} {plan[2]:7.1f}"
        f" {pyo.value(trained.expected_profit):12,.0f} {evaluated:11,.0f} {se:13,.0f}"
    )

N train   wheat    corn   beets   train mean   eval mean   eval SE ($)

     50   129.9    86.4   283.7      105,724     111,611           460

    200   136.1    85.7   278.2      110,180     111,713           476

   1000   135.0    84.6   280.4      110,314     111,694           471

Optional depth: discrete yields and continuous yields¶

Birge and Louveaux §1.1b, pp. 6–8, uses three equally likely, correlated yield vectors. Section 1.1d, pp. 11–15, uses independent uniforms on 80–120% of each mean yield. Their means match, but their marginal distributions do not. The optimal plans differ: (170, 80, 250) for the three-point model versus about (135.83, 85.07, 279.10) for the continuous model.

For a separate, controlled dependence comparison, the samplers below both use uniform marginals. One draws crops independently; the other uses a shared continuous weather factor. The correlated-uniform model is an instructor extension, not the three-point model in §1.1b.

Hold marginals fixed to isolate dependence¶

The farmer’s recourse separates across crops. Therefore the expected recourse is a sum of crop-wise expectations, each determined by that crop’s marginal (Birge and Louveaux, §1.1d, Eq. (1.6), p. 11).

Independent and correlated uniform yields consequently give the same expected-profit function and risk-neutral optimal plan. The three-point model has different marginals, so this argument does not equate its solution with the continuous model. Equal means alone are insufficient.

Sampling the continuous model, and checking it without a solver¶

Section 1.1d has no scenario list: the yields are continuous. Approximate the expectation by sampling, one pyo.Block per draw, and seed every draw so the cell reproduces on its own.

The check is the interesting part. Birge and Louveaux solve the continuous model analytically: pp. 12--14 give closed forms for $Q_1$ , $Q_2$ and $Q_3$ , and p. 15 solves the resulting Karush--Kuhn--Tucker conditions to $\lambda = 275.00$ , $\mathbf x = (135.83,\ 85.07,\ 279.10)$ . So this model can be checked by a route that uses no solver at all — a luxury worth taking whenever it is available.

# Two samplers, seeded. The only difference is whether one draw is shared by all
# three crops (an instructor extension) or each crop is drawn on its own
# (Section 1.1d's independence). Both give each crop the SAME marginal:
# Uniform(0.8 * mean, 1.2 * mean).
YIELD_BOUNDS = {c: (0.8 * MEAN_YIELD[c], 1.2 * MEAN_YIELD[c]) for c in CROPS}


def sample_independent(n, seed):
    """Section 1.1d, p. 11: each crop's yield drawn independently."""
    rng = np.random.default_rng(seed)
    return [
        {"yield": {c: rng.uniform(*YIELD_BOUNDS[c]) for c in CROPS}} for _ in range(n)
    ]


def sample_correlated(n, seed):
    """Instructor extension: one common continuous weather factor."""
    rng = np.random.default_rng(seed)
    return [correlated_scenario(f) for f in rng.uniform(0.8, 1.2, size=n)]


# --- the solver-free reference: B&L's closed-form expected recourse, pp. 13-14
def Q_wheat(x):
    if x <= 200 / 3:
        return 47600 - 595 * x
    if x <= 100:
        return 119 * (200 - 2 * x) ** 2 / x - 85 * (200 - 3 * x) ** 2 / x
    return 34000 - 425 * x


def Q_corn(x):
    if x <= 200 / 3:
        return 50400 - 630 * x
    if x <= 100:
        return 87.5 * (240 - 2.4 * x) ** 2 / x - 62.5 * (240 - 3.6 * x) ** 2 / x
    return 36000 - 450 * x


def Q_beets(x):
    """Equations (1.10)-(1.12), pp. 12-13, with l3 = 16 and u3 = 24 T/acre."""
    lo, hi = 16.0, 24.0
    if x <= 0:
        return 0.0
    if lo * x > 6000:  # quota exceeded for every yield
        return -156000 - 10 * MEAN_YIELD["BEETS"] * x
    if hi * x < 6000:  # quota never reached
        return -36 * MEAN_YIELD["BEETS"] * x
    return -36 * MEAN_YIELD["BEETS"] * x + 13 * (hi * x - 6000) ** 2 / (x * (hi - lo))


def analytic_expected_profit(plan):
    """Exact E[profit] of a plan under continuous uniform yields -- no solver."""
    planting = sum(PLANT_COST[c] * a for c, a in zip(CROPS, plan))
    return -(planting + Q_wheat(plan[0]) + Q_corn(plan[1]) + Q_beets(plan[2]))


book_plan = np.array([135.83, 85.07, 279.10])
print(
    f"B&L KKT solution, p. 15:  ({book_plan[0]}, {book_plan[1]}, {book_plan[2]})"
    f"   analytic E[profit] = ${analytic_expected_profit(book_plan):,.2f}"
)
print()
print(
    f"{'sampler':>12s} {'N':>6s} {'wheat':>8s} {'corn':>8s} {'beets':>8s}"
    f" {'sampled E':>12s} {'analytic E':>12s}"
)
for sampler, label in [
    (sample_independent, "independent"),
    (sample_correlated, "correlated"),
]:
    for n, seed in [(200, 1), (1000, 2), (4000, 3)]:
        m = solve_farmer(build_farmer(sampler(n, seed)))
        p = plan_of(m)
        print(
            f"{label:>12s} {n:6d} {p[0]:8.2f} {p[1]:8.2f} {p[2]:8.2f}"
            f" {pyo.value(m.expected_profit):12,.0f} {analytic_expected_profit(p):12,.0f}"
        )

B&L KKT solution, p. 15:  (135.83, 85.07, 279.1)   analytic E[profit] = $111,237.44

     sampler      N    wheat     corn    beets    sampled E   analytic E

 independent    200   134.31    84.21   281.49      110,390      111,220

 independent   1000   137.32    84.96   277.71      111,730      111,232

 independent   4000   135.64    84.76   279.60      110,604      111,237

  correlated    200   137.09    84.47   278.44      112,259      111,235

  correlated   1000   135.27    85.15   279.58      111,452      111,237

  correlated   4000   136.23    85.29   278.48      110,886      111,236

Both samplers target the continuous-uniform optimum; finite samples yield different plans. The analytic column evaluates a sampled plan under the true uniform marginals. Neither a particular sample objective nor a sequence of sampled decisions must improve monotonically with sample size. Training objectives are optimistic in expectation, not on every run.

# Hold the plan FIXED at the book's solution and look at the distribution of
# profit it produces under each assumption. Same decision, same mean -- the only
# difference is whether the three crops fail together.
N_EVAL, EVAL_SEED = 20_000, 20260921


def tail_summary(profits, alpha):
    """Empirical mean, spread and lower-tail risk of a profit sample."""
    s = np.sort(profits)
    mass = (1 - alpha) * len(s)
    k = min(len(s), int(np.floor(mass + 1e-10)))
    fractional = max(0.0, mass - k)
    tail_mean = (s[:k].sum() + (fractional * s[k] if k < len(s) else 0)) / mass
    cutoff = max(0, int(np.ceil(mass - 1e-10)) - 1)
    return {
        "mean": s.mean(),
        "sd": s.std(ddof=1),
        "min": s[0],
        f"VaR_{alpha}": s[cutoff],
        f"CVaR_{alpha}": tail_mean,
    }


print(
    f"Plan held fixed at ({book_plan[0]}, {book_plan[1]}, {book_plan[2]}) acres;"
    f" {N_EVAL:,} seeded draws each."
)
print(
    f"True E[profit] for BOTH, analytically: ${analytic_expected_profit(book_plan):,.2f}"
)
print()
print(
    f"{'assumption':>22s} {'E (sampled)':>13s} {'std dev':>10s} {'min':>10s}"
    f" {'VaR_0.95':>10s} {'CVaR_0.95':>10s}"
)
summaries = {}
for label, sampler in [
    ("independent (1.1d)", sample_independent),
    ("correlated uniform", sample_correlated),
]:
    _, profits = grade_plan(book_plan, sampler(N_EVAL, EVAL_SEED))
    s = tail_summary(profits, 0.95)
    summaries[label] = s
    print(
        f"{label:>22s} {s['mean']:13,.0f} {s['sd']:10,.0f} {s['min']:10,.0f}"
        f" {s['VaR_0.95']:10,.0f} {s['CVaR_0.95']:10,.0f}"
    )

gap = (
    summaries["independent (1.1d)"]["CVaR_0.95"]
    - summaries["correlated uniform"]["CVaR_0.95"]
)
print(f"\nSame plan, same true mean, but CVaR_0.95 differs by ${gap:,.0f}.")
print("The two sampled means differ only by Monte Carlo noise; the spread does not.")

Plan held fixed at (135.83, 85.07, 279.1) acres; 20,000 seeded draws each.
True E[profit] for BOTH, analytically: $111,237.44

            assumption   E (sampled)    std dev        min   VaR_0.95  CVaR_0.95

    independent (1.1d)       111,106     21,294     56,525     74,911     69,317

    correlated uniform       110,667     31,035     52,928     59,198     56,059

Same plan, same true mean, but CVaR_0.95 differs by $13,258.
The two sampled means differ only by Monte Carlo noise; the spread does not.

At this fixed plan, the independent-uniform and correlated-uniform models have the same true expected profit but different downside risk. Their sampled means need not coincide. This comparison isolates dependence because the marginals are held fixed; it does not compare §1.1b directly with §1.1d.

Risk aversion: the farmer’s efficient frontier¶

Birge and Louveaux pose this as Exercise 7 of §1.1, pp. 18--19:

“Economic theory tells us that, like many other people, the farmer would normally act as a risk-averse person. … One simple way is to plan for the worst case.”

They note that the worst situation is Scenario 3 (below-average yields), and that planning for it makes the Table 4 solution optimal — 100 acres of wheat, 25 of corn, 375 of sugar beets.

(a) compute the loss in expected profit if that plan is taken;
(b) instead require the worst-case profit not to fall below $58,000, and maximize expected profit subject to that;
(c) repeat at $56,000, $54,000, $52,000, $50,000 and $48,000, and graph the curve of expected profit loss.

Everything needed is already in hand. No new data, no new setting.

# Exercise 7(a): grade the worst-case plan of Table 4, p. 7. Recourse is still
# optimized in every scenario -- only the acreage is fixed.
table_4_plan = np.array([100.0, 25.0, 375.0])
e_worst_case, profits_worst_case = grade_plan(table_4_plan, THREE)

print("Table 4 plan (plan for the worst case): 100 / 25 / 375 acres")
print(
    "  scenario profits  above ${:,.0f} | average ${:,.0f} | below ${:,.0f}".format(
        *profits_worst_case
    )
)
print(f"  expected profit = ${e_worst_case:,.0f}")
print(
    f"\nExercise 7(a): loss in expected profit = "
    f"${rp_expected:,.0f} - ${e_worst_case:,.0f} = ${rp_expected - e_worst_case:,.0f}"
)

# The below-average profit is printed in the book: Table 4, p. 7, gives $59,950.
assert abs(profits_worst_case.min() - 59950.0) < 1.0

Table 4 plan (plan for the worst case): 100 / 25 / 375 acres
  scenario profits  above $113,250 | average $86,600 | below $59,950
  expected profit = $86,600

Exercise 7(a): loss in expected profit = $108,390 - $86,600 = $21,790

protected = solve_farmer(build_farmer(THREE, risk={"floor": 56_000}))
protected_mean, protected_profits = grade_plan(plan_of(protected), THREE)

# Exercises 7(b) and 7(c): raise the floor and watch what the insurance costs.
print(
    f"{'floor':>10s} {'wheat':>8s} {'corn':>8s} {'beets':>8s}"
    f" {'E[profit]':>12s} {'loss vs RP':>12s}"
)
frontier = []
for floor in [48_000, 50_000, 52_000, 54_000, 56_000, 58_000]:
    m = solve_farmer(build_farmer(THREE, risk={"floor": floor}))
    p, e = plan_of(m), pyo.value(m.expected_profit)
    frontier.append((floor, p, e))
    print(
        f"{floor:10,d} {p[0]:8.1f} {p[1]:8.1f} {p[2]:8.1f} {e:12,.0f}"
        f" {rp_expected - e:12,.0f}"
    )

print("\nMarginal cost of each extra $2,000 of floor:")
for (f0, _, e0), (f1, _, e1) in zip(frontier, frontier[1:]):
    print(f"  ${f0:,} -> ${f1:,}:  ${e0 - e1:>8,.0f}")

     floor    wheat     corn    beets    E[profit]   loss vs RP

    48,000    170.0     80.0    250.0      108,390            0

    50,000    156.0     94.0    250.0      108,292           98

    52,000    138.1    100.0    261.9      107,976          414

    54,000    122.2    100.0    277.8      107,611          779

    56,000    106.3    100.0    293.7      107,246        1,144

    58,000    100.0     71.4    328.6      101,176        7,214

Marginal cost of each extra $2,000 of floor:
  $48,000 -> $50,000:  $      98
  $50,000 -> $52,000:  $     315
  $52,000 -> $54,000:  $     365
  $54,000 -> $56,000:  $     365
  $56,000 -> $58,000:  $   6,070

The curve has a kink, and the kink is the discussion.

From $50,000 to $56,000 each extra $2,000 of guaranteed profit costs a near-constant $365 or so. The last step to $58,000 costs $6,070 — sixteen times the marginal rate.

The acreage column explains it. Through the linear stretch, corn sits pinned at its upper limit of 100 acres and the adjustment is made by trading wheat for beets. At $58,000 corn falls off that bound, dropping to 71.4 acres while beets jump to 328.6, and the plan collapses toward the worst-case allocation of 100 / 25 / 375.

Insurance is cheap until it isn’t.

From a worst-case floor to CVaR¶

Exercise 7’s floor turns out to be a conditional value at risk constraint in disguise, and seeing why makes the general CVaR machinery arrive as a generalization of something already understood rather than as new construction.

Write CVaR on the lower tail of profit — the mirror of the loss convention in Birge and Louveaux (9.7), p. 85:

\mathrm{CVaR}_\alpha(\pi)=\max_{\nu}\left\{\nu-\frac{1}{1-\alpha}\, \mathbb E\big[(\nu-\pi)_+\big]\right\},

(1)

the average profit over the worst $1-\alpha$ of the probability mass.

With $S$ equally likely scenarios and tail mass $1-\alpha \le 1/S$ , the worst scenario alone carries the entire tail, so $\mathrm{CVaR}_\alpha(\pi) = \min_s \pi_s$ and the CVaR floor is the worst-case floor. For $S = 3$ that means every $\alpha \ge 2/3$ .

# Verify the equivalence rather than asserting it: solve BOTH models and diff.
print(
    f"{'alpha':>8s} {'floor':>9s} {'E (floor model)':>17s} {'E (CVaR model)':>16s}"
    f" {'dE':>8s} {'max|d acres|':>13s}"
)
for alpha in [2 / 3, 0.75, 0.90, 0.95]:
    for floor, plan_floor, e_floor in frontier:
        if floor < 52_000:
            continue
        m = solve_farmer(build_farmer(THREE, risk={"cvar_floor": (alpha, floor)}))
        e_cvar, plan_cvar = pyo.value(m.expected_profit), plan_of(m)
        d_acres = np.abs(plan_cvar - plan_floor).max()
        print(
            f"{alpha:8.3f} {floor:9,d} {e_floor:17,.2f} {e_cvar:16,.2f}"
            f" {e_cvar - e_floor:8.2f} {d_acres:13.2e}"
        )
        assert abs(e_cvar - e_floor) < 1e-4 and d_acres < 1e-6

   alpha     floor   E (floor model)   E (CVaR model)       dE  max|d acres|

   0.667    52,000        107,976.19       107,976.19     0.00      0.00e+00

   0.667    54,000        107,611.11       107,611.11     0.00      0.00e+00

   0.667    56,000        107,246.03       107,246.03     0.00      0.00e+00

   0.667    58,000        101,176.19       101,176.19     0.00      0.00e+00

   0.750    52,000        107,976.19       107,976.19     0.00      0.00e+00

   0.750    54,000        107,611.11       107,611.11     0.00      0.00e+00

   0.750    56,000        107,246.03       107,246.03     0.00      0.00e+00

   0.750    58,000        101,176.19       101,176.19     0.00      0.00e+00

   0.900    52,000        107,976.19       107,976.19     0.00      0.00e+00

   0.900    54,000        107,611.11       107,611.11     0.00      0.00e+00

   0.900    56,000        107,246.03       107,246.03     0.00      0.00e+00

   0.900    58,000        101,176.19       101,176.19     0.00      0.00e+00

   0.950    52,000        107,976.19       107,976.19     0.00      0.00e+00

   0.950    54,000        107,611.11       107,611.11     0.00      0.00e+00

   0.950    56,000        107,246.03       107,246.03     0.00      0.00e+00

   0.950    58,000        101,176.19       101,176.19     0.00      0.00e+00

Partial mass at the cutoff, on a model you already own¶

The equivalence above is a degeneracy of three equally likely scenarios, not a general fact — and the quickest way to see CVaR do something the minimum cannot is to drop $\alpha$ below $2/3$ so that the tail splits an atom.

At the risk-neutral plan, the three equally likely profits are $48,820, $109,350 and $167,000. At $\alpha = 0.5$ the worst half of the probability mass is all of the below-average scenario ( $1/3$ ) plus half of the average one ( $1/6$ ). A conditional expectation cannot do this: it must take whole scenarios, so it averages $2/3$ of the mass when $1/2$ was asked for.

# Hand arithmetic first, then the LP epigraph, then check they agree.
sorted_profits = np.sort(rp_profits)
var_50 = sorted_profits[1]  # smallest t with P(profit <= t) >= 0.5
cvar_50 = (1 / 3 * sorted_profits[0] + 1 / 6 * sorted_profits[1]) / 0.5
conditional_mean = sorted_profits[sorted_profits <= var_50].mean()

print(
    "profits at the risk-neutral plan: "
    + " | ".join(f"${p:,.0f}" for p in sorted_profits)
)
print(f"  VaR_0.5  (lower tail)          = ${var_50:,.2f}")
print(
    f"  CVaR_0.5 = (1/3 x {sorted_profits[0]:,.0f} + 1/6 x {sorted_profits[1]:,.0f}) / 0.5"
    f" = ${cvar_50:,.2f}"
)
print(
    f"  E[profit | profit <= VaR_0.5]  = ${conditional_mean:,.2f}"
    "   <- averages 2/3 of the mass, not 1/2"
)

# Second route: the Rockafellar-Uryasev epigraph solved at the same fixed plan.
m = build_farmer(THREE, risk={"cvar": (0.5, 1.0)})
for crop, acres in zip(CROPS, rp_plan):
    m.acreage[crop].fix(acres)
solve_farmer(m)
print(f"\n  LP epigraph at the same plan   = ${pyo.value(m.cvar):,.2f}")
assert abs(pyo.value(m.cvar) - cvar_50) < 1e-6

profits at the risk-neutral plan: $48,820 | $109,350 | $167,000
  VaR_0.5  (lower tail)          = $109,350.00
  CVaR_0.5 = (1/3 x 48,820 + 1/6 x 109,350) / 0.5 = $68,996.67
  E[profit | profit <= VaR_0.5]  = $79,085.00   <- averages 2/3 of the mass, not 1/2


  LP epigraph at the same plan   = $68,996.67

When CVaR is the expected value, and when it is the worst case¶

CVaR prices the worst $1-\alpha$ of the probability mass, so the two ends of the $\alpha$ range are not really CVaR problems at all:

limit	tail mass	CVaR becomes	efficient formulation
$\alpha \to 0$	everything	$\mathbb E[\pi]$	drop $\nu$ and every shortfall variable; take the mean
$\alpha \to 1$	the single worst outcome	$\min_s \pi_s$	minimax: one variable $t$ with $t \le \pi_s$ for all $s$

Both limits are reached at finite $\alpha$ when the scenarios are finitely many and equally likely: with $S$ scenarios, any $\alpha \ge 1 - 1/S$ already gives exactly the worst case, so pushing $\alpha$ higher buys nothing.

⚠ And pushing it higher actively costs something. The epigraph carries a $1/(1-\alpha)$ coefficient, which grows without bound as $\alpha \to 1$ : at $\alpha = 0.999$ it is 1,000, and at $\alpha = 0.9999$ it is 10,000. Near the worst-case end the general formulation is therefore both larger and worse conditioned than the minimax model it is approximating. Switching formulation is the right move; raising $\alpha$ is not.

# The limits, verified rather than asserted. Read the PLAN off the risk solve,
# then grade it with a separate expectation solve (the next section explains why).
print(
    f"{'alpha':>8s} {'CVaR*':>12s} {'wheat':>7s} {'corn':>7s} {'beets':>7s}"
    f" {'E[profit]':>12s} {'worst':>10s}"
)
for alpha in [0.0, 0.1, 1 / 3, 0.5, 2 / 3, 0.9, 0.99]:
    m = solve_farmer(build_farmer(THREE, risk={"cvar": (alpha, 1.0)}))
    p = plan_of(m)
    e, scenario_profits = grade_plan(p, THREE)
    print(
        f"{alpha:8.3f} {pyo.value(m.cvar):12,.2f} {p[0]:7.1f} {p[1]:7.1f} {p[2]:7.1f}"
        f" {e:12,.0f} {scenario_profits.min():10,.0f}"
    )

print(
    f"\nalpha = 0    : CVaR* = ${rp_expected:,.2f} = the risk-neutral expected profit"
)
print(
    f"alpha >= 2/3 : CVaR* = ${profits_worst_case.min():,.2f} = the worst scenario profit"
)

   alpha        CVaR*   wheat    corn   beets    E[profit]      worst

   0.000   108,390.00   170.0    80.0   250.0      108,390     48,820

   0.100   102,711.11   120.0    80.0   300.0      107,240     55,120

   0.333    87,150.00   100.0   100.0   300.0      107,100     56,800

   0.500    77,033.33   100.0   100.0   300.0      107,100     56,800

   0.667    59,950.00   100.0    25.0   375.0       86,600     59,950

   0.900    59,950.00   100.0    25.0   375.0       86,600     59,950

   0.990    59,950.00   100.0    25.0   375.0       86,600     59,950

alpha = 0    : CVaR* = $108,390.00 = the risk-neutral expected profit
alpha >= 2/3 : CVaR* = $59,950.00 = the worst scenario profit

# What the purpose-built models cost, against the general epigraph.
def size_of(m):
    """Active variables and constraints, scenario blocks included."""
    return (
        sum(1 for _ in m.component_data_objects(pyo.Var, active=True)),
        sum(1 for _ in m.component_data_objects(pyo.Constraint, active=True)),
    )


expectation_model = solve_farmer(build_farmer(THREE))
epigraph_at_0 = solve_farmer(build_farmer(THREE, risk={"cvar": (0.0, 1.0)}))
minimax_model = solve_farmer(build_farmer(THREE, risk={"minimax": True}))
epigraph_at_high = solve_farmer(build_farmer(THREE, risk={"cvar": (0.99, 1.0)}))

rows = [
    (
        "expectation (alpha -> 0 limit)",
        expectation_model,
        pyo.value(expectation_model.expected_profit),
    ),
    ("epigraph at alpha = 0", epigraph_at_0, pyo.value(epigraph_at_0.cvar)),
    ("minimax (alpha -> 1 limit)", minimax_model, pyo.value(minimax_model.worst)),
    ("epigraph at alpha = 0.99", epigraph_at_high, pyo.value(epigraph_at_high.cvar)),
]
print(f"{'model':>32s} {'objective':>12s} {'vars':>6s} {'cons':>6s}")
for label, m, value in rows:
    v, c = size_of(m)
    print(f"{label:>32s} {value:12,.2f} {v:6d} {c:6d}")

print("\nSame two answers, reached twice. The epigraph adds one threshold variable")
print("and one shortfall variable per scenario, plus one constraint per scenario,")
print("to compute numbers that a mean and a minimax deliver directly.")
print("\n1/(1-alpha), the epigraph coefficient:")
for alpha in [0.5, 0.9, 0.99, 0.999, 0.9999]:
    print(f"  alpha = {alpha:<8.4f}  1/(1-alpha) = {1 / (1 - alpha):>10,.1f}")

                           model    objective   vars   cons
  expectation (alpha -> 0 limit)   108,390.00     21     10
           epigraph at alpha = 0   108,390.00     25     13
      minimax (alpha -> 1 limit)    59,950.00     22     13
        epigraph at alpha = 0.99    59,950.00     25     13

Same two answers, reached twice. The epigraph adds one threshold variable
and one shortfall variable per scenario, plus one constraint per scenario,
to compute numbers that a mean and a minimax deliver directly.

1/(1-alpha), the epigraph coefficient:
  alpha = 0.5000    1/(1-alpha) =        2.0
  alpha = 0.9000    1/(1-alpha) =       10.0
  alpha = 0.9900    1/(1-alpha) =      100.0
  alpha = 0.9990    1/(1-alpha) =    1,000.0
  alpha = 0.9999    1/(1-alpha) =   10,000.0

⚠ A pure-CVaR objective leaves the non-tail scenarios undetermined¶

This is a different problem from the limits above, and it is easy to conflate them. The limits are about what CVaR equals. This is about what the solver reports.

A pure-CVaR objective ( $\lambda = 1$ ) references only the tail scenarios. Every other scenario’s recourse variables appear nowhere in the objective, so the solver is free to return any feasible values for them. The model is solved correctly; the per-scenario numbers you read off it are simply not pinned down.

The two ideas share one root — CVaR only prices the tail — and the $\alpha \to 1$ end is where they meet: exactly one scenario is priced, so the recourse in all the others is unpriced and free.

# Demonstrate it. Read the scenario profits straight off a pure-CVaR solve, then
# read them again after fixing the plan and re-optimizing recourse.
print(
    f"{'lambda':>7s} {'wheat':>7s} {'corn':>7s} {'beets':>7s}"
    f" {'E (read off)':>14s} {'E (re-optimized)':>18s} {'gap':>10s}"
)
for lam in [0.0, 0.25, 0.5, 0.75, 0.9, 1.0]:
    m = solve_farmer(build_farmer(THREE, risk={"cvar": (2 / 3, lam)}))
    p = plan_of(m)
    read_off = np.array([pyo.value(m.scenarios[i].profit) for i in m.SCENARIOS])
    e_true, _ = grade_plan(p, THREE)
    print(
        f"{lam:7.2f} {p[0]:7.1f} {p[1]:7.1f} {p[2]:7.1f}"
        f" {read_off.mean():14,.0f} {e_true:18,.0f} {e_true - read_off.mean():10,.0f}"
    )

# The same trap in a model that contains no CVaR at all.
mm = solve_farmer(build_farmer(THREE, risk={"minimax": True}))
mm_read_off = np.array([pyo.value(mm.scenarios[i].profit) for i in mm.SCENARIOS])
mm_true, mm_profits = grade_plan(plan_of(mm), THREE)
print(f"\nminimax (B&L 9.15): worst = ${pyo.value(mm.worst):,.0f}")
print(
    f"  scenario profits read off the solve : "
    + " | ".join(f"${p:,.0f}" for p in mm_read_off)
)
print(
    f"  same plan, recourse re-optimized    : "
    + " | ".join(f"${p:,.0f}" for p in mm_profits)
)
print(f"  error if the first line were reported: ${mm_true - mm_read_off.mean():,.0f}")

 lambda   wheat    corn   beets   E (read off)   E (re-optimized)        gap

   0.00   170.0    80.0   250.0        108,390            108,390          0

   0.25   100.0   100.0   300.0        107,100            107,100          0

   0.50   100.0   100.0   300.0        107,100            107,100          0

   0.75   100.0   100.0   300.0        107,100            107,100          0

   0.90   100.0    25.0   375.0         86,600             86,600          0

   1.00   100.0    25.0   375.0         59,950             86,600     26,650


minimax (B&L 9.15): worst = $59,950
  scenario profits read off the solve : $59,950 | $59,950 | $59,950
  same plan, recourse re-optimized    : $113,250 | $86,600 | $59,950
  error if the first line were reported: $26,650

Every row with $\lambda < 1$ agrees, because a positive weight on the expectation prices all three scenarios. At $\lambda = 1$ the two columns differ by tens of thousands of dollars, and nothing errors — the solve terminates optimally and the reported numbers are feasible.

The minimax model shows the same failure without any CVaR in it, which is the proof that this is about which scenarios the objective touches, not about CVaR.

The rule: if you are going to report per-scenario recourse from a risk solve, pin it down first. Three ways, in increasing order of effort:

Do not report it. Report the plan, which is determined, and the risk measure, which is the objective.
Keep a small positive weight on the expectation ( $\lambda < 1$ ). This is what the sweep above does, and it is usually what you wanted anyway.
Grade the plan in a second solve with the acreage fixed and expected profit as the objective — what grade_plan does throughout this notebook.

The general statement covers more than CVaR: any objective that references only part of the scenario set leaves the rest of the recourse undetermined.

Compare expected profit and tail protection¶

The table selects acreage with weighted expected profit and lower-tail profit CVaR, then uses grade_plan to optimize recourse in every scenario. This matters especially at weight 1: pure CVaR need not determine non-tail trades. With three equally likely scenarios and alpha = 2/3, the tail is the worst scenario. More general alpha values can use only part of a scenario’s mass.

selected = solve_farmer(build_farmer(THREE, risk={"cvar": (2/3, 0.25)}))
plan = plan_of(selected)
mean, profits = grade_plan(plan, THREE)

print("weight   wheat   corn   beets   mean profit   worst profit ($)")
for weight in [0.0, 0.25, 0.5, 1.0]:
    selected = solve_farmer(build_farmer(THREE, risk={"cvar": (2 / 3, weight)}))
    plan = plan_of(selected)
    mean, profits = grade_plan(plan, THREE)
    print(
        f"{weight:6.2f} {plan[0]:7.1f} {plan[1]:6.1f} {plan[2]:7.1f}"
        f" {mean:13,.0f} {profits.min():16,.0f}"
    )

weight   wheat   corn   beets   mean profit   worst profit ($)

  0.00   170.0   80.0   250.0       108,390           48,820

  0.25   100.0  100.0   300.0       107,100           56,800

  0.50   100.0  100.0   300.0       107,100           56,800

  1.00   100.0   25.0   375.0        86,600           59,950

The other risk measures, and why they disagree¶

Birge and Louveaux §2.9, pp. 84--87, defines the measures this course uses: value at risk (9.3), the four coherence axioms of Artzner et al. (Definition 2.1, p. 85), CVaR (9.5)--(9.7) with its linear-programming form (9.8)--(9.10), the CVaR constraint (9.11)--(9.13), and robust optimization (9.15). Section 2.5, p. 68, adds expected downside risk from Eppen, Martin and Schrage:

u(\xi(\omega)) \ge g - \mathbf{q}^{\mathsf T}(\omega)\mathbf{y}(\omega), \qquad u \ge 0, \qquad \mathbb E_\xi\,u(\xi) \le l. \tag{5.1--5.3}

(2)

That is the CVaR epigraph with the $1/(1-\alpha)$ scaling removed and the threshold fixed at a target $g$ instead of optimized over $\nu$ — which is the clearest way to see what $\nu$ is doing in the CVaR form.

⚠ Neither section contains a worked numerical example. The book supplies the definitions; the numbers below are this course’s, computed on those definitions.

# Four plans, seven measures. Every plan is graded with recourse re-optimized.
def risk_profile(plan, target=100_000.0):
    _, p = grade_plan(plan, THREE)
    p = np.sort(p)
    out = {
        "E": p.mean(),
        "worst": p[0],
        "std dev": p.std(ddof=1),
        "MAD": np.abs(p - p.mean()).mean(),
    }
    for alpha in (2 / 3, 0.9):
        k = max(1, int(np.floor((1 - alpha) * len(p))))
        out[f"CVaR_{alpha:.2f}"] = p[:k].mean()
    out[f"EDR_{int(target / 1000)}k"] = np.maximum(target - p, 0).mean()
    return out


plans = [
    ("risk neutral (RP)", rp_plan),
    ("floor $56,000", [p for f, p, _ in frontier if f == 56_000][0]),
    ("Table 4 / worst case", table_4_plan),
    ("robust minimax (9.15)", plan_of(mm)),
]
header = None
for label, plan in plans:
    row = risk_profile(plan)
    if header is None:
        header = list(row)
        print(f"{'plan':>22s}" + "".join(f"{h:>12s}" for h in header))
    print(f"{label:>22s}" + "".join(f"{row[h]:12,.0f}" for h in header))

print("\nThe robust model returns exactly the Table 4 plan, which is what Birge and")
print("Louveaux assert on p. 19 without solving it.")

                  plan           E       worst     std dev         MAD   CVaR_0.67   CVaR_0.90    EDR_100k
     risk neutral (RP)     108,390      48,820      59,096      39,713      48,820      48,820      17,060

         floor $56,000     107,246      56,000      47,364      34,164      56,000      56,000      14,667

  Table 4 / worst case      86,600      59,950      26,650      17,767      59,950      59,950      17,817

 robust minimax (9.15)      86,600      59,950      26,650      17,767      59,950      59,950      17,817

The robust model returns exactly the Table 4 plan, which is what Birge and
Louveaux assert on p. 19 without solving it.

Read the last two columns against each other.

The worst-case plan has the best worst case, the best standard deviation, the best mean absolute deviation and the best CVaR at every confidence level — and the worst expected downside risk against a $100,000 target of all four plans, worse even than the plan that ignores risk entirely. It gives up so much upside that it falls short of $100,000 on average more often than the risk-neutral plan does.

“Less risky” is not a total order. Which plan is safest depends on which question you asked, and four numbers on the running example make that concrete in a way an axiom list does not.

# The expected-downside-risk frontier, B&L (5.1)-(5.3), p. 68, at g = $100,000.
TARGET = 100_000.0

# What is the smallest achievable expected downside risk? Ask, rather than guess:
# a risk CONSTRAINT can be infeasible, where a risk OBJECTIVE never is.
m = build_farmer(THREE)
m.u = pyo.Var(m.SCENARIOS, domain=pyo.NonNegativeReals, units=u.USD)
m.u_def = pyo.Constraint(
    m.SCENARIOS, rule=lambda m, i: m.u[i] >= TARGET * u.USD - m.scenarios[i].profit
)
m.objective.deactivate()
m.min_edr = pyo.Objective(expr=sum(m.u[i] for i in m.SCENARIOS) / 3, sense=pyo.minimize)
solve_farmer(m)
edr_floor = pyo.value(m.min_edr)
print(f"smallest attainable expected downside risk = ${edr_floor:,.0f}")
print(f"any limit below that is infeasible\n")

print(
    f"{'limit l':>10s} {'wheat':>8s} {'corn':>8s} {'beets':>8s}"
    f" {'E[profit]':>12s} {'achieved EDR':>14s} {'worst':>10s}"
)
for level in [20_000, 17_000, 15_000, 14_000]:
    m = solve_farmer(build_farmer(THREE, risk={"downside": (TARGET, level)}))
    p = plan_of(m)
    e, scenario_profits = grade_plan(p, THREE)
    achieved = np.maximum(TARGET - scenario_profits, 0).mean()
    print(
        f"{level:10,d} {p[0]:8.1f} {p[1]:8.1f} {p[2]:8.1f} {e:12,.0f}"
        f" {achieved:14,.0f} {scenario_profits.min():10,.0f}"
    )

smallest attainable expected downside risk = $13,758
any limit below that is infeasible

   limit l    wheat     corn    beets    E[profit]   achieved EDR      worst

    20,000    170.0     80.0    250.0      108,390         17,060     48,820

    17,000    167.9     82.1    250.0      108,375         17,000     49,000

    15,000    114.3    100.0    285.7      107,429         15,000     55,000

    14,000    100.0     71.4    328.6      101,176         14,000     58,000

Uncertain prices, not just uncertain yields¶

Everything so far has held prices fixed and let yields vary. Birge and Louveaux flag the other direction in the body text on p. 6 — “The influence of price fluctuations, illustrated by the dramatic price increases in 2007, is discussed in Exercise 8” — and then pose it as Exercise 8, “Data fluctuations”, p. 19, with its own data table:

	Wheat	Corn	Sugar Beets
Yield (T/acre)	2.5	3	20
Planting cost ($/acre)	180	280	310
Selling price ($/T)	300	170	41 under 6000 T, 11 above

“Consider a model with a random selling price of wheat being 300 or 220 with equal probability. Purchase prices are as before 40% higher than selling prices. … How much would a farmer be willing to pay for a perfect forecast on the selling price of wheat?”

Nothing in the builder needs to change: a scenario already carries its own prices alongside its own yields.

# Exercise 8 data, p. 19. Purchase prices are 40% above selling prices.
EX8_PLANT_COST = {"WHEAT": 180.0, "CORN": 280.0, "BEETS": 310.0}
EX8_SELL_PRICE = {
    "WHEAT": 300.0,
    "CORN": 170.0,
    "BEETS_FAVORABLE": 41.0,
    "BEETS_UNFAVORABLE": 11.0,
}


def price_scenario(wheat_price, yield_multiplier=1.0):
    """One scenario of Exercise 8: a wheat selling price, optionally with a yield."""
    sell = dict(EX8_SELL_PRICE, WHEAT=wheat_price)
    return {
        "yield": {c: MEAN_YIELD[c] * yield_multiplier for c in CROPS},
        "sell": sell,
        "buy": {c: 1.4 * sell[c] for c in FEED_CROPS},
    }


# Wait-and-see: what the farmer would do knowing next year's wheat price.
wait_and_see = {}
for price in (300.0, 220.0):
    m = solve_farmer(build_farmer([price_scenario(price)], plant_cost=EX8_PLANT_COST))
    wait_and_see[price] = pyo.value(m.expected_profit)
    p = plan_of(m) + 0.0  # normalize any -0.0 the LP returns
    print(
        f"perfect forecast, wheat at ${price:.0f}/T:"
        f"  acres {p[0]:5.1f} / {p[1]:5.1f} / {p[2]:5.1f}"
        f"   profit ${wait_and_see[price]:,.0f}"
    )

# Here-and-now: one plan, both prices, recourse after the price is revealed.
price_sp = solve_farmer(
    build_farmer(
        [price_scenario(300.0), price_scenario(220.0)], plant_cost=EX8_PLANT_COST
    )
)
price_plan = plan_of(price_sp)
rp_price = pyo.value(price_sp.expected_profit)
ws_price = np.mean(list(wait_and_see.values()))
print(
    f"\nstochastic plan: acres {price_plan[0]:5.1f} / {price_plan[1]:5.1f}"
    f" / {price_plan[2]:5.1f}   E[profit] ${rp_price:,.0f}"
)
print(f"WS   = ${ws_price:,.0f}")
print(f"RP   = ${rp_price:,.0f}")
print(f"EVPI = ${ws_price - rp_price:,.0f}   <- Exercise 8's question")

# EEV: plan at the mean price, then live with both.
mean_price_model = solve_farmer(
    build_farmer([price_scenario(260.0)], plant_cost=EX8_PLANT_COST)
)
eev, _ = grade_plan(
    plan_of(mean_price_model),
    [price_scenario(300.0), price_scenario(220.0)],
    plant_cost=EX8_PLANT_COST,
)
print(f"EEV  = ${eev:,.0f}")
print(f"VSS  = ${rp_price - eev:,.0f}")

perfect forecast, wheat at $300/T:  acres 500.0 /   0.0 /   0.0   profit $167,880

perfect forecast, wheat at $220/T:  acres 120.0 /  80.0 / 300.0   profit $131,000


stochastic plan: acres 200.0 /   0.0 / 300.0   E[profit] $137,880
WS   = $149,440
RP   = $137,880
EVPI = $11,560   <- Exercise 8's question

EEV  = $137,880
VSS  = $0

EVPI is $11,560 and VSS is $0 — at the same time¶

Lecture 10 makes the point that EVPI and VSS measure different things and do not order each other. Exercise 8 is a worked instance on the running example: a perfect price forecast is worth $11,560 a year, and the stochastic model itself is worth nothing over simply plugging in the mean price.

The reason is checkable, which is what makes it teachable.

# Why VSS = 0 here: does the recourse actually move with the price?
for i in price_sp.SCENARIOS:
    b = price_sp.scenarios[i]
    sales = {s: pyo.value(b.sales[s]) for s in SALE_TYPES}
    buys = {c: pyo.value(b.purchases[c]) for c in FEED_CROPS}
    label = "wheat $300" if i == 0 else "wheat $220"
    print(
        f"{label}:  sales "
        + ", ".join(f"{k} {v:,.0f}" for k, v in sales.items())
        + "   purchases "
        + ", ".join(f"{k} {v:,.0f}" for k, v in buys.items())
    )

identical = all(
    abs(
        pyo.value(price_sp.scenarios[0].sales[s])
        - pyo.value(price_sp.scenarios[1].sales[s])
    )
    < 1e-6
    for s in SALE_TYPES
) and all(
    abs(
        pyo.value(price_sp.scenarios[0].purchases[c])
        - pyo.value(price_sp.scenarios[1].purchases[c])
    )
    < 1e-6
    for c in FEED_CROPS
)
print(f"\nrecourse identical in both price scenarios: {identical}")

wheat $300:  sales WHEAT 300, CORN 0, BEETS_FAVORABLE 6,000, BEETS_UNFAVORABLE 0   purchases WHEAT 0, CORN 240
wheat $220:  sales WHEAT 300, CORN 0, BEETS_FAVORABLE 6,000, BEETS_UNFAVORABLE 0   purchases WHEAT 0, CORN 240

recourse identical in both price scenarios: True

The second-stage decision does not move. Both prices leave the farmer a wheat seller, so no decision flips, and $\mathbf{q}(\xi)^{\mathsf T}\mathbf{y}$ is then affine in the price. An affine function commutes with the expectation, so $\mathbb E_\xi\,Q(\mathbf x,\xi) = Q(\mathbf x, \mathbb E\,\xi)$ and the mean-price model is exact — hence $VSS = 0$ .

Birge and Louveaux ask for exactly this argument in Exercise 9, p. 20: show that for the news vendor with independent random prices and demands the solution is the one obtained with $\mathbf{q}$ and $r$ replaced by their expected values, and “indicate under which conditions the same proposition is true for the farmer’s problem.” The condition is the one just demonstrated — the recourse must not change with the price. It is a property of this instance, not a theorem: a wheat price low enough to make the farmer a buyer would break it.

⚠ Note what did not happen: making prices random did not make the problem nonlinear.

# Prices AND yields: two prices x three yields = six scenarios.
joint = [
    price_scenario(price, factor)
    for price in (300.0, 220.0)
    for factor in (1.2, 1.0, 0.8)
]
joint_model = solve_farmer(build_farmer(joint, plant_cost=EX8_PLANT_COST))
p = plan_of(joint_model)
print(
    f"acres {p[0]:.1f} / {p[1]:.1f} / {p[2]:.1f}"
    f"   E[profit] ${pyo.value(joint_model.expected_profit):,.0f}"
)
print(
    "scenario profits  "
    + " | ".join(
        f"${pyo.value(joint_model.scenarios[i].profit):,.0f}"
        for i in joint_model.SCENARIOS
    )
)
print(
    f"\nobjective polynomial degree = "
    f"{joint_model.objective.expr.polynomial_degree()}   (1 means linear)"
)

acres 250.0 / 0.0 / 250.0   E[profit] $135,880
scenario profits  $231,380 | $152,880 | $74,380 | $187,380 | $118,880 | $50,380

objective polynomial degree = 1   (1 means linear)

Each scenario block owns both its price data and its sales variables, so the revenue term $\mathbf{q}_s^{\mathsf T}\mathbf y_s$ is a constant times a variable — data times decision, never decision times decision. The extensive form stays a linear program.

The real costs of adding price uncertainty are different ones, and both are worth naming:

the scenario count multiplies, here $2\times3 = 6$ ; and
the crop-wise separability of §1.1d is lost, so the closed-form $Q_i$ used earlier in this notebook no longer supply a solver-free reference answer.

Modeling qualifications¶

A continuous random vector has uncountably many outcomes but can have finite dimension. A probability density can exceed one; its integral is one.
These historical farmer data are Birge and Louveaux (2011), 2nd ed., §1.1a–c, pp. 4–11. They are not current crop prices or policy.
The acreage heuristic applies to this data set. Higher purchase prices alone do not imply that buying feed is always suboptimal; the low-yield optimum buys corn.
Information values use §§4.1–4.3, pp. 163–166. Stage decisions are distinguished by information timing, not expense or irreversibility.