Showing posts with label Sudoku. Show all posts
Showing posts with label Sudoku. Show all posts

Monday, May 25, 2020

The Miracle Sudoku

Introduction


In [1] another strange Sudoku variant is presented.



The givens are just:



Obviously, more than just the standard Sudoku constraints are needed to have just one unique feasible solution. So, for this problem we have the following rules:


  1. Standard Sudoku rules apply.  So, for each row, column, and sub-block (a.k.a. nonet [3]), we have uniqueness constraints (each cell in a region has a unique number between 1 and 9).
  2. Borrowing from chess: any two cells corresponding to a knight's move, must have different values.
  3. Similarly, any two cells that form a king's move, must also be different.
  4. Orthogonally adjacent cells must be non-consecutive.

Standard Sudoku setup


When we want to solve Sudokus, the easiest approach is to define the following binary decision variables [4]: \[x_{i,j,k} = \begin{cases} 1 & \text{if cell $(i,j)$ has value $k$} \\ 0 & \text{otherwise}\end{cases}\] Here \(k \in \{1,\dots,9\}\). We have 27 areas we need to check for unique values: rows, columns and nonets. We can organize this as a set:  \[u_{a,i,j}\>\text{exists if and only if area $a$ contains cell $(i,j)$}\] This is data. We also have these two given cells, i.e. fixed variables for these cells. The resulting MIP model can look like [4]:

Mixed-Integer Programming Model for standard Sudoku
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\ &\color{darkred}x_{i,j,k} \in \{0,1\} \end{align}\]

During post-processing, we can calculate the completed grid using the optimal values of \(x\): \[v_{i,j} = \sum_k k \cdot x^*_{i,j,k}\] or we can include them as "accounting rows" to the model. Accounting rows are constraints that have the function to calculate some quantities for reporting purposes.

The implementation in GAMS uses a set u(a,i,j)that is populated as:
set
  a 'areas' /r1*r9,c1*c9,s1*s9/
  i(a) 'rows' /r1*r9/
  j(a) 'columns' /c1*c9/
  s(a) 'squares' /s1*s9/
  u(a,i,j) 'areas with unique values'
;

* populate u(a,i,j): u(a,i,j)=YES if area a contains cell (i,j)
* see https://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html
u(i,i,j) = yes;
u(j,i,j) = yes;
u(s,i,j)$(ord(s)=3*ceil(ord(i)/3)+ceil(ord(j)/3)-3) = yes;

The model equations look like:

binary variable x(i,j,k);
* v0(i,j) are the givens
x.fx(i,j,k)$(v0(i,j)=ord(k)) = 1;
variable z 'dummy objective';
integer variable v(i,j);
v.lo(i,j) = 1;
v.up(i,j) = 9;


equation
   dummy 'objective'
   unique(a,k) 'all-different'
   values(i,j) 'one value per cell'
   eqv(i,j)  'calculate values'
;

dummy.. z =e= 0;
unique(a,k).. sum(u(a,i,j), x(i,j,k)) =e= 1;
values(i,j).. sum(k, x(i,j,k)) =e= 1;
eqv(i,j).. v(i,j) =e= sum(k, ord(k)*x(i,j,k));

model sudoku /all/;
solve sudoku minimizing z using mip;

We can test this approach with the above data. We will actually find a feasible solution, but of course it will not be unique. I believe there are likely millions of different solutions (or more). Counting them is not that easy. I tried to give it a shot with a solution pool approach: after setting a limit if 1,000,000 solutions, the solver stopped after reaching this limit. So, there are more than 1 million solutions for this (partial) problem.

Knight's jumps


To model how knight's moves are covering cells, we can have a look at [5].

To place knights we set up a set \(\mathit{jump}_{i,j,i',j'}\) indicating if we can jump from cell \((i,j)\) to cell \((i',j')\). If we can jump from \((i,j) \rightarrow (i',j')\) we can also jump from \((i',j') \rightarrow (i,j)\). We don't want to check both cases. To prevent this double check, we only need to look forward. So for each \((i,j)\) we need to consider just four cases:
\(j\)\(j+1\)\(j+2\)
\(i-2\)\(x_{i-2,j+1}\)
\(i-1\)\(x_{i-1,j+2}\)
\(i\)\(x_{i,j}\)
\(i+1\)\(x_{i+1,j+2}\)
\(i+2\)\(x_{i+2,j+1}\)

Note that near the border we may have fewer than four cases. In GAMS we can populate the set \(\mathit{jump}\) in a straightforward manner:

alias (i,ii),(j,jj);
set jump(i,j,ii,jj);

jump(i,j,i-2,j+1) = yes;
Jump(i,j,i-1,j+2) = yes;
Jump(i,j,i+1,j+2) = yes;
Jump(i,j,i+2,j+1) = yes;


There are 224 elements in the set \(\mathit{jump}\).

Now comes the more complicated part: how to make sure that the values of cell \((i,j)\) and \((i',j')\) are not the same. There are basically two approaches:

  • Compare the values \(v_{i,j}\) and \(v_{i',j'}\). The constraint becomes \[|v_{i,j}-v_{i',j'}| \ge 1 \>\> \forall i,j,i',j'|\mathit{jump}(i,j,i',j')\] This can be linearized using \[v_{i,j}-v_{i',j'} \ge 1 \textbf{ or } v_{i',j'}-v_{i,j} \ge 1\] or \[\begin{align}& v_{i,j}-v_{i',j'} \ge 1 - M\cdot \delta_{i,j,i',j'} \\ & v_{i',j'}-v_{i,j} \ge 1- M (1-\delta_{i,j,i',j'}) \\ &\delta_{i,j,i',j'} \in \{0,1\}\end{align}\] Here \(M\) is a large enough constant (\(M=10\) should suffice).
  • Directly work with \(x_{i,j,k}\) and \(x_{i',j',k}\). The constraint should be \[x_{i,j,k} x_{i',j',k} = 0 \>\> \forall k,  \forall i,j,i',j'|\mathit{jump}(i,j,i',j') \] This non-linear constraint can be linearized as \[x_{i,j,k}+x_{i',j',k}\le 1\] This clearly is the easier route. 

In GAMS this can look like:

equation Different(i,j,ii,jj,k) 'cells forming a knights move should be different';

Different(jump(i,j,ii,jj),k)..
     x(i,j,k) + x(ii,jj,k) =l= 1;


Kings's moves


Here we want to implement a similar restriction as for the knight's move. As the constraints are the same, we just have to augment the set  \(\mathit{jump}\) a bit.


alias (i,ii),(j,jj);
set jump(i,j,ii,jj) 'knights or kings move';

* knight
jump(i,j,i-2,j+1) = yes;
jump(i,j,i-1,j+2) = yes;
jump(i,j,i+1,j+2) = yes;
jump(i,j,i+2,j+1) = yes;
* king
jump(i,j,i+1,j+1) = yes;
jump(i,j,i-1,j+1) = yes;


After this, the set \(\mathit{jump}\)  has grown from 224 to 352 elements.

Orthogonal neighbors can't be consecutive


First we need to build a set that models "orthogonal neighbors".  Similar to our set \(\mathit{jump}\), we want to prevent comparing pairs twice. So given a cell \((i,j)\), we only need to consider two neighbors.
\(j\)\(j+1\)
\(i-1\)\(x_{i-1,j}\)
\(i\)\(x_{i,j}\)\(x_{i,j+1}\)
\(i+1\)


This is coded in GAMS as:

set nb(i,j,ii,jj) 'orthogonal neighbors';
nb(i,j,i-1,j) = yes;
nb(i,j,i,j+1) = yes;

This set has 144 elements.

The constraints are: \[\begin{align}&x_{i,j,k} + x_{i',j',k+1} \le 1\\&x_{i,j,k} + x_{i',j',k-1} \le 1\end{align}\>\> \forall i,j,i',j'|nb(i,j,i',j'), \forall k\]

The GAMS representation is:

equations
   NonConsecutive1(i,j,ii,jj,k) 'orthogonal neighbors constraint'
   NonConsecutive2(i,j,ii,jj,k) 'orthogonal neighbors constraint'
;

NonConsecutive1(nb(i,j,ii,jj),k+1)..
     x(i,j,k) + x(ii,jj,k+1) =l= 1;
NonConsecutive2(nb(i,j,ii,jj),k-1)..
     x(i,j,k) + x(ii,jj,k-1) =l= 1;

Detail: The last index of the equation definition has a \(k+1\) or \(k-1\). This tells GAMS not to generate all possile constraints for all \(k\) but rather one less. This is not strictly needed: in GAMS addressing outside the domain results in a zero. But this trick will prevent singleton constraints \(x_{i,j,k} + 0 \le 1\) to be generated. 

Results


The complete model has 811 variables (808 of them binary or integer). The number of constraints is 5,878. The results look like:


----    111 VARIABLE v.L  

            c1          c2          c3          c4          c5          c6          c7          c8          c9

r1           4           8           3           7           2           6           1           5           9
r2           7           2           6           1           5           9           4           8           3
r3           1           5           9           4           8           3           7           2           6
r4           8           3           7           2           6           1           5           9           4
r5           2           6           1           5           9           4           8           3           7
r6           5           9           4           8           3           7           2           6           1
r7           3           7           2           6           1           5           9           4           8
r8           6           1           5           9           4           8           3           7           2
r9           9           4           8           3           7           2           6           1           5


The problem solves in the presolve phase, i.e. 0 iterations and 0 nodes.

Uniqueness of the solution


To prove the uniqueness of the solution, we add a cut that forbids the current solution \(x^*\): \[\sum_{i,j,k} x^*_{i,j,k} x_{i,j,k} \le 9^2-1\] The resulting model is infeasible. This means that the solution is indeed unique.


Conclusion


This Sudoku variant combines standard Sudoku rules with additional constraints, some of them borrowed from chess. This makes the modeling at bit trickier. In the development of the model, I have placed much of the logic in sets. This has the advantage that the constraints are fairly simple. In general, sets are easier to debug than constraints: we can display and debug sets before we have a working model.

Both the mathematical model and the GAMS representation is interesting. The mathematical formulation shows we can model this without extra (discrete) variables. The GAMS implementation details assembling sets that simplify constraints.

Sudoku models typically solve extremely fast: they are solved by the presolver. This model is no exception.

References


Posted by at No comments:
Labels: ,

Thursday, March 26, 2020

Special Sudoku (2)

In [1], a variant of Sudoku was discussed, where a side constraint is:

Orthogonally neigboring cells differ in value by at most 5

This max 5 difference is a convex constraint, so it can be expressed quite easily: \[-5 \le v_{i,j} - v_{i',j'} \le 5\] for neighboring cells \((i,j)\) and \((i',j')\).

A more difficult restriction is: neighboring cells must differ by at least 2. Here is such a problem, again designed by Aad van de Wetering:


Problem data

The rules are:
  • Standard Sudoku-rules apply. cells in each row, column and sub-block have unique values from 1 through 9.
  • In addition, neighboring (horizontal and vertical) cells have values that differ by 2 or more.
  • Also, values 1 and 9 are considered to have a difference of 1. So neighboring cells, cannot have values 1 and 9.

I will try to solve this problem in two different ways. The first is an extension of what we did in [1]. The second one is a bit different.

Approach I


The constraint: neighbors must be different by at least two, can be formalized as \[\left| v_{i,j} - v_{i',j'} \right| \ge 2\] for neighboring cells \((i,j)\) and \((i',j')\). The non-convexity of this constraint  \[\begin{align} & v_{i,j} - v_{i',j'} \ge 2 \\ & \mathbf{or} \\&v_{i,'j'} - v_{i,j} \ge 2 \end{align} \] requires extra binary variables and big-M constraints: \[\begin{align} & v_{i,j} - v_{i',j'} \ge 2 - M\cdot \delta_{i,j,i',j'} \\ & v_{i,'j'} - v_{i,j} \ge 2 - M\cdot (1-\delta_{i,j,i',j'})\\ &  \delta_{i,j,i',j'} \in \{0,1\}\end{align}\] So, how large should \(M\) be? The maximum difference between two different values is 8. So we need \(2-M\le -8\) or \(M=10\).

We ignored that values 1 and 9 also have a difference of 1. As all other combinations are allowed, we can handle this by \[\left|v_{i,j}-v_{i',j'}\right| \le 7\] This is a constraint we already know how to handle: \[-7 \le v_{i,j}-v_{i',j'} \le 7\] Interestingly, this addition restriction also can help us in tightening our big-M value. We can now use \(M = 9\). Note that with this smaller big-M we still need this constraint: the big-M constraint only covers one of the needed two sides of this constraint.

The complete model can look like:


Mixed-Integer Programming Model I
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ &\color{darkred}v_{i,j} =\sum_k \color{darkblue}k\cdot\color{darkred}x_{i,j,k} && \forall i,j &&\text{Value of cell}\\ & \color{darkred}v_{i,j} - \color{darkred}v_{i',j'} \ge 2 - \color{darkblue}M\cdot \color{darkred}\delta_{i,j,i',j'}&&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} &&\text{Min difference of 2} \\ & \color{darkred}v_{i,'j'} - \color{darkred}v_{i,j} \ge 2 - \color{darkblue}M\cdot (1-\color{darkred}\delta_{i,j,i',j'}) &&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'}\\ & -7 \le \color{darkred}v_{i,j} - \color{darkred}v_{i',j'} \le 7 &&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} && \text{Neighbors differ by max 7} \\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\&\color{darkred}x_{i,j,k} \in \{0,1\} \\ &\color{darkred}v_{i,j} \in \{1,\dots,9\} \\ & \color{darkred}\delta_{i,j,i',j'} \in \{0,1\} \\ & \color{darkblue}M = 9\end{align}\]


Notes:

  • This model solves quickly: the presolve can remove all rows and columns.
  • The solution is unique (see [1]).
  • Using the above data set we actually can drop the "max 7" constraint. I think that is just a property of this data set. In general, we should include this constraint.


Approach II


A simpler approach is to consider the pairs of values that cannot be neighbors: \[(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,1)\] We also need to add their symmetrical counterparts \[(2,1),(3,2),\dots\] Detail: this is due to the way we organized the neighbors. The neighbor set did not have both: i.e. cells \(\text{$(1,2)$-$(1,3)$}\) are neighbors, so this is stored. But, I did not also store the combination \(\text{$(1,3)$-$(1,2)$}\).

After automatically generating the combinations \((k,k')\) that we need to prevent, the forbid data can look like:



----     55 SET forbid  combinations of values we want to exclude

            v1          v2          v3          v4          v5          v6          v7          v8          v9

v1                     YES                                                                                 YES
v2         YES                     YES
v3                     YES                     YES
v4                                 YES                     YES
v5                                             YES                     YES
v6                                                         YES                     YES
v7                                                                     YES                     YES
v8                                                                                 YES                     YES
v9         YES                                                                                 YES


This table contains elements \(\mathit{forbid}_{k,k'}\) we want to forbid. Note that we can exclude the diagonal from this: neighboring cells already cannot have the same value. For this set of value-pairs, we want to add the constraint:\[x_{i,j,k} + x_{i',j',k'} \le 1 \] for all neighboring cells \((i,j)\) and \((i',j')\) and for all forbidden value combinations \((k,k')\). We have 18 of these combinations and 144 neighbors to check, yielding 2592 of these constraints.

The full model can look like:

Mixed-Integer Programming Model II
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ & \color{darkred}x_{i,j,k} + \color{darkred}x_{i',j',k'} \le 1 && \begin{aligned}& \forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} \\ & \forall k,k'|\color{darkblue}{\mathit{forbid}}_{k,k'} \end{aligned}&&\text{Forbid combinations} \\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\&\color{darkred}x_{i,j,k} \in \{0,1\} \end{align}\]


When we run this model, we will see the solution:


----     84 PARAMETER v  values (postprocessing)

            c1          c2          c3          c4          c5          c6          c7          c8          c9

r1           9           2           6           4           1           7           3           5           8
r2           3           5           8           2           6           9           7           1           4
r3           1           7           4           8           3           5           2           6           9
r4           7           4           1           6           9           2           5           8           3
r5           5           8           3           1           7           4           9           2           6
r6           2           6           9           3           5           8           4           7           1
r7           6           9           2           5           8           3           1           4           7
r8           8           3           5           7           4           1           6           9           2
r9           4           1           7           9           2           6           8           3           5

This solution is unique. The presolver completely solves the model, so no actual iterations are performed.

References



  1. Special Sudoku, http://yetanothermathprogrammingconsultant.blogspot.com/2020/03/special-sudoku.html
Posted by at No comments:
Labels: , , ,

Sunday, March 22, 2020

Special Sudoku

Here is a difficult Sudoku variant designed by Aad van de Wetering. Aad is a former colleague of mine.

Problem
The standard Sudoku rules apply:

  • Each cell has an integer value between 1 and 9 
  • In each row, column and sub-block cells must be unique
  • Some cells have a predefined value

In addition, we have the following restriction:

  • The maximum difference between the value of two neighboring cells (orthogonal, that is horizontal or vertical) is 5

To solve this you can follow this video [1]:



Of course, we want to try to solve this as a Mixed-Integer Programming model.

Review: a linear integer programming model for standard Sudoku


When we want to solve Sudokus, the easiest approach is to define the following binary decision variables [2]: \[x_{i,j,k} = \begin{cases} 1 & \text{if cell $(i,j)$ has value $k$} \\ 0 & \text{otherwise}\end{cases}\] Here \(k \in \{1,\dots,9\}\). We have 27 areas we need to check for unique values: rows, columns and sub-blocks. We can organize this as a set:  \[u_{a,i,j}\>\text{exists if and only if area $a$ contains cell $(i,j)$}\] This is data. We also have a set of given cells, which we can fix. The resultimg MIP model can look like [2]:

Mixed-Integer Programming Model for standard Sudoku
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\ &\color{darkred}x_{i,j,k} \in \{0,1\} \end{align}\]

During post-processing, we can calculate the completed grid using the optimal values of \(x\): \[v_{i,j} = \sum_k k \cdot x^*_{i,j,k}\]

Linear integer programming model for Sudoku variant


We need to introduce the concept of "neighbors". We can define the set: \[nb_{i,j,i',j'} \> \text{exists iff cells $(i,j)$ and $(i',j')$ are neighbors}\] I excluded symmetric cases. To be precise: in GAMS, I populated this set with:

alias (i,ii),(j,jj);
set
  nb(i,j,ii,jj) 'neighbors'
;

nb(i,j,i-1,j) = yes;
nb(i,j,i,j-1) = yes;


We can introduce the post-processing step directly into the MIP and then enforce our difference constraint on the variables \(v_{i,j}\). This can look like \[\begin{align} & v_{i,j} = \sum_k k \cdot x_{i,j,k} && \text{cell values as constraint} \\ & -5 \le v_{i,j} - v_{i',j'} \le 5 && \text{for neigboring cells $(i,j)$ and $(i',j')$}\end{align}\]
With this, we can formulate:


Mixed-Integer Programming Model for Sudoku Variant
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ &\color{darkred}v_{i,j} =\sum_k \color{darkblue}k\cdot\color{darkred}x_{i,j,k} && \forall i,j &&\text{Value of cell}\\ & -5 \le \color{darkred}v_{i,j} - \color{darkred}v_{i',j'} \le 5 &&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} && \text{Neighbors differ by max 5} \\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\&\color{darkred}x_{i,j,k} \in \{0,1\} \\ &\color{darkred}v_{i,j} \in \{1,\dots,9\} \end{align}\]

This model solves very easily. The presolver can solve this model completely, so we don't even have to start the branch-and-bound phase. As a result, we can make a much shorter movie:




The solution of this model looks like:


----     80 VARIABLE v.L  value of each cell

            c1          c2          c3          c4          c5          c6          c7          c8          c9

r1           2           4           9           6           7           5           8           3           1
r2           3           8           5           1           4           9           6           7           2
r3           7           6           1           2           3           8           9           5           4
r4           8           9           6           3           1           4           5           2           7
r5           5           7           3           8           6           2           1           4           9
r6           1           2           4           5           9           7           3           8           6
r7           6           1           2           4           8           3           7           9           5
r8           4           3           7           9           5           1           2           6           8
r9           9           5           8           7           2           6           4           1           3


Unique solution


A proper Sudoku puzzle has one unique solution. We can verify this by adding a cut to the problem after solving it. The cut should forbid the current solution (and allow all other solutions). In this case, the cut can look like \[\sum_{i,j,k} x^*_{i,j,k} x_{i,j,k}\le 9^2-1\] where \(x^*\) is the previously found solution. Indeed, this puzzle has a single, unique solution: after adding the cut, the problem becomes (integer) infeasible.

References


  1. A Sudoku With Just 11 Digits = Dutch Genius!, https://www.youtube.com/watch?v=ZU5fSDHJq8k
  2. MIP Modeling: from Sudoku to KenKen via Logarithms, https://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html

Posted by at No comments:
Labels: , , ,
Subscribe to: Posts (Atom)