diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index f7e7d7a4b9d..0ac105fc7ba 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -16,14 +16,16 @@ from pyomo.common import deprecated from pyomo.common.collections import ComponentMap, ComponentSet, DefaultComponentMap from pyomo.common.modeling import unique_component_name +from pyomo.common.errors import DeveloperError from pyomo.core.expr.numvalue import ZeroConstant import pyomo.core.expr as EXPR -from pyomo.core.base import TransformationFactory +from pyomo.core.base import TransformationFactory, SortComponents from pyomo.core import ( Block, BooleanVar, Connector, Constraint, + ConstraintList, Param, Set, SetOf, @@ -38,6 +40,9 @@ value, NonNegativeIntegers, Binary, + ConcreteModel, + Objective, + Reference, ) from pyomo.gdp import Disjunct, Disjunction, GDP_Error from pyomo.gdp.disjunct import DisjunctData @@ -48,9 +53,21 @@ is_child_of, _warn_for_active_disjunct, ) +from pyomo.gdp.plugins.multiple_bigm import Solver from pyomo.core.util import target_list +from pyomo.core.expr.visitor import ( + IdentifyVariableVisitor, + StreamBasedExpressionVisitor, +) +from pyomo.repn.linear import LinearRepnVisitor +from pyomo.repn.util import VarRecorder from pyomo.util.vars_from_expressions import get_vars_from_components +from pyomo.opt.base.solvers import SolverFactory +from pyomo.opt.results.solver import TerminationCondition +from pyomo.opt.solver import SolverStatus from weakref import ref as weakref_ref +import math +import itertools logger = logging.getLogger('pyomo.gdp.hull') @@ -61,6 +78,7 @@ class _HullTransformationData(AutoSlots.Mixin): 'original_var_map', 'bigm_constraint_map', 'disaggregation_constraint_map', + 'well_defined_points_map', ) def __init__(self): @@ -68,6 +86,7 @@ def __init__(self): self.original_var_map = ComponentMap() self.bigm_constraint_map = DefaultComponentMap(ComponentMap) self.disaggregation_constraint_map = DefaultComponentMap(ComponentMap) + self.well_defined_points_map = ComponentMap() Block.register_private_data_initializer(_HullTransformationData) @@ -129,7 +148,7 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): domain=cfg.In(['FurmanSawayaGrossmann', 'LeeGrossmann', 'GrossmannLee']), description='perspective function used for variable disaggregation', doc=""" - The perspective function used for variable disaggregation + The perspective function used to transform nonlinear functions "LeeGrossmann" is the original NL convex hull from Lee & Grossmann (2000) [1]_, which substitutes nonlinear constraints @@ -155,6 +174,20 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): ((1-eps)*y_ik + eps) * h_ik( nu_ik/((1-eps)*y_ik + eps) ) \ - eps * h_ki(0) * ( 1-y_ik ) <= 0 + + The default, "FurmanSawayaGrossmann", is strongly + recommended. When "FurmanSawayaGrossmann" is used, any value of + epsilon in (0, 1) leads to an exact reformulation, and + decreasing epsilon improves the quality of the continuous + relaxation (see [3]_). When "GrossmannLee" is used, epsilon + should be set very small, as the formulation is only correct in + the epsilon -> 0 limit. In particular, it should be small + enough to put spurious O(eps) sized constraint violations + within solver tolerances. Both "GrossmannLee" (when epsilon is + small enough) and the original "LeeGrossmann" have numerical + and feasibility issues. + + References ---------- .. [1] Lee, S., & Grossmann, I. E. (2000). New algorithms for @@ -178,6 +211,7 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): default=1e-4, domain=cfg.PositiveFloat, description="Epsilon value to use in perspective function", + doc="See the doc for 'perspective function' for discussion.", ), ) CONFIG.declare( @@ -197,6 +231,45 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): """, ), ) + CONFIG.declare( + 'well_defined_points', + cfg.ConfigValue( + default=ComponentMap(), + domain=ComponentMap, + description="Distinguished points at which constraints with restricted " + "domain are well-defined. This will be used as a center point for " + "transformed constraints.", + doc=""" + Dict-like mapping Disjunctions to ComponentMaps + mapping Vars appearing on the Disjuncts of that Disjunction to + float values, such that each constraint function appearing on + those disjuncts is well-defined (no division by zero, logarithm + of negative, etc.) when those vars are set to those values. The + outer map need not contain every Disjunction in the model as a + key, but the inner ComponentMaps (if any) should have as keys + every variable appearing on those disjuncts (including fixed + vars unless assume_fixed_vars_permanent is set to True). When + this is not provided for a disjunction, as it usually need not + be, we first try the point with all variables zero, then we + make a best effort to find a nonzero point through a subsolver + call, then we raise GDP_Error if neither attempt was successful. + """, + ), + ) + CONFIG.declare( + 'well_defined_points_heuristic_solver', + cfg.ConfigValue( + default='gurobi_direct_minlp', + domain=Solver, + description="Solver used for the base points heuristic", + doc=""" + Solver used for the base points heuristic. This must be a + nonconvex NLP solver in general. Pass as a solver object + supporting the V1 solver API version, or a corresponding string + for SolverFactory. + """, + ), + ) transformation_name = 'hull' def __init__(self): @@ -292,6 +365,114 @@ def _add_transformation_block(self, to_block): return transBlock, True + # From a test expression test_expr containing exactly the variables + # regular_vars and fallback_vars, except possibly containing some + # additional fixed vars, get a point at which test_expr is + # well-defined according to the following process: + # (1) try the origin + # (2) try fixing fallback_vars at zero, and allow a solver to + # change regular_vars + # (3) try allowing a solver to change all the vars + # + # If a point is found, return a ComponentMap x0_map from Var to the + # found numeric value, and a ComponentSet used_vars of all Vars that + # were given a nonzero value as part of that process. If no such + # point can be found, raise a GDP_Error. + # + # If we have assume_fixed_vars_permanent=True, then the fixed vars + # will not appear in regular_vars or fallback_vars, and we will not + # try to change them for the purpose of finding a well-defined + # point, nor will they appear in the x0_map. Otherwise, i.e., + # whenever they are actually passed in, they will be treated as + # normal variables. We will restore the original variable values and + # fixed statuses when we return successfully. + def _get_well_defined_point( + self, test_expr, regular_vars, fallback_vars, disj_name + ): + # First, see if test_expr is well-defined at the origin. + x0_map = ComponentMap() + orig_values = ComponentMap() + orig_fixed = ComponentMap() + for x in itertools.chain(regular_vars, fallback_vars): + x0_map[x] = 0 # ZeroConstant? + # also do some setup here + orig_values[x] = value(x, exception=False) + orig_fixed[x] = x.fixed + try: + val = value( + EXPR.ExpressionReplacementVisitor(substitute=x0_map).walk_expression( + test_expr + ) + ) + if math.isfinite(val): + return x0_map, ComponentSet() + except ValueError: # ('math domain error') + pass + except ZeroDivisionError: + pass + except Exception as e: # can anything else be thrown here? + logger.error( + "While trying to evaluate an expression, got unexpected exception type " + f"{e.__class__.__name__} (was prepared for success or a ValueError)." + ) + raise + # Second, try making it well-defined by editing only the regular vars + for x in fallback_vars: + x.fix(0) + for x in regular_vars: + x.set_value(0) + x.unfix() + test_model = ConcreteModel() + test_model.test_expr = Expression(expr=test_expr) + test_model.obj = Objective(expr=0) + # In case the solver can't deal with Vars it doesn't know about + for x in itertools.chain(regular_vars, fallback_vars): + test_model.add_component( + unique_component_name(test_model, x.name), Reference(x) + ) + test_model.well_defined_cons = ConstraintList() + _WellDefinedConstraintGenerator( + cons_list=test_model.well_defined_cons + ).walk_expression(test_expr) + feasible = self._solve_for_first_feasible_solution(test_model) + # Third, try again, but edit all the vars + if not feasible: + for x in fallback_vars: + x.unfix() + feasible = self._solve_for_first_feasible_solution(test_model) + if not feasible: + raise GDP_Error( + f"Unable to find a well-defined point on disjunction {disj_name}. " + "To carry out the hull transformation, each disjunction must have a " + "point at which every constraint function appearing in its " + "disjuncts is well-defined and finite. Please ensure such a point " + "actually exists, then if we still cannot find it, override our " + "search process using the `well_defined_points` option." + ) + # Found a point + x0_map = ComponentMap() + used_vars = ComponentSet() + for x in itertools.chain(regular_vars, fallback_vars): + x0_map[x] = value(x) + if x0_map[x] != 0: + used_vars.add(x) + x.set_value(orig_values[x]) + x.fixed = orig_fixed[x] + return x0_map, used_vars + + # Use gurobi_direct_minlp for the heuristic for now. It needs to be + # a nonlinear solver. + def _solve_for_first_feasible_solution(self, test_model): + results = self._config.well_defined_points_heuristic_solver.solve( + test_model, load_solutions=False + ) + if results.solver.termination_condition is TerminationCondition.infeasible: + return False + if results.solver.status is not SolverStatus.ok: + raise GDP_Error(f"Unexpected solver status {results.solver.status}.") + test_model.solutions.load_from(results) + return True + def _transform_disjunctionData( self, obj, index, parent_disjunct, local_vars_by_disjunct ): @@ -329,20 +510,20 @@ def _transform_disjunctionData( # which diaggregated variables we need. var_order = ComponentSet() disjuncts_var_appears_in = ComponentMap() + active_constraints = set() # For each disjunct in the disjunction, we will store a list of Vars # that need a disaggregated counterpart in that disjunct. disjunct_disaggregated_var_map = {} for disjunct in active_disjuncts: # create the key for each disjunct now disjunct_disaggregated_var_map[disjunct] = ComponentMap() - for var in get_vars_from_components( - disjunct, + for con in disjunct.component_data_objects( Constraint, - include_fixed=not self._config.assume_fixed_vars_permanent, active=True, sort=SortComponents.deterministic, descend_into=Block, ): + active_constraints.add(con) # [ESJ 02/14/2020] By default, we disaggregate fixed variables # on the philosophy that fixing is not a promise for the future # and we are mathematically wrong if we don't transform these @@ -354,11 +535,18 @@ def _transform_disjunctionData( # Note that, because ComponentSets are ordered, we will # eventually disaggregate the vars in a deterministic order # (the order that we found them) - if var not in var_order: - var_order.add(var) - disjuncts_var_appears_in[var] = ComponentSet([disjunct]) - else: - disjuncts_var_appears_in[var].add(disjunct) + seen_vars = set() + for var in IdentifyVariableVisitor( + include_fixed=not self._config.assume_fixed_vars_permanent + ).walk_expression(con.expr): + if id(var) in seen_vars: + continue + seen_vars.add(id(var)) + if var not in var_order: + var_order.add(var) + disjuncts_var_appears_in[var] = ComponentSet([disjunct]) + else: + disjuncts_var_appears_in[var].add(disjunct) # Now, we will disaggregate all variables that are not explicitly # declared as being local. If we are moving up in a nested tree, we have @@ -371,28 +559,56 @@ def _transform_disjunctionData( # do this, we will explicitly collect the set of local_vars in this # loop. local_vars = defaultdict(ComponentSet) + all_local_vars = ComponentSet() + # This set contains vars that potentially must be disaggregated, + # but do not need to worry about global constraints, and can + # safely be added to the local vars of any parent disjunct. This + # is a superset of all_local_vars, and the extra elements are + # members of all_vars_to_disaggregate. + generalized_local_vars = ComponentSet() for var in var_order: disjuncts = disjuncts_var_appears_in[var] - # clearly not local if used in more than one disjunct - if len(disjuncts) > 1: - if self._generate_debug_messages: - logger.debug( - "Assuming '%s' is not a local var since it is" - "used in multiple disjuncts." % var.name - ) + for disj in disjuncts: + if disj in local_vars_by_disjunct: + if var in local_vars_by_disjunct[disj]: + if len(disjuncts) == 1: + local_vars[disj].add(var) + all_local_vars.add(var) + else: + vars_to_disaggregate[disjunct].add(var) + all_vars_to_disaggregate.add(var) + generalized_local_vars.add(var) + if var not in generalized_local_vars: for disj in disjuncts: + # Not a local var, so we must disaggregate, even if + # it's only on one disjunct. vars_to_disaggregate[disj].add(var) - all_vars_to_disaggregate.add(var) - else: # var only appears in one disjunct - disjunct = next(iter(disjuncts)) - # We check if the user declared it as local - if disjunct in local_vars_by_disjunct: - if var in local_vars_by_disjunct[disjunct]: - local_vars[disjunct].add(var) - continue - # It's not declared local to this Disjunct, so we - # disaggregate - vars_to_disaggregate[disjunct].add(var) + all_vars_to_disaggregate.add(var) + + # Find a well-defined point x_0. We need every constraint body + # to successfully evaluate to something. + if obj in self._config.well_defined_points: + x0_map = self._config.well_defined_points[obj] + offset_vars = ComponentSet() + for x, val in x0_map.items(): + if val != 0: + offset_vars.add(x) + else: + x0_map, offset_vars = self._get_well_defined_point( + test_expr=sum(con.body for con in active_constraints), + regular_vars=all_vars_to_disaggregate, + fallback_vars=all_local_vars, + disj_name=obj.name, + ) + transBlock.parent_block().private_data().well_defined_points_map[obj] = x0_map + # Any var that got an offset cannot be local anymore, but it can + # still be generalized local + for var in offset_vars: + if var in all_local_vars: + var_disjunct = next(iter(disjuncts_var_appears_in[var])) + local_vars[var_disjunct].remove(var) + all_local_vars.remove(var) + vars_to_disaggregate[var_disjunct].add(var) all_vars_to_disaggregate.add(var) # Now that we know who we need to disaggregate, we will do it @@ -413,7 +629,13 @@ def _transform_disjunctionData( parent_local_var_suffix=parent_local_var_list, parent_disjunct_local_vars=local_vars_by_disjunct[parent_disjunct], disjunct_disaggregated_var_map=disjunct_disaggregated_var_map, + x0_map=x0_map, + offset_vars=offset_vars, ) + # The parent disjunct's local vars should include all of this + # disjunction's local or generalized local vars. + local_vars_by_disjunct[parent_disjunct].update(generalized_local_vars) + xorConstraint.add(index, (or_expr, 1)) # map the DisjunctionData to its XOR constraint to mark it as # transformed @@ -421,13 +643,24 @@ def _transform_disjunctionData( # Now add the reaggregation constraints for var in all_vars_to_disaggregate: - # There are two cases here: Either the var appeared in every - # disjunct in the disjunction, or it didn't. If it did, there's - # nothing special to do: All of the disaggregated variables have - # been created, and we can just proceed and make this constraint. If - # it didn't, we need one more disaggregated variable, correctly - # defined. And then we can make the constraint. - if len(disjuncts_var_appears_in[var]) < len(active_disjuncts): + # If a var did not appear in every disjunct of the + # disjunction, then we (intentionally) did not create a + # complete set of disaggregated vars and corresponding + # bounds constraints for it. This would cause the variable + # to be forced to zero when no disjunct containing it is + # selected. If the var were local, this would not matter, + # but unless we were able to put it in + # generalized_local_vars earlier, it is possible that it + # could appear in other parts of the model. It is therefore + # necessary that it be unconstrained when no disjunct + # containing it is selected. We implement this by adding one + # more disaggregated variable which becomes active if none + # of the disjuncts containing the original var were + # selected. Its only constraints are the bounds constraints. + if ( + len(disjuncts_var_appears_in[var]) < len(active_disjuncts) + and var not in generalized_local_vars + ): # create one more disaggregated var idx = len(disaggregatedVars) disaggregated_var = disaggregatedVars[idx] @@ -446,6 +679,7 @@ def _transform_disjunctionData( disjunct=obj, bigmConstraint=disaggregated_var_bounds, var_free_indicator=var_free, + x0_map=x0_map, var_idx=idx, ) original_var_info = var.parent_block().private_data() @@ -476,7 +710,9 @@ def _transform_disjunctionData( # constraint will be transformed again. (And if it turns out # everything in it is local, then that transformation won't actually # change the mathematical expression, so it's okay. - disaggregationConstraint.add(cons_idx, var == disaggregatedExpr) + disaggregationConstraint.add( + cons_idx, var - x0_map[var] == disaggregatedExpr + ) # and update the map so that we can find this later. We index by # variable and the particular disjunction because there is a # different one for each disjunction @@ -494,6 +730,8 @@ def _transform_disjunct( parent_local_var_suffix, parent_disjunct_local_vars, disjunct_disaggregated_var_map, + x0_map, + offset_vars, ): relaxationBlock = self._get_disjunct_transformation_block(obj, transBlock) @@ -535,6 +773,7 @@ def _transform_disjunct( disjunct=obj, bigmConstraint=bigmConstraint, var_free_indicator=obj.indicator_var.get_associated_binary(), + x0_map=x0_map, ) # update the bigm constraint mappings data_dict = disaggregatedVar.parent_block().private_data() @@ -554,14 +793,13 @@ def _transform_disjunct( bigmConstraint = Constraint(transBlock.lbub) relaxationBlock.add_component(conName, bigmConstraint) - parent_block = var.parent_block() - self._declare_disaggregated_var_bounds( original_var=var, disaggregatedVar=var, disjunct=obj, bigmConstraint=bigmConstraint, var_free_indicator=obj.indicator_var.get_associated_binary(), + x0_map=x0_map, # trivial in this case ) # update the bigm constraint mappings data_dict = var.parent_block().private_data() @@ -571,19 +809,17 @@ def _transform_disjunct( var_substitute_map = dict( (id(v), newV) for v, newV in disjunct_disaggregated_var_map[obj].items() ) - zero_substitute_map = dict( - (id(v), ZeroConstant) + x0_substitute_map = dict( + (id(v), x0_map[v]) for v, newV in disjunct_disaggregated_var_map[obj].items() ) - # Transform each component within this disjunct + # Transform each component within this disjunct. In particular, + # call _transform_constraint on each constraint. self._transform_block_components( - obj, obj, var_substitute_map, zero_substitute_map + obj, obj, var_substitute_map, x0_substitute_map ) - # Anything that was local to this Disjunct is also local to the parent, - # and just got "promoted" up there, so to speak. - parent_disjunct_local_vars.update(local_vars) # deactivate disjunct so writers can be happy obj._deactivate_without_fixing_indicator() @@ -594,6 +830,7 @@ def _declare_disaggregated_var_bounds( disjunct, bigmConstraint, var_free_indicator, + x0_map, var_idx=None, ): # For updating mappings: @@ -603,14 +840,14 @@ def _declare_disaggregated_var_bounds( disaggregated_var_info.bigm_constraint_map[disaggregatedVar][disjunct] = {} - lb = original_var.lb - ub = original_var.ub - if lb is None or ub is None: + if original_var.lb is None or original_var.ub is None: raise GDP_Error( "Variables that appear in disjuncts must be " "bounded in order to use the hull " "transformation! Missing bound for %s." % (original_var.name) ) + lb = original_var.lb - x0_map[original_var] + ub = original_var.ub - x0_map[original_var] disaggregatedVar.setlb(min(0, lb)) disaggregatedVar.setub(max(0, ub)) @@ -656,7 +893,7 @@ def _get_local_var_list(self, parent_disjunct): return local_var_list def _transform_constraint( - self, obj, disjunct, var_substitute_map, zero_substitute_map + self, obj, disjunct, var_substitute_map, x0_substitute_map ): # we will put a new transformed constraint on the relaxation block. relaxationBlock = disjunct._transformation_block() @@ -676,7 +913,23 @@ def _transform_constraint( unique = len(newConstraint) name = c.local_name + "_%s" % unique - NL = c.body.polynomial_degree() not in (0, 1) + # HACK: If I'm not treating fixed vars as permanent, I need + # the LinearRepnVisitor to pick them up, but it really + # doesn't want to. So I'm going to unfix and refix, which + # involves walking the expression twice, plus it's generally + # stupid + to_refix = ComponentSet() + if not self._config.assume_fixed_vars_permanent: + for var in EXPR.identify_variables(c.body, include_fixed=True): + if var.fixed: + to_refix.add(var) + var.unfix() + linear_repn = LinearRepnVisitor( + {}, var_recorder=VarRecorder({}, SortComponents.deterministic) + ).walk_expression(c.body) + for var in to_refix: + var.fix() + NL = linear_repn.nonlinear is not None EPS = self._config.EPS mode = self._config.perspective_function @@ -684,8 +937,12 @@ def _transform_constraint( # we substitute the expression variables with the # disaggregated variables if not NL or mode == "FurmanSawayaGrossmann": - h_0 = clone_without_expression_components( - c.body, substitute=zero_substitute_map + # Only permanently fixed vars are not being substituted + # by the x0_substitute_map + h_x0 = value( + clone_without_expression_components( + c.body, substitute=x0_substitute_map + ) ) y = disjunct.binary_indicator_var @@ -694,7 +951,8 @@ def _transform_constraint( sub_expr = clone_without_expression_components( c.body, substitute=dict( - (var, subs / y) for var, subs in var_substitute_map.items() + (var, (subs / y) + x0_substitute_map[var]) + for var, subs in var_substitute_map.items() ), ) expr = sub_expr * y @@ -702,7 +960,7 @@ def _transform_constraint( sub_expr = clone_without_expression_components( c.body, substitute=dict( - (var, subs / (y + EPS)) + (var, (subs / (y + EPS)) + x0_substitute_map[var]) for var, subs in var_substitute_map.items() ), ) @@ -711,42 +969,77 @@ def _transform_constraint( sub_expr = clone_without_expression_components( c.body, substitute=dict( - (var, subs / ((1 - EPS) * y + EPS)) + ( + var, + (subs / ((1 - EPS) * y + EPS)) + x0_substitute_map[var], + ) for var, subs in var_substitute_map.items() ), ) - expr = ((1 - EPS) * y + EPS) * sub_expr - EPS * h_0 * (1 - y) + expr = ((1 - EPS) * y + EPS) * sub_expr - EPS * h_x0 * (1 - y) else: raise RuntimeError("Unknown NL Hull mode") else: - expr = clone_without_expression_components( - c.body, substitute=var_substitute_map + # For a linear constraint that looks like a^Tx + b <= c, + # the transformed constraint will be a^Tv <= lambda (c - + # b - a^Tx_0). Get the a^Tv here and note that b + + # a^Tx_0 is exactly h_x0 from earlier, so we will have + # it when we need it. + # + # Note: linear_repn.multiplier is always 1 when obtained + # from LinearRepnVisitor.walk_expression so we do not + # need to read it + expr = sum( + coef * var_substitute_map[var] + for var, coef in linear_repn.linear.items() + if coef != 0 ) if c.equality: if NL: - # ESJ TODO: This can't happen right? This is the only - # obvious case where someone has messed up, but this has to - # be nonconvex, right? Shouldn't we tell them? + # NOTE: This nonlinear equality constraint is + # probably nonconvex, depending on your definition + # of nonconvex (it is never in the standard form for + # a convex problem unless the constraint body is a + # disguised affine function, but the feasible region + # may still be convex under weaker conditions, + # e.g. quasilinearity). But even so, this + # reformulation is still correct, so we will not + # complain to the user. newConsExpr = expr == c.lower * y else: - v = list(EXPR.identify_variables(expr)) - if len(v) == 1 and not c.lower: - # Setting a variable to 0 in a disjunct is - # *very* common. We should recognize that in - # that structure, the disaggregated variable - # will also be fixed to 0. - v[0].fix(0) - # ESJ: If you ask where the transformed constraint is, - # the answer is nowhere. Really, it is in the bounds of - # this variable, so I'm going to return - # it. Alternatively we could return an empty list, but I - # think I like this better. - constraint_map.transformed_constraints[c].append(v[0]) - # Reverse map also (this is strange) - constraint_map.src_constraint[v[0]] = c - continue - newConsExpr = expr - (1 - y) * h_0 == c.lower * y + if len(linear_repn.linear) == 1: + var, coef = next(iter(linear_repn.linear.items())) + # Second clause of this condition happens iff + # the constraint implies x = x0 + if ( + coef != 0 + and (c.lower - linear_repn.constant) / coef + == x0_substitute_map[var] + ): + v = var_substitute_map[var] + # Setting a variable to 0 in a disjunct is + # *very* common. We should recognize that in + # that structure, the disaggregated variable + # will also be fixed to 0. In the + # general-offset case, this happens when the + # equality constraint is of the form x = + # x_0. We're unlikely to hit this unless x_0 + # is the origin or was passed manually, but + # that's fine - nonzero x_0 is already a + # rare special case. + v.fix(0) + # ESJ: If you ask where the transformed + # constraint is, the answer is + # nowhere. Really, it is in the bounds of + # this variable, so I'm going to return + # it. Alternatively we could return an empty + # list, but I think I like this better. + constraint_map.transformed_constraints[c].append(v) + # Reverse map also (this is strange) + constraint_map.src_constraint[v] = c + continue + newConsExpr = expr == (c.lower - h_x0) * y if obj.is_indexed(): newConstraint.add((name, i, 'eq'), newConsExpr) @@ -778,7 +1071,7 @@ def _transform_constraint( if NL: newConsExpr = expr >= c.lower * y else: - newConsExpr = expr - (1 - y) * h_0 >= c.lower * y + newConsExpr = expr >= (c.lower - h_x0) * y if obj.is_indexed(): newConstraint.add((name, i, 'lb'), newConsExpr) @@ -800,7 +1093,7 @@ def _transform_constraint( if NL: newConsExpr = expr <= c.upper * y else: - newConsExpr = expr - (1 - y) * h_0 <= c.upper * y + newConsExpr = expr <= (c.upper - h_x0) * y if obj.is_indexed(): newConstraint.add((name, i, 'ub'), newConsExpr) @@ -913,7 +1206,7 @@ def get_disaggregation_constraint( .private_data() .disaggregation_constraint_map[original_var][disjunction] ) - except: + except Exception: if raise_exception: logger.error( "It doesn't appear that '%s' is a variable that was " @@ -972,6 +1265,18 @@ def get_transformed_constraints(self, cons): cons = transformed_cons return cons + def get_well_defined_points_map(self, b): + """ + Retrieve the well-defined points originally used to transform + a Block. Format is a ComponentMap of ComponentMaps identical to + that of the parameter well_defined_points. + + Parameters + ---------- + b: a Block that was transformed by gdp.hull + """ + return b.private_data().well_defined_points_map + @TransformationFactory.register( 'gdp.chull', @@ -986,3 +1291,106 @@ def get_transformed_constraints(self, cons): class _Deprecated_Name_Hull(Hull_Reformulation): def __init__(self): super(_Deprecated_Name_Hull, self).__init__() + + +# Walk the expression and, whenever encountering an expression that +# could fail to be well-defined, create a constraint that keeps it +# well-defined. This is done at the Pyomo level to make this capability +# more generic instead of needing to specially set up the options for +# each solver (plus, some solvers get rather buggy when used for this +# task). +class _WellDefinedConstraintGenerator(StreamBasedExpressionVisitor): + def __init__(self, cons_list, **kwds): + self.cons_list = cons_list + super().__init__(**kwds) + + # Whenever we exit a node (entering would also be fine) check if it + # has a restricted domain, and if it does add a corresponding + # constraint + def exitNode(self, node, data): + if node.__class__ in _handlers: + for con in _handlers[node.__class__](node): + # note: con should never be a boolean True here, such + # cases should have been filtered out during the handler + # call + self.cons_list.add(con) + + +# Epsilon for handling function domains with strict inequalities. This +# is a heuristic so it's not important for this to be tight. +EPS = 1e-4 + + +def _handlePowExpression(node): + base, exp = node.args + # if base is not variable, nothing for us to do + if base.__class__ in EXPR.native_types or not base.is_potentially_variable(): + return () + + # If exp is a NPV nonnegative integer, there are no restrictions on + # base. If exp is a NPV negative integer, base should not be + # zero. If exp is a NPV nonnegative fraction, base should not be + # negative. Otherwise, base should be strictly positive (as we can't + # be sure that exp could not be negative or fractional). + + # Note: this is problematic for LP, but I don't want to potentially + # invoke a MIP solve here, so replace "x is nonzero" with "x is >= + # eps". It's a heuristic so this is not critical. + if exp.__class__ in EXPR.native_types or not exp.is_potentially_variable(): + val = value(exp) + if round(val) == val: + if val >= 0: + return () + else: + # return base != 0 + return (base >= EPS,) + elif val >= 0: + return (base >= 0,) + return (base >= EPS,) + + +def _handleDivisionExpression(node): + # No division by zero. Dividing by an NPV is always allowed. + arg = node.args[1] + if arg.__class__ in EXPR.native_types or not arg.is_potentially_variable(): + return () + # Same LP vs MIP problem as before + return (arg >= EPS,) + + +def _handleUnaryFunctionExpression(node): + arg = node.args[0] + if arg.__class__ in EXPR.native_types or not arg.is_potentially_variable(): + return () + if node.name == 'log': + return (arg >= EPS,) + if node.name == 'log10': + return (arg >= EPS,) + if node.name == 'sqrt': + return (arg >= 0,) + if node.name == 'asin': + return (arg >= -1, arg <= 1) + if node.name == 'acos': + return (arg >= -1, arg <= 1) + # It can't be exactly pi/2 plus a multiple of pi. Rather difficult + # to enforce, so make a conservative effort by instead keeping it in + # (-pi/2, pi/2). + if node.name == 'tan': + return (arg >= -(math.pi / 2) + EPS, arg <= (math.pi / 2) - EPS) + if node.name == 'acosh': + return (arg >= 1,) + if node.name == 'atanh': + return (arg >= -1 + EPS, arg <= 1 - EPS) + # Other Pyomo unary functions should be well-defined. + return () + + +# All expression types that can potentially be +# ill-defined: +_handlers = { + # You are on your own here + # EXPR.ExternalFunctionExpression, + EXPR.PowExpression: _handlePowExpression, + EXPR.DivisionExpression: _handleDivisionExpression, + EXPR.UnaryFunctionExpression: _handleUnaryFunctionExpression, +} diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 8c59ef0bb55..cbaf40bf073 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -11,6 +11,7 @@ import sys import random from io import StringIO +import math import pyomo.common.unittest as unittest @@ -28,6 +29,14 @@ ComponentMap, value, log, + log10, + sqrt, + asin, + acos, + acosh, + atanh, + tan, + sin, ConcreteModel, Any, Suffix, @@ -36,6 +45,8 @@ Param, Objective, TerminationCondition, + SortComponents, + ConstraintList, ) from pyomo.core.expr.compare import ( assertExpressionsEqual, @@ -49,6 +60,7 @@ from pyomo.gdp import Disjunct, Disjunction, GDP_Error import pyomo.gdp.tests.models as models import pyomo.gdp.tests.common_tests as ct +from pyomo.gdp.plugins.hull import _WellDefinedConstraintGenerator currdir = this_file_dir() @@ -2253,6 +2265,7 @@ def test_infeasible_xor_because_all_disjuncts_deactivated(self): relaxed_xor = hull.get_transformed_constraints( m.disjunction_disjuncts[0].nestedDisjunction.algebraic_constraint ) + self.assertEqual(len(relaxed_xor), 1) relaxed_xor = relaxed_xor[0] repn = generate_standard_repn(relaxed_xor.body) @@ -2659,10 +2672,7 @@ def test_logical_constraints_transformed(self): c = cons[0] # hull transformation of z1 >= 1 assertExpressionsStructurallyEqual( - self, - c.expr, - dis_z1 - (1 - m.d[1].binary_indicator_var) * 0 - >= m.d[1].binary_indicator_var, + self, c.expr, dis_z1 >= m.d[1].binary_indicator_var ) # then d[4]: @@ -2832,10 +2842,7 @@ def test_logical_constraints_transformed(self): self.assertEqual(len(cons), 1) cons = cons[0] assertExpressionsStructurallyEqual( - self, - cons.expr, - 1 - z3d - (2 - (z1d + z2d)) - (1 - m.d[4].binary_indicator_var) * (-1) - <= 0 * m.d[4].binary_indicator_var, + self, cons.expr, -z3d + z1d + z2d <= m.d[4].binary_indicator_var ) # hull transformation of z3 >= 1 @@ -2845,9 +2852,7 @@ def test_logical_constraints_transformed(self): self.assertEqual(len(cons), 1) cons = cons[0] assertExpressionsStructurallyEqual( - self, - cons.expr, - z3d - (1 - m.d[4].binary_indicator_var) * 0 >= m.d[4].binary_indicator_var, + self, cons.expr, z3d >= m.d[4].binary_indicator_var ) self.assertFalse(m.bwahaha.active) @@ -2882,6 +2887,87 @@ def test_dill_pickle(self): sys.setrecursionlimit(rl) +class DomainRestrictionTest(unittest.TestCase): + def test_no_nonlinear(self): + m = models.makeTwoTermDisj() + TransformationFactory('gdp.hull').apply_to(m) + + m_pts = TransformationFactory('gdp.hull').get_well_defined_points_map(m) + # both points have zero offset + self.assertEqual(m_pts[m.disjunction][m.a], 0) + self.assertEqual(m_pts[m.disjunction][m.x], 0) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_simple_case(self): + m = models.makeTwoTermDisj() + m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) + m2 = m.clone() + TransformationFactory('gdp.hull').apply_to(m) + + m_pts = TransformationFactory('gdp.hull').get_well_defined_points_map(m) + # x gets a nonzero offset; a does not (but it might by chance) + self.assertNotEqual(m_pts[m.disjunction][m.x], 0) + # using the recorded map is the same thing as the process that + # created that map + TransformationFactory('gdp.hull').apply_to(m2, well_defined_points=m_pts) + # how to assert that two blocks have (value) identical contents? + for c1, c2 in zip( + m.component_data_objects(Constraint, sort=SortComponents.deterministic), + m2.component_data_objects(Constraint, sort=SortComponents.deterministic), + ): + self.assertEqual(str(c1.expr), str(c2.expr)) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_handle_fixed_disagg(self): + m = models.makeTwoTermDisj() + m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) + m.x.fix(2) + hull = TransformationFactory('gdp.hull') + hull.apply_to(m, assume_fixed_vars_permanent=False) + self.assertIn(m.x, hull.get_well_defined_points_map(m)[m.disjunction]) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_handle_fixed_no_disagg(self): + m = models.makeTwoTermDisj() + m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) + m.x.fix(2) + hull = TransformationFactory('gdp.hull') + hull.apply_to(m, assume_fixed_vars_permanent=True) + self.assertNotIn(m.x, hull.get_well_defined_points_map(m)[m.disjunction]) + # transformed nonlinear constraint is trivial here + assertExpressionsEqual( + self, + m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + .transformedConstraints['nonlinear_1', 'lb'] + .body, + 0.0 * m.d[0].binary_indicator_var, + ) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_no_good_point(self): + m = models.makeTwoTermDisj() + m.d[0].nonlinear = Constraint(expr=log(-((m.x) ** 2) - 1) >= 0) + self.assertRaisesRegex( + GDP_Error, + "Unable to find a well-defined point on disjunction .*", + TransformationFactory('gdp.hull').apply_to, + m, + ) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_various_nonlinear(self): + m = models.makeTwoTermDisj() + m.y = Var(bounds=(-5, 5)) + m.d[0].log = Constraint(expr=log(m.a + 1) >= 0) + m.d[0].pow = Constraint(expr=m.y ** (-0.5) >= 0) + m.d[0].div = Constraint(expr=1 / (m.x) >= 0) + hull = TransformationFactory('gdp.hull') + hull.apply_to(m) + m_pts = hull.get_well_defined_points_map(m) + self.assertNotEqual(m_pts[m.disjunction][m.x], 0) + self.assertNotEqual(m_pts[m.disjunction][m.y], 0) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") class NestedDisjunctsInFlatGDP(unittest.TestCase): """ @@ -2890,3 +2976,45 @@ class NestedDisjunctsInFlatGDP(unittest.TestCase): def test_declare_disjuncts_in_disjunction_rule(self): ct.check_nested_disjuncts_in_flat_gdp(self, 'hull') + + +class WellDefinedConstraintWalkerTest(unittest.TestCase): + def test_many_constraints(self): + m = models.makeTwoTermDisj() + m.hard = Constraint(expr=1 / (log(atanh(m.x))) >= 0) + m.several = Constraint( + expr=log(m.x + 1) + + log10(m.x + 2) + + sqrt(m.x + 3) + + asin(m.x + 4) + + acos(m.x + 5) + + tan(m.x + 6) + + acosh(m.x + 7) + + atanh(m.x + 8) + >= 0 + ) + m.cons = ConstraintList() + walker = _WellDefinedConstraintGenerator(cons_list=m.cons) + walker.walk_expression(m.hard.body) + walker.walk_expression(m.several.body) + walked_strings = [str(con.expr) for con in m.cons.values()] + walker_eps = 1e-4 + for con in [ + log(atanh(m.x)) >= walker_eps, + atanh(m.x) >= walker_eps, + m.x >= -1 + walker_eps, + m.x <= 1 - walker_eps, + m.x + 1 >= walker_eps, + m.x + 2 >= walker_eps, + m.x + 3 >= 0, + m.x + 4 >= -1, + m.x + 4 <= 1, + m.x + 5 >= -1, + m.x + 5 <= 1, + m.x + 6 <= (math.pi / 2) - walker_eps, + m.x + 6 >= -(math.pi / 2) + walker_eps, + m.x + 7 >= 1, + m.x + 8 <= 1 - walker_eps, + m.x + 8 >= -1 + walker_eps, + ]: + self.assertIn(str(con), walked_strings)