Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
Manifest.toml
.DS_Store
114 changes: 111 additions & 3 deletions src/MathOptLazy.jl
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@ struct _LazyData{F<:MOI.AbstractScalarFunction,S<:MOI.AbstractScalarSet}
data::Vector{Tuple{F,S}}
status::Vector{_LAZY_CONSTRAINT_STATUS}
index::Vector{MOI.ConstraintIndex{F,S}}
stale::Vector{Int}

function _LazyData{
F,
Expand All @@ -99,6 +100,7 @@ struct _LazyData{F<:MOI.AbstractScalarFunction,S<:MOI.AbstractScalarSet}
Tuple{F,S}[],
_LAZY_CONSTRAINT_STATUS[],
MOI.ConstraintIndex{F,S}[],
Int[],
)
end
end
Expand All @@ -122,7 +124,7 @@ struct Algorithm <: MOI.AbstractOptimizerAttribute end
abstract type AbstractAlgorithm end

"""
Iterative()
Iterative(; stale_solve_limit::Int = typemax(Int), prune_batch_size::Int = 1)

This algorithm iteratively solves a sequence of problems that iteratively add
violated lazy constraints to the main problem.
Expand All @@ -145,8 +147,30 @@ warm-start from the previous basis after new constraints are added, whereas
every iteration of the second phase solves a mixed-integer program from
scratch. The second phase is still needed because a lazy constraint can be
violated by an integer solution even if the relaxation satisfies it.

## Stale constraints

An active lazy constraint is stale if it has been non-binding for at least
`stale_solve_limit` consecutive solves (including internal iterative solves).
Stale constraints are pruned before solving if there are at least
`prune_batch_size` stale constraints.
"""
struct Iterative <: AbstractAlgorithm end
struct Iterative <: AbstractAlgorithm
stale_solve_limit::Int
prune_batch_size::Int

function Iterative(;
stale_solve_limit::Int = typemax(Int),
prune_batch_size::Int = typemax(Int),
)
if stale_solve_limit < 1
throw(ArgumentError("stale_solve_limit must be at least 1."))
elseif prune_batch_size < 1
throw(ArgumentError("prune_batch_size must be at least 1."))
end
return new(stale_solve_limit, prune_batch_size)
end
end

"""
Callback()
Expand Down Expand Up @@ -489,6 +513,7 @@ function MOI.add_constraint(
push!(data.data, (f, s.set))
push!(data.status, _kLAZY_CONSTRAINT_INACTIVE)
push!(data.index, MOI.ConstraintIndex{F,S}(0))
push!(data.stale, 0)
return MOI.ConstraintIndex{F,LazyScalarSet{S}}(length(data.data))
end

Expand Down Expand Up @@ -655,6 +680,87 @@ function MOI.get(
return count(==(_kLAZY_CONSTRAINT_ACTIVE), _data(model, F, S).status)
end

### MathOptLazy.prune_stale_constraints!

function _is_stale(algorithm::Iterative, data::_LazyData, i::Int)
return data.status[i] == _kLAZY_CONSTRAINT_ACTIVE &&
data.stale[i] >= algorithm.stale_solve_limit
end

function _number_stale(algorithm::Iterative, data::_LazyData)
return count(i -> _is_stale(algorithm, data, i), eachindex(data.status))
end

function _prune_stale!(
model::Optimizer,
algorithm::Iterative,
data::_LazyData{F,S},
) where {F,S}
constraints_pruned = 0
for i in reverse(eachindex(data.status))
if _is_stale(algorithm, data, i)
MOI.delete(model.inner, data.index[i])

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We're better to cache the vector of indices that need to be pruned, and then call one MOI.delete(model.inner, list_to_be_deleted), otherwise we might call MOI.delete a lot of times

data.index[i] = MOI.ConstraintIndex{F,S}(0)
data.status[i] = _kLAZY_CONSTRAINT_INACTIVE
data.stale[i] = 0
constraints_pruned += 1
end
end
return constraints_pruned
end

function _prune_stale!(model::Optimizer, algorithm::Iterative)
return sum(
d -> _prune_stale!(model, algorithm, d),
values(model.lazy);
init = 0,
)
end

function _maybe_prune_stale!(model::Optimizer, algorithm::Iterative)
n = sum(d -> _number_stale(algorithm, d), values(model.lazy); init = 0)
if n < algorithm.prune_batch_size
return
end
_prune_stale!(model, algorithm)
if !model.silent
println("[MathOptLazy] pruned $(n) stale constraints")
end
return
end

_is_binding(y, s::MOI.LessThan) = y >= s.upper - 1e-2
_is_binding(y, s::MOI.GreaterThan) = y <= s.lower + 1e-2
_is_binding(y, s::MOI.Interval) = y <= s.lower + 1e-2 || y >= s.upper - 1e-2

# For sets whose slackness we can't measure, assume they bind and don't prune.
_is_binding(y, ::MOI.AbstractScalarSet) = true

function _update_stale!(model::Optimizer, data::_LazyData)
for (i, (_, s)) in enumerate(data.data)
if data.status[i] != _kLAZY_CONSTRAINT_ACTIVE
continue
end
y = MOI.get(model.inner, MOI.ConstraintPrimal(), data.index[i])
if _is_binding(y, s)
data.stale[i] = 0
else
data.stale[i] += 1
end
end
return
end

function _update_stale!(model::Optimizer)
if MOI.get(model, MOI.PrimalStatus()) != MOI.FEASIBLE_POINT
return
end
for data in values(model.lazy)
_update_stale!(model, data)
end
return
end

### MOI.optimize!

function _try_get(model::Optimizer, attr::MOI.AbstractModelAttribute, default)
Expand Down Expand Up @@ -690,7 +796,7 @@ end

### MathOptLazy.Iterative

function _optimize!(model::Optimizer, ::Iterative)
function _optimize!(model::Optimizer, algorithm::Iterative)
if (undo = _relax_integrality(model.inner)) !== nothing
if !model.silent
println("[MathOptLazy] relaxing binary and integer variables")
Expand All @@ -701,6 +807,7 @@ function _optimize!(model::Optimizer, ::Iterative)
end
undo()
end
_maybe_prune_stale!(model, algorithm)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting question whether this should go at the entry to _optimize!, here, or in both places. I don't know the answer.

_iterate(model; start = true)
return
end
Expand All @@ -722,6 +829,7 @@ function _iterate(model::Optimizer; start::Bool)
println("[MathOptLazy] solving current subproblem\n")
end
_optimize_inner!(model)
_update_stale!(model)
if MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
# The problem is unbounded, but it might not be if we add more
# constraints.
Expand Down
134 changes: 134 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -855,6 +855,140 @@ function test_result_attributes_BarrierIterations()
return
end

function _stale_model()
model = MathOptLazy.Optimizer(HiGHS.Optimizer)
MOI.set(model, MOI.Silent(), true)
x = MOI.add_variables(model, 2)
MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
MOI.add_constraint(model, 1.0 * x[1] + 1.0 * x[2], MOI.LessThan(10.0))
for xi in x
set = MathOptLazy.LazyScalarSet(MOI.LessThan(1.0))
MOI.add_constraint(model, 1.0 * xi, set)
end
F, S = MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}
return model, x, MathOptLazy._data(model, F, S)
end

function _stale_solve(model, sense, f)
MOI.set(model, MOI.ObjectiveSense(), sense)
MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
MOI.optimize!(model)
@test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
return
end

function _set_iterative(model; kwargs...)
algorithm = MathOptLazy.Iterative(; kwargs...)
MOI.set(model, MathOptLazy.Algorithm(), algorithm)
return
end

function test_stale_iterative_kwargs()
algorithm = MathOptLazy.Iterative()
@test algorithm.stale_solve_limit == typemax(Int)
@test algorithm.prune_batch_size == 1
algorithm =
MathOptLazy.Iterative(; stale_solve_limit = 3, prune_batch_size = 4)
@test algorithm.stale_solve_limit == 3
@test algorithm.prune_batch_size == 4
@test_throws ArgumentError MathOptLazy.Iterative(; stale_solve_limit = 0)
@test_throws ArgumentError MathOptLazy.Iterative(; prune_batch_size = 0)
model = MathOptLazy.Optimizer(HiGHS.Optimizer)
MOI.set(model, MathOptLazy.Algorithm(), algorithm)
MOI.empty!(model)
@test MOI.get(model, MathOptLazy.Algorithm()) == algorithm
return
end

function test_is_binding()
s = MOI.LessThan(1.0)
@test MathOptLazy._is_binding(1.0, s)
@test !MathOptLazy._is_binding(0.99, s)
s = MOI.GreaterThan(1.0)
@test MathOptLazy._is_binding(1.0, s)
@test !MathOptLazy._is_binding(1.01, s)
s = MOI.Interval(1.0, 2.0)
@test MathOptLazy._is_binding(1.0, s)
@test MathOptLazy._is_binding(2.0, s)
@test !MathOptLazy._is_binding(1.5, s)
@test MathOptLazy._is_binding(1.5, MOI.EqualTo(1.0))
@test MathOptLazy._is_binding(0.5, MOI.ZeroOne())
return
end

function test_stale_counter()
model, x, data = _stale_model()
@test data.stale == [0, 0]
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] + 1.0 * x[2])
@test data.stale == [0, 0]
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] - 1.0 * x[2])
@test data.stale == [0, 1]
MOI.optimize!(model)
@test data.stale == [0, 2]
_stale_solve(model, MOI.MIN_SENSE, 1.0 * x[1] + 1.0 * x[2])
@test data.stale == [1, 3]
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] + 1.0 * x[2])
@test data.stale == [0, 0]
return
end

function test_stale_counter_no_feasible_point()
model, x, data = _stale_model()
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] + 1.0 * x[2])
_stale_solve(model, MOI.MIN_SENSE, 1.0 * x[1] + 1.0 * x[2])
@test data.stale == [1, 1]
MOI.add_constraint(model, 1.0 * x[1], MOI.GreaterThan(2.0))
MOI.optimize!(model)
@test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
@test data.stale == [1, 1]
return
end

function test_prune()
model, x, data = _stale_model()
_set_iterative(model; stale_solve_limit = 2, prune_batch_size = 2)
F, S = MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}
attr =
MathOptLazy.NumberOfConstraintsActive{F,MathOptLazy.LazyScalarSet{S}}()
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] + 1.0 * x[2])
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] - 1.0 * x[2])
MOI.optimize!(model)
@test data.stale == [0, 2]
MOI.optimize!(model)
@test data.stale == [0, 3]
@test MOI.get(model, attr) == 2
_stale_solve(model, MOI.MIN_SENSE, 1.0 * x[1] + 1.0 * x[2])
MOI.optimize!(model)
@test data.stale == [2, 5]
@test MOI.get(model, attr) == 2
@test MOI.get(model.inner, MOI.NumberOfConstraints{F,S}()) == 3
MOI.optimize!(model)
@test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
@test MOI.get(model, attr) == 0
@test MOI.get(model.inner, MOI.NumberOfConstraints{F,S}()) == 1
@test all(==(MathOptLazy._kLAZY_CONSTRAINT_INACTIVE), data.status)
@test all(==(MOI.ConstraintIndex{F,S}(0)), data.index)
@test data.stale == [0, 0]
return
end

function test_prune_readded()
model, x, data = _stale_model()
_set_iterative(model; stale_solve_limit = 1)
F, S = MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}
attr =
MathOptLazy.NumberOfConstraintsActive{F,MathOptLazy.LazyScalarSet{S}}()
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] + 1.0 * x[2])
_stale_solve(model, MOI.MIN_SENSE, 1.0 * x[1] + 1.0 * x[2])
@test data.stale == [1, 1]
@test MOI.get(model, attr) == 0
_stale_solve(model, MOI.MAX_SENSE, 1.0 * x[1] + 1.0 * x[2])
@test all(==(MathOptLazy._kLAZY_CONSTRAINT_ACTIVE), data.status)
@test data.stale == [0, 0]
@test MOI.get(model, attr) == 2
return
end

end # TestMathOptLazy

TestMathOptLazy.runtests()
Loading