Relaxed PH spoke

doc/src/extensions.rst

@@ -74,7 +74,7 @@ is an example of its use in ``examples.sizes.sizes_demo.py``
 fixer.py
 ^^^^^^^^
 
-This extension provides methods for fixing variables (usually integers) for
+This extension provides methods for fixing nonanticipative variables (usually integers) for
 which all scenarios have agreed for some number of iterations. There
 is an example of its use in ``examples.sizes.sizes_demo.py`` also
 in ``examples.sizes.uc_ama.py``. The ``uc_ama`` example illustrates
@@ -89,6 +89,36 @@ to be on the ``Config`` object so the amalgamator can find it.
    So if you don't want to fix a variable at iteration zero, provide a
    tolerance, but set all count values to ``None``.
 
+reduced_cost_fixer
+^^^^^^^^^^^^^^^^^^
+
+This extension provides methods for fixing nonanticipative variables based on their expected
+reduced cost as calculated by the ReducedCostSpoke. The aggressiveness of the
+fixing can be controled through the ``zero_rc_tol`` parameter (reduced costs
+with magnitude below this value will be considered 0 and not eligible for fixing)
+and the ``fix_fraction_target`` paramemters, which set a maximum fraction of
+nonanticipative variables to be fixed based on expected reduced costs. These two
+parameters iteract with each other -- the expected reduced costs are sorted by
+magnitude, and if the `fix_fraction_target`` percental is below ``zero_rc_tol``,
+then fewer than ``fix_fraction_target`` variables will be fixed. Further, to
+have a defined expected reduced cost, all nonant variable values *must be* at
+the same bound in the ReducedCostSpoke.
+
+Variables will be unfixed if they no longer meet the expected reduced cost
+criterion for fixing, e.g., the variable's expected reduced cost became too
+low or the variable was not at its bound in every subproblem in the ReducedCostSpoke.
+
+relaxed_ph_fixer
+^^^^^^^^^^^^^^^^
+
+This extension will fix nonanticipative variables at their bound if they are at
+their bound in the RelaxedPHSpoke for that subproblem. It will similarily unfix
+nonanticipative variables which are not at their bounds in the RelaxedPHSpoke.
+Because different nonanticipative variables are fixed in different suproblems,
+it will also unfix nonanticipative variables if their value is *not* at the the current
+consensus solution xbar (because the variable was not fixed in a different subproblem
+and therefore came off its bound).
+
 xhat
 ^^^^
 

doc/src/hubs.rst

@@ -47,6 +47,13 @@ this value (default 1e-1).
 If only the binary terms should be 
 approximated, the option `linearize_binary_proximal_terms` can be used. 
 
+PHNonant
+--------
+This is exactly like the PH hub, except it does not send W values.
+It can be useful when some other cylinder is providing W values
+(e.g., RelaxedPHSpoke).
+
+
 lshaped
 -------
 
@@ -72,6 +79,15 @@ also supplies x and/or W values at every iteration, and is largely based
 on the PH implementation. It utilizes a constant step size rule based on
 `rho` unless modified by an extension.
 
+FWPH
+----
+
+The Frank-Wolfe progressive hedging hub can be used with most spokes
+because it supplies x and/or W values as part of its solution process.
+While FWPH is not known to converge to a primal solution for a SMIP, it
+often discovers excellent incumbent values along the way (when paired with
+an xhat spoke).
+
 Hub Convergers
 --------------
 

doc/src/spokes.rst

@@ -92,8 +92,9 @@ Reduced Costs
 
 The reduced cost spoke is equivalent to the Lagrangian spoke, except that it relaxes all 
 integrality contraints in the subproblems. This enables the computation of reduced costs
-for the first stage variables, which can be used for bound tightening or heuristic fixing
-in the hub.
+for the nonanticipative variables, which can be used for provable bound tightening, which
+is shared to all cylinders. The reduced cost can also be used for heuristic fixing in the
+hub via the reduced_cost_fixer extension.
 
 
 Inner Bounds
@@ -171,4 +172,17 @@ General
 cross scenario
 ^^^^^^^^^^^^^^
 
-Passes cross scenario cuts.
+Computes and passes cross scenario cuts.
+
+relaxed_ph
+^^^^^^^^^^
+
+For S-MIPs, runs progressive hedging on the linear programming relaxation
+of the scenario subproblems. Provides Ws and "relaxed nonants", the former
+of which can be utilized for Lagrangian to compute lower bounds, and the later
+of which can be utilized to inform fixings for the hub via the relaxed_ph_fixer
+extention. This method can be effective when the subproblem linear programming
+relaxation is "strong".
+
+The option ``relaxed_ph_rescale_rho_factor``, which defaults to 1, rescales
+rho for this spoke between iteration 0 and iteration 1.

mpisppy/cylinders/fwph_spoke.py

@@ -15,9 +15,6 @@ class FrankWolfeOuterBound(mpisppy.cylinders.spoke.OuterBoundSpoke):
     def main(self):
         self.opt.fwph_main()
 
-    def is_converged(self):
-        return self.got_kill_signal()
-
     def sync(self):
         # The FWPH spoke can call "sync" before it
         # even starts doing anything, so its possible

mpisppy/cylinders/hub.py

@@ -225,11 +225,21 @@ def sync_bounds(self):
         self.send_boundsout()
 
 
-class PHHub(Hub):
-
-    send_fields = (*Hub.send_fields, Field.NONANT, Field.DUALS)
+class PHNonantHub(Hub):
+    """
+    Like PHHub, but only sends nonants and omits Ws. To be used
+    when another cylinder is supplying Ws (like RelaxedPHSpoke).
+    Could be removed when mpi-sppy supports pointing consuming
+    spokes like Lagrangian to a specific dual (W) buffer.
+    """
+
+    send_fields = (*Hub.send_fields, Field.NONANT, )
     receive_fields = (*Hub.receive_fields,)
 
+    @property
+    def nonant_field(self):
+        return Field.NONANT
+
     def setup_hub(self):
         ## Generate some warnings if nothing is giving bounds
         if not self.receive_field_spcomms[Field.OBJECTIVE_OUTER_BOUND]:
@@ -316,17 +326,27 @@ def send_nonants(self):
         """ Gather nonants and send them to the appropriate spokes
         """
         ci = 0  ## index to self.nonant_send_buffer
-        nonant_send_buffer = self.send_buffers[Field.NONANT]
+        nonant_send_buffer = self.send_buffers[self.nonant_field]
         for k, s in self.opt.local_scenarios.items():
             for xvar in s._mpisppy_data.nonant_indices.values():
                 nonant_send_buffer[ci] = xvar._value
                 ci += 1
         logging.debug("hub is sending X nonants={}".format(nonant_send_buffer))
 
-        self.put_send_buffer(nonant_send_buffer, Field.NONANT)
+        self.put_send_buffer(nonant_send_buffer, self.nonant_field)
 
         return
 
+    def send_ws(self):
+        """ Nonant hub; do not send Ws
+        """
+        pass
+
+
+class PHHub(PHNonantHub):
+    send_fields = (*PHNonantHub.send_fields, Field.DUALS, )
+    receive_fields = (*PHNonantHub.receive_fields,)
+
     def send_ws(self):
         """ Send dual weights to the appropriate spokes
         """
@@ -339,6 +359,7 @@ def send_ws(self):
 
         return
 
+
 class LShapedHub(Hub):
 
     send_fields = (*Hub.send_fields, Field.NONANT,)

mpisppy/cylinders/relaxed_ph_spoke.py

@@ -0,0 +1,57 @@
+###############################################################################
+# 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.
+###############################################################################
+
+from mpisppy.cylinders.spwindow import Field
+from mpisppy.cylinders.spoke import Spoke
+from mpisppy.cylinders.hub import PHHub
+
+import pyomo.environ as pyo
+
+class RelaxedPHSpoke(Spoke, PHHub):
+
+    send_fields = (*Spoke.send_fields, Field.DUALS, Field.RELAXED_NONANT, )
+    receive_fields = (*Spoke.receive_fields, )
+
+    @property
+    def nonant_field(self):
+        return Field.RELAXED_NONANT
+
+    def send_boundsout(self):
+        # overwrite PHHub.send_boundsout (not a hub)
+        return
+
+    def update_rho(self):
+        rho_factor = self.opt.options.get("relaxed_ph_rho_factor", 1.0)
+        if rho_factor == 1.0:
+            return
+        for s in self.opt.local_scenarios.values():
+            for rho in s._mpisppy_model.rho.values():
+                rho._value = rho_factor * rho._value
+
+    def main(self):
+        # relax the integers
+        integer_relaxer = pyo.TransformationFactory("core.relax_integer_vars")
+        for s in self.opt.local_scenarios.values():
+            integer_relaxer.apply_to(s)
+
+        # setup, PH Iter0
+        smoothed = self.options.get('smoothed', 0)
+        attach_prox = True
+        self.opt.PH_Prep(attach_prox=attach_prox, attach_smooth = smoothed)
+        trivial_bound = self.opt.Iter0()
+        if self.opt._can_update_best_bound():
+            self.opt.best_bound_obj_val = trivial_bound
+
+        # update the rho
+        self.update_rho()
+
+        # rest of PH
+        self.opt.iterk_loop()
+
+        return self.opt.conv, None, trivial_bound

mpisppy/cylinders/spoke.py

@@ -29,6 +29,12 @@ def got_kill_signal(self):
         self.get_receive_buffer(shutdown_buf, Field.SHUTDOWN, 0, synchronize=False)
         return self.allreduce_or(shutdown_buf[0] == 1.0)
 
+    def is_converged(self):
+        """ Alias for got_kill_signal; useful for algorithms working as both
+            hub and spoke
+        """
+        return self.got_kill_signal()
+
     @abc.abstractmethod
     def main(self):
         """

mpisppy/cylinders/spwindow.py

@@ -20,6 +20,7 @@ class Field(enum.IntEnum):
     SHUTDOWN=-1000
     NONANT=1
     DUALS=2
+    RELAXED_NONANT=3
     BEST_OBJECTIVE_BOUNDS=100 # Both inner and outer bounds from the hub. Layout: [OUTER INNER ID]
     OBJECTIVE_INNER_BOUND=101
     OBJECTIVE_OUTER_BOUND=102
@@ -42,6 +43,7 @@ class Field(enum.IntEnum):
         Field.SHUTDOWN : 1,
         Field.NONANT : _field_length_components.local_nonant_length,
         Field.DUALS : _field_length_components.local_nonant_length,
+        Field.RELAXED_NONANT : _field_length_components.local_nonant_length,
         Field.BEST_OBJECTIVE_BOUNDS : 2,
         Field.OBJECTIVE_INNER_BOUND : 1,
         Field.OBJECTIVE_OUTER_BOUND : 1,

mpisppy/extensions/relaxed_ph_fixer.py

@@ -0,0 +1,119 @@
+###############################################################################
+# 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.
+###############################################################################
+from mpisppy.extensions.extension import Extension
+from mpisppy.utils.sputils import is_persistent
+
+from mpisppy.cylinders.spwindow import Field
+
+class RelaxedPHFixer(Extension):
+
+    def __init__(self, spobj):
+        super().__init__(spobj)
+
+        ph_options = spobj.options
+        ph_fixer_options = ph_options.get("relaxed_ph_fixer_options", {})
+        self.bound_tol = ph_fixer_options.get("bound_tol", 1e-4)
+        self.verbose = ph_fixer_options.get("verbose", True)
+        self.debug = ph_fixer_options.get("debug", False)
+
+        self._heuristic_fixed_vars = {}
+        self._current_relaxed_nonants = None
+
+    def pre_iter0(self):
+        self._modeler_fixed_nonants = set()
+        self.nonant_length = self.opt.nonant_length
+        for k,s in self.opt.local_scenarios.items():
+            for ndn_i, xvar in s._mpisppy_data.nonant_indices.items():
+                if xvar.fixed:
+                    self._modeler_fixed_nonants.add(ndn_i)
+
+        for k,sub in self.opt.local_subproblems.items():
+            self._heuristic_fixed_vars[k] = 0
+
+    def iter0_post_solver_creation(self):
+        # wait for relaxed iter0:
+        if self.relaxed_nonant_buf.id() == 0:
+            while not self.opt.spcomm.get_receive_buffer(self.relaxed_nonant_buf, Field.RELAXED_NONANT, self.relaxed_ph_spoke_index):
+                continue
+        self.relaxed_ph_fixing(self.relaxed_nonant_buf.value_array(), pre_iter0=True)
+
+    def register_receive_fields(self):
+        spcomm = self.opt.spcomm
+        relaxed_ph_ranks = spcomm.fields_to_ranks[Field.RELAXED_NONANT]
+        assert len(relaxed_ph_ranks) == 1
+        index = relaxed_ph_ranks[0]
+
+        self.relaxed_ph_spoke_index = index
+
+        self.relaxed_nonant_buf = spcomm.register_recv_field(
+            Field.RELAXED_NONANT,
+            self.relaxed_ph_spoke_index,
+        )
+
+        return
+
+    def miditer(self):
+        self.opt.spcomm.get_receive_buffer(
+            self.relaxed_nonant_buf,
+            Field.RELAXED_NONANT,
+            self.relaxed_ph_spoke_index,
+        )
+        self.relaxed_ph_fixing(self.relaxed_nonant_buf.value_array(), pre_iter0=False)
+        return
+
+    def relaxed_ph_fixing(self, relaxed_solution, pre_iter0 = False):
+
+        for k, sub in self.opt.local_subproblems.items():
+            raw_fixed_this_iter = 0
+            persistent_solver = is_persistent(sub._solver_plugin)
+            for sn in sub.scen_list:
+                s = self.opt.local_scenarios[sn]
+                for ci, (ndn_i, xvar) in enumerate(s._mpisppy_data.nonant_indices.items()):
+                    if ndn_i in self._modeler_fixed_nonants:
+                        continue
+                    if xvar in s._mpisppy_data.all_surrogate_nonants:
+                        continue
+                    relaxed_val = relaxed_solution[ci]
+                    xvar_value = xvar._value
+                    update_var = False
+                    if not pre_iter0 and xvar.fixed:
+                        if (relaxed_val - xvar.lb) > self.bound_tol and (xvar.ub - relaxed_val) > self.bound_tol:
+                            xvar.unfix()
+                            update_var = True
+                            raw_fixed_this_iter -= 1
+                            if self.debug and self.opt.cylinder_rank == 0:
+                                print(f"{k}: unfixing var {xvar.name}; {relaxed_val=} is off bounds {(xvar.lb, xvar.ub)}")
+                        # in case somebody else unfixs a variable in another rank...
+                        xb = s._mpisppy_model.xbars[ndn_i]._value
+                        if abs(xb - xvar_value) > self.bound_tol:
+                            xvar.unfix()
+                            update_var = True
+                            raw_fixed_this_iter -= 1
+                            if self.debug and self.opt.cylinder_rank == 0:
+                                print(f"{k}: unfixing var {xvar.name}; xbar {xb} differs from the fixed value {xvar_value}")
+                    elif (relaxed_val - xvar.lb <= self.bound_tol) and (pre_iter0 or (xvar_value - xvar.lb <= self.bound_tol)):
+                        xvar.fix(xvar.lb)
+                        if self.debug and self.opt.cylinder_rank == 0:
+                            print(f"{k}: fixing var {xvar.name} to lb {xvar.lb}; {relaxed_val=}, var value is {xvar_value}")
+                        update_var = True
+                        raw_fixed_this_iter += 1
+                    elif (xvar.ub - relaxed_val <= self.bound_tol) and (pre_iter0 or (xvar.ub - xvar_value <= self.bound_tol)):
+                        xvar.fix(xvar.ub)
+                        if self.debug and self.opt.cylinder_rank == 0:
+                            print(f"{k}: fixing var {xvar.name} to ub {xvar.ub}; {relaxed_val=}, var value is {xvar_value}")
+                        update_var = True
+                        raw_fixed_this_iter += 1
+
+                    if update_var and persistent_solver:
+                        sub._solver_plugin.update_var(xvar)
+
+            # Note: might count incorrectly with bundling?
+            self._heuristic_fixed_vars[k] += raw_fixed_this_iter
+            if self.verbose:
+                print(f"{k}: total unique vars fixed by heuristic: {int(round(self._heuristic_fixed_vars[k]))}/{self.nonant_length}")