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:
- Each cell can contain just one value. \[\sum_k \color{darkred}x_{i,j,k} = 1 \>\> \forall i,j\]
- Each value should occur exactly once. \[\sum_{i,j}\color{darkred}x_{i,j,k} = 1 \>\> \forall k\]
- 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}\]
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}\] |
sets |
---- 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
- Numbrix puzzles by Marilyn vos Savant, https://parade.com/numbrix/
- Sara Jensen, What is a Numbrix Puzzle? https://www.youtube.com/watch?v=FypFESO6rVo
- 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.
Thank you for sharing. It is always a joy reading your posts.
ReplyDeleteReally cool post! It would also be interesting to see how this compares to an MTZ-like model. And how well the two models scale to larger grids.
ReplyDeleteI have something working. First renumber the cells to node number. The network is sparse (only links to neighbors) and directed. Then add extra node n0 that serves as start/end. Add links from n0 to all nodes (and back). The (continuous) u variables in MTZ then give the ordering in the TSP and path in the puzzle. I used no objective (or all arc costs are 0). Cplex needed 827 nodes and 5 seconds. I'll try to clean up the model a bit.
Delete