diff --git a/README.md b/README.md index 177fc1c..a79e7be 100644 --- a/README.md +++ b/README.md @@ -54,6 +54,7 @@ Each algorithm supports only a subset of the attributes. Consult the algorithm's docstring for details on which attributes it supports, and how it uses them in the solution process. + * `MOA.EpsilonConstraintStep()` * `MOA.ObjectiveAbsoluteTolerance(index::Int)` * `MOA.ObjectivePriority(index::Int)` * `MOA.ObjectiveRelativeTolerance(index::Int)` diff --git a/src/MultiObjectiveAlgorithms.jl b/src/MultiObjectiveAlgorithms.jl index 0ed3748..9015b85 100644 --- a/src/MultiObjectiveAlgorithms.jl +++ b/src/MultiObjectiveAlgorithms.jl @@ -249,6 +249,17 @@ end default(::ObjectiveAbsoluteTolerance) = 0.0 +""" + EpsilonConstraintStep <: AbstractAlgorithmAttribute -> Float64 + +The step `ε` to use in epsilon-constraint methods. + +Defaults to `1.0`. +""" +struct EpsilonConstraintStep <: AbstractAlgorithmAttribute end + +default(::EpsilonConstraintStep) = 1.0 + ### RawOptimizerAttribute function MOI.supports(model::Optimizer, attr::MOI.RawOptimizerAttribute) diff --git a/src/algorithms/EpsilonConstraint.jl b/src/algorithms/EpsilonConstraint.jl index 04c2ac3..ddd14c8 100644 --- a/src/algorithms/EpsilonConstraint.jl +++ b/src/algorithms/EpsilonConstraint.jl @@ -11,19 +11,17 @@ bi-objective programs. ## Supported optimizer attributes - * `MOA.ObjectiveAbsoluteTolerance(1)`: if set, `EpsilonConstraint` uses this - tolerance as the epsilon by which it partitions the first-objective's space. - If the objective is a pure integer program, set the tolerance to `1` to - enumerate all non-dominated solutions. Note that you can set only the - tolerance for the first objective index; the tolerances for other objective - indices are ignored. - - * `MOA.SolutionLimit()`: if `MOA.ObjectiveAbsoluteTolerance(1)` is not set - then, with a slight abuse of notation, `EpsilonConstraint` divides the width - of the first-objective's domain in objective space by `SolutionLimit` to - obtain the epsilon to use when iterating. Thus, there can be at most - `SolutionLimit` solutions returned, but there may be fewer. If no value is - set, the default is `100`, instead of the typical `default(::SolutionLimit)`. + * `MOA.EpsilonConstraintStep()`: `EpsilonConstraint` uses this value + as the epsilon by which it partitions the first-objective's space. The + default is `1`, so that for a pure integer program this algorithm will + enumerate all non-dominated solutions. + + * `MOA.SolutionLimit()`: if this attribute is set then, instead of using the + `MOA.EpsilonConstraintStep`, with a slight abuse of notation, + `EpsilonConstraint` divides the width of the first-objective's domain in + objective space by `SolutionLimit` to obtain the epsilon to use when + iterating. Thus, there can be at most `SolutionLimit` solutions returned, but + there may be fewer. """ mutable struct EpsilonConstraint <: AbstractAlgorithm solution_limit::Union{Nothing,Int} @@ -32,8 +30,6 @@ mutable struct EpsilonConstraint <: AbstractAlgorithm EpsilonConstraint() = new(nothing, nothing) end -default(::EpsilonConstraint, ::SolutionLimit) = 100 - MOI.supports(::EpsilonConstraint, ::SolutionLimit) = true function MOI.set(alg::EpsilonConstraint, ::SolutionLimit, value) @@ -45,24 +41,30 @@ function MOI.get(alg::EpsilonConstraint, attr::SolutionLimit) return something(alg.solution_limit, default(alg, attr)) end -MOI.supports(::EpsilonConstraint, ::ObjectiveAbsoluteTolerance) = true +MOI.supports(::EpsilonConstraint, ::EpsilonConstraintStep) = true -function MOI.set( - alg::EpsilonConstraint, - attr::ObjectiveAbsoluteTolerance, - value, -) - if attr.index == 1 - alg.atol = value - end +function MOI.set(alg::EpsilonConstraint, ::EpsilonConstraintStep, value) + alg.atol = value return end -function MOI.get(alg::EpsilonConstraint, attr::ObjectiveAbsoluteTolerance) - @assert attr.index == 1 +function MOI.get(alg::EpsilonConstraint, attr::EpsilonConstraintStep) return something(alg.atol, default(alg, attr)) end +MOI.supports(::EpsilonConstraint, ::ObjectiveAbsoluteTolerance) = true + +function MOI.set(alg::EpsilonConstraint, ::ObjectiveAbsoluteTolerance, value) + @warn("This attribute is deprecated. Use `EpsilonConstraintStep` instead.") + MOI.set(alg, EpsilonConstraintStep(), value) + return +end + +function MOI.get(alg::EpsilonConstraint, ::ObjectiveAbsoluteTolerance) + @warn("This attribute is deprecated. Use `EpsilonConstraintStep` instead.") + return MOI.get(alg, EpsilonConstraintStep()) +end + function optimize_multiobjective!( algorithm::EpsilonConstraint, model::Optimizer, @@ -85,9 +87,9 @@ function optimize_multiobjective!( a, b = solution_1[1].y[1], solution_2[1].y[1] left, right = min(a, b), max(a, b) # Compute the epsilon that we will be incrementing by each iteration - ε = MOI.get(algorithm, ObjectiveAbsoluteTolerance(1)) - if iszero(ε) - n_points = MOI.get(algorithm, SolutionLimit()) + ε = MOI.get(algorithm, EpsilonConstraintStep()) + n_points = MOI.get(algorithm, SolutionLimit()) + if n_points != default(algorithm, SolutionLimit()) ε = abs(right - left) / (n_points - 1) end solutions = SolutionPoint[] @@ -95,19 +97,18 @@ function optimize_multiobjective!( MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f2)}(), f2) # Add epsilon constraint sense = MOI.get(model.inner, MOI.ObjectiveSense()) - SetType = ifelse( - sense == MOI.MIN_SENSE, - MOI.LessThan{Float64}, - MOI.GreaterThan{Float64}, - ) - ci = MOI.add_constraint(model, f1, SetType(left)) variables = MOI.get(model.inner, MOI.ListOfVariableIndices()) - rhs = left - while rhs <= right + ε / 2 - MOI.set(model, MOI.ConstraintSet(), ci, SetType(rhs)) + SetType, bound, direction = if sense == MOI.MIN_SENSE + MOI.LessThan{Float64}, right, -1.0 + else + MOI.GreaterThan{Float64}, left, 1.0 + end + ci = MOI.add_constraint(model, f1, SetType(bound)) + while true + MOI.set(model, MOI.ConstraintSet(), ci, SetType(bound)) MOI.optimize!(model.inner) if MOI.get(model.inner, MOI.TerminationStatus()) != MOI.OPTIMAL - return MOI.OTHER_ERROR, nothing + break end X = Dict{MOI.VariableIndex,Float64}( x => MOI.get(model.inner, MOI.VariablePrimal(), x) for @@ -117,7 +118,7 @@ function optimize_multiobjective!( if isempty(solutions) || !(Y ≈ solutions[end].y) push!(solutions, SolutionPoint(X, Y)) end - rhs += ε + bound = Y[1] + direction * ε end MOI.delete(model, ci) return MOI.OPTIMAL, filter_nondominated(sense, solutions) diff --git a/test/algorithms/EpsilonConstraint.jl b/test/algorithms/EpsilonConstraint.jl index fdce2d6..4123add 100644 --- a/test/algorithms/EpsilonConstraint.jl +++ b/test/algorithms/EpsilonConstraint.jl @@ -29,6 +29,7 @@ function test_biobjective_knapsack() w = [80, 87, 68, 72, 66, 77, 99, 85, 70, 93, 98, 72, 100, 89, 67, 86, 91] model = MOA.Optimizer(HiGHS.Optimizer) MOI.set(model, MOA.Algorithm(), MOA.EpsilonConstraint()) + MOI.set(model, MOA.SolutionLimit(), 100) MOI.set(model, MOI.Silent(), true) x = MOI.add_variables(model, length(w)) MOI.add_constraint.(model, x, MOI.ZeroOne()) @@ -72,7 +73,6 @@ function test_biobjective_knapsack_atol() w = [80, 87, 68, 72, 66, 77, 99, 85, 70, 93, 98, 72, 100, 89, 67, 86, 91] model = MOA.Optimizer(HiGHS.Optimizer) MOI.set(model, MOA.Algorithm(), MOA.EpsilonConstraint()) - MOI.set(model, MOA.ObjectiveAbsoluteTolerance(1), 1.0) MOI.set(model, MOI.Silent(), true) x = MOI.add_variables(model, length(w)) MOI.add_constraint.(model, x, MOI.ZeroOne()) @@ -110,6 +110,45 @@ function test_biobjective_knapsack_atol() return end +function test_biobjective_knapsack_atol_large() + p1 = [77, 94, 71, 63, 96, 82, 85, 75, 72, 91, 99, 63, 84, 87, 79, 94, 90] + p2 = [65, 90, 90, 77, 95, 84, 70, 94, 66, 92, 74, 97, 60, 60, 65, 97, 93] + w = [80, 87, 68, 72, 66, 77, 99, 85, 70, 93, 98, 72, 100, 89, 67, 86, 91] + model = MOA.Optimizer(HiGHS.Optimizer) + MOI.set(model, MOA.Algorithm(), MOA.EpsilonConstraint()) + @test MOI.supports(model, MOA.EpsilonConstraintStep()) + MOI.set(model, MOA.EpsilonConstraintStep(), 10.0) + MOI.set(model, MOI.Silent(), true) + x = MOI.add_variables(model, length(w)) + MOI.add_constraint.(model, x, MOI.ZeroOne()) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + f = MOI.Utilities.operate( + vcat, + Float64, + [sum(1.0 * p[i] * x[i] for i in 1:length(w)) for p in [p1, p2]]..., + ) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.add_constraint( + model, + sum(1.0 * w[i] * x[i] for i in 1:length(w)), + MOI.LessThan(900.0), + ) + MOI.optimize!(model) + results = Dict( + [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17], + [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17], + [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17], + ) + @test MOI.get(model, MOI.ResultCount()) == 3 + for i in 1:MOI.get(model, MOI.ResultCount()) + x_sol = MOI.get(model, MOI.VariablePrimal(i), x) + X = findall(elt -> elt > 0.9, x_sol) + Y = MOI.get(model, MOI.ObjectiveValue(i)) + @test results[round.(Int, Y)] == X + end + return +end + function test_biobjective_knapsack_min() p1 = [77, 94, 71, 63, 96, 82, 85, 75, 72, 91, 99, 63, 84, 87, 79, 94, 90] p2 = [65, 90, 90, 77, 95, 84, 70, 94, 66, 92, 74, 97, 60, 60, 65, 97, 93] @@ -160,6 +199,7 @@ function test_biobjective_knapsack_min_solution_limit() w = [80, 87, 68, 72, 66, 77, 99, 85, 70, 93, 98, 72, 100, 89, 67, 86, 91] model = MOA.Optimizer(HiGHS.Optimizer) MOI.set(model, MOA.Algorithm(), MOA.EpsilonConstraint()) + @test MOI.supports(model, MOA.SolutionLimit()) MOI.set(model, MOA.SolutionLimit(), 3) MOI.set(model, MOI.Silent(), true) x = MOI.add_variables(model, length(w)) @@ -178,17 +218,10 @@ function test_biobjective_knapsack_min_solution_limit() ) MOI.optimize!(model) results = Dict( - [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17], - [949, 915] => [1, 2, 5, 6, 8, 9, 10, 11, 15, 16, 17], - [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17], [943, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17], - [936, 942] => [1, 2, 3, 5, 6, 10, 11, 12, 15, 16, 17], - [935, 947] => [2, 5, 6, 8, 9, 10, 11, 12, 15, 16, 17], - [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17], - [927, 972] => [2, 3, 5, 6, 8, 9, 10, 11, 12, 16, 17], [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17], ) - @test MOI.get(model, MOI.ResultCount()) == 3 + @test MOI.get(model, MOI.ResultCount()) == 2 for i in 1:MOI.get(model, MOI.ResultCount()) x_sol = MOI.get(model, MOI.VariablePrimal(i), x) X = findall(elt -> elt > 0.9, x_sol) @@ -230,6 +263,15 @@ function test_unbounded() return end +function test_deprecated() + model = MOA.Optimizer(HiGHS.Optimizer) + MOI.set(model, MOA.Algorithm(), MOA.EpsilonConstraint()) + @test MOI.supports(model, MOA.ObjectiveAbsoluteTolerance(1)) + @test_logs (:warn,) MOI.set(model, MOA.ObjectiveAbsoluteTolerance(1), 1.0) + @test_logs (:warn,) MOI.get(model, MOA.ObjectiveAbsoluteTolerance(1)) + return +end + end TestEpsilonConstraint.run_tests()