Efficient multi-objective LP implementation

Specific Domains Optimization (Mathematical)
1

Hi, I want to solve the following problem (maximum flow with losses):

For a graph with source node s and target node t let \forall v \notin \{s,t\}: \sum_{u} f_{u, v} = \sum_u \alpha_{v,u} f_{v, u} where \alpha_{v,u} = \alpha_{u,v} > 1 represents the loss coefficient (such that there is always a clear path of least resistance).

Also 0 \leq f_{u,v} \leq c_{u,v} (capacity constraint) holds.

Define F := \max_{f_{u,v}} \sum_u f_{u,t} (maximum output flow) and find the unique (f_{u,v})_{u,v\in V} which attains \min_{f_{u,v}} \sum_v f(s, v) s.t. \sum_u f_{u,t} = F (minimum input flow aka least waste of energy by transportation).

For that, I use MOA.Hierarchical() from MultiObjectiveAlgorithms.jl with HiGHS.jl:

function create_solver(model, resistance_matrix, source, target)
    solver = Model(() -> MOA.Optimizer(HiGHS.Optimizer))
    set_silent(solver)
    n = size(resistance_matrix, 1)
    magic_number = 1000000
    @variable(solver, f[i = 1:n, j = 1:n; near(model, i, j)] >= 0)
    # Capacity constraints
    @constraint(solver, upper[i = 1:n, j = 1:n; near(model, i, j)], f[i, j] <= magic_number)
    # Flow conservation constraints with losses
    @constraint(solver, [i = 1:n; i != source && i != target], sum(f[:, i]) == sum(resistance_matrix[j, i] * f[i, j] for j in 1:n if near(model, i, j)))
    # Multi-objective optimization
    @objective(solver, Max, [sum(f[:, 2]), -sum(f[1, :])])
    set_attribute(solver, MOA.Algorithm(), MOA.Hierarchical())
    set_attribute.(solver, MOA.ObjectivePriority.(1:2), [2, 1])
    return solver
end

However, this is very slow (~20 times slower than the single objective LP). Any ideas on how to speed it up or rewrite it in a smarter way?

2

Do you have a reproducible example? If you have two objectives, it should be about 2x slower than solving a single LP.

1 Like
3

To be frank, I can’t reproduce it anymore. It is indeed around 2x slower.

Nevertheless, I wonder about speedup opportunities. This solver is executed around 8760 times for each function call of a joint cost function in (independent snapshots). Since I call costs multiple times, I create 8760 solvers and execute them with multiple threads. For each function evaluation, I need to update the constraint upper. Is there a smarter approach?

On a related note: Is DiffOpt.jl compatible with MOA? @gdalle

4

I don’t think so, but the right person to ask is @mbesancon. The optimality conditions for lexicographic programming are a bit more involved, and DiffOpt.jl probably cannot deal with them. Ironically, part of my PhD thesis was about trying to differentiate through lexicographic LPs, but it led nowhere (and the maths of that chapter are wrong).

1 Like
5

DiffOpt does not work with MOA.

You’ll need to iterate through setting scalar objectives to get the gradient for each objective.

1 Like
6

Fyi this is my attempt at differentiating this problem manually @gdalle @odow @mbesancon

function create_solver(model, distance_matrix, source, target)
    solver1 = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
    solver2 = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
    set_silent.([solver1, solver2])
    n = size(distance_matrix, 1)
    magic_number = 1000000

    function set_common_constraints!(solver)
        # Variable definition
        @variable(solver, f[i = 1:n, j = 1:n; near(model, i, j)] ≥ 0)
        # Capacity constraints are modified later
        @constraint(solver, upper[i = 1:n, j = 1:n; near(model, i, j)], f[i, j] ≤ magic_number)
        # Flow conservation constraints with losses
        @constraint(solver, [i = 1:n; i ≠ source && i ≠ target], sum(f[:, i]) == sum(distance_matrix[j, i] * f[i, j] for j in 1:n if near(model, i, j)))
    end

    set_common_constraints!.([solver1, solver2])

    # Maximize used energy
    @objective(solver1, Max, sum(solver1[:f][:, 2]))
    
    # Minimize wasted energy
    @constraint(solver2, result, sum(solver2[:f][:, 2]) ≥ magic_number)
    @objective(solver2, Min, sum(solver2[:f][1, :]))

    return (solver1, solver2)
end

The returned tuple of solvers delivers the same result as MOA via

JuMP.optimize!(solver[1])
set_normalized_rhs(solver[2][:result], JuMP.objective_value(solver[1]))
JuMP.optimize!(solver[2])
value.(solver[2][:f]) # result

Differentiation is done by:

dflow = 0.1 # in reality this is a SparseAxisArray
MOI.set.(solver[2], DiffOpt.ReverseVariablePrimal(), solver[2][:f], dflow)
DiffOpt.reverse_differentiate!(solver[2])
obj_exp = MOI.get.(solver[2], DiffOpt.ReverseConstraintFunction(), solver[2][:result])
grad_F = JuMP.constant.(obj_exp)

# Backpropagate the flow matrix derivatives
obj_exp = MOI.get.(solver[2], DiffOpt.ReverseConstraintFunction(), solver[2][:upper])
grad_solver2 = JuMP.constant.(obj_exp)

# Backpropagate the achieved max flow
MOI.set.(solver[1], DiffOpt.ReverseVariablePrimal(), solver[1][:f][:, 2], grad_F)
DiffOpt.reverse_differentiate!(solver[1])
obj_exp = MOI.get.(solver[1], DiffOpt.ReverseConstraintFunction(), solver[1][:upper])
grad_solver1 = JuMP.constant.(obj_exp)
return (grad_solver1, grad_solver2) # the gradient is the sum of these

I am not sure at all that this does what it is supposed to do though…