Vale, la solución a este enigma está escrita en el blog que encontraste, pero ya que estás pidiendo ideas sobre cómo abordarlas, te daré la mía.
Encontré esta solución (usando GLPK) $$ \begin {array}{cccc} x_1 = 13 & x_2 = 8 & x_3 = 7 & x_4 = 6 \\ x_5 = 12 & x_6 = 2 & x_7 = 14 & x_8 = 11 \\ x_9 = 1 & x_{10} = 5 & x_{11} = 16 & x_{12} = 3 \\ x_{13} = 4 & x_{14} = 10 & x_{15} = 15 & x_{16} = 9 \end {array} $$ donde estoy usando la misma notación que la de arriba. También traté de entender lo que decía, pero obtuve un grupo de letras con sentido común, lo cual supongo que se debe a que existen, como lo notaste, al menos 8 soluciones simétricas. Vi en el blog que algunas personas resolvieron el enigma con sólo mirar las letras, pero yo soy más de números, así que me saltaré esa parte del "dicho".
La parte de álgebra lineal
Mi enfoque consiste en utilizar las 8 condiciones que ya has escrito, pero añadiendo la condición de que cada número entre el 1 y el 16 debe aparecer exactamente una vez. También tenemos que enfrentar el problema de que nuestras variables son números enteros, no números reales. Ni siquiera sé si es posible obtener tales condiciones usando sólo ecuaciones de álgebra lineal, un comienzo sería agregar esto $$ \sum_ {i=1}^{16} x_i = 136 $$ pero todavía nos faltan ecuaciones y no sé cómo proceder. Así que he usado la programación lineal, y GLPK para resolverlo dadas estas condiciones.
Programación lineal y enfoque utilizado (un poco largo)
No sé si está familiarizado con la programación lineal, pero le explicaré el procedimiento de manera que pueda ser entendido, con suerte.
Básicamente definimos un nuevo conjunto de variables $x_{ij}$ que son binarios (es decir $0$ o $1$ ) y se relacionan con sus variables de esta manera $$ x_i = \sum_ {j=1}^{16} x_{ij} \cdot j $$ Por supuesto, entonces sólo permitimos uno de estos $x_{ij}$ que sea diferente de cero para cada $i$ porque "todo lo desconocido es sólo un número", así que añadimos la condición $$ \sum_ {j=1}^{16} x_{ij} = 1 $$ De la misma manera, también queremos que "cada número aparezca sólo una vez", lo que da una condición muy similar (note el cambio en el índice de la suma) $$ \sum_ {i=1}^{16} x_{ij} = 1 $$ Finalmente, sustituimos este nuevo $x_{ij}$ en las 8 ecuaciones que teníamos al principio y ejecutar GLPK. Trata de encontrar una solución factible, y lo logra, dando la que escribí arriba. Supongo que en diferentes máquinas el programa podría encontrar diferentes soluciones, ya que ya sabemos que no es único.
Extra (casi encuentro el dicho)
Estaba un poco decepcionado por no entender el dicho, así que decidí hacer un poco de trampa, eché un vistazo al blog, el texto es " hacer un buen giro a diario ". Su solución implica que la Y sea la última letra, lo que significa $x_{12}=16$ . Así que puse esta restricción en mi programa y... todavía tengo una solución diferente. $$ \begin {array}{cccc} x_1 = 8 & x_2 = 4 & x_3 = 15 & x_4 = 7 \\ x_5 = 6 & x_6 = 12 & x_7 = 9 & x_8 = 10 \\ x_9 = 11 & x_{10} = 14 & x_{11} = 3 & x_{12} = 16 \\ x_{13} = 5 & x_{14} = 13 & x_{15} = 2 & x_{16} = 1 \end {array} $$ No estoy seguro de cuántas soluciones tiene esta molesta estrella :)