## 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:

 sets   r 'rows'    /r1*r9/   c 'columns' /c1*c9/   i0 'nodes with start/end node'  /n0*n81/   i(i0) 'nodes w/o start/end node' /n1*n81/   map(r,c,i) '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
r9.c9.n81


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   r1   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   r9 ; parameter given(i) 'this is in terms of nodes'; given(i) = sum(map(r,c,i),board(r,c)); 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) = yes; nb(r,c,r+1,c) = yes; nb(r,c,r,c-1) = yes; nb(r,c,r,c+1) = yes; * 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); arcs('n0',i) = yes; arcs(i,'n0') = yes; option arcs:0:0:8; display arcs; scalar narcs 'number of arcs'; narcs = card(arcs); 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:  alias(i0,j0); binary variable x(i0,j0); positive variable u(i0) 'ordering'; u.up(i0) = card(i); variable dummy; equations assign1(i0) 'assignment constraint' assign2(j0) 'assignment constraint' sequence(i0,j) 'subtour elimination constraints' obj 'dummy objective' ; assign1(i0).. sum(arcs(i0,j0), x(i0,j0)) =e= 1; assign2(j0).. sum(arcs(i0,j0), x(i0,j0)) =e= 1; scalar M; M = card(i0); u.fx('n0') = 0; sequence(arcs(i0,j)).. 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 = card(i0); u.fx('n0') = 0; u.lo(i) = 1; u.up(i) = n-1; sequence(arcs(i,j))..      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.

#### Performance

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
Equations6,643534453
Variables6,562533532
Discrete variables6,561450450
Time0.23.60.2
Nodes05380
Iterations0270050

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.

Notes:
• 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.

#### References

1. Numbrix puzzle via MIP, https://yetanothermathprogrammingconsultant.blogspot.com/2020/12/numbrix-puzzle-via-mip.html
2. Traveling Salesman Problem, https://en.wikipedia.org/wiki/Travelling_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.