Integer Programming for Sudoku

In this section we will use a Sudoku game to illustrate how to use integer and multi-dimensional arrays in RSOME. Sudoku is a popular number puzzle. The goal is to place the digits in [1,9] on a nine-by-nine grid, with some of the digits already filled in. The solution must satisfy the following four rules:

1. Each cell contains an integer in [1,9].
2. Each row must contain each of the integers in [1,9].
3. Each column must contain each of the integers in [1,9].
4. Each of the nine 3x3 squares with bold outlines must contain each of the integers in [1,9].

The Sudoku game can be considered as a feasibility optimization problem with the objective to be zero and constraints used to fulfill above rules. Consider a binary variable \(x_{ijk}\in \{0, 1\}\), with \(i \in [0, 8]\), \(j \in [0, 8]\), and \(k \in [0, 8]\). It equals to one if an integer \(k+1\) is placed in a cell at the \(i\)th row and \(j\)th column. Let \(a_{ij}\) be the known number at the \(i\)th row and \(j\)th column, with \((i, j)\in\mathcal{I}\times\mathcal{J}\) as \(\mathcal{I}\) and \(\mathcal{J}\) are the row and column indices of numbers with known values, the Sudoku game can be thus written as the following integer program:

In the following code, we are using RSOME to implement such a model.

import rsome as rso
import numpy as np
from rsome import ro
from rsome import grb_solver as grb

# Sudoku puzzle
# zeros represent unknown numbers
puzzle = np.array([[5, 3, 0, 0, 7, 0, 0, 0, 2],
                   [6, 0, 0, 1, 9, 5, 0, 0, 0],
                   [0, 9, 8, 0, 0, 0, 0, 6, 0],
                   [8, 0, 0, 0, 6, 0, 0, 0, 3],
                   [4, 0, 0, 8, 0, 3, 0, 0, 1],
                   [7, 0, 0, 0, 2, 0, 0, 0, 6],
                   [0, 6, 0, 0, 0, 0, 2, 8, 0],
                   [0, 0, 0, 4, 1, 9, 0, 0, 5],
                   [0, 0, 0, 0, 8, 0, 0, 7, 9]])

# create model and binary decision variables
model = ro.Model()
x = model.dvar((9, 9, 9), vtype='B')

# objective is set to be zero
model.min(0)

# constraints 1 to 3
model.st(x.sum(axis=0) == 1,
         x.sum(axis=1) == 1,
         x.sum(axis=2) == 1)

# constraints 4
i, j = np.where(puzzle)
model.st(x[i, j, puzzle[i, j]-1] == 1)

# constraints 5
for i in range(0, 9, 3):
    for j in range(0, 9, 3):
        model.st(x[i: i+3, j: j+3, :].sum(axis=(0, 1)) == 1)

# solve the integer programming problem
model.solve(grb)

Being solved by Gurobi...
Solution status: 2
Running time: 0.0017s

The binary variable \(x_{ijk}\) is defined to be a three-dimensional array x with the shape to be (9, 9, 9). Based on the decision variable x, each set of constraints can be formulated as the array form by using the sum() method. The method sum() in RSOME is consistent with that in NumPy, where you may use the axis argument to specify along which axis the sum is performed.

The Sudoku puzzle and the its solution are presented below.

print(puzzle)                   # display the Sudoku puzzle

[[5 3 0 0 7 0 0 0 2]
 [6 0 0 1 9 5 0 0 0]
 [0 9 8 0 0 0 0 6 0]
 [8 0 0 0 6 0 0 0 3]
 [4 0 0 8 0 3 0 0 1]
 [7 0 0 0 2 0 0 0 6]
 [0 6 0 0 0 0 2 8 0]
 [0 0 0 4 1 9 0 0 5]
 [0 0 0 0 8 0 0 7 9]]

x_sol = x.get().astype('int')   # retrieve the solution as an array of integers

print((x_sol * np.arange(1, 10).reshape((1, 1, 9))).sum(axis=2))

[[5 3 4 6 7 8 9 1 2]
 [6 7 2 1 9 5 3 4 8]
 [1 9 8 3 4 2 5 6 7]
 [8 5 9 7 6 1 4 2 3]
 [4 2 6 8 5 3 7 9 1]
 [7 1 3 9 2 4 8 5 6]
 [9 6 1 5 3 7 2 8 4]
 [2 8 7 4 1 9 6 3 5]
 [3 4 5 2 8 6 1 7 9]]

Note that in defining “constraints 4”, variables i and j represent the row and column indices of the fixed elements, which can be retrieved by the np.where() function. An alternative approach is to use the boolean indexing of arrays, as the code below.

# an alternative approach for constraints 4
is_fixed = puzzle > 0
model.st(x[is_fixed, puzzle[is_fixed]-1] == 1)

The variable is_fixed is an array with elements to be True if the numbers are fixed and False if the numbers are unknown. Such an array with boolean data can also be used as indices, thus defining the same constraints.