diff --git a/Project.toml b/Project.toml
index bdec9be4..a10eb102 100644
--- a/Project.toml
+++ b/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,11 +17,14 @@ 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"
@@ -28,6 +32,8 @@ 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"]
diff --git a/README.md b/README.md
index 543d1f7a..1e9a96ef 100644
--- a/README.md
+++ b/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
diff --git a/docs/src/index.md b/docs/src/index.md
index aa1902e3..989148de 100644
--- a/docs/src/index.md
+++ b/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).
diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl
index 5f77dce7..cc036b66 100644
--- a/ext/InfiniteDisjunctiveProgramming.jl
+++ b/ext/InfiniteDisjunctiveProgramming.jl
@@ -148,6 +148,22 @@ function DP.get_constant(
     return constant + param_expr
 end
 
+# Parameters act as constants, so their perspective scales with the
+# indicator (e.g. a bare parameter row of a conic constraint)
+function DP.disaggregate_expression(
+    model::InfiniteOpt.InfiniteModel,
+    vref::InfiniteOpt.GeneralVariableRef,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::DP._Hull
+)
+    _is_parameter(vref) && return vref * bvref
+    if JuMP.is_binary(vref) ||
+        !haskey(method.disjunct_variables, (vref, bvref))
+        return vref
+    end
+    return method.disjunct_variables[vref, bvref]
+end
+
 function DP.disaggregate_expression(
     model::M,
     aff::JuMP.GenericAffExpr,
diff --git a/src/DisjunctiveProgramming.jl b/src/DisjunctiveProgramming.jl
index 85f03455..89ecda1f 100644
--- a/src/DisjunctiveProgramming.jl
+++ b/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
diff --git a/src/bigm.jl b/src/bigm.jl
index 81386b3f..4e2228a4 100644
--- a/src/bigm.jl
+++ b/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
diff --git a/src/constraints.jl b/src/constraints.jl
index a25270d4..2bcecc72 100644
--- a/src/constraints.jl
+++ b/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 """
diff --git a/src/datatypes.jl b/src/datatypes.jl
index bdcd4c2e..3b7a2a5b 100644
--- a/src/datatypes.jl
+++ b/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
diff --git a/src/hull.jl b/src/hull.jl
index 779369e9..e493850f 100644
--- a/src/hull.jl
+++ b/src/hull.jl
@@ -135,13 +135,12 @@ function disaggregate_expression(
     bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr}, 
     method::_Hull
     )
-    #get affine part
-    new_expr = disaggregate_expression(model, quad.aff, bvref, method)
-    #get quadratic part
+    quad_part, affine_part = _split_quad_terms(quad)
+    new_expr = _disaggregate_affine_terms(model, affine_part, bvref, method)
     ϵ = method.value
-    for (pair, coeff) in quad.terms
-        da_ref = method.disjunct_variables[pair.a, bvref]
-        db_ref = method.disjunct_variables[pair.b, bvref]
+    for (pair, coeff) in quad_part.terms
+        da_ref = disaggregate_expression(model, pair.a, bvref, method)
+        db_ref = disaggregate_expression(model, pair.b, bvref, method)
         new_expr += coeff * da_ref * db_ref / ((1-ϵ)*bvref+ϵ)
     end
     return new_expr
@@ -254,7 +253,8 @@ function reformulate_disjunction(model::JuMP.AbstractModel, disj::Disjunction, m
     return ref_cons
 end
 function reformulate_disjunction(model::JuMP.AbstractModel, disj::Disjunction, method::_Hull)
-    return reformulate_disjunction(model, disj, Hull(method.value))
+    hull = Hull(method.value; quadratic = method.quadratic)
+    return reformulate_disjunction(model, disj, hull)
 end
 
 function reformulate_disjunct_constraint(
@@ -269,15 +269,41 @@ function reformulate_disjunct_constraint(
     reform_con = JuMP.build_constraint(error, new_func, S(0))
     return [reform_con]
 end
+
+# Promote disaggregated vector rows to one concrete expression type
+# (rows mix affine and quadratic when parameters scale the indicator)
+function _combine_expressions(exprs::Vector)
+    T = mapreduce(typeof, promote_type, exprs)
+    return convert.(T, exprs)
+end
+
 function reformulate_disjunct_constraint(
-    model::JuMP.AbstractModel, 
-    con::JuMP.VectorConstraint{T, S, R}, 
-    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr}, 
+    model::JuMP.AbstractModel,
+    con::JuMP.VectorConstraint{T, S, R},
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
     method::_Hull
 ) where {T <: JuMP.AbstractJuMPScalar, S <: Union{_MOI.Nonpositives, _MOI.Nonnegatives, _MOI.Zeros}, R}
-    new_func = JuMP.@expression(model, [i=1:con.set.dimension],
-        disaggregate_expression(model, con.func[i], bvref, method)
-    )
+    new_func = _combine_expressions([
+        disaggregate_expression(model, func, bvref, method)
+        for func in con.func
+    ])
+    reform_con = JuMP.build_constraint(error, new_func, con.set)
+    return [reform_con]
+end
+
+function reformulate_disjunct_constraint(
+    model::JuMP.AbstractModel,
+    con::JuMP.VectorConstraint{T, S, R},
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+) where {
+    T <: Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    S <: _ConicSets, R
+}
+    new_func = _combine_expressions([
+        disaggregate_expression(model, func, bvref, method)
+        for func in con.func
+    ])
     reform_con = JuMP.build_constraint(error, new_func, con.set)
     return [reform_con]
 end
@@ -332,8 +358,8 @@ function reformulate_disjunct_constraint(
     return [reform_con_gt, reform_con_lt]
 end
 function reformulate_disjunct_constraint(
-    model::JuMP.AbstractModel, 
-    con::JuMP.ScalarConstraint{T, S}, 
+    model::JuMP.AbstractModel,
+    con::JuMP.ScalarConstraint{T, S},
     bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
     method::_Hull
 ) where {T <: JuMP.GenericNonlinearExpr, S <: _MOI.Interval}
@@ -350,3 +376,304 @@ function reformulate_disjunct_constraint(
     reform_con_lt = JuMP.build_constraint(error, new_func_lt, _MOI.LessThan(0))
     return [reform_con_gt, reform_con_lt]
 end
+
+################################################################################
+#                       EXACT QUADRATIC HULL (GEHR / CEHR)
+################################################################################
+# Exact hull reformulations for quadratic disjunct constraints
+# (Gusev & Bernal Neira 2025, arXiv:2508.16093), replacing the
+# ε-approximated perspective. Both are exact for the full relaxation
+# y ∈ [0, 1] given finite variable bounds (the disaggregated variable
+# bounds force ν = 0 when y = 0):
+# - GEHR (Eq. 13): multiply cl h̃(ν, y) ≤ 0 through by y to get
+#   ν'Qν + (a'ν)y + d*y² ≤ 0. Valid for any Q (nonconvex included)
+#   and for equalities, but the function is nonconvex in (ν, y).
+# - CEHR (Eq. 22): for convex constraints (Q ⪰ 0), epigraph variable
+#   t ≥ 0 with ν'Qν - t*y ≤ 0 and t + a'ν + d*y ≤ 0. The cone
+#   constraint is left unfactorized so conic-aware solvers recognize
+#   the rotated SOC in presolve; at y = 0 the link forces t = 0.
+
+# Assemble the symmetric coefficient matrix of the quadratic terms and
+# the variables indexing its rows, in order of first appearance
+function _quad_coefficient_matrix(
+    quad::JuMP.GenericQuadExpr{C, V}
+    ) where {C, V}
+    index = Dict{V, Int}()
+    vars = V[]
+    for (pair, _) in quad.terms, v in (pair.a, pair.b)
+        haskey(index, v) || (push!(vars, v); index[v] = length(vars))
+    end
+    quad_matrix = zeros(float(C), length(vars), length(vars))
+    for (pair, coeff) in quad.terms
+        i, j = index[pair.a], index[pair.b]
+        if i == j
+            quad_matrix[i, i] += coeff
+        else
+            quad_matrix[i, j] += coeff / 2
+            quad_matrix[j, i] += coeff / 2
+        end
+    end
+    return quad_matrix, vars
+end
+
+# Check whether the quadratic part is convex (Q ⪰ 0 up to a tolerance)
+function _is_convex_quad(quad::JuMP.GenericQuadExpr)
+    quad_matrix, _ = _quad_coefficient_matrix(quad)
+    tol = 1e-9 * max(one(eltype(quad_matrix)), maximum(abs, quad_matrix))
+    return LinearAlgebra.eigmin(LinearAlgebra.Symmetric(quad_matrix)) >= -tol
+end
+
+# Split the quadratic terms into `quad_part` (both factors get
+# disaggregated) and `affine_part` (coefficient-only factors such as
+# infinite parameters, plus the affine part), which is affine in the
+# disaggregated variables
+function _split_quad_terms(quad::JuMP.GenericQuadExpr)
+    quad_part = zero(typeof(quad))
+    affine_part = typeof(quad)(copy(quad.aff))
+    for (pair, coeff) in quad.terms
+        if requires_disaggregation(pair.a) &&
+            requires_disaggregation(pair.b)
+            JuMP.add_to_expression!(quad_part, coeff, pair.a, pair.b)
+        else
+            JuMP.add_to_expression!(affine_part, coeff, pair.a, pair.b)
+        end
+    end
+    return quad_part, affine_part
+end
+
+# Disaggregate the remainder terms: coefficient-only factors scale
+# with the indicator like constants do
+function _disaggregate_affine_terms(
+    model::JuMP.AbstractModel,
+    affine_part::JuMP.GenericQuadExpr,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+    )
+    new_expr = disaggregate_expression(model, affine_part.aff, bvref, method)
+    for (pair, coeff) in affine_part.terms
+        if requires_disaggregation(pair.a)
+            dref = disaggregate_expression(model, pair.a, bvref, method)
+            new_expr = JuMP.@expression(model,
+                new_expr + coeff*dref*pair.b)
+        elseif requires_disaggregation(pair.b)
+            dref = disaggregate_expression(model, pair.b, bvref, method)
+            new_expr = JuMP.@expression(model,
+                new_expr + coeff*pair.a*dref)
+        else
+            new_expr = JuMP.@expression(model,
+                new_expr + coeff*pair.a*pair.b*bvref)
+        end
+    end
+    return new_expr
+end
+
+# Disaggregate the quadratic terms without the ε-perspective division
+function _disaggregate_quad_terms(
+    model::JuMP.AbstractModel,
+    quad::JuMP.GenericQuadExpr,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+    )
+    new_expr = zero(typeof(quad))
+    for (pair, coeff) in quad.terms
+        da_ref = disaggregate_expression(model, pair.a, bvref, method)
+        db_ref = disaggregate_expression(model, pair.b, bvref, method)
+        JuMP.add_to_expression!(new_expr, coeff, da_ref, db_ref)
+    end
+    return new_expr
+end
+
+# GEHR expression (Eq. 13) for a constraint normalized to h(x) vs 0:
+# ν'Qν + (a'ν)*y + d*y²
+function _gehr_expression(
+    model::JuMP.AbstractModel,
+    h::JuMP.GenericQuadExpr,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+    )
+    quad_part, affine_part = _split_quad_terms(h)
+    quad_part = _disaggregate_quad_terms(model, quad_part, bvref, method)
+    affine_part = _disaggregate_affine_terms(model, affine_part, bvref, method)
+    return JuMP.@expression(model, quad_part + affine_part*bvref)
+end
+
+# Add the CEHR epigraph variable t >= 0 for one quadratic constraint;
+# its type is inferred from the constraint function so extensions can
+# match it to the constraint (e.g. infinite over its parameters)
+function _add_cehr_epigraph_variable(
+    model::JuMP.AbstractModel,
+    h::JuMP.GenericQuadExpr,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr}
+    )
+    base = "t_cehr_$(bvref)"
+    n = count(
+        v -> startswith(JuMP.name(v), base),
+        _reformulation_variables(model)
+    )
+    epigraph_vref = create_variable(model, VariableProperties(h))
+    JuMP.set_lower_bound(epigraph_vref, 0)
+    JuMP.set_name(epigraph_vref, n == 0 ? base : "$(base)_$(n + 1)")
+    push!(_reformulation_variables(model), epigraph_vref)
+    return epigraph_vref
+end
+
+# CEHR for solvers that accept cones: the quadratic is handed over as
+# an explicit rotated second-order cone
+function _cehr_conic_constraints(
+    model::JuMP.AbstractModel,
+    h::JuMP.GenericQuadExpr,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+    )
+    quad_part, affine_part = _split_quad_terms(h)
+    epigraph_vref = _add_cehr_epigraph_variable(model, h, bvref)
+    quad_matrix, vars = _quad_coefficient_matrix(quad_part)
+    factors = LinearAlgebra.eigen(LinearAlgebra.Symmetric(quad_matrix))
+    tol = 1e-9 * max(one(eltype(quad_matrix)), maximum(abs, quad_matrix))
+    dvars = [disaggregate_expression(model, v, bvref, method)
+             for v in vars]
+    rows = [0.5 * epigraph_vref, zero(0.5 * epigraph_vref) + bvref]
+    for k in findall(>(tol), factors.values) # drop zero eigenvalues
+        push!(rows, sqrt(factors.values[k]) * sum(
+            factors.vectors[j, k] * dvars[j] for j in eachindex(dvars)))
+    end
+    affine_part = _disaggregate_affine_terms(model, affine_part, bvref, method)
+    epigraph_link = JuMP.@expression(model, epigraph_vref + affine_part)
+    return [
+        JuMP.build_constraint(error, rows,
+            _MOI.RotatedSecondOrderCone(length(rows))),
+        JuMP.build_constraint(error, epigraph_link, _MOI.LessThan(0))
+    ]
+end
+
+# Reformulate a quadratic constraint normalized to h(x) ≤ 0 exactly:
+# CEHR when the quadratic part is convex (unless `:gehr` is forced),
+# GEHR otherwise. CEHR stays algebraic for solvers that do not accept
+# cones; `:cehr_conic` writes the cone out for solvers that do.
+function _exact_quad_hull(
+    model::JuMP.AbstractModel,
+    h::JuMP.GenericQuadExpr,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+    )
+    quad_part, affine_part = _split_quad_terms(h)
+    if isempty(quad_part.terms) # affine in the disaggregated variables
+        affine_part = _disaggregate_affine_terms(
+            model, affine_part, bvref, method)
+        return [JuMP.build_constraint(error, affine_part, _MOI.LessThan(0))]
+    elseif _is_convex_quad(quad_part) && method.quadratic == :cehr_conic
+        return _cehr_conic_constraints(model, h, bvref, method)
+    elseif _is_convex_quad(quad_part) && method.quadratic != :gehr
+        epigraph_vref = _add_cehr_epigraph_variable(model, h, bvref)
+        quad_part = _disaggregate_quad_terms(model, quad_part, bvref, method)
+        affine_part = _disaggregate_affine_terms(
+            model, affine_part, bvref, method)
+        cone_func = JuMP.@expression(model, quad_part - epigraph_vref*bvref)
+        epigraph_link = JuMP.@expression(model, epigraph_vref + affine_part)
+        return [
+            JuMP.build_constraint(error, cone_func, _MOI.LessThan(0)),
+            JuMP.build_constraint(error, epigraph_link, _MOI.LessThan(0))
+        ]
+    elseif method.quadratic in (:cehr, :cehr_conic)
+        error("`Hull(quadratic = :$(method.quadratic))` requires convex " *
+              "quadratic disjunct constraints. Use `quadratic = :exact` " *
+              "or `quadratic = :gehr` for nonconvex quadratic constraints.")
+    else
+        gehr_func = _gehr_expression(model, h, bvref, method)
+        return [JuMP.build_constraint(error, gehr_func, _MOI.LessThan(0))]
+    end
+end
+
+# Reformulate a quadratic equality normalized to h(x) = 0 exactly
+# (GEHR; CEHR does not apply since quadratic equalities are nonconvex)
+function _exact_quad_hull_eq(
+    model::JuMP.AbstractModel,
+    h::JuMP.GenericQuadExpr,
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+    )
+    quad_part, affine_part = _split_quad_terms(h)
+    if isempty(quad_part.terms) # affine in the disaggregated variables
+        affine_part = _disaggregate_affine_terms(
+            model, affine_part, bvref, method)
+        return [JuMP.build_constraint(error, affine_part, _MOI.EqualTo(0))]
+    elseif method.quadratic in (:cehr, :cehr_conic)
+        error("`Hull(quadratic = :$(method.quadratic))` does not support " *
+              "quadratic equality constraints. Use `quadratic = :exact` " *
+              "or `quadratic = :gehr` instead.")
+    end
+    gehr_func = _gehr_expression(model, h, bvref, method)
+    return [JuMP.build_constraint(error, gehr_func, _MOI.EqualTo(0))]
+end
+
+# scalar quadratic constraint
+function reformulate_disjunct_constraint(
+    model::JuMP.AbstractModel,
+    con::JuMP.ScalarConstraint{T, S},
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+) where {
+    T <: JuMP.GenericQuadExpr,
+    S <: Union{_MOI.LessThan, _MOI.GreaterThan, _MOI.EqualTo}
+}
+    set_value = _set_value(con.set)
+    if method.quadratic == :epsilon # ε-approximated perspective
+        new_func = disaggregate_expression(model, con.func, bvref, method)
+        new_func -= set_value*bvref
+        return [JuMP.build_constraint(error, new_func, S(0))]
+    elseif S <: _MOI.EqualTo
+        return _exact_quad_hull_eq(model, con.func - set_value, bvref, method)
+    else # normalize to h(x) ≤ 0 (flip GreaterThan constraints)
+        h = S <: _MOI.LessThan ? con.func - set_value : set_value - con.func
+        return _exact_quad_hull(model, h, bvref, method)
+    end
+end
+# scalar quadratic constraint in an interval
+function reformulate_disjunct_constraint(
+    model::JuMP.AbstractModel,
+    con::JuMP.ScalarConstraint{T, S},
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+) where {T <: JuMP.GenericQuadExpr, S <: _MOI.Interval}
+    if method.quadratic == :epsilon # ε-approximated perspective
+        new_func = disaggregate_expression(model, con.func, bvref, method)
+        new_func_gt = JuMP.@expression(model, new_func - con.set.lower*bvref)
+        new_func_lt = JuMP.@expression(model, new_func - con.set.upper*bvref)
+        return [
+            JuMP.build_constraint(error, new_func_gt, _MOI.GreaterThan(0)),
+            JuMP.build_constraint(error, new_func_lt, _MOI.LessThan(0))
+        ]
+    end
+    return vcat(
+        _exact_quad_hull(model, con.func - con.set.upper, bvref, method),
+        _exact_quad_hull(model, con.set.lower - con.func, bvref, method)
+    )
+end
+# vector quadratic constraint
+function reformulate_disjunct_constraint(
+    model::JuMP.AbstractModel,
+    con::JuMP.VectorConstraint{T, S, R},
+    bvref::Union{JuMP.AbstractVariableRef, JuMP.GenericAffExpr},
+    method::_Hull
+) where {
+    T <: JuMP.GenericQuadExpr,
+    S <: Union{_MOI.Nonpositives, _MOI.Nonnegatives, _MOI.Zeros}, R
+}
+    if method.quadratic == :epsilon # ε-approximated perspective
+        new_func = _combine_expressions([
+            disaggregate_expression(model, func, bvref, method)
+            for func in con.func
+        ])
+        return [JuMP.build_constraint(error, new_func, con.set)]
+    end
+    reform_cons = Vector{AbstractConstraint}()
+    for func in con.func # normalize each entry to h(x) ≤ 0 or h(x) = 0
+        h = S <: _MOI.Nonnegatives ? -func : func
+        if S <: _MOI.Zeros
+            append!(reform_cons, _exact_quad_hull_eq(model, h, bvref, method))
+        else
+            append!(reform_cons, _exact_quad_hull(model, h, bvref, method))
+        end
+    end
+    return reform_cons
+end
diff --git a/test/constraints/bigm.jl b/test/constraints/bigm.jl
index ac1281c9..16a387ab 100644
--- a/test/constraints/bigm.jl
+++ b/test/constraints/bigm.jl
@@ -313,6 +313,67 @@ function test_extension_bigm()
     @test refcons[4].set == MOI.GreaterThan(10.0 - 110)
 end
 
+function test_soc_bigm()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, t)
+    @variable(model, y, Logical)
+    @constraint(model, con, [t, x] in SecondOrderCone(), Disjunct(y))
+    bvref = binary_variable(y)
+    ref = reformulate_disjunct_constraint(model, constraint_object(con), bvref, BigM(100, false))
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func[1], t + 100*(1 - bvref))
+    @test isequal_canonical(ref[1].func[2], 1.0*x)
+    @test ref[1].set == MOI.SecondOrderCone(2)
+end
+
+function test_rsoc_bigm()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, t)
+    @variable(model, y, Logical)
+    @constraint(model, con, [0.5, t, x] in RotatedSecondOrderCone(), Disjunct(y))
+    bvref = binary_variable(y)
+    ref = reformulate_disjunct_constraint(model, constraint_object(con), bvref, BigM(100, false))
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func[1], 0.5 + 100*(1 - bvref))
+    @test isequal_canonical(ref[1].func[2], t + 100*(1 - bvref))
+    @test isequal_canonical(ref[1].func[3], 1.0*x)
+    @test ref[1].set == MOI.RotatedSecondOrderCone(3)
+end
+
+function test_exp_bigm()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, t)
+    @variable(model, w)
+    @variable(model, y, Logical)
+    @constraint(model, con, [x, t, w] in MOI.ExponentialCone(), Disjunct(y))
+    bvref = binary_variable(y)
+    ref = reformulate_disjunct_constraint(model, constraint_object(con), bvref, BigM(100, false))
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func[1], 1.0*x)
+    @test isequal_canonical(ref[1].func[2], t + 100*(1 - bvref))
+    @test isequal_canonical(ref[1].func[3], w + 200*(1 - bvref))
+    @test ref[1].set == MOI.ExponentialCone()
+end
+
+function test_power_bigm()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, t)
+    @variable(model, w)
+    @variable(model, y, Logical)
+    @constraint(model, con, [x, t, w] in MOI.PowerCone(0.5), Disjunct(y))
+    bvref = binary_variable(y)
+    ref = reformulate_disjunct_constraint(model, constraint_object(con), bvref, BigM(100, false))
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func[1], x + 100*(1 - bvref))
+    @test isequal_canonical(ref[1].func[2], t + 100*(1 - bvref))
+    @test isequal_canonical(ref[1].func[3], 1.0*w)
+    @test ref[1].set == MOI.PowerCone(0.5)
+end
+
 @testset "BigM Reformulation" begin
     test_default_bigm()
     test_default_tighten_bigm()
@@ -333,4 +394,8 @@ end
     test_zeros_bigm()
     test_nested_bigm()
     test_extension_bigm()
+    test_soc_bigm()
+    test_rsoc_bigm()
+    test_exp_bigm()
+    test_power_bigm()
 end
\ No newline at end of file
diff --git a/test/constraints/hull.jl b/test/constraints/hull.jl
index d01b30e2..d5739a1f 100644
--- a/test/constraints/hull.jl
+++ b/test/constraints/hull.jl
@@ -652,6 +652,327 @@ function test_extension_hull()
     # TODO add more tests
 end
 
+function test_vector_soc_hull()
+    model = GDPModel()
+    @variable(model, 10 <= x <= 100)
+    @variable(model, 10 <= t <= 100)
+    @variable(model, z, Logical)
+    @constraint(model, con, [t - 5, x - 5] in SecondOrderCone(), Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(1e-3), Set([x, t]))
+    prep_bounds([x, t], model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x, t]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    t_z = variable_by_name(model, "t_z")
+    ref = reformulate_disjunct_constraint(model, constraint_object(con), zbin, method)
+    @test length(ref) == 1
+    @test ref[1].func == [t_z - 5*zbin, x_z - 5*zbin]
+    @test ref[1].set == MOI.SecondOrderCone(2)
+end
+
+function test_vector_rsoc_hull()
+    model = GDPModel()
+    @variable(model, 10 <= x <= 100)
+    @variable(model, 10 <= t <= 100)
+    @variable(model, z, Logical)
+    @constraint(model, con, [0.5, t - 2, x - 3] in RotatedSecondOrderCone(), Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(1e-3), Set([x, t]))
+    prep_bounds([x, t], model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x, t]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    t_z = variable_by_name(model, "t_z")
+    ref = reformulate_disjunct_constraint(model, constraint_object(con), zbin, method)
+    @test length(ref) == 1
+    @test ref[1].func == [0.5*zbin, t_z - 2*zbin, x_z - 3*zbin]
+    @test ref[1].set == MOI.RotatedSecondOrderCone(3)
+end
+
+function test_vector_exp_hull()
+    model = GDPModel()
+    @variable(model, 10 <= x <= 100)
+    @variable(model, 10 <= t <= 100)
+    @variable(model, 10 <= w <= 100)
+    @variable(model, z, Logical)
+    @constraint(model, con, [x - 1, t - 1, w - 1] in MOI.ExponentialCone(), Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(1e-3), Set([x, t, w]))
+    prep_bounds([x, t, w], model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x, t, w]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    t_z = variable_by_name(model, "t_z")
+    w_z = variable_by_name(model, "w_z")
+    ref = reformulate_disjunct_constraint(model, constraint_object(con), zbin, method)
+    @test length(ref) == 1
+    @test ref[1].func == [x_z - zbin, t_z - zbin, w_z - zbin]
+    @test ref[1].set == MOI.ExponentialCone()
+end
+
+function test_hull_quadratic_option_error()
+    @test_throws ErrorException Hull(quadratic = :bad_option)
+    @test Hull().quadratic == :epsilon
+    @test Hull(1e-3, quadratic = :exact).quadratic == :exact
+end
+#less than, greater than, equalto with GEHR forced
+function test_scalar_gehr_hull_1sided(moiset)
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, x^2 + 3x in moiset(5), Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :gehr), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test length(ref) == 1
+    gehr = x_z^2 + 3*x_z*zbin - 5*zbin^2
+    if moiset == MOI.GreaterThan # flipped to h(x) ≤ 0
+        @test isequal_canonical(ref[1].func, -gehr)
+        @test ref[1].set == MOI.LessThan(0.0)
+    else
+        @test isequal_canonical(ref[1].func, gehr)
+        expected = moiset == MOI.EqualTo ? MOI.EqualTo(0.0) :
+            MOI.LessThan(0.0)
+        @test ref[1].set == expected
+    end
+end
+#convex quadratic routed to CEHR under :exact
+function test_scalar_cehr_hull()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, x^2 + 3x <= 5, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :exact), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test length(ref) == 2
+    t = variable_by_name(model, "t_cehr_z")
+    @test lower_bound(t) == 0
+    @test isequal_canonical(ref[1].func, x_z^2 - t*zbin)
+    @test ref[1].set == MOI.LessThan(0.0)
+    @test isequal_canonical(ref[2].func, t + 3*x_z - 5*zbin)
+    @test ref[2].set == MOI.LessThan(0.0)
+end
+#concave greater-than flips to a convex constraint and routes to CEHR
+function test_scalar_cehr_hull_concave_greater()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, -x^2 + 3x >= 1, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :exact), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test length(ref) == 2
+    t = variable_by_name(model, "t_cehr_z")
+    @test isequal_canonical(ref[1].func, x_z^2 - t*zbin)
+    @test isequal_canonical(ref[2].func, t - 3*x_z + 1*zbin)
+end
+#cehr_conic writes the cone out as an explicit rotated SOC
+function test_scalar_cehr_conic_hull()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, x^2 + 3x <= 5, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :cehr_conic), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test length(ref) == 2
+    t = variable_by_name(model, "t_cehr_z")
+    @test lower_bound(t) == 0
+    @test ref[1].set == MOI.RotatedSecondOrderCone(3)
+    @test isequal_canonical(ref[1].func[1], 0.5 * t)
+    @test isequal_canonical(ref[1].func[2], 1.0 * zbin)
+    # eigenvector sign is implementation-defined
+    @test abs(coefficient(ref[1].func[3], x_z)) ≈ 1.0
+    @test isequal_canonical(ref[2].func, t + 3*x_z - 5*zbin)
+    @test ref[2].set == MOI.LessThan(0.0)
+end
+#singular Q drops its zero eigenvalue rows from the cone
+function test_scalar_cehr_conic_rank_deficient()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, -2 <= w <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, (x + w)^2 <= 4, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :cehr_conic), Set([x, w]))
+    prep_bounds([x, w], model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x, w]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    w_z = variable_by_name(model, "w_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test ref[1].set == MOI.RotatedSecondOrderCone(3) # rank 1, not 2
+    row = ref[1].func[3]
+    s = sign(coefficient(row, x_z))
+    @test coefficient(row, x_z) ≈ s * 1.0
+    @test coefficient(row, w_z) ≈ s * 1.0
+end
+#cehr_conic rejects nonconvex and equality constraints like cehr
+function test_scalar_cehr_conic_errors()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, -2 <= w <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, x*w <= 5, Disjunct(z))
+    @constraint(model, con_eq, x^2 == 4, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :cehr_conic), Set([x, w]))
+    prep_bounds([x, w], model, Hull())
+    DP._disaggregate_variables(model, z, Set([x, w]), method)
+    @test_throws ErrorException reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test_throws ErrorException reformulate_disjunct_constraint(
+        model, constraint_object(con_eq), zbin, method)
+end
+#all quadratic terms are pure on normal models, rest keeps the affine
+function test_split_quad_terms()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, w)
+    quad = @expression(model, x^2 + 2*x*w + 3*x + 2)
+    quad_part, affine_part = DP._split_quad_terms(quad)
+    @test isequal_canonical(quad_part, @expression(model, x^2 + 2*x*w))
+    @test isempty(affine_part.terms)
+    @test isequal_canonical(affine_part.aff, @expression(model, 3*x + 2))
+end
+#nonconvex quadratic routed to GEHR under :exact, error under :cehr
+function test_scalar_exact_hull_nonconvex()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, -2 <= w <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, x*w <= 5, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :exact), Set([x, w]))
+    prep_bounds([x, w], model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x, w]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    w_z = variable_by_name(model, "w_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func, x_z*w_z - 5*zbin^2)
+    @test ref[1].set == MOI.LessThan(0.0)
+    forced = DP._Hull(Hull(quadratic = :cehr), Set([x, w]))
+    forced.disjunct_variables = method.disjunct_variables
+    @test_throws ErrorException reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, forced)
+end
+#quadratic equality reformulated by GEHR, error under :cehr
+function test_scalar_exact_hull_equality()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, x^2 == 5, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :exact), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func, x_z^2 - 5*zbin^2)
+    @test ref[1].set == MOI.EqualTo(0.0)
+    forced = DP._Hull(Hull(quadratic = :cehr), Set([x]))
+    forced.disjunct_variables = method.disjunct_variables
+    @test_throws ErrorException reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, forced)
+end
+#interval: CEHR on the upper side, GEHR on the (nonconvex) lower side
+function test_scalar_exact_hull_2sided()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, 2 <= x^2 <= 5, Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :exact), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    @test length(ref) == 3
+    t = variable_by_name(model, "t_cehr_z")
+    @test isequal_canonical(ref[1].func, x_z^2 - t*zbin)
+    @test ref[1].set == MOI.LessThan(0.0)
+    @test isequal_canonical(ref[2].func, t - 5*zbin)
+    @test ref[2].set == MOI.LessThan(0.0)
+    @test isequal_canonical(ref[3].func, -x_z^2 + 2*zbin^2)
+    @test ref[3].set == MOI.LessThan(0.0)
+end
+#nonpositives, nonnegatives, zeros with exact reformulations
+function test_vector_exact_hull_1sided(moiset)
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, [x^2 - 5; x^2 - 5] in moiset(2), Disjunct(z))
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :exact), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), zbin, method)
+    if moiset == MOI.Nonpositives # convex entries: CEHR per entry
+        @test length(ref) == 4
+        @test all(r.set == MOI.LessThan(0.0) for r in ref)
+        tvars = filter(v -> startswith(name(v), "t_cehr"), all_variables(model))
+        @test length(tvars) == 2
+        @test all(lower_bound(t) == 0 for t in tvars)
+    elseif moiset == MOI.Nonnegatives # flipped: nonconvex, GEHR
+        @test length(ref) == 2
+        for i in 1:2
+            @test isequal_canonical(ref[i].func, -x_z^2 + 5*zbin^2)
+            @test ref[i].set == MOI.LessThan(0.0)
+        end
+    else # Zeros: equalities, GEHR
+        @test length(ref) == 2
+        for i in 1:2
+            @test isequal_canonical(ref[i].func, x_z^2 - 5*zbin^2)
+            @test ref[i].set == MOI.EqualTo(0.0)
+        end
+    end
+end
+#quadratic expression without quadratic terms falls back to affine hull
+function test_exact_hull_affine_fallback()
+    model = GDPModel()
+    @variable(model, -2 <= x <= 3)
+    @variable(model, z, Logical)
+    zbin = variable_by_name(model, "z")
+    method = DP._Hull(Hull(quadratic = :exact), Set([x]))
+    prep_bounds(x, model, Hull())
+    @test DP._disaggregate_variables(model, z, Set([x]), method) isa Nothing
+    x_z = variable_by_name(model, "x_z")
+    con = JuMP.build_constraint(
+        error, convert(QuadExpr, 2.0x + 1.0), MOI.LessThan(5.0))
+    ref = reformulate_disjunct_constraint(model, con, zbin, method)
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func, 2*x_z + 1*zbin - 5*zbin)
+    @test ref[1].set == MOI.LessThan(0.0)
+    con = JuMP.build_constraint(
+        error, convert(QuadExpr, 2.0x + 1.0), MOI.EqualTo(5.0))
+    ref = reformulate_disjunct_constraint(model, con, zbin, method)
+    @test length(ref) == 1
+    @test isequal_canonical(ref[1].func, 2*x_z + 1*zbin - 5*zbin)
+    @test ref[1].set == MOI.EqualTo(0.0)
+end
+
 @testset "Hull Reformulation" begin
     test_default_hull()
     test_set_hull()
@@ -676,6 +997,7 @@ end
         test_scalar_var_hull_1sided(s)
         test_scalar_affine_hull_1sided(s)
         test_scalar_quadratic_hull_1sided(s)
+        test_scalar_gehr_hull_1sided(s)
         test_scalar_nonlinear_hull_1sided(s)
     end
     test_scalar_nonlinear_hull_1sided_error()
@@ -683,6 +1005,7 @@ end
         test_vector_var_hull_1sided(s)
         test_vector_affine_hull_1sided(s)
         test_vector_quadratic_hull_1sided(s)
+        test_vector_exact_hull_1sided(s)
         test_vector_nonlinear_hull_1sided(s)
     end
     test_vector_nonlinear_hull_1sided_error()
@@ -691,6 +1014,20 @@ end
     test_scalar_quadratic_hull_2sided()
     test_scalar_nonlinear_hull_2sided()
     test_scalar_nonlinear_hull_2sided_error()
+    test_hull_quadratic_option_error()
+    test_split_quad_terms()
+    test_scalar_cehr_hull()
+    test_scalar_cehr_hull_concave_greater()
+    test_scalar_cehr_conic_hull()
+    test_scalar_cehr_conic_rank_deficient()
+    test_scalar_cehr_conic_errors()
+    test_scalar_exact_hull_nonconvex()
+    test_scalar_exact_hull_equality()
+    test_scalar_exact_hull_2sided()
+    test_exact_hull_affine_fallback()
     test_exactly1_error()
     test_extension_hull()
+    test_vector_soc_hull()
+    test_vector_rsoc_hull()
+    test_vector_exp_hull()
 end
diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl
index 47bb884b..9201375a 100644
--- a/test/extensions/InfiniteDisjunctiveProgramming.jl
+++ b/test/extensions/InfiniteDisjunctiveProgramming.jl
@@ -1,4 +1,4 @@
-using InfiniteOpt, HiGHS, Ipopt, Juniper
+using InfiniteOpt, HiGHS, Ipopt, Juniper, Pajarito, Hypatia
 import DisjunctiveProgramming as DP
 
 # Helper to access internal function
@@ -121,6 +121,38 @@ function test_requires_disaggregation()
     @test DP.requires_disaggregation(y) == true
 end
 
+# Parameter refs in disjunct constraints answer the standard JuMP
+# bound queries (InfiniteOpt bound support): finite parameters their
+# value, infinite parameters their domain bounds, parameter functions
+# any user-declared bounds. The base set_variable_bound_info methods
+# handle them like variables, so no extension overrides are needed.
+function test_parameter_bound_info()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t in [0, 1])
+    @finite_parameter(model, p == 2.0)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, y, Infinite(t))
+    @parameter_function(model, pf == t -> 2t - 1)
+    set_lower_bound(pf, -1)
+    set_upper_bound(pf, 1)
+    @parameter_function(model, pf2 == t -> sin(t))
+    @test DP.set_variable_bound_info(pf, BigM()) == (-1.0, 1.0)
+    @test DP.set_variable_bound_info(p, BigM()) == (2.0, 2.0)
+    @test DP.set_variable_bound_info(t, BigM()) == (0.0, 1.0)
+    @test DP.set_variable_bound_info(x, BigM()) == (0.0, 10.0)
+    @test DP.set_variable_bound_info(y, BigM()) == (-Inf, Inf)
+    # parameter functions without declared bounds fall back
+    @test DP.set_variable_bound_info(pf2, BigM()) == (-Inf, Inf)
+    # Hull and PSplit clamp the bounds to include 0
+    @test DP.set_variable_bound_info(pf, Hull()) == (-1.0, 1.0)
+    @test DP.set_variable_bound_info(t, Hull()) == (0.0, 1.0)
+    @test DP.set_variable_bound_info(p, Hull()) == (0.0, 2.0)
+    @test DP.set_variable_bound_info(x, Hull()) == (0.0, 10.0)
+    @test DP.set_variable_bound_info(p, PSplit([[x]])) == (0.0, 2.0)
+    @test_throws ErrorException DP.set_variable_bound_info(y, Hull())
+    @test_throws ErrorException DP.set_variable_bound_info(pf2, Hull())
+end
+
 function test_all_variables_infiniteopt()
     model = InfiniteGDPModel()
     @infinite_parameter(model, t ∈ [0, 1])
@@ -197,6 +229,275 @@ function test_disaggregate_expression_infiniteopt()
     @test haskey(result_not_disagg.terms, y)
 end
 
+function test_disaggregate_expression_parameter()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @finite_parameter(model, p == 2.0)
+    @parameter_function(model, pf == t -> 2t)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, w, Bin)
+    @variable(model, z, InfiniteLogical(t))
+    bvref = DP._indicator_to_binary(model)[z]
+    method = DP._Hull(Hull(), Set([x]))
+    DP._variable_bounds(model)[x] = DP.set_variable_bound_info(x, Hull())
+    DP._disaggregate_variables(model, z, Set([x]), method)
+    dvref = method.disjunct_variables[x, bvref]
+    # parameters scale with the indicator, variables map through
+    @test isequal_canonical(
+        DP.disaggregate_expression(model, t, bvref, method), t * bvref)
+    @test isequal_canonical(
+        DP.disaggregate_expression(model, p, bvref, method), p * bvref)
+    @test isequal_canonical(
+        DP.disaggregate_expression(model, pf, bvref, method), pf * bvref)
+    @test isequal(DP.disaggregate_expression(model, x, bvref, method), dvref)
+    @test isequal(DP.disaggregate_expression(model, w, bvref, method), w)
+end
+
+function test_split_quad_terms_infinite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1])
+    @parameter_function(model, pf == t -> 2t)
+    @variable(model, x, Infinite(t))
+    quad = @expression(model, x^2 + t*x + t*pf + 3x + 2)
+    quad_part, affine_part = DP._split_quad_terms(quad)
+    @test isequal_canonical(quad_part, @expression(model, x^2))
+    @test isequal_canonical(
+        affine_part, @expression(model, t*x + t*pf + 3x + 2))
+end
+
+function test_epsilon_quad_hull_parameter()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @parameter_function(model, pf == t -> 2t)
+    @variable(model, -2 <= x <= 3, Infinite(t))
+    @variable(model, z, InfiniteLogical(t))
+    @constraint(model, con, x^2 + pf*x <= 1, Disjunct(z))
+    bvref = DP._indicator_to_binary(model)[z]
+    ϵ = 1e-3
+    method = DP._Hull(Hull(ϵ), Set([x]))
+    DP._variable_bounds(model)[x] = DP.set_variable_bound_info(x, Hull())
+    DP._disaggregate_variables(model, z, Set([x]), method)
+    x_z = method.disjunct_variables[x, bvref]
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), bvref, method)
+    @test length(ref) == 1
+    # pf*ν stays undivided, ν² is ε-divided, the set value scales with y
+    vals1 = Dict(x_z => 2.0, pf => 0.5, bvref => 1.0)
+    @test _eval_at(v -> vals1[v], ref[1].func) ≈ 4.0 + 1.0 - 1.0
+    vals2 = Dict(x_z => 1.0, pf => 0.5, bvref => 0.5)
+    expected2 = 1.0 / ((1 - ϵ)*0.5 + ϵ) + 0.5 - 0.5
+    @test _eval_at(v -> vals2[v], ref[1].func) ≈ expected2
+end
+
+# evaluate an expression at a point (InfiniteOpt lacks value(f, vref))
+_eval_at(f, x::Number) = x
+_eval_at(f, v::GeneralVariableRef) = f(v)
+function _eval_at(f, e::JuMP.GenericAffExpr)
+    return e.constant + sum(c * f(v) for (v, c) in e.terms; init = 0.0)
+end
+function _eval_at(f, e::JuMP.GenericQuadExpr)
+    return _eval_at(f, e.aff) +
+        sum(c * f(p.a) * f(p.b) for (p, c) in e.terms; init = 0.0)
+end
+function _eval_at(f, e::JuMP.GenericNonlinearExpr)
+    return getfield(Base, e.head)((_eval_at(f, a) for a in e.args)...)
+end
+
+function test_gehr_infinite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @variable(model, -2 <= x <= 3, Infinite(t))
+    @variable(model, z, InfiniteLogical(t))
+    @constraint(model, con, x^2 + t*x <= 1, Disjunct(z))
+    bvref = DP._indicator_to_binary(model)[z]
+    method = DP._Hull(Hull(quadratic = :gehr), Set([x]))
+    DP._variable_bounds(model)[x] = DP.set_variable_bound_info(x, Hull())
+    DP._disaggregate_variables(model, z, Set([x]), method)
+    x_z = method.disjunct_variables[x, bvref]
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), bvref, method)
+    @test length(ref) == 1
+    @test ref[1].set == MOI.LessThan(0)
+    # GEHR must give ν² + t*ν*y - y² (parameter terms scale once with y)
+    vals1 = Dict(x_z => 2.0, t => 0.5, bvref => 1.0)
+    @test _eval_at(v -> vals1[v], ref[1].func) ≈ 4.0 + 1.0 - 1.0
+    vals2 = Dict(x_z => 1.0, t => 1.0, bvref => 0.5)
+    @test _eval_at(v -> vals2[v], ref[1].func) ≈ 1.0 + 0.5 - 0.25
+end
+
+function test_cehr_infinite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @parameter_function(model, pf == t -> 2t)
+    @variable(model, -2 <= x <= 3, Infinite(t))
+    @variable(model, z, InfiniteLogical(t))
+    @constraint(model, con, x^2 + 2x - pf <= 0, Disjunct(z))
+    bvref = DP._indicator_to_binary(model)[z]
+    method = DP._Hull(Hull(quadratic = :exact), Set([x]))
+    DP._variable_bounds(model)[x] = DP.set_variable_bound_info(x, Hull())
+    DP._disaggregate_variables(model, z, Set([x]), method)
+    x_z = method.disjunct_variables[x, bvref]
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), bvref, method)
+    @test length(ref) == 2
+    tvars = filter(v -> startswith(name(v), "t_cehr"),
+                   DP._reformulation_variables(model))
+    @test length(tvars) == 1
+    tvref = only(tvars)
+    # the epigraph variable is infinite over the constraint's parameter
+    @test isequal(parameter_refs(tvref), (t,))
+    @test lower_bound(tvref) == 0
+    @test isequal_canonical(ref[1].func, x_z^2 - tvref*bvref)
+    @test ref[1].set == MOI.LessThan(0.0)
+    @test isequal_canonical(ref[2].func, tvref + 2*x_z - pf*bvref)
+    @test ref[2].set == MOI.LessThan(0.0)
+end
+
+function test_cehr_conic_infinite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @parameter_function(model, pf == t -> 2t)
+    @variable(model, -2 <= x <= 3, Infinite(t))
+    @variable(model, z, InfiniteLogical(t))
+    @constraint(model, con, x^2 + 2x - pf <= 0, Disjunct(z))
+    bvref = DP._indicator_to_binary(model)[z]
+    method = DP._Hull(Hull(quadratic = :cehr_conic), Set([x]))
+    DP._variable_bounds(model)[x] = DP.set_variable_bound_info(x, Hull())
+    DP._disaggregate_variables(model, z, Set([x]), method)
+    x_z = method.disjunct_variables[x, bvref]
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), bvref, method)
+    @test length(ref) == 2
+    tvref = only(filter(v -> startswith(name(v), "t_cehr"),
+                        DP._reformulation_variables(model)))
+    @test isequal(parameter_refs(tvref), (t,))
+    @test ref[1].set == MOI.RotatedSecondOrderCone(3)
+    @test isequal_canonical(ref[1].func[1], 0.5 * tvref)
+    @test isequal_canonical(ref[1].func[2], 1.0 * bvref)
+    @test abs(coefficient(ref[1].func[3], x_z)) ≈ 1.0
+    @test isequal_canonical(ref[2].func, tvref + 2*x_z - pf*bvref)
+    @test ref[2].set == MOI.LessThan(0.0)
+end
+
+function test_cehr_epigraph_finite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @variable(model, -2 <= v <= 3)
+    @variable(model, z, Logical)
+    @constraint(model, con, v^2 + 2v <= 1, Disjunct(z))
+    bvref = DP._indicator_to_binary(model)[z]
+    method = DP._Hull(Hull(quadratic = :exact), Set([v]))
+    DP._variable_bounds(model)[v] = DP.set_variable_bound_info(v, Hull())
+    DP._disaggregate_variables(model, z, Set([v]), method)
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), bvref, method)
+    @test length(ref) == 2
+    tvref = only(filter(v -> startswith(name(v), "t_cehr"),
+                        DP._reformulation_variables(model)))
+    # a purely finite constraint gets a finite epigraph variable
+    @test isempty(parameter_refs(tvref))
+end
+
+function test_conic_hull_infinite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @parameter_function(model, pf == t -> 5t)
+    @variable(model, -10 <= x <= 10, Infinite(t))
+    @variable(model, -10 <= w <= 10, Infinite(t))
+    @variable(model, z, InfiniteLogical(t))
+    @constraint(model, con, [1, x - pf, w] in SecondOrderCone(), Disjunct(z))
+    @constraint(model, con2, [pf, x, w] in SecondOrderCone(), Disjunct(z))
+    bvref = DP._indicator_to_binary(model)[z]
+    method = DP._Hull(Hull(), Set([x, w]))
+    for v in (x, w)
+        DP._variable_bounds(model)[v] = DP.set_variable_bound_info(v, Hull())
+    end
+    DP._disaggregate_variables(model, z, Set([x, w]), method)
+    x_z = method.disjunct_variables[x, bvref]
+    w_z = method.disjunct_variables[w, bvref]
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), bvref, method)
+    @test length(ref) == 1
+    @test ref[1].set == MOI.SecondOrderCone(3)
+    T = eltype(ref[1].func)
+    @test isequal_canonical(ref[1].func[1], convert(T, 1 * bvref))
+    @test isequal_canonical(ref[1].func[2], convert(T, x_z - pf*bvref))
+    @test isequal_canonical(ref[1].func[3], convert(T, 1 * w_z))
+    # bare parameter rows scale with the indicator too
+    ref2 = reformulate_disjunct_constraint(
+        model, constraint_object(con2), bvref, method)
+    T2 = eltype(ref2[1].func)
+    @test isequal_canonical(ref2[1].func[1], convert(T2, pf * bvref))
+    @test isequal_canonical(ref2[1].func[2], convert(T2, 1 * x_z))
+end
+
+function test_conic_all_cones_infinite()
+    cones = [
+        (x -> [1, x[1], x[2]], SecondOrderCone(), MOI.SecondOrderCone(3)),
+        (x -> [1, x[1], x[2]], RotatedSecondOrderCone(),
+         MOI.RotatedSecondOrderCone(3)),
+        (x -> [x[1], 1, x[2]], MOI.ExponentialCone(), MOI.ExponentialCone()),
+        (x -> [x[1], 1, x[2]], MOI.PowerCone(0.5), MOI.PowerCone(0.5)),
+    ]
+    for (rows, set, moi_set) in cones, meth in (Hull(), BigM(100))
+        model = InfiniteGDPModel()
+        @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+        @variable(model, -10 <= x <= 10, Infinite(t))
+        @variable(model, -10 <= w <= 10, Infinite(t))
+        @parameter_function(model, pf == t -> 1 + t)
+        set_lower_bound(pf, 1)
+        set_upper_bound(pf, 2)
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @constraint(model, rows([x - pf, w]) in set, Disjunct(Y[1]))
+        @constraint(model, rows([x, w]) in set, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        DP.reformulate_model(model, meth)
+        sets = [constraint_object(c).set
+                for c in DP._reformulation_constraints(model)]
+        @test count(==(moi_set), sets) == 2
+        InfiniteOpt.build_transformation_backend!(model)
+        @test num_variables(InfiniteOpt.transformation_model(model)) > 0
+    end
+end
+
+function test_conic_constant_row_infinite()
+    # a constant row in a zero-direction slot keeps no variables under
+    # BigM and transcribes to a number; InfiniteOpt must promote the
+    # row vector to a common type (fork transcription fix)
+    for meth in (Hull(), BigM(100))
+        model = InfiniteGDPModel()
+        @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+        @variable(model, -10 <= x <= 10, Infinite(t))
+        @variable(model, 0 <= u <= 10, Infinite(t))
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @constraint(model, [u, x, 3] in SecondOrderCone(), Disjunct(Y[1]))
+        @constraint(model, [u, x - 5, 3] in SecondOrderCone(),
+                    Disjunct(Y[2]))
+        @disjunction(model, Y)
+        DP.reformulate_model(model, meth)
+        InfiniteOpt.build_transformation_backend!(model)
+        @test num_variables(InfiniteOpt.transformation_model(model)) > 0
+    end
+end
+
+function test_conic_bigm_infinite()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @parameter_function(model, pf == t -> 5t)
+    @variable(model, -10 <= x <= 10, Infinite(t))
+    @variable(model, -10 <= w <= 10, Infinite(t))
+    @variable(model, z, InfiniteLogical(t))
+    @constraint(model, con, [1, x - pf, w] in SecondOrderCone(), Disjunct(z))
+    bvref = DP._indicator_to_binary(model)[z]
+    ref = reformulate_disjunct_constraint(
+        model, constraint_object(con), bvref, BigM(100))
+    @test length(ref) == 1
+    @test ref[1].set == MOI.SecondOrderCone(3)
+    @test isequal_canonical(ref[1].func[1], 1 + 100*(1 - bvref))
+    @test isequal_canonical(ref[1].func[2], x - pf + 0*bvref)
+    @test isequal_canonical(ref[1].func[3], 1.0*w + 0*bvref)
+end
+
 function test_variable_properties_infiniteopt()
     model = InfiniteGDPModel()
     @infinite_parameter(model, t ∈ [0, 1])
@@ -643,6 +944,136 @@ function test_methods()
     @test value(z) ≈ expected_z atol=tol
 end
 
+function test_conic_methods_infinite()
+    # Infinite version of the circle disjunction from test/solve.jl:
+    # unit disk at the origin or at (pf(t), 0) with pf = 5t; minimizing
+    # the integral of x picks the origin disk with x(t) = -1.
+    oa = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true)
+    cs = optimizer_with_attributes(Hypatia.Optimizer, MOI.Silent() => true)
+    paj = optimizer_with_attributes(Pajarito.Optimizer,
+        "oa_solver" => oa, "conic_solver" => cs, "verbose" => false)
+    for meth in (BigM(100), Hull())
+        model = InfiniteGDPModel(paj)
+        set_attribute(model, MOI.Silent(), true)
+        @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+        @variable(model, -10 <= x <= 10, Infinite(t))
+        @variable(model, -10 <= w <= 10, Infinite(t))
+        @parameter_function(model, pf == t -> 5t)
+        set_lower_bound(pf, 0)
+        set_upper_bound(pf, 5)
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @objective(model, Min, ∫(x, t))
+        @constraint(model, [1, x, w] in SecondOrderCone(), Disjunct(Y[1]))
+        @constraint(model, [1, x - pf, w] in SecondOrderCone(), Disjunct(Y[2]))
+        @disjunction(model, Y)
+        @test optimize!(model, gdp_method = meth) isa Nothing
+        @test termination_status(model) == MOI.OPTIMAL
+        @test objective_value(model) ≈ -1 atol = 1e-3
+        @test all(isapprox.(value(x), -1, atol = 1e-3))
+    end
+end
+
+function test_exact_quad_methods_infinite()
+    # Disk at (pf(t), 0) with pf = 4t - 2, or the unit disk at the
+    # origin: the optimal x(t) = min(4t - 3, -1) switches disjunct at
+    # t = 0.5. All methods must agree on the objective.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt)
+    objectives = Float64[]
+    for meth in (BigM(100), Hull(), Hull(quadratic = :exact),
+                 Hull(quadratic = :gehr))
+        model = InfiniteGDPModel(juniper)
+        set_attribute(model, MOI.Silent(), true)
+        @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+        @variable(model, -10 <= x <= 10, Infinite(t))
+        @variable(model, -10 <= w <= 10, Infinite(t))
+        @parameter_function(model, pf == t -> 4t - 2)
+        set_lower_bound(pf, -2)
+        set_upper_bound(pf, 2)
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @objective(model, Min, ∫(x, t))
+        @constraint(model, (x - pf)^2 + w^2 <= 1, Disjunct(Y[1]))
+        @constraint(model, x^2 + w^2 <= 1, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        @test optimize!(model, gdp_method = meth) isa Nothing
+        @test termination_status(model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+        push!(objectives, objective_value(model))
+        @test value(x) ≈ [-3.0, -1.0, -1.0] atol = 1e-3
+    end
+    @test all(isapprox.(objectives, objectives[1], atol = 1e-3))
+end
+
+function test_cehr_conic_methods_infinite()
+    # Same switching geometry without a pf² term so the CEHR link stays
+    # quadratic and the transcription is a true MISOCP for Pajarito:
+    # x² - 2pf⋅x + w² <= 1 is the disk at (pf, 0) with radius √(1+pf²).
+    oa = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true)
+    cs = optimizer_with_attributes(Hypatia.Optimizer, MOI.Silent() => true)
+    paj = optimizer_with_attributes(Pajarito.Optimizer,
+        "oa_solver" => oa, "conic_solver" => cs, "verbose" => false)
+    model = InfiniteGDPModel(paj)
+    set_attribute(model, MOI.Silent(), true)
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @variable(model, -10 <= x <= 10, Infinite(t))
+    @variable(model, -10 <= w <= 10, Infinite(t))
+    @parameter_function(model, pf == t -> 4t - 2)
+    set_lower_bound(pf, -2)
+    set_upper_bound(pf, 2)
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @objective(model, Min, ∫(x, t))
+    @constraint(model, x^2 - 2*pf*x + w^2 <= 1, Disjunct(Y[1]))
+    @constraint(model, x^2 + w^2 <= 1, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @test optimize!(model, gdp_method = Hull(quadratic = :cehr_conic)) isa
+        Nothing
+    @test termination_status(model) == MOI.OPTIMAL
+    expected = [min(p - sqrt(1 + p^2), -1.0) for p in (-2.0, 0.0, 2.0)]
+    @test value(x) ≈ expected atol = 1e-3
+end
+
+function test_exact_quad_ground_truth()
+    # Reformulate-then-transcribe must match transcribing first: fix pf
+    # at each support and solve the finite GDP with the same method.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt)
+    supports = [0.0, 0.5, 1.0]
+    for meth in (Hull(quadratic = :exact), Hull(quadratic = :gehr))
+        model = InfiniteGDPModel(juniper)
+        set_attribute(model, MOI.Silent(), true)
+        @infinite_parameter(model, t ∈ [0, 1], supports = supports)
+        @variable(model, -10 <= x <= 10, Infinite(t))
+        @variable(model, -10 <= w <= 10, Infinite(t))
+        @parameter_function(model, pf == t -> 4t - 2)
+        set_lower_bound(pf, -2)
+        set_upper_bound(pf, 2)
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @objective(model, Min, ∫(x, t))
+        @constraint(model, (x - pf)^2 + w^2 <= 1, Disjunct(Y[1]))
+        @constraint(model, x^2 + w^2 <= 1, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        optimize!(model, gdp_method = meth)
+        xvals = value(x)
+        for (k, s) in enumerate(supports)
+            m = GDPModel(juniper)
+            set_attribute(m, MOI.Silent(), true)
+            pfv = 4s - 2
+            @variable(m, -10 <= xs <= 10)
+            @variable(m, -10 <= ws <= 10)
+            @variable(m, Ys[1:2], Logical)
+            @objective(m, Min, xs)
+            @constraint(m, (xs - pfv)^2 + ws^2 <= 1, Disjunct(Ys[1]))
+            @constraint(m, xs^2 + ws^2 <= 1, Disjunct(Ys[2]))
+            @disjunction(m, Ys)
+            optimize!(m, gdp_method = meth)
+            @test xvals[k] ≈ value(xs) atol = 1e-3
+        end
+    end
+end
+
 function test_mbm_with_derivatives()
     model = InfiniteGDPModel(HiGHS.Optimizer)
     set_silent(model)
@@ -806,6 +1237,7 @@ end
         test_is_parameter()
         test_is_parameter_concrete_dispatches()
         test_requires_disaggregation()
+        test_parameter_bound_info()
         test_variable_properties_infiniteopt()
         test_variable_properties_from_expr()
         test_variable_properties_from_quad_expr()
@@ -828,6 +1260,23 @@ end
     @testset "Methods" begin
         test_get_constant()
         test_disaggregate_expression_infiniteopt()
+        test_disaggregate_expression_parameter()
+    end
+
+    @testset "Conic" begin
+        test_conic_hull_infinite()
+        test_conic_bigm_infinite()
+        test_conic_all_cones_infinite()
+        test_conic_constant_row_infinite()
+    end
+
+    @testset "Exact Quadratic Hull" begin
+        test_split_quad_terms_infinite()
+        test_epsilon_quad_hull_parameter()
+        test_gehr_infinite()
+        test_cehr_infinite()
+        test_cehr_conic_infinite()
+        test_cehr_epigraph_finite()
     end
 
     @testset "MBM" begin
@@ -844,6 +1293,10 @@ end
     @testset "Integration" begin
         test_infiniteopt_extension()
         test_methods()
+        test_conic_methods_infinite()
+        test_exact_quad_methods_infinite()
+        test_cehr_conic_methods_infinite()
+        test_exact_quad_ground_truth()
     end
 
     @testset "Cutting Planes" begin
diff --git a/test/solve.jl b/test/solve.jl
index c3ca9441..61bd640c 100644
--- a/test/solve.jl
+++ b/test/solve.jl
@@ -1,4 +1,4 @@
-using HiGHS, Ipopt, Juniper
+using HiGHS, Ipopt, Juniper, Pajarito, Hypatia
 function test_linear_gdp_example(m, use_complements = false)
     set_attribute(m, MOI.Silent(), true)
     @variable(m, 1 ≤ x[1:2] ≤ 9)
@@ -168,3 +168,96 @@ end
     test_quadratic_gdp_example()
     test_generic_model(GDPModel{Float32}(mockoptimizer))
 end
+
+function test_conic_gdp_example()
+    # Circle-containment idea from Bernal Neira & Grossmann (2021),
+    # Eq. 4.3: the point (x, y) must lie in the unit circle around
+    # (0, 0) OR around (5, 0), expressed as second-order cones. The
+    # disjunction picks one circle; minimizing x selects the first and
+    # lands on its leftmost point, (x, y) = (-1, 0).
+    oa = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true)
+    cs = optimizer_with_attributes(Hypatia.Optimizer, MOI.Silent() => true)
+    paj = optimizer_with_attributes(Pajarito.Optimizer,
+        "oa_solver" => oa, "conic_solver" => cs, "verbose" => false)
+    for meth in (BigM(100), Hull())
+        m = GDPModel(paj)
+        set_attribute(m, MOI.Silent(), true)
+        @variable(m, -10 <= x <= 10)
+        @variable(m, -10 <= y <= 10)
+        @variable(m, Y[1:2], Logical)
+        @objective(m, Min, x)
+        @constraint(m, c1, [1, x, y] in SecondOrderCone(), Disjunct(Y[1]))
+        @constraint(m, c2, [1, x - 5, y] in SecondOrderCone(), Disjunct(Y[2]))
+        @disjunction(m, Y)
+        @test optimize!(m, gdp_method = meth) isa Nothing
+        @test termination_status(m) == MOI.OPTIMAL
+        @test isapprox(objective_value(m), -1, atol = 1e-4)
+        @test isapprox(value(x), -1, atol = 1e-4)
+        @test isapprox(value(y), 0, atol = 1e-4)
+        @test value(Y[1])
+        @test !value(Y[2])
+    end
+end
+
+@testset "Solve Conic GDP" begin
+    test_conic_gdp_example()
+end
+
+function test_exact_quadratic_gdp_example()
+    # Quadratic version of the circle disjunction above: the point must
+    # lie in the unit disk around (0, 0) or around (5, 0); minimizing x
+    # selects the first disk at (x, y) = (-1, 0). The exact hull
+    # reformulations (CEHR/GEHR, Gusev & Bernal Neira 2025) must agree
+    # with BigM and the ε-approximated hull.
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,
+        "print_level" => 0, "sb" => "yes")
+    juniper = optimizer_with_attributes(Juniper.Optimizer,
+        "nl_solver" => ipopt)
+    for meth in (BigM(100), Hull(), Hull(quadratic = :exact),
+                 Hull(quadratic = :gehr))
+        m = GDPModel(juniper)
+        set_attribute(m, MOI.Silent(), true)
+        @variable(m, -10 <= x <= 10)
+        @variable(m, -10 <= y <= 10)
+        @variable(m, Y[1:2], Logical)
+        @objective(m, Min, x)
+        @constraint(m, x^2 + y^2 <= 1, Disjunct(Y[1]))
+        @constraint(m, (x - 5)^2 + y^2 <= 1, Disjunct(Y[2]))
+        @disjunction(m, Y)
+        @test optimize!(m, gdp_method = meth) isa Nothing
+        @test termination_status(m) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+        @test isapprox(objective_value(m), -1, atol = 1e-4)
+        @test isapprox(value(x), -1, atol = 1e-4)
+        @test value(Y[1])
+        @test !value(Y[2])
+    end
+end
+
+function test_cehr_conic_gdp_example()
+    # Same disks as above, but reformulated to explicit rotated SOCs
+    # (quadratic = :cehr_conic) and solved as a true MICP.
+    oa = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true)
+    cs = optimizer_with_attributes(Hypatia.Optimizer, MOI.Silent() => true)
+    paj = optimizer_with_attributes(Pajarito.Optimizer,
+        "oa_solver" => oa, "conic_solver" => cs, "verbose" => false)
+    m = GDPModel(paj)
+    set_attribute(m, MOI.Silent(), true)
+    @variable(m, -10 <= x <= 10)
+    @variable(m, -10 <= y <= 10)
+    @variable(m, Y[1:2], Logical)
+    @objective(m, Min, x)
+    @constraint(m, x^2 + y^2 <= 1, Disjunct(Y[1]))
+    @constraint(m, (x - 5)^2 + y^2 <= 1, Disjunct(Y[2]))
+    @disjunction(m, Y)
+    @test optimize!(m, gdp_method = Hull(quadratic = :cehr_conic)) isa Nothing
+    @test termination_status(m) == MOI.OPTIMAL
+    @test isapprox(objective_value(m), -1, atol = 1e-4)
+    @test isapprox(value(x), -1, atol = 1e-4)
+    @test value(Y[1])
+    @test !value(Y[2])
+end
+
+@testset "Solve Exact Quadratic GDP" begin
+    test_exact_quadratic_gdp_example()
+    test_cehr_conic_gdp_example()
+end