How to Develop a Sudoku Solver Using Mixed Integer Linear Programming

Posted on 8. December 20238. December 2023 by Fabian Leuthold

In today’s blog post, let’s dive into the fascinating realm of solving Sudoku puzzles using Mixed Integer Linear Programming (MILP). Check this previous post to see how a sudoku solver can be developed in python using dynamic programming and backtracking.

Sudoku, a classic logic-based number placement puzzle, has long captivated minds around the globe with its deceptively simple rules and mind-bending complexity. As we delve once more into the art of developing a Sudoku solver, this time we’ll leverage the elegance and efficiency of MILP.

The Problem

Sudoku is a number puzzle game played on a 9 x 9 grid, divided into 3×3 subgrids or “regions.” The objective is to fill the entire grid with numbers from 1 to 9, ensuring that each row, each column, and each 3×3 subgrid contains every digit exactly once. The puzzle begins with some cells pre-filled with numbers, providing initial clues for the player.

With a clear understanding of the problem at hand, we can now move on to formulating the mathematical model.

The Model

We want to find a model which solves our problem described. In the following sections I describe how the model is built. Before we start modeling it often helps to think of the answer we want to get from the model: We want to know, which cell in which row and column needs to contain which value.

Decision Variable(s)

In the pursuit of crafting an effective Sudoku solver, the choice of decision variables plays a pivotal role in formulating a mathematical model that accurately represents the puzzle-solving process. For our Mixed Integer Linear Programming (MILP) approach, the decision variables are typically binary indicators. Each variable represents a specific combination of a digit and a cell position within the Sudoku grid. By employing binary variables, we elegantly capture the essence of decision-making: whether a particular digit is assigned to a specific cell or not:

$X_{ijk}$ : Binary decision variable indicating whether in row $i$ and column $j$ the digit $k$ is placed (=1) or not (=0).

Indices and Sets

The following indices and sets provide the means to further characterize the variables and parameters of the MILP:

$i \in I$ : The $i$ -th row of the cell from the set of rows $I=\{1, 2, ..., 9\}$ .
$j \in J$ : The $j$ -th column of the cell from the set of columns $J=\{1, 2, ..., 9\}$ .
$k \in K$ : A specific digit $k$ from the set of digits $K=\{1, 2, ..., 9\}$ .
$p \in P$ : The $p$ -th subgrid-row from the set of subgrid-rows $P=\{1, 2, 3\}$ .
$q \in Q$ : The $q$ -th subgrid-column from the set of subgrid-columns $Q=\{1, 2, 3\}$ .

Parameter(s)

To find out, what parameters we need, we usually just start modeling the constraints and repeatedly add the parameters we need to write down the constraints step by step. However, here the situation is easy, we just need one parameter:

$c_{ijk}$ : Binary parameter that specifies whether in row $i$ and column $j$ the digit $k$ was provided as a clue (=1) or not (=0).

Constraints

The constraints ensuring that each row, column, and 3×3 subgrid contains every digit exactly once can be formulated as:

C0: Each cell must contain one digit
$\forall (i, j): \sum_k x_{ijk} = 1$
This constraint ensures that each cell must contain exactly one digit.
C1: Each row must contain each digit exactly once
$\forall (i, k): \sum_j x_{ijk} = 1$
This constraint ensures that each digit from 1-9 must appear exactly once in each row.
C2: Each column must contain each digit exactly once
$\forall (j, k): \sum_i x_{ijk} = 1$
This constraint ensures that each digit from 1-9 must appear exactly once in each column.
C3: Each subgrid must contain each digit exactly once
$\forall (k, p, q) \in (K \times \{1,2,3\} \times \{1,2,3\}): \sum_{i=3p-2}^{3p}\sum_{j=3q-2}^{3q} x_{ijk} = 1$
This constraint ensures that each digit from 1-9 must appear exactly once in each subgrid.

Additionally, the variables have to respect the initial clues given in the Sudoku puzzle, ensuring that the assigned digits remain unchanged:

C4: Comply with the clues provided
$\forall (i, j, k): \sum_{(i,j,k)} x_{ijk} \ge c_{(i, j, k)}$
This constraint ensures that each digit provided in the puzzle must remain unchanged.

Please note the first constraint: At the beginning, I had forgotten this constraint and therefore received incorrect solutions. Only after closer examination did I realize that the model without this constraint allows several values to be active in a cell at the same time!

Objective Function

Unlike traditional optimization problems where we aim to maximize or minimize a certain quantity, Sudoku solving is primarily about satisfying a set of constraints to reach a feasible solution that adheres to the rules of the game. The goal is to find a valid assignment of digits to cells that meets all the specified conditions. For this reason, we do not need an objective function here and therefore we just provide a pseudo objective function:

$min\;\; 0$

The Source

In the following I present the solution modeled using python together with Google OR-Tools configured to use the ‘SCIP’ solver, available with an Apache License 2.0. It is an exceptional strong and free open source solver.

The main part of the program lies in the solve() function, where the MILP is constructed and solved. The two additional functions convert_sudoku_to_binary_array() and convert_x_to_sudoku() are only used to generate a 2d-Sudoku solution from the 3d result matrix and the 3d parameter from the 2d-Sudoku specification. The source code shown beneath can also be obtained here from my github account.

from ortools.linear_solver import pywraplp
import numpy as np

def solve(sudoku):
    """
    Solve the sudoku passed
    :param sudoku: a 9x9 sudoku field
    :return: the completed 9x9 sudoku field
    """
    # Create the model.
    solver_name = 'SCIP'
    model = pywraplp.Solver.CreateSolver(solver_name)
    model.EnableOutput()

    # create clues parameter
    c = convert_sudoku_to_binary_array(sudoku)

    # create index sets
    I, J, K, P, Q = list(range(9)), list(range(9)), list(range(9)), list(range(3)), list(range(3))

    # create decision vars
    X = {}
    for i in I:
        for j in J:
            for k in K:
                X[(i, j, k)] = model.BoolVar('task_i%ij%ik%i' % (i, j, k))

    # constraint 0: Each cell must contain a digit
    for i in I:
        for j in J:
            model.Add(sum(X[(i, j, k)] for k in K) == 1)

    # constraint 1: Each row must contain each digit exactly once
    for i in I:
        for k in K:
            model.Add(sum(X[(i, j, k)] for j in J) == 1)

    # constraint 2: Each column must contain each digit exactly once
    for j in J:
        for k in K:
            model.Add(sum(X[(i, j, k)] for i in I) == 1)

    # constraint 3: Each subgrid must contain each digit exactly once
    for p in P:
        # calculate i_subset from p for sum-formula
        i_subset = list(range(3 * p, 3 * p + 3))
        for q in Q:
            # calculate j_subset from q for sum-formula
            j_subset = list(range(3 * q, 3 * q + 3))
            for k in K:
                model.Add(sum(X[(i, j, k)] for i in i_subset for j in j_subset) == 1)

    # constraint 4: Comply with the clues provided
    for i in I:
        for j in J:
            for k in K:
                model.Add(X[(i, j, k)] >= c[(i, j, k)])

    # Zielfunktion
    model.Minimize(0)

    status = model.Solve()
    if status == pywraplp.Solver.OPTIMAL:
        model_solved = True
        solution = convert_x_to_sudoku(X)
    else:
        model_solved = False
        solution = sudoku

    return model_solved, solution


def convert_sudoku_to_binary_array(sudoku):
    """
    Takes a 2x2 sudoku np-array and derives a 3-d clue param by adding a dimension
    of binary flags for each digit as required by the model

    :param sudoku: the sudoku np-array
    :return: the clue parameter for the model
    """

    binary_sudoku = np.zeros((9, 9, 9), dtype=int)
    for i in range(9):
        for j in range(9):
            if sudoku[i, j] != 0:
                binary_sudoku[i, j, sudoku[i, j] - 1] = 1

    return binary_sudoku


def convert_x_to_sudoku(X):
    """
    Convert the solver solution back to a sudoku 2x2 field

    :param X: the solution variable x[i,j,k]
    :return: a sudoku 2x2 field
    """
    sudoku = np.zeros(shape=(9, 9), dtype=np.int8)
    for i in range(0, 9):
        for j in range(0, 9):
            for k in range(0, 9):
                if X[(i, j, k)].solution_value() == 1:
                    sudoku[i, j] = k + 1
    return sudoku


if __name__ == '__main__':
    # Define a sudoku; 0 = undefined, 1-9 are clues
    sudoku = np.array([[0, 7, 2, 4, 8, 5, 0, 0, 0],
                       [4, 0, 8, 2, 0, 0, 0, 0, 0],
                       [5, 0, 0, 0, 0, 9, 4, 0, 0],
                       [0, 0, 5, 0, 0, 1, 0, 0, 8],
                       [0, 0, 0, 0, 6, 0, 0, 0, 0],
                       [1, 0, 0, 5, 0, 0, 9, 0, 0],
                       [0, 0, 4, 1, 0, 0, 0, 0, 5],
                       [0, 0, 0, 0, 0, 4, 3, 0, 7],
                       [0, 0, 0, 7, 3, 8, 2, 1, 0]])

    solved, solution = solve(sudoku)
    if solved:
        print(solution)
    else:
        print("No solution found. Unfeasible instance provided.")