In the continuation of this comment and this post, I was wondering about the implementation of a multi-objective MINLP problem in JuMP using MOA.

Since @Arpit_Bhatia used `MOA`

for a nonlinear multi-objective portfolio optimization problem (by Ipopt solver) and in the JuMP’s main example for multi-objective optimization, it has been applied to a MIP (by HiGHS solver), now how we can combine them to use `MOA`

for MINLP (e.g. by Juniper)?

Do you have an example that doesn’t work? It should be fine.

```
julia> using JuMP
julia> import Ipopt
julia> import Juniper
julia> import MultiObjectiveAlgorithms as MOA
julia> ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
MathOptInterface.OptimizerWithAttributes(Ipopt.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("print_level") => 0])
julia> juniper = optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt)
MathOptInterface.OptimizerWithAttributes(Juniper.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("nl_solver") => MathOptInterface.OptimizerWithAttributes(Ipopt.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("print_level") => 0])])
julia> model = Model(() -> MOA.Optimizer(juniper))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: MOA[algorithm=MultiObjectiveAlgorithms.Lexicographic, optimizer=Juniper]
julia> set_silent(model)
julia> set_attribute(model, MOA.Algorithm(), MOA.Dichotomy())
julia> set_attribute(model, MOA.SolutionLimit(), 10)
julia> @variable(model, 1 <= x[1:5] <= 5, Int)
5-element Vector{VariableRef}:
x[1]
x[2]
x[3]
x[4]
x[5]
julia> @variable(model, obj[1:2])
2-element Vector{VariableRef}:
obj[1]
obj[2]
julia> @objective(model, Max, obj)
2-element Vector{VariableRef}:
obj[1]
obj[2]
julia> @NLconstraint(model, obj[1] <= sum(log(x[i]) for i in 1:5))
(obj[1] - (log(x[1]) + log(x[2]) + log(x[3]) + log(x[4]) + log(x[5]))) - 0.0 ≤ 0
julia> c = rand(5)
5-element Vector{Float64}:
0.9564024316877537
0.3761897861649184
0.9402157284498638
0.5496119699788227
0.2809383737680249
julia> @NLconstraint(model, obj[2] <= sum(c[i] * x[i] for i in 1:5))
(obj[2] - (0.9564024316877537 * x[1] + 0.3761897861649184 * x[2] + 0.9402157284498638 * x[3] + 0.5496119699788227 * x[4] + 0.2809383737680249 * x[5])) - 0.0 ≤ 0
julia> @constraint(model, sum(x) <= 20)
x[1] + x[2] + x[3] + x[4] + x[5] ≤ 20
julia> optimize!(model)
julia> solution_summary(model)
* Solver : MOA[algorithm=MultiObjectiveAlgorithms.Dichotomy, optimizer=Juniper]
* Status
Result count : 10
Termination status : LOCALLY_SOLVED
Message from the solver:
"Solve complete. Found 10 solution(s)"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Objective value : [6.93147e+00,-8.06555e+03]
Objective bound : [6.93147e+00,1.40168e+01]
Relative gap : 0.00000e+00
* Work counters
Solve time (sec) : 4.35057e-01
```

Thanks for the response. Consider the following code:

```
A = rand(1,10)
B = rand(1,10)
T=length(A)
M=10^6
c_S = 60
c_R = 100
using JuMP
import Juniper
import HiGHS
import MultiObjectiveAlgorithms as MOA
model = JuMP.Model(() -> MOA.Optimizer(HiGHS.Optimizer))
set_attribute(model, MOA.Algorithm(), MOA.Dichotomy())
set_attribute(model, MOA.SolutionLimit(), 5)
@variable(model, S >= 0);
@variable(model, 1 <= R <= 3);
@variable(model, y1[t = 1:T], Bin);
@variable(model, y2[t = 1:T], Bin);
@variable(model, Y[t = 1:T] >= 0);
@variable(model, O[t = 1:T] >= 0);
@variable(model, V[t = 0:T] >= 0);
@variable(model, Q[t = 1:T] >= 0);
@constraint(model, V[0] == 0);
@constraint(model, [t in 1:T], Q[t] == B[t] * R);
@constraint(model, [t in 1:T], Y[t] <= A[t]);
@constraint(model, [t in 1:T], Y[t] <= V[t-1]);
@constraint(model, [t in 1:T], Y[t] >= A[t] - M * (1 - y1[t]));
@constraint(model, [t in 1:T], Y[t] >= V[t-1] - M * y1[t]);
@constraint(model, [t in 1:T], V[t] <= V[t-1] + Q[t] - O[t] - Y[t]);
@constraint(model, [t in 1:T], V[t] <= S - Y[t]);
@constraint(model, [t in 1:T], V[t] >= V[t-1] + Q[t] - O[t] - Y[t] - M * (1 - y2[t]));
@constraint(model, [t in 1:T], V[t] >= S - Y[t] - M * y2[t]);
@objective(model, Max, [-c_S * S - c_R * R, sum(Y[t] for t in 1:T)])
optimize!(model)
solution_summary(model)
println(objective_value(model))
println("S = ", value.(S), "\n","A = ",value.(R))
```

Now I want to have a solution when I am changing the first objective function to the quadratic term `-c_S * S^2 - c_R * R`

(similar to the approach of @Arpit_Bhatia in NL MO portfolio optimization).

You need to use a solver that supports mixed integer quadratic programs. Use Juniper like i did above, or use Gurobi.

When I just replaced `HiGHS`

with `Ipopt`

and `Juniper`

, without changing the objective function into a quadratic form, the code takes lots of time to run. It’s still running! What do you think about?

I just replaced

```
import HiGHS
import MultiObjectiveAlgorithms as MOA
model = JuMP.Model(() -> MOA.Optimizer(HiGHS.Optimizer))
set_attribute(model, MOA.Algorithm(), MOA.Dichotomy())
set_attribute(model, MOA.SolutionLimit(), 5)
```

with

```
import Ipopt
import Juniper
import MultiObjectiveAlgorithms as MOA
ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
juniper = optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt)
model = Model(() -> MOA.Optimizer(juniper))
set_silent(model)
set_attribute(model, MOA.Algorithm(), MOA.Dichotomy())
set_attribute(model, MOA.SolutionLimit(), 10)
```

Do you have any idea about difference between running time of these two set of solvers for this problem?

Gurobi will be faster

When I used the `Gurobi`

, I got this error:

MathOptInterface.SetAttributeNotAllowed{MathOptInterface.TimeLimitSec}: Setting attribute MathOptInterface.TimeLimitSec() cannot be performed. You may want to use a

`CachingOptimizer`

in`AUTOMATIC`

mode or you may need to call`reset_optimizer`

before doing this operation if the`CachingOptimizer`

is in`MANUAL`

mode.

This looks like a bug in Gurobi.jl. Let me fix it: Fix support for MOI.TimeLimitSec by odow · Pull Request #518 · jump-dev/Gurobi.jl · GitHub

Fixed if you update Gurobi.jl to the newest version (v1.0.3)

Here’s how I would write your model:

```
using JuMP
import Gurobi
import MultiObjectiveAlgorithms as MOA
A = rand(1,10)
B = rand(1,10)
T = length(A)
M = 10^6
c_S = 60
c_R = 100
model = Model(() -> MOA.Optimizer(Gurobi.Optimizer));
set_silent(model)
set_attribute(model, MOA.Algorithm(), MOA.Dichotomy())
set_attribute(model, MOA.SolutionLimit(), 5)
@variables(model, begin
0 <= S
1 <= R <= 3
y1[1:T], Bin
y2[1:T], Bin
0 <= Y[t = 1:T] <= A[t]
0 <= O[1:T]
0 <= V[0:T]
end);
fix(V[0], 0; force = true)
@expressions(model, begin
Q[t in 1:T], B[t] * R
end)
@constraints(model, begin
[t in 1:T], Y[t] <= V[t-1]
[t in 1:T], Y[t] >= A[t] - M * (1 - y1[t])
[t in 1:T], Y[t] >= V[t-1] - M * y1[t]
[t in 1:T], V[t] <= V[t-1] + Q[t] - O[t] - Y[t]
[t in 1:T], V[t] <= S - Y[t]
[t in 1:T], V[t] >= V[t-1] + Q[t] - O[t] - Y[t] - M * (1 - y2[t])
[t in 1:T], V[t] >= S - Y[t] - M * y2[t]
end);
@objective(model, Max, [-c_S * S - c_R * R, sum(Y)])
optimize!(model)
solution_summary(model)
```

It’s completely solved. I really appreciate you taking the time to rewrite my code into a nice structure.

My problem was solved, but since `Ipopt`

+ `Juniper`

got an uncommon time (for me infinite) for running the code, It might be troubleshooting to have a look at it.

Infinite time?

Yes, when I ran the above code with `Juniper`

and `Ipopt`

as the NL solver, I did not get response (my kernel was running for a long time …)

Could u please share the code.

Using the code written by @odow solved the issue!

**But still**, for my real case with vector size `A`

of 365*10, the solution time is infinite.

Just one question about quadratic term in the multi-objective optimization. As it is mentioned in this post “HiGHS supports quadratic objectives”, but I got `BoundsError: attempt to access 0-element Vector{MultiObjectiveAlgorithms.SolutionPoint} at index [1] `

for this type of objectives in the multi-objective mode! Can you explain the reason?

Your problem didn’t find a solution. What does `solution_summary(model)`

say?