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.