Adding (some) conic constraints.

[dnguyen227]

Adding some conic constraints that can be handled with BigM and Hull. Some things I found odd that maybe you can comment on or correct me.

1. Hull reformulation is standard throughout all cones, but cones in MOI inherit from AbstractVectorSet, so I can't just ask for any cone. Right now I'm stuck naming each type of cone rather than a generic MOI.Cone datatype.

2. BigM reformulation looks like it needs to be adapted per cone type, and right now I'm kind of having to reason out the correct direction based on the table provided in the paper. Snippet from the paper below.

Overall a good introduction to cones but I just can't get past these two things while I was trying to generalize the code.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.54%. Comparing base (a4ebe8a) to head (6e87fa2).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #136      +/-   ##
==========================================
+ Coverage   99.51%   99.54%   +0.02%     
==========================================
  Files          17       17              
  Lines        2061     2183     +122     
==========================================
+ Hits         2051     2173     +122     
  Misses         10       10              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[pulsipher]

Ya, we'll have to simply enumerate the conic sets. You could could make your own internal Union type though:

_ConicSets =Union{ExponentialCone, GeometricMeanCone, ...}

[pulsipher]

Thoughts @bernalde?

[bernalde]

Folks, this is exciting. I’ll do a more in depth review once you consider it useful, but one observation from my paper was the fact that the same cone appears in the BigM and HR reformulations, just involving the binary variables. If you already have access to the type of cone and how it was built, then you have everything to complete the recipe. Happy to talk more about this and provide code that we used (the repo in the paper should have all). Cheers!

[dnguyen227]

I had a look at using MOI for the hessian evaluation and in my opinion because it's quadratic we actually have simpler version since I can read off the coefficents.

What we do now. I can just read off whats given.

Q = zeros(n, n)
for (pair, c) in quad.terms
    i, j = index[pair.a], index[pair.b]
    i == j ? (Q[i, i] += c) : (Q[i, j] += c / 2; Q[j, i] += c / 2)
end
eigmin(Symmetric(Q)) >= -tol

For MOI, I end up having to build something, then evaluate (this version works locally, so I can swap to this if we'd prefer).

nl = MOI.Nonlinear.Model()
MOI.Nonlinear.add_constraint(nl, expr, MOI.LessThan(0.0))  # expr rebuilt in MOI's index space
ev = MOI.Nonlinear.Evaluator(nl, MOI.Nonlinear.SparseReverseMode(),  MOI.VariableIndex.(1:n))
MOI.initialize(ev, [:Hess])
rows = MOI.hessian_lagrangian_structure(ev)
H = zeros(length(rows))
MOI.eval_hessian_lagrangian(ev, H, zeros(n), 0.0, [1.0])
Q = zeros(n, n)
for (k, (i, j)) in enumerate(rows)
    Q[i, j] += H[k] / 2
    i == j || (Q[j, i] += H[k] / 2)
end
eigmin(Symmetric(Q)) >= -tol   # same test either way

[pulsipher]

I had a look at using MOI for the hessian evaluation and in my opinion because it's quadratic we actually have simpler version since I can read off the coefficents.

What we do now. I can just read off whats given.

Q = zeros(n, n)
for (pair, c) in quad.terms
    i, j = index[pair.a], index[pair.b]
    i == j ? (Q[i, i] += c) : (Q[i, j] += c / 2; Q[j, i] += c / 2)
end
eigmin(Symmetric(Q)) >= -tol

For MOI, I end up having to build something, then evaluate (this version works locally, so I can swap to this if we'd prefer).

nl = MOI.Nonlinear.Model()
MOI.Nonlinear.add_constraint(nl, expr, MOI.LessThan(0.0))  # expr rebuilt in MOI's index space
ev = MOI.Nonlinear.Evaluator(nl, MOI.Nonlinear.SparseReverseMode(),  MOI.VariableIndex.(1:n))
MOI.initialize(ev, [:Hess])
rows = MOI.hessian_lagrangian_structure(ev)
H = zeros(length(rows))
MOI.eval_hessian_lagrangian(ev, H, zeros(n), 0.0, [1.0])
Q = zeros(n, n)
for (k, (i, j)) in enumerate(rows)
    Q[i, j] += H[k] / 2
    i == j || (Q[j, i] += H[k] / 2)
end
eigmin(Symmetric(Q)) >= -tol   # same test either way

I would say that it is preferred to directly compute it if you can easily do so based on the type. If it made since, we could make the AD approach a fallback, but I believe you have to enumerate each conic type anyway in MOI.

[dnguyen227]

@pulsipher ready for review

[pulsipher]

Looks good to me, I just would like to see the docs reflect that these capabilities have been added.

@bernalde: Do you approve?

[bernalde]

Review of 6e87fa2 (7 commits on top of a4ebe8a; local diff file set matches GitHub's 11 files).

The mathematics checks out. The conic hull A·ν - b·y ∈ K is exact by homogeneity of the cone; GEHR (ν'Qν + (a'ν)y + d·y² ≤ 0) and CEHR (ν'Qν ≤ t·y, t + a'ν + d·y ≤ 0, t ≥ 0) reduce to the perspective for y ∈ (0, 1] and to ν = t = 0 at y = 0; the :cehr_conic rows (t/2, y, √λ_k v_k'ν) reproduce ν'Qν ≤ t·y with zero eigenvalues dropped; and the four Big-M interior directions produce no violations at M = 100 on a [-10, 10]³ grid. The docs Josh asked for are in at this head (README, docs/src/index.md, Hull docstring).

The one blocking item is on the InfiniteGDPModel side: the CEHR epigraph variable is finite even when the constraint is infinite, which over-restricts the model (counterexample inline). Details are in the inline comments.

Findings

Inline (6):

• Blocking — src/hull.jl: finite CEHR epigraph variable on infinite models cuts off feasible points.
• Nonblocking — src/bigm.jl: conic Big-M skips the finite-M guard.
• Nonblocking — src/bigm.jl: Float64 cone directions promote coefficients on Float32 models.
• Nonblocking — README.md: document that M is used as given for cones and that the curved cones need M above the variable range.
• Nonblocking — src/constraints.jl: Indicator() on a conic disjunct fails with a raw MethodError.
• Nonblocking — test/solve.jl: no solve-level test for RSOC / exp / power cones or the README example.

Review body (1):

Nonblocking — the description and commit history cover half the change.
The title and body describe conic sets, but the PR also adds the exact quadratic hull reformulations (GEHR / CEHR / :cehr_conic) behind a new public Hull(quadratic = …) keyword with five options, a LinearAlgebra dependency, and Hypatia / Pajarito test dependencies; the commit messages are Initial, cehr_conic tag, Fix tests. Please update the body to describe the quadratic option (with the Gusev & Bernal Neira reference the docs already cite) and squash-merge with a message naming both features — that text is what the release notes and git blame will point at.

Coordination with open PRs

git merge-tree reports textual conflicts with #139 and #140 (Project.toml [compat] / [targets]) and with #141 (README.md, src/constraints.jl). #141 also adds _set_support(::Union{Nonnegatives, Nonpositives, Zeros}) = true to gate which vector sets basic steps accept; after both land that whitelist rejects the conic sets added here, so whichever merges second should decide explicitly whether basic steps accept _ConicSets.

Tests run

• PR head, local: Pkg.test() at 6e87fa2 on Julia 1.10.11 / Linux — passed (32 testsets, all Pass).
• PR head, GitHub: Julia 1.10 on ubuntu / macOS / windows, docs build, codecov/patch, codecov/project — all success for 6e87fa2.
• PR head, scratch environment: cone-direction grid check; GDPModel{Float32} reformulation; Indicator / MBM / PSplit / CuttingPlanes on a conic disjunct; the infinite-model counterexample above under BigM, Hull(), :gehr, :exact.
• Merge result: merge-tree against master is clean; conflicts against #139 / #140 / #141 as listed (no tests run on merge results).

Summary

Blocking: 1. Nonblocking: 6. Questions: 0. Total: 7 (6 inline, 1 in this body).

I would not merge this until the blocking issues above are addressed.

[bernalde]

Blocking — On an InfiniteGDPModel the CEHR epigraph variable is created with a plain @variable, so it is a finite variable even when the disjunct constraint — and therefore ν'Qν - t·y ≤ 0 and t + a'ν + d·y ≤ 0 — is infinite in τ. One t then has to satisfy max_τ ν(τ)'Qν(τ)/y ≤ t ≤ min_τ(-a'ν(τ) - d·y), which is strictly stronger than the per-τ constraint, so the reformulation silently cuts off feasible points. This is reachable with optimize!(m, gdp_method = Hull(quadratic = :exact)) on any infinite model with a convex quadratic disjunct constraint (:exact routes those here), and with :cehr / :cehr_conic.

Counterexample (supports {0, 1}, Juniper + Ipopt):

m = InfiniteGDPModel(juniper)
@infinite_parameter(m, τ in [0, 1], num_supports = 2)
@variable(m, -2 <= x <= 2, Infinite(τ))
@variable(m, 0 <= w <= 10, Infinite(τ))
@variable(m, Y[1:2], Logical)
@constraint(m, x(0) == 1); @constraint(m, x(1) == 0)
@constraint(m, x^2 - 2w <= 0, Disjunct(Y[1]))   # w(τ) ≥ x(τ)²/2
@constraint(m, w >= 5, Disjunct(Y[2]))
@disjunction(m, Y)
@objective(m, Min, support_sum(w, τ))

BigM(100, false), Hull() and Hull(quadratic = :gehr) all return objective 0.5 with w = [0.5, 0.0]; Hull(quadratic = :exact) returns 1.0 with w = [0.5, 0.5]. The reformulated model shows why: x(τ)_Y[1](τ)² - t_cehr_Y[1]*Y[1] ≤ 0 with t_cehr_Y[1] a FiniteVariableIndex.

Fix options:

• minimal: make the epigraph creation a hook that the extension overrides for InfiniteOpt.InfiniteModel to error(...) in the :cehr / :cehr_conic paths, and have :exact fall back to GEHR on infinite models (GEHR needs no auxiliary variable and is exact per τ — it gives 0.5 above);
• full: in the extension, create t through the same VariableProperties / create_variable route _disaggregate_variable uses, with parameter dependence taken from InfiniteOpt.parameter_refs(h) (plus those of bvref when the indicator is an InfiniteLogical).

Either way, add a regression test in test/extensions/InfiniteDisjunctiveProgramming.jl for Hull(quadratic = :exact) on an infinite quadratic disjunct constraint; the model above is small enough for the suite.

[dnguyen227]

@pulsipher Thoughts? I'd have to investigate some more understand the nuance here. Should we error for now when someone tries to build an InfiniteGDP with cones or the GEHR/CEHR workflow?

[pulsipher]

If possible, I would like to make sure this is compatible with InfiniteOpt. I believe that should be relatively straightforward

[dnguyen227]

Will do

[bernalde]

Nonblocking — M = method.value bypasses _get_M_value. Skipping tighten is consistent with the nonlinear scalar path (which also ends at method.value), but the finiteness guard is skipped too: BigM(Inf, false) on a conic disjunct builds t - Inf y + Inf without complaint, whereas the scalar path raises A finite Big-M value must be used. Add the same isinf(M) && error(...) check here (or route through _get_M).

[bernalde]

Nonblocking — _conic_bigm_direction always returns Vector{Float64}, so M*(1 - bvref)*d[i] promotes the coefficients: on a GDPModel{Float32} with BigM(100f0, false) the stored conic reformulation constraint has GenericAffExpr{Float64, GenericVariableRef{Float32}} entries while the scalar disjunct in the same model stays GenericAffExpr{Float32, …}. Build d in JuMP.value_type(typeof(model)) (as cuttingplanes.jl does) — e.g. pass T in and use one(T) / zero(T). test_generic_model in test/solve.jl could gain a conic disjunct to pin this.

[bernalde]

Nonblocking — Two things a user picking M for a cone needs to know that this paragraph does not say:

1. For conic sets the BigM value is used as given — no bound-based tightening — so it has to be chosen by the user.
2. "Large enough" is not the variable range for the curved cones. With every entry in [-10, 10], (x, y + M, z + 2M) ∈ ExponentialCone still has violations at M = 30 (the slack enters as (y+M)·exp(x/(y+M)), so roughly M ≳ x + y - z + x²/(2(y+M)) is needed), while the SOC direction is safe as soon as M ≥ max‖x‖ - min t. All four directions are valid at M = 100 on that grid.

One sentence for each is enough here; the same caveat belongs in the BigM docstring or the _conic_bigm_direction comments.

[bernalde]

Nonblocking — Now that _ConicSets is accepted at build time, every reformulation method can receive a conic disjunct constraint. MBM and PSplit fail with clear "not supported" messages, but Indicator() dies with MethodError: no method matching _vec_to_scalar_set(::MathOptInterface.SecondOrderCone) because _vec_to_scalar_set (lines 12–14 of this file) has no conic method. Add _vec_to_scalar_set(::_ConicSets) = error("Indicator reformulation does not support conic disjunct constraints; use BigM or Hull.") (or an Indicator-specific reformulate_disjunct_constraint method for _ConicSets) so users get the guidance the docs give. CuttingPlanes hits AssertionError: dim >= 3 on the same model — that looks like a solver/cone-dimension issue rather than this PR's code, so just noting it.

[bernalde]

Nonblocking — End-to-end coverage is SOC-only. RotatedSecondOrderCone, ExponentialCone and PowerCone are checked only at the reformulation-shape level, and the README example (SOC + ExponentialCone) is never executed. Pajarito + Hypatia handle all four cones, so adding the README example (or one exp-cone disk) here would exercise each _conic_bigm_direction choice and the conic hull against a solver.