Thursday, December 24, 2020

Numbrix puzzle via TSP

 In [1], I demonstrated a small and rather elegant MIP model to solve Numbrix puzzles. Here is the comment by Austin Buchanan:

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.

The Miller-Tucker-Zemlin [2,3] formulation for the Traveling Salesman Problem (TSP) starts with a standard assignment problem, and then adds some sub-tour elimination constraints. From [3]:

Note:  here \(\color{darkblue}p\ge \color{darkblue}n\). This 1960 paper has some interesting computational tests. Many runs ended with: "the machine was stopped after over 250 pivot steps had failed to produce the solution."

Basically, the \(\color{darkred}u_i\) variables indicate the order in which the cities are visited. I like to write this constraint as: \[\color{darkred}u_j \ge \color{darkred}u_i  + 1 - \color{darkblue}M\cdot (1-\color{darkred}x_{i,j})\] which is obviously equivalent to \[\color{darkred}x_{i,j}=1 \Rightarrow \color{darkred}u_j = \color{darkred}u_i  + 1\] if the bounds on \(\color{darkred}u\) are chosen correctly. 

These \(\color{darkred}u\) variables look like they can form our path in the Numbrix puzzle.

To create a TSP from the puzzle we need to look at it as a network. We can consider each cell as a node, and have arcs indicating neighboring cells. As in:

Of course, we have an \(9 \times 9\) grid with 81 nodes. We don't know where the node with value 1 and with 81 will end up, so the concept of a starting city is a bit difficult. One way to deal with this is to add an extra node that can function as the starting (and final) city. Just connect to all nodes. The picture becomes a bit messy:

The only thing we need to add is the collection of "given cells". They just correspond to fixing a few \(\color{darkred}u_i\). Interestingly, we don't have cost and thus no objective. Just a sparse network for which we need a feasible solution. A TSP solution would visit all nodes exactly once, and that is what we are looking for. 

Developing the network

We want to map between the grid board(r,c) and the network nodes \(i\) and the arcs. This code for the mapping can look like this:

'rows'    /r1*r9/
'columns' /c1*c9/

'nodes with start/end node'  /n0*n81/
'nodes w/o start/end node' /n1*n81/

'mapping (r,c) <-> i'

* mapping between (r,c) and i
map(r,c,i) = (
ord(r)-1)*card(r) + ord(c) = ord(i);
option map:0:0:8;
display map;

The constructed mapping is:

----     14 SET map  mapping (r,c) <-> i

r1.c1.n1 ,    r1.c2.n2 ,    r1.c3.n3 ,    r1.c4.n4 ,    r1.c5.n5 ,    r1.c6.n6 ,    r1.c7.n7 ,    r1.c8.n8 
r1.c9.n9 ,    r2.c1.n10,    r2.c2.n11,    r2.c3.n12,    r2.c4.n13,    r2.c5.n14,    r2.c6.n15,    r2.c7.n16
r2.c8.n17,    r2.c9.n18,    r3.c1.n19,    r3.c2.n20,    r3.c3.n21,    r3.c4.n22,    r3.c5.n23,    r3.c6.n24
r3.c7.n25,    r3.c8.n26,    r3.c9.n27,    r4.c1.n28,    r4.c2.n29,    r4.c3.n30,    r4.c4.n31,    r4.c5.n32
r4.c6.n33,    r4.c7.n34,    r4.c8.n35,    r4.c9.n36,    r5.c1.n37,    r5.c2.n38,    r5.c3.n39,    r5.c4.n40
r5.c5.n41,    r5.c6.n42,    r5.c7.n43,    r5.c8.n44,    r5.c9.n45,    r6.c1.n46,    r6.c2.n47,    r6.c3.n48
r6.c4.n49,    r6.c5.n50,    r6.c6.n51,    r6.c7.n52,    r6.c8.n53,    r6.c9.n54,    r7.c1.n55,    r7.c2.n56
r7.c3.n57,    r7.c4.n58,    r7.c5.n59,    r7.c6.n60,    r7.c7.n61,    r7.c8.n62,    r7.c9.n63,    r8.c1.n64
r8.c2.n65,    r8.c3.n66,    r8.c4.n67,    r8.c5.n68,    r8.c6.n69,    r8.c7.n70,    r8.c8.n71,    r8.c9.n72
r9.c1.n73,    r9.c2.n74,    r9.c3.n75,    r9.c4.n76,    r9.c5.n77,    r9.c6.n78,    r9.c7.n79,    r9.c8.n80

We need to translate the grid with givens from \((r,c)\) to \(i\). This we do with the almighty sum:

table board(r,c)
c1 c2 c3 c4 c5 c6 c7 c8 c9
r2       25 20 13 14 15  8 41
r3       26                40
r4        1                39
r5       30                38
r6       67                37
r7       70                50
r8       73 72 79 60 61 52 53

parameter given(i) 'this is in terms of nodes';
given(i) =
option given:0;
display given;

The display shows:

----     33 PARAMETER given  this is in terms of nodes

n11 25,    n12 20,    n13 13,    n14 14,    n15 15,    n16  8,    n17 41,    n20 26,    n26 40,    n29  1,    n35 39
n38 30,    n44 38,    n47 67,    n53 37,    n56 70,    n62 50,    n65 73,    n66 72,    n67 79,    n68 60,    n69 61
n70 52,    n71 53
The network setup is as follows. First, we establish the neighbors on the grid. This is in \((r,c)\) space. Next, we move this into the node space. Finally, we add the extra node \(n_0\).

set nb(r,c,r,c) 'neighbors';
nb(r,c,r-1,c) =
nb(r,c,r+1,c) =
nb(r,c,r,c-1) =
nb(r,c,r,c+1) =

* set up our network
alias (r,rr),(c,cc),(i,j);

set arcs(i0,i0) 'network topology';
arcs(i,j) =
sum((r,c,rr,cc)$(map(r,c,i) and map(rr,cc,j) and nb(r,c,rr,cc)),1);
'n0',i) = yes;
'n0') = yes;
option arcs:0:0:8;
display arcs;

scalar narcs 'number of arcs';
narcs =
option narcs:0;
display narcs;

The output is:

----     50 SET arcs  network topology

n0 .n1 ,    n0 .n2 ,    n0 .n3 ,    n0 .n4 ,    n0 .n5 ,    n0 .n6 ,    n0 .n7 ,    n0 .n8 
n0 .n9 ,    n0 .n10,    n0 .n11,    n0 .n12,    n0 .n13,    n0 .n14,    n0 .n15,    n0 .n16
n0 .n17,    n0 .n18,    n0 .n19,    n0 .n20,    n0 .n21,    n0 .n22,    n0 .n23,    n0 .n24
n0 .n25,    n0 .n26,    n0 .n27,    n0 .n28,    n0 .n29,    n0 .n30,    n0 .n31,    n0 .n32
n0 .n33,    n0 .n34,    n0 .n35,    n0 .n36,    n0 .n37,    n0 .n38,    n0 .n39,    n0 .n40
n0 .n41,    n0 .n42,    n0 .n43,    n0 .n44,    n0 .n45,    n0 .n46,    n0 .n47,    n0 .n48
n0 .n49,    n0 .n50,    n0 .n51,    n0 .n52,    n0 .n53,    n0 .n54,    n0 .n55,    n0 .n56
n0 .n57,    n0 .n58,    n0 .n59,    n0 .n60,    n0 .n61,    n0 .n62,    n0 .n63,    n0 .n64
n0 .n65,    n0 .n66,    n0 .n67,    n0 .n68,    n0 .n69,    n0 .n70,    n0 .n71,    n0 .n72
n0 .n73,    n0 .n74,    n0 .n75,    n0 .n76,    n0 .n77,    n0 .n78,    n0 .n79,    n0 .n80
n0 .n81,    n1 .n0 ,    n1 .n2 ,    n1 .n10,    n2 .n0 ,    n2 .n1 ,    n2 .n3 ,    n2 .n11
n3 .n0 ,    n3 .n2 ,    n3 .n4 ,    n3 .n12,    n4 .n0 ,    n4 .n3 ,    n4 .n5 ,    n4 .n13
n5 .n0 ,    n5 .n4 ,    n5 .n6 ,    n5 .n14,    n6 .n0 ,    n6 .n5 ,    n6 .n7 ,    n6 .n15
n7 .n0 ,    n7 .n6 ,    n7 .n8 ,    n7 .n16,    n8 .n0 ,    n8 .n7 ,    n8 .n9 ,    n8 .n17
n9 .n0 ,    n9 .n8 ,    n9 .n18,    n10.n0 ,    n10.n1 ,    n10.n11,    n10.n19,    n11.n0 
n11.n2 ,    n11.n10,    n11.n12,    n11.n20,    n12.n0 ,    n12.n3 ,    n12.n11,    n12.n13
n12.n21,    n13.n0 ,    n13.n4 ,    n13.n12,    n13.n14,    n13.n22,    n14.n0 ,    n14.n5 
n14.n13,    n14.n15,    n14.n23,    n15.n0 ,    n15.n6 ,    n15.n14,    n15.n16,    n15.n24
n16.n0 ,    n16.n7 ,    n16.n15,    n16.n17,    n16.n25,    n17.n0 ,    n17.n8 ,    n17.n16
n17.n18,    n17.n26,    n18.n0 ,    n18.n9 ,    n18.n17,    n18.n27,    n19.n0 ,    n19.n10
n19.n20,    n19.n28,    n20.n0 ,    n20.n11,    n20.n19,    n20.n21,    n20.n29,    n21.n0 
n21.n12,    n21.n20,    n21.n22,    n21.n30,    n22.n0 ,    n22.n13,    n22.n21,    n22.n23
n22.n31,    n23.n0 ,    n23.n14,    n23.n22,    n23.n24,    n23.n32,    n24.n0 ,    n24.n15
n24.n23,    n24.n25,    n24.n33,    n25.n0 ,    n25.n16,    n25.n24,    n25.n26,    n25.n34
n26.n0 ,    n26.n17,    n26.n25,    n26.n27,    n26.n35,    n27.n0 ,    n27.n18,    n27.n26
n27.n36,    n28.n0 ,    n28.n19,    n28.n29,    n28.n37,    n29.n0 ,    n29.n20,    n29.n28
n29.n30,    n29.n38,    n30.n0 ,    n30.n21,    n30.n29,    n30.n31,    n30.n39,    n31.n0 
n31.n22,    n31.n30,    n31.n32,    n31.n40,    n32.n0 ,    n32.n23,    n32.n31,    n32.n33
n32.n41,    n33.n0 ,    n33.n24,    n33.n32,    n33.n34,    n33.n42,    n34.n0 ,    n34.n25
n34.n33,    n34.n35,    n34.n43,    n35.n0 ,    n35.n26,    n35.n34,    n35.n36,    n35.n44
n36.n0 ,    n36.n27,    n36.n35,    n36.n45,    n37.n0 ,    n37.n28,    n37.n38,    n37.n46
n38.n0 ,    n38.n29,    n38.n37,    n38.n39,    n38.n47,    n39.n0 ,    n39.n30,    n39.n38
n39.n40,    n39.n48,    n40.n0 ,    n40.n31,    n40.n39,    n40.n41,    n40.n49,    n41.n0 
n41.n32,    n41.n40,    n41.n42,    n41.n50,    n42.n0 ,    n42.n33,    n42.n41,    n42.n43
n42.n51,    n43.n0 ,    n43.n34,    n43.n42,    n43.n44,    n43.n52,    n44.n0 ,    n44.n35
n44.n43,    n44.n45,    n44.n53,    n45.n0 ,    n45.n36,    n45.n44,    n45.n54,    n46.n0 
n46.n37,    n46.n47,    n46.n55,    n47.n0 ,    n47.n38,    n47.n46,    n47.n48,    n47.n56
n48.n0 ,    n48.n39,    n48.n47,    n48.n49,    n48.n57,    n49.n0 ,    n49.n40,    n49.n48
n49.n50,    n49.n58,    n50.n0 ,    n50.n41,    n50.n49,    n50.n51,    n50.n59,    n51.n0 
n51.n42,    n51.n50,    n51.n52,    n51.n60,    n52.n0 ,    n52.n43,    n52.n51,    n52.n53
n52.n61,    n53.n0 ,    n53.n44,    n53.n52,    n53.n54,    n53.n62,    n54.n0 ,    n54.n45
n54.n53,    n54.n63,    n55.n0 ,    n55.n46,    n55.n56,    n55.n64,    n56.n0 ,    n56.n47
n56.n55,    n56.n57,    n56.n65,    n57.n0 ,    n57.n48,    n57.n56,    n57.n58,    n57.n66
n58.n0 ,    n58.n49,    n58.n57,    n58.n59,    n58.n67,    n59.n0 ,    n59.n50,    n59.n58
n59.n60,    n59.n68,    n60.n0 ,    n60.n51,    n60.n59,    n60.n61,    n60.n69,    n61.n0 
n61.n52,    n61.n60,    n61.n62,    n61.n70,    n62.n0 ,    n62.n53,    n62.n61,    n62.n63
n62.n71,    n63.n0 ,    n63.n54,    n63.n62,    n63.n72,    n64.n0 ,    n64.n55,    n64.n65
n64.n73,    n65.n0 ,    n65.n56,    n65.n64,    n65.n66,    n65.n74,    n66.n0 ,    n66.n57
n66.n65,    n66.n67,    n66.n75,    n67.n0 ,    n67.n58,    n67.n66,    n67.n68,    n67.n76
n68.n0 ,    n68.n59,    n68.n67,    n68.n69,    n68.n77,    n69.n0 ,    n69.n60,    n69.n68
n69.n70,    n69.n78,    n70.n0 ,    n70.n61,    n70.n69,    n70.n71,    n70.n79,    n71.n0 
n71.n62,    n71.n70,    n71.n72,    n71.n80,    n72.n0 ,    n72.n63,    n72.n71,    n72.n81
n73.n0 ,    n73.n64,    n73.n74,    n74.n0 ,    n74.n65,    n74.n73,    n74.n75,    n75.n0 
n75.n66,    n75.n74,    n75.n76,    n76.n0 ,    n76.n67,    n76.n75,    n76.n77,    n77.n0 
n77.n68,    n77.n76,    n77.n78,    n78.n0 ,    n78.n69,    n78.n77,    n78.n79,    n79.n0 
n79.n70,    n79.n78,    n79.n80,    n80.n0 ,    n80.n71,    n80.n79,    n80.n81,    n81.n0 
n81.n72,    n81.n80

----     55 PARAMETER narcs                =          450  number of arcs

Note that this is a directed graph, so we have both \((n_1,n_2)\) and \((n_2,n_1)\).

With this under our belt, we can start writing down the optimization model:


binary variable x(i0,j0);
positive variable u(i0) 'ordering';
u.up(i0) =
variable dummy;

'assignment constraint'
'assignment constraint'
   sequence(i0,j)  'subtour elimination constraints'
   obj   'dummy objective'
sum(arcs(i0,j0), x(i0,j0)) =e= 1;
sum(arcs(i0,j0), x(i0,j0)) =e= 1;

scalar M;
M =

'n0') = 0;
     u(j) =g= u(i0) + 1 - M*(1-x(i0,j));

obj.. dummy =e= 0;

u.fx(i)$given(i) = given(i);

model mtz /all/;
solve mtz using mip minimizing dummy;

option u:0;
display u.l;

In normal TSP models we are interested in \(\color{darkred}x_{i,j}\), but here we want to know \(\color{darkred}u_i\).

----     89 VARIABLE u.L  ordering

n1  23,    n2  22,    n3  21,    n4  12,    n5  11,    n6  10,    n7   9,    n8  42,    n9  43,    n10 24,    n11 25
n12 20,    n13 13,    n14 14,    n15 15,    n16  8,    n17 41,    n18 44,    n19 27,    n20 26,    n21 19,    n22 18
n23 17,    n24 16,    n25  7,    n26 40,    n27 45,    n28 28,    n29  1,    n30  2,    n31  3,    n32  4,    n33  5
n34  6,    n35 39,    n36 46,    n37 29,    n38 30,    n39 31,    n40 32,    n41 33,    n42 34,    n43 35,    n44 38
n45 47,    n46 68,    n47 67,    n48 66,    n49 65,    n50 64,    n51 63,    n52 36,    n53 37,    n54 48,    n55 69
n56 70,    n57 71,    n58 80,    n59 81,    n60 62,    n61 51,    n62 50,    n63 49,    n64 74,    n65 73,    n66 72
n67 79,    n68 60,    n69 61,    n70 52,    n71 53,    n72 54,    n73 75,    n74 76,    n75 77,    n76 78,    n77 59
n78 58,    n79 57,    n80 56,    n81 55
Note that \(\color{darkred}u_{0}=0\) is not printed (only nonzero values are printed).

To present the solution in a more meaningful way, we need to map it back to an \((r,c)\) representation. This will do that:

parameter sol(r,c) 'solution';
sol(r,c) =
sum(map(r,c,i), round(u.l(i)));
option sol:0;
display sol;

The display output looks like this:

----     94 PARAMETER sol  solution

            c1          c2          c3          c4          c5          c6          c7          c8          c9

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

MTZ variants

The original paper [3] has: \[\begin{align} &\color{darkred}u_i -\color{darkred}u_j +\color{darkblue}n \cdot \color{darkred}x_{i,j} \le \color{darkblue}n-1 && 1 \le  i\ne  j \le \color{darkblue}n\\ &  \color{darkred}u_i \in [0,\color{darkblue}n-1]\end{align} \] where \(\color{darkblue}n\) is the total number of nodes. Most of the times we see these written as: \[\begin{align} &\color{darkred}u_i -\color{darkred}u_j +(\color{darkblue}n-1) \cdot \color{darkred}x_{i,j} \le \color{darkblue}n-2&& 2 \le  i\ne  j \le \color{darkblue}n\\ & \color{darkred}u_i \in [1,\color{darkblue}n-1]\end{align} \] This implicitly fixes the starting  node \(\color{darkred}u_1=0\).  In [4] a tighter formulation is given by a technique called "lifting": \[\begin{align} &\color{darkred}u_i -\color{darkred}u_j +(\color{darkblue}n-1) \cdot \color{darkred}x_{i,j}+(\color{darkblue}n-3)\cdot \color{darkred}x_{j,i}  \le \color{darkblue}n-2 &&2 \le  i\ne  j \le \color{darkblue}n\\ & \color{darkred}u_i \in [1,\color{darkblue}n-1]\end{align} \]

This particular formulation can be written in GAMS as:

scalar n;
n =

'n0') = 0;
u.lo(i) = 1;
u.up(i) = n-1;

     u(i)-u(j) + (n-1)*x(i,j) + (n-3)*x(j,i) =l= n-2;

This is essentially the same as what Rob Pratt suggested in the comments. The difference is that I added that we only consider nodes 2 through \(n\). That reduces the size of the problem a bit.


Let's compare how the solver does on these three models: the original model from [1] and the original MTZ formulation and the lifted MTZ formulation here:

Original ModelMTZ Modellifted MTZ Model
Discrete variables6,561450450

What we see is that the original model is much bigger, but solves during preprocessing (heuristics before starting on the relaxation). The number of constraints and variables is not a good measure of the difficulty of a MIP model. The MTZ model is not as tight, but we can save it by using the listed version as suggested by Rob Pratt. This version solves completely in the presolve: all rows and columns were eliminated.

  • This is just one data set. So not much of a "computational study".
  • The MTZ formulation is actually quite fast for an 82 node problem. Reasons are: (a) we are only looking for a feasible solution and (b) the network is sparse (which we exploited in the model formulation).
  • There is a bit more data manipulation going on here. Mapping from one representation to another and back is a common technique in optimization modeling. In GAMS that often means: using multi-dimensional mapping sets. Note that the statement to map from \((r,c)\) to \(i\) is essentially the same as for mapping back (we use the same mapping set).
  • This is an interesting model.


  1. Numbrix puzzle via MIP, 
  2. Traveling Salesman Problem,
  3. Miller, C., A. Tucker, R. Zemlin. 1960. Integer programming formulation of traveling salesman problems. J. ACM 7 326-329
  4. Desrochers, M., Laporte, G., Improvements and extensions to the Miller-Tucker-Zemlin subtour elimination constraints, Operations Research Letters, 10 (1991) 27-36.


  1. Might also be interesting to see if the following lifted MTZ constraints perform better than the classical ones:
    u(j) =g= u(i0) + 1 - M*(1-x(i0,j)) + (M-2)*x(j,i0);

    1. Looks like with this the presolve has some fun: "All rows and columns eliminated".