diff --git a/Manifest.toml b/Manifest.toml index 5aa38b6..058a1a3 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -126,9 +126,9 @@ version = "1.1.1" [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677" +git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.9" +version = "0.18.10" [[DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" diff --git a/Project.toml b/Project.toml index aed7353..bc90f42 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d" Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" @@ -18,7 +19,9 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] diff --git a/src/MIPLearn.jl b/src/MIPLearn.jl index 61d164d..d721aa3 100644 --- a/src/MIPLearn.jl +++ b/src/MIPLearn.jl @@ -33,6 +33,8 @@ include("solvers/macros.jl") include("utils/benchmark.jl") include("utils/parse.jl") +include("bb/BB.jl") + function __init__() copy!(miplearn, pyimport("miplearn")) copy!(traceback, pyimport("traceback")) diff --git a/src/bb/BB.jl b/src/bb/BB.jl new file mode 100644 index 0000000..4c84463 --- /dev/null +++ b/src/bb/BB.jl @@ -0,0 +1,22 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +module BB + +frac(x) = x - floor(x) + +include("structs.jl") + +include("nodepool.jl") +include("optimize.jl") +include("log.jl") +include("lp.jl") +include("varbranch/hybrid.jl") +include("varbranch/infeasibility.jl") +include("varbranch/pseudocost.jl") +include("varbranch/random.jl") +include("varbranch/reliability.jl") +include("varbranch/strong.jl") + +end # module diff --git a/src/bb/log.jl b/src/bb/log.jl new file mode 100644 index 0000000..a3cb2f2 --- /dev/null +++ b/src/bb/log.jl @@ -0,0 +1,79 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using Printf + +function print_progress_header(; detailed_output::Bool) + @printf( + "%8s %9s %9s %13s %13s %13s %9s %8s", + "time", + "visited", + "pending", + "obj", + "primal-bound", + "dual-bound", + "gap", + "lp-iter" + ) + if detailed_output + @printf( + " %6s %6s %-24s %6s %6s %6s", + "node", + "parent", + "branch-var", + "b-val", + "depth", + "iinfes" + ) + end + println() + flush(stdout) +end + +function print_progress( + pool::NodePool, + node::Node; + time_elapsed::Float64, + print_interval::Int, + detailed_output::Bool, + primal_update::Bool, +)::Nothing + prefix = primal_update ? "*" : " " + if (pool.processed % print_interval == 0) || isempty(pool.pending) || primal_update + @printf( + "%8.2f %1s%9d %9d %13.6e %13.6e %13.6e %9.2e %8d", + 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.parent === nothing ? "---" : @sprintf("%d", node.parent.index), + branch_var_name, + branch_value, + length(node.branch_variables), + length(node.fractional_variables) + ) + end + println() + flush(stdout) + end +end diff --git a/src/bb/lp.jl b/src/bb/lp.jl new file mode 100644 index 0000000..7ed1a36 --- /dev/null +++ b/src/bb/lp.jl @@ -0,0 +1,212 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +import Base: values, convert +using Base.Threads +import Base.Threads: @threads, nthreads, threadid +using JuMP +using MathOptInterface +const MOI = MathOptInterface + +function init(constructor)::MIP + return MIP(constructor, Any[nothing for t = 1:nthreads()], Variable[], 1.0, 0) +end + +function read!(mip::MIP, filename::AbstractString)::Nothing + @threads for t = 1:nthreads() + model = read_from_file(filename) + mip.optimizers[t] = backend(model) + _replace_zero_one!(mip.optimizers[t]) + 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]) + end + _relax_integrality!(mip.optimizers[t]) + set_optimizer(model, mip.constructor) + set_silent(model) + end + return +end + +function _assert_supported(optimizer::MOI.AbstractOptimizer)::Nothing + types = MOI.get(optimizer, MOI.ListOfConstraints()) + for (F, S) in types + _assert_supported(F, S) + end +end + +function _assert_supported(F::DataType, S::DataType)::Nothing + if F in [MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable] && S in [ + MOI.LessThan{Float64}, + MOI.GreaterThan{Float64}, + MOI.EqualTo{Float64}, + MOI.Interval{Float64}, + ] + return + end + if F in [MOI.SingleVariable] && S in [MOI.Integer, MOI.ZeroOne] + return + end + error("MOI constraint not supported: $F in $S") +end + +function _get_objective_sense(optimizer::MOI.AbstractOptimizer)::Float64 + sense = MOI.get(optimizer, MOI.ObjectiveSense()) + if sense == MOI.MIN_SENSE + return 1.0 + elseif sense == MOI.MAX_SENSE + return -1.0 + else + error("objective sense not supported: $sense") + end +end + +_bounds_constraint(v::Variable) = + MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{Float64}}(v.index) + +function _replace_zero_one!(optimizer::MOI.AbstractOptimizer)::Nothing + constrs_to_delete = MOI.ConstraintIndex[] + funcs = MOI.SingleVariable[] + sets = Union{MOI.Interval,MOI.Integer}[] + for ci in + MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.ZeroOne}()) + func = MOI.get(optimizer, MOI.ConstraintFunction(), ci) + var = func.variable + push!(constrs_to_delete, ci) + push!(funcs, MOI.SingleVariable(var)) + push!(funcs, MOI.SingleVariable(var)) + push!(sets, MOI.Interval{Float64}(0.0, 1.0)) + push!(sets, MOI.Integer()) + end + MOI.delete(optimizer, constrs_to_delete) + MOI.add_constraints(optimizer, funcs, sets) + return +end + +function _get_binary_variables(optimizer::MOI.AbstractOptimizer)::Vector{Variable} + vars = Variable[] + for ci in + MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer}()) + func = MOI.get(optimizer, MOI.ConstraintFunction(), ci) + var = Variable(func.variable.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") + + push!(vars, var) + end + return vars +end + +function _relax_integrality!(optimizer::MOI.AbstractOptimizer)::Nothing + indices = + MOI.get(optimizer, MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer}()) + MOI.delete(optimizer, indices) +end + +""" + solve_relaxation(mip::MIP)::Tuple{Symbol, Float64} + +Solve the linear relaxation of `mip` and returns a tuple containing the +solution status (either :Optimal or :Infeasible) and the optimal objective +value. If the problem is infeasible, the optimal value is Inf for minimization +problems and -Inf for maximization problems.. +""" +function solve_relaxation!(mip::MIP)::Tuple{Symbol,Float64} + t = threadid() + MOI.optimize!(mip.optimizers[t]) + try + mip.lp_iterations += MOI.get(mip.optimizers[t], MOI.SimplexIterations()) + catch + # ignore + end + status = MOI.get(mip.optimizers[t], MOI.TerminationStatus()) + if status == MOI.OPTIMAL + obj = MOI.get(mip.optimizers[t], MOI.ObjectiveValue()) + return :Optimal, obj * mip.sense + elseif status in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED] + return :Infeasible, Inf * mip.sense + end + error("unknown status: $status") +end + +""" + values(mip::MIP, vars::Vector{Variable})::Array{Float64} + +Returns a vector `vals` which describes the current primal values for the +decision variables `vars`. More specifically, `vals[j]` is the current +primal value of `vars[j]`. +""" +function values(mip::MIP, vars::Vector{Variable})::Array{Float64} + return MOI.get( + mip.optimizers[threadid()], + MOI.VariablePrimal(), + convert.(MOI.VariableIndex, vars), + ) +end + +values(mip::MIP) = + values(mip, MOI.get(mip.optimizers[threadid()], MOI.ListOfVariableIndices())) + +""" + set_bounds!(mip::MIP, + vars::Vector{Variable}, + lb::Array{Float64}, + ub::Array{Float64}) + +Modify the bounds of the given variables. More specifically, sets +upper and lower bounds of `vars[j]` to `ub[j]` and `lb[j]`, respectively. +""" +function set_bounds!( + mip::MIP, + vars::Vector{Variable}, + lb::Array{Float64}, + ub::Array{Float64}, +)::Nothing + t = threadid() + MOI.delete(mip.optimizers[t], _bounds_constraint.(vars)) + funcs = MOI.SingleVariable[] + sets = MOI.Interval[] + for j = 1:length(vars) + push!(funcs, MOI.SingleVariable(vars[j])) + push!(sets, MOI.Interval(lb[j], ub[j])) + end + MOI.add_constraints(mip.optimizers[t], funcs, sets) + return +end + +""" + name(mip::MIP, var::Variable)::String + +Return the name of the decision variable `var`. +""" +function name(mip::MIP, var::Variable)::String + t = threadid() + return MOI.get(mip.optimizers[t], MOI.VariableName(), convert(MOI.VariableIndex, var)) +end + +convert(::Type{MOI.VariableIndex}, v::Variable) = MOI.VariableIndex(v.index) + +""" + probe(mip::MIP, var)::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. +""" +function probe(mip::MIP, var)::Tuple{Float64,Float64} + set_bounds!(mip, [var], [1.0], [1.0]) + status_up, obj_up = solve_relaxation!(mip) + set_bounds!(mip, [var], [0.0], [0.0]) + status_down, obj_down = solve_relaxation!(mip) + set_bounds!(mip, [var], [0.0], [1.0]) + return obj_up * mip.sense, obj_down * mip.sense +end diff --git a/src/bb/nodepool.jl b/src/bb/nodepool.jl new file mode 100644 index 0000000..dae435f --- /dev/null +++ b/src/bb/nodepool.jl @@ -0,0 +1,188 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using Statistics +using DataStructures +import Base.Threads: threadid + +function take( + pool::NodePool; + suggestions::Array{Node} = [], + time_remaining::Float64, + gap_limit::Float64, + node_limit::Int, +)::Union{Symbol,Node} + t = threadid() + lock(pool.lock) do + n_processing = length(pool.processing) + if ( + (pool.gap < gap_limit) || + (n_processing + pool.processed >= node_limit) || + (time_remaining < 0) + ) + return :END + end + if isempty(pool.pending) + if isempty(pool.processing) + return :END + else + return :WAIT + end + else + # If one of the suggested nodes is still pending, return it. + # This is known in the literature as plunging. + for s in suggestions + if s in keys(pool.pending) + delete!(pool.pending, s) + pool.processing[s] = s.obj + return s + end + end + + # If all suggestions have already been processed + # or pruned, find another node based on best bound. + node = dequeue!(pool.pending) + pool.processing[node] = node.obj + return node + end + end +end + +function offer( + pool::NodePool; + parent_node::Union{Nothing,Node}, + child_nodes::Vector{Node}, + time_elapsed::Float64 = 0.0, + print_interval::Int = 100, + detailed_output::Bool = false, +)::Nothing + lock(pool.lock) do + primal_update = false + + # Update node.processing and node.processed + pool.processed += 1 + if parent_node !== nothing + delete!(pool.processing, parent_node) + end + + # Queue child nodes + for node in child_nodes + if node.status == :Infeasible + continue + end + if node.obj >= pool.primal_bound + continue + end + if isempty(node.fractional_variables) + pool.primal_bound = min(pool.primal_bound, node.obj) + primal_update = true + continue + end + pool.pending[node] = node.obj + end + + # Update dual bound + pool.dual_bound = pool.primal_bound + if !isempty(pool.pending) + pool.dual_bound = min(pool.dual_bound, peek(pool.pending)[2]) + end + if !isempty(pool.processing) + pool.dual_bound = min(pool.dual_bound, peek(pool.processing)[2]) + end + + # Update gap + if pool.primal_bound == pool.dual_bound + pool.gap = 0 + else + pool.gap = abs((pool.primal_bound - pool.dual_bound) / pool.primal_bound) + end + + if parent_node !== nothing + # Update branching variable history + branch_var = child_nodes[1].branch_variables[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 + obj_change_down = child_nodes[2].obj - parent_node.obj + _update_var_history( + pool = pool, + var = branch_var, + x = x, + 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)) + + # Print progress + print_progress( + pool, + parent_node, + time_elapsed = time_elapsed, + print_interval = print_interval, + detailed_output = detailed_output, + primal_update = primal_update, + ) + end + end + return +end + +function _update_var_history(; + pool::NodePool, + var::Variable, + x::Float64, + obj_change_down::Float64, + obj_change_up::Float64, +)::Nothing + # Create new history entry + if var ∉ keys(pool.var_history) + pool.var_history[var] = VariableHistory() + end + varhist = pool.var_history[var] + + # Push fractional value + push!(varhist.fractional_values, x) + + # Push objective value changes + push!(varhist.obj_change_up, obj_change_up) + push!(varhist.obj_change_down, obj_change_down) + + # Push objective change ratios + f_up = x - floor(x) + f_down = ceil(x) - x + if isfinite(obj_change_up) + push!(varhist.obj_ratio_up, obj_change_up / f_up) + end + if isfinite(obj_change_down) + push!(varhist.obj_ratio_down, obj_change_down / f_down) + end + + # Update variable pseudocosts + varhist.pseudocost_up = 0.0 + varhist.pseudocost_down = 0.0 + if !isempty(varhist.obj_ratio_up) + varhist.pseudocost_up = sum(varhist.obj_ratio_up) / length(varhist.obj_ratio_up) + end + if !isempty(varhist.obj_ratio_down) + varhist.pseudocost_down = + sum(varhist.obj_ratio_down) / length(varhist.obj_ratio_down) + end + return +end + +function generate_indices(pool::NodePool, n::Int)::Vector{Int} + lock(pool.lock) do + result = Int[] + for i = 1:n + push!(result, pool.next_index) + pool.next_index += 1 + end + return result + end +end diff --git a/src/bb/optimize.jl b/src/bb/optimize.jl new file mode 100644 index 0000000..ed5cbd2 --- /dev/null +++ b/src/bb/optimize.jl @@ -0,0 +1,125 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using Printf +using Base.Threads +import Base.Threads: @threads, nthreads, threadid + +function solve!( + mip::MIP; + time_limit::Float64 = Inf, + node_limit::Int = typemax(Int), + gap_limit::Float64 = 1e-4, + print_interval::Int = 5, + initial_primal_bound::Float64 = Inf, + detailed_output::Bool = false, + branch_rule::VariableBranchingRule = HybridBranching(), +)::NodePool + time_initial = time() + 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) + println("root relaxation is integer feasible") + pool.dual_bound = root_node.obj + pool.primal_bound = root_node.obj + return pool + end + + offer(pool, parent_node = nothing, child_nodes = [root_node]) + @threads for t = 1:nthreads() + child_one, child_zero, suggestions = nothing, nothing, Node[] + while true + time_elapsed = time() - time_initial + if child_one !== nothing + suggestions = Node[child_one, child_zero] + end + node = take( + pool, + suggestions = suggestions, + time_remaining = time_limit - time_elapsed, + node_limit = node_limit, + gap_limit = gap_limit, + ) + if node == :END + break + elseif node == :WAIT + sleep(0.1) + continue + else + ids = generate_indices(pool, 2) + branch_variable = find_branching_var(branch_rule, node, pool) + child_zero = _create_node( + mip, + index = ids[1], + parent = node, + branch_variable = branch_variable, + branch_value = 0.0, + ) + child_one = _create_node( + mip, + index = ids[2], + parent = node, + branch_variable = branch_variable, + branch_value = 1.0, + ) + offer( + pool, + parent_node = node, + child_nodes = [child_one, child_zero], + time_elapsed = time_elapsed, + print_interval = print_interval, + detailed_output = detailed_output, + ) + end + end + end + return pool +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, +)::Node + if parent === nothing + branch_variables = Variable[] + branch_values = Float64[] + depth = 1 + else + branch_variables = [parent.branch_variables; branch_variable] + branch_values = [parent.branch_values; branch_value] + depth = parent.depth + 1 + end + set_bounds!(mip, branch_variables, branch_values, branch_values) + status, obj = solve_relaxation!(mip) + if status == :Optimal + vals = values(mip, mip.binary_variables) + fractional_indices = + [j for j in 1:length(mip.binary_variables) if 1e-6 < vals[j] < 1 - 1e-6] + fractional_values = vals[fractional_indices] + fractional_variables = mip.binary_variables[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)) + return Node( + mip, + index, + depth, + obj, + status, + branch_variables, + branch_values, + fractional_variables, + fractional_values, + parent, + ) +end diff --git a/src/bb/structs.jl b/src/bb/structs.jl new file mode 100644 index 0000000..36cc158 --- /dev/null +++ b/src/bb/structs.jl @@ -0,0 +1,70 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using DataStructures +abstract type VariableBranchingRule end + +struct Variable + index::Any +end + +mutable struct MIP + constructor + optimizers::Vector + binary_variables::Vector{Variable} + sense::Float64 + lp_iterations::Int64 +end + +struct Node + mip::MIP + index::Int + depth::Int + obj::Float64 + status::Symbol + branch_variables::Array{Variable} + branch_values::Array{Float64} + fractional_variables::Array{Variable} + fractional_values::Array{Float64} + parent::Union{Nothing,Node} +end + +Base.@kwdef mutable struct History + avg_pseudocost_up::Float64 = 1.0 + avg_pseudocost_down::Float64 = 1.0 +end + +mutable struct VariableHistory + fractional_values::Array{Float64} + obj_change_up::Array{Float64} + obj_change_down::Array{Float64} + obj_ratio_up::Array{Float64} + obj_ratio_down::Array{Float64} + pseudocost_up::Float64 + pseudocost_down::Float64 + + VariableHistory() = new( + Float64[], # fractional_values + Float64[], # obj_change_up + Float64[], # obj_change_down + Float64[], # obj_ratio_up + Float64[], # obj_ratio_up + 0.0, # pseudocost_up + 0.0, # pseudocost_down + ) +end + +Base.@kwdef mutable struct NodePool + mip::MIP + pending::PriorityQueue{Node,Float64} = PriorityQueue{Node,Float64}() + processing::PriorityQueue{Node,Float64} = PriorityQueue{Node,Float64}() + processed::Int = 0 + next_index::Int = 1 + lock::ReentrantLock = ReentrantLock() + primal_bound::Float64 = Inf + dual_bound::Float64 = Inf + gap::Float64 = Inf + history::History = History() + var_history::Dict{Variable,VariableHistory} = Dict() +end diff --git a/src/bb/varbranch/hybrid.jl b/src/bb/varbranch/hybrid.jl new file mode 100644 index 0000000..fa5d959 --- /dev/null +++ b/src/bb/varbranch/hybrid.jl @@ -0,0 +1,31 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +""" + HybridBranching(depth_cutoff::Int, + shallow_rule::VariableBranchingRule, + deep_rule::::VariableBranchingRule) + +Branching strategy that switches between two variable branching strategies, +according to the depth of the node. + +More specifically, if `node.depth <= depth_cutoff`, then `shallow_rule` is +applied. Otherwise, `deep_rule` is applied. +""" +mutable struct HybridBranching <: VariableBranchingRule + depth_cutoff::Int + shallow_rule::VariableBranchingRule + deep_rule::VariableBranchingRule +end + +HybridBranching(depth_cutoff::Int = 10) = + HybridBranching(depth_cutoff, StrongBranching(), PseudocostBranching()) + +function find_branching_var(rule::HybridBranching, node::Node, pool::NodePool)::Variable + if node.depth <= rule.depth_cutoff + return find_branching_var(rule.shallow_rule, node, pool) + else + return find_branching_var(rule.deep_rule, node, pool) + end +end diff --git a/src/bb/varbranch/infeasibility.jl b/src/bb/varbranch/infeasibility.jl new file mode 100644 index 0000000..2c99ed8 --- /dev/null +++ b/src/bb/varbranch/infeasibility.jl @@ -0,0 +1,54 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +""" + FirstInfeasibleBranching() + +Branching rule that always selects the first fractional variable. +""" +struct FirstInfeasibleBranching <: VariableBranchingRule end + +function find_branching_var( + rule::FirstInfeasibleBranching, + node::Node, + pool::NodePool, +)::Variable + return node.fractional_variables[1] +end + +""" + LeastInfeasibleBranching() + +Branching strategy that select the fractional variable whose value is the closest +to an integral value. +""" +struct LeastInfeasibleBranching <: VariableBranchingRule end + +function find_branching_var( + rule::LeastInfeasibleBranching, + node::Node, + pool::NodePool, +)::Variable + scores = [max(v - floor(v), ceil(v) - v) for v in node.fractional_values] + _, max_offset = findmax(scores) + return node.fractional_variables[max_offset] +end + +""" + MostInfeasibleBranching() + +Branching strategy that selects the fractional variable whose value is closest +to 1/2. +""" +struct MostInfeasibleBranching <: VariableBranchingRule end + +function find_branching_var( + rule::MostInfeasibleBranching, + node::Node, + pool::NodePool, +)::Variable + scores = [min(v - floor(v), ceil(v) - v) for v in node.fractional_values] + _, max_offset = findmax(scores) + return node.fractional_variables[max_offset] +end diff --git a/src/bb/varbranch/pseudocost.jl b/src/bb/varbranch/pseudocost.jl new file mode 100644 index 0000000..8b3cff0 --- /dev/null +++ b/src/bb/varbranch/pseudocost.jl @@ -0,0 +1,45 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +""" + PseudocostBranching() + +Branching strategy that uses historical changes in objective value to estimate +strong branching scores at lower computational cost. +""" +struct PseudocostBranching <: VariableBranchingRule end + +function find_branching_var(rule::PseudocostBranching, node::Node, pool::NodePool)::Variable + scores = [ + _pseudocost_score( + node, + pool, + node.fractional_variables[j], + node.fractional_values[j], + ) for j = 1:length(node.fractional_variables) + ] + _, max_offset = findmax(scores) + return node.fractional_variables[max_offset] +end + +function _pseudocost_score( + node::Node, + pool::NodePool, + var::Variable, + x::Float64, +)::Tuple{Float64,Int} + f_up = x - floor(x) + f_down = ceil(x) - x + pseudo_up = pool.history.avg_pseudocost_up * f_up + pseudo_down = pool.history.avg_pseudocost_down * f_down + if var in keys(pool.var_history) + if isfinite(pool.var_history[var].pseudocost_up) + pseudo_up = pool.var_history[var].pseudocost_up * f_up + end + if isfinite(pool.var_history[var].pseudocost_down) + pseudo_down = pool.var_history[var].pseudocost_down * f_down + end + end + return (pseudo_up * f_up * pseudo_down * f_down, var.index) +end diff --git a/src/bb/varbranch/random.jl b/src/bb/varbranch/random.jl new file mode 100644 index 0000000..fd09c30 --- /dev/null +++ b/src/bb/varbranch/random.jl @@ -0,0 +1,17 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using Random + + +""" + RandomBranching() + +Branching strategy that picks a fractional variable randomly. +""" +struct RandomBranching <: VariableBranchingRule end + +function find_branching_var(rule::RandomBranching, node::Node, pool::NodePool)::Variable + return shuffle(node.fractional_variables)[1] +end diff --git a/src/bb/varbranch/reliability.jl b/src/bb/varbranch/reliability.jl new file mode 100644 index 0000000..a664459 --- /dev/null +++ b/src/bb/varbranch/reliability.jl @@ -0,0 +1,78 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +""" + ReliabilityBranching + +Branching strategy that uses pseudocosts if there are enough observations +to make an accurate prediction of strong branching scores, or runs the +actual strong branching routine if not enough observations are available. +""" +Base.@kwdef mutable struct ReliabilityBranching <: VariableBranchingRule + min_samples::Int = 8 + max_sb_calls::Int = 100 + look_ahead::Int = 10 + n_sb_calls::Int = 0 + side_effect::Bool = true +end + +function find_branching_var( + rule::ReliabilityBranching, + node::Node, + pool::NodePool, +)::Variable + nfrac = length(node.fractional_variables) + pseudocost_scores = [ + _pseudocost_score( + node, + pool, + node.fractional_variables[j], + node.fractional_values[j], + ) for j = 1:nfrac + ] + σ = sortperm(pseudocost_scores, rev = true) + sorted_vars = node.fractional_variables[σ] + _strong_branch_start(node) + no_improv_count, n_sb_calls = 0, 0 + max_score, max_var = pseudocost_scores[σ[1]], sorted_vars[1] + for (i, var) in enumerate(sorted_vars) + use_strong_branch = true + if n_sb_calls >= rule.max_sb_calls + use_strong_branch = false + else + if var in keys(pool.var_history) + varhist = pool.var_history[var] + hlength = min( + length(varhist.obj_ratio_up), + length(varhist.obj_ratio_down), + ) + if hlength >= rule.min_samples + use_strong_branch = false + end + end + end + if use_strong_branch + n_sb_calls += 1 + rule.n_sb_calls += 1 + score = _strong_branch_score( + node = node, + pool = pool, + var = var, + x = node.fractional_values[σ[i]], + side_effect = rule.side_effect, + ) + else + score = pseudocost_scores[σ[i]] + end + if score > max_score + max_score, max_var = score, var + no_improv_count = 0 + else + no_improv_count += 1 + end + no_improv_count <= rule.look_ahead || break + end + _strong_branch_end(node) + return max_var +end diff --git a/src/bb/varbranch/strong.jl b/src/bb/varbranch/strong.jl new file mode 100644 index 0000000..f4fae31 --- /dev/null +++ b/src/bb/varbranch/strong.jl @@ -0,0 +1,93 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using Random + +""" + StrongBranching(look_ahead::Int, max_calls::Int) + +Branching strategy that selects a subset of fractional variables +as candidates (according to pseudocosts) the solves two linear +programming problems for each candidate. + + +""" +Base.@kwdef struct StrongBranching <: VariableBranchingRule + look_ahead::Int = 10 + max_calls::Int = 100 + side_effect::Bool = true +end + +function find_branching_var(rule::StrongBranching, node::Node, pool::NodePool)::Variable + nfrac = length(node.fractional_variables) + pseudocost_scores = [ + _pseudocost_score( + node, + pool, + node.fractional_variables[j], + node.fractional_values[j], + ) for j = 1:nfrac + ] + σ = sortperm(pseudocost_scores, rev = true) + 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] + for (i, var) in enumerate(sorted_vars) + call_count += 1 + score = _strong_branch_score( + node = node, + pool = pool, + var = var, + x = node.fractional_values[σ[i]], + side_effect = rule.side_effect, + ) + if score > max_score + max_score, max_var = score, var + no_improv_count = 0 + else + no_improv_count += 1 + end + no_improv_count <= rule.look_ahead || break + call_count <= rule.max_calls || break + end + _strong_branch_end(node) + return max_var +end + +function _strong_branch_score(; + node::Node, + pool::NodePool, + var::Variable, + x::Float64, + side_effect::Bool, +)::Tuple{Float64,Int} + obj_up, obj_down = 0, 0 + try + obj_up, obj_down = probe(node.mip, var) + catch + @warn "strong branch error" var=var + end + obj_change_up = obj_up - node.obj + obj_change_down = obj_down - node.obj + if side_effect + _update_var_history( + pool = pool, + var = var, + x = x, + obj_change_down = obj_change_down, + obj_change_up = obj_change_up, + ) + end + return (obj_change_up * obj_change_down, var.index) +end + +function _strong_branch_start(node::Node) + set_bounds!(node.mip, node.branch_variables, node.branch_values, node.branch_values) +end + +function _strong_branch_end(node::Node) + nbranch = length(node.branch_variables) + set_bounds!(node.mip, node.branch_variables, zeros(nbranch), ones(nbranch)) +end diff --git a/test/bb/lp_test.jl b/test/bb/lp_test.jl new file mode 100644 index 0000000..304e11b --- /dev/null +++ b/test/bb/lp_test.jl @@ -0,0 +1,115 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using Clp +using JuMP +using Test +using MIPLearn.BB + +basepath = @__DIR__ + +function runtests(optimizer_name, optimizer; large = true) + @testset "Solve ($optimizer_name)" begin + @testset "interface" begin + filename = "$basepath/../fixtures/danoint.mps.gz" + + mip = BB.init(optimizer) + BB.read!(mip, filename) + + @test mip.sense == 1.0 + @test length(mip.binary_variables) == 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" + + vals = BB.values(mip, mip.binary_variables) + @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]) + @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]) + 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]) + 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], + ) + status, obj = BB.solve_relaxation!(mip) + @test status == :Optimal + probe_up, probe_down = BB.probe(mip, mip.binary_variables[3]) + @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)) + 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)) + status, obj = BB.solve_relaxation!(mip) + @test status == :Optimal + @test round(obj, digits = 6) == 62.637280 + end + + @testset "varbranch" begin + branch_rules = [ + BB.RandomBranching(), + BB.FirstInfeasibleBranching(), + BB.LeastInfeasibleBranching(), + BB.MostInfeasibleBranching(), + BB.PseudocostBranching(), + BB.StrongBranching(), + BB.ReliabilityBranching(), + BB.HybridBranching(), + ] + for branch_rule in branch_rules + filename = "$basepath/../fixtures/vpm2.mps.gz" + mip = BB.init(optimizer) + BB.read!(mip, filename) + @info optimizer_name, branch_rule + @time BB.solve!( + mip, + initial_primal_bound = 13.75, + print_interval = 10, + node_limit = 100, + branch_rule = branch_rule, + ) + end + end + end +end + +@testset "BB" begin + @time runtests("Clp", Clp.Optimizer) + if is_gurobi_available + using Gurobi + @time runtests("Gurobi", Gurobi.Optimizer) + end +end diff --git a/test/fixtures/danoint.mps.gz b/test/fixtures/danoint.mps.gz new file mode 100644 index 0000000..1222df5 Binary files /dev/null and b/test/fixtures/danoint.mps.gz differ diff --git a/test/fixtures/vpm2.mps.gz b/test/fixtures/vpm2.mps.gz new file mode 100644 index 0000000..deb1c34 Binary files /dev/null and b/test/fixtures/vpm2.mps.gz differ diff --git a/test/runtests.jl b/test/runtests.jl index b41c71a..5122003 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using Test using MIPLearn MIPLearn.setup_logger() +const is_gurobi_available = ("GUROBI_HOME" in keys(ENV)) @testset "MIPLearn" begin include("fixtures/knapsack.jl") @@ -15,4 +16,5 @@ MIPLearn.setup_logger() include("solvers/learning_solver_test.jl") # include("utils/benchmark_test.jl") include("utils/parse_test.jl") + include("bb/lp_test.jl") end diff --git a/test/solvers/jump_solver_test.jl b/test/solvers/jump_solver_test.jl index d6bb22d..c353760 100644 --- a/test/solvers/jump_solver_test.jl +++ b/test/solvers/jump_solver_test.jl @@ -8,7 +8,6 @@ using MIPLearn using PyCall using Test -const is_gurobi_available = ("GUROBI_HOME" in keys(ENV)) if is_gurobi_available using Gurobi end