Add Sandwiching

Author: kofgokhan

Project.toml

@@ -7,6 +7,12 @@ version = "1.4.3"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 
+[weakdeps]
+Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
+
+[extensions]
+PolyhedraExt = "Polyhedra"
+
 [compat]
 Combinatorics = "1"
 HiGHS = "1"
@@ -21,6 +27,7 @@ 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"]

ext/PolyhedraExt.jl

@@ -0,0 +1,12 @@
+module PolyhedraExt
+
+using MultiObjectiveAlgorithms
+using Polyhedra
+
+function MultiObjectiveAlgorithms._halfspaces(IPS::Vector{Vector{Float64}})
+    V = vrep(IPS)
+    H = halfspaces(doubledescription(V))
+    return [(-H_i.a, -H_i.β) for H_i in H]
+end
+
+end

src/algorithms/Sandwiching.jl

@@ -0,0 +1,145 @@
+function _halfspaces(IPS::Vector{Vector{Float64}})
+    error("MOA.Sandwiching requires Polyhedra.jl to be loaded")
+end
+
+mutable struct Sandwiching <: AbstractAlgorithm
+    precision::Float64
+end
+
+tol = 1e-3
+
+function _compute_anchors(model::Optimizer)
+    anchors = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}()
+    n = MOI.output_dimension(model.f)
+    scalars = MOI.Utilities.scalarize(model.f)
+    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
+    yI, yUB = zeros(n), zeros(n)
+    for (i, f_i) in enumerate(scalars)
+        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i) # ρ * sum(f_j for (j, f_j) in enumerate(model.f) if j != i)
+        MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+        MOI.optimize!(model.inner)
+        # status check
+        X, Y = _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 check
+        _, Y = _compute_point(model, variables, f_i)
+        yUB[i] = Y
+    end
+    MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    return yI, yUB, anchors
+end
+
+@enum DistanceMeasure SUB CUR SOL
+
+function _distance(w̄, b̄, OPS, model)
+    n = MOI.output_dimension(model.f)
+    optimizer = typeof(model.inner.optimizer)
+    δ_optimizer = optimizer()
+    MOI.set(δ_optimizer, MOI.Silent(), true)
+    x = MOI.add_variables(δ_optimizer, n)
+    for (w, b) in OPS
+        MOI.add_constraint(
+            δ_optimizer,
+            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, x), 0.0),
+            MOI.GreaterThan(b),
+        )
+    end
+    MOI.set(
+        δ_optimizer,
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w̄, x), 0.0),
+    )
+    MOI.set(δ_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.optimize!(δ_optimizer)
+    δ = b̄ - MOI.get(δ_optimizer, MOI.ObjectiveValue())
+    return δ
+end
+
+function _select_next_halfspace(H, OPS, model)
+    distances = [_distance(w, b, OPS, model) for (w, b) in H]
+    @info "Distances: $(Dict(zip(H, distances)))"
+    index = argmax(distances)
+    w, b = H[index]
+    return distances[index], w, b
+end
+
+function optimize_multiobjective!(algorithm::Sandwiching, model::Optimizer)
+    @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+    ε = algorithm.precision
+    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)
+    yI, yUB, anchors = _compute_anchors(model)
+    merge!(solutions, anchors)
+    @info "yI: $(yI)"
+    @info "yUB: $(yUB)"
+    IPS = [yUB, keys(anchors)...]
+    OPS = Tuple{Vector{Float64}, Float64}[]
+    for i in 1:n
+        e_i = Float64.(1:n .== i)
+        push!(OPS, (e_i, yI[i])) # e_i' * y >= yI_i
+        push!(OPS, (-e_i, -yUB[i])) # -e_i' * y >= -yUB_i ⟹ e_i' * y <= yUB_i
+    end
+    @info "IPS: $(IPS)"
+    @info "OPS: $(OPS)"
+    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)
+        count += 1
+        @info "-- Iteration #$(count) --"
+        @info "HalfSpaces: $(H)"
+        δ, w, b = _select_next_halfspace(H, OPS, model)
+        @info "Selected halfspace: w: $(w), b: $(b)"
+        @info "δ: $(δ)"
+        if δ - tol > ε # added some convergence tolerance
+                       # 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)
+            β̄ = MOI.get(model.inner, MOI.ObjectiveValue())
+            @info "β̄: $(β̄)"
+            X, Y = _compute_point(model, variables, model.f)
+            @info "Y: $(Y)"
+            solutions[Y] = X
+            push!(OPS, (w, β̄))
+            @info "Added halfspace w: $(w), b: $(β̄) to OPS"
+            IPS = push!(IPS, Y)
+        else
+            break
+        end
+        @info "IPS: $(IPS)"
+        @info "OPS: $(OPS)"
+        H = _halfspaces(IPS)
+        if count == 10
+            break
+        end
+    end
+    MOI.delete.(model.inner, f_constraints)
+    MOI.delete.(model.inner, u_constraints)
+    MOI.delete.(model.inner, u)
+    return MOI.OPTIMAL, [SolutionPoint(X, Y) for (Y, X) in solutions]
+end

test/algorithms/Sandwiching.jl

@@ -0,0 +1,84 @@
+module TestSandwiching
+
+using Test
+
+import HiGHS
+using Polyhedra
+import MultiObjectiveAlgorithms as MOA
+import MultiObjectiveAlgorithms: 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
+
+
+# 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()
+    C = Float64[
+        3 1
+        -1 -2
+    ]
+    p = size(C, 1)
+    A = Float64[
+        0 1
+        3 -1
+    ]
+    m, n = size(A)
+    b = Float64[3, 6]
+    model = MOI.instantiate(; with_bridge_type = Float64) do
+        return MOA.Optimizer(HiGHS.Optimizer)
+    end
+    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(), MOI.MIN_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    N = MOI.get(model, MOI.ResultCount())
+    solutions = reverse([MOI.get(model, MOI.VariablePrimal(i), x) => MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N])
+    results = reverse([
+        [0., 0.] => [0., 0.],
+        [0., 3.] => [3., -6.],
+        [3., 3.] => [12., -9.],
+    ])
+    @test length(solutions) == 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
+end
+
+end
+
+TestSandwiching.run_tests()