From 3d98cd5c28a7a0dbd50cd209554e8248a7ffb0e5 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 22 Jul 2025 08:40:42 -0600 Subject: [PATCH 01/20] hull fix: initial --- pyomo/gdp/plugins/hull.py | 438 ++++++++++++++++++++++++++++++++------ 1 file changed, 371 insertions(+), 67 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 53ff1080940..5820ddee849 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -18,6 +18,7 @@ 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 @@ -40,6 +41,9 @@ value, NonNegativeIntegers, Binary, + ConcreteModel, + Objective, + Reference, ) from pyomo.gdp import Disjunct, Disjunction, GDP_Error from pyomo.gdp.disjunct import DisjunctData @@ -51,8 +55,14 @@ _warn_for_active_disjunct, ) from pyomo.core.util import target_list +from pyomo.core.expr.visitor import IdentifyVariableVisitor 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') @@ -131,7 +141,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 @@ -157,6 +167,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 serious + numerical issues. + + References ---------- .. [1] Lee, S., & Grossmann, I. E. (2000). New algorithms for @@ -180,6 +204,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( @@ -199,6 +224,32 @@ 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 finally we raise GDP_Error if neither attempt was + successful. + """, + ), + ) transformation_name = 'hull' def __init__(self): @@ -294,6 +345,111 @@ 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, 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. + # + # For the purposes of this function, fixed vars are treated as + # legitimate variables and may be changed to find the + # point. However, this function will restore the original variable + # values and fixed statuses when it returns successfully. + # + # TODO: treat vars like params if they're fixed and we have + # assume_fixed_vars_permanent set? Or not? + def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): + # 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: + # TODO: value() is only needed here because there can be + # fixed variables with values that did not appear in + # regular_vars or fallback_vars, so they do not appear in + # x0_map. This is probably an error. See test + # test_do_not_disaggregate_fixed_variables + # + # TODO actually, maybe that solves my problem for me? If I + # don't get the fixed variables when we have + # assume_fixed_vars_permanent=True, then I'll never put them + # in x0_map and therefore never need to offset + # them... Investigate this. + val = value( + EXPR.ExpressionReplacementVisitor(substitute=x0_map).walk_expression( + test_expr + ) + ) + if math.isfinite(val): + return x0_map, ComponentSet() + except Exception as e: # what types can we throw here? + breakpoint() + pass + # 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.unfix() + test_model = ConcreteModel() + test_model.obj = Objective(expr=test_expr) + # 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) + ) + 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( + # TODO finish error message + "Unable to find a well-defined point on disjunction {disjunction}. " + "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 {TODO_PARAM_NAME} option." + ) + # We have a point. + if not math.isfinite(value(test_expr)): + raise DeveloperError("Theoretically unreachable") + 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 + + # This one is going to get a little janky... + # Let's start with just baron for the time being. + def _solve_for_first_feasible_solution(self, test_model): + results = SolverFactory("baron").solve( + test_model, options={'NumSol': 1, 'FirstFeas': 1} + ) + 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}.") + return True + def _transform_disjunctionData( self, obj, index, parent_disjunct, local_vars_by_disjunct ): @@ -331,20 +487,28 @@ 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 var in get_vars_from_components( + # disjunct, + # Constraint, + # include_fixed=not self._config.assume_fixed_vars_permanent, + # active=True, + # sort=SortComponents.deterministic, + # descend_into=Block, + # ): + 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 @@ -356,11 +520,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 @@ -373,28 +544,75 @@ 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: - 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]: + # # 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: + # 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) + # all_local_vars.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) + for disj in disjuncts: + if disj in local_vars_by_disjunct[disjunct]: + if len(disjuncts) == 1: + # was this a noop before? local_vars[disjunct].add(var) - continue - # It's not declared local to this Disjunct, so we - # disaggregate - vars_to_disaggregate[disjunct].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) + else: + # Not a local var, so we must disaggregate, even if + # it's only on one disjunct. + 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, + ) + # 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 @@ -415,6 +633,8 @@ 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, ) xorConstraint.add(index, (or_expr, 1)) # map the DisjunctionData to its XOR constraint to mark it as @@ -429,7 +649,25 @@ def _transform_disjunctionData( # 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): + # TODO new comment: + # 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] @@ -448,6 +686,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() @@ -471,6 +710,7 @@ def _transform_disjunctionData( else: disaggregatedExpr = 0 for disjunct in disjuncts_var_appears_in[var]: + breakpoint() disaggregatedExpr += disjunct_disaggregated_var_map[disjunct][var] cons_idx = len(disaggregationConstraint) @@ -478,7 +718,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 @@ -496,6 +738,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) @@ -537,6 +781,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() @@ -556,7 +801,8 @@ def _transform_disjunct( bigmConstraint = Constraint(transBlock.lbub) relaxationBlock.add_component(conName, bigmConstraint) - parent_block = var.parent_block() + # TODO unused code? + # parent_block = var.parent_block() self._declare_disaggregated_var_bounds( original_var=var, @@ -564,6 +810,7 @@ def _transform_disjunct( 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() @@ -573,18 +820,22 @@ 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. + # + # TODO: this needs to add generalized local vars too. Why is + # this in transform_disjunct, can't I just do this centrally??? parent_disjunct_local_vars.update(local_vars) # deactivate disjunct so writers can be happy obj._deactivate_without_fixing_indicator() @@ -596,6 +847,7 @@ def _declare_disaggregated_var_bounds( disjunct, bigmConstraint, var_free_indicator, + x0_map, var_idx=None, ): # For updating mappings: @@ -605,8 +857,8 @@ def _declare_disaggregated_var_bounds( disaggregated_var_info.bigm_constraint_map[disaggregatedVar][disjunct] = {} - lb = original_var.lb - ub = original_var.ub + lb = original_var.lb - x0_map[original_var] + ub = original_var.ub - x0_map[original_var] if lb is None or ub is None: raise GDP_Error( "Variables that appear in disjuncts must be " @@ -658,7 +910,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() @@ -678,6 +930,7 @@ def _transform_constraint( unique = len(newConstraint) name = c.local_name + "_%s" % unique + # TODO: should this be a walker call? NL = c.body.polynomial_degree() not in (0, 1) EPS = self._config.EPS mode = self._config.perspective_function @@ -686,8 +939,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 + h_x0 = clone_without_expression_components( + # TODO: is this properly handling fixed vars when we + # consider them permanent? Do we need a value() + # call? + c.body, + substitute=x0_substitute_map, ) y = disjunct.binary_indicator_var @@ -696,7 +953,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 @@ -704,7 +962,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() ), ) @@ -713,42 +971,88 @@ 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 + c.body, + substitute=dict( + # It is possible to do this substituting v + # instead of v + x_0 here, but I think it would + # require picking apart the affine function into + # linear part and offset later... really that's + # just evaluating it at zero, which should be no + # problem... TODO consider this further. + (var, subs + x0_substitute_map[var]) + for var, subs in var_substitute_map.items() + ), ) 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 that we didn't catch, + # but the feasible region may still be convex under + # weaker conditions, e.g. quasilinearity). But even + # so, this reformulation is still correct (though + # some of its purpose evaporates), 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(v) == 1: + # TODO FIXME FIXME FIXME: this condition is + # wrong; it triggers in a case like 2x=x0, but + # that is not correct for x0 nonzero. + # + # TODO ask Emma about constraint + # canonicalization. The old version of this + # condition was clearly operating under the + # assumption that c.body has no affine part (and + # so is this edited version), but other parts of + # the code (eg setting up newConsExpr) appear + # not to be making that assumption, otherwise + # they would be simpler. + # + # TODO: also this is wrong for another reason: + # we already substituted so the id is not the + # same as before. I should be identifying + # variables on c.body or something. + # + # if c.lower == x0_substitute_map[id(v[0])]: + if False: + # 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[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_x0 == c.lower * y if obj.is_indexed(): newConstraint.add((name, i, 'eq'), newConsExpr) @@ -780,7 +1084,7 @@ def _transform_constraint( if NL: newConsExpr = expr >= c.lower * y else: - newConsExpr = expr - (1 - y) * h_0 >= c.lower * y + newConsExpr = expr - (1 - y) * h_x0 >= c.lower * y if obj.is_indexed(): newConstraint.add((name, i, 'lb'), newConsExpr) @@ -802,7 +1106,7 @@ def _transform_constraint( if NL: newConsExpr = expr <= c.upper * y else: - newConsExpr = expr - (1 - y) * h_0 <= c.upper * y + newConsExpr = expr - (1 - y) * h_x0 <= c.upper * y if obj.is_indexed(): newConstraint.add((name, i, 'ub'), newConsExpr) @@ -915,7 +1219,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 " From 2ea51345c1219d7abc0e13bf46209b419a186309 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 30 Jul 2025 09:47:30 -0600 Subject: [PATCH 02/20] tests passing --- pyomo/gdp/plugins/hull.py | 54 ++++++++++++++++++++++-------------- pyomo/gdp/tests/test_hull.py | 12 ++++++++ 2 files changed, 45 insertions(+), 21 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 5820ddee849..7d1d19fdacb 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -21,7 +21,7 @@ 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, @@ -56,6 +56,8 @@ ) from pyomo.core.util import target_list from pyomo.core.expr.visitor import IdentifyVariableVisitor +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 @@ -394,9 +396,10 @@ def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): ) if math.isfinite(val): return x0_map, ComponentSet() + except ValueError: # ('math domain error') + pass except Exception as e: # what types can we throw here? breakpoint() - pass # Second, try making it well-defined by editing only the regular vars for x in fallback_vars: x.fix(0) @@ -410,6 +413,7 @@ def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): unique_component_name(test_model, x.name), Reference(x) ) feasible = self._solve_for_first_feasible_solution(test_model) + # breakpoint() # Third, try again, but edit all the vars if not feasible: for x in fallback_vars: @@ -576,20 +580,22 @@ def _transform_disjunctionData( # vars_to_disaggregate[disjunct].add(var) # all_vars_to_disaggregate.add(var) for disj in disjuncts: - if disj in local_vars_by_disjunct[disjunct]: - if len(disjuncts) == 1: - # was this a noop before? - local_vars[disjunct].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) - else: + if disj in local_vars_by_disjunct: + if var in local_vars_by_disjunct[disj]: + if len(disjuncts) == 1: + # was this a noop before? + 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[disjunct].add(var) - all_vars_to_disaggregate.add(var) + vars_to_disaggregate[disj].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. @@ -607,6 +613,7 @@ def _transform_disjunctionData( ) # Any var that got an offset cannot be local anymore, but it can # still be generalized local + # breakpoint() for var in offset_vars: if var in all_local_vars: var_disjunct = next(iter(disjuncts_var_appears_in[var])) @@ -710,7 +717,7 @@ def _transform_disjunctionData( else: disaggregatedExpr = 0 for disjunct in disjuncts_var_appears_in[var]: - breakpoint() + # breakpoint() disaggregatedExpr += disjunct_disaggregated_var_map[disjunct][var] cons_idx = len(disaggregationConstraint) @@ -741,6 +748,7 @@ def _transform_disjunct( x0_map, offset_vars, ): + # breakpoint() relaxationBlock = self._get_disjunct_transformation_block(obj, transBlock) # Put the disaggregated variables all on their own block so that we can @@ -857,14 +865,14 @@ def _declare_disaggregated_var_bounds( disaggregated_var_info.bigm_constraint_map[disaggregatedVar][disjunct] = {} - lb = original_var.lb - x0_map[original_var] - ub = original_var.ub - x0_map[original_var] - 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)) @@ -930,8 +938,12 @@ def _transform_constraint( unique = len(newConstraint) name = c.local_name + "_%s" % unique - # TODO: should this be a walker call? - NL = c.body.polynomial_degree() not in (0, 1) + # vid -> var + repn_var_map = {} + linear_repn = LinearRepnVisitor( + {}, var_recorder=VarRecorder(repn_var_map, SortComponents.deterministic) + ).walk_expression(c.body) + NL = linear_repn.nonlinear is not None EPS = self._config.EPS mode = self._config.perspective_function @@ -1031,7 +1043,7 @@ def _transform_constraint( # variables on c.body or something. # # if c.lower == x0_substitute_map[id(v[0])]: - if False: + if c.lower - linear_repn.constant * linear_repn.multiplier == 0: # Setting a variable to 0 in a disjunct is # *very* common. We should recognize that in # that structure, the disaggregated variable diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 204bcfea4c0..e462c3b7fb9 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2884,6 +2884,18 @@ def test_dill_pickle(self): finally: sys.setrecursionlimit(rl) + +class DomainRestrictionTest(unittest.TestCase): + def test_simple_case(self): + m = models.makeTwoTermDisj() + m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) + TransformationFactory('gdp.hull').apply_to(m) + m.pprint() + # did not throw + # + # TODO: keep on the model somewhere what x0 was, and assert it + # has a nonzero element + @unittest.skipUnless(gurobi_available, "Gurobi is not available") class NestedDisjunctsInFlatGDP(unittest.TestCase): From 84b582ebae779178c4f74660a74dfc00436ad092 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 30 Jul 2025 14:10:44 -0600 Subject: [PATCH 03/20] hull tests: remove 0 * x terms from expected expressions --- pyomo/gdp/tests/test_hull.py | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index e462c3b7fb9..5dd79e83640 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2256,6 +2256,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) @@ -2662,10 +2663,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]: @@ -2835,10 +2833,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 @@ -2848,9 +2843,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) @@ -2884,7 +2877,7 @@ def test_dill_pickle(self): finally: sys.setrecursionlimit(rl) - + class DomainRestrictionTest(unittest.TestCase): def test_simple_case(self): m = models.makeTwoTermDisj() @@ -2895,7 +2888,7 @@ def test_simple_case(self): # # TODO: keep on the model somewhere what x0 was, and assert it # has a nonzero element - + @unittest.skipUnless(gurobi_available, "Gurobi is not available") class NestedDisjunctsInFlatGDP(unittest.TestCase): From 656de7e818dcd0c3ce21da73a5999e55119fcdf3 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 30 Jul 2025 14:12:10 -0600 Subject: [PATCH 04/20] simplify linear constraint handling and fix a bug --- pyomo/gdp/plugins/hull.py | 94 +++++++++++++++++++-------------------- 1 file changed, 47 insertions(+), 47 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 7d1d19fdacb..302983bf099 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -938,11 +938,22 @@ def _transform_constraint( unique = len(newConstraint) name = c.local_name + "_%s" % unique - # vid -> var - repn_var_map = {} + # HACK: If I'm not treating fixed vars as permanent, I need + # the LinearRepnVisitor to pick them up, but it very + # strongly 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(repn_var_map, SortComponents.deterministic) + {}, 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 @@ -994,18 +1005,15 @@ def _transform_constraint( else: raise RuntimeError("Unknown NL Hull mode") else: - expr = clone_without_expression_components( - c.body, - substitute=dict( - # It is possible to do this substituting v - # instead of v + x_0 here, but I think it would - # require picking apart the affine function into - # linear part and offset later... really that's - # just evaluating it at zero, which should be no - # problem... TODO consider this further. - (var, subs + x0_substitute_map[var]) - for var, subs in var_substitute_map.items() - ), + # 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. + expr = linear_repn.multiplier * sum( + coef * var_substitute_map[var] + for var, coef in linear_repn.linear.items() + if coef != 0 ) if c.equality: @@ -1022,28 +1030,19 @@ def _transform_constraint( # complain to the user. newConsExpr = expr == c.lower * y else: - v = list(EXPR.identify_variables(expr)) - if len(v) == 1: - # TODO FIXME FIXME FIXME: this condition is - # wrong; it triggers in a case like 2x=x0, but - # that is not correct for x0 nonzero. - # - # TODO ask Emma about constraint - # canonicalization. The old version of this - # condition was clearly operating under the - # assumption that c.body has no affine part (and - # so is this edited version), but other parts of - # the code (eg setting up newConsExpr) appear - # not to be making that assumption, otherwise - # they would be simpler. - # - # TODO: also this is wrong for another reason: - # we already substituted so the id is not the - # same as before. I should be identifying - # variables on c.body or something. - # - # if c.lower == x0_substitute_map[id(v[0])]: - if c.lower - linear_repn.constant * linear_repn.multiplier == 0: + 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 * linear_repn.multiplier + ) + / (coef * linear_repn.multiplier) + == 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 @@ -1054,17 +1053,18 @@ def _transform_constraint( # is the origin or was passed manually, but # that's fine - nonzero x_0 is already a # rare special case. - v[0].fix(0) - # ESJ: If you ask where the transformed constraint is, - # the answer is nowhere. Really, it is in the bounds of + 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[0]) + # 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[0]] = c + constraint_map.src_constraint[v] = c continue - newConsExpr = expr - (1 - y) * h_x0 == c.lower * y + newConsExpr = expr == (c.lower - h_x0) * y if obj.is_indexed(): newConstraint.add((name, i, 'eq'), newConsExpr) @@ -1096,7 +1096,7 @@ def _transform_constraint( if NL: newConsExpr = expr >= c.lower * y else: - newConsExpr = expr - (1 - y) * h_x0 >= c.lower * y + newConsExpr = expr >= (c.lower - h_x0) * y if obj.is_indexed(): newConstraint.add((name, i, 'lb'), newConsExpr) @@ -1118,7 +1118,7 @@ def _transform_constraint( if NL: newConsExpr = expr <= c.upper * y else: - newConsExpr = expr - (1 - y) * h_x0 <= c.upper * y + newConsExpr = expr <= (c.upper - h_x0) * y if obj.is_indexed(): newConstraint.add((name, i, 'ub'), newConsExpr) From 480ccdf27770faaf0bf9f034d595a8572c6247e7 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Thu, 31 Jul 2025 08:30:20 -0600 Subject: [PATCH 05/20] some comments and cleanup --- pyomo/gdp/plugins/hull.py | 94 ++++++++++++++---------------------- pyomo/gdp/tests/test_hull.py | 36 +++++++++++++- 2 files changed, 70 insertions(+), 60 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 302983bf099..4fb1aeae023 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -348,7 +348,8 @@ 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, get a point at which test_expr is + # 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 @@ -360,14 +361,16 @@ def _add_transformation_block(self, to_block): # were given a nonzero value as part of that process. If no such # point can be found, raise a GDP_Error. # - # For the purposes of this function, fixed vars are treated as - # legitimate variables and may be changed to find the - # point. However, this function will restore the original variable - # values and fixed statuses when it returns successfully. - # - # TODO: treat vars like params if they're fixed and we have - # assume_fixed_vars_permanent set? Or not? - def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): + # 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() @@ -378,17 +381,6 @@ def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): orig_values[x] = value(x, exception=False) orig_fixed[x] = x.fixed try: - # TODO: value() is only needed here because there can be - # fixed variables with values that did not appear in - # regular_vars or fallback_vars, so they do not appear in - # x0_map. This is probably an error. See test - # test_do_not_disaggregate_fixed_variables - # - # TODO actually, maybe that solves my problem for me? If I - # don't get the fixed variables when we have - # assume_fixed_vars_permanent=True, then I'll never put them - # in x0_map and therefore never need to offset - # them... Investigate this. val = value( EXPR.ExpressionReplacementVisitor(substitute=x0_map).walk_expression( test_expr @@ -398,8 +390,12 @@ def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): return x0_map, ComponentSet() except ValueError: # ('math domain error') pass - except Exception as e: # what types can we throw here? - breakpoint() + 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) @@ -413,7 +409,6 @@ def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): unique_component_name(test_model, x.name), Reference(x) ) feasible = self._solve_for_first_feasible_solution(test_model) - # breakpoint() # Third, try again, but edit all the vars if not feasible: for x in fallback_vars: @@ -421,13 +416,12 @@ def _get_well_defined_point(self, test_expr, regular_vars, fallback_vars): feasible = self._solve_for_first_feasible_solution(test_model) if not feasible: raise GDP_Error( - # TODO finish error message - "Unable to find a well-defined point on disjunction {disjunction}. " + 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 {TODO_PARAM_NAME} option." + "search process using the `well_defined_points` option." ) # We have a point. if not math.isfinite(value(test_expr)): @@ -498,14 +492,6 @@ def _transform_disjunctionData( 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, - # Constraint, - # include_fixed=not self._config.assume_fixed_vars_permanent, - # active=True, - # sort=SortComponents.deterministic, - # descend_into=Block, - # ): for con in disjunct.component_data_objects( Constraint, active=True, @@ -610,6 +596,7 @@ def _transform_disjunctionData( 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, ) # Any var that got an offset cannot be local anymore, but it can # still be generalized local @@ -650,13 +637,6 @@ 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. - # TODO new comment: # 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 @@ -809,9 +789,6 @@ def _transform_disjunct( bigmConstraint = Constraint(transBlock.lbub) relaxationBlock.add_component(conName, bigmConstraint) - # TODO unused code? - # parent_block = var.parent_block() - self._declare_disaggregated_var_bounds( original_var=var, disaggregatedVar=var, @@ -939,10 +916,10 @@ def _transform_constraint( name = c.local_name + "_%s" % unique # HACK: If I'm not treating fixed vars as permanent, I need - # the LinearRepnVisitor to pick them up, but it very - # strongly doesn't want to. So I'm going to unfix and refix - # which involves walking the expression twice, plus it's - # generally stupid + # 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): @@ -962,12 +939,12 @@ def _transform_constraint( # we substitute the expression variables with the # disaggregated variables if not NL or mode == "FurmanSawayaGrossmann": - h_x0 = clone_without_expression_components( - # TODO: is this properly handling fixed vars when we - # consider them permanent? Do we need a value() - # call? - c.body, - substitute=x0_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 @@ -1022,11 +999,10 @@ def _transform_constraint( # 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 that we didn't catch, - # but the feasible region may still be convex under - # weaker conditions, e.g. quasilinearity). But even - # so, this reformulation is still correct (though - # some of its purpose evaporates), so we will not + # 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: diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 5dd79e83640..550fa63872f 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2883,12 +2883,46 @@ def test_simple_case(self): m = models.makeTwoTermDisj() m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) TransformationFactory('gdp.hull').apply_to(m) - m.pprint() # did not throw # # TODO: keep on the model somewhere what x0 was, and assert it # has a nonzero element + def test_handle_fixed_disagg(self): + m = models.makeTwoTermDisj() + m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) + m.x.fix(2) + TransformationFactory('gdp.hull').apply_to(m, assume_fixed_vars_permanent=False) + # TODO: assert x0 has a key for m.x + + 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) + TransformationFactory('gdp.hull').apply_to(m, assume_fixed_vars_permanent=True) + # TODO: assert x0 does not have a key for m.x, also the + # transformed constraint here should be trivial + + def test_well_defined_points_arg(self): + m = models.makeTwoTermDisj() + m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) + TransformationFactory('gdp.hull').apply_to(m) + # did not throw + # + # TODO: check that well_defined_points is what is should be then + # redo by passing it + + 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") class NestedDisjunctsInFlatGDP(unittest.TestCase): From 25b6c559ae9cf41a9168a1847a4679d069586a0e Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 30 Sep 2025 18:51:09 -0600 Subject: [PATCH 06/20] Move LocalVars promotion and include "generalized local" vars --- pyomo/gdp/plugins/hull.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 4fb1aeae023..ccee262250c 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -630,6 +630,14 @@ def _transform_disjunctionData( 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. + # NOTE not doing this in _transform_disjunct anymore; verify this still works. + # Compared to the previous version, this could lead to adding variables on + # disactivated disjuncts to the parent's local vars, but this should not + # change anything. + 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 @@ -816,12 +824,9 @@ def _transform_disjunct( 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. - # - # TODO: this needs to add generalized local vars too. Why is - # this in transform_disjunct, can't I just do this centrally??? - parent_disjunct_local_vars.update(local_vars) + # NOTE: LocalVars promotion moved from here to caller; ensure + # nothing has gone horribly wrong + # deactivate disjunct so writers can be happy obj._deactivate_without_fixing_indicator() From ce00ea5ca1dc46b0c3b9b44a57cccf143d0e291e Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 30 Sep 2025 21:32:12 -0600 Subject: [PATCH 07/20] Introduce method for getting well-defined points with less dependence on solver particulars. Hopefully also less buggy. --- pyomo/gdp/plugins/hull.py | 120 +++++++++++++++++++++++++++++++++-- pyomo/gdp/tests/test_hull.py | 11 +++- 2 files changed, 123 insertions(+), 8 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index ccee262250c..b9efed923e6 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -27,6 +27,7 @@ BooleanVar, Connector, Constraint, + ConstraintList, Param, Set, SetOf, @@ -55,7 +56,10 @@ _warn_for_active_disjunct, ) from pyomo.core.util import target_list -from pyomo.core.expr.visitor import IdentifyVariableVisitor +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 @@ -390,6 +394,8 @@ def _get_well_defined_point( 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 " @@ -400,14 +406,20 @@ def _get_well_defined_point( for x in fallback_vars: x.fix(0) for x in regular_vars: + x.set_value(0) x.unfix() test_model = ConcreteModel() - test_model.obj = Objective(expr=test_expr) + 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: @@ -439,9 +451,7 @@ def _get_well_defined_point( # This one is going to get a little janky... # Let's start with just baron for the time being. def _solve_for_first_feasible_solution(self, test_model): - results = SolverFactory("baron").solve( - test_model, options={'NumSol': 1, 'FirstFeas': 1} - ) + results = SolverFactory("baron").solve(test_model, load_solutions=True) if results.solver.termination_condition is TerminationCondition.infeasible: return False if results.solver.status is not SolverStatus.ok: @@ -637,7 +647,7 @@ def _transform_disjunctionData( # disactivated disjuncts to the parent's local vars, but this should not # change anything. 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 @@ -826,7 +836,7 @@ def _transform_disjunct( # NOTE: LocalVars promotion moved from here to caller; ensure # nothing has gone horribly wrong - + # deactivate disjunct so writers can be happy obj._deactivate_without_fixing_indicator() @@ -1285,3 +1295,99 @@ 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): + self.cons_list.add(con) + + +# epsilon for handling function domains with strict inequalities +EPS = 1e-4 + +def _handlePowExpression(node): + (base, exp) = node.args + # 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). + + # Observe that this is problematic for LP and it needs to be a MIP! + # For now I'll just make it positive; this doesn't need to be perfect. + if exp.__class__ in EXPR.native_types or not exp.is_potentially_variable(): + val = value(exp) + rounded = round(val) + if rounded == 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 problem as before, this needs to be a MIP + 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 pi/2 plus a multiple of pi. Not much we can do about + # this one. + # if node.name == 'tan': + 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 = { + # Note: we will skip all the NPV expression types, since if one + # of those fails to be well-defined, changing variable values is + # not going to fix it + + # You're 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 550fa63872f..0b117c6d2e0 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2921,7 +2921,16 @@ def test_no_good_point(self): TransformationFactory('gdp.hull').apply_to, m, ) - + + def test_various_nonlinear(self): + m = models.makeTwoTermDisj() + m.y = Var(bounds=(-5, 5)) + m.d[0].log = Constraint(expr=log(m.y + 1) >= 0) + m.d[0].pow = Constraint(expr=m.y ** (-0.5) >= 0) + m.d[0].div = Constraint(expr=1 / (1 - m.y) >= 0) + TransformationFactory('gdp.hull').apply_to(m) + # did not throw + # TODO check that I properly handled the ill-defined stuff @unittest.skipUnless(gurobi_available, "Gurobi is not available") From 84e7b05bf961cb79301683ff45bb936d5f499aff Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Wed, 15 Oct 2025 08:55:15 -0600 Subject: [PATCH 08/20] fix an edge case --- pyomo/gdp/plugins/hull.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index b9efed923e6..4440fb33e13 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1314,14 +1314,22 @@ def __init__(self, cons_list, **kwds): 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 were supposed to have been filtered out during + # the handler call self.cons_list.add(con) # epsilon for handling function domains with strict inequalities 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 @@ -1384,7 +1392,6 @@ def _handleUnaryFunctionExpression(node): # Note: we will skip all the NPV expression types, since if one # of those fails to be well-defined, changing variable values is # not going to fix it - # You're on your own here # EXPR.ExternalFunctionExpression, EXPR.PowExpression: _handlePowExpression, From 91ac88084865fc71a17de416e768330ba1bcd526 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 30 Jan 2026 13:27:17 -0700 Subject: [PATCH 09/20] minor edits and made the walker class private --- pyomo/gdp/plugins/hull.py | 56 ++++++++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 21 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 4440fb33e13..df99db2cc51 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -79,6 +79,7 @@ class _HullTransformationData(AutoSlots.Mixin): 'original_var_map', 'bigm_constraint_map', 'disaggregation_constraint_map', + 'well_defined_points_map', ) def __init__(self): @@ -86,6 +87,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) @@ -183,8 +185,8 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): 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 serious - numerical issues. + small enough) and the original "LeeGrossmann" have numerical + and feasibility issues. References @@ -251,8 +253,7 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): 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 finally we raise GDP_Error if neither attempt was - successful. + call, then we raise GDP_Error if neither attempt was successful. """, ), ) @@ -417,7 +418,7 @@ def _get_well_defined_point( unique_component_name(test_model, x.name), Reference(x) ) test_model.well_defined_cons = ConstraintList() - WellDefinedConstraintGenerator( + _WellDefinedConstraintGenerator( cons_list=test_model.well_defined_cons ).walk_expression(test_expr) feasible = self._solve_for_first_feasible_solution(test_model) @@ -608,6 +609,7 @@ def _transform_disjunctionData( fallback_vars=all_local_vars, disj_name=obj.name, ) + transBlock.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 # breakpoint() @@ -1281,6 +1283,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._pyomo_gdp_hull_reformulation.private_data().well_defined_points_map + @TransformationFactory.register( 'gdp.chull', @@ -1303,7 +1317,7 @@ def __init__(self): # 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): +class _WellDefinedConstraintGenerator(StreamBasedExpressionVisitor): def __init__(self, cons_list, **kwds): self.cons_list = cons_list super().__init__(**kwds) @@ -1315,12 +1329,13 @@ 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 were supposed to have been filtered out during - # the handler call + # cases should have been filtered out during the handler + # call self.cons_list.add(con) -# epsilon for handling function domains with strict inequalities +# 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 @@ -1336,12 +1351,12 @@ def _handlePowExpression(node): # negative. Otherwise, base should be strictly positive (as we can't # be sure that exp could not be negative or fractional). - # Observe that this is problematic for LP and it needs to be a MIP! - # For now I'll just make it positive; this doesn't need to be perfect. + # 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" if exp.__class__ in EXPR.native_types or not exp.is_potentially_variable(): val = value(exp) - rounded = round(val) - if rounded == val: + if round(val) == val: if val >= 0: return () else: @@ -1357,7 +1372,7 @@ def _handleDivisionExpression(node): arg = node.args[1] if arg.__class__ in EXPR.native_types or not arg.is_potentially_variable(): return () - # Same problem as before, this needs to be a MIP + # Same LP vs MIP problem as before return (arg >= EPS,) @@ -1375,9 +1390,11 @@ def _handleUnaryFunctionExpression(node): return (arg >= -1, arg <= 1) if node.name == 'acos': return (arg >= -1, arg <= 1) - # It can't be pi/2 plus a multiple of pi. Not much we can do about - # this one. - # if node.name == 'tan': + # 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': @@ -1389,10 +1406,7 @@ def _handleUnaryFunctionExpression(node): # All expression types that can potentially be # ill-defined: _handlers = { - # Note: we will skip all the NPV expression types, since if one - # of those fails to be well-defined, changing variable values is - # not going to fix it - # You're on your own here + # You are on your own here # EXPR.ExternalFunctionExpression, EXPR.PowExpression: _handlePowExpression, EXPR.DivisionExpression: _handleDivisionExpression, From 7e3603ee6832ce4782e8d4f265d487c1510b60ab Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 30 Jan 2026 13:27:34 -0700 Subject: [PATCH 10/20] more properly test the hull changes --- pyomo/gdp/tests/test_hull.py | 126 ++++++++++++++++++++++++++++------- 1 file changed, 101 insertions(+), 25 deletions(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index fc1b998cfd3..371dc88a619 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -13,6 +13,7 @@ import sys import random from io import StringIO +import math import pyomo.common.unittest as unittest @@ -30,6 +31,14 @@ ComponentMap, value, log, + log10, + sqrt, + asin, + acos, + acosh, + atanh, + tan, + sin, ConcreteModel, Any, Suffix, @@ -38,6 +47,8 @@ Param, Objective, TerminationCondition, + SortComponents, + ConstraintList, ) from pyomo.core.expr.compare import ( assertExpressionsEqual, @@ -51,6 +62,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() @@ -2878,42 +2890,61 @@ def test_dill_pickle(self): 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) + 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) - # did not throw - # - # TODO: keep on the model somewhere what x0 was, and assert it - # has a nonzero element + + 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)) def test_handle_fixed_disagg(self): m = models.makeTwoTermDisj() m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) m.x.fix(2) - TransformationFactory('gdp.hull').apply_to(m, assume_fixed_vars_permanent=False) - # TODO: assert x0 has a key for m.x + 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]) 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) - TransformationFactory('gdp.hull').apply_to(m, assume_fixed_vars_permanent=True) - # TODO: assert x0 does not have a key for m.x, also the - # transformed constraint here should be trivial - - def test_well_defined_points_arg(self): - m = models.makeTwoTermDisj() - m.d[0].nonlinear = Constraint(expr=log(m.x - 1) >= 0) - TransformationFactory('gdp.hull').apply_to(m) - # did not throw - # - # TODO: check that well_defined_points is what is should be then - # redo by passing it + 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, + ) def test_no_good_point(self): m = models.makeTwoTermDisj() - m.d[0].nonlinear = Constraint(expr=log(-(m.x)**2 - 1) >= 0) + 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 .*", @@ -2924,13 +2955,15 @@ def test_no_good_point(self): def test_various_nonlinear(self): m = models.makeTwoTermDisj() m.y = Var(bounds=(-5, 5)) - m.d[0].log = Constraint(expr=log(m.y + 1) >= 0) + 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 / (1 - m.y) >= 0) - TransformationFactory('gdp.hull').apply_to(m) - # did not throw - # TODO check that I properly handled the ill-defined stuff - + 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): @@ -2940,3 +2973,46 @@ 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) + From 73c5af1c381c2f11415d42f0824071a368a079a1 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Fri, 30 Jan 2026 13:28:25 -0700 Subject: [PATCH 11/20] black --- pyomo/gdp/tests/test_hull.py | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 371dc88a619..be6a17fa74b 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2997,22 +2997,21 @@ def test_many_constraints(self): 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, + 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) - From 9abca2ef3f343f1cbdb3d4792c5e17212d82ca23 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 23 Mar 2026 12:54:01 -0600 Subject: [PATCH 12/20] remove a note --- pyomo/gdp/plugins/hull.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index df99db2cc51..6df05fb4786 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -644,10 +644,6 @@ def _transform_disjunctionData( ) # The parent disjunct's local vars should include all of this # disjunction's local or generalized local vars. - # NOTE not doing this in _transform_disjunct anymore; verify this still works. - # Compared to the previous version, this could lead to adding variables on - # disactivated disjuncts to the parent's local vars, but this should not - # change anything. local_vars_by_disjunct[parent_disjunct].update(generalized_local_vars) xorConstraint.add(index, (or_expr, 1)) From 08c0432b391c9e18c98f873b8908be71ec9ffdac Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 23 Mar 2026 12:56:06 -0600 Subject: [PATCH 13/20] cleanup --- pyomo/gdp/plugins/hull.py | 27 +-------------------------- 1 file changed, 1 insertion(+), 26 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 6df05fb4786..0e062f87988 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -554,28 +554,6 @@ def _transform_disjunctionData( 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: - # 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) - # all_local_vars.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) for disj in disjuncts: if disj in local_vars_by_disjunct: if var in local_vars_by_disjunct[disj]: @@ -832,9 +810,6 @@ def _transform_disjunct( obj, obj, var_substitute_map, x0_substitute_map ) - # NOTE: LocalVars promotion moved from here to caller; ensure - # nothing has gone horribly wrong - # deactivate disjunct so writers can be happy obj._deactivate_without_fixing_indicator() @@ -1349,7 +1324,7 @@ def _handlePowExpression(node): # 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" + # 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: From b7fa95f91139f9340ef2beb99797d55bd6a90bd9 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 23 Mar 2026 13:14:37 -0600 Subject: [PATCH 14/20] Switch to using gurobi for the well-defined point heuristic, since it's probably more likely to be installed than baron --- pyomo/gdp/plugins/hull.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 21243a64162..83202565453 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -434,9 +434,9 @@ def _get_well_defined_point( "actually exists, then if we still cannot find it, override our " "search process using the `well_defined_points` option." ) - # We have a point. + # Found a point if not math.isfinite(value(test_expr)): - raise DeveloperError("Theoretically unreachable") + raise DeveloperError("unreachable") x0_map = ComponentMap() used_vars = ComponentSet() for x in itertools.chain(regular_vars, fallback_vars): @@ -447,14 +447,17 @@ def _get_well_defined_point( x.fixed = orig_fixed[x] return x0_map, used_vars - # This one is going to get a little janky... - # Let's start with just baron for the time being. + # 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 = SolverFactory("baron").solve(test_model, load_solutions=True) + results = SolverFactory("gurobi_direct_minlp").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( From 7be9e9860f9da8ac801544b75e0e2c30c1f461b7 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Mon, 23 Mar 2026 15:54:46 -0600 Subject: [PATCH 15/20] check gurobi availability in tests --- pyomo/gdp/plugins/hull.py | 17 +++++++++++++++-- pyomo/gdp/tests/test_hull.py | 5 +++++ 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 83202565453..3c2a6b7a62b 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -53,6 +53,7 @@ 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, @@ -255,6 +256,18 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): """, ), ) + 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 + nonlinear solver in general. + """, + ), + ) transformation_name = 'hull' def __init__(self): @@ -450,7 +463,7 @@ def _get_well_defined_point( # 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 = SolverFactory("gurobi_direct_minlp").solve( + results = self._config.well_defined_points_heuristic_solver.solve( test_model, load_solutions=False ) if results.solver.termination_condition is TerminationCondition.infeasible: @@ -1312,7 +1325,7 @@ def exitNode(self, node, data): def _handlePowExpression(node): - (base, exp) = node.args + 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 () diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 906031dfa6b..cbaf40bf073 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2897,6 +2897,7 @@ def test_no_nonlinear(self): 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) @@ -2916,6 +2917,7 @@ def test_simple_case(self): ): 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) @@ -2924,6 +2926,7 @@ def test_handle_fixed_disagg(self): 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) @@ -2940,6 +2943,7 @@ def test_handle_fixed_no_disagg(self): 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) @@ -2950,6 +2954,7 @@ def test_no_good_point(self): m, ) + @unittest.skipUnless(gurobi_available, "Gurobi is not available") def test_various_nonlinear(self): m = models.makeTwoTermDisj() m.y = Var(bounds=(-5, 5)) From c22212912e942436e3e810ef637480154e333e03 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Mar 2026 16:59:03 -0600 Subject: [PATCH 16/20] cleanup --- pyomo/gdp/plugins/hull.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 3c2a6b7a62b..7f353069a01 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -448,8 +448,6 @@ def _get_well_defined_point( "search process using the `well_defined_points` option." ) # Found a point - if not math.isfinite(value(test_expr)): - raise DeveloperError("unreachable") x0_map = ComponentMap() used_vars = ComponentSet() for x in itertools.chain(regular_vars, fallback_vars): @@ -604,7 +602,6 @@ def _transform_disjunctionData( transBlock.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 - # breakpoint() for var in offset_vars: if var in all_local_vars: var_disjunct = next(iter(disjuncts_var_appears_in[var])) @@ -705,7 +702,6 @@ def _transform_disjunctionData( else: disaggregatedExpr = 0 for disjunct in disjuncts_var_appears_in[var]: - # breakpoint() disaggregatedExpr += disjunct_disaggregated_var_map[disjunct][var] cons_idx = len(disaggregationConstraint) @@ -736,7 +732,6 @@ def _transform_disjunct( x0_map, offset_vars, ): - # breakpoint() relaxationBlock = self._get_disjunct_transformation_block(obj, transBlock) # Put the disaggregated variables all on their own block so that we can From b8575974de97ce1cdc373e9f38f1863b36e46dbb Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Mar 2026 17:03:58 -0600 Subject: [PATCH 17/20] improve doc for the solver parameter --- pyomo/gdp/plugins/hull.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 7f353069a01..08dba2bc4eb 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -264,7 +264,9 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): description="Solver used for the base points heuristic", doc=""" Solver used for the base points heuristic. This must be a - nonlinear solver in general. + nonconvex NLP solver in general. Pass as a solver object + supporting the V1 solver API version, or a corresponding string + for SolverFactory. """, ), ) From b386530887378ad039672a7426ef43ad491c3c75 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Mar 2026 17:10:54 -0600 Subject: [PATCH 18/20] remove commentary on history of codebase --- pyomo/gdp/plugins/hull.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 08dba2bc4eb..3e5d3757860 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -572,7 +572,6 @@ def _transform_disjunctionData( if disj in local_vars_by_disjunct: if var in local_vars_by_disjunct[disj]: if len(disjuncts) == 1: - # was this a noop before? local_vars[disj].add(var) all_local_vars.add(var) else: From 5ce5434acc8dc4cc04b2e50ac05793c97c0109dc Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Mar 2026 18:05:10 -0600 Subject: [PATCH 19/20] multiplier is always 1 here --- pyomo/gdp/plugins/hull.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 3e5d3757860..9bb894ea288 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -985,7 +985,11 @@ def _transform_constraint( # 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. - expr = linear_repn.multiplier * sum( + # + # 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 @@ -1010,10 +1014,7 @@ def _transform_constraint( # the constraint implies x = x0 if ( coef != 0 - and ( - c.lower - linear_repn.constant * linear_repn.multiplier - ) - / (coef * linear_repn.multiplier) + and (c.lower - linear_repn.constant) / coef == x0_substitute_map[var] ): v = var_substitute_map[var] From b074f67e90e98b356040df1a8825db4cfde0a741 Mon Sep 17 00:00:00 2001 From: Soren Davis Date: Tue, 24 Mar 2026 18:29:39 -0600 Subject: [PATCH 20/20] move well_defined_points_map to a more suitable location --- pyomo/gdp/plugins/hull.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 9bb894ea288..0ac105fc7ba 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -600,7 +600,7 @@ def _transform_disjunctionData( fallback_vars=all_local_vars, disj_name=obj.name, ) - transBlock.private_data().well_defined_points_map[obj] = x0_map + 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: @@ -1275,7 +1275,7 @@ def get_well_defined_points_map(self, b): ---------- b: a Block that was transformed by gdp.hull """ - return b._pyomo_gdp_hull_reformulation.private_data().well_defined_points_map + return b.private_data().well_defined_points_map @TransformationFactory.register(