Create JuMP model with multi-threading in Julia

Question

I have an optimization model which turns to be very difficult to build. This model has many if-else conditions and many loops as well. So I was thinking of using multi-threading for building this single JuMP model object.

A very simplified version of one loop of the code looks like this:

Threads.@threads for g in sets["A"]

    Array_1 = [gg for gg in [sets["B"];sets["A"]] if data2[gg] == g]
    Array_2 = [gg for gg in sets["B"] if data[gg] == g]

    for t in STAGES
        Array_3 = [gg for gg in [sets["B"];sets["A"]] if data2[gg] == g && (gg, t) in sets["C"] ]
        for b in BLOCKS
            name = @constraint( model, ((g, t, b) in sets["C"] ? X1[(g,t,b)] : 0)
            - sum(X1[(gg,t,b)] for gg in Array_3 )
            + X2[(g,t,b)] - sum(X2[(gg,t,b)] for gg in Array_1)
            - sum(data3[gg] for gg in Array_2) == data4[(g, t, b)])
        end
    end

    a=string("con_",g,"_",t,"_",b)
    JuMP.set_name(name,a)
end

I have many of those loops with many if-else conditions inside. So I added @Threads.threads before the first for g in sets["A"] aiming to reduce the time of building the model.

The problem is that I obtain an ERROR: LoadError: TaskFailedException: UndefRefError: access to undefined reference when renaming the constraint. Is there any problem about my approach? If I don't put the Threads.@threads there isn't any problem at all, it just works very slow.

Some information about the infrastructure:

julia> versioninfo()
Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2660 v3 @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, haswell)
Environment:
  JULIA_NUM_THREADS = 40

and packages:

(@v1.4) pkg> status
Status `~/.julia/environments/v1.4/Project.toml`
  [c7e460c6] ArgParse v1.1.0
  [a076750e] CPLEX v0.6.6
  [336ed68f] CSV v0.7.7
  [e2554f3b] Clp v0.8.1
  [a93c6f00] DataFrames v0.21.7
  [5789e2e9] FileIO v1.4.3
  [2e9cd046] Gurobi v0.8.1
  [033835bb] JLD2 v0.2.1
  [4076af6c] JuMP v0.21.5
  [438e738f] PyCall v1.91.4
  [2913bbd2] StatsBase v0.33.1
  [bd369af6] Tables v1.0.5
  [6dd1b50a] Tulip v0.6.2
  [1a1011a3] SharedArrays
  [10745b16] Statistics

Thanks in advance!

Full stacktrace:

ERROR: LoadError: TaskFailedException:
UndefRefError: access to undefined reference
Stacktrace:
 [1] getindex at ./array.jl:788 [inlined]
 [2] ht_keyindex2!(::Dict{MathOptInterface.ConstraintIndex,String}, ::MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}}) at ./dict.jl:326
 [3] setindex!(::Dict{MathOptInterface.ConstraintIndex,String}, ::String, ::MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}}) at ./dict.jl:381
 [4] set at /home/user/.julia/packages/MathOptInterface/k7UUH/src/Utilities/model.jl:349 [inlined]
 [5] set at /home/user/.julia/packages/MathOptInterface/k7UUH/src/Utilities/universalfallback.jl:354 [inlined]
 [6] set(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, ::MathOptInterface.ConstraintName, ::MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}}, ::String) at /home/user/.julia/packages/MathOptInterface/k7UUH/src/Utilities/cachingoptimizer.jl:646
 [7] set(::Model, ::MathOptInterface.ConstraintName, ::ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},ScalarShape}, ::String) at /home/user/.julia/packages/JuMP/qhoVb/src/JuMP.jl:903
 [8] set_name(::ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},ScalarShape}, ::String) at /home/user/.julia/packages/JuMP/qhoVb/src/constraints.jl:68
 [9] macro expansion at /home/user/code/model_formulation.jl:117 [inlined]
 [10] (::var"#20#threadsfor_fun#255"{Dict{Any,Any},Dict{Any,Any},JuMP.Containers.DenseAxisArray{VariableRef,1,Tuple{Array{Tuple{String,Int64,Int64},1}},Tuple{Dict{Tuple{String,Int64,Int64},Int64}}},JuMP.Containers.DenseAxisArray{VariableRef,1,Tuple{Array{Tuple{String,Int64,Int64},1}},Tuple{Dict{Tuple{String,Int64,Int64},Int64}}},Array{String,1}})(::Bool) at ./threadingconstructs.jl:61
 [11] (::var"#20#threadsfor_fun#255"{Dict{Any,Any},Dict{Any,Any},JuMP.Containers.DenseAxisArray{VariableRef,1,Tuple{Array{Tuple{String,Int64,Int64},1}},Tuple{Dict{Tuple{String,Int64,Int64},Int64}}},JuMP.Containers.DenseAxisArray{VariableRef,1,Tuple{Array{Tuple{String,Int64,Int64},1}},Tuple{Dict{Tuple{String,Int64,Int64},Int64}}},Array{String,1}})() at ./threadingconstructs.jl:28
Stacktrace:
 [1] wait(::Task) at ./task.jl:267
 [2] macro expansion at ./threadingconstructs.jl:69 [inlined]
 [3] model_formulation(::Dict{Any,Any}, ::Dict{Any,Any}, ::Dict{Any,Any}, ::Dict{String,Bool}, ::String) at /home/user/code/model_formulation.jl:102
 [4] functionA(::Dict{Any,Any}, ::Dict{Any,Any}, ::Dict{Any,Any}, ::String, ::Dict{String,Bool}) at /home/user/code/functionA.jl:178
 [5] top-level scope at /home/user/code/main.jl:81
 [6] include(::Module, ::String) at ./Base.jl:377
 [7] exec_options(::Base.JLOptions) at ./client.jl:288
 [8] _start() at ./client.jl:484
in expression starting at /home/user/code/main.jl:81

Answer 1

You have two options to parallelize a JuMP optimization model

1. Run a multi-threaded version of the Solver (provided that the solver supports it) - in that case the parallelism is fully handled by the external solver library and your Julia process remains single-threaded.

2. Run several single-threaded solver processes in parallel threads controlled by Julia. In this case several copies of the model need to be separately created which you can try to send to the solver at the same time.

#1:

Solvers support parameters including multi-threading control (on the other hand they might be simply using all available threads by default). Here is an example with Gurobi:

using JuMP, Gurobi
m = Model(optimizer_with_attributes(Gurobi.Optimizer,  "Threads" => 2))
@variable(m, 0 <= x <= 2)
@variable(m, 0 <= y <= 30)
@objective(m, Max, 5x + 3 * y)
@constraint(m, con, 1x + 5y <= 3)
optimize!(m)  # the model will be optimized using 2 threads

#2:

Running many solver copies in parallel you need to have separate model copies. In my code they differ by the range for x parameter:

Threads.@threads for z in 1:4
    m = Model(optimizer_with_attributes(Gurobi.Optimizer,  "Threads" => 1))
    @variable(m, 0 <= x <= z)
    @variable(m, 0 <= y <= 30)
    @objective(m, Max, 5x + 3 * y)
    @constraint(m, con, 1x + 5y <= 3)
    optimize!(m) 
    #todo collect results
end

These are two separate approaches and you cannot mix them. If you parallelize execution each thread needs to get a separate model copy because JuMP mutates the Model object.