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 removedExplored 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.