Finding multiple feasible solutions

This tutorial was generated using Literate.jl. Download the source as a .jl file.

Author: James Foster (@jd-foster)

This tutorial demonstrates how to formulate and solve a combinatorial problem with multiple feasible solutions. In fact, we will see how to find all feasible solutions to our problem. We will also see how to enforce an "all-different" constraint on a set of integer variables.

Required packages

This tutorial uses the following packages:

using JuMP
import Gurobi
import Test

Warning

This tutorial uses Gurobi.jl as the solver because it supports returning multiple feasible solutions, something that open-source MIP solvers such as HiGHS do not currently support. Gurobi is a commercial solver and requires a paid license. However, there are free licenses available for academic and student users. See Gurobi.jl for more details.

Symmetric number squares

Symmetric number squares and their sums often arise in recreational mathematics. Here are a few examples:

   1 5 2 9       2 3 1 8        5 2 1 9
   5 8 3 7       3 7 9 0        2 3 8 4
+  2 3 4 0     + 1 9 5 6      + 1 8 6 7
=  9 7 0 6     = 8 0 6 4      = 9 4 7 0

Notice how all the digits 0 to 9 are used at least once, the first three rows sum to the last row, the columns in each are the same as the corresponding rows (forming a symmetric matrix), and 0 does not appear in the first column.

We will answer the question: how many such squares are there?

JuMP model

We now encode the symmetric number square as a JuMP model. First, we need a symmetric matrix of decision variables between 0 and 9 to represent each number:

n = 4
model = Model()
set_silent(model)
@variable(model, 0 <= x_digits[row in 1:n, col in 1:n] <= 9, Int, Symmetric)

4×4 LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}:
 x_digits[1,1]  x_digits[1,2]  x_digits[1,3]  x_digits[1,4]
 x_digits[1,2]  x_digits[2,2]  x_digits[2,3]  x_digits[2,4]
 x_digits[1,3]  x_digits[2,3]  x_digits[3,3]  x_digits[3,4]
 x_digits[1,4]  x_digits[2,4]  x_digits[3,4]  x_digits[4,4]

We modify the lower bound to ensure that the first column cannot contain 0:

set_lower_bound.(x_digits[:, 1], 1)

4-element Vector{Nothing}:
 nothing
 nothing
 nothing
 nothing

Then, we need a constraint that the sum of the first three rows equals the last row:

@expression(model, x_base_10, x_digits * [1_000, 100, 10, 1]);
@constraint(model, sum(x_base_10[i] for i in 1:n-1) == x_base_10[n])

\[ 1000 x\_digits_{1,1} + 1100 x\_digits_{1,2} + 100 x\_digits_{2,2} + 1010 x\_digits_{1,3} + 110 x\_digits_{2,3} + 10 x\_digits_{3,3} - 999 x\_digits_{1,4} - 99 x\_digits_{2,4} - 9 x\_digits_{3,4} - x\_digits_{4,4} = 0 \]

And we use MOI.AllDifferent to ensure that each digit is used exactly once in the upper triangle matrix of x_digits:

x_digits_upper = [x_digits[i, j] for j in 1:n for i in 1:j]
@constraint(model, x_digits_upper in MOI.AllDifferent(length(x_digits_upper)));

If we optimize this model, we find that Gurobi has returned one solution:

set_optimizer(model, Gurobi.Optimizer)
optimize!(model)
Test.@test is_solved_and_feasible(model)
Test.@test result_count(model) == 1
solution_summary(model)

* Solver : Gurobi

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "Model was solved to optimality (subject to tolerances), and an optimal solution is available."

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : 0.00000e+00
  Objective bound    : 0.00000e+00
  Relative gap       : 0.00000e+00
  Dual objective value : 0.00000e+00

* Work counters
  Solve time (sec)   : 4.50032e-02
  Simplex iterations : 1587
  Barrier iterations : 0
  Node count         : 255

To return multiple solutions, we need to set Gurobi-specific parameters to enable the solution pool. Moreover, there is a bug in Gurobi that means the solution pool is not activated if we have already solved the model once. To work around the bug, we need to reset the optimizer. If you turn the solution pool options on before the first solve you do not need to reset the optimizer.

set_optimizer(model, Gurobi.Optimizer)

Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 722777
Set parameter GURO_PAR_SPECIAL
WLS license 722777 - registered to JuMP Development

The first option turns on the exhaustive search mode for multiple solutions:

set_attribute(model, "PoolSearchMode", 2)

The second option sets a limit for the number of solutions found:

set_attribute(model, "PoolSolutions", 100)

Here the value 100 is an "arbitrary but large enough" whole number for our particular model (and in general will depend on the application).

We can then call optimize! and view the results.

optimize!(model)
Test.@test is_solved_and_feasible(model)
solution_summary(model)

* Solver : Gurobi

* Status
  Result count       : 20
  Termination status : OPTIMAL
  Message from the solver:
  "Model was solved to optimality (subject to tolerances), and an optimal solution is available."

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : 0.00000e+00
  Objective bound    : 0.00000e+00
  Relative gap       : 0.00000e+00
  Dual objective value : 0.00000e+00

* Work counters
  Solve time (sec)   : 3.73312e-01
  Simplex iterations : 19526
  Barrier iterations : 0
  Node count         : 4661

Now Gurobi has found 20 solutions:

Test.@test result_count(model) == 20

Test Passed

Viewing the Results

Access the various feasible solutions by using the value function with the result keyword:

solutions =
    [round.(Int, value.(x_digits; result = i)) for i in 1:result_count(model)];

Here we have converted the solution to an integer after rounding off very small numerical tolerances.

An example of one feasible solution is:

solutions[1]

4×4 Matrix{Int64}:
 1  5  2  9
 5  8  3  7
 2  3  4  0
 9  7  0  6

and we can nicely print out all the feasible solutions with

function solution_string(x::Matrix)
    header = [" ", " ", "+", "="]
    return join([join(vcat(header[i], x[i, :]), " ") for i in 1:4], "\n")
end

for i in 1:result_count(model)
    println("Solution $i: \n", solution_string(solutions[i]), "\n")
end

Solution 1:
  1 5 2 9
  5 8 3 7
+ 2 3 4 0
= 9 7 0 6

Solution 2:
  1 5 2 9
  5 7 4 6
+ 2 4 0 8
= 9 6 8 3

Solution 3:
  1 3 2 7
  3 6 5 4
+ 2 5 0 8
= 7 4 8 9

Solution 4:
  2 3 1 7
  3 5 6 4
+ 1 6 0 8
= 7 4 8 9

Solution 5:
  2 1 4 8
  1 9 6 7
+ 4 6 3 5
= 8 7 5 0

Solution 6:
  3 2 1 7
  2 9 4 5
+ 1 4 0 6
= 7 5 6 8

Solution 7:
  1 2 3 7
  2 9 6 8
+ 3 6 4 5
= 7 8 5 0

Solution 8:
  2 3 1 8
  3 7 9 0
+ 1 9 5 6
= 8 0 6 4

Solution 9:
  5 1 3 9
  1 0 4 6
+ 3 4 8 7
= 9 6 7 2

Solution 10:
  5 2 1 9
  2 6 8 7
+ 1 8 3 4
= 9 7 4 0

Solution 11:
  5 2 1 9
  2 3 8 4
+ 1 8 6 7
= 9 4 7 0

Solution 12:
  2 1 6 9
  1 3 0 5
+ 6 0 7 4
= 9 5 4 8

Solution 13:
  1 2 5 9
  2 6 4 3
+ 5 4 7 8
= 9 3 8 0

Solution 14:
  1 2 5 9
  2 4 3 0
+ 5 3 8 7
= 9 0 7 6

Solution 15:
  2 1 4 8
  1 5 6 3
+ 4 6 7 9
= 8 3 9 0

Solution 16:
  5 2 1 9
  2 7 4 3
+ 1 4 0 6
= 9 3 6 8

Solution 17:
  5 1 3 9
  1 4 0 6
+ 3 0 8 2
= 9 6 2 7

Solution 18:
  1 2 3 7
  2 5 6 4
+ 3 6 8 9
= 7 4 9 0

Solution 19:
  1 4 2 8
  4 7 5 6
+ 2 5 0 9
= 8 6 9 3

Solution 20:
  3 2 1 6
  2 0 4 7
+ 1 4 9 5
= 6 7 5 8

The result is the full list of feasible solutions. So the answer to "how many such squares are there?" turns out to be 20.