Adding (some) conic constraints.

Project.toml

@@ -21,13 +21,17 @@ 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"]

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/hull.jl

@@ -281,6 +281,22 @@ function reformulate_disjunct_constraint(
     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 = JuMP.@expression(model, [i=1:_MOI.dimension(con.set)],
+        disaggregate_expression(model, con.func[i], bvref, method)
+    )
+    reform_con = JuMP.build_constraint(error, new_func, con.set)
+    return [reform_con]
+end
 function reformulate_disjunct_constraint(
     model::JuMP.AbstractModel, 
     con::JuMP.ScalarConstraint{T, S}, 

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

test/constraints/hull.jl

@@ -662,6 +662,62 @@ 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
+
 @testset "Hull Reformulation" begin
     test_default_hull()
     test_set_hull()
@@ -703,4 +759,7 @@ end
     test_scalar_nonlinear_hull_2sided_error()
     test_exactly1_error()
     test_extension_hull()
+    test_vector_soc_hull()
+    test_vector_rsoc_hull()
+    test_vector_exp_hull()
 end

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,37 @@ 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