|
| 1 | +# Copyright 2019, Oscar Dowson and contributors |
| 2 | +# This Source Code Form is subject to the terms of the Mozilla Public License, |
| 3 | +# v.2.0. If a copy of the MPL was not distributed with this file, You can |
| 4 | +# obtain one at http://mozilla.org/MPL/2.0/. |
| 5 | + |
| 6 | +function _halfspaces(IPS::Vector{Vector{Float64}}) |
| 7 | + V = Polyhedra.vrep(IPS) |
| 8 | + H = Polyhedra.halfspaces(Polyhedra.doubledescription(V)) |
| 9 | + return [(-H_i.a, -H_i.β) for H_i in H] |
| 10 | +end |
| 11 | + |
| 12 | +function _distance(w̄, b̄, δ_OPS_optimizer) |
| 13 | + y = MOI.get(δ_OPS_optimizer, MOI.ListOfVariableIndices()) |
| 14 | + MOI.set( |
| 15 | + δ_OPS_optimizer, |
| 16 | + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), |
| 17 | + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w̄, y), 0.0), |
| 18 | + ) |
| 19 | + MOI.set(δ_OPS_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) |
| 20 | + MOI.optimize!(δ_OPS_optimizer) |
| 21 | + return b̄ - MOI.get(δ_OPS_optimizer, MOI.ObjectiveValue()) |
| 22 | +end |
| 23 | + |
| 24 | +function _select_next_halfspace(H, δ_OPS_optimizer) |
| 25 | + distances = [_distance(w, b, δ_OPS_optimizer) for (w, b) in H] |
| 26 | + index = argmax(distances) |
| 27 | + w, b = H[index] |
| 28 | + return distances[index], w, b |
| 29 | +end |
| 30 | + |
| 31 | +function MOA.minimize_multiobjective!( |
| 32 | + algorithm::MOA.Sandwiching, |
| 33 | + model::MOA.Optimizer, |
| 34 | + inner::MOI.ModelLike, |
| 35 | + f::MOI.AbstractVectorFunction, |
| 36 | +) |
| 37 | + @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE |
| 38 | + solutions = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}() |
| 39 | + variables = MOI.get(inner, MOI.ListOfVariableIndices()) |
| 40 | + n = MOI.output_dimension(f) |
| 41 | + scalars = MOI.Utilities.scalarize(f) |
| 42 | + status = MOI.OPTIMAL |
| 43 | + δ_OPS_optimizer = MOI.instantiate(model.optimizer_factory) |
| 44 | + if MOI.supports(δ_OPS_optimizer, MOI.Silent()) |
| 45 | + MOI.set(δ_OPS_optimizer, MOI.Silent(), true) |
| 46 | + end |
| 47 | + y = MOI.add_variables(δ_OPS_optimizer, n) |
| 48 | + anchors = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}() |
| 49 | + yI, yUB = zeros(n), zeros(n) |
| 50 | + for (i, f_i) in enumerate(scalars) |
| 51 | + MOI.set(inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i) |
| 52 | + MOA.optimize_inner!(model) |
| 53 | + status = MOI.get(inner, MOI.TerminationStatus()) |
| 54 | + if !MOA._is_scalar_status_optimal(model) |
| 55 | + return status, nothing |
| 56 | + end |
| 57 | + X, Y = MOA._compute_point(model, variables, f) |
| 58 | + model.ideal_point[i] = Y[i] |
| 59 | + yI[i] = Y[i] |
| 60 | + anchors[Y] = X |
| 61 | + MOI.set(inner, MOI.ObjectiveSense(), MOI.MAX_SENSE) |
| 62 | + MOA.optimize_inner!(model) |
| 63 | + status = MOI.get(inner, MOI.TerminationStatus()) |
| 64 | + if !MOA._is_scalar_status_optimal(model) |
| 65 | + MOA._warn_on_nonfinite_anti_ideal(algorithm, MOI.MIN_SENSE, i) |
| 66 | + return status, nothing |
| 67 | + end |
| 68 | + _, Y = MOA._compute_point(model, variables, f_i) |
| 69 | + yUB[i] = Y |
| 70 | + MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE) |
| 71 | + e_i = Float64.(1:n .== i) |
| 72 | + MOI.add_constraint( |
| 73 | + δ_OPS_optimizer, |
| 74 | + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(e_i, y), 0.0), |
| 75 | + MOI.GreaterThan(yI[i]), |
| 76 | + ) |
| 77 | + MOI.add_constraint( |
| 78 | + δ_OPS_optimizer, |
| 79 | + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(e_i, y), 0.0), |
| 80 | + MOI.LessThan(yUB[i]), |
| 81 | + ) |
| 82 | + end |
| 83 | + IPS = [yUB, keys(anchors)...] |
| 84 | + merge!(solutions, anchors) |
| 85 | + u = MOI.add_variables(inner, n) |
| 86 | + u_constraints = [ # u_i >= 0 for all i = 1:n |
| 87 | + MOI.add_constraint(inner, u_i, MOI.GreaterThan{Float64}(0)) for u_i in u |
| 88 | + ] |
| 89 | + f_constraints = [ # f_i + u_i <= yUB_i for all i = 1:n |
| 90 | + MOI.Utilities.normalize_and_add_constraint( |
| 91 | + inner, |
| 92 | + scalars[i] + u[i], |
| 93 | + MOI.LessThan(yUB[i]), |
| 94 | + ) for i in 1:n |
| 95 | + ] |
| 96 | + H = _halfspaces(IPS) |
| 97 | + count = 0 |
| 98 | + while !isempty(H) |
| 99 | + ret = MOA._check_premature_termination(model) |
| 100 | + if ret !== nothing |
| 101 | + status = ret |
| 102 | + break |
| 103 | + end |
| 104 | + count += 1 |
| 105 | + δ, w, b = _select_next_halfspace(H, δ_OPS_optimizer) |
| 106 | + if δ - 1e-3 <= algorithm.precision # added some convergence tolerance |
| 107 | + break |
| 108 | + end |
| 109 | + # would not terminate when precision is set to 0 |
| 110 | + new_f = sum(w[i] * (scalars[i] + u[i]) for i in 1:n) # w' * (f(x) + u) |
| 111 | + MOI.set(inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f) |
| 112 | + MOA.optimize_inner!(model) |
| 113 | + status = MOI.get(inner, MOI.TerminationStatus()) |
| 114 | + if !MOA._is_scalar_status_optimal(model) |
| 115 | + return status, nothing |
| 116 | + end |
| 117 | + β̄ = MOI.get(inner, MOI.ObjectiveValue()) |
| 118 | + X, Y = MOA._compute_point(model, variables, f) |
| 119 | + solutions[Y] = X |
| 120 | + MOI.add_constraint( |
| 121 | + δ_OPS_optimizer, |
| 122 | + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, y), 0.0), |
| 123 | + MOI.GreaterThan(β̄), |
| 124 | + ) |
| 125 | + IPS = push!(IPS, Y) |
| 126 | + H = _halfspaces(IPS) |
| 127 | + end |
| 128 | + MOI.delete.(inner, f_constraints) |
| 129 | + MOI.delete.(inner, u_constraints) |
| 130 | + MOI.delete.(inner, u) |
| 131 | + return status, [MOA.SolutionPoint(X, Y) for (Y, X) in solutions] |
| 132 | +end |
0 commit comments