Source code for pyomo.contrib.preprocessing.plugins.induced_linearity

#  ___________________________________________________________________________
#  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.
#  ___________________________________________________________________________

"""Transformation to reformulate nonlinear models with linearity induced from
discrete variables.

Ref: Grossmann, IE; Voudouris, VT; Ghattas, O. Mixed integer linear
reformulations for some nonlinear discrete design optimization problems.


from __future__ import division

import logging
import textwrap
from math import fabs

from pyomo.common.collections import ComponentMap, ComponentSet
from pyomo.common.config import (ConfigBlock, ConfigValue, NonNegativeFloat,
from pyomo.common.modeling import unique_component_name
from pyomo.contrib.preprocessing.util import SuppressConstantObjectiveWarning
from pyomo.core import (Binary, Block, Constraint, Objective, Set,
                        TransformationFactory, Var, summation, value)
from pyomo.core.plugins.transform.hierarchy import IsomorphicTransformation
from pyomo.gdp import Disjunct
from pyomo.opt import TerminationCondition as tc
from pyomo.opt import SolverFactory
from pyomo.repn import generate_standard_repn

logger = logging.getLogger('pyomo.contrib.preprocessing')

[docs]@TransformationFactory.register('contrib.induced_linearity', doc="Reformulate nonlinear constraints with induced linearity.") class InducedLinearity(IsomorphicTransformation): """Reformulate nonlinear constraints with induced linearity. Finds continuous variables :math:`v` where :math:`v = d_1 + d_2 + d_3`, where :math:`d`'s are discrete variables. These continuous variables may participate nonlinearly in other expressions, which may then be induced to be linear. The overall algorithm flow can be summarized as: 1. Detect effectively discrete variables and the constraints that imply discreteness. 2. Determine the set of valid values for each effectively discrete variable 3. Find nonlinear expressions in which effectively discrete variables participate. 4. Reformulate nonlinear expressions appropriately. .. note:: Tasks 1 & 2 must incorporate scoping considerations (Disjuncts) Keyword arguments below are specified for the ``apply_to`` and ``create_using`` functions. """ CONFIG = ConfigBlock("contrib.induced_linearity") CONFIG.declare('equality_tolerance', ConfigValue( default=1E-6, domain=NonNegativeFloat, description="Tolerance on equality constraints." )) CONFIG.declare('pruning_solver', ConfigValue( default='glpk', description="Solver to use when pruning possible values." )) __doc__ = add_docstring_list(__doc__, CONFIG) def _apply_to(self, model, **kwds): """Apply the transformation to the given model.""" config = self.CONFIG(kwds.pop('options', {})) config.set_value(kwds) _process_container(model, config) _process_subcontainers(model, config)
def _process_subcontainers(blk, config): for disj in blk.component_data_objects( Disjunct, active=True, descend_into=True): _process_container(disj, config) _process_subcontainers(disj, config) def _process_container(blk, config): if not hasattr(blk, '_induced_linearity_info'): blk._induced_linearity_info = Block() else: assert blk._induced_linearity_info.ctype == Block eff_discr_vars = detect_effectively_discrete_vars( blk, config.equality_tolerance) # TODO will need to go through this for each disjunct, since it does # not (should not) descend into Disjuncts. # Determine the valid values for the effectively discrete variables possible_var_values = determine_valid_values(blk, eff_discr_vars, config) # Collect find bilinear expressions that can be reformulated using # knowledge of effectively discrete variables bilinear_map = _bilinear_expressions(blk) # Relevant constraints are those with bilinear terms that involve # effectively_discrete_vars processed_pairs = ComponentSet() for v1, var_values in possible_var_values.items(): v1_pairs = bilinear_map.get(v1, ()) for v2, bilinear_constrs in v1_pairs.items(): if (v1, v2) in processed_pairs: continue _process_bilinear_constraints( blk, v1, v2, var_values, bilinear_constrs) processed_pairs.add((v2, v1)) # processed_pairs.add((v1, v2)) # TODO is this necessary? def determine_valid_values(block, discr_var_to_constrs_map, config): """Calculate valid values for each effectively discrete variable. We need the set of possible values for the effectively discrete variable in order to do the reformulations. Right now, we select a naive approach where we look for variables in the discreteness-inducing constraints. We then adjust their values and see if things are stil feasible. Based on their coefficient values, we can infer a set of allowable values for the effectively discrete variable. Args: block: The model or a disjunct on the model. """ possible_values = ComponentMap() for eff_discr_var, constrs in discr_var_to_constrs_map.items(): # get the superset of possible values by looking through the # constraints for constr in constrs: repn = generate_standard_repn(constr.body) var_coef = sum(coef for i, coef in enumerate(repn.linear_coefs) if repn.linear_vars[i] is eff_discr_var) const = -(repn.constant - constr.upper) / var_coef possible_vals = set((const,)) for i, var in enumerate(repn.linear_vars): if var is eff_discr_var: continue coef = -repn.linear_coefs[i] / var_coef if var.is_binary(): var_values = (0, coef) elif var.is_integer(): var_values = [v * coef for v in range(, var.ub + 1)] else: raise ValueError( '%s has unacceptable variable domain: %s' % (, var.domain)) possible_vals = set( (v1 + v2 for v1 in possible_vals for v2 in var_values)) old_possible_vals = possible_values.get(eff_discr_var, None) if old_possible_vals is not None: possible_values[eff_discr_var] = old_possible_vals & possible_vals else: possible_values[eff_discr_var] = possible_vals possible_values = prune_possible_values(block, possible_values, config) return possible_values def prune_possible_values(block_scope, possible_values, config): # Prune the set of possible values by solving a series of feasibility # problems top_level_scope = block_scope.model() tmp_name = unique_component_name( top_level_scope, '_induced_linearity_prune_data') tmp_orig_blk = Block() setattr(top_level_scope, tmp_name, tmp_orig_blk) tmp_orig_blk._possible_values = possible_values tmp_orig_blk._possible_value_vars = list(v for v in possible_values) tmp_orig_blk._tmp_block_scope = (block_scope,) model = top_level_scope.clone() tmp_clone_blk = getattr(model, tmp_name) for obj in model.component_data_objects(Objective, active=True): obj.deactivate() for constr in model.component_data_objects( Constraint, active=True, descend_into=(Block, Disjunct)): if constr.body.polynomial_degree() not in (1, 0): constr.deactivate() if block_scope.ctype == Disjunct: disj = tmp_clone_blk._tmp_block_scope[0] disj.indicator_var.fix(1) TransformationFactory('gdp.bigm').apply_to(model) tmp_clone_blk.test_feasible = Constraint() tmp_clone_blk._obj = Objective(expr=1) for eff_discr_var, vals in tmp_clone_blk._possible_values.items(): val_feasible = {} for val in vals: tmp_clone_blk.test_feasible.set_value(eff_discr_var == val) with SuppressConstantObjectiveWarning(): res = SolverFactory(config.pruning_solver).solve(model) if res.solver.termination_condition is tc.infeasible: val_feasible[val] = False tmp_clone_blk._possible_values[eff_discr_var] = set( v for v in tmp_clone_blk._possible_values[eff_discr_var] if val_feasible.get(v, True)) for i, var in enumerate(tmp_orig_blk._possible_value_vars): possible_values[var] = tmp_clone_blk._possible_values[ tmp_clone_blk._possible_value_vars[i]] return possible_values def _process_bilinear_constraints(block, v1, v2, var_values, bilinear_constrs): # TODO check that the appropriate variable bounds exist. if not (v2.has_lb() and v2.has_ub()): logger.warning(textwrap.dedent("""\ Attempting to transform bilinear term {v1} * {v2} using effectively discrete variable {v1}, but {v2} is missing a lower or upper bound: ({v2lb}, {v2ub}). """.format(v1=v1, v2=v2,, v2ub=v2.ub)).strip()) return False blk = Block() unique_name = unique_component_name( block, ("%s_%s_bilinear" % (v1.local_name, v2.local_name)) .replace('[', '').replace(']', '')) block._induced_linearity_info.add_component(unique_name, blk) # TODO think about not using floats as indices in a set blk.valid_values = Set(initialize=sorted(var_values)) blk.x_active = Var(blk.valid_values, domain=Binary, initialize=1) blk.v_increment = Var( blk.valid_values, domain=v2.domain, bounds=(, v2.ub), initialize=v2.value) blk.v_defn = Constraint(expr=v2 == summation(blk.v_increment)) @blk.Constraint(blk.valid_values) def v_lb(blk, val): return * blk.x_active[val] <= blk.v_increment[val] @blk.Constraint(blk.valid_values) def v_ub(blk, val): return blk.v_increment[val] <= v2.ub * blk.x_active[val] blk.select_one_value = Constraint(expr=summation(blk.x_active) == 1) # Categorize as case 1 or case 2 for bilinear_constr in bilinear_constrs: # repn = generate_standard_repn(bilinear_constr.body) # Case 1: no other variables besides bilinear term in constraint. v1 # (effectively discrete variable) is positive. # if (len(repn.quadratic_vars) == 1 and len(repn.linear_vars) == 0 # and repn.nonlinear_expr is None): # _reformulate_case_1(v1, v2, discrete_constr, bilinear_constr) # NOTE: Case 1 is left unimplemented for now, because it involves some # messier logic with respect to how the transformation needs to happen. # Case 2: this is everything else, but do we want to have a special # case if there are nonlinear expressions involved with the constraint? pass _reformulate_case_2(blk, v1, v2, bilinear_constr) pass def _reformulate_case_2(blk, v1, v2, bilinear_constr): repn = generate_standard_repn(bilinear_constr.body) replace_index = next( i for i, var_tup in enumerate(repn.quadratic_vars) if (var_tup[0] is v1 and var_tup[1] is v2) or (var_tup[0] is v2 and var_tup[1] is v1)) bilinear_constr.set_value(( bilinear_constr.lower, sum(coef * repn.linear_vars[i] for i, coef in enumerate(repn.linear_coefs)) + repn.quadratic_coefs[replace_index] * sum( val * blk.v_increment[val] for val in blk.valid_values) + sum(repn.quadratic_coefs[i] * var_tup[0] * var_tup[1] for i, var_tup in enumerate(repn.quadratic_vars) if not i == replace_index) + repn.constant + zero_if_None(repn.nonlinear_expr), bilinear_constr.upper )) def zero_if_None(val): return 0 if val is None else val def _bilinear_expressions(model): # TODO for now, we look for only expressions where the bilinearities are # exposed on the root level SumExpression, and thus accessible via # generate_standard_repn. This will not detect exp(x*y). We require a # factorization transformation to be applied beforehand in order to pick # these constraints up. pass # Bilinear map will be stored in the format: # x --> (y --> [constr1, constr2, ...], z --> [constr2, constr3]) bilinear_map = ComponentMap() for constr in model.component_data_objects( Constraint, active=True, descend_into=(Block, Disjunct)): if constr.body.polynomial_degree() in (1, 0): continue # Skip trivial and linear constraints repn = generate_standard_repn(constr.body) for pair in repn.quadratic_vars: v1, v2 = pair v1_pairs = bilinear_map.get(v1, ComponentMap()) if v2 in v1_pairs: # bilinear term has been found before. Simply add constraint to # the set associated with the bilinear term. v1_pairs[v2].add(constr) else: # We encounter the bilinear term for the first time. bilinear_map[v1] = v1_pairs bilinear_map[v2] = bilinear_map.get(v2, ComponentMap()) constraints_with_bilinear_pair = ComponentSet([constr]) bilinear_map[v1][v2] = constraints_with_bilinear_pair bilinear_map[v2][v1] = constraints_with_bilinear_pair return bilinear_map def detect_effectively_discrete_vars(block, equality_tolerance): """Detect effectively discrete variables. These continuous variables are the sum of discrete variables. """ # Map of effectively_discrete var --> inducing constraints effectively_discrete = ComponentMap() for constr in block.component_data_objects(Constraint, active=True): if constr.lower is None or constr.upper is None: continue # skip inequality constraints if fabs(value(constr.lower) - value(constr.upper) ) > equality_tolerance: continue # not equality constriant. Skip. if constr.body.polynomial_degree() not in (1, 0): continue # skip nonlinear expressions repn = generate_standard_repn(constr.body) if len(repn.linear_vars) < 2: # TODO should this be < 2 or < 1? # TODO we should make sure that trivial equality relations are # preprocessed before this, or we will end up reformulating # expressions that we do not need to here. continue non_discrete_vars = list(v for v in repn.linear_vars if v.is_continuous()) if len(non_discrete_vars) == 1: # We know that this is an effectively discrete continuous # variable. Add it to our identified variable list. var = non_discrete_vars[0] inducing_constraints = effectively_discrete.get(var, []) inducing_constraints.append(constr) effectively_discrete[var] = inducing_constraints # TODO we should eventually also look at cases where all other # non_discrete_vars are effectively_discrete_vars return effectively_discrete