# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2025
# 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')
[docs]
def do_not_tighten(m):
return m
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))