diff --git a/Project.toml b/Project.toml index 727c1bb..b2f4fd0 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,8 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +KLU = "ef3ab10e-7fda-4108-b977-705223b18434" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" diff --git a/src/MIPLearn.jl b/src/MIPLearn.jl index 99f1434..97c6627 100644 --- a/src/MIPLearn.jl +++ b/src/MIPLearn.jl @@ -35,6 +35,7 @@ include("utils/benchmark.jl") include("utils/parse.jl") include("bb/BB.jl") +include("cuts/Cuts.jl") function __init__() copy!(miplearn, pyimport("miplearn")) diff --git a/src/cuts/Cuts.jl b/src/cuts/Cuts.jl new file mode 100644 index 0000000..b5826dc --- /dev/null +++ b/src/cuts/Cuts.jl @@ -0,0 +1,14 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +module Cuts + +include("structs.jl") +include("moi.jl") +include("transform.jl") +include("tableau.jl") +include("gmi.jl") +include("collect.jl") + +end # module diff --git a/src/cuts/collect.jl b/src/cuts/collect.jl new file mode 100644 index 0000000..0957f4d --- /dev/null +++ b/src/cuts/collect.jl @@ -0,0 +1,201 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +import ..Hdf5Sample + +using OrderedCollections + +function collect_gmi( + mps_filename; + optimizer, + max_rounds=10, + max_cuts_per_round=100, +) + @info mps_filename + reset_timer!() + + # Open HDF5 file + h5_filename = replace(mps_filename, ".mps.gz" => ".h5") + h5 = Hdf5Sample(h5_filename) + + # Read optimal solution + sol_opt_dict = Dict( + zip( + h5.get_array("static_var_names"), + convert(Array{Float64}, h5.get_array("mip_var_values")), + ) + ) + + # Read optimal value + obj_mip = h5.get_scalar("mip_lower_bound") + if obj_mip === nothing + obj_mip = h5.get_scalar("mip_obj_value") + end + obj_lp = nothing + h5.file.close() + + # Define relative MIP gap + gap(v) = 100 * abs(obj_mip - v) / abs(v) + + # Initialize stats + stats_obj = [] + stats_gap = [] + stats_ncuts = [] + stats_time_convert = 0 + stats_time_solve = 0 + stats_time_select = 0 + stats_time_tableau = 0 + stats_time_gmi = 0 + all_cuts = nothing + + # Read problem + model = read_from_file(mps_filename) + + for round in 1:max_rounds + @info "Round $(round)..." + + stats_time_convert = @elapsed 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 all(data.constr_lb .- 1e-3 .<= data.constr_lhs * sol_opt) + @assert all(data.constr_lhs * sol_opt .<= data.constr_ub .+ 1e-3) + + # Convert to standard form + data_s, transforms = convert_to_standard_form(data) + model_s = to_model(data_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 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 + push!(stats_obj, obj) + push!(stats_gap, gap(obj)) + push!(stats_ncuts, 0) + end + if termination_status(model_s) != MOI.OPTIMAL + return + end + + # Select tableau rows + basis = get_basis(model_s) + 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, + ) + end + + # Compute selected tableau rows + stats_time_tableau += @elapsed begin + 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 + end + + # Compute GMI cuts + stats_time_gmi += @elapsed begin + 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 + + # Abort if no cuts are left + if length(cuts_s.lb) == 0 + @info "No cuts generated. Aborting." + break + end + end + + # Add GMI cuts to original problem + cuts = backwards(transforms, cuts_s) + assert_does_not_cut_off(cuts, sol_opt) + constrs = add_constraint_set(model, cuts) + + # Optimize original form + set_optimizer(model, optimizer) + undo_relax = relax_integrality(model) + optimize!(model) + obj = objective_value(model) + push!(stats_obj, obj) + 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 + ] + drop = findall(useful .== false) + keep = findall(useful .== true) + delete.(model, constrs[drop]) + if all_cuts === nothing + all_cuts = cuts + else + all_cuts.lhs = [all_cuts.lhs; cuts.lhs[keep, :]] + all_cuts.lb = [all_cuts.lb; cuts.lb[keep]] + all_cuts.lb = [all_cuts.lb; cuts.lb[keep]] + end + push!(stats_ncuts, length(all_cuts.lb)) + + undo_relax() + end + + # Store cuts + if all_cuts !== nothing + @info "Storing $(length(all_cuts.ub)) GMI cuts..." + h5 = Hdf5Sample(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.file.close() + end + + return OrderedDict( + "instance" => mps_filename, + "max_rounds" => max_rounds, + "rounds" => length(stats_obj) - 1, + "time_convert" => stats_time_convert, + "time_solve" => stats_time_solve, + "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, + ) +end + +export collect_gmi \ No newline at end of file diff --git a/src/cuts/gmi.jl b/src/cuts/gmi.jl new file mode 100644 index 0000000..cdfe1a8 --- /dev/null +++ b/src/cuts/gmi.jl @@ -0,0 +1,91 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using SparseArrays + +@inline frac(x::Float64) = x - floor(x) + +function select_gmi_rows(data, basis, x; max_rows=10, atol=0.001) + candidate_rows = [ + r + for r in 1:length(basis.var_basic) + if (data.var_types[basis.var_basic[r]] != 'C') && (frac(x[basis.var_basic[r]]) > atol) + ] + candidate_vals = frac.(x[basis.var_basic[candidate_rows]]) + score = abs.(candidate_vals .- 0.5) + perm = sortperm(score) + return [candidate_rows[perm[i]] for i in 1:min(length(perm), max_rows)] +end + +function compute_gmi( + data::ProblemData, + tableau::Tableau, + tol=1e-8, +)::ConstraintSet + nrows, ncols = size(tableau.lhs) + ub = Float64[Inf for _ in 1:nrows] + lb = Float64[0.999 for _ in 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 in 1:nnz(tableau.lhs) + i::Int = tableau_I[k] + v::Float64 = 0.0 + alpha_j = frac(tableau_V[k]) + beta = frac(tableau.rhs[i]) + if data.var_types[i] == "C" + if alpha_j >= 0 + v = alpha_j / beta + else + v = alpha_j / (1 - beta) + end + else + if alpha_j <= beta + v = alpha_j / beta + else + v = (1 - alpha_j) / (1 - beta) + end + end + if abs(v) > tol + push!(lhs_I, i) + push!(lhs_J, tableau_J[k]) + push!(lhs_V, v) + end + end + lhs = sparse(lhs_I, lhs_J, lhs_V, nrows, ncols) + end + return ConstraintSet(; lhs, ub, lb) +end + +function assert_cuts_off( + cuts::ConstraintSet, + x::Vector{Float64}, + tol=1e-6 +) + for i in 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 in 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 \ No newline at end of file diff --git a/src/cuts/moi.jl b/src/cuts/moi.jl new file mode 100644 index 0000000..34c6eb1 --- /dev/null +++ b/src/cuts/moi.jl @@ -0,0 +1,177 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using JuMP + +function ProblemData(model::Model)::ProblemData + vars = all_variables(model) + + # Objective function + obj = objective_function(model) + obj = [v ∈ keys(obj.terms) ? obj.terms[v] : 0.0 for v in vars] + + # Variable types, lower bounds and upper bounds + var_lb = [is_binary(v) ? 0.0 : has_lower_bound(v) ? lower_bound(v) : -Inf for v in vars] + var_ub = [is_binary(v) ? 1.0 : has_upper_bound(v) ? upper_bound(v) : Inf for v in vars] + var_types = [is_binary(v) || is_integer(v) ? 'I' : 'C' for v in vars] + var_names = [name(v) for v in vars] + + # Constraints + constr_lb = Float64[] + constr_ub = Float64[] + constr_lhs_rows = Int[] + constr_lhs_cols = Int[] + constr_lhs_values = Float64[] + constr_index = 1 + for (ftype, stype) in list_of_constraint_types(model) + for constr in all_constraints(model, ftype, stype) + cset = MOI.get(constr.model.moi_backend, MOI.ConstraintSet(), constr.index) + cf = MOI.get( + constr.model.moi_backend, + MOI.ConstraintFunction(), + constr.index, + ) + if ftype == VariableRef + var_idx = cf.value + if stype == MOI.Integer || stype == MOI.ZeroOne + # nop + elseif stype == MOI.EqualTo{Float64} + var_lb[var_idx] = max(var_lb[var_idx], cset.value) + var_ub[var_idx] = min(var_ub[var_idx], cset.value) + elseif stype == MOI.LessThan{Float64} + var_ub[var_idx] = min(var_ub[var_idx], cset.upper) + elseif stype == MOI.GreaterThan{Float64} + var_lb[var_idx] = max(var_lb[var_idx], cset.lower) + elseif stype == MOI.Interval{Float64} + var_lb[var_idx] = max(var_lb[var_idx], cset.lower) + var_ub[var_idx] = min(var_ub[var_idx], cset.upper) + else + error("Unsupported set: $stype") + end + elseif ftype == AffExpr + if stype == MOI.EqualTo{Float64} + push!(constr_lb, cset.value) + push!(constr_ub, cset.value) + elseif stype == MOI.LessThan{Float64} + push!(constr_lb, -Inf) + push!(constr_ub, cset.upper) + elseif stype == MOI.GreaterThan{Float64} + push!(constr_lb, cset.lower) + push!(constr_ub, Inf) + elseif stype == MOI.Interval{Float64} + push!(constr_lb, cset.lower) + push!(constr_ub, cset.upper) + else + error("Unsupported set: $stype") + end + for term in cf.terms + push!(constr_lhs_cols, term.variable.value) + push!(constr_lhs_rows, constr_index) + push!(constr_lhs_values, term.coefficient) + end + constr_index += 1 + else + error("Unsupported constraint type: ($ftype, $stype)") + end + end + end + + n = length(vars) + m = constr_index - 1 + constr_lhs = sparse( + constr_lhs_rows, + constr_lhs_cols, + constr_lhs_values, + m, + n, + ) + + @assert length(obj) == n + @assert length(var_lb) == n + @assert length(var_ub) == n + @assert length(var_types) == n + @assert length(var_names) == n + @assert length(constr_lb) == m + @assert length(constr_ub) == m + + return ProblemData( + obj_offset=0.0; + obj, + constr_lb, + constr_ub, + constr_lhs, + var_lb, + var_ub, + var_types, + var_names + ) +end + +function to_model(data::ProblemData, tol=1e-6)::Model + model = Model() + + # Variables + nvars = length(data.obj) + @variable(model, x[1:nvars]) + for i = 1:nvars + set_name(x[i], data.var_names[i]) + if data.var_types[i] == 'B' + set_binary(x[i]) + elseif data.var_types[i] == 'I' + set_integer(x[i]) + end + if isfinite(data.var_lb[i]) + set_lower_bound(x[i], data.var_lb[i]) + end + if isfinite(data.var_ub[i]) + set_upper_bound(x[i], data.var_ub[i]) + end + set_objective_coefficient(model, x[i], data.obj[i]) + end + + # Constraints + lhs = data.constr_lhs * x + for (j, lhs_expr) in enumerate(lhs) + lb = data.constr_lb[j] + ub = data.constr_ub[j] + if abs(lb - ub) < tol + @constraint(model, lb == lhs_expr) + elseif isfinite(lb) && !isfinite(ub) + @constraint(model, lb <= lhs_expr) + elseif !isfinite(lb) && isfinite(ub) + @constraint(model, lhs_expr <= ub) + else + @constraint(model, lb <= lhs_expr <= ub) + end + end + + return model +end + +function add_constraint_set(model::JuMP.Model, cs::ConstraintSet) + vars = all_variables(model) + nrows, _ = size(cs.lhs) + constrs = [] + for i in 1:nrows + c = nothing + if isinf(cs.ub[i]) + c = @constraint(model, cs.lb[i] <= dot(cs.lhs[i, :], vars)) + elseif isinf(cs.lb[i]) + c = @constraint(model, dot(cs.lhs[i, :], vars) <= cs.ub[i]) + else + c = @constraint(model, cs.lb[i] <= dot(cs.lhs[i, :], vars) <= cs.ub[i]) + end + push!(constrs, c) + end + return constrs +end + +function set_warm_start(model::JuMP.Model, x::Vector{Float64}) + vars = all_variables(model) + for (i, xi) in enumerate(x) + set_start_value(vars[i], xi) + end +end + +export to_model, ProblemData, add_constraint_set, set_warm_start \ No newline at end of file diff --git a/src/cuts/structs.jl b/src/cuts/structs.jl new file mode 100644 index 0000000..55200ad --- /dev/null +++ b/src/cuts/structs.jl @@ -0,0 +1,39 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using SparseArrays + +Base.@kwdef mutable struct ProblemData + obj::Vector{Float64} + obj_offset::Float64 + constr_lb::Vector{Float64} + constr_ub::Vector{Float64} + constr_lhs::SparseMatrixCSC + var_lb::Vector{Float64} + var_ub::Vector{Float64} + var_types::Vector{Char} + var_names::Vector{String} +end + +Base.@kwdef mutable struct Tableau + obj + lhs + rhs + z +end + +Base.@kwdef mutable struct Basis + var_basic + var_nonbasic + constr_basic + constr_nonbasic +end + +Base.@kwdef mutable struct ConstraintSet + lhs::SparseMatrixCSC + ub::Vector{Float64} + lb::Vector{Float64} +end + +export ProblemData, Tableau, Basis, ConstraintSet diff --git a/src/cuts/tableau.jl b/src/cuts/tableau.jl new file mode 100644 index 0000000..6498342 --- /dev/null +++ b/src/cuts/tableau.jl @@ -0,0 +1,130 @@ +# MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization +# Copyright (C) 2020-2022, UChicago Argonne, LLC. All rights reserved. +# Released under the modified BSD license. See COPYING.md for more details. + +using KLU +using TimerOutputs + +function get_basis(model::JuMP.Model)::Basis + var_basic = Int[] + var_nonbasic = Int[] + constr_basic = Int[] + constr_nonbasic = Int[] + + # Variables + for (i, var) in enumerate(all_variables(model)) + bstatus = MOI.get(model, MOI.VariableBasisStatus(), var) + if bstatus == MOI.BASIC + push!(var_basic, i) + elseif bstatus == MOI.NONBASIC_AT_LOWER + push!(var_nonbasic, i) + else + error("Unknown basis status: $bstatus") + end + end + + # Constraints + constr_index = 1 + for (ftype, stype) in list_of_constraint_types(model) + for constr in all_constraints(model, ftype, stype) + if ftype == VariableRef + # nop + elseif ftype == AffExpr + bstatus = MOI.get(model, MOI.ConstraintBasisStatus(), constr) + if bstatus == MOI.BASIC + push!(constr_basic, constr_index) + elseif bstatus == MOI.NONBASIC + push!(constr_nonbasic, constr_index) + else + error("Unknown basis status: $bstatus") + end + constr_index += 1 + else + error("Unsupported constraint type: ($ftype, $stype)") + end + end + end + + return Basis(; var_basic, var_nonbasic, constr_basic, constr_nonbasic) +end + +function get_x(model::JuMP.Model) + return JuMP.value.(all_variables(model)) +end + +function compute_tableau( + data::ProblemData, + basis::Basis, + x::Vector{Float64}; + rows::Union{Vector{Int},Nothing}=nothing, + tol=1e-8 +)::Tableau + @timeit "Split data" begin + nrows, ncols = size(data.constr_lhs) + lhs_slacks = sparse(I, nrows, nrows) + lhs_b = [data.constr_lhs[:, basis.var_basic] lhs_slacks[:, basis.constr_basic]] + obj_b = [data.obj[basis.var_basic]; zeros(length(basis.constr_basic))] + if rows === nothing + rows = 1:nrows + end + end + + @timeit "Factorize basis matrix" begin + factor = klu(sparse(lhs_b')) + end + + @timeit "Compute tableau LHS" begin + tableau_lhs_I = Int[] + tableau_lhs_J = Int[] + tableau_lhs_V = Float64[] + for k in 1:length(rows) + @timeit "Prepare inputs" begin + i = rows[k] + e = zeros(nrows) + e[i] = 1.0 + end + @timeit "Solve" begin + sol = factor \ e + end + @timeit "Multiply" begin + row = sol' * data.constr_lhs + end + @timeit "Sparsify & copy" begin + for (j, v) in enumerate(row) + if abs(v) < tol + continue + end + push!(tableau_lhs_I, k) + push!(tableau_lhs_J, j) + push!(tableau_lhs_V, v) + end + end + end + tableau_lhs = sparse( + tableau_lhs_I, + tableau_lhs_J, + tableau_lhs_V, + length(rows), + ncols, + ) + end + + @timeit "Compute tableau RHS" begin + tableau_rhs = [x[basis.var_basic]; zeros(length(basis.constr_basic))][rows] + end + + @timeit "Compute tableau objective row" begin + sol = factor \ obj_b + tableau_obj = -data.obj' + sol' * data.constr_lhs + tableau_obj[abs.(tableau_obj).