diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py index f7e7d7a4b9d..3daf31df982 100644 --- a/pyomo/gdp/plugins/hull.py +++ b/pyomo/gdp/plugins/hull.py @@ -11,6 +11,8 @@ from collections import defaultdict +from pyomo.common.dependencies import numpy as np, numpy_available + from pyomo.common.autoslots import AutoSlots import pyomo.common.config as cfg from pyomo.common import deprecated @@ -19,6 +21,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 +38,7 @@ Any, RangeSet, Reals, + NonNegativeReals, value, NonNegativeIntegers, Binary, @@ -101,6 +105,27 @@ 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, 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 + 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}` + (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] @@ -169,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 """, ), ) @@ -180,6 +210,64 @@ 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, + following Gusev & Bernal Neira (2025), arXiv:2508.16093. + + 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( + default=1e-10, + domain=cfg.NonNegativeFloat, + 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 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. + """, + ), + ) CONFIG.declare( 'assume_fixed_vars_permanent', cfg.ConfigValue( @@ -658,6 +746,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,20 +783,42 @@ 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 # 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 ) 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, @@ -724,7 +853,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 +904,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 +926,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 +949,403 @@ 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. + + 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). + + 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``). + """ + # 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(): + raise GDP_Error( + "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 + + 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() + 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 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): + 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 = 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) + + # Determine which bounds can use the conic formulation + use_conic_upper = False + use_conic_lower = False + negate_for_conic = False + + if Q_is_psd and Q_is_nsd: + # 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. 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: + 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: + 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 + 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 + # 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 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) + 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. diff --git a/pyomo/gdp/tests/models.py b/pyomo/gdp/tests/models.py index dde0b007c31..c34bcfbea9b 100644 --- a/pyomo/gdp/tests/models.py +++ b/pyomo/gdp/tests/models.py @@ -1262,3 +1262,311 @@ 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 + + +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 + + +def makeTwoTermDisj_ConvexQuadCrossProduct(): + """Two-term disjunction with a convex quadratic constraint containing + cross-product (off-diagonal) terms. + + Disjunct d1 has ``x**2 + x*y + y**2 <= 4``. The Q matrix is + ``[[1, 0.5], [0.5, 1]]`` (eigenvalues 0.5 and 1.5, both positive), + so Q is PSD and the conic exact hull reformulation should be used. + """ + 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.x * m.y + 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_NonconvexQuadCrossProduct(): + """Two-term disjunction with a non-convex quadratic constraint containing + cross-product (off-diagonal) terms. + + Disjunct d1 has ``x**2 + 3*x*y - y**2 <= 1``. The Q matrix is + ``[[1, 1.5], [1.5, -1]]`` (eigenvalues approximately -1.803 and 1.803), + so Q is indefinite and the general exact hull reformulation should be used. + """ + 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 * m.y - 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_QuadMutableParam(): + """Two-term disjunction with a quadratic constraint containing a mutable + Param as a coefficient. + + Disjunct d1 has ``p*x**2 + y**2 <= 4`` where ``p`` is a mutable Param + (default value 1). The exact hull reformulation should issue a warning + about the mutable parameter being evaluated to its current numeric value. + """ + m = ConcreteModel() + m.x = Var(bounds=(-2, 2)) + m.y = Var(bounds=(-2, 2)) + m.p = Param(initialize=1, mutable=True) + m.d1 = Disjunct() + m.d1.c = Constraint(expr=m.p * 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_IndexedQuadEq(): + """Two-term disjunction with an IndexedConstraint containing quadratic + equality constraints. + + Disjunct d1 has ``x[s]**2 + y[s]**2 == 4`` for s in {1, 2}. + Used to test the ``obj.is_indexed()`` branch for equality constraints + in the exact hull reformulation. + """ + m = ConcreteModel() + m.s = Set(initialize=[1, 2]) + m.x = Var(m.s, bounds=(-2, 2)) + m.y = Var(m.s, bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(m.s, rule=lambda _, s: m.x[s] ** 2 + m.y[s] ** 2 == 4) + m.d2 = Disjunct() + m.d2.c = Constraint(m.s, rule=lambda _, s: m.x[s] + m.y[s] <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_IndexedQuadUB(): + """Two-term disjunction with an IndexedConstraint containing convex + quadratic upper-bound constraints. + + Disjunct d1 has ``x[s]**2 + y[s]**2 <= 4`` for s in {1, 2}. + Used to test the ``obj.is_indexed()`` branch for upper-bound constraints + in the exact hull reformulation. + """ + m = ConcreteModel() + m.s = Set(initialize=[1, 2]) + m.x = Var(m.s, bounds=(-2, 2)) + m.y = Var(m.s, bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(m.s, rule=lambda _, s: m.x[s] ** 2 + m.y[s] ** 2 <= 4) + m.d2 = Disjunct() + m.d2.c = Constraint(m.s, rule=lambda _, s: m.x[s] + m.y[s] <= 1) + m.disj = Disjunction(expr=[m.d1, m.d2]) + return m + + +def makeTwoTermDisj_IndexedQuadLB(): + """Two-term disjunction with an IndexedConstraint containing convex + quadratic lower-bound constraints. + + Disjunct d1 has ``-x[s]**2 - y[s]**2 >= -4`` for s in {1, 2}. + Used to test the ``obj.is_indexed()`` branch for lower-bound constraints + in the exact hull reformulation. + """ + m = ConcreteModel() + m.s = Set(initialize=[1, 2]) + m.x = Var(m.s, bounds=(-2, 2)) + m.y = Var(m.s, bounds=(-2, 2)) + m.d1 = Disjunct() + m.d1.c = Constraint(m.s, rule=lambda _, s: -m.x[s] ** 2 - m.y[s] ** 2 >= -4) + m.d2 = Disjunct() + m.d2.c = Constraint(m.s, rule=lambda _, s: m.x[s] + m.y[s] <= 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 8c59ef0bb55..bc9b430d59a 100644 --- a/pyomo/gdp/tests/test_hull.py +++ b/pyomo/gdp/tests/test_hull.py @@ -8,13 +8,15 @@ # ____________________________________________________________________________________ import logging +import re import sys import random 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 +38,7 @@ Param, Objective, TerminationCondition, + NonNegativeReals, ) from pyomo.core.expr.compare import ( assertExpressionsEqual, @@ -47,6 +50,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 +2894,845 @@ class NestedDisjunctsInFlatGDP(unittest.TestCase): def test_declare_disjuncts_in_disjunction_rule(self): ct.check_nested_disjuncts_in_flat_gdp(self, 'hull') + + +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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_ConvexQuadUB() + + 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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_ConvexQuadLB() + + 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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_NonconvexQuad() + + 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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_QuadEquality() + + 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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_QuadRange() + + 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 = models.makeTwoTermDisj_ConvexQuadUB() + + 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 + # ------------------------------------------------------------------ + @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 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-9 + ) + 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 + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + 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 = models.makeTwoTermDisj_NearPSDQuad() + + 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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_ConvexQuadUB() + + 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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_ConvexQuad_LinearTerms() + + 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 + # ------------------------------------------------------------------ + @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. + + 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 = models.makeTwoTermDisj_NonconvexQuad_LinearTerms() + + 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) + + # ------------------------------------------------------------------ + # 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): + """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 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). + """ + m = models.makeTwoTermDisj_AllZeroEigenvalueQuad() + + output = StringIO() + with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): + self.hull.apply_to(m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10) + + warning_text = output.getvalue() + self.assertIn("quadratic terms", warning_text) + self.assertIn("eigenvalues", 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'). + 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. + """ + m = models.makeTwoTermDisj_AllZeroEigenvalueQuadRange() + + output = StringIO() + with LoggingIntercept(output, 'pyomo.gdp.hull', logging.WARNING): + self.hull.apply_to(m, exact_hull_quadratic=True, eigenvalue_tolerance=1e-10) + + # Warning must be issued. + self.assertIn("general exact hull", 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}", + ) + + # ------------------------------------------------------------------ + # 14. Cross-product terms in conic reformulation (PSD Q with off-diagonal) + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_conic_reformulation_with_cross_product_terms(self): + """PSD Q with off-diagonal terms should produce a conic reformulation + that includes cross-product terms in the rotated SOC constraint. + + For ``x**2 + x*y + y**2 <= 4`` the Q matrix is + ``[[1, 0.5], [0.5, 1]]`` (eigenvalues 0.5 and 1.5, PSD). + The conic constraint should contain the cross-product ``v_x * v_y`` + with coefficient 1 (the original x*y coefficient). + """ + m = models.makeTwoTermDisj_ConvexQuadCrossProduct() + + 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 + + t_var = relaxBlock.component('_conic_aux_t_c') + self.assertIsNotNone( + t_var, "Expected conic 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 = next( + c + for c in trans_cons + if generate_standard_repn(c.body, compute_values=False).is_quadratic() + ) + linear_cons = next( + c + for c in trans_cons + if generate_standard_repn(c.body, compute_values=False).is_linear() + ) + + # --- Conic constraint: v_x**2 + v_x*v_y + 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) + cross_coef = quad_pairs.get( + (id(v_x), id(v_y)), quad_pairs.get((id(v_y), id(v_x)), None) + ) + self.assertIsNotNone(cross_coef, "Missing cross-product term v_x*v_y") + self.assertAlmostEqual(cross_coef, 1) + + # --- 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) + + # ------------------------------------------------------------------ + # 15. Cross-product terms in general exact hull (indefinite Q with off-diagonal) + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_general_exact_hull_with_cross_product_terms(self): + """Indefinite Q with off-diagonal terms should produce a general exact + hull expression that includes cross-product terms. + + For ``x**2 + 3*x*y - y**2 <= 1`` the Q matrix is + ``[[1, 1.5], [1.5, -1]]`` (indefinite). The general exact hull + constraint should be: + ``v_x**2 + 3*v_x*v_y - v_y**2 - y_ind**2 <= 0`` + """ + m = models.makeTwoTermDisj_NonconvexQuadCrossProduct() + + 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 + + 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*v_y cross-product term + cross_coef = quad_pairs.get( + (id(v_x), id(v_y)), quad_pairs.get((id(v_y), id(v_x)), None) + ) + self.assertIsNotNone(cross_coef, "Missing cross-product term v_x*v_y") + self.assertAlmostEqual(cross_coef, 3) + # -1*y_ind**2 (from rhs: 1 * y_ind**2 moved to LHS) + self.assertAlmostEqual(quad_pairs.get((id(y_ind), id(y_ind)), 0), -1) + + # ------------------------------------------------------------------ + # 16. Mutable parameter warning in exact hull quadratic + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_mutable_param_warning(self): + """A quadratic constraint with a mutable Param should trigger a warning. + + For ``p*x**2 + y**2 <= 4`` where ``p`` is a mutable Param, the exact + hull reformulation evaluates ``p`` to its current numeric value for the + eigenvalue-based convexity check. A warning must be emitted mentioning + the constraint name and the mutable parameter name. The transformation + should still succeed (using the evaluated value). + """ + m = models.makeTwoTermDisj_QuadMutableParam() + + 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("mutable parameters", warning_text) + self.assertIn("d1.c", warning_text) + self.assertIn("p", warning_text) + self.assertIn("will not update", warning_text) + + # The transformation should still produce a valid result. + # With p=1, Q = I (PSD), so the conic path should be taken. + relaxBlock = m._pyomo_gdp_hull_reformulation.relaxedDisjuncts[0] + t_var = relaxBlock.component('_conic_aux_t_c') + self.assertIsNotNone( + t_var, "Expected conic auxiliary variable (p=1 makes Q=I, PSD)" + ) + + # ------------------------------------------------------------------ + # 17. Indexed quadratic equality -> is_indexed() branch + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_indexed_quadratic_equality(self): + """An IndexedConstraint with quadratic equalities should use the + ``obj.is_indexed()`` branch in the exact hull equality path. + + For ``x[s]**2 + y[s]**2 == 4`` (s in {1,2}), each ConstraintData + should produce one general exact hull equality constraint indexed + with the ``(name, i, 'eq')`` tuple key. + """ + m = models.makeTwoTermDisj_IndexedQuadEq() + + self.hull.apply_to(m, exact_hull_quadratic=True) + + for s in m.s: + trans_cons = self.hull.get_transformed_constraints(m.d1.c[s]) + self.assertEqual(len(trans_cons), 1) + eq_cons = trans_cons[0] + self.assertEqual(value(eq_cons.lower), 0) + self.assertEqual(value(eq_cons.upper), 0) + self.assertIs(self.hull.get_src_constraint(eq_cons), m.d1.c[s]) + + # ------------------------------------------------------------------ + # 18. Indexed quadratic upper-bound -> is_indexed() branch + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_indexed_quadratic_upper_bound(self): + """An IndexedConstraint with convex quadratic upper bounds should use + the ``obj.is_indexed()`` branch in the conic upper-bound path. + + For ``x[s]**2 + y[s]**2 <= 4`` (s in {1,2}, Q = I, PSD), each + ConstraintData should produce a conic SOC constraint and a linear + bound constraint, both indexed with ``(name, i, ...)`` tuple keys. + """ + m = models.makeTwoTermDisj_IndexedQuadUB() + + self.hull.apply_to(m, exact_hull_quadratic=True) + + for s in m.s: + trans_cons = self.hull.get_transformed_constraints(m.d1.c[s]) + self.assertEqual(len(trans_cons), 2) + for tc in trans_cons: + self.assertIs(self.hull.get_src_constraint(tc), m.d1.c[s]) + + # ------------------------------------------------------------------ + # 19. Indexed quadratic lower-bound -> is_indexed() branch + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_indexed_quadratic_lower_bound(self): + """An IndexedConstraint with convex quadratic lower bounds should use + the ``obj.is_indexed()`` branch in the conic lower-bound path. + + For ``-x[s]**2 - y[s]**2 >= -4`` (s in {1,2}, Q = -I, NSD), each + ConstraintData should produce a conic SOC constraint and a linear + bound constraint via negation, both indexed with ``(name, i, ...)`` + tuple keys. + """ + m = models.makeTwoTermDisj_IndexedQuadLB() + + self.hull.apply_to(m, exact_hull_quadratic=True) + + for s in m.s: + trans_cons = self.hull.get_transformed_constraints(m.d1.c[s]) + self.assertEqual(len(trans_cons), 2) + for tc in trans_cons: + self.assertIs(self.hull.get_src_constraint(tc), m.d1.c[s]) + + # ------------------------------------------------------------------ + # 20. Debug logging for exact hull quadratic bounds + # ------------------------------------------------------------------ + @unittest.skipUnless(numpy_available, "NumPy is not available") + def test_debug_logging_exact_hull_quadratic(self): + """When the logger is at DEBUG level, the exact hull reformulation + should emit debug messages for lower-bound and upper-bound quadratic + constraints. + + Uses a range constraint ``1 <= x**2 + y**2 <= 4`` so that both the + lower-bound and upper-bound debug paths are exercised. + """ + m = models.makeTwoTermDisj_QuadRange() + + output = StringIO() + with LoggingIntercept(output, 'pyomo.gdp.hull', logging.DEBUG): + self.hull.apply_to(m, exact_hull_quadratic=True) + + debug_text = output.getvalue() + self.assertIn("GDP(Hull): Transforming constraint", debug_text) + self.assertIn("d1.c", debug_text)