Merge pull request Pyomo#2988 from ZedongPeng/add_grey_box

Add the support of GreyBox in MindtPy

pyomo/contrib/mindtpy/algorithm_base_class.py

@@ -84,6 +84,9 @@
 
 single_tree, single_tree_available = attempt_import('pyomo.contrib.mindtpy.single_tree')
 tabu_list, tabu_list_available = attempt_import('pyomo.contrib.mindtpy.tabu_list')
+egb, egb_available = attempt_import(
+    'pyomo.contrib.pynumero.interfaces.external_grey_box'
+)
 
 
 class _MindtPyAlgorithm(object):
@@ -140,6 +143,8 @@ def __init__(self, **kwds):
         self.last_iter_cuts = False
         # Store the OA cuts generated in the mip_start_process.
         self.mip_start_lazy_oa_cuts = []
+        # Whether to load solutions in solve() function
+        self.load_solutions = True
 
     # Support use as a context manager under current solver API
     def __enter__(self):
@@ -289,7 +294,7 @@ def model_is_valid(self):
                 results = self.mip_opt.solve(
                     self.original_model,
                     tee=config.mip_solver_tee,
-                    load_solutions=False,
+                    load_solutions=self.load_solutions,
                     **config.mip_solver_args,
                 )
                 if len(results.solution) > 0:
@@ -323,6 +328,14 @@ def build_ordered_component_lists(self, model):
                 ctype=Constraint, active=True, descend_into=(Block)
             )
         )
+        if egb_available:
+            util_block.grey_box_list = list(
+                model.component_data_objects(
+                    ctype=egb.ExternalGreyBoxBlock, active=True, descend_into=(Block)
+                )
+            )
+        else:
+            util_block.grey_box_list = []
         util_block.linear_constraint_list = list(
             c
             for c in util_block.constraint_list
@@ -350,11 +363,20 @@ def build_ordered_component_lists(self, model):
 
         # We use component_data_objects rather than list(var_set) in order to
         # preserve a deterministic ordering.
-        util_block.variable_list = list(
-            v
-            for v in model.component_data_objects(ctype=Var, descend_into=(Block))
-            if v in var_set
-        )
+        if egb_available:
+            util_block.variable_list = list(
+                v
+                for v in model.component_data_objects(
+                    ctype=Var, descend_into=(Block, egb.ExternalGreyBoxBlock)
+                )
+                if v in var_set
+            )
+        else:
+            util_block.variable_list = list(
+                v
+                for v in model.component_data_objects(ctype=Var, descend_into=(Block))
+                if v in var_set
+            )
         util_block.discrete_variable_list = list(
             v for v in util_block.variable_list if v in var_set and v.is_integer()
         )
@@ -802,18 +824,21 @@ def init_rNLP(self, add_oa_cuts=True):
             MindtPy unable to handle the termination condition of the relaxed NLP.
         """
         config = self.config
-        m = self.working_model.clone()
+        self.rnlp = self.working_model.clone()
         config.logger.debug('Relaxed NLP: Solve relaxed integrality')
-        MindtPy = m.MindtPy_utils
-        TransformationFactory('core.relax_integer_vars').apply_to(m)
+        MindtPy = self.rnlp.MindtPy_utils
+        TransformationFactory('core.relax_integer_vars').apply_to(self.rnlp)
         nlp_args = dict(config.nlp_solver_args)
         update_solver_timelimit(self.nlp_opt, config.nlp_solver, self.timing, config)
         with SuppressInfeasibleWarning():
             results = self.nlp_opt.solve(
-                m, tee=config.nlp_solver_tee, load_solutions=False, **nlp_args
+                self.rnlp,
+                tee=config.nlp_solver_tee,
+                load_solutions=self.load_solutions,
+                **nlp_args,
             )
             if len(results.solution) > 0:
-                m.solutions.load_from(results)
+                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]
@@ -841,24 +866,24 @@ def init_rNLP(self, add_oa_cuts=True):
                 ):
                     # TODO: recover the opposite dual when cyipopt issue #2831 is solved.
                     dual_values = (
-                        list(-1 * m.dual[c] for c in MindtPy.constraint_list)
+                        list(-1 * self.rnlp.dual[c] for c in MindtPy.constraint_list)
                         if config.calculate_dual_at_solution
                         else None
                     )
                 else:
                     dual_values = (
-                        list(m.dual[c] for c in MindtPy.constraint_list)
+                        list(self.rnlp.dual[c] for c in MindtPy.constraint_list)
                         if config.calculate_dual_at_solution
                         else None
                     )
                 copy_var_list_values(
-                    m.MindtPy_utils.variable_list,
+                    self.rnlp.MindtPy_utils.variable_list,
                     self.mip.MindtPy_utils.variable_list,
                     config,
                 )
                 if config.init_strategy == 'FP':
                     copy_var_list_values(
-                        m.MindtPy_utils.variable_list,
+                        self.rnlp.MindtPy_utils.variable_list,
                         self.working_model.MindtPy_utils.variable_list,
                         config,
                     )
@@ -867,6 +892,7 @@ def init_rNLP(self, add_oa_cuts=True):
                     linearize_active=True,
                     linearize_violated=True,
                     cb_opt=None,
+                    nlp=self.rnlp,
                 )
                 for var in self.mip.MindtPy_utils.discrete_variable_list:
                     # We don't want to trigger the reset of the global stale
@@ -936,7 +962,7 @@ def init_max_binaries(self):
         mip_args = dict(config.mip_solver_args)
         update_solver_timelimit(self.mip_opt, config.mip_solver, self.timing, config)
         results = self.mip_opt.solve(
-            m, tee=config.mip_solver_tee, load_solutions=False, **mip_args
+            m, tee=config.mip_solver_tee, load_solutions=self.load_solutions, **mip_args
         )
         if len(results.solution) > 0:
             m.solutions.load_from(results)
@@ -1055,7 +1081,7 @@ def solve_subproblem(self):
                 results = self.nlp_opt.solve(
                     self.fixed_nlp,
                     tee=config.nlp_solver_tee,
-                    load_solutions=False,
+                    load_solutions=self.load_solutions,
                     **nlp_args,
                 )
                 if len(results.solution) > 0:
@@ -1153,6 +1179,7 @@ def handle_subproblem_optimal(self, fixed_nlp, cb_opt=None, fp=False):
             linearize_active=True,
             linearize_violated=True,
             cb_opt=cb_opt,
+            nlp=self.fixed_nlp,
         )
 
         var_values = list(v.value for v in fixed_nlp.MindtPy_utils.variable_list)
@@ -1229,6 +1256,7 @@ def handle_subproblem_infeasible(self, fixed_nlp, cb_opt=None):
             linearize_active=True,
             linearize_violated=True,
             cb_opt=cb_opt,
+            nlp=feas_subproblem,
         )
         # Add a no-good cut to exclude this discrete option
         var_values = list(v.value for v in fixed_nlp.MindtPy_utils.variable_list)
@@ -1311,6 +1339,12 @@ def solve_feasibility_subproblem(self):
         update_solver_timelimit(
             self.feasibility_nlp_opt, config.nlp_solver, self.timing, config
         )
+        TransformationFactory('contrib.deactivate_trivial_constraints').apply_to(
+            feas_subproblem,
+            tmp=True,
+            ignore_infeasible=False,
+            tolerance=config.constraint_tolerance,
+        )
         with SuppressInfeasibleWarning():
             try:
                 with time_code(self.timing, 'feasibility subproblem'):
@@ -1346,6 +1380,9 @@ def solve_feasibility_subproblem(self):
             constr.activate()
         active_obj.activate()
         MindtPy.feas_obj.deactivate()
+        TransformationFactory('contrib.deactivate_trivial_constraints').revert(
+            feas_subproblem
+        )
         return feas_subproblem, feas_soln
 
     def handle_feasibility_subproblem_tc(self, subprob_terminate_cond, MindtPy):
@@ -1480,7 +1517,10 @@ def fix_dual_bound(self, last_iter_cuts):
                 self.mip_opt, config.mip_solver, self.timing, config
             )
             main_mip_results = self.mip_opt.solve(
-                self.mip, tee=config.mip_solver_tee, load_solutions=False, **mip_args
+                self.mip,
+                tee=config.mip_solver_tee,
+                load_solutions=self.load_solutions,
+                **mip_args,
             )
             if len(main_mip_results.solution) > 0:
                 self.mip.solutions.load_from(main_mip_results)
@@ -1564,7 +1604,10 @@ def solve_main(self):
 
         try:
             main_mip_results = self.mip_opt.solve(
-                self.mip, tee=config.mip_solver_tee, load_solutions=False, **mip_args
+                self.mip,
+                tee=config.mip_solver_tee,
+                load_solutions=self.load_solutions,
+                **mip_args,
             )
             # update_attributes should be before load_from(main_mip_results), since load_from(main_mip_results) may fail.
             if len(main_mip_results.solution) > 0:
@@ -1617,7 +1660,10 @@ def solve_fp_main(self):
         mip_args = self.set_up_mip_solver()
 
         main_mip_results = self.mip_opt.solve(
-            self.mip, tee=config.mip_solver_tee, load_solutions=False, **mip_args
+            self.mip,
+            tee=config.mip_solver_tee,
+            load_solutions=self.load_solutions,
+            **mip_args,
         )
         # update_attributes should be before load_from(main_mip_results), since load_from(main_mip_results) may fail.
         # if config.single_tree or config.use_tabu_list:
@@ -1659,7 +1705,7 @@ def solve_regularization_main(self):
         main_mip_results = self.regularization_mip_opt.solve(
             self.mip,
             tee=config.mip_solver_tee,
-            load_solutions=False,
+            load_solutions=self.load_solutions,
             **dict(config.mip_solver_args),
         )
         if len(main_mip_results.solution) > 0:
@@ -1871,7 +1917,7 @@ def handle_main_unbounded(self, main_mip):
             main_mip_results = self.mip_opt.solve(
                 main_mip,
                 tee=config.mip_solver_tee,
-                load_solutions=False,
+                load_solutions=self.load_solutions,
                 **config.mip_solver_args,
             )
             if len(main_mip_results.solution) > 0:
@@ -2200,6 +2246,17 @@ def check_config(self):
                     config.logger.info("Solution pool does not support APPSI solver.")
                 config.mip_solver = 'cplex_persistent'
 
+        # related to https://github.com/Pyomo/pyomo/issues/2363
+        if (
+            'appsi' in config.mip_solver
+            or 'appsi' in config.nlp_solver
+            or (
+                config.mip_regularization_solver is not None
+                and 'appsi' in config.mip_regularization_solver
+            )
+        ):
+            self.load_solutions = False
+
     ################################################################################################################################
     # Feasibility Pump
 
@@ -2263,7 +2320,10 @@ def solve_fp_subproblem(self):
         with SuppressInfeasibleWarning():
             with time_code(self.timing, 'fp subproblem'):
                 results = self.nlp_opt.solve(
-                    fp_nlp, tee=config.nlp_solver_tee, load_solutions=False, **nlp_args
+                    fp_nlp,
+                    tee=config.nlp_solver_tee,
+                    load_solutions=self.load_solutions,
+                    **nlp_args,
                 )
                 if len(results.solution) > 0:
                     fp_nlp.solutions.load_from(results)
@@ -2482,6 +2542,9 @@ def initialize_mip_problem(self):
             getattr(self.mip, 'ipopt_zU_out', _DoNothing()).deactivate()
 
         MindtPy = self.mip.MindtPy_utils
+        if len(MindtPy.grey_box_list) > 0:
+            for grey_box in MindtPy.grey_box_list:
+                grey_box.deactivate()
 
         if config.init_strategy == 'FP':
             MindtPy.cuts.fp_orthogonality_cuts = ConstraintList(

pyomo/contrib/mindtpy/cut_generation.py

@@ -181,6 +181,55 @@ def add_oa_cuts(
                         )
 
 
+def add_oa_cuts_for_grey_box(
+    target_model, jacobians_model, config, objective_sense, mip_iter, cb_opt=None
+):
+    sign_adjust = -1 if objective_sense == minimize else 1
+    if config.add_slack:
+        slack_var = target_model.MindtPy_utils.cuts.slack_vars.add()
+    for target_model_grey_box, jacobian_model_grey_box in zip(
+        target_model.MindtPy_utils.grey_box_list,
+        jacobians_model.MindtPy_utils.grey_box_list,
+    ):
+        jacobian_matrix = (
+            jacobian_model_grey_box.get_external_model()
+            .evaluate_jacobian_outputs()
+            .toarray()
+        )
+        # Enumerate over values works well now. However, it might be stable if the values() method changes.
+        for index, output in enumerate(target_model_grey_box.outputs.values()):
+            dual_value = jacobians_model.dual[jacobian_model_grey_box][
+                output.name.replace("outputs", "output_constraints")
+            ]
+            target_model.MindtPy_utils.cuts.oa_cuts.add(
+                expr=copysign(1, sign_adjust * dual_value)
+                * (
+                    sum(
+                        jacobian_matrix[index][var_index] * (var - value(var))
+                        for var_index, var in enumerate(
+                            target_model_grey_box.inputs.values()
+                        )
+                    )
+                )
+                - (output - value(output))
+                - (slack_var if config.add_slack else 0)
+                <= 0
+            )
+    # TODO: gurobi_persistent currently does not support greybox model.
+    # https://github.com/Pyomo/pyomo/issues/3000
+    # if (
+    #     config.single_tree
+    #     and config.mip_solver == 'gurobi_persistent'
+    #     and mip_iter > 0
+    #     and cb_opt is not None
+    # ):
+    #     cb_opt.cbLazy(
+    #         target_model.MindtPy_utils.cuts.oa_cuts[
+    #             len(target_model.MindtPy_utils.cuts.oa_cuts)
+    #         ]
+    #     )
+
+
 def add_ecp_cuts(
     target_model,
     jacobians,

pyomo/contrib/mindtpy/feasibility_pump.py

@@ -52,7 +52,12 @@ def initialize_mip_problem(self):
         )
 
     def add_cuts(
-        self, dual_values, linearize_active=True, linearize_violated=True, cb_opt=None
+        self,
+        dual_values,
+        linearize_active=True,
+        linearize_violated=True,
+        cb_opt=None,
+        nlp=None,
     ):
         add_oa_cuts(
             self.mip,

pyomo/contrib/mindtpy/global_outer_approximation.py

@@ -95,6 +95,7 @@ def add_cuts(
         linearize_active=True,
         linearize_violated=True,
         cb_opt=None,
+        nlp=None,
     ):
         add_affine_cuts(self.mip, self.config, self.timing)
 

pyomo/contrib/mindtpy/outer_approximation.py

@@ -16,7 +16,7 @@
 from pyomo.opt import SolverFactory
 from pyomo.contrib.mindtpy.config_options import _get_MindtPy_OA_config
 from pyomo.contrib.mindtpy.algorithm_base_class import _MindtPyAlgorithm
-from pyomo.contrib.mindtpy.cut_generation import add_oa_cuts
+from pyomo.contrib.mindtpy.cut_generation import add_oa_cuts, add_oa_cuts_for_grey_box
 
 
 @SolverFactory.register(
@@ -102,7 +102,12 @@ def initialize_mip_problem(self):
         )
 
     def add_cuts(
-        self, dual_values, linearize_active=True, linearize_violated=True, cb_opt=None
+        self,
+        dual_values,
+        linearize_active=True,
+        linearize_violated=True,
+        cb_opt=None,
+        nlp=None,
     ):
         add_oa_cuts(
             self.mip,
@@ -117,6 +122,10 @@ def add_cuts(
             linearize_active,
             linearize_violated,
         )
+        if len(self.mip.MindtPy_utils.grey_box_list) > 0:
+            add_oa_cuts_for_grey_box(
+                self.mip, nlp, self.config, self.objective_sense, self.mip_iter, cb_opt
+            )
 
     def deactivate_no_good_cuts_when_fixing_bound(self, no_good_cuts):
         # Only deactivate the last OA cuts may not be correct.