Solving a Sudoku puzzle with Integer Programming

From Wikipedia:

Sudoku, originally called Number Place, is a logic-based, combinatorial number-placement puzzle. The objective is to fill a 9×9 grid with digits so that each column, each row, and each of the nine 3×3 sub-grids that compose the grid (also called "boxes", "blocks", "regions", or "sub-squares") contains all of the digits from 1 to 9. The puzzle setter provides a partially completed grid, which for a well-posed puzzle has a unique solution.

In order to solve the problem we'll use Gurobi:


In [ ]:
import gurobipy as gurobi

The first step is to define decision variables. One possible choice is to use the binary variables $x_{ijk}$ which is equal to $1$ if, and only if, number $k$ is placed on cell $(i,j)$, with $i,j \in \{0,1,2,3,4,5,6,7,8\}$ and $k \in \{1,2,3,4,5,6,7,8,9\}$. I'm using $i$ as row index and $j$ as column index.

Note that I'm numbering the cells from 0 to 8 instead of from 1 to 9.

So, let's add them to a Gurobi model:


In [ ]:
model = gurobi.Model("SUDOKU")

x = {}
rows = 9
cols = 9
numbers = [1,2,3,4,5,6,7,8,9]

# Adding all the variables:
for i in range(rows):
    for j in range(cols):
        for k in numbers:
            varName = 'x_' +str(i) + str(j) + str(k) 
            x[i,j,k] = model.addVar(vtype = gurobi.GRB.BINARY, name = varName)  # Adding the binary variables

model.update()
model

Executing the previous cell, it is possible to see that the model, so far, has zero contraints and 729 variables.

Now it is time to start modelling the rules of the game. The first one is: the grid must be full, meaning that each cell of the 9x9 grid must have a number from 1 to 9 assigned to it. This is modeled by the following set of constraints:

$$\sum_{k=1}^9 x_{ijk} = 1, \qquad i \in \{0, ...,8\} \quad and \quad j \in \{0, ..., 8\}$$

Time to add it to our model:


In [ ]:
for i in range(rows):
    for j in range(cols):
        model.addConstr(gurobi.quicksum(x[i,j,k] for k in numbers) == 1) 
        
model.update()
model

Now our model has 81 constraints: one for each cell.

The next constraint we are going to model is each column contains all of the digits from 1 to 9 or, alternatively, each column must not contain repeated numbers. This give us:

$$\sum_{i=0}^8 x_{ijk} = 1, \qquad j \in \{0, ..., 8\} \quad and \quad k \in \{1, ..., 9\}.$$

Translating to Python:


In [ ]:
for j in range(cols): # No repeating numbers in the columns
    for k in numbers:
        model.addConstr(gurobi.quicksum(x[i,j,k] for i in range(rows)) == 1)

model.update()
model

Which give us 81 more constraints, one for each combination of column and possible number.

The next constraint is each row must not contain repetead numbers. Which is somewhat similar to the previous one:

$$\sum_{j=0}^8 x_{ijk} = 1, \qquad i \in \{0, ..., 8\} \quad and \quad k \in \{1, ..., 9\}.$$

Adding this to the model:


In [ ]:
for i in range(rows):
    for k in numbers:
        model.addConstr(gurobi.quicksum(x[i,j,k] for j in range(cols)) == 1)       
        
model.update()
model

As expected, we got 81 more constraints.

The last constraint is each of the nine 3x3 subgrid must not contain repeated numbers. This set of constraints is a little tricker and we are going to use some extra index control to help us move along the grid.

$$\sum_{i=0}^2 \sum_{j=0}^2 x_{i+3m, j+3n, k} = 1, \qquad k \in \{1, ..., 9\} \quad m \in \{0,1,2\} \quad n \in \{0,1,2\}$$

Note that the summations in $i$ and $j$, are related to a 3x3 subgrid and the indexes $m$ and $n$ are used to tell which of the nine subgrids is being considered.

Now, putting this into our model:


In [ ]:
for k in numbers:
    for m in range(3): 
        for n in range(3):
            model.addConstr(gurobi.quicksum(x[i + 3*m ,j + 3*n ,k] for i in range(3) for j in range(3)) == 1)
        
model.update()
model

And, with a total of 324 constraints and 729 variables we can solve our model and hope everything will be all right.


In [ ]:
model.optimize()

And a very lazy way to print the solution...


In [ ]:
print 'Optimal solution? ' + str(model.status == gurobi.GRB.status.OPTIMAL)

toPrint = ''
for i in range(rows):
    if i % 3 == 0:
        toPrint += '\n'
    
    line = ''
    for j in range(cols):
        if j % 3 == 0:
            line += '\t'
        
        for k in numbers:
            if x[i,j,k].x > 0.5:
                line += str(k) + ' '
        
    toPrint += line + '\n'

print toPrint

Note that the model only contemplates the general rules of the game. As a homework you can add constraints to solve a particular boards, forcing some numbers to be in particular positions.