BB: Support general int vars

master
Alinson S. Xavier 3 years ago
parent 372a92d5d2
commit 9b3a0da5f7
Signed by: isoron
GPG Key ID: 0DA8E4B9E1109DCA

@ -6,27 +6,22 @@ using Printf
function print_progress_header(; detailed_output::Bool)
@printf(
"%8s %9s %9s %13s %13s %13s %9s %8s",
"%8s %9s %9s %13s %13s %9s %6s %13s %6s %-24s %9s %9s %6s %6s",
"time",
"visited",
"processed",
"pending",
"obj",
"primal-bound",
"dual-bound",
"gap",
"lp-iter"
)
if detailed_output
@printf(
" %6s %6s %-24s %6s %6s %6s",
"node",
"obj",
"parent",
"branch-var",
"b-val",
"branch-lb",
"branch-ub",
"depth",
"iinfes"
)
end
println()
flush(stdout)
end
@ -39,40 +34,35 @@ function print_progress(
detailed_output::Bool,
primal_update::Bool,
)::Nothing
prefix = primal_update ? "*" : " "
if (pool.processed % print_interval == 0) || isempty(pool.pending) || primal_update
if isempty(node.branch_vars)
branch_var_name = "---"
branch_lb = "---"
branch_ub = "---"
else
branch_var_name = name(node.mip, last(node.branch_vars))
L = min(24, length(branch_var_name))
branch_var_name = branch_var_name[1:L]
branch_lb = @sprintf("%9.2f", last(node.branch_lb))
branch_ub = @sprintf("%9.2f", last(node.branch_ub))
end
@printf(
"%8.2f %1s%9d %9d %13.6e %13.6e %13.6e %9.2e %8d",
"%8.2f %9d %9d %13.6e %13.6e %9.2e %6d %13.6e %6s %-24s %9s %9s %6d %6d",
time_elapsed,
prefix,
pool.processed,
length(pool.processing) + length(pool.pending),
node.obj * node.mip.sense,
pool.primal_bound * node.mip.sense,
pool.dual_bound * node.mip.sense,
pool.gap,
pool.mip.lp_iterations,
)
if detailed_output
if isempty(node.branch_variables)
branch_var_name = "---"
branch_value = "---"
else
branch_var_name = name(node.mip, last(node.branch_variables))
L = min(24, length(branch_var_name))
branch_var_name = branch_var_name[1:L]
branch_value = @sprintf("%.2f", last(node.branch_values))
end
@printf(
" %6d %6s %-24s %6s %6d %6d",
node.index,
node.obj * node.mip.sense,
node.parent === nothing ? "---" : @sprintf("%d", node.parent.index),
branch_var_name,
branch_value,
length(node.branch_variables),
branch_lb,
branch_ub,
length(node.branch_vars),
length(node.fractional_variables)
)
end
println()
flush(stdout)
end

@ -10,21 +10,36 @@ using MathOptInterface
const MOI = MathOptInterface
function init(constructor)::MIP
return MIP(constructor, Any[nothing for t = 1:nthreads()], Variable[], 1.0, 0)
return MIP(
constructor,
Any[nothing for t = 1:nthreads()],
Variable[],
Float64[],
Float64[],
1.0,
0,
)
end
function read!(mip::MIP, filename::AbstractString)::Nothing
load!(mip, read_from_file(filename))
return
end
function load!(mip::MIP, prototype::JuMP.Model)
@threads for t = 1:nthreads()
model = read_from_file(filename)
set_optimizer(model, mip.constructor)
mip.optimizers[t] = backend(model)
_replace_zero_one!(mip.optimizers[t])
model = Model()
MOI.copy_to(model, backend(prototype))
_replace_zero_one!(backend(model))
if t == 1
_assert_supported(mip.optimizers[t])
mip.binary_variables = _get_binary_variables(mip.optimizers[t])
mip.sense = _get_objective_sense(mip.optimizers[t])
_assert_supported(backend(model))
mip.int_vars, mip.int_vars_lb, mip.int_vars_ub =
_get_int_variables(backend(model))
mip.sense = _get_objective_sense(backend(model))
end
_relax_integrality!(mip.optimizers[t])
_relax_integrality!(backend(model))
set_optimizer(model, mip.constructor)
mip.optimizers[t] = backend(model)
set_silent(model)
end
return
@ -63,9 +78,16 @@ function _get_objective_sense(optimizer::MOI.AbstractOptimizer)::Float64
end
end
_bounds_constraint(v::Variable) =
_interval_index(v::Variable) =
MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Float64}}(v.index)
_lower_bound_index(v::Variable) =
MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(v.index)
_upper_bound_index(v::Variable) =
MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(v.index)
function _replace_zero_one!(optimizer::MOI.AbstractOptimizer)::Nothing
constrs_to_delete = MOI.ConstraintIndex[]
funcs = MOI.VariableIndex[]
@ -85,22 +107,46 @@ function _replace_zero_one!(optimizer::MOI.AbstractOptimizer)::Nothing
return
end
function _get_binary_variables(optimizer::MOI.AbstractOptimizer)::Vector{Variable}
function _get_int_variables(
optimizer::MOI.AbstractOptimizer,
)::Tuple{Vector{Variable},Vector{Float64},Vector{Float64}}
vars = Variable[]
lb = Float64[]
ub = Float64[]
for ci in
MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer}())
func = MOI.get(optimizer, MOI.ConstraintFunction(), ci)
var = Variable(func.value)
MOI.is_valid(optimizer, _bounds_constraint(var)) ||
error("$var is not interval-constrained")
interval = MOI.get(optimizer, MOI.ConstraintSet(), _bounds_constraint(var))
interval.lower == 0.0 || error("$var has lb != 0")
interval.upper == 1.0 || error("$var has ub != 1")
var_lb, var_ub = -Inf, Inf
if MOI.is_valid(optimizer, _interval_index(var))
constr = MOI.get(optimizer, MOI.ConstraintSet(), _interval_index(var))
var_ub = constr.upper
var_lb = constr.lower
else
# If interval constraint is not found, query individual lower/upper bound
# constraints and replace them by an interval constraint.
if MOI.is_valid(optimizer, _lower_bound_index(var))
constr = MOI.get(optimizer, MOI.ConstraintSet(), _lower_bound_index(var))
var_lb = constr.lower
MOI.delete(optimizer, _lower_bound_index(var))
end
if MOI.is_valid(optimizer, _upper_bound_index(var))
constr = MOI.get(optimizer, MOI.ConstraintSet(), _upper_bound_index(var))
var_ub = constr.upper
MOI.delete(optimizer, _upper_bound_index(var))
end
MOI.add_constraint(
optimizer,
var,
MOI.Interval(var_lb, var_ub),
)
end
push!(vars, var)
push!(lb, var_lb)
push!(ub, var_ub)
end
return vars
return vars, lb, ub
end
function _relax_integrality!(optimizer::MOI.AbstractOptimizer)::Nothing
@ -169,14 +215,14 @@ function set_bounds!(
ub::Array{Float64},
)::Nothing
t = threadid()
MOI.delete(mip.optimizers[t], _bounds_constraint.(vars))
funcs = MOI.VariableIndex[]
sets = MOI.Interval[]
for j = 1:length(vars)
push!(funcs, MOI.VariableIndex(vars[j].index))
push!(sets, MOI.Interval(lb[j], ub[j]))
MOI.delete(mip.optimizers[t], _interval_index(vars[j]))
MOI.add_constraint(
mip.optimizers[t],
MOI.VariableIndex(vars[j].index),
MOI.Interval(lb[j], ub[j]),
)
end
MOI.add_constraints(mip.optimizers[t], funcs, sets)
return
end
@ -190,23 +236,23 @@ function name(mip::MIP, var::Variable)::String
return MOI.get(mip.optimizers[t], MOI.VariableName(), MOI.VariableIndex(var.index))
end
convert(::Type{MOI.VariableIndex}, v::Variable) = MOI.VariableIndex(v.index)
# convert(::Type{MOI.VariableIndex}, v::Variable) = MOI.VariableIndex(v.index)
"""
probe(mip::MIP, var)::Tuple{Float64, Float64}
probe(mip::MIP, var, frac, lb, ub)::Tuple{Float64, Float64}
Suppose that the LP relaxation of `mip` has been solved and that `var` holds
a fractional value `f`. This function returns two numbers corresponding,
respectively, to the the optimal values of the LP relaxations having the
constraints var=1 and var=0 enforced. If any branch is infeasible, the optimal
value for that branch is Inf for minimization problems and -Inf for maximization
problems.
constraints `ceil(frac) <= var <= ub` and `lb <= var <= floor(frac)` enforced.
If any branch is infeasible, the optimal value for that branch is Inf for
minimization problems and -Inf for maximization problems.
"""
function probe(mip::MIP, var)::Tuple{Float64,Float64}
set_bounds!(mip, [var], [1.0], [1.0])
function probe(mip::MIP, var::Variable, frac::Float64, lb::Float64, ub::Float64)::Tuple{Float64,Float64}
set_bounds!(mip, [var], [ceil(frac)], [ub])
status_up, obj_up = solve_relaxation!(mip)
set_bounds!(mip, [var], [0.0], [0.0])
set_bounds!(mip, [var], [lb], [floor(frac)])
status_down, obj_down = solve_relaxation!(mip)
set_bounds!(mip, [var], [0.0], [1.0])
set_bounds!(mip, [var], [lb], [ub])
return obj_up * mip.sense, obj_down * mip.sense
end

@ -61,8 +61,8 @@ function offer(
primal_update = false
# Update node.processing and node.processed
pool.processed += 1
if parent_node !== nothing
pool.processed += 1
delete!(pool.processing, parent_node)
end
@ -71,7 +71,7 @@ function offer(
if node.status == :Infeasible
continue
end
if node.obj >= pool.primal_bound
if node.obj >= pool.primal_bound - 1e-6
continue
end
if isempty(node.fractional_variables)
@ -100,7 +100,7 @@ function offer(
if parent_node !== nothing
# Update branching variable history
branch_var = child_nodes[1].branch_variables[end]
branch_var = child_nodes[1].branch_vars[end]
offset = findfirst(isequal(branch_var), parent_node.fractional_variables)
x = parent_node.fractional_values[offset]
obj_change_up = child_nodes[1].obj - parent_node.obj
@ -112,21 +112,21 @@ function offer(
obj_change_down = obj_change_down,
obj_change_up = obj_change_up,
)
# Update global history
pool.history.avg_pseudocost_up =
mean(vh.pseudocost_up for vh in values(pool.var_history))
pool.history.avg_pseudocost_down =
mean(vh.pseudocost_down for vh in values(pool.var_history))
end
# Print progress
for node in child_nodes
print_progress(
pool,
parent_node,
node,
time_elapsed = time_elapsed,
print_interval = print_interval,
detailed_output = detailed_output,
primal_update = primal_update,
primal_update = isfinite(node.obj) && isempty(node.fractional_variables),
)
end
end

@ -6,6 +6,8 @@ using Printf
using Base.Threads
import Base.Threads: @threads, nthreads, threadid
import ..load_data, ..Hdf5Sample
function solve!(
mip::MIP;
time_limit::Float64 = Inf,
@ -14,13 +16,12 @@ function solve!(
print_interval::Int = 5,
initial_primal_bound::Float64 = Inf,
detailed_output::Bool = false,
branch_rule::VariableBranchingRule = HybridBranching(),
branch_rule::VariableBranchingRule = ReliabilityBranching(),
enable_plunging = true,
)::NodePool
time_initial = time()
pool = NodePool(mip = mip)
pool = NodePool(mip=mip)
pool.primal_bound = initial_primal_bound
print_progress_header(detailed_output = detailed_output)
root_node = _create_node(mip)
if isempty(root_node.fractional_variables)
@ -28,9 +29,17 @@ function solve!(
pool.dual_bound = root_node.obj
pool.primal_bound = root_node.obj
return pool
else
print_progress_header(detailed_output=detailed_output)
end
offer(pool, parent_node = nothing, child_nodes = [root_node])
offer(
pool,
parent_node=nothing,
child_nodes=[root_node],
print_interval=print_interval,
detailed_output=detailed_output,
)
@threads for t = 1:nthreads()
child_one, child_zero, suggestions = nothing, nothing, Node[]
while true
@ -40,10 +49,10 @@ function solve!(
end
node = take(
pool,
suggestions = suggestions,
time_remaining = time_limit - time_elapsed,
node_limit = node_limit,
gap_limit = gap_limit,
suggestions=suggestions,
time_remaining=time_limit - time_elapsed,
node_limit=node_limit,
gap_limit=gap_limit,
)
if node == :END
break
@ -52,28 +61,46 @@ function solve!(
continue
else
ids = generate_indices(pool, 2)
branch_variable = find_branching_var(branch_rule, node, pool)
branch_var = find_branching_var(branch_rule, node, pool)
# Find current variable lower and upper bounds
offset = findfirst(isequal(branch_var), mip.int_vars)
var_lb = mip.int_vars_lb[offset]
var_ub = mip.int_vars_ub[offset]
for (offset, v) in enumerate(node.branch_vars)
if v == branch_var
var_lb = max(var_lb, node.branch_lb[offset])
var_ub = min(var_ub, node.branch_ub[offset])
end
end
# Query current fractional value
offset = findfirst(isequal(branch_var), node.fractional_variables)
var_value = node.fractional_values[offset]
child_zero = _create_node(
mip,
index = ids[1],
parent = node,
branch_variable = branch_variable,
branch_value = 0.0,
index=ids[2],
parent=node,
branch_var=branch_var,
branch_var_lb=var_lb,
branch_var_ub=floor(var_value),
)
child_one = _create_node(
mip,
index = ids[2],
parent = node,
branch_variable = branch_variable,
branch_value = 1.0,
index=ids[1],
parent=node,
branch_var=branch_var,
branch_var_lb=ceil(var_value),
branch_var_ub=var_ub,
)
offer(
pool,
parent_node = node,
child_nodes = [child_one, child_zero],
time_elapsed = time_elapsed,
print_interval = print_interval,
detailed_output = detailed_output,
parent_node=node,
child_nodes=[child_one, child_zero],
time_elapsed=time_elapsed,
print_interval=print_interval,
detailed_output=detailed_output,
)
end
end
@ -83,44 +110,95 @@ end
function _create_node(
mip;
index::Int = 0,
parent::Union{Nothing,Node} = nothing,
branch_variable::Union{Nothing,Variable} = nothing,
branch_value::Union{Nothing,Float64} = nothing,
index::Int=0,
parent::Union{Nothing,Node}=nothing,
branch_var::Union{Nothing,Variable}=nothing,
branch_var_lb::Union{Nothing,Float64}=nothing,
branch_var_ub::Union{Nothing,Float64}=nothing
)::Node
if parent === nothing
branch_variables = Variable[]
branch_values = Float64[]
branch_vars = Variable[]
branch_lb = Float64[]
branch_ub = Float64[]
depth = 1
else
branch_variables = [parent.branch_variables; branch_variable]
branch_values = [parent.branch_values; branch_value]
branch_vars = [parent.branch_vars; branch_var]
branch_lb = [parent.branch_lb; branch_var_lb]
branch_ub = [parent.branch_ub; branch_var_ub]
depth = parent.depth + 1
end
set_bounds!(mip, branch_variables, branch_values, branch_values)
set_bounds!(mip, branch_vars, branch_lb, branch_ub)
status, obj = solve_relaxation!(mip)
if status == :Optimal
vals = values(mip, mip.binary_variables)
vals = values(mip, mip.int_vars)
fractional_indices =
[j for j in 1:length(mip.binary_variables) if 1e-6 < vals[j] < 1 - 1e-6]
[j for j in 1:length(mip.int_vars) if 1e-6 < vals[j] - floor(vals[j]) < 1 - 1e-6]
fractional_values = vals[fractional_indices]
fractional_variables = mip.binary_variables[fractional_indices]
fractional_variables = mip.int_vars[fractional_indices]
else
fractional_variables = Variable[]
fractional_values = Float64[]
end
n_branch = length(branch_variables)
set_bounds!(mip, branch_variables, zeros(n_branch), ones(n_branch))
n_branch = length(branch_vars)
set_bounds!(mip, branch_vars, zeros(n_branch), ones(n_branch))
return Node(
mip,
index,
depth,
obj,
status,
branch_variables,
branch_values,
branch_vars,
branch_lb,
branch_ub,
fractional_variables,
fractional_values,
parent,
)
end
function solve!(
optimizer,
filename::String;
time_limit::Float64=Inf,
node_limit::Int=typemax(Int),
gap_limit::Float64=1e-4,
print_interval::Int=5,
detailed_output::Bool=false,
branch_rule::VariableBranchingRule=ReliabilityBranching()
)::NodePool
model = read_from_file("$filename.mps.gz")
mip = init(optimizer)
load!(mip, model)
h5 = Hdf5Sample("$filename.h5")
primal_bound = h5.get_scalar("mip_obj_value")
nvars = length(h5.get_array("static_var_names"))
pool = solve!(
mip;
initial_primal_bound=primal_bound,
time_limit,
node_limit,
gap_limit,
print_interval,
detailed_output,
branch_rule
)
pseudocost_up = [NaN for _ = 1:nvars]
pseudocost_down = [NaN for _ = 1:nvars]
priorities = [0.0 for _ in 1:nvars]
for (var, var_hist) in pool.var_history
pseudocost_up[var.index] = var_hist.pseudocost_up
pseudocost_down[var.index] = var_hist.pseudocost_down
x = mean(var_hist.fractional_values)
f_up = x - floor(x)
f_down = ceil(x) - x
priorities[var.index] = var_hist.pseudocost_up * f_up * var_hist.pseudocost_down * f_down
end
h5.put_array("bb_var_pseudocost_up", pseudocost_up)
h5.put_array("bb_var_pseudocost_down", pseudocost_down)
h5.put_array("bb_var_priority", priorities)
return pool
end

@ -12,7 +12,9 @@ end
mutable struct MIP
constructor::Any
optimizers::Vector
binary_variables::Vector{Variable}
int_vars::Vector{Variable}
int_vars_lb::Vector{Float64}
int_vars_ub::Vector{Float64}
sense::Float64
lp_iterations::Int64
end
@ -23,8 +25,9 @@ struct Node
depth::Int
obj::Float64
status::Symbol
branch_variables::Array{Variable}
branch_values::Array{Float64}
branch_vars::Array{Variable}
branch_lb::Array{Float64}
branch_ub::Array{Float64}
fractional_variables::Array{Variable}
fractional_values::Array{Float64}
parent::Union{Nothing,Node}

@ -33,7 +33,7 @@ function find_branching_var(rule::StrongBranching, node::Node, pool::NodePool)::
sorted_vars = node.fractional_variables[σ]
_strong_branch_start(node)
no_improv_count, call_count = 0, 0
max_score, max_var = pseudocost_scores[σ[1]], sorted_vars[1]
max_score, max_var = (-Inf, -Inf), sorted_vars[1]
for (i, var) in enumerate(sorted_vars)
call_count += 1
score = _strong_branch_score(
@ -43,6 +43,7 @@ function find_branching_var(rule::StrongBranching, node::Node, pool::NodePool)::
x = node.fractional_values[σ[i]],
side_effect = rule.side_effect,
)
# @show name(node.mip, var), round(score[1], digits=2)
if score > max_score
max_score, max_var = score, var
no_improv_count = 0
@ -63,9 +64,21 @@ function _strong_branch_score(;
x::Float64,
side_effect::Bool,
)::Tuple{Float64,Int}
# Find current variable lower and upper bounds
offset = findfirst(isequal(var), node.mip.int_vars)
var_lb = node.mip.int_vars_lb[offset]
var_ub = node.mip.int_vars_ub[offset]
for (offset, v) in enumerate(node.branch_vars)
if v == var
var_lb = max(var_lb, node.branch_lb[offset])
var_ub = min(var_ub, node.branch_ub[offset])
end
end
obj_up, obj_down = 0, 0
try
obj_up, obj_down = probe(node.mip, var)
obj_up, obj_down = probe(node.mip, var, x, var_lb, var_ub)
catch
@warn "strong branch error" var = var
end
@ -84,10 +97,9 @@ function _strong_branch_score(;
end
function _strong_branch_start(node::Node)
set_bounds!(node.mip, node.branch_variables, node.branch_values, node.branch_values)
set_bounds!(node.mip, node.branch_vars, node.branch_lb, node.branch_ub)
end
function _strong_branch_end(node::Node)
nbranch = length(node.branch_variables)
set_bounds!(node.mip, node.branch_variables, zeros(nbranch), ones(nbranch))
set_bounds!(node.mip, node.mip.int_vars, node.mip.int_vars_lb, node.mip.int_vars_ub)
end

@ -18,56 +18,59 @@ function runtests(optimizer_name, optimizer; large = true)
BB.read!(mip, filename)
@test mip.sense == 1.0
@test length(mip.binary_variables) == 56
@test length(mip.int_vars) == 56
status, obj = BB.solve_relaxation!(mip)
@test status == :Optimal
@test round(obj, digits = 6) == 62.637280
@test BB.name(mip, mip.binary_variables[1]) == "xab"
@test BB.name(mip, mip.binary_variables[2]) == "xac"
@test BB.name(mip, mip.binary_variables[3]) == "xad"
@test BB.name(mip, mip.int_vars[1]) == "xab"
@test BB.name(mip, mip.int_vars[2]) == "xac"
@test BB.name(mip, mip.int_vars[3]) == "xad"
vals = BB.values(mip, mip.binary_variables)
@test mip.int_vars_lb[1] == 0.0
@test mip.int_vars_ub[1] == 1.0
vals = BB.values(mip, mip.int_vars)
@test round(vals[1], digits = 6) == 0.046933
@test round(vals[2], digits = 6) == 0.000841
@test round(vals[3], digits = 6) == 0.248696
# Probe (up and down are feasible)
probe_up, probe_down = BB.probe(mip, mip.binary_variables[1])
probe_up, probe_down = BB.probe(mip, mip.int_vars[1], 0.5, 0.0, 1.0)
@test round(probe_down, digits = 6) == 62.690000
@test round(probe_up, digits = 6) == 62.714100
# Fix one variable to zero
BB.set_bounds!(mip, mip.binary_variables[1:1], [0.0], [0.0])
BB.set_bounds!(mip, mip.int_vars[1:1], [0.0], [0.0])
status, obj = BB.solve_relaxation!(mip)
@test status == :Optimal
@test round(obj, digits = 6) == 62.690000
# Fix one variable to one and another variable variable to zero
BB.set_bounds!(mip, mip.binary_variables[1:2], [1.0, 0.0], [1.0, 0.0])
BB.set_bounds!(mip, mip.int_vars[1:2], [1.0, 0.0], [1.0, 0.0])
status, obj = BB.solve_relaxation!(mip)
@test status == :Optimal
@test round(obj, digits = 6) == 62.714777
# Probe (up is infeasible, down is feasible)
BB.set_bounds!(mip, mip.binary_variables[1:3], [1.0, 1.0, 0.0], [1.0, 1.0, 1.0])
BB.set_bounds!(mip, mip.int_vars[1:3], [1.0, 1.0, 0.0], [1.0, 1.0, 1.0])
status, obj = BB.solve_relaxation!(mip)
@test status == :Optimal
probe_up, probe_down = BB.probe(mip, mip.binary_variables[3])
probe_up, probe_down = BB.probe(mip, mip.int_vars[3], 0.5, 0.0, 1.0)
@test round(probe_up, digits = 6) == Inf
@test round(probe_down, digits = 6) == 63.073992
# Fix all binary variables to one, making problem infeasible
N = length(mip.binary_variables)
BB.set_bounds!(mip, mip.binary_variables, ones(N), ones(N))
N = length(mip.int_vars)
BB.set_bounds!(mip, mip.int_vars, ones(N), ones(N))
status, obj = BB.solve_relaxation!(mip)
@test status == :Infeasible
@test obj == Inf
# Restore original problem
N = length(mip.binary_variables)
BB.set_bounds!(mip, mip.binary_variables, zeros(N), ones(N))
N = length(mip.int_vars)
BB.set_bounds!(mip, mip.int_vars, zeros(N), ones(N))
status, obj = BB.solve_relaxation!(mip)
@test status == :Optimal
@test round(obj, digits = 6) == 62.637280

Loading…
Cancel
Save