Solver-specific callback for custom termination rule for LP Barrier using Gurobi and JuMP

Hi all,

I am implementing a Benders-like bundle method algorithm to solve a large-scale stochastic linear program. I want to use on-demand inexact cuts, i.e. I would like to solve the subproblems accurately only if they meet a certain criterion. To that end, I would like to provide the subproblem minimisation a threshold value corresponding to the upper bound of the overall algorithm. Solving the subproblem with Barrier, I can trace the dual bound and terminate the run early if it surpasses the threshold value for then I know that the subproblem will not yield a better result than the current upper bound.

I have implemented this using a Gurobi-specific callback. Please refer to the reproducible example below:

#%% -----------------------------
# Reproducible example custom stopping rule for LP barrier based on dual bound

using JuMP
using Gurobi
using JSON

data = JSON.parse("""
    "plants": {
        "Seattle": {"capacity": 350},
        "San-Diego": {"capacity": 600}
    "markets": {
        "New-York": {"demand": 300},
        "Chicago": {"demand": 300},
        "Topeka": {"demand": 300}
    "distances": {
        "Seattle => New-York": 2.5,
        "Seattle => Chicago": 1.7,
        "Seattle => Topeka": 1.8,
        "San-Diego => New-York": 2.5,
        "San-Diego => Chicago": 1.8,
        "San-Diego => Topeka": 1.4
P = keys(data["plants"])
M = keys(data["markets"])
distance(p::String, m::String) = data["distances"]["$(p) => $(m)"]
model = direct_model(Gurobi.Optimizer())
@variable(model, x[P, M] >= 0)
@constraint(model, [p in P], sum(x[p, :]) <= data["plants"][p]["capacity"])
@constraint(model, [m in M], sum(x[:, m]) >= data["markets"][m]["demand"])
@objective(model, Min, sum(distance(p, m) * x[p, m] for p in P, m in M));

# enforcing Barrier 
# No cross over

# Define dual bound threshold
thr = 1e3

# Define callback
c = Float64[]
function absolute_stopping_criterion(cb_data,cb_where)
    objbound =Ref{Cdouble}()
    objbound = objbound[]
    if maximum(c) > thr
    return obj 

# Attach callback function to model

# Optimize model

# retrieve termination_status

# Obtain function value and primal solution

While the algorithm terminates once the dual bound surpasses the threshold of 1000, as intended, the termination status is INTERRUPTED and I cannot access the (primal) objective value or the variable solutions of the prematurely terminated model. Is there a way around this, e.g. an alternative stopping command that treats the stopping criterion like the Gurobi Parameter BarConvTol. After all, using a IPM, I should be inside the convex hull defined by the problem constraints.

Thank you very much.



cc @vasyfa works for Gurobi support and might be able to explain what’s going on and offer a suggestion.

But I assume that at the point when you terminate, Gurobi does not have a feasible primal-dual solution to return.

Hi Oscar,

Thank you very much for your answer!



You may be interested in GitHub - odow/SDDP.jl: Stochastic Dual Dynamic Programming in Julia,

I am! I have been experimenting with SDDP.jl for another project!

1 Like

Hi! I think the attributes are there just not accessible. e.g. you can run:

ret = GRBgetintattr(unsafe_backend(model), "SolCount", valueP)
println("solcount = $(valuePdouble[]) ret = $(ret)")
ret = GRBgetdblattrelement(unsafe_backend(model), "RC", 0, valuePdouble)
println("RC = $(valuePdouble[]) ret = $(ret)")
ret = GRBgetdblattrelement(unsafe_backend(model), "X", 0, valuePdouble)
println("X = $(valuePdouble[]) ret = $(ret)")


solcount = 0 ret = 0
RC = 0.02203038263388657 ret = 0
X = 140.39787499441337 ret = 0

The issue lies in that the SolCount attribute is 0 which is used in the logic for Gurobi.jl MOI_wrapper. I will open an issue there.

1 Like

Hi David,
thanks for getting back to me! I see, so I reckon there is no way to do this at the moment?



There’s a hack:

# Hack variable values out
for p in P
    for m in M
        col = unsafe_backend(model).variable_info[x[p, m].index].column
        colc = Cint(col - 1)
        ret = GRBgetdblattrelement(unsafe_backend(model), "X", colc, valuePdouble)
        println("x = $(x[p,m]) has value $(valuePdouble[])")

But perhaps someone more knowledgeable can suggest a nicer way.


thanks so much! I will give this a try.



Something to bear in mind (see GitHub) is that the solution might not be primal feasible. You can check some quality attributes such as MaxVio or ComplVio.

1 Like

Thank you very much for looking into this. I will give the solution on Github a shot.

Does terminating early make a big difference to the run-time? If you don’t have a guarantee on the solution, I don’t really see the motivation for doing this.

It can potentially make a big difference, as each subproblem corresponds to solving an economic dispatch problem for a European energy system model with 8,760 time steps. Using on-demand accuracy oracles (e.g. can substantially speed up the algorithm especially if the subproblems are computationally cumbersome and I had hoped that this might be a simple way of implementing it. But I agree, it does not make much sense if the early stopping almost always returns infeasible solutions. I will run some tests.

I guess this is really a feature request for Gurobi to implement for LPs as well.

1 Like