An idea to explain a drastic difference in computation times for optimizing a JuMP model, and the unexpected ranking of solvers used?

I have to solve a basic parametric MILP.:

\begin{equation} \begin{split} \mbox{maximize} & \quad \alpha \\ \mbox{s.t.} & \quad \alpha \leq - \lambda_j(\psi)(r_{*j} - z_j(x)), \,\, j=1,\dots,d \\ & \quad x\in \mathbb X \subset\mathbb R^{n}, \ \alpha \ge 0 \end{split} \end{equation}

The parameter is a weight value \lambda_j(\psi) of dimension d from 2 to 6.
The number of weights to consider varies between 100 and 10000.
I am repeating the optimization 20 times.

To avoid to rebuild the model from scratch, I modify it at each iteration, in deleting and adding d constraints concerned by a weight.

I implemented it, and I solved instances respectively with GLPK, HiGHS, Gurobi, and CPLEX. The results collected are valid (for the 4 MIP solvers) but the time consumed are drastically differents… in favor of GLPK.
Ok, my comparison is unfair; the preprocessing is not disabled for the others solvers…but the elapsed times are very amazing:

first run (with 1000 weights) to compile:
GLPK.Optimizer | time(s)= 2.8956220149993896

with 1000 weights:
GLPK.Optimizer | time(s)= 0.2676060199737549
GLPK.Optimizer | time(s)= 0.260451078414917
GLPK.Optimizer | time(s)= 0.2592289447784424

HiGHS.Optimizer | time(s)= 15.928501844406128
HiGHS.Optimizer | time(s)= 15.65535306930542
HiGHS.Optimizer | time(s)= 15.687047004699707

Gurobi.Optimizer | time(s)= 1.675018072128296
Gurobi.Optimizer | time(s)= 1.4513440132141113
Gurobi.Optimizer | time(s)= 1.4410531520843506

CPLEX.Optimizer | time(s)= 105.86109781265259
CPLEX.Optimizer | time(s)= 118.24663805961609
CPLEX.Optimizer | time(s)= 116.49558091163635

with 10000 weights:
GLPK.Optimizer | time(s)= 12.330639839172363
HiGHS.Optimizer | time(s)= 1553.121855020523
Gurobi.Optimizer | time(s)= 154.71333193778992
too long with CPLEX

I am wondering about the cause of these differences; is it

  • a consequence of by my approach for modifying the model?
  • inhered to the MIP solver activity?
  • something internal to the packages bridging MOI and the corresponding solver?

Any opinion/explanation/advice?

Here the code (extracted from my implementation and simplified to illustrate this) used to report here the elapsed times:

using JuMP, GLPK, HiGHS, Gurobi, CPLEX

function L(  solver::DataType, 
             p::Matrix{Int64}, w::Vector{Int64}, c::Int64,
             rp::Vector{Int64}, λ_ψ::Vector{Vector{Float64}}
          )

    L = (Float64)[]  # a series of L(S,r*,ψ) values
    
    d,n = size(p)       # number of objectives and number of variables
    
    # define a JuMP model
    model = Model(solver)  
    set_silent(model)

    # setup X according (1)
    @variable(model, x[1:n], Bin)                                    # the n binary variables of 01UKP
    @constraint(model, sum(w[i] * x[i] for i in 1:n) ≤ c)            # the constraint of 01UKP

    # declare the objective functions
    @expression(model, z[k=1:d], sum(p[k,j] * x[j] for j in 1:n))    # the d objectives of 01UKP

    # ------------------------------------------------------------------------------
    # Concerning the Chebychev part of the model

    @variable(model, α ≥ 0)        
    @objective(model, Max, α)

    for i in 1:length(λ_ψ)

        # add the d constraints for a given λ_ψ to the model
        @constraint(model, con[k=1:d], α ≤ - λ_ψ[i][k] * (rp[k] -  z[k]) )
        
        optimize!(model)
        @assert is_solved_and_feasible(model) "Error: optimal solution not found"

        push!(L,objective_value(model))
        
        for k=1:d
            delete(model, con[k])
        end
        unregister(model, :con)

    end
    
    return L
end


p = [ 13 10  3 16 12 11  1  9 19 13 ;   1 10  3 13 12 19 16 13 11  9  ]   
w = [ 4, 4, 3, 5, 5, 3, 2, 3, 5, 4  ]   
c = 19                                   

rp = [40,40]
λ_ψ = fill([1.0238499551062132, 4.66018735258353],1000)
solver = GLPK.Optimizer
#solver = HiGHS.Optimizer
#solver = Gurobi.Optimizer
#solver = CPLEX.Optimizer

start = time()
for _ in 1:20
    L(solver, p, w, c, rp, λ_ψ) # => 25.596248877655327 for this weight
end
t_elapsed = time() - start
println(solver, " | time(s)= ", t_elapsed)

My configuration is the following:

  • macBook pro M2
  • julia 1.10.5

and packages are:

  • JuMP v1.26.0
  • GLPK v1.2.1
  • HiGHS v1.18.1
  • Gurobi v1.7.4
  • CPLEX v1.1.1

Why can’t rp - z be modeled as a fixed expression?
Why is your λ_ψ unaltered in the 1000 times?
Is α >= 0 necessary in a Max context? Can you drop it?


You can use anonymous constraint syntax, in which case you obviate the need to unregister.
You can use assert_is_solved_and_feasible.
You can leverage import LinearAlgebra.⋅ as ⋅ to shorten you expression, e.g.
JuMP.@constraint(model, w ⋅ x ≤ c)
JuMP.@expression(model, z[k = 1:d], view(p, k, :) ⋅ x)

I think this is just solver variability:

julia> using JuMP

julia> import CPLEX

julia> import GLPK

julia> import HiGHS

julia> import Gurobi

julia> function L(
           solver,
           p::Matrix{Int},
           w::Vector{Int},
           c::Int,
           rp::Vector{Int},
           λ_ψ::Vector{Vector{Float64}},
       )
           obj_vals = Float64[]
           d, n = size(p)
           model = Model(solver)
           set_silent(model)
           @variable(model, x[1:n], Bin)
           @constraint(model, w' * x <= c)
           @expression(model, z, p * x)
           @variable(model, α >= 0)
           @objective(model, Max, α)
           solve_time_sec = 0.0
           for i in 1:length(λ_ψ)
               con = @constraint(model, [k in 1:d], α <= - λ_ψ[i][k] * (rp[k] -  z[k]))
               optimize!(model)
               assert_is_solved_and_feasible(model)
               push!(obj_vals, objective_value(model))
               solve_time_sec += solve_time(model)
               delete.(model, con)
           end
           return obj_vals, solve_time_sec
       end
L (generic function with 1 method)

julia> function main(solver, n, rep)
           p = [
               13 10  3 16 12 11  1  9 19 13
               1 10  3 13 12 19 16 13 11  9
           ]
           w = [4, 4, 3, 5, 5, 3, 2, 3, 5, 4]
           c = 19
           rp = [40, 40]
           λ_ψ = fill([1.0238499551062132, 4.66018735258353], n)
           start = time()
           solve_time_sec = 0.0
           obj_vals = Float64[]
           for _ in 1:rep
               obj_val, t = L(solver, p, w, c, rp, λ_ψ)
               solve_time_sec += t
               append!(obj_vals, obj_val)
           end
           t_elapsed = time() - start
           println(solver, " | time(s)      = ", t_elapsed)
           println(solver, " | solve_time(s)= ", solve_time_sec)
           println(solver, " | overhead(s)  = ", t_elapsed - solve_time_sec)
           return sum(obj_vals)
       end
main (generic function with 2 methods)

julia> obj = main(GLPK.Optimizer, 100, 2)
GLPK.Optimizer | time(s)      = 0.03159523010253906
GLPK.Optimizer | solve_time(s)= 0.019971847534179688
GLPK.Optimizer | overhead(s)  = 0.011623382568359375
5119.249775531065

julia> obj = main(HiGHS.Optimizer, 100, 2)
HiGHS.Optimizer | time(s)      = 7.252434968948364
HiGHS.Optimizer | solve_time(s)= 7.196737084072083
HiGHS.Optimizer | overhead(s)  = 0.05569788487628102
5119.249775531065

julia> obj = main(Gurobi.Optimizer, 100, 2)
Set parameter LicenseID to value 890341
Set parameter LicenseID to value 890341
Gurobi.Optimizer | time(s)      = 1.224890947341919
Gurobi.Optimizer | solve_time(s)= 1.1032283306121826
Gurobi.Optimizer | overhead(s)  = 0.12166261672973633
5119.249775531065

julia> obj = main(CPLEX.Optimizer, 100, 2)
CPLEX.Optimizer | time(s)      = 4.104698896408081
CPLEX.Optimizer | solve_time(s)= 4.063469409942627
CPLEX.Optimizer | overhead(s)  = 0.0412294864654541
5119.249775531065

If I had to guess, it’s that GLPK just cracks on and solves the problem naively. Whereas the others are typically engineered for much larger problems, so they spend more time doing stuff that works well on large problems, but isn’t required for a model like this.

If you change to the “modify constraint coefficient” approach: no change in GLPK, HiGHS is much faster, Gurobi is faster, CPLEX is slower.

julia> using JuMP

julia> import CPLEX

julia> import GLPK

julia> import HiGHS

julia> import Gurobi

julia> function L(
           solver,
           p::Matrix{Int},
           w::Vector{Int},
           c::Int,
           rp::Vector{Int},
           λ_ψ::Vector{Vector{Float64}},
       )
           obj_vals = Float64[]
           d, n = size(p)
           model = Model(solver)
           set_silent(model)
           @variable(model, x[1:n], Bin)
           @variable(model, rhs[1:d])
           @variable(model, α)
           @constraint(model, w' * x <= c)
           @constraint(model, rhs .== rp .- p * x)
           @objective(model, Max, 1.0 * α)
           @constraint(model, con[k in 1:d], α + 1.0 * rhs[k] <= 0)
           solve_time_sec = 0.0
           for λi in λ_ψ
               set_normalized_coefficient.(con, rhs, λi)
               optimize!(model)
               assert_is_solved_and_feasible(model)
               push!(obj_vals, objective_value(model))
               solve_time_sec += solve_time(model)
           end
           return obj_vals, solve_time_sec
       end
L (generic function with 1 method)

julia> function main(solver, n, rep)
           p = [
               13 10  3 16 12 11  1  9 19 13
               1 10  3 13 12 19 16 13 11  9
           ]
           w = [4, 4, 3, 5, 5, 3, 2, 3, 5, 4]
           c = 19
           rp = [40, 40]
           λ_ψ = fill([1.0238499551062132, 4.66018735258353], n)
           start = time()
           solve_time_sec = 0.0
           obj_vals = Float64[]
           for _ in 1:rep
               obj_val, t = L(solver, p, w, c, rp, λ_ψ)
               solve_time_sec += t
               append!(obj_vals, obj_val)
           end
           t_elapsed = time() - start
           println(solver, " | time(s)      = ", t_elapsed)
           println(solver, " | solve_time(s)= ", solve_time_sec)
           println(solver, " | overhead(s)  = ", t_elapsed - solve_time_sec)
           return sum(obj_vals)
       end
main (generic function with 2 methods)

julia> obj = main(GLPK.Optimizer, 100, 2)
GLPK.Optimizer | time(s)      = 0.10955095291137695
GLPK.Optimizer | solve_time(s)= 0.04365062713623047
GLPK.Optimizer | overhead(s)  = 0.06590032577514648
5119.249775531065

julia> obj = main(HiGHS.Optimizer, 100, 2)
HiGHS.Optimizer | time(s)      = 1.5504350662231445
HiGHS.Optimizer | solve_time(s)= 1.4861632268875837
HiGHS.Optimizer | overhead(s)  = 0.0642718393355608
5119.249775531056

julia> obj = main(Gurobi.Optimizer, 100, 2)
Set parameter LicenseID to value 890341
Set parameter LicenseID to value 890341
Gurobi.Optimizer | time(s)      = 0.6475129127502441
Gurobi.Optimizer | solve_time(s)= 0.5366668701171875
Gurobi.Optimizer | overhead(s)  = 0.11084604263305664
5119.249775531065

julia> obj = main(CPLEX.Optimizer, 100, 2)
CPLEX.Optimizer | time(s)      = 8.036567211151123
CPLEX.Optimizer | solve_time(s)= 7.830169677734375
CPLEX.Optimizer | overhead(s)  = 0.20639753341674805
5119.249775531065

julia> obj = main(GLPK.Optimizer, 100, 2)
GLPK.Optimizer | time(s)      = 0.0437159538269043
GLPK.Optimizer | solve_time(s)= 0.03648018836975098
GLPK.Optimizer | overhead(s)  = 0.00723576545715332
5119.249775531065

julia> obj = main(HiGHS.Optimizer, 100, 2)
HiGHS.Optimizer | time(s)      = 1.1625919342041016
HiGHS.Optimizer | solve_time(s)= 1.1198832280933857
HiGHS.Optimizer | overhead(s)  = 0.042708706110715866
5119.249775531056

julia> obj = main(Gurobi.Optimizer, 100, 2)
Set parameter LicenseID to value 890341
Set parameter LicenseID to value 890341
Gurobi.Optimizer | time(s)      = 0.3628709316253662
Gurobi.Optimizer | solve_time(s)= 0.30057835578918457
Gurobi.Optimizer | overhead(s)  = 0.06229257583618164
5119.249775531065

julia> obj = main(CPLEX.Optimizer, 100, 2)
CPLEX.Optimizer | time(s)      = 9.893035888671875
CPLEX.Optimizer | solve_time(s)= 9.836957931518555
CPLEX.Optimizer | overhead(s)  = 0.05607795715332031
5119.249775531065

Thank you very much for your feedbacks, and your suggestions to improve the code.

To clarify two points noticed in your comments, the weights are computed to estimate a particular measure for an optimization problem with d objectives. The number of weights is supposed to be very large (e.g. 10 000). Also, I provided this snipped of code using the knapsack as decision space X. But X may be any IP (and indeed, GLPK will probably struggles faced to a less basic problem), explaining this test with different solvers.

1 Like

Ill flag this up with the HiGHS folks next week. They may have some insight. Ill also see if I can get a better profile.

1 Like