Sunday, April 26, 2009

Alphametics (3)

In http://www-cs-faculty.stanford.edu/~knuth/fasc2b.ps.gz we find the huge alphametic (I like the term cryptarithm mentioned as alternative even better):

A N
A C C E L E R A T I N G
I N F E R E N T I A L
E N G I N E E R I N G
T A L E
E L I T E
G R A N T
F E E
E T
C E T E R A

A R T I F I C I A L
I N T E L L I G E N C E

The simple formulation from this posting leads to problems even for Cplex. The stable formulation from this posting is a winner here: it yields a quick solution that is correct.

Simple formulation

$ontext

Alphametics, Simple formulation
This does not work very well.

Erwin Kalvelagen

AN
+ ACCELERATING
+ INFERENTIAL
+ ENGINEERING
+ TALE
+ ELITE
+ GRANT
+ FEE
+ ET
+ CETERA =
-------------
ARTIFICIAL
+ INTELLIGENCE



$offtext

set i /A,N,C,E,L,R,T,I,G,F/;
alias (i,j);

set leading(i) /A,I,E,T,G,F,C/;

abort$(card(i) <> 10) "set i should have 10 elements";

parameter v(i); v(i) = ord(i)-1;

variables
y(i) 'decision variables'
x(i,j) 'auxiliary variables for uniqueness'
z 'dummy objective variable'
;
binary variable x;

equation
addition 'the actual problem'
ydef(i) 'calculate y'
xrow(i) 'row sums'
xcol(j) 'column sums'
dummy_objective
;

addition..
10*y('A')+y('N')+
1e11*y('A')+1e10*y('C')+1e9*y('C')+1e8*y('E')+1e7*y('L')+1e6*y('E')+1e5*y('R')+1e4*y('A')+1e3*y('T')+100*y('I')+10*y('N')+y('G')+
1e10*y('I')+1e9*y('N')+1e8*y('F')+1e7*y('E')+1e6*y('R')+1e5*y('E')+1e4*y('N')+1e3*y('T')+100*y('I')+10*y('A')+y('L')+
1e10*y('E')+1e9*y('N')+1e8*y('G')+1e7*y('I')+1e6*y('N')+1e5*y('E')+1e4*y('E')+1e3*y('R')+100*y('I')+10*y('N')+y('G')+
1e3*y('T')+100*y('A')+10*y('L')+y('E')+
1e4*y('E')+1e3*y('L')+100*y('I')+10*y('T')+y('E')+
1e4*y('G')+1e3*y('R')+100*y('A')+10*y('N')+y('T')+
100*y('F')+10*y('E')+y('E')+
10*y('E')+y('T')+
1e5*y('C')+1e4*y('E')+1e3*y('T')+100*y('E')+10*y('R')+y('A') =e=
1e9*y('A')+1e8*y('R')+1e7*y('T')+1e6*y('I')+1e5*y('F')+1e4*y('I')+1e3*y('C')+100*y('I')+10*y('A')+y('L')+
1e11*y('I')+1e10*y('N')+1e9*y('T')+1e8*y('E')+1e7*y('L')+1e6*y('L')+1e5*y('I')+1e4*y('G')+1e3*y('E')+100*y('N')+10*y('C')+y('E');


ydef(i).. y(i) =e= sum(j, x(i,j)*v(j));
xrow(i).. sum(j, x(i,j)) =e= 1;
xcol(j).. sum(i, x(i,j)) =e= 1;

*
* bounds on leading digits: they can not be 0
*
y.lo(leading) = 1;

dummy_objective.. z =e= sum(i, y(i));

model m /all/;
m.optfile=1;
solve m using mip minimizing z;

option y:0;
display y.l;


The numbers in the addition equation are really become too big for comfort. It is surprising that some solvers such as Gurobi can deliver a good solution for this model, but as expected some others can not. The poor scaling can create severe havoc. The model can be downloaded from here.



Stable formulation






$ontext

Alphametics, stable formulation

AN
+ ACCELERATING
+ INFERENTIAL
+ ENGINEERING
+ TALE
+ ELITE
+ GRANT
+ FEE
+ ET
+ CETERA =
-------------
ARTIFICIAL
+ INTELLIGENCE

Erwin Kalvelagen, 2009

$offtext

set
i 'letters' /A,N,C,E,L,R,T,I,G,F/
w 'words' /an,accelerating,inferential,engineering,tale,elite,grant,fee,et,cetera,
artificial,intelligence/
w1(w) /an,accelerating,inferential,engineering,tale,elite,grant,fee,et,cetera/
w2(w) /artificial,intelligence/
;

alias (i,j);
abort$(card(i) <> 10) "set i should have 10 elements";

set
k 'slices (slice1 if rightmost slice)' /slice1*slice12/
lead(i) 'leading letters' /A,I,E,T,G,F,C/
;

set letters(w,i,k) /
an.(a.slice2, n.slice1)
accelerating.(a.slice12, c.slice11, c.slice10, e.slice9, l.slice8, e.slice7,
r.slice6, a.slice5, t.slice4, i.slice3, n.slice2, g.slice1)
inferential.(i.slice11, n.slice10, f.slice9, e.slice8, r.slice7, e.slice6,
n.slice5, t.slice4, i.slice3, a.slice2, l.slice1)
engineering.(e.slice11, n.slice10, g.slice9, i.slice8, n.slice7, e.slice6,
e.slice5, r.slice4, i.slice3, n.slice2, g.slice1)
tale.(t.slice4, a.slice3, l.slice2, e.slice1)
elite.(e.slice5, l.slice4, i.slice3, t.slice2, e.slice1)
grant.(g.slice5, r.slice4, a.slice3, n.slice2, t.slice1)
fee.(f.slice3, e.slice2, e.slice1)
et.(e.slice2, t.slice1)
cetera.(c.slice6, e.slice5, t.slice4, e.slice3, r.slice2, a.slice1)
artificial.(a.slice10, r.slice9, t.slice8, i.slice7,
f.slice6, i.slice5, c.slice4, i.slice3, a.slice2, l.slice1)
intelligence.(i.slice12, n.slice11, t.slice10, e.slice9, l.slice8, l.slice7,
i.slice6, g.slice5, e.slice4, n.slice3, c.slice2, e.slice1)

/;

parameter lhs(k,i), rhs(k,i);
lhs(k,i) = sum(letters(w1,i,k),1);
rhs(k,i) = sum(letters(w2,i,k),1);

parameter v(i); v(i) = ord(i)-1;

variables
y(i) 'decision variables'
x(i,j) 'auxiliary variables for uniqueness'
carry(k)
z 'dummy objective variable'
;
binary variable x;
integer variable carry;

equation
addition(k) 'the actual problem'
ydef(i) 'calculate y'
xrow(i) 'row sums'
xcol(j) 'column sums'
dummy_objective

;


addition(k)..
carry(k-1) + sum(i, lhs(k,i)*y(i)) =e=
sum(i, rhs(k,i)*y(i)) + 10*carry(k);

carry.fx('slice12')=0;

ydef(i).. y(i) =e= sum(j, x(i,j)*v(j));
xrow(i).. sum(j, x(i,j)) =e= 1;
xcol(j).. sum(i, x(i,j)) =e= 1;

*
* bounds on leading digits: they can not be 0
*
y.lo(lead) = 1;

dummy_objective.. z =e= sum(i, y(i));


model m /all/;
solve m using mip minimizing z;
option y:0;
display y.l;

The solution I get is:






----    106 VARIABLE y.L  decision variables




A 5,    N 9,    C 7,    E 4,    R 2,    T 1,    I 6,    G 8,    F 3



where y(‘L’)=0 is implied. This formulation is stable: it yields correct solutions for all MIP solvers. The model can be downloaded here.