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.
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.
ReplyDeleteJuro Bystricky has a free compiler ... K1compiler .. that has automatic backtracking and logic constraint ... that addresses magic square problems quite efficently