# Reusing a JuMP \`Model\` for several, different optimization problems

**URL:** <https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031>\
**Category:** Optimization (Mathematical)\
**Created:** [October 31, 2024, 8:12am UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031 "2024-10-31T08:12:48Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [October 31, 2024, 8:12am UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/1 "2024-10-31T08:12:48Z")

</div>

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:

```julia
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:

```julia
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
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](https://en.wikipedia.org/wiki/Subset_sum_problem), implemented as a binary integer problem.
> 
> ````julia
> """
> 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
> 
> ````

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 31, 2024, 8:55pm UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/2 "2024-10-31T20:55:10Z")

</div>

Can you provide example input data?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 31, 2024, 8:59pm UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/3 "2024-10-31T20:59:46Z")

</div>

> [@tchr](#):
>
> `model = Model(optimizer)`

A quick win is probably changing this to:

```julia
    inner = MOI.instantiate(optimizer)
    model = direct_model(inner)
    set_string_names_on_creation(model, false)

```

> This is not the way

I said this because `Model(optimizer)` should not be the bottleneck. Trying to be clever by using `empty!` just feels like the wrong approach.

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [October 31, 2024, 9:03pm UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/4 "2024-10-31T21:03:05Z")

</div>

Sure (it’s in the doc string example):

```julia
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=false);

```

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [October 31, 2024, 9:06pm UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/5 "2024-10-31T21:06:58Z")

</div>

> [@odow](#):
>
> I said this because `Model(optimizer)` should not be the bottleneck. Trying to be clever by using `empty!` just feels like the wrong approach.

I guess the main thing I’m not understanding is why `Model(opt)` is taking on the order of ms?

Your suggestion is indeed much faster (~43 μs). Why is it so much faster?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 31, 2024, 9:14pm UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/6 "2024-10-31T21:14:49Z")

</div>

> Why is it so much faster?

Because it does so many more things. It is equivalent to:

```julia
using JuMP, HiGHS
inner = MOI.Utilities.CachingOptimizer(
    MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    MOI.Bridges.full_bridge_optimizer(
        HiGHS.Optimizer(),
        Float64,
    ),
)
MOI.set(inner, MOI.RawOptimizerAttribute("presolve"), "off")
model = direct_model(inner)

```

You end up with two copies of the problem, one in the cache and one in HiGHS, and there is also a bridging layer (that ends up not being used).

All these extras are on by default, because we assume that most people do not care about 3ms of overhead and want the convenience of swapping solvers, bridge reformulations, etc. 😄

But you can get very close to the solver with minimal overhead by using `direct_model`, and there’s no need for `empty!`.

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [October 31, 2024, 9:24pm UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/7 "2024-10-31T21:24:45Z")

</div>

Thanks a bunch for the explanations. Is there a rule of thumb for when I might need the `CachingOptimizer` and bridges etc.? In this case, why don’t I need it?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 31, 2024, 9:49pm UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/8 "2024-10-31T21:49:59Z")

</div>

You need a cache if you:

- want to change the solver
- want to implement a modification that the solver doesn’t support (it’s a bit moot with HiGHS, because I think I’ve wrapped everything. But some solvers like SCS doesn’t support adding constraints after solve, for example)
- want an efficient `copy_to` method to begin with, at the expense of having an additional copy of the model

You need a bridge if you:

- Use constraint types that are not natively supported by the solver. For HiGHS, this is things like indicators, SOS1 and SOS2, and constraint programming sets like AllDifferent.

If you use _only_ HiGHS, and you are solving simple LP or MIPs, then direct mode is a good choice. But it isn’t the default because most people don’t need or care about it.

See these parts of the documentation:

- [Models · JuMP](https://jump.dev/JuMP.jl/stable/manual/models/#Backends)
- [Models · JuMP](https://jump.dev/JuMP.jl/stable/manual/models/#Direct-mode)
- [Performance tips · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/getting_started/performance_tips/#Disable-string-names)

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [November 1, 2024, 10:25am UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/9 "2024-11-01T10:25:09Z")

</div>

Thanks a lot for taking the time to explain all this in detail: I understand much better now (and the defaults make perfect sense, yeah).  
Thanks also for JuMP.jl - it’s such a powerful enabler, and allows me to explore things that otherwise would’ve been way out of scope.

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [November 1, 2024, 10:29am UTC](https://discourse.julialang.org/t/reusing-a-jump-model-for-several-different-optimization-problems/122031/10 "2024-11-01T10:29:53Z")

</div>

> [@odow](#):
>
> See these parts of the documentation:
> 
> - [Models · JuMP](https://jump.dev/JuMP.jl/stable/manual/models/#Backends)
> - [Models · JuMP](https://jump.dev/JuMP.jl/stable/manual/models/#Direct-mode)
> - [Performance tips · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/getting_started/performance_tips/#Disable-string-names)

And the docs explain it really well too; thanks for those pointers also.