diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8369946..57ef58a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -26,10 +26,10 @@ jobs: - version: '1' # The latest point-release (Windows) os: windows-latest arch: x64 - - version: '1.6' # 1.6 LTS (64-bit Linux) + - version: '1.10' # 1.10 LTS (64-bit Linux) os: ubuntu-latest arch: x64 - - version: '1.6' # 1.6 LTS (32-bit Linux) + - version: '1.10' # 1.10 LTS (32-bit Linux) os: ubuntu-latest arch: x86 steps: diff --git a/Project.toml b/Project.toml index 95437f5..23c221b 100644 --- a/Project.toml +++ b/Project.toml @@ -7,20 +7,28 @@ version = "1.4.3" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +[weakdeps] +Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" + +[extensions] +MultiObjectiveAlgorithmsPolyhedraExt = "Polyhedra" + [compat] Combinatorics = "1" HiGHS = "1" Ipopt = "1" JSON = "0.21" MathOptInterface = "1.19" -Test = "<0.0.1, 1.6" -julia = "1.6" +Polyhedra = "0.8" +Test = "1" +julia = "1.10" [extras] HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" [targets] -test = ["HiGHS", "Ipopt", "JSON", "Test"] +test = ["HiGHS", "Ipopt", "JSON", "Test", "Polyhedra"] diff --git a/README.md b/README.md index 826974c..818662e 100644 --- a/README.md +++ b/README.md @@ -63,6 +63,7 @@ The value must be one of the algorithms supported by MOA: * `MOA.KirlikSayin()` * `MOA.Lexicographic()` [default] * `MOA.RandomWeighting()` + * `MOA.Sandwiching()` * `MOA.TambyVanderpooten()` Consult their docstrings for details. diff --git a/ext/MultiObjectiveAlgorithmsPolyhedraExt.jl b/ext/MultiObjectiveAlgorithmsPolyhedraExt.jl new file mode 100644 index 0000000..b68034f --- /dev/null +++ b/ext/MultiObjectiveAlgorithmsPolyhedraExt.jl @@ -0,0 +1,138 @@ +# 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 MultiObjectiveAlgorithmsPolyhedraExt + +import MathOptInterface as MOI +import MultiObjectiveAlgorithms as MOA +import Polyhedra + +function _halfspaces(IPS::Vector{Vector{Float64}}) + V = Polyhedra.vrep(IPS) + H = Polyhedra.halfspaces(Polyhedra.doubledescription(V)) + return [(-H_i.a, -H_i.β) for H_i in H] +end + +function _distance(w̄, b̄, δ_OPS_optimizer) + y = MOI.get(δ_OPS_optimizer, MOI.ListOfVariableIndices()) + MOI.set( + δ_OPS_optimizer, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w̄, y), 0.0), + ) + MOI.set(δ_OPS_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.optimize!(δ_OPS_optimizer) + return b̄ - MOI.get(δ_OPS_optimizer, MOI.ObjectiveValue()) +end + +function _select_next_halfspace(H, δ_OPS_optimizer) + distances = [_distance(w, b, δ_OPS_optimizer) for (w, b) in H] + index = argmax(distances) + w, b = H[index] + return distances[index], w, b +end + +function MOA.minimize_multiobjective!( + algorithm::MOA.Sandwiching, + model::MOA.Optimizer, +) + @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE + start_time = time() + solutions = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}() + variables = MOI.get(model.inner, MOI.ListOfVariableIndices()) + n = MOI.output_dimension(model.f) + scalars = MOI.Utilities.scalarize(model.f) + status = MOI.OPTIMAL + optimizer = typeof(model.inner.optimizer) + δ_OPS_optimizer = optimizer() + MOI.set(δ_OPS_optimizer, MOI.Silent(), true) + y = MOI.add_variables(δ_OPS_optimizer, n) + anchors = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}() + yI, yUB = zeros(n), zeros(n) + for (i, f_i) in enumerate(scalars) + MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i) + MOI.optimize!(model.inner) + status = MOI.get(model.inner, MOI.TerminationStatus()) + if !MOA._is_scalar_status_optimal(model) + return status, nothing + end + X, Y = MOA._compute_point(model, variables, model.f) + model.ideal_point[i] = Y[i] + yI[i] = Y[i] + anchors[Y] = X + MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.optimize!(model.inner) + status = MOI.get(model.inner, MOI.TerminationStatus()) + if !MOA._is_scalar_status_optimal(model) + MOA._warn_on_nonfinite_anti_ideal(algorithm, MOI.MIN_SENSE, i) + return status, nothing + end + _, Y = MOA._compute_point(model, variables, f_i) + yUB[i] = Y + MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE) + e_i = Float64.(1:n .== i) + MOI.add_constraint( + δ_OPS_optimizer, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(e_i, y), 0.0), + MOI.GreaterThan(yI[i]), + ) + MOI.add_constraint( + δ_OPS_optimizer, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(e_i, y), 0.0), + MOI.LessThan(yUB[i]), + ) + end + IPS = [yUB, keys(anchors)...] + merge!(solutions, anchors) + u = MOI.add_variables(model.inner, n) + u_constraints = [ # u_i >= 0 for all i = 1:n + MOI.add_constraint(model.inner, u_i, MOI.GreaterThan{Float64}(0)) + for u_i in u + ] + f_constraints = [ # f_i + u_i <= yUB_i for all i = 1:n + MOI.Utilities.normalize_and_add_constraint( + model.inner, + scalars[i] + u[i], + MOI.LessThan(yUB[i]), + ) for i in 1:n + ] + H = _halfspaces(IPS) + count = 0 + while !isempty(H) + if MOA._time_limit_exceeded(model, start_time) + status = MOI.TIME_LIMIT + break + end + count += 1 + δ, w, b = _select_next_halfspace(H, δ_OPS_optimizer) + if δ - 1e-3 <= algorithm.precision # added some convergence tolerance + break + end + # would not terminate when precision is set to 0 + new_f = sum(w[i] * (scalars[i] + u[i]) for i in 1:n) # w' * (f(x) + u) + MOI.set(model.inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f) + MOI.optimize!(model.inner) + status = MOI.get(model.inner, MOI.TerminationStatus()) + if !MOA._is_scalar_status_optimal(model) + return status, nothing + end + β̄ = MOI.get(model.inner, MOI.ObjectiveValue()) + X, Y = MOA._compute_point(model, variables, model.f) + solutions[Y] = X + MOI.add_constraint( + δ_OPS_optimizer, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, y), 0.0), + MOI.GreaterThan(β̄), + ) + IPS = push!(IPS, Y) + H = _halfspaces(IPS) + end + MOI.delete.(model.inner, f_constraints) + MOI.delete.(model.inner, u_constraints) + MOI.delete.(model.inner, u) + return status, [MOA.SolutionPoint(X, Y) for (Y, X) in solutions] +end + +end # module MultiObjectiveAlgorithmsPolyhedraExt diff --git a/src/algorithms/Sandwiching.jl b/src/algorithms/Sandwiching.jl new file mode 100644 index 0000000..a53af81 --- /dev/null +++ b/src/algorithms/Sandwiching.jl @@ -0,0 +1,23 @@ +# 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/. + +""" + Sandwiching(precision::Float64) + +An algorithm that implemennts the paper described in Koenen, M., Balvert, M., & Fleuren, H. A. (2023). A Renewed Take on Weighted Sum in Sandwich Algorithms: Modification of the Criterion Space. (Center Discussion Paper; Vol. 2023-012). CentER, Center for Economic Research. + +## Compat + +To use this algorithm you MUST first load the Polyhedra.jl Julia package: + +```julia +import MultiObjectiveAlgorithms as MOA +import Polyhedra +algorithm = MOA.Sandwiching(0.0) +``` +""" +mutable struct Sandwiching <: AbstractAlgorithm + precision::Float64 +end diff --git a/test/algorithms/Sandwiching.jl b/test/algorithms/Sandwiching.jl new file mode 100644 index 0000000..5f6f64b --- /dev/null +++ b/test/algorithms/Sandwiching.jl @@ -0,0 +1,194 @@ +# 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 TestSandwiching + +using Test + +import HiGHS +import MultiObjectiveAlgorithms as MOA +import MultiObjectiveAlgorithms: MOI +import Polyhedra + +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_molp(C, A, b, results, sense) + p = size(C, 1) + m, n = size(A) + model = MOA.Optimizer(HiGHS.Optimizer) + MOI.set(model, MOA.Algorithm(), MOA.Sandwiching(0.0)) + MOI.set(model, MOI.Silent(), true) + x = MOI.add_variables(model, n) + MOI.add_constraint.(model, x, MOI.GreaterThan(0.0)) + for i in 1:m + MOI.add_constraint( + model, + MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(A[i, j], x[j]) for j in 1:n], + 0.0, + ), + MOI.LessThan(b[i]), + ) + end + f = MOI.VectorAffineFunction( + [ + MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for + i in 1:p for j in 1:n + ], + zeros(p), + ) + MOI.set(model, MOI.ObjectiveSense(), sense) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.optimize!(model) + N = MOI.get(model, MOI.ResultCount()) + solutions = sort([ + MOI.get(model, MOI.VariablePrimal(i), x) => + MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N + ]) + @test N == length(results) + for (sol, res) in zip(solutions, results) + x_sol, y_sol = sol + x_res, y_res = res + @test ≈(x_sol, x_res; atol = 1e-6) + @test ≈(y_sol, y_res; atol = 1e-6) + end + return +end + +# From International Doctoral School Algorithmic Decision Theory: MCDA and MOO +# Lecture 2: Multiobjective Linear Programming +# Matthias Ehrgott +# Department of Engineering Science, The University of Auckland, New Zealand +# Laboratoire d’Informatique de Nantes Atlantique, CNRS, Universit´e de Nantes, France +function test_molp_1() + C = Float64[3 1; -1 -2] + A = Float64[0 1; 3 -1] + b = Float64[3, 6] + results = sort([ + [0.0, 0.0] => [0.0, 0.0], + [0.0, 3.0] => [3.0, -6.0], + [3.0, 3.0] => [12.0, -9.0], + ]) + sense = MOI.MIN_SENSE + return _test_molp(C, A, b, results, sense) +end + +# From Civil and Environmental Systems Engineering +# Chapter 5 Exercise 5.A.3 A graphical Interpretation of Noninferiority +function test_molp_2() + C = Float64[3 -2; -1 2] + A = Float64[-4 -8; 3 -6; 4 -2; 1 0; -1 3; -2 4; -6 3] + b = Float64[-8, 6, 14, 6, 15, 18, 9] + results = sort([ + [1.0, 5.0] => [-7.0, 9.0], # not sure about this + [3.0, 6.0] => [-3.0, 9.0], + [4.0, 1.0] => [10.0, -2.0], + [6.0, 5.0] => [8.0, 4.0], + [6.0, 7.0] => [4.0, 8.0], + ]) + sense = MOI.MAX_SENSE + return _test_molp(C, A, b, results, sense) +end + +function test_infeasible() + model = MOA.Optimizer(HiGHS.Optimizer) + MOI.set(model, MOA.Algorithm(), MOA.Sandwiching(0.0)) + MOI.set(model, MOI.Silent(), true) + x = MOI.add_variables(model, 2) + MOI.add_constraint.(model, x, MOI.GreaterThan(0.0)) + MOI.add_constraint(model, 1.0 * x[1] + 1.0 * x[2], MOI.LessThan(-1.0)) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + f = MOI.Utilities.operate(vcat, Float64, 1.0 .* x...) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE + @test MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION + @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION + return +end + +function test_unbounded() + model = MOA.Optimizer(HiGHS.Optimizer) + MOI.set(model, MOA.Algorithm(), MOA.Sandwiching(0.0)) + MOI.set(model, MOI.Silent(), true) + x = MOI.add_variables(model, 2) + MOI.add_constraint.(model, x, MOI.GreaterThan(0.0)) + f = MOI.Utilities.operate(vcat, Float64, 1.0 .* x...) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE + @test MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION + @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION + return +end + +function test_no_bounding_box() + model = MOA.Optimizer(HiGHS.Optimizer) + MOI.set(model, MOA.Algorithm(), MOA.Sandwiching(0.0)) + MOI.set(model, MOI.Silent(), true) + x = MOI.add_variables(model, 2) + MOI.add_constraint.(model, x, MOI.GreaterThan(0.0)) + f = MOI.Utilities.operate(vcat, Float64, 1.0 .* x...) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + @test_logs (:warn,) MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE + @test MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION + @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION + return +end + +function test_time_limit() + p = 3 + n = 10 + W = 2137.0 + C = Float64[ + 566 611 506 180 817 184 585 423 26 317 + 62 84 977 979 874 54 269 93 881 563 + 664 982 962 140 224 215 12 869 332 537 + ] + w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272] + model = MOA.Optimizer(HiGHS.Optimizer) + MOI.set(model, MOA.Algorithm(), MOA.Sandwiching(0.0)) + MOI.set(model, MOI.TimeLimitSec(), 0.0) + MOI.set(model, MOI.Silent(), true) + x = MOI.add_variables(model, n) + MOI.add_constraint.(model, x, MOI.ZeroOne()) + MOI.add_constraint( + model, + MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n], + 0.0, + ), + MOI.LessThan(W), + ) + f = MOI.VectorAffineFunction( + [ + MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j])) + for i in 1:p for j in 1:n + ], + fill(0.0, p), + ) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.TIME_LIMIT + @test MOI.get(model, MOI.ResultCount()) == size(C, 1) # anchor points are already computed when the time limit is checked + return +end + +end # TestSandwiching + +TestSandwiching.run_tests()