Invalid Gurobi Licence when instance size increases using lazy constraints to generate cuts

Specific Domains Optimization (Mathematical)
1

Hello,
I’m using CPLEX and Gurobi to solve a benders decomposition scheme for the mining problem.
I defined the master problem using CPLEX solver and the cut generation subproblems using Gurobi.
Whenever I increase the instance size too much (number of periods or scenarios) I get the following error (which is reallly weird because LinearQuadraticModel function is in the log of the error, which I’m sure is linear):

Invalid Gurobi license
Gurobi.Env() at grb_env.jl:13
#GurobiMathProgModel#10(::Array{Any,1}, ::Type{T} where T, ::Void) at GurobiSolverInterface.jl:24
(::Core.#kw#Type)(::Array{Any,1}, ::Type{Gurobi.GurobiMathProgModel}, ::Void) at <missing>:0
LinearQuadraticModel(::Gurobi.GurobiSolver) at GurobiSolverInterface.jl:68
#build#119(::Bool, ::Bool, ::JuMP.ProblemTraits, ::Function, ::JuMP.Model) at solvers.jl:356
(::JuMP.#kw##build)(::Array{Any,1}, ::JuMP.#build, ::JuMP.Model) at <missing>:0
#solve#116(::Bool, ::Bool, ::Bool, ::Array{Any,1}, ::Function, ::JuMP.Model) at solvers.jl:168
solveDMPbend(::Array{Any,1}, ::Bool) at adp.jl:221
SolveADP(::Int64, ::Array{Int64,2}, ::Float64, ::Array{Float64,1}, ::Array{Int32,1}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,2}, ::UnitRange{Int64}, ::Int64, ::Int64, ::Array{Any,1}, ::Array{Float64,2}, ::UnitRange{Int64}, ::Array{Int32,1}, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Any,1}, ::Array{Int64,2}, ::Array{Array{Int64,2},1}, ::Array{Int64,2}, ::Array{Any,1}, ::Int64, ::Int64, ::Int64, ::Float64, ::Bool, ::Float64, ::Float64, ::Bool, ::Bool) at adp.jl:505
Simulate(::Float64, ::Array{Float64,1}, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Float64, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Array{Array{Int64,2},1}, ::Array{Float64,1}, ::Array{Int64,2}, ::Array{Int64,2}, ::Array{Int64,2}, ::Array{Any,1}, ::Bool, ::Int64, ::Bool, ::Bool, ::Bool) at sim.jl:242
macro expansion at sim.jl:677 [inlined]
anonymous at <missing>:?
include_string(::String, ::String) at loading.jl:515
include_string(::String, ::String, ::Int64) at eval.jl:30
include_string(::Module, ::String, ::String, ::Int64, ::Vararg{Int64,N} where N) at eval.jl:34
(::Atom.##100#105{String,Int64,String})() at eval.jl:75
withpath(::Atom.##100#105{String,Int64,String}, ::String) at utils.jl:30
withpath(::Function, ::String) at eval.jl:38
hideprompt(::Atom.##99#104{String,Int64,String}) at repl.jl:61
macro expansion at eval.jl:73 [inlined]
(::Atom.##98#103{Dict{String,Any}})() at task.jl:80

Here the implementation of the callback and the model:

@everywhere function definemaster(ʖ,ℓ,T,τ,Edges,Ωext,Captotalextraccion)
  m = Model(solver = mastersolver);
  @variable(m, x[j=1:length(ʖ),t=1:(T-τ+1)] ,Bin);
  @variable(m, z)
  @objective(m,:Max,z);
  @constraint(m,
    ext_once_every_cluster[i=1:length(ʖ)] ,
    sum(x[i,ϕ] for ϕ=1:(T-τ+1)) <= 1
  );

  # print("Precedence constraints\n");
  @constraint(m,
    precedences[ee=1:length(Edges[1]), t=1:(T-τ+1)],
    sum(x[ℓ[Edges[1][ee]],s] for s=1:t) >= sum(x[ℓ[Edges[2][ee]],s] for s=1:t)
  );

  # print("KP constraints\n")
  @constraint(m,
    KP[t=1:(T-τ+1)],
    sum(x[i,t]*Ωext[ʖ[i]] for i=1:length(ʖ)) <= Captotalextraccion
  );

  return m,x,z,ext_once_every_cluster,precedences,KP

end


@everywhere function solvesubproblem(sol,ω,blockcluster,blocks,ℓ,weights,Captotalproceso1)
  sb=Model(solver=subproblemsolver);
  @variable(sb,μ>=0);
  @variable(sb,λ[b in blocks]);
  @objective(sb,:Min,
    Captotalproceso1*μ +
    sum(λ[b]* (max(sol[ℓ[blockcluster[b]]],0)) for b in blocks )
  );
  @constraint(sb,
    diet1[b in blocks],
    λ[b] + weights[b]*μ  >= ω[b,2]
  );
  @constraint(sb,
    diet2[b in blocks],
    λ[b] >= ω[b,1]
  );
  sta=solve(sb);
  if sta==:Unbounded
    print(sb)
    writeLP(sb,string(wdir,"out"))
  end
  return getobjectivevalue(sb),getvalue(λ),getvalue(μ)
end




@everywhere function solveDMPbend(arg,givesolution=false)
  Captotalextraccion,Captotalproceso1,Ωext,ν,Ω,edges,blockcluster,bci,clusterblocks,wx,Ξ,τ,T,D,r,
  angulocut,U=arg
  scen,lawscb=ν
  NoEscenarios=length(Ξ)

  scenarioslaws=linearTprc(scen,lawscb,Ω,blockcluster,NoEscenarios,bci)
  𝝹=length(clusterblocks)

  Edges = [Int64[],Int64[]]
  if τ==1
    Edges=edges
  else
    for e in 1:length(edges[1])
      bol=true
      if sum(wx[edges[1][e],ϕ] for ϕ = 1:(τ-1)) > 0 || sum(wx[edges[2][e],ϕ] for ϕ = 1:(τ-1)) > 0
        bol=false
      end
      if bol
        push!(Edges[1],edges[1][e])
        push!(Edges[2],edges[2][e])
      end
    end
  end

  ʖ = Int64[]
  ℓ = Int64[]
  blocks=Int32[]
  j=1
  for i=1:𝝹
    if τ==1 || sum(wx[i,1:(τ-1)])<0.99
      push!(ʖ,i)
      push!(ℓ,j)
      blocks=union(blocks,clusterblocks[i])
      j=j+1
    else
      push!(ℓ,0)
    end
  end


  if length(ʖ)==0
    if givesolution
      return zeros(𝝹,T-τ+1)
    else
      return 0.0
    end
  end

  blocks = setdiff(bci,setdiff(bci,blocks))

  m,x,z,ext_once_every_cluster,precedences,KP=definemaster(ʖ,ℓ,T,τ,Edges,Ωext,Captotalextraccion);



  function cut(cb)
    masterval=getvalue(z);
    sol=getvalue(x)
    if angulocut aexpr=zero(AffExpr) end
    bexpr=zero(AffExpr)
    solval=0;
    solτ=[]
    cf=zeros(length(ʖ),T-τ+1)
    for t = 0:(T-τ)
      push!(solτ,Float64[])
      for j in ʖ
        if ℓ[j]>0 push!(solτ[t+1],sol[ℓ[j],t+1]) end
      end
      for ϕ=1:length(Ξ)
        ξ=Ξ[ϕ]
        sub,λ,μ = solvesubproblem(solτ[t+1], scenarioslaws[ : , (D*ξ - 1):(D*ξ) ] * (1/(1+r))^t /length(Ξ) ,blockcluster,blocks,ℓ,Ω,Captotalproceso1)
        solval=solval + sub  ;
        bexpr +=  μ*Captotalproceso1
        for c =1:length(ʖ)
          for b in clusterblocks[ʖ[c]]
            cf[c,t+1]+=λ[b]
          end
        end
      end

      for c =1:length(ʖ)
        push!(bexpr,cf[c,t+1],x[c,t+1])
      end
   
    if masterval-0.0001<=solval<=masterval+0.0001
      return
    end

    # println("Adding sample average scenarios elimination cut")
    # println("----")
    if solval!=-Inf
      #bexpr=bexpr-z
      push!(bexpr,-1.0,z)
    end
    @lazyconstraint(cb, bexpr >= 0)
  end
  addlazycallback(m, cut)

  TT = STDOUT # save original STDOUT stream
  redirect_stdout()#"ttttttttt.txt")
  solve(m)
  redirect_stdout(TT) # restore STDOUT

  if givesolution
    return getvalue(x)
  else
    return getobjectivevalue(m)* (1/(1+r))^(τ-1)
  end

end
2

Without a proper MWE, it’s hard to know what is going on. However from the sounds of it, you’re running into an issue running to many concurrent Gurobi models for your license.

This part of the README explains what to do: GitHub - jump-dev/Gurobi.jl: Julia interface for Gurobi Optimizer

which is reallly weird because LinearQuadraticModel function is in the log of the error, which I’m sure is linear

LinearQuadraticModel is used to store pure linear models as well, so that’s not the problem.

3

Hello,
Thank you very much.

From your idea I see that the way to solve the problem is to define and solve the subproblems (the ones in Gurobi) within a function, in order to get the solutions values without the model (so this last is dropped from memory).
I’ve already done that:

@everywhere function solvesubproblem(sol,ω,blockcluster,blocks,ℓ,weights,Captotalproceso1)
sb=Model(solver=subproblemsolver);
@variable(sb,μ>=0);
@variable(sb,λ[b in blocks]);
@objective(sb,:Min,
Captotalproceso1μ +
sum(λ[b] (max(sol[ℓ[blockcluster[b]]],0)) for b in blocks )
);
@constraint(sb,
diet1[b in blocks],
λ[b] + weights[b]*μ >= ω[b,2]
);
@constraint(sb,
diet2[b in blocks],
λ[b] >= ω[b,1]
);
sta=solve(sb);
if sta==:Unbounded
print(sb)
writeLP(sb,string(wdir,“out”))
end
return getobjectivevalue(sb),getvalue(λ),getvalue(μ)
end

Note that the function is defined using @everywhere, but any way is solved in serial.

Thanks again.

4

Thank you!
It worked!