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