Simple multi-objective examples

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial presents three worked examples of multi-objective linear programs —a bi-objective LP, a bi-objective assignment problem, and a bi-objective shortest path problem—drawn from the operations research literature. It introduces the core JuMP and MultiObjectiveAlgorithms.jl workflow for enumerating Pareto-optimal solutions.

Learning intentions:

  • Formulate a multi-objective program in JuMP by passing a vector of objective expressions to @objective
  • Select a multi-objective algorithm from MultiObjectiveAlgorithms.jl and understand the trade-offs between the Lexicographic and EpsilonConstraint approaches
  • Iterate over the Pareto-optimal solutions returned by the solver using result_count and the result keyword in objective_value

Required packages

This tutorial uses the following packages:

using JuMP
import HiGHS
import MultiObjectiveAlgorithms as MOA

Bi-objective linear problem

This example is taken from Example 6.3 (from Steuer, 1985), page 154 of Ehrgott, M. (2005). Multicriteria Optimization. Springer, Berlin. The code was adapted from an example in vOptGeneric by @xgandibleux.

model = Model(() -> MOA.Optimizer(HiGHS.Optimizer))
set_attribute(model, MOA.Algorithm(), MOA.Lexicographic())
@variable(model, x1 >= 0)
@variable(model, 0 <= x2 <= 3)
@objective(model, Min, [3x1 + x2, -x1 - 2x2])
@constraint(model, 3x1 - x2 <= 6)
optimize!(model)
----------------------------------------------
        MultiObjectiveAlgorithms.jl
----------------------------------------------
Algorithm: Lexicographic
----------------------------------------------
solve #     Obj. 1       Obj. 2       Time
----------------------------------------------
    1    0.00000e+00  0.00000e+00  1.32040e+00
    2    0.00000e+00  0.00000e+00  1.32082e+00
    3    1.20000e+01 -9.00000e+00  1.40796e+00
    4    1.20000e+01 -9.00000e+00  1.40819e+00
----------------------------------------------
termination_status: OPTIMAL
result_count: 2

Total solve time:          1.51544e+00
Time spent in subproblems: 8.01662e-01 (53%)
Number of subproblems:     6
----------------------------------------------
solution_summary(model)
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.Lexicographic, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 2
│ ├ raw_status         : Solve complete. Found 2 solution(s)
│ └ objective_bound    : [0.00000e+00,-9.00000e+00]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [0.00000e+00,0.00000e+00]
└ Work counters
  └ solve_time (sec)   : 1.51544e+00
for i in 1:result_count(model)
    assert_is_solved_and_feasible(model; result = i)
    print(i, ": z = ", objective_value(model; result = i), " | ")
    println("x = ", value.([x1, x2]; result = i))
end
1: z = [0.0, 0.0] | x = [0.0, -0.0]
2: z = [12.0, -9.0] | x = [3.0, 3.0]

Bi-objective linear assignment problem

This example is taken from Example 9.38 (from Ulungu and Teghem, 1994), page 255 of Ehrgott, M. (2005). Multicriteria Optimization. Springer, Berlin. The code was adapted from an example in vOptGeneric by @xgandibleux.

C1 = [5 1 4 7; 6 2 2 6; 2 8 4 4; 3 5 7 1]
C2 = [3 6 4 2; 1 3 8 3; 5 2 2 3; 4 2 3 5]
n = size(C2, 1)
model = Model(() -> MOA.Optimizer(HiGHS.Optimizer))
set_attribute(model, MOA.Algorithm(), MOA.EpsilonConstraint())
@variable(model, x[1:n, 1:n], Bin)
@objective(model, Min, [sum(C1 .* x), sum(C2 .* x)])
@constraint(model, [i in 1:n], sum(x[i, :]) == 1)
@constraint(model, [j in 1:n], sum(x[:, j]) == 1)
optimize!(model)
----------------------------------------------
        MultiObjectiveAlgorithms.jl
----------------------------------------------
Algorithm: EpsilonConstraint
----------------------------------------------
solve #     Obj. 1       Obj. 2       Time
----------------------------------------------
    1    6.00000e+00  2.40000e+01  7.97205e-01
    2    6.00000e+00  2.40000e+01  7.97817e-01
    3    2.20000e+01  7.00000e+00  8.04817e-01
    4    2.20000e+01  7.00000e+00  8.05339e-01
    6    1.90000e+01  1.00000e+01  8.13406e-01
    8    1.60000e+01  1.10000e+01  8.21765e-01
   10    1.20000e+01  1.30000e+01  8.29917e-01
   12    9.00000e+00  1.70000e+01  8.30952e-01
   14    6.00000e+00  2.40000e+01  8.31771e-01
----------------------------------------------
termination_status: OPTIMAL
result_count: 6

Total solve time:          8.32327e-01
Time spent in subproblems: 4.11470e-02 (5%)
Number of subproblems:     15
----------------------------------------------
solution_summary(model)
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.EpsilonConstraint, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 6
│ ├ raw_status         : Solve complete. Found 6 solution(s)
│ └ objective_bound    : [6.00000e+00,7.00000e+00]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [6.00000e+00,2.40000e+01]
└ Work counters
  └ solve_time (sec)   : 8.32327e-01
for i in 1:result_count(model)
    assert_is_solved_and_feasible(model; result = i)
    print(i, ": z = ", round.(Int, objective_value(model; result = i)), " | ")
    println("x = ", round.(Int, value.(x; result = i)))
end
1: z = [6, 24] | x = [0 1 0 0; 0 0 1 0; 1 0 0 0; 0 0 0 1]
2: z = [9, 17] | x = [0 0 1 0; 0 1 0 0; 1 0 0 0; 0 0 0 1]
3: z = [12, 13] | x = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
4: z = [16, 11] | x = [0 0 0 1; 0 1 0 0; 0 0 1 0; 1 0 0 0]
5: z = [19, 10] | x = [0 0 1 0; 1 0 0 0; 0 0 0 1; 0 1 0 0]
6: z = [22, 7] | x = [0 0 0 1; 1 0 0 0; 0 0 1 0; 0 1 0 0]

Bi-objective shortest path problem

This example is taken from Exercise 9.5 page 269 of Ehrgott, M. (2005). Multicriteria Optimization. Springer, Berlin. The code was adapted from an example in vOptGeneric by @xgandibleux.

M = 50
C1 = [
    M 4 5 M M M
    M M 2 1 2 7
    M M M 5 2 M
    M M 5 M M 3
    M M M M M 4
    M M M M M M
]
C2 = [
    M 3 1 M M M
    M M 1 4 2 2
    M M M 1 7 M
    M M 1 M M 2
    M M M M M 2
    M M M M M M
]
n = size(C2, 1)
model = Model(() -> MOA.Optimizer(HiGHS.Optimizer))
set_attribute(model, MOA.Algorithm(), MOA.EpsilonConstraint())
@variable(model, x[1:n, 1:n], Bin)
@objective(model, Min, [sum(C1 .* x), sum(C2 .* x)])
@constraint(model, sum(x[1, :]) == 1)
@constraint(model, sum(x[:, n]) == 1)
@constraint(model, [i in 2:(n-1)], sum(x[i, :]) - sum(x[:, i]) == 0)
optimize!(model)
----------------------------------------------
        MultiObjectiveAlgorithms.jl
----------------------------------------------
Algorithm: EpsilonConstraint
----------------------------------------------
solve #     Obj. 1       Obj. 2       Time
----------------------------------------------
    1    8.00000e+00  9.00000e+00  7.17807e-03
    2    8.00000e+00  9.00000e+00  7.75599e-03
    3    1.30000e+01  4.00000e+00  1.43430e-02
    4    1.30000e+01  4.00000e+00  1.49052e-02
    6    1.10000e+01  5.00000e+00  1.60201e-02
    8    1.00000e+01  7.00000e+00  1.69952e-02
   10    8.00000e+00  9.00000e+00  1.79200e-02
----------------------------------------------
termination_status: OPTIMAL
result_count: 4

Total solve time:          1.85461e-02
Time spent in subproblems: 1.71220e-02 (92%)
Number of subproblems:     11
----------------------------------------------
solution_summary(model)
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.EpsilonConstraint, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 4
│ ├ raw_status         : Solve complete. Found 4 solution(s)
│ └ objective_bound    : [8.00000e+00,4.00000e+00]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [8.00000e+00,9.00000e+00]
└ Work counters
  └ solve_time (sec)   : 1.85461e-02
for i in 1:result_count(model)
    assert_is_solved_and_feasible(model; result = i)
    print(i, ": z = ", round.(Int, objective_value(model; result = i)), " | ")
    X = round.(Int, value.(x; result = i))
    print("Path:")
    for ind in findall(val -> val ≈ 1, X)
        i, j = ind.I
        print(" $i->$j")
    end
    println()
end
1: z = [8, 9] | Path: 1->2 2->4 4->6
2: z = [10, 7] | Path: 1->2 2->5 5->6
3: z = [11, 5] | Path: 1->2 2->6
4: z = [13, 4] | Path: 1->3 3->4 4->6