Merge pull request #3034 from ZedongPeng/benchmark

Fix MindtPy bugs

pyomo/contrib/mindtpy/MindtPy.py

@@ -50,6 +50,14 @@
 - Add single-tree implementation.
 - Add support for cplex_persistent solver.
 - Fix bug in OA cut expression in cut_generation.py.
+
+24.1.11 changes:
+- fix gurobi single tree termination check bug
+- fix Gurobi single tree cycle handling
+- fix bug in feasibility pump method
+- add special handling for infeasible relaxed NLP
+- update the log format of infeasible fixed NLP subproblems
+- create a new copy_var_list_values function
 """
 
 from pyomo.contrib.mindtpy import __version__

pyomo/contrib/mindtpy/__init__.py

@@ -1 +1 @@
-__version__ = (0, 1, 0)
+__version__ = (1, 0, 0)

pyomo/contrib/mindtpy/cut_generation.py

@@ -210,8 +210,8 @@ def add_oa_cuts_for_grey_box(
                             target_model_grey_box.inputs.values()
                         )
                     )
+                    - (output - value(output))
                 )
-                - (output - value(output))
                 - (slack_var if config.add_slack else 0)
                 <= 0
             )
@@ -271,8 +271,9 @@ def add_ecp_cuts(
                 try:
                     upper_slack = constr.uslack()
                 except (ValueError, OverflowError) as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.error(
-                        str(e) + '\nConstraint {} has caused either a '
+                        'Constraint {} has caused either a '
                         'ValueError or OverflowError.'
                         '\n'.format(constr)
                     )
@@ -300,8 +301,9 @@ def add_ecp_cuts(
                 try:
                     lower_slack = constr.lslack()
                 except (ValueError, OverflowError) as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.error(
-                        str(e) + '\nConstraint {} has caused either a '
+                        'Constraint {} has caused either a '
                         'ValueError or OverflowError.'
                         '\n'.format(constr)
                     )
@@ -424,9 +426,9 @@ def add_affine_cuts(target_model, config, timing):
             try:
                 mc_eqn = mc(constr.body)
             except MCPP_Error as e:
+                config.logger.error(e, exc_info=True)
                 config.logger.error(
-                    '\nSkipping constraint %s due to MCPP error %s'
-                    % (constr.name, str(e))
+                    'Skipping constraint %s due to MCPP error' % (constr.name)
                 )
                 continue  # skip to the next constraint
 

pyomo/contrib/mindtpy/extended_cutting_plane.py

@@ -66,12 +66,6 @@ def MindtPy_iteration_loop(self):
 
             add_ecp_cuts(self.mip, self.jacobians, self.config, self.timing)
 
-        # if add_no_good_cuts is True, the bound obtained in the last iteration is no reliable.
-        # we correct it after the iteration.
-        if (
-            self.config.add_no_good_cuts or self.config.use_tabu_list
-        ) and not self.should_terminate:
-            self.fix_dual_bound(self.last_iter_cuts)
         self.config.logger.info(
             ' ==============================================================================================='
         )
@@ -84,9 +78,12 @@ def check_config(self):
         super().check_config()
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
-        self.jacobians = calc_jacobians(self.mip, self.config)  # preload jacobians
+        self.jacobians = calc_jacobians(
+            self.mip.MindtPy_utils.nonlinear_constraint_list,
+            self.config.differentiate_mode,
+        )  # preload jacobians
         self.mip.MindtPy_utils.cuts.ecp_cuts = ConstraintList(
             doc='Extended Cutting Planes'
         )
@@ -140,7 +137,7 @@ def all_nonlinear_constraint_satisfied(self):
                     lower_slack = nlc.lslack()
                 except (ValueError, OverflowError) as e:
                     # Set lower_slack (upper_slack below) less than -config.ecp_tolerance in this case.
-                    config.logger.error(e)
+                    config.logger.error(e, exc_info=True)
                     lower_slack = -10 * config.ecp_tolerance
                 if lower_slack < -config.ecp_tolerance:
                     config.logger.debug(
@@ -153,7 +150,7 @@ def all_nonlinear_constraint_satisfied(self):
                 try:
                     upper_slack = nlc.uslack()
                 except (ValueError, OverflowError) as e:
-                    config.logger.error(e)
+                    config.logger.error(e, exc_info=True)
                     upper_slack = -10 * config.ecp_tolerance
                 if upper_slack < -config.ecp_tolerance:
                     config.logger.debug(

pyomo/contrib/mindtpy/feasibility_pump.py

@@ -44,9 +44,12 @@ def check_config(self):
         super().check_config()
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
-        self.jacobians = calc_jacobians(self.mip, self.config)  # preload jacobians
+        self.jacobians = calc_jacobians(
+            self.mip.MindtPy_utils.nonlinear_constraint_list,
+            self.config.differentiate_mode,
+        )  # preload jacobians
         self.mip.MindtPy_utils.cuts.oa_cuts = ConstraintList(
             doc='Outer approximation cuts'
         )

pyomo/contrib/mindtpy/global_outer_approximation.py

@@ -67,7 +67,7 @@ def check_config(self):
         super().check_config()
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
         self.mip.MindtPy_utils.cuts.aff_cuts = ConstraintList(doc='Affine cuts')
 
@@ -108,4 +108,5 @@ def deactivate_no_good_cuts_when_fixing_bound(self, no_good_cuts):
             if self.config.use_tabu_list:
                 self.integer_list = self.integer_list[:valid_no_good_cuts_num]
         except KeyError as e:
-            self.config.logger.error(str(e) + '\nDeactivating no-good cuts failed.')
+            self.config.logger.error(e, exc_info=True)
+            self.config.logger.error('Deactivating no-good cuts failed.')

pyomo/contrib/mindtpy/outer_approximation.py

@@ -94,9 +94,12 @@ def check_config(self):
         _MindtPyAlgorithm.check_config(self)
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
-        self.jacobians = calc_jacobians(self.mip, self.config)  # preload jacobians
+        self.jacobians = calc_jacobians(
+            self.mip.MindtPy_utils.nonlinear_constraint_list,
+            self.config.differentiate_mode,
+        )  # preload jacobians
         self.mip.MindtPy_utils.cuts.oa_cuts = ConstraintList(
             doc='Outer approximation cuts'
         )

pyomo/contrib/mindtpy/single_tree.py

@@ -16,12 +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 (
+from pyomo.contrib.mindtpy.util import (
+    get_integer_solution,
     copy_var_list_values,
-    get_main_elapsed_time,
-    time_code,
+    set_var_valid_value,
 )
+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
@@ -35,17 +35,9 @@ 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.
-
         Rounds to Binary/Integer if necessary.
         Sets to zero for NonNegativeReals if necessary.
 
@@ -54,61 +46,29 @@ def copy_lazy_var_list_values(
         opt : SolverFactory
             The cplex_persistent solver.
         from_list : list
-            The variables that provides the values to copy from.
+            The variable list that provides the values to copy from.
         to_list : list
-            The variables that need to set value.
+            The variable list that needs to set value.
         config : ConfigBlock
             The specific configurations for MindtPy.
         skip_stale : bool, optional
             Whether to skip the stale variables, by default False.
         skip_fixed : bool, optional
             Whether to skip the fixed variables, by default True.
-        ignore_integrality : bool, optional
-            Whether to ignore the integrality of integer variables, by default False.
         """
         for v_from, v_to in zip(from_list, to_list):
             if skip_stale and v_from.stale:
                 continue  # Skip stale variable 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
+            set_var_valid_value(
+                v_to,
+                v_val,
+                config.integer_tolerance,
+                config.zero_tolerance,
+                ignore_integrality=False,
+            )
 
     def add_lazy_oa_cuts(
         self,
@@ -309,12 +269,11 @@ def add_lazy_affine_cuts(self, mindtpy_solver, config, opt):
                 try:
                     mc_eqn = mc(constr.body)
                 except MCPP_Error as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.debug(
-                        'Skipping constraint %s due to MCPP error %s'
-                        % (constr.name, str(e))
+                        'Skipping constraint %s due to MCPP error' % (constr.name)
                     )
                     continue  # skip to the next constraint
-                # TODO: check if the value of ccSlope and cvSlope is not Nan or inf. If so, we skip this.
                 ccSlope = mc_eqn.subcc()
                 cvSlope = mc_eqn.subcv()
                 ccStart = mc_eqn.concave()
@@ -705,10 +664,11 @@ def __call__(self):
         main_mip = self.main_mip
         mindtpy_solver = self.mindtpy_solver
 
+        # The lazy constraint callback may be invoked during MIP start processing. In that case get_solution_source returns mip_start_solution.
         # 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()
             )
         )
@@ -717,6 +677,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
@@ -727,6 +688,7 @@ def __call__(self):
             mindtpy_solver.mip_start_lazy_oa_cuts = []
 
         if mindtpy_solver.should_terminate:
+            # TODO: check the performance difference if we don't use self.abort() and let cplex terminate by itself.
             self.abort()
             return
         self.handle_lazy_main_feasible_solution(main_mip, mindtpy_solver, config, opt)
@@ -744,9 +706,9 @@ def __call__(self):
                         mindtpy_solver.mip, None, mindtpy_solver, config, opt
                     )
                 except ValueError as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.error(
-                        str(e)
-                        + "\nUsually this error is caused by the MIP start solution causing a math domain error. "
+                        "Usually this error is caused by the MIP start solution causing a math domain error. "
                         "We will skip it."
                     )
                     return
@@ -782,6 +744,7 @@ def __call__(self):
                 )
             )
             mindtpy_solver.results.solver.termination_condition = tc.optimal
+            # TODO: check the performance difference if we don't use self.abort() and let cplex terminate by itself.
             self.abort()
             return
 
@@ -909,19 +872,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
 
@@ -952,15 +903,31 @@ 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.
+                begin_index, end_index = mindtpy_solver.integer_solution_to_cuts_index[
+                    mindtpy_solver.curr_int_sol
+                ]
+                for ind in range(begin_index, end_index + 1):
+                    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.integer_solution_to_cuts_index[
+                mindtpy_solver.curr_int_sol
+            ] = [cut_ind + 1, len(mindtpy_solver.mip.MindtPy_utils.cuts.oa_cuts)]
 
 
 def handle_lazy_main_feasible_solution_gurobi(cb_m, cb_opt, mindtpy_solver, config):

pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py

@@ -114,15 +114,15 @@ def evaluate_jacobian_equality_constraints(self):
             """Evaluate the Jacobian of the equality constraints."""
             return None
 
-        '''
+        """
         def _extract_and_assemble_fim(self):
             M = np.zeros((self.n_parameters, self.n_parameters))
             for i in range(self.n_parameters):
                 for k in range(self.n_parameters):                
                     M[i,k] = self._input_values[self.ele_to_order[(i,k)]]
 
             return M
-        '''
+        """
 
         def evaluate_jacobian_outputs(self):
             """Evaluate the Jacobian of the outputs."""