Presolve Refactor and Distributed OBBT

doc/src/generic_cylinders.rst

@@ -118,3 +118,40 @@ warmstart-subproblems
 Loosely speaking, this option causes subproblem solves to be given the
 previous iteration solution as a warm-start. This is particularly important
 when using an option to lineraize proximal terms.
+
+presolve (FBBT and OBBT)
+-------------------------
+
+The ``presolve`` option enables variable bounds tightening steps before solving,
+which can improve solver performance and numerical stability. The presolve process
+includes two main techniques:
+
+**Feasibility-Based Bounds Tightening (FBBT)**: Uses constraint propagation to
+tighten variable bounds based on the feasible region defined by the constraints.
+This is a fast, lightweight technique that is always applied when presolve is enabled.
+
+**Optimization-Based Bounds Tightening (OBBT)**: Uses optimization to 
+find the tightest possible bounds by solving auxiliary min and max optimization
+problems. This is more computationally expensive but can achieve significantly
+tighter bounds.
+
+To enable basic presolve with FBBT only:
+
+.. code-block:: bash
+
+   python -m mpisppy.generic_cylinders farmer.farmer --num-scens 3 --presolve
+
+
+To enable OBBT in addition to FBBT:
+
+.. code-block:: bash
+
+   python -m mpisppy.generic_cylinders farmer.farmer --num-scens 3 --presolve --obbt
+
+.. Warning::
+   OBBT may not be compatible with all solver interfaces, particularly persistent
+   solvers. If you encounter issues, try using a different solver for OBBT operations.
+
+The OBBT implementation uses Pyomo's ``obbt_analysis`` function and supports
+various configuration options through the ``obbt_options`` dictionary. For
+more details on the implementation, see ``mpisppy.opt.presolve.py``.

examples/netdes/netdes_cylinders.py

@@ -30,6 +30,7 @@ def _parse_args():
     cfg.slammax_args()
     cfg.cross_scenario_cuts_args()
     cfg.reduced_costs_args()
+    cfg.presolve_args()
     cfg.add_to_config("instance_name",
                         description="netdes instance name (e.g., network-10-20-L-01)",
                         domain=str,

examples/run_all.py

@@ -273,10 +273,12 @@ def do_one_mmw(dirname, modname, runefstring, npyfile, mmwargstring):
        f"--num-scens 3 --crops-multiplier=1  --EF-solver-name={solver_name} "
        "--BPL-c0 25 --BPL-eps 100 --confidence-level 0.95 --BM-vs-BPL BPL")
 
+# NOTE: Pyomo OBBT does not support persistent solvers as of Aug 2025
+direct_solver_name = solver_name.replace("_persistent", "_direct") if "_persistent" in solver_name else solver_name
 do_one("netdes", "netdes_cylinders.py", 4,
        "--max-iterations=3 --instance-name=network-10-20-L-01 "
-       "--solver-name={} --rel-gap=0.0 --default-rho=10000 --presolve "
-       "--slammax --subgradient-hub --xhatshuffle --cross-scenario-cuts --max-solver-threads=2".format(solver_name))
+       "--solver-name={} --rel-gap=0.0 --default-rho=10000 --presolve --obbt --obbt-solver={} "
+       "--slammax --subgradient-hub --xhatshuffle --cross-scenario-cuts --max-solver-threads=2".format(solver_name, direct_solver_name))
 
 # sizes is slow for xpress so try linearizing the proximal term.
 do_one("sizes",

mpisppy/opt/presolve.py

@@ -23,8 +23,10 @@
 
 from pyomo.common.errors import InfeasibleConstraintException, DeferredImportError
 from pyomo.contrib.appsi.fbbt import IntervalTightener
+from pyomo.contrib.alternative_solutions.obbt import obbt_analysis
+import pyomo.environ as pyo
 
-from mpisppy import MPI
+from mpisppy import MPI, global_toc
 
 _INF = 1e100
 
@@ -59,47 +61,36 @@ def opt(self, value):
             raise RuntimeError("SPPresolve.opt should only be set once")
 
 
-class SPIntervalTightener(_SPPresolver):
-    """Interval Tightener (feasibility-based bounds tightening)
-    TODO: enable options
+class _SPIntervalTightenerBase(_SPPresolver):
+    """Interval Tightener Base Class (feasibility-based bounds tightening OR optimization-based bounds tightening)
     """
 
     def __init__(self, spbase, verbose=False):
         super().__init__(spbase, verbose)
 
-        self.subproblem_tighteners = {}
-        for k, s in self.opt.local_subproblems.items():
-            try:
-                self.subproblem_tighteners[k] = it = IntervalTightener()
-            except DeferredImportError as e:
-                # User may not have extension built
-                raise ImportError(f"presolve needs the APPSI extensions for pyomo: {e}")
-
-            # ideally, we'd be able to share the `_cmodel`
-            # here between interfaces, etc.
-            try:
-                it.set_instance(s)
-            except Exception as e:
-                # TODO: IntervalTightener won't handle
-                # every Pyomo model smoothly, see:
-                # https://github.com/Pyomo/pyomo/issues/3002
-                # https://github.com/Pyomo/pyomo/issues/3184
-                # https://github.com/Pyomo/pyomo/issues/1864#issuecomment-1989164335
-                raise Exception(
-                    f"Issue with IntervalTightener; cannot apply presolve to this problem. Error: {e}"
-                )
-
         self._lower_bound_cache = {}
         self._upper_bound_cache = {}
 
+        self.printed_warning = False
+
+    @abc.abstractmethod
+    def _tighten_intervals(self) -> bool:
+        """ return True if we should stop, otherwise return False """
+        raise NotImplementedError
+
+    @property
+    def subproblem_tighteners(self):
+        """ return the dictionary of subproblem tighteners """
+        raise NotImplementedError
+
     def presolve(self):
-        """Run the interval tightener (FBBT):
-        1. FBBT on each subproblem
+        """Run the interval tightener (FBBT/OBBT):
+        1. FBBT/OBBT on each subproblem
         2. Narrow bounds on the nonants across all subproblems
         3. If the bounds are updated, go to (1)
         """
+        global_toc(f"Start {self.__class__.__name__}")
 
-        printed_warning = False
         update = False
 
         same_nonant_bounds, global_lower_bounds, global_upper_bounds = (
@@ -125,17 +116,17 @@ def presolve(self):
                                 msg = f"Nonant {var.name} has lower bound greater than upper bound; lb: {lb}, ub: {ub}"
                                 raise InfeasibleConstraintException(msg)
                             if (
-                                (not printed_warning or self.verbose)
+                                (not self.printed_warning or self.verbose)
                                 and (lb, ub) != var.bounds
                                 and (node_comm.Get_rank() == 0)
                             ):
-                                if not printed_warning:
-                                    msg = "WARNING: SPIntervalTightener found different bounds on nonanticipative variables from different scenarios."
+                                if not self.printed_warning:
+                                    msg = f"WARNING: {self.__class__.__name__} found different bounds on nonanticipative variables from different scenarios."
                                     if self.verbose:
                                         print(msg + " See below.")
                                     else:
                                         print(msg + " Use verbose=True to see details.")
-                                    printed_warning = True
+                                    self.printed_warning = True
                                 if self.verbose:
                                     print(
                                         f"Tightening bounds on nonant {var.name} in scenario {k} from {var.bounds} to {(lb, ub)} based on global bound information."
@@ -146,33 +137,27 @@ def presolve(self):
                                 ub = None
                             var.bounds = (lb, ub)
 
-            # Now do FBBT
-            big_iters = 0.0
-            for k, it in self.subproblem_tighteners.items():
-                n_iters = it.perform_fbbt(self.opt.local_subproblems[k])
-                # get the number of constraints after we do
-                # FBBT so we get any updates on the subproblem
-                big_iters = max(big_iters, n_iters / len(it._cmodel.constraints))
-
-            update_this_pass = big_iters > 1.0
-            update_this_pass = self.opt.allreduce_or(update_this_pass)
+            # Now do FBBT/OBBT 
+            stop = self._tighten_intervals()
 
-            if not update_this_pass:
+            if stop:
                 break
-
-            update = True
-
+
             same_nonant_bounds, global_lower_bounds, global_upper_bounds = (
                 self._compare_nonant_bounds()
             )
 
             if same_nonant_bounds:
                 # The nonant bounds did not change, so it is unlikely
-                # that further rounds of FBBT will do any good.
+                # that further rounds of FBBT/OBBT will do any good.
                 break
 
+            update = True
+
+        bounds_tightened = 0
         if update:
-            self._print_bound_movement()
+            bounds_tightened = self._print_bound_movement()
+        global_toc(f"{self.__class__.__name__} tightend {bounds_tightened} nonanticipative variable bound(s).")
 
         return update
 
@@ -261,6 +246,21 @@ def _compare_nonant_bounds(self):
 
         return same_nonant_bounds, global_lower_bounds, global_upper_bounds
 
+    def _check_bounds(self, model):
+        for var in model.component_data_objects(pyo.Var, active=True):
+                lb = var.lb
+                ub = var.ub
+
+                if lb is not None and ub is not None:
+                    if lb > ub:
+                        delta = lb - ub
+                        lb = lb - delta
+                        ub = ub + delta
+
+                        # update revised bounds
+                        var.lb = lb
+                        var.ub = ub
+
     def _print_bound_movement(self):
         lower_bound_movement = {}
         upper_bound_movement = {}
@@ -369,12 +369,107 @@ def _print_bound_movement(self):
                             )
 
                 printed_nodes.add(ndn)
+        return bounds_tightened
+
+
+class SPFBBT(_SPIntervalTightenerBase):
 
-        if (
-            bounds_tightened > 0
-            and self.opt.cylinder_rank == 0
-        ):
-            print(f"SPIntervalTightener tightend {bounds_tightened} bounds.")
+    def __init__(self, spbase, verbose=False):
+
+        super().__init__(spbase, verbose)
+
+        self._subproblem_tighteners = {}
+        for k, s in self.opt.local_subproblems.items():
+            try:
+                self._subproblem_tighteners[k] = it = IntervalTightener()
+            except DeferredImportError as e:
+                # User may not have extension built
+                raise ImportError(f"presolve needs the APPSI extensions for pyomo: {e}")
+
+            # ideally, we'd be able to share the `_cmodel`
+            # here between interfaces, etc.
+            try:
+                it.set_instance(s)
+            except Exception as e:
+                # TODO: IntervalTightener won't handle
+                # every Pyomo model smoothly, see:
+                # https://github.com/Pyomo/pyomo/issues/3002
+                # https://github.com/Pyomo/pyomo/issues/3184
+                # https://github.com/Pyomo/pyomo/issues/1864#issuecomment-1989164335
+                raise Exception(
+                    f"Issue with IntervalTightener; cannot apply presolve to this problem. Error: {e}"
+                )
+
+    @property
+    def subproblem_tighteners(self):
+        return self._subproblem_tighteners
+
+    def _tighten_intervals(self):
+        big_iters = 0.0
+        for k, it in self.subproblem_tighteners.items():
+            n_iters = it.perform_fbbt(self.opt.local_subproblems[k])
+            self._check_bounds(self.opt.local_subproblems[k])
+            # get the number of constraints after we do
+            # FBBT so we get any updates on the subproblem
+            big_iters = max(big_iters, n_iters / len(it._cmodel.constraints))
+        update_this_pass = big_iters > 1.0
+        update_this_pass = self.opt.allreduce_or(update_this_pass)
+
+        return not update_this_pass
+
+
+class SPOBBT(_SPIntervalTightenerBase):
+
+    def __init__(self, spbase, fbbt, verbose=False, obbt_options=None):
+        super().__init__(spbase, verbose)
+
+        self.fbbt = weakref.ref(fbbt)
+        # we expect this in OBBT, so don't warn
+        self.printed_warning = True
+
+        if obbt_options is None:
+            obbt_options = {}
+
+        # TODO: this will create the solver twice, once here and again
+        #       before the first solve ... need to resolve this issue
+        self.solver = obbt_options.get("solver_name", self.opt.options["solver_name"])
+        self.solver_options = obbt_options.get("solver_options")
+        if self.solver_options is None:
+            self.solver_options = {}
+        self.nonant_variables = obbt_options.get("nonant_variables_only", True)
+
+    @property
+    def subproblem_tighteners(self):
+        return self.fbbt().subproblem_tighteners
+
+    def _tighten_intervals(self):
+
+        for k, it in self.fbbt().subproblem_tighteners.items():
+            s = self.opt.local_subproblems[k]
+            it.perform_fbbt(self.opt.local_subproblems[k])
+            self._check_bounds(self.opt.local_subproblems[k])
+
+            if self.nonant_variables:
+                bounds = obbt_analysis(s, variables=s._mpisppy_data.nonant_indices.values(), solver=self.solver,
+                                            solver_options=self.solver_options, warmstart=False)
+
+                for var in s.component_data_objects(pyo.Var, active=True):
+                    if var in bounds.keys():
+                        lb, ub = bounds[var]
+                        var.bounds = (lb, ub)
+
+                it.perform_fbbt(self.opt.local_subproblems[k])
+                self._check_bounds(self.opt.local_subproblems[k])
+
+            else:
+                bounds = obbt_analysis(s, variables=None, solver=self.solver, solver_options=self.solver_options)
+
+                for var in s.component_data_objects(pyo.Var, active=True):
+                    if var.name in bounds.keys():
+                        lb, ub = bounds[var.name]
+                        var.bounds = (lb, ub)
+
+        return False
 
 
 def _lb_generator(var_iterable):
@@ -400,10 +495,19 @@ class SPPresolve(_SPPresolver):
         spbase (SPBase): an SPBase object
     """
 
-    def __init__(self, spbase, verbose=False):
+    def __init__(self, spbase, presolve_options=None, verbose=False):
         super().__init__(spbase, verbose)
+        if presolve_options is None:
+            presolve_options = {}
 
-        self.interval_tightener = SPIntervalTightener(spbase, verbose)
+        self.fbbt = SPFBBT(spbase, verbose)
+        if presolve_options.get("obbt", False):
+            self.obbt = SPOBBT(spbase, self.fbbt, verbose, presolve_options.get("obbt_options", None))
+        else:
+            self.obbt = None
 
     def presolve(self):
-        return self.interval_tightener.presolve()
+        changes_made = self.fbbt.presolve()
+        if self.obbt is not None:
+            changes_made = self.obbt.presolve() or changes_made
+        return changes_made

mpisppy/spopt.py

@@ -69,7 +69,7 @@ def __init__(
             #       models, it is imperative we allow this
             #       object to get garbage collected to
             #       free the memory the C++ model uses.
-            SPPresolve(self).presolve()
+            SPPresolve(self, options.get("presolve_options",None)).presolve()
         self._create_fixed_nonant_cache()
         self.current_solver_options = None
         self.extensions = extensions

mpisppy/utils/cfg_vanilla.py

@@ -92,6 +92,15 @@ def shared_options(cfg):
         shoptions["rc_bound_tol"] = cfg.rc_bound_tol
     if _hasit(cfg, "solver_log_dir"):
         shoptions["solver_log_dir"] = cfg.solver_log_dir
+    if _hasit(cfg, "obbt"):
+        shoptions["presolve_options"] = {
+            "obbt" : cfg.obbt,
+            "obbt_options" : {
+                "nonant_variables_only" : not cfg.full_obbt,
+                "solver_name": cfg.solver_name if cfg.obbt_solver is None else cfg.obbt_solver,
+                "solver_options" : sputils.option_string_to_dict(cfg.obbt_solver_options)
+            },
+        }
 
     return shoptions