# Efficient multi-objective LP implementation
**URL:** https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615
**Category:** Optimization (Mathematical)
**Created:** [May 23, 2024, 8:49am UTC](https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615 "2024-05-23T08:49:30Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@PaterPen](https://discourse.julialang.org/u/PaterPen)
#### Post date: [May 23, 2024, 8:49am UTC](https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615/1 "2024-05-23T08:49:31Z")
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:
```julia
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?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 23, 2024, 8:59am UTC](https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615/2 "2024-05-23T08:59:53Z")
Do you have a reproducible example? If you have two objectives, it should be about 2x slower than solving a single LP.
---
### Author:  [@PaterPen](https://discourse.julialang.org/u/PaterPen)
#### Post date: [May 23, 2024, 8:36pm UTC](https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615/3 "2024-05-23T20:36:19Z")
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
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [May 24, 2024, 6:31am UTC](https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615/4 "2024-05-24T06:31:56Z")
> [@PaterPen](#):
>
> On a related note: Is [DiffOpt.jl](https://juliahub.com/ui/Packages/General/DiffOpt) compatible with MOA? @gdalle
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).
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 24, 2024, 7:12am UTC](https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615/5 "2024-05-24T07:12:20Z")
DiffOpt does not work with MOA.
You’ll need to iterate through setting scalar objectives to get the gradient for each objective.
---
### Author:  [@PaterPen](https://discourse.julialang.org/u/PaterPen)
#### Post date: [May 31, 2024, 1:22pm UTC](https://discourse.julialang.org/t/efficient-multi-objective-lp-implementation/114615/6 "2024-05-31T13:22:45Z")
Fyi this is my attempt at differentiating this problem manually @gdalle @odow @mbesancon
```julia
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
```julia
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:
```julia
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…