first working version of obbt in MindtPy

pyomo/contrib/coramin/relaxations/hessian.py

@@ -3,6 +3,18 @@
 from pyomo.common.collections import ComponentSet
 import enum
 import pyomo.environ as pe
+from pyomo.core import (
+    minimize,
+    maximize,
+    ComponentMap,
+    ConcreteModel,
+    Var,
+    Objective,
+    Reference,
+    Constraint,
+    value,
+)
+from pyomo.opt.base import SolverFactory
 from pyomo.core.expr.numvalue import is_fixed
 import math
 from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr
@@ -56,9 +68,7 @@ def __init__(
         self._constant_hessian_max_eig = None
         self._expr = expr
         self._var_list = list(identify_variables(expr=expr, include_fixed=False))
-        self._ndx_map = pe.ComponentMap(
-            (v, ndx) for ndx, v in enumerate(self._var_list)
-        )
+        self._ndx_map = ComponentMap((v, ndx) for ndx, v in enumerate(self._var_list))
         self._hessian = self.compute_symbolic_hessian()
         self._eigenvalue_problem: Optional[_BlockData] = None
         self._eigenvalue_relaxation: Optional[_BlockData] = None
@@ -69,7 +79,7 @@ def __init__(
     def variables(self):
         return tuple(self._var_list)
 
-    def formulate_eigenvalue_problem(self, sense=pe.minimize):
+    def formulate_eigenvalue_problem(self, sense=minimize):
         if self._eigenvalue_problem is not None:
             min_eig, max_eig = self.bound_eigenvalues_from_interval_hessian()
             if min_eig > self._eigenvalue_problem.eig.lb:
@@ -79,11 +89,11 @@ def formulate_eigenvalue_problem(self, sense=pe.minimize):
             self._eigenvalue_problem.obj.sense = sense
             return self._eigenvalue_problem
         min_eig, max_eig = self.bound_eigenvalues_from_interval_hessian()
-        m = pe.ConcreteModel()
-        m.eig = pe.Var(bounds=(min_eig, max_eig))
-        m.obj = pe.Objective(expr=m.eig, sense=sense)
+        m = ConcreteModel()
+        m.eig = Var(bounds=(min_eig, max_eig))
+        m.obj = Objective(expr=m.eig, sense=sense)
         for v in self._var_list:
-            m.add_component(v.name, pe.Reference(v))
+            m.add_component(v.name, Reference(v))
 
         n = len(self._var_list)
         np_hess = np.empty((n, n), dtype=object)
@@ -99,11 +109,11 @@ def formulate_eigenvalue_problem(self, sense=pe.minimize):
                     np_hess[ndx1, ndx2] -= m.eig
                 else:
                     np_hess[ndx2, ndx1] = np_hess[ndx1, ndx2]
-        m.det_con = pe.Constraint(expr=_determinant(np_hess) == 0)
+        m.det_con = Constraint(expr=_determinant(np_hess) == 0)
         self._eigenvalue_problem = m
         return m
 
-    def formulate_eigenvalue_relaxation(self, sense=pe.minimize):
+    def formulate_eigenvalue_relaxation(self, sense=minimize):
         if self._eigenvalue_relaxation is not None:
             for orig_v, rel_v in self._orig_to_relaxation_vars.items():
                 orig_lb, orig_ub = orig_v.bounds
@@ -123,9 +133,7 @@ def formulate_eigenvalue_relaxation(self, sense=pe.minimize):
             self._eigenvalue_relaxation.obj.sense = sense
             return self._eigenvalue_relaxation
         m = self.formulate_eigenvalue_problem(sense=sense)
-        all_vars = list(
-            ComponentSet(m.component_data_objects(pe.Var, descend_into=True))
-        )
+        all_vars = list(ComponentSet(m.component_data_objects(Var, descend_into=True)))
         from .auto_relax import relax
 
         relaxation = relax(m)
@@ -168,10 +176,10 @@ def get_maximum_eigenvalue(self):
         elif self.method <= EigenValueBounder.GershgorinWithSimplification:
             res = self.bound_eigenvalues_from_interval_hessian()[1]
         elif self.method == EigenValueBounder.LinearProgram:
-            m = self.formulate_eigenvalue_relaxation(sense=pe.maximize)
+            m = self.formulate_eigenvalue_relaxation(sense=maximize)
             res = self.opt.solve(m).best_objective_bound
         else:
-            m = self.formulate_eigenvalue_problem(sense=pe.maximize)
+            m = self.formulate_eigenvalue_problem(sense=maximize)
             res = self.opt.solve(m).best_objective_bound
         return res
 
@@ -198,13 +206,13 @@ def bound_eigenvalues_from_interval_hessian(self):
 
     def compute_symbolic_hessian(self):
         ders = reverse_sd(self._expr)
-        ders2 = pe.ComponentMap()
+        ders2 = ComponentMap()
         for v in self._var_list:
             ders2[v] = reverse_sd(ders[v])
 
-        res = pe.ComponentMap()
+        res = ComponentMap()
         for v in self._var_list:
-            res[v] = pe.ComponentMap()
+            res[v] = ComponentMap()
 
         n = len(self._var_list)
         for v1 in self._var_list:
@@ -215,7 +223,7 @@ def compute_symbolic_hessian(self):
                     continue
                 der = ders2[v1][v2]
                 if is_fixed(der):
-                    val = pe.value(der)
+                    val = value(der)
                     res[v1][v2] = val
                 else:
                     if self.method >= EigenValueBounder.GershgorinWithSimplification:
@@ -229,16 +237,16 @@ def compute_symbolic_hessian(self):
         for v1, dd in res.items():
             for v2, d2 in dd.items():
                 if is_fixed(d2):
-                    res[v1][v2] = pe.value(d2)
+                    res[v1][v2] = value(d2)
                 else:
                     self._constant_hessian = False
 
         return res
 
     def compute_interval_hessian(self):
-        res = pe.ComponentMap()
+        res = ComponentMap()
         for v in self._var_list:
-            res[v] = pe.ComponentMap()
+            res[v] = ComponentMap()
 
         n = len(self._var_list)
         for ndx1, v1 in enumerate(self._var_list):

pyomo/contrib/mindtpy/algorithm_base_class.py

@@ -82,6 +82,10 @@
     update_solver_timelimit,
     copy_var_list_values,
 )
+from pyomo.contrib import appsi
+from pyomo.contrib import coramin
+
+# coramin, coramin_available = attempt_import('pyomo.contrib.coramin')
 
 single_tree, single_tree_available = attempt_import('pyomo.contrib.mindtpy.single_tree')
 tabu_list, tabu_list_available = attempt_import('pyomo.contrib.mindtpy.tabu_list')
@@ -154,6 +158,7 @@ def __init__(self, **kwds):
         self.mip_start_lazy_oa_cuts = []
         # Whether to load solutions in solve() function
         self.load_solutions = True
+        self.obbt_performed = False
 
     # Support use as a context manager under current solver API
     def __enter__(self):
@@ -753,15 +758,16 @@ def set_up_solve_data(self, model):
             )
             config.use_dual_bound = False
 
+        self.original_model = model
+        self.working_model = model.clone()
+
         if config.use_fbbt:
-            IntervalTightener().perform_fbbt(model)
+            with time_code(self.timing, 'presolve - fbbt'):
+                IntervalTightener().perform_fbbt(self.working_model)
             # fbbt(model)
             # TODO: logging_level is not logging.INFO here
             config.logger.info('Use the fbbt to tighten the bounds of variables')
 
-        self.original_model = model
-        self.working_model = model.clone()
-
         # set up bounds
         if obj.sense == minimize:
             self.primal_bound = float('inf')
@@ -1208,6 +1214,58 @@ def handle_subproblem_optimal(self, fixed_nlp, cb_opt=None, fp=False):
         if self.primal_bound_improved:
             self.best_solution_found = fixed_nlp.clone()
             self.best_solution_found_time = get_main_elapsed_time(self.timing)
+            if config.use_obbt and not self.obbt_performed:
+                with time_code(self.timing, 'obbt'):
+                    vars_to_tighten = ComponentSet()
+                    vars_to_tighten.update(self.mip.MindtPy_utils.variable_list)
+                    vars_to_tighten.update(self.mip.coramin_relaxation.aux_vars[:])
+                    obbt_opt = appsi.solvers.Gurobi()
+
+                    full_space_lb_vars, full_space_ub_vars = (
+                        coramin.domain_reduction.aggressive_filter(
+                            candidate_variables=vars_to_tighten,
+                            relaxation=self.mip.coramin_relaxation,
+                            solver=obbt_opt,
+                            tolerance=1e-4,
+                            objective_bound=self.primal_bound,
+                        )
+                    )
+                    obbt_opt = appsi.solvers.Gurobi()
+                    elapsed = get_main_elapsed_time(self.timing)
+                    remaining = math.ceil(max(config.time_limit - elapsed, 1))
+                    obbt_opt = appsi.solvers.Gurobi()
+                    new_lbs, new_ubs = coramin.domain_reduction.perform_obbt(
+                        model=self.mip.coramin_relaxation,
+                        solver=obbt_opt,
+                        varlist=full_space_lb_vars,
+                        objective_bound=self.primal_bound,
+                        update_bounds=True,
+                        direction='lbs',
+                        time_limit=remaining,
+                        parallel=False,
+                    )
+                    for var in full_space_lb_vars:
+                        if "coramin_relaxation.aux_vars" not in var.name:
+                            self.fixed_nlp.find_component(var.name).setlb(var.lb)
+
+                    obbt_opt = appsi.solvers.Gurobi()
+                    new_lbs, new_ubs = coramin.domain_reduction.perform_obbt(
+                        model=self.mip.coramin_relaxation,
+                        solver=obbt_opt,
+                        varlist=full_space_ub_vars,
+                        objective_bound=self.primal_bound,
+                        update_bounds=True,
+                        direction='ubs',
+                        time_limit=remaining,
+                        parallel=False,
+                    )
+                    for var in full_space_ub_vars:
+                        if (
+                            "coramin_relaxation.aux_vars" not in var.name
+                        ) and not var.is_integer():
+                            self.fixed_nlp.find_component(var.name).setub(var.ub)
+                self.obbt_performed = True
+
         # Add the linear cut
         copy_var_list_values(
             fixed_nlp.MindtPy_utils.variable_list,
@@ -1665,6 +1723,7 @@ def solve_main(self):
         update_solver_timelimit(self.mip_opt, config.mip_solver, self.timing, config)
 
         try:
+            self.mip.coramin_relaxation.relaxations.rel2.pprint(verbose=True)
             main_mip_results = self.mip_opt.solve(
                 self.mip,
                 tee=config.mip_solver_tee,
@@ -1868,12 +1927,19 @@ def handle_main_optimal(self, main_mip, update_bound=True):
         )
 
         if update_bound:
-            self.update_dual_bound(value(MindtPy.mip_obj.expr))
+            mip_obj = next(
+                self.mip.component_data_objects(
+                    ctype=Objective, active=True, descend_into=(Block)
+                )
+            )
+            self.update_dual_bound(value(mip_obj.expr))
+            # self.update_dual_bound(value(MindtPy.mip_obj.expr))
             self.config.logger.info(
                 self.log_formatter.format(
                     self.mip_iter,
                     'MILP',
-                    value(MindtPy.mip_obj.expr),
+                    value(mip_obj.expr),
+                    # value(MindtPy.mip_obj.expr),
                     self.primal_bound,
                     self.dual_bound,
                     self.rel_gap,
@@ -2063,14 +2129,17 @@ def setup_main(self):
         config = self.config
         MindtPy = self.mip.MindtPy_utils
 
+        # TODO: this can be moved to creating the MIP model.
         for c in MindtPy.constraint_list:
+            # if config.use_obbt:
+            #     c.deactivate()
             if c.body.polynomial_degree() not in self.mip_constraint_polynomial_degree:
                 c.deactivate()
 
         MindtPy.cuts.activate()
 
         sign_adjust = 1 if self.objective_sense == minimize else -1
-        MindtPy.del_component('mip_obj')
+
         if config.add_regularization is not None and config.add_no_good_cuts:
             MindtPy.cuts.no_good_cuts.deactivate()
 
@@ -2083,11 +2152,13 @@ def setup_main(self):
                 * sum(v for v in MindtPy.cuts.slack_vars.values())
             )
         main_objective = MindtPy.objective_list[-1]
-        MindtPy.mip_obj = Objective(
-            expr=main_objective.expr
-            + (MindtPy.aug_penalty_expr if config.add_slack else 0),
-            sense=self.objective_sense,
-        )
+        if not config.use_obbt:
+            MindtPy.del_component('mip_obj')
+            MindtPy.mip_obj = Objective(
+                expr=main_objective.expr
+                + (MindtPy.aug_penalty_expr if config.add_slack else 0),
+                sense=self.objective_sense,
+            )
 
         if config.use_dual_bound:
             # Delete previously added dual bound constraint
@@ -2606,6 +2677,48 @@ def initialize_mip_problem(self):
             add_var_bound(self.working_model, config)
 
         self.mip = self.working_model.clone()
+
+        # OBBT
+        if config.use_obbt:
+            with time_code(self.timing, 'presolve - convex relaxation'):
+                self.mip.coramin_relaxation = coramin.relaxations.relax(self.mip)
+                # TODO: only needed for avm cuts
+                # for rel in coramin.relaxations.relaxation_data_objects(
+                #     self.mip.coramin_relaxation
+                # ):
+                #     if rel.relaxation_side == RelaxationSide.BOTH:
+                #         rel._original_constraint = Constraint(
+                #             expr=rel.get_aux_var() == rel.get_rhs_expr()
+                #         )
+                #     elif rel.relaxation_side == RelaxationSide.UNDER:
+                #         rel._original_constraint = Constraint(
+                #             expr=rel.get_aux_var() >= rel.get_rhs_expr()
+                #         )
+                #     else:
+                #         rel._original_constraint = Constraint(
+                #             expr=rel.get_aux_var() <= rel.get_rhs_expr()
+                #         )
+                #     rel._original_constraint.deactivate()
+
+        # TODO: if we use avm cuts for MIP master. We should also apply avm to fixed NLP. Otherwise, the auxiliary variables has no values.
+        # if config.use_avm_cuts:
+        #     # self.mip.MindtPy_utils.constraint_list = []
+        #     self.mip.MindtPy_utils.linear_constraint_list = [
+        #         c for c in self.mip.coramin_relaxation.linear.cons[:]
+        #     ]
+        #     print(self.mip.MindtPy_utils.linear_constraint_list)
+        #     self.mip.MindtPy_utils.nonlinear_constraint_list = [
+        #         rel._original_constraint
+        #         for rel in coramin.relaxations.relaxation_data_objects(
+        #             self.mip.coramin_relaxation
+        #         )
+        #     ]
+        #     print(self.mip.MindtPy_utils.nonlinear_constraint_list)
+        #     self.mip.MindtPy_utils.constraint_list = (
+        #         self.mip.MindtPy_utils.linear_constraint_list
+        #         + self.mip.MindtPy_utils.nonlinear_constraint_list
+        #     )
+
         next(self.mip.component_data_objects(Objective, active=True)).deactivate()
         if hasattr(self.mip, 'dual') and isinstance(self.mip.dual, Suffix):
             self.mip.del_component('dual')
@@ -2778,14 +2891,14 @@ def solve(self, model, **kwds):
         new_logging_level = logging.INFO if config.tee else None
         with lower_logger_level_to(config.logger, new_logging_level):
             self.check_config()
-
+        # TODO: fbbt applied here, might be able to move behind
         self.set_up_solve_data(model)
 
         if config.integer_to_binary:
             TransformationFactory('contrib.integer_to_binary').apply_to(
                 self.working_model
             )
-
+        # constraint_list built here
         self.create_utility_block(self.working_model, 'MindtPy_utils')
         with time_code(self.timing, 'total', is_main_timer=True), lower_logger_level_to(
             config.logger, new_logging_level

pyomo/contrib/mindtpy/config_options.py

@@ -457,11 +457,25 @@ def _add_common_configs(CONFIG):
     CONFIG.declare(
         'use_fbbt',
         ConfigValue(
-            default=False,
+            default=True,
             description='Use fbbt to tighten the feasible region of the problem.',
             domain=bool,
         ),
     )
+    CONFIG.declare(
+        'use_obbt',
+        ConfigValue(
+            default=True,
+            description='Use obbt to tighten the feasible region of the problem.',
+            domain=bool,
+        ),
+    )
+    # CONFIG.declare(
+    #     'use_avm_cuts',
+    #     ConfigValue(
+    #         default=True, description='Use AVM to generate the OA cuts.', domain=bool
+    #     ),
+    # )
     CONFIG.declare(
         'use_dual_bound',
         ConfigValue(