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

Specific Domains Optimization (Mathematical)
,
1

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

2

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.

3

Hi Oscar,

Thank you very much for your answer!

Best

Felix

4

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

5

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

1 Like
6

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
7

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

8

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.

9

Hi,

thanks so much! I will give this a try.

Best

Felix

10

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
11

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

12

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.

13

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.

14

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