Fix MindtPy bugs

pyomo/contrib/mindtpy/algorithm_base_class.py

@@ -55,7 +55,6 @@
     SuppressInfeasibleWarning,
     _DoNothing,
     lower_logger_level_to,
-    copy_var_list_values,
     get_main_elapsed_time,
     time_code,
 )
@@ -80,6 +79,7 @@
     set_solver_mipgap,
     set_solver_constraint_violation_tolerance,
     update_solver_timelimit,
+    copy_var_list_values,
 )
 
 single_tree, single_tree_available = attempt_import('pyomo.contrib.mindtpy.single_tree')
@@ -108,6 +108,8 @@ def __init__(self, **kwds):
         self.curr_int_sol = []
         self.should_terminate = False
         self.integer_list = []
+        # dictionary {integer solution (list): cuts index (list)}
+        self.int_sol_2_cuts_ind = dict()
 
         # Set up iteration counters
         self.nlp_iter = 0
@@ -126,6 +128,9 @@ def __init__(self, **kwds):
         self.fixed_nlp_log_formatter = (
             '{:1}{:>9}   {:>15}   {:>15g}   {:>12g}   {:>12g}   {:>7.2%}   {:>7.2f}'
         )
+        self.infeasible_fixed_nlp_log_formatter = (
+            '{:1}{:>9}   {:>15}   {:>15}   {:>12g}   {:>12g}   {:>7.2%}   {:>7.2f}'
+        )
         self.log_note_formatter = ' {:>9}   {:>15}   {:>15}'
 
         # Flag indicating whether the solution improved in the past
@@ -805,6 +810,9 @@ def MindtPy_initialization(self):
             self.integer_list.append(self.curr_int_sol)
             fixed_nlp, fixed_nlp_result = self.solve_subproblem()
             self.handle_nlp_subproblem_tc(fixed_nlp, fixed_nlp_result)
+            self.int_sol_2_cuts_ind[self.curr_int_sol] = list(
+                range(1, len(self.mip.MindtPy_utils.cuts.oa_cuts) + 1)
+            )
         elif config.init_strategy == 'FP':
             self.init_rNLP()
             self.fp_loop()
@@ -840,6 +848,29 @@ def init_rNLP(self, add_oa_cuts=True):
             if len(results.solution) > 0:
                 self.rnlp.solutions.load_from(results)
         subprob_terminate_cond = results.solver.termination_condition
+
+        # Sometimes, the NLP solver might be trapped in a infeasible solution if the objective function is nonlinear and partition_obj_nonlinear_terms is True. If this happens, we will use the original objective function instead.
+        if (
+            subprob_terminate_cond == tc.infeasible
+            and config.partition_obj_nonlinear_terms
+            and self.rnlp.MindtPy_utils.objective_list[0].expr.polynomial_degree()
+            not in self.mip_objective_polynomial_degree
+        ):
+            config.logger.info(
+                'Initial relaxed NLP problem is infeasible. This might be related to partition_obj_nonlinear_terms. Try to solve it again without partitioning nonlinear objective function.'
+            )
+            self.rnlp.MindtPy_utils.objective.deactivate()
+            self.rnlp.MindtPy_utils.objective_list[0].activate()
+            results = self.nlp_opt.solve(
+                self.rnlp,
+                tee=config.nlp_solver_tee,
+                load_solutions=self.load_solutions,
+                **nlp_args,
+            )
+            if len(results.solution) > 0:
+                self.rnlp.solutions.load_from(results)
+            subprob_terminate_cond = results.solver.termination_condition
+
         if subprob_terminate_cond in {tc.optimal, tc.feasible, tc.locallyOptimal}:
             main_objective = MindtPy.objective_list[-1]
             if subprob_terminate_cond == tc.optimal:
@@ -880,12 +911,14 @@ def init_rNLP(self, add_oa_cuts=True):
                     self.rnlp.MindtPy_utils.variable_list,
                     self.mip.MindtPy_utils.variable_list,
                     config,
+                    ignore_integrality=True,
                 )
                 if config.init_strategy == 'FP':
                     copy_var_list_values(
                         self.rnlp.MindtPy_utils.variable_list,
                         self.working_model.MindtPy_utils.variable_list,
                         config,
+                        ignore_integrality=True,
                     )
                 self.add_cuts(
                     dual_values=dual_values,
@@ -1219,7 +1252,18 @@ def handle_subproblem_infeasible(self, fixed_nlp, cb_opt=None):
         # TODO try something else? Reinitialize with different initial
         # value?
         config = self.config
-        config.logger.info('NLP subproblem was locally infeasible.')
+        config.logger.info(
+            self.infeasible_fixed_nlp_log_formatter.format(
+                ' ',
+                self.nlp_iter,
+                'Fixed NLP',
+                'Infeasible',
+                self.primal_bound,
+                self.dual_bound,
+                self.rel_gap,
+                get_main_elapsed_time(self.timing),
+            )
+        )
         self.nlp_infeasible_counter += 1
         if config.calculate_dual_at_solution:
             for c in fixed_nlp.MindtPy_utils.constraint_list:
@@ -1241,7 +1285,7 @@ def handle_subproblem_infeasible(self, fixed_nlp, cb_opt=None):
         #         elif var.has_lb() and abs(value(var) - var.lb) < config.absolute_bound_tolerance:
         #             fixed_nlp.ipopt_zU_out[var] = -1
 
-        config.logger.info('Solving feasibility problem')
+        # config.logger.info('Solving feasibility problem')
         feas_subproblem, feas_subproblem_results = self.solve_feasibility_subproblem()
         # TODO: do we really need this?
         if self.should_terminate:
@@ -1375,6 +1419,18 @@ def solve_feasibility_subproblem(self):
         self.handle_feasibility_subproblem_tc(
             feas_soln.solver.termination_condition, MindtPy
         )
+        config.logger.info(
+            self.fixed_nlp_log_formatter.format(
+                ' ',
+                self.nlp_iter,
+                'Feasibility NLP',
+                value(feas_subproblem.MindtPy_utils.feas_obj),
+                self.primal_bound,
+                self.dual_bound,
+                self.rel_gap,
+                get_main_elapsed_time(self.timing),
+            )
+        )
         MindtPy.feas_opt.deactivate()
         for constr in MindtPy.nonlinear_constraint_list:
             constr.activate()
@@ -1601,6 +1657,7 @@ def solve_main(self):
         # setup main problem
         self.setup_main()
         mip_args = self.set_up_mip_solver()
+        update_solver_timelimit(self.mip_opt, config.mip_solver, self.timing, config)
 
         try:
             main_mip_results = self.mip_opt.solve(
@@ -1626,7 +1683,11 @@ def solve_main(self):
                         "No-good cuts are added and GOA algorithm doesn't converge within the time limit. "
                         'No integer solution is found, so the CPLEX solver will report an error status. '
                     )
-            return None, None
+            # Value error will be raised if the MIP problem is unbounded and appsi solver is used when loading solutions. Although the problem is unbounded, a valid result is provided and we do not return None to let the algorithm continue.
+            if 'main_mip_results' in dir():
+                return self.mip, main_mip_results
+            else:
+                return None, None
         if config.solution_pool:
             main_mip_results._solver_model = self.mip_opt._solver_model
             main_mip_results._pyomo_var_to_solver_var_map = (
@@ -1658,6 +1719,7 @@ def solve_fp_main(self):
         config = self.config
         self.setup_fp_main()
         mip_args = self.set_up_mip_solver()
+        update_solver_timelimit(self.mip_opt, config.mip_solver, self.timing, config)
 
         main_mip_results = self.mip_opt.solve(
             self.mip,
@@ -2342,6 +2404,7 @@ def handle_fp_subproblem_optimal(self, fp_nlp):
             fp_nlp.MindtPy_utils.variable_list,
             self.working_model.MindtPy_utils.variable_list,
             self.config,
+            ignore_integrality=True,
         )
         add_orthogonality_cuts(self.working_model, self.mip, self.config)
 
@@ -2585,7 +2648,7 @@ def initialize_subsolvers(self):
             self.nlp_opt, config.nlp_solver, config
         )
         set_solver_constraint_violation_tolerance(
-            self.feasibility_nlp_opt, config.nlp_solver, config
+            self.feasibility_nlp_opt, config.nlp_solver, config, warm_start=False
         )
 
         self.set_appsi_solver_update_config()

pyomo/contrib/mindtpy/single_tree.py

@@ -16,15 +16,12 @@
 from pyomo.repn import generate_standard_repn
 import pyomo.core.expr as EXPR
 from math import copysign
-from pyomo.contrib.mindtpy.util import get_integer_solution
-from pyomo.contrib.gdpopt.util import (
-    copy_var_list_values,
-    get_main_elapsed_time,
-    time_code,
-)
+from pyomo.contrib.mindtpy.util import get_integer_solution, copy_var_list_values
+from pyomo.contrib.gdpopt.util import get_main_elapsed_time, time_code
 from pyomo.opt import TerminationCondition as tc
 from pyomo.core import minimize, value
 from pyomo.core.expr import identify_variables
+import math
 
 cplex, cplex_available = attempt_import('cplex')
 
@@ -35,14 +32,7 @@ class LazyOACallback_cplex(
     """Inherent class in CPLEX to call Lazy callback."""
 
     def copy_lazy_var_list_values(
-        self,
-        opt,
-        from_list,
-        to_list,
-        config,
-        skip_stale=False,
-        skip_fixed=True,
-        ignore_integrality=False,
+        self, opt, from_list, to_list, config, skip_stale=False, skip_fixed=True
     ):
         """This function copies variable values from one list to another.
 
@@ -72,43 +62,48 @@ def copy_lazy_var_list_values(
             if skip_fixed and v_to.is_fixed():
                 continue  # Skip fixed variables.
             v_val = self.get_values(opt._pyomo_var_to_solver_var_map[v_from])
-            try:
-                # We don't want to trigger the reset of the global stale
-                # indicator, so we will set this variable to be "stale",
-                # knowing that set_value will switch it back to "not
-                # stale"
-                v_to.stale = True
-                # NOTE: PEP 2180 changes the var behavior so that domain
-                # / bounds violations no longer generate exceptions (and
-                # instead log warnings).  This means that the following
-                # will always succeed and the ValueError should never be
-                # raised.
-                v_to.set_value(v_val, skip_validation=True)
-            except ValueError as e:
-                # Snap the value to the bounds
-                config.logger.error(e)
-                if (
-                    v_to.has_lb()
-                    and v_val < v_to.lb
-                    and v_to.lb - v_val <= config.variable_tolerance
-                ):
-                    v_to.set_value(v_to.lb, skip_validation=True)
-                elif (
-                    v_to.has_ub()
-                    and v_val > v_to.ub
-                    and v_val - v_to.ub <= config.variable_tolerance
-                ):
-                    v_to.set_value(v_to.ub, skip_validation=True)
-                # ... or the nearest integer
-                elif v_to.is_integer():
-                    rounded_val = int(round(v_val))
-                    if (
-                        ignore_integrality
-                        or abs(v_val - rounded_val) <= config.integer_tolerance
-                    ) and rounded_val in v_to.domain:
-                        v_to.set_value(rounded_val, skip_validation=True)
-                else:
-                    raise
+            rounded_val = int(round(v_val))
+            # We don't want to trigger the reset of the global stale
+            # indicator, so we will set this variable to be "stale",
+            # knowing that set_value will switch it back to "not
+            # stale"
+            v_to.stale = True
+            # NOTE: PEP 2180 changes the var behavior so that domain
+            # / bounds violations no longer generate exceptions (and
+            # instead log warnings).  This means that the following
+            # will always succeed and the ValueError should never be
+            # raised.
+            if (
+                v_val in v_to.domain
+                and not ((v_to.has_lb() and v_val < v_to.lb))
+                and not ((v_to.has_ub() and v_val > v_to.ub))
+            ):
+                v_to.set_value(v_val)
+            # Snap the value to the bounds
+            # TODO: check the performance of
+            # v_to.lb - v_val <= config.variable_tolerance
+            elif (
+                v_to.has_lb()
+                and v_val < v_to.lb
+                # and v_to.lb - v_val <= config.variable_tolerance
+            ):
+                v_to.set_value(v_to.lb)
+            elif (
+                v_to.has_ub()
+                and v_val > v_to.ub
+                # and v_val - v_to.ub <= config.variable_tolerance
+            ):
+                v_to.set_value(v_to.ub)
+            # ... or the nearest integer
+            elif (
+                v_to.is_integer()
+                and math.fabs(v_val - rounded_val) <= config.integer_tolerance
+            ):  # and rounded_val in v_to.domain:
+                v_to.set_value(rounded_val)
+            elif abs(v_val) <= config.zero_tolerance and 0 in v_to.domain:
+                v_to.set_value(0)
+            else:
+                raise ValueError('copy_lazy_var_list_values failed.')
 
     def add_lazy_oa_cuts(
         self,
@@ -708,8 +703,8 @@ def __call__(self):
 
         # Reference: https://www.ibm.com/docs/en/icos/22.1.1?topic=SSSA5P_22.1.1/ilog.odms.cplex.help/refpythoncplex/html/cplex.callbacks.SolutionSource-class.htm
         # Another solution source is user_solution = 118, but it will not be encountered in LazyConstraintCallback.
-        config.logger.debug(
-            "Solution source: %s (111 node_solution, 117 heuristic_solution, 119 mipstart_solution)".format(
+        config.logger.info(
+            "Solution source: {} (111 node_solution, 117 heuristic_solution, 119 mipstart_solution)".format(
                 self.get_solution_source()
             )
         )
@@ -718,6 +713,7 @@ def __call__(self):
         # Lazy constraints separated when processing a MIP start will be discarded after that MIP start has been processed.
         # This means that the callback may have to separate the same constraint again for the next MIP start or for a solution that is found later in the solution process.
         # https://www.ibm.com/docs/en/icos/22.1.1?topic=SSSA5P_22.1.1/ilog.odms.cplex.help/refpythoncplex/html/cplex.callbacks.LazyConstraintCallback-class.htm
+        # For the MINLP3_simple example, all the solutions are obtained from mip_start (solution source). Therefore, it will not go to a branch and bound process.Cause an error output.
         if (
             self.get_solution_source()
             != cplex.callbacks.SolutionSource.mipstart_solution
@@ -910,19 +906,7 @@ def LazyOACallback_gurobi(cb_m, cb_opt, cb_where, mindtpy_solver, config):
             if mindtpy_solver.dual_bound != mindtpy_solver.dual_bound_progress[0]:
                 mindtpy_solver.add_regularization()
 
-        if (
-            abs(mindtpy_solver.primal_bound - mindtpy_solver.dual_bound)
-            <= config.absolute_bound_tolerance
-        ):
-            config.logger.info(
-                'MindtPy exiting on bound convergence. '
-                '|Primal Bound: {} - Dual Bound: {}| <= (absolute tolerance {})  \n'.format(
-                    mindtpy_solver.primal_bound,
-                    mindtpy_solver.dual_bound,
-                    config.absolute_bound_tolerance,
-                )
-            )
-            mindtpy_solver.results.solver.termination_condition = tc.optimal
+        if mindtpy_solver.bounds_converged() or mindtpy_solver.reached_time_limit():
             cb_opt._solver_model.terminate()
             return
 
@@ -953,15 +937,32 @@ def LazyOACallback_gurobi(cb_m, cb_opt, cb_where, mindtpy_solver, config):
                     )
                 return
             elif config.strategy == 'OA':
+                # Refer to the official document of GUROBI.
+                # Your callback should be prepared to cut off solutions that violate any of your lazy constraints, including those that have already been added. Node solutions will usually respect previously added lazy constraints, but not always.
+                # https://www.gurobi.com/documentation/current/refman/cs_cb_addlazy.html
+                # If this happens, MindtPy will look for the index of corresponding cuts, instead of solving the fixed-NLP again.
+                for ind in mindtpy_solver.int_sol_2_cuts_ind[
+                    mindtpy_solver.curr_int_sol
+                ]:
+                    cb_opt.cbLazy(mindtpy_solver.mip.MindtPy_utils.cuts.oa_cuts[ind])
                 return
         else:
             mindtpy_solver.integer_list.append(mindtpy_solver.curr_int_sol)
+            if config.strategy == 'OA':
+                cut_ind = len(mindtpy_solver.mip.MindtPy_utils.cuts.oa_cuts)
 
         # solve subproblem
         # The constraint linearization happens in the handlers
         fixed_nlp, fixed_nlp_result = mindtpy_solver.solve_subproblem()
 
         mindtpy_solver.handle_nlp_subproblem_tc(fixed_nlp, fixed_nlp_result, cb_opt)
+        if config.strategy == 'OA':
+            # store the cut index corresponding to current integer solution.
+            mindtpy_solver.int_sol_2_cuts_ind[mindtpy_solver.curr_int_sol] = list(
+                range(
+                    cut_ind + 1, len(mindtpy_solver.mip.MindtPy_utils.cuts.oa_cuts) + 1
+                )
+            )
 
 
 def handle_lazy_main_feasible_solution_gurobi(cb_m, cb_opt, mindtpy_solver, config):