inconsistent outputs with 2 methods for problems with 3 objectives

[xgandibleux]

In trying to set-up a basic 3-obj problem, I observed (if I am not wrong) two cases of inconsistent outputs with the KirlikSayin and the TambyVanderpooten implementations.

The code is given below.

1. result_count(model) returns 0 for KirlikSayin. DominguezRios and TambyVanderpooten return 1 non-dominated point and the corresponding solution is valid.

----- KirlikSayin -----
Statut: OPTIMAL
number of non-dominated points: 0

----- DominguezRios -----
Statut: OPTIMAL
number of non-dominated points: 1

Variables (number of trucks) :
x[1,1] = 5.0
x[1,2] = 0.0
x[1,3] = 5.0
x[2,1] = 0.0
x[2,2] = 8.0
x[2,3] = -0.0
Objectives :
 - Total cost  = 2380.0
 - Total time  = 46.0
 - Total truck = 18.0

----- TambyVanderpooten -----
Statut: OPTIMAL
number of non-dominated points: 1

Variables (number of trucks) :
x[1,1] = 5.0
x[1,2] = 0.0
x[1,3] = 5.0
x[2,1] = 0.0
x[2,2] = 8.0
x[2,3] = -0.0
Objectives :
 - Total cost  = 2380.0
 - Total time  = 46.0
 - Total truck = 18.0

2. I modified 1 value in the matrix time as follow [ 8 4 5; 3 2 4 ] in order to get a problem with more non-dominated points:

• KirlikSayin: result_count(model) returns 0
• DominguezRios: result_count(model) returns 6; all solutions are valid
• TambyVanderpooten: result_count(model) returns 1; the solution is valid (and provided also by DominguezRios)

----- KirlikSayin -----
Statut: OPTIMAL
number of non-dominated points: 0

----- DominguezRios -----
Statut: OPTIMAL
number of non-dominated points: 6

Variables (number of trucks) :
x[1,1] = 3.0
x[1,2] = 2.0
x[1,3] = 5.0
x[2,1] = 2.0
x[2,2] = 6.0
x[2,3] = 0.0
Objectives :
 - Total cost  = 2500.0
 - Total time  = 75.0
 - Total truck = 18.0

Variables (number of trucks) :
x[1,1] = 2.0
x[1,2] = 3.0
x[1,3] = 5.0
x[2,1] = 3.0
x[2,2] = 5.0
x[2,3] = 0.0
Objectives :
 - Total cost  = 2560.0
 - Total time  = 72.0
 - Total truck = 18.0

Variables (number of trucks) :
x[1,1] = 4.0
x[1,2] = 1.0
x[1,3] = 5.0
x[2,1] = 1.0
x[2,2] = 7.0
x[2,3] = 0.0
Objectives :
 - Total cost  = 2440.0
 - Total time  = 78.0
 - Total truck = 18.0

Variables (number of trucks) :
x[1,1] = 1.0
x[1,2] = 4.0
x[1,3] = 5.0
x[2,1] = 4.0
x[2,2] = 4.0
x[2,3] = 0.0
Objectives :
 - Total cost  = 2620.0
 - Total time  = 69.0
 - Total truck = 18.0

Variables (number of trucks) :
x[1,1] = -8.881784197001252e-16
x[1,2] = 5.000000000000001
x[1,3] = 5.0
x[2,1] = 5.000000000000001
x[2,2] = 2.999999999999999
x[2,3] = 0.0
Objectives :
 - Total cost  = 2680.0
 - Total time  = 66.0
 - Total truck = 18.0

Variables (number of trucks) :
x[1,1] = 5.0
x[1,2] = 0.0
x[1,3] = 5.0
x[2,1] = -0.0
x[2,2] = 8.0
x[2,3] = 0.0
Objectives :
 - Total cost  = 2380.0
 - Total time  = 81.0
 - Total truck = 18.0

----- TambyVanderpooten -----
Statut: OPTIMAL
number of non-dominated points: 1

Variables (number of trucks) :
x[1,1] = 5.0
x[1,2] = 0.0
x[1,3] = 5.0
x[2,1] = 0.0
x[2,2] = 8.0
x[2,3] = -0.0
Objectives :
 - Total cost  = 2380.0
 - Total time  = 81.0
 - Total truck = 18.0

Here the JuMP+MOA model and the 3 optimizations:

using JuMP, HiGHS
import MultiObjectiveAlgorithms as MOA

# Data
cost = [ 100 150 200;    # warehouse 1 to customers 1, 2, 3
         120 110 170     # warehouse 2 to customers 1, 2, 3
       ]

time = [ 1 4 5;
         3 2 4
       ]

capacity = [100, 80]     # capacity of warehouse  (in units)
demand   = [50, 80, 50]  # demand of customers (in units)
truck_capacity = 10      # capacity a truck

nWarehouses, nCustomers = size(cost)

model = Model()
set_silent(model)
set_optimizer(model,()->MOA.Optimizer(HiGHS.Optimizer))

# variables: x[i,j] = nb of trucks from i to j
@variable(model, x[1:nWarehouses, 1:nCustomers] >= 0, Int)  

# Constraints : capacity of warehouses
@constraint(model, [i=1:nWarehouses], 10 * sum(x[i, j] for j in 1:nCustomers) <= capacity[i])
# Constraints : demands of customers
@constraint(model, [j=1:nCustomers], 10 * sum(x[i,j] for i in 1:nWarehouses) == demand[j])

# Objectives 
@expression(model, total_cost, sum(cost[i,j] * x[i,j] for i in 1:nWarehouses, j in 1:nCustomers))
@expression(model, total_time, sum(time[i,j] * x[i,j] for i in 1:nWarehouses, j in 1:nCustomers))
@expression(model, total_trucks, sum(x[i,j] for i in 1:nWarehouses, j in 1:nCustomers))
@objective(model, Min, [total_cost, total_time, total_trucks] )


# --------- KirlikSayin ---------
set_attribute(model, MOA.Algorithm(), MOA.KirlikSayin())
optimize!(model)

# Results
println("\n----- KirlikSayin -----")
println("Statut: ", termination_status(model))
println("number of non-dominated points: ", result_count(model))

for s in 1:result_count(model)
    println("\nVariables (number of trucks) :")
    for i in 1:2, j in 1:3
        println("x[$i,$j] = ", value(x[i,j]; result = s))
    end

    println("Objectives :")
    println(" - Total cost  = ", value(total_cost; result = s))
    println(" - Total time  = ", value(total_time; result = s))
    println(" - Total truck = ", value(total_trucks; result = s))
end


# --------- DominguezRios ---------
set_attribute(model, MOA.Algorithm(), MOA.DominguezRios())
optimize!(model)

# Results
println("\n----- DominguezRios -----")
println("Statut: ", termination_status(model))
println("number of non-dominated points: ", result_count(model))

for s in 1:result_count(model)
    println("\nVariables (number of trucks) :")
    for i in 1:2, j in 1:3
        println("x[$i,$j] = ", value(x[i,j]; result = s))
    end

    println("Objectives :")
    println(" - Total cost  = ", value(total_cost; result = s))
    println(" - Total time  = ", value(total_time; result = s))
    println(" - Total truck = ", value(total_trucks; result = s))
end

# --------- TambyVanderpooten ---------
set_attribute(model, MOA.Algorithm(), MOA.TambyVanderpooten())
optimize!(model)

# Results
println("\n----- TambyVanderpooten -----")
println("Statut: ", termination_status(model))
println("number of non-dominated points: ", result_count(model))

for s in 1:result_count(model)
    println("\nVariables (number of trucks) :")
    for i in 1:2, j in 1:3
        println("x[$i,$j] = ", value(x[i,j]; result = s))
    end

    println("Objectives :")
    println(" - Total cost  = ", value(total_cost; result = s))
    println(" - Total time  = ", value(total_time; result = s))
    println(" - Total truck = ", value(total_trucks; result = s))
end

My config:

• HiGHS v1.17.0
• JuMP v1.25.0
• MultiObjectiveAlgorithms v1.4.2

[kofgokhan]

Some time ago, we have increased to value of the nadir point by one in order to not miss any nondominated points with DominguezRios when the objective functions are vector of variables (there was a special edge case). I incremented the ideal point instead of the nadir point. That's part of the problem.

Another issue is that Tamby & Vanderpooten is trying to solve an infeasible problem but that algorithm is not supposed to try solve any infeasible problems. That terminates the algorithm immaturely and it returns nothing. Also when it terminates immaturely it modifies the problem with epsilon constraints so the original model is modified.

I updated the TV for the nadir point computation and KS for the ideal point computation. Thats what I am getting now with the changes:

for alg in [MOA.KirlikSayin(), MOA.TambyVanderpooten(), MOA.DominguezRios()]
    set_optimizer_attribute(model, MOA.Algorithm(), alg)
    optimize!(model)
    @info alg, result_count(model)
end

[ Info: (MultiObjectiveAlgorithms.KirlikSayin(), 6)
[ Info: (MultiObjectiveAlgorithms.TambyVanderpooten(), 6)
[ Info: (MultiObjectiveAlgorithms.DominguezRios(), 6)

[kofgokhan]

Some time ago, we have increased to value of the nadir point by one in order to not miss any nondominated points with DominguezRios when the objective functions are vector of variables (there was a special edge case). I incremented the ideal point instead of the nadir point. That's part of the problem.

Another issue is that Tamby & Vanderpooten is trying to solve an infeasible problem but that algorithm is not supposed to try solve any infeasible problems. That terminates the algorithm immaturely and it returns nothing. Also when it terminates immaturely it modifies the problem with epsilon constraints so the original model is modified.

I updated the TV for the nadir point computation and KS for the ideal point computation. Thats what I am getting now with the changes:

for alg in [MOA.KirlikSayin(), MOA.TambyVanderpooten(), MOA.DominguezRios()]
    set_optimizer_attribute(model, MOA.Algorithm(), alg)
    optimize!(model)
    @info alg, result_count(model)
end

[ Info: (MultiObjectiveAlgorithms.KirlikSayin(), 6)
[ Info: (MultiObjectiveAlgorithms.TambyVanderpooten(), 6)
[ Info: (MultiObjectiveAlgorithms.DominguezRios(), 6)

This did not work for the other tests so I am gonna have to look further.

[odow]

I just fixed a bug in DominguezRios: #104. The algorithms on paper don't seem to take into consideration the tolerances of the solver used for the subproblems.

[odow]

I tweaked the example a bit to make it more convenient for testing. I'll take a look:

julia> using JuMP, HiGHS

julia> import MultiObjectiveAlgorithms as MOA

julia> function main(algorithm)
           cost, time = [100 150 200; 120 110 170], [8 4 5; 3 2 4]
           capacity, demand = [10, 8], [5, 8, 5]
           m, n = size(cost)
           model = Model(() -> MOA.Optimizer(HiGHS.Optimizer))
           set_silent(model)
           @variable(model, x[1:m, 1:n] >= 0, Int)  
           @constraint(model, [i in 1:m], sum(x[i,:]) <= capacity[i])
           @constraint(model, [j in 1:n], sum(x[:,j]) == demand[j])
           @objective(model, Min, [sum(cost .* x), sum(time .* x), sum(x)])
           set_attribute(model, MOA.Algorithm(), algorithm)
           optimize!(model)
           return solution_summary(model)
       end
main (generic function with 1 method)

julia> main(MOA.KirlikSayin())
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.KirlikSayin, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 0
│ ├ raw_status         : Solve complete. Found 0 solution(s)
│ └ objective_bound    : [2.38000e+03,6.60000e+01,1.80000e+01]
├ Solution (result = 1)
│ ├ primal_status        : NO_SOLUTION
│ └ dual_status          : NO_SOLUTION
└ Work counters
  └ solve_time (sec)   : 1.99218e-02

julia> main(MOA.DominguezRios())

solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.DominguezRios, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 6
│ ├ raw_status         : Solve complete. Found 6 solution(s)
│ └ objective_bound    : [2.38000e+03,6.60000e+01,1.80000e+01]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [2.50000e+03,7.50000e+01,1.80000e+01]
└ Work counters
  └ solve_time (sec)   : 7.67150e-02

julia> main(MOA.TambyVanderpooten())
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.TambyVanderpooten, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : Solve complete. Found 1 solution(s)
│ └ objective_bound    : [2.38000e+03,6.60000e+01,1.80000e+01]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [2.38000e+03,8.10000e+01,1.80000e+01]
└ Work counters
  └ solve_time (sec)   : 4.19478e-02

[kofgokhan]

I just fixed a bug in DominguezRios: #104. The algorithms on paper don't seem to take into consideration the tolerances of the solver used for the subproblems.

Can this be a better formulation for the weighted Tchebychev scalarization? I understand the scaling needs to remain but this seems cleaner:

from: https://doi.org/10.1007/s00186-023-00841-0

[kofgokhan]

I have updated selection procedure for TV. I think playing around with the nadir point changed the hypervolume computations for the selection. There is an inconsistency in the algorithm where the algorithm stops when the upper bounds are depleted but there is also a selection rule when there are no upper bounds ( k = 1, u = (M, ..., M) ).

I think we should have more problems for testing. Especially different types of problems other than knapsack, assignment, etc.

julia> using JuMP, HiGHS

julia> import MultiObjectiveAlgorithms as MOA
Precompiling MultiObjectiveAlgorithms finished.
  1 dependency successfully precompiled in 1 seconds. 42 already precompiled.

julia> function main(algorithm)
                  cost, time = [100 150 200; 120 110 170], [8 4 5; 3 2 4]
                  capacity, demand = [10, 8], [5, 8, 5]
                  m, n = size(cost)
                  model = Model(() -> MOA.Optimizer(HiGHS.Optimizer))
                  set_silent(model)
                  @variable(model, x[1:m, 1:n] >= 0, Int)  
                  @constraint(model, [i in 1:m], sum(x[i,:]) <= capacity[i])
                  @constraint(model, [j in 1:n], sum(x[:,j]) == demand[j])
                  @objective(model, Min, [sum(cost .* x), sum(time .* x), sum(x)])
                  set_attribute(model, MOA.Algorithm(), algorithm)
                  optimize!(model)
                  return solution_summary(model)
              end
main (generic function with 1 method)

julia> main(MOA.KirlikSayin())
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.KirlikSayin, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 6
│ ├ raw_status         : Solve complete. Found 6 solution(s)
│ └ objective_bound    : [2.38000e+03,6.60000e+01,1.80000e+01]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [2.38000e+03,8.10000e+01,1.80000e+01]
└ Work counters
  └ solve_time (sec)   : 9.36304e-01

julia> main(MOA.DominguezRios())
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.DominguezRios, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 6
│ ├ raw_status         : Solve complete. Found 6 solution(s)
│ └ objective_bound    : [2.38000e+03,6.60000e+01,1.80000e+01]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [2.50000e+03,7.50000e+01,1.80000e+01]
└ Work counters
  └ solve_time (sec)   : 9.63204e-01

julia> main(MOA.TambyVanderpooten())
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.TambyVanderpooten, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 6
│ ├ raw_status         : Solve complete. Found 6 solution(s)
│ └ objective_bound    : [2.38000e+03,6.60000e+01,1.80000e+01]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [2.62000e+03,6.90000e+01,1.80000e+01]
└ Work counters
  └ solve_time (sec)   : 3.98926e-01

[odow]

Fixed by #108, #109, #110, and #111.

I think we should have more problems for testing

I'm very open to ideas for how we could automate this. Are there libraries of small/trivial instances with solutions? I don't particularly want to have to solve anything large in CI.

[kofgokhan]

Fixed by #108, #109, #110, and #111.

I think we should have more problems for testing

I'm very open to ideas for how we could automate this. Are there libraries of small/trivial instances with solutions? I don't particularly want to have to solve anything large in CI.

https://github.com/vOptSolver/vOptLib had some problems but I do not know if solutions are available.

[odow]

We'd need them to be in a standardized file format, like https://github.com/jump-dev/MathOptFormat

Multi objective is also hard because we want tests like "did you find all the solutions". Otherwise it's easy to make a mistake like this one above where we return only a subset of the efficient solutions.

[kofgokhan]

To simplify tests, maybe we can put different problems into different folders with a text file for each problem and have something like:

struct MOKP
...
end

function MOKP(filename)
...
end

function build_model(mokp::MOKP, optimizer)
...
end

[odow]

I'm going to close this issue. I'll open a new one with some thoughts.