Add EpsilonConstraint algorithm

README.md

@@ -36,8 +36,9 @@ the choice of solution algorithm.
 
 There are a number of algorithms supported by the algorithms in MOO.
 
- * `MOO.Lexicographic()` [default]
+ * `MOO.EpsilonConstraint()`
  * `MOO.Hierarchical()`
+ * `MOO.Lexicographic()` [default]
  * `MOO.NISE()`
 
 Consult their docstrings for details.
@@ -50,7 +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.SolutionLimit()`
  * `MOO.ObjectivePriority(index::Int)`
- * `MOO.ObjectiveWeight(index::Int)`
  * `MOO.ObjectiveRelativeTolerance(index::Int)`
+ * `MOO.ObjectiveWeight(index::Int)`
+ * `MOO.SolutionLimit()`

src/MOO.jl

@@ -124,6 +124,8 @@ A super-type for MOO-specific optimizer attributes.
 """
 abstract type AbstractAlgorithmAttribute <: MOI.AbstractOptimizerAttribute end
 
+default(::AbstractAlgorithm, attr::AbstractAlgorithmAttribute) = default(attr)
+
 function MOI.supports(model::Optimizer, attr::AbstractAlgorithmAttribute)
     return MOI.supports(model.algorithm, attr)
 end
@@ -192,13 +194,6 @@ end
 
 default(::ObjectiveRelativeTolerance) = 0.0
 
-function _append_default(attr::AbstractAlgorithmAttribute, x::Vector)
-    for _ in (1+length(x)):attr.index
-        push!(x, default(attr))
-    end
-    return
-end
-
 ### RawOptimizerAttribute
 
 function MOI.supports(model::Optimizer, attr::MOI.RawOptimizerAttribute)

src/algorithms/EpsilonConstraint.jl

@@ -0,0 +1,89 @@
+#  Copyright 2019, Oscar Dowson and contributors
+#  This Source Code Form is subject to the terms of the Mozilla Public License,
+#  v.2.0. If a copy of the MPL was not distributed with this file, You can
+#  obtain one at http://mozilla.org/MPL/2.0/.
+
+"""
+    EpsilonConstraint()
+
+`EpsilonConstraint` implements the epsilon-constraint algorithm for
+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)`.
+"""
+mutable struct EpsilonConstraint <: AbstractAlgorithm
+    solution_limit::Union{Nothing,Int}
+
+    EpsilonConstraint() = new(nothing)
+end
+
+default(::EpsilonConstraint, ::SolutionLimit) = 100
+
+MOI.supports(::EpsilonConstraint, ::SolutionLimit) = true
+
+function MOI.set(alg::EpsilonConstraint, ::SolutionLimit, value)
+    alg.solution_limit = value
+    return
+end
+
+function MOI.get(alg::EpsilonConstraint, attr::SolutionLimit)
+    return something(alg.solution_limit, default(alg, attr))
+end
+
+function optimize_multiobjective!(
+    algorithm::EpsilonConstraint,
+    model::Optimizer,
+)
+    if MOI.output_dimension(model.f) != 2
+        error("EpsilonConstraint requires exactly two objectives")
+    end
+    # Compute the bounding box ofthe objectives using Hierarchical().
+    alg = Hierarchical()
+    MOI.set.(Ref(alg), ObjectivePriority.(1:2), [1, 0])
+    status, solution_1 = optimize_multiobjective!(alg, model)
+    @assert status == MOI.OPTIMAL
+    MOI.set(alg, ObjectivePriority(2), 2)
+    status, solution_2 = optimize_multiobjective!(alg, model)
+    @assert status == MOI.OPTIMAL
+    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)
+    solutions = ParetoSolution[]
+    f1, f2 = MOI.Utilities.eachscalar(model.f)
+    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
+    # Add epsilon constraint
+    SetType = ifelse(
+        MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE,
+        MOI.LessThan{Float64},
+        MOI.GreaterThan{Float64},
+    )
+    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) * ε
+        MOI.set(model, MOI.ConstraintSet(), ci, SetType(rhs))
+        MOI.optimize!(model.inner)
+        if MOI.get(model.inner, MOI.TerminationStatus()) != MOI.OPTIMAL
+            return MOI.OTHER_ERROR, nothing
+        end
+        X = Dict{MOI.VariableIndex,Float64}(
+            x => MOI.get(model.inner, MOI.VariablePrimal(), x) for
+            x in variables
+        )
+        Y = MOI.Utilities.eval_variables(x -> X[x], model.f)
+        if isempty(solutions) || !(Y ≈ solutions[end].y)
+            push!(solutions, ParetoSolution(X, Y))
+        end
+    end
+    MOI.delete(model, ci)
+    return MOI.OPTIMAL, unique(solutions)
+end

src/algorithms/Hierarchical.jl

@@ -36,35 +36,46 @@ end
 MOI.supports(::Hierarchical, ::ObjectivePriority) = true
 
 function MOI.get(alg::Hierarchical, attr::ObjectivePriority)
-    return get(alg.priorities, attr.index, default(attr))
+    return get(alg.priorities, attr.index, default(alg, attr))
+end
+
+function _append_default(
+    alg::Hierarchical,
+    attr::AbstractAlgorithmAttribute,
+    x::Vector,
+)
+    for _ in (1+length(x)):attr.index
+        push!(x, default(alg, attr))
+    end
+    return
 end
 
 function MOI.set(alg::Hierarchical, attr::ObjectivePriority, value)
-    _append_default(attr, alg.priorities)
+    _append_default(alg, attr, alg.priorities)
     alg.priorities[attr.index] = value
     return
 end
 
 MOI.supports(::Hierarchical, ::ObjectiveWeight) = true
 
 function MOI.get(alg::Hierarchical, attr::ObjectiveWeight)
-    return get(alg.weights, attr.index, default(attr))
+    return get(alg.weights, attr.index, default(alg, attr))
 end
 
 function MOI.set(alg::Hierarchical, attr::ObjectiveWeight, value)
-    _append_default(attr, alg.weights)
+    _append_default(alg, attr, alg.weights)
     alg.weights[attr.index] = value
     return
 end
 
 MOI.supports(::Hierarchical, ::ObjectiveRelativeTolerance) = true
 
 function MOI.get(alg::Hierarchical, attr::ObjectiveRelativeTolerance)
-    return get(alg.rtol, attr.index, default(attr))
+    return get(alg.rtol, attr.index, default(alg, attr))
 end
 
 function MOI.set(alg::Hierarchical, attr::ObjectiveRelativeTolerance, value)
-    _append_default(attr, alg.rtol)
+    _append_default(alg, attr, alg.rtol)
     alg.rtol[attr.index] = value
     return
 end
@@ -80,13 +91,13 @@ function optimize_multiobjective!(algorithm::Hierarchical, model::Optimizer)
     variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
     # Find list of objectives with same priority
     constraints = Any[]
-    objective_subsets = _sorted_priorities([
-        MOI.get(algorithm, ObjectivePriority(i)) for i in 1:N
-    ])
+    priorities = [MOI.get(algorithm, ObjectivePriority(i)) for i in 1:N]
+    weights = [MOI.get(algorithm, ObjectiveWeight(i)) for i in 1:N]
+    objective_subsets = _sorted_priorities(priorities)
     for (round, indices) in enumerate(objective_subsets)
         # Solve weighted sum
         new_vector_f = objectives[indices]
-        new_f = _scalarise(new_vector_f, algorithm.weights[indices])
+        new_f = _scalarise(new_vector_f, weights[indices])
         MOI.set(model.inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f)
         MOI.optimize!(model.inner)
         if MOI.get(model.inner, MOI.TerminationStatus()) != MOI.OPTIMAL

src/algorithms/Lexicographic.jl

@@ -22,11 +22,13 @@ end
 MOI.supports(::Lexicographic, ::ObjectiveRelativeTolerance) = true
 
 function MOI.get(alg::Lexicographic, attr::ObjectiveRelativeTolerance)
-    return get(alg.rtol, attr.index, default(attr))
+    return get(alg.rtol, attr.index, default(alg, attr))
 end
 
 function MOI.set(alg::Lexicographic, attr::ObjectiveRelativeTolerance, value)
-    _append_default(attr, alg.rtol)
+    for _ in (1+length(alg.rtol)):attr.index
+        push!(alg.rtol, default(alg, attr))
+    end
     alg.rtol[attr.index] = value
     return
 end

src/algorithms/NISE.jl

@@ -29,7 +29,9 @@ function MOI.set(alg::NISE, ::SolutionLimit, value)
     return
 end
 
-MOI.get(alg::NISE, ::SolutionLimit) = alg.solution_limit
+function MOI.get(alg::NISE, attr::SolutionLimit)
+    return something(alg.solution_limit, default(alg, attr))
+end
 
 function _solve_weighted_sum(model::Optimizer, alg::NISE, weight::Float64)
     return _solve_weighted_sum(model, alg, [weight, 1 - weight])
@@ -68,7 +70,7 @@ function optimize_multiobjective!(algorithm::NISE, model::Optimizer)
     if !(solutions[0.0] ≈ solutions[1.0])
         push!(queue, (0.0, 1.0))
     end
-    limit = something(algorithm.solution_limit, default(SolutionLimit()))
+    limit = MOI.get(algorithm, SolutionLimit())
     while length(queue) > 0 && length(solutions) < limit
         (a, b) = popfirst!(queue)
         y_d = solutions[a].y .- solutions[b].y