Skip to content

Commit

Permalink
Adding Node Finalizer (#568)
Browse files Browse the repository at this point in the history
  • Loading branch information
artalvpes authored Oct 2, 2021
1 parent 6edb71d commit 5870353
Show file tree
Hide file tree
Showing 5 changed files with 254 additions and 6 deletions.
4 changes: 3 additions & 1 deletion src/Algorithm/branching/branchingalgo.jl
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,7 @@ function perform_strong_branching_with_phases!(
getmaster(reform), getoptstate(input), exploitsprimalsolutions, false
)

children_generated = false
for (phase_index, current_phase) in enumerate(algo.phases)
nb_candidates_for_next_phase::Int64 = 1
if phase_index < length(algo.phases)
Expand Down Expand Up @@ -127,7 +128,7 @@ function perform_strong_branching_with_phases!(

for (group_index,group) in enumerate(groups)
#TO DO: verify if time limit is reached
if phase_index == 1
if !children_generated
generate_children!(group, env, reform, parent)
else
regenerate_children!(group, parent)
Expand Down Expand Up @@ -178,6 +179,7 @@ function perform_strong_branching_with_phases!(
end
print_bounds_and_score(group, phase_index, max_descr_length)
end
children_generated = true

sort!(groups, rev = true, by = x -> (x.isconquered, x.score))

Expand Down
5 changes: 3 additions & 2 deletions src/Algorithm/colgen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -584,8 +584,6 @@ function move_convexity_constrs_dual_values!(
for (spuid, spinfo) in spinfos
spinfo.lb_dual = dualsol[spinfo.lb_constr_id]
spinfo.ub_dual = dualsol[spinfo.ub_constr_id]
dualsol[spinfo.lb_constr_id] = zero(0.0)
dualsol[spinfo.ub_constr_id] = zero(0.0)
newbound -= (spinfo.lb_dual * spinfo.lb + spinfo.ub_dual * spinfo.ub)
end
constrids = Vector{ConstrId}()
Expand Down Expand Up @@ -681,6 +679,9 @@ function cg_main_loop!(
# this is needed due to convention that MOI uses for signs of duals in the maximization case
change_values_sign!(lp_dual_sol)
end
if lp_dual_sol !== nothing
set_lp_dual_sol!(cg_optstate, lp_dual_sol)
end
lp_dual_sol = move_convexity_constrs_dual_values!(spinfos, lp_dual_sol)

TO.@timeit Coluna._to "Getting primal solution" begin
Expand Down
71 changes: 68 additions & 3 deletions src/Algorithm/conquer.jl
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,17 @@ ParamRestrictedMasterHeuristic() =
1.0, 1.0, 1, 1000, "Restricted Master IP"
)

####################################################################
# NodeFinalizer
####################################################################

struct NodeFinalizer
algorithm::AbstractOptimizationAlgorithm
frequency::Integer
min_depth::Integer
name::String
end

####################################################################
# BendersConquer
####################################################################
Expand Down Expand Up @@ -149,6 +160,7 @@ Parameters :
@with_kw struct ColCutGenConquer <: AbstractConquerAlgorithm
stages::Vector{ColumnGeneration} = [ColumnGeneration()]
primal_heuristics::Vector{ParameterisedHeuristic} = [ParamRestrictedMasterHeuristic()]
node_finalizer::Union{Nothing, NodeFinalizer} = nothing
preprocess = PreprocessAlgorithm()
cutgen = CutCallbacks()
max_nb_cut_rounds::Int = 3 # TODO : tailing-off ?
Expand Down Expand Up @@ -184,7 +196,7 @@ function run!(algo::ColCutGenConquer, env::Env, reform::Reformulation, input::Co
restore_from_records!(input)
node = getnode(input)
nodestate = getoptstate(node)
if algo.run_preprocessing && isinfeasible(run!(algo.preprocess, reform, EmptyInput()))
if algo.run_preprocessing && isinfeasible(run!(algo.preprocess, env, reform, EmptyInput()))
setterminationstatus!(nodestate, INFEASIBLE)
return
end
Expand Down Expand Up @@ -245,7 +257,7 @@ function run!(algo::ColCutGenConquer, env::Env, reform::Reformulation, input::Co
end
sort!(heuristics_to_run, by = x -> last(x), rev=true)

for (heur_algorithm, name, priority) in heuristics_to_run
for (heur_algorithm, name, _) in heuristics_to_run
if ip_gap_closed(nodestate, atol = algo.opt_atol, rtol = algo.opt_rtol)
break
end
Expand All @@ -258,7 +270,6 @@ function run!(algo::ColCutGenConquer, env::Env, reform::Reformulation, input::Co
heur_output = run!(heur_algorithm, env, reform, OptimizationInput(nodestate))
status = getterminationstatus(getoptstate(heur_output))
status == TIME_LIMIT && setterminationstatus!(nodestate, status)

ip_primal_sols = get_ip_primal_sols(getoptstate(heur_output))
if ip_primal_sols !== nothing && length(ip_primal_sols) > 0
# we start with worst solution to add all improving solutions
Expand All @@ -280,6 +291,60 @@ function run!(algo::ColCutGenConquer, env::Env, reform::Reformulation, input::Co
break
end
end

# if the gap is still unclosed, try to run the node finalizer if any
run_node_finalizer = (algo.node_finalizer !== nothing)
run_node_finalizer = run_node_finalizer && getterminationstatus(nodestate) != TIME_LIMIT
run_node_finalizer =
run_node_finalizer && !ip_gap_closed(nodestate, atol = algo.opt_atol, rtol = algo.opt_rtol)
run_node_finalizer = run_node_finalizer && getdepth(node) >= algo.node_finalizer.min_depth
run_node_finalizer =
run_node_finalizer && mod(get_tree_order(node) - 1, algo.node_finalizer.frequency) == 0

if run_node_finalizer
# get the algorithm info
nodefinalizer = algo.node_finalizer.algorithm
name = algo.node_finalizer.name

@info "Running $name node finalizer"
if ismanager(nodefinalizer)
recordids = store_records!(reform)
end

nf_output = run!(nodefinalizer, env, reform, OptimizationInput(nodestate))
status = getterminationstatus(getoptstate(nf_output))
status == TIME_LIMIT && setterminationstatus!(nodestate, status)
ip_primal_sols = get_ip_primal_sols(getoptstate(nf_output))

# if the node has been conquered by the node finalizer
if status in (OPTIMAL, INFEASIBLE)
# set the ip solutions found without checking the cuts and finish
if ip_primal_sols !== nothing && length(ip_primal_sols) > 0
for sol in sort(ip_primal_sols)
update_ip_primal_sol!(nodestate, sol)
end
end

# make sure that the gap is closed for the current node
dual_bound = DualBound(reform, getvalue(get_ip_primal_bound(nodestate)))
update_ip_dual_bound!(nodestate, dual_bound)
else
if ip_primal_sols !== nothing && length(ip_primal_sols) > 0
# we start with worst solution to add all improving solutions
for sol in sort(ip_primal_sols)
cutgen = CutCallbacks(call_robust_facultative = false)
# TO DO : Node finalizer should ensure itselves that the returned solution is feasible
cutcb_output = run!(cutgen, env, getmaster(reform), CutCallbacksInput(sol))
if cutcb_output.nb_cuts_added == 0
update_ip_primal_sol!(nodestate, sol)
end
end
end
if ismanager(nodefinalizer)
ColunaBase.restore_from_records!(input.units_to_restore, recordids)
end
end
end
end

if ip_gap_closed(nodestate, atol = algo.opt_atol, rtol = algo.opt_rtol)
Expand Down
173 changes: 173 additions & 0 deletions test/node_finalizer_tests.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
# This file implements a toy bin packing model for Node Finalizer. It solves an instance with
# three items where any two of them fits into a bin but the three together do not. Pricing is
# solved by inspection onn the set of six possible solutions (three singletons and three pairs)
# which gives a fractional solution at the root node. Then node finalizer function
# "enumerative_finalizer" is called to find the optimal solution still at the root node and
# avoid branching (which would fail because maxnumnodes is set to 1).
# If "heuristic_finalizer" is true, then it allows branching and assumes that the solution found
# is not necessarily optimal.
CL.@with_kw struct EnumerativeFinalizer <: ClA.AbstractOptimizationAlgorithm
optimizer::Function
end

function ClA.run!(
algo::EnumerativeFinalizer, env::CL.Env, reform::ClMP.Reformulation, input::ClA.OptimizationInput
)::ClA.OptimizationOutput
masterform = ClMP.getmaster(reform)
_, spform = first(ClMP.get_dw_pricing_sps(reform))
cbdata = ClMP.PricingCallbackData(spform, 1)
isopt, primal_sol = algo.optimizer(masterform, cbdata)
result = ClA.OptimizationState(
masterform,
ip_primal_bound = ClA.get_ip_primal_bound(ClA.getoptstate(input)),
termination_status = isopt ? CL.OPTIMAL : CL.OTHER_LIMIT
)
if primal_sol !== nothing
ClA.add_ip_primal_sol!(result, primal_sol)
end
return ClA.OptimizationOutput(result)
end

function node_finalizer_tests(heuristic_finalizer)

function build_toy_model(optimizer)
toy = BlockModel(optimizer, direct_model = true)
I = [1, 2, 3]
@axis(B, [1])
@variable(toy, y[b in B] >= 0, Int)
@variable(toy, x[b in B, i in I], Bin)
@constraint(toy, sp[i in I], sum(x[b,i] for b in B) == 1)
@objective(toy, Min, sum(y[b] for b in B))
@dantzig_wolfe_decomposition(toy, dec, B)

return toy, x, y, dec, B
end

@testset "Optimization algorithms that may conquer a node" begin

call_enumerative_finalizer(masterform, cbdata) = enumerative_finalizer(masterform, cbdata)

coluna = JuMP.optimizer_with_attributes(
CL.Optimizer,
"default_optimizer" => GLPK.Optimizer,
"params" => CL.Params(
solver = ClA.TreeSearchAlgorithm(
conqueralg = ClA.ColCutGenConquer(
stages = [ClA.ColumnGeneration(
pricing_prob_solve_alg = ClA.SolveIpForm(
optimizer_id = 1
))
],
primal_heuristics = [],
node_finalizer = ClA.NodeFinalizer(
EnumerativeFinalizer(optimizer = call_enumerative_finalizer),
1, 0, "Enumerative"
)
),
maxnumnodes = heuristic_finalizer ? 10000 : 1
)
)
)

model, x, y, dec, B = build_toy_model(coluna)

function enumerative_pricing(cbdata)
# Get the reduced costs of the original variables
I = [1, 2, 3]
b = BlockDecomposition.callback_spid(cbdata, model)
rc_y = BD.callback_reduced_cost(cbdata, y[b])
rc_x = [BD.callback_reduced_cost(cbdata, x[b, i]) for i in I]

# check all possible solutions
sols = [[1], [2], [3], [1, 2], [1, 3], [2, 3]]
best_s = Int[]
best_rc = Inf
for s in sols
rc_s = rc_y + sum(rc_x[i] for i in s)
if rc_s < best_rc
best_rc = rc_s
best_s = s
end
end

# build the best one and submit
solcost = best_rc
solvars = JuMP.VariableRef[]
solvarvals = Float64[]
for i in best_s
push!(solvars, x[b, i])
push!(solvarvals, 1.0)
end
push!(solvars, y[b])
push!(solvarvals, 1.0)

# Submit the solution
MOI.submit(
model, BD.PricingSolution(cbdata), solcost, solvars, solvarvals
)
return
end
subproblems = BD.getsubproblems(dec)
BD.specify!.(
subproblems, lower_multiplicity = 0, upper_multiplicity = 3,
solver = enumerative_pricing
)

function enumerative_finalizer(masterform, cbdata)
# Get the reduced costs of the original variables
I = [1, 2, 3]
b = BlockDecomposition.callback_spid(cbdata, model)
rc_y = BD.callback_reduced_cost(cbdata, y[b])
rc_x = [BD.callback_reduced_cost(cbdata, x[b, i]) for i in I]
@test (rc_y, rc_x) == (1.0, [-0.5, -0.5, -0.5])

# Add the columns that are possibly missing for the solution [[1], [2,3]] in the master problem
# [1]
opt = JuMP.backend(model)
vars = [y[b], x[b, 1]]
varids = [CL._get_orig_varid_in_form(opt, cbdata.form, v) for v in JuMP.index.(vars)]
push!(varids, cbdata.form.duty_data.setup_var)
sol = ClMP.PrimalSolution(cbdata.form, varids, [1.0, 1.0, 1.0], 1.0, CL.FEASIBLE_SOL)
var_was_inserted, sol_id = ClMP.setprimalsol!(cbdata.form, sol)
if var_was_inserted
mc_1 = ClMP.setcol_from_sp_primalsol!(
masterform, cbdata.form, sol_id, string("MC_", ClA.getsortuid(sol_id)), ClMP.MasterCol
)
else
mc_1 = ClMP.getvar(masterform, sol_id)
end
# [2, 3]
vars = [y[b], x[b, 2], x[b, 3]]
varids = [CL._get_orig_varid_in_form(opt, cbdata.form, v) for v in JuMP.index.(vars)]
push!(varids, cbdata.form.duty_data.setup_var)
sol = ClMP.PrimalSolution(cbdata.form, varids, [1.0, 1.0, 1.0, 1.0], 1.0, CL.FEASIBLE_SOL)
var_was_inserted, sol_id = ClMP.setprimalsol!(cbdata.form, sol)
if var_was_inserted
mc_2_3 = ClMP.setcol_from_sp_primalsol!(
masterform, cbdata.form, sol_id, string("MC_", ClA.getsortuid(sol_id)), ClMP.MasterCol
)
else
mc_2_3 = ClMP.getvar(masterform, sol_id)
end

# add the solution to the master problem
varids = [ClMP.getid(mc_1), ClMP.getid(mc_2_3)]
primal_sol = ClMP.PrimalSolution(masterform, varids, [1.0, 1.0], 2.0, CL.FEASIBLE_SOL)
return !heuristic_finalizer, primal_sol
end

JuMP.optimize!(model)
@show JuMP.objective_value(model)
@test JuMP.termination_status(model) == MOI.OPTIMAL
for b in B
sets = BD.getsolutions(model, b)
for s in sets
@test BD.value(s) == 1.0 # value of the master column variable
@test BD.value(s, x[b, 1]) != BD.value(s, x[b, 2]) # only x[1,1] in its set
@test BD.value(s, x[b, 1]) != BD.value(s, x[b, 3]) # only x[1,1] in its set
@test BD.value(s, x[b, 2]) == BD.value(s, x[b, 3]) # x[1,2] and x[1,3] in the same set
end
end
end

end
7 changes: 7 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ include("optimizer_with_attributes_test.jl")
include("subproblem_solvers_tests.jl")
include("custom_var_cuts_tests.jl")
include("sol_disaggregation_tests.jl")
include("node_finalizer_tests.jl")

rng = MersenneTwister(1234123)

Expand Down Expand Up @@ -86,3 +87,9 @@ end
@testset "Solution Disaggregation" begin
sol_disaggregation_tests()
end

@testset "Node Finalizer" begin
node_finalizer_tests(false) # exact node finalizer

node_finalizer_tests(true) # heuristic node finalizer
end

0 comments on commit 5870353

Please sign in to comment.