Merge pull request #468 from bknueven/int_relax_and_enforce

Initial implementation of IntegerRelaxThenEnforce

doc/src/extensions.rst

@@ -88,6 +88,21 @@ xhat
 Most of the xhat methods can be used as an extension instead of being used
 as a spoke, when that is desired (e.g. for serial applications).
 
+integer_relax_then_enforce
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This extension is for problems with integer variables. The scenario subproblems
+have the integrality restrictions initially relaxed, and then at a later point
+the subproblem integrality restrictions are re-enabled. The parameter ``ratio``
+(default = 0.5) controls how much of the progressive hedging algorithm, either
+in the iteration or time limit, is used for relaxed progressive hedging iterations.
+The extension will also re-enforce the integrality restrictions if the convergence
+threshold is within 10\%  of the convergence tolerance.
+
+This extension can be especially effective if (1) solving the relaxation
+is much easier than solving the problem with integrality constraints or (2) the
+relaxation is reasonably "tight".
+
 WXBarWriter and WXBarReader
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

examples/run_all.py

@@ -275,9 +275,12 @@ def do_one_mmw(dirname, runefstring, npyfile, mmwargstring):
 do_one("sizes", "sizes_pysp.py", 1, "3 {}".format(solver_name))
 do_one("sslp",
        "sslp_cylinders.py",
-       3,
+       4,
        "--instance-name=sslp_15_45_10 --bundles-per-rank=0 "
-       "--max-iterations=5 --default-rho=1 "
+       "--integer-relax-then-enforce "
+       "--integer-relax-then-enforce-ratio=0.95 "
+       "--lagrangian "
+       "--max-iterations=100 --default-rho=1 "
        "--reduced-costs --rc-fixer --xhatshuffle "
        "--linearize-proximal-terms "
        "--rel-gap=0.0 "

examples/sslp/sslp_cylinders.py

@@ -33,6 +33,7 @@ def _parse_args():
     cfg.subgradient_args()
     cfg.reduced_costs_args()
     cfg.coeff_rho_args()
+    cfg.integer_relax_then_enforce_args()
     cfg.parse_command_line("sslp_cylinders")
     return cfg
 
@@ -91,6 +92,9 @@ def main():
     if cfg.coeff_rho:
         vanilla.add_coeff_rho(hub_dict, cfg)
 
+    if cfg.integer_relax_then_enforce:
+        vanilla.add_integer_relax_then_enforce(hub_dict, cfg)
+
     # FWPH spoke
     if fwph:
         fw_spoke = vanilla.fwph_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs)

mpisppy/cylinders/xhatlooper_bounder.py

@@ -73,7 +73,9 @@ def main(self):
                 logger.debug(f'   *localnonants={str(self.localnonants)}')
 
                 self.opt._put_nonant_cache(self.localnonants)
-                self.opt._restore_nonants()
+                # just for sending the values to other scenarios
+                # so we don't need to tell persistent solvers
+                self.opt._restore_nonants(update_persistent=False)
                 upperbound, srcsname = xhatter.xhat_looper(scen_limit=scen_limit, restore_nonants=False)
 
                 # send a bound to the opt companion

mpisppy/cylinders/xhatshufflelooper_bounder.py

@@ -139,7 +139,9 @@ def _vb(msg):
 
                 # update the caches
                 self.opt._put_nonant_cache(self.localnonants)
-                self.opt._restore_nonants()
+                # just for sending the values to other scenarios
+                # so we don't need to tell persistent solvers
+                self.opt._restore_nonants(update_persistent=False)
 
                 scenario_cycler.begin_epoch()
 

mpisppy/cylinders/xhatspecific_bounder.py

@@ -92,7 +92,10 @@ def main(self):
                 logging.debug('  localnonants={}'.format(str(self.localnonants)))
 
                 self.opt._put_nonant_cache(self.localnonants)  # don't really need all caches
-                self.opt._restore_nonants()
+                # just for sending the values to other scenarios
+                # so we don't need to tell persistent solvers
+                self.opt._restore_nonants(update_persistent=False)
+
                 innerbound = xhatter.xhat_tryit(xhat_scenario_dict, restore_nonants=False)
 
                 self.update_if_improving(innerbound)

mpisppy/cylinders/xhatxbar_bounder.py

@@ -96,7 +96,9 @@ def main(self):
                               format(global_rank))
                 logging.debug('  localnonants={}'.format(str(self.localnonants)))
                 self.opt._put_nonant_cache(self.localnonants)  # don't really need all caches
-                self.opt._restore_nonants()
+                # just for sending the values to other scenarios
+                # so we don't need to tell persistent solvers
+                self.opt._restore_nonants(update_persistent=False)
                 innerbound = xhatter.xhat_tryit(restore_nonants=False)
 
                 self.update_if_improving(innerbound)

mpisppy/extensions/integer_relax_then_enforce.py

@@ -0,0 +1,63 @@
+###############################################################################
+# 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.
+###############################################################################
+
+import time
+import pyomo.environ as pyo
+import mpisppy.extensions.extension
+from mpisppy.utils.sputils import is_persistent
+from mpisppy import global_toc
+
+class IntegerRelaxThenEnforce(mpisppy.extensions.extension.Extension):
+    """ Class for relaxing integer variables, running PH, and then
+        enforcing the integality constraints after some condition.
+    """
+
+    def __init__(self, opt):
+        super().__init__(opt)
+        self.integer_relaxer = pyo.TransformationFactory('core.relax_integer_vars')
+        options = opt.options.get("integer_relax_then_enforce_options", {})
+        # fraction of iterations or time to spend in relaxed mode
+        self.ratio = options.get("ratio", 0.5)
+
+
+    def pre_iter0(self):
+        global_toc(f"{self.__class__.__name__}: relaxing integrality constraints", self.opt.cylinder_rank == 0)
+        for s in self.opt.local_scenarios.values():
+            self.integer_relaxer.apply_to(s) 
+        self._integers_relaxed = True
+
+    def _unrelax_integers(self):
+        for sub in self.opt.local_subproblems.values():
+            for sn in sub.scen_list:
+                s = self.opt.local_scenarios[sn]
+                subproblem_solver = sub._solver_plugin
+                vlist = None
+                if is_persistent(subproblem_solver):
+                    vlist = list(v for v,d in s._relaxed_integer_vars[None].values())
+                self.integer_relaxer.apply_to(s, options={"undo":True}) 
+                if is_persistent(subproblem_solver):
+                    for v in vlist:
+                        subproblem_solver.update_var(v)
+        self._integers_relaxed = False
+
+    def miditer(self):
+        if not self._integers_relaxed:
+            return
+        # time is running out
+        if self.opt.options["time_limit"] is not None and ( time.perf_counter() - self.opt.start_time ) > (self.opt.options["time_limit"] * self.ratio):
+            global_toc(f"{self.__class__.__name__}: enforcing integrality constraints, ran so far for more than {self.opt.options['time_limit']*self.ratio} seconds", self.opt.cylinder_rank == 0)
+            self._unrelax_integers()
+        # iterations are running out
+        if self.opt._PHIter > self.opt.options["PHIterLimit"] * self.ratio:
+            global_toc(f"{self.__class__.__name__}: enforcing integrality constraints, ran so far for {self.opt._PHIter - 1} iterations", self.opt.cylinder_rank == 0)
+            self._unrelax_integers()
+        # nearly converged
+        if self.opt.conv < (self.opt.options["convthresh"] * 1.1):
+            global_toc(f"{self.__class__.__name__}: Enforcing integrality constraints, PH is nearly converged", self.opt.cylinder_rank == 0)
+            self._unrelax_integers()

mpisppy/extensions/xhatbase.py

@@ -192,8 +192,9 @@ def _try_one(self, snamedict, solver_options=None, verbose=False,
                     sputils.reactivate_objs(s)
                 # if you hit infeas, return None
                 if not pyo.check_optimal_termination(results):
-                   self.opt._restore_nonants()
-                   return None
+                    if restore_nonants:
+                        self.opt._restore_nonants()
+                    return None
 
             # feasible xhat found, so finish up 2EF part and return
             if verbose and src_rank == self.cylinder_rank:
@@ -219,9 +220,8 @@ def _try_one(self, snamedict, solver_options=None, verbose=False,
 
         infeasP = self.opt.infeas_prob()
         if infeasP != 0.:
-            # restoring does no harm
-            # if this solution is infeasible
-            self.opt._restore_nonants()
+            if restore_nonants:
+                self.opt._restore_nonants()
             return None
         else:
             if verbose and src_rank == self.cylinder_rank:

mpisppy/extensions/xhatxbar.py

@@ -94,9 +94,8 @@ def _vb(msg):
 
         infeasP = self.opt.infeas_prob()
         if infeasP != 0.:
-            # restoring does no harm
-            # if this solution is infeasible
-            self.opt._restore_nonants()
+            if restore_nonants:
+                self.opt._restore_nonants()
             return None
         else:
             if verbose and self.cylinder_rank == 0:

mpisppy/generic_cylinders.py

@@ -57,6 +57,7 @@ def _parse_args(m):
     cfg.ph_args()
     cfg.aph_args()
     cfg.fixer_args()    
+    cfg.integer_relax_then_enforce_args()
     cfg.gapper_args()    
     cfg.fwph_args()
     cfg.lagrangian_args()
@@ -171,6 +172,9 @@ def _do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, scenario_
     if cfg.rc_fixer:
         vanilla.add_reduced_costs_fixer(hub_dict, cfg)
 
+    if cfg.integer_relax_then_enforce:
+        vanilla.add_integer_relax_then_enforce(hub_dict, cfg)
+
     if cfg.grad_rho:
         ext_classes.append(Gradient_extension)
         hub_dict['opt_kwargs']['options']['gradient_extension_options'] = {'cfg': cfg}        

mpisppy/spopt.py

@@ -59,6 +59,7 @@ def __init__(
         self._subproblem_creation(options.get("verbose", False))
         if options.get("presolve", False):
             self._presolver = SPPresolve(self)
+            self._presolver.presolve()
         else:
             self._presolver = None
         self.current_solver_options = None
@@ -672,7 +673,7 @@ def _fix_root_nonants(self,root_cache):
 
 
 
-    def _restore_nonants(self):
+    def _restore_nonants(self, update_persistent=True):
         """ Restore nonanticipative variables to their original values.
 
         This function works in conjunction with _save_nonants.
@@ -688,7 +689,7 @@ def _restore_nonants(self):
         for k,s in self.local_scenarios.items():
 
             persistent_solver = None
-            if (sputils.is_persistent(s._solver_plugin)):
+            if (update_persistent and sputils.is_persistent(s._solver_plugin)):
                 persistent_solver = s._solver_plugin
 
             for ci, vardata in enumerate(s._mpisppy_data.nonant_indices.values()):
@@ -875,9 +876,6 @@ def _subproblem_creation(self, verbose=False):
 
     def _create_solvers(self, presolve=True):
 
-        if self._presolver is not None and presolve:
-            self._presolver.presolve()
-
         dtiming = ("display_timing" in self.options) and self.options["display_timing"]
         local_sit = [] # Local set instance time for time tracking
         for sname, s in self.local_subproblems.items(): # solver creation

mpisppy/utils/cfg_vanilla.py

@@ -38,6 +38,7 @@
 from mpisppy.cylinders.hub import APHHub
 from mpisppy.extensions.extension import MultiExtension
 from mpisppy.extensions.fixer import Fixer
+from mpisppy.extensions.integer_relax_then_enforce import IntegerRelaxThenEnforce
 from mpisppy.extensions.cross_scen_extension import CrossScenarioExtension
 from mpisppy.extensions.reduced_costs_fixer import ReducedCostsFixer
 from mpisppy.extensions.reduced_costs_rho import ReducedCostsRho
@@ -204,6 +205,13 @@ def add_fixer(hub_dict,
                                               "id_fix_list_fct": cfg.id_fix_list_fct}
     return hub_dict
 
+def add_integer_relax_then_enforce(hub_dict,
+              cfg,
+              ):
+    hub_dict = extension_adder(hub_dict,IntegerRelaxThenEnforce)
+    hub_dict["opt_kwargs"]["options"]["integer_relax_then_enforce_options"] = {"ratio":cfg.integer_relax_then_enforce_ratio}
+    return hub_dict
+
 def add_reduced_costs_rho(hub_dict, cfg):
     hub_dict = extension_adder(hub_dict,ReducedCostsRho)
     hub_dict["opt_kwargs"]["options"]["reduced_costs_rho_options"] = {"multiplier" : cfg.reduced_costs_rho_multiplier, "cfg": cfg}

mpisppy/utils/config.py

@@ -433,6 +433,18 @@ def fixer_args(self):
                            domain=float,
                            default=1e-2)
 
+    def integer_relax_then_enforce_args(self):
+        self.add_to_config('integer_relax_then_enforce',
+                           description="have an integer relax then enforce extensions",
+                           domain=bool,
+                           default=False)
+
+        self.add_to_config('integer_relax_then_enforce_ratio',
+                           description="fraction of time limit or iterations (whichever is sooner) "
+                                       "to spend with relaxed integers",
+                           domain=float,
+                           default=0.5)
+
     def reduced_costs_rho_args(self):
         self.add_to_config("reduced_costs_rho",
                            description="have a ReducedCostsRho extension",