diff --git a/src/algorithms/TambyVanderpooten.jl b/src/algorithms/TambyVanderpooten.jl index 87c2040..bae276b 100644 --- a/src/algorithms/TambyVanderpooten.jl +++ b/src/algorithms/TambyVanderpooten.jl @@ -33,32 +33,42 @@ This algorithm is restricted to problems with: struct TambyVanderpooten <: AbstractAlgorithm end """ - const _TAMBY_VANDERPOOTEN_ATOL = 1e-3 + _round_objective(y::Float64, constant::Float64)::Int -This constant is important. +Return `y - constant`, rounded to the nearest integer. -When we have a solution, we need to decide if the `Y` vector is at the edge of -one of the bounds. The paper assumes exact arithmetic. Solvers have tolerances. +The paper assumes exact arithmetic, and the algorithm relies on exact +comparisons, like whether `y[k] == u[k]` or `y[k] < u[k]`, to decide if a point +is on the edge of a search zone or strictly inside it. Solvers have tolerances, +so the same nondominated point can be returned with different noise, for +example, as `1234.00005` and `1233.99998`. If we compare these values directly, +the algorithm can revisit the same point and cycle. -In theory the algorithm works only for problems with integer-valued objectives. -So we could make this as big as `0.5` or something. We previously used `1e-6`, -but this caused problems when the solver's feasibility tolerance meant that the -objective value was something like `1234.00005`. +The algorithm works only for problems with integer-valued objectives (we use an +ε of `1.0` for the strict inequalities), so we round the non-constant part of +every objective value to the nearest integer and store it as an `Int`. This +ensures that all comparisons are exact. -In future, we could consider how to rewrite the algorithm to make it more robust -to numerical issues like this. +All objective values in the algorithm (`yI`, `yN`, `U_N`, `V`, and the keys of +`solutions`) are stored without the objective constant. We still return the +original floating point objective values to the user. This does not change +the algorithm because the comparisons are elementwise and `_h` uses only +differences between points. It also means that the right-hand side of the +ε-constraints, which have their constant moved to the set, is the stored value. """ -const _TAMBY_VANDERPOOTEN_ATOL = 1e-3 +function _round_objective(y::Float64, constant::Float64) + return round(Int, y - constant) +end # This is Algorithm 1 from the paper. function _update_search_region( - U_N::Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}, - y::Vector{Float64}, - yN::Vector{Float64}, + U_N::Dict{Vector{Int},Vector{Vector{Vector{Int}}}}, + y::Vector{Int}, + yN::Vector{Int}, ) p = length(y) - bounds_to_remove = Vector{Float64}[] - bounds_to_add = Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}() + bounds_to_remove = Vector{Int}[] + bounds_to_add = Dict{Vector{Int},Vector{Vector{Vector{Int}}}}() for (u, u_n) in U_N if _is_less(y, u, 0) # k=0 here because we want to check every element push!(bounds_to_remove, u) @@ -102,19 +112,20 @@ function _is_less(y, u, k) return true end -function _get_child(u::Vector{Float64}, y::Vector{Float64}, k::Int) +function _get_child(u::Vector{Int}, y::Vector{Int}, k::Int) @assert length(u) == length(y) child = copy(u) child[k] = y[k] return child end -# This is function h(u, k) in the paper. -function _h(k::Int, u::Vector{Float64}, yI::Vector{Float64}) +# This is function h(u, k) in the paper. We use Float64 because the product can +# overflow an Int. +function _h(k::Int, u::Vector{Int}, yI::Vector{Int}) h = 1.0 for (i, (u_i, yI_i)) in enumerate(zip(u, yI)) if i != k - h *= u_i - yI_i + h *= Float64(u_i - yI_i) end end return h @@ -122,13 +133,13 @@ end # This is Problem (6) in the paper. The upper bound M is our nadir point yN. function _select_search_zone( - U_N::Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}, - yI::Vector{Float64}, - yN::Vector{Float64}, + U_N::Dict{Vector{Int},Vector{Vector{Vector{Int}}}}, + yI::Vector{Int}, + yN::Vector{Int}, ) k_star, u_star, v_star = 1, yN, -Inf for k in 1:length(yI), u in keys(U_N) - if !isapprox(u[k], yN[k]; atol = _TAMBY_VANDERPOOTEN_ATOL) + if u[k] != yN[k] v = _h(k, u, yI) if v > v_star k_star, u_star, v_star = k, u, v @@ -144,7 +155,7 @@ function minimize_multiobjective!( inner::MOI.ModelLike, f::MOI.AbstractVectorFunction, ) - solutions = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}() + solutions = Dict{Vector{Int},SolutionPoint}() status = _minimize_multiobjective!( algorithm, model, @@ -152,20 +163,21 @@ function minimize_multiobjective!( f, solutions, )::MOI.TerminationStatusCode - return status, SolutionPoint[SolutionPoint(X, Y) for (Y, X) in solutions] + return status, collect(values(solutions)) end # To ensure the type stability of `.inner` and `.f`, we use a function barrier. # # We also pass in `solutions` to simplify the function and to ensure that we can # return at any point from the this inner function with a partial list of primal -# solutions. +# solutions. The keys of `solutions` are the rounded objective values that we +# use in the algorithm, and the values are the original solutions. function _minimize_multiobjective!( algorithm::TambyVanderpooten, model::Optimizer, inner::MOI.ModelLike, f::MOI.AbstractVectorFunction, - solutions::Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}, + solutions::Dict{Vector{Int},SolutionPoint}, ) @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE warm_start_supported = @@ -173,8 +185,9 @@ function _minimize_multiobjective!( variables = MOI.get(inner, MOI.ListOfVariableIndices()) n = MOI.output_dimension(f) # Compute the ideal (yI) and nadir (yN) points - yI, yN = zeros(n), zeros(n) + yI, yN = zeros(Int, n), zeros(Int, n) scalars = MOI.Utilities.scalarize(f) + obj_constants = MOI.constant.(scalars) for (i, f_i) in enumerate(scalars) MOI.set(inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i) MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE) @@ -185,7 +198,8 @@ function _minimize_multiobjective!( end _, y_i = _compute_point(model, variables, f_i) _log_subproblem_solve(model, variables) - model.ideal_point[i] = yI[i] = y_i + model.ideal_point[i] = y_i + yI[i] = _round_objective(y_i, obj_constants[i]) MOI.set(inner, MOI.ObjectiveSense(), MOI.MAX_SENSE) optimize_inner!(model) status = MOI.get(inner, MOI.TerminationStatus()) @@ -197,7 +211,7 @@ function _minimize_multiobjective!( _log_subproblem_solve(model, variables) # This is not really the nadir point, it's an upper bound on the # nadir point because we later add epsilon constraints with a -1. - yN[i] = y_i + 1 + yN[i] = _round_objective(y_i, obj_constants[i]) + 1 end MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE) # Instead of adding and deleting constraints during the algorithm, we add @@ -206,15 +220,14 @@ function _minimize_multiobjective!( MOI.Utilities.normalize_and_add_constraint( inner, f_i, - MOI.LessThan(yN[i]), + MOI.LessThan(yN[i] + obj_constants[i]), ) for (i, f_i) in enumerate(scalars) ] - obj_constants = MOI.constant.(scalars) - U_N = Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}( + U_N = Dict{Vector{Int},Vector{Vector{Vector{Int}}}}( # The nadir point, except for the ideal point in position k yN => [[_get_child(yN, yI, k)] for k in 1:n], ) - V = [Tuple{Vector{Float64},Vector{Float64}}[] for k in 1:n] + V = [Tuple{Vector{Int},Vector{Int}}[] for k in 1:n] status = MOI.OPTIMAL sum_f = MOI.Utilities.operate(+, Float64, scalars...) while !isempty(U_N) @@ -230,11 +243,11 @@ function _minimize_multiobjective!( # computing the nadir point. for i in 1:n u_i = ifelse(i == k, yN[i], u[i]) - set = MOI.LessThan(u_i - 1.0 - obj_constants[i]) + set = MOI.LessThan(Float64(u_i - 1)) MOI.set(inner, MOI.ConstraintSet(), ε_constraints[i], set) end if warm_start_supported && !isempty(solutions) - for (x_i, start_i) in solutions[last(U_N[u][k])] + for (x_i, start_i) in solutions[last(U_N[u][k])].x MOI.set(inner, MOI.VariablePrimalStart(), x_i, start_i) end end @@ -244,8 +257,9 @@ function _minimize_multiobjective!( break end # Now solve problem Π²(k, u) - y_k = MOI.get(inner, MOI.ObjectiveValue())::Float64 - set = MOI.LessThan(y_k - obj_constants[k]) + obj = MOI.get(inner, MOI.ObjectiveValue())::Float64 + y_k = _round_objective(obj, obj_constants[k]) + set = MOI.LessThan(Float64(y_k)) MOI.set(inner, MOI.ConstraintSet(), ε_constraints[k], set) MOI.set(inner, MOI.ObjectiveFunction{typeof(sum_f)}(), sum_f) optimize_inner!(model) @@ -253,14 +267,13 @@ function _minimize_multiobjective!( if !_is_scalar_status_optimal(status) break end - X, Y = _compute_point(model, variables, f) - _log_subproblem_solve(model, Y) + X, Y_float = _compute_point(model, variables, f) + _log_subproblem_solve(model, Y_float) + Y = _round_objective.(Y_float, obj_constants) push!(V[k], (u, Y)) - # We want `if !(Y in U_N[u][k])` but this tests exact equality. We - # want an approximate comparison. - if all(!isapprox(Y; atol = _TAMBY_VANDERPOOTEN_ATOL), U_N[u][k]) + if !(Y in U_N[u][k]) _update_search_region(U_N, Y, yN) - solutions[Y] = X + solutions[Y] = SolutionPoint(X, Y_float) end _clean_search_region(U_N, yI, V, k) end @@ -288,7 +301,7 @@ end function _clean_search_region_inner(u′, U_N, yI, V, k) for (u′_k, yI_k) in zip(u′, yI) - if isapprox(u′_k, yI_k; atol = _TAMBY_VANDERPOOTEN_ATOL) + if u′_k == yI_k return true end end @@ -303,7 +316,7 @@ function _clean_search_region_inner(u′, U_N, yI, V, k) end function _comparison_line_16(u′, u, y_k, k) - if !≈(y_k[k], u′[k]; atol = _TAMBY_VANDERPOOTEN_ATOL) + if y_k[k] != u′[k] return false end for (i, (u′_i, u_i)) in enumerate(zip(u′, u)) diff --git a/test/algorithms/test_TambyVanderpooten.jl b/test/algorithms/test_TambyVanderpooten.jl index 675b63d..bb36a75 100644 --- a/test/algorithms/test_TambyVanderpooten.jl +++ b/test/algorithms/test_TambyVanderpooten.jl @@ -173,6 +173,68 @@ function test_solve_failures() return end +function _solve_noisy_mock(mock, counter::Ref{Int}) + highs = HiGHS.Optimizer() + MOI.set(highs, MOI.Silent(), true) + index_map = MOI.copy_to(highs, mock) + MOI.optimize!(highs) + x = [index_map[xi] for xi in MOI.get(mock, MOI.ListOfVariableIndices())] + # Perturb the solution by a deterministic amount that is within the + # tolerances that a typical MIP solver might return. + noise() = 1e-6 * sin(counter[] += 1) + x_primal = [xi + noise() for xi in MOI.get(highs, MOI.VariablePrimal(), x)] + status = MOI.get(highs, MOI.TerminationStatus()) + MOI.Utilities.mock_optimize!(mock, status, x_primal) + obj = MOI.get(highs, MOI.ObjectiveValue()) + MOI.set(mock, MOI.ObjectiveValue(), obj + noise()) + return +end + +function test_noisy_solutions() + C = [ + [3 13 3 20 9; 4 13 3 4 9; 18 9 3 9 17; 1 16 16 11 3; 18 4 7 13 11], + [1 12 4 6 20; 18 5 13 20 7; 9 3 3 2 15; 11 18 1 7 13; 6 7 3 15 18], + [13 13 13 14 15; 16 3 20 7 13; 12 17 9 8 3; 6 10 4 18 5; 4 9 7 14 11], + ] + n = size(C[1], 1) + solutions = Dict{Bool,Set{Vector{Int}}}() + for noisy in (false, true) + counter = Ref(0) + model = MOA.Optimizer() do + mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}()) + MOI.Utilities.set_mock_optimize!( + mock, + m -> noisy ? _solve_noisy_mock(m, counter) : _solve_mock(m), + ) + return mock + end + MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten()) + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.TimeLimitSec(), 60.0) + x = reshape(MOI.add_variables(model, n * n), n, n) + MOI.add_constraint.(model, x, MOI.ZeroOne()) + for i in 1:n + MOI.add_constraint(model, sum(1.0 .* x[i, :]), MOI.EqualTo(1.0)) + MOI.add_constraint(model, sum(1.0 .* x[:, i]), MOI.EqualTo(1.0)) + end + f_k = [sum(1.0 * C[k] .* x) for k in 1:3] + f = MOI.Utilities.operate(vcat, Float64, f_k...) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL + N = MOI.get(model, MOI.ResultCount()) + Y_float = [MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N] + # The algorithm rounds internally, but we return the solver's values. + @test any(y -> any(!isinteger, y), Y_float) == noisy + Y = [round.(Int, y) for y in Y_float] + @test allunique(Y) + solutions[noisy] = Set(Y) + end + @test solutions[false] == solutions[true] + return +end + end # module TestTambyVanderpooten TestTambyVanderpooten.run_tests()