From 87bc7212b207c079118242341143ff02c859b4ed Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Thu, 9 Jul 2026 11:12:06 -0700 Subject: [PATCH 1/6] docs: design for Value of the Stochastic Solution (--vss report) Add doc/designs/vss_design.md: a post-run --vss report in generic_cylinders that computes VSS = EEV - RP. Reuses shipped machinery (average_scenario_creator for the mean-value solve; Xhat_Eval.evaluate for the honest cross-scenario EEV; wheel/EF for RP), mirroring the do_mmw post-run report structure. Covers the minimization/maximization sign convention, the EF-exact vs decomposition-bracketed RP, the infeasible-mean-value-solution case, and a prominent warning that EEV re-solves every scenario and can roughly double run time. Two-stage first; multistage plumbing (ef_xhat_nonants) noted as future. Closes the pysp_but_not_mpisppy.md A7 gap. Co-Authored-By: Claude Opus 4.8 --- doc/designs/vss_design.md | 437 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 437 insertions(+) create mode 100644 doc/designs/vss_design.md diff --git a/doc/designs/vss_design.md b/doc/designs/vss_design.md new file mode 100644 index 000000000..da60c3684 --- /dev/null +++ b/doc/designs/vss_design.md @@ -0,0 +1,437 @@ +# Value of the Stochastic Solution (VSS) design (two-stage first) + +Status: design approved by DLW (all §10 decisions resolved). +Implementation pending. + +Related: `doc/designs/pysp_but_not_mpisppy.md` §A7 (PySP `ef_vss.py` — +`create_expected_value_instance` + `fix_ef_first_stage_variables` — has no +mpi-sppy equivalent). This design closes that gap by *reusing* machinery +that already shipped, rather than porting PySP's code. + +--- + +## 0. What VSS is (and why anyone cares) + +The **Value of the Stochastic Solution** measures how much you gain by +actually modeling the uncertainty instead of collapsing it to a single +"average" future and solving that deterministic problem. + +The comparison is between two *first-stage* (here-and-now) decisions: + +- the decision you get from the full stochastic program, and +- the decision you get from the **mean-value problem** — replace every + random parameter by its expected value, solve the resulting + deterministic model, and take its first-stage solution. + +The mean-value decision is cheap and intuitive, and practitioners reach +for it constantly ("just plan for the average demand"). VSS answers the +question *"what does that shortcut cost me, in expectation, once the real +uncertainty shows up?"* A large VSS is the empirical justification for +having built a stochastic model at all; a near-zero VSS says the +deterministic shortcut was almost as good and the modeling effort bought +little. + +### 0.1 The three numbers (minimization convention) + +Write the two-stage stochastic program with first-stage vector `x`, +random data `ξ`, and second-stage value function `Q(x, ξ)`: + +``` +RP = min_x { c·x + E_ξ[ Q(x, ξ) ] } +``` + +- **RP** — *Recourse Problem* value: the optimal objective of the full + stochastic program. This is the here-and-now solution; **the run already + computes it** (see §3 on where it comes from and how exact it is). + +- **EV** — *Expected-Value (mean-value) problem*: replace `ξ` by its mean + `ξ̄` and solve the deterministic model + `EV = min_x { c·x + Q(x, ξ̄) }`. + Call its optimal first-stage solution `x̄`. Note `EV` is an *objective + value*; `x̄` is the *decision* we actually reuse. In the SAA world + mpi-sppy lives in, `ξ̄` is the average of the data over the explicit + scenario set the run uses — exactly what `average_scenario_creator` + already builds for Jensen's bound. + +- **EEV** — *Expected result of the EV solution*: pin the first stage at + `x̄`, then honestly evaluate expected cost across **every real + scenario**: + `EEV = c·x̄ + E_ξ[ Q(x̄, ξ) ]`. + +Then + +``` +VSS = EEV − RP (minimization; VSS ≥ 0 always) +``` + +`VSS ≥ 0` because `x̄` is feasible for the stochastic program but not +necessarily optimal, so evaluating it can only do worse than `RP`. + +For a **maximization** model every inequality flips and the definition is +`VSS = RP − EEV` (still `≥ 0`). The implementation reads +`is_minimizing` off the scenarios and picks the sign; VSS is always +reported as a non-negative "cost of using the average." + +### 0.2 Not to be confused with EVPI + +VSS is frequently confused with **EVPI** (Expected Value of Perfect +Information), `EVPI = RP − WS`, where `WS` ("wait-and-see") is the +expected value of solving each scenario *with perfect foresight*. EVPI +measures the value of *knowing the future*; VSS measures the value of +*modeling that you don't*. They are different numbers and answer +different questions. This design implements **VSS only**; EVPI is a +natural sibling (it reuses the same evaluation plumbing) and is noted as +future work in §8. + +--- + +## 1. Front-and-center warning: VSS can be expensive + +**This is a required, prominent note in the user-facing docs and in the +`--vss` help string.** + +Computing VSS is *not* free relative to the run that produced `RP`. The +three pieces cost very differently: + +| Quantity | Work | Typical cost | +| :--- | :--- | :--- | +| `RP` | already done by the run | free | +| `EV` / `x̄` | one deterministic solve of the average scenario | cheap | +| **`EEV`** | **fix `x̄`, solve the recourse subproblem for every scenario** | **can rival the original solve** | + +`EEV` is the expensive piece. It is a *full second pass over all +scenarios* — `N` second-stage solves (or one first-stage-fixed extensive +form). For a hard MIP with many scenarios, `--vss` can roughly **double** +the wall-clock of the job. On easy LPs it is negligible; on large +integer recourse it is not. The docs must say this plainly so nobody is +surprised when `--vss` turns a 10-minute run into a 20-minute run. + +A second, subtler cost: if `RP` came from a decomposition run that did +**not** close the optimality gap, then `RP` is only known to lie in a +bracket, and so is VSS (§3.2). Getting a tight point value for VSS may +require solving the extensive form for an exact `RP`, which can be far +more expensive than the decomposition run itself. `--vss` never does that +automatically; it reports the bracket instead. + +--- + +## 2. User-facing surface: the `--vss` flag + +The request calls it `-vss`; in mpi-sppy's Pyomo-`ConfigDict` CLI +convention that is the boolean flag **`--vss`** (single-dash `-vss` would +argparse-split into `-v -s -s`). Default `False`. + +Behavior: after the main algorithm finishes and `RP` is known, if +`cfg.vss` is set, `generic_cylinders` computes `EV`, `x̄`, `EEV`, and +`VSS` and prints a **VSS report** at the end of the run. It is +**report-only** — it never changes the solution, the bounds, or any +written solution file. It is an optional add-on to *either* an `--EF` +run or a decomposition (cylinders) run. + +Help string (drafted, includes the cost warning): + +``` +--vss After the run, report the Value of the Stochastic Solution + (VSS = EEV - RP). Requires the scenario module to define + average_scenario_creator. WARNING: computing EEV re-solves every + scenario with the first stage fixed and can roughly double run + time on large/integer models. Two-stage only (see docs). +``` + +### 2.1 Module contract and hard-fail + +`--vss` reuses the **same contract as Jensen's bound**: the scenario +module must define + +```python +def average_scenario_creator(scenario_name, **kwargs): + """Return a single deterministic scenario built from the sample-mean + of the random data (probability 1.0).""" +``` + +This already exists in `examples/farmer/farmer.py`, +`examples/sizes/sizes.py`, and `examples/netdes/netdes.py`, and is +discovered with the established `getattr(module, "average_scenario_creator", +None)` pattern. + +If `--vss` is set and the module has no `average_scenario_creator`, fail +**early and loudly at setup** (not after the whole run) with a message +that names the missing function and points at the Jensen's docs — the +same fail-fast posture used by the xhat-feasibility-cuts and Jensen's +paths. Wasting a long solve only to discover at the report step that VSS +cannot be computed is the failure mode to avoid. + +### 2.2 V1 scope restrictions (explicit, hard-error) + +To keep V1 honest and small, `--vss` refuses to combine with transforms +that change what "the objective" means, mirroring how `--cvar` already +refuses proper bundles / ADMM (`generic_cylinders.py`): + +- **Two-stage only.** Multistage is out of scope for V1 (§7). +- **No proper bundles, no ADMM, no CVaR** in V1. Each of these rewrites + the scenario objective or the scenario/first-stage structure, so `RP` + and `EEV` would have to be defined against the *transformed* problem to + be comparable. That is meaningful but subtle (what is "the average + scenario" of a risk measure?) and is deferred. Combining `--vss` with + any of them is a clear setup-time error in V1. + +Each restriction is a one-line guard with an explanatory message, not a +silent no-op. + +--- + +## 3. Where `RP` comes from, and how exact it is + +`RP` is the objective of the stochastic program, and its exactness +depends on which driver path produced it. + +### 3.1 `--EF` path — exact `RP` + +`do_EF` already solves the extensive form and has +`ef.get_objective_value()`. That *is* `RP`, exactly (up to the solver's +own gap). VSS from an EF run is a clean point value. + +### 3.2 Decomposition path — bracketed `RP` + +`do_decomp` returns the `WheelSpinner`, which exposes +`wheel.BestInnerBound` and `wheel.BestOuterBound`. For a minimization +run the true `RP` satisfies + +``` +BestOuterBound ≤ RP ≤ BestInnerBound +``` + +so, with `EEV` fixed, + +``` +VSS = EEV − RP ∈ [ EEV − BestInnerBound , EEV − BestOuterBound ]. +``` + +**Resolved (§10.2):** on a decomposition run the report therefore always +does both: + +- prints the **point value** `VSS = EEV − BestInnerBound` (using the + incumbent — the honest here-and-now value the run actually achieved), + clearly labeled as *conservative* (it underestimates the true VSS, + since the incumbent is `≥ RP`), and +- prints the **bracket** `[EEV − BestInnerBound, EEV − BestOuterBound]` + whenever the run's gap is nonzero, so the reader sees exactly how much + of the uncertainty in VSS is just the unclosed optimality gap. + +We report the incumbent-based point value rather than refusing one, +because the incumbent is the decision the run actually stands behind; the +always-printed bracket keeps that honest. The bracket falls out for free +from bounds the wheel already carries; no extra solve. The maximization +case swaps inner/outer roles accordingly. + +--- + +## 4. Where `EV`, `x̄`, and `EEV` come from — reuse, don't reinvent + +All three map onto functions that already exist and are tested. + +### 4.1 `EV` and `x̄` — the mean-value solve + +`mpisppy/utils/xhat_helpers.py::average_xhat_nonants` already: + +1. calls `average_scenario_creator(...)`, +2. asserts two-stage, +3. solves it, and +4. returns the ROOT first-stage values as a 1-D `np.ndarray`. + +That array **is** `x̄`. For VSS we additionally want the *objective* of +that solve to report `EV`, which `average_xhat_nonants` does not return. +**Resolved (§10.1):** `do_vss` builds and solves the average scenario +**inline** — a few lines that call `average_scenario_creator`, solve, and +read both `pyo.value()` (→ `EV`) and the root nonants (→ `x̄`). +This keeps `xhat_helpers.average_xhat_nonants` untouched; we do not add an +objective-returning sibling unless a second caller ever needs one. + +`EV` is reported for context (and to sanity-check `EEV ≥ EV`), but it is +**not** part of the VSS arithmetic — only `x̄` is. + +### 4.2 `EEV` — honest cross-scenario evaluation + +This is exactly what `Xhat_Eval.evaluate` does +(`mpisppy/utils/xhat_eval.py:257`): fix a nonant cache across all local +scenarios, `solve_loop`, and return the probability-weighted +`Eobjective`. Construct an `Xhat_Eval` over **all** scenarios (the +canonical pattern is in `confidence_intervals/ciutils.py:403` and +`zhat4xhat.py:94`), pack `x̄` as `{"ROOT": x̄}`, and call `evaluate`: + +```python +ev = Xhat_Eval(options, all_scenario_names, scenario_creator, + scenario_creator_kwargs=scenario_creator_kwargs, + all_nodenames=all_nodenames) +EEV = ev.evaluate({"ROOT": xbar}) # xbar is x̄ from §4.1 +``` + +`Xhat_Eval` already distributes scenarios across the MPI comm and reduces +the expected objective, so `EEV` is computed in parallel with no new +communication code. Crucially, `evaluate` uses the **same +`scenario_creator`** the run used, so `EEV` and `RP` are the same +objective measured two ways — apples to apples. + +### 4.3 Infeasibility of the mean-value solution is a *real answer* + +If `x̄` (built for the average scenario) is infeasible when fixed into +some real scenario — i.e. the model lacks relatively complete recourse — +then `Q(x̄, ξ) = +∞` for that scenario, `EEV = +∞`, and `VSS = +∞`. This +is not a bug; it is the strongest possible statement that the +deterministic shortcut is *unusable*. The report must therefore: + +- detect per-scenario infeasibility from the evaluate pass rather than + crashing, +- list which scenarios were infeasible, and +- print `EEV = +inf`, `VSS = +inf` with a one-line explanation. + +This is the same "fix a candidate, it may not be feasible everywhere" +concern documented for feasible-xhat (`doc/src/feasible_xhat.rst`); VSS +does not attempt any model-specific repair — a repaired `x̄` would no +longer be *the mean-value decision*, so repairing it would silently +change the quantity being measured. + +--- + +## 5. The report + +Printed once, on `cylinder_rank == 0` (per the rank-gating convention), +via `global_toc`, after the run's normal output. Sketch: + +``` +================= VSS report ================= + RP (stochastic solution, here-and-now) : 108900.0 [EF, exact] + EV (mean-value problem objective) : 107240.0 + EEV (EV first stage over all scenarios) : 115405.6 + VSS = EEV - RP : 6505.6 (5.98% of |RP|) +============================================= +``` + +Decomposition variant adds the bracket: + +``` + RP (stochastic solution, here-and-now) : 108900.0 [decomposition incumbent] + optimality bracket [outer, inner] : [108640.0, 108900.0] + ... + VSS = EEV - RP (point, conservative) : 6505.6 + VSS bracket [EEV-inner, EEV-outer] : [ 6505.6, 6765.6] +``` + +Infeasible variant: + +``` + EEV : +inf (EV solution infeasible in scenarios: Scenario7, Scenario12) + VSS : +inf (mean-value first stage is not usable across all scenarios) +``` + +The VSS percentage is `VSS / |RP|` guarded against `RP == 0`. + +--- + +## 6. Code layout + +Mirror the MMW post-run report exactly. + +- **New `mpisppy/generic/vss.py`** with `do_vss(module, cfg, ef=None, + wheel=None, scenario_creator=..., scenario_creator_kwargs=...)`, + paralleling `mpisppy/generic/mmw.py::do_mmw`. It sources `RP` from `ef` + or `wheel`, computes `EV`/`x̄`/`EEV`, and prints the report. A + `vss_requested(cfg)` predicate mirrors `mmw_requested`. +- **`mpisppy/generic_cylinders.py`**: after `do_EF` (pass the returned + `ef`) and after `do_decomp` (pass the returned `wheel`), call `do_vss` + when `cfg.get("vss")`. Two call sites, symmetric with the existing + `do_mmw` calls (`generic_cylinders.py:141,146`). +- **`mpisppy/utils/config.py`**: one `add_to_config('vss', ...)` boolean, + in whichever driver-level args group owns post-run report options + (alongside the MMW flags), so it is available to `generic_cylinders`. +- **`mpisppy/utils/xhat_helpers.py`**: untouched. `do_vss` solves the + average scenario inline to get both `EV` and `x̄` (§4.1); no + objective-returning sibling is added in V1. + +No spoke edits, no hub edits, no changes to the algorithms — VSS is +strictly a post-processing consumer of results the drivers already +return, exactly like MMW. + +--- + +## 7. Multistage (future, but the plumbing is already here) + +Two-stage is V1. The multistage generalization is well-defined and the +building blocks exist: + +- The mean-value analogue is the **expected-value tree** — the real + branching structure with the random data at each node pinned to its + conditional mean. Solving that EF yields a first-stage-through-last + decision at *every* non-leaf node. +- `mpisppy/utils/xhat_helpers.py::ef_xhat_nonants` already solves an EF + and returns the whole `{node_name: np.ndarray}` nonant tree — the exact + cache form `Xhat_Eval._fix_nonants` / `fix_nonants_upto_stage` + consume. So `EEV` for the multistage case is "fix the whole EV tree, + evaluate across real scenarios" and needs no new evaluation code. + +What is missing for multistage is a **module contract for building the +expected-value tree** (an `ev_tree` creator, or a documented convention +for driving `average_scenario_creator` per node). That contract, plus +the choice among the several textbook multistage-VSS variants, is why +multistage is deferred rather than shipped in V1. The `--vss` two-stage +guard (§2.2) is the placeholder that later lifts. + +--- + +## 8. Testing plan + +New `mpisppy/tests/test_vss.py` (and wire it into `run_coverage.bash` +**and** `test_pr_and_main.yml` in the same commit, per the coverage-harness +convention): + +1. **Farmer EF, exact RP.** Run `--EF --vss`. Independently, in the test: + solve the average scenario for `x̄`, `Xhat_Eval.evaluate({"ROOT": x̄})` + for `EEV`, and `ef.get_objective_value()` for `RP`. Assert the + reported VSS equals `EEV − RP` to tolerance, and that `VSS ≥ 0`. +2. **Farmer decomposition, bracket.** Run cylinders `--vss`; assert the + point VSS uses `BestInnerBound` and the printed bracket endpoints + match `EEV − BestInnerBound` / `EEV − BestOuterBound`. +3. **Maximization sign.** Sign-flipped farmer variant; assert + `VSS = RP − EEV ≥ 0` and that the report labels it correctly (reuse + the max-support test fixtures). +4. **Missing `average_scenario_creator`.** A module without it + `--vss` + ⇒ `RuntimeError` at setup, before any long solve. +5. **Infeasible mean-value solution.** A small model with no relatively + complete recourse where `x̄` is infeasible in ≥1 scenario ⇒ report + shows `EEV = +inf`, `VSS = +inf`, names the offending scenarios, and + does not crash. +6. **V1 restriction guards.** `--vss` with `--EF`... plus `--cvar` (or a + proper-bundle / ADMM config) ⇒ clear setup-time error. + +--- + +## 9. Documentation + +- New `doc/src/vss.rst`: the §0 "what VSS is" explanation (RP/EV/EEV/VSS, + the minimization convention, and the EVPI distinction), the §1 cost + warning **at the top**, a worked farmer example, and the decomposition + bracket + infeasibility semantics. Cross-link `doc/src/jensens.rst` + (shared `average_scenario_creator` contract) and + `doc/src/feasible_xhat.rst` (shared fix-a-candidate feasibility caveat). +- Add `--vss` to the generic_cylinders option reference. +- One line in `doc/designs/pysp_but_not_mpisppy.md` §A7 noting the gap is + addressed (two-stage) once this ships. + +--- + +## 10. Resolved decisions + +All resolved; implementation can proceed. + +1. **`EV` objective plumbing (§4.1/§6):** solve the average scenario + **inline** in `do_vss` and read both its objective (`EV`) and root + nonants (`x̄`). `xhat_helpers.average_xhat_nonants` stays untouched; no + objective-returning sibling until a second caller needs one. +2. **Decomposition point value (§3.2):** report the **incumbent-based + point value** `VSS = EEV − BestInnerBound`, labeled conservative, and + **always** print the `[EEV − BestInnerBound, EEV − BestOuterBound]` + bracket when the gap is open. Do not refuse a point value. +3. **EVPI (§0.2):** **VSS only** in V1, for a review-sized PR. EVPI is a + natural follow-on (it reuses the same evaluate plumbing, with WS = + per-scenario perfect-foresight solves) and is left to a later PR. +``` From a7e23a84817976be7c6e7b125de87af06a0b0170 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Thu, 9 Jul 2026 13:55:04 -0700 Subject: [PATCH 2/6] feat(vss): add --vss Value of the Stochastic Solution report (V1, two-stage) generic_cylinders --vss prints RP / EV / EEV / VSS after an --EF or a decomposition run. Reuses shipped machinery: average_scenario_creator for the mean-value solve (x_bar) and Xhat_Eval.evaluate for the honest cross-scenario EEV; RP is exact from the EF or the incumbent (with a BestOuterBound/BestInnerBound bracket) from the wheel. Structure mirrors do_mmw: new mpisppy/generic/vss.py (vss_prep + do_vss), one config flag, two call sites in generic_cylinders. - Sign-aware: VSS = EEV-RP (min) or RP-EEV (max), always reported >= 0. - Decomposition run reports a conservative point value plus the VSS bracket implied by the unclosed optimality gap. - Infeasible mean-value first stage => EEV = VSS = +inf, naming the offending scenarios (MPI-safe: solves with need_solution=False so a clean per-scenario infeasibility does not deadlock the collectives). - Two-stage only; fails fast at setup (before the solve) on multistage, proper bundles, ADMM, --cvar, or a missing average_scenario_creator. Tests: mpisppy/tests/test_vss.py (guards, sign/bracket helpers, EF min/max e2e, decomposition bracket, infeasible EEV), wired into run_coverage.bash and test_pr_and_main.yml. Docs: doc/src/vss.rst (+ index toctree); pysp_but_not_mpisppy.md A7 marked addressed for two-stage. Design in doc/designs/vss_design.md. Co-Authored-By: Claude Opus 4.8 --- .github/workflows/test_pr_and_main.yml | 1 + doc/designs/pysp_but_not_mpisppy.md | 3 + doc/designs/vss_design.md | 5 +- doc/src/index.rst | 1 + doc/src/vss.rst | 120 ++++++++++ mpisppy/generic/parsing.py | 1 + mpisppy/generic/vss.py | 316 +++++++++++++++++++++++++ mpisppy/generic_cylinders.py | 18 +- mpisppy/tests/test_vss.py | 200 ++++++++++++++++ mpisppy/utils/config.py | 12 + run_coverage.bash | 3 + 11 files changed, 676 insertions(+), 4 deletions(-) create mode 100644 doc/src/vss.rst create mode 100644 mpisppy/generic/vss.py create mode 100644 mpisppy/tests/test_vss.py diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index 7e90c591b..b01340052 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -1192,6 +1192,7 @@ jobs: mpisppy/tests/test_solver_spec.py \ mpisppy/tests/test_extensions.py \ mpisppy/tests/test_jensens.py \ + mpisppy/tests/test_vss.py \ mpisppy/tests/test_feasible_xhat.py \ mpisppy/tests/test_proper_bundler.py \ mpisppy/tests/test_incumbent_writing.py \ diff --git a/doc/designs/pysp_but_not_mpisppy.md b/doc/designs/pysp_but_not_mpisppy.md index a031cbe3e..70cc3cdd9 100644 --- a/doc/designs/pysp_but_not_mpisppy.md +++ b/doc/designs/pysp_but_not_mpisppy.md @@ -85,6 +85,9 @@ dual-bound extension itself is not reproduced. computes the VSS / evaluates the expected-value solution across scenarios. mpi-sppy has `average_scenario_creator` (an EV building block) but no dedicated VSS/EEV computation utility. +**Addressed (two-stage):** `generic_cylinders --vss` now reports RP/EV/EEV/VSS +after a run (`mpisppy/generic/vss.py`; see `doc/src/vss.rst` and +`doc/designs/vss_design.md`). Multistage VSS is still a gap. ## B. PH algorithmic knobs that are missing diff --git a/doc/designs/vss_design.md b/doc/designs/vss_design.md index da60c3684..a5d6dfabd 100644 --- a/doc/designs/vss_design.md +++ b/doc/designs/vss_design.md @@ -1,7 +1,8 @@ # Value of the Stochastic Solution (VSS) design (two-stage first) -Status: design approved by DLW (all §10 decisions resolved). -Implementation pending. +Status: design approved by DLW (all §10 decisions resolved). V1 implemented +(`mpisppy/generic/vss.py`, `--vss` flag, `mpisppy/tests/test_vss.py`, +`doc/src/vss.rst`). Related: `doc/designs/pysp_but_not_mpisppy.md` §A7 (PySP `ef_vss.py` — `create_expected_value_instance` + `fix_ef_first_stage_variables` — has no diff --git a/doc/src/index.rst b/doc/src/index.rst index 2f577ac23..6a4b63f95 100644 --- a/doc/src/index.rst +++ b/doc/src/index.rst @@ -58,6 +58,7 @@ MPI is used. properbundles.rst pickling.rst jensens.rst + vss.rst feasible_xhat.rst xhat_from_file.rst iis.rst diff --git a/doc/src/vss.rst b/doc/src/vss.rst new file mode 100644 index 000000000..1c0e58f91 --- /dev/null +++ b/doc/src/vss.rst @@ -0,0 +1,120 @@ +.. _vss: + +Value of the Stochastic Solution (VSS) +====================================== + +The **Value of the Stochastic Solution** measures how much is gained by +modeling the uncertainty instead of collapsing it to a single "average" +future and solving that deterministic problem. It compares two +*here-and-now* (first-stage) decisions: the one from the full stochastic +program, and the one from the **mean-value problem** (replace every random +parameter by its expected value and solve the resulting deterministic +model). VSS is the expected extra cost of using the mean-value decision +once the real uncertainty shows up. A large VSS justifies building a +stochastic model; a near-zero VSS says the deterministic shortcut was +almost as good. + +.. warning:: + + **Computing VSS can be expensive.** The mean-value solve is cheap, but + the ``EEV`` term (below) re-solves the recourse problem for *every* + scenario with the first stage fixed — a full second pass over all + scenarios. On large or integer models, ``--vss`` can roughly **double** + the run time. It is off by default. + +The three numbers (minimization) +-------------------------------- + +Write the two-stage program with first-stage vector :math:`x`, random data +:math:`\xi`, and second-stage value function :math:`Q(x,\xi)`: + +.. math:: + + RP = \min_x \; \big\{\, c\cdot x + \mathbb{E}_\xi[\,Q(x,\xi)\,] \,\big\}. + +- **RP** — the *Recourse Problem* optimum: the here-and-now solution the run + already computes. It is **exact** from an ``--EF`` run and the + **incumbent** (with a bracket) from a decomposition run. +- **EV** — the *mean-value problem* + :math:`\min_x \{ c\cdot x + Q(x,\bar\xi)\}` with :math:`\bar\xi` the + average of the data over the scenario set. Its first-stage solution is + :math:`\bar x`. +- **EEV** — the *Expected result of the EV solution*: pin the first stage at + :math:`\bar x` and evaluate honestly across every scenario, + :math:`EEV = c\cdot\bar x + \mathbb{E}_\xi[\,Q(\bar x,\xi)\,]`. + +Then + +.. math:: + + VSS = EEV - RP \;\; (\ge 0). + +For a **maximization** model the sign flips: :math:`VSS = RP - EEV` +(still :math:`\ge 0`). ``mpi-sppy`` reads the model sense and reports VSS as +a non-negative "cost of using the average." + +VSS is not EVPI (Expected Value of Perfect Information, +:math:`RP - WS`), which measures the value of *knowing* the future rather +than *modeling* that you do not. ``--vss`` reports VSS only. + +Usage +----- + +Add ``--vss`` to a ``generic_cylinders`` run. The scenario module must +define ``average_scenario_creator`` (the same contract used by +:ref:`jensens`); it builds the mean-value scenario whose first-stage +solution VSS evaluates. The report is printed at the end and does not change +the solution. + +.. code-block:: bash + + cd examples/farmer + python ../../mpisppy/generic_cylinders.py --module-name farmer \ + --num-scens 3 --EF --EF-solver-name gurobi --solver-name gurobi --vss + +produces:: + + ================= VSS report ================= + RP (stochastic solution, here-and-now): -108390 [EF, exact] + EV (mean-value problem objective) : -118600 + EEV (EV first stage over scenarios) : -107240 + VSS = EEV - RP (sense=min) : 1150 (1.06% of |RP|) + ============================================= + +Exact vs. bracketed RP +~~~~~~~~~~~~~~~~~~~~~~~~ + +From an ``--EF`` run, ``RP`` is exact. From a decomposition (cylinders) run +that did not close the optimality gap, ``RP`` is only known to lie between +the outer and inner bounds, so VSS is reported both as a conservative point +value (against the incumbent) and as a bracket:: + + RP (stochastic solution, here-and-now): -108389 [decomposition incumbent] + optimality bracket [outer, inner] : [-108508, -108389] + ... + VSS = EEV - RP (sense=min) : 1149.33 (1.06% of |RP|) + VSS bracket from RP bracket : [1149.33, 1268.39] + +The width of the VSS bracket is exactly the unclosed optimality gap; tighten +the run (or use ``--EF``) for a sharper VSS. + +Infeasible mean-value solution +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +If the mean-value first stage :math:`\bar x` is infeasible when fixed in +some scenario (the model lacks relatively complete recourse), then that +scenario's recourse cost is :math:`+\infty`, so ``EEV`` and ``VSS`` are +``+inf``. This is a meaningful result — the deterministic shortcut is +unusable — and the report names the offending scenarios. See +:ref:`feasible_xhat` for the related "fix a candidate, repair if needed" +concern; VSS deliberately does **not** repair :math:`\bar x`, since a +repaired point would no longer be the mean-value decision. + +Limitations +----------- + +This version is **two-stage only** and cannot be combined with proper +bundles, ADMM, or ``--cvar`` (each rewrites the objective or the +scenario/first-stage structure, so ``RP`` and ``EEV`` would no longer be +comparable). ``--vss`` fails fast at setup if the model or configuration is +unsupported, so a long run is never wasted. diff --git a/mpisppy/generic/parsing.py b/mpisppy/generic/parsing.py index b336fdc74..cd90ea941 100644 --- a/mpisppy/generic/parsing.py +++ b/mpisppy/generic/parsing.py @@ -187,6 +187,7 @@ def parse_args(m): # TBD - think about adding directory for json options files cfg.mmw_args() + cfg.vss_args() from mpisppy.generic.admm import admm_args admm_args(cfg) diff --git a/mpisppy/generic/vss.py b/mpisppy/generic/vss.py new file mode 100644 index 000000000..0b07db1aa --- /dev/null +++ b/mpisppy/generic/vss.py @@ -0,0 +1,316 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +"""Value of the Stochastic Solution (VSS) report for generic_cylinders. + +VSS = EEV - RP (minimization; VSS = RP - EEV for maximization), where + + RP = the stochastic-program optimum -- the here-and-now solution the run + already computed: exact from an --EF run, or the decomposition + incumbent (BestInnerBound), with a bracket from BestOuterBound. + EV = the mean-value (expected-value) problem: replace the random data by + its average and solve. x_bar is its first-stage solution. + EEV = x_bar fixed in the first stage and evaluated honestly across every + scenario. + +Two-stage only in this version. See doc/designs/vss_design.md. + +WARNING: computing EEV re-solves every scenario with the first stage fixed; +on large or integer models this can rival the cost of the original run. +""" + +import math + +import numpy as np +import pyomo.environ as pyo + +import mpisppy.utils.sputils as sputils +from mpisppy import global_toc, MPI +from mpisppy.utils import xhat_eval +from mpisppy.generic.parsing import name_lists, proper_bundles + + +def vss_prep(module, cfg): + """Validate that a VSS report can be produced, and fail fast -- before + the main solve -- if not. Called early by generic_cylinders when --vss + is set, so a long run is never wasted only to fail at the report step. + + V1 is two-stage only and does not (yet) support the objective-rewriting + transform (CVaR) or the scenario/first-stage-restructuring paths (proper + bundles, ADMM), because RP and EEV would then have to be defined against + the transformed problem to be comparable. + """ + if getattr(module, "average_scenario_creator", None) is None: + raise RuntimeError( + "--vss requires the scenario module to define " + "average_scenario_creator (the same contract used by Jensen's " + "bound; see doc/src/jensens.rst). It builds the mean-value " + "scenario whose first-stage solution VSS evaluates." + ) + if cfg.get("branching_factors") is not None: + raise RuntimeError( + "--vss is two-stage only in this version; it cannot be used with " + "multistage runs (--branching-factors)." + ) + if proper_bundles(cfg): + raise RuntimeError("--vss cannot (yet) be combined with proper bundles.") + if cfg.get("admm", ifmissing=False) or cfg.get("stoch_admm", ifmissing=False): + raise RuntimeError("--vss cannot (yet) be combined with ADMM.") + if cfg.get("cvar", ifmissing=False): + raise RuntimeError("--vss cannot (yet) be combined with --cvar.") + + +def do_vss(module, cfg, scenario_creator, scenario_creator_kwargs, + scenario_denouement, ef=None, wheel=None): + """Compute and print a VSS report after the main algorithm. + + Exactly one of ``ef`` (from an --EF run) or ``wheel`` (from a + decomposition run) supplies RP. MPI-collective: EEV is evaluated across + all ranks, so every rank must call this together. + + Args: + module: the model module (must define average_scenario_creator) + cfg (Config): parsed options + scenario_creator (function): the run's scenario creator + scenario_creator_kwargs (dict): kwargs for the scenario creator + scenario_denouement (function): unused for the quiet EEV pass + ef (ExtensiveForm or None): source of an exact RP + wheel (WheelSpinner or None): source of an incumbent RP + bracket + """ + if (ef is None) == (wheel is None): + raise ValueError("do_vss: pass exactly one of ef= or wheel=") + comm = MPI.COMM_WORLD + + all_scenario_names, _ = name_lists(module, cfg) + + # EV problem: solve the mean-value scenario for its objective and its + # first-stage solution x_bar. + ev_obj, x_bar, is_min = _solve_average_scenario( + module, cfg, scenario_creator_kwargs) + + # EEV: fix x_bar in the first stage, evaluate across all scenarios. + eev, infeasible_names = _compute_eev( + cfg, scenario_creator, scenario_creator_kwargs, + all_scenario_names, x_bar) + + # RP: exact from EF, or incumbent (+ bracket) from the wheel. + if ef is not None: + rp_point = comm.bcast(ef.get_objective_value() if comm.Get_rank() == 0 + else None, root=0) + rp_source = "EF, exact" + inner = outer = None + else: + inner, outer = _reduce_bounds(comm, wheel.BestInnerBound, + wheel.BestOuterBound, is_min) + rp_point = inner + rp_source = "decomposition incumbent" + + vss = _vss_value(is_min, rp_point, eev) + + # Bracket: with RP only known to lie in [outer, inner], VSS lies in a + # corresponding interval. Only meaningful when the run left a gap and + # EEV is finite. + have_bracket = (inner is not None and outer is not None + and math.isfinite(inner) and math.isfinite(outer) + and inner != outer) + vss_bracket = None + if have_bracket and math.isfinite(eev): + if is_min: + vss_bracket = (eev - inner, eev - outer) + else: + vss_bracket = (inner - eev, outer - eev) + + result = {"RP": rp_point, "EV": ev_obj, "EEV": eev, "VSS": vss, + "is_min": is_min, "rp_source": rp_source, + "inner": inner, "outer": outer, "vss_bracket": vss_bracket, + "infeasible_scenarios": infeasible_names} + _print_report(result) + return result + + +def _solve_average_scenario(module, cfg, scenario_creator_kwargs): + """Build and solve the mean-value scenario. Return + (EV_objective, x_bar, is_minimizing) where x_bar is the ROOT first-stage + solution as a 1-D np.ndarray in nonant_vardata_list order. + + Two-stage only (asserted). Every rank solves this identical deterministic + scenario independently -- no collective communication. + """ + avg = module.average_scenario_creator( + "AverageScenario", **(scenario_creator_kwargs or {})) + if not hasattr(avg, "_mpisppy_node_list") or len(avg._mpisppy_node_list) != 1: + raise RuntimeError( + "--vss is two-stage only; average_scenario_creator must return a " + "model with exactly one tree node (ROOT)." + ) + solver = pyo.SolverFactory(cfg.solver_name) + if sputils.is_persistent(solver): + solver.set_instance(avg) + results = solver.solve(tee=False) + else: + results = solver.solve(avg, tee=False) + if not pyo.check_optimal_termination(results): + raise RuntimeError( + "--vss: the mean-value (EV) problem did not solve to optimality " + f"(termination_condition={results.solver.termination_condition}). " + "VSS cannot be computed." + ) + obj = sputils.find_active_objective(avg) + ev_obj = pyo.value(obj) + is_min = (obj.sense == pyo.minimize) + root = avg._mpisppy_node_list[0] + x_bar = np.array([pyo.value(v) for v in root.nonant_vardata_list], dtype="d") + return ev_obj, x_bar, is_min + + +def _compute_eev(cfg, scenario_creator, scenario_creator_kwargs, + all_scenario_names, x_bar): + """Fix x_bar in the first stage and evaluate expected cost across all + scenarios. Return (EEV, infeasible_scenario_names). EEV is math.inf if + the mean-value first stage is infeasible in any scenario. + + MPI-collective: every rank builds its share of scenarios, so all ranks + must call this together. + """ + options = { + "iter0_solver_options": None, + "iterk_solver_options": None, + "display_timing": False, + "solver_name": cfg.solver_name, + "verbose": False, + "solver_options": None, + } + ev = xhat_eval.Xhat_Eval( + options, + all_scenario_names, + scenario_creator, + scenario_denouement=None, + scenario_creator_kwargs=scenario_creator_kwargs, + all_nodenames=None, + mpicomm=MPI.COMM_WORLD, + ) + ev._lazy_create_solvers() + ev._fix_nonants({"ROOT": x_bar}) + # compute_val_at_nonant=False => need_solution=False in solve_one, so a + # clean per-scenario infeasibility does NOT raise. That matters for MPI: + # a raise on only some ranks would deadlock the collectives below. + ev.solve_loop(solver_options=None, gripe=True, tee=False, + compute_val_at_nonant=False) + + local_infeasible = [ + k for k, s in ev.local_scenarios.items() + if not getattr(s._mpisppy_data, "solution_available", False) + ] + n_local = np.array([float(len(local_infeasible))]) + n_global = np.zeros(1) + ev.mpicomm.Allreduce(n_local, n_global, op=MPI.SUM) + + if n_global[0] > 0: + gathered = ev.mpicomm.gather(local_infeasible, root=0) + names = [] + if ev.mpicomm.Get_rank() == 0: + for lst in gathered: + names.extend(lst) + return math.inf, names + + # All feasible: Eobjective is collective, so every rank calls it. + return ev.Eobjective(), [] + + +def _reduce_bounds(comm, best_inner, best_outer, is_min): + """Make the wheel's best inner/outer bounds available on every rank. + + BestInnerBound / BestOuterBound are set only on the hub rank(s) and are + None elsewhere. Reduce with the sense-appropriate op, treating a missing + value as the identity (so the hub's real value wins). + """ + if is_min: + inner = comm.allreduce(best_inner if best_inner is not None + else math.inf, op=MPI.MIN) + outer = comm.allreduce(best_outer if best_outer is not None + else -math.inf, op=MPI.MAX) + else: + inner = comm.allreduce(best_inner if best_inner is not None + else -math.inf, op=MPI.MAX) + outer = comm.allreduce(best_outer if best_outer is not None + else math.inf, op=MPI.MIN) + return inner, outer + + +def _vss_value(is_min, rp, eev): + """VSS with the sign convention; math.inf if EEV is infinite.""" + if math.isinf(eev): + return math.inf + return (eev - rp) if is_min else (rp - eev) + + +def _print_report(result): + """Print the VSS report once, on global rank 0, from the do_vss dict.""" + if MPI.COMM_WORLD.Get_rank() != 0: + return + + def _f(x): + if x is None: + return "n/a" + if math.isinf(x): + return "+inf" if x > 0 else "-inf" + return f"{x:.6g}" + + def _row(label, value, suffix=""): + return f" {label:<39}: {value:>14}{suffix}" + + is_min = result["is_min"] + eev = result["EEV"] + inner, outer = result["inner"], result["outer"] + sense = "min" if is_min else "max" + formula = "EEV - RP" if is_min else "RP - EEV" + + lines = ["", "================= VSS report =================", + _row("RP (stochastic solution, here-and-now)", + _f(result["RP"]), f" [{result['rp_source']}]")] + + have_bracket = (inner is not None and outer is not None + and math.isfinite(inner) and math.isfinite(outer) + and inner != outer) + if have_bracket: + lo, hi = (outer, inner) if is_min else (inner, outer) + lines.append(_row(" optimality bracket [outer, inner]", + f"[{_f(lo)}, {_f(hi)}]")) + + lines.append(_row("EV (mean-value problem objective)", _f(result["EV"]))) + + if math.isinf(eev): + names = result["infeasible_scenarios"] + shown = names[:8] + more = "" if len(names) <= 8 else f" (+{len(names) - 8} more)" + lines.append(_row("EEV (EV first stage over scenarios)", "+inf", + f" infeasible in: {', '.join(shown)}{more}")) + lines.append(_row(f"VSS = {formula} (sense={sense})", "+inf", + " (mean-value first stage not usable everywhere)")) + lines.append("=============================================") + print("\n".join(lines)) + global_toc("VSS = +inf (mean-value first stage infeasible " + "in some scenarios)") + return + + lines.append(_row("EEV (EV first stage over scenarios)", _f(eev))) + + vss, rp = result["VSS"], result["RP"] + pct = "" + if rp != 0 and math.isfinite(rp): + pct = f" ({100.0 * vss / abs(rp):.2f}% of |RP|)" + lines.append(_row(f"VSS = {formula} (sense={sense})", _f(vss), pct)) + + if result["vss_bracket"] is not None: + vlo, vhi = result["vss_bracket"] + lines.append(_row(" VSS bracket from RP bracket", + f"[{_f(vlo)}, {_f(vhi)}]")) + + lines.append("=============================================") + print("\n".join(lines)) + global_toc(f"VSS = {_f(vss)} ({formula})") diff --git a/mpisppy/generic_cylinders.py b/mpisppy/generic_cylinders.py index d8d8ab832..dc5100883 100644 --- a/mpisppy/generic_cylinders.py +++ b/mpisppy/generic_cylinders.py @@ -45,6 +45,12 @@ if hasattr(module, "get_mpisppy_helper_object"): module = module.get_mpisppy_helper_object(cfg) + # Fail fast (before the main solve) if a --vss report was requested but + # cannot be produced for this model/config. + if cfg.get("vss", ifmissing=False): + from mpisppy.generic.vss import vss_prep + vss_prep(module, cfg) + bundle_wrapper = None # the default if proper_bundles(cfg): # Nonant name validation will fail with proper bundles because @@ -136,12 +142,20 @@ def scenario_denouement(rank, sname, s): bundle_wrapper=bundle_wrapper, ) elif cfg.EF: - do_EF(module, cfg, scenario_creator, scenario_creator_kwargs, - scenario_denouement, bundle_wrapper=bundle_wrapper) + ef = do_EF(module, cfg, scenario_creator, scenario_creator_kwargs, + scenario_denouement, bundle_wrapper=bundle_wrapper) + if cfg.get("vss", ifmissing=False): + from mpisppy.generic.vss import do_vss + do_vss(module, cfg, scenario_creator, scenario_creator_kwargs, + scenario_denouement, ef=ef) if mmw_requested(cfg) and cfg.get("mmw_xhat_input_file_name") is not None: do_mmw(fname, cfg) else: wheel = do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_denouement, bundle_wrapper=bundle_wrapper) + if cfg.get("vss", ifmissing=False): + from mpisppy.generic.vss import do_vss + do_vss(module, cfg, scenario_creator, scenario_creator_kwargs, + scenario_denouement, wheel=wheel) if mmw_requested(cfg): do_mmw(fname, cfg, wheel=wheel) diff --git a/mpisppy/tests/test_vss.py b/mpisppy/tests/test_vss.py new file mode 100644 index 000000000..f4aa57e17 --- /dev/null +++ b/mpisppy/tests/test_vss.py @@ -0,0 +1,200 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +"""Tests for the Value of the Stochastic Solution (--vss) report. + +Covers vss_prep's fail-fast guards, the sign/bracket helpers, and end-to-end +EF runs (minimize and maximize) plus the infeasible-mean-value-solution path, +all on the two-stage farmer example. +""" + +import math +import os +import sys +import types +import unittest + +import numpy as np +import pyomo.environ as pyo + +import mpisppy.opt.ef +import mpisppy.utils.sputils as sputils +import mpisppy.generic.parsing as parsing +from mpisppy.generic import vss +from mpisppy import MPI +from mpisppy.tests.utils import get_solver + +# make the farmer example importable (it defines average_scenario_creator) +_EXAMPLES_DIR = os.path.join( + os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))), + "examples") +sys.path.insert(0, os.path.join(_EXAMPLES_DIR, "farmer")) +import farmer # noqa: E402 + +solver_available, solver_name, _, _ = get_solver() + + +class _FakeCfg: + """Minimal cfg exposing only .get(key, ifmissing=...) for vss_prep.""" + def __init__(self, d): + self._d = d + + def get(self, key, ifmissing=None): + return self._d.get(key, ifmissing) + + +def _make_cfg(extra_argv): + """Build a real, fully-parsed Config for the farmer module.""" + argv = ["prog", "--module-name", "farmer", "--num-scens", "3", + "--solver-name", solver_name or "gurobi"] + extra_argv + old = sys.argv + sys.argv = argv + try: + return parsing.parse_args(farmer) + finally: + sys.argv = old + + +class TestVssPrep(unittest.TestCase): + """vss_prep should fail fast on unsupported configurations.""" + + def test_clean_two_stage_passes(self): + vss.vss_prep(farmer, _FakeCfg({})) # farmer has average_scenario_creator + + def test_missing_average_scenario_creator(self): + dummy = types.SimpleNamespace() # no average_scenario_creator + with self.assertRaises(RuntimeError) as ctx: + vss.vss_prep(dummy, _FakeCfg({})) + self.assertIn("average_scenario_creator", str(ctx.exception)) + + def test_multistage_rejected(self): + with self.assertRaises(RuntimeError) as ctx: + vss.vss_prep(farmer, _FakeCfg({"branching_factors": [2, 3]})) + self.assertIn("two-stage", str(ctx.exception)) + + def test_cvar_rejected(self): + with self.assertRaises(RuntimeError) as ctx: + vss.vss_prep(farmer, _FakeCfg({"cvar": True})) + self.assertIn("cvar", str(ctx.exception).lower()) + + def test_admm_rejected(self): + with self.assertRaises(RuntimeError) as ctx: + vss.vss_prep(farmer, _FakeCfg({"admm": True})) + self.assertIn("ADMM", str(ctx.exception)) + + def test_proper_bundles_rejected(self): + with self.assertRaises(RuntimeError) as ctx: + vss.vss_prep(farmer, _FakeCfg({"scenarios_per_bundle": 2})) + self.assertIn("bundles", str(ctx.exception)) + + +class TestVssHelpers(unittest.TestCase): + """Sign convention and bound-reduction helpers (no solver needed).""" + + def test_vss_value_min(self): + # minimize: VSS = EEV - RP + self.assertAlmostEqual(vss._vss_value(True, rp=-100.0, eev=-90.0), 10.0) + + def test_vss_value_max(self): + # maximize: VSS = RP - EEV + self.assertAlmostEqual(vss._vss_value(False, rp=100.0, eev=90.0), 10.0) + + def test_vss_value_infinite_eev(self): + self.assertTrue(math.isinf(vss._vss_value(True, rp=-100.0, eev=math.inf))) + + def test_reduce_bounds_min(self): + # min: inner is the tightest (lowest) incumbent, outer the highest bound + inner, outer = vss._reduce_bounds(MPI.COMM_WORLD, 5.0, 3.0, is_min=True) + self.assertEqual((inner, outer), (5.0, 3.0)) + + def test_reduce_bounds_max(self): + inner, outer = vss._reduce_bounds(MPI.COMM_WORLD, 5.0, 7.0, is_min=False) + self.assertEqual((inner, outer), (5.0, 7.0)) + + +@unittest.skipIf(not solver_available, "no solver is available") +class TestVssEndToEnd(unittest.TestCase): + """End-to-end VSS on farmer (3 scenarios).""" + + def _ef(self, kwargs): + snames = farmer.scenario_names_creator(3) + ef = mpisppy.opt.ef.ExtensiveForm( + {"solver": solver_name}, snames, farmer.scenario_creator, + scenario_creator_kwargs=kwargs) + ef.solve_extensive_form(tee=False) + return ef + + def test_ef_minimize(self): + cfg = _make_cfg(["--EF", "--EF-solver-name", solver_name]) + kwargs = farmer.kw_creator(cfg) + ef = self._ef(kwargs) + res = vss.do_vss(farmer, cfg, farmer.scenario_creator, kwargs, + farmer.scenario_denouement, ef=ef) + self.assertTrue(res["is_min"]) + self.assertAlmostEqual(res["RP"], ef.get_objective_value(), places=3) + # VSS = EEV - RP, and VSS >= 0 (the mean-value solution can't beat RP) + self.assertAlmostEqual(res["VSS"], res["EEV"] - res["RP"], places=3) + self.assertGreaterEqual(res["VSS"], -1e-6) + # regression on the known farmer-3 values + self.assertAlmostEqual(res["RP"], -108390.0, delta=1.0) + self.assertAlmostEqual(res["EV"], -118600.0, delta=1.0) + self.assertAlmostEqual(res["EEV"], -107240.0, delta=1.0) + self.assertAlmostEqual(res["VSS"], 1150.0, delta=1.0) + + def test_ef_maximize(self): + cfg = _make_cfg(["--EF", "--EF-solver-name", solver_name]) + kwargs = dict(farmer.kw_creator(cfg), sense=pyo.maximize) + ef = self._ef(kwargs) + res = vss.do_vss(farmer, cfg, farmer.scenario_creator, kwargs, + farmer.scenario_denouement, ef=ef) + self.assertFalse(res["is_min"]) + # maximize: VSS = RP - EEV, still >= 0 + self.assertAlmostEqual(res["VSS"], res["RP"] - res["EEV"], places=3) + self.assertGreaterEqual(res["VSS"], -1e-6) + + def test_decomposition_bracket(self): + # A stub wheel exercises the incumbent+bracket branch without cylinders. + # Bounds chosen to bracket the true RP (~-108390) for minimization. + cfg = _make_cfg([]) + kwargs = farmer.kw_creator(cfg) + fake_wheel = types.SimpleNamespace(BestInnerBound=-108389.0, + BestOuterBound=-108508.0) + res = vss.do_vss(farmer, cfg, farmer.scenario_creator, kwargs, + farmer.scenario_denouement, wheel=fake_wheel) + self.assertEqual(res["RP"], -108389.0) # incumbent + self.assertEqual(res["inner"], -108389.0) + self.assertEqual(res["outer"], -108508.0) + self.assertIsNotNone(res["vss_bracket"]) + lo, hi = res["vss_bracket"] + self.assertAlmostEqual(lo, res["EEV"] - res["inner"], places=4) + self.assertAlmostEqual(hi, res["EEV"] - res["outer"], places=4) + self.assertLessEqual(lo, hi) + + def test_eev_infeasible_gives_inf(self): + # Fixing the first stage x=5 is infeasible in the scenario whose + # demand is 50, so EEV (hence VSS) is +inf. + def _infeas_creator(sname, **kw): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 100)) + demand = 5.0 if sputils.extract_num(sname) == 0 else 50.0 + m.meet = pyo.Constraint(expr=m.x >= demand) + m.obj = pyo.Objective(expr=m.x, sense=pyo.minimize) + m._mpisppy_probability = 0.5 + sputils.attach_root_node(m, m.x, [m.x]) + return m + + cfg = types.SimpleNamespace(solver_name=solver_name) + eev, names = vss._compute_eev( + cfg, _infeas_creator, {}, ["Scenario0", "Scenario1"], + np.array([5.0])) + self.assertTrue(math.isinf(eev)) + self.assertIn("Scenario1", names) + + +if __name__ == "__main__": + unittest.main() diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index f0badce7f..a5b8eb90b 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -1697,6 +1697,18 @@ def mmw_args(self): default=None, ) + def vss_args(self): + self.add_to_config( + "vss", + description="After the run, report the Value of the Stochastic " + "Solution (VSS = EEV - RP). Requires the scenario module to define " + "average_scenario_creator. WARNING: computing EEV re-solves every " + "scenario with the first stage fixed and can roughly double run " + "time on large or integer models. Two-stage only.", + domain=bool, + default=False, + ) + #================ def create_parser(self,progname=None): # seldom used diff --git a/run_coverage.bash b/run_coverage.bash index 794919895..51cf0b3b0 100755 --- a/run_coverage.bash +++ b/run_coverage.bash @@ -137,6 +137,9 @@ run_phase "test_w_oscillation (serial)" \ run_phase "test_jensens (serial)" \ coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_jensens.py -v +run_phase "test_vss (serial)" \ + coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_vss.py -v + run_phase "test_feasible_xhat (serial)" \ coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_feasible_xhat.py -v From 79f85df1c6ad9677cbcddf2db5eea2514a8efd35 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Thu, 9 Jul 2026 13:58:43 -0700 Subject: [PATCH 3/6] test(vss): run test_vss in the solver-equipped regression job The unit-tests job installs no solver, so test_vss's solver-gated e2e cases (the ones asserting the actual RP/EV/EEV/VSS numbers on farmer) would skip there. Add a dedicated step in the cplex/xpress regression job -- mirroring test_feasible_xhat -- so the numerical results are actually checked in CI. The farmer-3 cases are LPs with a unique, solver-independent optimum, well within community-cplex size limits. Co-Authored-By: Claude Opus 4.8 --- .github/workflows/test_pr_and_main.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index b01340052..ad9b605ad 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -148,6 +148,10 @@ jobs: run: | coverage run $COV_ARGS -m pytest mpisppy/tests/test_feasible_xhat.py -v + - name: Test vss + run: | + coverage run $COV_ARGS -m pytest mpisppy/tests/test_vss.py -v + - name: Test prox_approx end-to-end run: | coverage run $COV_ARGS -m pytest mpisppy/tests/test_prox_approx_e2e.py -v From aca6cf5b9d7776ece3449a51b1af913428a6fd9b Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Thu, 9 Jul 2026 14:00:40 -0700 Subject: [PATCH 4/6] test(jensens): drive-by -- run test_jensens in the solver-equipped regression job Like test_vss, test_jensens was only in the no-solver unit-tests job, so its solver-gated end-to-end cases (farmer-3 average-scenario bound, sizes-3 integer-guard paths) skipped in CI. Add a dedicated regression-job step so they actually run under cplex/xpress. Both models are small and within community-cplex limits; the file is green locally. Co-Authored-By: Claude Opus 4.8 --- .github/workflows/test_pr_and_main.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index ad9b605ad..9dc9462bc 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -148,6 +148,10 @@ jobs: run: | coverage run $COV_ARGS -m pytest mpisppy/tests/test_feasible_xhat.py -v + - name: Test jensens + run: | + coverage run $COV_ARGS -m pytest mpisppy/tests/test_jensens.py -v + - name: Test vss run: | coverage run $COV_ARGS -m pytest mpisppy/tests/test_vss.py -v From 95d929cd52ef78622b44fe767ebbd3a25377c21e Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Thu, 9 Jul 2026 14:50:00 -0700 Subject: [PATCH 5/6] fix(vss): thread the run's solver options (incl. mipgap) into the VSS solves The EV and EEV solves were passing no solver options, so they ran at the solver's default mipgap -- inconsistent with the run that produced RP, and uncontrollable. Resolve solver name + options via solver_specification (EF_solver_options after --EF, else solver_options) and apply them to both the mean-value solve and the Xhat_Eval EEV pass. Also lets --vss work on an --EF run that sets only --EF-solver-name. Because VSS = EEV - RP is a difference of two optimized values, mipgap matters; document that the VSS solves inherit the run's solver-options string and that a loose gap makes a small VSS unreliable. Also reconcile the cost docs: drop the "roughly doubles run time on large/ integer models" claim (backwards -- fixing the first stage decouples the scenarios and makes each EEV solve easier; for decomposition EEV is ~one iteration's worth). State honestly that the added time is model-dependent and can be significant for some problems. Adds test_ef_solver_options_threaded; updates _compute_eev's signature and its infeasible-path test caller. Co-Authored-By: Claude Opus 4.8 --- doc/designs/vss_design.md | 37 ++++++++++++++++++++++++++-------- doc/src/vss.rst | 30 ++++++++++++++++++++++------ mpisppy/generic/vss.py | 42 ++++++++++++++++++++++++++++----------- mpisppy/tests/test_vss.py | 17 +++++++++++++--- mpisppy/utils/config.py | 8 +++++--- 5 files changed, 102 insertions(+), 32 deletions(-) diff --git a/doc/designs/vss_design.md b/doc/designs/vss_design.md index a5d6dfabd..4babd73e3 100644 --- a/doc/designs/vss_design.md +++ b/doc/designs/vss_design.md @@ -98,14 +98,35 @@ three pieces cost very differently: | :--- | :--- | :--- | | `RP` | already done by the run | free | | `EV` / `x̄` | one deterministic solve of the average scenario | cheap | -| **`EEV`** | **fix `x̄`, solve the recourse subproblem for every scenario** | **can rival the original solve** | - -`EEV` is the expensive piece. It is a *full second pass over all -scenarios* — `N` second-stage solves (or one first-stage-fixed extensive -form). For a hard MIP with many scenarios, `--vss` can roughly **double** -the wall-clock of the job. On easy LPs it is negligible; on large -integer recourse it is not. The docs must say this plainly so nobody is -surprised when `--vss` turns a 10-minute run into a 20-minute run. +| **`EEV`** | **fix `x̄`, solve the recourse subproblem for every scenario** | **an extra pass over all scenarios** | + +`EEV` is the added piece: `N` second-stage solves with the first stage +fixed. How much wall-clock that adds is genuinely hard to state in general, +and the naive "it's a second solve, so ~2x" intuition is usually wrong. +Fixing `x̄` *decouples* the scenarios and removes the here-and-now decisions, +so each `EEV` subproblem is smaller and easier than the original coupled +solve: + +- vs. an **EF** run, `EEV` solves `N` decoupled fixed-first-stage + subproblems instead of one big coupled model — typically much cheaper, + and dramatically so for a MIP (fixing integer first-stage vars is exactly + what makes the recourse easy). +- vs. a **decomposition** run, `EEV` is roughly *one iteration's* worth of + subproblem solves, against the run's many iterations — a small fraction. + +So `--vss` adds real but usually sub-linear work; it becomes noticeable +mainly with very many scenarios or genuinely expensive recourse. The docs +should say this honestly rather than promise a fixed multiplier. + +**Solver options / mipgap.** The `EV` and `EEV` solves reuse the run's +solver and `solver_options` (via `solver_specification`) — `EF_solver_options` +after an `--EF` run, `solver_options` after a decomposition run — so all +three numbers are solved consistently, and any `mipgap` in that option +string applies to the VSS solves too. This matters because `VSS = EEV - RP` +is a *difference* of two optimized values: with a loose gap a small VSS is +gap-noise. V1 threads only the solver-options *string*, not the separate +`EF_mipgap` / `*_iter*_mipgap` knobs or a `*_solver_options_file` (a +possible later refinement). A second, subtler cost: if `RP` came from a decomposition run that did **not** close the optimality gap, then `RP` is only known to lie in a diff --git a/doc/src/vss.rst b/doc/src/vss.rst index 1c0e58f91..2c00279de 100644 --- a/doc/src/vss.rst +++ b/doc/src/vss.rst @@ -14,13 +14,15 @@ once the real uncertainty shows up. A large VSS justifies building a stochastic model; a near-zero VSS says the deterministic shortcut was almost as good. -.. warning:: +.. note:: - **Computing VSS can be expensive.** The mean-value solve is cheap, but - the ``EEV`` term (below) re-solves the recourse problem for *every* - scenario with the first stage fixed — a full second pass over all - scenarios. On large or integer models, ``--vss`` can roughly **double** - the run time. It is off by default. + **Computing VSS does extra work.** ``RP`` comes for free from the run and + the mean-value solve is cheap, but ``EEV`` re-solves the recourse problem + for *every* scenario once with the first stage fixed. Fixing the first + stage decouples the scenarios and usually makes each solve easier than the + original, so how much wall-clock this adds is hard to predict: for many + problems it is a modest fraction of the run, but with very many scenarios + or expensive recourse it can be significant. ``--vss`` is off by default. The three numbers (minimization) -------------------------------- @@ -81,6 +83,22 @@ produces:: VSS = EEV - RP (sense=min) : 1150 (1.06% of |RP|) ============================================= +Solver options and mipgap +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The EV and EEV solves **reuse the run's solver and solver options** so that +all three numbers are solved the same way. After an ``--EF`` run they use +``--EF-solver-options`` (falling back to ``--solver-options`` if the former +is unset); after a decomposition run they use ``--solver-options``. Any +``mipgap`` in that option string therefore applies to the VSS solves too. + +Because ``VSS = EEV - RP`` is a *difference* of two optimized values, the +mipgap matters: if you solve with a loose gap, a small VSS can be dominated +by that gap and is not trustworthy. When VSS is small relative to ``RP``, +solve tightly, for example ``--EF-solver-options "mipgap=1e-7"``. (VSS +inherits only the solver-options *string*; it does not read the separate +``--EF-mipgap`` / ``--*-iter*-mipgap`` knobs or a ``--*-solver-options-file``.) + Exact vs. bracketed RP ~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/mpisppy/generic/vss.py b/mpisppy/generic/vss.py index 0b07db1aa..7865f22b6 100644 --- a/mpisppy/generic/vss.py +++ b/mpisppy/generic/vss.py @@ -20,8 +20,14 @@ Two-stage only in this version. See doc/designs/vss_design.md. -WARNING: computing EEV re-solves every scenario with the first stage fixed; -on large or integer models this can rival the cost of the original run. +Cost: computing EEV re-solves every scenario once with the first stage +fixed. Fixing the first stage decouples the scenarios, so each solve is +usually easier than the original; how much wall-clock this adds depends on +the model and can be significant for some problems (very many scenarios or +expensive recourse). The EV and EEV solves reuse the run's solver and +solver_options (including any mipgap) via solver_specification, so all three +numbers are solved consistently -- but VSS is a difference of two optimized +values, so a loose mipgap makes a small VSS unreliable. """ import math @@ -32,6 +38,7 @@ import mpisppy.utils.sputils as sputils from mpisppy import global_toc, MPI from mpisppy.utils import xhat_eval +from mpisppy.utils import solver_spec from mpisppy.generic.parsing import name_lists, proper_bundles @@ -88,15 +95,23 @@ def do_vss(module, cfg, scenario_creator, scenario_creator_kwargs, all_scenario_names, _ = name_lists(module, cfg) + # The VSS solves (EV and EEV) reuse the SAME solver and solver_options as + # the run that produced RP -- including any mipgap in that option string -- + # so the three numbers are solved consistently. After an EF run prefer the + # EF spec (falling back to the default); after a decomposition run use the + # default (subproblem) spec. + prefix = ["EF", ""] if ef is not None else "" + _, solver_name, solver_options = solver_spec.solver_specification(cfg, prefix) + # EV problem: solve the mean-value scenario for its objective and its # first-stage solution x_bar. ev_obj, x_bar, is_min = _solve_average_scenario( - module, cfg, scenario_creator_kwargs) + module, solver_name, solver_options, scenario_creator_kwargs) # EEV: fix x_bar in the first stage, evaluate across all scenarios. eev, infeasible_names = _compute_eev( - cfg, scenario_creator, scenario_creator_kwargs, - all_scenario_names, x_bar) + solver_name, solver_options, scenario_creator, + scenario_creator_kwargs, all_scenario_names, x_bar) # RP: exact from EF, or incumbent (+ bracket) from the wheel. if ef is not None: @@ -133,7 +148,8 @@ def do_vss(module, cfg, scenario_creator, scenario_creator_kwargs, return result -def _solve_average_scenario(module, cfg, scenario_creator_kwargs): +def _solve_average_scenario(module, solver_name, solver_options, + scenario_creator_kwargs): """Build and solve the mean-value scenario. Return (EV_objective, x_bar, is_minimizing) where x_bar is the ROOT first-stage solution as a 1-D np.ndarray in nonant_vardata_list order. @@ -148,7 +164,9 @@ def _solve_average_scenario(module, cfg, scenario_creator_kwargs): "--vss is two-stage only; average_scenario_creator must return a " "model with exactly one tree node (ROOT)." ) - solver = pyo.SolverFactory(cfg.solver_name) + solver = pyo.SolverFactory(solver_name) + for k, v in (solver_options or {}).items(): + solver.options[k] = v if sputils.is_persistent(solver): solver.set_instance(avg) results = solver.solve(tee=False) @@ -168,8 +186,8 @@ def _solve_average_scenario(module, cfg, scenario_creator_kwargs): return ev_obj, x_bar, is_min -def _compute_eev(cfg, scenario_creator, scenario_creator_kwargs, - all_scenario_names, x_bar): +def _compute_eev(solver_name, solver_options, scenario_creator, + scenario_creator_kwargs, all_scenario_names, x_bar): """Fix x_bar in the first stage and evaluate expected cost across all scenarios. Return (EEV, infeasible_scenario_names). EEV is math.inf if the mean-value first stage is infeasible in any scenario. @@ -181,9 +199,9 @@ def _compute_eev(cfg, scenario_creator, scenario_creator_kwargs, "iter0_solver_options": None, "iterk_solver_options": None, "display_timing": False, - "solver_name": cfg.solver_name, + "solver_name": solver_name, "verbose": False, - "solver_options": None, + "solver_options": solver_options, } ev = xhat_eval.Xhat_Eval( options, @@ -199,7 +217,7 @@ def _compute_eev(cfg, scenario_creator, scenario_creator_kwargs, # compute_val_at_nonant=False => need_solution=False in solve_one, so a # clean per-scenario infeasibility does NOT raise. That matters for MPI: # a raise on only some ranks would deadlock the collectives below. - ev.solve_loop(solver_options=None, gripe=True, tee=False, + ev.solve_loop(solver_options=solver_options, gripe=True, tee=False, compute_val_at_nonant=False) local_infeasible = [ diff --git a/mpisppy/tests/test_vss.py b/mpisppy/tests/test_vss.py index f4aa57e17..3c16c9f03 100644 --- a/mpisppy/tests/test_vss.py +++ b/mpisppy/tests/test_vss.py @@ -146,6 +146,18 @@ def test_ef_minimize(self): self.assertAlmostEqual(res["EEV"], -107240.0, delta=1.0) self.assertAlmostEqual(res["VSS"], 1150.0, delta=1.0) + def test_ef_solver_options_threaded(self): + # A mipgap in --EF-solver-options must flow into the VSS (EV and EEV) + # solves without error. farmer is an LP, so the value doesn't change + # the result; this locks the threading, not the numeric effect. + cfg = _make_cfg(["--EF", "--EF-solver-name", solver_name, + "--EF-solver-options", "mipgap=1e-9"]) + kwargs = farmer.kw_creator(cfg) + ef = self._ef(kwargs) + res = vss.do_vss(farmer, cfg, farmer.scenario_creator, kwargs, + farmer.scenario_denouement, ef=ef) + self.assertAlmostEqual(res["VSS"], 1150.0, delta=1.0) + def test_ef_maximize(self): cfg = _make_cfg(["--EF", "--EF-solver-name", solver_name]) kwargs = dict(farmer.kw_creator(cfg), sense=pyo.maximize) @@ -188,10 +200,9 @@ def _infeas_creator(sname, **kw): sputils.attach_root_node(m, m.x, [m.x]) return m - cfg = types.SimpleNamespace(solver_name=solver_name) eev, names = vss._compute_eev( - cfg, _infeas_creator, {}, ["Scenario0", "Scenario1"], - np.array([5.0])) + solver_name, None, _infeas_creator, {}, + ["Scenario0", "Scenario1"], np.array([5.0])) self.assertTrue(math.isinf(eev)) self.assertIn("Scenario1", names) diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index a5b8eb90b..5a8375117 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -1702,9 +1702,11 @@ def vss_args(self): "vss", description="After the run, report the Value of the Stochastic " "Solution (VSS = EEV - RP). Requires the scenario module to define " - "average_scenario_creator. WARNING: computing EEV re-solves every " - "scenario with the first stage fixed and can roughly double run " - "time on large or integer models. Two-stage only.", + "average_scenario_creator. Computing EEV re-solves every scenario " + "once with the first stage fixed; the added time depends on the " + "model and can be significant for some problems. The EV/EEV solves " + "reuse the run's solver options (--EF-solver-options for --EF, else " + "--solver-options), including mipgap. Two-stage only.", domain=bool, default=False, ) From e3664d0c9ce6165d1d60d875c65eb662281afcc9 Mon Sep 17 00:00:00 2001 From: David L Woodruff Date: Tue, 14 Jul 2026 14:07:20 -0700 Subject: [PATCH 6/6] fix(vss): don't claim RP is "exact" when the EF has a nonzero mipgap get_objective_value() returns the EF solver's incumbent, which for a MIP solved to a nonzero gap is only known to within that gap -- but do_vss unconditionally labelled RP "[EF, exact]", contradicting the doc's own mipgap warning. The dual bound existed at solve time (results.problem lower/upper_bound) but do_EF discarded it. do_EF now stashes ef.ef_objective_bounds; do_vss reads the sense- appropriate dual bound and, when the relative gap exceeds 1e-9, relabels RP "EF incumbent, gap=..." and brackets VSS with the existing bracket machinery (same inner/outer semantics as the decomposition case). LP, zero-gap, or no-bound-reported keeps "EF, exact" with no bracket. Docs (vss.rst + module docstring) updated for the exact-vs-bracketed distinction. Adds helper unit tests plus solver-gated end-to-end cases for the no-bound, nonzero-gap, and zero-gap paths. Co-Authored-By: Claude Opus 4.8 (1M context) --- doc/src/vss.rst | 25 +++++++++------ mpisppy/generic/ef.py | 10 ++++++ mpisppy/generic/vss.py | 56 +++++++++++++++++++++++++++++----- mpisppy/tests/test_vss.py | 64 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 139 insertions(+), 16 deletions(-) diff --git a/doc/src/vss.rst b/doc/src/vss.rst index 2c00279de..cbf90c34a 100644 --- a/doc/src/vss.rst +++ b/doc/src/vss.rst @@ -35,8 +35,10 @@ Write the two-stage program with first-stage vector :math:`x`, random data RP = \min_x \; \big\{\, c\cdot x + \mathbb{E}_\xi[\,Q(x,\xi)\,] \,\big\}. - **RP** — the *Recourse Problem* optimum: the here-and-now solution the run - already computes. It is **exact** from an ``--EF`` run and the - **incumbent** (with a bracket) from a decomposition run. + already computes. From an ``--EF`` run it is **exact** when the EF is solved + to a zero gap — always the case for an LP — and otherwise the **incumbent** + bracketed by the solver's dual bound (see *Exact vs. bracketed RP* below). + From a decomposition run it is the **incumbent** (with a bracket). - **EV** — the *mean-value problem* :math:`\min_x \{ c\cdot x + Q(x,\bar\xi)\}` with :math:`\bar\xi` the average of the data over the scenario set. Its first-stage solution is @@ -102,19 +104,24 @@ inherits only the solver-options *string*; it does not read the separate Exact vs. bracketed RP ~~~~~~~~~~~~~~~~~~~~~~~~ -From an ``--EF`` run, ``RP`` is exact. From a decomposition (cylinders) run -that did not close the optimality gap, ``RP`` is only known to lie between -the outer and inner bounds, so VSS is reported both as a conservative point -value (against the incumbent) and as a bracket:: +From an ``--EF`` run that reaches a zero gap, ``RP`` is exact. But an EF MIP +solved to a *nonzero* mipgap returns only an incumbent, known to lie between +the solver's dual bound and that incumbent — so ``RP`` is then reported the +same way as the decomposition case: a conservative point value (the +incumbent) plus a bracket, with the label naming the gap:: - RP (stochastic solution, here-and-now): -108389 [decomposition incumbent] + RP (stochastic solution, here-and-now): -108389 [EF incumbent, gap=0.001] optimality bracket [outer, inner] : [-108508, -108389] ... VSS = EEV - RP (sense=min) : 1149.33 (1.06% of |RP|) VSS bracket from RP bracket : [1149.33, 1268.39] -The width of the VSS bracket is exactly the unclosed optimality gap; tighten -the run (or use ``--EF``) for a sharper VSS. +A decomposition (cylinders) run that did not close its optimality gap +brackets ``RP`` the same way (labelled ``[decomposition incumbent]``). In +both cases the width of the VSS bracket is exactly the unclosed optimality +gap; solve the EF to a tighter mipgap (or an LP) for a sharper VSS. When the +solver reports no dual bound, ``RP`` is labelled ``[EF, exact]`` and no +bracket is shown. Infeasible mean-value solution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/mpisppy/generic/ef.py b/mpisppy/generic/ef.py index 6590ffd92..3d57fd241 100644 --- a/mpisppy/generic/ef.py +++ b/mpisppy/generic/ef.py @@ -72,6 +72,16 @@ def do_EF(module, cfg, scenario_creator, scenario_creator_kwargs, if not pyo.check_optimal_termination(results): print("Warning: non-optimal solver termination") + # Stash the solver's objective (dual) bound alongside the incumbent so a + # VSS report can tell whether RP is exact or only known to within the MIP + # gap -- get_objective_value() returns only the incumbent. For an LP the + # bound equals the objective; None if the solver reported no bound. + try: + ef.ef_objective_bounds = (results.problem.lower_bound, + results.problem.upper_bound) + except (AttributeError, KeyError): + ef.ef_objective_bounds = None + global_toc(f"EF objective: {ef.get_objective_value()}") if ef.extensions is not None: diff --git a/mpisppy/generic/vss.py b/mpisppy/generic/vss.py index 7865f22b6..f960091a0 100644 --- a/mpisppy/generic/vss.py +++ b/mpisppy/generic/vss.py @@ -11,8 +11,10 @@ VSS = EEV - RP (minimization; VSS = RP - EEV for maximization), where RP = the stochastic-program optimum -- the here-and-now solution the run - already computed: exact from an --EF run, or the decomposition - incumbent (BestInnerBound), with a bracket from BestOuterBound. + already computed: from an --EF run it is exact when the EF is solved + to zero gap (always so for an LP), otherwise the incumbent bracketed + by the solver's dual bound; from a decomposition run it is the + incumbent (BestInnerBound) bracketed by BestOuterBound. EV = the mean-value (expected-value) problem: replace the random data by its average and solve. x_bar is its first-stage solution. EEV = x_bar fixed in the first stage and evaluated honestly across every @@ -113,12 +115,27 @@ def do_vss(module, cfg, scenario_creator, scenario_creator_kwargs, solver_name, solver_options, scenario_creator, scenario_creator_kwargs, all_scenario_names, x_bar) - # RP: exact from EF, or incumbent (+ bracket) from the wheel. + # RP: from EF (exact when solved to zero gap; otherwise the incumbent with + # a bracket from the solver's dual bound), or the incumbent (+ bracket) + # from the wheel. if ef is not None: - rp_point = comm.bcast(ef.get_objective_value() if comm.Get_rank() == 0 - else None, root=0) - rp_source = "EF, exact" - inner = outer = None + if comm.Get_rank() == 0: + payload = (ef.get_objective_value(), _ef_dual_bound(ef, is_min)) + else: + payload = None + incumbent, dual = comm.bcast(payload, root=0) + rp_point = incumbent + # A MIP left with a nonzero gap makes the incumbent only known to within + # [dual, incumbent]; relabel and bracket so the "exact" claim is not + # overstated. An LP (or a MIP solved to zero gap) reports dual==incumbent. + if (dual is not None and incumbent is not None + and math.isfinite(incumbent) + and _rel_gap(dual, incumbent) > _EF_EXACT_RELGAP): + inner, outer = incumbent, dual + rp_source = f"EF incumbent, gap={_rel_gap(dual, incumbent):.2g}" + else: + rp_source = "EF, exact" + inner = outer = None else: inner, outer = _reduce_bounds(comm, wheel.BestInnerBound, wheel.BestOuterBound, is_min) @@ -240,6 +257,31 @@ def _compute_eev(solver_name, solver_options, scenario_creator, return ev.Eobjective(), [] +# Relative gap at or below which an EF RP is reported as exact (an LP, or a +# MIP solved to a negligible gap, reports dual bound == incumbent). +_EF_EXACT_RELGAP = 1e-9 + + +def _ef_dual_bound(ef, is_min): + """The EF solver's best objective (dual) bound for the run's sense, or + None if unavailable. do_EF stashes (lower_bound, upper_bound) as + ef.ef_objective_bounds; the sense-appropriate one brackets the incumbent. + """ + bounds = getattr(ef, "ef_objective_bounds", None) + if bounds is None: + return None + val = bounds[0] if is_min else bounds[1] + if val is None or not math.isfinite(val): + return None + return val + + +def _rel_gap(dual, incumbent): + """Relative gap between the dual bound and the incumbent objective.""" + denom = abs(incumbent) if incumbent != 0 else 1.0 + return abs(incumbent - dual) / denom + + def _reduce_bounds(comm, best_inner, best_outer, is_min): """Make the wheel's best inner/outer bounds available on every rank. diff --git a/mpisppy/tests/test_vss.py b/mpisppy/tests/test_vss.py index 3c16c9f03..ba2ecb107 100644 --- a/mpisppy/tests/test_vss.py +++ b/mpisppy/tests/test_vss.py @@ -116,6 +116,24 @@ def test_reduce_bounds_max(self): inner, outer = vss._reduce_bounds(MPI.COMM_WORLD, 5.0, 7.0, is_min=False) self.assertEqual((inner, outer), (5.0, 7.0)) + def test_ef_dual_bound_picks_by_sense(self): + ef = types.SimpleNamespace(ef_objective_bounds=(3.0, 7.0)) + self.assertEqual(vss._ef_dual_bound(ef, is_min=True), 3.0) # lower + self.assertEqual(vss._ef_dual_bound(ef, is_min=False), 7.0) # upper + + def test_ef_dual_bound_missing_or_infinite(self): + self.assertIsNone(vss._ef_dual_bound(types.SimpleNamespace(), is_min=True)) + ef_none = types.SimpleNamespace(ef_objective_bounds=None) + self.assertIsNone(vss._ef_dual_bound(ef_none, is_min=True)) + ef_inf = types.SimpleNamespace(ef_objective_bounds=(-math.inf, math.inf)) + self.assertIsNone(vss._ef_dual_bound(ef_inf, is_min=True)) + self.assertIsNone(vss._ef_dual_bound(ef_inf, is_min=False)) + + def test_rel_gap(self): + self.assertAlmostEqual(vss._rel_gap(-108500.0, -108390.0), + 110.0 / 108390.0) + self.assertEqual(vss._rel_gap(0.0, 0.0), 0.0) # denom guard, no div0 + @unittest.skipIf(not solver_available, "no solver is available") class TestVssEndToEnd(unittest.TestCase): @@ -169,6 +187,52 @@ def test_ef_maximize(self): self.assertAlmostEqual(res["VSS"], res["RP"] - res["EEV"], places=3) self.assertGreaterEqual(res["VSS"], -1e-6) + def test_ef_exact_when_no_bound_reported(self): + # An EF object without ef_objective_bounds (or with none reported) is + # labelled exact and carries no bracket -- the back-compat default. + cfg = _make_cfg(["--EF", "--EF-solver-name", solver_name]) + kwargs = farmer.kw_creator(cfg) + ef = self._ef(kwargs) # solve_extensive_form doesn't set the bounds + res = vss.do_vss(farmer, cfg, farmer.scenario_creator, kwargs, + farmer.scenario_denouement, ef=ef) + self.assertEqual(res["rp_source"], "EF, exact") + self.assertIsNone(res["inner"]) + self.assertIsNone(res["outer"]) + self.assertIsNone(res["vss_bracket"]) + + def test_ef_nonzero_gap_brackets_rp(self): + # A MIP left with a nonzero gap: the incumbent is bracketed by the + # solver's dual bound, so RP is relabelled and a VSS bracket appears. + cfg = _make_cfg(["--EF", "--EF-solver-name", solver_name]) + kwargs = farmer.kw_creator(cfg) + ef = self._ef(kwargs) + incumbent = ef.get_objective_value() # ~ -108390 (min) + dual = incumbent - 110.0 # dual bound below it + ef.ef_objective_bounds = (dual, incumbent) # (lower, upper) + res = vss.do_vss(farmer, cfg, farmer.scenario_creator, kwargs, + farmer.scenario_denouement, ef=ef) + self.assertIn("gap=", res["rp_source"]) + self.assertEqual(res["RP"], incumbent) # point value = incumbent + self.assertEqual(res["inner"], incumbent) + self.assertEqual(res["outer"], dual) + self.assertIsNotNone(res["vss_bracket"]) + lo, hi = res["vss_bracket"] + self.assertAlmostEqual(lo, res["EEV"] - res["inner"], places=4) + self.assertAlmostEqual(hi, res["EEV"] - res["outer"], places=4) + self.assertLessEqual(lo, hi) + + def test_ef_zero_gap_stays_exact(self): + # dual bound == incumbent (an LP, or a MIP solved tight) => exact label. + cfg = _make_cfg(["--EF", "--EF-solver-name", solver_name]) + kwargs = farmer.kw_creator(cfg) + ef = self._ef(kwargs) + incumbent = ef.get_objective_value() + ef.ef_objective_bounds = (incumbent, incumbent) + res = vss.do_vss(farmer, cfg, farmer.scenario_creator, kwargs, + farmer.scenario_denouement, ef=ef) + self.assertEqual(res["rp_source"], "EF, exact") + self.assertIsNone(res["vss_bracket"]) + def test_decomposition_bracket(self): # A stub wheel exercises the incumbent+bracket branch without cylinders. # Bounds chosen to bracket the true RP (~-108390) for minimization.