# Right solver for JuMP to find every solution of a linear system of equations with integer solutions

**URL:** <https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [August 10, 2020, 9:44pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709 "2020-08-10T21:44:52Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [August 10, 2020, 9:44pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/1 "2020-08-10T21:44:52Z")

</div>

I have been struggling to find an equivalent Julia function  
to `Mathematica` `FindInstance`.  
The reason behind this is to find every solution of a system  
of equations with positive integer solutions. I provide a concrete example below.

I have tried using `JuMP` (using solvers like `gurobi`) with no avail.  
For instance, I adapted one example from `JuMP` documentation  
to find try to find every solution of the following system of equations (there  
are only two solutions)

x + y = 3  
u + v = 1  
x + u = 2  
y + v = 2

> function example\_little\_sum(; verbose = true)  
> #model = Model(GLPK.Optimizer)  
> #model = Model()  
> #model = Model(with\_optimizer(Gurobi.Optimizer))  
> #model = Model(Cbc.Optimizer)  
> #m = Model(solver=IpoptSolver(print\_level=0))  
> #model = Model(solver=IpoptSolver(print\_level=0))  
> #model = Model(IpoptSolver(print\_level=0))  
> #model = IpoptSolver(print\_level=0)  
> #model = Model(IpoptSolver)  
> #model = Model(IpoptSolver)  
> model = Model(Gurobi.Optimizer)
> 
> ```
> @variable(model, 0 <= x <= 4, Int)
> @variable(model, 0 <= y <= 4, Int)
> @variable(model, 0 <= u <= 4, Int)
> @variable(model, 0 <= v <= 4, Int)
> @constraint(model, x + y == 3)
> @constraint(model, u + v == 1)
> @constraint(model, x + u == 2)
> @constraint(model, y + v == 2)
> 
> if verbose
> print(model)
> end
> 
> JuMP.optimize!(model)
> 
> x_value = JuMP.value(x)
> y_value = JuMP.value(y)
> u_value = JuMP.value(u)
> v_value = JuMP.value(v)
> num_results = result_count(model)
> if verbose
> #println("Objective value: ", obj_value)
> println("Num of results: ", num_results)
> println("x = ", x_value)
> println("y = ", y_value)
> println("u = ", u_value)
> println("v = ", v_value)
> end
> 
> ```
> 
> end

Output of the previous code is

> Academic license - for non-commercial use only  
> Feasibility  
> Subject to  
> x + y = 3.0  
> u + v = 1.0  
> x + u = 2.0  
> y + v = 2.0  
> x ≥ 0.0  
> y ≥ 0.0  
> u ≥ 0.0  
> v ≥ 0.0  
> x ≤ 4.0  
> y ≤ 4.0  
> u ≤ 4.0  
> v ≤ 4.0  
> x integer  
> y integer  
> u integer  
> v integer  
> Academic license - for non-commercial use only  
> Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (linux64)  
> Optimize a model with 4 rows, 4 columns and 8 nonzeros  
> Model fingerprint: 0xd24d95ae  
> Variable types: 0 continuous, 4 integer (0 binary)  
> Coefficient statistics:  
> Matrix range [1e+00, 1e+00]  
> Objective range [0e+00, 0e+00]  
> Bounds range [4e+00, 4e+00]  
> RHS range [1e+00, 3e+00]  
> Presolve removed 4 rows and 4 columns  
> Presolve time: 0.04s  
> Presolve: All rows and columns removed
> 
> Explored 0 nodes (0 simplex iterations) in 0.08 seconds  
> Thread count was 1 (of 4 available processors)
> 
> Solution count 1: 0
> 
> Optimal solution found (tolerance 1.00e-04)  
> Best objective 0.000000000000e+00, best bound 0.000000000000e+00, gap 0.0000%  
> Num of results: 1  
> x = 2.0  
> y = 1.0  
> u = -0.0  
> v = 1.0

As you can see, JuMP —via gurobi— only finds one solution.  
On the other hand, Mathematica (correctly) find two solutions

x = 1, y = 2, u = 1, v = 0  
x = 2, y = 1, u = 0, v = 1

So I have two questions: is there any solver available (free or free for  
Academics if possible) that can return every solution of the system mentioned  
above? Another option I can  
think of is using Interval Arithmetic. Unfortunately, I have been unable to  
have a working version of this same problem.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 10, 2020, 10:00pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/2 "2020-08-10T22:00:45Z")

</div>

You can use `IntervalConstraintProgramming.jl` to eliminate pieces of the domain:

```julia

julia> using IntervalConstraintProgramming, ModelingToolkit

julia> vars = @variables x, y, u, v
(x, y, u, v)

julia> C1 = Separator(vars, x + y == 3);

julia> C2 = Separator(vars, u + v == 1);

julia> C3 = Separator(vars, x + u == 2);

julia> C4 = Separator(vars, y + v == 2)

julia> X = IntervalBox(0..4, 4) # domain
[0, 4] × [0, 4] × [0, 4] × [0, 4]

julia> X = C1(X)[1]
[0, 3] × [0, 3] × [0, 4] × [0, 4]

julia> X = C2(X)[1]
[0, 3] × [0, 3] × [0, 1] × [0, 1]

julia> X = C3(X)[1]
[1, 2] × [0, 3] × [0, 1] × [0, 1]

julia> X = C4(X)[1]
[1, 2] × [1, 2] × [0, 1] × [0, 1]

```

This at least vastly reduces the number of integers you need to try.  
Currently it is not possible to use constraints that restrict to being close to an integer, though that  
would definitely be nice.

In general the JuMP solvers return only one solution.  
A “cheating” method I have seen used at some point to find more solutions is to add an explicit constraints of not being equal to any previously-found solution.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 10, 2020, 10:07pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/3 "2020-08-10T22:07:53Z")

</div>

In fact you can almost solve the problem already:

```julia
julia> p = pave(C1 ∩ C2 ∩ C3 ∩ C4, X, 0.01)
Paving:
- tolerance ϵ = 0.01
- inner approx. of length 0
- boundary approx. of length 374

julia> p.boundary
374-element Array{IntervalBox{4,Float64},1}:
 [1.99175, 2] × [1, 1.00825] × [0, 0.0082499] × [0.99175, 1]
 [1.99175, 2] × [1.00824, 1.00825] × [0.00824989, 0.0082499] × [0.99175, 0.991751]
 [1.98362, 1.99176] × [1.00824, 1.00825] × [0.00824989, 0.0082499] × [0.99175, 0.991751]
 [1.98362, 1.99176] × [1.00824, 1.01638] × [0.00824989, 0.0163719] × [0.983628, 0.991751]
 [1.98362, 1.98363] × [1.01637, 1.01638] × [0.0163718, 0.0163719] × [0.983628, 0.983629]
...
[snipped]

julia> unique([round.(X) for X in p.boundary])
4-element Array{IntervalBox{4,Float64},1}:
 [2, 2] × [1, 1] × [0, 0] × [1, 1]
 [1, 2] × [1, 1] × [0, 0] × [1, 1]
 [1, 2] × [1, 2] × [0, 1] × [0, 1]
 [1, 1] × [2, 2] × [1, 1] × [0, 0]

```

I’m not sure why we get those intervals still.

---

<div class="post-metadata">

**Author:** ![mtanneau](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mtanneau/32/17787_2.png) [@mtanneau](https://discourse.julialang.org/u/mtanneau)\
**Post date:** [August 10, 2020, 11:14pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/4 "2020-08-10T23:14:22Z")

</div>

> find every solution of a system of equations with positive integer solutions

This is not what MIP solvers have been designed to do, so you might have to hack into things a bit.  
You might have better luck with Constraint Programming solvers, though they are not as well interfaced in Julia.

Not 100% sure, but [SCIP](https://www.scipopt.org/) might have this functionality.  
I don’t think it will be exposed at the JuMP level.

Otherwise, if all your variables are binary, you can use so-called “no-good cuts” to eliminate integer-feasible solutions as they are found, thereby enumerating all feasible solutions.

> A “cheating” method I have seen used at some point to find more solutions is to add an explicit constraints of not being equal to any previously-found solution.

==\> That’s the one 🙂

Assuming that all your variables are binaries, and given \bar{x} \in \{0, 1\}^{n}, the inequaliy

\sum\_{j | \bar{x}\_{j} = 1} (1 - x\_{j}) + \sum\_{j | \bar{x}\_{j} = 0} x\_{j} \geq 1

cuts off exactly \bar{x} from the binary feasible set.

You add those as lazy constraints within a callback, and you will end up enumerating all feasible solutions.  
Just be careful about computing times.

If you have general integer variables, the above method does not work, because removing a single integer point would yield a non-convex set. You can always use binary representation of bounded integers though.

---

<div class="post-metadata">

**Author:** ![mtanneau](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mtanneau/32/17787_2.png) [@mtanneau](https://discourse.julialang.org/u/mtanneau)\
**Post date:** [August 10, 2020, 11:17pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/5 "2020-08-10T23:17:52Z")

</div>

Actually, Gurobi has a [parameter for finding multiple solutions](https://www.gurobi.com/documentation/9.0/refman/finding_multiple_solutions.html).

Never tried it myself.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 10, 2020, 11:28pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/6 "2020-08-10T23:28:45Z")

</div>

> [@dpsanders](#):
>
> ```julia
> 4
> 
> ```

@mtanneau: I gave a pure-Julia constraint-programming solution one post above 😉

To finish it off, define

```julia
function integerize(X, i)
    Xi = X[i]
    a = ceil(Xi.lo)
    b = floor(Xi.hi)
    
    if a > b
        return emptyinterval(X)
    end

    return setindex(X, interval(a, b), i)
end 

function integerize(X)
    for i in 1:length(X)
        X = integerize(X, i)
    end
    
    return X
end

```

Then

```julia
julia> p = pave(C1 ∩ C2 ∩ C3 ∩ C4, X, 1e-2);

julia> unique(integerize.(p.boundary))
3-element Array{IntervalBox{4,Float64},1}:
 [2, 2] × [1, 1] × [0, 0] × [1, 1]
 ∅ × ∅ × ∅ × ∅
 [1, 1] × [2, 2] × [1, 1] × [0, 0]

```

give the desired solutions!

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [August 10, 2020, 11:39pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/7 "2020-08-10T23:39:18Z")

</div>

Thanks for your responses @dpsanders @mtanneau. I’ll see how does @dpsanders’s solution works for my original purpose and provide feedback. Which in case you’re curious, it refers to Lemma 5.4 with Corollary 5.1 in [here.](http://www.ma.rhul.ac.uk/~uvah099/Maths/doubleRevised2.pdf)

---

<div class="post-metadata">

**Author:** ![leethargo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leethargo/32/6004_2.png) [@leethargo](https://discourse.julialang.org/u/leethargo)\
**Post date:** [August 11, 2020, 8:00am UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/8 "2020-08-11T08:00:24Z")

</div>

> [@mtanneau](#):
>
> Not 100% sure, but [SCIP](https://www.scipopt.org/) might have this functionality.  
> I don’t think it will be exposed at the JuMP level.

That’s exactly right. SCIP has a feature to [count solutions](https://www.scipopt.org/doc/html/COUNTER.php) (and can also store the solution values).  
But this function is not exposed in SCIP.jl.

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [August 11, 2020, 4:18pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/9 "2020-08-11T16:18:14Z")

</div>

I’m not positive now.

I tried @dpsanders’s solution on a larger —but still on topic—  
problem and found a huge use of resources. Example below.

```julia
μ = [0,2,3,1]
ν = [2,2,1,1]

components = 4 # μ |> length
vars = @variables y[1:components,1:components]

constriction_list = Separator[]
for i in 1:components
    push!(constriction_list, Separator(vars..., sum(vars[1][i,1:components]) == μ[i]))
end
for i in 1:components
    push!(constriction_list, Separator(vars..., sum(vars[1][1:components,i]) == ν[i]))
end

X = IntervalBox(0..(components*components), (components*components)) # domain
for c1 in constriction_list
    X = c1(X)[1]
end

ini = ∩(constriction_list[1], constriction_list[2])
for ci in constriction_list[3:end]
    ini = ∩(ini, ci)
end
p = pave(ini, X, 1e-2)

unique(integerize.(p.boundary))

```

So, is this (large memory usage) an issue with the way I extended the solution or indeed `IntervalConstraintProgramming`  
is not apt to this new problem?

ps1. I was unable to use a function for the previous code. I had a problem similar to `The applicable method may be too new: running in world age 28047, while current world is 28055.`  
ps2. In case someone wants to copy-paste code from this post, you need to add to @dpsanders’s code `using IntervalArithmetic`

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 11, 2020, 4:29pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/10 "2020-08-11T16:29:30Z")

</div>

So that’s in 16 dimensions? That is indeed rather high. I can’t look at this right now but will try later

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 11, 2020, 6:42pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/11 "2020-08-11T18:42:40Z")

</div>

Here’s a gist with a better approach:

[https://gist.github.com/dpsanders/5cf6448b3ef3d280b6d70a8fefe34963](https://gist.github.com/dpsanders/5cf6448b3ef3d280b6d70a8fefe34963)

Basically `integerize` needs to be mixed in throughout the process, and we use Contractors instead of Separators.

(Separators are used when you want to find inner and outer approximations. We want only an outer approximation here.)

This new version produces the following results in a few seconds:

```julia
24-element Array{StaticArrays.SArray{Tuple{16},Int64,1,16},1}:
 [0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1]
 [0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 1, 0, 0, 1, 0]
 [0, 2, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0]
 [0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1]
 [0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0]
 ⋮
 [0, 0, 1, 1, 0, 2, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0]
 [0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0]
 [0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0]
 [0, 0, 2, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0]
 [0, 0, 1, 1, 0, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 0]

```

Do you know what the correct solutions are in this case?

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 11, 2020, 6:45pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/12 "2020-08-11T18:45:29Z")

</div>

By the way, after applying each contractor once to the initial box `[0..16]^16`, it already reduces to

```julia
[0, 0] × [0, 2] × [0, 2] × [0, 1] × [0, 0] × [0, 2] × [0, 2] × [0, 1] × [0, 0] × [0, 1] × [0, 1] × [0, 1] × [0, 0] × [0, 1] × [0, 1] × [0, 1]

```

This really shows the power of constraint propagation!

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 11, 2020, 6:52pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/13 "2020-08-11T18:52:51Z")

</div>

I forgot to say that the world-age error comes from code generation inside `Separator` / `Contractor`.

If you want to get around this while still being inside a function, you need something like

`Base.invokelatest(pave, C, X)`

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [August 11, 2020, 7:24pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/14 "2020-08-11T19:24:26Z")

</div>

There are the 24 solutions (returned by Mathematica):

```julia
0,0,0,0,0,0,1,1,1,2,0,0,1,0,0,0
0,0,0,0,0,0,1,1,2,1,0,0,0,1,0,0
0,0,0,0,0,1,0,1,1,1,1,0,1,0,0,0
0,0,0,0,0,1,0,1,2,0,1,0,0,1,0,0
0,0,0,0,0,1,0,1,2,1,0,0,0,0,1,0
0,0,0,0,0,1,1,0,1,1,0,1,1,0,0,0
0,0,0,0,0,1,1,0,2,0,0,1,0,1,0,0
0,0,0,0,0,1,1,0,2,1,0,0,0,0,0,1
0,0,0,0,0,2,0,0,1,0,1,1,1,0,0,0
0,0,0,0,0,2,0,0,2,0,0,1,0,0,1,0
0,0,0,0,0,2,0,0,2,0,1,0,0,0,0,1
0,0,0,0,1,0,0,1,0,2,1,0,1,0,0,0
0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0
0,0,0,0,1,0,0,1,1,2,0,0,0,0,1,0
0,0,0,0,1,0,1,0,0,2,0,1,1,0,0,0
0,0,0,0,1,0,1,0,1,1,0,1,0,1,0,0
0,0,0,0,1,0,1,0,1,2,0,0,0,0,0,1
0,0,0,0,1,1,0,0,0,1,1,1,1,0,0,0
0,0,0,0,1,1,0,0,1,0,1,1,0,1,0,0
0,0,0,0,1,1,0,0,1,1,0,1,0,0,1,0
0,0,0,0,1,1,0,0,1,1,1,0,0,0,0,1
0,0,0,0,2,0,0,0,0,1,1,1,0,1,0,0
0,0,0,0,2,0,0,0,0,2,0,1,0,0,1,0
0,0,0,0,2,0,0,0,0,2,1,0,0,0,0,1

```

I tried checking both solutions (yours and Mathematica) using `sort` with no avail. I’ll send the results soon, but I’m confident both are equal.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 11, 2020, 7:30pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/15 "2020-08-11T19:30:25Z")

</div>

It looks like the variables are stored in a different order, which is making things difficult to compare?

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [August 11, 2020, 7:55pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/16 "2020-08-11T19:55:58Z")

</div>

Both solutions are equal.  
I exported your solution into an csv and compared it against `Mathematica` solution. Unfortunately, I have been unable to compare it using Julia.  
I tried the following  
`sortslices(mma_output,dims=1)-sortslices(hcat(collect(output)...)',dims=1)`, where `mma_output` was obtained by `mma_output = readdlm("/home/david/example.csv",',',Int)`.  
I was expecting to find a zero matrix but, as you can check, it is not the case. example.csv is a file containing the list in [post 14](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/14).

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 11, 2020, 8:33pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/17 "2020-08-11T20:33:05Z")

</div>

What order is Mathematica using for the variables?  
I tried just transposing the matrix of variables in the Julia code, but that didn’t help.

I would just apply a permutation matrix that maps the order to the correct one.  
I’m not sure that just sorting is enough.

---

<div class="post-metadata">

**Author:** ![damaro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/damaro/32/17002_2.png) [@damaro](https://discourse.julialang.org/u/damaro)\
**Post date:** [August 11, 2020, 9:19pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/18 "2020-08-11T21:19:38Z")

</div>

The orden I gave to Mathematica is

```julia
# non Julia code
{x[1, 1], x[2, 1], x[3, 1], x[4, 1], x[1, 2], x[2, 2], x[3, 2], 
 x[4, 2], x[1, 3], x[2, 3], x[3, 3], x[4, 3], x[1, 4], x[2, 4], 
 x[3, 4], x[4, 4]}

```

Inspecting `vars`, it looks like this is the same order as in Julia. Still, I am unable to match both solutions in Julia.

Edit. Silly mistake. I interchanged `\nu` and `\mu`. Now running `sortslices(mma_output,dims=1)-sortslices(hcat(collect(output)...)',dims=1)` results in the zero matrix.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 11, 2020, 9:36pm UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/19 "2020-08-11T21:36:38Z")

</div>

Here’s a way to force Julia to use exactly that order of the variables. Then the output looks the same:

```julia
vars = @variables y[1:n, 1:n]
vars = vars[1]
vars = permutedims(vars, (1, 2))

contractors = []

for i in 1:n
    push!(contractors, Contractor(vec(vars), sum(vars[i, 1:n]) - μ[i]))
end

for i in 1:n
    push!(contractors, Contractor(vec(vars), sum(vars[1:n, i]) - ν[i]))
end

```

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [May 19, 2021, 2:00am UTC](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709/20 "2021-05-19T02:00:46Z")

</div>

Hi @damaro, here is how you do this in Gurobi (I’m using Gurobi 0.9.12, JuMP 0.21.7)

```julia
## Packages
using JuMP
using Gurobi

## This function will return an array of named tuples in the variables with objective value.
function example_little_sum(; verbose::Bool = true)
    ## Model
    model = Model(Gurobi.Optimizer)

    @variable(model, 0 <= x <= 4, Int)
    @variable(model, 0 <= y <= 4, Int)
    @variable(model, 0 <= u <= 4, Int)
    @variable(model, 0 <= v <= 4, Int)
    @constraint(model, x + y == 3)
    @constraint(model, u + v == 1)
    @constraint(model, x + u == 2)
    @constraint(model, y + v == 2)
    if verbose; print(model); end

    ## Gurobi parameters - see https://www.gurobi.com/documentation/9.0/refman/finding_multiple_solutions.html
    JuMP.set_optimizer_attribute(model, "PoolSearchMode", 2) 
    JuMP.set_optimizer_attribute(model, "PoolSolutions", 100) # 100 is an arbitrary (large enough) whole number

    if !verbose; JuMP.set_optimizer_attribute(model, "OutputFlag", 0); end

    ## Optimize:
    JuMP.optimize!(model)

    ## Results
    num_results = result_count(model)
    if verbose; println("Num of results: ", num_results); end

    Results = Array{NamedTuple{(:x,:y,:u,:v,:obj),NTuple{5,Int64}},1}() ## Note the conversion to Int64
    for i in 1:num_results
        sol = (x=JuMP.value(x; result=i), y=JuMP.value(y; result=i),
               u=JuMP.value(u; result=i), v=JuMP.value(v; result=i),
               obj=JuMP.objective_value(model; result = i))
        push!(Results,sol)
        if verbose
            println(Results[i])
        end
    end

    return Results
end

```

[Next page](https://discourse.julialang.org/t/right-solver-for-jump-to-find-every-solution-of-a-linear-system-of-equations-with-integer-solutions/44709.md?page=2)