# 10.7 SUDOKU¶

SUDOKU is a famous simple yet mind-blowing game. 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. For more information see http://en.wikipedia.org/wiki/Sudoku. Here is a simple example:

In a more general setting we are given a grid of dimension $$n\times n$$, with $$n=m^2, m\in \natnumb$$. Each cell $$(i,j)$$ must be filled with an integer $$y_{ij} \in [1,n]$$. Along each row and each column there must be no repetitions. No repetitions are allowed also in each sub-grid with corners $$\{(mt,ml),\ (m(t+1)-1,m(l+1)-1)\}$$, for $$t,l=0,\ldots,m-1$$ (we index cells from $$(0,0)$$).

In general, each SUDOKU instance comes with a set $$F$$ of predetermined values which:

• reduce the complexity of the game by removing symmetries and guiding the initial moves of the player;
• ensure that there will be a unique solution.

We represent the set $$F$$ as list of triplets $$(i,j,v)$$, meaning that the cell $$(i,j)$$ contains the value $$v$$.

Note that SUDOKU is a feasibility problem. A typical Integer Programming formulation is straightforward: let $$x_{ijk}$$ be a binary variable that takes value $$1$$ if $$k$$ is written in cell $$(i,j)$$. Then we look for a feasible solution of a system of constraints given below.

SUDOKU is a typical assignment problem. Its constraints are commonly found in optimization problems concerning scheduling or resource allocation. SUDOKU has also been a nice problem to fiddle with for many researchers in the optimization community. Indeed, its simple structure and the easy way in which the results can be tested make it a perfect test problem.

We will approach SUDOKU as a standard integer linear program, and we will show how easily and elegantly it can be implemented in Fusion.

Mathematical Formulation

In this section we formulate SUDOKU as a mixed-integer linear optimization problem. Let’s introduce a binary variable $$x_{ijk}$$ that takes value $$1$$ if $$k$$ is written in the cell $$(i,j)$$, or $$0$$ otherwise. We first ask that for each cell exactly one digit is selected:

(1)$\sum_{k=0}^{n-1} x_{ijk} = 1, \qquad i,j=0,\ldots,n-1.$

Similar constraints can be used to force each digit to appear only once in each row or column:

(2)$\begin{split}\begin{array}{lr} \sum_{i=0}^{n-1} x_{ijk} = 1, & j,k=0,\ldots,n-1,\\ \sum_{j=0}^{n-1} x_{ijk} = 1, & i,k=0,\ldots,n-1. \end{array}\end{split}$

To force a digit to appear only once in each sub-grid we can use the following

(3)$\sum_{i=0}^{m-1} \sum_{j=0}^{m-1} x_{(i+tm)(j+tl)k} = 1 \qquad k=0,\ldots,n-1 \text{ and } t,l=0,\ldots,m-1$

If a cell $$(i,j)$$ has a predetermined value, i.e. $$(i,j,k)\in F$$ then we set

$x_{ijk} = 1.$

Summarizing, and considering that there is no objective function to minimize, the optimization model for the SUDOKU problem takes the form

(4)$\begin{split}\begin{array}{lr} \min 0 &\\ \mbox{s.t.} &\\ \sum_{i=0}^{n-1} x_{ijk} = 1, & j,k=0,\ldots,n-1,\\ \sum_{j=0}^{n-1} x_{ijk} = 1, & i,k=0,\ldots,n-1,\\ \sum_{k=0}^{n-1} x_{ijk} = 1, & i,j=0,\ldots,n-1,\\ \sum_{i=0}^{m-1} \sum_{j=0}^{m-1} x_{(i+tm)(j+tl)k} = 1, & k=0,\ldots,n-1 \text{ and } \\ & t,l=0,\ldots,m-1,\\ &\\ x_{ijk} = 1, & \forall (i,j,k)\in F. \end{array}\end{split}$

Implementation with Fusion

The implementation in Fusion is straightforward. First, we represent the variable $$x$$ using a three dimensional Fusion variable:

x= M.variable([n,n,n], Domain.binary());


Then we can define constraints (1) and (2) simply using the Expr.sum operator, that allows to sum the elements of an expression (in this case of the variable itself) along arbitrary dimensions. The code reads:

%each value only once per dim
for d = 1:m
M.constraint(Expr.sum(x,d-1), Domain.equalsTo(1.));
end


The last set of constraints (3) , i.e. the sum over block, needs a little more effort: we must loop over all blocks and select the proper slice:

%each number must appears only once in a block
for k = 1:n
for i = 1:m
for j = 1:m
M.constraint( Expr.sum( x.slice([1+(i-1)*m,1+(j-1)*m,k], [1+i*m, 1+j*m, k+1]) ), ...
Domain.equalsTo(1.) );
end
end
end


To set the triplets given in the set $$F$$ we can use the Variable.pick method that returns a one dimensional view of an arbitrary set of elements of the variable.

M.constraint( x.pick(hr_fixed), Domain.equalsTo(1.0) );


SUDOKU: the complete example code.

The complete code for the SUDOKU problem is shown in Listing 10.17.

Listing 10.17 Fusion implementation to solve SUDOKU. Click here to download.
function sudoku()

import mosek.fusion.*;

m= 3;
n= m*m;

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

M= Model('SUDOKU');

x= M.variable([n,n,n], Domain.binary());

%each value only once per dim
for d = 1:m
M.constraint(Expr.sum(x,d-1), Domain.equalsTo(1.));
end

%each number must appears only once in a block
for k = 1:n
for i = 1:m
for j = 1:m
M.constraint( Expr.sum( x.slice([1+(i-1)*m,1+(j-1)*m,k], [1+i*m, 1+j*m, k+1]) ), ...
Domain.equalsTo(1.) );
end
end
end

M.constraint( x.pick(hr_fixed), Domain.equalsTo(1.0) );

M.setLogHandler(java.io.PrintWriter(java.lang.System.out));
M.solve();

%print the solution, if any...
if M.getPrimalSolutionStatus() == SolutionStatus.Optimal || ...
M.getPrimalSolutionStatus() == SolutionStatus.NearOptimal

fprintf('\n');
for i = 1:n
fprintf(' |');
for j = 1:n

for k = 1:n
if x.index([i,j,k]).level()>0.5
fprintf(' %d', k)
break;
end
end
if mod(j,m) == 0
fprintf(' |');
end

end
fprintf('\n');
if mod(i,m) == 0
fprintf('\n');
end

end
else
fprintf('No solution found!\n');
end

M.dispose()


The problem instance corresponding to Fig. 10.6 is hard-coded for the sake of simplicity. It will produce the following output

Problem
Name                   : SUDOKU
Objective sense        : min
Type                   : LO (linear optimization problem)
Constraints            : 350
Cones                  : 0
Scalar variables       : 1000
Matrix variables       : 0
Integer variables      : 729

Optimizer started.
Mixed integer optimizer started.
Presolve started.
Presolve terminated. Time = 0.00
Presolved problem: 0 variables, 0 constraints, 0 non-zeros
Presolved problem: 0 general integer, 0 binary, 0 continuous
Clique table size: 0
BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME
0        1        0        0        0.0000000000e+00     0.0000000000e+00     0.00e+00    0.0
An optimal solution satisfying the relative gap tolerance of 1.00e-02(%) has been located.
The relative gap is 0.00e+00(%).
An optimal solution satisfying the absolute gap tolerance of 0.00e+00 has been located.
The absolute gap is 0.00e+00.

Objective of best integer solution : 0.000000000000e+00
Best objective bound               : -0.000000000000e+00
Construct solution objective       : Not employed
Construct solution # roundings     : 0
User objective cut value           : 0
Number of cuts generated           : 0
Number of branches                 : 0
Number of relaxations solved       : 1
Number of interior point iterations: 0
Number of simplex iterations       : 0
Time spend presolving the root     : 0.00
Time spend in the heuristic        : 0.00
Time spend in the sub optimizers   : 0.00
Time spend optimizing the root   : 0.00
Mixed integer optimizer terminated. Time: 0.02

Optimizer terminated. Time: 0.02

| 3 2 9 | 5 4 7 | 6 1 8 |
| 6 5 8 | 9 1 3 | 4 2 7 |
| 4 1 7 | 2 8 6 | 9 5 3 |

| 9 7 3 | 1 5 2 | 8 4 6 |
| 5 6 2 | 8 7 4 | 3 9 1 |
| 8 4 1 | 6 3 9 | 2 7 5 |

| 1 9 4 | 3 6 5 | 7 8 2 |
| 7 3 5 | 4 2 8 | 1 6 9 |
| 2 8 6 | 7 9 1 | 5 3 4 |