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:
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 |
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) |
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'; |
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); |
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'; |
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; |
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 Model | MTZ Model | lifted MTZ Model | |
|---|---|---|---|
| Equations | 6,643 | 534 | 453 |
| Variables | 6,562 | 533 | 532 |
| Discrete variables | 6,561 | 450 | 450 |
| Time | 0.2 | 3.6 | 0.2 |
| Nodes | 0 | 538 | 0 |
| Iterations | 0 | 27005 | 0 |
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
- Numbrix puzzle via MIP, https://yetanothermathprogrammingconsultant.blogspot.com/2020/12/numbrix-puzzle-via-mip.html
- Traveling Salesman Problem, https://en.wikipedia.org/wiki/Travelling_salesman_problem
- Miller, C., A. Tucker, R. Zemlin. 1960. Integer programming formulation of traveling salesman problems. J. ACM 7 326-329
- Desrochers, M., Laporte, G., Improvements and extensions to the Miller-Tucker-Zemlin subtour elimination constraints, Operations Research Letters, 10 (1991) 27-36.
Might also be interesting to see if the following lifted MTZ constraints perform better than the classical ones:
ReplyDeleteu(j) =g= u(i0) + 1 - M*(1-x(i0,j)) + (M-2)*x(j,i0);
Looks like with this the presolve has some fun: "All rows and columns eliminated".
Delete