I’m writing the Benders’ decomposition for an MIP and encountering a very weird behavior of Gurobi.
There is an MIP that I need to modify its obj and repeatedly re-solve (Gurobi can (according to its Logging) Loaded MIP start from previous solve with objective). Some obj coefficient might leads to a harder model (surely).
But I discover a model that once become very very slow to resolve. Therefore I set a TimeLimit and intercept its status which at that moment is questionable. Then I write it as an MPS file, and opened a new Julia REPL to directly solve that—it suddenly becomes very fast (at the supposed speed).
The code to generate the questionable model is here
import JuMP.value as ı
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.⋅ as ⋅
import Random
seed = 4599192817953253712 # may need to cherry-pick
Random.seed!(seed)
J = 9; # number of resources
W = rand(9:49, J); # the sole width type of each resource
o = rand(30:0.01:70, J); # occupied cost
u = 7 * rand(0.75:0.01:1.25, J)/30 .* W; # unit price
L = rand(6:17, J); # purchasing num limit
I = 10; # number of customer
w = rand(7:0.7:24, I); # width needed
S = 2; # number of scenes
d = rand(6:35, I, S); # number needed
function Min_s_prim_model(s) # evaluate the 2nd stage exactly ⚠️ Feasible Region is NOT fixed
Min_s = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
JuMP.@variable(Min_s, m2[1:J]) # JuMP.fix.(m2, m) # 🟡 [Jit] Linking
JuMP.@variable(Min_s, 0 ≤ n[1:I, 1:J], Int)
JuMP.@constraint(Min_s, [i = 1:I], sum(view(n, i, :)) ≥ d[i, s])
JuMP.@constraint(Min_s, [j = 1:J], m2[j]W[j] ≥ w ⋅ view(n, :, j))
JuMP.@objective(Min_s, Min, 0/S)
return Min_s
end; function Rhs_model_s(s) # 🟢 Feasible Region is fixed
Rhs = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
JuMP.@variable(Rhs, 0 ≤ b[1:J] ≤ 1, Int);
JuMP.@variable(Rhs, 0 ≤ m[1:J], Int);
JuMP.@constraint(Rhs, [j = 1:J], L[j]b[j] ≥ m[j]);
JuMP.@objective(Rhs, Min, 0) # to be filled
JuMP.@variable(Rhs, 0 ≤ n[1:I, 1:J], Int)
JuMP.@constraint(Rhs, [i = 1:I], sum(view(n, i, :)) ≥ d[i, s])
JuMP.@constraint(Rhs, [j = 1:J], m[j]W[j] ≥ w ⋅ view(n, :, j));
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
# Rhs_1 = Rhs_model_vec[1];
# μ = -0.1 * o
# π = 0.1 * u
# JuMP.set_objective_coefficient.(Rhs_1, Rhs_1[:b], μ)
# JuMP.set_objective_coefficient.(Rhs_1, Rhs_1[:m], π)
# JuMP.optimize!(Rhs_1); JuMP.assert_is_solved_and_feasible(Rhs_1, allow_local = false)
# JuMP.objective_value(Rhs_1)
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
return Rhs
end; function Rhs_surrogate_plus_s(s) # 🔵 θ-Cuts can be reused, whereas `varycon` needs to be modified
Sur = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
JuMP.@variable(Sur, ϕ ≤ 500.0); # ⚠️ initial upper bound need to tune
JuMP.@objective(Sur, Max, ϕ)
JuMP.@variable(Sur, θ); # This is the `Rhs_surrogate`, And the follows are `plus`
JuMP.@variable(Sur, μ[1:J]);
JuMP.@variable(Sur, π[1:J]);
# θ-Cuts to be added, which reads `θ ≤ b ⋅ μ + m ⋅ π` in which (b, m) are vertices of the 🟢 Feasible Region of Rhs_model_s(s)
JuMP.@constraint(Sur, varycon, ϕ - θ ≤ false); # TBD
return Sur
end
prim_model_vec = [Min_s_prim_model(s) for s = 1:S]; # ✅ evaluate the 2nd stage ObjVals
Rhs_model_vec = [Rhs_model_s(s) for s = 1:S]; # ✅ evaluate the RHS of Lag cut
Sur_model_vec = [Rhs_surrogate_plus_s(s) for s = 1:S]; # ✅ used to generate trial coefficients of Lag cut
JuMP.set_silent.(prim_model_vec); JuMP.set_silent.(Rhs_model_vec); JuMP.set_silent.(Sur_model_vec);
Mst = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # ✅ used to generate trial 1st-stage solutions
JuMP.@variable(Mst, 0 ≤ b[1:J] ≤ 1, Int);
JuMP.@variable(Mst, 0 ≤ m[1:J], Int);
JuMP.@constraint(Mst, [j = 1:J], L[j]b[j] ≥ m[j]);
JuMP.@objective(Mst, Min, o ⋅ b + u ⋅ m);
JuMP.set_silent.(Mst);
JuMP.unset_silent(Rhs_model_vec[1]) # This might be hard
JuMP.set_attribute(Rhs_model_vec[1], "TimeLimit", 60)
suspect = nothing
JuMP.optimize!(Mst); JuMP.assert_is_solved_and_feasible(Mst; allow_local = false)
B = ı.(b); # to be globally ✂️ off
M = ı.(m); # to be globally ✂️ off
while true # iteratively find a violating cut for the Master problem
infeasible_block = 0
for s = 1:S
Min_s = prim_model_vec[s]; JuMP.fix.(Min_s[:m2], M);
JuMP.optimize!(Min_s); termination = JuMP.termination_status(Min_s)
if termination == JuMP.INFEASIBLE
infeasible_block = s
break # [early break]
elseif termination != JuMP.OPTIMAL
error("unknown termination status")
end
end
s = infeasible_block;
s == 0 && error("We have now derived a 1st-stage FEASIBLE solution!")
Sur = Sur_model_vec[s]; varycon = Sur[:varycon];
θ, μ, π = Sur[:θ], Sur[:μ], Sur[:π];
JuMP.set_normalized_coefficient.(varycon, μ, B);
JuMP.set_normalized_coefficient.(varycon, π, M);
JuMP.optimize!(Sur); JuMP.assert_is_solved_and_feasible(Sur; allow_local = false); ub = JuMP.objective_bound(Sur)
Μ, Π = ı.(μ), ı.(π); # to be locally ✂️ off
Θ = ı(θ); # [🪜Upper] to be locally ✂️ off
rhs = NaN; # [🪜Lower] pre-allocate
while true # try generate the coefficient of Lag cut until success
Rhs_s = Rhs_model_vec[s]
local b, m = Rhs_s[:b], Rhs_s[:m]
JuMP.set_objective_coefficient.(Rhs_s, b, Μ)
JuMP.set_objective_coefficient.(Rhs_s, m, Π)
JuMP.optimize!(Rhs_s); # JuMP.assert_is_solved_and_feasible(Rhs_s, allow_local = false)
if JuMP.termination_status(Rhs_s) == JuMP.TIME_LIMIT
suspect = Rhs_s
return # 🍅🍅🍅
end
rhs = JuMP.objective_value(Rhs_s)
lb = rhs - Μ ⋅ B - Π ⋅ M
if lb > 4e-3 # 🟣 global COT
@info "Find a cut for master with vio = $lb"
break # go back to add cut for the outer Master
end
local_B, local_M = ı.(b), ı.(m) # vertices
vertical_distance = Θ - (local_B ⋅ Μ + local_M ⋅ Π) # 🪜 having this since you approximate `Rhs_model_s`
if vertical_distance > 8e-4 # 🟣 local COT
JuMP.@constraint(Sur, θ ≤ local_B ⋅ μ + local_M ⋅ π)
else
error("weird 568435686586442")
end
JuMP.optimize!(Sur); JuMP.assert_is_solved_and_feasible(Sur; allow_local = false); ub = JuMP.objective_bound(Sur)
Μ, Π = ı.(μ), ı.(π); # to be locally ✂️ off
Θ = ı(θ); # [🪜Upper] to be locally ✂️ off
end
JuMP.@constraint(Mst, Μ ⋅ b + Π ⋅ m ≥ rhs)
JuMP.optimize!(Mst); JuMP.assert_is_solved_and_feasible(Mst; allow_local = false)
M = ı.(m); # to be globally ✂️ off
B = ı.(b); # to be globally ✂️ off
end
Gurobi.GRBwrite(JuMP.unsafe_backend(suspect), "data/suspect.mps")
On my machine it branches enormously (If I don’t set a TIME_LIMIT, it will branch more and non-stop in a reasonable time):
Gurobi Optimizer version 12.0.2 build v12.0.2rc0 (win64 - Windows 11.0 (26100.2))
CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Non-default parameters:
TimeLimit 60
Optimize a model with 28 rows, 108 columns and 207 nonzeros
Model fingerprint: 0x8bbf69b7
Variable types: 0 continuous, 108 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 5e+01]
Objective range [2e+01, 3e+04]
Bounds range [1e+00, 1e+00]
RHS range [6e+00, 4e+01]
Loaded MIP start from previous solve with objective -1.36424e-12
Presolve removed 7 rows and 7 columns
Presolve time: 0.00s
Presolved: 21 rows, 101 columns, 193 nonzeros
Variable types: 0 continuous, 101 integer (2 binary)
Root relaxation: objective -3.505153e+02, 53 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -350.51533 0 11 -0.00000 -350.51533 - - 0s
H 0 0 -51.7171719 -350.51533 578% - 0s
0 0 -339.13976 0 20 -51.71717 -339.13976 556% - 0s
H 0 0 -53.4817681 -338.94230 534% - 0s
0 0 -327.61043 0 20 -53.48177 -327.61043 513% - 0s
0 0 -327.61043 0 20 -53.48177 -327.61043 513% - 0s
0 0 -327.61043 0 19 -53.48177 -327.61043 513% - 0s
0 0 -327.61043 0 15 -53.48177 -327.61043 513% - 0s
0 2 -327.57472 0 13 -53.48177 -327.57472 512% - 0s
H 136 322 -199.4028229 -327.57472 64.3% 3.0 0s
H 264 322 -226.5814028 -327.57472 44.6% 3.2 0s
H 1138 806 -300.0734385 -327.57472 9.16% 2.9 0s
H 1578 442 -302.8250805 -326.42894 7.79% 2.9 0s
H33784 1717 -302.8251207 -315.62847 4.23% 2.3 2s
H33791 1636 -311.2578474 -315.62847 1.40% 2.3 2s
66994 5793 infeasible 91 -311.25785 -315.62847 1.40% 2.4 5s
184282 9413 infeasible 100 -311.25785 -315.62847 1.40% 2.1 10s
390579 24443 -315.62847 107 14 -311.25785 -315.62847 1.40% 1.6 15s
621856 32762 -315.62847 100 9 -311.25785 -315.62847 1.40% 1.4 20s
820322 32676 infeasible 90 -311.25785 -315.62847 1.40% 1.3 25s
1054675 32766 -315.62847 95 8 -311.25785 -315.62847 1.40% 1.2 30s
1271044 32907 infeasible 100 -311.25785 -315.62847 1.40% 1.2 35s
1483400 33327 -315.62847 115 9 -311.25785 -315.62847 1.40% 1.1 40s
1691494 32747 -315.62847 96 7 -311.25785 -315.62847 1.40% 1.1 45s
1876933 33645 infeasible 99 -311.25785 -315.62847 1.40% 1.1 50s
2030884 33891 infeasible 117 -311.25785 -315.62847 1.40% 1.1 55s
Cutting planes:
Gomory: 2
MIR: 38
StrongCG: 18
Inf proof: 2
Explored 2206702 nodes (2405517 simplex iterations) in 60.01 seconds (23.12 work units)
Thread count was 8 (of 8 available processors)
Solution count 9: -311.258 -302.825 -302.825 ... -1.36424e-12
Time limit reached
Best objective -3.112578474459e+02, best bound -3.156284748041e+02, gap 1.4042%
User-callback calls 4426781, time in user-callback 1.42 sec
julia>
julia> Gurobi.GRBwrite(JuMP.unsafe_backend(suspect), "data/suspect.mps")
0
After that, I open a new julia REPL and solve it directly, it is very quick.
julia> import Gurobi
julia> env = Gurobi.Env();
Set parameter Username
Set parameter LicenseID to value 2663413
Academic license - for non-commercial use only - expires 2026-05-11
julia> pModel = Ref{Ptr{Cvoid}}(C_NULL);
julia> Gurobi.GRBreadmodel(env, "data/suspect.mps", pModel)
Read MPS format model from file data/suspect.mps
Reading time = 0.00 seconds
: 28 rows, 108 columns, 207 nonzeros
0
julia> Gurobi.GRBoptimize(pModel[])
Gurobi Optimizer version 12.0.2 build v12.0.2rc0 (win64 - Windows 11.0 (26100.2))
CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 28 rows, 108 columns and 207 nonzeros
Model fingerprint: 0x8bbf69b7
Variable types: 0 continuous, 108 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 5e+01]
Objective range [2e+01, 3e+04]
Bounds range [1e+00, 1e+00]
RHS range [6e+00, 4e+01]
Presolve removed 7 rows and 7 columns
Presolve time: 0.00s
Presolved: 21 rows, 101 columns, 193 nonzeros
Variable types: 0 continuous, 101 integer (2 binary)
Found heuristic solution: objective 6470.8950980
Root relaxation: objective -3.505153e+02, 53 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -350.51533 0 11 6470.89510 -350.51533 105% - 0s
H 0 0 1617.7237745 -350.51533 122% - 0s
H 0 0 262.8039698 -350.51533 233% - 0s
H 0 0 46.5482749 -350.51533 853% - 0s
H 0 0 -53.4817681 -350.51533 555% - 0s
0 0 -339.13976 0 20 -53.48177 -339.13976 534% - 0s
0 0 -327.61043 0 20 -53.48177 -327.61043 513% - 0s
H 0 0 -99.3727799 -327.61043 230% - 0s
0 0 -327.61043 0 20 -99.37278 -327.61043 230% - 0s
0 0 -327.61043 0 19 -99.37278 -327.61043 230% - 0s
H 0 0 -179.8996033 -327.57472 82.1% - 0s
0 0 -327.57472 0 11 -179.89960 -327.57472 82.1% - 0s
0 0 -327.57472 0 15 -179.89960 -327.57472 82.1% - 0s
0 0 -327.57472 0 15 -179.89960 -327.57472 82.1% - 0s
0 0 -327.57472 0 17 -179.89960 -327.57472 82.1% - 0s
H 0 0 -226.5814028 -327.57472 44.6% - 0s
0 0 -327.57472 0 25 -226.58140 -327.57472 44.6% - 0s
0 0 -327.57472 0 17 -226.58140 -327.57472 44.6% - 0s
H 0 0 -251.5349570 -327.57472 30.2% - 0s
0 0 -327.57472 0 18 -251.53496 -327.57472 30.2% - 0s
0 0 -327.57472 0 18 -251.53496 -327.57472 30.2% - 0s
0 0 -327.57472 0 22 -251.53496 -327.57472 30.2% - 0s
0 0 -327.57472 0 22 -251.53496 -327.57472 30.2% - 0s
0 2 -327.57472 0 22 -251.53496 -327.57472 30.2% - 0s
H 71 152 -299.8524323 -327.57472 9.25% 1.4 0s
H 1687 920 -303.1600430 -327.57472 8.05% 1.8 0s
H 2084 969 -303.1600841 -324.45429 7.02% 1.8 0s
* 5703 486 104 -307.7487989 -315.62847 2.56% 2.2 0s
*16609 1257 128 -311.2578474 -311.25785 0.00% 2.6 1s
Cutting planes:
Gomory: 5
MIR: 15
StrongCG: 11
Inf proof: 1
Explored 17226 nodes (44703 simplex iterations) in 1.11 seconds (0.46 work units)
Thread count was 8 (of 8 available processors)
Solution count 10: -311.258 -307.749 -303.16 ... -53.4818
No other solutions better than -311.258
Optimal solution found (tolerance 1.00e-04)
Best objective -3.112578474459e+02, best bound -3.112578474459e+02, gap 0.0000%
0
Where is wrong?
.