MultiObjectiveAlgorithms.EpsilonConstraint runs infinitely

Specific Domains Optimization (Mathematical)
1

I am trying to solve an optimization problem (obtain the efficient frontier of optimal portfolios by maximizing expected return per given level of sharpe ratio). However, the solution runs infinitely (though not always, depending on the data).

Sharpe ratio is defined as (expected_return - Rf) / sqrt(variance) and would like to maximize expected_return for a given level of sharpe ratio.

Any help is much appreciated.

# packages
using JuMP # v1.9.0
using Ipopt # v1.2.0
using MultiObjectiveAlgorithms # v0.1.4
using Statistics

To reproduce the example, I provided the contents of the μ and Q below.

# the content of μ and Q
julia> μ = Vector{Float64}(vec(Statistics.mean(R; dims = 1)))
2-element Vector{Float64}:
  0.006898463772627643
 -0.02972609131603086

julia> Q = Statistics.cov(R)
2×2 Matrix{Float64}:
 0.030446    0.00393731
 0.00393731  0.00713285

# Optimization
μ = Vector{Float64}(vec(Statistics.mean(R; dims = 1)))
Q = Statistics.cov(R)
Rf = 0

model = Model(() -> MultiObjectiveAlgorithms.Optimizer(Ipopt.Optimizer))
# set_silent(model)
set_optimizer_attribute(model, MultiObjectiveAlgorithms.Algorithm(), MultiObjectiveAlgorithms.EpsilonConstraint())
set_optimizer_attribute(model, MultiObjectiveAlgorithms.SolutionLimit(), 25)

@variable(model, 0 <= w[1:2] <= 1)
@constraint(model, sum(w) == 1)
@expression(model, variance, w' * Q * w)
@expression(model, expected_return, w' * μ)

@variable(model, sharpe)

@NLconstraint(model, sharpe == (expected_return - Rf) / sqrt(variance)) 
@objective(model, Max, [expected_return, sharpe]) 
optimize!(model) # => This infinitely runs
2

I don’t think the code runs infinitely, it just tries to find 25 solutions as stated in the SolutionLimit. Probably EpsilonConstraint() requires some tuning to be able to find the required number of solutions or there is some external parameter to set how many runs can happen.

Some code that I tried - If I comment out the line with SolutionLimit and run

using JuMP # v1.9.0
using Ipopt # v1.2.0
using MultiObjectiveAlgorithms # v0.1.4
using Statistics


μ = [  0.006898463772627643
-0.02972609131603086]
Q = [ 0.030446    0.00393731
0.00393731  0.00713285]

Rf = 0

model = Model(() -> MultiObjectiveAlgorithms.Optimizer(Ipopt.Optimizer))
# set_silent(model)
set_optimizer_attribute(model, MultiObjectiveAlgorithms.Algorithm(), MultiObjectiveAlgorithms.EpsilonConstraint())
# set_optimizer_attribute(model, MultiObjectiveAlgorithms.SolutionLimit(), 25)

@variable(model, 0 <= w[1:2] <= 1)
@constraint(model, sum(w) == 1)
@expression(model, variance, w' * Q * w)
@expression(model, expected_return, w' * μ)

@variable(model, sharpe)

@NLconstraint(model, sharpe == (expected_return - Rf) / sqrt(variance)) 
@objective(model, Max, [expected_return, sharpe]) 
optimize!(model) # => This infinitely runs

I get:

* Solver : MOA[algorithm=MultiObjectiveAlgorithms.EpsilonConstraint, optimizer=Ipopt]

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "Solve complete. Found 1 solution(s)"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : [6.89846e-03,3.95355e-02]
  Objective bound    : [6.89846e-03,3.95355e-02]
  Dual objective value : 3.18483e-09

* Work counters
  Solve time (sec)   : 5.39999e-02

Conversely, if I change the algorithm to Dichotomy (for no particular reason, this algorithm was the first on the github documentation):

using JuMP # v1.9.0
using Ipopt # v1.2.0
using MultiObjectiveAlgorithms # v0.1.4
using Statistics


μ = [  0.006898463772627643
-0.02972609131603086]
Q = [ 0.030446    0.00393731
0.00393731  0.00713285]

Rf = 0

model = Model(() -> MultiObjectiveAlgorithms.Optimizer(Ipopt.Optimizer))
# set_silent(model)
set_optimizer_attribute(model, MultiObjectiveAlgorithms.Algorithm(), MultiObjectiveAlgorithms.Dichotomy())
set_optimizer_attribute(model, MultiObjectiveAlgorithms.SolutionLimit(), 25)

@variable(model, 0 <= w[1:2] <= 1)
@constraint(model, sum(w) == 1)
@expression(model, variance, w' * Q * w)
@expression(model, expected_return, w' * μ)

@variable(model, sharpe)

@NLconstraint(model, sharpe == (expected_return - Rf) / sqrt(variance)) 
@objective(model, Max, [expected_return, sharpe]) 
optimize!(model) # => This infinitely runs

I get 25 solutions:

* Solver : MOA[algorithm=MultiObjectiveAlgorithms.Dichotomy, optimizer=Ipopt]

* Status
  Result count       : 25
  Termination status : OPTIMAL
  Message from the solver:
  "Solve complete. Found 25 solution(s)"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : [-2.97261e-02,-3.51970e-01]
  Objective bound    : [6.89846e-03,3.95355e-02]
  Dual objective value : 3.18483e-09

* Work counters
  Solve time (sec)   : 5.39999e-02
3

Thanks very much @blob for the answer. This does not look like a proper solution. I don’t know how EpsilonConstraint differs from Dichotomy, but the solution looks strange.

image

4

I’ve opened an issue to take a look at this: EpsilonConstraint doesn't terminate · Issue #50 · jump-dev/MultiObjectiveAlgorithms.jl · GitHub

But the SolutionLimit option is a little weird for EpsilonConstraint:

help?> MOA.EpsilonConstraint
  EpsilonConstraint()

  EpsilonConstraint implements the epsilon-constraint algorithm for bi-objective programs.

  Supported optimizer attributes
  ================================

    •  MOA.EpsilonConstraintStep(): EpsilonConstraint uses this value as the epsilon by which
       it partitions the first-objective's space. The default is 1, so that for a pure integer
       program this algorithm will enumerate all non-dominated solutions.

    •  MOA.SolutionLimit(): if this attribute is set then, instead of using the
       MOA.EpsilonConstraintStep, with a slight abuse of notation, EpsilonConstraint divides
       the width of the first-objective's domain in objective space by SolutionLimit to obtain
       the epsilon to use when iterating. Thus, there can be at most SolutionLimit solutions
       returned, but there may be fewer.
1 Like