From 25e4914a65a9325cd77f85a73b9a2d6ea65d4f06 Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 6 Feb 2023 12:38:59 +1300 Subject: [PATCH 1/2] Add EpsilonConstraint algorithm --- README.md | 7 +- src/MOO.jl | 9 +- src/algorithms/EpsilonConstraint.jl | 89 +++++++++++++++ src/algorithms/Hierarchical.jl | 31 ++++-- src/algorithms/Lexicographic.jl | 6 +- src/algorithms/NISE.jl | 2 +- test/algorithms/EpsilonConstraint.jl | 159 +++++++++++++++++++++++++++ 7 files changed, 280 insertions(+), 23 deletions(-) create mode 100644 src/algorithms/EpsilonConstraint.jl create mode 100644 test/algorithms/EpsilonConstraint.jl diff --git a/README.md b/README.md index b860fde..8b4f784 100644 --- a/README.md +++ b/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()` diff --git a/src/MOO.jl b/src/MOO.jl index 14be38f..c5932ab 100644 --- a/src/MOO.jl +++ b/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) diff --git a/src/algorithms/EpsilonConstraint.jl b/src/algorithms/EpsilonConstraint.jl new file mode 100644 index 0000000..502ef6b --- /dev/null +++ b/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 diff --git a/src/algorithms/Hierarchical.jl b/src/algorithms/Hierarchical.jl index a9d0869..0754d84 100644 --- a/src/algorithms/Hierarchical.jl +++ b/src/algorithms/Hierarchical.jl @@ -36,11 +36,22 @@ 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 @@ -48,11 +59,11 @@ 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 @@ -60,11 +71,11 @@ 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 diff --git a/src/algorithms/Lexicographic.jl b/src/algorithms/Lexicographic.jl index 1e9c636..493dcc1 100644 --- a/src/algorithms/Lexicographic.jl +++ b/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 diff --git a/src/algorithms/NISE.jl b/src/algorithms/NISE.jl index 88c9476..1a5b9eb 100644 --- a/src/algorithms/NISE.jl +++ b/src/algorithms/NISE.jl @@ -68,7 +68,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 = something(algorithm.solution_limit, default(alg, SolutionLimit())) while length(queue) > 0 && length(solutions) < limit (a, b) = popfirst!(queue) y_d = solutions[a].y .- solutions[b].y diff --git a/test/algorithms/EpsilonConstraint.jl b/test/algorithms/EpsilonConstraint.jl new file mode 100644 index 0000000..256bd41 --- /dev/null +++ b/test/algorithms/EpsilonConstraint.jl @@ -0,0 +1,159 @@ +# 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/. + +module TestEpsilonConstraint + +using Test + +import HiGHS +import MOO + +const MOI = MOO.MOI + +function run_tests() + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$name" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +function test_biobjective_knapsack() + 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, 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] + 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.SolutionLimit(), 100) + 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.MIN_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_solution_limit() + 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.SolutionLimit(), 3) + 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()) == 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 + +end + +TestEpsilonConstraint.run_tests() From 7d4a9afed3eb1ee8308a26e63eebf4bdb505a76c Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 6 Feb 2023 13:07:59 +1300 Subject: [PATCH 2/2] Fix tests --- src/algorithms/NISE.jl | 6 ++++-- test/algorithms/NISE.jl | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/algorithms/NISE.jl b/src/algorithms/NISE.jl index 1a5b9eb..7ef5670 100644 --- a/src/algorithms/NISE.jl +++ b/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(alg, 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 diff --git a/test/algorithms/NISE.jl b/test/algorithms/NISE.jl index b8b6ea7..775d11f 100644 --- a/test/algorithms/NISE.jl +++ b/test/algorithms/NISE.jl @@ -26,7 +26,8 @@ end function test_NISE_SolutionLimit() model = MOO.Optimizer(HiGHS.Optimizer) MOI.set(model, MOO.Algorithm(), MOO.NISE()) - @test MOI.get(model, MOO.SolutionLimit()) == nothing + @test MOI.get(model, MOO.SolutionLimit()) == + MOO.default(MOO.SolutionLimit()) MOI.set(model, MOO.SolutionLimit(), 1) @test MOI.get(model, MOO.SolutionLimit()) == 1 return