From 05e7d1619ceca0dc4b83d99820d0a013304dfba1 Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Fri, 1 Aug 2025 09:43:09 -0500 Subject: [PATCH] Make dual GMI cuts stronger --- src/Cuts/tableau/gmi.jl | 4 +- src/Cuts/tableau/gmi_dual.jl | 341 +++++++++++++++++++++++++++++++++- src/Cuts/tableau/transform.jl | 11 +- 3 files changed, 350 insertions(+), 6 deletions(-) diff --git a/src/Cuts/tableau/gmi.jl b/src/Cuts/tableau/gmi.jl index 506ad2e..bdde291 100644 --- a/src/Cuts/tableau/gmi.jl +++ b/src/Cuts/tableau/gmi.jl @@ -185,7 +185,7 @@ function collect_gmi( ) end -function select_gmi_rows(data, basis, x; max_rows = 10, atol = 1e-4) +function select_gmi_rows(data, basis, x; max_rows = 10, atol = 0.001) candidate_rows = [ r for r = 1:length(basis.var_basic) if ( (data.var_types[basis.var_basic[r]] != 'C') && @@ -202,7 +202,7 @@ end function compute_gmi(data::ProblemData, tableau::Tableau)::ConstraintSet nrows, ncols = size(tableau.lhs) ub = Float64[Inf for _ = 1:nrows] - lb = Float64[0.9999 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[] diff --git a/src/Cuts/tableau/gmi_dual.jl b/src/Cuts/tableau/gmi_dual.jl index aae4bb7..d6535eb 100644 --- a/src/Cuts/tableau/gmi_dual.jl +++ b/src/Cuts/tableau/gmi_dual.jl @@ -7,6 +7,7 @@ using JuMP using HiGHS using Random using DataStructures +using Statistics import ..H5FieldsExtractor @@ -25,7 +26,7 @@ function collect_gmi_dual( mps_filename; optimizer, max_rounds = 10, - max_cuts_per_round = 500, + max_cuts_per_round = 1_000_000, time_limit = 3_600, ) reset_timer!() @@ -263,6 +264,342 @@ function collect_gmi_dual( ) end +function collect_gmi_FisSal2011( + mps_filename; + optimizer, + max_rounds = 10_000, + max_cuts_per_round = 1_000_000, + time_limit = 30, + interval_print=1, + max_cut_age=10, +) + reset_timer!() + initial_time = time() + + @timeit "Read H5" begin + h5_filename = replace(mps_filename, ".mps.gz" => ".h5") + h5 = H5File(h5_filename, "r") + sol_opt_dict = Dict( + zip( + h5.get_array("static_var_names"), + convert(Array{Float64}, h5.get_array("mip_var_values")), + ), + ) + obj_mip = h5.get_scalar("mip_obj_value") + h5.file.close() + 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 + gapcl_best = 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) + p = (1 - ε) .+ 2 * ε * rand(length(v)) + return v .* p + end + + @timeit "Read problem" begin + model = read_from_file(mps_filename) + set_optimizer(model, optimizer) + end + + @timeit "Convert model to standard form" begin + # Extract problem data + data = ProblemData(model) + + # Construct optimal solution vector (with correct variable sequence) + sol_opt = [sol_opt_dict[n] for n in data.var_names] + + # Assert optimal solution is feasible for the original problem + 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) + model_s = to_model(data_s) + vars_s = all_variables(model_s) + orig_obj_s = objective_function(model_s) + set_optimizer(model_s, optimizer) + relax_integrality(model_s) + + # Convert optimal solution to standard form + sol_opt_s = forward(transforms, sol_opt) + + # Assert converted solution is feasible for standard form problem + assert_eq(data_s.constr_lhs * sol_opt_s, data_s.constr_lb) + end + + 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)) + for offset in 1:nnz(v) + var_idx = v.nzind[offset] + add_to_expression!( + lagr_term, + vars_s[var_idx], + - v.nzval[offset], + ) + end + + # Update objective + set_objective_function( + model_s, + orig_obj_s + lagr_term, + ) + end + 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 + end + + if obj_curr >= obj_best + obj_best = obj_curr + gapcl_best = gapcl(obj_best) + lambda_best = lambda_curr + noimprov_count = 0 + else + noimprov_count += 1 + end + + if noimprov_count > 10 + lambda_curr = lambda_best + backtrack_count += 1 + noimprov_count = 0 + continue + end + + push!(stats_obj, obj_curr) + push!(stats_gap, gap(obj_curr)) + if round == 1 + push!(stats_ncuts, 0) + else + push!(stats_ncuts, length(pool.lb)) + end + 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 + 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) + + # 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 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 + 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)] + 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 = [] + + @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 "Find most violated cut" begin + violations = pool.lb .- pool.lhs * sol_frac + σ = sortperm(violations, rev=true) + end + + # Stop if all cuts are satisfied + 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] + 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 + 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 "Delete all cut constraints" begin + delete.(model_s, selected_contrs) + end + end + + push!(stats_ncuts, length(pool.lb)) + + elapsed_time = time() - initial_time + if elapsed_time > time_limit + @info "Time limit exceeded. Stopping." + break + end + + if round == 1 + @printf( + "%8s %12s %12s %12s %8s %8s %9s\n", + "round", + "obj", + "gapcl_curr", + "gapcl_best", + "active", + "pool", + "backtrack", + ) + end + + if time() - last_print_time > interval_print + last_print_time = time() + @printf( + "%8d %12.6e %12.2f %12.2f %8d %8d %9d\n", + round, + obj_curr, + gapcl(obj_curr), + gapcl_best, + length(selected_idx), + length(pool.ub), + backtrack_count, + ) + end + end + + # @timeit "Store cuts in H5 file" begin + # if all_cuts !== nothing + # ncuts = length(all_cuts_rows) + # total = + # length(original_basis.var_basic) + + # length(original_basis.var_nonbasic) + + # length(original_basis.constr_basic) + + # length(original_basis.constr_nonbasic) + # all_cuts_basis_sizes = Array{Int64,2}(undef, ncuts, 4) + # all_cuts_basis_vars = Array{Int64,2}(undef, ncuts, total) + # for i = 1:ncuts + # vb = all_cuts_bases[i].var_basic + # vn = all_cuts_bases[i].var_nonbasic + # cb = all_cuts_bases[i].constr_basic + # cn = all_cuts_bases[i].constr_nonbasic + # all_cuts_basis_sizes[i, :] = [length(vb) length(vn) length(cb) length(cn)] + # all_cuts_basis_vars[i, :] = [vb' vn' cb' cn'] + # end + # @info "Storing $(length(all_cuts.ub)) GMI cuts..." + # h5 = H5File(h5_filename) + # h5.put_sparse("cuts_lhs", all_cuts.lhs) + # h5.put_array("cuts_lb", all_cuts.lb) + # h5.put_array("cuts_ub", all_cuts.ub) + # h5.put_array("cuts_basis_vars", all_cuts_basis_vars) + # h5.put_array("cuts_basis_sizes", all_cuts_basis_sizes) + # h5.put_array("cuts_rows", all_cuts_rows) + # h5.file.close() + # end + # end + + to = TimerOutputs.get_defaulttimer() + stats_time = TimerOutputs.tottime(to) / 1e9 + print_timer() + + return OrderedDict( + "instance" => mps_filename, + "max_rounds" => max_rounds, + "rounds" => length(stats_obj) - 1, + "obj_mip" => obj_mip, + "stats_obj" => stats_obj, + "stats_ncuts" => stats_ncuts, + "stats_time" => stats_time, + ) +end + function add_constraint_set_dual_v2(model::JuMP.Model, cs::ConstraintSet) vars = all_variables(model) nrows, ncols = size(cs.lhs) @@ -576,5 +913,5 @@ function __init_gmi_dual__() copy!(ExpertDualGmiComponent, ExpertDualGmiComponentPy) end -export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent +export collect_gmi_dual, expert_gmi_dual, ExpertDualGmiComponent, KnnDualGmiComponent, collect_gmi_FisSal2011 diff --git a/src/Cuts/tableau/transform.jl b/src/Cuts/tableau/transform.jl index 2cfcf5f..bce5a1e 100644 --- a/src/Cuts/tableau/transform.jl +++ b/src/Cuts/tableau/transform.jl @@ -102,7 +102,14 @@ function forward!(t::AddSlackVariables, data::ProblemData) ge = [i for i = 1:nrows if isfinite(data.constr_lb[i]) && !isequality[i]] le = [i for i = 1:nrows if isfinite(data.constr_ub[i]) && !isequality[i]] EQ, GE, LE = length(eq), length(ge), length(le) - + is_integral(row_idx, rhs) = ( + abs(rhs - round(rhs)) <= 1e-6 && + all(j -> data.var_types[j] ∈ ['I', 'B'], findnz(data.constr_lhs[row_idx, :])[1]) + ) + slack_types = [ + [is_integral(ge[i], data.constr_lb[ge[i]]) ? 'I' : 'C' for i = 1:GE]; + [is_integral(le[i], data.constr_ub[le[i]]) ? 'I' : 'C' for i = 1:LE] + ] t.M1 = [ I spzeros(ncols, GE + LE) data.constr_lhs[ge, :] spzeros(GE, GE + LE) @@ -129,7 +136,7 @@ function forward!(t::AddSlackVariables, data::ProblemData) data.var_lb = [data.var_lb; zeros(GE + LE)] data.var_ub = [data.var_ub; [Inf for _ = 1:(GE+LE)]] data.var_names = [data.var_names; ["__s$i" for i = 1:(GE+LE)]] - data.var_types = [data.var_types; ['C' for _ = 1:(GE+LE)]] + data.var_types = [data.var_types; slack_types] data.constr_lb = [ data.constr_lb[eq] data.constr_lb[ge]