I am using JuMP to solve many (~tens to hundreds of thousands) small integer-feasibility problems (variants of a subset sum problem). Each individual problem is parameterized by different integers but the structure is quite similar.
My initial implementation (pseudo-code here; full implementation in details below) ran something like this:
function subset_sum_variant(parameters, opt=optimizer_with_attributes(HiGHS.Optimizer, "presolve"=>"off"))
# Step 1. instantiate a JuMP model to hold the problem
m = Model(opt)
# Step 2. Set up model, using `parameters` via `@variable` and `@constraint`
...
# Step 3. check if feasible
optimize!(m)
termination_status(m) == INFEASIBLE
end
To my surprise, I found that maybe ~40% of the time was spent in step 1, which I had thought to be near-free.
Since I’m solving these problems repeatedly and Step 1 is the same each time, I was wondering whether there’s a way to set up a Model outside to reuse in the function?
On Slack, I asked whether something like the following would be the intended way to do that:
const m = Model(optimizer_with_attributes(HiGHS.Optimizer, "presolve"=>"off"))
function solve_problem(inputs)
empty!(m)
# Step 2
...
# Step 3
optimize!(m)
# ...
end
But, @odow let me know that “This is not the way” and directed me to ask here again. Hence: what is the way? (Pragmatically, I tried the above version, which seems to work fine and indeed cuts away the instantiation time entirely?)
As more context, the crux of the issue is that initializing a Model is comparatively slow, relative to setting up and solving the problem. In absolute numbers:
julia> using Chairmarks, JuMP, HiGHS
julia> opt = optimizer_with_attributes(HiGHS.Optimizer, "presolve"=>"off")
julia> @b Model($opt)
1.459 ms (2162 allocs: 52.797 KiB)
In relative terms, this is ~45% of the total time spent for my function (3.2 ms).
Full problem details & implementation
Below is the problem I actually solve: a variant of the subset sum problem, implemented as a binary integer problem.
"""
solve_subset_sum_variant(Sʲs, T, R;
verbose::Bool = false,
optimizer = optimizer_with_attributes(HiGHS.Optimizer, "presolve" => "off"))
Solves the following variant of the subset sum problem, using binary integer programming.
## Problem formulation
Given a target list `T` = (T₁, …, T_N), a set of multisets `Sʲs` = {S⁽ʲ⁾}_{j=1}^{J} with
S⁽ʲ⁾ = {s⁽ʲ⁾₁, …, s⁽ʲ⁾_{M⁽ʲ⁾}}, as well as a multiset `R` = {r₁, …, r_L}, determine if there
exist disjoint subset selections {S⁽ʲ⁾ₙ} and {Rₙ} with S⁽ʲ⁾ₙ ⊂ S⁽ʲ⁾ and Rₙ ⊂ R such that:
∑_{s ∈ S⁽ʲ⁾ₙ} s = tₙ + ∑_{r ∈ Rₙ} r
for all j = 1, …, J and n = 1, …, N and with ⋃ₙ S⁽ʲ⁾ₙ = S⁽ʲ⁾} and ⋃ₙ Rₙ = R (and
S⁽ʲ⁾_n ∩ S⁽ʲ⁾_{n′} = ∅ if n ≠ n′ and similarly for R_n ∩ R_{n′}).
## Output
A two-tuple `output` with:
- `output[1] :: Bool`: whether the problem is feasible (i.e., solvable).
- `output[2] :: Model`: the solution in binary variables `x` (into `Sʲs`) and `y` (into
`R`).
## Example
```jl
julia> Sʲs = [[1, 3, 2, 1], [3, 2, 2], [1, 1, 2, 2, 1]];
julia> T = [3, 3];
julia> R = [1];
julia> solve_subset_sum_variant(Sʲs, T, R; verbose=true);
Feasible
T = [3, 3]
R = [1]
Sʲs = [[1, 3, 2, 1], [3, 2, 2], [1, 1, 2, 2, 1]]
j = 1:
n = 1: S[1] + S[2] = T[1] + (R[1])
1 + 3 = 3 + (1)
n = 2: S[3] + S[4] = T[2]
2 + 1 = 3
j = 2:
n = 1: S[2] + S[3] = T[1] + (R[1])
2 + 2 = 3 + (1)
n = 2: S[1] = T[2]
3 = 3
j = 3:
n = 1: S[1] + S[2] + S[3] = T[1] + (R[1])
1 + 1 + 2 = 3 + (1)
n = 2: S[4] + S[5] = T[2]
2 + 1 = 3
```
"""
function solve_subset_sum_variant(
Sʲs :: AbstractVector{<:AbstractVector{<:Integer}},
T :: AbstractVector{<:Integer},
R :: AbstractVector{<:Integer};
verbose :: Bool = false,
optimizer = optimizer_with_attributes(HiGHS.Optimizer, "presolve" => "off")
)
J = length(Sʲs)
N, Mʲs, L = length(T), length.(Sʲs), length(R)
model = Model(optimizer)
verbose || MOI.set(model, MOI.Silent(), true)
@variable(model, x[1:sum(Mʲs), 1:N], Bin) # binary "indexing" into `Sʲs`
if L ≠ 0
@variable(model, y[1:L, 1:N], Bin) # binary "indexing" into `R`
end
# constraints on `x` and `y`
for n in 1:N
RHS = L ≠ 0 ? T[n] + sum(y[l,n]*R[l] for l in 1:L) : AffExpr(T[n])
for j in 1:J
mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j])
@constraint(model, sum(x[mᵍ,n]*Sʲs[j][m]
for (m,mᵍ) in enumerate(mʲ_global_idxs)) == RHS)
end
end
for j in 1:J
mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j])
@constraint(model, sum(x[mᵍ,n] for mᵍ in mʲ_global_idxs for n in 1:N) == Mʲs[j])
for mᵍ in mʲ_global_idxs
@constraint(model, sum(x[mᵍ,n] for n in 1:N) == 1)
end
end
# constraints only on `y`
if L ≠ 0
@constraint(model, sum(y[l,n] for l in 1:L for n in 1:N) == L)
for l in 1:L
@constraint(model, sum(y[l,n] for n in 1:N) == 1)
end
end
optimize!(model)
feasible = is_solved_and_feasible(model)
if !(feasible || termination_status(model) == INFEASIBLE)
error(lazy"unexpected termination status: $(termination_status(m))")
end
verbose && _print_subset_sum_result(model, Sʲs, T, R)
return feasible, model
end
function _print_subset_sum_result(io :: IO, model :: Model, Sʲs, T, R)
J = length(Sʲs)
N, Mʲs, L = length(T), length.(Sʲs), length(R)
if is_solved_and_feasible(model)
printstyled(io, "Feasible\n"; color=:green)
printstyled(io, "T = ", T, "\nR = ", R, "\nSʲs = ", Sʲs, "\n\n";
color=:light_black)
x = object_dictionary(model)[:x] # into `Sʲs`
xv = value.(x)
yv = if L ≠ 0
y = object_dictionary(model)[:y] # into `R`
value.(y)
else
Matrix{Float64}(undef, 0, 0)
end
for j in 1:J
mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j])
printstyled(io, "j = ", j, ":\n", color=:light_black)
for n in 1:N
printstyled(io, " n = ", n, ": ", color=:light_black)
s = ""
v = ""
for (m, mᵍ) in enumerate(mʲ_global_idxs)
if xv[mᵍ,n] > 0.5
s = s * (isempty(s) ? "" : " + ") * "S[$m]"
v = v * (isempty(v) ? "" : " + ") * string(Sʲs[j][m])
end
end
s = *(s, " = T[", string(n), "]")
v = *(v, " = ", string(T[n]))
if L ≠ 0
for l in 1:L
first = true
if yv[l,n] > 0.5
if first
s *= " + ("
v *= " + ("
first = false
else
s *= " + "
v *= " + "
end
s *= "R[$l]"
v *= string(R[l])
end
first || (s *= ")"; v *= ")")
end
end
println(io, s, "\n ", v)
end
end
elseif termination_status(model) == INFEASIBLE
printstyled(io, "Infeasible\n"; color=:red)
end
end
function _print_subset_sum_result(model :: Model, Sʲs, T, R)
_print_subset_sum_result(stdout, model, Sʲs, T, R)
end
