Another Rook Tour on the Chessboard

In [1] the following problem is described:

1. The chessboard shown above has \(8\times 8 = 64\) cells.
2. Place numbers 1 through 64 in the cells such that the next number is either directly left, right, below, or above the current number. 
3. The dark cells must be occupied by prime numbers. Note there are 18 dark cells, and 18 primes in the sequence \(1,\dots,64\). These prime numbers are: 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61.

This problem looks very much like the Numbrix puzzle [3]. We don't have the given cells that are provided in Numbrix, but instead, we have some cells that need to be covered by prime numbers.

The path formed by the consecutive values looks like how a rook moves on a chessboard. Which explains the title of this post (and of the original post). If we require a proper tour we need values 1 and 64 to be neighbors. Otherwise, if this is not required, I think a better term is: Hamiltonian path. I'll discuss both cases, starting with the Hamiltonian path problem.

Mathematical Model

As usual for this type of problems, we use the following decision variables: \[\color{darkred}x_{i,j,k} = \begin{cases} 1 & \text{if cell $(i,j)$ has a value $k$} \\ 0 & \text{otherwise}\end{cases}\] With this, our model can look like: 

Hamiltonian Path Model

\[\begin{align} & \sum_k \color{darkred}x_{i,j,k} = 1 && \forall i,j &&\text{One value in a cell}\\ & \sum_{i,j} \color{darkred}x_{i,j,k} = 1 && \forall k &&\text{Unique values} \\ & \color{darkred}x_{i-1,j,k+1}+\color{darkred}x_{i+1,j,k+1}+\color{darkred}x_{i,j-1,k+1}+\color{darkred}x_{i,j+1,k+1} \ge \color{darkred}x_{i,j,k} && \forall i,j,k\lt 64 && \text{Next value is a neighbor} \\ & \sum_p \color{darkred}x_{i,j,p} = 1 && \text{for dark cells $(i,j)$} &&\text{Prime values} \\ & \color{darkred}x_{i,j,k} \in \{0,1\} \end{align} \]

Here \(p\) is a subset of \(k\): all the prime values. The next-value constraint is explained in more detail in [3]. This is a fascinating constraint. The only thing we are saying here is: the next value is found among the neighbors of \((i,j)\). Or to be more precise: \[\color{darkred}x_{i,j,k}=1\Rightarrow \color{darkred}x_{i-1,j,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i+1,j,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i,j-1,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i,j+1,k+1}=1\]There is no concept of a path here. Rather, the formation of a path is a natural consequence. I assume in this constraint that when we address outside the board, the corresponding variable is zero. I think this model is quite elegant. It precisely states the conditions we must obey for a valid solution. 

Notes:

• This model has no objective.
• We don't require a tour here. There is no restriction that says that values 1 and 64 must be neighbors. That version of the problem is looked at further down.
• This is not a very small model. The number of \(\color{darkred}x_{i,j,k}\) variables is \(8 \times 8 \times 64 = 4096\). The model size (after adding a dummy objective) is:

MODEL STATISTICS

BLOCKS OF EQUATIONS           6     SINGLE EQUATIONS        4,179
BLOCKS OF VARIABLES           2     SINGLE VARIABLES        4,097
NON ZERO ELEMENTS        26,661     DISCRETE VARIABLES      4,096

A final question is: is the solution unique or are there multiple solutions? To forbid a previously found solution, we can add the cut: \[\sum_{i,j,k} \color{darkblue}x^*_{i,j,k} \cdot \color{darkred}x_{i,j,k} \le  {\bf{card}}(k)-1\] where \(\color{darkblue}x^*_{i,j,k}\) is such a previously found solution. This constraint says: at least one cell should change. Note that we never can only change a 0 value to a 1: the total number of ones is given by \({\bf{card}}(k)=64\). So we only need to keep track of the ones. And that is what this constraint is doing: force to flip at least one 1 to a 0. We repeat solving and adding a cut to the model until the model becomes infeasible. An alternative approach would be to use the solution pool to multiple optimal solutions. The solution pool is a technology found in some advanced solvers. As the number of solutions is small, I use here the DIY approach using the "no good" cuts.

The feasible paths we find are:

----    106 PARAMETER result  negative numbers indicate prime cells

INDEX 1 = solution1

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          51          50         -47          -2          -5           8           9          16
i6          52          49          48           1           6          -7          10          15
i7         -53          56          57          58         -59          60         -11          14
i8          54          55          64          63          62         -61          12         -13

INDEX 1 = solution2

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          51          50         -47          -2          -5           8           9          16
i6          52          49          48           1           6          -7          10          15
i7         -53          64          63          62         -61          60         -11          14
i8          54          55          56          57          58         -59          12         -13

INDEX 1 = solution3

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          49          48         -47          -2          -5           8           9          16
i6          50          51          64           1           6          -7          10          15
i7         -53          52          63          62         -61          60         -11          14
i8          54          55          56          57          58         -59          12         -13

INDEX 1 = solution4

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -43          44          45           6          -7           8           9          10
i2          42         -47          46          -5          64          25          24         -11
i3         -41          48          49           4          63          26         -23          12
i4          40          51          50          -3          62          27          22         -13
i5          39          52         -53          -2         -61          28          21          14
i6          38          55          54           1          60         -29          20          15
i7         -37          56          57          58         -59          30         -19          16
i8          36          35          34          33          32         -31          18         -17

Or plotted using R/ggplot:

All feasible paths

Finding feasible tours

Looking at these solutions, we see that just one solution forms a proper tour. The obvious question is: can we add a constraint to the model that just delivers this tour. Actually, this is rather simple. We need to change the equation next. The mathematical formulation looks like: \[\begin{cases}\color{darkred}x_{i-1,j,k+1}+\color{darkred}x_{i+1,j,k+1}+\color{darkred}x_{i,j-1,k+1}+\color{darkred}x_{i,j+1,k+1} \ge \color{darkred}x_{i,j,k} & \forall i,j,k\lt 64 \\ \color{darkred}x_{i-1,j,1}+\color{darkred}x_{i+1,j,1}+\color{darkred}x_{i,j-1,1}+\color{darkred}x_{i,j+1,1} \ge \color{darkred}x_{i,j,64} & \forall i,j\end{cases}\]

Here is how we can do this a bit more easily in GAMS. The original version had (see the appendix below for the complete model):

next(i,j,k+1)..   
   x(i-1,j,k+1)+x(i+1,j,k+1)+x(i,j-1,k+1)+x(i,j+1,k+1) =g= x(i,j,k);

This version has k+1 as equation index. That means we generate this equation for all but the last k. If we want to require that value 1 and 64 are neighbors, we can simply write:

next(i,j,k)..   
   x(i-1,j,k++1)+x(i+1,j,k++1)+x(i,j-1,k++1)+x(i,j+1,k++1) =g= x(i,j,k);

Here the ++ operator is used to implement a circular lead. Also, note that the equation next is now indexed by all k. The net effect is that we add 64 rows to the model (one for each cell). If we run the model with this equation we get a single solution:

----    106 PARAMETER result  negative numbers indicate prime cells

INDEX 1 = solution1

            j1          j2          j3          j4          j5          j6          j7          j8

i1         -41          40          39          30         -29          28          21          20
i2          42         -37          38         -31          32          27          22         -19
i3         -43          36          35          34          33          26         -23          18
i4          44          45          46          -3           4          25          24         -17
i5          49          48         -47          -2          -5           8           9          16
i6          50          51          64           1           6          -7          10          15
i7         -53          52          63          62         -61          60         -11          14
i8          54          55          56          57          58         -59          12         -13

This is the third path from our earlier picture.

Conclusions

The problem of finding a Hamiltonian path while obeying the prime-value-only cells, is a fairly close variant of the Numbrix model discussed in [3]. Finding a tour can be handled by changing a single equation (but in an interesting way).

This was an intriguing modeling exercise. Both the underlying mathematical model and the implementation are not very complicated but somewhat surprising. Some puzzles can be elegantly stated as mathematical models and simply solved with off-the-shelf solvers. Probably this was not intended, but I often find doing these modeling exercises just as rewarding as solving these puzzles by hand. It forces you to think about the problem in a different way, totally different from solving it by hand.

References

1. Another Rook's Tour of the Chessboard, https://puzzling.stackexchange.com/questions/111818/another-rooks-tour-of-the-chessboard
2. Chessboard, https://en.wikipedia.org/wiki/Chessboard
3. Numbrix puzzle via MIP, https://yetanothermathprogrammingconsultant.blogspot.com/2020/12/numbrix-puzzle-via-mip.html 
4. Prime numbers in GAMS, https://www.gams.com/34/gamslib_ml/libhtml/gamslib_prime.html. Unfortunately, it skips the prime number 2.

Appendix: GAMS model

$ontext    A Numbrix-like puzzle.    Place the numbers 1..64 on an 8x8 board such that:       1. next values are neighbors. A neighbor is defined as          below, left, right, or above.       2. certain cells need to contain a prime number    https://puzzling.stackexchange.com/questions/111818/another-rooks-tour-of-the-chessboard $offtext *---------------------------------------------------------------------- * data *---------------------------------------------------------------------- sets   i 'rows'    /i1*i8/   j 'columns' /j1*j8/   k 'values'  /1*64/   kprime(k) 'primes' /              2,3,5,7,11,13,17,19,23,29,              31,37,41,43,47,53,59,61     / ; acronym p; table board(i,j) 'locations with a p should get a prime value'         j1 j2 j3 j4 j5 j6 j7 j8   i1     p           p   i2        p     p           p   i3     p                 p   i4              p           p   i5           p  p  p   i6                    p   i7     p           p     p   i8                    p     p ; *---------------------------------------------------------------------- * model *---------------------------------------------------------------------- set sol /solution1*solution10/; set solfound(sol) 'solutions found'; parameter solution(sol,i,j,k) 'solutions found so far'; solfound(sol) = no; solution(sol,i,j,k) = 0; binary variable x(i,j,k); variable dummy; equations    onevalue(i,j) 'each cell has one value (between 1 and 81)'    unique(k)     'each value must be in exactly one cell'    next(i,j,k)   'x(i,j,k)=1 => a neighbor is k+1'    isprime(i,j)  'we need a prime on (i,j)'    cut(sol)      'cut off previous solutions'    obj           'dummy' ; obj.. dummy =e= 0; onevalue(i,j).. sum(k, x(i,j,k)) =e= 1; unique(k)..     sum((i,j),x(i,j,k)) =e= 1; next(i,j,k+1)..     x(i-1,j,k+1)+x(i+1,j,k+1)+x(i,j-1,k+1)+x(i,j+1,k+1) =g= x(i,j,k); isprime(i,j)$(board(i,j)=p)..     sum(kprime, x(i,j,kprime)) =e= 1; cut(solfound)..     sum((i,j,k),solution(solfound,i,j,k)*x(i,j,k)) =l= card(k)-1; model m /all/; option mip=cplex,threads=8; *---------------------------------------------------------------------- * solve loop *---------------------------------------------------------------------- loop(sol,     solve m minimizing dummy using mip;     break$(m.modelstat <> %modelStat.optimal%            and m.modelstat <> %modelStat.integerSolution%);     solfound(sol) = yes;     solution(sol,i,j,k) = round(x.l(i,j,k));     display x.l; ); *---------------------------------------------------------------------- * reporting *---------------------------------------------------------------------- parameter result(sol,i,j) 'negative numbers indicate prime cells'; result(sol,i,j) = sum(k, k.val*solution(sol,i,j,k)); result(sol,i,j)$(board(i,j)=p) = -result(sol,i,j); option result:0:1:1; display result;

Notes:

1. An acronym allows a string to be used in a numerical table (like inf, eps etc.).
2. The equation next is indexed by \(k+1\). This means the equation is generated for all but the last k.
3. If \(i-1\),\(i+1\),\(j-1\),\(j+1\) is outside the domain (i.e. outside the board), GAMS will return zero for the variable (or in other words: ignore that entry). That is exactly what we need here. So we don't need to add additional protection to remain inside the board.
4. This version of the model generates all feasible Hamiltonian paths.
5. We could have reduced the number of variables a bit by only generating valid \(\color{darkred}x_{i,j,k}\): for prime cells, non-prime \(k\)'s are not allowed. That would also eliminate the isprime equation.