diff --git a/src/Cuts/tableau/gmi_dual.jl b/src/Cuts/tableau/gmi_dual.jl index 71c964f..edd7e60 100644 --- a/src/Cuts/tableau/gmi_dual.jl +++ b/src/Cuts/tableau/gmi_dual.jl @@ -266,13 +266,13 @@ end function collect_gmi_FisSal2011( mps_filename; - optimizer, - max_cuts_per_round = 1_000_000, - time_limit = 3_600, interval_print_sec=0.1, - silent_solver=true, + max_cuts_per_round = 1_000_000, max_pool_size_mb = 1024, - variant = :fast, + optimizer, + silent_solver=true, + time_limit = 900, + variant = :faster, ) variant in [:subg, :hybr, :fast, :faster] || error("unknown variant: $variant") if variant == :subg @@ -292,6 +292,7 @@ function collect_gmi_FisSal2011( interval_large_lp = 50 interval_read_tableau = 1 end + gapcl_best_patience = 2 * interval_large_lp + 1 reset_timer!() initial_time = time() @@ -314,13 +315,14 @@ function collect_gmi_FisSal2011( count_backtrack = 0 count_deterioration = 0 gapcl_best = 0 + gapcl_best_history = CircularBuffer{Float64}(gapcl_best_patience) gapcl_curr = 0 last_print_time = 0 multipliers_best = nothing multipliers_curr = nothing obj_best = nothing obj_curr = nothing - obj_hist = [] + obj_hist = CircularBuffer{Float64}(100) obj_initial = nothing pool = nothing pool_cut_age = nothing @@ -347,6 +349,11 @@ function collect_gmi_FisSal2011( # 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) + for (var_idx, var_type) in enumerate(data.var_types) + if var_type in ['B', 'I'] + assert_int(sol_opt[var_idx]) + end + end # Convert to standard form data_s, transforms = convert_to_standard_form(data) @@ -363,13 +370,19 @@ function collect_gmi_FisSal2011( sol_opt_s = forward(transforms, sol_opt) # Assert converted solution is feasible for standard form problem + for (var_idx, var_type) in enumerate(data_s.var_types) + if var_type in ['B', 'I'] + assert_int(sol_opt_s[var_idx]) + end + end 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" for round = 1:max_rounds - should_print_log = false + log_prefix = ' ' + log_should_print = false if round > 1 @timeit "Update objective function" begin @@ -395,15 +408,17 @@ function collect_gmi_FisSal2011( @timeit "Optimize LP (lagrangian)" begin optimize!(model_s) + status = termination_status(model_s) + if status != MOI.OPTIMAL + error("Non-optimal termination status: $status") + end sol_frac = get_x(model_s) obj_curr = objective_value(model_s) push!(obj_hist, obj_curr) - if length(obj_hist) > 100 - popfirst!(obj_hist) - end if obj_best === nothing || obj_curr > obj_best + log_prefix = '*' obj_best = obj_curr multipliers_best = multipliers_curr end @@ -412,9 +427,7 @@ function collect_gmi_FisSal2011( end gapcl_curr = gapcl(obj_curr) gapcl_best = gapcl(obj_best) - if termination_status(model_s) != MOI.OPTIMAL - error("Non-optimal termination status") - end + push!(gapcl_best_history, gapcl_best) @timeit "Update μ" begin if variant in [:subg, :hybr] @@ -508,10 +521,38 @@ function collect_gmi_FisSal2011( end end end + + @timeit "Prune the 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_curr[idx_keep] + multipliers_best = multipliers_best[idx_keep] + pool_cut_age = pool_cut_age[idx_keep] + end + pool_size_mb = Base.summarysize(pool) / 1024^2 + end + end end if mod(round - 1, interval_large_lp) == 0 || round == max_rounds - should_print_log = true + log_should_print = true @timeit "Update multipliers (large LP)" begin selected_idx = [] selected_contrs = [] @@ -519,6 +560,10 @@ function collect_gmi_FisSal2011( @timeit "Optimize LP (extended)" begin set_objective_function(model_s, orig_obj_s) optimize!(model_s) + status = termination_status(model_s) + if status != MOI.OPTIMAL + error("Non-optimal termination status: $status") + end obj_curr = objective_value(model_s) sol_frac = get_x(model_s) end @@ -568,38 +613,13 @@ function collect_gmi_FisSal2011( @timeit "Update best" begin if obj_curr > obj_best + log_prefix = '*' obj_best = obj_curr multipliers_best = multipliers_curr end gapcl_curr = gapcl(obj_curr) gapcl_best = gapcl(obj_best) - end - - @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 + push!(gapcl_best_history, gapcl_best) end @timeit "Delete all cut constraints" begin @@ -609,14 +629,19 @@ function collect_gmi_FisSal2011( else @timeit "Update multipliers (subgradient)" begin subgrad = (pool.lb' - (sol_frac' * pool.lhs))' - λ = μ * (obj_mip - obj_curr) / norm(subgrad)^2 + subgrad_norm_sq = norm(subgrad)^2 + if subgrad_norm_sq < 1e-10 + λ = 0 + else + λ = μ * (obj_mip - obj_curr) / subgrad_norm_sq + end multipliers_curr = max.(0, multipliers_curr .+ λ * subgrad) end end if round == 1 @printf( - "%8s %10s %9s %9s %9s %9s %9s %4s %8s %8s %8s\n", + " %8s %10s %9s %9s %9s %9s %9s %4s %8s %8s %8s\n", "round", "obj", "cl_curr", @@ -632,13 +657,14 @@ function collect_gmi_FisSal2011( end if time() - last_print_time > interval_print_sec - should_print_log = true + log_should_print = true end - if should_print_log + if log_should_print last_print_time = time() @printf( - "%8d %10.3e %9.2e %9.2e %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n", + "%c %8d %10.3e %9.2e %9.2e %9d %9.2f %9d %4d %8.2e %8.2e %8.2e\n", + log_prefix, round, obj_curr, gapcl_curr, @@ -653,6 +679,13 @@ function collect_gmi_FisSal2011( ) end + if length(gapcl_best_history) >= gapcl_best_patience + if gapcl_best <= gapcl_best_history[1] + @info "No gap closure improvement. Stopping." + break + end + end + elapsed_time = time() - initial_time if elapsed_time > time_limit @info "Time limit exceeded. Stopping." diff --git a/src/Cuts/tableau/numerics.jl b/src/Cuts/tableau/numerics.jl index 15c38e4..da2da15 100644 --- a/src/Cuts/tableau/numerics.jl +++ b/src/Cuts/tableau/numerics.jl @@ -41,6 +41,8 @@ function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol = ub = cuts.ub[i] lb = cuts.lb[i] if (val >= ub) || (val <= lb) + @show cuts.lhs[i, :] + @show x[cuts.lhs[i, :].nzind] throw( ErrorException( "inequality $i cuts off integer solution ($lb <= $val <= $ub)", @@ -49,3 +51,10 @@ function assert_does_not_cut_off(cuts::ConstraintSet, x::Vector{Float64}; tol = end end end + +function assert_int(x::Float64, tol=1e-6) + fx = frac(x) + if min(fx, 1 - fx) >= tol + throw(ErrorException("Number must be integer: $x")) + end +end