## Thursday, October 23, 2008

### Magic squares in GAMS

In http://yetanothermathprogrammingconsultant.blogspot.com/2008/10/magic-squares-with-ms-csp-solver.html we describe a CSP formulation using the Microsoft Solver Foundation tools. Here are some GAMS formulations using a MIP.

Instead of x(i,j) we use a single index: x(i). That makes it easier to model the "all-different" constraints, at the expense of making the row, column and diagonal sums more difficult. The latter complication can be alleviated by using a mapping set that maps between the two representations.

This version uses a permutation matrix paradigm to model the all-different constraint:

``\$ontextMagic Squares Take consecutive numbers from 1 to n^2 andarrange them in a n x n array such thatevery row, every column and the two diagonalsall add up to the same total. See also: http://forum.swarthmore.edu/alejandre/magic.square.html Erwin Kalvelagen, nov 1999\$offtext  set k "board size" /k1*k5/;alias (k,kk);set i "numbering of cells" /i1*i25/;alias(i,j); abort\$(sqr(card(k)) <> card(i)) "card(i) should be card(k)^k"; \$ontext  cells are numbered:       1  2  3  4      5  6  7  8      9 10 11 12     13 14 15 16 \$offtext   variable x(i)          "value of board cell i";binary variable y(i,j) "permutation matrix";variable s             "sum";  parameter p(i);p(i)=ord(i);  equations   unique(i) "make sure each x(i) is unique"   yrow(i)   "row sum of y"   ycol(j)   "col sum of y"   row(k)    "row of board"   col(k)    "col of board"   diag1     "first diagonal"   diag2     "second diagonal";   set map(k,kk,i);scalar cnt /0/;loop((k,kk),    cnt = cnt + 1;    map(k,kk,i)\$(ord(i)=cnt)=yes;);display map; set revdiag(i);loop(map(k,kk,i)\$((ord(k)+ord(kk))=card(k)+1),    revdiag(i) = yes;); display revdiag; unique(i)..x(i)=e=sum(j,y(i,j)*p(j));yrow(i)..  sum(j,y(i,j))=e=1;ycol(j)..  sum(i,y(i,j))=e=1; row(k).. sum(map(k,kk,i), x(i)) =e= s;col(k).. sum(map(kk,k,i), x(i)) =e= s;diag1..  sum(map(k,k,i), x(i)) =e= s;diag2..  sum(revdiag, x(revdiag)) =e= s;  equations   extra1    "extra constraint for performance"   extra2    "another redundant constraint for performance" ; * additional constraints belowscalar n  "size of set i";n = card(i);extra1.. sum(i,x(i)) =e= n*(n+1)/2;extra2.. sum((i,j),y(i,j)) =e= n; * fix s with its known valuescalar magicsum;magicsum = card(k)*(n+1)/2;s.fx = magicsum;  model m1 /all/;solve m1 using mip minimizing s; display x.l; parameter board(k,kk);board(k,kk) = sum(map(k,kk,i), x.l(i));display board;``

This version is using or-constraints to model the all-different condition.

``\$ontextMagic Squares Take consecutive numbers from 1 to n^2 andarrange them in a n x n array such thatevery row, every column and the two diagonalsall add up to the same total. See also: http://forum.swarthmore.edu/alejandre/magic.square.html Erwin Kalvelagen, nov 1999\$offtext set k "board size" /k1*k5/;alias (k,kk);set i "numbering of cells" /i1*i25/;alias(i,j); abort\$(sqr(card(k)) <> card(i)) "card(i) should be card(k)^k"; \$ontext  cells are numbered:       1  2  3  4      5  6  7  8      9 10 11 12     13 14 15 16 \$offtext   variable x(i)          "value of board cell i";binary variable y(i,j) "permutation matrix";variable s             "sum";  parameter p(i);p(i)=ord(i);  equations   unique(i) "make sure each x(i) is unique"   yrow(i)   "row sum of y"   ycol(j)   "col sum of y"   row(k)    "row of board"   col(k)    "col of board"   diag1     "first diagonal"   diag2     "second diagonal";   set map(k,kk,i);scalar cnt /0/;loop((k,kk),    cnt = cnt + 1;    map(k,kk,i)\$(ord(i)=cnt)=yes;);display map; set revdiag(i);loop(map(k,kk,i)\$((ord(k)+ord(kk))=card(k)+1),    revdiag(i) = yes;); display revdiag; unique(i)..x(i)=e=sum(j,y(i,j)*p(j));yrow(i)..  sum(j,y(i,j))=e=1;ycol(j)..  sum(i,y(i,j))=e=1; row(k).. sum(map(k,kk,i), x(i)) =e= s;col(k).. sum(map(kk,k,i), x(i)) =e= s;diag1..  sum(map(k,k,i), x(i)) =e= s;diag2..  sum(revdiag, x(revdiag)) =e= s;  equations   extra1    "extra constraint for performance"   extra2    "another redundant constraint for performance" ; * additional constraints belowscalar n  "size of set i";n = card(i);extra1.. sum(i,x(i)) =e= n*(n+1)/2;extra2.. sum((i,j),y(i,j)) =e= n; * fix s with its known valuescalar magicsum;magicsum = card(k)*(n+1)/2;s.fx = magicsum;  model m1 /all/;solve m1 using mip minimizing s; display x.l; parameter board(k,kk);board(k,kk) = sum(map(k,kk,i), x.l(i));display board;``

The mapping set looks like:

``----     60 SET map                  i1          i2          i3          i4          i5          i6          i7          i8          i9 k1.k1         YESk1.k2                     YESk1.k3                                 YESk1.k4                                             YESk1.k5                                                         YESk2.k1                                                                     YESk2.k2                                                                                 YESk2.k3                                                                                             YESk2.k4                                                                                                         YES    +         i10         i11         i12         i13         i14         i15         i16         i17         i18k2.k5         YESk3.k1                     YESk3.k2                                 YESk3.k3                                             YESk3.k4                                                         YESk3.k5                                                                     YESk4.k1                                                                                 YESk4.k2                                                                                             YESk4.k3                                                                                                         YES     +         i19         i20         i21         i22         i23         i24         i25 k4.k4         YESk4.k5                     YESk5.k1                                 YESk5.k2                                             YESk5.k3                                                         YESk5.k4                                                                     YESk5.k5                                                                                 YES``

The solution is displayed as:

``----     99 VARIABLE x.L  value of board cell i i1   6.000,    i2  22.000,    i3   7.000,    i4  13.000,    i5  17.000,    i6  20.000,    i7  24.000,    i8   2.000i9  18.000,    i10  1.000,    i11 19.000,    i12 12.000,    i13 16.000,    i14  8.000,    i15 10.000,    i16  9.000i17  3.000,    i18 25.000,    i19  5.000,    i20 23.000,    i21 11.000,    i22  4.000,    i23 15.000,    i24 21.000i25 14.000  ----    103 PARAMETER board               k1          k2          k3          k4          k5 k1       6.000      22.000       7.000      13.000      17.000k2      20.000      24.000       2.000      18.000       1.000k3      19.000      12.000      16.000       8.000      10.000k4       9.000       3.000      25.000       5.000      23.000k5      11.000       4.000      15.000      21.000      14.000``

The N=5 problem is already very difficult to solve.

#### 1 comment:

1. There are 275 million different order 5 magic squares. Walter Trump ... on his web site lists a program written in German Basic .. that will produce the solutions ... ? in about a hour.

Juro Bystricky has a free compiler ... K1compiler .. that has automatic backtracking and logic constraint ... that addresses magic square problems quite efficently