Source code for pyomo.gdp.plugins.cuttingplane
# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2024
# 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.
# ___________________________________________________________________________
"""
Cutting plane-based GDP reformulation.
Implements a general cutting plane-based reformulation for linear and
convex GDPs.
"""
from pyomo.common.config import (
ConfigBlock,
ConfigValue,
NonNegativeFloat,
PositiveInt,
In,
)
from pyomo.common.modeling import unique_component_name
from pyomo.core import (
Any,
Block,
Constraint,
Objective,
Param,
Var,
SortComponents,
Transformation,
TransformationFactory,
value,
NonNegativeIntegers,
Reals,
NonNegativeReals,
Suffix,
ComponentMap,
)
from pyomo.core.expr import differentiate
from pyomo.common.collections import ComponentSet
from pyomo.opt import SolverFactory
from pyomo.repn import generate_standard_repn
from pyomo.gdp import Disjunct, Disjunction, GDP_Error
from pyomo.gdp.util import (
verify_successful_solve,
NORMAL,
clone_without_expression_components,
)
from pyomo.contrib.fme.fourier_motzkin_elimination import (
Fourier_Motzkin_Elimination_Transformation,
)
import logging
logger = logging.getLogger('pyomo.gdp.cuttingplane')
def _get_constraint_exprs(constraints, hull_to_bigm_map):
"""Returns a list of expressions which are constrain.expr translated
into the bigm space, for each constraint in constraints.
"""
cuts = []
for cons in constraints:
cuts.append(
clone_without_expression_components(cons.expr, substitute=hull_to_bigm_map)
)
return cuts
def _constraint_tight(model, constraint, TOL):
"""
Returns a list [a,b] where a is -1 if the lower bound is tight or
slightly violated, b is 1 if the upper bound is tight of slightly
violated, and [a,b]=[-1,1] if we have an exactly satisfied (or
slightly violated) equality.
"""
val = value(constraint.body)
ans = [0, 0]
if constraint.lower is not None:
if val - value(constraint.lower) <= TOL:
# tight or in violation of LB
ans[0] -= 1
if constraint.upper is not None:
if value(constraint.upper) - val <= TOL:
# tight or in violation of UB
ans[1] += 1
return ans
def _get_linear_approximation_expr(normal_vec, point):
"""Returns constraint linearly approximating constraint normal to normal_vec
at point"""
body = 0
for coef, v in zip(point, normal_vec):
body -= coef * v
return body >= -sum(normal_vec[idx] * v.value for (idx, v) in enumerate(point))
def _precompute_potentially_useful_constraints(transBlock_rHull, disaggregated_vars):
instance_rHull = transBlock_rHull.model()
constraints = transBlock_rHull.constraints_for_FME = []
for constraint in instance_rHull.component_data_objects(
Constraint, active=True, descend_into=Block, sort=SortComponents.deterministic
):
# we don't care about anything that does not involve at least one
# disaggregated variable.
repn = generate_standard_repn(constraint.body)
for v in repn.linear_vars + repn.quadratic_vars + repn.nonlinear_vars:
# ESJ: This is why disaggregated_vars is a ComponentSet
if v in disaggregated_vars:
constraints.append(constraint)
break
[docs]
def create_cuts_fme(
transBlock_rHull,
var_info,
hull_to_bigm_map,
rBigM_linear_constraints,
rHull_vars,
disaggregated_vars,
norm,
cut_threshold,
zero_tolerance,
integer_arithmetic,
constraint_tolerance,
):
"""Returns a cut which removes x* from the relaxed bigm feasible region.
Finds all the constraints which are tight at xhat (assumed to be the
solution currently in instance_rHull), and calculates a composite normal
vector by summing the vectors normal to each of these constraints. Then
Fourier-Motzkin elimination is used to project the disaggregated variables
out of the polyhedron formed by the composite normal and the collection
of tight constraints. This results in multiple cuts, of which we select
one that cuts of x* by the greatest margin, as long as that margin is
more than cut_threshold. If no cut satisfies the margin specified by
cut_threshold, we return None.
Parameters
-----------
transBlock_rHull: transformation block on relaxed hull instance
var_info: List of tuples (rBigM_var, rHull_var, xstar_param)
hull_to_bigm_map: For expression substitution, maps id(hull_var) to
corresponding bigm var
rBigM_linear_constraints: list of linear constraints in relaxed bigM
rHull_vars: list of all variables in relaxed hull
disaggregated_vars: ComponentSet of disaggregated variables in hull
reformulation
norm: norm used in the separation problem
cut_threshold: Amount x* needs to be infeasible in generated cut in order
to consider the cut for addition to the bigM model.
zero_tolerance: Tolerance at which a float will be treated as 0 during
Fourier-Motzkin elimination
integer_arithmetic: boolean, whether or not to require Fourier-Motzkin
Elimination does integer arithmetic. Only possible
when all data is integer.
constraint_tolerance: Tolerance at which we will consider a constraint
tight.
"""
instance_rHull = transBlock_rHull.model()
# In the first iteration, we will compute a list of constraints that could
# ever be interesting: Everything that involves at least one disaggregated
# variable.
if transBlock_rHull.component("constraints_for_FME") is None:
_precompute_potentially_useful_constraints(transBlock_rHull, disaggregated_vars)
tight_constraints = Block()
conslist = tight_constraints.constraints = Constraint(NonNegativeIntegers)
conslist.construct()
something_interesting = False
for constraint in transBlock_rHull.constraints_for_FME:
multipliers = _constraint_tight(
instance_rHull, constraint, constraint_tolerance
)
for multiplier in multipliers:
if multiplier:
something_interesting = True
f = constraint.body
firstDerivs = differentiate(f, wrt_list=rHull_vars)
normal_vec = [multiplier * value(_) for _ in firstDerivs]
# check if constraint is linear
if f.polynomial_degree() == 1:
conslist[len(conslist)] = constraint.expr
else:
# we will use the linear approximation of this constraint at
# x_hat
conslist[len(conslist)] = _get_linear_approximation_expr(
normal_vec, rHull_vars
)
# NOTE: we now have all the tight Constraints (in the pyomo sense of the
# word "Constraint"), but we are missing some variable bounds. The ones for
# the disaggregated variables will be added by FME
# It is possible that the separation problem returned a point in the
# interior of the convex hull. It is also possible that the only active
# constraints do not involve the disaggregated variables. In these
# situations, there are not constraints from which to create a valid cut.
if not something_interesting:
return None
tight_constraints.construct()
logger.info(
"Calling FME transformation on %s constraints to eliminate"
" %s variables" % (len(tight_constraints.constraints), len(disaggregated_vars))
)
TransformationFactory('contrib.fourier_motzkin_elimination').apply_to(
tight_constraints,
vars_to_eliminate=disaggregated_vars,
zero_tolerance=zero_tolerance,
do_integer_arithmetic=integer_arithmetic,
projected_constraints_name="fme_constraints",
)
fme_results = tight_constraints.fme_constraints
projected_constraints = [cons for i, cons in fme_results.items()]
# we created these constraints with the variables from rHull. We
# actually need constraints for BigM and rBigM now!
cuts = _get_constraint_exprs(projected_constraints, hull_to_bigm_map)
# We likely have some cuts that duplicate other constraints now. We will
# filter them to make sure that they do in fact cut off x*. If that's the
# case, we know they are not already in the BigM relaxation. Because they
# came from FME, they are very likely redundant, so we'll keep the best one
# we find
best = 0
best_cut = None
cuts_to_keep = []
for i, cut in enumerate(cuts):
# x* is still in rBigM, so we can just remove this constraint if it
# is satisfied at x*
logger.info("FME: Post-processing cut %s" % cut)
if value(cut):
logger.info("FME:\t Doesn't cut off x*")
continue
# we have found a constraint which cuts of x* by some convincing amount
# and is not already in rBigM.
cuts_to_keep.append(i)
# We know cut is lb <= expr and that it's violated
assert len(cut.args) == 2
cut_off = value(cut.args[0]) - value(cut.args[1])
if cut_off > cut_threshold and cut_off > best:
best = cut_off
best_cut = cut
logger.info("FME:\t New best cut: Cuts off x* by %s." % best)
# NOTE: this is not used right now, but it's not hard to imagine a world in
# which we would want to keep multiple cuts from FME, so leaving it in for
# now.
cuts = [cuts[i] for i in cuts_to_keep]
if best_cut is not None:
return [best_cut]
return None
[docs]
def create_cuts_normal_vector(
transBlock_rHull,
var_info,
hull_to_bigm_map,
rBigM_linear_constraints,
rHull_vars,
disaggregated_vars,
norm,
cut_threshold,
zero_tolerance,
integer_arithmetic,
constraint_tolerance,
):
"""Returns a cut which removes x* from the relaxed bigm feasible region.
Ignores all parameters except var_info and cut_threshold, and constructs
a cut at x_hat, the projection of the relaxed bigM solution x* onto the hull,
which is perpendicular to the vector from x* to x_hat.
Note that this method will often lead to numerical difficulties since both
x* and x_hat are solutions to optimization problems. To mitigate this,
use some method of backing off the cut to make it a bit more conservative.
Parameters
-----------
transBlock_rHull: transformation block on relaxed hull instance. Ignored by
this callback.
var_info: List of tuples (rBigM_var, rHull_var, xstar_param)
hull_to_bigm_map: For expression substitution, maps id(hull_var) to
corresponding bigm var. Ignored by this callback
rBigM_linear_constraints: list of linear constraints in relaxed bigM.
Ignored by this callback.
rHull_vars: list of all variables in relaxed hull. Ignored by this callback.
disaggregated_vars: ComponentSet of disaggregated variables in hull
reformulation. Ignored by this callback
norm: The norm used in the separation problem, will be used to calculate
the subgradient used to generate the cut
cut_threshold: Amount x* needs to be infeasible in generated cut in order
to consider the cut for addition to the bigM model.
zero_tolerance: Tolerance at which a float will be treated as 0 during
Fourier-Motzkin elimination. Ignored by this callback
integer_arithmetic: Ignored by this callback (specifies FME use integer
arithmetic)
constraint_tolerance: Ignored by this callback (specifies when constraints
are considered tight in FME)
"""
cutexpr = 0
if norm == 2:
for x_rbigm, x_hull, x_star in var_info:
cutexpr += (x_hull.value - x_star.value) * (x_rbigm - x_hull.value)
elif norm == float('inf'):
duals = transBlock_rHull.model().dual
if len(duals) == 0:
raise GDP_Error(
"No dual information in the separation problem! "
"To use the infinity norm and the "
"create_cuts_normal_vector method, you must use "
"a solver which provides dual information."
)
i = 0
for x_rbigm, x_hull, x_star in var_info:
# ESJ: We wrote this so duals will be nonnegative
mu_plus = value(duals[transBlock_rHull.inf_norm_linearization[i]])
mu_minus = value(duals[transBlock_rHull.inf_norm_linearization[i + 1]])
assert mu_plus >= 0
assert mu_minus >= 0
cutexpr += (mu_plus - mu_minus) * (x_rbigm - x_hull.value)
i += 2
# make sure we're cutting off x* by enough.
if value(cutexpr) < -cut_threshold:
return [cutexpr >= 0]
logger.warning(
"Generated cut did not remove relaxed BigM solution by more "
"than the specified threshold of %s. Stopping cut "
"generation." % cut_threshold
)
return None
[docs]
def back_off_constraint_with_calculated_cut_violation(
cut, transBlock_rHull, bigm_to_hull_map, opt, stream_solver, TOL
):
"""Calculates the maximum violation of cut subject to the relaxed hull
constraints. Increases this violation by TOL (to account for optimality
tolerance in solving the problem), and, if it finds that cut can be violated
up to this tolerance, makes it more conservative such that it no longer can.
Parameters
----------
cut: The cut to be made more conservative, a Constraint
transBlock_rHull: the relaxed hull model's transformation Block
bigm_to_hull_map: Dictionary mapping ids of bigM variables to the
corresponding variables on the relaxed hull instance
opt: SolverFactory object for solving the maximum violation problem
stream_solver: Whether or not to set tee=True while solving the maximum
violation problem.
TOL: An absolute tolerance to be added to the calculated cut violation,
to account for optimality tolerance in the maximum violation problem
solve.
"""
instance_rHull = transBlock_rHull.model()
logger.info("Post-processing cut: %s" % cut.expr)
# Take a constraint. We will solve a problem maximizing its violation
# subject to rHull. We will add some user-specified tolerance to that
# violation, and then add that much padding to it if it can be violated.
transBlock_rHull.separation_objective.deactivate()
transBlock_rHull.infeasibility_objective = Objective(
expr=clone_without_expression_components(cut.body, substitute=bigm_to_hull_map)
)
results = opt.solve(instance_rHull, tee=stream_solver, load_solutions=False)
if verify_successful_solve(results) is not NORMAL:
logger.warning(
"Problem to determine how much to "
"back off the new cut "
"did not solve normally. Leaving the constraint as is, "
"which could lead to numerical trouble%s" % (results,)
)
# restore the objective
transBlock_rHull.del_component(transBlock_rHull.infeasibility_objective)
transBlock_rHull.separation_objective.activate()
return
instance_rHull.solutions.load_from(results)
# we're minimizing, val is <= 0
val = value(transBlock_rHull.infeasibility_objective) - TOL
if val <= 0:
logger.info("\tBacking off cut by %s" % val)
lb, body, ub = cut.to_bounded_expression()
cut.set_value((lb, body + abs(val), ub))
# else there is nothing to do: restore the objective
transBlock_rHull.del_component(transBlock_rHull.infeasibility_objective)
transBlock_rHull.separation_objective.activate()
[docs]
def back_off_constraint_by_fixed_tolerance(
cut, transBlock_rHull, bigm_to_hull_map, opt, stream_solver, TOL
):
"""Makes cut more conservative by absolute tolerance TOL
Parameters
----------
cut: the cut to be made more conservative, a Constraint
transBlock_rHull: the relaxed hull model's transformation Block. Ignored by
this callback
bigm_to_hull_map: Dictionary mapping ids of bigM variables to the
corresponding variables on the relaxed hull instance.
Ignored by this callback.
opt: SolverFactory object. Ignored by this callback
stream_solver: Whether or not to set tee=True while solving. Ignored by
this callback
TOL: An absolute tolerance to be added to make cut more conservative.
"""
lb, body, ub = cut.to_bounded_expression()
cut.set_value((lb, body + TOL, ub))
[docs]
@TransformationFactory.register(
'gdp.cuttingplane',
doc="Relaxes a linear disjunctive model by "
"adding cuts from convex hull to Big-M "
"reformulation.",
)
class CuttingPlane_Transformation(Transformation):
"""Relax convex disjunctive model by forming the bigm relaxation and then
iteratively adding cuts from the hull relaxation (or the hull relaxation
after some basic steps) in order to strengthen the formulation.
Note that gdp.cuttingplane is not a structural transformation: If variables
on the model are fixed, they will be treated as data, and unfixing them
after transformation will very likely result in an invalid model.
This transformation accepts the following keyword arguments:
Parameters
----------
solver : Solver name (as string) to use to solve relaxed BigM and separation
problems
solver_options : dictionary of options to pass to the solver
stream_solver : Whether or not to display solver output
verbose : Enable verbose output from cuttingplanes algorithm
cuts_name : Optional name for the IndexedConstraint containing the projected
cuts (must be a unique name with respect to the instance)
minimum_improvement_threshold : Stopping criterion based on improvement in
Big-M relaxation. This is the minimum
difference in relaxed BigM objective
values between consecutive iterations
separation_objective_threshold : Stopping criterion based on separation
objective. If separation objective is not
at least this large, cut generation will
terminate.
cut_filtering_threshold : Stopping criterion based on effectiveness of the
generated cut: This is the amount by which
a cut must be violated at the relaxed bigM
solution in order to be added to the bigM model
max_number_of_cuts : The maximum number of cuts to add to the big-M model
norm : norm to use in the objective of the separation problem
tighten_relaxation : callback to modify the GDP model before the hull
relaxation is taken (e.g. could be used to perform
basic steps)
create_cuts : callback to create cuts using the solved relaxed bigM and hull
problems
post_process_cut : callback to perform post-processing on created cuts
back_off_problem_tolerance : tolerance to use while post-processing
zero_tolerance : Tolerance at which a float will be considered 0 when
using Fourier-Motzkin elimination to create cuts.
do_integer_arithmetic : Whether or not to require Fourier-Motzkin elimination
to do integer arithmetic. Only possible when all
data is integer.
tight_constraint_tolerance : Tolerance at which a constraint is considered
tight for the Fourier-Motzkin cut generation
procedure
By default, the callbacks will be set such that the algorithm performed is
as presented in [1], but with an additional post-processing procedure to
reduce numerical error, which calculates the maximum violation of the cut
subject to the relaxed hull constraints, and then pads the constraint by
this violation plus an additional user-specified tolerance.
In addition, the create_cuts_fme function provides an (exponential time)
method of generating cuts which reduces numerical error (and can eliminate
it if all data is integer). It collects the hull constraints which are
tight at the solution of the separation problem. It creates a cut in the
extended space perpendicular to a composite normal vector created by
summing the directions normal to these constraints. It then performs
fourier-motzkin elimination on the collection of constraints and the cut
to project out the disaggregated variables. The resulting constraint which
is most violated by the relaxed bigM solution is then returned.
References
----------
[1] Sawaya, N. W., Grossmann, I. E. (2005). A cutting plane method for
solving linear generalized disjunctive programming problems. Computers
and Chemical Engineering, 29, 1891-1913
"""
CONFIG = ConfigBlock("gdp.cuttingplane")
CONFIG.declare(
'solver',
ConfigValue(
default='ipopt',
domain=str,
description="""Solver to use for relaxed BigM problem and the separation
problem""",
doc="""
This specifies the solver which will be used to solve LP relaxation
of the BigM problem and the separation problem. Note that this solver
must be able to handle a quadratic objective because of the separation
problem.
""",
),
)
CONFIG.declare(
'minimum_improvement_threshold',
ConfigValue(
default=0.01,
domain=NonNegativeFloat,
description="Threshold value for difference in relaxed bigM problem "
"objectives used to decide when to stop adding cuts",
doc="""
If the difference between the objectives in two consecutive iterations is
less than this value, the algorithm terminates without adding the cut
generated in the last iteration.
""",
),
)
CONFIG.declare(
'separation_objective_threshold',
ConfigValue(
default=0.01,
domain=NonNegativeFloat,
description="Threshold value used to decide when to stop adding cuts: "
"If separation problem objective is not at least this quantity, cut "
"generation will terminate.",
doc="""
If the separation problem objective (distance between relaxed bigM
solution and its projection onto the relaxed hull feasible region)
does not exceed this threshold, the algorithm will terminate.
""",
),
)
CONFIG.declare(
'max_number_of_cuts',
ConfigValue(
default=100,
domain=PositiveInt,
description="The maximum number of cuts to add before the algorithm "
"terminates.",
doc="""
If the algorithm does not terminate due to another criterion first,
cut generation will stop after adding this many cuts.
""",
),
)
CONFIG.declare(
'norm',
ConfigValue(
default=2,
domain=In([2, float('inf')]),
description="Norm to use in the separation problem: 2, or float('inf')",
doc="""
Norm used to calculate distance in the objective of the separation
problem which finds the nearest point on the hull relaxation region
to the current solution of the relaxed bigm problem.
Supported norms are the Euclidean norm (specify 2) and the infinity
norm (specify float('inf')). Note that the first makes the separation
problem objective quadratic and the latter makes it linear.
""",
),
)
CONFIG.declare(
'verbose',
ConfigValue(
default=False,
domain=bool,
description="Flag to enable verbose output",
doc="""
If True, prints subproblem solutions, as well as potential and added cuts
during algorithm.
If False, only the relaxed BigM objective and minimal information about
cuts is logged.
""",
),
)
CONFIG.declare(
'stream_solver',
ConfigValue(
default=False,
domain=bool,
description="""If true, sets tee=True for every solve performed over
"the course of the algorithm""",
),
)
CONFIG.declare(
'solver_options',
ConfigBlock(
implicit=True,
description="Dictionary of solver options",
doc="""
Dictionary of solver options that will be set for the solver for both the
relaxed BigM and separation problem solves.
""",
),
)
CONFIG.declare(
'tighten_relaxation',
ConfigValue(
default=do_not_tighten,
description="Callback which takes the GDP formulation and returns a "
"GDP formulation with a tighter hull relaxation",
doc="""
Function which accepts the GDP formulation of the problem and returns
a GDP formulation which the transformation will then take the hull
reformulation of.
Most typically, this callback would be used to apply basic steps before
taking the hull reformulation, but anything which tightens the GDP can
be performed here.
""",
),
)
CONFIG.declare(
'create_cuts',
ConfigValue(
default=create_cuts_normal_vector,
description="Callback which generates a list of cuts, given the solved "
"relaxed bigM and relaxed hull solutions. If no cuts can be "
"generated, returns None",
doc="""
Callback to generate cuts to be added to the bigM problem based on
solutions to the relaxed bigM problem and the separation problem.
Arguments
---------
transBlock_rBigm: transformation block on relaxed bigM instance
transBlock_rHull: transformation block on relaxed hull instance
var_info: List of tuples (rBigM_var, rHull_var, xstar_param)
hull_to_bigm_map: For expression substitution, maps id(hull_var) to
corresponding bigm var
rBigM_linear_constraints: list of linear constraints in relaxed bigM
rHull_vars: list of all variables in relaxed hull
disaggregated_vars: ComponentSet of disaggregated variables in hull
reformulation
cut_threshold: Amount x* needs to be infeasible in generated cut in order
to consider the cut for addition to the bigM model.
zero_tolerance: Tolerance at which a float will be treated as 0
Returns
-------
list of cuts to be added to bigM problem (and relaxed bigM problem),
represented as expressions using variables from the bigM model
""",
),
)
CONFIG.declare(
'post_process_cut',
ConfigValue(
default=back_off_constraint_with_calculated_cut_violation,
description="Callback which takes a generated cut and post processes "
"it, presumably to back it off in the case of numerical error. Set to "
"None if not post-processing is desired.",
doc="""
Callback to adjust a cut returned from create_cuts before adding it to
the model, presumably to make it more conservative in case of numerical
error.
Arguments
---------
cut: the cut to be made more conservative, a Constraint
transBlock_rHull: the relaxed hull model's transformation Block.
bigm_to_hull_map: Dictionary mapping ids of bigM variables to the
corresponding variables on the relaxed hull instance.
opt: SolverFactory object for subproblem solves in this procedure
stream_solver: Whether or not to set tee=True while solving.
TOL: A tolerance
Returns
-------
None, modifies the cut in place
""",
),
)
# back off problem tolerance (on top of the solver's (sometimes))
CONFIG.declare(
'back_off_problem_tolerance',
ConfigValue(
default=1e-8,
domain=NonNegativeFloat,
description="Tolerance to pass to the post_process_cut callback.",
doc="""
Tolerance passed to the post_process_cut callback.
Depending on the callback, different values could make sense, but
something on the order of the solver's optimality or constraint
tolerances is appropriate.
""",
),
)
CONFIG.declare(
'cut_filtering_threshold',
ConfigValue(
default=0.001,
domain=NonNegativeFloat,
description="Tolerance used to decide if a cut removes x* from the "
"relaxed BigM problem by enough to be added to the bigM problem.",
doc="""
Absolute tolerance used to decide whether to keep a cut. We require
that, when evaluated at x* (the relaxed BigM optimal solution), the
cut be infeasible by at least this tolerance.
""",
),
)
CONFIG.declare(
'zero_tolerance',
ConfigValue(
default=1e-9,
domain=NonNegativeFloat,
description="Tolerance at which floats are assumed to be 0 while "
"performing Fourier-Motzkin elimination",
doc="""
Only relevant when create_cuts=create_cuts_fme, this sets the
zero_tolerance option for the Fourier-Motzkin elimination transformation.
""",
),
)
CONFIG.declare(
'do_integer_arithmetic',
ConfigValue(
default=False,
domain=bool,
description="Only relevant if using Fourier-Motzkin Elimination (FME) "
"and if all problem data is integer, requires FME transformation to "
"perform integer arithmetic to eliminate numerical error.",
doc="""
Only relevant when create_cuts=create_cuts_fme and if all problem data
is integer, this sets the do_integer_arithmetic flag to true for the
FME transformation, meaning that the projection to the big-M space
can be done with exact precision.
""",
),
)
CONFIG.declare(
'cuts_name',
ConfigValue(
default=None,
domain=str,
description="Optional name for the IndexedConstraint containing the "
"projected cuts. Must be a unique name with respect to the "
"instance.",
doc="""
Optional name for the IndexedConstraint containing the projected
constraints. If not specified, the cuts will be stored on a
private block created by the transformation, so if you want access
to them after the transformation, use this argument.
Must be a string which is a unique component name with respect to the
Block on which the transformation is called.
""",
),
)
CONFIG.declare(
'tight_constraint_tolerance',
ConfigValue(
default=1e-6, # Gurobi constraint tolerance
domain=NonNegativeFloat,
description="Tolerance at which a constraint is considered tight for "
"the Fourier-Motzkin cut generation procedure.",
doc="""
For a constraint a^Tx <= b, the Fourier-Motzkin cut generation procedure
will consider the constraint tight (and add it to the set of constraints
being projected) when a^Tx - b is less than this tolerance.
It is recommended to set this tolerance to the constraint tolerance of
the solver being used.
""",
),
)
def _apply_to(self, instance, bigM=None, **kwds):
original_log_level = logger.level
log_level = logger.getEffectiveLevel()
try:
self._config = self.CONFIG(kwds.pop('options', {}))
self._config.set_value(kwds)
if self._config.verbose and log_level > logging.INFO:
logger.setLevel(logging.INFO)
self.verbose = True
elif log_level <= logging.INFO:
self.verbose = True
else:
self.verbose = False
(instance_rBigM, cuts_obj, instance_rHull, var_info, transBlockName) = (
self._setup_subproblems(instance, bigM, self._config.tighten_relaxation)
)
self._generate_cuttingplanes(
instance_rBigM, cuts_obj, instance_rHull, var_info, transBlockName
)
# restore integrality
TransformationFactory('core.relax_integer_vars').apply_to(
instance, undo=True
)
finally:
del self._config
del self.verbose
# restore logging level
logger.setLevel(original_log_level)
def _setup_subproblems(self, instance, bigM, tighten_relaxation_callback):
# create transformation block
transBlockName, transBlock = self._add_transformation_block(instance)
# We store a list of all vars so that we can efficiently
# generate maps among the subproblems
transBlock.all_vars = list(
v
for v in instance.component_data_objects(
Var, descend_into=(Block, Disjunct), sort=SortComponents.deterministic
)
if not v.is_fixed()
)
# we'll store all the cuts we add together
nm = self._config.cuts_name
if nm is None:
cuts_obj = transBlock.cuts = Constraint(NonNegativeIntegers)
else:
# check that this really is an available name
if instance.component(nm) is not None:
raise GDP_Error(
"cuts_name was specified as '%s', but this is "
"already a component on the instance! Please "
"specify a unique name." % nm
)
instance.add_component(nm, Constraint(NonNegativeIntegers))
cuts_obj = instance.component(nm)
# get bigM and hull relaxations
bigMRelaxation = TransformationFactory('gdp.bigm')
hullRelaxation = TransformationFactory('gdp.hull')
relaxIntegrality = TransformationFactory('core.relax_integer_vars')
#
# Generate the Hull relaxation (used for the separation
# problem to generate cutting planes)
#
tighter_instance = tighten_relaxation_callback(instance)
instance_rHull = hullRelaxation.create_using(tighter_instance)
relaxIntegrality.apply_to(instance_rHull, transform_deactivated_blocks=True)
#
# Reformulate the instance using the BigM relaxation (this will
# be the final instance returned to the user)
#
bigMRelaxation.apply_to(instance, bigM=bigM)
#
# Generate the continuous relaxation of the BigM transformation. We'll
# restore it at the end.
#
relaxIntegrality.apply_to(instance, transform_deactivated_blocks=True)
#
# Add the xstar parameter for the Hull problem
#
transBlock_rHull = instance_rHull.component(transBlockName)
#
# this will hold the solution to rbigm each time we solve it. We
# add it to the transformation block so that we don't have to
# worry about name conflicts.
transBlock_rHull.xstar = Param(
range(len(transBlock.all_vars)), mutable=True, default=0, within=Reals
)
# we will add a block that we will deactivate to use to store the
# extended space cuts. We never need to solve these, but we need them to
# be constructed for the sake of Fourier-Motzkin Elimination
extendedSpaceCuts = transBlock_rHull.extendedSpaceCuts = Block()
extendedSpaceCuts.deactivate()
extendedSpaceCuts.cuts = Constraint(Any)
#
# Generate the mapping between the variables on all the
# instances and the xstar parameter.
#
var_info = [
(
v, # this is the bigM variable
transBlock_rHull.all_vars[i],
transBlock_rHull.xstar[i],
)
for i, v in enumerate(transBlock.all_vars)
]
# NOTE: we wait to add the separation objective to the rHull problem
# because it is best to do it in the first iteration, so that we can
# skip stale variables.
return (instance, cuts_obj, instance_rHull, var_info, transBlockName)
# this is the map that I need to translate my projected cuts and add
# them to bigM
def _create_hull_to_bigm_substitution_map(self, var_info):
return dict((id(var_info[i][1]), var_info[i][0]) for i in range(len(var_info)))
# this map is needed to solve the back-off problem for post-processing
def _create_bigm_to_hull_substitution_map(self, var_info):
return dict((id(var_info[i][0]), var_info[i][1]) for i in range(len(var_info)))
def _get_disaggregated_vars(self, hull):
disaggregatedVars = ComponentSet()
for disjunction in hull.component_data_objects(
Disjunction, descend_into=(Disjunct, Block)
):
for disjunct in disjunction.disjuncts:
transBlock = disjunct.transformation_block
if transBlock is not None:
for v in transBlock.disaggregatedVars.component_data_objects(Var):
disaggregatedVars.add(v)
return disaggregatedVars
def _get_rBigM_obj_and_constraints(self, instance_rBigM):
# We try to grab the first active objective. If there is more
# than one, the writer will yell when we try to solve below. If
# there are 0, we will yell here.
rBigM_obj = next(
instance_rBigM.component_data_objects(Objective, active=True), None
)
if rBigM_obj is None:
raise GDP_Error(
"Cannot apply cutting planes transformation "
"without an active objective in the model!"
)
#
# Collect all of the linear constraints that are in the rBigM
# instance. We will need these so that we can compare what we get from
# FME to them and make sure we aren't adding redundant constraints to
# the model. For convenience, we will make sure they are all in the form
# lb <= expr (so we will break equality constraints)
#
fme = TransformationFactory('contrib.fourier_motzkin_elimination')
rBigM_linear_constraints = []
for cons in instance_rBigM.component_data_objects(
Constraint,
descend_into=Block,
sort=SortComponents.deterministic,
active=True,
):
body = cons.body
if body.polynomial_degree() != 1:
# We will never get a nonlinear constraint out of FME, so we
# don't risk it being identical to this one.
continue
# TODO: Guess this shouldn't have been private...
rBigM_linear_constraints.extend(fme._process_constraint(cons))
# [ESJ Aug 13 2020] NOTE: We actually don't need to worry about variable
# bounds here because the FME transformation will take care of them
# (i.e. convert those of the disaggregated variables to constraints for
# the purposes of the projection.)
return rBigM_obj, rBigM_linear_constraints
def _generate_cuttingplanes(
self, instance_rBigM, cuts_obj, instance_rHull, var_info, transBlockName
):
opt = SolverFactory(self._config.solver)
stream_solver = self._config.stream_solver
opt.options = dict(self._config.solver_options)
improving = True
prev_obj = None
epsilon = self._config.minimum_improvement_threshold
cuts = None
transBlock_rHull = instance_rHull.component(transBlockName)
rBigM_obj, rBigM_linear_constraints = self._get_rBigM_obj_and_constraints(
instance_rBigM
)
# Get list of all variables in the rHull model which we will use when
# calculating the composite normal vector.
rHull_vars = [
i
for i in instance_rHull.component_data_objects(
Var, descend_into=Block, sort=SortComponents.deterministic
)
]
# collect a list of disaggregated variables.
disaggregated_vars = self._get_disaggregated_vars(instance_rHull)
hull_to_bigm_map = self._create_hull_to_bigm_substitution_map(var_info)
bigm_to_hull_map = self._create_bigm_to_hull_substitution_map(var_info)
xhat = ComponentMap()
while improving:
# solve rBigM, solution is xstar
results = opt.solve(instance_rBigM, tee=stream_solver, load_solutions=False)
if verify_successful_solve(results) is not NORMAL:
logger.warning(
"Relaxed BigM subproblem "
"did not solve normally. Stopping cutting "
"plane generation.\n\n%s" % (results,)
)
return
instance_rBigM.solutions.load_from(results)
rBigM_objVal = value(rBigM_obj)
logger.warning("rBigM objective = %s" % (rBigM_objVal,))
#
# Add the separation objective to the hull subproblem if it's not
# there already (so in the first iteration). We're waiting until now
# to avoid it including variables that came back stale from the
# rbigm solve.
#
if transBlock_rHull.component("separation_objective") is None:
self._add_separation_objective(var_info, transBlock_rHull)
# copy over xstar
logger.info("x* is:")
for x_rbigm, x_hull, x_star in var_info:
if not x_rbigm.stale:
x_star.set_value(x_rbigm.value)
# initialize the X values
x_hull.set_value(x_rbigm.value, skip_validation=True)
if self.verbose:
logger.info(
"\t%s = %s"
% (x_rbigm.getname(fully_qualified=True), x_rbigm.value)
)
# compare objectives: check absolute difference close to 0, relative
# difference further from 0.
if prev_obj is None:
improving = True
else:
obj_diff = prev_obj - rBigM_objVal
improving = (
abs(obj_diff) > epsilon
if abs(rBigM_objVal) < 1
else abs(obj_diff / prev_obj) > epsilon
)
# solve separation problem to get xhat.
results = opt.solve(instance_rHull, tee=stream_solver, load_solutions=False)
if verify_successful_solve(results) is not NORMAL:
logger.warning(
"Hull separation subproblem "
"did not solve normally. Stopping cutting "
"plane generation.\n\n%s" % (results,)
)
return
instance_rHull.solutions.load_from(results)
logger.warning(
"separation problem objective value: %s"
% value(transBlock_rHull.separation_objective)
)
# save xhat to initialize rBigM with in the next iteration
if self.verbose:
logger.info("xhat is: ")
for x_rbigm, x_hull, x_star in var_info:
xhat[x_rbigm] = value(x_hull)
if self.verbose:
logger.info(
"\t%s = %s"
% (x_hull.getname(fully_qualified=True), x_hull.value)
)
# [JDS 19 Dec 18] Note: we check that the separation objective was
# significantly nonzero. If it is too close to zero, either the
# rBigM solution was in the convex hull, or the separation vector is
# so close to zero that the resulting cut is likely to have
# numerical issues.
if (
value(transBlock_rHull.separation_objective)
< self._config.separation_objective_threshold
):
logger.warning(
"Separation problem objective below threshold of"
" %s: Stopping cut generation."
% self._config.separation_objective_threshold
)
break
cuts = self._config.create_cuts(
transBlock_rHull,
var_info,
hull_to_bigm_map,
rBigM_linear_constraints,
rHull_vars,
disaggregated_vars,
self._config.norm,
self._config.cut_filtering_threshold,
self._config.zero_tolerance,
self._config.do_integer_arithmetic,
self._config.tight_constraint_tolerance,
)
# We are done if the cut generator couldn't return a valid cut
if cuts is None:
logger.warning("Did not generate a valid cut, stopping cut generation.")
break
if not improving:
logger.warning(
"Difference in relaxed BigM problem objective "
"values from past two iterations is below "
"threshold of %s: Stopping cut generation." % epsilon
)
break
for cut in cuts:
# we add the cut to the model and then post-process it in place.
cut_number = len(cuts_obj)
logger.warning("Adding cut %s to BigM model." % (cut_number,))
cuts_obj.add(cut_number, cut)
if self._config.post_process_cut is not None:
self._config.post_process_cut(
cuts_obj[cut_number],
transBlock_rHull,
bigm_to_hull_map,
opt,
stream_solver,
self._config.back_off_problem_tolerance,
)
if cut_number + 1 == self._config.max_number_of_cuts:
logger.warning("Reached maximum number of cuts.")
break
prev_obj = rBigM_objVal
# Initialize rbigm with xhat (for the next iteration)
for x_rbigm, x_hull, x_star in var_info:
x_rbigm.set_value(xhat[x_rbigm], skip_validation=True)
def _add_transformation_block(self, instance):
# creates transformation block with a unique name based on name, adds it
# to instance, and returns it.
transBlockName = unique_component_name(
instance, '_pyomo_gdp_cuttingplane_transformation'
)
transBlock = Block()
instance.add_component(transBlockName, transBlock)
return transBlockName, transBlock
def _add_separation_objective(self, var_info, transBlock_rHull):
# creates the separation objective. That is just adding an objective for
# Euclidean norm, it means adding an auxiliary variable to linearize the
# L-infinity norm. We do this assuming that rBigM has been solved, and
# if any variables come back stale, we leave them out of the separation
# problem, as they aren't doing anything and they could cause numerical
# issues later.
# Deactivate any/all other objectives
for o in transBlock_rHull.model().component_data_objects(Objective):
o.deactivate()
norm = self._config.norm
to_delete = []
if norm == 2:
obj_expr = 0
for i, (x_rbigm, x_hull, x_star) in enumerate(var_info):
if not x_rbigm.stale:
obj_expr += (x_hull - x_star) ** 2
else:
if self.verbose:
logger.info(
"The variable %s will not be included in "
"the separation problem: It was stale in "
"the rBigM solve." % x_rbigm.getname(fully_qualified=True)
)
to_delete.append(i)
elif norm == float('inf'):
u = transBlock_rHull.u = Var(domain=NonNegativeReals)
inf_cons = transBlock_rHull.inf_norm_linearization = Constraint(
NonNegativeIntegers
)
i = 0
for j, (x_rbigm, x_hull, x_star) in enumerate(var_info):
if not x_rbigm.stale:
# NOTE: these are written as >= constraints so that we know
# the duals will come back nonnegative.
inf_cons[i] = u - x_hull >= -x_star
inf_cons[i + 1] = u + x_hull >= x_star
i += 2
else:
if self.verbose:
logger.info(
"The variable %s will not be included in "
"the separation problem: It was stale in "
"the rBigM solve." % x_rbigm.getname(fully_qualified=True)
)
to_delete.append(j)
# we'll need the duals of these to get the subgradient
self._add_dual_suffix(transBlock_rHull.model())
obj_expr = u
# delete the unneeded x_stars so that we don't add cuts involving
# useless variables later.
for i in sorted(to_delete, reverse=True):
del var_info[i]
# add separation objective to transformation block
transBlock_rHull.separation_objective = Objective(expr=obj_expr)
def _add_dual_suffix(self, rHull):
# rHull is our model and we aren't giving it back (unless in the future
# we we add a callback to do basic steps to it...), so we just check if
# dual is there. If it's a Suffix, we'll borrow it. If it's something
# else we'll rename it and add the Suffix.
dual = rHull.component("dual")
if dual is None:
rHull.dual = Suffix(direction=Suffix.IMPORT)
else:
if dual.ctype is Suffix:
return
rHull.del_component(dual)
rHull.dual = Suffix(direction=Suffix.IMPORT)
rHull.add_component(unique_component_name(rHull, "dual"), dual)