Numbrix puzzle via MIP

The puzzle of numbrix [1,2] is simple to state.

• We have a 9x9 board with 81 cells. 
• We need to fill each cell with a unique integer value between 1 and 81.
• There are some prefilled cells.
• We should be able to start with the cell containing value 1 and then following a path (only going one up, down, left or right) to the next values and thus eventually arriving at the cell with value 81.  

Example

Below is an example of an input (the "given" cells) and a solution. The colors may help to see that we have a path from the cell with value 1 to the cell with value 81.

Of course, we want to see if we can create a MIP model for this.

Like in Sudoku models we use a binary variable \(\color{darkred}x_{i,j,k} \in \{0,1\}\) defined by \[\color{darkred}x_{i,j,k} = \begin{cases} 1 & \text{if cell $(i,j)$ has value $k$} \\ 0 & \text{otherwise}\end{cases}\] This introduces a large number of binary variables (\(9 \times 9 \times 81=6,561\)), but it makes the problem much simpler when formulating as a MIP model. This is really the big decision: once you have the main variables defined, things follow from there.

We need to formulate just three constraints:

1. Each cell can contain just one value. \[\sum_k \color{darkred}x_{i,j,k} = 1 \>\> \forall i,j\]
2. Each value should occur exactly once. \[\sum_{i,j}\color{darkred}x_{i,j,k} = 1 \>\> \forall k\]
3. We need: if \(\color{darkred}x_{i,j,k} = 1\) then one of its neighbors should contain value \(k+1\). A neighbor is defined as going one up, down, left or right. This can be formulated as an implication \[\color{darkred}x_{i,j,k} = 1 \Rightarrow  \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} = 1\] Or, as inequality: \[\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}\]

This is all! 

Note that I don't explicitly form a path in the model. Just the constraint about "next k should be a neighbor" is enough. So, our complete model looks like:

MIP Model

\[\begin{align} & \sum_k \color{darkred}x_{i,j,k} = 1 && \forall i,j \\ & \sum_{i,j}\color{darkred}x_{i,j,k} = 1 && \forall k \\ & \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, \forall k\lt 81 \\ & \color{darkred}x_{i,j,k}=1 && \text{for given cells}\\ & \color{darkred}x_{i,j,k} \in \{0,1\} \end{align}\]

In the neighbor constraint, I am careful about \(k\) being possibly outside the board, but less so about \(i\) and \(j\). The reason is, that with GAMS, when addressing using leads/lags outside of the domain, a value of zero is returned. That is exactly what we need here, for \(i\) and \(j\). The \(k\) index needs protection. The GAMS model can look like:

sets   i 'rows'    /i1*i9/   j 'columns' /j1*j9/   k 'values'  /1*81/ ; table board(i,j)         j1 j2 j3 j4 j5 j6 j7 j8 j9   i1   i2       25 20 13 14 15  8 41   i3       26                40   i4        1                39   i5       30                38   i6       67                37   i7       70                50   i8       73 72 79 60 61 52 53   i9 ; 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'    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); x.fx(i,j,k)$(board(i,j)=k.val) = 1; model m /all/; option mip=cplex; solve m minimizing dummy using mip; display x.l; parameter result(i,j); result(i,j) = sum(k, k.val*x.l(i,j,k)); option result:0; display result;

This model is solved completely in preprocessing. Cplex needs 0 iterations and 0 nodes. The output looks like:

----     46 PARAMETER result  

            j1          j2          j3          j4          j5          j6          j7          j8          j9

i1          23          22          21          12          11          10           9          42          43
i2          24          25          20          13          14          15           8          41          44
i3          27          26          19          18          17          16           7          40          45
i4          28           1           2           3           4           5           6          39          46
i5          29          30          31          32          33          34          35          38          47
i6          68          67          66          65          64          63          36          37          48
i7          69          70          71          80          81          62          51          50          49
i8          74          73          72          79          60          61          52          53          54
i9          75          76          77          78          59          58          57          56          55

 Notes:

• We don't need an optcr (gap) setting: we are looking for a feasible solution.
• GAMS has a rule: "return zero when a lead/lag is outside the domain". We need to know this when looking at equation next.
• The k+1 in next(i,j,k+1) means that we skip the last k when generating this equation. This protects us against generating \(0 \ge \color{darkred}x_{i,j,81}\). That would make the model infeasible.
• We were lucky with addressing \(i\pm1, j\pm 1\). GAMS returning zero there when outside the domain is exactly what we want. 
• We could have defined equation next in terms of \(k-1\) and \(k\). Then you would not need the k+1 in next(i,j,k+1). I.e., next(i,j,k).. x(i-1,j,k)+x(i+1,j,k)+x(i,j-1,k)+x(i,j+1,k) =g= x(i,j,k-1); This formulation does not need any protection, as the GAMS rule: "return zero when a lead/lag is outside the domain" is good in all cases. I still like the original formulation a little bit better as it is closer to what we stated in English.
• Interestingly when solving as an RMIP (Relaxed MIP), Cplex will need to do a few iterations. It will still find the integer solution. I suspect this is just by accident.
• A well-formed Numbrix puzzle has one unique solution. To verify this, we can add a cut to the model that forbids the current solution \(x^*\): \[\sum_{i,j,k} \color{darkblue}x^*_{i,j,k}\cdot \color{darkred} x_{i,j,k}\le 80\] and resolve. This must be infeasible. This is indeed the case for this input. In GAMS this test can be stated as:
  

  equation cut 'cut off current solution x.l'; cut.. sum((i,j,k)$(x.l(i,j,k)>0.5), x(i,j,k)) =l= card(k)-1; model m2 /all/; solve m2 minimizing dummy using mip; abort$(m2.modelstat <> %modelStat.infeasible% and        m2.modelstat <> %modelStat.integerInfeasible%) "Should be infeasible";

• Q: can we drop some given cells while making sure there is still one unique solution? A: yes
• Q: what is the minimum of given cells that we need to keep so the solution is unique? A: don't know. Looks like a difficult question to answer. Maybe write a model that tries to find two different solutions with given cells represented by binary variables.
• I think this model looks quite nice!

References

1. Numbrix puzzles by Marilyn vos Savant, https://parade.com/numbrix/
2. Sara Jensen, What is a Numbrix Puzzle? https://www.youtube.com/watch?v=FypFESO6rVo
3. Numbrix puzzle via TSP, http://yetanothermathprogrammingconsultant.blogspot.com/2020/12/numbrix-puzzle-via-tsp.html. Follow up post. Here we model the problem as a TSP problem. The MTZ formulation will give us the path as a byproduct. Although much smaller in size, the TSP problem is not as cheap to solve.