Solving Killer Sudoku as a Mixed Integer Program
In killer sudoku you have a 9x9 grid of numbers. One has to place numbers between 1 and 9 in each grid cell such that the rules of Sudoku are satisfied. Additionally, each square is part of a contiguous group of cells around it with a number, and the sum of the cells in that group must add up to this number.
This notebook-turned-blog-post formulates killer sudoku as a mixed integer program, then solves the following example using Gurobi in Python.
Methodology
We will use the official Gurobi Python bindings.
import gurobipy as gp
Killer sudoku data
Each line contains a list of cells (x,y) and what that group must sum to.
data_killer = [
([(1,1),(1,2)],15),
([(2,1),(2,2)],8),
([(3,1),(3,2)],3),
([(4,1),(5,1)],16),
([(6,1),(7,1)],5),
([(8,1),(9,1),(8,2),(9,2)],24),
([(4,2),(3,3),(3,4),(4,3)],18),
([(5,2),(5,3),(5,4),(5,5),(5,6)],25),
([(6,2),(6,3)],5),
([(7,2),(7,3)],7),
([(1,3),(1,4),(1,5),(1,6)],22),
([(2,3),(2,4)],5),
([(8,3),(9,3)],12),
([(4,4),(4,5)],9),
([(6,4),(7,4),(8,4)],18),
([(9,4),(9,5)],14),
([(2,5),(3,5),(3,6),(4,6)],24),
([(6,5),(6,6)],16),
([(7,5),(7,6),(8,5),(8,6)],14),
([(2,6),(1,7),(2,7),(1,8),(1,9)],21),
([(9,6),(9,7)],10),
([(3,7),(4,7)],17),
([(5,7),(6,7),(7,7),(8,7)],11),
([(2,8),(2,9)],11),
([(3,8),(3,9)],11),
([(4,8),(4,9)],9),
([(5,8),(5,9),(6,9),(7,9)],22),
([(6,8),(7,8)],15),
([(8,8),(9,8)],9),
([(8,9),(9,9)],9),
]
Check
Visually check we’ve got the right groups and sums:
for squares, sum_ in data_killer[6:9]:
print(f"Sum of cells marked below must add up to {sum_}")
out = ""
for x in range(1,10):
for y in range(1,10):
out += "x" if (y,x) in squares else "."
out += "\n"
print(out)
print()
Sum of cells marked below must add up to 18
.........
...x.....
..xx.....
..x......
.........
.........
.........
.........
.........
Sum of cells marked below must add up to 25
.........
....x....
....x....
....x....
....x....
....x....
.........
.........
.........
Sum of cells marked below must add up to 5
.........
.....x...
.....x...
.........
.........
.........
.........
.........
.........
Also check each cell is in exactly one group.
sqs = set()
for squares, sum_ in data_killer:
for s in squares:
assert s not in sqs, "Repeated killer square"
sqs.add(s)
Sudoku rules
Each of these groups must contain the numbers from 1 to 9 exactly once.
data_sudoku = []
for x in range(1,10,3):
for y in range(1,10,3):
data_sudoku.append([(x+i,y+j) for i in range(3) for j in range(3)])
for q in range(1,10):
data_sudoku.append([(q,w) for w in range(10)])
data_sudoku.append([(w,q) for w in range(10)])
Check
for squares in data_sudoku[:2]:
out = ""
for x in range(1,10):
for y in range(1,10):
out += "x" if (y,x) in squares else "."
out += "\n"
print(out)
print()
xxx......
xxx......
xxx......
.........
.........
.........
.........
.........
.........
.........
.........
.........
xxx......
xxx......
xxx......
.........
.........
.........
Build a MIP model
m = gp.Model()
Set parameter Username
Academic license - for non-commercial use only - expires 2023-11-08
Create binary variables for each grid square and each number from 1 to 9, so grid[4-1][2-1][7-1] is one if and only if cell (4,2) has value 7 (note the minus one for zero-based indexing of the arrays).
grid = [[[m.addVar(vtype="b", name=f"{x},{y}={i}") for i in range(1,10)] for x in range(1,10)] for y in range(1,10)]
Now add constraints to ensure each grid cell has exactly one value.
for cell_x in range(1,10):
for cell_y in range(1,10):
m.addConstr(sum(grid[cell_x-1][cell_y-1]) == 1.)
The following constraint ensures each sudoku group (line across or down, or 3x3 sudoku square) contains each number from 1 to 9 at least once.
for cells in data_sudoku:
for i in range(1,10):
m.addConstr(sum(grid[x-1][y-1][i-1] for (x,y) in cells) >= 1.)
The following constraint ensures that the sum of the cells in each killer sudoku square adds up to the number indicated.
for cells, sum_ in data_killer:
m.addConstr(sum(i*grid[x-1][y-1][i-1] for i in range(1,10) for (x,y) in cells) == sum_)
Now call Gurobi to find a solution!
m.optimize()
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (linux64)
CPU model: AMD Ryzen 7 5800X 8-Core Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 354 rows, 729 columns and 3645 nonzeros
Model fingerprint: 0x4a587d32
Variable types: 0 continuous, 729 integer (729 binary)
Coefficient statistics:
Matrix range [1e+00, 9e+00]
Objective range [0e+00, 0e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+01]
Presolve removed 16 rows and 140 columns
Presolve time: 0.01s
Presolved: 338 rows, 589 columns, 2848 nonzeros
Variable types: 0 continuous, 589 integer (589 binary)
Root relaxation: objective 0.000000e+00, 1179 iterations, 0.03 seconds (0.04 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 0.00000 0 201 - 0.00000 - - 0s
0 0 0.00000 0 216 - 0.00000 - - 0s
* 0 0 0 0.0000000 0.00000 0.00% - 0s
Cutting planes:
Gomory: 19
Cover: 15
Clique: 26
MIR: 3
StrongCG: 2
GUB cover: 16
Zero half: 13
RLT: 15
Explored 1 nodes (4242 simplex iterations) in 0.17 seconds (0.23 work units)
Thread count was 16 (of 16 available processors)
Solution count 1: 0
Optimal solution found (tolerance 1.00e-04)
Best objective 0.000000000000e+00, best bound 0.000000000000e+00, gap 0.0000%
Print out the solution
out = ""
for y in range(1,10):
for x in range(1,10):
for i in range(1,10):
if grid[x-1][y-1][i-1].X >= 0.5:
out += " " + str(i)
if x in [3, 6]: out += " |"
if y in [3, 6]: out += "\n-------+-------+-------"
out += "\n"
print(out)
6 5 1 | 7 9 3 | 2 4 8
9 3 2 | 6 8 4 | 1 5 7
8 4 7 | 2 5 1 | 6 9 3
-------+-------+-------
2 1 3 | 8 4 5 | 7 6 9
7 8 4 | 1 6 9 | 3 2 5
5 6 9 | 3 2 7 | 8 1 4
-------+-------+-------
4 7 8 | 9 1 2 | 5 3 6
3 2 5 | 4 7 6 | 9 8 1
1 9 6 | 5 3 8 | 4 7 2
Caveats
A MIP solver is not really the right tool to solve sudoku (or killer sudoku). There are purpose-built constraint-solvers that are better and faster at solving binary programs like the one formulated here. I ran into some slight numerical issues with Gurobi and had to tweaks some constraints ever so slightly to get it to solve this instance, but in the end it ran surprisingly fast and of course gave a valid solution!