parent
db6456dbaa
commit
d8b80f00ca
GPG Key ID:
0DA8E4B9E1109DCA
@ -0,0 +1,14 @@
|
||||
module Cuts
|
||||
|
||||
import ..to_str_array
|
||||
|
||||
include("tableau/structs.jl")
|
||||
|
||||
include("blackbox/cplex.jl")
|
||||
include("tableau/collect.jl")
|
||||
include("tableau/gmi.jl")
|
||||
include("tableau/moi.jl")
|
||||
include("tableau/tableau.jl")
|
||||
include("tableau/transform.jl")
|
||||
|
||||
end # module
|
||||
@ -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 ..H5File
|
||||
|
||||
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 = H5File(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 = 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.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
|
||||
@ -0,0 +1,92 @@
|
||||
# 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
|
||||
using TimerOutputs
|
||||
|
||||
@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
|
||||
@ -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
|
||||
@ -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
|
||||
@ -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).<tol] .= 0
|
||||
end
|
||||
|
||||
return Tableau(
|
||||
obj=tableau_obj,
|
||||
lhs=tableau_lhs,
|
||||
rhs=tableau_rhs,
|
||||
z=dot(data.obj, x),
|
||||
)
|
||||
end
|
||||
|
||||
export get_basis, get_x, compute_tableau
|
||||
@ -0,0 +1,314 @@
|
||||
# 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 LinearAlgebra
|
||||
using TimerOutputs
|
||||
|
||||
abstract type Transform end
|
||||
|
||||
function _isbounded(x)
|
||||
isfinite(x) || return false
|
||||
abs(x) < 1e15 || return false
|
||||
return true
|
||||
end
|
||||
|
||||
function backwards!(transforms::Vector{Transform}, m::ConstraintSet; tol=1e-8)
|
||||
for t in reverse(transforms)
|
||||
backwards!(t, m)
|
||||
end
|
||||
for (idx, val) in enumerate(m.lhs.nzval)
|
||||
if abs(val) < tol
|
||||
m.lhs.nzval[idx] = 0
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
function backwards(transforms::Vector{Transform}, m::ConstraintSet; tol=1e-8)
|
||||
m2 = deepcopy(m)
|
||||
backwards!(transforms, m2; tol)
|
||||
return m2
|
||||
end
|
||||
|
||||
function forward(transforms::Vector{Transform}, p::Vector{Float64})::Vector{Float64}
|
||||
for t in transforms
|
||||
p = forward(t, p)
|
||||
end
|
||||
return p
|
||||
end
|
||||
|
||||
# -----------------------------------------------------------------------------
|
||||
|
||||
Base.@kwdef mutable struct ShiftVarLowerBoundsToZero <: Transform
|
||||
lb::Vector{Float64} = []
|
||||
end
|
||||
|
||||
function forward!(t::ShiftVarLowerBoundsToZero, data::ProblemData)
|
||||
t.lb = copy(data.var_lb)
|
||||
data.obj_offset += dot(data.obj, t.lb)
|
||||
data.var_ub -= t.lb
|
||||
data.var_lb -= t.lb
|
||||
data.constr_lb -= data.constr_lhs * t.lb
|
||||
data.constr_ub -= data.constr_lhs * t.lb
|
||||
end
|
||||
|
||||
function backwards!(t::ShiftVarLowerBoundsToZero, c::ConstraintSet)
|
||||
c.lb += c.lhs * t.lb
|
||||
c.ub += c.lhs * t.lb
|
||||
end
|
||||
|
||||
function forward(t::ShiftVarLowerBoundsToZero, p::Vector{Float64})::Vector{Float64}
|
||||
return p - t.lb
|
||||
end
|
||||
|
||||
# -----------------------------------------------------------------------------
|
||||
|
||||
Base.@kwdef mutable struct MoveVarUpperBoundsToConstrs <: Transform end
|
||||
|
||||
function forward!(t::MoveVarUpperBoundsToConstrs, data::ProblemData)
|
||||
_, ncols = size(data.constr_lhs)
|
||||
data.constr_lhs = [data.constr_lhs; I]
|
||||
data.constr_lb = [data.constr_lb; [-Inf for _ in 1:ncols]]
|
||||
data.constr_ub = [data.constr_ub; data.var_ub]
|
||||
data.var_ub .= Inf
|
||||
end
|
||||
|
||||
function backwards!(::MoveVarUpperBoundsToConstrs, ::ConstraintSet)
|
||||
# nop
|
||||
end
|
||||
|
||||
function forward(t::MoveVarUpperBoundsToConstrs, p::Vector{Float64})::Vector{Float64}
|
||||
return p
|
||||
end
|
||||
|
||||
# -----------------------------------------------------------------------------
|
||||
|
||||
Base.@kwdef mutable struct AddSlackVariables <: Transform
|
||||
M1::SparseMatrixCSC = spzeros(0)
|
||||
M2::Vector{Float64} = []
|
||||
ncols_orig::Int = 0
|
||||
GE::Int = 0
|
||||
LE::Int = 0
|
||||
lhs_ge::SparseMatrixCSC = spzeros(0)
|
||||
lhs_le::SparseMatrixCSC = spzeros(0)
|
||||
rhs_le::Vector{Float64} = []
|
||||
rhs_ge::Vector{Float64} = []
|
||||
end
|
||||
|
||||
function forward!(t::AddSlackVariables, data::ProblemData)
|
||||
nrows, ncols = size(data.constr_lhs)
|
||||
isequality = abs.(data.constr_ub .- data.constr_lb) .< 1e-6
|
||||
eq = [i for i in 1:nrows if isequality[i]]
|
||||
ge = [i for i in 1:nrows if isfinite(data.constr_lb[i]) && !isequality[i]]
|
||||
le = [i for i in 1:nrows if isfinite(data.constr_ub[i]) && !isequality[i]]
|
||||
EQ, GE, LE = length(eq), length(ge), length(le)
|
||||
|
||||
t.M1 = [
|
||||
I spzeros(ncols, GE + LE)
|
||||
data.constr_lhs[ge, :] spzeros(GE, GE + LE)
|
||||
-data.constr_lhs[le, :] spzeros(LE, GE + LE)
|
||||
]
|
||||
t.M2 = [
|
||||
zeros(ncols)
|
||||
data.constr_lb[ge]
|
||||
-data.constr_ub[le]
|
||||
]
|
||||
t.ncols_orig = ncols
|
||||
t.GE, t.LE = GE, LE
|
||||
t.lhs_ge = data.constr_lhs[ge, :]
|
||||
t.lhs_le = data.constr_lhs[le, :]
|
||||
t.rhs_ge = data.constr_lb[ge]
|
||||
t.rhs_le = data.constr_ub[le]
|
||||
|
||||
data.constr_lhs = [
|
||||
data.constr_lhs[eq, :] spzeros(EQ, GE) spzeros(EQ, LE)
|
||||
data.constr_lhs[ge, :] -I spzeros(GE, LE)
|
||||
data.constr_lhs[le, :] spzeros(LE, GE) I
|
||||
]
|
||||
data.obj = [data.obj; zeros(GE + LE)]
|
||||
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 in 1:(GE+LE)]]
|
||||
data.var_types = [data.var_types; ['C' for _ in 1:(GE+LE)]]
|
||||
data.constr_lb = [
|
||||
data.constr_lb[eq]
|
||||
data.constr_lb[ge]
|
||||
data.constr_ub[le]
|
||||
]
|
||||
data.constr_ub = copy(data.constr_lb)
|
||||
end
|
||||
|
||||
function backwards!(t::AddSlackVariables, c::ConstraintSet)
|
||||
c.lb += c.lhs * t.M2
|
||||
c.ub += c.lhs * t.M2
|
||||
c.lhs = (c.lhs*t.M1)[:, 1:t.ncols_orig]
|
||||
end
|
||||
|
||||
function forward(t::AddSlackVariables, x::Vector{Float64})::Vector{Float64}
|
||||
return [
|
||||
x
|
||||
t.lhs_ge * x - t.rhs_ge
|
||||
t.rhs_le - t.lhs_le * x
|
||||
]
|
||||
end
|
||||
|
||||
# -----------------------------------------------------------------------------
|
||||
|
||||
Base.@kwdef mutable struct SplitFreeVars <: Transform
|
||||
F::Int = 0
|
||||
B::Int = 0
|
||||
free::Vector{Int}=[]
|
||||
others::Vector{Int}=[]
|
||||
end
|
||||
|
||||
function forward!(t::SplitFreeVars, data::ProblemData)
|
||||
lhs = data.constr_lhs
|
||||
_, ncols = size(lhs)
|
||||
free = [i for i in 1:ncols if !isfinite(data.var_lb[i]) && !isfinite(data.var_ub[i])]
|
||||
others = [i for i in 1:ncols if isfinite(data.var_lb[i]) || isfinite(data.var_ub[i])]
|
||||
t.F = length(free)
|
||||
t.B = length(others)
|
||||
t.free, t.others = free, others
|
||||
data.obj = [
|
||||
data.obj[others]
|
||||
data.obj[free]
|
||||
-data.obj[free]
|
||||
]
|
||||
data.constr_lhs = [lhs[:, others] lhs[:, free] -lhs[:, free]]
|
||||
data.var_lb = [
|
||||
data.var_lb[others]
|
||||
[0.0 for _ in free]
|
||||
[0.0 for _ in free]
|
||||
]
|
||||
data.var_ub = [
|
||||
data.var_ub[others]
|
||||
[Inf for _ in free]
|
||||
[Inf for _ in free]
|
||||
]
|
||||
data.var_types = [
|
||||
data.var_types[others]
|
||||
data.var_types[free]
|
||||
data.var_types[free]
|
||||
]
|
||||
data.var_names = [
|
||||
data.var_names[others]
|
||||
["$(v)_p" for v in data.var_names[free]]
|
||||
["$(v)_m" for v in data.var_names[free]]
|
||||
]
|
||||
end
|
||||
|
||||
function backwards!(t::SplitFreeVars, c::ConstraintSet)
|
||||
# Convert GE constraints into LE
|
||||
nrows, _ = size(c.lhs)
|
||||
ge = [i for i in 1:nrows if isfinite(c.lb[i])]
|
||||
c.ub[ge], c.lb[ge] = -c.lb[ge], -c.ub[ge]
|
||||
c.lhs[ge, :] *= -1
|
||||
|
||||
# Assert only LE constraints are left (EQ constraints are not supported)
|
||||
@assert all(c.lb .== -Inf)
|
||||
|
||||
# Take minimum (weakest) coefficient
|
||||
B, F = t.B, t.F
|
||||
for i in 1:F
|
||||
c.lhs[:, B + i] = min.(c.lhs[:, B + i], -c.lhs[:, B + F + i])
|
||||
end
|
||||
c.lhs = c.lhs[:, 1:(B+F)]
|
||||
end
|
||||
|
||||
function forward(t::SplitFreeVars, p::Vector{Float64})::Vector{Float64}
|
||||
return [
|
||||
p[t.others]
|
||||
max.(p[t.free], 0)
|
||||
max.(-p[t.free], 0)
|
||||
]
|
||||
end
|
||||
|
||||
# -----------------------------------------------------------------------------
|
||||
|
||||
Base.@kwdef mutable struct FlipUnboundedBelowVars <: Transform
|
||||
flip_idx::Vector{Int} = []
|
||||
end
|
||||
|
||||
function forward!(t::FlipUnboundedBelowVars, data::ProblemData)
|
||||
_, ncols = size(data.constr_lhs)
|
||||
for i in 1:ncols
|
||||
if isfinite(data.var_lb[i])
|
||||
continue
|
||||
end
|
||||
data.obj[i] *= -1
|
||||
data.constr_lhs[:, i] *= -1
|
||||
data.var_lb[i], data.var_ub[i] = -data.var_ub[i], -data.var_lb[i]
|
||||
push!(t.flip_idx, i)
|
||||
end
|
||||
end
|
||||
|
||||
function backwards!(t::FlipUnboundedBelowVars, c::ConstraintSet)
|
||||
for i in t.flip_idx
|
||||
c.lhs[:, i] *= -1
|
||||
end
|
||||
end
|
||||
|
||||
function forward(t::FlipUnboundedBelowVars, p::Vector{Float64})::Vector{Float64}
|
||||
p2 = copy(p)
|
||||
p2[t.flip_idx] *= -1
|
||||
return p2
|
||||
end
|
||||
|
||||
# -----------------------------------------------------------------------------
|
||||
|
||||
function _assert_standard_form(data::ProblemData)
|
||||
# Check sizes
|
||||
nrows, ncols = size(data.constr_lhs)
|
||||
@assert length(data.constr_lb) == nrows
|
||||
@assert length(data.constr_ub) == nrows
|
||||
@assert length(data.obj) == ncols
|
||||
@assert length(data.var_lb) == ncols
|
||||
@assert length(data.var_ub) == ncols
|
||||
@assert length(data.var_names) == ncols
|
||||
@assert length(data.var_types) == ncols
|
||||
|
||||
# Check standard form
|
||||
@assert all(data.var_lb .== 0.0)
|
||||
@assert all(data.var_ub .== Inf)
|
||||
@assert all(data.constr_lb .== data.constr_ub)
|
||||
end
|
||||
|
||||
function convert_to_standard_form!(data::ProblemData)::Vector{Transform}
|
||||
transforms = []
|
||||
function _apply!(t)
|
||||
push!(transforms, t)
|
||||
forward!(t, data)
|
||||
end
|
||||
@timeit "Split free vars" begin
|
||||
_apply!(SplitFreeVars())
|
||||
end
|
||||
@timeit "Flip unbounded-below vars" begin
|
||||
_apply!(FlipUnboundedBelowVars())
|
||||
end
|
||||
@timeit "Shift var lower bounds to zero" begin
|
||||
_apply!(ShiftVarLowerBoundsToZero())
|
||||
end
|
||||
@timeit "Move var upper bounds to constrs" begin
|
||||
_apply!(MoveVarUpperBoundsToConstrs())
|
||||
end
|
||||
@timeit "Add slack vars" begin
|
||||
_apply!(AddSlackVariables())
|
||||
end
|
||||
_assert_standard_form(data)
|
||||
return transforms
|
||||
end
|
||||
|
||||
function convert_to_standard_form(data::ProblemData)::Tuple{ProblemData,Vector{Transform}}
|
||||
data2 = deepcopy(data)
|
||||
transforms = convert_to_standard_form!(data2)
|
||||
return (data2, transforms)
|
||||
end
|
||||
|
||||
export convert_to_standard_form!,
|
||||
convert_to_standard_form,
|
||||
forward!,
|
||||
backwards!,
|
||||
backwards,
|
||||
AddSlackVariables,
|
||||
SplitFreeVars,
|
||||
forward
|
||||
Binary file not shown.
Loading…
Reference in new issue