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

damaro · August 10, 2020, 9:44pm

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.

dpsanders · August 10, 2020, 10:00pm

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


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.

dpsanders · August 10, 2020, 10:07pm

In fact you can almost solve the problem already:

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.

mtanneau · August 10, 2020, 11:14pm

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

mtanneau · August 10, 2020, 11:17pm

Actually, Gurobi has a parameter for finding multiple solutions.

Never tried it myself.

dpsanders · August 10, 2020, 11:28pm

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

To finish it off, define

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> 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!

damaro · August 10, 2020, 11:39pm

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.

leethargo · August 11, 2020, 8:00am

That’s exactly right. SCIP has a feature to count solutions (and can also store the solution values).
But this function is not exposed in SCIP.jl.

damaro · August 11, 2020, 4:18pm

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.

μ = [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

dpsanders · August 11, 2020, 4:29pm

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

dpsanders · August 11, 2020, 6:42pm

Here’s a gist with a better approach:

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:

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?

dpsanders · August 11, 2020, 6:45pm

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

[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!

dpsanders · August 11, 2020, 6:52pm

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)

damaro · August 11, 2020, 7:24pm

There are the 24 solutions (returned by Mathematica):

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.

dpsanders · August 11, 2020, 7:30pm

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

damaro · August 11, 2020, 7:55pm

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.

dpsanders · August 11, 2020, 8:33pm

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.

damaro · August 11, 2020, 9:19pm

The orden I gave to Mathematica is

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

dpsanders · August 11, 2020, 9:36pm

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

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

jd-foster · May 19, 2021, 2:00am

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

## 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