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:

$ontext
Magic Squares

Take consecutive numbers from 1 to n^2 and
arrange them in a n x n array such that
every row, every column and the two diagonals
all 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 below
scalar 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 value
scalar 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.

$ontext
Magic Squares

Take consecutive numbers from 1 to n^2 and
arrange them in a n x n array such that
every row, every column and the two diagonals
all 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 below
scalar 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 value
scalar 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 YES
k1.k2 YES
k1.k3 YES
k1.k4 YES
k1.k5 YES
k2.k1 YES
k2.k2 YES
k2.k3 YES
k2.k4 YES

+ i10 i11 i12 i13 i14 i15 i16 i17 i18

k2.k5 YES
k3.k1 YES
k3.k2 YES
k3.k3 YES
k3.k4 YES
k3.k5 YES
k4.k1 YES
k4.k2 YES
k4.k3 YES

+ i19 i20 i21 i22 i23 i24 i25

k4.k4 YES
k4.k5 YES
k5.k1 YES
k5.k2 YES
k5.k3 YES
k5.k4 YES
k5.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.000
i9 18.000, i10 1.000, i11 19.000, i12 12.000, i13 16.000, i14 8.000, i15 10.000, i16 9.000
i17 3.000, i18 25.000, i19 5.000, i20 23.000, i21 11.000, i22 4.000, i23 15.000, i24 21.000
i25 14.000


---- 103 PARAMETER board

k1 k2 k3 k4 k5

k1 6.000 22.000 7.000 13.000 17.000
k2 20.000 24.000 2.000 18.000 1.000
k3 19.000 12.000 16.000 8.000 10.000
k4 9.000 3.000 25.000 5.000 23.000
k5 11.000 4.000 15.000 21.000 14.000


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