|
| 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 | +""" |
| 7 | + TambyVanderpooten() |
| 8 | +
|
| 9 | +`TambyVanderpooten` implements the algorithm of: |
| 10 | +
|
| 11 | +Satya Tamby, Daniel Vanderpooten (2021) Enumeration of the Nondominated Set |
| 12 | +of Multiobjective Discrete Optimization Problems. INFORMS Journal on |
| 13 | +Computing 33(1):72-85. |
| 14 | +
|
| 15 | +This is an algorithm to generate all nondominated solutions for multi-objective |
| 16 | +discrete optimization problems. The algorithm maintains upper bounds (for |
| 17 | +minimization problems) and their associated defining points. At each iteration, |
| 18 | +one of the objectives and an upper bound is picked and the single objective |
| 19 | +reformulation is solved using one of the defining points as a starting solution. |
| 20 | +
|
| 21 | +## Supported optimizer attributes |
| 22 | +
|
| 23 | + * `MOI.TimeLimitSec()`: terminate if the time limit is exceeded and return the |
| 24 | + list of current solutions. |
| 25 | +""" |
| 26 | +mutable struct TambyVanderpooten <: AbstractAlgorithm end |
| 27 | + |
| 28 | +function _project(x::Vector{Float64}, axis::Int) |
| 29 | + return [x[i] for i in 1:length(x) if i != axis] |
| 30 | +end |
| 31 | + |
| 32 | +_volume(r::_Rectangle, l::Vector{Float64}) = prod(r.u - l) |
| 33 | + |
| 34 | +function _update_search_region( |
| 35 | + U_N::Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}, |
| 36 | + y::Vector{Float64}, |
| 37 | + yN::Vector{Float64}, |
| 38 | +) |
| 39 | + bounds_to_remove = Vector{Float64}[] |
| 40 | + p = length(y) |
| 41 | + for u in keys(U_N) |
| 42 | + if all(y .< u) |
| 43 | + push!(bounds_to_remove, u) |
| 44 | + for l in 1:p |
| 45 | + u_l = _get_child(u, y, l) |
| 46 | + N = [ |
| 47 | + k != l ? [yi for yi in U_N[u][k] if yi[l] < y[l]] : [y] |
| 48 | + for k in 1:p |
| 49 | + ] |
| 50 | + if all(!isempty(N[k]) for k in 1:p if u_l[k] ≠ yN[k]) |
| 51 | + U_N[u_l] = N |
| 52 | + end |
| 53 | + end |
| 54 | + else |
| 55 | + for k in 1:p |
| 56 | + if (y[k] == u[k]) && all(_project(y, k) .< _project(u, k)) |
| 57 | + push!(U_N[u][k], y) |
| 58 | + end |
| 59 | + end |
| 60 | + end |
| 61 | + end |
| 62 | + for bound_to_remove in bounds_to_remove |
| 63 | + delete!(U_N, bound_to_remove) |
| 64 | + end |
| 65 | + return |
| 66 | +end |
| 67 | + |
| 68 | +function _get_child(u::Vector{Float64}, y::Vector{Float64}, k::Int) |
| 69 | + @assert length(u) == length(y) |
| 70 | + return vcat(u[1:k-1], y[k], u[k+1:length(y)]) |
| 71 | +end |
| 72 | + |
| 73 | +function _select_search_zone( |
| 74 | + U_N::Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}, |
| 75 | + yI::Vector{Float64}, |
| 76 | +) |
| 77 | + i, j = |
| 78 | + argmax([ |
| 79 | + prod(_project(u, k) - _project(yI, k)) for k in 1:length(yI), |
| 80 | + u in keys(U_N) |
| 81 | + ]).I |
| 82 | + return i, collect(keys(U_N))[j] |
| 83 | +end |
| 84 | + |
| 85 | +function optimize_multiobjective!( |
| 86 | + algorithm::TambyVanderpooten, |
| 87 | + model::Optimizer, |
| 88 | +) |
| 89 | + start_time = time() |
| 90 | + sense = MOI.get(model.inner, MOI.ObjectiveSense()) |
| 91 | + if sense == MOI.MAX_SENSE |
| 92 | + old_obj, neg_obj = copy(model.f), -model.f |
| 93 | + MOI.set(model, MOI.ObjectiveFunction{typeof(neg_obj)}(), neg_obj) |
| 94 | + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) |
| 95 | + status, solutions = optimize_multiobjective!(algorithm, model) |
| 96 | + MOI.set(model, MOI.ObjectiveFunction{typeof(old_obj)}(), old_obj) |
| 97 | + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) |
| 98 | + if solutions !== nothing |
| 99 | + solutions = [SolutionPoint(s.x, -s.y) for s in solutions] |
| 100 | + end |
| 101 | + return status, solutions |
| 102 | + end |
| 103 | + warm_start_supported = false |
| 104 | + if MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex) |
| 105 | + warm_start_supported = true |
| 106 | + end |
| 107 | + solutions = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}() |
| 108 | + YN = Vector{Float64}[] |
| 109 | + variables = MOI.get(model.inner, MOI.ListOfVariableIndices()) |
| 110 | + n = MOI.output_dimension(model.f) |
| 111 | + yI, yN = zeros(n), zeros(n) |
| 112 | + scalars = MOI.Utilities.scalarize(model.f) |
| 113 | + for (i, f_i) in enumerate(scalars) |
| 114 | + MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i) |
| 115 | + MOI.set(model.inner, MOI.ObjectiveSense(), sense) |
| 116 | + MOI.optimize!(model.inner) |
| 117 | + status = MOI.get(model.inner, MOI.TerminationStatus()) |
| 118 | + if !_is_scalar_status_optimal(status) |
| 119 | + return status, nothing |
| 120 | + end |
| 121 | + _, Y = _compute_point(model, variables, f_i) |
| 122 | + yI[i] = Y |
| 123 | + MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MAX_SENSE) |
| 124 | + MOI.optimize!(model.inner) |
| 125 | + status = MOI.get(model.inner, MOI.TerminationStatus()) |
| 126 | + if !_is_scalar_status_optimal(status) |
| 127 | + _warn_on_nonfinite_anti_ideal(algorithm, sense, i) |
| 128 | + return status, nothing |
| 129 | + end |
| 130 | + _, Y = _compute_point(model, variables, f_i) |
| 131 | + yN[i] = Y |
| 132 | + end |
| 133 | + MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE) |
| 134 | + U_N = Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}() |
| 135 | + V = [Tuple{Vector{Float64},Vector{Float64}}[] for k in 1:n] |
| 136 | + U_N[yN] = [[_get_child(yN, yI, k)] for k in 1:n] |
| 137 | + status = MOI.OPTIMAL |
| 138 | + while !isempty(U_N) |
| 139 | + if _time_limit_exceeded(model, start_time) |
| 140 | + status = MOI.TIME_LIMIT |
| 141 | + break |
| 142 | + end |
| 143 | + k, u = _select_search_zone(U_N, yI) |
| 144 | + MOI.set( |
| 145 | + model.inner, |
| 146 | + MOI.ObjectiveFunction{typeof(scalars[k])}(), |
| 147 | + scalars[k], |
| 148 | + ) |
| 149 | + ε_constraints = Any[] |
| 150 | + for (i, f_i) in enumerate(scalars) |
| 151 | + if i != k |
| 152 | + ci = MOI.add_constraint( |
| 153 | + model.inner, |
| 154 | + f_i, |
| 155 | + MOI.LessThan{Float64}(u[i] - 1), |
| 156 | + ) |
| 157 | + push!(ε_constraints, ci) |
| 158 | + end |
| 159 | + end |
| 160 | + if u[k] ≠ yN[k] |
| 161 | + if warm_start_supported |
| 162 | + variables_start = solutions[first(U_N[u][k])] |
| 163 | + for x_i in variables |
| 164 | + MOI.set( |
| 165 | + model.inner, |
| 166 | + MOI.VariablePrimalStart(), |
| 167 | + x_i, |
| 168 | + variables_start[x_i], |
| 169 | + ) |
| 170 | + end |
| 171 | + end |
| 172 | + end |
| 173 | + MOI.optimize!(model.inner) |
| 174 | + if !_is_scalar_status_optimal(model) |
| 175 | + return status, nothing |
| 176 | + end |
| 177 | + y_k = MOI.get(model.inner, MOI.ObjectiveValue()) |
| 178 | + sum_f = sum(1.0 * s for s in scalars) |
| 179 | + MOI.set(model.inner, MOI.ObjectiveFunction{typeof(sum_f)}(), sum_f) |
| 180 | + y_k_constraint = |
| 181 | + MOI.add_constraint(model.inner, scalars[k], MOI.EqualTo(y_k)) |
| 182 | + MOI.optimize!(model.inner) |
| 183 | + if !_is_scalar_status_optimal(model) |
| 184 | + return status, nothing |
| 185 | + end |
| 186 | + X, Y = _compute_point(model, variables, model.f) |
| 187 | + MOI.delete.(model, ε_constraints) |
| 188 | + MOI.delete(model, y_k_constraint) |
| 189 | + push!(V[k], (u, Y)) |
| 190 | + if Y ∉ U_N[u][k] |
| 191 | + _update_search_region(U_N, Y, yN) |
| 192 | + solutions[Y] = X |
| 193 | + end |
| 194 | + bounds_to_remove = Vector{Float64}[] |
| 195 | + for u_i in keys(U_N) |
| 196 | + for k in 1:n |
| 197 | + if u_i[k] == yI[k] |
| 198 | + push!(bounds_to_remove, u_i) |
| 199 | + else |
| 200 | + for (u_j, y_j) in V[k] |
| 201 | + if all(_project(u_i, k) .<= _project(u_j, k)) && |
| 202 | + (y_j[k] == u_i[k]) |
| 203 | + push!(bounds_to_remove, u_i) |
| 204 | + end |
| 205 | + end |
| 206 | + end |
| 207 | + end |
| 208 | + end |
| 209 | + if !isempty(bounds_to_remove) |
| 210 | + for bound_to_remove in bounds_to_remove |
| 211 | + delete!(U_N, bound_to_remove) |
| 212 | + end |
| 213 | + end |
| 214 | + end |
| 215 | + solutions = [SolutionPoint(X, Y) for (Y, X) in solutions] |
| 216 | + return status, solutions |
| 217 | +end |
0 commit comments