From 651253bc69e92ae1c6f8a28cf1b29a7e48d2a5e4 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Fri, 6 Mar 2026 15:19:00 -0500 Subject: [PATCH 01/39] added exact hull reformulations --- pyomo/gdp/plugins/hull.py | 430 +++++++++++++++++++++++++++++++++++++- 1 file changed, 424 insertions(+), 6 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index f7e7d7a4b9d..a5b80620095 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -7,10 +7,29 @@ # software. This software is distributed under the 3-clause BSD License. # ____________________________________________________________________________________ +"""Hull reformulation for Generalized Disjunctive Programming. + +This module implements the hull reformulation (also known as the convex hull +relaxation) for disjunctive models in Pyomo's GDP framework. It transforms +disjunctive constraints into algebraic form using variable disaggregation and +perspective functions. + +When the ``exact_hull_quadratic`` configuration option is enabled, quadratic +constraints receive an exact hull treatment: + +* **Conic exact hull** -- for convex quadratics, an auxiliary variable and a + rotated second-order cone constraint replace the standard perspective + function, yielding a tighter relaxation without a Cholesky factorisation. +* **General exact hull** -- for non-convex (or equality) quadratics, the + disaggregated quadratic form ``v'Qv + c'v y + d y**2`` is used directly. +""" + import logging from collections import defaultdict +import numpy as np + from pyomo.common.autoslots import AutoSlots import pyomo.common.config as cfg from pyomo.common import deprecated @@ -19,6 +38,7 @@ from pyomo.core.expr.numvalue import ZeroConstant import pyomo.core.expr as EXPR from pyomo.core.base import TransformationFactory +from pyomo.repn import generate_standard_repn from pyomo.core import ( Block, BooleanVar, @@ -35,6 +55,7 @@ Any, RangeSet, Reals, + NonNegativeReals, value, NonNegativeIntegers, Binary, @@ -197,6 +218,39 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): """, ), ) + CONFIG.declare( + 'exact_hull_quadratic', + cfg.ConfigValue( + default=False, + domain=bool, + description="Use exact hull reformulation for quadratic constraints.", + doc=""" + If True, quadratic constraints (polynomial degree 2) are reformulated + using the exact hull instead of the standard perspective function. + + For a quadratic constraint of the form + + x'Qx + c'x + d <= 0 + + the reformulation depends on convexity: + + **Conic exact hull** (convex quadratics): An auxiliary variable *t* + and a rotated second-order cone constraint ``v'Qv <= t * y`` are + introduced, and the original bound becomes ``t + c'v + d*y <= 0``. + Convexity is determined via eigenvalue decomposition of the Hessian + matrix *Q*: the quadratic is convex for an upper-bound constraint + when *Q* is positive semi-definite, and for a lower-bound constraint + when *Q* is negative semi-definite. + + **General exact hull** (non-convex quadratics and equalities): The + constraint is reformulated as ``v'Qv + c'v*y + d*y**2``, where *v* + are the disaggregated variables and *y* is the binary indicator. + + Default is False, which uses the standard perspective function for + all nonlinear constraints. + """, + ), + ) transformation_name = 'hull' def __init__(self): @@ -658,6 +712,25 @@ def _get_local_var_list(self, parent_disjunct): def _transform_constraint( self, obj, disjunct, var_substitute_map, zero_substitute_map ): + """Transform a single Constraint on a Disjunct. + + Applies the appropriate hull reformulation to each + ``ConstraintData`` in *obj*. When ``exact_hull_quadratic`` is + enabled and the constraint body has polynomial degree 2, an exact + hull formulation is used instead of the perspective function. + + Parameters + ---------- + obj : Constraint + The Constraint component to transform. + disjunct : Disjunct + The Disjunct that owns *obj*. + var_substitute_map : dict + Mapping from ``id(original_var)`` to its disaggregated + counterpart. + zero_substitute_map : dict + Mapping from ``id(original_var)`` to ``ZeroConstant``. + """ # we will put a new transformed constraint on the relaxation block. relaxationBlock = disjunct._transformation_block() constraint_map = relaxationBlock.private_data('pyomo.gdp') @@ -676,9 +749,17 @@ def _transform_constraint( unique = len(newConstraint) name = c.local_name + "_%s" % unique - NL = c.body.polynomial_degree() not in (0, 1) + polynomial_degree = c.body.polynomial_degree() EPS = self._config.EPS mode = self._config.perspective_function + exact_quad = self._config.exact_hull_quadratic + + use_exact_quad = ( + exact_quad and polynomial_degree == 2 + ) + + NL = polynomial_degree not in (0, 1) + general_NL = NL and not use_exact_quad # We need to evaluate the expression at the origin *before* # we substitute the expression variables with the @@ -689,12 +770,21 @@ def _transform_constraint( ) y = disjunct.binary_indicator_var - if NL: + + if use_exact_quad: + self._build_exact_quadratic_hull( + c, y, disjunct, relaxationBlock, constraint_map, + var_substitute_map, newConstraint, name, i, obj, + ) + continue + + if general_NL: if mode == "LeeGrossmann": sub_expr = clone_without_expression_components( c.body, substitute=dict( - (var, subs / y) for var, subs in var_substitute_map.items() + (var, subs / y) + for var, subs in var_substitute_map.items() ), ) expr = sub_expr * y @@ -724,7 +814,7 @@ def _transform_constraint( ) if c.equality: - if NL: + if general_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? @@ -775,7 +865,7 @@ def _transform_constraint( if self._generate_debug_messages: _name = c.getname(fully_qualified=True) logger.debug("GDP(Hull): Transforming constraint " + "'%s'", _name) - if NL: + if general_NL: newConsExpr = expr >= c.lower * y else: newConsExpr = expr - (1 - y) * h_0 >= c.lower * y @@ -797,7 +887,7 @@ def _transform_constraint( if self._generate_debug_messages: _name = c.getname(fully_qualified=True) logger.debug("GDP(Hull): Transforming constraint " + "'%s'", _name) - if NL: + if general_NL: newConsExpr = expr <= c.upper * y else: newConsExpr = expr - (1 - y) * h_0 <= c.upper * y @@ -820,6 +910,334 @@ def _transform_constraint( # deactivate now that we have transformed obj.deactivate() + def _build_exact_quadratic_hull( + self, c, y, disjunct, relaxationBlock, constraint_map, + var_substitute_map, newConstraint, name, i, obj, + ): + """Build the exact hull reformulation for a single quadratic constraint. + + For a constraint whose body is a quadratic of the form + ``x'Qx + c'x + d``, this method constructs either the *conic exact + hull* (when the quadratic is convex with respect to the bound + direction) or the *general exact hull* (otherwise). + + **Conic exact hull** (convex case): introduces an auxiliary variable + *t >= 0* and a rotated second-order cone constraint + ``v'Q_psd v <= t * y``, then replaces the original bound with a + linear constraint on ``t + c'v + d*y``. + + **General exact hull** (non-convex / equality case): directly + substitutes the quadratic form to ``v'Qv + c'v*y + d*y**2``. + + Parameters + ---------- + c : ConstraintData + The individual constraint data object being transformed. + y : Var + The binary indicator variable for the parent disjunct. + disjunct : Disjunct + The Disjunct that owns the constraint. + relaxationBlock : Block + The transformation block for this disjunct. + constraint_map : object + Private data object tracking constraint mappings. + var_substitute_map : dict + Mapping from ``id(original_var)`` to disaggregated variable. + newConstraint : Constraint + The indexed Constraint container for transformed constraints. + name : str + Base name for the transformed constraint indices. + i : object + The index of the constraint in its parent component. + obj : Constraint + The parent Constraint component (needed for ``is_indexed``). + """ + repn = generate_standard_repn(c.body) + + if not repn.is_quadratic(): + raise RuntimeError( + "Constraint '%s' has polynomial degree 2 but its standard " + "representation is not quadratic." % c.getname(fully_qualified=True) + ) + + const_term = repn.constant if repn.constant is not None else 0 + + # --- Build the symmetric Q matrix and determine convexity --- + all_vars = [] + var_ids_seen = set() + for var_i, var_j in repn.quadratic_vars: + if id(var_i) not in var_ids_seen: + all_vars.append(var_i) + var_ids_seen.add(id(var_i)) + if id(var_j) not in var_ids_seen: + all_vars.append(var_j) + var_ids_seen.add(id(var_j)) + + n_vars = len(all_vars) + var_to_idx = {id(var): idx for idx, var in enumerate(all_vars)} + Q = np.zeros((n_vars, n_vars)) + + for coef, (var_i, var_j) in zip( + repn.quadratic_coefs, repn.quadratic_vars + ): + idx_i = var_to_idx[id(var_i)] + idx_j = var_to_idx[id(var_j)] + if var_i is var_j: + Q[idx_i, idx_i] += coef + else: + Q[idx_i, idx_j] += 0.5 * coef + Q[idx_j, idx_i] += 0.5 * coef + + numerical_tolerance = 1e-10 + eigenvalues, _ = np.linalg.eigh(Q) + Q_is_psd = not np.any(eigenvalues < -numerical_tolerance) + Q_is_nsd = not np.any(eigenvalues > numerical_tolerance) + + # Determine which bounds can use the conic formulation + use_conic_upper = False + use_conic_lower = False + negate_for_conic = False + + if c.upper is not None and not c.equality: + if Q_is_psd: + use_conic_upper = True + if c.lower is not None and not c.equality: + if Q_is_nsd: + use_conic_lower = True + negate_for_conic = True + + if negate_for_conic: + Q_for_conic = -Q + else: + Q_for_conic = Q + + # --- Decide which expression forms are needed --- + need_non_convex = False + if c.equality: + need_non_convex = True + if c.upper is not None and not use_conic_upper: + need_non_convex = True + if c.lower is not None and not use_conic_lower: + need_non_convex = True + + non_conv_expr = None + conic_expr_linear = None + + if need_non_convex: + non_conv_expr = self._build_general_exact_hull_expr( + repn, var_substitute_map, y, const_term, + ) + + if use_conic_upper or use_conic_lower: + conic_expr_linear = self._build_conic_exact_hull_expr( + c, y, disjunct, relaxationBlock, constraint_map, + repn, var_substitute_map, const_term, + negate_for_conic, + ) + + # --- Equality constraints always use general exact hull --- + if c.equality: + newConsExpr = non_conv_expr == c.lower * y**2 + + if obj.is_indexed(): + newConstraint.add((name, i, 'eq'), newConsExpr) + constraint_map.transformed_constraints[c].append( + newConstraint[name, i, 'eq'] + ) + constraint_map.src_constraint[newConstraint[name, i, 'eq']] = c + else: + newConstraint.add((name, 'eq'), newConsExpr) + constraint_map.transformed_constraints[c].append( + newConstraint[name, 'eq'] + ) + constraint_map.src_constraint[newConstraint[name, 'eq']] = c + return + + # --- Lower bound --- + if c.lower is not None: + if self._generate_debug_messages: + _name = c.getname(fully_qualified=True) + logger.debug("GDP(Hull): Transforming constraint '%s'", _name) + + if use_conic_lower: + newConsExpr = conic_expr_linear <= -c.lower * y + else: + newConsExpr = non_conv_expr >= c.lower * y**2 + + if obj.is_indexed(): + newConstraint.add((name, i, 'lb'), newConsExpr) + constraint_map.transformed_constraints[c].append( + newConstraint[name, i, 'lb'] + ) + constraint_map.src_constraint[newConstraint[name, i, 'lb']] = c + else: + newConstraint.add((name, 'lb'), newConsExpr) + constraint_map.transformed_constraints[c].append( + newConstraint[name, 'lb'] + ) + constraint_map.src_constraint[newConstraint[name, 'lb']] = c + + # --- Upper bound --- + if c.upper is not None: + if self._generate_debug_messages: + _name = c.getname(fully_qualified=True) + logger.debug("GDP(Hull): Transforming constraint '%s'", _name) + + if use_conic_upper: + newConsExpr = conic_expr_linear <= c.upper * y + else: + newConsExpr = non_conv_expr <= c.upper * y**2 + + if obj.is_indexed(): + newConstraint.add((name, i, 'ub'), newConsExpr) + constraint_map.transformed_constraints[c].append( + newConstraint[name, i, 'ub'] + ) + constraint_map.src_constraint[newConstraint[name, i, 'ub']] = c + else: + newConstraint.add((name, 'ub'), newConsExpr) + constraint_map.transformed_constraints[c].append( + newConstraint[name, 'ub'] + ) + constraint_map.src_constraint[newConstraint[name, 'ub']] = c + + def _build_general_exact_hull_expr( + self, repn, var_substitute_map, y, const_term + ): + """Build the general exact hull expression for a quadratic constraint. + + Constructs the expression ``v'Qv + c'v*y + d*y**2`` where *v* are + disaggregated variables, *y* is the indicator, *c* are linear + coefficients, and *d* is the constant term. + + Parameters + ---------- + repn : StandardRepn + The standard representation of the constraint body. + var_substitute_map : dict + Mapping from ``id(original_var)`` to disaggregated variable. + y : Var + The binary indicator variable. + const_term : float + The constant term from the standard representation. + + Returns + ------- + expression + The general exact hull Pyomo expression. + """ + expr = 0 + + for coef, (var_i, var_j) in zip( + repn.quadratic_coefs, repn.quadratic_vars + ): + v_i = var_substitute_map.get(id(var_i), var_i) + v_j = var_substitute_map.get(id(var_j), var_j) + if var_i is var_j: + expr += coef * v_i**2 + else: + expr += coef * v_i * v_j + + lin_coefs = repn.linear_coefs or [] + lin_vars = repn.linear_vars or [] + for coef, var in zip(lin_coefs, lin_vars): + v = var_substitute_map.get(id(var), var) + expr += coef * v * y + + if const_term: + expr += const_term * y**2 + + return expr + + def _build_conic_exact_hull_expr( + self, c, y, disjunct, relaxationBlock, constraint_map, + repn, var_substitute_map, const_term, negate_for_conic, + ): + """Build the conic exact hull expression for a convex quadratic. + + Creates an auxiliary variable ``t >= 0`` and a rotated second-order + cone constraint ``v'Q_psd v <= t * y``, then returns the linear + expression ``t + c'v + d*y`` (with signs adjusted for lower-bound + constraints that required negation). + + Parameters + ---------- + c : ConstraintData + The constraint data being transformed. + y : Var + The binary indicator variable. + disjunct : Disjunct + The parent Disjunct. + relaxationBlock : Block + The transformation block for the disjunct. + constraint_map : object + Private data tracking constraint mappings. + repn : StandardRepn + The standard representation of the constraint body. + var_substitute_map : dict + Mapping from ``id(original_var)`` to disaggregated variable. + const_term : float + The constant term from the standard representation. + negate_for_conic : bool + If ``True``, coefficients are negated (used when the + lower-bound constraint is reformulated by negation to obtain a + PSD form). + + Returns + ------- + expression + The linear Pyomo expression ``t + c'v + d*y`` (or its negated + variant) to be bounded by the constraint's RHS. + """ + t = Var(domain=NonNegativeReals) + t_name = unique_component_name( + relaxationBlock, + '_conic_aux_t_%s' % c.getname( + fully_qualified=True, relative_to=disjunct + ), + ) + relaxationBlock.add_component(t_name, t) + + linear_expr = t + + lin_coefs = repn.linear_coefs or [] + lin_vars = repn.linear_vars or [] + for coef, var in zip(lin_coefs, lin_vars): + v = var_substitute_map.get(id(var), var) + actual_coef = -coef if negate_for_conic else coef + linear_expr += actual_coef * v + + if const_term: + actual_const = -const_term if negate_for_conic else const_term + linear_expr += actual_const * y + + # Build rotated SOC: v'Q_psd v <= t * y + quadratic_form = 0 + for coef, (var_i, var_j) in zip( + repn.quadratic_coefs, repn.quadratic_vars + ): + v_i = var_substitute_map.get(id(var_i), var_i) + v_j = var_substitute_map.get(id(var_j), var_j) + actual_coef = -coef if negate_for_conic else coef + if var_i is var_j: + quadratic_form += actual_coef * v_i**2 + else: + quadratic_form += actual_coef * v_i * v_j + + conic_constraint_name = unique_component_name( + relaxationBlock, + '_conic_constraint_%s' % c.getname( + fully_qualified=True, relative_to=disjunct + ), + ) + conic_constraint = Constraint(expr=quadratic_form <= t * y) + relaxationBlock.add_component(conic_constraint_name, conic_constraint) + + constraint_map.transformed_constraints[c].append(conic_constraint) + constraint_map.src_constraint[conic_constraint] = c + + return linear_expr + def _get_local_var_suffix(self, disjunct): # If the Suffix is there, we will borrow it. If not, we make it. If it's # something else, we complain. From 4b53b335c84ca3c3ce41f7b15e7a28ccb9fda5f5 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Fri, 6 Mar 2026 21:19:23 -0500 Subject: [PATCH 02/39] changed comments --- pyomo/gdp/plugins/hull.py | 20 +------------------- 1 file changed, 1 insertion(+), 19 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index a5b80620095..678459f6460 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -7,23 +7,6 @@ # software. This software is distributed under the 3-clause BSD License. # ____________________________________________________________________________________ -"""Hull reformulation for Generalized Disjunctive Programming. - -This module implements the hull reformulation (also known as the convex hull -relaxation) for disjunctive models in Pyomo's GDP framework. It transforms -disjunctive constraints into algebraic form using variable disaggregation and -perspective functions. - -When the ``exact_hull_quadratic`` configuration option is enabled, quadratic -constraints receive an exact hull treatment: - -* **Conic exact hull** -- for convex quadratics, an auxiliary variable and a - rotated second-order cone constraint replace the standard perspective - function, yielding a tighter relaxation without a Cholesky factorisation. -* **General exact hull** -- for non-convex (or equality) quadratics, the - disaggregated quadratic form ``v'Qv + c'v y + d y**2`` is used directly. -""" - import logging from collections import defaultdict @@ -783,8 +766,7 @@ 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) for var, subs in var_substitute_map.items() ), ) expr = sub_expr * y From ab99b79fd38df1089b9e3ced904432b78c12cd14 Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Fri, 6 Mar 2026 23:25:19 -0500 Subject: [PATCH 03/39] Update pyomo/gdp/plugins/hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- 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 678459f6460..076e593fab3 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -988,11 +988,6 @@ def _build_exact_quadratic_hull( use_conic_lower = True negate_for_conic = True - if negate_for_conic: - Q_for_conic = -Q - else: - Q_for_conic = Q - # --- Decide which expression forms are needed --- need_non_convex = False if c.equality: From 37b9d811a5763a19d22c921580194a2054db0d22 Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Fri, 6 Mar 2026 23:34:53 -0500 Subject: [PATCH 04/39] Update pyomo/gdp/plugins/hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 076e593fab3..aa1aca91bd4 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -11,7 +11,7 @@ from collections import defaultdict -import numpy as np +from pyomo.common.dependencies import numpy as np, numpy_available from pyomo.common.autoslots import AutoSlots import pyomo.common.config as cfg From da67052dcfe0c34462dfe6ba450d86246cb46340 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Fri, 6 Mar 2026 23:40:29 -0500 Subject: [PATCH 05/39] added np available check --- pyomo/gdp/plugins/hull.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index aa1aca91bd4..8e1c8124f36 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -944,6 +944,12 @@ def _build_exact_quadratic_hull( const_term = repn.constant if repn.constant is not None else 0 + if not numpy_available: + raise GDP_Error( + "exact_hull_quadratic requires NumPy for convexity checks. " + "NumPy is not available in this environment." + ) + # --- Build the symmetric Q matrix and determine convexity --- all_vars = [] var_ids_seen = set() From 471a020155026c82a9d72e2a6fd520aea7c13043 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 7 Mar 2026 04:44:24 +0000 Subject: [PATCH 06/39] Initial plan From 14d76428e1cf25c4a49da5c66df49b49e1ba7008 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 7 Mar 2026 04:46:22 +0000 Subject: [PATCH 07/39] Make eigenvalue_tolerance a configurable parameter in hull.py Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 8e1c8124f36..da9095a6175 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -105,6 +105,15 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): 'LeeGrossmann', or 'GrossmannLee' EPS : float The value to use for epsilon [default: 1e-4] + eigenvalue_tolerance : float + Numerical tolerance for eigenvalue-based positive/negative + semi-definite checks when using the exact hull reformulation for + quadratic constraints (``exact_hull_quadratic=True``). An + eigenvalue is considered non-negative (non-positive) if it is + greater than (less than) ``-eigenvalue_tolerance`` + (``eigenvalue_tolerance``). Increasing this value makes the + convexity check more conservative; decreasing it makes it more + permissive. [default: 1e-10] targets : block, disjunction, or list of those types The targets to transform. This can be a block, disjunction, or a list of blocks and Disjunctions [default: the instance] @@ -184,6 +193,29 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): description="Epsilon value to use in perspective function", ), ) + CONFIG.declare( + 'eigenvalue_tolerance', + cfg.ConfigValue( + default=1e-10, + domain=cfg.PositiveFloat, + description="Numerical tolerance for eigenvalue-based PSD/NSD checks " + "in exact hull quadratic reformulations", + doc=""" + Numerical tolerance used when determining positive semi-definiteness + (PSD) or negative semi-definiteness (NSD) of the Hessian matrix Q in + the exact hull reformulation for quadratic constraints + (``exact_hull_quadratic=True``). + + An eigenvalue ``lam`` is treated as non-negative if + ``lam >= -eigenvalue_tolerance``, and non-positive if + ``lam <= eigenvalue_tolerance``. Increasing this value makes the + convexity classification more conservative (i.e., Q must have + eigenvalues further from zero to be considered PSD/NSD); decreasing + it makes the check more permissive. For ill-conditioned Q matrices a + larger tolerance may be appropriate. + """, + ), + ) CONFIG.declare( 'assume_fixed_vars_permanent', cfg.ConfigValue( @@ -976,7 +1008,7 @@ def _build_exact_quadratic_hull( Q[idx_i, idx_j] += 0.5 * coef Q[idx_j, idx_i] += 0.5 * coef - numerical_tolerance = 1e-10 + numerical_tolerance = self._config.eigenvalue_tolerance eigenvalues, _ = np.linalg.eigh(Q) Q_is_psd = not np.any(eigenvalues < -numerical_tolerance) Q_is_nsd = not np.any(eigenvalues > numerical_tolerance) From 0c7e4460b43a8d5ce9d635470d73691e6c2417f1 Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Sat, 7 Mar 2026 00:18:37 -0500 Subject: [PATCH 08/39] Update pyomo/gdp/plugins/hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index da9095a6175..b6334aab4ad 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -209,10 +209,11 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): An eigenvalue ``lam`` is treated as non-negative if ``lam >= -eigenvalue_tolerance``, and non-positive if ``lam <= eigenvalue_tolerance``. Increasing this value makes the - convexity classification more conservative (i.e., Q must have - eigenvalues further from zero to be considered PSD/NSD); decreasing - it makes the check more permissive. For ill-conditioned Q matrices a - larger tolerance may be appropriate. + convexity classification more permissive (i.e., a wider band around + zero is treated as numerically zero, so more eigenvalues are accepted + as PSD/NSD); decreasing it makes the check more conservative (i.e., + eigenvalues must be further from zero). For ill-conditioned Q matrices + a larger tolerance may be appropriate. """, ), ) From a885997031905ec1611814e0a86c678b96a1ad3b Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Sat, 7 Mar 2026 00:20:02 -0500 Subject: [PATCH 09/39] Update pyomo/gdp/plugins/hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- 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 b6334aab4ad..4f99e6fc77b 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -109,11 +109,14 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): Numerical tolerance for eigenvalue-based positive/negative semi-definite checks when using the exact hull reformulation for quadratic constraints (``exact_hull_quadratic=True``). An - eigenvalue is considered non-negative (non-positive) if it is - greater than (less than) ``-eigenvalue_tolerance`` - (``eigenvalue_tolerance``). Increasing this value makes the - convexity check more conservative; decreasing it makes it more - permissive. [default: 1e-10] + eigenvalue :math:`\lambda` is treated as non-negative if + :math:`\lambda >= -\text{eigenvalue_tolerance}` and as + non-positive if :math:`\lambda <= \text{eigenvalue_tolerance}` + (i.e., eigenvalues in + ``[-eigenvalue_tolerance, eigenvalue_tolerance]`` are treated as + zero). Increasing this value makes the convexity check more + permissive; decreasing it makes it more conservative. + [default: 1e-10] targets : block, disjunction, or list of those types The targets to transform. This can be a block, disjunction, or a list of blocks and Disjunctions [default: the instance] From 567730fb34da5030803e0355e4c798d4020adff8 Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Sat, 7 Mar 2026 00:20:24 -0500 Subject: [PATCH 10/39] Update pyomo/gdp/plugins/hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 4f99e6fc77b..17bd022233f 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -200,7 +200,7 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): 'eigenvalue_tolerance', cfg.ConfigValue( default=1e-10, - domain=cfg.PositiveFloat, + domain=cfg.NonNegativeFloat, description="Numerical tolerance for eigenvalue-based PSD/NSD checks " "in exact hull quadratic reformulations", doc=""" From 50b3f5c422821848921950ea7d5bdd78aa431d33 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 05:03:31 +0000 Subject: [PATCH 11/39] Initial plan From 800af39513acbc25489ad9b50d492bb13695b486 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 05:10:40 +0000 Subject: [PATCH 12/39] Initial plan From 802112d689bc97dc2972fb7049cc003b4e5437ff Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 05:11:58 +0000 Subject: [PATCH 13/39] Add exact_hull_quadratic to class docstring before eigenvalue_tolerance Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 17bd022233f..adf5f4f81db 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -105,6 +105,14 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): 'LeeGrossmann', or 'GrossmannLee' EPS : float The value to use for epsilon [default: 1e-4] + exact_hull_quadratic : bool + If ``True``, quadratic constraints (polynomial degree 2) are + reformulated using the exact hull instead of the standard + perspective function. Convex quadratics are handled with a conic + reformulation (rotated second-order cone), while non-convex + quadratics and equalities use the general exact hull + reformulation. Convexity is determined via eigenvalue + decomposition of the Hessian matrix. [default: False] eigenvalue_tolerance : float Numerical tolerance for eigenvalue-based positive/negative semi-definite checks when using the exact hull reformulation for From 6db6ec90c95c3ce15c00e6b7e84a4fca7cd52986 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 05:14:39 +0000 Subject: [PATCH 14/39] Add comprehensive TestExactHullQuadratic tests for exact_hull_quadratic feature Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/tests/test_hull.py | 563 ++++++++++++++++++++++++++++++++++- 1 file changed, 562 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 8c59ef0bb55..3e60239a082 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -13,8 +13,9 @@ from io import StringIO import pyomo.common.unittest as unittest +import unittest.mock as mock -from pyomo.common.dependencies import dill_available +from pyomo.common.dependencies import dill_available, numpy_available from pyomo.common.log import LoggingIntercept from pyomo.common.fileutils import this_file_dir @@ -36,6 +37,7 @@ Param, Objective, TerminationCondition, + NonNegativeReals, ) from pyomo.core.expr.compare import ( assertExpressionsEqual, @@ -47,6 +49,7 @@ from pyomo.repn.linear import LinearRepnVisitor from pyomo.gdp import Disjunct, Disjunction, GDP_Error +import pyomo.gdp.plugins.hull as hull_module import pyomo.gdp.tests.models as models import pyomo.gdp.tests.common_tests as ct @@ -2890,3 +2893,561 @@ class NestedDisjunctsInFlatGDP(unittest.TestCase): def test_declare_disjuncts_in_disjunction_rule(self): ct.check_nested_disjuncts_in_flat_gdp(self, 'hull') + + +@unittest.skipUnless(numpy_available, "NumPy is not available") +class TestExactHullQuadratic(unittest.TestCase): + """Tests for the ``exact_hull_quadratic`` option of the hull transformation. + + Covers: + - Convex quadratic upper-bound (conic reformulation, PSD Q) + - Convex quadratic lower-bound (conic reformulation with negation, NSD Q) + - Non-convex quadratic (general exact hull, mixed eigenvalues) + - Equality constraint (always general exact hull) + - Range constraint with both bounds (mixed conic + general) + - ``numpy_available`` guard + - ``eigenvalue_tolerance`` configuration + - ``get_transformed_constraints`` mapping + - Linear and constant terms in both conic and general reformulations + """ + + def setUp(self): + self.hull = TransformationFactory('gdp.hull') + + # ------------------------------------------------------------------ + # 1. Convex quadratic upper-bound → conic reformulation + # ------------------------------------------------------------------ + def test_convex_quadratic_upper_bound_conic_reformulation(self): + """PSD Q with upper bound should produce a conic reformulation. + + For ``x**2 + y**2 <= 4`` (Q = I, PSD), the exact hull should: + - Introduce an auxiliary variable ``t >= 0``. + - Add a rotated SOC constraint ``v_x**2 + v_y**2 - t*y_ind <= 0``. + - Replace the original bound with the linear constraint + ``t - 4*y_ind <= 0``. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + v_y = relaxBlock.disaggregatedVars.y + y_ind = m.d1.binary_indicator_var + + # An auxiliary t variable should have been created. + t_var = relaxBlock.component('_conic_aux_t_c') + self.assertIsNotNone(t_var, "Expected auxiliary variable '_conic_aux_t_c'") + self.assertIs(t_var.domain, NonNegativeReals) + + # Two transformed constraints: conic SOC + linear bound. + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 2) + + conic_cons, linear_cons = trans_cons + + # --- Conic constraint: v_x**2 + v_y**2 - t*y_ind <= 0 --- + self.assertIsNone(conic_cons.lower) + self.assertEqual(value(conic_cons.upper), 0) + repn = generate_standard_repn(conic_cons.body, compute_values=False) + self.assertTrue(repn.is_quadratic()) + quad_pairs = dict( + zip( + [(id(a), id(b)) for a, b in repn.quadratic_vars], + repn.quadratic_coefs, + ) + ) + self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) + self.assertAlmostEqual(quad_pairs[(id(v_y), id(v_y))], 1) + self.assertAlmostEqual( + quad_pairs.get((id(t_var), id(y_ind)), quad_pairs.get((id(y_ind), id(t_var)), None)), + -1, + ) + self.assertEqual(repn.constant, 0) + self.assertEqual(len(repn.linear_vars), 0) + + # --- Linear bound constraint: t - 4*y_ind <= 0 --- + self.assertIsNone(linear_cons.lower) + self.assertEqual(value(linear_cons.upper), 0) + repn2 = generate_standard_repn(linear_cons.body, compute_values=False) + self.assertTrue(repn2.is_linear()) + linear_map = dict(zip([id(v) for v in repn2.linear_vars], repn2.linear_coefs)) + self.assertAlmostEqual(linear_map[id(t_var)], 1) + self.assertAlmostEqual(linear_map[id(y_ind)], -4) + + # ------------------------------------------------------------------ + # 2. Convex quadratic lower-bound → conic reformulation with negation + # ------------------------------------------------------------------ + def test_convex_quadratic_lower_bound_conic_reformulation(self): + """NSD Q with lower bound should produce a conic reformulation via negation. + + For ``-x**2 - y**2 >= -4`` (Q = -I, NSD), the exact hull negates Q + to obtain a PSD form and produces: + - Rotated SOC constraint ``v_x**2 + v_y**2 - t*y_ind <= 0``. + - Linear bound constraint ``t - 4*y_ind <= 0``. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=-m.x**2 - m.y**2 >= -4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + v_y = relaxBlock.disaggregatedVars.y + y_ind = m.d1.binary_indicator_var + + # Auxiliary t variable must exist. + t_var = relaxBlock.component('_conic_aux_t_c') + self.assertIsNotNone(t_var, "Expected auxiliary variable '_conic_aux_t_c'") + self.assertIs(t_var.domain, NonNegativeReals) + + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 2) + + conic_cons, linear_cons = trans_cons + + # --- Conic constraint (uses negated Q, i.e. +I): v_x**2 + v_y**2 - t*y_ind <= 0 --- + self.assertIsNone(conic_cons.lower) + self.assertEqual(value(conic_cons.upper), 0) + repn = generate_standard_repn(conic_cons.body, compute_values=False) + self.assertTrue(repn.is_quadratic()) + quad_pairs = dict( + zip( + [(id(a), id(b)) for a, b in repn.quadratic_vars], + repn.quadratic_coefs, + ) + ) + self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) + self.assertAlmostEqual(quad_pairs[(id(v_y), id(v_y))], 1) + + # --- Linear bound constraint: t - 4*y_ind <= 0 (from t <= -c.lower*y = 4y) --- + self.assertIsNone(linear_cons.lower) + self.assertEqual(value(linear_cons.upper), 0) + repn2 = generate_standard_repn(linear_cons.body, compute_values=False) + self.assertTrue(repn2.is_linear()) + linear_map = dict(zip([id(v) for v in repn2.linear_vars], repn2.linear_coefs)) + self.assertAlmostEqual(linear_map[id(t_var)], 1) + self.assertAlmostEqual(linear_map[id(y_ind)], -4) + + # ------------------------------------------------------------------ + # 3. Non-convex quadratic → general exact hull + # ------------------------------------------------------------------ + def test_nonconvex_quadratic_general_exact_hull(self): + """Mixed-eigenvalue Q should produce the general exact hull expression. + + For ``x**2 - y**2 <= 1`` (Q = diag(1,-1), indefinite), the exact + hull should produce a single constraint: + ``v_x**2 - v_y**2 - y_ind**2 <= 0`` + with no auxiliary variable or extra conic constraint. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 - m.y**2 <= 1) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + v_y = relaxBlock.disaggregatedVars.y + y_ind = m.d1.binary_indicator_var + + # No auxiliary t variable should exist. + self.assertIsNone(relaxBlock.component('_conic_aux_t_c')) + + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 1) + ub_cons = trans_cons[0] + + self.assertIsNone(ub_cons.lower) + self.assertEqual(value(ub_cons.upper), 0) + + repn = generate_standard_repn(ub_cons.body, compute_values=False) + self.assertTrue(repn.is_quadratic()) + self.assertEqual(repn.constant, 0) + self.assertEqual(len(repn.linear_vars), 0) + + quad_pairs = dict( + zip( + [(id(a), id(b)) for a, b in repn.quadratic_vars], + repn.quadratic_coefs, + ) + ) + # v_x**2 coefficient: +1 + self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) + # v_y**2 coefficient: -1 + self.assertAlmostEqual(quad_pairs[(id(v_y), id(v_y))], -1) + # y_ind**2 coefficient: -1 (from rhs: 1 * y_ind**2 moved to LHS) + self.assertAlmostEqual(quad_pairs[(id(y_ind), id(y_ind))], -1) + + # ------------------------------------------------------------------ + # 4. Equality constraint → always uses general exact hull + # ------------------------------------------------------------------ + def test_equality_constraint_general_exact_hull(self): + """Equality constraints should always use the general exact hull. + + Even with a PSD Q (``x**2 + y**2 == 4``), the exact hull must + use the general formulation: + ``v_x**2 + v_y**2 - 4*y_ind**2 == 0`` + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + m.y**2 == 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + v_y = relaxBlock.disaggregatedVars.y + y_ind = m.d1.binary_indicator_var + + # No auxiliary t variable should exist (no conic path for equality). + self.assertIsNone(relaxBlock.component('_conic_aux_t_c')) + + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 1) + eq_cons = trans_cons[0] + + # The constraint is an equality. + self.assertEqual(value(eq_cons.lower), 0) + self.assertEqual(value(eq_cons.upper), 0) + + repn = generate_standard_repn(eq_cons.body, compute_values=False) + self.assertTrue(repn.is_quadratic()) + self.assertEqual(repn.constant, 0) + self.assertEqual(len(repn.linear_vars), 0) + + quad_pairs = dict( + zip( + [(id(a), id(b)) for a, b in repn.quadratic_vars], + repn.quadratic_coefs, + ) + ) + self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) + self.assertAlmostEqual(quad_pairs[(id(v_y), id(v_y))], 1) + # rhs 4*y_ind**2 moves to LHS as -4*y_ind**2 + self.assertAlmostEqual(quad_pairs[(id(y_ind), id(y_ind))], -4) + + # ------------------------------------------------------------------ + # 5. Range constraint with both bounds → mixed reformulations + # ------------------------------------------------------------------ + def test_range_constraint_mixed_reformulations(self): + """Range constraint with PSD Q uses conic for upper, general for lower. + + For ``1 <= x**2 + y**2 <= 4`` (Q = I, PSD but not NSD): + - Upper bound uses the conic reformulation (``t - 4*y_ind <= 0``). + - Lower bound uses the general exact hull (``y_ind**2 - v_x**2 - v_y**2 <= 0``). + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=(1, m.x**2 + m.y**2, 4)) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + v_y = relaxBlock.disaggregatedVars.y + y_ind = m.d1.binary_indicator_var + + # Auxiliary t should exist (conic upper bound). + t_var = relaxBlock.component('_conic_aux_t_c') + self.assertIsNotNone(t_var) + + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + # The mapping returns: [conic_SOC, lb_general, ub_linear] + # (conic constraint first, then the bound constraints in order lb, ub). + self.assertEqual(len(trans_cons), 3) + + # Pre-compute all repns to avoid duplicate calls in the search loops. + repns = {id(c): generate_standard_repn(c.body, compute_values=False) + for c in trans_cons} + + # Identify each constraint by its standard representation. + # The linear one (involving t) is the conic upper bound. + # The quadratic one without t_var quad terms is the general lower bound. + ub_cons = next( + c for c in trans_cons + if repns[id(c)].is_linear() + ) + lb_cons = next( + c for c in trans_cons + if repns[id(c)].is_quadratic() + and id(t_var) not in { + id(v) + for pair in repns[id(c)].quadratic_vars + for v in pair + } + ) + + # --- Upper bound: linear conic bound ``t - 4*y_ind <= 0`` --- + self.assertIsNone(ub_cons.lower) + self.assertEqual(value(ub_cons.upper), 0) + repn_ub = repns[id(ub_cons)] + self.assertTrue(repn_ub.is_linear()) + linear_map = dict(zip([id(v) for v in repn_ub.linear_vars], repn_ub.linear_coefs)) + self.assertAlmostEqual(linear_map[id(t_var)], 1) + self.assertAlmostEqual(linear_map[id(y_ind)], -4) + + # --- Lower bound: general exact hull ``y_ind**2 - v_x**2 - v_y**2 <= 0`` --- + self.assertIsNone(lb_cons.lower) + self.assertEqual(value(lb_cons.upper), 0) + repn_lb = repns[id(lb_cons)] + self.assertTrue(repn_lb.is_quadratic()) + quad_map = dict( + zip( + [(id(a), id(b)) for a, b in repn_lb.quadratic_vars], + repn_lb.quadratic_coefs, + ) + ) + self.assertAlmostEqual(quad_map.get((id(v_x), id(v_x)), 0), -1) + self.assertAlmostEqual(quad_map.get((id(v_y), id(v_y)), 0), -1) + self.assertAlmostEqual(quad_map.get((id(y_ind), id(y_ind)), 0), 1) + + # ------------------------------------------------------------------ + # 6. numpy_available guard + # ------------------------------------------------------------------ + def test_numpy_unavailable_raises_error(self): + """When NumPy is unavailable, ``exact_hull_quadratic=True`` must raise. + + The transformation should raise a ``GDP_Error`` mentioning NumPy + before attempting any eigenvalue computation. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + with mock.patch.object(hull_module, 'numpy_available', False): + self.assertRaisesRegex( + GDP_Error, + "exact_hull_quadratic requires NumPy", + TransformationFactory('gdp.hull').apply_to, + m, + exact_hull_quadratic=True, + ) + + # ------------------------------------------------------------------ + # 7. eigenvalue_tolerance: larger tolerance accepts a near-PSD matrix + # ------------------------------------------------------------------ + def test_eigenvalue_tolerance_permissive(self): + """A larger ``eigenvalue_tolerance`` should accept a near-PSD Q. + + The matrix Q = diag(1, -1e-11) has a very small negative eigenvalue + (-1e-11). With the default tolerance of 1e-10 the eigenvalue is + within the tolerance band and Q is treated as PSD → conic + reformulation. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 - 1e-11 * m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + m_perm = self.hull.create_using( + m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10 + ) + relaxBlock = m_perm._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + # Conic path should have been taken → t variable present. + self.assertIsNotNone( + relaxBlock.component('_conic_aux_t_c'), + "Expected conic aux variable with permissive eigenvalue_tolerance", + ) + + # ------------------------------------------------------------------ + # 8. eigenvalue_tolerance: strict tolerance rejects a near-PSD matrix + # ------------------------------------------------------------------ + def test_eigenvalue_tolerance_strict(self): + """A strict ``eigenvalue_tolerance`` should reject a near-PSD Q. + + Same Q = diag(1, -1e-11) as above. With a tolerance of 1e-12 the + eigenvalue -1e-11 is outside the tolerance band → Q is not treated + as PSD → general exact hull (no t variable). + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 - 1e-11 * m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + m_strict = self.hull.create_using( + m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-12 + ) + relaxBlock = m_strict._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + # General path should have been taken → no t variable. + self.assertIsNone( + relaxBlock.component('_conic_aux_t_c'), + "Expected no conic aux variable with strict eigenvalue_tolerance", + ) + + # ------------------------------------------------------------------ + # 9. Constraint mapping: get_transformed_constraints + # ------------------------------------------------------------------ + def test_constraint_mappings_preserved(self): + """``get_transformed_constraints`` must return all constraints for a + conic-reformulated quadratic constraint. + + For a PSD upper-bound constraint the return list should contain both + the conic SOC constraint and the linear bound constraint. + ``get_src_constraint`` must map back to the original constraint. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 2) + + for tc in trans_cons: + self.assertIs(self.hull.get_src_constraint(tc), m.d1.c) + + # ------------------------------------------------------------------ + # 10. Linear and constant terms in conic reformulation + # ------------------------------------------------------------------ + def test_conic_reformulation_with_linear_and_constant_terms(self): + """Conic reformulation correctly incorporates linear/constant terms. + + For ``x**2 + 3*x + 2 <= 0`` (PSD, with a linear and a constant term) + the linear bound constraint must be: + ``t + 3*v_x + 2*y_ind <= 0`` + (where the RHS of 0 is absorbed because upper=0). + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + 3 * m.x + 2 <= 0) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + y_ind = m.d1.binary_indicator_var + t_var = relaxBlock.component('_conic_aux_t_c') + self.assertIsNotNone(t_var) + + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 2) + + # The linear bound constraint is the one involving t. + linear_cons = next( + c for c in trans_cons + if generate_standard_repn(c.body, compute_values=False).is_linear() + ) + repn = generate_standard_repn(linear_cons.body, compute_values=False) + self.assertTrue(repn.is_linear()) + linear_map = dict(zip([id(v) for v in repn.linear_vars], repn.linear_coefs)) + + # t coefficient: 1 + self.assertAlmostEqual(linear_map[id(t_var)], 1) + # 3*v_x term + self.assertAlmostEqual(linear_map[id(v_x)], 3) + # 2*y_ind constant term (constant=2 * y) + self.assertAlmostEqual(linear_map[id(y_ind)], 2) + + # ------------------------------------------------------------------ + # 11. Linear and constant terms in general exact hull + # ------------------------------------------------------------------ + def test_general_exact_hull_with_linear_and_constant_terms(self): + """General exact hull correctly incorporates linear/constant terms. + + For ``x**2 - y**2 + 3*x - 2 <= 0`` (indefinite Q, with linear/const + terms) the single transformed constraint should be: + ``v_x**2 - v_y**2 + 3*v_x*y_ind - 2*y_ind**2 <= 0`` + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 - m.y**2 + 3 * m.x - 2 <= 0) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + + self.hull.apply_to(m, exact_hull_quadratic=True) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + v_y = relaxBlock.disaggregatedVars.y + y_ind = m.d1.binary_indicator_var + + # No conic aux variable (indefinite Q → general path). + self.assertIsNone(relaxBlock.component('_conic_aux_t_c')) + + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 1) + ub_cons = trans_cons[0] + + self.assertIsNone(ub_cons.lower) + self.assertEqual(value(ub_cons.upper), 0) + + repn = generate_standard_repn(ub_cons.body, compute_values=False) + self.assertTrue(repn.is_quadratic()) + self.assertEqual(repn.constant, 0) + self.assertEqual(len(repn.linear_vars), 0) + + quad_pairs = dict( + zip( + [(id(a), id(b)) for a, b in repn.quadratic_vars], + repn.quadratic_coefs, + ) + ) + # v_x**2: +1 + self.assertAlmostEqual(quad_pairs.get((id(v_x), id(v_x)), 0), 1) + # v_y**2: -1 + self.assertAlmostEqual(quad_pairs.get((id(v_y), id(v_y)), 0), -1) + # 3*v_x*y_ind cross term + cross_coef = quad_pairs.get( + (id(v_x), id(y_ind)), + quad_pairs.get((id(y_ind), id(v_x)), None), + ) + self.assertIsNotNone(cross_coef, "Missing cross term v_x*y_ind") + self.assertAlmostEqual(cross_coef, 3) + # -2*y_ind**2 (constant term * y^2) + self.assertAlmostEqual(quad_pairs.get((id(y_ind), id(y_ind)), 0), -2) From e008716e003ab4ea6ed57547428bdb2f556aab21 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 05:43:48 +0000 Subject: [PATCH 15/39] Move TestExactHullQuadratic inline models into models.py as factory functions Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/tests/models.py | 148 +++++++++++++++++++++++++++++++++++ pyomo/gdp/tests/test_hull.py | 99 +++-------------------- 2 files changed, 159 insertions(+), 88 deletions(-) diff --git a/pyomo/gdp/tests/models.py b/pyomo/gdp/tests/models.py index dde0b007c31..6ad75496a4f 100644 --- a/pyomo/gdp/tests/models.py +++ b/pyomo/gdp/tests/models.py @@ -1262,3 +1262,151 @@ def disj_rule(m, t): m.obj = Objective(expr=m.x[1] + m.x[2], sense=minimize) return m + + +# --------------------------------------------------------------------------- +# Models for exact hull quadratic tests +# --------------------------------------------------------------------------- + + +def makeTwoTermDisj_ConvexQuadUB(): + """Two-term disjunction with a convex quadratic upper-bound constraint. + + Disjunct d1 has ``x**2 + y**2 <= 4`` (Q = I, positive semi-definite), + which should be handled by the conic exact hull reformulation. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_ConvexQuadLB(): + """Two-term disjunction with a convex quadratic lower-bound constraint. + + Disjunct d1 has ``-x**2 - y**2 >= -4`` (Q = -I, negative semi-definite), + which should be handled by the conic exact hull reformulation via negation. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=-m.x**2 - m.y**2 >= -4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_NonconvexQuad(): + """Two-term disjunction with a non-convex (indefinite) quadratic constraint. + + Disjunct d1 has ``x**2 - y**2 <= 1`` (Q = diag(1, -1), indefinite), + which should be handled by the general exact hull reformulation. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 - m.y**2 <= 1) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_QuadEquality(): + """Two-term disjunction with a quadratic equality constraint. + + Disjunct d1 has ``x**2 + y**2 == 4`` (Q = I, PSD), which should always + be handled by the general exact hull because it is an equality. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + m.y**2 == 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_QuadRange(): + """Two-term disjunction with a quadratic range (double-bounded) constraint. + + Disjunct d1 has ``1 <= x**2 + y**2 <= 4`` (Q = I, PSD but not NSD). + The upper bound should use the conic reformulation and the lower bound + should use the general exact hull. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=(1, m.x**2 + m.y**2, 4)) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_NearPSDQuad(): + """Two-term disjunction with a near-PSD quadratic upper-bound constraint. + + Disjunct d1 has ``x**2 - 1e-11*y**2 <= 4``. The Q matrix + ``diag(1, -1e-11)`` has a very small negative eigenvalue (-1e-11), so + its classification as PSD or indefinite depends on the + ``eigenvalue_tolerance`` setting. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 - 1e-11 * m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_ConvexQuad_LinearTerms(): + """Two-term disjunction with a convex quadratic constraint that has linear + and constant terms. + + Disjunct d1 has ``x**2 + 3*x + 2 <= 0`` (Q = [[1]], PSD). Tests that + linear coefficients and the constant term are correctly incorporated into + the conic exact hull reformulation. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 + 3 * m.x + 2 <= 0) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_NonconvexQuad_LinearTerms(): + """Two-term disjunction with an indefinite quadratic constraint that has + linear and constant terms. + + Disjunct d1 has ``x**2 - y**2 + 3*x - 2 <= 0`` (Q = diag(1, -1), + indefinite). Tests that linear and constant terms are correctly + incorporated into the general exact hull reformulation. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.x**2 - m.y**2 + 3 * m.x - 2 <= 0) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 3e60239a082..b444de3385c 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2926,14 +2926,7 @@ def test_convex_quadratic_upper_bound_conic_reformulation(self): - Replace the original bound with the linear constraint ``t - 4*y_ind <= 0``. """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 + m.y**2 <= 4) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_ConvexQuadUB() self.hull.apply_to(m, exact_hull_quadratic=True) @@ -2993,14 +2986,7 @@ def test_convex_quadratic_lower_bound_conic_reformulation(self): - Rotated SOC constraint ``v_x**2 + v_y**2 - t*y_ind <= 0``. - Linear bound constraint ``t - 4*y_ind <= 0``. """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=-m.x**2 - m.y**2 >= -4) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_ConvexQuadLB() self.hull.apply_to(m, exact_hull_quadratic=True) @@ -3053,14 +3039,7 @@ def test_nonconvex_quadratic_general_exact_hull(self): ``v_x**2 - v_y**2 - y_ind**2 <= 0`` with no auxiliary variable or extra conic constraint. """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 - m.y**2 <= 1) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_NonconvexQuad() self.hull.apply_to(m, exact_hull_quadratic=True) @@ -3107,14 +3086,7 @@ def test_equality_constraint_general_exact_hull(self): use the general formulation: ``v_x**2 + v_y**2 - 4*y_ind**2 == 0`` """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 + m.y**2 == 4) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_QuadEquality() self.hull.apply_to(m, exact_hull_quadratic=True) @@ -3160,14 +3132,7 @@ def test_range_constraint_mixed_reformulations(self): - Upper bound uses the conic reformulation (``t - 4*y_ind <= 0``). - Lower bound uses the general exact hull (``y_ind**2 - v_x**2 - v_y**2 <= 0``). """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=(1, m.x**2 + m.y**2, 4)) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_QuadRange() self.hull.apply_to(m, exact_hull_quadratic=True) @@ -3239,14 +3204,7 @@ def test_numpy_unavailable_raises_error(self): The transformation should raise a ``GDP_Error`` mentioning NumPy before attempting any eigenvalue computation. """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 + m.y**2 <= 4) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_ConvexQuadUB() with mock.patch.object(hull_module, 'numpy_available', False): self.assertRaisesRegex( @@ -3268,14 +3226,7 @@ def test_eigenvalue_tolerance_permissive(self): within the tolerance band and Q is treated as PSD → conic reformulation. """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 - 1e-11 * m.y**2 <= 4) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_NearPSDQuad() m_perm = self.hull.create_using( m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10 @@ -3297,14 +3248,7 @@ def test_eigenvalue_tolerance_strict(self): eigenvalue -1e-11 is outside the tolerance band → Q is not treated as PSD → general exact hull (no t variable). """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 - 1e-11 * m.y**2 <= 4) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_NearPSDQuad() m_strict = self.hull.create_using( m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-12 @@ -3327,14 +3271,7 @@ def test_constraint_mappings_preserved(self): the conic SOC constraint and the linear bound constraint. ``get_src_constraint`` must map back to the original constraint. """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 + m.y**2 <= 4) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_ConvexQuadUB() self.hull.apply_to(m, exact_hull_quadratic=True) @@ -3355,14 +3292,7 @@ def test_conic_reformulation_with_linear_and_constant_terms(self): ``t + 3*v_x + 2*y_ind <= 0`` (where the RHS of 0 is absorbed because upper=0). """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 + 3 * m.x + 2 <= 0) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_ConvexQuad_LinearTerms() self.hull.apply_to(m, exact_hull_quadratic=True) @@ -3401,14 +3331,7 @@ def test_general_exact_hull_with_linear_and_constant_terms(self): terms) the single transformed constraint should be: ``v_x**2 - v_y**2 + 3*v_x*y_ind - 2*y_ind**2 <= 0`` """ - m = ConcreteModel() - m.x = Var(bounds=(-2, 2)) - m.y = Var(bounds=(-2, 2)) - m.d1 = Disjunct() - m.d1.c = Constraint(expr=m.x**2 - m.y**2 + 3 * m.x - 2 <= 0) - m.d2 = Disjunct() - m.d2.c = Constraint(expr=m.x + m.y <= 1) - m.disj = Disjunction(expr=[m.d1, m.d2]) + m = models.makeTwoTermDisj_NonconvexQuad_LinearTerms() self.hull.apply_to(m, exact_hull_quadratic=True) From 40c0077197bd91dd60f712f862b60022d4ee17f5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 06:08:13 +0000 Subject: [PATCH 16/39] Apply review feedback: per-test skipUnless decorators and updated eigenvalue_tolerance in permissive test Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/tests/test_hull.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index b444de3385c..cc3c6bdd914 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2895,7 +2895,6 @@ def test_declare_disjuncts_in_disjunction_rule(self): ct.check_nested_disjuncts_in_flat_gdp(self, 'hull') -@unittest.skipUnless(numpy_available, "NumPy is not available") class TestExactHullQuadratic(unittest.TestCase): """Tests for the ``exact_hull_quadratic`` option of the hull transformation. @@ -2917,6 +2916,7 @@ def setUp(self): # ------------------------------------------------------------------ # 1. Convex quadratic upper-bound → conic reformulation # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_convex_quadratic_upper_bound_conic_reformulation(self): """PSD Q with upper bound should produce a conic reformulation. @@ -2978,6 +2978,7 @@ def test_convex_quadratic_upper_bound_conic_reformulation(self): # ------------------------------------------------------------------ # 2. Convex quadratic lower-bound → conic reformulation with negation # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_convex_quadratic_lower_bound_conic_reformulation(self): """NSD Q with lower bound should produce a conic reformulation via negation. @@ -3031,6 +3032,7 @@ def test_convex_quadratic_lower_bound_conic_reformulation(self): # ------------------------------------------------------------------ # 3. Non-convex quadratic → general exact hull # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_nonconvex_quadratic_general_exact_hull(self): """Mixed-eigenvalue Q should produce the general exact hull expression. @@ -3079,6 +3081,7 @@ def test_nonconvex_quadratic_general_exact_hull(self): # ------------------------------------------------------------------ # 4. Equality constraint → always uses general exact hull # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_equality_constraint_general_exact_hull(self): """Equality constraints should always use the general exact hull. @@ -3125,6 +3128,7 @@ def test_equality_constraint_general_exact_hull(self): # ------------------------------------------------------------------ # 5. Range constraint with both bounds → mixed reformulations # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_range_constraint_mixed_reformulations(self): """Range constraint with PSD Q uses conic for upper, general for lower. @@ -3218,18 +3222,19 @@ def test_numpy_unavailable_raises_error(self): # ------------------------------------------------------------------ # 7. eigenvalue_tolerance: larger tolerance accepts a near-PSD matrix # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_eigenvalue_tolerance_permissive(self): """A larger ``eigenvalue_tolerance`` should accept a near-PSD Q. The matrix Q = diag(1, -1e-11) has a very small negative eigenvalue - (-1e-11). With the default tolerance of 1e-10 the eigenvalue is + (-1e-11). With a permissive tolerance of 1e-9 the eigenvalue is within the tolerance band and Q is treated as PSD → conic reformulation. """ m = models.makeTwoTermDisj_NearPSDQuad() m_perm = self.hull.create_using( - m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10 + m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-9 ) relaxBlock = m_perm._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] # Conic path should have been taken → t variable present. @@ -3241,6 +3246,7 @@ def test_eigenvalue_tolerance_permissive(self): # ------------------------------------------------------------------ # 8. eigenvalue_tolerance: strict tolerance rejects a near-PSD matrix # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_eigenvalue_tolerance_strict(self): """A strict ``eigenvalue_tolerance`` should reject a near-PSD Q. @@ -3263,6 +3269,7 @@ def test_eigenvalue_tolerance_strict(self): # ------------------------------------------------------------------ # 9. Constraint mapping: get_transformed_constraints # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_constraint_mappings_preserved(self): """``get_transformed_constraints`` must return all constraints for a conic-reformulated quadratic constraint. @@ -3284,6 +3291,7 @@ def test_constraint_mappings_preserved(self): # ------------------------------------------------------------------ # 10. Linear and constant terms in conic reformulation # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_conic_reformulation_with_linear_and_constant_terms(self): """Conic reformulation correctly incorporates linear/constant terms. @@ -3324,6 +3332,7 @@ def test_conic_reformulation_with_linear_and_constant_terms(self): # ------------------------------------------------------------------ # 11. Linear and constant terms in general exact hull # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") def test_general_exact_hull_with_linear_and_constant_terms(self): """General exact hull correctly incorporates linear/constant terms. From ce7ed66e7f0254a0e874ba5e6d58211d48dc7593 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 06:28:06 +0000 Subject: [PATCH 17/39] Initial plan From 928b3e19f1981673053959953d65921a726e250a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 06:34:45 +0000 Subject: [PATCH 18/39] Fix all-zero-eigenvalue bug: warn and fall back to general hull path Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 30 ++++++-- pyomo/gdp/tests/models.py | 42 +++++++++++ pyomo/gdp/tests/test_hull.py | 131 +++++++++++++++++++++++++++++++++++ 3 files changed, 196 insertions(+), 7 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index adf5f4f81db..1957d4ee45f 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1030,13 +1030,29 @@ def _build_exact_quadratic_hull( use_conic_lower = False negate_for_conic = False - if c.upper is not None and not c.equality: - if Q_is_psd: - use_conic_upper = True - if c.lower is not None and not c.equality: - if Q_is_nsd: - use_conic_lower = True - negate_for_conic = True + if Q_is_psd and Q_is_nsd: + # All eigenvalues lie within the tolerance band [-tol, tol], + # so the quadratic part is numerically zero and the constraint + # is treated as linear. Warn the user in case this is unexpected. + max_abs_eigenvalue = float(np.max(np.abs(eigenvalues))) + logger.warning( + "GDP(Hull): Constraint '%s' has quadratic terms, but all " + "eigenvalues of the Q matrix are within the " + "eigenvalue_tolerance band (largest eigenvalue by modulus: " + "%g). The constraint will be treated as linear. If this is " + "not the expected behavior, consider using a tighter (smaller) " + "eigenvalue_tolerance.", + c.getname(fully_qualified=True), + max_abs_eigenvalue, + ) + else: + if c.upper is not None and not c.equality: + if Q_is_psd: + use_conic_upper = True + if c.lower is not None and not c.equality: + if Q_is_nsd: + use_conic_lower = True + negate_for_conic = True # --- Decide which expression forms are needed --- need_non_convex = False diff --git a/pyomo/gdp/tests/models.py b/pyomo/gdp/tests/models.py index 6ad75496a4f..8b9e33b8442 100644 --- a/pyomo/gdp/tests/models.py +++ b/pyomo/gdp/tests/models.py @@ -1410,3 +1410,45 @@ def makeTwoTermDisj_NonconvexQuad_LinearTerms(): m.d2.c = Constraint(expr=m.x + m.y <= 1) m.disj = Disjunction(expr=[m.d1, m.d2]) return m + + +def makeTwoTermDisj_AllZeroEigenvalueQuad(): + """Two-term disjunction with a quadratic constraint whose Q matrix has + all eigenvalues within the default eigenvalue_tolerance band. + + Disjunct d1 has ``1e-11*x**2 + 1e-11*y**2 <= 4``. The Q matrix + ``diag(1e-11, 1e-11)`` has all eigenvalues equal to 1e-11, which is + within the default tolerance of 1e-10, so Q is both PSD and NSD. + The transformation should issue a warning and use the general exact hull + (no conic path) for this constraint. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=1e-11 * m.x**2 + 1e-11 * m.y**2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_AllZeroEigenvalueQuadRange(): + """Two-term disjunction with a quadratic range constraint whose Q matrix + has all eigenvalues within the default eigenvalue_tolerance band. + + Disjunct d1 has ``(1, 1e-11*x**2 + 1e-11*y**2, 4)``. The Q matrix + ``diag(1e-11, 1e-11)`` is both PSD and NSD under the default tolerance. + This exercises the two-sided bound bug: without the fix, a single + ``conic_expr_linear`` would be built with ``negate_for_conic=True`` + and then incorrectly reused for the upper bound. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=(1, 1e-11 * m.x**2 + 1e-11 * m.y**2, 4)) + m.d2 = Disjunct() + m.d2.c = Constraint(expr=m.x + m.y <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index cc3c6bdd914..b1c17520ca5 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -3383,3 +3383,134 @@ def test_general_exact_hull_with_linear_and_constant_terms(self): self.assertAlmostEqual(cross_coef, 3) # -2*y_ind**2 (constant term * y^2) self.assertAlmostEqual(quad_pairs.get((id(y_ind), id(y_ind)), 0), -2) + + # ------------------------------------------------------------------ + # 12. All-zero eigenvalues: warning issued, treated as linear + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_all_zero_eigenvalues_issues_warning(self): + """When all Q eigenvalues are within the tolerance band, a warning + should be issued and the general exact hull path should be taken. + + For ``1e-11*x**2 + 1e-11*y**2 <= 4`` with the default tolerance of + 1e-10, every eigenvalue of Q = diag(1e-11, 1e-11) lies in the band + [-1e-10, 1e-10]. The transformation must: + - Emit a warning mentioning the constraint name, that it is treated + as linear, the largest eigenvalue by modulus, and the suggestion + to use a tighter tolerance. + - Fall back to the general exact hull (no conic auxiliary variable). + """ + from pyomo.common.log import LoggingIntercept + from io import StringIO + + m = models.makeTwoTermDisj_AllZeroEigenvalueQuad() + + output = StringIO() + with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): + self.hull.apply_to(m, exact_hull_quadratic=True) + + warning_text = output.getvalue() + self.assertIn("quadratic terms", warning_text) + self.assertIn("eigenvalues", warning_text) + self.assertIn("treated as linear", warning_text) + self.assertIn("eigenvalue_tolerance", warning_text) + # The largest eigenvalue by modulus (1e-11) should appear in the message. + # Accept any floating-point rendering of 1e-11 (e.g. '1e-11', '1.0e-11'). + import re + self.assertTrue( + re.search(r'1[.,]?0*e-0*11', warning_text), + f"Expected the largest eigenvalue (1e-11) in the warning message, " + f"got: {warning_text!r}", + ) + + # General path should have been taken: no conic auxiliary variable. + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + self.assertIsNone( + relaxBlock.component('_conic_aux_t_c'), + "Expected no conic aux variable when all eigenvalues are near zero", + ) + + # There should be exactly one transformed constraint (general hull). + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 1) + + # ------------------------------------------------------------------ + # 13. All-zero eigenvalues with range constraint: no incorrect negation + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_all_zero_eigenvalues_range_constraint_no_negation_bug(self): + """Range constraint with all-zero Q eigenvalues must not produce + an incorrectly negated upper-bound reformulation. + + Without the fix, when Q is both PSD and NSD a single + ``conic_expr_linear`` is built with ``negate_for_conic=True`` + (from the lower-bound branch) and then reused for the upper bound, + producing sign-flipped coefficients. + + After the fix the general exact hull path is taken for both bounds, + resulting in two quadratic transformed constraints with correct signs. + """ + from pyomo.common.log import LoggingIntercept + from io import StringIO + + m = models.makeTwoTermDisj_AllZeroEigenvalueQuadRange() + + output = StringIO() + with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): + self.hull.apply_to(m, exact_hull_quadratic=True) + + # Warning must be issued. + self.assertIn("treated as linear", output.getvalue()) + + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + v_x = relaxBlock.disaggregatedVars.x + v_y = relaxBlock.disaggregatedVars.y + y_ind = m.d1.binary_indicator_var + + # No conic auxiliary variable. + self.assertIsNone(relaxBlock.component('_conic_aux_t_c')) + + # Both bounds should produce general hull constraints. + trans_cons = self.hull.get_transformed_constraints(m.d1.c) + self.assertEqual(len(trans_cons), 2) + + repns = [ + generate_standard_repn(tc.body, compute_values=False) + for tc in trans_cons + ] + + # Both should be quadratic (general hull includes the tiny quad terms). + for repn in repns: + self.assertTrue(repn.is_quadratic()) + + # Collect all v_x**2 coefficients across both transformed constraints. + # The upper-bound constraint body is + # ``1e-11*v_x**2 + 1e-11*v_y**2 - 4*y_ind**2 <= 0`` + # (coefficient of v_x**2 is +1e-11), while the lower-bound body is + # ``1*y_ind**2 - 1e-11*v_x**2 - 1e-11*v_y**2 <= 0`` + # (coefficient of v_x**2 is -1e-11). + # Without the fix, only a single conic expression with + # ``negate_for_conic=True`` would be built and reused, producing two + # negated constraints instead of one positive and one negative. + # We therefore assert that at least one constraint has a positive + # v_x**2 coefficient and at least one has a negative v_x**2 coefficient. + coefs_vx = [] + for repn in repns: + quad_map = dict( + zip( + [(id(a), id(b)) for a, b in repn.quadratic_vars], + repn.quadratic_coefs, + ) + ) + coefs_vx.append(quad_map.get((id(v_x), id(v_x)), 0)) + + self.assertTrue( + any(c > 0 for c in coefs_vx), + f"Expected at least one constraint with positive v_x**2 coefficient " + f"(upper-bound general hull), got: {coefs_vx}", + ) + self.assertTrue( + any(c < 0 for c in coefs_vx), + f"Expected at least one constraint with negative v_x**2 coefficient " + f"(lower-bound general hull), got: {coefs_vx}", + ) From 1981a70b2905f64559c54611ac3e9f21cd0477b5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 06:45:06 +0000 Subject: [PATCH 19/39] Fix warning message: says general exact hull, not linear Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 15 +++++++++------ pyomo/gdp/tests/test_hull.py | 4 ++-- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 1957d4ee45f..ff8e5d4bfeb 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1031,17 +1031,20 @@ def _build_exact_quadratic_hull( negate_for_conic = False if Q_is_psd and Q_is_nsd: - # All eigenvalues lie within the tolerance band [-tol, tol], - # so the quadratic part is numerically zero and the constraint - # is treated as linear. Warn the user in case this is unexpected. + # All eigenvalues lie within the tolerance band [-tol, tol]. + # The conic path cannot be used for either bound because PSD and + # NSD simultaneously would require a single negated expression to + # serve both bounds (causing sign-flip bugs). Fall back to the + # general exact hull reformulation instead. max_abs_eigenvalue = float(np.max(np.abs(eigenvalues))) logger.warning( "GDP(Hull): Constraint '%s' has quadratic terms, but all " "eigenvalues of the Q matrix are within the " "eigenvalue_tolerance band (largest eigenvalue by modulus: " - "%g). The constraint will be treated as linear. If this is " - "not the expected behavior, consider using a tighter (smaller) " - "eigenvalue_tolerance.", + "%g). The conic reformulation cannot be applied; the " + "constraint will be handled by the general exact hull " + "reformulation instead. If this is not the expected behavior, " + "consider using a tighter (smaller) eigenvalue_tolerance.", c.getname(fully_qualified=True), max_abs_eigenvalue, ) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index b1c17520ca5..47ea8eed598 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -3412,7 +3412,7 @@ def test_all_zero_eigenvalues_issues_warning(self): warning_text = output.getvalue() self.assertIn("quadratic terms", warning_text) self.assertIn("eigenvalues", warning_text) - self.assertIn("treated as linear", warning_text) + self.assertIn("general exact hull", warning_text) self.assertIn("eigenvalue_tolerance", warning_text) # The largest eigenvalue by modulus (1e-11) should appear in the message. # Accept any floating-point rendering of 1e-11 (e.g. '1e-11', '1.0e-11'). @@ -3460,7 +3460,7 @@ def test_all_zero_eigenvalues_range_constraint_no_negation_bug(self): self.hull.apply_to(m, exact_hull_quadratic=True) # Warning must be issued. - self.assertIn("treated as linear", output.getvalue()) + self.assertIn("general exact hull", output.getvalue()) relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] v_x = relaxBlock.disaggregatedVars.x From 97caf8932db910c315acfb61cc781254e5a38606 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 06:51:23 +0000 Subject: [PATCH 20/39] Remove last "treated as linear" references from comment and docstring Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/tests/test_hull.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 47ea8eed598..b71f80ab549 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -3385,7 +3385,7 @@ def test_general_exact_hull_with_linear_and_constant_terms(self): self.assertAlmostEqual(quad_pairs.get((id(y_ind), id(y_ind)), 0), -2) # ------------------------------------------------------------------ - # 12. All-zero eigenvalues: warning issued, treated as linear + # 12. All-zero eigenvalues: warning issued, falls back to general hull # ------------------------------------------------------------------ @unittest.skipUnless(numpy_available, "NumPy is not available") def test_all_zero_eigenvalues_issues_warning(self): @@ -3395,9 +3395,9 @@ def test_all_zero_eigenvalues_issues_warning(self): For ``1e-11*x**2 + 1e-11*y**2 <= 4`` with the default tolerance of 1e-10, every eigenvalue of Q = diag(1e-11, 1e-11) lies in the band [-1e-10, 1e-10]. The transformation must: - - Emit a warning mentioning the constraint name, that it is treated - as linear, the largest eigenvalue by modulus, and the suggestion - to use a tighter tolerance. + - Emit a warning mentioning the constraint name, that the conic + reformulation cannot be applied, the largest eigenvalue by modulus, + and the suggestion to use a tighter tolerance. - Fall back to the general exact hull (no conic auxiliary variable). """ from pyomo.common.log import LoggingIntercept From 5032059963108cc8c280fad4e085fda359245a7f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 06:55:55 +0000 Subject: [PATCH 21/39] Move inline imports to module level in test_hull.py Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/tests/test_hull.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index b71f80ab549..10c284443b3 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -8,6 +8,7 @@ # ____________________________________________________________________________________ import logging +import re import sys import random from io import StringIO @@ -3400,9 +3401,6 @@ def test_all_zero_eigenvalues_issues_warning(self): and the suggestion to use a tighter tolerance. - Fall back to the general exact hull (no conic auxiliary variable). """ - from pyomo.common.log import LoggingIntercept - from io import StringIO - m = models.makeTwoTermDisj_AllZeroEigenvalueQuad() output = StringIO() @@ -3416,7 +3414,6 @@ def test_all_zero_eigenvalues_issues_warning(self): self.assertIn("eigenvalue_tolerance", warning_text) # The largest eigenvalue by modulus (1e-11) should appear in the message. # Accept any floating-point rendering of 1e-11 (e.g. '1e-11', '1.0e-11'). - import re self.assertTrue( re.search(r'1[.,]?0*e-0*11', warning_text), f"Expected the largest eigenvalue (1e-11) in the warning message, " @@ -3450,9 +3447,6 @@ def test_all_zero_eigenvalues_range_constraint_no_negation_bug(self): After the fix the general exact hull path is taken for both bounds, resulting in two quadratic transformed constraints with correct signs. """ - from pyomo.common.log import LoggingIntercept - from io import StringIO - m = models.makeTwoTermDisj_AllZeroEigenvalueQuadRange() output = StringIO() From c512f6e06facd1249b46a024b2c8d00e0127ba7e Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Sun, 8 Mar 2026 03:09:29 -0400 Subject: [PATCH 22/39] Update pyomo/gdp/tests/test_hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyomo/gdp/tests/test_hull.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 10c284443b3..1ae29216406 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -3405,7 +3405,9 @@ def test_all_zero_eigenvalues_issues_warning(self): output = StringIO() with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): - self.hull.apply_to(m, exact_hull_quadratic=True) + self.hull.apply_to( + m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10 + ) warning_text = output.getvalue() self.assertIn("quadratic terms", warning_text) From 68c350083ba7bc512ff40a25dc18586df0402230 Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Sun, 8 Mar 2026 03:09:45 -0400 Subject: [PATCH 23/39] Update pyomo/gdp/tests/test_hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyomo/gdp/tests/test_hull.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 1ae29216406..629624fdc9f 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -3453,7 +3453,9 @@ def test_all_zero_eigenvalues_range_constraint_no_negation_bug(self): output = StringIO() with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): - self.hull.apply_to(m, exact_hull_quadratic=True) + self.hull.apply_to( + m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10 + ) # Warning must be issued. self.assertIn("general exact hull", output.getvalue()) From 0adfba778c4c66c70f12c6ff2952cb95730a9351 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 07:10:32 +0000 Subject: [PATCH 24/39] Gate all-zero-eigenvalue fallback warning on not c.equality Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index ff8e5d4bfeb..26a8f32d7e8 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1036,18 +1036,21 @@ def _build_exact_quadratic_hull( # NSD simultaneously would require a single negated expression to # serve both bounds (causing sign-flip bugs). Fall back to the # general exact hull reformulation instead. - max_abs_eigenvalue = float(np.max(np.abs(eigenvalues))) - logger.warning( - "GDP(Hull): Constraint '%s' has quadratic terms, but all " - "eigenvalues of the Q matrix are within the " - "eigenvalue_tolerance band (largest eigenvalue by modulus: " - "%g). The conic reformulation cannot be applied; the " - "constraint will be handled by the general exact hull " - "reformulation instead. If this is not the expected behavior, " - "consider using a tighter (smaller) eigenvalue_tolerance.", - c.getname(fully_qualified=True), - max_abs_eigenvalue, - ) + # Only warn for inequality constraints; equality constraints always + # use the general exact hull path, so no fallback warning is needed. + if not c.equality: + max_abs_eigenvalue = float(np.max(np.abs(eigenvalues))) + logger.warning( + "GDP(Hull): Constraint '%s' has quadratic terms, but all " + "eigenvalues of the Q matrix are within the " + "eigenvalue_tolerance band (largest eigenvalue by modulus: " + "%g). The conic reformulation cannot be applied; the " + "constraint will be handled by the general exact hull " + "reformulation instead. If this is not the expected behavior, " + "consider using a tighter (smaller) eigenvalue_tolerance.", + c.getname(fully_qualified=True), + max_abs_eigenvalue, + ) else: if c.upper is not None and not c.equality: if Q_is_psd: From 55cee770307de341c17273fcf9838f5a878e1f3d Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Sun, 8 Mar 2026 03:17:19 -0400 Subject: [PATCH 25/39] Update pyomo/gdp/plugins/hull.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 26a8f32d7e8..80e73edf1fb 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1032,10 +1032,12 @@ def _build_exact_quadratic_hull( if Q_is_psd and Q_is_nsd: # All eigenvalues lie within the tolerance band [-tol, tol]. - # The conic path cannot be used for either bound because PSD and - # NSD simultaneously would require a single negated expression to - # serve both bounds (causing sign-flip bugs). Fall back to the - # general exact hull reformulation instead. + # In this numerically ambiguous case, we conservatively avoid the + # conic reformulation and fall back to the general exact hull + # reformulation. The underlying sign-flip bug motivating this + # choice occurs for two-sided (range) constraints, where a single + # conic expression would otherwise be reused (possibly negated) + # for both bounds. # Only warn for inequality constraints; equality constraints always # use the general exact hull path, so no fallback warning is needed. if not c.equality: From 52cf0a76b2866cf66b19d91880aa699bd23787f8 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 07:32:03 +0000 Subject: [PATCH 26/39] Initial plan From 87f07d3d9261d98f22a1fae8f89715bddbc0d779 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 8 Mar 2026 07:33:53 +0000 Subject: [PATCH 27/39] Reorder CONFIG: declare exact_hull_quadratic before eigenvalue_tolerance Co-authored-by: sergey-gusev94 <101810399+sergey-gusev94@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 66 +++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 80e73edf1fb..7676c7f3334 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -204,6 +204,39 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): description="Epsilon value to use in perspective function", ), ) + CONFIG.declare( + 'exact_hull_quadratic', + cfg.ConfigValue( + default=False, + domain=bool, + description="Use exact hull reformulation for quadratic constraints.", + doc=""" + If True, quadratic constraints (polynomial degree 2) are reformulated + using the exact hull instead of the standard perspective function. + + For a quadratic constraint of the form + + x'Qx + c'x + d <= 0 + + the reformulation depends on convexity: + + **Conic exact hull** (convex quadratics): An auxiliary variable *t* + and a rotated second-order cone constraint ``v'Qv <= t * y`` are + introduced, and the original bound becomes ``t + c'v + d*y <= 0``. + Convexity is determined via eigenvalue decomposition of the Hessian + matrix *Q*: the quadratic is convex for an upper-bound constraint + when *Q* is positive semi-definite, and for a lower-bound constraint + when *Q* is negative semi-definite. + + **General exact hull** (non-convex quadratics and equalities): The + constraint is reformulated as ``v'Qv + c'v*y + d*y**2``, where *v* + are the disaggregated variables and *y* is the binary indicator. + + Default is False, which uses the standard perspective function for + all nonlinear constraints. + """, + ), + ) CONFIG.declare( 'eigenvalue_tolerance', cfg.ConfigValue( @@ -245,39 +278,6 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): """, ), ) - CONFIG.declare( - 'exact_hull_quadratic', - cfg.ConfigValue( - default=False, - domain=bool, - description="Use exact hull reformulation for quadratic constraints.", - doc=""" - If True, quadratic constraints (polynomial degree 2) are reformulated - using the exact hull instead of the standard perspective function. - - For a quadratic constraint of the form - - x'Qx + c'x + d <= 0 - - the reformulation depends on convexity: - - **Conic exact hull** (convex quadratics): An auxiliary variable *t* - and a rotated second-order cone constraint ``v'Qv <= t * y`` are - introduced, and the original bound becomes ``t + c'v + d*y <= 0``. - Convexity is determined via eigenvalue decomposition of the Hessian - matrix *Q*: the quadratic is convex for an upper-bound constraint - when *Q* is positive semi-definite, and for a lower-bound constraint - when *Q* is negative semi-definite. - - **General exact hull** (non-convex quadratics and equalities): The - constraint is reformulated as ``v'Qv + c'v*y + d*y**2``, where *v* - are the disaggregated variables and *y* is the binary indicator. - - Default is False, which uses the standard perspective function for - all nonlinear constraints. - """, - ), - ) transformation_name = 'hull' def __init__(self): From ed73232f307dd887043a76035d91e66a20e8c110 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Tue, 10 Mar 2026 23:57:41 -0400 Subject: [PATCH 28/39] run black --- pyomo/gdp/plugins/hull.py | 78 ++++++++++++++++++++++-------------- pyomo/gdp/tests/test_hull.py | 71 ++++++++++++-------------------- 2 files changed, 74 insertions(+), 75 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 7676c7f3334..ce6d4a900c8 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -781,9 +781,7 @@ def _transform_constraint( mode = self._config.perspective_function exact_quad = self._config.exact_hull_quadratic - use_exact_quad = ( - exact_quad and polynomial_degree == 2 - ) + use_exact_quad = exact_quad and polynomial_degree == 2 NL = polynomial_degree not in (0, 1) general_NL = NL and not use_exact_quad @@ -800,8 +798,16 @@ def _transform_constraint( if use_exact_quad: self._build_exact_quadratic_hull( - c, y, disjunct, relaxationBlock, constraint_map, - var_substitute_map, newConstraint, name, i, obj, + c, + y, + disjunct, + relaxationBlock, + constraint_map, + var_substitute_map, + newConstraint, + name, + i, + obj, ) continue @@ -937,8 +943,17 @@ def _transform_constraint( obj.deactivate() def _build_exact_quadratic_hull( - self, c, y, disjunct, relaxationBlock, constraint_map, - var_substitute_map, newConstraint, name, i, obj, + self, + c, + y, + disjunct, + relaxationBlock, + constraint_map, + var_substitute_map, + newConstraint, + name, + i, + obj, ): """Build the exact hull reformulation for a single quadratic constraint. @@ -1009,9 +1024,7 @@ def _build_exact_quadratic_hull( var_to_idx = {id(var): idx for idx, var in enumerate(all_vars)} Q = np.zeros((n_vars, n_vars)) - for coef, (var_i, var_j) in zip( - repn.quadratic_coefs, repn.quadratic_vars - ): + for coef, (var_i, var_j) in zip(repn.quadratic_coefs, repn.quadratic_vars): idx_i = var_to_idx[id(var_i)] idx_j = var_to_idx[id(var_j)] if var_i is var_j: @@ -1076,13 +1089,19 @@ def _build_exact_quadratic_hull( if need_non_convex: non_conv_expr = self._build_general_exact_hull_expr( - repn, var_substitute_map, y, const_term, + repn, var_substitute_map, y, const_term ) if use_conic_upper or use_conic_lower: conic_expr_linear = self._build_conic_exact_hull_expr( - c, y, disjunct, relaxationBlock, constraint_map, - repn, var_substitute_map, const_term, + c, + y, + disjunct, + relaxationBlock, + constraint_map, + repn, + var_substitute_map, + const_term, negate_for_conic, ) @@ -1152,9 +1171,7 @@ def _build_exact_quadratic_hull( ) constraint_map.src_constraint[newConstraint[name, 'ub']] = c - def _build_general_exact_hull_expr( - self, repn, var_substitute_map, y, const_term - ): + def _build_general_exact_hull_expr(self, repn, var_substitute_map, y, const_term): """Build the general exact hull expression for a quadratic constraint. Constructs the expression ``v'Qv + c'v*y + d*y**2`` where *v* are @@ -1179,9 +1196,7 @@ def _build_general_exact_hull_expr( """ expr = 0 - for coef, (var_i, var_j) in zip( - repn.quadratic_coefs, repn.quadratic_vars - ): + for coef, (var_i, var_j) in zip(repn.quadratic_coefs, repn.quadratic_vars): v_i = var_substitute_map.get(id(var_i), var_i) v_j = var_substitute_map.get(id(var_j), var_j) if var_i is var_j: @@ -1201,8 +1216,16 @@ def _build_general_exact_hull_expr( return expr def _build_conic_exact_hull_expr( - self, c, y, disjunct, relaxationBlock, constraint_map, - repn, var_substitute_map, const_term, negate_for_conic, + self, + c, + y, + disjunct, + relaxationBlock, + constraint_map, + repn, + var_substitute_map, + const_term, + negate_for_conic, ): """Build the conic exact hull expression for a convex quadratic. @@ -1243,9 +1266,7 @@ def _build_conic_exact_hull_expr( t = Var(domain=NonNegativeReals) t_name = unique_component_name( relaxationBlock, - '_conic_aux_t_%s' % c.getname( - fully_qualified=True, relative_to=disjunct - ), + '_conic_aux_t_%s' % c.getname(fully_qualified=True, relative_to=disjunct), ) relaxationBlock.add_component(t_name, t) @@ -1264,9 +1285,7 @@ def _build_conic_exact_hull_expr( # Build rotated SOC: v'Q_psd v <= t * y quadratic_form = 0 - for coef, (var_i, var_j) in zip( - repn.quadratic_coefs, repn.quadratic_vars - ): + for coef, (var_i, var_j) in zip(repn.quadratic_coefs, repn.quadratic_vars): v_i = var_substitute_map.get(id(var_i), var_i) v_j = var_substitute_map.get(id(var_j), var_j) actual_coef = -coef if negate_for_conic else coef @@ -1277,9 +1296,8 @@ def _build_conic_exact_hull_expr( conic_constraint_name = unique_component_name( relaxationBlock, - '_conic_constraint_%s' % c.getname( - fully_qualified=True, relative_to=disjunct - ), + '_conic_constraint_%s' + % c.getname(fully_qualified=True, relative_to=disjunct), ) conic_constraint = Constraint(expr=quadratic_form <= t * y) relaxationBlock.add_component(conic_constraint_name, conic_constraint) diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index 629624fdc9f..f7b0509be9d 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2953,15 +2953,14 @@ def test_convex_quadratic_upper_bound_conic_reformulation(self): repn = generate_standard_repn(conic_cons.body, compute_values=False) self.assertTrue(repn.is_quadratic()) quad_pairs = dict( - zip( - [(id(a), id(b)) for a, b in repn.quadratic_vars], - repn.quadratic_coefs, - ) + zip([(id(a), id(b)) for a, b in repn.quadratic_vars], repn.quadratic_coefs) ) self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) self.assertAlmostEqual(quad_pairs[(id(v_y), id(v_y))], 1) self.assertAlmostEqual( - quad_pairs.get((id(t_var), id(y_ind)), quad_pairs.get((id(y_ind), id(t_var)), None)), + quad_pairs.get( + (id(t_var), id(y_ind)), quad_pairs.get((id(y_ind), id(t_var)), None) + ), -1, ) self.assertEqual(repn.constant, 0) @@ -3013,10 +3012,7 @@ def test_convex_quadratic_lower_bound_conic_reformulation(self): repn = generate_standard_repn(conic_cons.body, compute_values=False) self.assertTrue(repn.is_quadratic()) quad_pairs = dict( - zip( - [(id(a), id(b)) for a, b in repn.quadratic_vars], - repn.quadratic_coefs, - ) + zip([(id(a), id(b)) for a, b in repn.quadratic_vars], repn.quadratic_coefs) ) self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) self.assertAlmostEqual(quad_pairs[(id(v_y), id(v_y))], 1) @@ -3067,10 +3063,7 @@ def test_nonconvex_quadratic_general_exact_hull(self): self.assertEqual(len(repn.linear_vars), 0) quad_pairs = dict( - zip( - [(id(a), id(b)) for a, b in repn.quadratic_vars], - repn.quadratic_coefs, - ) + zip([(id(a), id(b)) for a, b in repn.quadratic_vars], repn.quadratic_coefs) ) # v_x**2 coefficient: +1 self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) @@ -3116,10 +3109,7 @@ def test_equality_constraint_general_exact_hull(self): self.assertEqual(len(repn.linear_vars), 0) quad_pairs = dict( - zip( - [(id(a), id(b)) for a, b in repn.quadratic_vars], - repn.quadratic_coefs, - ) + zip([(id(a), id(b)) for a, b in repn.quadratic_vars], repn.quadratic_coefs) ) self.assertAlmostEqual(quad_pairs[(id(v_x), id(v_x))], 1) self.assertAlmostEqual(quad_pairs[(id(v_y), id(v_y))], 1) @@ -3156,24 +3146,21 @@ def test_range_constraint_mixed_reformulations(self): self.assertEqual(len(trans_cons), 3) # Pre-compute all repns to avoid duplicate calls in the search loops. - repns = {id(c): generate_standard_repn(c.body, compute_values=False) - for c in trans_cons} + repns = { + id(c): generate_standard_repn(c.body, compute_values=False) + for c in trans_cons + } # Identify each constraint by its standard representation. # The linear one (involving t) is the conic upper bound. # The quadratic one without t_var quad terms is the general lower bound. - ub_cons = next( - c for c in trans_cons - if repns[id(c)].is_linear() - ) + ub_cons = next(c for c in trans_cons if repns[id(c)].is_linear()) lb_cons = next( - c for c in trans_cons + c + for c in trans_cons if repns[id(c)].is_quadratic() - and id(t_var) not in { - id(v) - for pair in repns[id(c)].quadratic_vars - for v in pair - } + and id(t_var) + not in {id(v) for pair in repns[id(c)].quadratic_vars for v in pair} ) # --- Upper bound: linear conic bound ``t - 4*y_ind <= 0`` --- @@ -3181,7 +3168,9 @@ def test_range_constraint_mixed_reformulations(self): self.assertEqual(value(ub_cons.upper), 0) repn_ub = repns[id(ub_cons)] self.assertTrue(repn_ub.is_linear()) - linear_map = dict(zip([id(v) for v in repn_ub.linear_vars], repn_ub.linear_coefs)) + linear_map = dict( + zip([id(v) for v in repn_ub.linear_vars], repn_ub.linear_coefs) + ) self.assertAlmostEqual(linear_map[id(t_var)], 1) self.assertAlmostEqual(linear_map[id(y_ind)], -4) @@ -3316,7 +3305,8 @@ def test_conic_reformulation_with_linear_and_constant_terms(self): # The linear bound constraint is the one involving t. linear_cons = next( - c for c in trans_cons + c + for c in trans_cons if generate_standard_repn(c.body, compute_values=False).is_linear() ) repn = generate_standard_repn(linear_cons.body, compute_values=False) @@ -3366,10 +3356,7 @@ def test_general_exact_hull_with_linear_and_constant_terms(self): self.assertEqual(len(repn.linear_vars), 0) quad_pairs = dict( - zip( - [(id(a), id(b)) for a, b in repn.quadratic_vars], - repn.quadratic_coefs, - ) + zip([(id(a), id(b)) for a, b in repn.quadratic_vars], repn.quadratic_coefs) ) # v_x**2: +1 self.assertAlmostEqual(quad_pairs.get((id(v_x), id(v_x)), 0), 1) @@ -3377,8 +3364,7 @@ def test_general_exact_hull_with_linear_and_constant_terms(self): self.assertAlmostEqual(quad_pairs.get((id(v_y), id(v_y)), 0), -1) # 3*v_x*y_ind cross term cross_coef = quad_pairs.get( - (id(v_x), id(y_ind)), - quad_pairs.get((id(y_ind), id(v_x)), None), + (id(v_x), id(y_ind)), quad_pairs.get((id(y_ind), id(v_x)), None) ) self.assertIsNotNone(cross_coef, "Missing cross term v_x*y_ind") self.assertAlmostEqual(cross_coef, 3) @@ -3405,9 +3391,7 @@ def test_all_zero_eigenvalues_issues_warning(self): output = StringIO() with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): - self.hull.apply_to( - m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10 - ) + self.hull.apply_to(m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10) warning_text = output.getvalue() self.assertIn("quadratic terms", warning_text) @@ -3453,9 +3437,7 @@ def test_all_zero_eigenvalues_range_constraint_no_negation_bug(self): output = StringIO() with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): - self.hull.apply_to( - m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10 - ) + self.hull.apply_to(m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10) # Warning must be issued. self.assertIn("general exact hull", output.getvalue()) @@ -3473,8 +3455,7 @@ def test_all_zero_eigenvalues_range_constraint_no_negation_bug(self): self.assertEqual(len(trans_cons), 2) repns = [ - generate_standard_repn(tc.body, compute_values=False) - for tc in trans_cons + generate_standard_repn(tc.body, compute_values=False) for tc in trans_cons ] # Both should be quadratic (general hull includes the tiny quad terms). From 239e0382ae00f7a0a947e06ed9e2f3a53d625ed2 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Wed, 11 Mar 2026 00:14:11 -0400 Subject: [PATCH 29/39] fixed comments and formating --- coverage.xml | 520 +++++++++++++++++++++++++++++++++++ pyomo/gdp/plugins/hull.py | 10 +- pyomo/gdp/tests/test_hull.py | 22 +- 3 files changed, 537 insertions(+), 15 deletions(-) create mode 100644 coverage.xml diff --git a/coverage.xml b/coverage.xml new file mode 100644 index 00000000000..3e1e21acc4b --- /dev/null +++ b/coverage.xml @@ -0,0 +1,520 @@ + + + + + + /home/sgusev/repo/pyomo + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index ce6d4a900c8..d060fab3da1 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1047,10 +1047,12 @@ def _build_exact_quadratic_hull( # All eigenvalues lie within the tolerance band [-tol, tol]. # In this numerically ambiguous case, we conservatively avoid the # conic reformulation and fall back to the general exact hull - # reformulation. The underlying sign-flip bug motivating this - # choice occurs for two-sided (range) constraints, where a single - # conic expression would otherwise be reused (possibly negated) - # for both bounds. + # reformulation. When Q is simultaneously PSD and NSD (within + # tolerance), both use_conic_upper and use_conic_lower would be + # set for two-sided (range) constraints, but a single conic + # expression built with negate_for_conic=True cannot correctly + # serve both bounds. Falling back to the general path avoids + # this issue entirely. # Only warn for inequality constraints; equality constraints always # use the general exact hull path, so no fallback warning is needed. if not c.equality: diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py index f7b0509be9d..d9d2167ddf5 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -2915,7 +2915,7 @@ def setUp(self): self.hull = TransformationFactory('gdp.hull') # ------------------------------------------------------------------ - # 1. Convex quadratic upper-bound → conic reformulation + # 1. Convex quadratic upper-bound -> conic reformulation # ------------------------------------------------------------------ @unittest.skipUnless(numpy_available, "NumPy is not available") def test_convex_quadratic_upper_bound_conic_reformulation(self): @@ -2976,7 +2976,7 @@ def test_convex_quadratic_upper_bound_conic_reformulation(self): self.assertAlmostEqual(linear_map[id(y_ind)], -4) # ------------------------------------------------------------------ - # 2. Convex quadratic lower-bound → conic reformulation with negation + # 2. Convex quadratic lower-bound -> conic reformulation with negation # ------------------------------------------------------------------ @unittest.skipUnless(numpy_available, "NumPy is not available") def test_convex_quadratic_lower_bound_conic_reformulation(self): @@ -3027,7 +3027,7 @@ def test_convex_quadratic_lower_bound_conic_reformulation(self): self.assertAlmostEqual(linear_map[id(y_ind)], -4) # ------------------------------------------------------------------ - # 3. Non-convex quadratic → general exact hull + # 3. Non-convex quadratic -> general exact hull # ------------------------------------------------------------------ @unittest.skipUnless(numpy_available, "NumPy is not available") def test_nonconvex_quadratic_general_exact_hull(self): @@ -3073,7 +3073,7 @@ def test_nonconvex_quadratic_general_exact_hull(self): self.assertAlmostEqual(quad_pairs[(id(y_ind), id(y_ind))], -1) # ------------------------------------------------------------------ - # 4. Equality constraint → always uses general exact hull + # 4. Equality constraint -> always uses general exact hull # ------------------------------------------------------------------ @unittest.skipUnless(numpy_available, "NumPy is not available") def test_equality_constraint_general_exact_hull(self): @@ -3117,7 +3117,7 @@ def test_equality_constraint_general_exact_hull(self): self.assertAlmostEqual(quad_pairs[(id(y_ind), id(y_ind))], -4) # ------------------------------------------------------------------ - # 5. Range constraint with both bounds → mixed reformulations + # 5. Range constraint with both bounds -> mixed reformulations # ------------------------------------------------------------------ @unittest.skipUnless(numpy_available, "NumPy is not available") def test_range_constraint_mixed_reformulations(self): @@ -3218,7 +3218,7 @@ def test_eigenvalue_tolerance_permissive(self): The matrix Q = diag(1, -1e-11) has a very small negative eigenvalue (-1e-11). With a permissive tolerance of 1e-9 the eigenvalue is - within the tolerance band and Q is treated as PSD → conic + within the tolerance band and Q is treated as PSD -> conic reformulation. """ m = models.makeTwoTermDisj_NearPSDQuad() @@ -3227,7 +3227,7 @@ def test_eigenvalue_tolerance_permissive(self): m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-9 ) relaxBlock = m_perm._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] - # Conic path should have been taken → t variable present. + # Conic path should have been taken -> t variable present. self.assertIsNotNone( relaxBlock.component('_conic_aux_t_c'), "Expected conic aux variable with permissive eigenvalue_tolerance", @@ -3241,8 +3241,8 @@ def test_eigenvalue_tolerance_strict(self): """A strict ``eigenvalue_tolerance`` should reject a near-PSD Q. Same Q = diag(1, -1e-11) as above. With a tolerance of 1e-12 the - eigenvalue -1e-11 is outside the tolerance band → Q is not treated - as PSD → general exact hull (no t variable). + eigenvalue -1e-11 is outside the tolerance band -> Q is not treated + as PSD -> general exact hull (no t variable). """ m = models.makeTwoTermDisj_NearPSDQuad() @@ -3250,7 +3250,7 @@ def test_eigenvalue_tolerance_strict(self): m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-12 ) relaxBlock = m_strict._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] - # General path should have been taken → no t variable. + # General path should have been taken -> no t variable. self.assertIsNone( relaxBlock.component('_conic_aux_t_c'), "Expected no conic aux variable with strict eigenvalue_tolerance", @@ -3340,7 +3340,7 @@ def test_general_exact_hull_with_linear_and_constant_terms(self): v_y = relaxBlock.disaggregatedVars.y y_ind = m.d1.binary_indicator_var - # No conic aux variable (indefinite Q → general path). + # No conic aux variable (indefinite Q -> general path). self.assertIsNone(relaxBlock.component('_conic_aux_t_c')) trans_cons = self.hull.get_transformed_constraints(m.d1.c) From d61c54b13bb659d3ce64d5854ed7299b62b230cb Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Wed, 11 Mar 2026 00:15:45 -0400 Subject: [PATCH 30/39] removed coverage file --- coverage.xml | 520 --------------------------------------------------- 1 file changed, 520 deletions(-) delete mode 100644 coverage.xml diff --git a/coverage.xml b/coverage.xml deleted file mode 100644 index 3e1e21acc4b..00000000000 --- a/coverage.xml +++ /dev/null @@ -1,520 +0,0 @@ - - - - - - /home/sgusev/repo/pyomo - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - From dcf29062e7cce627c1067993588b716a52bc5306 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Wed, 11 Mar 2026 00:27:52 -0400 Subject: [PATCH 31/39] fixed formating of comments and unssesary h0 calculation --- pyomo/gdp/plugins/hull.py | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index d060fab3da1..296cdb714d4 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -220,17 +220,17 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): the reformulation depends on convexity: - **Conic exact hull** (convex quadratics): An auxiliary variable *t* + Conic exact hull (convex quadratics): An auxiliary variable ``t`` and a rotated second-order cone constraint ``v'Qv <= t * y`` are introduced, and the original bound becomes ``t + c'v + d*y <= 0``. Convexity is determined via eigenvalue decomposition of the Hessian - matrix *Q*: the quadratic is convex for an upper-bound constraint - when *Q* is positive semi-definite, and for a lower-bound constraint - when *Q* is negative semi-definite. + matrix ``Q``: the quadratic is convex for an upper-bound constraint + when ``Q`` is positive semi-definite, and for a lower-bound constraint + when ``Q`` is negative semi-definite. - **General exact hull** (non-convex quadratics and equalities): The - constraint is reformulated as ``v'Qv + c'v*y + d*y**2``, where *v* - are the disaggregated variables and *y* is the binary indicator. + General exact hull (non-convex quadratics and equalities): The + constraint is reformulated as ``v'Qv + c'v*y + d*y**2``, where ``v`` + are the disaggregated variables and ``y`` is the binary indicator. Default is False, which uses the standard perspective function for all nonlinear constraints. @@ -742,7 +742,7 @@ def _transform_constraint( """Transform a single Constraint on a Disjunct. Applies the appropriate hull reformulation to each - ``ConstraintData`` in *obj*. When ``exact_hull_quadratic`` is + ``ConstraintData`` in ``obj``. When ``exact_hull_quadratic`` is enabled and the constraint body has polynomial degree 2, an exact hull formulation is used instead of the perspective function. @@ -751,7 +751,7 @@ def _transform_constraint( obj : Constraint The Constraint component to transform. disjunct : Disjunct - The Disjunct that owns *obj*. + The Disjunct that owns ``obj``. var_substitute_map : dict Mapping from ``id(original_var)`` to its disaggregated counterpart. @@ -789,7 +789,7 @@ def _transform_constraint( # We need to evaluate the expression at the origin *before* # we substitute the expression variables with the # disaggregated variables - if not NL or mode == "FurmanSawayaGrossmann": + if not use_exact_quad and (not NL or mode == "FurmanSawayaGrossmann"): h_0 = clone_without_expression_components( c.body, substitute=zero_substitute_map ) @@ -958,16 +958,16 @@ def _build_exact_quadratic_hull( """Build the exact hull reformulation for a single quadratic constraint. For a constraint whose body is a quadratic of the form - ``x'Qx + c'x + d``, this method constructs either the *conic exact - hull* (when the quadratic is convex with respect to the bound - direction) or the *general exact hull* (otherwise). + ``x'Qx + c'x + d``, this method constructs either the conic exact + hull (when the quadratic is convex with respect to the bound + direction) or the general exact hull (otherwise). - **Conic exact hull** (convex case): introduces an auxiliary variable - *t >= 0* and a rotated second-order cone constraint + Conic exact hull (convex case): introduces an auxiliary variable + ``t >= 0`` and a rotated second-order cone constraint ``v'Q_psd v <= t * y``, then replaces the original bound with a linear constraint on ``t + c'v + d*y``. - **General exact hull** (non-convex / equality case): directly + General exact hull (non-convex / equality case): directly substitutes the quadratic form to ``v'Qv + c'v*y + d*y**2``. Parameters @@ -1176,9 +1176,9 @@ def _build_exact_quadratic_hull( def _build_general_exact_hull_expr(self, repn, var_substitute_map, y, const_term): """Build the general exact hull expression for a quadratic constraint. - Constructs the expression ``v'Qv + c'v*y + d*y**2`` where *v* are - disaggregated variables, *y* is the indicator, *c* are linear - coefficients, and *d* is the constant term. + Constructs the expression ``v'Qv + c'v*y + d*y**2`` where ``v`` are + disaggregated variables, ``y`` is the indicator, ``c`` are linear + coefficients, and ``d`` is the constant term. Parameters ---------- From 8bc211a822808eb4033e24c20fb22d0e516de5a7 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Wed, 11 Mar 2026 15:24:50 -0400 Subject: [PATCH 32/39] trigger CI From 27744c296b1c701033fcb500f1a2fa7a0f5ee14a Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Thu, 12 Mar 2026 23:26:08 -0400 Subject: [PATCH 33/39] added comments and reference to the paper --- pyomo/gdp/plugins/hull.py | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 296cdb714d4..b5f8ae45c74 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -108,11 +108,12 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): exact_hull_quadratic : bool If ``True``, quadratic constraints (polynomial degree 2) are reformulated using the exact hull instead of the standard - perspective function. Convex quadratics are handled with a conic - reformulation (rotated second-order cone), while non-convex - quadratics and equalities use the general exact hull - reformulation. Convexity is determined via eigenvalue - decomposition of the Hessian matrix. [default: False] + perspective function, following Gusev & Bernal Neira (2025) [4]_. + Convex quadratics are handled with a conic reformulation + (rotated second-order cone), while non-convex quadratics and + equalities use the general exact hull reformulation. Convexity + is determined via eigenvalue decomposition of the Hessian matrix. + [default: False] eigenvalue_tolerance : float Numerical tolerance for eigenvalue-based positive/negative semi-definite checks when using the exact hull reformulation for @@ -193,6 +194,11 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): useful algebraic representation of nonlinear disjunctive convex sets using the perspective function. Optimization Online (2016). http://www.optimization-online.org/DB_HTML/2016/07/5544.html. + + .. [4] Gusev, S., & Bernal Neira, D. E. (2025). Exact Hull + Reformulation for Quadratically Constrained Generalized + Disjunctive Programs. arXiv preprint arXiv:2508.16093. + https://arxiv.org/abs/2508.16093 """, ), ) @@ -212,7 +218,8 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): description="Use exact hull reformulation for quadratic constraints.", doc=""" If True, quadratic constraints (polynomial degree 2) are reformulated - using the exact hull instead of the standard perspective function. + using the exact hull instead of the standard perspective function, + following Gusev & Bernal Neira (2025), arXiv:2508.16093. For a quadratic constraint of the form @@ -957,8 +964,9 @@ def _build_exact_quadratic_hull( ): """Build the exact hull reformulation for a single quadratic constraint. - For a constraint whose body is a quadratic of the form - ``x'Qx + c'x + d``, this method constructs either the conic exact + Implements the reformulation from Gusev & Bernal Neira (2025), + arXiv:2508.16093. For a constraint whose body is a quadratic of the + form ``x'Qx + c'x + d``, this method constructs either the conic exact hull (when the quadratic is convex with respect to the bound direction) or the general exact hull (otherwise). @@ -1022,6 +1030,9 @@ def _build_exact_quadratic_hull( n_vars = len(all_vars) var_to_idx = {id(var): idx for idx, var in enumerate(all_vars)} + # Q is built symmetric (off-diagonal coefs split equally between + # Q[i,j] and Q[j,i]) because np.linalg.eigh assumes symmetry and + # x'Qx depends only on the symmetric part of Q. Q = np.zeros((n_vars, n_vars)) for coef, (var_i, var_j) in zip(repn.quadratic_coefs, repn.quadratic_vars): From e583d5ea1cd85044532463fe9f9317ecd652b040 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Sun, 15 Mar 2026 02:42:30 -0400 Subject: [PATCH 34/39] warning for mutable parameters --- pyomo/gdp/plugins/hull.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index b5f8ae45c74..1cb74961294 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1001,6 +1001,24 @@ def _build_exact_quadratic_hull( obj : Constraint The parent Constraint component (needed for ``is_indexed``). """ + # generate_standard_repn below uses compute_values=True by default, + # which evaluates mutable Params to numeric constants. The transformed + # expressions therefore won't track later Param updates. This is + # acceptable because the eigenvalue-based convexity check requires + # numeric coefficients, and a Param change could invalidate the + # PSD/NSD classification anyway. + mutable_params = list(EXPR.identify_mutable_parameters(c.body)) + if mutable_params: + logger.warning( + "GDP(Hull): Constraint '%s' contains mutable parameters %s. " + "The exact hull reformulation evaluates these to their current " + "numeric values for the eigenvalue-based convexity check; the " + "transformed constraint will not update if these parameters " + "change later.", + c.getname(fully_qualified=True), + [p.name for p in mutable_params], + ) + repn = generate_standard_repn(c.body) if not repn.is_quadratic(): From a667a187037b9eeee22f5232acbccdead3e2c433 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Sun, 15 Mar 2026 02:46:20 -0400 Subject: [PATCH 35/39] fixed syntax issue --- pyomo/gdp/plugins/hull.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 1cb74961294..7756fca36bd 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -118,9 +118,9 @@ class Hull_Reformulation(GDP_to_MIP_Transformation): Numerical tolerance for eigenvalue-based positive/negative semi-definite checks when using the exact hull reformulation for quadratic constraints (``exact_hull_quadratic=True``). An - eigenvalue :math:`\lambda` is treated as non-negative if - :math:`\lambda >= -\text{eigenvalue_tolerance}` and as - non-positive if :math:`\lambda <= \text{eigenvalue_tolerance}` + eigenvalue :math:`\\lambda` is treated as non-negative if + :math:`\\lambda >= -\\text{eigenvalue_tolerance}` and as + non-positive if :math:`\\lambda <= \\text{eigenvalue_tolerance}` (i.e., eigenvalues in ``[-eigenvalue_tolerance, eigenvalue_tolerance]`` are treated as zero). Increasing this value makes the convexity check more From 4d8e7dde1757436c7874d418b875516101fa438e Mon Sep 17 00:00:00 2001 From: Sergey Gusev <101810399+sergey-gusev94@users.noreply.github.com> Date: Sun, 15 Mar 2026 04:30:12 -0400 Subject: [PATCH 36/39] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- pyomo/gdp/plugins/hull.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index 7756fca36bd..c3713d9a837 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1022,7 +1022,7 @@ def _build_exact_quadratic_hull( repn = generate_standard_repn(c.body) if not repn.is_quadratic(): - raise RuntimeError( + raise GDP_Error( "Constraint '%s' has polynomial degree 2 but its standard " "representation is not quadratic." % c.getname(fully_qualified=True) ) From b654e88f50e684a3623a062f805f4efe87eab645 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Mon, 16 Mar 2026 19:27:29 -0400 Subject: [PATCH 37/39] added comment on conic formulation --- pyomo/gdp/plugins/hull.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index c3713d9a837..a20eecac139 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1315,6 +1315,14 @@ def _build_conic_exact_hull_expr( linear_expr += actual_const * y # Build rotated SOC: v'Q_psd v <= t * y + # NOTE: For a general PSD matrix Q this is not in canonical rotated + # second-order cone form (sum-of-squares <= product). Some solvers or + # writer interfaces may therefore treat the bilinear term t*y on the + # right-hand side as a generic nonconvex product rather than + # recognizing a rotated quadratic cone. This formulation was chosen + # deliberately based on computational testing with Gurobi and SCIP, + # where the generic quadratic form performed at least as well as (and + # sometimes better than) an explicit conic representation. quadratic_form = 0 for coef, (var_i, var_j) in zip(repn.quadratic_coefs, repn.quadratic_vars): v_i = var_substitute_map.get(id(var_i), var_i) From 71beedf581c69bccd485996c6f994f6081269db1 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Mon, 16 Mar 2026 19:28:15 -0400 Subject: [PATCH 38/39] added comment on conic formulation --- 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 a20eecac139..ed9690ea0df 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1321,8 +1321,8 @@ def _build_conic_exact_hull_expr( # right-hand side as a generic nonconvex product rather than # recognizing a rotated quadratic cone. This formulation was chosen # deliberately based on computational testing with Gurobi and SCIP, - # where the generic quadratic form performed at least as well as (and - # sometimes better than) an explicit conic representation. + # where the generic quadratic form overall performed better than + # other formulations tested. quadratic_form = 0 for coef, (var_i, var_j) in zip(repn.quadratic_coefs, repn.quadratic_vars): v_i = var_substitute_map.get(id(var_i), var_i) From f626263fcc2275f03304d7be65af024abbbe6fa9 Mon Sep 17 00:00:00 2001 From: Sergey Gusev Date: Mon, 16 Mar 2026 19:37:02 -0400 Subject: [PATCH 39/39] fixed formating with black --- pyomo/gdp/plugins/hull.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index ed9690ea0df..3daf31df982 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -1321,7 +1321,7 @@ def _build_conic_exact_hull_expr( # right-hand side as a generic nonconvex product rather than # recognizing a rotated quadratic cone. This formulation was chosen # deliberately based on computational testing with Gurobi and SCIP, - # where the generic quadratic form overall performed better than + # where the generic quadratic form overall performed better than # other formulations tested. quadratic_form = 0 for coef, (var_i, var_j) in zip(repn.quadratic_coefs, repn.quadratic_vars):