[Dichotomy] fix partitioning when weight is outside domain (#192)

Author: odow

src/algorithms/Dichotomy.jl

@@ -117,6 +117,7 @@ function optimize_multiobjective!(
     end
     limit = MOI.get(algorithm, SolutionLimit())
     status = MOI.OPTIMAL
+    atol = 1e-6
     while length(queue) > 0 && length(solutions) < limit
         if (ret = _check_premature_termination(model)) !== nothing
             status = ret
@@ -125,15 +126,24 @@ function optimize_multiobjective!(
         (a, b) = popfirst!(queue)
         y_d = solutions[a].y .- solutions[b].y
         w = y_d[2] / (y_d[2] - y_d[1])
+        if !(a <= w <= b)
+            # The weight is outside the domain. This can happen if the solver
+            # finds a solution "outside" the domain by some small numerical
+            # value, either because of constraint feasibility, or, more likely,
+            # because of the MIP gap. See MultiObjectiveAlgorithms.jl#191.
+            # We don't need to search this weight
+            continue
+        end
         status, solution =
             _solve_weighted_sum(algorithm, model, inner, f, [w, 1.0 - w])
         if !_is_scalar_status_optimal(status)
             break # Exit the solve with some error.
-        elseif solution ≈ solutions[a] || solution ≈ solutions[b]
+        end
+        if ≈(solution, solutions[a]; atol) || ≈(solution, solutions[b]; atol)
             # We have found an existing solution. We're free to prune (a, b)
             # from the search space.
         else
-            # Solution is identical to a and b, so search the domain (a, w) and
+            # Solution is different to a and b, so search the domain (a, w) and
             # (w, b), and add solution as a new Pareto-optimal solution!
             push!(queue, (a, w))
             push!(queue, (w, b))

test/algorithms/test_Dichotomy.jl

@@ -461,6 +461,30 @@ function test_scalar_time_limit()
     return
 end
 
+function test_dichotomy_issue_191()
+    #!format:off
+    C = Float64[
+        31 42 34 29 28 24 93 81 60 65 70 26 73 78 70 6 53 17 75 53 37 5 68 84 18 78 14 14 33 13 77 60 34 12 87 47 48 47 22 38 22 46 16 12 40 82 42 51 51 10 100 40 81 86 29 95 97 72 89 87 77 86 40 80 2 64 66 76 24 36 82 68 9 75 94 87 35 36 57 34 10 12 57 46 42 53 22 55 8 38 53 73 77 91 61 98 52 75 47 39 10 30 65 68 65 24 57 70 57 11 77 49 76 14 56 81 77 90 97 70 84 22 41 5 3 50 98 20 23 94 46 39 25 44 63 91 64 8 75 15 11 14 64 72 56 80 79 61 75 83 14 21 92 58 77 92 74 82 96 54 19 15 23 12 89 44 37 29 15 41 61 73 16 86 81 47 74 83 4 20 67 87 9 52 87 21 4 52 16 87 86 35 21 3 97 65 73 8 65 75
+        38 45 53 27 44 71 84 35 66 100 69 95 13 24 90 5 95 28 13 76 65 68 16 42 41 58 31 53 74 91 24 83 43 46 44 26 22 4 47 41 49 56 20 44 23 96 35 94 84 1 20 84 18 99 52 7 64 27 81 70 87 1 72 45 96 35 31 25 42 80 90 70 3 6 96 76 80 67 78 99 66 68 10 91 57 81 65 67 18 49 13 17 7 64 33 16 94 52 55 28 58 74 85 43 8 77 33 57 64 67 66 15 65 8 14 38 12 87 68 8 67 12 2 37 94 32 53 99 63 8 90 25 38 27 68 47 79 43 40 75 4 25 10 58 6 79 18 78 71 11 45 96 77 69 20 57 58 81 6 63 17 82 97 63 21 10 39 71 83 57 50 59 83 47 31 81 33 42 12 46 30 60 41 54 78 92 24 96 74 46 19 92 68 91 39 100 38 2 41 59
+    ]
+    w = Float64[2, 9, 32, 24, 39, 8, 57, 44, 98, 35, 54, 9, 37, 60, 52, 40, 94, 92, 8, 52, 4, 91, 10, 57, 75, 20, 32, 19, 53, 59, 96, 84, 32, 86, 13, 84, 47, 70, 37, 33, 40, 40, 63, 59, 27, 51, 51, 67, 42, 50, 40, 73, 7, 74, 94, 30, 40, 29, 85, 37, 62, 51, 76, 32, 72, 53, 11, 61, 17, 84, 44, 76, 41, 70, 71, 53, 7, 65, 12, 84, 93, 47, 57, 83, 97, 23, 79, 94, 55, 84, 94, 28, 54, 29, 26, 82, 70, 12, 34, 87, 96, 77, 55, 57, 70, 65, 54, 8, 47, 77, 94, 59, 82, 31, 80, 69, 13, 56, 14, 75, 52, 89, 11, 28, 10, 82, 17, 41, 67, 51, 97, 75, 29, 54, 98, 34, 96, 79, 78, 42, 98, 26, 89, 22, 40, 63, 73, 23, 72, 23, 35, 26, 24, 68, 93, 89, 84, 39, 3, 84, 33, 39, 12, 85, 11, 5, 91, 26, 10, 71, 25, 30, 66, 36, 67, 94, 74, 9, 81, 34, 41, 86, 7, 3, 75, 94, 84, 95, 34, 47, 26, 90, 17, 43, 88, 94, 35, 10, 34, 13.0]
+    #!format:on
+    model = MOA.Optimizer(HiGHS.Optimizer)
+    MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOA.Algorithm(), MOA.Dichotomy())
+    x = MOI.add_variables(model, 200)
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    MOI.add_constraint(model, w' * x, MOI.LessThan(5207.5))
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    f = MOI.Utilities.vectorize(C * x)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
+    @test 38 <= MOI.get(model, MOI.ResultCount()) <= 39
+    @test 70 <= MOI.get(model, MOA.SubproblemCount()) <= 80
+    return
+end
+
 end  # module TestDichotomy
 
 TestDichotomy.run_tests()