|
| 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 | + DominguezRios() |
| 8 | +
|
| 9 | +`DominguezRios` implements the algorithm of: |
| 10 | +
|
| 11 | +Dominguez-Rios, M.A. & Chicano, F., & Alba, E. (2021). Effective anytime |
| 12 | +algorithm for multiobjective combinatorial optimization problems. Information |
| 13 | +Sciences, 565(7), 210-228. |
| 14 | +""" |
| 15 | +mutable struct DominguezRios <: AbstractAlgorithm end |
| 16 | + |
| 17 | +mutable struct _DominguezRiosBox |
| 18 | + l::Vector{Float64} |
| 19 | + u::Vector{Float64} |
| 20 | + priority::Float64 |
| 21 | + function _DominguezRiosBox( |
| 22 | + l::Vector{Float64}, |
| 23 | + u::Vector{Float64}, |
| 24 | + p::Float64 = 0.0, |
| 25 | + ) |
| 26 | + @assert length(l) == length(u) "Dimension mismatch between l and u" |
| 27 | + return new(l, u, p) |
| 28 | + end |
| 29 | +end |
| 30 | + |
| 31 | +function _reduced_scaled_priority( |
| 32 | + l::Vector{Float64}, |
| 33 | + u::Vector{Float64}, |
| 34 | + i::Int, |
| 35 | + z::Vector{Float64}, |
| 36 | + yI::Vector{Float64}, |
| 37 | + yN::Vector{Float64}, |
| 38 | +) |
| 39 | + ret = prod((u - l) ./ (yN - yI)) |
| 40 | + if i != length(z) |
| 41 | + return ret |
| 42 | + end |
| 43 | + return ret - prod((z - l) ./ (yN - yI)) |
| 44 | +end |
| 45 | + |
| 46 | +function _p_partition( |
| 47 | + B::_DominguezRiosBox, |
| 48 | + z::Vector{Float64}, |
| 49 | + yI::Vector{Float64}, |
| 50 | + yN::Vector{Float64}, |
| 51 | +) |
| 52 | + ẑ = max.(z, B.l) |
| 53 | + ret = _DominguezRiosBox[] |
| 54 | + for i in 1:length(z) |
| 55 | + new_l = vcat(B.l[1:i], ẑ[i+1:end]) |
| 56 | + new_u = vcat(B.u[1:i-1], ẑ[i], B.u[i+1:end]) |
| 57 | + new_priority = _reduced_scaled_priority(new_l, new_u, i, ẑ, yI, yN) |
| 58 | + push!(ret, _DominguezRiosBox(new_l, new_u, new_priority)) |
| 59 | + end |
| 60 | + return ret |
| 61 | +end |
| 62 | + |
| 63 | +function _select_next_box(L::Vector{Vector{_DominguezRiosBox}}, k::Int) |
| 64 | + p = length(L) |
| 65 | + if any(.!isempty.(L)) |
| 66 | + k = k % p + 1 |
| 67 | + while isempty(L[k]) |
| 68 | + k = k % p + 1 |
| 69 | + end |
| 70 | + i = argmax([B.priority for B in L[k]]) |
| 71 | + end |
| 72 | + return i, k |
| 73 | +end |
| 74 | + |
| 75 | +function _join( |
| 76 | + A::_DominguezRiosBox, |
| 77 | + B::_DominguezRiosBox, |
| 78 | + i::Int, |
| 79 | + z::Vector{Float64}, |
| 80 | + yI::Vector{Float64}, |
| 81 | + yN::Vector{Float64}, |
| 82 | +) |
| 83 | + lᵃ, uᵃ, lᵇ, uᵇ = A.l, A.u, B.l, B.u |
| 84 | + @assert all(uᵃ .<= uᵇ) "`join` operation not valid. (uᵃ ≰ uᵇ)" |
| 85 | + lᶜ, uᶜ = min.(lᵃ, lᵇ), uᵇ |
| 86 | + ẑ = max.(z, lᶜ) |
| 87 | + priority = _reduced_scaled_priority(lᶜ, uᶜ, i, ẑ, yI, yN) |
| 88 | + return _DominguezRiosBox(lᶜ, uᶜ, priority) |
| 89 | +end |
| 90 | + |
| 91 | +function Base.isempty(B::_DominguezRiosBox) |
| 92 | + return any(isapprox(B.l[i], B.u[i]) for i in 1:length(B.u)) |
| 93 | +end |
| 94 | + |
| 95 | +function _update!( |
| 96 | + L::Vector{Vector{_DominguezRiosBox}}, |
| 97 | + z::Vector{Float64}, |
| 98 | + yI::Vector{Float64}, |
| 99 | + yN::Vector{Float64}, |
| 100 | +) |
| 101 | + T = [_DominguezRiosBox[] for _ in 1:length(L)] |
| 102 | + for j in 1:length(L) |
| 103 | + for B in L[j] |
| 104 | + if all(z .< B.u) |
| 105 | + for (i, Bᵢ) in enumerate(_p_partition(B, z, yI, yN)) |
| 106 | + if !isempty(Bᵢ) |
| 107 | + push!(T[i], Bᵢ) |
| 108 | + end |
| 109 | + end |
| 110 | + else |
| 111 | + push!(T[j], B) |
| 112 | + end |
| 113 | + end |
| 114 | + end |
| 115 | + L .= T |
| 116 | + for k in 1:length(L) |
| 117 | + i = 1 |
| 118 | + N = length(L[k]) |
| 119 | + while i < N |
| 120 | + index_to_remove = Int[] |
| 121 | + for j in i:N |
| 122 | + if i != j |
| 123 | + if all(L[k][i].u .<= L[k][j].u) |
| 124 | + L[k][i] = _join(L[k][i], L[k][j], k, z, yI, yN) |
| 125 | + push!(index_to_remove, j) |
| 126 | + elseif all(L[k][i].u .>= L[k][j].u) |
| 127 | + L[k][i] = _join(L[k][j], L[k][i], k, z, yI, yN) |
| 128 | + push!(index_to_remove, j) |
| 129 | + end |
| 130 | + end |
| 131 | + end |
| 132 | + i += 1 |
| 133 | + N -= length(index_to_remove) |
| 134 | + deleteat!(L[k], index_to_remove) |
| 135 | + end |
| 136 | + end |
| 137 | + return |
| 138 | +end |
| 139 | + |
| 140 | +function optimize_multiobjective!(algorithm::DominguezRios, model::Optimizer) |
| 141 | + sense = MOI.get(model.inner, MOI.ObjectiveSense()) |
| 142 | + if sense == MOI.MAX_SENSE |
| 143 | + old_obj, neg_obj = copy(model.f), -model.f |
| 144 | + MOI.set(model, MOI.ObjectiveFunction{typeof(neg_obj)}(), neg_obj) |
| 145 | + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) |
| 146 | + status, solutions = optimize_multiobjective!(algorithm, model) |
| 147 | + MOI.set(model, MOI.ObjectiveFunction{typeof(old_obj)}(), old_obj) |
| 148 | + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) |
| 149 | + if solutions !== nothing |
| 150 | + solutions = [SolutionPoint(s.x, -s.y) for s in solutions] |
| 151 | + end |
| 152 | + return status, solutions |
| 153 | + end |
| 154 | + n = MOI.output_dimension(model.f) |
| 155 | + L = [_DominguezRiosBox[] for i in 1:n] |
| 156 | + scalars = MOI.Utilities.scalarize(model.f) |
| 157 | + variables = MOI.get(model.inner, MOI.ListOfVariableIndices()) |
| 158 | + yI, yN = zeros(n), zeros(n) |
| 159 | + # Ideal and Nadir point estimation |
| 160 | + for (i, f_i) in enumerate(scalars) |
| 161 | + MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i) |
| 162 | + MOI.set(model.inner, MOI.ObjectiveSense(), sense) |
| 163 | + MOI.optimize!(model.inner) |
| 164 | + status = MOI.get(model.inner, MOI.TerminationStatus()) |
| 165 | + if !_is_scalar_status_optimal(status) |
| 166 | + return status, nothing |
| 167 | + end |
| 168 | + _, Y = _compute_point(model, variables, f_i) |
| 169 | + yI[i] = Y |
| 170 | + rev_sense = sense == MOI.MIN_SENSE ? MOI.MAX_SENSE : MOI.MIN_SENSE |
| 171 | + MOI.set(model.inner, MOI.ObjectiveSense(), rev_sense) |
| 172 | + MOI.optimize!(model.inner) |
| 173 | + status = MOI.get(model.inner, MOI.TerminationStatus()) |
| 174 | + if !_is_scalar_status_optimal(status) |
| 175 | + @warn( |
| 176 | + "Unable to solve problem using `DominguezRios()` because " * |
| 177 | + "objective $i does not have a finite domain.", |
| 178 | + ) |
| 179 | + return status, nothing |
| 180 | + end |
| 181 | + _, Y = _compute_point(model, variables, f_i) |
| 182 | + yN[i] = Y |
| 183 | + end |
| 184 | + MOI.set(model.inner, MOI.ObjectiveSense(), sense) |
| 185 | + ϵ = 1 / (2 * n * (maximum(yN - yI) - 1)) |
| 186 | + push!(L[1], _DominguezRiosBox(yI, yN, 0.0)) |
| 187 | + t_max = MOI.add_variable(model.inner) |
| 188 | + solutions = SolutionPoint[] |
| 189 | + k = 0 |
| 190 | + while any(!isempty(l) for l in L) |
| 191 | + i, k = _select_next_box(L, k) |
| 192 | + B = L[k][i] |
| 193 | + w = 1 ./ max.(1, B.u - yI) |
| 194 | + constraints = [ |
| 195 | + MOI.Utilities.normalize_and_add_constraint( |
| 196 | + model.inner, |
| 197 | + t_max - (w[i] * (scalars[i] - yI[i])), |
| 198 | + MOI.GreaterThan(0.0), |
| 199 | + ) for i in 1:n |
| 200 | + ] |
| 201 | + new_f = t_max + ϵ * sum(w[i] * (scalars[i] - yI[i]) for i in 1:n) |
| 202 | + MOI.set(model.inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f) |
| 203 | + MOI.optimize!(model.inner) |
| 204 | + if MOI.get(model.inner, MOI.TerminationStatus()) == MOI.OPTIMAL |
| 205 | + X, Y = _compute_point(model, variables, model.f) |
| 206 | + obj = MOI.get(model.inner, MOI.ObjectiveValue()) |
| 207 | + if (obj < 1) && all(yI .< B.u) |
| 208 | + push!(solutions, SolutionPoint(X, Y)) |
| 209 | + _update!(L, Y, yI, yN) |
| 210 | + else |
| 211 | + deleteat!(L[k], i) |
| 212 | + end |
| 213 | + end |
| 214 | + MOI.delete.(model.inner, constraints) |
| 215 | + end |
| 216 | + MOI.delete(model.inner, t_max) |
| 217 | + return MOI.OPTIMAL, solutions |
| 218 | +end |
0 commit comments