11.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:

_images/sudoku-ex.png

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:

(11.27)\[\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:

(11.28)\[\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

(11.29)\[\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

(11.30)\[\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:

        Variable x = M.Variable(new int[] { n, n, n }, Domain.Binary());

Then we can define constraints (11.27) and (11.28) 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 dimension
        for (int d = 0; d < m; d++)
          M.Constraint(Expr.Sum(x, d), Domain.EqualsTo(1.0));

The last set of constraints (11.29) , 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 (int k = 0; k < n; k++)
          for (int i = 0; i < m; i++)
            for (int j = 0; j < m; j++)
            {
              Variable block = x.Slice(new int[] { i * m, j * m, k },
                                       new int[] { (i + 1)*m, (j + 1)*m, k + 1 });
              M.Constraint(Expr.Sum(block), Domain.EqualsTo(1.0));
            }

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(fixed_cells), Domain.EqualsTo(1.0));

SUDOKU: the complete example code.

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

Listing 11.16 Fusion implementation to solve SUDOKU. Click here to download.
namespace mosek.fusion.example
{
  public class sudoku
  {
    public static void Main(string[] args)
    {
      int m = 3;
      int n = m * m;

      //fixed cells in human readable (i.e. 1-based) format
      int[,] 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}
      };
      int nf = hr_fixed.Length / 3;
      int[,] fixed_cells = new int[nf, 3];
      for (int i = 0; i < nf; i++)
        for (int d = 0; d < m; d++)
          fixed_cells[i, d] = hr_fixed[i, d] - 1;

      using (Model M = new Model("SUDOKU")) {
        M.SetLogHandler(Console.Out);
        Variable x = M.Variable(new int[] { n, n, n }, Domain.Binary());

        //each value only once per dimension
        for (int d = 0; d < m; d++)
          M.Constraint(Expr.Sum(x, d), Domain.EqualsTo(1.0));

        //each number must appears only once in a block
        for (int k = 0; k < n; k++)
          for (int i = 0; i < m; i++)
            for (int j = 0; j < m; j++)
            {
              Variable block = x.Slice(new int[] { i * m, j * m, k },
                                       new int[] { (i + 1)*m, (j + 1)*m, k + 1 });
              M.Constraint(Expr.Sum(block), Domain.EqualsTo(1.0));
            }

        M.Constraint(x.Pick(fixed_cells), Domain.EqualsTo(1.0));

        M.Solve();

        //print the solution, if any...
        if (M.GetPrimalSolutionStatus() == SolutionStatus.Optimal)
          print_solution(m, x);
        else
          Console.WriteLine("No solution found!\n");
      }
    }

    public static void print_solution(int m, Variable x)
    {
      int n = m * m;
      Console.WriteLine("\n");
      try
      {
        for (int i = 0; i < n; i++)
        {
          String s = "";
          for (int j = 0; j < n; j++)
          {
            if (j % m == 0) s += " |";
            for (int k = 0; k < n; k++)
            {
              double[] sol = x.Index(new int[] { i, j, k }).Level();
              if (sol[0] > 0.5)
              {
                s += " " + (k + 1);
                break;
              }
            }
          }
          Console.WriteLine(s + " |");
          if ((i + 1) % m == 0)
            Console.WriteLine("\n");
        }
      }
      catch (SolutionError se) { }
    }
  }
}

The problem instance corresponding to Fig. 11.8 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.
Threads used: 2
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 |