Crie um modelo JuMP com multi-threading em Julia

Question

Tenho um modelo de otimização que se torna muito difícil de construir. Este modelo tem muitas condições if-else e também muitos loops. Então, eu estava pensando em usar multi-threading para construir este único objeto de modelo JuMP.

Uma versão muito simplificada de um loop do código tem a seguinte aparência:

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

Eu tenho muitos desses loops com muitas condições if-else dentro. Então acrescentei @Threads.threadsantes do primeiro for g in sets["A"]com o objetivo de reduzir o tempo de construção do modelo.

O problema é que obtenho um ERROR: LoadError: TaskFailedException: UndefRefError: access to undefined referenceao renomear a restrição. Existe algum problema com a minha abordagem? Se eu não colocar não Threads.@threadshá problema nenhum, só funciona muito devagar.

Algumas informações sobre a infraestrutura:

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

e pacotes:

(@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

Desde já, obrigado!

Stacktrace completo:

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

Você tem duas opções para paralelizar um modelo de otimização JuMP

1. Execute uma versão multiencadeada do Solver (desde que seja compatível com o solucionador) - nesse caso, o paralelismo é totalmente controlado pela biblioteca externa do solver e seu processo Julia permanece com uma única thread.

2. Execute vários processos de solucionador de thread único em threads paralelos controlados por Julia. Nesse caso, várias cópias do modelo precisam ser criadas separadamente, que você pode tentar enviar ao solucionador ao mesmo tempo.

# 1:

Os solucionadores suportam parâmetros, incluindo controle de multi-threading (por outro lado, eles podem simplesmente usar todos os threads disponíveis por padrão). Aqui está um exemplo com 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:

Executando muitas cópias do solver em paralelo, você precisa ter cópias separadas do modelo. No meu código, eles diferem pelo intervalo do xparâmetro:

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

Essas são duas abordagens separadas e você não pode misturá-las. Se você paralelizar a execução, cada thread precisa obter uma cópia do modelo separada porque o JuMP modifica o Modelobjeto.