ANL-CEEESA
/
MIPLearn.jl
mirror of https://github.com/ANL-CEEESA/MIPLearn.jl.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

FisSal2011: partial implementation

Browse Source
fs11_01
Alinson S. Xavier 2 months ago
parent 05e7d1619c
commit 0b5ec4740e
4 changed files with 328 additions and 158 deletions
Show Stats Download Patch File Download Diff File

109
src/Cuts/tableau/gmi.jl
Unescape View File

@ -199,23 +199,81 @@ function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001)
return [candidate_rows[perm[i]] for i = 1:min(length(perm), max_rows)]
end
# function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
# @timeit "Initialization" begin
# nrows, ncols = size(tableau.lhs)
# ub = Float64[Inf for _ = 1:nrows]
# lb = Float64[0.999 for _ = 1:nrows]
# tableau_I, tableau_J, tableau_V = findnz(tableau.lhs)
# lhs_I = Int[]
# lhs_J = Int[]
# lhs_V = Float64[]
# end
# @timeit "Compute coefficients" begin
# for k = 1:nnz(tableau.lhs)
# i::Int = tableau_I[k]
# j::Int = tableau_J[k]
# v::Float64 = 0.0
# frac_alpha_j = frac(tableau_V[k])
# alpha_j = tableau_V[k]
# beta = frac(tableau.rhs[i])
# if data.var_types[j] == 'C'
# if alpha_j >= 0
# v = alpha_j / beta
# else
# v = -alpha_j / (1 - beta)
# end
# else
# if frac_alpha_j < beta
# v = frac_alpha_j / beta
# else
# v = (1 - frac_alpha_j) / (1 - beta)
# end
# end
# if abs(v) > 1e-8
# push!(lhs_I, i)
# push!(lhs_J, tableau_J[k])
# push!(lhs_V, v)
# end
# end
# end
# @timeit "Convert to ConstraintSet" begin
# lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols)
# cs = ConstraintSet(; lhs, ub, lb)
# end
# return cs
# end
function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
nrows, ncols = size(tableau.lhs)
ub = Float64[Inf for _ = 1:nrows]
lb = Float64[0.999 for _ = 1:nrows]
tableau_I, tableau_J, tableau_V = findnz(tableau.lhs)
lhs_I = Int[]
lhs_J = Int[]
lhs_V = Float64[]
@timeit "Initialization" begin
nrows::Int, ncols::Int = size(tableau.lhs)
var_types::Vector{Char} = data.var_types
tableau_rhs::Vector{Float64} = tableau.rhs
tableau_I::Vector{Int}, tableau_J::Vector{Int}, tableau_V::Vector{Float64} = findnz(tableau.lhs)
end
@timeit "Pre-allocation" begin
ub::Vector{Float64} = fill(Inf, nrows)
lb::Vector{Float64} = fill(0.999, nrows)
nnz_tableau::Int = length(tableau_I)
lhs_I::Vector{Int} = Vector{Int}(undef, nnz_tableau)
lhs_J::Vector{Int} = Vector{Int}(undef, nnz_tableau)
lhs_V::Vector{Float64} = Vector{Float64}(undef, nnz_tableau)
nnz_count::Int = 0
end
@timeit "Compute coefficients" begin
for k = 1:nnz(tableau.lhs)
@inbounds for k = 1:nnz_tableau
i::Int = tableau_I[k]
j::Int = tableau_J[k]
v::Float64 = 0.0
frac_alpha_j = frac(tableau_V[k])
alpha_j = tableau_V[k]
beta = frac(tableau.rhs[i])
if data.var_types[j] == 'C'
alpha_j::Float64 = tableau_V[k]
frac_alpha_j::Float64 = alpha_j - floor(alpha_j)
beta_i::Float64 = tableau_rhs[i]
beta::Float64 = beta_i - floor(beta_i)
v::Float64 = 0
# Compute coefficient
if var_types[j] == 'C'
if alpha_j >= 0
v = alpha_j / beta
else
@ -228,16 +286,31 @@ function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet
v = (1 - frac_alpha_j) / (1 - beta)
end
end
# Store if significant
if abs(v) > 1e-8
push!(lhs_I, i)
push!(lhs_J, tableau_J[k])
push!(lhs_V, v)
nnz_count += 1
lhs_I[nnz_count] = i
lhs_J[nnz_count] = j
lhs_V[nnz_count] = v
end
end
lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols)
end
return ConstraintSet(; lhs, ub, lb)
@timeit "Resize arrays to actual size" begin
resize!(lhs_I, nnz_count)
resize!(lhs_J, nnz_count)
resize!(lhs_V, nnz_count)
end
@timeit "Convert to ConstraintSet" begin
lhs::SparseMatrixCSC = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols)
cs::ConstraintSet = ConstraintSet(; lhs, ub, lb)
end
return cs
end
export compute_gmi,
frac, select_gmi_rows, assert_cuts_off, assert_does_not_cut_off, collect_gmi

368
src/Cuts/tableau/gmi_dual.jl
Unescape View File

@ -269,9 +269,10 @@ function collect_gmi_FisSal2011(
optimizer,
max_rounds = 10_000,
max_cuts_per_round = 1_000_000,
time_limit = 30,
time_limit = 1_000_000,
interval_print=1,
max_cut_age=10,
silent_solver=true,
max_pool_size_mb = 1024,
)
reset_timer!()
initial_time = time()
@ -290,26 +291,37 @@ function collect_gmi_FisSal2011(
end
@timeit "Initialize" begin
stats_obj = []
stats_gap = []
stats_ncuts = []
pool = nothing
cut_age = nothing
lambda_curr = nothing
lambda_best = nothing
last_print_time = 0
obj_initial = nothing
obj_curr = 0
obj_best = 0
noimprov_count = 0
backtrack_count = 0
deterioration_count = 0
# ε = 0.01
basis_curr = nothing
basis_prev = nothing
basis_seen = Set{UInt64}()
gapcl_best = 0
gapcl_curr = 0
multipliers_curr = nothing
multipliers_best = nothing
last_print_time = 0
obj_curr = 0
obj_initial = nothing
obj_best = nothing
obj_hist = []
pool = nothing
pool_active = 0
pool_cut_age = nothing
pool_cut_hashes = Set{UInt64}()
pool_size_mb = 0
stats_gap = []
stats_ncuts = []
stats_obj = []
μ = 10
λ = 0
end
gap(v) = 100 * abs(obj_mip - v) / abs(obj_mip)
gapcl(v) = 100 * (v - obj_initial) / (obj_mip - obj_initial)
function perturb(v, ε=0.001)
function perturb(v, ε)
p = (1 - ε) .+ 2 * ε * rand(length(v))
return v .* p
end
@ -333,6 +345,9 @@ function collect_gmi_FisSal2011(
# Convert to standard form
data_s, transforms = convert_to_standard_form(data)
model_s = to_model(data_s)
if silent_solver
set_silent(model_s)
end
vars_s = all_variables(model_s)
orig_obj_s = objective_function(model_s)
set_optimizer(model_s, optimizer)
@ -345,13 +360,16 @@ function collect_gmi_FisSal2011(
assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb)
end
@info "Standard form model has $(length(data.var_lb)) vars, $(length(data.constr_lb)) constrs"
@show obj_mip
for round = 1:max_rounds
if round > 1
@timeit "Update objective function" begin
# Build Lagrangian term
lambda_perturbed = perturb(lambda_curr)
v = sparse(pool.lhs' * lambda_perturbed)
lagr_term = AffExpr(dot(lambda_perturbed, pool.lb))
# lambda_perturbed = perturb(multipliers, 0.1)
v = sparse(pool.lhs * multipliers_curr)
lagr_term = AffExpr(dot(multipliers_curr, pool.lb))
for offset in 1:nnz(v)
var_idx = v.nzind[offset]
add_to_expression!(
@ -370,32 +388,25 @@ function collect_gmi_FisSal2011(
end
@timeit "Optimize LP (lagrangian)" begin
set_silent(model_s)
optimize!(model_s)
sol_frac = get_x(model_s)
obj_curr = objective_value(model_s)
obj_curr <= obj_mip || error("LP value higher than MIP value: $(obj_curr) > $(obj_mip)")
if round == 1
obj_initial = obj_curr
push!(obj_hist, obj_curr)
if length(obj_hist) > 100
popfirst!(obj_hist)
end
if obj_curr >= obj_best
if obj_best === nothing || obj_curr > obj_best
obj_best = obj_curr
gapcl_best = gapcl(obj_best)
lambda_best = lambda_curr
noimprov_count = 0
else
noimprov_count += 1
multipliers_best = multipliers_curr
end
if noimprov_count > 10
lambda_curr = lambda_best
backtrack_count += 1
noimprov_count = 0
continue
if round == 1
obj_initial = obj_curr
end
gapcl_curr = gapcl(obj_curr)
gapcl_best = gapcl(obj_best)
push!(stats_obj, obj_curr)
push!(stats_gap, gap(obj_curr))
if round == 1
@ -406,121 +417,191 @@ function collect_gmi_FisSal2011(
if termination_status(model_s) != MOI.OPTIMAL
error("Non-optimal termination status")
end
end
@timeit "Select tableau rows" begin
basis = get_basis(model_s)
if round == 1
original_basis = basis
Δ = obj_mip - obj_best
if obj_curr < obj_best - Δ
deterioration_count += 1
else
deterioration_count = 0
end
if deterioration_count >= 10
μ *= 0.5
multipliers_curr = multipliers_best
deterioration_count = 0
backtrack_count += 1
elseif length(obj_hist) >= 100
obj_hist_avg = mean(obj_hist)
improv = obj_best - obj_hist[1]
if improv < 0.01 * Δ
if obj_best - obj_hist_avg < 0.001 * Δ
μ = 10 * μ
elseif obj_best - obj_hist_avg < 0.01 * Δ
μ = 2 * μ
else
μ = 0.5 * μ
end
end
end
selected_rows =
select_gmi_rows(data_s, basis, sol_frac, max_rows = max_cuts_per_round)
end
@timeit "Compute tableau rows" begin
tableau = compute_tableau(data_s, basis, x = sol_frac, rows = selected_rows)
if mod(round, 10) == 1
@timeit "Select tableau rows" begin
basis_prev = basis_curr
basis_curr = get_basis(model_s)
basis_hash = hash(basis_curr)
push!(basis_seen, basis_hash)
selected_rows =
select_gmi_rows(data_s, basis_curr, sol_frac, max_rows = max_cuts_per_round)
end
# Assert tableau rows have been computed correctly
assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3)
assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3)
end
@timeit "Compute tableau rows" begin
tableau = compute_tableau(data_s, basis_curr, x = sol_frac, rows = selected_rows)
@timeit "Compute GMI cuts" begin
cuts_s = compute_gmi(data_s, tableau)
assert_cuts_off(cuts_s, sol_frac)
assert_does_not_cut_off(cuts_s, sol_opt_s)
ncuts = length(cuts_s.lb)
end
# Assert tableau rows have been computed correctly
assert_eq(tableau.lhs * sol_frac, tableau.rhs, atol=1e-3)
assert_eq(tableau.lhs * sol_opt_s, tableau.rhs, atol=1e-3)
end
@timeit "Add new cuts to the pool" begin
if round == 1
pool = cuts_s
lambda_curr = zeros(ncuts)
lambda_best = zeros(ncuts)
cut_age = zeros(ncuts)
else
pool.lhs = [pool.lhs; cuts_s.lhs]
pool.lb = [pool.lb; cuts_s.lb]
pool.ub = [pool.ub; cuts_s.ub]
lambda_curr = [lambda_curr; zeros(ncuts)]
lambda_best = [lambda_best; zeros(ncuts)]
cut_age = [cut_age; zeros(ncuts)]
@timeit "Compute GMI cuts" begin
cuts_s = compute_gmi(data_s, tableau)
assert_cuts_off(cuts_s, sol_frac)
assert_does_not_cut_off(cuts_s, sol_opt_s)
ncuts = length(cuts_s.lb)
end
@timeit "Add new cuts to the pool" begin
@timeit "Compute cut hashses" begin
unique_indices = Int[]
for i in 1:ncuts
cut_data = (cuts_s.lhs[i, :], cuts_s.lb[i], cuts_s.ub[i])
cut_hash = hash(cut_data)
if !(cut_hash in pool_cut_hashes)
push!(pool_cut_hashes, cut_hash)
push!(unique_indices, i)
end
end
end
@timeit "Append unique cuts" begin
if round == 1
pool = ConstraintSet(
lhs = sparse(cuts_s.lhs[unique_indices, :]'),
lb = cuts_s.lb[unique_indices],
ub = cuts_s.ub[unique_indices]
)
ncuts_unique = length(unique_indices)
multipliers_curr = zeros(ncuts_unique)
multipliers_best = zeros(ncuts_unique)
pool_cut_age = zeros(ncuts_unique)
else
if !isempty(unique_indices)
ncuts_unique = length(unique_indices)
pool.lhs = [pool.lhs sparse(cuts_s.lhs[unique_indices, :]')]
pool.lb = [pool.lb; cuts_s.lb[unique_indices]]
pool.ub = [pool.ub; cuts_s.ub[unique_indices]]
multipliers_curr = [multipliers_curr; zeros(ncuts_unique)]
multipliers_best = [multipliers_best; zeros(ncuts_unique)]
pool_cut_age = [pool_cut_age; zeros(ncuts_unique)]
end
end
end
end
end
# @timeit "Update multipliers (subgradient)" begin
# subgrad = pool.lb .- pool.lhs * sol_frac
# lambda = max.(0, lambda .+ 0.01 * subgrad)
# end
selected_idx = []
selected_contrs = []
if round == 1 || round == max_rounds
@timeit "Update multipliers (large LP)" begin
while true
@timeit "Optimize LP (extended)" begin
set_objective_function(model_s, orig_obj_s)
optimize!(model_s)
sol_frac = get_x(model_s)
end
@timeit "Update multipliers (large LP)" begin
while true
@timeit "Optimize LP (extended)" begin
set_objective_function(model_s, orig_obj_s)
optimize!(model_s)
sol_frac = get_x(model_s)
end
@timeit "Computing cut violations" begin
violations = (pool.lb' - (sol_frac' * pool.lhs))'
end
@timeit "Find most violated cut" begin
violations = pool.lb .- pool.lhs * sol_frac
σ = sortperm(violations, rev=true)
end
@timeit "Sorting cut violations" begin
σ = sortperm(violations, rev=true)
end
# Stop if all cuts are satisfied
if violations[σ[1]] <= 1e-6
break
end
if violations[σ[1]] <= 1e-6
break
end
@timeit "Add constraint to the model" begin
push!(selected_idx, σ[1])
cut_lhs = pool.lhs[σ[1], :]
cut_lhs_value = 0.0
cut_lb = pool.lb[σ[1]]
cut_expr = AffExpr()
for offset in 1:nnz(cut_lhs)
var_idx = cut_lhs.nzind[offset]
add_to_expression!(
cut_expr,
vars_s[var_idx],
cut_lhs.nzval[offset],
)
cut_lhs_value += sol_frac[var_idx] * cut_lhs.nzval[offset]
@timeit "Add constraint to the model" begin
push!(selected_idx, σ[1])
cut_lhs = pool.lhs[:, σ[1]]
cut_lhs_value = 0.0
cut_lb = pool.lb[σ[1]]
cut_expr = AffExpr()
for offset in 1:nnz(cut_lhs)
var_idx = cut_lhs.nzind[offset]
add_to_expression!(
cut_expr,
vars_s[var_idx],
cut_lhs.nzval[offset],
)
cut_lhs_value += sol_frac[var_idx] * cut_lhs.nzval[offset]
end
cut_constr = @constraint(model_s, cut_expr >= cut_lb)
push!(selected_contrs, cut_constr)
end
cut_constr = @constraint(model_s, cut_expr >= cut_lb)
push!(selected_contrs, cut_constr)
end
end
@timeit "Find dual values for all selected cuts" begin
lambda_curr .= 0
cut_age .+= 1
for (offset, idx) in enumerate(selected_idx)
lambda_curr[idx] = -shadow_price(selected_contrs[offset])
if lambda_curr[idx] > 1e-5
cut_age[idx] = 0
pool_active = length(selected_idx)
@timeit "Find dual values for all selected cuts" begin
multipliers_curr .= 0
pool_cut_age .+= 1
for (offset, idx) in enumerate(selected_idx)
multipliers_curr[idx] = -shadow_price(selected_contrs[offset])
if multipliers_curr[idx] > 1e-5
pool_cut_age[idx] = 0
end
end
end
end
# Filter cut pool
keep = findall(cut_age .< max_cut_age)
pool.lhs = pool.lhs[keep, :]
pool.ub = pool.ub[keep]
pool.lb = pool.lb[keep]
lambda_curr = lambda_curr[keep]
lambda_best = lambda_best[keep]
cut_age = cut_age[keep]
@timeit "Prune cut pool" begin
pool_size_mb = Base.summarysize(pool) / 1024^2
while pool_size_mb >= max_pool_size_mb
@timeit "Identify cuts to remove" begin
σ = sortperm(pool_cut_age, rev=true)
pool_size = length(pool.ub)
n_keep = Int(floor(pool_size * 0.8))
idx_keep = 1:n_keep
idx_remove = σ[(n_keep+1):end]
end
@timeit "Update cut hashes" begin
for idx in idx_remove
cut_data = (pool.lhs[:, idx], pool.lb[idx], pool.ub[idx])
delete!(pool_cut_hashes, hash(cut_data))
end
end
@timeit "Update cut pool" begin
pool.lhs = pool.lhs[:, idx_keep]
pool.ub = pool.ub[idx_keep]
pool.lb = pool.lb[idx_keep]
multipliers_curr = multipliers[idx_keep]
pool_cut_age = pool_cut_age[idx_keep]
end
pool_size_mb = Base.summarysize(pool) / 1024^2
end
end
@timeit "Delete all cut constraints" begin
delete.(model_s, selected_contrs)
@timeit "Delete all cut constraints" begin
delete.(model_s, selected_contrs)
end
end
else
@timeit "Update multipliers (subgradient)" begin
subgrad = (pool.lb' - (sol_frac' * pool.lhs))'
λ = μ * (obj_mip - obj_curr) / norm(subgrad)^2
multipliers_curr = max.(0, multipliers_curr .+ λ * subgrad)
end
end
push!(stats_ncuts, length(pool.lb))
elapsed_time = time() - initial_time
if elapsed_time > time_limit
@info "Time limit exceeded. Stopping."
@ -529,28 +610,38 @@ function collect_gmi_FisSal2011(
if round == 1
@printf(
"%8s %12s %12s %12s %8s %8s %9s\n",
"%8s %9s %7s %7s %8s %9s %9s %9s %4s %8s %8s %8s\n",
"round",
"obj",
"gapcl_curr",
"gapcl_best",
"cl_curr",
"cl_best",
"active",
"pool",
"backtrack",
"pool_cuts",
"pool_mb",
"bases",
"bktk",
"Δ",
"μ",
"λ",
)
end
if time() - last_print_time > interval_print
last_print_time = time()
@printf(
"%8d %12.6e %12.2f %12.2f %8d %8d %9d\n",
"%8d %9.3e %7.2f %7.2f %8d %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n",
round,
obj_curr,
gapcl(obj_curr),
gapcl_curr,
gapcl_best,
length(selected_idx),
length(pool.ub),
pool_size_mb,
length(basis_seen),
backtrack_count,
Δ,
μ,
λ,
)
end
end
@ -585,18 +676,23 @@ function collect_gmi_FisSal2011(
# end
# end
@info "Best gap closure: $(gapcl_best)"
to = TimerOutputs.get_defaulttimer()
stats_time = TimerOutputs.tottime(to) / 1e9
print_timer()
return OrderedDict(
"gapcl_best" => gapcl_best,
"gapcl_curr" => gapcl_curr,
"instance" => mps_filename,
"max_rounds" => max_rounds,
"rounds" => length(stats_obj) - 1,
"obj_final" => obj_curr,
"obj_initial" => obj_initial,
"obj_mip" => obj_mip,
"stats_obj" => stats_obj,
"stats_ncuts" => stats_ncuts,
"stats_time" => stats_time,
"pool_active" => pool_active,
"pool_size_mb" => pool_size_mb,
"pool_total" => length(pool.lb),
"time" => stats_time,
)
end

8
src/Cuts/tableau/structs.jl
Unescape View File

@ -18,10 +18,10 @@ Base.@kwdef mutable struct ProblemData
end
Base.@kwdef mutable struct Tableau
obj::Any
lhs::Any
rhs::Any
z::Any
obj::Vector{Float64}
lhs::SparseMatrixCSC
rhs::Vector{Float64}
z::Float64
end
Base.@kwdef mutable struct Basis

1
src/Cuts/tableau/tableau.jl
Unescape View File

@ -102,6 +102,7 @@ function compute_tableau(
sol = factor \ obj_b
tableau_obj = -data.obj' + sol' * data.constr_lhs
tableau_obj[abs.(tableau_obj).<tol] .= 0
tableau_obj = Array(tableau_obj')
end
# Compute z if solution is provided

Write Preview
Loading…
Cancel
Save