Reoptimization with GLPK in Julia/JuMP

Question

I've been looking everywhere, and I can't seem to find any information on whether GLPK can do successive optimizations on the same model, on which we add a constraint at each iteration. More concretely, I want to implement a Branch and Bound algorithm on a 01ILP, so at each iteration, I'll fix a variable either to 0 or 1 (depending on the branch) and solve the relaxed LP, but I don't want to reoptimize from scratch each time.

Do you know if such a thing is possible using GLPK (and Julia/JuMP as modeling language)? And if not, do you see a solver that would do that?

Answer 1

Yes, you can do such a thing in JuMP. If you use the so-called direct mode, then incremental changes to the model are passed directly to the solver. If your solver supports incremental changes and re-solve, then it will do so. For instance, if you are solving an LP, adding a constraint, and re-solving, then the second solve will (likely) be a lot faster than the first.

All mainstream solvers I know of, e.g., CPLEX, Gurobi, GLPK, HiGHS, etc., support this. BTW: nowadays, I would recommend using HiGHS instead of GLPK. It's much much faster and more numerically robust.

Here is a small example based on a knapsack problem

using JuMP
using GLPK

# Generate instance data
n = 8        # number of items
w = rand(n)  # item weights
c = rand(n)  # item values
W = sum(w) / 2  # total knapsack capacity

# Build the optimization model
# (we only give the LP relaxation to GLPK)
model = JuMP.Model(GLPK.Optimizer)

# Create x variables in [0, 1]
# they are not marked as binary to do branching manually later
@variable(model, x[1:n])
@constraint(model, bounds, 0 .<= x .<= 1)

@constraint(model, capacity, w'x <= W)  # knapsack capacity
@objective(model, Max, c'x)  # maximize value of selected items

# Main loop
optimize!(model)
z_root = objective_value(model)
x_ = value.(x)

# branch on first fractional variable
# (here we assume that `x1` is fractional)
# Left (x1 = 0) node
JuMP.fix(x[1], 0.0)
optimize!(model)
z0 = objective_value(model)
# Right (x1 = 1) node
JuMP.unfix(x[1])
JuMP.fix(x[1], 1.0)
optimize!(model)
z1 = objective_value(model)

If you have more technical questions, I encourage you to post them on the Julia discourse forum on mathematical optimization. You'll get more traction there.