Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add the support of GreyBox in MindtPy #2988

Merged
merged 38 commits into from
Nov 27, 2023
Merged

Add the support of GreyBox in MindtPy #2988

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
38 commits
Select commit Hold shift + click to select a range
bb03a18
fix load_state_into_pyomo bug
ZedongPeng Aug 21, 2023
1d508e7
Merge branch 'main' into add_grey_box
ZedongPeng Aug 25, 2023
ca54afc
add the support of greybox in mindtpy
ZedongPeng Sep 8, 2023
ec05a3d
add grey box test in mindtpy
ZedongPeng Sep 8, 2023
5048b42
Merge branch 'main' into add_grey_box
ZedongPeng Sep 8, 2023
c2877a3
black format
ZedongPeng Sep 8, 2023
ccff416
black format
ZedongPeng Sep 8, 2023
e461e93
add comment
ZedongPeng Sep 11, 2023
13a7213
remove the support of greybox in LP/NLP gurobi
ZedongPeng Sep 19, 2023
20ccbdc
black format
ZedongPeng Sep 19, 2023
4288e33
update import source
ZedongPeng Sep 20, 2023
8d7f6c5
remove redundant import
ZedongPeng Sep 20, 2023
807e4f5
add config check for load_solutions
ZedongPeng Sep 20, 2023
30771b3
recover numpy and scipy import
ZedongPeng Sep 20, 2023
231abc6
Merge branch 'main' into add_grey_box
ZedongPeng Sep 20, 2023
72b142b
change numpy and scipy import
ZedongPeng Sep 20, 2023
ff401df
change scipy import
ZedongPeng Sep 21, 2023
d82dcde
fix import error
ZedongPeng Sep 21, 2023
cb1c2a9
black format
ZedongPeng Sep 21, 2023
e80c6dc
update mindtpy import in pynumero
ZedongPeng Sep 21, 2023
8959666
remove redundant scipy import
ZedongPeng Sep 21, 2023
e0c245b
add skip
ZedongPeng Sep 25, 2023
5ffaeb1
black format
ZedongPeng Sep 25, 2023
260b2d6
move numpy and scipy check forward
ZedongPeng Sep 25, 2023
c10e212
Merge branch 'main' into add_grey_box
ZedongPeng Sep 25, 2023
46cf8a5
fix
ZedongPeng Oct 31, 2023
ba135b6
black format
ZedongPeng Oct 31, 2023
8eacf0b
fix import bug
ZedongPeng Oct 31, 2023
bd12664
fix import bug
ZedongPeng Nov 1, 2023
ebe91a6
black format
ZedongPeng Nov 1, 2023
2b45756
fix import bug
ZedongPeng Nov 1, 2023
96efbd4
Merge branch 'main' into add_grey_box
emma58 Nov 14, 2023
73f2e16
change load_solutions to attribute
ZedongPeng Nov 14, 2023
a08fd36
black format
ZedongPeng Nov 14, 2023
ae84558
remove load_solutions configuration
ZedongPeng Nov 15, 2023
7bf0268
add comments to enumerate over values()
ZedongPeng Nov 15, 2023
993290f
black format
ZedongPeng Nov 15, 2023
3cf0fc2
Merge branch 'main' into add_grey_box
ZedongPeng Nov 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
109 changes: 86 additions & 23 deletions pyomo/contrib/mindtpy/algorithm_base_class.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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()
)
Expand Down Expand Up @@ -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()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you really want to store this publicly on the solver object? If your users can find it and it doesn't start with underscore, they might start to rely on it. That's not necessarily bad--just something to think through.

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]
Expand Down Expand Up @@ -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,
)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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,
)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use of this transformation is new. Is it motivated by our new tests? Is the transformation well-supported? @emma58

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After MindtPy rewrite, we repeatedly use the fixed-nlp subproblem. Therefore, we will call 'contrib.deactivate_trivial_constraints' transformation and revert it after solving. I missed the feasibility subproblem and this might be a bug.

with SuppressInfeasibleWarning():
try:
with time_code(self.timing, 'feasibility subproblem'):
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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(
Expand Down
49 changes: 49 additions & 0 deletions pyomo/contrib/mindtpy/cut_generation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
Comment on lines +200 to +217
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure that this is always going to create a correct cut, even if greybox changes its behavior? The part that looks potentially scary to me is calling enumerate over getting the values() from dictionaries coming from greybox. You will get deterministic ordering from values() as long as the dictionary is constructed in the same order every time, but if that changes, is this still right? Or does greybox make a promise about consistency with that order and the jacobian? (This is me wondering out of ignorance, but it may be worth adding a comment even if it is fine.)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we don't call enumerate over values, any suggestions about how to implement it? Directly enumerate over variable?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't really know--you would know better than me. It's completely possible that this is fine. My concern is that you will end up with the wrong coefficient (jacobian_matrix[index][var_index]) paired with the wrong var - value(var) term if the target_model_greyboox.inputs dictionary is every constructed in a different order that doesn't correspond to the indices in jacobian_matrix. I don't know greybox well at all, so maybe it makes a promise that what I just said will never happen, in which case your code is totally fine. I'm just asking if you know if that's the case or not. It's always worth thinking carefully when you are relying on the order things come up when iterating over a dictionary...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Currently, I cannot come up with a good way to improve this since the evaluate_jacobian_outputs returns a scipy.sparse._coo.coo_matrix.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe just put in a comment that you are relying on the order in which items are added to the target_model_greybox.inputs dictionary, so that future you (or someone else) will have a hint if this ever doesn't work.

# TODO: gurobi_persistent currently does not support greybox model.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this become an issue?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gurobi will fail even when we deactivate the greybox block. Therefore, we cannot use Gurobi as the mip solver.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this a bug in gurobi_persistent? Why does it behave differently?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will submit an issue related to this.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@emma58 Could you provide me a reference example?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

from pyomo.environ import *

m = ConcreteModel()
m.b = Block()
m.b.x = Var(bounds=(2, 5))
m.y = Var(domain=NonNegativeReals)
m.c = Constraint(expr=m.b.x - m.y <= 4)
m.obj = Objective(expr=m.b.x + m.y)

m.b.deactivate()

m.somewhere_safe_to_keep_things = Block()
m.somewhere_safe_to_keep_things.x = Reference(m.b.x)

opt = SolverFactory('gurobi_persistent')
opt.set_instance(m)

opt.solve(m, tee=True)

This is the main idea. You'll want to figure out how to do your own accounting in terms of finding the Vars you need to reference and naming the references, etc.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

#3000 Issue opened

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@emma58 . I don't think Reference is an elegant way to resolve this issue. I will comment the code of adding Gurobi lazy constraints and only support greybox in multi-tree implementation this time.
Btw, since persistent solver will be totally replaced by appsi_solver, does appsi_solver support both CPLEX and Gurobi callbacks as persistent solver now? @michaelbynum

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

appsi_gurobi does support callbacks. appsi_cplex will, but does not currently.

# 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)
# ]
# )


ZedongPeng marked this conversation as resolved.
Show resolved Hide resolved
def add_ecp_cuts(
target_model,
jacobians,
Expand Down
7 changes: 6 additions & 1 deletion pyomo/contrib/mindtpy/feasibility_pump.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
1 change: 1 addition & 0 deletions pyomo/contrib/mindtpy/global_outer_approximation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
13 changes: 11 additions & 2 deletions pyomo/contrib/mindtpy/outer_approximation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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,
Expand All @@ -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.
Expand Down
Loading