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