diff --git a/pyomo/contrib/mindtpy/algorithm_base_class.py b/pyomo/contrib/mindtpy/algorithm_base_class.py index 7def1dcaab3..c9169ab8c62 100644 --- a/pyomo/contrib/mindtpy/algorithm_base_class.py +++ b/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( diff --git a/pyomo/contrib/mindtpy/cut_generation.py b/pyomo/contrib/mindtpy/cut_generation.py index c0449054baa..28d302104a3 100644 --- a/pyomo/contrib/mindtpy/cut_generation.py +++ b/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, diff --git a/pyomo/contrib/mindtpy/feasibility_pump.py b/pyomo/contrib/mindtpy/feasibility_pump.py index 5716400598a..990f56b8f93 100644 --- a/pyomo/contrib/mindtpy/feasibility_pump.py +++ b/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, diff --git a/pyomo/contrib/mindtpy/global_outer_approximation.py b/pyomo/contrib/mindtpy/global_outer_approximation.py index ee3ffb62f55..dfb7ef54630 100644 --- a/pyomo/contrib/mindtpy/global_outer_approximation.py +++ b/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) diff --git a/pyomo/contrib/mindtpy/outer_approximation.py b/pyomo/contrib/mindtpy/outer_approximation.py index 99d9cea1bd4..6cf0b26cb37 100644 --- a/pyomo/contrib/mindtpy/outer_approximation.py +++ b/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. diff --git a/pyomo/contrib/mindtpy/tests/MINLP_simple.py b/pyomo/contrib/mindtpy/tests/MINLP_simple.py index 5663c93af8b..7454b595986 100644 --- a/pyomo/contrib/mindtpy/tests/MINLP_simple.py +++ b/pyomo/contrib/mindtpy/tests/MINLP_simple.py @@ -35,16 +35,25 @@ RangeSet, Var, minimize, + Block, ) from pyomo.common.collections import ComponentMap +from pyomo.contrib.mindtpy.tests.MINLP_simple_grey_box import ( + GreyBoxModel, + build_model_external, +) class SimpleMINLP(ConcreteModel): """Convex MINLP problem Assignment 6 APSE.""" - def __init__(self, *args, **kwargs): + def __init__(self, grey_box=False, *args, **kwargs): """Create the problem.""" kwargs.setdefault('name', 'SimpleMINLP') + if grey_box and GreyBoxModel is None: + m = None + return + super(SimpleMINLP, self).__init__(*args, **kwargs) m = self @@ -83,6 +92,38 @@ def __init__(self, *args, **kwargs): m.objective = Objective( expr=Y[1] + 1.5 * Y[2] + 0.5 * Y[3] + X[1] ** 2 + X[2] ** 2, sense=minimize ) + + if not grey_box: + m.objective = Objective( + expr=Y[1] + 1.5 * Y[2] + 0.5 * Y[3] + X[1] ** 2 + X[2] ** 2, + sense=minimize, + ) + else: + + def _model_i(b): + build_model_external(b) + + m.my_block = Block(rule=_model_i) + + for i in m.I: + + def eq_inputX(m): + return m.X[i] == m.my_block.egb.inputs["X" + str(i)] + + con_name = "con_X_" + str(i) + m.add_component(con_name, Constraint(expr=eq_inputX)) + + for j in m.J: + + def eq_inputY(m): + return m.Y[j] == m.my_block.egb.inputs["Y" + str(j)] + + con_name = "con_Y_" + str(j) + m.add_component(con_name, Constraint(expr=eq_inputY)) + + # add objective + m.objective = Objective(expr=m.my_block.egb.outputs['z'], sense=minimize) + """Bound definitions""" # x (continuous) upper bounds x_ubs = {1: 4, 2: 4} diff --git a/pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py b/pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py new file mode 100644 index 00000000000..547efc0a74c --- /dev/null +++ b/pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py @@ -0,0 +1,151 @@ +from pyomo.common.dependencies import numpy as np +import pyomo.common.dependencies.scipy.sparse as scipy_sparse +from pyomo.common.dependencies import attempt_import + +egb, egb_available = attempt_import( + 'pyomo.contrib.pynumero.interfaces.external_grey_box' +) + +if egb_available: + + class GreyBoxModel(egb.ExternalGreyBoxModel): + """Greybox model to compute the example objective function.""" + + def __init__(self, initial, use_exact_derivatives=True, verbose=False): + """ + Parameters + + use_exact_derivatives: bool + If True, the exact derivatives are used. If False, the finite difference + approximation is used. + verbose: bool + If True, print information about the model. + """ + self._use_exact_derivatives = use_exact_derivatives + self.verbose = verbose + self.initial = initial + + # For use with exact Hessian + self._output_con_mult_values = np.zeros(1) + + if not use_exact_derivatives: + raise NotImplementedError( + "use_exact_derivatives == False not supported" + ) + + def input_names(self): + """Return the names of the inputs.""" + self.input_name_list = ["X1", "X2", "Y1", "Y2", "Y3"] + + return self.input_name_list + + def equality_constraint_names(self): + """Return the names of the equality constraints.""" + # no equality constraints + return [] + + def output_names(self): + """Return the names of the outputs.""" + return ['z'] + + def set_output_constraint_multipliers(self, output_con_multiplier_values): + """Set the values of the output constraint multipliers.""" + # because we only have one output constraint + assert len(output_con_multiplier_values) == 1 + np.copyto(self._output_con_mult_values, output_con_multiplier_values) + + def finalize_block_construction(self, pyomo_block): + """Finalize the construction of the ExternalGreyBoxBlock.""" + if self.initial is not None: + print("initialized") + pyomo_block.inputs["X1"].value = self.initial["X1"] + pyomo_block.inputs["X2"].value = self.initial["X2"] + pyomo_block.inputs["Y1"].value = self.initial["Y1"] + pyomo_block.inputs["Y2"].value = self.initial["Y2"] + pyomo_block.inputs["Y3"].value = self.initial["Y3"] + + else: + print("uninitialized") + for n in self.input_name_list: + pyomo_block.inputs[n].value = 1 + + pyomo_block.inputs["X1"].setub(4) + pyomo_block.inputs["X1"].setlb(0) + + pyomo_block.inputs["X2"].setub(4) + pyomo_block.inputs["X2"].setlb(0) + + pyomo_block.inputs["Y1"].setub(1) + pyomo_block.inputs["Y1"].setlb(0) + + pyomo_block.inputs["Y2"].setub(1) + pyomo_block.inputs["Y2"].setlb(0) + + pyomo_block.inputs["Y3"].setub(1) + pyomo_block.inputs["Y3"].setlb(0) + + def set_input_values(self, input_values): + """Set the values of the inputs.""" + self._input_values = list(input_values) + + def evaluate_equality_constraints(self): + """Evaluate the equality constraints.""" + return None + + def evaluate_outputs(self): + """Evaluate the output of the model.""" + # form matrix as a list of lists + # M = self._extract_and_assemble_fim() + x1 = self._input_values[0] + x2 = self._input_values[1] + y1 = self._input_values[2] + y2 = self._input_values[3] + y3 = self._input_values[4] + # z + z = x1**2 + x2**2 + y1 + 1.5 * y2 + 0.5 * y3 + + if self.verbose: + print("\n Consider inputs [x1,x2,y1,y2,y3] =\n", x1, x2, y1, y2, y3) + print(" z = ", z, "\n") + + return np.asarray([z], dtype=np.float64) + + 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.""" + if self._use_exact_derivatives: + # compute gradient of log determinant + row = np.zeros(5) # to store row index + col = np.zeros(5) # to store column index + data = np.zeros(5) # to store data + + row[0], col[0], data[0] = (0, 0, 2 * self._input_values[0]) # x1 + row[0], col[1], data[1] = (0, 1, 2 * self._input_values[1]) # x2 + row[0], col[2], data[2] = (0, 2, 1) # y1 + row[0], col[3], data[3] = (0, 3, 1.5) # y2 + row[0], col[4], data[4] = (0, 4, 0.5) # y3 + + # sparse matrix + return scipy_sparse.coo_matrix((data, (row, col)), shape=(1, 5)) + + def build_model_external(m): + ex_model = GreyBoxModel(initial={"X1": 0, "X2": 0, "Y1": 0, "Y2": 1, "Y3": 1}) + m.egb = egb.ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + +else: + GreyBoxModel = None + build_model_external = None diff --git a/pyomo/contrib/mindtpy/tests/test_mindtpy_grey_box.py b/pyomo/contrib/mindtpy/tests/test_mindtpy_grey_box.py new file mode 100644 index 00000000000..f84136ca6bf --- /dev/null +++ b/pyomo/contrib/mindtpy/tests/test_mindtpy_grey_box.py @@ -0,0 +1,76 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +"""Tests for the MindtPy solver.""" +from pyomo.core.expr.calculus.diff_with_sympy import differentiate_available +import pyomo.common.unittest as unittest +from pyomo.environ import SolverFactory, value, maximize +from pyomo.opt import TerminationCondition +from pyomo.common.dependencies import numpy_available, scipy_available +from pyomo.contrib.mindtpy.tests.MINLP_simple import SimpleMINLP as SimpleMINLP + +model_list = [SimpleMINLP(grey_box=True)] +required_solvers = ('cyipopt', 'glpk') +if all(SolverFactory(s).available(exception_flag=False) for s in required_solvers): + subsolvers_available = True +else: + subsolvers_available = False + + +@unittest.skipIf(model_list[0] is None, 'Unable to generate the Grey Box model.') +@unittest.skipIf( + not subsolvers_available, + 'Required subsolvers %s are not available' % (required_solvers,), +) +@unittest.skipIf( + not differentiate_available, 'Symbolic differentiation is not available' +) +class TestMindtPy(unittest.TestCase): + """Tests for the MindtPy solver plugin.""" + + def check_optimal_solution(self, model, places=1): + for var in model.optimal_solution: + self.assertAlmostEqual( + var.value, model.optimal_solution[var], places=places + ) + + def test_OA_rNLP(self): + """Test the outer approximation decomposition algorithm.""" + with SolverFactory('mindtpy') as opt: + for model in model_list: + model = model.clone() + results = opt.solve( + model, + strategy='OA', + init_strategy='rNLP', + mip_solver=required_solvers[1], + nlp_solver=required_solvers[0], + calculate_dual_at_solution=True, + nlp_solver_args={ + 'options': { + 'hessian_approximation': 'limited-memory', + 'linear_solver': 'mumps', + } + }, + ) + + self.assertIn( + results.solver.termination_condition, + [TerminationCondition.optimal, TerminationCondition.feasible], + ) + self.assertAlmostEqual( + value(model.objective.expr), model.optimal_value, places=1 + ) + self.check_optimal_solution(model) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/pynumero/interfaces/external_grey_box.py b/pyomo/contrib/pynumero/interfaces/external_grey_box.py index 8fd728a7c9b..642fd3bf310 100644 --- a/pyomo/contrib/pynumero/interfaces/external_grey_box.py +++ b/pyomo/contrib/pynumero/interfaces/external_grey_box.py @@ -11,8 +11,8 @@ import abc import logging -import numpy as np from scipy.sparse import coo_matrix +from pyomo.common.dependencies import numpy as np from pyomo.common.deprecation import RenamedClass from pyomo.common.log import is_debug_set diff --git a/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py b/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py index ea51e38983d..945e9a05f51 100644 --- a/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py +++ b/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py @@ -426,7 +426,13 @@ def load_state_into_pyomo(self, bound_multipliers=None): # not we are minimizing or maximizing - this is done in the ASL interface # for ipopt, but does not appear to be done in cyipopt. obj_sign = 1.0 - objs = list(m.component_data_objects(ctype=pyo.Objective, descend_into=True)) + # since we will assert the number of objective functions, + # we only focus on active objective function. + objs = list( + m.component_data_objects( + ctype=pyo.Objective, active=True, descend_into=True + ) + ) assert len(objs) == 1 if objs[0].sense == pyo.maximize: obj_sign = -1.0