Are you are adding user cuts and not lazy cuts on purpose? In the link you shared there is a distinction:
The main difference between user cuts and lazy constraints is that the former are not allowed to cut off integer-feasible solutions.
This works if Lazy attribute is set to 1, or Presolve is turned off.
You comment intrigued me, and I find something new. In this example, unlike my previous post, the user cuts are used by Gurobi as normal cuts.
import JuMP, Gurobi; GRB_ENV = Gurobi.Env()
function set_user_cut(model, c)
JuMP.MOI.set(JuMP.backend(model), Gurobi.ConstraintAttribute("Lazy"), JuMP.index(c), -1)
end
user_cut_vec = JuMP.ConstraintRef[]
model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.@variable(model, -0.8 <= x <= 1.8)
JuMP.@objective(model, Min, x)
for b in [0.1, 0.2]
push!(user_cut_vec, JuMP.@constraint(model, x >= b))
end
JuMP.optimize!(model); JuMP.objective_value(model) # 0.2
JuMP.set_binary(x)
JuMP.set_attribute(model, "LazyConstraints", 1)
set_user_cut.(model, user_cut_vec)
JuMP.optimize!(model)
JuMP.solution_summary(model; verbose = true) # x = 1.0 now, as opposed to previous `0`
Yes, purposefully. As far as Iโm concerned, user cuts embodies our a priori experience, whereas lazy cuts are those on-the-fly generated cuts that aims to cut off the current trial solution.
I donโt know if there is any relation to Presolve, I canโt find it in Gurobiโs doc.
I adapt your example to an MIQP problem
Minimize (x + 0.54)^2 - 1 over the interval [-1, 1] and x being Int.
I intend to use Gurobiโs MIP callback-lazy cuts.
But I have a novel ideaโadding a Phase I before we conventionally start.
The main question is:
๐
line is advantageous (e.g. has the potential to speed up the overall solving procedure)?The secondary question is: see the last for-loop (๐ )
import JuMP, Gurobi; GRB_ENV = Gurobi.Env()
import JuMP.value as v
COT = 1e-6 # cut-off tolerance
a = 0.54 # data of the problem
f(x) = (x + a)^2 - 1
fห(x) = 2 * (x + a)
function set_user_cut(model, c) # via Gurobi's API
JuMP.MOI.set(JuMP.backend(model), Gurobi.ConstraintAttribute("Lazy"), JuMP.index(c), -1)
end
function solve_to_optimality(model)
JuMP.optimize!(model)
JuMP.assert_is_solved_and_feasible(model; allow_local = false)
end
function my_callback_function(cb_data, cb_where::Cint)
cb_where != Gurobi.GRB_CB_MIPSOL && return
Gurobi.load_callback_variable_primal(cb_data, cb_where)
xv = JuMP.callback_value(cb_data, x)
yv = JuMP.callback_value(cb_data, y)
@info "current trial at x = $xv, y = $yv"
if yv < f(xv) - COT # โ
generate cuts on-the-fly
con = JuMP.@build_constraint(y >= f(xv) + fห(xv) * (x - xv))
@info "push & add lazy constr" func=con.func set=con.set
push!(from_MIP_cut_vec, con)
JuMP.MOI.submit(model, JuMP.MOI.LazyConstraint(cb_data), con)
end
end
# ๐ฃ Phase 1: Relax all "Int" constraints in the `model`
# Phase 1 is supposed to be fast, due to its continuous nature
model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.set_silent(model)
JuMP.@variable(model, -1.98 <= x <= 1.98)
JuMP.@variable(model, y >= -7) # the lower bound -7 is an artificial one to avoid initial unboundedness
# we intend to emulate the following Quad constraint via linear cuts
# JuMP.@constraint(model, y >= f(x))
JuMP.@objective(model, Min, y)
solve_to_optimality(model)
from_LP_cut_vec = JuMP.ConstraintRef[]
while true
if v(y) < f(v(x)) - COT # โ
generate cuts on-the-fly
push!(from_LP_cut_vec, JuMP.@constraint(model, y >= f(v(x)) + fห(v(x)) * (x - v(x))))
solve_to_optimality(model)
else
@info "After adding $(length(from_LP_cut_vec)) cuts, (x, y) converges to ($(v(x)), $(v(y)))"
break
end
end
# ๐ต Phase 2: This is the `model` with `Int` constrs that we intend to solve
JuMP.unset_silent(model)
set_user_cut.(model, from_LP_cut_vec) # ๐
JuMP.set_integer(x)
JuMP.set_attribute(model, "LazyConstraints", 1) # this setting is exclusively for the callback purpose
JuMP.MOI.set(model, Gurobi.CallbackFunction(), my_callback_function)
from_MIP_cut_vec = JuMP.ScalarConstraint[]
solve_to_optimality(model) # ๐ก From the Logging we observe that user cut is indeed adopted by Gurobi
JuMP.solution_summary(model; verbose = true) # This final result is correct
for i in from_MIP_cut_vec[1:2] # ๐ [question] why are they the same ??
@info i
end
julia> a = from_MIP_cut_vec[1]
JuMP.ScalarConstraint{JuMP.AffExpr, MathOptInterface.GreaterThan{Float64}}(y + 0.9199999999999999 x, MathOptInterface.GreaterThan{Float64}(-1.7084))
julia> b = from_MIP_cut_vec[2]
JuMP.ScalarConstraint{JuMP.AffExpr, MathOptInterface.GreaterThan{Float64}}(y + 0.9199999999999999 x, MathOptInterface.GreaterThan{Float64}(-1.7084))
julia> a.func == b.func && a.set == b.set
true
julia> a == b
false
Therefore a related question is: how to remove this sort of redundancy, such that one cut occurs only once in the collection from_MIP_cut_vec?
Your code is incorrect. If you declare that a constraint is a user-cut, it cannot remove any integer feasible solutions.
The Lazy attribute is likely ignored in the LP.
When you solve the MIP, your constraints cut off the x = 0 solution.
You may want to do some research on the difference between lazy constraints and user cuts. For example, read OR in an OB World: User Cuts versus Lazy Constraints
That post was outdated. It was a nascent one. The latest one is more meaningful. Please take a look.
As with many things MIP, this is almost certainly problem dependent. Try it out on your model and see. If you spend a lot of time generating useless cuts, it probably isnโt useful. If you can generate a few good cuts, it probably is. My expectation for your current model is โnoโ.
The secondary question is: see the last for-loop (
๐)
[question] why are they the same ??
If you add a print to observe JuMP.callback_value(cb_data, y), youโll see that they are generated at different (x, y) points. But your constraint does not depend one the value of y. Gurobi doesnโt know this, so you end up generating the same constraint twice.
why MOI has this behavior?
Because we havenโt added a method for == between ScalarConstraint objects. I also donโt really think we need to. This is a pretty uncommon use-case.
Iโve revised my code above, according to your advice.
Now the only issue is that Gurobiโs 2nd trial point violates the 1st lazy constraint, as shown in
[ Info: current trial at x = -1.0, y = -0.9662027806043625
โ Info: push & add lazy constr # ๐ 1st lazy constraint
โ func = y + 0.9199999999999999 x
โ set = MathOptInterface.GreaterThan{Float64}(-1.7084)
Presolve removed 13 rows and 0 columns
Presolve time: 0.00s
Presolved: 0 rows, 2 columns, 0 nonzeros
Crushed 13 out of 13 user cuts to presolved model
Variable types: 1 continuous, 1 integer (0 binary)
[ Info: current trial at x = -1.0, y = -1.0000200521257165 # ๐ 2nd trial point
But whatever, this is a minor issue, as long as the final result is correct. You answer is fantastic. 
This should be due to the nature of this exemplary problem.
I here furnish a large-scale case. I hope that given this new case, the user cut will be advantageous.
import JuMP, Gurobi; import JuMP.value as ๐; GRB_ENV = Gurobi.Env()
using Logging; global_logger(ConsoleLogger(Logging.Info))
import Random
Random.seed!(1);
COT = 1e-4
# Test case 1๏ธโฃ modest scale
# M, N = 9, 10; โญ = N/2
# ๐บ, C_y, C_z, ๐, ๐ฑ, ๐
= let
# ๐ = [2.75 12.33 8.58 0.67 13.58 7.33 9.42 5.67 9.83 4.0; 1.08 5.67 6.5 14.42 0.67 4.0 12.75 6.08 4.0 9.42; 4.42 6.92 12.75 4.42 11.5 7.75 6.92 7.33 13.17 11.08; 4.0 11.08 4.0 13.58 12.33 5.25 4.83 4.0 11.5 14.0; 3.17 1.08 11.92 8.58 11.5 1.08 7.75 1.5 5.67 3.17; 3.58 2.75 11.5 2.75 2.33 14.42 4.83 5.67 3.17 10.25; 11.92 9.83 11.92 2.33 2.75 13.17 1.92 2.33 12.75 8.17; 4.0 14.0 4.0 12.33 1.08 9.83 1.92 10.67 5.67 3.58; 14.0 9.83 12.33 14.0 7.33 11.92 8.58 2.75 12.33 10.67]
# C_z = 2 * [4.0, 2.5, 4.0, 4.0, 2.5, 2.5, 3.0, 3.5, 3.0]
# C_y = [10, 12, 10, 9, 12, 9, 12, 9, 11] .* C_z
# ๐บ = [224, 224, 221, 223, 221, 226, 224, 219, 218]
# ๐ฑ = [394, 70, 238, 6, 225, 383, 65, 136, 13, 83]
# ๐
= [40, 40, 40, 30, 10, 10, 40, 10, 30, 30]
# ๐บ, C_y, C_z, ๐, ๐ฑ, ๐
# end; ๐ = 20 * C_z;
# Test case 2๏ธโฃ large scale
M, N = 69, 100; โญ = N/2
๐บ, C_y, C_z, ๐, ๐ฑ, ๐
= let
๐ = R_x = rand(0.3:.00117:8.7, M, N)
C_z = rand(2:0.00117:5, M)
C_y = rand(9:0.017:15, M) .* C_z
๐บ = [2.42, 2.42, 1.78, 2.7, 2.06, 1.78, 1.14, 3.5, 3.1, 2.02, 1.82, 3.42, 2.42, 0.9, 3.06, 1.18, 1.7, 3.46, 1.46, 1.62, 1.26, 1.3, 3.46, 2.66, 1.82, 2.34, 1.54, 1.82, 2.58, 2.9, 2.66, 3.14, 2.5, 1.62, 2.78, 3.3, 1.38, 1.26, 2.74, 1.82, 3.14, 1.38, 2.74, 0.94, 3.06, 3.22, 2.62, 2.42, 3.02, 3.18, 1.58, 2.7, 3.06, 1.7, 3.22, 2.74, 2.54, 3.26, 1.3, 2.38, 3.3, 1.02, 2.42, 2.38, 3.1, 1.3, 2.02, 0.94, 3.46];
๐ฑ = [0.11, 0.57, 0.77, 0.35, 0.49, 0.25, 0.25, 0.45, 0.05, 0.77, 0.39, 0.83, 0.47, 0.55, 0.43, 0.47, 0.29, 0.49, 0.47, 0.75, 0.09, 0.59, 0.11, 0.09, 0.11, 0.49, 0.29, 0.61, 0.23, 0.53, 0.15, 0.85, 0.17, 0.23, 0.25, 0.57, 0.65, 0.47, 0.65, 0.11, 0.09, 0.13, 0.57, 0.29, 0.19, 0.39, 0.31, 0.63, 0.61, 0.51, 0.75, 0.41, 0.57, 0.59, 0.47, 0.05, 0.29, 0.49, 0.53, 0.79, 0.21, 0.83, 0.13, 0.31, 0.63, 0.43, 0.05, 0.33, 0.77, 0.43, 0.11, 0.41, 0.19, 0.61, 0.89, 0.43, 0.15, 0.35, 0.43, 0.63, 0.59, 0.89, 0.23, 0.23, 0.33, 0.63, 0.69, 0.35, 0.79, 0.31, 0.69, 0.81, 0.15, 0.19, 0.27, 0.29, 0.17, 0.83, 0.85, 0.45];
๐
= [0.56, 0.26, 0.5, 0.5, 0.71, 0.77, 0.32, 0.53, 0.38, 0.2, 0.2, 0.29, 0.68, 0.59, 0.38, 0.29, 0.38, 0.41, 0.26, 0.35, 0.56, 0.68, 0.8, 0.8, 0.62, 0.44, 0.47, 0.53, 0.68, 0.44, 0.26, 0.35, 0.44, 0.59, 0.59, 0.41, 0.53, 0.8, 0.29, 0.62, 0.29, 0.59, 0.23, 0.56, 0.68, 0.2, 0.44, 0.74, 0.65, 0.62, 0.8, 0.71, 0.35, 0.56, 0.23, 0.44, 0.65, 0.35, 0.38, 0.38, 0.44, 0.47, 0.5, 0.38, 0.74, 0.29, 0.2, 0.29, 0.62, 0.23, 0.41, 0.8, 0.77, 0.2, 0.35, 0.65, 0.32, 0.23, 0.74, 0.38, 0.68, 0.56, 0.59, 0.32, 0.8, 0.74, 0.38, 0.2, 0.5, 0.68, 0.74, 0.8, 0.68, 0.53, 0.44, 0.59, 0.26, 0.68, 0.44, 0.38];
๐บ, C_y, C_z, ๐, ๐ฑ, ๐
end; ๐ = 20 * C_z;
function sm(p, x) return sum(p .* x) end;
function optimise(m)
JuMP.optimize!(m)
return (
JuMP.termination_status(m),
JuMP.primal_status(m),
JuMP.dual_status(m)
)
end;
macro set_objective_function(m, f) return esc(:(JuMP.set_objective_function($m, JuMP.@expression($m, $f)))) end;
function Model(name)
m = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.MOI.set(m, JuMP.MOI.Name(), name)
s = name[end-1]
s == 'm' && JuMP.set_objective_sense(m, JuMP.MIN_SENSE)
s == 'M' && JuMP.set_objective_sense(m, JuMP.MAX_SENSE)
last(name) == 's' && JuMP.set_silent(m)
m
end;
function solve_to_normality(m)
t, p, d = optimise(m)
name = JuMP.name(m)
t == JuMP.OPTIMAL || error("$name: $(string(t))")
p == JuMP.FEASIBLE_POINT || error("$name: (primal) $(string(p))")
(name[end-2] == 'l' && d != p) && error("$name: (dual) $(string(d))")
return JuMP.result_count(m)
end;
function stage2_problem_primal(z, d)
st2 = Model("st2_primal_lms")
JuMP.@variable(st2, x[i = 1:M, j = 1:N] โฅ 0)
JuMP.@variable(st2, ฮถ[i = 1:M] โฅ 0) # an expensive substitute for `z`
JuMP.@constraint(st2, ๐ธ[i = 1:M], z[i] + ฮถ[i] โฅ sum(x[i, :]))
JuMP.@constraint(st2, ๐น[j = 1:N], sum(x[:, j]) โฅ d[j] )
@set_objective_function(st2, sm(๐, x) + sm(๐, ฮถ))
solve_to_normality(st2)
return JuMP.objective_value(st2)
end;
function g2d(g) return ๐ฑ .+ ๐
.* g end;
function set_๐ป_obj(z) return JuMP.set_objective_coefficient.(๐ป, ฤฑ, -z) end # a โ
Partial โ
setting
function set_๐ณ_obj(z, g) return JuMP.set_objective_coefficient.(๐ณ, [๐น; ๐ธ], [g2d(g); -z]) end
function set_๐ถ_obj(z, I, J) return @set_objective_function(๐ถ, sm(J, g2d(๐)) - sm(I, z)) end # also involve a constant
function set_user_cut(m, c) return JuMP.MOI.set(JuMP.backend(m), Gurobi.ConstraintAttribute("Lazy"), JuMP.index(c), -1) end
function mip_cb_function(cb_data, cb_where::Cint) # only used in Phase 2
cb_where != Gurobi.GRB_CB_MIPSOL && return
Gurobi.load_callback_variable_primal(cb_data, cb_where)
common_value = JuMP.callback_value(cb_data, common_expr)
o = JuMP.callback_value(cb_data, ๐)
z = JuMP.callback_value.(cb_data, ๐ฃ)
# then we seek a hazardous scene given the 1st stage trial solution `z`
set_๐ป_obj(z); Lหs = (ul = [Inf, -Inf], g = Vector{Float64}(undef, N)); solve_to_normality(๐ป);
Lหs.ul[1], Lหs.g[:] = JuMP.objective_bound(๐ป), ๐.(๐ ) # setting upper bound is one-off, get a heuristic primal solution
ub = common_value + Lหs.ul[1]
if ub < ๐ผหs.ul[1]
๐ผหs.ul[1], ๐ผหs.z[:] = ub, z
@info "update (global)ub to $ub โช"
end
while true
set_๐ณ_obj(z, Lหs.g); solve_to_normality(๐ณ);
set_๐ถ_obj(z, ๐.(๐ธ), ๐.(๐น)); solve_to_normality(๐ถ);
lb = JuMP.objective_value(๐ถ)
if lb > Lหs.ul[2]
Lหs.ul[2], Lหs.g[:] = lb, ๐.(๐)
else
@debug "After RLT_BCA, (sub)lb = $(Lหs.ul[2]) < $(Lหs.ul[1]) = (sub)ub"
break
end
end
set_๐ณ_obj(z, Lหs.g); solve_to_normality(๐ณ); cn, pz = sm(๐.(๐น), g2d(Lหs.g)), -๐.(๐ธ)
if o < cn + sm(pz, z) - COT
con = JuMP.@build_constraint(๐ >= cn + sm(pz, ๐ฃ))
JuMP.MOI.submit(๐ผ, JuMP.MOI.LazyConstraint(cb_data), con)
end
end
๐ผ = Model("st1_bms"); JuMP.@variable(๐ผ, ๐);
JuMP.@variables(๐ผ, begin 0 โค ๐ข[1:M] โค 1, Int; ๐ฃ[1:M] โฅ 0 end); JuMP.@constraint(๐ผ, ๐ฃ .โค ๐บ .* ๐ข); JuMP.set_objective_coefficient.(๐ผ, [๐ข; ๐ฃ], [C_y; C_z]);
๐ณ = Model("st2_dual_lMs");
JuMP.@variable(๐ณ, 0 โค ๐ธ[i = 1:M] โค ๐[i]);
JuMP.@variable(๐ณ, 0 โค ๐น[j = 1:N]); JuMP.@constraint(๐ณ, [j = 1:N, i = 1:M],
๐น[j] โค ๐[i, j] + ๐ธ[i]
); ๐ถ = Model("opt_g_lMs"); JuMP.@variable(๐ถ, 0 โค ๐[n = 1:N] โค 1); JuMP.@constraint(๐ถ, sum(๐) โค โญ);
๐ป = Model("lpr_lMs"); JuMP.@variable(๐ป, 0 โค ๐ [n = 1:N] โค 1); JuMP.@constraint(๐ป, sum(๐ ) โค โญ);
JuMP.@variable(๐ป, 0 โค ฤฑ[i = 1:M] โค ๐[i]);
JuMP.@variable(๐ป, 0 โค ศท[j = 1:N]); JuMP.@constraint(๐ป, ๐[i = 1:M, j = 1:N],
ศท[j] โค ๐[i, j] + ฤฑ[i]
);
JuMP.@variable(๐ป, 0 โค ๐ ศท[nj = 1:N]) # only diagonal, as it occurs in Obj
JuMP.@variable(๐ป, 0 โค ๐ Xฤฑ[n = 1:N, i = 1:M]) # needed to specify the upper_bound of `๐ ศท`
JuMP.set_objective_coefficient.(๐ป, [ศท; ๐ ศท], [๐ฑ; ๐
]) # Obj: the โ
Fixed โ
part
JuMP.@constraint(๐ป, ๐ ศท .โค ศท) # (๐ โค 1) โจ (0 โค ศท) [a part of 1/10] โ
JuMP.@constraint(๐ป, [n = 1:N, i = 1:M], ๐ ศท[n] โค ๐ [n] * ๐[i, n] + ๐ Xฤฑ[n, i]) # (g โฅ 0) โจ ๐ [a part of 1/10] โ
JuMP.@constraint(๐ป, [n = 1:N, i = 1:M], ๐ Xฤฑ[n, i] โค ๐[i] * ๐ [n]) # (g โฅ 0) โจ (๐ธ โค ๐) [1/10] ๐ฃ
JuMP.@constraint(๐ป, [n = 1:N, i = 1:M], ๐ Xฤฑ[n, i] โค ฤฑ[i]) # (๐ โค 1) โจ (0 โค ๐ธ) [1/10] ๐ฃ
JuMP.@constraint(๐ป, [i = 1:M], sum(๐ Xฤฑ[:, i]) โค โญ * ฤฑ[i]) # (sum(๐ ) โค โญ) โจ (0 โค ๐ธ) [1/10]
@info "๐ฃ Start Phase 0: adding an initial cut"
JuMP.unset_integer.(๐ข); ๐ผหs = (ul = [Inf, -Inf], z = Vector{Float64}(undef, M))
solve_to_normality(๐ผ); z = ๐.(๐ฃ); # โ๏ธ the initial solve of `st1_bms`, No cut yet, and lb is invalid
set_๐ป_obj(z); Lหs = (ul = [Inf, -Inf], g = Vector{Float64}(undef, N)); solve_to_normality(๐ป);
Lหs.ul[1], Lหs.g[:] = JuMP.objective_bound(๐ป), ๐.(๐ ) # setting upper bound is one-off, get a heuristic primal solution
while true
set_๐ณ_obj(z, Lหs.g); solve_to_normality(๐ณ)
set_๐ถ_obj(z, ๐.(๐ธ), ๐.(๐น)); solve_to_normality(๐ถ);
lb = JuMP.objective_value(๐ถ)
if lb > Lหs.ul[2]
Lหs.ul[2], Lหs.g[:] = lb, ๐.(๐)
else
@debug "After RLT_BCA, (sub)lb = $(Lหs.ul[2]) < $(Lหs.ul[1]) = (sub)ub"
break
end
end
๐ผหs.ul[1], ๐ผหs.z[:] = JuMP.objective_value(๐ผ) + Lหs.ul[1], z # The feasible value and solution at initialization
JuMP.@expression(๐ผ, common_expr, JuMP.objective_function(๐ผ)); JuMP.set_objective_coefficient(๐ผ, ๐, 1); # โ
a one-shot turn on
set_๐ณ_obj(z, Lหs.g); solve_to_normality(๐ณ); the_first_cut = JuMP.@constraint(๐ผ, ๐ >= sm(๐.(๐น), g2d(Lหs.g)) - sm(๐.(๐ธ), ๐ฃ))
@info "๐ฃ Start Phase 1: adding some cuts for the continuous relaxation of the master problem"
from_LP_cut_vec = JuMP.ConstraintRef[]; push!(from_LP_cut_vec, the_first_cut)
solve_to_normality(๐ผ); # โ
This is indeed a regular solve
๐ผหs.ul[2], z = JuMP.objective_bound(๐ผ), ๐.(๐ฃ);
@info "The 1st (global)lb = $(๐ผหs.ul[2]) < $(๐ผหs.ul[1]) = (global)ub, then we start the main loop of Phase 1"
t0 = time()
for ite = 1:typemax(Int)
global z, Lหs, t0
set_๐ป_obj(z); Lหs = (ul = [Inf, -Inf], g = Vector{Float64}(undef, N)); solve_to_normality(๐ป);
Lหs.ul[1], Lหs.g[:] = JuMP.objective_bound(๐ป), ๐.(๐ ) # setting upper bound is one-off, get a heuristic primal solution
ub = ๐(common_expr) + Lหs.ul[1]
if ub < ๐ผหs.ul[1]
๐ผหs.ul[1], ๐ผหs.z[:] = ub, z
@info "ite = $ite โถ $(๐ผหs.ul[2]) < $(๐ผหs.ul[1]) โช"
else
@info "ite = $ite โถ $(๐ผหs.ul[2]) < $(๐ผหs.ul[1]) = ubs | ub = $ub"
end
while true
set_๐ณ_obj(z, Lหs.g); solve_to_normality(๐ณ);
set_๐ถ_obj(z, ๐.(๐ธ), ๐.(๐น)); solve_to_normality(๐ถ);
lb = JuMP.objective_value(๐ถ)
if lb > Lหs.ul[2]
Lหs.ul[2], Lหs.g[:] = lb, ๐.(๐)
else
@debug "After RLT_BCA, (sub)lb = $(Lหs.ul[2]) < $(Lหs.ul[1]) = (sub)ub"
break
end
end
set_๐ณ_obj(z, Lหs.g); solve_to_normality(๐ณ); cn, pz = sm(๐.(๐น), g2d(Lหs.g)), -๐.(๐ธ)
if ๐(๐) < cn + sm(pz, z) - COT
push!(from_LP_cut_vec, JuMP.@constraint(๐ผ, ๐ >= cn + sm(pz, ๐ฃ)))
else
@info "โถ Phase 1: After $(time() - t0) seconds, we find $(length(from_LP_cut_vec)) valid cuts"
break
end
solve_to_normality(๐ผ); # โ
This is indeed a regular solve
๐ผหs.ul[2], z = JuMP.objective_bound(๐ผ), ๐.(๐ฃ);
end
@info "๐ฃ Start Phase 2: employ the callback+lazy_cut method for MIP"
JuMP.unset_silent(๐ผ)
set_user_cut.(๐ผ, from_LP_cut_vec)
JuMP.set_binary.(๐ข); ๐ผหs = (ul = [Inf, -Inf], z = Vector{Float64}(undef, M))
JuMP.set_attribute(๐ผ, "LazyConstraints", 1)
JuMP.MOI.set(๐ผ, Gurobi.CallbackFunction(), mip_cb_function)
solve_to_normality(๐ผ) # MIP with callback_lazyCut
@info "โถ Phase 2:After $(JuMP.solve_time(๐ผ)) seconds, (global)lb = $(JuMP.objective_bound(๐ผ)) < $(๐ผหs.ul[1]) = (global)ub"
set_user_cut line, and delete the code block of ๐ฃ Start Phase 1. This way it takes โฅ 1 hour to return to REPL.Yet the convergence quality (in terms of lb and ub) of the 2 methods are almost at the same level. (I am not expecting exact convergence, since the application is a 2SARO, which has a NP hard subproblem in it. I adopted an approximate algorithmโRLT-BCAโto generate Bendersโ cut.)
Update: I find that the constraint attribute Lazy can be modified from -1 (user cut) to 1/2/3 (lazy cut). It will be faster at 1, and even faster at mode 2 and 3.
I have a question on how to use multithread (i.e. parallel computing) in Julia.
i.e., let my computer execute the ๐
for-loop below in parallel (theoretically achieving N times faster than the current style).
# start a new julia REPL
# load Line 1โ128 of my code in my last post
# then follows
function do_BCA_at(g)
stable_lb = -Inf
while true
set_๐ณ_obj(z, g); solve_to_normality(๐ณ)
set_๐ถ_obj(z, ๐.(๐ธ), ๐.(๐น)); solve_to_normality(๐ถ);
lb = JuMP.objective_value(๐ถ)
if lb > stable_lb
stable_lb, g = lb, ๐.(๐)
else
return stable_lb, g
end
end
end
# I need to call the above function N times. Note that each call is independent of the other (N-1) calls
# therefore it's beneficial to carry this out in parallel
# now suppose N = 2 for demo purpose
g1 = zeros(N)
g2 = [1; zeros(N-1)]
g_vec = [g1, g2]
g_stable_vec = Vector{Vector{Float64}}(undef, 2)
stable_lb_vec = Vector{Float64}(undef, 2)
for i in 1:2 # ๐
stable_lb_vec[i], g_stable_vec[i] = do_BCA_at(g_vec[i])
end
@info "see" stable_lb_vec
I find this
# This usage is incorrect
import JuMP, Gurobi
model = JuMP.Model(Gurobi.Optimizer)
JuMP.@variable(model, x >= 0)
JuMP.set_objective_sense(model, JuMP.MOI.MIN_SENSE)
Threads.@threads for i in 1:4
JuMP.set_objective_coefficient(model, x, i)
JuMP.optimize!(model)
end
# This usage might be correct, but it entails model building
import JuMP, Gurobi
v = ones(10)
Threads.@threads for i in eachindex(v)
model = JuMP.Model(Gurobi.Optimizer)
JuMP.@variable(model, x >= 0)
JuMP.@objective(model, Min, i * (x + 0.5))
JuMP.optimize!(model)
v[i] = JuMP.objective_value(model)
end
I find an intriguing behavior of Gurobiโit does NOT memorize lazy cuts.
The model is a MIP, I solve it with Gurobiโs callback + lazy cuts.
It is reported that Gurobi add at least 3 lazy cuts, which can make model bounded below, even without the presence of lower_bound(y).
Then I delete_lower_bound(y) and solve it again, finding that Gurobi reports unbounded.
Please run the following code in REPL to see what happened (I haved made detailed illustrations using @info to facilitate comprehension).
import JuMP, Gurobi; GRB_ENV = Gurobi.Env()
COT = 1e-6 # cut-off tolerance
a = 0.54 # data of the problem
f(x) = (x + a)^2 - 1
fห(x) = 2 * (x + a)
function cbf(cb_data, cb_where::Cint)
cb_where != Gurobi.GRB_CB_MIPSOL && return
Gurobi.load_callback_variable_primal(cb_data, cb_where)
xv, yv = JuMP.callback_value.(cb_data, [x, y])
@info "current trial at P = (x = $xv; y = $yv)"
if yv < f(xv) - COT
con = JuMP.@build_constraint(y >= f(xv) + fห(xv) * (x - xv))
@info "P violates 'y โฅ f(x)',thus we build this constr: LHS โฅ 0, in which LHS = $(con.func) - $(con.set.lower)"
slope_y, slope_x = con.func.terms[y], con.func.terms[x]
@info "the LHS evaluated at P is $slope_y * $yv + $slope_x * $xv - $(con.set.lower) = $(slope_y * yv + slope_x * xv - con.set.lower)"
push!(from_MIP_cut_vec, con)
JuMP.MOI.submit(model, JuMP.MOI.LazyConstraint(cb_data), con)
end
end
model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.@variable(model, -1.98 <= x <= 1.98, Int)
JuMP.@variable(model, y >= -7)
JuMP.@objective(model, Min, y)
JuMP.set_attribute(model, "LazyConstraints", 1); JuMP.MOI.set(model, Gurobi.CallbackFunction(), cbf)
from_MIP_cut_vec = JuMP.ScalarConstraint[]
JuMP.optimize!(model);
@info "The initial solve is normal" JuMP.solution_summary(model; verbose = true)
JuMP.delete_lower_bound(y) # We made a modification to `model`
@error "The following 2nd solve is problematic, which might indicates that Gurobi doesn't memorize lazy cuts"
JuMP.optimize!(model); # โ Gurobi: Model is unbounded
To remedy this, we might need to introduce an if-else block as the ๐
follows
import JuMP, Gurobi; GRB_ENV = Gurobi.Env()
COT = 1e-6 # cut-off tolerance
a = 0.54 # data of the problem
f(x) = (x + a)^2 - 1
fห(x) = 2 * (x + a)
function cbf(cb_data, cb_where::Cint)
cb_where != Gurobi.GRB_CB_MIPSOL && return
Gurobi.load_callback_variable_primal(cb_data, cb_where)
xv, yv = JuMP.callback_value.(cb_data, [x, y])
@info "current trial at P = (x = $xv; y = $yv)"
if yv < f(xv) - COT
con = JuMP.@build_constraint(y >= f(xv) + fห(xv) * (x - xv))
xv != last(from_MIP_xs) && push!(from_MIP_xs, xv) # ๐
JuMP.MOI.submit(model, JuMP.MOI.LazyConstraint(cb_data), con)
end
end
model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.@variable(model, -1.98 <= x <= 1.98, Int)
JuMP.@variable(model, y >= -7)
JuMP.@objective(model, Min, y)
JuMP.set_attribute(model, "LazyConstraints", 1); JuMP.MOI.set(model, Gurobi.CallbackFunction(), cbf)
from_MIP_xs = [-2.0] # ๐
arbitrarily introduce a value which is outside of the feasible region
JuMP.optimize!(model) # the 1st solve
# ๐
the following if-else makes provision for the 2nd solve
if length(from_MIP_xs) <= 1
error("No lazy cuts are added")
else
popfirst!(from_MIP_xs)
for xv in from_MIP_xs
c = JuMP.@constraint(model, y >= f(xv) + fห(xv) * (x - xv))
JuMP.MOI.set(JuMP.backend(model), Gurobi.ConstraintAttribute("Lazy"), JuMP.index(c), 1)
end
end
JuMP.delete_lower_bound(y) # We made a modification to `model`
JuMP.optimize!(model) # the 2nd solve
JuMP.solution_summary(model; verbose = true)
See the following code. Does this constitute a bug of Gurobi?
import JuMP, Gurobi; GRB_ENV = Gurobi.Env()
COT = 1e-6 # cut-off tolerance
a = 0.54 # data of the problem
f(x) = (x + a)^2 - 1
fห(x) = 2 * (x + a)
function cbf(cb_data, cb_where::Cint)
cb_where != Gurobi.GRB_CB_MIPSOL && return
Gurobi.load_callback_variable_primal(cb_data, cb_where)
xv, yv = JuMP.callback_value.(cb_data, [x, y])
@info "current trial at P = (x = $xv; y = $yv)"
if yv < f(xv) - COT # the Quad constraint is violated at the current trial point
con = JuMP.@build_constraint(y โฅ f(xv) + fห(xv) * (x - xv))
neg_con_set_lower = -con.set.lower
@info "P violates 'y โฅ f(x)',thus we build this constr: LHS โฅ 0, in which LHS = $(con.func) + $neg_con_set_lower"
slope_y, slope_x = con.func.terms[y], con.func.terms[x]
@info "the LHS evaluated at P is $slope_y * $yv + $slope_x * $xv + $neg_con_set_lower = $(slope_y * yv + slope_x * xv + neg_con_set_lower)"
JuMP.MOI.submit(model, JuMP.MOI.LazyConstraint(cb_data), con)
end
end
model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.@variable(model, -1.98 <= x <= 1.98, Int)
JuMP.@variable(model, y)
JuMP.@objective(model, Min, y)
JuMP.set_attribute(model, "LazyConstraints", 1); JuMP.MOI.set(model, Gurobi.CallbackFunction(), cbf)
JuMP.optimize!(model);
error("model: ", JuMP.termination_status(model))
From the Loggings (both Gurobiโs and Juliaโs) we know that weโve summitted these 2 constraints
[ Info: LHS โฅ 0, in which LHS = y + 0.9199999999999999 x + 1.7084
[ Info: LHS โฅ 0, in which LHS = y - 3.08 x + 1.7084000000000001
Therefore y is certainly bounded below. Yet Gurobi reports unbounded
julia> JuMP.optimize!(model);
Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (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:
LazyConstraints 1
Optimize a model with 0 rows, 2 columns and 0 nonzeros
Model fingerprint: 0xb7a8a0a7
Variable types: 1 continuous, 1 integer (0 binary)
Coefficient statistics:
Matrix range [0e+00, 0e+00]
Objective range [1e+00, 1e+00]
Bounds range [2e+00, 2e+00]
RHS range [0e+00, 0e+00]
[ Info: current trial at P = (x = -1.0; y = -1.0e30)
[ Info: P violates 'y โฅ f(x)',thus we build this constr: LHS โฅ 0, in which LHS = y + 0.9199999999999999 x + 1.7084
[ Info: the LHS evaluated at P is 1.0 * -1.0e30 + 0.9199999999999999 * -1.0 + 1.7084 = -1.0e30
[ Info: current trial at P = (x = 1.0; y = -1.0e30)
[ Info: P violates 'y โฅ f(x)',thus we build this constr: LHS โฅ 0, in which LHS = y - 3.08 x + 1.7084000000000001
[ Info: the LHS evaluated at P is 1.0 * -1.0e30 + -3.08 * 1.0 + 1.7084000000000001 = -1.0e30
[ Info: current trial at P = (x = 1.0; y = -1.0e30)
[ Info: P violates 'y โฅ f(x)',thus we build this constr: LHS โฅ 0, in which LHS = y - 3.08 x + 1.7084000000000001
[ Info: the LHS evaluated at P is 1.0 * -1.0e30 + -3.08 * 1.0 + 1.7084000000000001 = -1.0e30
[ Info: current trial at P = (x = 1.0; y = -1.0e30)
[ Info: P violates 'y โฅ f(x)',thus we build this constr: LHS โฅ 0, in which LHS = y - 3.08 x + 1.7084000000000001
[ Info: the LHS evaluated at P is 1.0 * -1.0e30 + -3.08 * 1.0 + 1.7084000000000001 = -1.0e30
Presolve time: 0.00s
Presolved: 0 rows, 2 columns, 0 nonzeros
Variable types: 1 continuous, 1 integer (0 binary)
Root relaxation: unbounded, 0 iterations, 0.00 seconds (0.00 work units)
Explored 1 nodes (0 simplex iterations) in 0.97 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)
Solution count 0
Model is unbounded
Best objective -, best bound -, gap -
User-callback calls 76, time in user-callback 1.11 sec
This is the most concise case that explain things, the above posts might be skipped
import JuMP, Gurobi; GRB_ENV = Gurobi.Env()
function cbf(cb_data, cb_where::Cint)
cb_where != Gurobi.GRB_CB_MIPSOL && return
Gurobi.load_callback_variable_primal(cb_data, cb_where)
xv, yv = JuMP.callback_value.(cb_data, [x, y])
if xv == -1.0
@info "current trial at P = (x = $xv; y = $yv)"
con = JuMP.@build_constraint(y โฅ 0.6 - x)
JuMP.MOI.submit(model, JuMP.MOI.LazyConstraint(cb_data), con)
elseif xv == 1.0 || xv == 2.0
@info "current trial at P = (x = $xv; y = $yv)"
con = JuMP.@build_constraint(y โฅ x - 0.6)
JuMP.MOI.submit(model, JuMP.MOI.LazyConstraint(cb_data), con)
else
@warn "This is an other trial point at (x = $xv, y = $yv)"
end
end
# The normal case
model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.set_attribute(model, "LazyConstraints", 1)
JuMP.MOI.set(model, Gurobi.CallbackFunction(), cbf)
JuMP.@variable(model, -1.98 <= x <= 2.98, Int)
JuMP.@variable(model, y >= -0.12)
JuMP.@objective(model, Min, y)
JuMP.optimize!(model);
JuMP.solution_summary(model; verbose = true)
# abnormal case
model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
JuMP.set_attribute(model, "LazyConstraints", 1)
JuMP.MOI.set(model, Gurobi.CallbackFunction(), cbf)
JuMP.@variable(model, -1.98 <= x <= 2.98, Int)
JuMP.@variable(model, y)
JuMP.@objective(model, Min, y)
JuMP.optimize!(model); # The first solve
c_left = JuMP.@constraint(model, y โฅ 0.6 - x)
JuMP.MOI.set(JuMP.backend(model), Gurobi.ConstraintAttribute("Lazy"), JuMP.index(c_left), 1)
JuMP.optimize!(model); # The second solve
JuMP.solution_summary(model; verbose = true)
You are adding your callback only at MIPSOL nodes. Try adding them at MIPNODE as well.
julia> using JuMP, Gurobi
julia> function my_callback_fn(cb_data, cb_where::Cint)
if cb_where == Gurobi.GRB_CB_MIPSOL || cb_where == Gurobi.GRB_CB_MIPNODE
Gurobi.load_callback_variable_primal(cb_data, cb_where)
xv, yv = JuMP.callback_value.(cb_data, [x, y])
if !(yv >= 0.6 - xv - 1e-6)
con = JuMP.@build_constraint(y >= 0.6 - x)
MOI.submit(model, MOI.LazyConstraint(cb_data), con)
elseif !(yv >= xv - 0.6 - 1e-6)
con = JuMP.@build_constraint(y >= x - 0.6)
MOI.submit(model, MOI.LazyConstraint(cb_data), con)
end
end
return
end
my_callback_fn (generic function with 1 method)
julia> model = direct_model(Gurobi.Optimizer())
Set parameter LicenseID to value 890341
A JuMP Model
โ mode: DIRECT
โ solver: Gurobi
โ objective_sense: FEASIBILITY_SENSE
โ num_variables: 0
โ num_constraints: 0
โ Names registered in the model: none
julia> set_attribute(model, "LazyConstraints", 1)
Set parameter LazyConstraints to value 1
julia> set_attribute(model, Gurobi.CallbackFunction(), my_callback_fn)
julia> @variable(model, -1.98 <= x <= 1.98, Int)
x
julia> @variable(model, y)
y
julia> @objective(model, Min, y)
y
julia> optimize!(model)
Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (mac64[x86] - Darwin 24.1.0 24B83)
CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Non-default parameters:
LazyConstraints 1
Optimize a model with 0 rows, 2 columns and 0 nonzeros
Model fingerprint: 0xb7a8a0a7
Variable types: 1 continuous, 1 integer (0 binary)
Coefficient statistics:
Matrix range [0e+00, 0e+00]
Objective range [1e+00, 1e+00]
Bounds range [2e+00, 2e+00]
RHS range [0e+00, 0e+00]
Presolve time: 0.00s
Presolved: 0 rows, 2 columns, 0 nonzeros
Variable types: 1 continuous, 1 integer (0 binary)
Root relaxation: unbounded, 0 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 postponed 0 - - - - 0s
0 0 -0.40000 0 - - -0.40000 - - 0s
0 0 0.04444 0 1 - 0.04444 - - 0s
H 0 0 0.4000000 0.04444 88.9% - 0s
0 0 0.04444 0 1 0.40000 0.04444 88.9% - 0s
Cutting planes:
Lazy constraints: 2
Explored 1 nodes (2 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)
Solution count 1: 0.4
Optimal solution found (tolerance 1.00e-04)
Best objective 4.000000000000e-01, best bound 4.000000000000e-01, gap 0.0000%
User-callback calls 90, time in user-callback 0.03 sec
is this the sole standard style to add lazy cuts in callback?
(or ask, is my styleโonly MIPSOLโincorrect)?
Maybe both are okay? I find this python example that only with MIPNODE And this python TSP problem also.
def __call__(self, model, where):
"""Callback entry point: call lazy constraints routine when new
solutions are found. Stop the optimization if there is an exception in
user code."""
if where == GRB.Callback.MIPSOL:
try:
self.eliminate_subtours(model)
except Exception:
logging.exception("Exception occurred in MIPSOL callback")
model.terminate()
I donโt fully understand what this means (from Gurobi Doc)
MIP solutions may be generated outside of a MIP node. Thus, generating lazy constraints is optional when the
wherevalue in the callback function equalsGRB_CB_MIPNODE. To avoid this, we recommend to always check when thewherevalue equalsGRB_CB_MIPSOL.
Maybe he wants to say:
(am I right?)
Update: I spend quite a few time to read Gurobi docs and other books, concluding that the usage in my 23th post is a disciplined one (i.e., not error-prone, to the largest extent acceptable). The features mainly include
MIPSOL.
If there is anyone who want to discuss this, please contact me.
The details of how callbacks work are subtle. I assume that the underlying issue is that your problem is unbounded without the callback. You could argue that it should call a MIPSOL with a feasible point, find out that the problem is bounded, and continue. But they might also argue this is just a documentation issue and that you must also check MIPNODE.
@torressa may want to chime in here.