From 9c61b98cb9a3223cdcff811fd0bfecd2409c31b4 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Tue, 12 Mar 2024 13:56:34 -0500 Subject: [PATCH] Make GMI cuts more stable --- src/Cuts/Cuts.jl | 1 + src/Cuts/tableau/collect.jl | 40 ++++++++++++------------- src/Cuts/tableau/gmi.jl | 57 ++++++++++++------------------------ src/Cuts/tableau/numerics.jl | 51 ++++++++++++++++++++++++++++++++ 4 files changed, 89 insertions(+), 60 deletions(-) create mode 100644 src/Cuts/tableau/numerics.jl diff --git a/src/Cuts/Cuts.jl b/src/Cuts/Cuts.jl index 77738d6..b2a4254 100644 --- a/src/Cuts/Cuts.jl +++ b/src/Cuts/Cuts.jl @@ -9,6 +9,7 @@ import ..to_str_array include("tableau/structs.jl") # include("blackbox/cplex.jl") +include("tableau/numerics.jl") include("tableau/collect.jl") include("tableau/gmi.jl") include("tableau/moi.jl") diff --git a/src/Cuts/tableau/collect.jl b/src/Cuts/tableau/collect.jl index 01e84ff..0246449 100644 --- a/src/Cuts/tableau/collect.jl +++ b/src/Cuts/tableau/collect.jl @@ -5,8 +5,10 @@ import ..H5File using OrderedCollections +using Statistics -function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_round = 100) + +function collect_gmi(mps_filename; optimizer, max_rounds=10, max_cuts_per_round=100, atol=1e-4) @info mps_filename reset_timer!() @@ -27,7 +29,7 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun if obj_mip === nothing obj_mip = h5.get_scalar("mip_obj_value") end - obj_lp = nothing + obj_lp = h5.get_scalar("lp_obj_value") h5.file.close() # Define relative MIP gap @@ -58,8 +60,8 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun sol_opt = [sol_opt_dict[n] for n in data.var_names] # Assert optimal solution is feasible for the original problem - @assert all(data.constr_lb .- 1e-3 .<= data.constr_lhs * sol_opt) - @assert all(data.constr_lhs * sol_opt .<= data.constr_ub .+ 1e-3) + assert_leq(data.constr_lb, data.constr_lhs * sol_opt) + assert_leq(data.constr_lhs * sol_opt, data.constr_ub) # Convert to standard form data_s, transforms = convert_to_standard_form(data) @@ -71,15 +73,17 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun sol_opt_s = forward(transforms, sol_opt) # Assert converted solution is feasible for standard form problem - @assert data_s.constr_lhs * sol_opt_s ≈ data_s.constr_lb + assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb) end # Optimize standard form optimize!(model_s) stats_time_solve += solve_time(model_s) obj = objective_value(model_s) + data_s.obj_offset - if obj_lp === nothing - obj_lp = obj + + if round == 1 + # Assert standard form problem has same value as original + assert_eq(obj, obj_lp) push!(stats_obj, obj) push!(stats_gap, gap(obj)) push!(stats_ncuts, 0) @@ -93,16 +97,16 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun sol_frac = get_x(model_s) stats_time_select += @elapsed begin selected_rows = - select_gmi_rows(data_s, basis, sol_frac, max_rows = max_cuts_per_round) + select_gmi_rows(data_s, basis, sol_frac, max_rows=max_cuts_per_round) end # Compute selected tableau rows stats_time_tableau += @elapsed begin - tableau = compute_tableau(data_s, basis, sol_frac, rows = selected_rows) + tableau = compute_tableau(data_s, basis, sol_frac, rows=selected_rows) # Assert tableau rows have been computed correctly - @assert tableau.lhs * sol_frac ≈ tableau.rhs - @assert tableau.lhs * sol_opt_s ≈ tableau.rhs + assert_eq(tableau.lhs * sol_frac, tableau.rhs) + assert_eq(tableau.lhs * sol_opt_s, tableau.rhs) end # Compute GMI cuts @@ -110,17 +114,12 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun cuts_s = compute_gmi(data_s, tableau) # Assert cuts have been generated correctly - try - assert_cuts_off(cuts_s, sol_frac) - assert_does_not_cut_off(cuts_s, sol_opt_s) - catch - @warn "Invalid cuts detected. Discarding round $round cuts and aborting." - break - end + assert_cuts_off(cuts_s, sol_frac) + assert_does_not_cut_off(cuts_s, sol_opt_s) # Abort if no cuts are left if length(cuts_s.lb) == 0 - @info "No cuts generated. Aborting." + @info "No cuts generated. Stopping." break end end @@ -139,7 +138,7 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun push!(stats_gap, gap(obj)) # Store useful cuts; drop useless ones from the problem - useful = [abs(shadow_price(c)) > 1e-3 for c in constrs] + useful = [abs(shadow_price(c)) > atol for c in constrs] drop = findall(useful .== false) keep = findall(useful .== true) delete.(model, constrs[drop]) @@ -174,7 +173,6 @@ function collect_gmi(mps_filename; optimizer, max_rounds = 10, max_cuts_per_roun "time_tableau" => stats_time_tableau, "time_gmi" => stats_time_gmi, "obj_mip" => obj_mip, - "obj_lp" => obj_lp, "stats_obj" => stats_obj, "stats_gap" => stats_gap, "stats_ncuts" => stats_ncuts, diff --git a/src/Cuts/tableau/gmi.jl b/src/Cuts/tableau/gmi.jl index 29e65e1..f4d578f 100644 --- a/src/Cuts/tableau/gmi.jl +++ b/src/Cuts/tableau/gmi.jl @@ -5,13 +5,14 @@ using SparseArrays using TimerOutputs -@inline frac(x::Float64) = x - floor(x) - -function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001) +function select_gmi_rows(data, basis, x; max_rows=10, atol=1e-4) candidate_rows = [ r for - r = 1:length(basis.var_basic) if (data.var_types[basis.var_basic[r]] != 'C') && - (frac(x[basis.var_basic[r]]) > atol) + r in 1:length(basis.var_basic) if ( + (data.var_types[basis.var_basic[r]] != 'C') && + (frac(x[basis.var_basic[r]]) > atol) && + (frac2(x[basis.var_basic[r]]) > atol) + ) ] candidate_vals = frac.(x[basis.var_basic[candidate_rows]]) score = abs.(candidate_vals .- 0.5) @@ -19,34 +20,36 @@ 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, tol = 1e-8)::ConstraintSet +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] + lb = Float64[0.9999 for _ = 1:nrows] tableau_I, tableau_J, tableau_V = findnz(tableau.lhs) lhs_I = Int[] lhs_J = Int[] lhs_V = Float64[] @timeit "Compute coefficients" begin - for k = 1:nnz(tableau.lhs) + for k in 1:nnz(tableau.lhs) i::Int = tableau_I[k] + j::Int = tableau_J[k] v::Float64 = 0.0 - alpha_j = frac(tableau_V[k]) + frac_alpha_j = frac(tableau_V[k]) + alpha_j = tableau_V[k] beta = frac(tableau.rhs[i]) - if data.var_types[i] == "C" + if data.var_types[j] == 'C' if alpha_j >= 0 v = alpha_j / beta else - v = alpha_j / (1 - beta) + v = -alpha_j / (1 - beta) end else - if alpha_j <= beta - v = alpha_j / beta + if frac_alpha_j < beta + v = frac_alpha_j / beta else - v = (1 - alpha_j) / (1 - beta) + v = (1 - frac_alpha_j) / (1 - beta) end end - if abs(v) > tol + if abs(v) > 1e-8 push!(lhs_I, i) push!(lhs_J, tableau_J[k]) push!(lhs_V, v) @@ -57,28 +60,4 @@ function compute_gmi(data::ProblemData, tableau::Tableau, tol = 1e-8)::Constrain return ConstraintSet(; lhs, ub, lb) end -function assert_cuts_off(cuts::ConstraintSet, x::Vector{Float64}, tol = 1e-6) - for i = 1:length(cuts.lb) - val = cuts.lhs[i, :]' * x - if (val <= cuts.ub[i] - tol) && (val >= cuts.lb[i] + tol) - throw(ErrorException("inequality fails to cut off fractional solution")) - end - end -end - -function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol = 1e-6) - for i = 1:length(cuts.lb) - val = cuts.lhs[i, :]' * x - ub = cuts.ub[i] - lb = cuts.lb[i] - if (val >= ub) || (val <= lb) - throw( - ErrorException( - "inequality $i cuts off integer solution ($lb <= $val <= $ub)", - ), - ) - end - end -end - export compute_gmi, frac, select_gmi_rows, assert_cuts_off, assert_does_not_cut_off diff --git a/src/Cuts/tableau/numerics.jl b/src/Cuts/tableau/numerics.jl new file mode 100644 index 0000000..2a77547 --- /dev/null +++ b/src/Cuts/tableau/numerics.jl @@ -0,0 +1,51 @@ +@inline frac(x::Float64) = x - floor(x) + +@inline frac2(x::Float64) = ceil(x) - x + +function assert_leq(a, b; atol=0.01) + if !all(a .<= b .+ atol) + delta = a .- b + for i in eachindex(delta) + if delta[i] > atol + @info "Assertion failed: a[$i] = $(a[i]) <= $(b[i]) = b[$i]" + end + end + error("assert_leq failed") + end +end + +function assert_eq(a, b; atol=1e-4) + if !all(abs.(a .- b) .<= atol) + delta = abs.(a .- b) + for i in eachindex(delta) + if delta[i] > atol + @info "Assertion failed: a[$i] = $(a[i]) == $(b[i]) = b[$i]" + end + end + error("assert_eq failed") + end +end + +function assert_cuts_off(cuts::ConstraintSet, x::Vector{Float64}, tol=1e-6) + for i = 1:length(cuts.lb) + val = cuts.lhs[i, :]' * x + if (val <= cuts.ub[i] - tol) && (val >= cuts.lb[i] + tol) + throw(ErrorException("inequality fails to cut off fractional solution")) + end + end +end + +function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol=1e-6) + for i = 1:length(cuts.lb) + val = cuts.lhs[i, :]' * x + ub = cuts.ub[i] + lb = cuts.lb[i] + if (val >= ub) || (val <= lb) + throw( + ErrorException( + "inequality $i cuts off integer solution ($lb <= $val <= $ub)", + ), + ) + end + end +end