Adding (some) conic constraints.

Project.toml

@@ -5,6 +5,7 @@ version = "0.6.1"
 
 [deps]
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [weakdeps]
@@ -16,18 +17,23 @@ InfiniteDisjunctiveProgramming = "InfiniteOpt"
 [compat]
 Aqua = "0.8"
 JuMP = "1.18"
+LinearAlgebra = "1"
 Reexport = "1"
 julia = "1.10"
 Juniper = "0.9.3"
 Ipopt = "1.9.0"
 InfiniteOpt = "0.6.3"
+Hypatia = "0.10"
+Pajarito = "0.8"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
+Hypatia = "b99e6be6-89ff-11e8-14f8-45c827f4f8f2"
+Pajarito = "2f354839-79df-5901-9f0a-cdb2aac6fe30"
 
 [targets]
-test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt"]
+test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt", "Hypatia", "Pajarito"]

README.md

@@ -169,6 +169,7 @@ The following reformulation methods are currently supported:
 2. [Hull](https://optimization.cbe.cornell.edu/index.php?title=Disjunctive_inequalities#Convex-Hull_Reformulation[1][2]): The `Hull` struct is created with the following optional arguments:
 
     - `value`: `ϵ` value to use when reformulating quadratic or nonlinear constraints via the perspective function proposed by [Furman, et al. [2020]](https://link.springer.com/article/10.1007/s10589-020-00176-0). Default: `1e-6`. `ϵ` values are currently global to the model. Constraint specific tolerances can be supported in future releases.
+    - `quadratic`: Reformulation to use for quadratic disjunct constraints. Default: `:epsilon` (the ε-approximated perspective function above). The exact hull reformulations of [Gusev & Bernal Neira [2025]](https://arxiv.org/abs/2508.16093) are available as `:exact` (automatically uses CEHR for constraints with a convex quadratic part and GEHR otherwise), `:gehr`, `:cehr`, and `:cehr_conic` (CEHR with the cone written out explicitly as a rotated second-order cone for solvers that use cones natively).
 
 3. [Indicator](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints): This method reformulates each disjunct constraint into an indicator constraint with the Boolean reformulation counterpart of the Logical variable used to define the disjunct constraint.
 
@@ -192,6 +193,22 @@ The following reformulation methods are currently supported:
     - `final_reform_method`: Reformulation method to apply after cutting plane iterations. Default: `BigM()`.
     - `M_value`: Big-M value to use in the relaxed Big-M reformulation during iterations. Default: `1e9`.
 
+## Conic Constraints
+
+Disjunct constraints can also be defined with the conic sets `SecondOrderCone`, `RotatedSecondOrderCone`, `MOI.ExponentialCone`, and `MOI.PowerCone`. Big-M reformulates these by relaxing the constraint function along a fixed interior direction of the cone, and Hull produces the exact conic hull of [Bernal Neira & Grossmann [2021]](https://arxiv.org/abs/2109.09657).
+
+```julia
+using DisjunctiveProgramming
+
+model = GDPModel()
+@variable(model, -5 <= x[1:2] <= 5)
+@variable(model, 0 <= t <= 10)
+@variable(model, Y[1:2], Logical)
+@constraint(model, [t; x] in SecondOrderCone(), Disjunct(Y[1]))          # t >= ‖x‖
+@constraint(model, [x[1], 1, t] in MOI.ExponentialCone(), Disjunct(Y[2])) # t >= exp(x[1])
+@disjunction(model, Y)
+```
+
 ## Infinite-Dimensional GDP
 To model disjunctions, logical variables, and logical constraints with infinite-dimensional optimization problems (e.g., dynamic and stochastic optimization), DisjunctiveProgramming is also compatible with [InfiniteOpt.jl](https://github.com/infiniteopt/InfiniteOpt.jl). For this, the syntax is largely the same, users simply need to import `InfiniteOpt` and use `InfiniteGDPModel`. They also can use `InfiniteLogical` to declare infinite logical variables as shown below:
 ```julia

docs/src/index.md

@@ -167,6 +167,8 @@ The following reformulation methods are currently supported:
 
 3. [Indicator](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints): This method reformulates each disjunct constraint into an indicator constraint with the Boolean reformulation counterpart of the Logical variable used to define the disjunct constraint. This is invoked with [`Indicator`](@ref).
 
+Disjunct constraints can also be defined with the conic sets `SecondOrderCone`, `RotatedSecondOrderCone`, `MOI.ExponentialCone`, and `MOI.PowerCone`; these are supported by the [`BigM`](@ref) and [`Hull`](@ref) reformulations, where the hull of a conic constraint is exact. Exact hull reformulations for quadratic disjunct constraints (GEHR/CEHR) are available via the `quadratic` keyword argument of [`Hull`](@ref).
+
 ## Release Notes
 
 Prior to `v0.4.0`, the package did not leverage the JuMP extension capabilities and was not as robust. For these earlier releases, refer to [Perez, Joshi, and Grossmann, 2023](https://arxiv.org/abs/2304.10492v1) and the following [JuliaCon 2022 Talk](https://www.youtube.com/watch?v=AMIrgTTfUkI).

src/DisjunctiveProgramming.jl

@@ -6,6 +6,8 @@ Reexport.@reexport using JuMP
 
 # Use Meta for metaprogramming
 using Base.Meta
+# Convexity checks for the exact quadratic hull reformulations
+import LinearAlgebra
 # Create aliases
 import JuMP.MOI as _MOI
 import JuMP.MOIU.CleverDicts as _MOIUC

src/bigm.jl

@@ -251,3 +251,45 @@ function reformulate_disjunct_constraint(
     reform_con_np = JuMP.build_constraint(error, new_func_np, _MOI.Nonpositives(con.set.dimension))
     return [reform_con_nn, reform_con_np]
 end
+
+################################################################################
+#                       BIG-M FOR CONIC CONSTRAINTS
+################################################################################
+# Big-M for `func in K`: add slack `M*(1 - y)*d` along a fixed interior
+# direction `d` of the cone, so `func + M*d in K`.
+
+# SOC: (t, x...) with t >= ||x||.  d = e1 = (1, 0, ...); interior since
+# 1 > ||0|| = 0.  Bumping t alone makes (t + M) >= ||x|| hold for big M.
+_conic_bigm_direction(set::_MOI.SecondOrderCone) =
+    [i == 1 ? 1.0 : 0.0 for i in 1:_MOI.dimension(set)]
+# Rotated SOC: (t, u, x...) with 2*t*u >= ||x||^2, t,u >= 0.  d =
+# (1,1,0,...); interior since 2*1*1 = 2 > 0.  Both t,u must grow so
+# 2(t+M)(u+M) ~ 2*M^2 dominates ||x||^2; bumping one leaves 2*t*0 = 0.
+_conic_bigm_direction(set::_MOI.RotatedSecondOrderCone) =
+    [i <= 2 ? 1.0 : 0.0 for i in 1:_MOI.dimension(set)]
+# Exp cone: (x, y, z) with z >= y*exp(x/y), y >= 0.  d = (0, 1, 2);
+# interior since 2 > 1*exp(0) = 1.  As M grows, exp(x/(y+M)) -> 1 so the
+# RHS ~ y + M ~ M while z grows as 2M; the factor 2 keeps z above it.
+_conic_bigm_direction(::_MOI.ExponentialCone) = [0.0, 1.0, 2.0]
+# Power cone: (x, y, z) with x^a * y^(1-a) >= |z|, x,y >= 0, a the cone
+# exponent.  d = (1,1,0); interior since 1^a * 1^(1-a) = 1 > 0.  Bumping
+# x,y makes (x+M)^a (y+M)^(1-a) ~ M dominate the fixed |z|; z untouched.
+_conic_bigm_direction(::_MOI.PowerCone) = [1.0, 1.0, 0.0]
+
+function reformulate_disjunct_constraint(
+    model::JuMP.AbstractModel,
+    con::JuMP.VectorConstraint{T, S, R},
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::BigM
+) where {
+    T <: Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    S <: _ConicSets, R
+}
+    M = method.value
+    d = _conic_bigm_direction(con.set)
+    new_func = JuMP.@expression(model, [i=1:_MOI.dimension(con.set)],
+        con.func[i] + M*(1 - bvref)*d[i]
+    )
+    reform_con = JuMP.build_constraint(error, new_func, con.set)
+    return [reform_con]
+end

src/constraints.jl

@@ -207,10 +207,23 @@ function JuMP.build_constraint(
     return _DisjunctConstraint(constr, tag.indicator)
 end
 
+# MOI conic sets handled by the BigM and Hull conic reformulations
+# (Bernal Neira & Grossmann 2021). The affine map A*x - b sits inside
+# the cone; the nonlinearity is carried by the cone itself.
+const _ConicSets = Union{
+    _MOI.SecondOrderCone,
+    _MOI.RotatedSecondOrderCone,
+    _MOI.ExponentialCone,
+    _MOI.PowerCone,
+}
+
 # Allows for building DisjunctConstraints for VectorConstraints since these get parsed differently by JuMP (JuMP changes the set to a MOI.AbstractScalarSet)
 for SetType in (
     JuMP.Nonnegatives, JuMP.Nonpositives, JuMP.Zeros,
-    _MOI.Nonnegatives, _MOI.Nonpositives, _MOI.Zeros
+    _MOI.Nonnegatives, _MOI.Nonpositives, _MOI.Zeros,
+    JuMP.SecondOrderCone, JuMP.RotatedSecondOrderCone,
+    _MOI.SecondOrderCone, _MOI.RotatedSecondOrderCone,
+    _MOI.ExponentialCone, _MOI.PowerCone
 )
     @eval begin
         @doc """

src/datatypes.jl

@@ -409,28 +409,51 @@ end
 """
     Hull{T} <: AbstractReformulationMethod
 
-A type for using the convex hull reformulation approach for disjunctive 
+A type for using the convex hull reformulation approach for disjunctive
 constraints.
 
 **Fields**
 - `value::T`: epsilon value for nonlinear hull reformulations (default = `1e-6`).
+- `quadratic::Symbol`: reformulation used for quadratic disjunct
+  constraints (default = `:epsilon`). Options are:
+  - `:epsilon`: ε-approximated perspective (Furman, Sawaya & Grossmann
+    2020).
+  - `:exact`: exact hull (Gusev & Bernal Neira 2025); routes each
+    constraint to CEHR when its quadratic part is convex and to GEHR
+    otherwise (equality constraints always use GEHR since they are
+    nonconvex).
+  - `:gehr`: always use the General Exact Hull Reformulation.
+  - `:cehr`: always use the Conic Exact Hull Reformulation (errors on
+    nonconvex quadratic constraints).
+  - `:cehr_conic`: CEHR with the cone written out explicitly as a
+    rotated second-order cone via a spectral factorization of the
+    quadratic part, for solvers that consume cones natively (errors on
+    nonconvex quadratic constraints).
 """
 struct Hull{T} <: AbstractReformulationMethod
     value::T
-    function Hull(ϵ::T = 1e-6) where {T}
-        new{T}(ϵ)
+    quadratic::Symbol
+    function Hull(ϵ::T = 1e-6; quadratic::Symbol = :epsilon) where {T}
+        if !(quadratic in (:epsilon, :exact, :gehr, :cehr, :cehr_conic))
+            error("Invalid `quadratic` option `:$(quadratic)`. Choose " *
+                  "from `:epsilon`, `:exact`, `:gehr`, `:cehr`, or " *
+                  "`:cehr_conic`.")
+        end
+        new{T}(ϵ, quadratic)
     end
 end
 
 # temp struct to store variable disaggregations (reset for each disjunction)
 mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationMethod
     value::T
+    quadratic::Symbol
     disjunction_variables::Dict{V, Vector{V}}
     disjunct_variables::Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V}
     function _Hull(method::Hull{T}, vrefs::Set{V}) where {T, V <: JuMP.AbstractVariableRef}
         new{V, T}(
             method.value,
-            Dict{V, Vector{V}}(vref => V[] for vref in vrefs), 
+            method.quadratic,
+            Dict{V, Vector{V}}(vref => V[] for vref in vrefs),
             Dict{Tuple{V, Union{V, JuMP.GenericAffExpr{T, V}}}, V}()
         )
     end