Add complementarity constraints

docs/src/apireference.md

@@ -273,6 +273,7 @@ SettingSingleVariableFunctionNotAllowed
 List of recognized functions.
 ```@docs
 AbstractFunction
+AbstractVectorFunction
 SingleVariable
 VectorOfVariables
 ScalarAffineTerm
@@ -349,6 +350,7 @@ DualPowerCone
 SOS1
 SOS2
 IndicatorSet
+Complements
 ```
 
 ### Matrix sets

src/Test/UnitTests/basic_constraint_tests.jl

@@ -64,7 +64,7 @@ const BasicConstraintTests = Dict(
     (MOI.VectorAffineFunction{Float64}, MOI.Zeros)        => ( dummy_vector_affine, 2, MOI.Zeros(2) ),
     (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) => ( dummy_vector_affine, 2, MOI.Nonpositives(2) ),
     (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => ( dummy_vector_affine, 2, MOI.Nonnegatives(2) ),
-
+    (MOI.VectorAffineFunction{Float64}, MOI.Complements) => (dummy_vector_affine, 2, MOI.Complements(1)),
     (MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone)       => ( dummy_vector_affine, 3, MOI.NormInfinityCone(3) ),
     (MOI.VectorAffineFunction{Float64}, MOI.NormOneCone)            => ( dummy_vector_affine, 3, MOI.NormOneCone(3) ),
     (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone)        => ( dummy_vector_affine, 3, MOI.SecondOrderCone(3) ),

src/Test/nlp.jl

@@ -337,3 +337,156 @@ const nlptests = Dict("hs071" => hs071_test,
                       )
 
 @moitestset nlp
+
+"""
+    test_linear_mixed_complementarity(model::MOI.ModelLike, config::TestConfig)
+
+Test the solution of the linear mixed-complementarity problem:
+`F(x) complements x`, where `F(x) = M * x .+ q` and `0 <= x <= 10`.
+"""
+function test_linear_mixed_complementarity(model::MOI.ModelLike, config::TestConfig)
+    MOI.empty!(model)
+    x = MOI.add_variables(model, 4)
+    MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.0))
+    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
+    M = Float64[0  0 -1 -1; 0  0  1 -2; 1 -1  2 -2; 1  2 -2  4]
+    q = [2; 2; -2; -6]
+    terms = MOI.VectorAffineTerm{Float64}[]
+    for i = 1:4
+        push!(
+            terms,
+            MOI.VectorAffineTerm(4 + i, MOI.ScalarAffineTerm(1.0, x[i]))
+        )
+        for j = 1:4
+            iszero(M[i, j]) && continue
+            push!(
+                terms,
+                MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(M[i, j], x[j]))
+            )
+        end
+    end
+    MOI.add_constraint(
+        model,
+        MOI.VectorAffineFunction(terms, [q; 0.0; 0.0; 0.0; 0.0]),
+        MOI.Complements(4)
+    )
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
+    if config.solve
+        MOI.optimize!(model)
+        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
+        x_val = MOI.get.(model, MOI.VariablePrimal(), x)
+        @test isapprox(
+            x_val, [2.8, 0.0, 0.8, 1.2], atol = config.atol, rtol = config.rtol
+        )
+    end
+end
+
+const mixed_complementaritytests = Dict(
+    "test_linear_mixed_complementarity" => test_linear_mixed_complementarity,
+)
+
+@moitestset mixed_complementarity
+
+"""
+    test_qp_complementarity_constraint(model::MOI.ModelLike, config::TestConfig)
+
+Test the solution of the quadratic program with complementarity constraints:
+
+```
+    min  (x0 - 5)^2 +(2 x1 + 1)^2
+    s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2
+        x2 complements(3 x0 - x1 - 3)
+        x3 complements(-x0 + 0.5 x1 + 4)
+        x4 complements(-x0 - x1 + 7)
+        x0, x1, x2, x3, x4 >= 0
+
+```
+which rewrites, with auxiliary variables
+
+```
+    min  (x0 - 5)^2 +(2 x1 + 1)^2
+    s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2  (cf1)
+        3 x0 - x1 - 3 - x5 = 0                 (cf2)
+        -x0 + 0.5 x1 + 4 - x6 = 0              (cf3)
+        -x0 - x1 + 7 - x7 = 0                  (cf4)
+        x2 complements x5
+        x3 complements x6
+        x4 complements x7
+        x0, x1, x2, x3, x4, x5, x6, x7 >= 0
+
+```
+"""
+function test_qp_complementarity_constraint(
+    model::MOI.ModelLike, config::TestConfig
+)
+    MOI.empty!(model)
+    x = MOI.add_variables(model, 8)
+    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
+    MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0))
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
+        MOI.ScalarQuadraticFunction(
+            MOI.ScalarAffineTerm.([-10.0, 4.0], x[[1, 2]]),
+            MOI.ScalarQuadraticTerm.([2.0, 8.0], x[1:2], x[1:2]),
+            26.0
+        )
+    )
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            MOI.ScalarAffineTerm.([-1.5, 2.0, 1.0, 0.5, 1.0], x[1:5]), 0.0
+        ),
+        MOI.EqualTo(2.0)
+   )
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            MOI.ScalarAffineTerm.([3.0, -1.0, -1.0], x[[1, 2, 6]]), 0.0
+        ),
+        MOI.EqualTo(3.0)
+   )
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            MOI.ScalarAffineTerm.([-1.0, 0.5, -1.0], x[[1, 2, 7]]), 0.0
+        ),
+        MOI.EqualTo(-4.0)
+   )
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+           MOI.ScalarAffineTerm.(-1.0, x[[1, 2, 8]]), 0.0
+        ),
+        MOI.EqualTo(-7.0)
+   )
+    MOI.add_constraint(
+        model,
+        MOI.VectorOfVariables([x[3], x[4], x[5], x[6], x[7], x[8]]),
+        MOI.Complements(3)
+    )
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
+    if config.solve
+        MOI.optimize!(model)
+        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
+        x_val = MOI.get.(model, MOI.VariablePrimal(), x)
+        @test isapprox(
+            x_val,
+            [1.0, 0.0, 3.5, 0.0, 0.0, 0.0, 3.0, 6.0],
+            atol = config.atol,
+            rtol = config.rtol
+        )
+        @test isapprox(
+            MOI.get(model, MOI.ObjectiveValue()),
+            17.0,
+            atol = config.atol,
+            rtol = config.rtol
+        )
+    end
+end
+
+const math_program_complementarity_constraintstests = Dict(
+    "test_qp_complementarity_constraint" => test_qp_complementarity_constraint,
+)
+
+@moitestset math_program_complementarity_constraints

src/Utilities/model.jl

@@ -1027,7 +1027,7 @@ const LessThanIndicatorSetZero{T} = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO, MOI.L
        (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval,
         MOI.Semicontinuous, MOI.Semiinteger),
        (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
-        MOI.NormInfinityCone, MOI.NormOneCone,
+        MOI.Complements, MOI.NormInfinityCone, MOI.NormOneCone,
         MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
         MOI.GeometricMeanCone, MOI.ExponentialCone, MOI.DualExponentialCone,
         MOI.PositiveSemidefiniteConeTriangle, MOI.PositiveSemidefiniteConeSquare,

src/sets.jl

@@ -669,18 +669,74 @@ end
 
 Base.:(==)(set1::IndicatorSet{A, S}, set2::IndicatorSet{A, S}) where {A, S} = set1.set == set2.set
 
+"""
+    Complements(dimension::Int)
+
+The set corresponding to a mixed complementarity constraint.
+
+Complementarity constraints should be specified with an
+[`AbstractVectorFunction`](@ref)-in-`Complements(dimension)` constraint.
+
+The dimension of the vector-valued function `F` must be `2 * dimension`. This
+defines a complementarity constraint between the scalar function `F[i]` and the
+variable in `F[i + dimension]`. Thus, `F[i + dimension]` must be interpretable
+as a single variable `x_i` (e.g., `1.0 * x + 0.0`).
+
+The mixed complementarity problem consists of finding `x_i` in the interval
+`[lb, ub]` (i.e., in the set `Interval(lb, ub)`), such that the following holds:
+
+  1. `F_i(x) == 0` if `lb_i < x_i < ub_i`
+  2. `F_i(x) >= 0` if `lb_i == x_i`
+  3. `F_i(x) <= 0` if `x_i == ub_i`
+
+Classically, the bounding set for `x_i` is `Interval(0, Inf)`, which recovers:
+`0 <= F_i(x) ⟂ x_i >= 0`, where the `⟂` operator implies `F_i(x) * x_i = 0`.
+
+### Examples
+
+The problem:
+
+    x -in- Interval(-1, 1)
+    [-4 * x - 3, x] -in- Complements(1)
+
+defines the mixed complementarity problem where the following holds:
+
+  1. `-4 * x - 3 == 0` if `-1 < x < 1`
+  2. `-4 * x - 3 >= 0` if `x == -1`
+  3. `-4 * x - 3 <= 0` if `x == 1`
+
+There are three solutions:
+
+  1. `x = -3/4` with `F(x) = 0`
+  2. `x = -1` with `F(x) = 1`
+  3. `x = 1` with `F(x) = -7`
+
+The function `F` can also be defined in terms of single variables. For example,
+the problem:
+
+    [x_3, x_4] -in- Nonnegatives(2)
+    [x_1, x_2, x_3, x_4] -in- Complements(2)
+
+defines the complementarity problem where `0 <= x_1 ⟂ x_3 >= 0` and
+`0 <= x_2 ⟂ x_4 >= 0`.
+"""
+struct Complements <: AbstractVectorSet
+    dimension::Int
+end
+
 # isbits types, nothing to copy
-function Base.copy(set::Union{Reals, Zeros, Nonnegatives, Nonpositives,
-                              GreaterThan, LessThan, EqualTo, Interval,
-                              NormInfinityCone, NormOneCone,
-                              SecondOrderCone, RotatedSecondOrderCone,
-                              GeometricMeanCone, ExponentialCone,
-                              DualExponentialCone, PowerCone, DualPowerCone,
-                              PositiveSemidefiniteConeTriangle,
-                              PositiveSemidefiniteConeSquare,
-                              LogDetConeTriangle, LogDetConeSquare,
-                              RootDetConeTriangle, RootDetConeSquare,
-                              Integer, ZeroOne, Semicontinuous, Semiinteger})
+function Base.copy(
+    set::Union{
+        Reals, Zeros, Nonnegatives, Nonpositives, GreaterThan, LessThan,
+        EqualTo, Interval, NormInfinityCone, NormOneCone, SecondOrderCone,
+        RotatedSecondOrderCone, GeometricMeanCone, ExponentialCone,
+        DualExponentialCone, PowerCone, DualPowerCone,
+        PositiveSemidefiniteConeTriangle, PositiveSemidefiniteConeSquare,
+        LogDetConeTriangle, LogDetConeSquare, RootDetConeTriangle,
+        RootDetConeSquare, Complements, Integer, ZeroOne, Semicontinuous,
+        Semiinteger
+    }
+)
     return set
 end
 Base.copy(set::S) where {S <: Union{SOS1, SOS2}} = S(copy(set.weights))

test/Test/Test.jl

@@ -1,23 +1,8 @@
 using Test
 
-@testset "Config" begin
-    include("config.jl")
-end
-@testset "Unit" begin
-    include("unit.jl")
-end
-@testset "Continuous Linear" begin
-    include("contlinear.jl")
-end
-@testset "Continuous Conic" begin
-    include("contconic.jl")
-end
-@testset "Continuous Quadratic" begin
-    include("contquadratic.jl")
-end
-@testset "Integer Linear" begin
-    include("intlinear.jl")
-end
-@testset "Integer Conic" begin
-    include("intconic.jl")
+@testset "$(file)" for file in readdir(@__DIR__)
+    if file == "Test.jl"
+        continue
+    end
+    include(file)
 end

test/Test/nlp.jl

@@ -0,0 +1,32 @@
+using Test
+
+import MathOptInterface
+const MOI = MathOptInterface
+
+@testset "mixed_complementarity" begin
+    mock = MOI.Utilities.MockOptimizer(
+        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
+    )
+    config = MOI.Test.TestConfig(optimal_status = MOI.LOCALLY_SOLVED)
+    MOI.Utilities.set_mock_optimize!(
+        mock,
+        (mock) -> MOI.Utilities.mock_optimize!(
+            mock, config.optimal_status, [2.8, 0.0, 0.8, 1.2]
+        )
+    )
+    MOI.Test.mixed_complementaritytest(mock, config)
+end
+
+@testset "math_program_complementarity_constraints" begin
+    mock = MOI.Utilities.MockOptimizer(
+        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
+    )
+    config = MOI.Test.TestConfig(optimal_status = MOI.LOCALLY_SOLVED)
+    MOI.Utilities.set_mock_optimize!(
+        mock,
+        (mock) -> MOI.Utilities.mock_optimize!(
+            mock, config.optimal_status, [1.0, 0.0, 3.5, 0.0, 0.0, 0.0, 3.0, 6.0]
+        )
+    )
+    MOI.Test.math_program_complementarity_constraintstest(mock, config)
+end