diff --git a/test/instances.jl b/test/instances.jl index 5fac208..5914456 100644 --- a/test/instances.jl +++ b/test/instances.jl @@ -3,7 +3,12 @@ # Use of this source code is governed by an MIT-style license that can be found # in the LICENSE.md file or at https://opensource.org/licenses/MIT. +module Instances + using JuMP +using Test + +import Ipopt # Solution reported by Sven Leyffer in MacMPEC: # https://wiki.mcs.anl.gov/leyffer/index.php/MacMPEC @@ -19,6 +24,29 @@ const MACMPEC_SOLUTIONS = Dict{Symbol,Float64}( :water_net_model => 930.6687, # Solution reported is slightly worse than in MacMPEC (929.169) ) +# Some instances are flaky in CI around the reference optimum, so we loosen +# the tolerance on a per-instance basis. +const MACMPEC_TOLERANCES = Dict(:water_net_model => (rtol = 1e-2, atol = 1e-2)) + +function runtests(make_opt) + is_test(s) = !startswith("$s", "#") && endswith("$s", "_model") + @testset "$name" for name in filter(is_test, names(@__MODULE__; all = true)) + model = getfield(@__MODULE__, name)() + set_optimizer(model, () -> make_opt(Ipopt.Optimizer())) + set_attribute(model, "bound_relax_factor", 0.0) + set_attribute(model, "mu_strategy", "adaptive") + set_attribute(model, "bound_push", 1e-1) + set_silent(model) + optimize!(model) + @test is_solved_and_feasible(model) + if haskey(MACMPEC_SOLUTIONS, name) + tol = get(MACMPEC_TOLERANCES, name, (rtol = 1e-4, atol = 1e-4)) + @test objective_value(model) ≈ MACMPEC_SOLUTIONS[name] rtol=tol.rtol atol=tol.atol + end + end + return +end + # Ex. (2.2) from "Local convergence of SQP methods for MPECs". # Solution is (0.5, 0.5), basic multipliers is (0, 1, 0) function fletcher_leyffer_ex1_model() @@ -439,3 +467,5 @@ function water_net_model() return model end + +end # module Instances diff --git a/test/runtests.jl b/test/runtests.jl index 73393c7..72e0982 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,858 +3,12 @@ # Use of this source code is governed by an MIT-style license that can be found # in the LICENSE.md file or at https://opensource.org/licenses/MIT. -module Instances - -include("instances.jl") - -end - using Test -using JuMP -using Ipopt -using HiGHS - -using MathOptComplements - -function fletcher_leyffer_ex1_nonlinear_model() - model = Model() - @variable(model, z[1:2]) - @variable(model, slack >= 0.0) - set_lower_bound(z[2], 0) - @objective(model, Min, (z[1] - 1)^2 + z[2]^2) - @constraint(model, z[2] - z[1] == slack) - @constraint(model, slack * z[2] <= 0) - return model -end - -function nonlinear_test_reformulated_model() - model = Model() - @variable(model, x >= 0.0) - @variable(model, y >= 0.0) - @variable(model, slack >= 0) - @objective(model, Min, x^2 + y^2 - 4*x*y) - # Build complementarity constraints with nonlinear expression - @constraint(model, sin(x) == slack) - @constraint(model, slack * y <= 0.0) - return model -end - -function simple_ncp() - model = Model() - @variable(model, y <= 0) - @constraint(model, y + 1 ⟂ y) - return model, [y], [-1.0] -end - -# Solve min_x x subject to 0 <= x <= 1 -function simple_lp_1() - model = Model() - @variable(model, 0 <= x <= 1) - @variable(model, μ) - @constraint(model, 1 - μ == 0.0) - @constraint(model, μ ⟂ x) - return model, [x, μ], [0.0, 1.0] -end - -# Solve min_x -x subject to 0 <= x <= 1 -function simple_lp_2() - model = Model() - @variable(model, 0 <= x <= 1) - @variable(model, μ) - @constraint(model, -1 - μ == 0.0) - @constraint(model, μ ⟂ x) - return model, [x, μ], [1.0, -1.0] -end - -# Solve min -x2 subject to (x1, x2) >= 0; x1 + x2 = 1 -function simple_lp_3() - model = Model() - @variable(model, 0.0 <= x[1:2]) - @variable(model, 0.0 <= z[1:2]) - @variable(model, y) - @constraint(model, -z[1] + y == 0.0) - @constraint(model, -1.0 - z[2] + y == 0.0) - @constraint(model, x[1] + x[2] == 1.0) - @constraint(model, z[1] ⟂ x[1]) - @constraint(model, z[2] ⟂ x[2]) - return model, [x; z; y], [0.0, 1.0, 1.0, 0.0, 1.0] -end - -function nonlinear_test_model() - model = Model() - @variable(model, x >= 0.0) - @variable(model, y >= 0.0) - @objective(model, Min, x^2 + y^2 - 4*x*y) - # Build complementarity constraints with nonlinear expression - @constraint(model, [sin(x), y] ∈ MOI.Complements(2)) - return model -end - -# Variable in the left-hand-side should not have two bounds -function test_vertical_mispecified_1() - model = Model() - @variable(model, 0.0 <= x <= 1.0) - @variable(model, 0.0 <= y) - @constraint(model, [x, y] ∈ MOI.Complements(2)) - return model -end - -# Variable in the right-hand-side should not be an expression -function test_vertical_mispecified_2() - model = Model() - @variable(model, 0.0 <= x) - @variable(model, 0.0 <= y) - @constraint(model, [x, 1.0*y + x] ∈ MOI.Complements(2)) - return model -end - -function test_vertical_formulation() - model = Model() - # Case 1: LHS is already a variable (do nothing) - @variable(model, x1) - @variable(model, 0.0 <= y1) - @constraint(model, [x1, y1] ∈ MOI.Complements(2)) - # Case 2: RHS is unbounded (convert LHS to equality) - @variable(model, x2) - @variable(model, y2) - @constraint(model, [1.0*x2, y2] ∈ MOI.Complements(2)) - # Case 3: LHS is a ScalarAffineFunction with a single variable - @variable(model, x3) - @variable(model, 0.0 <= y3) - @constraint(model, [1.0*x3, y3] ∈ MOI.Complements(2)) - return model -end - -function test_nonlinear_reformulation() - model = Model() - # Case 1: Complementarity defined as lower-bound on RHS - @variable(model, x1) - @variable(model, 0.0 <= y1) - @constraint(model, [x1, y1] ∈ MOI.Complements(2)) - # Case 2: Complementarity defined as upper-bound on RHS - @variable(model, x2) - @variable(model, y2 <= 1.0) - @constraint(model, [x2, y2] ∈ MOI.Complements(2)) - return model -end - -# LHS has a non-trivial lower-bound -function test_nonlinear_mispecified_1() - model = Model() - @variable(model, 1.0 <= x1) - @variable(model, 0.0 <= y1) - @constraint(model, [x1, y1] ∈ MOI.Complements(2)) - return model -end - -# LHS has a non-trivial upper-bound -function test_nonlinear_mispecified_2() - model = Model() - @variable(model, x1 <= 1.0) - @variable(model, y1 <= 0.0) - @constraint(model, [x1, y1] ∈ MOI.Complements(2)) - return model -end - -# RHS is unbounded -function test_nonlinear_mispecified_3() - model = Model() - @variable(model, x1) - @variable(model, y1) - @constraint(model, [x1, y1] ∈ MOI.Complements(2)) - return model -end - -expected_models = Dict( - Instances.fletcher_leyffer_ex1_model => - fletcher_leyffer_ex1_nonlinear_model, -) - -# The tests that do not set `DefaultComplementarityReformulation` should work -# both with `MathOptComplements.Optimizer` and with a `MOI.Bridges.LazyBridgeOptimizer` -# in which all `MathOptComplements` bridges have been registered via `add_all_bridges`. -const OPTIMIZER_FACTORIES = [ - ( - "MathOptComplements.Optimizer", - inner -> MathOptComplements.Optimizer(inner), - ), - ( - "LazyBridgeOptimizer + add_all_bridges", - inner -> begin - lazy = MOI.Bridges.full_bridge_optimizer(inner, Float64) - MathOptComplements.Bridges.add_all_bridges(lazy) - return lazy - end, - ), -] - -# Some instances are flaky in CI around the reference optimum, so we loosen -# the tolerance on a per-instance basis. -const MACMPEC_TOLERANCES = Dict(:water_net_model => (rtol = 1e-2, atol = 1e-2)) - -function test_model(model_func, make_opt) - model = model_func() - set_optimizer(model, () -> make_opt(Ipopt.Optimizer())) - JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0) - JuMP.set_optimizer_attribute(model, "mu_strategy", "adaptive") - JuMP.set_optimizer_attribute(model, "bound_push", 1e-1) - JuMP.set_silent(model) - JuMP.optimize!(model) - name = Symbol(model_func) - @test JuMP.is_solved_and_feasible(model) - if haskey(Instances.MACMPEC_SOLUTIONS, name) - tol = get(MACMPEC_TOLERANCES, name, (rtol = 1e-4, atol = 1e-4)) - @test JuMP.objective_value(model) ≈ Instances.MACMPEC_SOLUTIONS[name] rtol=tol.rtol atol=tol.atol - end -end - -function test_nonlinear_expr(original_model, reformulated_model) - model = original_model() - inner = MOI.Utilities.Model{Float64}() - set_optimizer(model, () -> MathOptComplements.Optimizer(inner)) - MOI.Utilities.attach_optimizer(model) - expected = reformulated_model() - return MOI.Bridges._test_structural_identical( - unsafe_backend(model).model, - backend(expected), - ) -end - -@testset "Test vertical formulation ($(opt_name))" for (opt_name, make_opt) in - OPTIMIZER_FACTORIES - - model = test_vertical_formulation() - set_optimizer(model, () -> make_opt(Ipopt.Optimizer())) - MOI.Utilities.attach_optimizer(model) - - model = test_vertical_mispecified_2() - set_optimizer(model, () -> make_opt(Ipopt.Optimizer())) - @test_throws Exception MOI.Utilities.attach_optimizer(model) -end - -instances = filter(names(Instances; all = true)) do name - # The function types start with `#` - s = String(name) - return endswith(s, "_model") && !startswith(s, "#") -end - -@testset "$(opt_name): $name" for (opt_name, make_opt) in OPTIMIZER_FACTORIES, - name in instances - - test_model(getfield(Instances, name), make_opt) -end - -@testset "Test reformulation for $original_model" for ( - original_model, - reformulated_model, -) in [ - (nonlinear_test_model, nonlinear_test_reformulated_model), - ( - Instances.fletcher_leyffer_ex1_model, - fletcher_leyffer_ex1_nonlinear_model, - ), -] - test_nonlinear_expr(original_model, reformulated_model) -end - -@testset "Relaxation method: $(relax)" for relax in [ - MathOptComplements.ScholtesRelaxation(0.0), - MathOptComplements.FischerBurmeisterRelaxation(1e-8), - MathOptComplements.LiuFukushimaRelaxation(1e-8), - MathOptComplements.KanzowSchwarzRelaxation(1e-8), -] - @testset "Test reformulation" begin - model = test_nonlinear_reformulation() - set_optimizer( - model, - () -> MathOptComplements.Optimizer(Ipopt.Optimizer()), - ) - MOI.set( - model, - MathOptComplements.DefaultComplementarityReformulation(), - relax, - ) - MOI.Utilities.attach_optimizer(model) - - for test_func in - (test_nonlinear_mispecified_1, test_nonlinear_mispecified_2) - model = test_func() - set_optimizer( - model, - () -> MathOptComplements.Optimizer(Ipopt.Optimizer()), - ) - MOI.set( - model, - MathOptComplements.DefaultComplementarityReformulation(), - relax, - ) - @test_throws Exception MOI.Utilities.attach_optimizer(model) - end - end - - @testset "Solve Fletcher-Leyffer Ex1 problem" begin - model = Instances.fletcher_leyffer_ex1_model() - JuMP.set_optimizer( - model, - () -> MathOptComplements.Optimizer(Ipopt.Optimizer()), - ) - MOI.set( - model, - MathOptComplements.DefaultComplementarityReformulation(), - relax, - ) - JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0) - JuMP.set_silent(model) - JuMP.optimize!(model) - - @test JuMP.is_solved_and_feasible(model) - @test JuMP.objective_value(model) ≈ 0.5 atol=1e-7 - @test JuMP.value.(model[:z]) ≈ [0.5, 0.5] atol=1e-7 - end - - @testset "Solve NCP problem $(func)" for func in [simple_ncp, simple_lp_3] - model, vars, sol = func() - JuMP.set_optimizer( - model, - () -> MathOptComplements.Optimizer(Ipopt.Optimizer()), - ) - MOI.set( - model, - MathOptComplements.DefaultComplementarityReformulation(), - relax, - ) - JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0) - JuMP.set_silent(model) - JuMP.optimize!(model) - - @test JuMP.is_solved_and_feasible(model) - @test JuMP.value.(vars) ≈ sol atol=1e-7 - end -end - -# N.B.: at the moment, mixed-complementarity problems are supported only -# with ScholtesRelaxation and with FischerBurmeisterRelaxation -@testset "Mixed-complementarity problem with $(relax)" for relax in [ - MathOptComplements.ScholtesRelaxation(0.0), - MathOptComplements.FischerBurmeisterRelaxation(1e-8), -] - @testset "Solve NCP problem $(func)" for func in [simple_lp_1] #, simple_lp_2] - model, vars, sol = func() - JuMP.set_optimizer( - model, - () -> MathOptComplements.Optimizer( - MOI.instantiate(Ipopt.Optimizer, with_cache_type = Float64), - ), - ) - MOI.set( - model, - MathOptComplements.DefaultComplementarityReformulation(), - relax, - ) - JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0) - JuMP.set_silent(model) - JuMP.optimize!(model) - - inner = backend(model).optimizer.model.model - - if relax isa MathOptComplements.ScholtesRelaxation - F = MOI.ScalarQuadraticFunction{Float64} - G = MOI.ScalarNonlinearFunction - else - G = MOI.ScalarQuadraticFunction{Float64} - F = MOI.ScalarNonlinearFunction - end - @test MOI.get( - inner, - MOI.NumberOfConstraints{F,MOI.LessThan{Float64}}(), - ) > 0 - @test MOI.get( - inner, - MOI.NumberOfConstraints{G,MOI.LessThan{Float64}}(), - ) == 0 - - @test JuMP.is_solved_and_feasible(model) - @test JuMP.value.(vars) ≈ sol atol=1e-7 - end -end - -@testset "Per-constraint reformulation" begin - model = Model() - @variable(model, x1) - @variable(model, 0.0 <= y1) - c1 = @constraint(model, x1 ⟂ y1) - @variable(model, x2) - @variable(model, 0.0 <= y2) - c2 = @constraint(model, x2 ⟂ y2) - @objective(model, Min, (x1 - 1)^2 + y1^2 + (x2 - 1)^2 + y2^2) - JuMP.set_optimizer( - model, - () -> MathOptComplements.Optimizer( - MOI.instantiate(Ipopt.Optimizer, with_cache_type = Float64), - ), - with_cache_type = Float64, - ) - # Default is Scholtes - MOI.set( - model, - MathOptComplements.DefaultComplementarityReformulation(), - MathOptComplements.ScholtesRelaxation(0.0), - ) - # Override c1 with FischerBurmeister - MOI.set( - model, - MathOptComplements.ComplementarityReformulation(), - c1, - MathOptComplements.FischerBurmeisterRelaxation(1e-8), - ) - @test MOI.supports( - JuMP.unsafe_backend(model), - MathOptComplements.DefaultComplementarityReformulation(), - ) - b = JuMP.unsafe_backend(model) - attr = MathOptComplements.ComplementarityReformulation() - F = MOI.VectorOfVariables - S = MOI.Complements - @test MOI.Bridges.is_bridged(b, S) - @test MOI.supports_add_constrained_variables(b, S) - @test !MOI.Bridges.is_variable_bridged(b, S) - bridge_type = MOI.Bridges.Constraint.concrete_bridge_type(b, F, S) - @test bridge_type == - MathOptComplements.Bridges.SpecifySetTypeBridge{Float64} - @test MOI.supports(b, attr, bridge_type) - @test MOI.supports( - JuMP.unsafe_backend(model), - MathOptComplements.ComplementarityReformulation(), - MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Complements}, - ) - @test MOI.supports( - JuMP.unsafe_backend(model), - MathOptComplements.ComplementarityReformulation(), - MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Complements}, - ) - @test MOI.supports( - JuMP.unsafe_backend(model), - MathOptComplements.ComplementarityReformulation(), - MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Complements}, - ) - @test MOI.get( - model, - MathOptComplements.ComplementarityReformulation(), - c1, - ) isa MathOptComplements.FischerBurmeisterRelaxation - JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0) - JuMP.set_silent(model) - JuMP.optimize!(model) - @test JuMP.is_solved_and_feasible(model) - # Test get through the LazyBridgeOptimizer - lazy = JuMP.backend(model).optimizer - @test !MOI.Bridges.is_bridged(lazy, S) - ci_mapped = first( - MOI.get( - lazy, - MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Complements}(), - ), - ) - @test MOI.get( - lazy, - MathOptComplements.ComplementarityReformulation(), - ci_mapped, - ) isa MathOptComplements.FischerBurmeisterRelaxation - @test MOI.get( - model, - MathOptComplements.ComplementarityReformulation(), - c1, - ) isa MathOptComplements.FischerBurmeisterRelaxation - @test isnothing( - MOI.get(model, MathOptComplements.ComplementarityReformulation(), c2), - ) -end - -@testset "Per-constraint reformulation after optimize!" begin - model = Model() - @variable(model, x) - @variable(model, 0.0 <= y) - c = @constraint(model, x ⟂ y) - @objective(model, Min, (x - 1)^2 + y^2) - JuMP.set_optimizer( - model, - () -> MathOptComplements.Optimizer( - MOI.instantiate(Ipopt.Optimizer, with_cache_type = Float64), - ), - with_cache_type = Float64, - ) - MOI.set( - model, - MathOptComplements.DefaultComplementarityReformulation(), - MathOptComplements.ScholtesRelaxation(0.0), - ) - JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0) - JuMP.set_silent(model) - JuMP.optimize!(model) - @test JuMP.is_solved_and_feasible(model) - # Scholtes produces a quadratic constraint - inner = backend(model).optimizer.model.optimizer.model - @test MOI.get( - inner, - MOI.NumberOfConstraints{ - MOI.ScalarQuadraticFunction{Float64}, - MOI.LessThan{Float64}, - }(), - ) == 1 - # Change reformulation after first optimize! (bridge.constraints is populated) - MOI.set( - model, - MathOptComplements.ComplementarityReformulation(), - c, - MathOptComplements.FischerBurmeisterRelaxation(1e-8), - ) - @test MOI.get( - model, - MathOptComplements.ComplementarityReformulation(), - c, - ) isa MathOptComplements.FischerBurmeisterRelaxation - JuMP.optimize!(model) - @test JuMP.is_solved_and_feasible(model) - inner = backend(model).optimizer.model.optimizer.model - # FB produces a nonlinear constraint - @test MOI.get( - inner, - MOI.NumberOfConstraints{ - MOI.ScalarNonlinearFunction, - MOI.LessThan{Float64}, - }(), - ) == 1 -end - -@testset "ComplementsWithSetType dimension" begin - s = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives}(4) - @test MOI.dimension(s) == 4 -end - -@testset "Optimizer bridge dispatch" begin - # NLP path: VectorAffineFunction in ComplementsWithSetType → NonlinearBridge - opt_nlp = MathOptComplements.Optimizer( - MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), - ) - S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} - @test MOI.Bridges.is_bridged(opt_nlp, S) - @test MOI.Bridges.supports_bridging_constrained_variable(opt_nlp, S) - @test MOI.Bridges.supports_bridging_constraint( - opt_nlp, - MOI.VectorAffineFunction{Float64}, - S, - ) - @test MOI.Bridges.bridge_type( - opt_nlp, - MOI.VectorAffineFunction{Float64}, - S, - ) == MathOptComplements.Bridges.NonlinearBridge{Float64,MOI.Nonnegatives} - # SOS1 path: VectorOfVariables in ComplementsWithSetType{Zeros} → ToSOS1Bridge - opt_sos1 = MathOptComplements.Optimizer( - MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64), - ) - S_zeros = MathOptComplements.ComplementsWithSetType{MOI.Zeros} - @test MOI.Bridges.bridge_type(opt_sos1, MOI.VectorOfVariables, S_zeros) == - MathOptComplements.Bridges.ToSOS1Bridge{Float64} -end - -@testset "ComplementarityReformulation supports on bridge types" begin - model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) - attr = MathOptComplements.ComplementarityReformulation() - @test MOI.supports( - model, - attr, - MathOptComplements.Bridges.SplitIntervalBridge, - ) - @test MOI.supports(model, attr, MathOptComplements.Bridges.FlipSignBridge) - @test MOI.supports( - model, - attr, - MathOptComplements.Bridges.ComplementsVectorizeBridge, - ) - @test MOI.supports(model, attr, MathOptComplements.Bridges.NonlinearBridge) -end - -@testset "ComplementarityReformulation through SplitInterval" begin - # Stack SplitInterval on top of a model with NonlinearBridge - inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) - nl_model = MOI.Bridges.Constraint.SingleBridgeOptimizer{ - MathOptComplements.Bridges.NonlinearBridge{Float64}, - }( - inner, - ) - model = MOI.Bridges.Constraint.SingleBridgeOptimizer{ - MathOptComplements.Bridges.SplitIntervalBridge{Float64}, - }( - nl_model, - ) - x = MOI.add_variable(model) - y, _ = MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0)) - ci = MOI.add_constraint( - model, - MOI.VectorOfVariables([x, y]), - MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}}(2), - ) - # Default Scholtes → quadratic constraints - MOI.Bridges.final_touch(model) - MOI.Bridges.final_touch(nl_model) - @test MOI.get( - inner, - MOI.NumberOfConstraints{ - MOI.ScalarQuadraticFunction{Float64}, - MOI.LessThan{Float64}, - }(), - ) > 0 - # Change to FB and re-reformulate → nonlinear constraints - attr = MathOptComplements.ComplementarityReformulation() - MOI.set( - model, - attr, - ci, - MathOptComplements.FischerBurmeisterRelaxation(1e-8), - ) - MOI.Bridges.final_touch(model) - MOI.Bridges.final_touch(nl_model) - @test MOI.get( - inner, - MOI.NumberOfConstraints{ - MOI.ScalarNonlinearFunction, - MOI.LessThan{Float64}, - }(), - ) > 0 - # Test SplitIntervalBridge metadata - b_split = MOI.Bridges.bridge(model, ci) - @test b_split isa MathOptComplements.Bridges.SplitIntervalBridge - @test length( - MOI.get( - b_split, - MOI.ListOfConstraintIndices{ - MOI.VariableIndex, - MOI.GreaterThan{Float64}, - }(), - ), - ) == 1 - @test length( - MOI.get( - b_split, - MOI.ListOfConstraintIndices{ - MOI.VariableIndex, - MOI.LessThan{Float64}, - }(), - ), - ) == 1 -end - -@testset "ComplementarityReformulation through FlipSign" begin - inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) - model = MOI.Bridges.Constraint.SingleBridgeOptimizer{ - MathOptComplements.Bridges.FlipSignBridge{Float64}, - }( - inner, - ) - x, _ = MOI.add_constrained_variable(model, MOI.LessThan(0.0)) - y, _ = MOI.add_constrained_variable(model, MOI.LessThan(0.0)) - ci = MOI.add_constraint( - model, - MOI.VectorOfVariables([x, y]), - MathOptComplements.ComplementsWithSetType{MOI.Nonpositives}(2), - ) - attr = MathOptComplements.ComplementarityReformulation() - relax = MathOptComplements.FischerBurmeisterRelaxation(1e-8) - MOI.set(model, attr, ci, relax) - # FlipSign flips Nonpositives → Nonnegatives, check child exists - S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} - @test MOI.get( - inner, - MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},S}(), - ) == 1 - # Check the attribute was propagated to the child constraint - b = MOI.Bridges.bridge(model, ci) - @test MOI.get(inner, attr, b.constraint) isa - MathOptComplements.FischerBurmeisterRelaxation -end - -@testset "ComplementarityReformulation through ComplementsVectorize" begin - inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) - model = MOI.Bridges.Constraint.SingleBridgeOptimizer{ - MathOptComplements.Bridges.ComplementsVectorizeBridge{Float64}, - }( - inner, - ) - x, _ = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0)) - y, _ = MOI.add_constrained_variable(model, MOI.GreaterThan(3.0)) - ci = MOI.add_constraint( - model, - MOI.VectorOfVariables([x, y]), - MathOptComplements.ComplementsWithSetType{MOI.GreaterThan{Float64}}(2), - ) - attr = MathOptComplements.ComplementarityReformulation() - relax = MathOptComplements.FischerBurmeisterRelaxation(1e-8) - MOI.set(model, attr, ci, relax) - # ComplementsVectorize shifts GreaterThan → Nonnegatives, check child exists - S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} - @test MOI.get( - inner, - MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},S}(), - ) == 1 - # Check the attribute was propagated to the child constraint - b = MOI.Bridges.bridge(model, ci) - @test MOI.get(inner, attr, b.constraint) isa - MathOptComplements.FischerBurmeisterRelaxation -end - -@testset "SpecifySetTypeBridge ($(opt_name))" for (opt_name, make_opt) in - OPTIMIZER_FACTORIES - - @testset "Lower bound (Nonnegatives)" begin - model = Model() - @variable(model, x) - @variable(model, 0.0 <= y) - @constraint(model, [x, y] ∈ MOI.Complements(2)) - inner = MOI.Utilities.Model{Float64}() - set_optimizer(model, () -> make_opt(inner)) - MOI.Utilities.attach_optimizer(model) - # ComplementsWithSetType is bridged further to nonlinear constraints - S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} - @test MOI.get( - inner, - MOI.NumberOfConstraints{MOI.VectorOfVariables,S}(), - ) == 0 - end - - @testset "Lower bound (GreaterThan)" begin - model = Model() - @variable(model, x) - @variable(model, 3.0 <= y) - @constraint(model, [x, y] ∈ MOI.Complements(2)) - inner = MOI.Utilities.Model{Float64}() - set_optimizer(model, () -> make_opt(inner)) - MOI.Utilities.attach_optimizer(model) - S = MathOptComplements.ComplementsWithSetType{MOI.GreaterThan{Float64}} - @test MOI.get( - inner, - MOI.NumberOfConstraints{MOI.VectorOfVariables,S}(), - ) == 0 - end - - @testset "Upper bound (LessThan)" begin - model = Model() - @variable(model, x) - @variable(model, y <= 1.0) - @constraint(model, [x, y] ∈ MOI.Complements(2)) - inner = MOI.Utilities.Model{Float64}() - set_optimizer(model, () -> make_opt(inner)) - MOI.Utilities.attach_optimizer(model) - S = MathOptComplements.ComplementsWithSetType{MOI.LessThan{Float64}} - @test MOI.get( - inner, - MOI.NumberOfConstraints{MOI.VectorOfVariables,S}(), - ) == 0 - end - - @testset "Range case (x1 unbounded → Interval)" begin - model = Model() - @variable(model, x) - @variable(model, 0.0 <= y <= 1.0) - @constraint(model, [x, y] ∈ MOI.Complements(2)) - inner = MOI.Utilities.Model{Float64}() - set_optimizer(model, () -> make_opt(inner)) - MOI.Utilities.attach_optimizer(model) - S = MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}} - @test MOI.get( - inner, - MOI.NumberOfConstraints{MOI.VectorOfVariables,S}(), - ) == 0 - end - - @testset "Range case (x1 bounded nonneg)" begin - model = Model() - @variable(model, 0.0 <= x <= 10.0) - @variable(model, 0.0 <= y <= 10.0) - @constraint(model, [x, y] ∈ MOI.Complements(2)) - inner = MOI.Utilities.Model{Float64}() - set_optimizer(model, () -> make_opt(inner)) - MOI.Utilities.attach_optimizer(model) - S = MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}} - @test MOI.get( - inner, - MOI.NumberOfConstraints{MOI.VectorOfVariables,S}(), - ) == 0 - end -end - -@testset "Per-constraint reformulation with VerticalBridge" begin - # Use an expression LHS so that the constraint goes through VerticalBridge - model = Model() - @variable(model, x >= 0.0) - @variable(model, y >= 0.0) - c = @constraint(model, [x + y, y] ∈ MOI.Complements(2)) - @objective(model, Min, x^2 + y^2) - JuMP.set_optimizer( - model, - () -> MathOptComplements.Optimizer( - MOI.instantiate(Ipopt.Optimizer, with_cache_type = Float64), - ), - ) - attr = MathOptComplements.ComplementarityReformulation() - reformulation = MathOptComplements.FischerBurmeisterRelaxation(1e-8) - MOI.set(model, attr, c, reformulation) - JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0) - JuMP.set_silent(model) - JuMP.optimize!(model) - @test MOI.supports(backend(model), attr, typeof(index(c))) - @test MOI.get(model, attr, c) == reformulation - @test JuMP.is_solved_and_feasible(model) -end - -@testset "Bridge chain: SpecifySetType → SplitInterval → FlipSign → ToSOS1" begin - # HiGHS does not support ScalarNonlinearFunction, so the Optimizer - # uses the SOS1 path instead of NonlinearBridge. - opt = MathOptComplements.Optimizer( - MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64), - ) - - # Create a model with an Interval complementarity constraint - x = MOI.add_variable(opt) - y, _ = MOI.add_constrained_variable(opt, MOI.Interval(0.0, 1.0)) - ci = MOI.add_constraint( - opt, - MOI.VectorOfVariables([x, y]), - MOI.Complements(2), - ) - MOI.Bridges.final_touch(opt) - - # Step 1: Complements → SpecifySetTypeBridge - bridge1 = MOI.Bridges.bridge(opt, ci) - @test bridge1 isa MathOptComplements.Bridges.SpecifySetTypeBridge - ci_interval = bridge1.constraints[1] - @test ci_interval isa MOI.ConstraintIndex{ - MOI.VectorOfVariables, - MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}}, - } - - # Step 2: Interval → SplitIntervalBridge (not NonlinearBridge!) - bridge2 = MOI.Bridges.bridge(opt, ci_interval) - @test bridge2 isa MathOptComplements.Bridges.SplitIntervalBridge - - # Step 3: LessThan part → FlipSignBridge - bridge3 = MOI.Bridges.bridge(opt, bridge2.upper) - @test bridge3 isa MathOptComplements.Bridges.FlipSignBridge -end - -@testset "add_all_bridges(::JuMP.GenericModel)" begin - model = Model(Ipopt.Optimizer) - MathOptComplements.Bridges.add_all_bridges(model) - set_silent(model) - @variable(model, y <= 0) - @constraint(model, y + 1 ⟂ y) - optimize!(model) - @test is_solved_and_feasible(model) - @test value(y) ≈ -1.0 atol = 1e-7 -end is_test(f) = startswith(f, "test_") && endswith(f, ".jl") -@testset "$file" for file in filter(is_test, readdir("Bridges")) - include(joinpath(@__DIR__, "Bridges", file)) +@testset "$root" for (root, dirs, files) in walkdir(@__DIR__) + @testset "$file" for file in filter(is_test, files) + include(joinpath(root, file)) + end end diff --git a/test/test_optimizer.jl b/test/test_optimizer.jl new file mode 100644 index 0000000..b5bf20c --- /dev/null +++ b/test/test_optimizer.jl @@ -0,0 +1,799 @@ +# Copyright (c) 2025 François Pacaud, Benoît Legat, and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +module TestOptimizer + +using JuMP +using Test + +import HiGHS +import Ipopt +import MathOptComplements +import MathOptInterface as MOI + +include(joinpath(@__DIR__, "instances.jl")) + +function runtests() + is_test(name) = startswith("$name", "test_") + @testset "$name" for name in filter(is_test, names(@__MODULE__; all = true)) + getfield(@__MODULE__, name)() + end + return +end + +function test_per_constraint_reformulation() + model = Model() + @variable(model, x1) + @variable(model, 0.0 <= y1) + c1 = @constraint(model, x1 ⟂ y1) + @variable(model, x2) + @variable(model, 0.0 <= y2) + c2 = @constraint(model, x2 ⟂ y2) + @objective(model, Min, (x1 - 1)^2 + y1^2 + (x2 - 1)^2 + y2^2) + set_optimizer( + model, + () -> MathOptComplements.Optimizer( + MOI.instantiate(Ipopt.Optimizer; with_cache_type = Float64), + ); + with_cache_type = Float64, + ) + # Default is Scholtes + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + MathOptComplements.ScholtesRelaxation(0.0), + ) + # Override c1 with FischerBurmeister + MOI.set( + model, + MathOptComplements.ComplementarityReformulation(), + c1, + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + ) + @test MOI.supports( + unsafe_backend(model), + MathOptComplements.DefaultComplementarityReformulation(), + ) + b = unsafe_backend(model) + attr = MathOptComplements.ComplementarityReformulation() + F = MOI.VectorOfVariables + S = MOI.Complements + @test MOI.Bridges.is_bridged(b, S) + @test MOI.supports_add_constrained_variables(b, S) + @test !MOI.Bridges.is_variable_bridged(b, S) + bridge_type = MOI.Bridges.Constraint.concrete_bridge_type(b, F, S) + @test bridge_type == + MathOptComplements.Bridges.SpecifySetTypeBridge{Float64} + @test MOI.supports(b, attr, bridge_type) + @test MOI.supports( + unsafe_backend(model), + MathOptComplements.ComplementarityReformulation(), + MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Complements}, + ) + @test MOI.supports( + unsafe_backend(model), + MathOptComplements.ComplementarityReformulation(), + MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Complements}, + ) + @test MOI.supports( + unsafe_backend(model), + MathOptComplements.ComplementarityReformulation(), + MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Complements}, + ) + @test MOI.get( + model, + MathOptComplements.ComplementarityReformulation(), + c1, + ) isa MathOptComplements.FischerBurmeisterRelaxation + set_attribute(model, "bound_relax_factor", 0.0) + set_silent(model) + optimize!(model) + @test is_solved_and_feasible(model) + # Test get through the LazyBridgeOptimizer + lazy = backend(model).optimizer + @test !MOI.Bridges.is_bridged(lazy, S) + ci_mapped = first( + MOI.get( + lazy, + MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Complements}(), + ), + ) + @test MOI.get( + lazy, + MathOptComplements.ComplementarityReformulation(), + ci_mapped, + ) isa MathOptComplements.FischerBurmeisterRelaxation + @test MOI.get( + model, + MathOptComplements.ComplementarityReformulation(), + c1, + ) isa MathOptComplements.FischerBurmeisterRelaxation + @test isnothing( + MOI.get(model, MathOptComplements.ComplementarityReformulation(), c2), + ) + return +end + +function test_per_constraint_reformulation_after_optimize!() + model = Model() + @variable(model, x) + @variable(model, 0.0 <= y) + c = @constraint(model, x ⟂ y) + @objective(model, Min, (x - 1)^2 + y^2) + set_optimizer( + model, + () -> MathOptComplements.Optimizer( + MOI.instantiate(Ipopt.Optimizer; with_cache_type = Float64), + ); + with_cache_type = Float64, + ) + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + MathOptComplements.ScholtesRelaxation(0.0), + ) + set_attribute(model, "bound_relax_factor", 0.0) + set_silent(model) + optimize!(model) + @test is_solved_and_feasible(model) + # Scholtes produces a quadratic constraint + inner = backend(model).optimizer.model.optimizer.model + @test MOI.get( + inner, + MOI.NumberOfConstraints{ + MOI.ScalarQuadraticFunction{Float64}, + MOI.LessThan{Float64}, + }(), + ) == 1 + # Change reformulation after first optimize! (bridge.constraints is populated) + MOI.set( + model, + MathOptComplements.ComplementarityReformulation(), + c, + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + ) + @test MOI.get( + model, + MathOptComplements.ComplementarityReformulation(), + c, + ) isa MathOptComplements.FischerBurmeisterRelaxation + optimize!(model) + @test is_solved_and_feasible(model) + inner = backend(model).optimizer.model.optimizer.model + # FB produces a nonlinear constraint + @test MOI.get( + inner, + MOI.NumberOfConstraints{ + MOI.ScalarNonlinearFunction, + MOI.LessThan{Float64}, + }(), + ) == 1 + return +end + +function test_ComplementsWithSetType_dimension() + s = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives}(4) + @test MOI.dimension(s) == 4 + return +end + +function test_Optimizer_bridge_dispatch() + # NLP path: VectorAffineFunction in ComplementsWithSetType → NonlinearBridge + opt_nlp = MathOptComplements.Optimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + ) + S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} + @test MOI.Bridges.is_bridged(opt_nlp, S) + @test MOI.Bridges.supports_bridging_constrained_variable(opt_nlp, S) + @test MOI.Bridges.supports_bridging_constraint( + opt_nlp, + MOI.VectorAffineFunction{Float64}, + S, + ) + @test MOI.Bridges.bridge_type( + opt_nlp, + MOI.VectorAffineFunction{Float64}, + S, + ) == MathOptComplements.Bridges.NonlinearBridge{Float64,MOI.Nonnegatives} + # SOS1 path: VectorOfVariables in ComplementsWithSetType{Zeros} → ToSOS1Bridge + opt_sos1 = MathOptComplements.Optimizer( + MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64), + ) + S_zeros = MathOptComplements.ComplementsWithSetType{MOI.Zeros} + @test MOI.Bridges.bridge_type(opt_sos1, MOI.VectorOfVariables, S_zeros) == + MathOptComplements.Bridges.ToSOS1Bridge{Float64} + return +end + +function test_ComplementarityReformulation_supports_on_bridge_types() + model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + attr = MathOptComplements.ComplementarityReformulation() + @test MOI.supports( + model, + attr, + MathOptComplements.Bridges.SplitIntervalBridge, + ) + @test MOI.supports(model, attr, MathOptComplements.Bridges.FlipSignBridge) + @test MOI.supports( + model, + attr, + MathOptComplements.Bridges.ComplementsVectorizeBridge, + ) + @test MOI.supports(model, attr, MathOptComplements.Bridges.NonlinearBridge) + return +end + +function test_ComplementarityReformulation_through_SplitInterval() + # Stack SplitInterval on top of a model with NonlinearBridge + inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + B = MathOptComplements.Bridges.NonlinearBridge{Float64} + nl_model = MOI.Bridges.Constraint.SingleBridgeOptimizer{B}(inner) + B = MathOptComplements.Bridges.SplitIntervalBridge{Float64} + model = MOI.Bridges.Constraint.SingleBridgeOptimizer{B}(nl_model) + x = MOI.add_variable(model) + y, _ = MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0)) + ci = MOI.add_constraint( + model, + MOI.VectorOfVariables([x, y]), + MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}}(2), + ) + # Default Scholtes → quadratic constraints + MOI.Bridges.final_touch(model) + MOI.Bridges.final_touch(nl_model) + @test MOI.get( + inner, + MOI.NumberOfConstraints{ + MOI.ScalarQuadraticFunction{Float64}, + MOI.LessThan{Float64}, + }(), + ) > 0 + # Change to FB and re-reformulate → nonlinear constraints + attr = MathOptComplements.ComplementarityReformulation() + MOI.set( + model, + attr, + ci, + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + ) + MOI.Bridges.final_touch(model) + MOI.Bridges.final_touch(nl_model) + @test MOI.get( + inner, + MOI.NumberOfConstraints{ + MOI.ScalarNonlinearFunction, + MOI.LessThan{Float64}, + }(), + ) > 0 + # Test SplitIntervalBridge metadata + b_split = MOI.Bridges.bridge(model, ci) + @test b_split isa MathOptComplements.Bridges.SplitIntervalBridge + @test length( + MOI.get( + b_split, + MOI.ListOfConstraintIndices{ + MOI.VariableIndex, + MOI.GreaterThan{Float64}, + }(), + ), + ) == 1 + @test length( + MOI.get( + b_split, + MOI.ListOfConstraintIndices{ + MOI.VariableIndex, + MOI.LessThan{Float64}, + }(), + ), + ) == 1 + return +end + +function test_ComplementarityReformulation_through_FlipSign() + inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + B = MathOptComplements.Bridges.FlipSignBridge{Float64} + model = MOI.Bridges.Constraint.SingleBridgeOptimizer{B}(inner) + x, _ = MOI.add_constrained_variable(model, MOI.LessThan(0.0)) + y, _ = MOI.add_constrained_variable(model, MOI.LessThan(0.0)) + ci = MOI.add_constraint( + model, + MOI.VectorOfVariables([x, y]), + MathOptComplements.ComplementsWithSetType{MOI.Nonpositives}(2), + ) + attr = MathOptComplements.ComplementarityReformulation() + relax = MathOptComplements.FischerBurmeisterRelaxation(1e-8) + MOI.set(model, attr, ci, relax) + # FlipSign flips Nonpositives → Nonnegatives, check child exists + S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} + @test MOI.get( + inner, + MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},S}(), + ) == 1 + # Check the attribute was propagated to the child constraint + b = MOI.Bridges.bridge(model, ci) + @test MOI.get(inner, attr, b.constraint) isa + MathOptComplements.FischerBurmeisterRelaxation + return +end + +function test_ComplementarityReformulation_through_ComplementsVectorize() + inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + B = MathOptComplements.Bridges.ComplementsVectorizeBridge{Float64} + model = MOI.Bridges.Constraint.SingleBridgeOptimizer{B}(inner) + x, _ = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0)) + y, _ = MOI.add_constrained_variable(model, MOI.GreaterThan(3.0)) + ci = MOI.add_constraint( + model, + MOI.VectorOfVariables([x, y]), + MathOptComplements.ComplementsWithSetType{MOI.GreaterThan{Float64}}(2), + ) + attr = MathOptComplements.ComplementarityReformulation() + relax = MathOptComplements.FischerBurmeisterRelaxation(1e-8) + MOI.set(model, attr, ci, relax) + # ComplementsVectorize shifts GreaterThan → Nonnegatives, check child exists + S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} + @test MOI.get( + inner, + MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},S}(), + ) == 1 + # Check the attribute was propagated to the child constraint + b = MOI.Bridges.bridge(model, ci) + @test MOI.get(inner, attr, b.constraint) isa + MathOptComplements.FischerBurmeisterRelaxation + return +end + +function test_per_constraint_reformulation_with_VerticalBridge() + # Use an expression LHS so that the constraint goes through VerticalBridge + model = Model() + @variable(model, x >= 0.0) + @variable(model, y >= 0.0) + c = @constraint(model, [x + y, y] ∈ MOI.Complements(2)) + @objective(model, Min, x^2 + y^2) + set_optimizer( + model, + () -> MathOptComplements.Optimizer( + MOI.instantiate(Ipopt.Optimizer; with_cache_type = Float64), + ), + ) + attr = MathOptComplements.ComplementarityReformulation() + reformulation = MathOptComplements.FischerBurmeisterRelaxation(1e-8) + MOI.set(model, attr, c, reformulation) + set_attribute(model, "bound_relax_factor", 0.0) + set_silent(model) + optimize!(model) + @test MOI.supports(backend(model), attr, typeof(index(c))) + @test MOI.get(model, attr, c) == reformulation + @test is_solved_and_feasible(model) + return +end + +# SpecifySetType → SplitInterval → FlipSign → ToSOS1 +function test_bridge_chain() + # HiGHS does not support ScalarNonlinearFunction, so the Optimizer + # uses the SOS1 path instead of NonlinearBridge. + opt = MathOptComplements.Optimizer( + MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64), + ) + # Create a model with an Interval complementarity constraint + x = MOI.add_variable(opt) + y, _ = MOI.add_constrained_variable(opt, MOI.Interval(0.0, 1.0)) + ci = MOI.add_constraint( + opt, + MOI.VectorOfVariables([x, y]), + MOI.Complements(2), + ) + MOI.Bridges.final_touch(opt) + # Step 1: Complements → SpecifySetTypeBridge + bridge1 = MOI.Bridges.bridge(opt, ci) + @test bridge1 isa MathOptComplements.Bridges.SpecifySetTypeBridge + ci_interval = bridge1.constraints[1] + @test ci_interval isa MOI.ConstraintIndex{ + MOI.VectorOfVariables, + MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}}, + } + # Step 2: Interval → SplitIntervalBridge (not NonlinearBridge!) + bridge2 = MOI.Bridges.bridge(opt, ci_interval) + @test bridge2 isa MathOptComplements.Bridges.SplitIntervalBridge + # Step 3: LessThan part → FlipSignBridge + bridge3 = MOI.Bridges.bridge(opt, bridge2.upper) + @test bridge3 isa MathOptComplements.Bridges.FlipSignBridge + return +end + +function test_add_all_bridges_JuMP_GenericModel() + model = Model(Ipopt.Optimizer) + MathOptComplements.Bridges.add_all_bridges(model) + set_silent(model) + @variable(model, y <= 0) + @constraint(model, y + 1 ⟂ y) + optimize!(model) + @test is_solved_and_feasible(model) + @test value(y) ≈ -1.0 atol = 1e-7 +end + +function test_opt_default() + is_test(f) = startswith("$f", "_test_opt_") + for name in filter(is_test, names(@__MODULE__; all = true)) + getfield(@__MODULE__, name)(MathOptComplements.Optimizer) + end + return +end + +function test_opt_add_all_bridges() + is_test(f) = startswith("$f", "_test_opt_") + function make_opt(inner) + lazy = MOI.Bridges.full_bridge_optimizer(inner, Float64) + MathOptComplements.Bridges.add_all_bridges(lazy) + return lazy + end + for name in filter(is_test, names(@__MODULE__; all = true)) + getfield(@__MODULE__, name)(make_opt) + end + return +end + +function _test_opt_lower_bound_nonnegatives(make_opt) + model = Model() + @variable(model, x) + @variable(model, 0.0 <= y) + @constraint(model, [x, y] ∈ MOI.Complements(2)) + inner = MOI.Utilities.Model{Float64}() + set_optimizer(model, () -> make_opt(inner)) + MOI.Utilities.attach_optimizer(model) + # ComplementsWithSetType is bridged further to nonlinear constraints + S = MathOptComplements.ComplementsWithSetType{MOI.Nonnegatives} + @test MOI.get(inner, MOI.NumberOfConstraints{MOI.VectorOfVariables,S}()) == + 0 + return +end + +function _test_opt_lower_bound_greater_than(make_opt) + model = Model() + @variable(model, x) + @variable(model, 3.0 <= y) + @constraint(model, [x, y] ∈ MOI.Complements(2)) + inner = MOI.Utilities.Model{Float64}() + set_optimizer(model, () -> make_opt(inner)) + MOI.Utilities.attach_optimizer(model) + S = MathOptComplements.ComplementsWithSetType{MOI.GreaterThan{Float64}} + @test MOI.get(inner, MOI.NumberOfConstraints{MOI.VectorOfVariables,S}()) == + 0 + return +end + +function _test_opt_upper_bound_less_than(make_opt) + model = Model() + @variable(model, x) + @variable(model, y <= 1.0) + @constraint(model, [x, y] ∈ MOI.Complements(2)) + inner = MOI.Utilities.Model{Float64}() + set_optimizer(model, () -> make_opt(inner)) + MOI.Utilities.attach_optimizer(model) + S = MathOptComplements.ComplementsWithSetType{MOI.LessThan{Float64}} + @test MOI.get(inner, MOI.NumberOfConstraints{MOI.VectorOfVariables,S}()) == + 0 + return +end + +function _test_opt_range_case_x1_bounded_interval(make_opt) + model = Model() + @variable(model, x) + @variable(model, 0.0 <= y <= 1.0) + @constraint(model, [x, y] ∈ MOI.Complements(2)) + inner = MOI.Utilities.Model{Float64}() + set_optimizer(model, () -> make_opt(inner)) + MOI.Utilities.attach_optimizer(model) + S = MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}} + @test MOI.get(inner, MOI.NumberOfConstraints{MOI.VectorOfVariables,S}()) == + 0 + return +end + +function _test_opt_range_case_x1_bounded_nonneg(make_opt) + model = Model() + @variable(model, 0.0 <= x <= 10.0) + @variable(model, 0.0 <= y <= 10.0) + @constraint(model, [x, y] ∈ MOI.Complements(2)) + inner = MOI.Utilities.Model{Float64}() + set_optimizer(model, () -> make_opt(inner)) + MOI.Utilities.attach_optimizer(model) + S = MathOptComplements.ComplementsWithSetType{MOI.Interval{Float64}} + @test MOI.get(inner, MOI.NumberOfConstraints{MOI.VectorOfVariables,S}()) == + 0 + return +end + +function _test_opt_Vertical(make_opt) + model = Model() + # Case 1: LHS is already a variable (do nothing) + @variable(model, x1) + @variable(model, 0.0 <= y1) + @constraint(model, [x1, y1] ∈ MOI.Complements(2)) + # Case 2: RHS is unbounded (convert LHS to equality) + @variable(model, x2) + @variable(model, y2) + @constraint(model, [1.0*x2, y2] ∈ MOI.Complements(2)) + # Case 3: LHS is a ScalarAffineFunction with a single variable + @variable(model, x3) + @variable(model, 0.0 <= y3) + @constraint(model, [1.0*x3, y3] ∈ MOI.Complements(2)) + set_optimizer(model, () -> make_opt(Ipopt.Optimizer())) + MOI.Utilities.attach_optimizer(model) + return +end + +function _test_opt_Vertical_errors(make_opt) + model = Model() + @variable(model, 0.0 <= x) + @variable(model, 0.0 <= y) + @constraint(model, [x, 1.0*y + x] ∈ MOI.Complements(2)) + set_optimizer(model, () -> make_opt(Ipopt.Optimizer())) + @test_throws Exception MOI.Utilities.attach_optimizer(model) + return +end + +_test_opt_Instances(make_opt) = Instances.runtests(make_opt) + +function _test_relaxation(relax) + model = Model() + @variable(model, 0 <= x <= 1) + @variable(model, μ) + @constraint(model, 1 - μ == 0.0) + @constraint(model, μ ⟂ x) + vars = [x, μ] + sol = [0.0, 1.0] + set_optimizer( + model, + () -> MathOptComplements.Optimizer( + MOI.instantiate(Ipopt.Optimizer; with_cache_type = Float64), + ), + ) + set_attribute( + model, + MathOptComplements.DefaultComplementarityReformulation(), + relax, + ) + set_attribute(model, "bound_relax_factor", 0.0) + set_silent(model) + optimize!(model) + inner = backend(model).optimizer.model.model + if relax isa MathOptComplements.ScholtesRelaxation + F = MOI.ScalarQuadraticFunction{Float64} + G = MOI.ScalarNonlinearFunction + else + G = MOI.ScalarQuadraticFunction{Float64} + F = MOI.ScalarNonlinearFunction + end + @test MOI.get(inner, MOI.NumberOfConstraints{F,MOI.LessThan{Float64}}()) > 0 + @test MOI.get(inner, MOI.NumberOfConstraints{G,MOI.LessThan{Float64}}()) == + 0 + @test is_solved_and_feasible(model) + @test value.(vars) ≈ sol atol=1e-7 + return +end + +function test_relaxation_scholtes() + return _test_relaxation(MathOptComplements.ScholtesRelaxation(0.0)) +end + +function test_relaxation_fisher_burmeister() + return _test_relaxation( + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + ) +end + +function test_simple_ncp() + for relax in [ + MathOptComplements.ScholtesRelaxation(0.0), + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + MathOptComplements.LiuFukushimaRelaxation(1e-8), + MathOptComplements.KanzowSchwarzRelaxation(1e-8), + ] + model = Model() + @variable(model, y <= 0) + @constraint(model, y + 1 ⟂ y) + set_optimizer( + model, + () -> MathOptComplements.Optimizer(Ipopt.Optimizer()), + ) + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + relax, + ) + set_attribute(model, "bound_relax_factor", 0.0) + set_silent(model) + optimize!(model) + @test is_solved_and_feasible(model) + @test value(y) ≈ -1.0 atol=1e-7 + end + return +end + +function test_simple_lp_3() + for relax in [ + MathOptComplements.ScholtesRelaxation(0.0), + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + MathOptComplements.LiuFukushimaRelaxation(1e-8), + MathOptComplements.KanzowSchwarzRelaxation(1e-8), + ] + model = Model() + @variable(model, 0.0 <= x[1:2]) + @variable(model, 0.0 <= z[1:2]) + @variable(model, y) + @constraint(model, -z[1] + y == 0.0) + @constraint(model, -1.0 - z[2] + y == 0.0) + @constraint(model, x[1] + x[2] == 1.0) + @constraint(model, z[1] ⟂ x[1]) + @constraint(model, z[2] ⟂ x[2]) + set_optimizer( + model, + () -> MathOptComplements.Optimizer(Ipopt.Optimizer()), + ) + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + relax, + ) + set_attribute(model, "bound_relax_factor", 0.0) + set_silent(model) + optimize!(model) + @test is_solved_and_feasible(model) + @test value([x; z; y]) ≈ [0.0, 1.0, 1.0, 0.0, 1.0] atol=1e-7 + end + return +end + +function test_fletcher_leyffer_ex1_model() + for relax in [ + MathOptComplements.ScholtesRelaxation(0.0), + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + MathOptComplements.LiuFukushimaRelaxation(1e-8), + MathOptComplements.KanzowSchwarzRelaxation(1e-8), + ] + model = Model() + @variable(model, z[1:2]) + set_lower_bound(z[2], 0) + @objective(model, Min, (z[1] - 1)^2 + z[2]^2) + @constraint(model, [z[2] - z[1], z[2]] ∈ MOI.Complements(2)) + set_optimizer( + model, + () -> MathOptComplements.Optimizer(Ipopt.Optimizer()), + ) + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + relax, + ) + set_attribute(model, "bound_relax_factor", 0.0) + set_silent(model) + optimize!(model) + @test is_solved_and_feasible(model) + @test objective_value(model) ≈ 0.5 atol=1e-7 + @test value.(model[:z]) ≈ [0.5, 0.5] atol=1e-7 + end + return +end + +function test_reformulation_doesnt_error() + for relax in [ + MathOptComplements.ScholtesRelaxation(0.0), + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + MathOptComplements.LiuFukushimaRelaxation(1e-8), + MathOptComplements.KanzowSchwarzRelaxation(1e-8), + ] + model = Model(() -> MathOptComplements.Optimizer(Ipopt.Optimizer())) + # Case 1: Complementarity defined as lower-bound on RHS + @variable(model, x1) + @variable(model, 0.0 <= y1) + @constraint(model, [x1, y1] ∈ MOI.Complements(2)) + # Case 2: Complementarity defined as upper-bound on RHS + @variable(model, x2) + @variable(model, y2 <= 1.0) + @constraint(model, [x2, y2] ∈ MOI.Complements(2)) + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + relax, + ) + MOI.Utilities.attach_optimizer(model) + end + return +end + +function test_reformulation_errors_1() + for relax in [ + MathOptComplements.ScholtesRelaxation(0.0), + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + MathOptComplements.LiuFukushimaRelaxation(1e-8), + MathOptComplements.KanzowSchwarzRelaxation(1e-8), + ] + model = Model(() -> MathOptComplements.Optimizer(Ipopt.Optimizer())) + @variable(model, 1.0 <= x1) + @variable(model, 0.0 <= y1) + @constraint(model, [x1, y1] ∈ MOI.Complements(2)) + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + relax, + ) + # FIXME(odow): why is this? + @test_throws Exception MOI.Utilities.attach_optimizer(model) + end + return +end + +function test_reformulation_errors_2() + for relax in [ + MathOptComplements.ScholtesRelaxation(0.0), + MathOptComplements.FischerBurmeisterRelaxation(1e-8), + MathOptComplements.LiuFukushimaRelaxation(1e-8), + MathOptComplements.KanzowSchwarzRelaxation(1e-8), + ] + model = Model(() -> MathOptComplements.Optimizer(Ipopt.Optimizer())) + @variable(model, x1 <= 1.0) + @variable(model, y1 <= 0.0) + @constraint(model, [x1, y1] ∈ MOI.Complements(2)) + MOI.set( + model, + MathOptComplements.DefaultComplementarityReformulation(), + relax, + ) + # FIXME(odow): why is this? + @test_throws Exception MOI.Utilities.attach_optimizer(model) + end + return +end + +function test_nonlinear_reformulation() + model = Model() + @variable(model, x >= 0.0) + @variable(model, y >= 0.0) + @objective(model, Min, x^2 + y^2 - 4*x*y) + @constraint(model, [sin(x), y] ∈ MOI.Complements(2)) + inner = MOI.Utilities.Model{Float64}() + set_optimizer(model, () -> MathOptComplements.Optimizer(inner)) + MOI.Utilities.attach_optimizer(model) + expected = Model() + @variable(expected, x >= 0.0) + @variable(expected, y >= 0.0) + @variable(expected, slack >= 0) + @objective(expected, Min, x^2 + y^2 - 4*x*y) + # Build complementarity constraints with nonlinear expression + @constraint(expected, sin(x) == slack) + @constraint(expected, slack * y <= 0.0) + MOI.Bridges._test_structural_identical( + unsafe_backend(model).model, + backend(expected), + ) + return +end + +function test_reformulation_fletcher_leyffer_ex1() + model = Model() + @variable(model, z[1:2]) + set_lower_bound(z[2], 0) + @objective(model, Min, (z[1] - 1)^2 + z[2]^2) + @constraint(model, [z[2] - z[1], z[2]] ∈ MOI.Complements(2)) + inner = MOI.Utilities.Model{Float64}() + set_optimizer(model, () -> MathOptComplements.Optimizer(inner)) + MOI.Utilities.attach_optimizer(model) + expected = Model() + @variable(expected, z[1:2]) + @variable(expected, slack >= 0.0) + set_lower_bound(z[2], 0) + @objective(expected, Min, (z[1] - 1)^2 + z[2]^2) + @constraint(expected, z[2] - z[1] == slack) + @constraint(expected, slack * z[2] <= 0) + MOI.Bridges._test_structural_identical( + unsafe_backend(model).model, + backend(expected), + ) + return +end + +end # TestOptimizer + +TestOptimizer.runtests()