Solve Sudoku Puzzles via Integer Programming: Problem-Based

Open Live Script

This example shows how to solve a Sudoku puzzle using binary integer programming. For the solver-based approach, see Solve Sudoku Puzzles via Integer Programming: Solver-Based.

You probably have seen Sudoku puzzles. A puzzle is to fill a 9-by-9 grid with integers from 1 through 9 so that each integer appears only once in each row, column, and major 3-by-3 square. The grid is partially populated with clues, and your task is to fill in the rest of the grid.

Initial Puzzle

Here is a data matrix B of clues. The first row, B(1,2,2), means row 1, column 2 has a clue 2. The second row, B(1,5,3), means row 1, column 5 has a clue 3. Here is the entire matrix B.

B = [1,2,2;
    1,5,3;
    1,8,4;
    2,1,6;
    2,9,3;
    3,3,4;
    3,7,5;
    4,4,8;
    4,6,6;
    5,1,8;
    5,5,1;
    5,9,6;
    6,4,7;
    6,6,5;
    7,3,7;
    7,7,6;
    8,1,4;
    8,9,8;
    9,2,3;
    9,5,4;
    9,8,2];

drawSudoku(B) % For the listing of this program, see the end of this example.

Figure contains an axes object. The hidden axes object contains 30 objects of type rectangle, text.

This puzzle, and an alternative MATLAB® solution technique, was featured in Cleve's Corner in 2009.

There are many approaches to solving Sudoku puzzles manually, as well as many programmatic approaches. This example shows a straightforward approach using binary integer programming.

This approach is particularly simple because you do not give a solution algorithm. Just express the rules of Sudoku, express the clues as constraints on the solution, and then MATLAB produces the solution.

Binary Integer Programming Approach

The key idea is to transform a puzzle from a square 9-by-9 grid to a cubic 9-by-9-by-9 array of binary values (0 or 1). Think of the cubic array as being 9 square grids stacked on top of each other, where each layer corresponds to an integer from 1 through 9. The top grid, a square layer of the array, has a 1 wherever the solution or clue has a 1. The second layer has a 1 wherever the solution or clue has a 2. The ninth layer has a 1 wherever the solution or clue has a 9.

This formulation is precisely suited for binary integer programming.

The objective function is not needed here, and might as well be a constant term 0. The problem is really just to find a feasible solution, meaning one that satisfies all the constraints. However, for tie breaking in the internals of the integer programming solver, giving increased solution speed, use a nonconstant objective function.

Express the Rules for Sudoku as Constraints

Suppose a solution $x$ is represented in a 9-by-9-by-9 binary array. What properties does $x$ have? First, each square in the 2-D grid (i,j) has exactly one value, so there is exactly one nonzero element among the 3-D array entries $x (i, j, 1), . . ., x (i, j, 9)$ . In other words, for every $i$ and $j$ ,

$\sum_{k = 1}^{9} x (i, j, k) = 1 .$

Similarly, in each row $i$ of the 2-D grid, there is exactly one value out of each of the digits from 1 to 9. In other words, for each $i$ and $k$ ,

$\sum_{j = 1}^{9} x (i, j, k) = 1 .$

And each column $j$ in the 2-D grid has the same property: for each $j$ and $k$ ,

$\sum_{i = 1}^{9} x (i, j, k) = 1 .$

The major 3-by-3 grids have a similar constraint. For the grid elements $1 \leq i \leq 3$ and $1 \leq j \leq 3$ , and for each $1 \leq k \leq 9$ ,

$\sum_{i = 1}^{3} \sum_{j = 1}^{3} x (i, j, k) = 1 .$

To represent all nine major grids, just add 3 or 6 to each $i$ and $j$ index:

$\sum_{i = 1}^{3} \sum_{j = 1}^{3} x (i + U, j + V, k) = 1,$ where $U, V ϵ {0, 3, 6} .$

Express Clues

Each initial value (clue) can be expressed as a constraint. Suppose that the $(i, j)$ clue is $m$ for some $1 \leq m \leq 9$ . Then $x (i, j, m) = 1$ . The constraint $\sum_{k = 1}^{9} x (i, j, k) = 1$ ensures that all other $x (i, j, k) = 0$ for $k \neq m$ .

Sudoku in Optimization Problem Form

Create an optimization variable x that is binary and of size 9-by-9-by-9.

x = optimvar('x',9,9,9,'Type','integer','LowerBound',0,'UpperBound',1);

Create an optimization problem with a rather arbitrary objective function. The objective function can help the solver by destroying the inherent symmetry of the problem.

sudpuzzle = optimproblem;
mul = ones(1,1,9);
mul = cumsum(mul,3);
sudpuzzle.Objective = sum(sum(sum(x,1),2).*mul);

Represent the constraints that the sums of x in each coordinate direction are one.

sudpuzzle.Constraints.consx = sum(x,1) == 1;
sudpuzzle.Constraints.consy = sum(x,2) == 1;
sudpuzzle.Constraints.consz = sum(x,3) == 1;

Create the constraints that the sums of the major grids sum to one as well.

majorg = optimconstr(3,3,9);

for u = 1:3
    for v = 1:3
        arr = x(3*(u-1)+1:3*(u-1)+3,3*(v-1)+1:3*(v-1)+3,:);
        majorg(u,v,:) = sum(sum(arr,1),2) == ones(1,1,9);
    end
end
sudpuzzle.Constraints.majorg = majorg;

Include the initial clues by setting lower bounds of 1 at the clue entries. This setting fixes the value of the corresponding entry to be 1, and so sets the solution at each clued value to be the clue entry.

for u = 1:size(B,1)
    x.LowerBound(B(u,1),B(u,2),B(u,3)) = 1;
end

Solve the Sudoku puzzle.

sudsoln = solve(sudpuzzle)

Solving problem using intlinprog.
Running HiGHS 1.7.0: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+00, 9e+00]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Presolving model
221 rows, 246 cols, 990 nonzeros  0s
42 rows, 28 cols, 192 nonzeros  0s
22 rows, 18 cols, 86 nonzeros  0s
Objective function is integral with scale 1

Solving MIP model with:
   22 rows
   18 cols (18 binary, 0 integer, 0 implied int., 0 continuous)
   86 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   372             inf                  inf        0      0      0         0     0.0s
         0       0         0   0.00%   405             inf                  inf        0      0      2         7     0.0s

Solving report
  Status            Optimal
  Primal bound      405
  Dual bound        405
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    405 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.01 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     7 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)

Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.

sudsoln = struct with fields:
    x: [9x9x9 double]

Round the solution to ensure that all entries are integers, and display the solution.

sudsoln.x = round(sudsoln.x);

y = ones(size(sudsoln.x));
for k = 2:9
    y(:,:,k) = k; % multiplier for each depth k
end
S = sudsoln.x.*y; % multiply each entry by its depth
S = sum(S,3); % S is 9-by-9 and holds the solved puzzle
drawSudoku(S)

Figure contains an axes object. The hidden axes object contains 90 objects of type rectangle, text.

You can easily check that the solution is correct.

Function to Draw the Sudoku Puzzle

type drawSudoku

function drawSudoku(B)
% Function for drawing the Sudoku board

%   Copyright 2014 The MathWorks, Inc. 


figure;hold on;axis off;axis equal % prepare to draw
rectangle('Position',[0 0 9 9],'LineWidth',3,'Clipping','off') % outside border
rectangle('Position',[3,0,3,9],'LineWidth',2) % heavy vertical lines
rectangle('Position',[0,3,9,3],'LineWidth',2) % heavy horizontal lines
rectangle('Position',[0,1,9,1],'LineWidth',1) % minor horizontal lines
rectangle('Position',[0,4,9,1],'LineWidth',1)
rectangle('Position',[0,7,9,1],'LineWidth',1)
rectangle('Position',[1,0,1,9],'LineWidth',1) % minor vertical lines
rectangle('Position',[4,0,1,9],'LineWidth',1)
rectangle('Position',[7,0,1,9],'LineWidth',1)

% Fill in the clues
%
% The rows of B are of the form (i,j,k) where i is the row counting from
% the top, j is the column, and k is the clue. To place the entries in the
% boxes, j is the horizontal distance, 10-i is the vertical distance, and
% we subtract 0.5 to center the clue in the box.
%
% If B is a 9-by-9 matrix, convert it to 3 columns first

if size(B,2) == 9 % 9 columns
    [SM,SN] = meshgrid(1:9); % make i,j entries
    B = [SN(:),SM(:),B(:)]; % i,j,k rows
end

for ii = 1:size(B,1)
    text(B(ii,2)-0.5,9.5-B(ii,1),num2str(B(ii,3)))
end

hold off

end