How to use the solver-dependent callback

I want to use Gurobi to implement Bender’s decomposition. I write it in Python:
LB = model.cbGet(GRB.Callback.MIPSOL_OBJBND)

How to call the cbGet(GRB.Callback.MIPSOL_OBJBND) in Julia? Thank you so much.

1 Like

Hi there,

Since this is your first post, take a read of Please read: make it easier to help you. In general, it’s easier to help if you add more content to your question. You can also put JuMP-related questions in the “Optimization (Mathematical)” section.

You need to use the C API for Gurobi: GitHub - jump-dev/Gurobi.jl: Julia interface for Gurobi Optimizer

So something like this should get you started

lb_pointer = Ref{Cdouble}(0.0)
GRBcbget(cb_data, ch_where, CRB_MIPSOL_OBJBND, lb_pointer)
lb = lp_pointer[]

There’s a JuMP tutorial for Benders decomposition with the solver-independent API, but it doesn’t use the objective bound:
https://jump.dev/JuMP.jl/stable/tutorials/algorithms/benders_decomposition/#Callback-method

Hi, thank you for you timely reply. I have viewed the page and know how to post a question more efficiently.

I have another question: how to retrieve the UnbdRay in Julia. Indeed, there is a quite similar question but no solution:

Thank you so much!!!

The easiest way is GitHub - jump-dev/Gurobi.jl: Julia interface for Gurobi Optimizer.

Otherwise, Gurobi.jl wraps the complete C API, so follow the examples here: https://www.gurobi.com/documentation/9.5/refman/c_attribute_examples.html

Using the C API is quite an advanced feature, so it takes some getting used to. You also have to account for the 1-based indexing in Julia and the 0-based indexing in C.

For example:

 # Get the full ray as a vector
num_cols = ... 
x = fill(0.0, num_cols) 
GRBgetdblattrarray(m, "UnbndRay", 0, num_cols, x)

or

 # Get the coefficient of the ray for the 5th variable
valueP = Ref{Cdouble}()
GRBgetdblattrelement(m, "UnbndRay", 4, valueP)  # Note that the element is index-1
ray = valueP[]

Also, if you solve an unbounded LP, the ray is available via value(x). Check primal_status(model) to see if it is INFEASIBILITY_CERTIFICATE.

Hi, Julia is recommended due to its high speed relative to Python. However, when I use the tsp.jl and tsp.py to solve the tsp problem with sub-tour elimination callback function under the same instances and the same solver(Gurobi). How should I explain it?

*** Python is always faster than Julia when node=10, 50, 100…)

Please provide the code. If you’re using JuMP and gurobipy, then python will likely be faster.

The python code: https://www.gurobi.com/documentation/9.5/examples/tsp_py.html
The Julia code: JuMP.jl/tsp_lazy_constraints.jl at master · jump-dev/JuMP.jl · GitHub
The computation instance: 100 customers
I always use the same IDE and same slover(Gurobi)。

Julia is recommended due to its high speed relative to Python

In these cases, if you only want to use Gurobi, there is little benefit to using Julia. Most of the time is spent in Gurobi. In addition, JuMP is solver-independent, so it has some overhead compared to a tool like gurobipy, although you can reduce this with some effort.

Also, how are you timing things?

1 Like

Hi, Mr. odow. To testify whether the code used Julia would run faster then that with python, I translate my code to Julia. However there is something wrong during the callback process. The question is that: I don’t know why the callback doesn’t end by my condition.

function BendersCut(cb_data, cb_where)
status = callback_node_status(cb_data, m_master)
if status != MOI.CALLBACK_NODE_STATUS_INTEGER
return # Only run at integer solutions
end

Gurobi.load_callback_variable_primal(cb_data, cb_where)
r_value = callback_value.(Ref(cb_data), r_ij)
r_value = convert(Matrix{Float64}, r_value)


z_value = callback_value.(Ref(cb_data), z_ijl)
z_value = convert(Array{Float64, 3}, z_value)

# Get the current best objbound MIPSOL_OBJBND
LB_pointer = Ref{Cdouble}()  
Gurobi.GRBcbget(cb_data, cb_where, Gurobi.GRB_CB_MIPSOL_OBJBND, LB_pointer)
LB = LB_pointer[]

# Set the objective function of m_sub(up to z_value)
@objective(m_sub, Max, sum(w_3[i, j] * path_y[i, j] for i in 1:N for j in B)
    - sum(w_4[i, j] * path_y[i, j] for i in 1:N for j in B)
    + sum(w_6[i] * path_D[i, 1] for i in 1:N)
    - sum(w_7[i] * path_D[i, 2] for i in 1:N)
    - K * sum(z_value[i+1, j+1, l+1] * u[i, j, l] for l in L for j in 0:N for i in 0:N)
    - K * sum(w_8[i] for i in 1:N))
optimize!(m_sub)
m_sub_status = primal_status(m_sub)
w_1_value = JuMP.value.(w_1)
w_2_value = JuMP.value.(w_2)
w_3_value = JuMP.value.(w_3)
w_4_value = JuMP.value.(w_4)
w_5_value = JuMP.value.(w_5)
w_8_value = JuMP.value.(w_8)
u_value = JuMP.value.(u)

if m_sub_status == NO_SOLUTION
    println("This is something wrong!!!")
    return

else
    lin_exp_2 = K * sum(
        u_value[i, j, l] * m_master[:z_ijl][i, j, l] for l in L
        for j in 0:N for i in 0:N)

    lin_exp = sum(
    w_3_value[i, j] * path_y[i, j] for j in B for i in 1:N) - sum(
    w_4_value[i, j] * path_y[i, j] for j in B for i in 1:N) - lin_exp_2 - K * sum(w_8_value[i] for i in 1:N)
    
        
    if m_sub_status == INFEASIBILITY_CERTIFICATE
        println("add the feasibility constraints")
        # Add the feasibility cut(Unbounded)
        UB = -float(Inf)
        # Add the lazy constraint
        con = @build_constraint(0 >= lin_exp)
        MOI.submit(m_master, MOI.LazyConstraint(cb_data), con)

    elseif m_sub_status == FEASIBLE_POINT
        println("add the optimality constraints")
        # Add the optimality cut
        UB = objective_value(m_sub)
        outer_value = sum(path_t[i, j] * r_value[i, j] for i in 1:N+1 for j in 1:N+1) + sum(l * z_value[i, j, l+1] for l in L for j in 1:N+1 for i in 1:N+1)
        UB += outer_value

        println("UB!!!!")
        println(UB)
        println(LB)


        if abs(UB - LB)<0.00001
            Gurobi.GRBterminate(backend(m_master).inner)
        else

            println("Continue")
                

            # Add the lazy constraint
            con = @build_constraint(m_master[:q] >= lin_exp)
            MOI.submit(m_master, MOI.LazyConstraint(cb_data), con)
            println("Continue")
        end
    end
end

end

Indeed, I have know the optimal value before hand(which is 16.6) and I have tested that the cuts added are all correct.

But when the callback function are called twice, the process ends.

Root relaxation: objective 4.000000e+00, 8 iterations, 0.00 seconds (0.00 work units)
Set parameter InfUnbdInfo to value 1
add the optimality constraints
UB!!!
(UB)16.6
(LB) 4.0
Continue
Continue
add the optimality constraints
UB!!!
(UB)16.6
(LB)10.600000000000001
Continue
Continue

Nodes    |    Current Node    |     Objective Bounds      |     Work

Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time

  • 0 0 0 10.6000000 10.60000 0.00% - 0s

Cutting planes:
Lazy constraints: 1

From the result we could see that the UB(the dual problem’s value) is right and the LB (the master problem’s value) is increasing. But why it would end at 10.6?