Skip to content

Commit

Permalink
Merge branch 'benchmark' into OA_convexification
Browse files Browse the repository at this point in the history
Zedong Peng committed Oct 24, 2023
2 parents 8d4a8d0 + 2a6f1d7 commit e14eb96
Showing 12 changed files with 507 additions and 43 deletions.
153 changes: 129 additions & 24 deletions pyomo/contrib/mindtpy/algorithm_base_class.py

Large diffs are not rendered by default.

8 changes: 8 additions & 0 deletions pyomo/contrib/mindtpy/config_options.py
Original file line number Diff line number Diff line change
@@ -494,6 +494,14 @@ def _add_common_configs(CONFIG):
domain=bool,
),
)
CONFIG.declare(
'load_solutions',
ConfigValue(
default=True,
description='Whether to load solutions in solve() function',
domain=bool,
),
)


def _add_subsolver_configs(CONFIG):
48 changes: 48 additions & 0 deletions pyomo/contrib/mindtpy/cut_generation.py
Original file line number Diff line number Diff line change
@@ -184,6 +184,54 @@ 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()
)
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,
7 changes: 6 additions & 1 deletion pyomo/contrib/mindtpy/feasibility_pump.py
Original file line number Diff line number Diff line change
@@ -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,
1 change: 1 addition & 0 deletions pyomo/contrib/mindtpy/global_outer_approximation.py
Original file line number Diff line number Diff line change
@@ -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)

13 changes: 11 additions & 2 deletions pyomo/contrib/mindtpy/outer_approximation.py
Original file line number Diff line number Diff line change
@@ -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.
8 changes: 4 additions & 4 deletions pyomo/contrib/mindtpy/single_tree.py
Original file line number Diff line number Diff line change
@@ -16,9 +16,8 @@
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.mindtpy.util import get_integer_solution, copy_var_list_values
from pyomo.contrib.gdpopt.util import (
copy_var_list_values,
get_main_elapsed_time,
time_code,
)
@@ -708,8 +707,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 +717,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
43 changes: 42 additions & 1 deletion pyomo/contrib/mindtpy/tests/MINLP_simple.py
Original file line number Diff line number Diff line change
@@ -35,14 +35,23 @@
RangeSet,
Var,
minimize,
Block,
)
from pyomo.common.collections import ComponentMap
from pyomo.contrib.mindtpy.tests.MINLP_simple_grey_box import GreyBoxModel
from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock


def build_model_external(m):
ex_model = GreyBoxModel(initial={"X1": 0, "X2": 0, "Y1": 0, "Y2": 1, "Y3": 1})
m.egb = ExternalGreyBoxBlock()
m.egb.set_external_model(ex_model)


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')
super(SimpleMINLP, self).__init__(*args, **kwargs)
@@ -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}
138 changes: 138 additions & 0 deletions pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
from pyomo.common.dependencies import numpy as np
import pyomo.common.dependencies.scipy.sparse as scipy_sparse
from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxModel
from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock


class GreyBoxModel(ExternalGreyBoxModel):
"""Greybox model to compute the example OF."""

def __init__(self, initial, use_exact_derivatives=True, verbose=True):
"""
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."""
# Not sure what this function should return with no 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:
pass
# 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))
Loading

0 comments on commit e14eb96

Please sign in to comment.