Multiobjective support in JuMP

[odow]

This issue is for a discussion on introducing multiobjective support to JuMP (and MOI).

Passing multiple objectives to the solver

There are two possibilities:

1. Extending @objective to support vector-valued functions, so the user would write:

   @objective(model, Min, C * x)
   

2. Adding an index keyword, so the user would write:

   @objective(model, Min, C[1, :] * x, index = 1)
   @objective(model, Min, C[2, :] * x, index = 2)
   

Option (1) is probably the easiest, since it fits quite well into the MOI framework, and JuMP already has support for parsing vector-valued functions.

M.O. solvers would just declare:

function MOI.supports(
    ::Optimizer, ::MOI.ObjectiveFunction{MOI.VectorAffineFunction{Float64}}
)
    return true
end

The biggest downside with (1) is that it would only support linear and quadratic objectives. We wouldn't be able to directly support nonlinear objectives. (You could, of course, introduce a dummy variable with a nonlinear equality constraint.)

Querying results from the solver

I think the easiest way for M.O. solvers to return the Pareto frontier is for them to sligthly abuse the notion of ResultCount.

We should implement

function result_count(model::Model)
    return MOI.get(model, MOI.ResultCount())
end

And then add keyword arguments to value, objective_value, so we have:

function value(x::VariableRef; result=1)
    return MOI.get(owner_model(x), MOI.VariablePrimal(result), index(x))
end

The only issue with objective_value (and objective_bound and dual_objective_bound) is that they expect a Float64 return type. This would be relaxed to depend on the type of the objective function.

Example

If MultiJuMP implemented a MOI.Optimizer solver, we could write:

using JuMP
using MultiJuMP

model = Model(() -> MultiJuMP.Optimizer(GLPK.Optimizer))
@variable(model, x[1:2] >= 0)
@objective(model, Min, [2x[1] + x[2]; x[1] + 2x[2]])
@constraint(model, sum(x) >= 1)
optimize!(model)

pareto_frontier = []
for i = 1:result_count(model)
    @assert primal_status(model; result = i) == MOI.OPTIMAL
    push!(pareto_frontier, (
        x_val = value.(x; result = i),
        obj = objective_value(model; result = i)
    ))
end

cc @anriseth, @matbesancon

[blegat]

Another option is to define

function MOI.set(model, ::ObjectiveFunction, func::AbstractVectorFunction)
   for i in 1:MOI.output_dimension(func)
       MOI.set(model, MultiObjectiveFunction(i), func[i]) 
    end
end

so that both syntax works for the user

[odow]

The more I've thought about this, the more sure I am that 1 is the correct way. You really are declaring a vector valued objective, instead of multiple different objectives. It causes a natural interpretation of the return type of ObjectiveValue and needs no change to MOI.

[anriseth]

Great to see a renewed and improved effort to bring vector-optimzation to JuMP :) I've left academia and, therefore, sadly stopped using Julia. I don't think I can provide much useful input, but look forward to trying these features at some point in the future:)

[xgandibleux]

Thanks @odow for having initiated this discussion, and I am happy to see it.

We are developping vOptSolver which is currently composed of two packages, vOptSpecific and vOptGeneric. vOptGeneric is devoted to multi-objective linear optimization problems with discrete variables, see https://github.com/vOptSolver (but we have also some beta versions of algorithms for problems with mixed integer variables and continuous variables).

The version online is compliant with julia 1.x and JuMP 0.18. You can give a look on the notebook (here : https://github.com/vOptSolver/vOptSolver-notebook) that I have prepared as support for a seminar held in November 2018. The version of vOptGeneric compliant with JuMP 0.19 is ready on a branch (https://github.com/vOptSolver/vOptGeneric.jl/branches). It is a matter of free time for me for updating all the documents (for some reasons I have been much less available since the last summer).

Here an example with vOptGeneric+JuMP0.19 for the bi-objective knapsack problem:

using vOptGeneric, JuMP, Cbc, LinearAlgebra
m = vModel(with_optimizer(Cbc.Optimizer, logLevel=0))

p1 = [77,94,71,63,96,82,85,75,72,91,99,63,84,87,79,94,90,60,69,62]
p2 = [65,90,90,77,95,84,70,94,66,92,74,97,60,60,65,97,93,60,69,74]
w = [80,87,68,72,66,77,99,85,70,93,98,72,100,89,67,86,91,79,71,99]
c = 900
size = length(p1)

@variable(m, x[1:size], Bin)
@addobjective(m, Max, dot(x, p1))
@addobjective(m, Max, dot(x, p2))
@constraint(m, dot(x, w) <= c)

vSolve(m, method=:dichotomy)

Y_N = getY_N(m)
for n = 1:length(Y_N)
    X = value.(x, n)
    print(findall(elt -> elt ≈ 1, X))
    println("| z = ",Y_N[n])
end

As you can see:

1. we have chosen to handle the objectives separately. Adding an index which denote a specific objective as suggested by @odow is our preferred option because we would like in a near future to remove one objective without redefining all the model.
2. with vSolve, we specify the solving method used for computing the exact solutions (X_E) and the corresponding non-dominated points (Y_N). Currently four methods are available (epsilon-constraint; dichotomy; Chalmet; lexicographic), and these methods are currently limited to the bi-objective case (except for the lexicographic ).
3. we have developed some functions for extracting Y_N and X_E.

Thus:

• a discussion about how to pass the resolution method to invoke (and possibly the required parameters) is needed.

• a discussion about how to retrieve Y_N and X_E is needed as soon as the problem is not discrete. Indeed, for the MO-MILP case, Y_N may be composed of points, segments (open or closed), facets (full or not; open or closed). We have not yet stated about a protocol in the vOpt team. One option may be to return each object separately, with its proprieties, and a description of their adjacency relations.

[blegat]

a discussion about how to pass the resolution method to invoke (and possibly the required parameters) is needed.

vOpt would can either create

• one MOI optimizer per resolution method, so the user has do to model = Model(with_optimizer(vOptSolver.Chalmet.Optimizer) or
• one single MOI optimizer and the resolution method is given as parameter, so the user ahs to do model = Model(with_optimizer(vOptSolver.Optimizer) and then set_parameter(model, "method", "Chalmet").

a discussion about how to retrieve Y_N and X_E is needed as soon as the problem is not discrete

You can freely define any interface, JuMP will not restrict you. You just need to define new MOI attributes, so that the user can do MOI.get(model, VOptSolver.Y_N()) and then you can define convenience function with the interface you like, e.g.: getY_N(model) = MOI.get(model, VOptSolver.Y_N()).

[odow]

we have chosen to handle the objectives separately. Adding an index which denote a specific objective as suggested by @odow is our preferred option because we would like in a near future to remove one objective without redefining all the model.

Is it okay to just redefine the objective? (E.g., below.) Or are you keeping the pre-computed frontier and projecting it onto the lower-dimensional space when you remove an objective?

@objective(model, Min, C * x)
# Then later
@objective(model, Min, C[1, :] * x)

If X_E and Y_N are points, then you will be able to query them from JuMP using value(x, result=i) and objective_value(model, result=i) as I outlined. If you have more solver-specific solutions, it's easy to define MOI attributes as @blegat mentioned. For example:

model = Model()
@variable(model, x[1:2])
@objective(model, Min, C * x)
Y_N = MOI.get(model, vOpt.EfficientFacets())

Gurobi.jl uses this to query the IIS, for example.

There is one very big change for vOptSolver: you should re-write it as a MOI solver, rather than a JuMP-extension. That means you don't have to worry about macros like @addobjective, and defining types like vModel which wrap JuMP models. You would have first-class support for everything JuMP offers :)

[joaquimg]

From the user/modeler perspective it seems that option 2 is better. Mode one is sort of not natural mainly if the two objectives are completely different.

[odow]

So my reasoning is this:

• MathOptInterface can support multiple objectives via vector-valued objective functions
• MO solvers should declare support for MOI.ObjectiveFunction{MOI.VectorOfVariables}, MOI.ObjectiveFunction{MOI.VectorAffineFunction} or MOI.ObjectiveFunction{MOI.VectorQuadraticFunction}.
• If they do this, then things like bridges should JustWork™
• MathOptInterface can interpret the Pareto frontier using the ResultCount attribute. We already have support for this, and all VariablePrimal etc are parameterized by the ResultCount.

So, at the MOI level, we don't need to change anything to support M.O.

One exception is if solvers want to return more complicated solution artifacts (e.g., facets). If so, they should define solver-specific MOI attributes.

So, given we don't need to change MOI, what is the closest mapping to JuMP? Vector-valued objectives.

If you want to give the objectives separately, you can go:

@expression(model, first_objective, ...)
@expression(model, second_objective, ...)
@objective(model, Min, [first_objective; second_objective])

The other approach, JuMP giving objectives separately, would require changes to MOI to parameterize the objective function based on the index. This raises some questions.

• What does ObjectiveValue return?
• What if users give ObjectiveFunction{1} and ObjectiveFunction{3}?
• What if one of the objectives is a vector-valued function?
• How do we delete an objective? MOI only has the concept of setting attributes, so we could set it an objective to zero, but then does the user expect the dimension of the returned ObjectiveValue to change? Or do they expect a 0 in that element?

[joaquimg]

expressions are reasonable for a quick fix, although doing separately seems a good future goal. Delete should delete and not set to zero, otherwise it could grow a lot...

[matbesancon]

The inconvenience I see to the @objective(model, Min, C * x) syntax is that people can get confused when accidentally providing a vector function, compared to having a multi-objective specific syntax

[odow]

The inconvenience I see to the @objective(model, Min, C * x) syntax is that people can get confused when accidentally providing a vector function, compared to having a multi-objective specific syntax

To which I counter: there is already precedent for vector-valued functions in constraints; most solvers would throw "VectorAffineFunction in objective not supported"; and even if they did solve it, objective_value would return a vector not a Float64.

People could get equally confused with

@objective(model, Min, x)
@objective(model, Min, y)

And how would we handle multiple objective senses?

Here is @xgandibleux's example with what I am proposing:

using vOptGeneric, JuMP, Cbc, LinearAlgebra

p1 = [77,94,71,63,96,82,85,75,72,91,99,63,84,87,79,94,90,60,69,62]
p2 = [65,90,90,77,95,84,70,94,66,92,74,97,60,60,65,97,93,60,69,74]
w = [80,87,68,72,66,77,99,85,70,93,98,72,100,89,67,86,91,79,71,99]
c = 900
size = length(p1)

model = Model(with_optimizer(
    vOptGeneric.Optimizer(
        Cbc.Optimizer(logLevel=0), method=:dichotomy
    )
))
@variable(model, x[1:size], Bin)
@objective(model, Max, [dot(x, p1), dot(x, p2)])
@constraint(m, dot(x, w) <= c)
optimize!(model)
for n = 1:result_count(model)
    X = value.(x; result = n)
    print(findall(elt -> elt ≈ 1, X))
    println("| z = ", objective_value(model; result = n))
end

[xgandibleux]

In writing

@objective(model, Max, [dot(x, p1), dot(x, p2)])

the user cannot (1) blend functions to minimize and to maximize in the model, and (2) point out a function to delete without rebuilding all the model.

• Concerning (1): of course, it is trivial to rewrite the function from min to max, but having the choice here is natural (same level of why not express the constraints only in <=).

• Concerning (2): MO-models are also solved within an interactive approach, where the user can adopt a "what-if" behavior with the objectives. This is a feedback that I have received from users of vOptGeneric.

I am working with my students on a case study (a car sequencing problem) with 2 or 3 objective functions. The (technical) constraints are static but the objectives may change between production plants. This is an other illustration where preferences of a decision-maker may suggest to delete an objective or replace an objective by an other (with the idea of comparing different optimization policy) in a (large) MIP model.

Is it technically possible to imagine of naming the objectives as the constraints:

@objective(model, obj1, Max, [dot(x, p1)])
@objective(model, obj2, Min, [dot(x, p2)])

or

@objective(model, cost[1], Min, [dot(x, p1)])
@objective(model, stability[2], Max, [dot(x, p2)])

or

@objective(model, [1], Min, [dot(x, p1)])
@objective(model, [2], Max, [dot(x, p2)])

or

@objective(model, obj[i = 1:3], [ dot(x, p[i]) ])

and then to refer to one objective by name (e.g. stability) or by its index (e.g. 2)?

[odow]

blend functions to minimize and to maximize in the model

This is a reasonable point to make. The question of "objective sense" isn't obvious. The hang-up we seem to be having is the standard form. Is it:

min: f_0(x)
s.t. f_i(x) in S_i

where f_0 can be vector-valued, or is it

sense_i: f_i(x), i=1,2,...
s.t. g_j(x) in S_j, j=1,2,...

where f_i must be scalar-valued.

I'm going to argue strongly that it's the first, since it is the simplest concept. More importantly, it's the design of MathOptInterface as it stands today. We wouldn't need to add any special features for multi objective problems. This is pretty important for actually implementing this, because I don't think we want to re-think the entire concept of objectives this close to the release of JuMP 1.0.

So at the expense of (i) providing a list of expressions not being the nicest syntax and (ii) not being able to specify different objective senses, perhaps we can table this idea until JuMP 2.0?

point out a function to delete without rebuilding all the model.

They can just reset the objective:

model = Model()
@variable(model, x[1:2])
@objective(model, Min, [x[1], x[2]])
optimize!(model)
@objective(model, Min, x[1] + x[2])
optimize!(model)

Is it technically possible to imagine of naming the objectives as the constraints:

This would be possible using @expression:

model = Model()
@variable(model, x[1:2])
@expression(model, obj[i=1:2], x[i])
@objective(model, Min, obj)
optimize!(model)
value(obj[1], result=1)

Relating to the changing objectives point:

model = Model()
@variable(model, x[1:3])
@expression(model, stability, x[1])
@expression(model, cost, x[2])
@expression(model, production, x[3])
@objective(model, Max, [stability, production, -cost])
optimize!(model)
@show value(stability; result=2)
@show value(cost; result=2)
@objective(model, Max, [stability, production])
optimize!(model)
# Note: we didn't have cost as an objective, but we can still query the value!
@show value(cost; result=2)

[odow]

Okay, so I started to implement this, and I found that one thing would have to change in MOI. Essentially, we just need to relax the AbstractScalarFunction aspect of this definition:

"""
    ObjectiveFunction{F<:AbstractScalarFunction}()

A model attribute for the objective function which has a type `F<:AbstractScalarFunction`.
`F` should be guaranteed to be equivalent but not necessarily identical to the function type provided by the user.
Throws an `InexactError` if the objective function cannot be converted to `F`,
e.g. the objective function is quadratic and `F` is `ScalarAffineFunction{Float64}` or
it has non-integer coefficient and `F` is `ScalarAffineFunction{Int}`.
"""
struct ObjectiveFunction{F<:AbstractScalarFunction} <: AbstractModelAttribute end

https://github.com/JuliaOpt/MathOptInterface.jl/blob/1d4f0d3fd0a545eb2c926731b7c428963092974d/src/attributes.jl#L804-L813

[odow]

So I talked to @mlubin, and we won't be pushing this into JuMP unless there is general agreement on the approach. We definitely don't want to force @xgandibleux et al. to implement our (edit) my approach when they are the subject matter experts.

However, since I require this for my personal work, I will move forward with an implementation. Once it's ready for review, we can have another discussion about the pro's and con's of each approach, this time with some actual working code. But, to re-iterate, we won't merge it unless there is general agreement.

[odow]

Okay, I now have a proof of concept multi-objective solver: https://github.com/odow/MOO.jl

The best place to look is the test, where we solve a trivial 2 variable BOLP: https://github.com/odow/MOO.jl/blob/4a1efcbfe04802ca477932b2c1d7fca55eaccb16/test/nise.jl

Arguments for min VectorFunction as opposed to Vector[min ScalarFunction1, max ScalarFunction2]:

• Aside from a 1 line change in MOI (the rest of the diff is for MOI.Utilities), everything JustWorks™. Note in particular that we can use loadfromstring! without modifications.
• The use of ResultCount to represent the solver returning a set of Pareto-optimal solutions
• The interpretation of ObjectiveBound as the utopian point
• It's nice that the entire implementation of the solver is only ~200 lines, and most of those are lines which forward methods onto the inner solver.
• The MOP format, the VLP format, Gurobi, and CPLEX all assume min/max F(x)

If we want to move forward with this, now that JuMP has support for accessing multiple solutions (https://github.com/JuliaOpt/JuMP.jl/pull/2100), the only missing piece is JuMP passing vector-valued functions to solvers as the objective function.

[freemin7]

I recently had a conversation with Kaisa Miettinen who leads the multiobjective optimization group at University of Jyväskylä. (She made me aware of @xgandibleux 's work using Julia) The reason for using multi objective optimization according to her is that the proper weighing of the objectives depends on querying a decision maker. There a 4 ways how that might be done:

• The a proper (method of) weighing is know ahead of time and the problem can be solved as a single objective optimization
• A surrogate model of the pareto front is required so that a decision maker can judge/explore the space computationally cheaply
• After some sampling the decision maker get queried to might judge which objectives/parts of the parameter space are interesting to be explored more in case the space of options is really huge
• A sampling of the pareto front is required such that a decision maker can choose between those samples, possibly with some software that helps him explore that space.

I see that 1&4 can be accomplished with the current frame work. 3 could be accomplished using callbacks. Although a standardized interface for user interaction would be really nice so that many Multi objective solvers can share one exploratory interface (or multiple interfaces could compete in a far away future). While i imagine that 2 could be accomplished using either callback or parsing the multiple results but i imagine that a solver might have more useful information for constructing surrogate than just the set of the pareto front. (I can think right now of local gradient information on the pareto manifold (think tangent space) or rejected points close to the pareto front so that the surrogate doesn't touch them.)

If consider those features non-essential but think they are low hanging fruits on the way of JuMP becoming a best in class multi objective framework. Although my knowledge on multi objective optimization is second hand only and i am open to be corrected.

[manuelbb-upb]

I would like to add my two cents to the discussion, mainly concerning way objectives are declared. (tldr: I go with @xgandibleux and am in favor of declaring them individually).

Background: I have recently worked on a trust region solver for nonlinear multiobjective problems. This solver is specifically aimed at expensive or heterogeneous problems where some objectives are considered blackbox and might not even be Julia functions. Hence, it is not necessarily something I would write a MOI interface for.

Despite this fact, the need for separate evaluation of the different objectives might arise in other approaches as well. • There is, of course, the Pascoletti-Serafini scalarization, a very powerful tool that can be used on its own or be integrated in iterative schemes. It requires an utopia point in objective space that results from the minimization of the individual objectives. • I guess for the ε-constraint method, accessibility of the individual objectives would be of use too. • Suppose, for a you want to calculate the multiobjective Newton direction. This requires positive definite Hessians which is why often times BFGS approximations are used. A sophisticated solver could allow for the user to specify (or automatically determine) how the Hessian should be calculated/approximated for each objective. • I know of investigations into the hierarchical structure of the Pareto Set/Front, where reduced problems were used, i.e., where one or several objectives were dropped.

Hope this helps in evaluating how to proceed :)

[odow]

@/all: Thanks for the input.

If anyone watching this thread (or your students) is interested in progressing this further, I've proposed this as a Google Summer of Code (GSOC) for this (Northern) summer: https://github.com/jump-dev/GSOC2021#multiobjective-optimization

You can find more about GSOC here:

• https://summerofcode.withgoogle.com
• https://github.com/numfocus/gsoc

If we get selected, student applications run March 30 - April 14.

[jacktang]

Hi @odow , how is the progress about MOO? or recommend to use https://github.com/anriseth/MultiJuMP.jl?

[odow]

No progress. MultiJuMP hasn't been updated to JuMP 1.0. Digging in to understand how it works and updating it to JuMP 1.0 would be a good project, if you have the time.

[xgandibleux]

vOptGeneric solves generic linear problems with 2 objectives with e.g. an epsilon constraint method (+ MIP solver available from JuMP). It is compliant with the last versions of julia and JuMP. A quick starter is presented in recent slides here. Several examples are provided here. Several extensions/new features are under developments, I hope to come back soon with news. Feedbacks are welcome.

[odow]

I'm back at this in #3176.

Here's a related project that we should look to for ideas: https://pymoo.org. It's a bit limited though. It supports only minimize of a vector of objectives.

[odow]

Single objective sense

• Gurobi, CPLEX, and Xpress support multiple objectives, but require weights (and/or priorities)
• CPLEX's extension to the LP file format and MPS file format (can give negative weights to flip)
• pymoo (minimize only)

Mixed objective sense

• https://github.com/vOptSolver
• LocalSolver has a lexicographic solver that lets you mix and match Min and Max: https://www.localsolver.com/docs/last/modelingfeatures/multiobjectiveresolution.html
• AMPL lets you formulate multiple objectives, but you must pick one prior to sending to a solver. It doesn't support multi-objective solvers
• Pyomo lets you create multiple objectives, but you need to deactivate N-1 before calling the solver

[rschwarz]

If I may comment on Gurobi's capability: You can specify both weights and priorities for the objectives. The objectives will then be grouped by the priority, and each group will be "blended" with a weighted sum.

In addition, the user can specify tolerances (both rel. and abs.?) to allow for some relaxation when going one level deeper into the hierarchy. This feature in particular is quite useful in practice, we found, as you'd otherwise quickly end up in extreme corners.

[odow]

You can specify both weights and priorities for the objectives

Yeah. I haven't tested, but this should work:

set_optimizer_attribute.(
    model, 
    Gurobi.MultiObjectivePriority.(1:3), 
    [2, 1, 1],
)
set_optimizer_attribute.(
    model, 
    Gurobi.MultiObjectiveWeight.(1:3), 
    [1.0, 0.25, 0.75],
)

I think these are attributes of the optimizer though. They're unrelated to how we model the objectives at the JuMP level.

[rschwarz]

I think these are attributes of the optimizer though. They're unrelated to how we model the objectives at the JuMP level.

Fair enough, I'm not sure about the scope of this issue. If this is just about defining attributes that some of the solvers may support, there's not much more to say about it.

But the algorithm for hierarchical objectives is easy enough to implement on top of any (LP/MIP) solver, so it could make sense to implement some kind of wrapper for it.

[odow]

But the algorithm for hierarchical objectives is easy enough to implement on top of any (LP/MIP) solver, so it could make sense to implement some kind of wrapper for it.

Yeah. That's my plan for: https://github.com/odow/MOO.jl

Currently it just has a single algorithm, but the idea is to collate a few simple ones.

[xgandibleux]

Before the integration of MOP features into CPLEX and Gurobi, we have integrated in vOptGeneric two methods aiming to solve weighted sum and the lexicographic LP/MIP with any solver interfaced to JuMP. The methods to invoke are

:dicho or :dichotomy :lex or :lexico

These last months, we have extended HiGHS to 2LP (expected by several colleagues for research needs) and developped the epsilon-contraint method for 3 objectives based on the Kirlik-and-Sayin algorithm (vOptGeneric is mainly used today for solving bi-objective problems with the epslion-contraint method with declared with JuMP). Both should be integrated to vOptGeneric shorthly (it is a matter of freetime for us).

BTW, FICO has also announced recently the integration of MOP features into XPRESS.

[odow]

These last months, we have extended HiGHS to 2LP

Nice :smile: Is there a link to the code? Does it support independent objective senses?