Problem
The 6snakes problem is a board game for one person. The board is a hexagon with 72 notches where the 72 pellets must be placed (see Figure 1). The pellets are divided into 6 groups of 12 pellets. In each group the pellets have the same color – hence there are 6 colors.
In each group, exactly one pellet is a little bit larger that is called the “head of the snake”. The inital situation of the game is to place the 6 heads of the snakes anywhere in a notch of the board. The goal is to place all other pellets such that the pellets with the same color form a “snake” on the board starting with the head of the snake, that is, each following pellets with the same color must be “neighbor” of the previous pellet with the same color and the snake cannot cross over another snake. A partial solution is given by Figure 2
The game is also called “Glastropfenspiel” or “Schlangengrube” in German.
Modeling Steps
The set of notches (or pellets (72 in total)) is defined as a set I = {1,…, 72} with the indexes i,j ∈ I. There are 6 groups (6 snakes) also introduced as a set H = {1,…, 6} with index h ∈ H.
First, the hexagonal board is mapped to a graphic and the x– and ycoordinates of all 72 notches are calculated. For this purpose, a grid is defined with the (0, 0) coordinate at the top left (see Figure 3, left side). There are 14 vertical and 12 horizontal lines in the grid. The vertical distances between the horizontal grid lines are fixed at 1 unit. Then the xcoordinates of the vertical grid lines are calculated as follows (to get the nice hexagonal arrangement of the notches inside the board):
Two notches i ∈ I and j ∈ I are defined to be “neighbors” if the Euclidean distance between the two notches is ≤ 2∕. If two notches are neighors then a line is drawn to link them (see Figure 3, right side).
To start the game, the 6 heads of the snakes (the six larger pellets) are placed anywhere on the board, for example at the notches (1, 29, 34, 39, 44, 50) as in Figure 4. A possible snake is also drawn in the same Figure. It consists of the head at notch 1, it continues to the neighbor notch 4, then to 8 and to 13, and so on, until the tail notch at 61. Hence, a “snake” is a path of length 11 (consisting of 12 pellets) starting at the head “jumping” only from neighbor to neighbor.
“The path distance of a notch from the head plus 1” is defined as position. Hence, the head has position 1 and its neighbor notch 4 has position 2, and the tail notch 61 has position 12. Let’s introduce an (ordered) set K = {1,…, 12} with index k ∈ K to denote the position of a notch inside a snake.
Now finally, we can introduce a binary variable z_{i,h,k} with the following meaning: z_{i,h,k} = 1 if and only if notch i is in snake h at the position k (otherwise z_{i,h,k} = 0).
Using this definition of the variables the model constraints are very straightforward:

Let a subset J ⊂ I with J = 6 be the head notches, then the six heads have position 1, that means:

Each notch is exactly in one snake and has a unique position:

The same snake and the same position is assigned to exactly one notch:

Finally, this is the key constraint: If a notch in a snake has position k then there must be at least one other neighbor that has position k + 1 (except the tail notch). One can formulate this as a logical constraint as follows (where e_{i,j} is a relation that is true if i and j are neighbors and K^{–} = K {12}
LPL translates this constraint automatically into a linear form as follows (which can be alternatively formulated).

We are looking only for a feasible solution, if one exists.
model snake "6Snakes Problem";
set i,j "set of notches";
e{i,j} "i and j are neighours";
h "set of snakes";
k "Length of snakes (ordered set of position)";
parameter J{h} "6 starting positions (heads of snake)";
binary variable z{i,h,k};
constraint
A{h}: z[J,h,1] = 1;
B{i}: sum{h,k} z = 1;
C{h,k}: sum{i} z = 1;
D{i,h,kk<#k}: sum{j in e} z[j,h,k+1] >= z[i,h,k]; //Sl(1,1000);
D{i,h,kk<#k}: z[i,h,k] > or{j in e} z[j,h,k+1];
solve;
end
Solution
A solution for the inital placements of the head as given above is given in Figure 5
Gurobi 9.0 takes 6secs to find this solution. The problem is not only very compact but also efficient. This underlines once again the importance to choose the “right” variable. The key in this variable is to integrate the position k inside the snake to get a connected path as is formulated in the last constraint.
Questions

For some initial placements of the head, there is no solution. Find one.

“Solve” the problem with the subset J in the previous answer.
Answers

A subset is J = {1, 2, 3, 70, 71, 72}. The solver returns an infeasibility. One can relax the last constraint using the LPL SL() function to find the best possible solution cutting one snake into two pieces.

The problem is infeasible, so no solution exists. However we can find a solution with a minimal “deviation” from a feasible solution. Replace the constraint
D{i,h,kk<#k}: sum{j in e} z[j,h,k+1] >= z[i,h,k]
with the constraint
D{i,h,kk<#k}: sum{j in e} z[j,h,k+1] >= z[i,h,k]  Sl(1,1000);
and solve it. The result is shown in Figure 6 (one snake is cut into two pieces).
References
[1] MatMod. Homepage for Learning Mathematical Modeling : https://matmod.ch.
[2] Hürlimann T. Reference Manual for the LPL Modeling Language, most recent version. https://matmod.ch/lpl/doc/manual.pdf.