Saturday, April 25, 2009

Alphametics (1)

Alphametics are a type of mathematical puzzle. Given a set of words, written down as a long-hand addition, find the mapping between letters and digits (values). For example:

G E O R G I A
O R E G O N
V E R M O N T

V I R G I N I A

The solution is:

8 4 6 5 8 0 2
6 5 4 8 6 3
1 4 5 9 6 3 7

1 0 5 8 0 3 0 2

A somewhat simpler example is:

  S E N D
  M O R E

M O N E Y
  9 5 6 7
  1 0 8 5

1 0 6 5 2

A simple way to find the digits is using the following equation:

1000ys + 100ye + 10yn + yd + 1000ym + 100yo + 10yr + ye = 10000ym + 1000yo + 100yn + 10ye + yy

with the variables yi∈{0,1,...,9}. This is not sufficient as a solution yi=0 would be allowed. In addition we need to add an all-different constraint, and we may also need to forbid leading zero’s.

Forbidding leading zero’s is easy. We add the bounds: ys≥1 and ym≥1. The other constraint is somewhat more complicated. The uniqueness constraint can be implemented using an extra binary variable xi,j with form a permutation matrix.

We need all 9 y's so we introduce ydummy1 and ydummy2. This now results in the MIP model:
$ontext

    Alphametics, Simple formulation
    Erwin Kalvelagen, 1998

$offtext

set i /s,e,n,d,m,o,r,y,dummy1,dummy2/;
alias (i,j);

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.. 1000*y(’s’) + 100*y(’e’) + 10*y(’n’) + y(’d’) + 1000*y(’m’) + 100*y(’o’) + 10*y(’r’) + y(’e’) =e= 10000*y(’m’) + 1000*y(’o’) + 100*y(’n’) + 10*y(’e’) + y(’y’);

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(’s’) = 1; y.lo(’m’) = 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 is:

----     47 VARIABLE y.L  decision variables

s      9,    e      5,    n      6,    d      7,    m      1,    r      8,    y      2,    dummy1 4,    dummy2 3

Note that GAMS does not report zero’s. So you will need to deduct that the missing y(‘O’) is zero. This model solves quickly with all GAMS solvers except BDMLP (the LP seems to cycle on one of the nodes).

The objective function is somewhat interesting: the objective of the relaxed problem is always identical to the objective of the MIP. We can see this from the following:

ydef

As a result the MIP solver can stop as soon as it finds the first integer solution: it is immediately optimal. Of course a simpler way of achieving the same would be to use the following objective:

dummy_objective.. z =e= 0;

The bigger, first example can now be easily coded using:

set i /G,E,O,R,I,A,N,V,M,T/;  
addition..
     1e6*y('G') + 1e5*y('E') + 1e4*y('O') + 1e3*y('R') + 100*y('G') + 10*y('I') + y('A') +
     1e5*y('O') + 1e4*y('R') + 1e3*y('E') + 100*y('G') + 10*y('O') + y('N') +
     1e6*y('V') + 1e5*y('E') + 1e4*y('R') + 1e3*y('M') + 100*y('O') + 10*y('N') + y('T') =e=
     1e7*y('V') + 1e6*y('I') + 1e5*y('R') + 1e4*y('G') + 1e3*y('I') + 100*y('N') + 10*y('I') + y('A');

This formulation gives more problems, which is not surprising given the large coefficients. Some solvers will not be able to handle this problem. E.g. COINGLPK gives:

**** SOLVER STATUS     1 NORMAL COMPLETION        
**** MODEL STATUS      8 INTEGER SOLUTION         
**** OBJECTIVE VALUE               45.0000


**** REPORT SUMMARY :       10     NONOPT ( NOPT)

----     60 VARIABLE y.L  decision variables

G 1,    V 1

The reporting of NONOPT’s is probably just a bug in the GAMS link: with the given solver and model status there cannot be any NONOPT. More importantly, the reported solution is completely wrong.

Not all solvers fail on this example. E.g. Cplex gives:

**** SOLVER STATUS     1 NORMAL COMPLETION        
**** MODEL STATUS      1 OPTIMAL                  
**** OBJECTIVE VALUE               45.0000

----     60 VARIABLE y.L  decision variables

G 8,    E 4,    O 6,    R 5,    A 2,    N 3,    V 1,    M 9,    T 7

We can conclude y(‘I’)=0.

I will show how the formulation can be improved in a follow-up posting.