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

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
optimizer = Gurobi.Optimizer
model = Model(optimizer)
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:

optimize!(model)
assert_is_solved_and_feasible(model)
Test.@test result_count(model) == 1
solution_summary(model)

solution_summary(; result = 1, verbose = false)
├ solver_name          : HiGHS
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : kHighsModelStatusOptimal
│ └ objective_bound    : 0.00000e+00
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 0.00000e+00
│ ├ dual_objective_value : NaN
│ └ relative_gap         : 0.00000e+00
└ Work counters
  ├ solve_time (sec)   : 1.96747e-01
  ├ simplex_iterations : 3389
  ├ barrier_iterations : -1
  └ node_count         : 209

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, optimizer)

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)
assert_is_solved_and_feasible(model)
solution_summary(model)

solution_summary(; result = 1, verbose = false)
├ solver_name          : HiGHS
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : kHighsModelStatusOptimal
│ └ objective_bound    : 0.00000e+00
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 0.00000e+00
│ ├ dual_objective_value : NaN
│ └ relative_gap         : 0.00000e+00
└ Work counters
  ├ solve_time (sec)   : 1.95732e-01
  ├ simplex_iterations : 3389
  ├ barrier_iterations : -1
  └ node_count         : 209

Now Gurobi has found 20 solutions:

result_count(model)

1

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}:
 5  2  1  9
 2  3  8  4
 1  8  6  7
 9  4  7  0

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:
  5 2 1 9
  2 3 8 4
+ 1 8 6 7
= 9 4 7 0

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