Add Sandwiching

.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:

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"]

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.

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

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