Multi-objective MINLP problem in JuMP

ashefa · August 23, 2023, 4:12pm

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)?

odow · August 23, 2023, 7:40pm

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

ashefa · August 24, 2023, 8:47am

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).

odow · August 24, 2023, 9:32am

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

ashefa · August 24, 2023, 10:39am

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?

ashefa · August 24, 2023, 10:43am

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)

ashefa · August 24, 2023, 1:01pm

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

odow · August 24, 2023, 5:23pm

Gurobi will be faster

ashefa · August 25, 2023, 3:45pm

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.

odow · August 25, 2023, 10:05pm

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

odow · August 25, 2023, 11:42pm

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)

ashefa · August 26, 2023, 9:18am

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

ashefa · August 26, 2023, 9:25am

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.

JackaChou · August 27, 2023, 1:19am

Infinite time?

ashefa · August 27, 2023, 12:58pm

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 …)

JackaChou · August 28, 2023, 3:08am

Could u please share the code.

ashefa · August 28, 2023, 10:49am

Yes, exactly this

ashefa · August 28, 2023, 11:04am

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.

ashefa · September 28, 2023, 2:50pm

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?

odow · September 28, 2023, 5:43pm

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