Sudoku Solver with Linopy

Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Sudoku Solver with Linopy#

The linopy package is good at tracking indices and writing optimization problems as vectorized functions (i.e., functions that operate on all rows of a column at once). This notebook adapts the Medium article Creating Sudoku Solver with Python and Pyomo in Easy Steps by Dhanalakota Mohan for linopy to demonstrate how indices (i.e., dimensions and coordinates) work.

[1]:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

import linopy

pd.options.display.float_format = "{:.0f}".format

The Puzzle#

Sudoku is a 9x9 grid where each of the digits 1 through 9 appears once and only once in each row, column, and 3x3 square. The puzzle starts with values already filled into some cells, which we call hints. Here is an example puzzle with its hints as a pandas DataFrame.

[2]:

# look at the hints for the puzzle
puzzle_hints = [
    (1, 7, 2),
    (2, 2, 8),
    (2, 6, 7),
    (2, 8, 9),
    (3, 1, 6),
    (3, 3, 2),
    (3, 7, 5),
    (4, 2, 7),
    (4, 5, 6),
    (5, 4, 9),
    (5, 6, 1),
    (6, 5, 2),
    (6, 8, 4),
    (7, 3, 5),
    (7, 7, 6),
    (7, 9, 3),
    (8, 2, 9),
    (8, 4, 4),
    (8, 8, 7),
    (9, 3, 6),
]
puzzle_hints = pd.DataFrame(puzzle_hints, columns=["row", "column", "digit"])
puzzle_hints_piv = puzzle_hints.pivot(
    index="row", columns="column", values="digit"
).replace(np.nan, "-")
puzzle_hints_piv

[2]:
column 1 2 3 4 5 6 7 8 9
row
1 - - - - - - 2 - -
2 - 8 - - - 7 - 9 -
3 6 - 2 - - - 5 - -
4 - 7 - - 6 - - - -
5 - - - 9 - 1 - - -
6 - - - - 2 - - 4 -
7 - - 5 - - - 6 - 3
8 - 9 - 4 - - - 7 -
9 - - 6 - - - - - -

Defining Dimensions and Coordinates in Xarray#

In Pandas parlance, the index are the row labels and columns are the column labels. When converting to Xarray, it is handy to name the index and columns with df.index.name = "row_name" and df.columns.name = "column_name".

In Xarray parlance, the Pandas df.index.name and df.columns.name are both called dimensions. Then the array of values of the index and columns are called coordinates. Think of latitude and longitude as dimensions with coordinate values. Xarray will infer the dimensions and coordinates from a Pandas DataFrame, provided columns and the index have names. The following cell shows how this looks.

[3]:

xr.DataArray(puzzle_hints_piv)

[3]:

<xarray.DataArray (row: 9, column: 9)> Size: 648B
array([['-', '-', '-', '-', '-', '-', 2.0, '-', '-'],
       ['-', 8.0, '-', '-', '-', 7.0, '-', 9.0, '-'],
       [6.0, '-', 2.0, '-', '-', '-', 5.0, '-', '-'],
       ['-', 7.0, '-', '-', 6.0, '-', '-', '-', '-'],
       ['-', '-', '-', 9.0, '-', 1.0, '-', '-', '-'],
       ['-', '-', '-', '-', 2.0, '-', '-', 4.0, '-'],
       ['-', '-', 5.0, '-', '-', '-', 6.0, '-', 3.0],
       ['-', 9.0, '-', 4.0, '-', '-', '-', 7.0, '-'],
       ['-', '-', 6.0, '-', '-', '-', '-', '-', '-']], dtype=object)
Coordinates:
  * row      (row) int64 72B 1 2 3 4 5 6 7 8 9
  * column   (column) int64 72B 1 2 3 4 5 6 7 8 9

In the Mohan article, the puzzle is solved with a clever reformulation. The obvious encoding — one integer variable per cell holding its digit — is hard to optimize: the rule “the nine cells of a row are all different” is not a linear constraint, and neither is “this cell equals 5”.

Instead we add a third dimension, digit, and make the variable a binary on/off switch: sudoku[row, column, digit] == 1 means “this cell holds this digit”. Every Sudoku rule then becomes a linear statement about sums of 0/1 values — for example, “exactly one digit per cell” is just “the nine switches for that cell sum to 1”. That linearity is what a linear solver needs, and it is the whole reason for the extra dimension.

Here is what an empty (all zeros) xarray DataArray of these switches looks like:

[4]:

range_coord = range(1, 10)
xr.DataArray(
    np.zeros(shape=(9, 9, 9)),
    dims=["row", "column", "digit"],
    coords={"row": range_coord, "column": range_coord, "digit": range_coord},
)

[4]:

<xarray.DataArray (row: 9, column: 9, digit: 9)> Size: 6kB
array([[[0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.]],

...

       [[0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0.]]])
Coordinates:
  * row      (row) int64 72B 1 2 3 4 5 6 7 8 9
  * column   (column) int64 72B 1 2 3 4 5 6 7 8 9
  * digit    (digit) int64 72B 1 2 3 4 5 6 7 8 9

And here is a graphical representation of the 9x9x9 cube of switches — each small cube is one (row, column, digit) variable:

[5]:

# Function to create a single cube at a specific location
def create_cube(x, y, z, size=1):
    # Define the vertices of a cube
    vertices = [
        [x, y, z],
        [x + size, y, z],
        [x + size, y + size, z],
        [x, y + size, z],
        [x, y, z + size],
        [x + size, y, z + size],
        [x + size, y + size, z + size],
        [x, y + size, z + size],
    ]
    # Define the 6 faces of the cube
    faces = [
        [vertices[j] for j in [0, 1, 2, 3]],
        [vertices[j] for j in [4, 5, 6, 7]],
        [vertices[j] for j in [0, 3, 7, 4]],
        [vertices[j] for j in [1, 2, 6, 5]],
        [vertices[j] for j in [0, 1, 5, 4]],
        [vertices[j] for j in [2, 3, 7, 6]],
    ]
    return faces


# Create a 3D plot
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111, projection="3d")

# Colors for the slices
colors = plt.cm.viridis(np.linspace(0, 1, 9))

# Plot each cube
for i in range(9):
    for j in range(9):
        for k in range(9):
            faces = create_cube(i, j, k)
            poly3d = Poly3DCollection(
                faces, color=colors[j], edgecolor="k", linewidth=0.1
            )  # Color by digit dimension
            ax.add_collection3d(poly3d)

# Set custom labels
ax.set_xlabel("Column")
ax.set_ylabel("Digit")
ax.set_zlabel("")
ax.set_title("Sudoku Dimensions and Coordinates")

# Customize ticks and labels
ticks = np.arange(9)
labels = np.arange(1, 10)  # Labels from 1 to 9
ax.set_xticks(ticks)
ax.set_yticks(ticks)
ax.set_zticks([])  # Remove z-axis ticks
ax.set_xticklabels(labels)
ax.set_yticklabels(labels)

# Set the aspect ratio to be equal and adjust limits to fit the plot
ax.set_box_aspect([1, 1, 1])  # aspect ratio is 1:1:1
ax.set_xlim([0, 9])
ax.set_ylim([0, 9])
ax.set_zlim([0, 9])

# Adjust layout to ensure labels are visible
plt.subplots_adjust(left=0.5, right=0.8, top=0.9, bottom=0.0)

# Manually add z-axis labels
for z in range(9):
    ax.text(x=-1, y=-1, z=z + 0.5, s=9 - z)
ax.text(x=-3, y=-1, z=4.5, s="Row", ha="center")

plt.tight_layout()

/tmp/ipykernel_3197/4151666252.py:72: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all Axes decorations.
  plt.tight_layout()
_images/sudoku_9_1.svg

Indices for the 3x3 Squares#

If you play Sudoku, you know each 3x3 square must also contain each digit once and only once. To express that, we first need to know which square every cell belongs to. We compute a square number for each (row, column) and store it as an xarray; later the model groups cells by this number to build the square constraints. Here is how the squares are numbered:

[6]:

# Define the function to get the 3x3 square index
def get_square_index(row, col):
    return (row // 3) * 3 + (col // 3)


# Create the grid with row, column, and square indices
square_index = pd.DataFrame(
    {
        "row": np.repeat(range_coord, 9),
        "column": list(range_coord) * 9,
        "square": [
            get_square_index(row - 1, col - 1) + 1
            for row in range_coord
            for col in range_coord
        ],
    }
)
square_index = square_index.pivot(index="row", columns="column", values="square")
square_index

[6]:
column 1 2 3 4 5 6 7 8 9
row
1 1 1 1 2 2 2 3 3 3
2 1 1 1 2 2 2 3 3 3
3 1 1 1 2 2 2 3 3 3
4 4 4 4 5 5 5 6 6 6
5 4 4 4 5 5 5 6 6 6
6 4 4 4 5 5 5 6 6 6
7 7 7 7 8 8 8 9 9 9
8 7 7 7 8 8 8 9 9 9
9 7 7 7 8 8 8 9 9 9

We put that into an xarray. Each cell now carries the number of the 3x3 square it belongs to:

[7]:

square_index = xr.DataArray(square_index).rename("square")
square_index

[7]:

<xarray.DataArray 'square' (row: 9, column: 9)> Size: 648B
array([[1, 1, 1, 2, 2, 2, 3, 3, 3],
       [1, 1, 1, 2, 2, 2, 3, 3, 3],
       [1, 1, 1, 2, 2, 2, 3, 3, 3],
       [4, 4, 4, 5, 5, 5, 6, 6, 6],
       [4, 4, 4, 5, 5, 5, 6, 6, 6],
       [4, 4, 4, 5, 5, 5, 6, 6, 6],
       [7, 7, 7, 8, 8, 8, 9, 9, 9],
       [7, 7, 7, 8, 8, 8, 9, 9, 9],
       [7, 7, 7, 8, 8, 8, 9, 9, 9]])
Coordinates:
  * row      (row) int64 72B 1 2 3 4 5 6 7 8 9
  * column   (column) int64 72B 1 2 3 4 5 6 7 8 9

Model#

Now we build the model. Initialize a linopy model:

[8]:

model = linopy.Model()

Make the binary sudoku variable: a 9x9x9 cube of on/off switches over the row, column, and digit dimensions.

[9]:

sudoku = model.add_variables(
    name="sudoku",
    coords=[range_coord, range_coord, range_coord],
    dims=["row", "column", "digit"],
    binary=True,
)

It is handy to look at what linopy is doing. The code defines a single variable, but it is indexed on a 9x9x9 cube, so it really holds 729 binary variables.

[10]:

sudoku

[10]:

Variable (row: 9, column: 9, digit: 9)
--------------------------------------
[1, 1, 1]: sudoku[1, 1, 1] ∈ {0, 1}
[1, 1, 2]: sudoku[1, 1, 2] ∈ {0, 1}
[1, 1, 3]: sudoku[1, 1, 3] ∈ {0, 1}
[1, 1, 4]: sudoku[1, 1, 4] ∈ {0, 1}
[1, 1, 5]: sudoku[1, 1, 5] ∈ {0, 1}
[1, 1, 6]: sudoku[1, 1, 6] ∈ {0, 1}
[1, 1, 7]: sudoku[1, 1, 7] ∈ {0, 1}
                ...
[9, 9, 3]: sudoku[9, 9, 3] ∈ {0, 1}
[9, 9, 4]: sudoku[9, 9, 4] ∈ {0, 1}
[9, 9, 5]: sudoku[9, 9, 5] ∈ {0, 1}
[9, 9, 6]: sudoku[9, 9, 6] ∈ {0, 1}
[9, 9, 7]: sudoku[9, 9, 7] ∈ {0, 1}
[9, 9, 8]: sudoku[9, 9, 8] ∈ {0, 1}
[9, 9, 9]: sudoku[9, 9, 9] ∈ {0, 1}

Now the constraints. Each one sums the on/off switches over a single dimension and sets that sum to 1, which forces exactly one switch to be on along that dimension. The three lines below read as:

  • summing over column → each digit appears once in every row,

  • summing over row → each digit appears once in every column,

  • summing over digit → each cell holds exactly one digit.

The last constraint is displayed by Jupyter, a handy way to see the structure linopy has built.

[11]:

model.add_constraints(sudoku.sum(dim=["column"]) == 1, name="row_digit_constraint")
model.add_constraints(sudoku.sum(dim=["row"]) == 1, name="column_digit_constraint")
model.add_constraints(sudoku.sum(dim=["digit"]) == 1, name="row_column_constraint")

[11]:

Constraint `row_column_constraint` [row: 9, column: 9]:
-------------------------------------------------------
[1, 1]: +1 sudoku[1, 1, 1] + 1 sudoku[1, 1, 2] + 1 sudoku[1, 1, 3] ... +1 sudoku[1, 1, 7] + 1 sudoku[1, 1, 8] + 1 sudoku[1, 1, 9] = 1.0
[1, 2]: +1 sudoku[1, 2, 1] + 1 sudoku[1, 2, 2] + 1 sudoku[1, 2, 3] ... +1 sudoku[1, 2, 7] + 1 sudoku[1, 2, 8] + 1 sudoku[1, 2, 9] = 1.0
[1, 3]: +1 sudoku[1, 3, 1] + 1 sudoku[1, 3, 2] + 1 sudoku[1, 3, 3] ... +1 sudoku[1, 3, 7] + 1 sudoku[1, 3, 8] + 1 sudoku[1, 3, 9] = 1.0
[1, 4]: +1 sudoku[1, 4, 1] + 1 sudoku[1, 4, 2] + 1 sudoku[1, 4, 3] ... +1 sudoku[1, 4, 7] + 1 sudoku[1, 4, 8] + 1 sudoku[1, 4, 9] = 1.0
[1, 5]: +1 sudoku[1, 5, 1] + 1 sudoku[1, 5, 2] + 1 sudoku[1, 5, 3] ... +1 sudoku[1, 5, 7] + 1 sudoku[1, 5, 8] + 1 sudoku[1, 5, 9] = 1.0
[1, 6]: +1 sudoku[1, 6, 1] + 1 sudoku[1, 6, 2] + 1 sudoku[1, 6, 3] ... +1 sudoku[1, 6, 7] + 1 sudoku[1, 6, 8] + 1 sudoku[1, 6, 9] = 1.0
[1, 7]: +1 sudoku[1, 7, 1] + 1 sudoku[1, 7, 2] + 1 sudoku[1, 7, 3] ... +1 sudoku[1, 7, 7] + 1 sudoku[1, 7, 8] + 1 sudoku[1, 7, 9] = 1.0
                ...
[9, 3]: +1 sudoku[9, 3, 1] + 1 sudoku[9, 3, 2] + 1 sudoku[9, 3, 3] ... +1 sudoku[9, 3, 7] + 1 sudoku[9, 3, 8] + 1 sudoku[9, 3, 9] = 1.0
[9, 4]: +1 sudoku[9, 4, 1] + 1 sudoku[9, 4, 2] + 1 sudoku[9, 4, 3] ... +1 sudoku[9, 4, 7] + 1 sudoku[9, 4, 8] + 1 sudoku[9, 4, 9] = 1.0
[9, 5]: +1 sudoku[9, 5, 1] + 1 sudoku[9, 5, 2] + 1 sudoku[9, 5, 3] ... +1 sudoku[9, 5, 7] + 1 sudoku[9, 5, 8] + 1 sudoku[9, 5, 9] = 1.0
[9, 6]: +1 sudoku[9, 6, 1] + 1 sudoku[9, 6, 2] + 1 sudoku[9, 6, 3] ... +1 sudoku[9, 6, 7] + 1 sudoku[9, 6, 8] + 1 sudoku[9, 6, 9] = 1.0
[9, 7]: +1 sudoku[9, 7, 1] + 1 sudoku[9, 7, 2] + 1 sudoku[9, 7, 3] ... +1 sudoku[9, 7, 7] + 1 sudoku[9, 7, 8] + 1 sudoku[9, 7, 9] = 1.0
[9, 8]: +1 sudoku[9, 8, 1] + 1 sudoku[9, 8, 2] + 1 sudoku[9, 8, 3] ... +1 sudoku[9, 8, 7] + 1 sudoku[9, 8, 8] + 1 sudoku[9, 8, 9] = 1.0
[9, 9]: +1 sudoku[9, 9, 1] + 1 sudoku[9, 9, 2] + 1 sudoku[9, 9, 3] ... +1 sudoku[9, 9, 7] + 1 sudoku[9, 9, 8] + 1 sudoku[9, 9, 9] = 1.0

Finally, the 3x3 square rule: each digit appears once and only once in every square. We groupby the cells by their square_index and sum within each group. Because the grouper is 2D (indexed by row and column) while the variable also carries digit, the grouped sum collapses row and column into a single square dimension and keeps digit — giving one constraint per (square, digit), all in one vectorized call.

[12]:

model.add_constraints(
    sudoku.groupby(square_index).sum() == 1,
    name="square_digit_constraint",
)

[12]:

Constraint `square_digit_constraint` [digit: 9, square: 9]:
-----------------------------------------------------------
[1, 1]: +1 sudoku[1, 1, 1] + 1 sudoku[1, 2, 1] + 1 sudoku[1, 3, 1] ... +1 sudoku[3, 1, 1] + 1 sudoku[3, 2, 1] + 1 sudoku[3, 3, 1] = 1.0
[1, 2]: +1 sudoku[1, 4, 1] + 1 sudoku[1, 5, 1] + 1 sudoku[1, 6, 1] ... +1 sudoku[3, 4, 1] + 1 sudoku[3, 5, 1] + 1 sudoku[3, 6, 1] = 1.0
[1, 3]: +1 sudoku[1, 7, 1] + 1 sudoku[1, 8, 1] + 1 sudoku[1, 9, 1] ... +1 sudoku[3, 7, 1] + 1 sudoku[3, 8, 1] + 1 sudoku[3, 9, 1] = 1.0
[1, 4]: +1 sudoku[4, 1, 1] + 1 sudoku[4, 2, 1] + 1 sudoku[4, 3, 1] ... +1 sudoku[6, 1, 1] + 1 sudoku[6, 2, 1] + 1 sudoku[6, 3, 1] = 1.0
[1, 5]: +1 sudoku[4, 4, 1] + 1 sudoku[4, 5, 1] + 1 sudoku[4, 6, 1] ... +1 sudoku[6, 4, 1] + 1 sudoku[6, 5, 1] + 1 sudoku[6, 6, 1] = 1.0
[1, 6]: +1 sudoku[4, 7, 1] + 1 sudoku[4, 8, 1] + 1 sudoku[4, 9, 1] ... +1 sudoku[6, 7, 1] + 1 sudoku[6, 8, 1] + 1 sudoku[6, 9, 1] = 1.0
[1, 7]: +1 sudoku[7, 1, 1] + 1 sudoku[7, 2, 1] + 1 sudoku[7, 3, 1] ... +1 sudoku[9, 1, 1] + 1 sudoku[9, 2, 1] + 1 sudoku[9, 3, 1] = 1.0
                ...
[9, 3]: +1 sudoku[1, 7, 9] + 1 sudoku[1, 8, 9] + 1 sudoku[1, 9, 9] ... +1 sudoku[3, 7, 9] + 1 sudoku[3, 8, 9] + 1 sudoku[3, 9, 9] = 1.0
[9, 4]: +1 sudoku[4, 1, 9] + 1 sudoku[4, 2, 9] + 1 sudoku[4, 3, 9] ... +1 sudoku[6, 1, 9] + 1 sudoku[6, 2, 9] + 1 sudoku[6, 3, 9] = 1.0
[9, 5]: +1 sudoku[4, 4, 9] + 1 sudoku[4, 5, 9] + 1 sudoku[4, 6, 9] ... +1 sudoku[6, 4, 9] + 1 sudoku[6, 5, 9] + 1 sudoku[6, 6, 9] = 1.0
[9, 6]: +1 sudoku[4, 7, 9] + 1 sudoku[4, 8, 9] + 1 sudoku[4, 9, 9] ... +1 sudoku[6, 7, 9] + 1 sudoku[6, 8, 9] + 1 sudoku[6, 9, 9] = 1.0
[9, 7]: +1 sudoku[7, 1, 9] + 1 sudoku[7, 2, 9] + 1 sudoku[7, 3, 9] ... +1 sudoku[9, 1, 9] + 1 sudoku[9, 2, 9] + 1 sudoku[9, 3, 9] = 1.0
[9, 8]: +1 sudoku[7, 4, 9] + 1 sudoku[7, 5, 9] + 1 sudoku[7, 6, 9] ... +1 sudoku[9, 4, 9] + 1 sudoku[9, 5, 9] + 1 sudoku[9, 6, 9] = 1.0
[9, 9]: +1 sudoku[7, 7, 9] + 1 sudoku[7, 8, 9] + 1 sudoku[7, 9, 9] ... +1 sudoku[9, 7, 9] + 1 sudoku[9, 8, 9] + 1 sudoku[9, 9, 9] = 1.0

A Sudoku puzzle starts with some cells already filled in — the hints. We fix these too. Each hint is a (row, column, digit) triple; comparing each hint coordinate against the matching grid axis broadcasts a True wherever they line up, and & keeps only the cells that match on all three. Reducing with .any("hint") unions the 20 hints into one boolean mask, which we hand to add_constraints to pin those switches to 1 — again in a single vectorized call.

[13]:

hints = puzzle_hints.to_xarray().rename(index="hint")
hint_mask = (
    (sudoku.coords["row"] == hints.row)
    & (sudoku.coords["column"] == hints.column)
    & (sudoku.coords["digit"] == hints.digit)
).any("hint")
model.add_constraints(sudoku == 1, mask=hint_mask, name="hints")

[13]:

Constraint `hints` [row: 9, column: 9, digit: 9] - 709 masked entries:
----------------------------------------------------------------------
[1, 1, 1]: None
[1, 1, 2]: None
[1, 1, 3]: None
[1, 1, 4]: None
[1, 1, 5]: None
[1, 1, 6]: None
[1, 1, 7]: None
                ...
[9, 9, 3]: None
[9, 9, 4]: None
[9, 9, 5]: None
[9, 9, 6]: None
[9, 9, 7]: None
[9, 9, 8]: None
[9, 9, 9]: None

This is really a feasibility problem: the constraints alone pin down the unique solution, so solving just means finding any point that satisfies them. An objective is optional here — we add a dummy one (maximize the sum of the switches) purely because it can nudge the solver; its value carries no meaning.

[14]:

model.add_objective(sudoku.sum(), sense="max")

Solve#

Solve the model with the default solver (HiGHS).

[15]:

model.solve()

Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-_o3te90z.lp
Reading time = 0.00 seconds
obj: 344 rows, 729 columns, 2936 nonzeros
Gurobi Optimizer version 13.0.2 build v13.0.2rc1 (linux64 - "Ubuntu 24.04 LTS")

CPU model: AMD EPYC 7R13 Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads

Optimize a model with 344 rows, 729 columns and 2936 nonzeros (Max)
Model fingerprint: 0x95a93691
Model has 729 linear objective coefficients
Variable types: 0 continuous, 729 integer (729 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]

Presolve removed 344 rows and 729 columns
Presolve time: 0.01s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 2 available processors)

Solution count 1: 81

Optimal solution found (tolerance 1.00e-04)
Best objective 8.100000000000e+01, best bound 8.100000000000e+01, gap 0.0000%

Dual values of MILP couldn't be parsed

[15]:

('ok', 'optimal')

The solution comes out as an xarray representing a cube of binary variables.

[16]:

solution = model.solution.sudoku
solution

[16]:

<xarray.DataArray 'sudoku' (row: 9, column: 9, digit: 9)> Size: 6kB
array([[[0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0.]],

...

       [[0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1.]]])
Coordinates:
  * row      (row) int64 72B 1 2 3 4 5 6 7 8 9
  * column   (column) int64 72B 1 2 3 4 5 6 7 8 9
  * digit    (digit) int64 72B 1 2 3 4 5 6 7 8 9

Finally, turn the cube of 0/1 switches back into a 9x9 grid of digits. Since exactly one switch is on per cell, multiplying each switch by its own digit coordinate and summing over digit recovers the chosen digit — the exact inverse of the on/off encoding we started with.

[17]:

final_sudoku_grid = (solution * solution.digit).sum("digit").astype(int)

result = pd.DataFrame(
    data=final_sudoku_grid.values,
    columns=puzzle_hints_piv.columns,
    index=puzzle_hints_piv.index,
)
result

[17]:
column 1 2 3 4 5 6 7 8 9
row
1 9 5 7 6 1 3 2 8 4
2 4 8 3 2 5 7 1 9 6
3 6 1 2 8 4 9 5 3 7
4 1 7 8 3 6 4 9 5 2
5 5 2 4 9 7 1 3 6 8
6 3 6 9 5 2 8 7 4 1
7 8 4 5 7 9 2 6 1 3
8 2 9 1 4 3 6 8 7 5
9 7 3 6 1 8 5 4 2 9

The grid above is the completed puzzle, filled in from just the 20 hints we started with.

Check the result#

As a cross-check, every row, column, and 3x3 square must contain each digit once, so each of their sums must equal 1 + 2 + … + 9 = 45.

[18]:

row_sums = result.sum(axis=1)
column_sums = result.sum(axis=0)
square_sums = result.stack().groupby(square_index.to_series()).sum()
sums = pd.DataFrame({"row": row_sums, "column": column_sums, "square": square_sums})
sums

[18]:
row column square
1 45 45 45
2 45 45 45
3 45 45 45
4 45 45 45
5 45 45 45
6 45 45 45
7 45 45 45
8 45 45 45
9 45 45 45 
[19]:

sums.eq(45).all().all()

[19]:

True

Conclusion#

By encoding each cell as a set of binary on/off switches, every Sudoku rule became a linear constraint on sums of 0/1 values, and linopy’s dimension-aware sum and groupby expressed all 729 variables and their constraints in just a handful of vectorized calls.