[GeneralDichotomy] allow any primal feasible point from subsolve

ext/Polyhedra/GeneralDichotomy.jl

@@ -24,37 +24,30 @@ function MOA.minimize_multiobjective!(
     # Storage we need for the algorithm.
     weights, solutions = Weight[], MOA.SolutionPoint[]
     n_obj = MOI.output_dimension(model.f)
+    existing_sol = Dict{Vector{Int},Int}()
     # First, search for an initial primal feasible point.
-    init_sol_idx = 0
-    status, solution = nothing, nothing
     for i in 1:n_obj
         w = zeros(Float64, n_obj)
         w[i] = 1.0
         status, solution = MOA._solve_weighted_sum(model, alg, w)
-        if solution !== nothing
-            init_sol_idx = i
-            push!(solutions, solution)
-            break
+        if solution === nothing
+            # One of the subproblems failed to solve. This means something went
+            # really wrong. A common reason is that the problem is unbounded.
+            return status, nothing
         end
+        push!(solutions, solution)
+        model.ideal_point[i] = solution.y[i]
     end
-    if length(solutions) == 0
-        return status, nothing
-    end
-    # Initialize the weights. There is one weight vector for each objective, and
-    # the weight is set to 1.0 for each objective. We use the current solution
-    # obtained by minimizing the 1st objective as the reference.
+    solution = solutions[1]
+    existing_sol[_round(solution.y; atol)] = 1
     for i in 1:n_obj
         w = zeros(Float64, n_obj)
         w[i] = 1.0
         z = w' * solution.y
         adj_bnd = Int[-j for j in 1:n_obj if j != i]
-        tested = i <= init_sol_idx
-        removed = i < init_sol_idx
-        push!(weights, Weight(w, z, adj_bnd, [1], tested, removed))
+        push!(weights, Weight(w, z, adj_bnd, [1], i == 1, false))
     end
-    # Prevent solution duplicates: existing_sol maps an rounded objective vector
-    # to its index in `solutions::Vector{MOA.SolutionPoint}`.
-    existing_sol = Dict(_round(solution.y; atol) => 1)
+    status = MOI.OPTIMAL
     n_removed = 0
     while length(solutions) < MOI.get(alg, MOA.SolutionLimit())
         if (ret = MOA._check_premature_termination(model)) !== nothing

src/algorithms/GeneralDichotomy.jl

@@ -59,13 +59,20 @@ function _solve_weighted_sum(
     f = _scalarise(model.f, weight)
     MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f)}(), f)
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
-    if !_is_scalar_status_optimal(status)
-        _log_subproblem_solve(model, "subproblem not optimal")
-        return status, nothing
+    term_status = MOI.get(model.inner, MOI.TerminationStatus())
+    primal_status = MOI.get(model.inner, MOI.PrimalStatus())
+    if !_is_scalar_status_feasible_point(primal_status)
+        _log_subproblem_solve(model, "subproblem failed to solve")
+        return term_status, nothing
+    elseif term_status == MOI.DUAL_INFEASIBLE
+        _log_subproblem_solve(model, "subproblem is unbounded")
+        return term_status, nothing
     end
     variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
     X, Y = _compute_point(model, variables, model.f)
+    if !_is_scalar_status_optimal(term_status)
+        _log_subproblem_solve(model, "subproblem not optimal")
+    end
     _log_subproblem_solve(model, Y)
-    return status, SolutionPoint(X, Y)
+    return term_status, SolutionPoint(X, Y)
 end

test/algorithms/test_GeneralDichotomy.jl

@@ -456,7 +456,7 @@ function test_quadratic()
     MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
     MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
     MOI.optimize!(model)
-    @test MOI.get(model, MOI.ResultCount()) == 10
+    @test 9 <= MOI.get(model, MOI.ResultCount()) <= 10
     for i in 1:MOI.get(model, MOI.ResultCount())
         w_sol = MOI.get(model, MOI.VariablePrimal(i), w)
         y = MOI.get(model, MOI.ObjectiveValue(i))
@@ -507,7 +507,7 @@ function test_solve_failures()
         end
         MOI.optimize!(model)
         @test MOI.get(model, MOI.TerminationStatus()) == MOI.NUMERICAL_ERROR
-        fails = [0, 1, 2, 3]
+        fails = [0, 0, 2, 2]
         @test MOI.get(model, MOI.ResultCount()) == fails[fail_after+1]
     end
     return
@@ -537,6 +537,34 @@ function test_dichotomy_issue_191()
     return
 end
 
+function test_general_dichotomy_highs_solution_limit()
+    p1 = [77, 94, 71, 63, 96, 82, 85, 75, 72, 91, 99, 63, 84, 87, 79, 94, 90]
+    p2 = [65, 90, 90, 77, 95, 84, 70, 94, 66, 92, 74, 97, 60, 60, 65, 97, 93]
+    w = [80, 87, 68, 72, 66, 77, 99, 85, 70, 93, 98, 72, 100, 89, 67, 86, 91]
+    model = MOA.Optimizer(HiGHS.Optimizer)
+    MOI.set(model, MOI.RawOptimizerAttribute("mip_max_improving_sols"), 1)
+    MOI.set(model, MOA.Algorithm(), MOA.GeneralDichotomy())
+    # MOI.set(model, MOI.Silent(), true)
+    x = MOI.add_variables(model, length(w))
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    f = MOI.Utilities.operate(
+        vcat,
+        Float64,
+        [sum(1.0 * p[i] * x[i] for i in 1:length(w)) for p in [p1, p2]]...,
+    )
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.add_constraint(
+        model,
+        sum(1.0 * w[i] * x[i] for i in 1:length(w)),
+        MOI.LessThan(900.0),
+    )
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
+    @test MOI.get(model, MOI.ResultCount()) == 3
+    return
+end
+
 end  # module TestGeneralDichotomy
 
 TestGeneralDichotomy.run_tests()