Mutable scenario probabilities (issue #797)

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::

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:

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()

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.