# 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
set_optimizer_attribute(model,"Method",2)
# No cross over
set_optimizer_attribute(model,"Crossover",0)

# Define dual bound threshold
thr = 1e3

# Define callback
c = Float64[]
function absolute_stopping_criterion(cb_data,cb_where)
objbound =Ref{Cdouble}()
Gurobi.GRBcbget(cb_data,cb_where,Gurobi.GRB_CB_BARRIER_DUALOBJ,objbound)
objbound = objbound[]
push!(c,objbound)
if maximum(c) > thr
GRBterminate(backend(model).inner)
end
return obj
end

# Attach callback function to model
MOI.set(model,Gurobi.CallbackFunction(),absolute_stopping_criterion)

# Optimize model
optimize!(model)

# retrieve termination_status
termination_status(model)

# Obtain function value and primal solution
objective_value(model)
value.(model[:x])
``````

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.

Best

Felix

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,

Best

Felix

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

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)")
``````

Gives:

``````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?

Best

Felix

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[])")
end
end
``````

But perhaps someone more knowledgeable can suggest a nicer way.

Hi,

thanks so much! I will give this a try.

Best

Felix

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. https://www.tandfonline.com/doi/abs/10.1080/10556788.2013.871282) 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 https://www.gurobi.com/documentation/current/refman/bestobjstop.html#parameter:BestObjStop for LPs as well.

1 Like