mixle.stochastic_opt module

Scenario-based CVaR selection under uncertainty (work-plan §7-H, H4).

Bridges a calibrated Posterior over an uncertain per-item value (IC-1, .samples(n, rng)) into a single-period binary selection decision. Rather than optimizing against one point-estimate value, two_stage_stochastic_plan draws k_scenarios realizations and solves the sample-average program: one shared selection decision x_b in {0, 1} per item, with per-scenario realized value v_k(x) = sum_b x_b * (price * g[k, b] - item_cost[b]). The objective trades expected value against downside risk via CVaR_alpha(-v_k(x)) (Rockafellar–Uryasev), so an item whose average realization looks profitable but whose scenario-conditional downside is large is priced correctly instead of naively included on its mean alone. CVaR-penalized scenario selection is a general finance/operations-research technique – portfolio selection under uncertain returns and project selection under uncertain payoffs are the same construction; this module’s worked instantiation is mine-planning block selection under grade uncertainty (an item is a block, price is ore price, and the downside risk being priced is grade uncertainty / ore-waste misclassification risk).

What “two-stage” does and does not mean here. The implemented program has exactly one decision vector, x, taken before any scenario is observed, plus the eta/u_k variables that linearize CVaR. There is no scenario-indexed operational decision, no recourse constraint or recourse cost, no second information stage, and therefore no nonanticipativity coupling to enforce – a scenario’s “recourse value” above is simply the payoff of that same x evaluated under scenario k. This is a single-stage sample-average risk-adjusted selection model (the deterministic equivalent of a chance-constrained/CVaR selection), not a two-stage stochastic program with recourse. The function name is kept for API compatibility with existing callers; a genuine second stage would require per-scenario decision variables y_k with their own constraints and costs, which nothing here builds. Read the outputs accordingly: expected_value and cvar describe the one committed selection, not the value of an adaptive policy.

cvar_epigraph is the reusable LP-representable epigraph of CVaR: given losses[k] = L_k(x) as a linear map of the (as yet undetermined) decision vector — an (K, n) coefficient matrix, not realized numbers — it emits the extra eta (Value-at-Risk) / u_k (excess-loss) variable block plus the a_ub rows enforcing u_k >= L_k(x) - eta, u_k >= 0, ready to be concatenated onto any existing MILP built on mixle.relations.branch_and_bound_milp().

Repo-boundary note (see the PR body for the full explanation): H1 (min_cost_flow et al., IC-9) and H3 (mixle.mine_planning, since renamed mixle.precedence_scheduling) had not landed on release/0.8.0 as of this PR, so neither this task’s frozen Public API nor its Algorithm section was written to call into either module directly — the scenario program here is built entirely on the already-landed mixle.reason.posterior_protocol (IC-1) and mixle.relations.branch_and_bound_milp(). Both H1 and H3 have since landed, but this module still imports neither mixle.relations’ flow surface nor mixle.precedence_scheduling: its scenario program never needed them, so the boundary held on its own merits rather than by landing-order accident.

J6 (the grand synthesis, work-plan §7-J) extends this scenario program with risk_adjusted_plan(): the same grade-uncertainty CVaR objective, but net of the priced environmental/health/carbon liabilities mixle.analysis.objective.priced_liabilities() assembles, and subject to the hard no-mine/exposure/ water constraints mixle.analysis.objective.hard_constraints() assembles — grade, cost, carbon, and enviro/health terms all trading off on one objective, one MILP. It is a distinct symbol from two_stage_stochastic_plan() above (never edited by J6), so the two land in different waves without a same-file write conflict.

class StochasticPlan(extract, expected_value, cvar, scenarios)[source]

Bases: NamedTuple

A scenario-optimal selection plan: which items, and its risk profile.

extract is the boolean per-item decision (named for this module’s mine-planning instantiation; kept as-is rather than renamed to something like selected since existing callers already depend on this field name); expected_value is E_k[v_k(extract)] over the scenarios the plan was optimized against; cvar is CVaR_alpha(-v_k(extract)) (a more negative value is safer — the tail is still profitable; a less negative or positive value is riskier); scenarios is the raw (K, n_items) value draws used for planning.

Both figures are recomputed from extract itself, not read out of the solver’s relaxation, so they always describe the selection actually being returned.

Parameters:
extract: ndarray

Alias for field number 0

expected_value: float

Alias for field number 1

cvar: float

Alias for field number 2

scenarios: ndarray

Alias for field number 3

cvar_epigraph(losses, alpha)[source]

Rockafellar–Uryasev LP epigraph of CVaR_alpha for a loss that is linear in the decision.

losses is a (K, n) matrix such that scenario k’s loss is losses[k] @ x for the not-yet-fixed length-n decision vector x — e.g. losses[k, b] = -(price * g[k, b] - cost[b]), the negated per-scenario recourse value. Returns the pieces of:

CVaR_alpha(L(x)) = min_{eta, u >= 0}  eta + (1 / ((1 - alpha) * K)) * sum_k u_k
                    s.t.  u_k >= losses[k] @ x - eta

embeddable alongside any existing constraints/variables on x:

  • c_add: length n + 1 + K objective row over [x, eta, u] giving the CVaR value itself (eta’s coefficient is 1, each u_k’s coefficient is 1 / ((1 - alpha) * K), zero on x) — combine with an expected-value objective as c_ev_padded - lam * c_add for a sense="max" solve of E[v] - lam * CVaR_alpha(-v).

  • a_ub_rows: (K, n + 1 + K) rows encoding losses[k] @ x - eta - u_k <= 0.

  • b_ub: length-K zeros, the right-hand side of a_ub_rows.

  • var_index: the column index of eta in the [x, eta, u] layout (= n); u_k sits at var_index + 1 + k. Give eta bounds (-inf, inf) and each u_k bounds (0, inf).

Parameters:
Return type:

tuple[ndarray, ndarray, ndarray, int]

risk_adjusted_plan(posterior, cost, price, liabilities, constraints, *, k_scenarios=50, alpha=0.9, rng)[source]

J6 — the risk-adjusted-NPV extension of two_stage_stochastic_plan().

Same two-stage scenario program (CVaR objective under scenario uncertainty), but:

  1. Per-scenario recourse value is net of liabilities["total"] — a length-n_items per-item dollar array (typically mixle.analysis.objective.priced_liabilities()’s output) folding in remediation, health, and carbon-price terms: v_k(x) = sum_b x_b * (price * g[k, b] - cost[b] - liabilities["total"][b]). An empty/absent liabilities ({} or no "total" key) means zero liability, i.e. identical to two_stage_stochastic_plan().

  2. constraints (typically mixle.analysis.objective.hard_constraints()’s output) adds hard constraints on top of the shared x_b in {0, 1} bounds:

    • "no_mine_mask": boolean array, True items are hard-fixed to x_b = 0 (this module’s mine-planning instantiation uses it for G9 no-mine/buffer zones; the mechanism itself is a general hard-exclusion mask) by tightening their variable bounds to (0, 0) — exact, not just penalized.

    • "caps": a list of {"coeffs": array, "bound": float} linear rows (already normalized to the solver’s <= convention), each added as an extra a_ub row coeffs @ x <= bound (K6 exposure budgets, L6 water budgets, or any other per-item activity cap).

Reuses cvar_epigraph() for the CVaR epigraph and mixle.relations.branch_and_bound_milp() for the extended MILP — the same solver two_stage_stochastic_plan uses, so value, cost, and any other priced terms all trade off against each other on one objective.

Parameters:
Return type:

StochasticPlan

two_stage_stochastic_plan(posterior, cost, price, *, k_scenarios=50, alpha=0.9, rng)[source]

Scenario program: select items to maximize E[v] - lambda * CVaR_alpha(-v).

Draws g = posterior.samples(k_scenarios, rng) (IC-1) as the calibrated value scenarios, forms the per-scenario realized value v_k(x) = sum_b x_b * (price * g[k, b] - cost[b]), and solves the joint MILP — binary x plus the free eta/u_k >= 0 block from cvar_epigraph() — via mixle.relations.branch_and_bound_milp().

Despite the name (kept for API compatibility), the program has a single here-and-now decision and no second-stage recourse variables; see the module docstring for exactly what is and is not modeled. Every scenario/economic input is validated as a finite, aligned schema first: a NaN price, an infinite cost, a (2, 3) cost matrix read as six items, or a fractional/Boolean scenario count is refused rather than turned into a different, plausible-looking problem.

Parameters:
Return type:

StochasticPlan