diff --git a/doc/designs/mutable_scenario_probabilities_design.md b/doc/designs/mutable_scenario_probabilities_design.md new file mode 100644 index 000000000..8ba522b1a --- /dev/null +++ b/doc/designs/mutable_scenario_probabilities_design.md @@ -0,0 +1,442 @@ +# Mutable scenario probabilities + +Status: phases 0–1 implemented and verified (addresses issue #797); phases +2–4 remain (§9). Covers the current behavior (§1), the use case and goals +(§2–3), the design (§4), the EF path in detail (§5), the PH/decomposition +path (§6), API and compatibility (§7), and open questions (§8). + +Implemented so far (EF path, the issue's request): +`sputils.has_persistent_solve_api` (recognizes APPSI / `pyomo.contrib.solver` +as persistent for the EF workflow); `mutable_probability` option on +`_create_EF_from_scen_dict` and `ExtensiveForm` (option-B Param objective); +`ExtensiveForm.set_scenario_probabilities`; and a `reuse_instance` argument to +`solve_extensive_form`. Verified end-to-end on the farmer example against a +rebuild oracle — machine-precision agreement across a probability sweep, with +`set_instance` called once — for **both** `appsi_highs` (auto-tracks the +change) and `gurobi_persistent` (requires the explicit `set_objective` +re-push). Tests: `Test_mutable_probability` in `mpisppy/tests/test_ef_ph.py`. + +--- + +## 0. Motivation (issue #797) + +A user solves a two-stage stochastic program with the HiGHS *persistent* +solver in a rolling-horizon loop. The persistent model is built once, and +between solves only data changes. For each roll they need to change the +**scenario probabilities** and re-solve. + +Today that is not possible without rebuilding: mpi-sppy bakes each +scenario's probability into the Extensive Form objective as a Python +**float constant** when the EF is constructed. Changing a probability +requires reconstructing the EF objective (and, for a persistent solver, +re-loading the instance), which defeats the point of a persistent model. + +The request: represent scenario probabilities as **mutable Pyomo +parameters** so they can be updated in place, and the (persistent) solver +re-solved cheaply. + +--- + +## 1. Current behavior + +### 1.1 Where probabilities live + +Each scenario model carries a scalar attribute `_mpisppy_probability` +(a Python `float`), set by the user's `scenario_creator` or defaulted to +uniform in `sputils.py`: + +```python +# sputils.py ~281-288 (_check_get_default_probability / attach) +if probs_specified and all(... == "uniform" ...): + ... + scen._mpisppy_probability = 1 / total_number_of_scenarios +``` + +For multistage problems, `ScenarioNode.cond_prob` holds the conditional +probability at each tree node, and `spbase._compute_unconditional_node_probabilities` +derives, per scenario and node, an unconditional coefficient +`prob_coeff`: + +```python +# spbase.py ~419-427 +root.uncond_prob = 1.0 +for parent, child in zip(node_list[:-1], node_list[1:]): + child.uncond_prob = parent.uncond_prob * child.cond_prob +... +s._mpisppy_data.prob_coeff[node.name] = (s._mpisppy_probability / node.uncond_prob) +s._mpisppy_data.prob0_mask[node.name] = 1.0 +``` + +So there are **two** representations of "how much this scenario counts": + +- `_mpisppy_probability` — the scenario's (leaf) probability; used to + weight the **objective**. +- `_mpisppy_data.prob_coeff[ndn]` — per-node conditional coefficient; + used to weight **xbar / W / rho / residuals** in PH-family algorithms. + +Both are plain floats (or numpy float arrays, for variable probability). + +### 1.2 Where the probability enters the EF objective + +`sputils._create_EF_from_scen_dict` (the function behind +`ExtensiveForm` and behind PH bundles) builds the EF objective by +folding each scenario's float probability directly into the expression: + +```python +# sputils.py 348-364 +EF_instance._mpisppy_probability = 0 +for (sname, scenario_instance) in scen_dict.items(): + EF_instance.add_component(sname, scenario_instance) + ... + obj_func = scenario_objs[0] + EF_instance.EF_Obj.expr += scenario_instance._mpisppy_probability * obj_func.expr + EF_instance._mpisppy_probability += scenario_instance._mpisppy_probability +# normalize (matters for bundles; a no-op weight-sum==1 for a full EF) +EF_instance.EF_Obj.expr /= EF_instance._mpisppy_probability +``` + +Because `_mpisppy_probability` is a `float`, the probabilities become +**constant coefficients** in the compiled objective. A persistent solver +that has already ingested this objective has no handle to update them; +the only recourse today is to rebuild the objective and re-set the +instance. + +### 1.3 Where the probability enters PH + +For PH/APH/subgradient/FWPH etc. the probability weighting is applied +numerically (not in a Pyomo expression) via `prob_coeff`: + +- `phbase.py:94,164` — xbar and weighted averages +- `convergers/norms_and_residuals.py`, `convergers/primal_dual_converger.py` +- `extensions/dyn_rho_base.py:158`, `extensions/primal_dual_rho.py:75`, + `extensions/grad_rho.py:113` +- `opt/aph.py:525` + +These read `s._mpisppy_data.prob_coeff[ndn]` fresh each iteration, so they +are *already* update-friendly **provided** the stored value is refreshed. +The proximal/objective terms attached to each subproblem (`attach_Ws_and_prox`) +do **not** multiply by probability — PH's per-scenario objective is the raw +scenario objective plus W and prox terms; probability weighting happens in +the aggregation, not in the subproblem objective. That is important: for the +PH path, "mutable probabilities" is mostly a matter of recomputing +`prob_coeff`, not of touching any Pyomo expression. + +--- + +## 2. Use case in scope + +Primary (issue #797): **ExtensiveForm** with a persistent solver, probabilities +changed between solves, no rebuild. Two-stage is the concrete ask; the design +should not preclude multistage. + +Secondary: keeping the PH/decomposition path consistent so that a user who +updates probabilities there gets correct xbar/W/rho weighting on the next +iteration. + +## 3. Goals and non-goals + +Goals: + +1. Allow scenario probabilities to be updated after model construction and + re-solved without rebuilding the Pyomo model. +2. For persistent solvers, make the update cheap (update param values + + re-push the objective, not `set_instance`). +3. Keep the default path (probabilities fixed for the life of the run) + behaving exactly as today, at no measurable overhead, and opt-in for the + mutable path. +4. One source of truth: updating a probability must keep `_mpisppy_probability`, + the EF objective, and `prob_coeff` mutually consistent. + +Non-goals: + +- Changing the **structure** of the scenario tree at runtime (adding/removing + scenarios or nodes). Only the probability *values* on an existing tree. +- Re-deriving sample-average or confidence-interval machinery on the fly. +- Automatic re-solve orchestration in the rolling-horizon loop — that stays + the user's driver code. + +--- + +## 4. Design overview + +Introduce an **opt-in mutable representation** of scenario probability, keyed +off a flag so the default stays a baked-in constant. + +- A new option, `mutable_probability` (bool, default `False`), threaded to + `ExtensiveForm` / `sputils._create_EF_from_scen_dict` (and available to + `SPBase` for the PH path). +- When enabled, each scenario's probability is stored as a Pyomo **mutable + `Param`** rather than folded in as a float, and the EF objective references + that Param. +- A single public method, `set_scenario_probabilities(mapping)`, updates the + Param values, refreshes the derived `prob_coeff`, and (for persistent + solvers / EF) re-pushes the objective. + +The flag keeps the common case free of any Param-indirection overhead and +avoids perturbing the many call sites that read `_mpisppy_probability`. + +## 5. EF path (the issue's actual request) + +### 5.1 Representation + +In `_create_EF_from_scen_dict`, when `mutable_probability` is set, attach a +mutable Param per scenario on the EF's `_mpisppy_model` block and build the +objective against it: + +```python +EF_instance._mpisppy_model.prob = pyo.Param( + EF_instance._ef_scenario_names, mutable=True, within=pyo.NonNegativeReals, + initialize={sname: scen._mpisppy_probability for sname, scen in scen_dict.items()}, +) +# objective term +EF_instance.EF_Obj.expr += EF_instance._mpisppy_model.prob[sname] * obj_func.expr +``` + +Normalization (`/= sum`) is the subtlety: the current code divides the whole +objective by the accumulated probability sum. That division only does real +work for **bundles**, where member probabilities sum to <1 and the divisor +rescales the bundle objective into a proper within-bundle conditional +expectation. For a standalone **full EF** — the issue's use case — the +probabilities already sum to 1, so the division is a no-op. + +The two considered options: + +- **(A) Normalize inside the expression with a mutable divisor.** Store an + unnormalized objective `sum_s prob[s] * obj_s` and a separate mutable Param + `prob_sum`; write the objective as `(sum_s prob[s]*obj_s) / prob_sum`, and + have `set_scenario_probabilities` update both `prob[*]` and `prob_sum`. + General (handles sum != 1) but carries a division node and an extra Param. +- **(B) Require normalized input and drop the divisor.** Impose that the + supplied probabilities sum to 1 whenever `mutable_probability` is set, and + build the objective as just `sum_s prob[s] * obj_s` with no division. + +**Recommendation: (B).** The requirement is cheap to impose *only* on the +mutable path because the mutable and needs-a-divisor regimes are disjoint in +practice: mutable probability is a full-EF feature, and bundles (the only +sum != 1 case) are built by this same function but never set the flag. So a +single branch at construction covers both without touching existing behavior: + +```python +if mutable_probability: + _validate_sum_to_one(scen_dict) # raise on violation (see §8.5) + EF_instance._mpisppy_model.prob = pyo.Param( + EF_instance._ef_scenario_names, mutable=True, within=pyo.NonNegativeReals, + initialize={sname: scen._mpisppy_probability for sname, scen in scen_dict.items()}, + ) + for (sname, scenario_instance) in scen_dict.items(): + ... + EF_instance.EF_Obj.expr += EF_instance._mpisppy_model.prob[sname] * obj_func.expr + # no /= sum +else: + # unchanged float-coefficient path, including the /= accumulated-sum + # normalization that bundles rely on + ... +``` + +The same `_validate_sum_to_one` runs inside `set_scenario_probabilities`, so +the "mutable ⇒ normalized" contract holds at build time and on every update. + +Why (B) is not just simpler but genuinely better here: + +- **No `prob_sum` Param and no division node.** This matters for the + persistent path (§5.2): re-pushing via `set_objective` becomes a plain + linear combination of Params, with nothing for the solver to re-derive + around a division. +- **Nothing existing changes.** The requirement is gated on the opt-in flag, + so the float path — and bundle scaling in particular — is untouched. + +Two guards close the one case where the disjointness assumption could be +violated: + +1. **Bad input on the mutable path.** Reject (raise) rather than silently + renormalizing — silent renormalization would secretly reintroduce the + divisor B is trying to avoid. See §8.5. +2. **Mutable probabilities requested for a bundle.** Disallow explicitly + (raise) rather than silently skipping the divisor a bundle needs. + +### 5.2 Update + re-solve + +```python +def set_scenario_probabilities(self, prob_map): + # prob_map: {scenario_name: float} (must keep the full set summing to 1) + _validate_sum_to_one(prob_map_applied_to_all_scenarios) # §5.1, §8.5 + for sname, p in prob_map.items(): + self.ef._mpisppy_model.prob[sname].value = p + getattr(self.ef, sname)._mpisppy_probability = p # §5.3 + if and hasattr(self.solver, "set_objective"): + self.solver.set_objective(self.ef.EF_Obj) # re-push objective coefficients +``` + +Under option (B) there is no `prob_sum` to update — just the per-scenario +Params. `set_objective` re-extracts the objective coefficients from the +(already-updated) Param values without touching constraints or variables — +far cheaper than `set_instance`. + +**Cross-solver behavior (verified by experiment, see §8.3).** Two persistent +families behave differently, and calling `set_objective` unconditionally on +the persistent path is correct for both: + +- *APPSI / `pyomo.contrib.solver` (e.g. `appsi_highs`, the issue's solver).* + The wrapper auto-tracks model changes on the next `solve()` — its + `update_config.update_params` and `update_objective` default to `True` — so + a mutable-Param objective change is picked up automatically even *without* + `set_objective`. Calling `set_objective` is redundant but harmless. +- *Legacy `PersistentSolver` (e.g. `gurobi_persistent`, `cplex_persistent`).* + These do **not** auto-track model changes; `set_objective` (or an explicit + coefficient update) is **required**, or the solver re-solves the stale + objective. + +**Reuse guard.** `solve_extensive_form` currently re-instances on every call +for solvers it considers persistent (ef.py:117-118). For the rolling-horizon +reuse case we want to **skip** re-instancing when the instance is already +loaded and only the objective changed. *Resolved (§8.2):* add an explicit +`reuse_instance=False` argument to `solve_extensive_form`; when `True`, skip +`set_instance` and rely on the loaded instance plus the re-pushed objective. + +**Detection caveat (found by experiment, see §8.3).** mpi-sppy does not currently +recognize `appsi_highs` as persistent at all: `ef.py:117` tests the substring +`"persistent" in solver_name` (false for `"appsi_highs"`), and +`sputils.is_persistent` tests `isinstance(solver, PersistentSolver)` +(also false — APPSI solvers are `LegacySolver` wrappers, not that base class). +So today `ExtensiveForm` treats `appsi_highs` like a non-persistent solver and +never takes the `set_instance` / `load_vars` path. This feature must therefore +extend persistence detection to the APPSI / `pyomo.contrib.solver` interfaces +(e.g. duck-type on `set_instance`/`set_objective`/`load_vars`, or check for the +contrib base classes) — otherwise the `reuse_instance` path is unreachable for +exactly the solver in the issue. This is a prerequisite for phase 1, not a +nice-to-have. + +### 5.3 Keeping `_mpisppy_probability` consistent + +`set_scenario_probabilities` must also write `scen._mpisppy_probability = p` +on the underlying scenario models, so any downstream reader (solution +reporting, `get_objective_value`, zhat evaluation) sees the same numbers the +objective used. + +## 6. PH / decomposition path + +The PH path does not bake probabilities into a Pyomo objective, so no Param is +needed there. What it needs is a way to **recompute `prob_coeff`** after a +change. Today `_compute_unconditional_node_probabilities` only computes +`prob_coeff` if it is missing (`if not hasattr(...)`), so it will not pick up +an updated `_mpisppy_probability` on its own. + +Design (implemented, phase 2): `SPBase.set_scenario_probabilities(prob_map, +check_sum=True, reset_ph_duals=True)` that + +1. updates `_mpisppy_probability` on each local scenario named in `prob_map` + (two-stage only for now; a multistage `_mpisppy_node_list` raises + `NotImplementedError` — updating `ScenarioNode.cond_prob` is a later phase), +2. forces recomputation of `uncond_prob` and `prob_coeff` via a new `force` + flag on `_compute_unconditional_node_probabilities` that bypasses the + `hasattr` compute-once short-circuit, +3. preserves any `has_variable_probability` overrides (re-applies + `_use_variable_probability_setter` after the refresh), +4. zeroes the PH multipliers `W` on each local scenario (`reset_ph_duals`, + default True; no-op for objects without PH `W` terms), and +5. optionally checks (default) with an MPI reduction that the resulting + probabilities sum to 1 (option B). + +Because xbar/W/rho all read `prob_coeff` fresh each iteration, refreshing that +dict is sufficient for the next PH iteration to be correct. + +**Warm-start caveat (motivates step 4).** At a converged PH solution every +scenario sits at the same nonanticipative point, so the probability-weighted +`xbar` equals that point *regardless of the weights*. Re-solving from there +with new probabilities but stale `W` leaves `xbar` unmoved and PH reports +immediate (false) convergence at the old solution — the exact +probability-sweep use case would silently return the wrong answer. Zeroing `W` +(step 4) breaks that consensus so the next `ph_main()` re-converges for the new +problem. Verified on farmer: fresh PH at a skewed vector and an in-place +`set_scenario_probabilities` + re-solve both reach the EF-oracle solution, +while `reset_ph_duals=False` after convergence stays stuck at the old optimum +(`test_ph_reuse_after_prob_change` / `..._without_dual_reset_stays_stuck` in +`test_ef_ph.py`). Note the EF's persistent-solver reuse (avoiding a monolithic +rebuild) is the real payoff of this feature; PH subproblems already persist +across iterations, so reusing a PH object across a sweep saves only scenario +construction — rebuilding per vector is also fine. + +This is a lower-priority companion to the EF work — the issue is specifically +about the EF — but it belongs in the same design so the two representations +stay in sync and there is one method name across both. + +## 7. API and backward compatibility + +- New option `mutable_probability` (default `False`) — no behavior change for + existing users; the objective is still a float-coefficient expression. +- New method `set_scenario_probabilities` on `ExtensiveForm` (and `SPBase`). + Purely additive. +- `_mpisppy_probability` remains a float attribute and is still the read API; + the Param is an internal, opt-in representation, not a replacement. +- CLI: a `--mutable-probability` `config.py` arg was considered but not added. + `generic_cylinders.py` has no EF-construction path that would consume it, and + a probability *sweep* is inherently a driver loop (change the vector, + re-solve) that a single CLI invocation cannot express. The CLI exposure of + the feature is instead the runnable driver script + `examples/farmer/farmer_prob_sensitivity.py`, which has its own argparse + (`--solver-name`, `--num-scens`). + +## 8. Open questions + +1. **Normalization** — design settles on option (B) (require sum-to-1 on the + mutable path, drop the divisor; see §5.1). Remaining to confirm: the exact + tolerance for the sum-to-1 check, and that no in-tree caller builds a + *bundle* with `mutable_probability` set (the guard in §5.1 should make that + an explicit error). +2. **Persistent re-instancing guard** — *resolved:* an explicit + `reuse_instance=False` argument to `solve_extensive_form`. When `True`, + skip the `set_instance` at ef.py:117 and rely on the already-loaded + instance plus the updated objective. Explicit over an invisible internal + flag so the reuse contract is visible at the call site. +3. **HiGHS `set_objective` support** — *resolved by experiment* + (appsi_highs / highspy 1.15.1). Findings: + (a) `set_objective(EF_Obj)` **does** re-extract objective coefficients from + mutable Params — re-solve moved the shared first-stage var to the new + optimum after only re-pushing the objective, no `set_instance`. + (b) Stronger: appsi_highs **auto-tracks** the change on the next `solve()` + even without `set_objective` (`update_config.update_params` / + `update_objective` default `True`), and retains its instance across + `.solve()` calls. So the mutable-Param objective "just works" for + appsi_highs; `set_objective` is redundant-but-safe there and *required* for + legacy persistent solvers (see §5.2). + (c) *New prerequisite surfaced:* mpi-sppy classifies `appsi_highs` as + **non-persistent** — `sputils.is_persistent` returns `False` and the + `"persistent" in solver_name` check at ef.py:117 is also `False` — so the + reuse path is currently unreachable for it. Persistence detection must be + extended to APPSI / `pyomo.contrib.solver` interfaces (§5.2). +4. **Scope of a probability change** — do we ever need to change the *number* + of scenarios or the tree structure between rolls? Declared out of scope + here; if the user's rolling horizon changes the scenario set, that is a + rebuild, not a probability update. +5. **Validation** — *resolved in the implementation.* Under option (B), + `set_scenario_probabilities` **raises** if the resulting full probability + vector does not sum to 1 (tolerance `1e-9`), with no silent renormalization + (which would reintroduce the divisor B removes). Partial-mapping updates + **are** allowed: scenarios omitted from the mapping keep their current + probability, as long as the resulting total still sums to 1. The check runs + *before* any value is written, so a rejected call leaves the model + unchanged (transactional). + +--- + +## 9. Suggested implementation phases + +0. **[done]** Prerequisite: extend persistence detection to APPSI / + `pyomo.contrib.solver` interfaces so `appsi_highs` (the issue's solver) is + recognized as persistent (§5.2, §8.3c). Implemented as + `sputils.has_persistent_solve_api`; `ExtensiveForm.solve_extensive_form` + uses it in place of the old `"persistent" in name` / `is_persistent` checks. +1. **[done]** EF-only, two-stage: `mutable_probability` flag, Param objective + with option-(B) normalization (require sum-to-1, no divisor), + `set_scenario_probabilities`, explicit `reuse_instance` argument to + `solve_extensive_form`. Closes the issue. Verified against a rebuild oracle + on farmer for `appsi_highs` and `gurobi_persistent`. +2. **[done]** PH path: `SPBase.set_scenario_probabilities` + `prob_coeff` + refresh (with `reset_ph_duals` for the consensus warm-start caveat, §6). + Verified on farmer against the EF oracle. +3. Multistage node probabilities and variable-probability interaction. +4. **[done]** Docs + a probability-sensitivity example under `examples/`. + `examples/farmer/farmer_prob_sensitivity.py` (sweeps the weight on one + scenario, reuses the persistent instance across the sweep) and + `doc/src/mutable_probability.rst` (linked from `index.rst` after `ef.rst`). + A dedicated `generic_cylinders` CLI flag was intentionally omitted — see §7. diff --git a/doc/src/index.rst b/doc/src/index.rst index 2f577ac23..531b84c57 100644 --- a/doc/src/index.rst +++ b/doc/src/index.rst @@ -27,6 +27,7 @@ MPI is used. generic_cylinders.rst examples.rst ef.rst + mutable_probability.rst chance_constraints.rst .. toctree:: diff --git a/doc/src/mutable_probability.rst b/doc/src/mutable_probability.rst new file mode 100644 index 000000000..4d7230cd5 --- /dev/null +++ b/doc/src/mutable_probability.rst @@ -0,0 +1,103 @@ +.. _mutable_probability: + +Mutable Scenario Probabilities +============================== + +By default, each scenario's probability is folded into the Extensive Form +(EF) objective as a floating-point constant when the model is built. Changing +a probability therefore requires rebuilding the objective and, for a +persistent solver, re-loading the instance. + +The ``mutable_probability`` option (issue #797) makes the probabilities +updatable in place. When it is set, each scenario's probability is stored as a +mutable Pyomo ``Param`` in the objective instead of a constant, so the +probability vector can be changed and the (persistent) solver re-solved +cheaply, without rebuilding the model. This is useful for probability +sensitivity studies and rolling-horizon loops that keep a fixed scenario set +and only re-weight it between solves. + +Requirements and scope +----------------------- + +- The supplied probabilities must sum to 1 (option "B" in the design). There + is no re-normalization: a vector that does not sum to 1 (within ``1e-9``) is + rejected. This keeps the objective a plain probability-weighted sum with no + division node. +- ``mutable_probability`` is a full-EF feature and is **not** supported for + scenario bundles (which rely on the normalization divisor). Requesting it + for a bundle raises an error. +- Updates are transactional: if a call is rejected, no probability is changed. + +Extensive Form +-------------- + +Build the EF with ``mutable_probability=True``, then call +``set_scenario_probabilities`` with a mapping and re-solve. Pass +``reuse_instance=True`` after the first solve so a persistent solver keeps its +loaded instance and only the objective coefficients are re-pushed: + +.. code-block:: python + + from mpisppy.opt.ef import ExtensiveForm + + ef = ExtensiveForm( + options={"solver": "gurobi_persistent"}, + all_scenario_names=scenario_names, + scenario_creator=scenario_creator, + scenario_creator_kwargs=scenario_creator_kwargs, + mutable_probability=True, + ) + + for i, prob_map in enumerate(probability_vectors): + ef.set_scenario_probabilities(prob_map) # must sum to 1 + ef.solve_extensive_form(reuse_instance=(i > 0)) + print(ef.get_objective_value(), ef.get_root_solution()) + +A partial mapping is allowed: scenarios omitted from ``prob_map`` keep their +current probability, as long as the resulting full vector still sums to 1. + +Persistent solvers, including the APPSI / ``pyomo.contrib.solver`` interface +(e.g. ``appsi_highs``, the solver in the issue) are recognized via +``sputils.has_persistent_solve_api`` so the ``reuse_instance`` path is taken. +Legacy persistent solvers (e.g. ``gurobi_persistent``) require the objective +to be re-pushed, which ``set_scenario_probabilities`` handles; APPSI solvers +auto-track the change on the next solve. + +A complete, runnable example is ``examples/farmer/farmer_prob_sensitivity.py``, +which sweeps the weight on one farmer scenario and reports how the optimal +first-stage planting decision responds. + +Progressive Hedging and other decomposition +-------------------------------------------- + +The PH family does not bake probabilities into a Pyomo objective; it weights +``xbar``/``W``/``rho`` numerically through ``prob_coeff``, read fresh each +iteration. ``SPBase.set_scenario_probabilities`` updates +``_mpisppy_probability`` on the scenarios and forces ``prob_coeff`` to be +recomputed so the next iteration uses the new weights: + +.. code-block:: python + + ph.set_scenario_probabilities(prob_map, reset_ph_duals=True) + +``reset_ph_duals=True`` (the default) zeroes the PH multipliers ``W``. This +matters when re-solving from a *converged* PH solution: there, every scenario +sits at the same nonanticipative point, so the probability-weighted ``xbar`` +is independent of the weights and PH would otherwise report immediate (false) +convergence at the old solution. Zeroing ``W`` breaks that consensus so PH +re-converges for the new probabilities. + +.. note:: + + ``SPBase.set_scenario_probabilities`` currently supports two-stage + problems. Multistage node probabilities (``ScenarioNode.cond_prob``) are a + later phase and raise ``NotImplementedError``. + +API +--- + +.. automethod:: mpisppy.opt.ef.ExtensiveForm.set_scenario_probabilities + :noindex: + +.. automethod:: mpisppy.spbase.SPBase.set_scenario_probabilities + :noindex: diff --git a/examples/farmer/farmer_prob_sensitivity.py b/examples/farmer/farmer_prob_sensitivity.py new file mode 100644 index 000000000..ab4500ff1 --- /dev/null +++ b/examples/farmer/farmer_prob_sensitivity.py @@ -0,0 +1,99 @@ +############################################################################### +# 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. +############################################################################### +# Probability sensitivity on the Extensive Form (issue #797). +# +# The farmer decides how many acres to plant of each crop before the yields +# are known. Here we fix the scenario set and sweep the probability vector: +# we vary how much weight is placed on the low-yield scenario and watch the +# optimal first-stage planting decision respond. +# +# The point of the example is that the Extensive Form is built ONCE with +# mutable_probability=True. Each new probability vector is pushed in with +# set_scenario_probabilities() and re-solved with reuse_instance=True, so a +# persistent solver keeps its loaded instance across the whole sweep instead +# of rebuilding the model for every vector. +# +# Run, e.g.: +# python farmer_prob_sensitivity.py --solver-name gurobi_persistent +# python farmer_prob_sensitivity.py --solver-name appsi_highs --num-scens 5 + +import argparse + +import farmer + +from mpisppy.opt.ef import ExtensiveForm + + +def make_probability_vector(scenario_names, low_yield_weight): + """Weight the first ("low yield") scenario by low_yield_weight and split + the remaining probability uniformly over the others. Returns a mapping + that sums to 1, as set_scenario_probabilities requires.""" + others = scenario_names[1:] + rest = (1.0 - low_yield_weight) / len(others) + pv = {scenario_names[0]: low_yield_weight} + for sn in others: + pv[sn] = rest + return pv + + +def main(): + parser = argparse.ArgumentParser( + description="Sweep scenario probabilities on a farmer EF without " + "rebuilding the model (issue #797).") + parser.add_argument("--solver-name", default="gurobi_persistent", + help="Pyomo solver name. A persistent solver " + "(e.g. gurobi_persistent, cplex_persistent) or an " + "APPSI solver (e.g. appsi_highs) keeps its loaded " + "instance across the sweep. Default: " + "gurobi_persistent.") + parser.add_argument("--num-scens", type=int, default=3, + help="Number of scenarios. Default: 3.") + args = parser.parse_args() + + scenario_names = farmer.scenario_names_creator(args.num_scens) + scenario_creator_kwargs = {"num_scens": args.num_scens} + + # Build the EF once. mutable_probability=True stores each scenario's + # probability as a mutable Pyomo Param in the objective instead of baking + # it in as a float constant, so the objective can be re-weighted in place. + ef = ExtensiveForm( + options={"solver": args.solver_name}, + all_scenario_names=scenario_names, + scenario_creator=farmer.scenario_creator, + scenario_creator_kwargs=scenario_creator_kwargs, + mutable_probability=True, + ) + + # Sweep the weight on the low-yield scenario. A larger weight makes the + # farmer more worried about a bad harvest and shifts the planting decision. + low_yield_weights = [1.0 / args.num_scens, 0.25, 0.5, 0.75, 0.9] + + print(f"Probability sensitivity for {args.num_scens} farmer scenarios " + f"using {args.solver_name}") + print(f"(weight is the probability placed on {scenario_names[0]})\n") + header = f"{'weight':>8} {'objective':>14} first-stage acreage" + print(header) + print("-" * len(header)) + + for i, w in enumerate(low_yield_weights): + pv = make_probability_vector(scenario_names, w) + ef.set_scenario_probabilities(pv) + # reuse_instance=True after the first solve keeps the persistent + # solver's instance loaded; only the objective coefficients change. + ef.solve_extensive_form(reuse_instance=(i > 0)) + + obj = ef.get_objective_value() + root = ef.get_root_solution() + acreage = " ".join(f"{name.split('[')[-1].rstrip(']')}={val:8.2f}" + for name, val in sorted(root.items())) + print(f"{w:8.3f} {obj:14.2f} {acreage}") + + +if __name__ == "__main__": + main() diff --git a/mpisppy/opt/ef.py b/mpisppy/opt/ef.py index abac84908..75ccd08c0 100644 --- a/mpisppy/opt/ef.py +++ b/mpisppy/opt/ef.py @@ -59,6 +59,7 @@ def __init__( suppress_warnings=False, extensions=None, extension_kwargs=None, + mutable_probability=None, ): """ Create the EF and associated solver. """ super().__init__( @@ -76,6 +77,16 @@ def __init__( self._options_check(required, self.options) self.solver = pyo.SolverFactory(self.options["solver"]) + # When True, scenario probabilities are stored as mutable Pyomo Params + # so they can be updated in place (see set_scenario_probabilities). + # Falls back to the "mutable_probability" option key for convenience. + if mutable_probability is None: + mutable_probability = self.options.get("mutable_probability", False) + self.mutable_probability = mutable_probability + # Tracks whether a persistent solver already has this EF loaded, so + # solve_extensive_form(reuse_instance=True) can skip set_instance. + self._instance_loaded = False + self.extensions = extensions self.extension_kwargs = extension_kwargs @@ -97,25 +108,40 @@ def __init__( raise FileExistsError(f"solver-log-dir={directory} already exists!") self.ef = sputils._create_EF_from_scen_dict(self.local_scenarios, - EF_name=model_name) + EF_name=model_name, + mutable_probability=self.mutable_probability) - def solve_extensive_form(self, solver_options=None, tee=False): + def solve_extensive_form(self, solver_options=None, tee=False, + reuse_instance=False): """ Solve the extensive form. - + Args: solver_options (dict, optional): Dictionary of solver-specific options (e.g. Gurobi options, CPLEX options, etc.). tee (bool, optional): If True, displays solver output. Default False. + reuse_instance (bool, optional): + If True and this EF has already been loaded into a + persistent solver, skip re-loading the instance + (``set_instance``) and re-solve the already-loaded model. + Use this to re-solve cheaply after updating mutable data + such as scenario probabilities (see + set_scenario_probabilities). Ignored for non-persistent + solvers. Default False. Returns: :class:`pyomo.opt.results.results_.SolverResults`: Result returned by the Pyomo solve method. - + """ - if "persistent" in self.options["solver"]: + # Recognizes both legacy PersistentSolver and APPSI / + # pyomo.contrib.solver interfaces (e.g. appsi_highs); see + # sputils.has_persistent_solve_api for why this is not is_persistent. + persistent = sputils.has_persistent_solve_api(self.solver) + if persistent and not (reuse_instance and self._instance_loaded): self.solver.set_instance(self.ef) + self._instance_loaded = True solve_keyword_args = dict() @@ -142,16 +168,68 @@ def solve_extensive_form(self, solver_options=None, tee=False): # this should catch infeasible and unbounded cases return results - if sputils.is_persistent(self.solver): + if persistent: self.solver.load_vars() else: self.ef.solutions.load_from(results) self.first_stage_solution_available = True self.tree_solution_available = True - + return results + def set_scenario_probabilities(self, prob_map): + """ Update scenario probabilities on a mutable-probability EF in place. + + Requires the EF to have been created with + ``mutable_probability=True``. The supplied probabilities replace the + current values; the full set of scenario probabilities must sum to + 1 after the update (this is a full EF). If a persistent solver has + already been loaded, the objective is re-pushed so a subsequent + ``solve_extensive_form(reuse_instance=True)`` re-solves with the new + probabilities without rebuilding the model. + + Args: + prob_map (dict): + Maps scenario name to its new probability. Names not present + are left unchanged. + + Raises: + RuntimeError: + If the EF was not built with ``mutable_probability=True``. + KeyError: + If ``prob_map`` contains an unknown scenario name. + ValueError: + If the resulting probabilities do not sum to 1. + """ + if not self.mutable_probability: + raise RuntimeError( + "set_scenario_probabilities requires the ExtensiveForm to be " + "created with mutable_probability=True.") + prob = self.ef._mpisppy_model.prob + # Validate before applying so a bad call leaves the model unchanged. + # Unmentioned scenarios keep their current probability (partial updates + # are allowed as long as the resulting full vector still sums to 1). + for sname in prob_map: + if sname not in prob: + raise KeyError(f"Unknown scenario name '{sname}' in prob_map.") + resulting = {sn: prob_map.get(sn, pyo.value(prob[sn])) + for sn in self.ef._ef_scenario_names} + total = sum(resulting.values()) + if abs(total - 1.0) > 1e-9: + raise ValueError( + f"scenario probabilities must sum to 1; got {total}.") + for sname, p in prob_map.items(): + prob[sname].value = p + # keep _mpisppy_probability consistent for downstream readers + getattr(self.ef, sname)._mpisppy_probability = p + # Re-push the objective so a persistent solver picks up the new + # coefficients. Required for legacy persistent solvers; harmless (and + # redundant with auto-tracking) for APPSI / pyomo.contrib.solver. + if self._instance_loaded and \ + sputils.has_persistent_solve_api(self.solver): + self.solver.set_objective(self.ef.EF_Obj) + def get_objective_value(self): """ Retrieve the objective value. diff --git a/mpisppy/spbase.py b/mpisppy/spbase.py index fcfcbd0dc..4f67e3fc1 100644 --- a/mpisppy/spbase.py +++ b/mpisppy/spbase.py @@ -411,15 +411,22 @@ def _create_communicators(self): raise RuntimeError(f"For the node {nodename}, the scenario {sname} has the rank {rank} from scenario_names_to_rank and {comm.Get_rank()} from its comm.") - def _compute_unconditional_node_probabilities(self): + def _compute_unconditional_node_probabilities(self, force=False): """ calculates unconditional node probabilities and prob_coeff - and prob0_mask is set to a scalar 1 (used by variable_probability)""" + and prob0_mask is set to a scalar 1 (used by variable_probability) + + Args: + force (bool): if True, rebuild prob_coeff/prob0_mask even when + they already exist. Used by set_scenario_probabilities to + pick up updated _mpisppy_probability values; the default + (False) keeps the compute-once behavior relied on at setup. + """ for k,s in self.local_scenarios.items(): root = s._mpisppy_node_list[0] root.uncond_prob = 1.0 for parent,child in zip(s._mpisppy_node_list[:-1],s._mpisppy_node_list[1:]): child.uncond_prob = parent.uncond_prob * child.cond_prob - if not hasattr(s._mpisppy_data, 'prob_coeff'): + if force or not hasattr(s._mpisppy_data, 'prob_coeff'): s._mpisppy_data.prob_coeff = dict() s._mpisppy_data.prob0_mask = dict() for node in s._mpisppy_node_list: @@ -427,6 +434,112 @@ def _compute_unconditional_node_probabilities(self): s._mpisppy_data.prob0_mask[node.name] = 1.0 # needs to be a float + def set_scenario_probabilities(self, prob_map, check_sum=True, + reset_ph_duals=True): + """ Update scenario probabilities in place for the PH / decomposition + path and refresh the derived ``prob_coeff`` so the next iteration + (xbar, W, rho) uses the new values. + + Unlike the EF, PH does not bake probabilities into a Pyomo + objective, so there is no Param to update and no persistent-solver + objective to re-push: the aggregation code reads + ``_mpisppy_data.prob_coeff`` fresh each iteration. This method + updates ``_mpisppy_probability`` on the local scenarios named in + ``prob_map`` and forces ``prob_coeff`` to be recomputed. + + Two-stage only for now (multistage node ``cond_prob`` handling is a + later phase; a multistage tree raises). Any variable-probability + overrides are re-applied after the refresh. + + Warm-start caveat (why ``reset_ph_duals`` defaults to True): at a + converged PH solution every scenario sits at the same + (nonanticipative) point, so the probability-weighted ``xbar`` equals + that point *regardless of the weights*. Re-solving from there with + new probabilities but stale W leaves ``xbar`` unmoved and PH reports + immediate (false) convergence at the old solution. Zeroing W breaks + that consensus so the next solve converges to the new problem. Pass + ``reset_ph_duals=False`` only to deliberately keep the current W + (e.g. a small mid-run perturbation before the run has converged). + + Args: + prob_map (dict): + Maps scenario name to its new probability. May be the full + set or a partial update; names not present keep their + current probability. Each rank applies only the scenarios + it owns, so the same full map may be passed on every rank. + check_sum (bool, optional): + If True (default), verify with an MPI reduction that the + resulting probabilities over all scenarios sum to 1 + (option B; see the design doc). Pass False to skip the + collective (e.g. when the caller has already validated). + reset_ph_duals (bool, optional): + If True (default), zero the PH multipliers (``W``) on each + local scenario so a subsequent solve re-converges for the + new probabilities (see the warm-start caveat). No-op for + objects without PH ``W`` terms (e.g. a plain EF). + + Raises: + KeyError: + If ``prob_map`` contains a name not in + ``all_scenario_names``. + NotImplementedError: + If any local scenario has more than the root and a single + leaf node (multistage is a later phase). + ValueError: + If ``check_sum`` and the resulting probabilities do not + sum to 1 within ``E1_tolerance``. + """ + unknown = [sn for sn in prob_map if sn not in self.all_scenario_names] + if unknown: + raise KeyError( + f"set_scenario_probabilities got unknown scenario name(s): " + f"{unknown}") + + for k, s in self.local_scenarios.items(): + # Two-stage only: root plus one leaf. Multistage would need the + # relevant ScenarioNode.cond_prob updated as well (later phase). + if len(s._mpisppy_node_list) > 1: + raise NotImplementedError( + "set_scenario_probabilities currently supports two-stage " + "problems only; multistage node probabilities are a later " + "phase.") + if k in prob_map: + s._mpisppy_probability = prob_map[k] + + # Force recomputation of uncond_prob and prob_coeff; the setup-time + # compute-once short-circuit would otherwise ignore the new values. + self._compute_unconditional_node_probabilities(force=True) + + # Rebuilding prob_coeff above reset any per-variable overrides to the + # scalar node coefficient, so re-apply them. + if self.variable_probability is not None: + self._use_variable_probability_setter() + + # Break a stale consensus so a re-solve tracks the new probabilities + # (see the warm-start caveat above). Guarded because SPBase subclasses + # without PH terms have no W. + if reset_ph_duals: + for s in self.local_scenarios.values(): + if hasattr(s, "_mpisppy_model") and \ + hasattr(s._mpisppy_model, "W"): + for idx in s._mpisppy_model.W: + s._mpisppy_model.W[idx]._value = 0.0 + + if check_sum: + localP = np.zeros(1, dtype='d') + for s in self.local_scenarios.values(): + localP[0] += s._mpisppy_probability + globalP = np.zeros(1, dtype='d') + self.mpicomm.Allreduce([localP, MPI.DOUBLE], + [globalP, MPI.DOUBLE], + op=MPI.SUM) + total = float(globalP[0]) + if abs(total - 1.0) > self.E1_tolerance: + raise ValueError( + f"scenario probabilities must sum to 1; got {total} " + f"(E1_tolerance={self.E1_tolerance}).") + + def _use_variable_probability_setter(self, verbose=False): """ set variable probability unconditional values using a function self.variable_probability that gives us a list of (id(vardata), probability)] diff --git a/mpisppy/tests/test_ef_ph.py b/mpisppy/tests/test_ef_ph.py index 36bc79942..fda4d1bcd 100644 --- a/mpisppy/tests/test_ef_ph.py +++ b/mpisppy/tests/test_ef_ph.py @@ -804,8 +804,214 @@ def test_ph_mult_rho_updater(self): obj2 = round_pos_sig(obj, 2) self.assertEqual(210, obj2) - + # MultRhoUpdater - + + +class Test_mutable_probability(unittest.TestCase): + """ Mutable scenario probabilities on the ExtensiveForm (issue #797). """ + + def setUp(self): + import mpisppy.tests.examples.farmer as farmer + from mpisppy.opt.ef import ExtensiveForm + self.farmer = farmer + self.ExtensiveForm = ExtensiveForm + self.snames = ["scen0", "scen1", "scen2"] + self.sck = {"num_scens": 3, "sense": pyo.minimize} + + def _pv(self, p_bad): + rest = (1.0 - p_bad) / 2.0 + return {"scen0": p_bad, "scen1": rest, "scen2": rest} + + def _make_ef(self, mutable_probability=False): + return self.ExtensiveForm( + options={"solver": solver_name}, + all_scenario_names=self.snames, + scenario_creator=self.farmer.scenario_creator, + scenario_creator_kwargs=self.sck, + mutable_probability=mutable_probability, + ) + + def _rebuild_obj(self, pv): + # baked-in EF at the given probabilities (correctness oracle) + scen_dict = {sn: self.farmer.scenario_creator(sn, **self.sck) + for sn in self.snames} + for sn, scen in scen_dict.items(): + scen._mpisppy_probability = pv[sn] + ef = sputils._create_EF_from_scen_dict(scen_dict, EF_name="oracle") + solver = pyo.SolverFactory(solver_name) + if '_persistent' in solver_name: + solver.set_instance(ef) + solver.solve(ef) + return pyo.value(ef.EF_Obj) + + def test_mutable_prob_builds_param(self): + ef = self._make_ef(mutable_probability=True) + self.assertTrue(ef.mutable_probability) + self.assertTrue(hasattr(ef.ef._mpisppy_model, "prob")) + # defaults to the scenario_creator probabilities (uniform here) + for sn in self.snames: + self.assertAlmostEqual(pyo.value(ef.ef._mpisppy_model.prob[sn]), + 1.0 / 3.0) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_matches_rebuild_across_sweep(self): + ef = self._make_ef(mutable_probability=True) + for i, p_bad in enumerate([0.0, 0.2, 0.5, 0.9]): + pv = self._pv(p_bad) + ef.set_scenario_probabilities(pv) + ef.solve_extensive_form(reuse_instance=(i > 0)) + self.assertAlmostEqual(ef.get_objective_value(), + self._rebuild_obj(pv), places=4) + + @unittest.skipIf(not persistent_available, + "no persistent solver is available") + def test_reuse_instance_loads_once(self): + options = {"solver": persistent_solver_name} + ef = self.ExtensiveForm( + options=options, all_scenario_names=self.snames, + scenario_creator=self.farmer.scenario_creator, + scenario_creator_kwargs=self.sck, mutable_probability=True) + calls = {"n": 0} + orig = ef.solver.set_instance + def counting(*a, **k): + calls["n"] += 1 + return orig(*a, **k) + ef.solver.set_instance = counting + for i, p_bad in enumerate([0.2, 0.5, 0.9]): + ef.set_scenario_probabilities(self._pv(p_bad)) + ef.solve_extensive_form(reuse_instance=(i > 0)) + self.assertEqual(calls["n"], 1) + + def test_guards(self): + ef = self._make_ef(mutable_probability=True) + # non-mutable EF rejects the setter + plain = self._make_ef(mutable_probability=False) + with self.assertRaises(RuntimeError): + plain.set_scenario_probabilities(self._pv(0.2)) + # unknown scenario name + with self.assertRaises(KeyError): + ef.set_scenario_probabilities({"nope": 0.5}) + # probabilities not summing to 1, and the failed call is transactional + before = pyo.value(ef.ef._mpisppy_model.prob["scen0"]) + with self.assertRaises(ValueError): + ef.set_scenario_probabilities({"scen0": before + 0.25}) + self.assertAlmostEqual(pyo.value(ef.ef._mpisppy_model.prob["scen0"]), + before) + + +class Test_mutable_probability_ph(unittest.TestCase): + """ Mutable scenario probabilities on the PH path (issue #797, phase 2). """ + + def setUp(self): + import mpisppy.tests.examples.farmer as farmer + self.farmer = farmer + self.snames = ["scen0", "scen1", "scen2"] + self.sck = {"num_scens": 3, "sense": pyo.minimize} + + def _pv(self, p0): + rest = (1.0 - p0) / 2.0 + return {"scen0": p0, "scen1": rest, "scen2": rest} + + def _denouement(self, *a, **k): + pass + + def _make_ph(self, iters=0): + options = _get_ph_base_options() + options["PHIterLimit"] = iters + return mpisppy.opt.ph.PH( + options, self.snames, self.farmer.scenario_creator, + self._denouement, scenario_creator_kwargs=self.sck) + + def _ef_first_stage(self, pv): + # solved EF at the given probabilities: the nonanticipative first-stage + # DevotedAcreage (independent oracle; the EF path is verified in + # Test_mutable_probability against a rebuild). + scen_dict = {sn: self.farmer.scenario_creator(sn, **self.sck) + for sn in self.snames} + for sn, scen in scen_dict.items(): + scen._mpisppy_probability = pv[sn] + ef = sputils._create_EF_from_scen_dict(scen_dict, EF_name="oracle") + solver = pyo.SolverFactory(solver_name) + if '_persistent' in solver_name: + solver.set_instance(ef) + solver.solve(ef) + block = getattr(ef, self.snames[0]) + return {c: pyo.value(block.DevotedAcreage[c]) + for c in block.DevotedAcreage} + + def _ph_first_stage(self, ph): + # xbar is the consensus first-stage value; read it off any local scenario + s = ph.local_scenarios[self.snames[0]] + return {c: pyo.value(v) for c, v in + zip(s.DevotedAcreage, s.DevotedAcreage.values())} + + def test_prob_coeff_refreshes(self): + # Two-stage: prob_coeff["ROOT"] == _mpisppy_probability after an update. + ph = self._make_ph(iters=0) + for s in ph.local_scenarios.values(): + self.assertAlmostEqual(s._mpisppy_data.prob_coeff["ROOT"], 1.0 / 3.0) + pv = self._pv(0.5) + ph.set_scenario_probabilities(pv) + for sname, s in ph.local_scenarios.items(): + self.assertAlmostEqual(s._mpisppy_probability, pv[sname]) + self.assertAlmostEqual(s._mpisppy_data.prob_coeff["ROOT"], pv[sname]) + + def test_guards(self): + ph = self._make_ph(iters=0) + with self.assertRaises(KeyError): + ph.set_scenario_probabilities({"nope": 0.5}) + with self.assertRaises(ValueError): + ph.set_scenario_probabilities({"scen0": 0.9}) # sum != 1 + + @unittest.skipIf(not solver_available, "no solver is available") + def test_ph_matches_ef_oracle(self): + # A fresh PH at each probability vector should converge to the EF + # solution at that vector, whether the vector was set at construction + # (uniform) or via set_scenario_probabilities before the first solve. + for p0 in (1.0 / 3.0, 0.6): + pv = self._pv(p0) + ph = self._make_ph(iters=100) + if abs(p0 - 1.0 / 3.0) > 1e-12: + ph.set_scenario_probabilities(pv) + ph.ph_main() + got = self._ph_first_stage(ph) + want = self._ef_first_stage(pv) + for c in want: + self.assertAlmostEqual(got[c], want[c], places=2) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_ph_reuse_after_prob_change(self): + # In-place probability change on an already-solved PH object: re-solving + # must track the new probabilities (the probability-sensitivity use + # case). The default reset_ph_duals=True breaks the stale consensus. + ph = self._make_ph(iters=100) + ph.ph_main() + pv = self._pv(0.6) + ph.set_scenario_probabilities(pv) + ph.ph_main() + got = self._ph_first_stage(ph) + want = self._ef_first_stage(pv) + for c in want: + self.assertAlmostEqual(got[c], want[c], places=2) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_ph_reuse_without_dual_reset_stays_stuck(self): + # Documents the warm-start caveat: keeping stale W (reset_ph_duals=False) + # after convergence leaves PH at the old consensus, wrong for new probs. + ph = self._make_ph(iters=100) + ph.ph_main() + uniform = self._ph_first_stage(ph) + pv = self._pv(0.6) + ph.set_scenario_probabilities(pv, reset_ph_duals=False) + ph.ph_main() + stuck = self._ph_first_stage(ph) + want = self._ef_first_stage(pv) + # still at the uniform solution, not the (different) skewed optimum + self.assertNotAlmostEqual(stuck["WHEAT0"], want["WHEAT0"], places=2) + for c in uniform: + self.assertAlmostEqual(stuck[c], uniform[c], places=2) + + if __name__ == '__main__': unittest.main() diff --git a/mpisppy/utils/sputils.py b/mpisppy/utils/sputils.py index 5c9cad942..08c75f14d 100644 --- a/mpisppy/utils/sputils.py +++ b/mpisppy/utils/sputils.py @@ -296,7 +296,8 @@ def create_EF(scenario_names, scenario_creator, scenario_creator_kwargs=None, return EF_instance def _create_EF_from_scen_dict(scen_dict, EF_name=None, - nonant_for_fixed_vars=True): + nonant_for_fixed_vars=True, + mutable_probability=False): """ Create a ConcreteModel of the extensive form from a scenario dictionary. @@ -329,6 +330,15 @@ def _create_EF_from_scen_dict(scen_dict, EF_name=None, enforced non-anticipativity for non-fixed vars, which is not always desirable in the context of bundling. This allows for more fine-grained control. + + If mutable_probability is True, the scenario probabilities are + represented as a mutable Pyomo Param (``_mpisppy_model.prob``, + indexed by scenario name) instead of being folded into the + objective as float constants, so they can be updated after the EF + is built without rebuilding it (see ExtensiveForm. + set_scenario_probabilities). This requires the probabilities to + sum to 1 (as a full EF does) and is therefore NOT supported for + bundles, whose member probabilities sum to less than 1. """ is_min, clear = _models_have_same_sense(scen_dict) if (not clear): @@ -346,6 +356,18 @@ def _create_EF_from_scen_dict(scen_dict, EF_name=None, EF_instance._ef_scenario_names = [] EF_instance._mpisppy_probability = 0 + if mutable_probability: + # Probabilities become mutable Params so they can be updated in place + # (option B of the design: require a normalized full EF, no divisor). + try: + prob_init = {sname: float(scen._mpisppy_probability) + for sname, scen in scen_dict.items()} + except (AttributeError, TypeError, ValueError) as e: + raise ValueError("mutable_probability requires every scenario to " + "have a numeric _mpisppy_probability.") from e + EF_instance._mpisppy_model.prob = pyo.Param( + list(scen_dict.keys()), mutable=True, within=pyo.NonNegativeReals, + initialize=prob_init) for (sname, scenario_instance) in scen_dict.items(): EF_instance.add_component(sname, scenario_instance) EF_instance._ef_scenario_names.append(sname) @@ -353,15 +375,29 @@ def _create_EF_from_scen_dict(scen_dict, EF_name=None, scenario_objs = deact_objs(scenario_instance) obj_func = scenario_objs[0] # Select the first objective try: - EF_instance.EF_Obj.expr += scenario_instance._mpisppy_probability * obj_func.expr + if mutable_probability: + EF_instance.EF_Obj.expr += \ + EF_instance._mpisppy_model.prob[sname] * obj_func.expr + else: + EF_instance.EF_Obj.expr += scenario_instance._mpisppy_probability * obj_func.expr EF_instance._mpisppy_probability += scenario_instance._mpisppy_probability except AttributeError as e: raise AttributeError("Scenario " + sname + " has no specified " "probability. Specify a value for the attribute " " _mpisppy_probability and try again.") from e - # Normalization does nothing when solving the full EF, but is required for - # appropriate scaling of EFs used as bundles. - EF_instance.EF_Obj.expr /= EF_instance._mpisppy_probability + if mutable_probability: + # Require a normalized full EF; a bundle (sum < 1) would land here and + # is intentionally rejected, since the divisor it needs cannot be a + # baked constant when probabilities are mutable. + if abs(EF_instance._mpisppy_probability - 1.0) > 1e-9: + raise RuntimeError( + "mutable_probability requires scenario probabilities summing " + f"to 1; got {EF_instance._mpisppy_probability}. " + "(mutable_probability is not supported for scenario bundles.)") + else: + # Normalization does nothing when solving the full EF, but is required + # for appropriate scaling of EFs used as bundles. + EF_instance.EF_Obj.expr /= EF_instance._mpisppy_probability # For each node in the scenario tree, we need to collect the # nonanticipative vars and create the constraints for them, @@ -513,6 +549,27 @@ def is_persistent(solver): pyo.pyomo.solvers.plugins.solvers.persistent_solver.PersistentSolver) +def has_persistent_solve_api(solver): + """Return True if the solver object supports the persistent EF workflow: + ``set_instance``, ``set_objective``, and ``load_vars``. + + This recognizes both the legacy ``PersistentSolver`` interface and the + APPSI / ``pyomo.contrib.solver`` interfaces (e.g. ``appsi_highs``), which + are not subclasses of the legacy base class and so are *not* reported by + :func:`is_persistent`. + + It is deliberately kept separate from :func:`is_persistent`. Many PH/FWPH + call sites gate ``update_var``/``add_var``/``add_constraint`` on + ``is_persistent``; the APPSI legacy wrapper does not expose those methods, + so broadening ``is_persistent`` itself would break those paths. The EF + solve path only needs the three methods checked here. + """ + if is_persistent(solver): + return True + return all(hasattr(solver, name) + for name in ("set_instance", "set_objective", "load_vars")) + + def solver_quadratic_objective_capability(solver_plugin): """Tri-state probe of whether a solver can handle a quadratic objective.