Add ObjectiveAbsoluteTolerance and use in EpsilonConstraint

README.md

@@ -51,6 +51,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.
 
+ * `MOO.ObjectiveAbsoluteTolerance(index::Int)`
  * `MOO.ObjectivePriority(index::Int)`
  * `MOO.ObjectiveRelativeTolerance(index::Int)`
  * `MOO.ObjectiveWeight(index::Int)`

src/MOO.jl

@@ -194,6 +194,21 @@ end
 
 default(::ObjectiveRelativeTolerance) = 0.0
 
+"""
+    ObjectiveAbsoluteTolerance(index::Int) <: AbstractAlgorithmAttribute -> Float64
+
+Assign a `Float64` tolerance to objective number `index`. This is most commonly
+used to constrain an objective to a range in absolute terms to the optimal
+objective value of that objective.
+
+Defaults to `0.0`.
+"""
+struct ObjectiveAbsoluteTolerance <: AbstractAlgorithmAttribute
+    index::Int
+end
+
+default(::ObjectiveAbsoluteTolerance) = 0.0
+
 ### RawOptimizerAttribute
 
 function MOI.supports(model::Optimizer, attr::MOI.RawOptimizerAttribute)

src/algorithms/EpsilonConstraint.jl

@@ -11,17 +11,25 @@ bi-objective programs.
 
 ## Supported optimizer attributes
 
- * `MOO.SolutionLimit()`: 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)`.
+ * `MOO.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.
+
+ * `MOO.SolutionLimit()`: if `MOO.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)`.
 """
 mutable struct EpsilonConstraint <: AbstractAlgorithm
     solution_limit::Union{Nothing,Int}
+    atol::Union{Nothing,Float64}
 
-    EpsilonConstraint() = new(nothing)
+    EpsilonConstraint() = new(nothing, nothing)
 end
 
 default(::EpsilonConstraint, ::SolutionLimit) = 100
@@ -37,6 +45,24 @@ function MOI.get(alg::EpsilonConstraint, attr::SolutionLimit)
     return something(alg.solution_limit, default(alg, attr))
 end
 
+MOI.supports(::EpsilonConstraint, ::ObjectiveAbsoluteTolerance) = true
+
+function MOI.set(
+    alg::EpsilonConstraint,
+    attr::ObjectiveAbsoluteTolerance,
+    value,
+)
+    if attr.index == 1
+        alg.atol = value
+    end
+    return
+end
+
+function MOI.get(alg::EpsilonConstraint, attr::ObjectiveAbsoluteTolerance)
+    @assert attr.index == 1
+    return something(alg.atol, default(alg, attr))
+end
+
 function optimize_multiobjective!(
     algorithm::EpsilonConstraint,
     model::Optimizer,
@@ -55,8 +81,11 @@ 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
-    n_points = MOI.get(algorithm, SolutionLimit())
-    ε = abs(right - left) / (n_points - 1)
+    ε = MOI.get(algorithm, ObjectiveAbsoluteTolerance(1))
+    if iszero(ε)
+        n_points = MOI.get(algorithm, SolutionLimit())
+        ε = abs(right - left) / (n_points - 1)
+    end
     solutions = SolutionPoint[]
     f1, f2 = MOI.Utilities.eachscalar(model.f)
     MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
@@ -66,10 +95,10 @@ function optimize_multiobjective!(
         MOI.LessThan{Float64},
         MOI.GreaterThan{Float64},
     )
-    ci = MOI.add_constraint(model, f1, SetType(left + ε))
+    ci = MOI.add_constraint(model, f1, SetType(left))
     variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
-    for i in 1:n_points
-        rhs = left + (i - 1) * ε
+    rhs = left
+    while rhs <= right + ε / 2
         MOI.set(model, MOI.ConstraintSet(), ci, SetType(rhs))
         MOI.optimize!(model.inner)
         if MOI.get(model.inner, MOI.TerminationStatus()) != MOI.OPTIMAL
@@ -83,6 +112,7 @@ function optimize_multiobjective!(
         if isempty(solutions) || !(Y ≈ solutions[end].y)
             push!(solutions, SolutionPoint(X, Y))
         end
+        rhs += ε
     end
     MOI.delete(model, ci)
     return MOI.OPTIMAL, unique(solutions)

test/algorithms/EpsilonConstraint.jl

@@ -66,6 +66,50 @@ function test_biobjective_knapsack()
     return
 end
 
+function test_biobjective_knapsack_atol()
+    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 = MOO.Optimizer(HiGHS.Optimizer)
+    MOI.set(model, MOO.Algorithm(), MOO.EpsilonConstraint())
+    MOI.set(model, MOO.ObjectiveAbsoluteTolerance(1), 1.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(
+        [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()) == 9
+    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]