Wednesday, April 29, 2009

Excel Solver Model conversion

I was involved in a small project where we converted a spreadsheet model using Solver to GAMS/CONOPT. I was unable to reproduce the spreadsheet results, and it turns out the spreadsheet had some errors in the formula’s. I have converted quite a number of Excel models, and almost all had some problems that were uncovered by moving the model to a modeling language environment. See table 1 in http://panko.shidler.hawaii.edu/SSR/Mypapers/whatknow.htm for more evidence of how difficult it is to create error-free spreadsheet models.

PS: after fixing the spreadsheet errors we were able to get similar results. This is actually not a bad approach to verify and repair spreadsheet models.

Ceiling in OML

Hello, I have to deal with an objective function of the following form: min [C1 * ceiling (X/m) plus C2 * (X-ceiling(X/m))];
where C1, C2, m are parameters and X is the decision variable. Ceiling is a mathematical operator that rounds a number up, to the nearest integer (see Excel's definition).
My first question is: is the Ceiling operator defined in OML?
My second question: is there a solver within the Microsoft Solver Foundation that can handle this type of objective function?
Thank you.

I believe you can do the following in OML:

Let Y be an integer variable with the constraint

X/m ≤ Y ≤ 1+X/m

then minimize

C1 * Y + C2 * (X-Y)

Hopefully the rest of the model is linear. In that case this is a much better formulation as it gets rid of a non-smooth function: solvers often have troubles with that.

See also http://yetanothermathprogrammingconsultant.blogspot.com/2009/02/floor-in-minlp-model.html and http://yetanothermathprogrammingconsultant.blogspot.com/2008/08/gams-documentation-endogenous-floor.html.

Sunday, April 26, 2009

Alphametics (4)

A CSP solver with the All-different constraint can make the modeling of alphametic problems even easier. For instance with Microsoft Solver Foundation:

Model[
   Decisions[Integers[0,9],s,e,n,d,m,o,r,y],
   Constraints[
     1e3*s + 100*e + 10*n + d + 1000*m + 100*o + 10*r + e == 10000*m + 1000*o + 100*n + 10*e + y,
     Unequal[s,e,n,d,m,o,r,y],
     s>=1,m>=1]
]

This solves just fine, but now lets move to the more difficult problem:
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
I encoded this as:

Model[
  Decisions[Integers[0,9],A,N,C,E,L,R,T,I,G,F],
  Constraints[
    10*A+N+
    1e11*A+1e10*C+1e9*C+1e8*E+1e7*L+1e6*E+1e5*R+1e4*A+1e3*T+100*I+10*N+G+
    1e10*I+1e9*N+1e8*F+1e7*E+1e6*R+1e5*E+1e4*N+1e3*T+100*I+10*A+L+
    1e10*E+1e9*N+1e8*G+1e7*I+1e6*N+1e5*E+1e4*E+1e3*R+100*I+10*N+G+
    1e3*T+100*A+10*L+E+
    1e4*E+1e3*L+100*I+10*T+E+
    1e4*G+1e3*R+100*A+10*N+T+
    100*F+10*E+E+
    10*E+T+
    1e5*C+1e4*E+1e3*T+100*E+10*R+A
       ==
    1e9*A+1e8*R+1e7*T+1e6*I+1e5*F+1e4*I+1e3*C+100*I+10*A+L+
    1e11*I+1e10*N+1e9*T+1e8*E+1e7*L+1e6*L+1e5*I+1e4*G+1e3*E+100*N+10*C+E,
    A>=1, I>=1,E>=1,T>=1,G>=1,F>=1,C>=1,
    Unequal[A,N,C,E,L,R,T,I,G,F]
  ]
]

This model produces Infeasible:

alphameticcsp2

I had hopes that the rational arithmetic of MSF could help here. Of course we can apply the more stable formulation we used in the MIP formulation here.

Model[
  Decisions[Integers[0,9],A,N,C,E,L,R,T,I,G,F],
  Decisions[Integers[0,Infinity],carry1,carry2,carry3,carry4,carry5,carry6,carry7,carry8,carry9,carry10,carry11],
  Constraints[
    N + G + L + G + E + E + T + E + T + A == L + E + 10*carry1,
    A + N + A + N + L + T + N + E + E + R + carry1 == A + C + 10*carry2,
    I + I + I + A + I + A + F + E + carry2 == I + N + 10*carry3,
    T + T + R + T + L + R + T + carry3 == C + E + 10*carry4,
    A + N + E + E + G + E + carry4 == I + G + 10*carry5,
    R + E + E + C + carry5 == F + I + 10*carry6,
    E + R + N + carry6 == I + L + 10*carry7,
    L + E + I + carry7 == T + L + 10*carry8,
    E + F + G + carry8 == R + E + 10*carry9,
    C + N + N + carry9 == A + T + 10*carry10,
    C + I + E + carry10 == N + 10*carry11,
    A + carry11 == I,
    A>=1, I>=1,E>=1,T>=1,G>=1,F>=1,C>=1,
    Unequal[A,N,C,E,L,R,T,I,G,F]
  ]
]

This now solves quickly:

alphameticcsp1

So why is the simple formulation not working. It seems related to the Unequal construct. When we feed the solution into the problem we get:

  1. Infeasible if Unequal statement is part of the model.
  2. Feasible is Unequal statement is removed.

This may indicate a bug or some intricate numerics suddenly playing a role.

Model[
  Decisions[Integers[0,9],A,N,C,E,L,R,T,I,G,F],
  Constraints[
    10*A+N+
    1e11*A+1e10*C+1e9*C+1e8*E+1e7*L+1e6*E+1e5*R+1e4*A+1e3*T+100*I+10*N+G+
    1e10*I+1e9*N+1e8*F+1e7*E+1e6*R+1e5*E+1e4*N+1e3*T+100*I+10*A+L+
    1e10*E+1e9*N+1e8*G+1e7*I+1e6*N+1e5*E+1e4*E+1e3*R+100*I+10*N+G+
    1e3*T+100*A+10*L+E+
    1e4*E+1e3*L+100*I+10*T+E+
    1e4*G+1e3*R+100*A+10*N+T+
    100*F+10*E+E+
    10*E+T+
    1e5*C+1e4*E+1e3*T+100*E+10*R+A
       ==
    1e9*A+1e8*R+1e7*T+1e6*I+1e5*F+1e4*I+1e3*C+100*I+10*A+L+
    1e11*I+1e10*N+1e9*T+1e8*E+1e7*L+1e6*L+1e5*I+1e4*G+1e3*E+100*N+10*C+E,
   A>=1, I>=1,E>=1,T>=1,G>=1,F>=1,C>=1,
  //  Unequal[A,N,C,E,L,R,T,I,G,F],
    A==5,N==9,C==7,E==4,L==0,R==2,T==1,I==6,G==8,F==3
  ]
]

Actually I would have expected problems in the big equality condition and not in the Unequal constraint.

Update: I understand now the behavior a little bit more. Without the Unequal the MIP solver Gurobi is called. This produces the correct conclusion: Feasible. With the Unequal statement enabled the CSP solver is used and that solver has problems with this model. Actually Gurobi uses floating point instead of rational arithmetic. So it is somewhat surprising that Gurobi is more accurate on this model than the CSP solver.

Update2: see expert explanation in comment. CSP wants all intermediate results (I think all terms used in any expression) also to fit in an integer. 

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.

Saturday, April 25, 2009

Alphametics (2)

The problem we are trying to solve is the alphametic puzzle:

alphametic3

In this posting we developed an equation

alphameticeq

that caused some problems.

A different formulation can be developed that circumvents this problem. We can use the addition algorithm taught to us in primary school: first add the right most digits, and save a possible carry. Then add up the digits in the second rightmost column, including the previous carry, etc. This scheme can actually be implemented in constraints as follows:

alphameticeq2

The complete model can look like:

$ontext

  Alphametics, stable formulation

     G E O R G I A        8 4 6 5 8 0 2
       O R E G O N   ->     6 5 4 8 6 3
     V E R M O N T        1 4 5 9 6 3 7
     -------------      ---------------
   V I R G I N I A      1 0 5 8 0 3 0 2

   Erwin Kalvelagen, 1998, updated 2009

$offtext

set
   i 'letters' /g,e,o,r,i,a,n,v,m,t/
   w 'words'   /georgia, oregon, vermont, virginia/
   wadd(w)     /georgia, oregon, vermont /
;

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

set
  k 'slices (slice1 if rightmost slice)' /slice1*slice8/
  lead(i) 'leading letters' /g,o,v/
;

set letters(w,i,k) /
georgia.(
    g.slice7
    e.slice6
    o.slice5
    r.slice4
    g.slice3
    i.slice2
    a.slice1
   )
oregon.(
    o.slice6
    r.slice5
    e.slice4
    g.slice3
    o.slice2
    n.slice1
)
vermont.(
    v.slice7
    e.slice6
    r.slice5
    m.slice4
    o.slice3
    n.slice2
    t.slice1
)
virginia.(
    v.slice8
    i.slice7
    r.slice6
    g.slice5
    i.slice4
    n.slice3
    i.slice2
    a.slice1
)
/;

parameter lhs(k,i), rhs(k,i);
lhs(k,i) = sum(letters(wadd,i,k),1);
rhs(k,i) = sum(letters('virginia',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('slice8')=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;

This model solves without a problem even with the GAMS version of GLPK:

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

----    113 VARIABLE y.L  decision variables

g 8,    e 4,    o 6,    r 5,    a 2,    n 3,    v 1,    m 9,    t 7

The only small bug here is that this solution is actually optimal, so the model status is reported incorrectly.

GLPK: gams vs glpsol

The GAMS version of GLPK seems numerically less stable than the standalone GLPK version (at least for this case). Here is a GAMS model, where we used an option file to not use a final LP to get duals:

$ontext

    Alphametics, Simple formulation
    Erwin Kalvelagen, 1998

     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

$offtext

set i /G,E,O,R,I,A,N,V,M,T/;
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..
              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');

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('G') = 1; y.lo('V') = 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;

$onecho > coinglpk.opt
solvefinal=0
$offecho

This gives:

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

---- VAR y  decision variables

         LOWER          LEVEL          UPPER         MARGINAL

G         1.0000         1.0000        +INF            EPS     NOPT
E        -INF             .            +INF            EPS     NOPT
O        -INF             .            +INF            EPS     NOPT
R        -INF             .            +INF            EPS     NOPT
I        -INF             .            +INF            EPS     NOPT
A        -INF             .            +INF            EPS     NOPT
N        -INF             .            +INF            EPS     NOPT
V         1.0000         1.0000        +INF            EPS     NOPT
M        -INF             .            +INF            EPS     NOPT
T        -INF             .            +INF            EPS     NOPT

This solution is not correct. I suspect the NOPTS are just the result of a bug in GAMS link. As explained in a previous posting this behavior may be because of scaling. The log file looks like:

--- Job alphametic2.gms Start 04/25/09 19:02:59 WEX-VIS 23.0.2 x86/MS Windows       
GAMS Rev 230  Copyright (C) 1987-2009 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen                               G090213/0001CV-WIN
          Amsterdam Optimization Modeling Group                      DC4572
--- Starting compilation
--- alphametic2.gms(63) 3 Mb
--- Starting execution: elapsed 0:00:00.004
--- alphametic2.gms(55) 4 Mb
--- Generating MIP model m
--- alphametic2.gms(56) 4 Mb
---   32 rows  111 columns  320 non-zeroes
---   100 discrete-columns
--- Executing COINGLPK: elapsed 0:00:00.009

GAMS/CoinGlpk LP/MIP Solver (Glpk Library 4.32 )
written by A. Makhorin
Reading parameter(s) from "C:\projects\test\coinglpk.opt"
>>  solvefinal=0
Finished reading from "C:\projects\test\coinglpk.opt"

Starting GLPK Branch and Bound...
ipp_basic_tech:  0 row(s) and 0 column(s) removed
ipp_reduce_bnds: 4 pass(es) made, 14 bound(s) reduced
ipp_basic_tech:  0 row(s) and 3 column(s) removed
ipp_reduce_coef: 1 pass(es) made, 0 coefficient(s) reduced
glp_intopt: presolved MIP has 31 rows, 107 columns, 300 non-zeros
glp_intopt: 97 integer columns, all of which are binary
Scaling...
A: min|aij| = 1.000e+000  max|aij| = 9.000e+006  ratio = 9.000e+006
GM: min|aij| = 1.304e-001  max|aij| = 7.666e+000  ratio = 5.877e+001
EQ: min|aij| = 1.707e-002  max|aij| = 1.000e+000  ratio = 5.858e+001
2N: min|aij| = 9.766e-003  max|aij| = 1.500e+000  ratio = 1.536e+002
Crashing...
Size of triangular part = 29
Solving LP relaxation...
      0: obj = -1.766865672e+002  infeas = 1.380e+003 (2)
*    37: obj =  4.500000000e+001  infeas = 3.052e-014 (1)
OPTIMAL SOLUTION FOUND
Integer optimization begins...
Gomory's cuts enabled
Cover cuts enabled
Clique cuts enabled
Creating the conflict graph...
The conflict graph has 2*97 vertices and 955 edges
+    37: mip =     not found yet >=              -inf        (1; 0)
+ 62738: mip =     not found yet >=  4.500000000e+001        (231; 4046)
+ 64138: >>>>>  4.500000000e+001 >=  4.500000000e+001 < 0.1 (236; 4157)
+ 64138: mip =  4.500000000e+001 >=  4.500000000e+001 < 0.1 (236; 4157)
RELATIVE MIP GAP TOLERANCE REACHED; SEARCH TERMINATED

Solving fixed problem skipped.

Integer Solution.
Best solution:                   45   (64138 iterations, 5.137 seconds)
Writing primal solution. Objective: 45 Time: 5.137 s
--- Restarting execution
--- alphametic2.gms(56) 0 Mb
--- Reading solution for model m
--- Executing after solve: elapsed 0:00:05.203
--- alphametic2.gms(59) 3 Mb
*** Status: Normal completion
--- Job alphametic2.gms Stop 04/25/09 19:03:04 elapsed 0:00:05.204

The same model formulated in AMPL or rather GNU Mathprog looks like:

set I;
var x{i in I, j in I}, binary;
var y{i in I}, >= 0, <= 9;
param v{i in I};
minimize obj: sum{i in I} y[i];
addition:
     1000000*y['G'] + 100000*y['E'] + 10000*y['O'] + 1000*y['R'] + 100*y['G'] + 10*y['I'] + y['A'] +
        100000*y['O'] + 10000*y['R'] + 1000*y['E'] + 100*y['G'] + 10*y['O'] + y['N'] +
        1000000*y['V'] + 100000*y['E'] + 10000*y['R'] + 1000*y['M'] + 100*y['O'] + 10*y['N'] + y['T'] =
     10000000*y['V'] + 1000000*y['I'] + 100000*y['R'] + 10000*y['G'] + 1000*y['I'] + 100*y['N'] + 10*y['I'] + y['A'];
ydef{i in I}: y[i] = sum{j in I} x[i,j]*v[j];
xrow{i in I}: sum{j in I} x[i,j] = 1;
xcol{j in I}: sum{i in I} x[i,j] = 1;
bound1: y['G']>=1;
bound2: y['V']>=1;
solve;
display y;
data;
set I := G E O R I A N V M T;
param v := G 0  E 1  O 2  R 3  I 4  A 5  N 6  V 7  M 8  T 9;
end;

The result looks much better:

>cmd /c path=%path%;C:\projects\glpk\glpk\gusek&& glpsol.exe --cover --clique --gomory --mir -m "alphametic.mod"    && if not '00'=='00' gusek.exe 
Reading model section from alphametic.mod...
Reading data section from alphametic.mod...
21 lines were read
Generating obj...
Generating addition...
Generating ydef...
Generating xrow...
Generating xcol...
Generating bound1...
Generating bound2...
Model has been successfully generated
ipp_basic_tech:  3 row(s) and 0 column(s) removed
ipp_reduce_bnds: 5 pass(es) made, 31 bound(s) reduced
ipp_basic_tech:  3 row(s) and 27 column(s) removed
ipp_reduce_coef: 1 pass(es) made, 0 coefficient(s) reduced
glp_intopt: presolved MIP has 28 rows, 83 columns, 230 non-zeros
glp_intopt: 74 integer columns, all of which are binary
Scaling...
A: min|aij| = 1.000e+000  max|aij| = 1.001e+006  ratio = 1.001e+006
GM: min|aij| = 1.619e-001  max|aij| = 6.178e+000  ratio = 3.817e+001
EQ: min|aij| = 2.673e-002  max|aij| = 1.000e+000  ratio = 3.741e+001
2N: min|aij| = 1.563e-002  max|aij| = 1.250e+000  ratio = 8.000e+001
Crashing...
Size of triangular part = 26
Solving LP relaxation...
      0: obj =  1.008907291e+001  infeas = 2.123e+002 (2)
*    33: obj =  4.500000000e+001  infeas = 4.764e-015 (1)
OPTIMAL SOLUTION FOUND
Integer optimization begins...
Gomory's cuts enabled
MIR cuts enabled
Cover cuts enabled
Clique cuts enabled
Creating the conflict graph...
The conflict graph has 2*74 vertices and 784 edges
+    33: mip =     not found yet >=              -inf        (1; 0)
+ 20878: >>>>>  4.500000000e+001 >=  4.500000000e+001   0.0% (34; 1002)
+ 20878: mip =  4.500000000e+001 >=     tree is empty   0.0% (0; 1181)
INTEGER OPTIMAL SOLUTION FOUND
Time used:   3.6 secs
Memory used: 1.0 Mb (1014952 bytes)
Display statement at line 17
y[G] = 8
y[E] = 4.00000000000001
y[O] = 5.99999999999998
y[R] = 4.99999999999998
y[I] = 0
y[A] = 2.00000000000001
y[N] = 3
y[V] = 1
y[M] = 9
y[T] = 7
Model has been successfully processed
>Exit code: 0

This solution is correct. At this moment I have no explanation why the GAMS version is doing worse. There may be issues with compilers, GAMS handling of the objective or other possible reasons. There is a difference in the log indicating the scaling is done differently. I believe the standalone version is 4.34 (Dec 2008) compared to 4.32 (built Feb 12 2009) for the GAMS version, so the GAMS system was little bit outdated. Also note that I used a 64 bit GAMS system and a 32 bit GLPSOL system. It remains very treacherous to just compare results from a numerical algorithm, but in this case we expect things to be more equal than they are now, and it may be possible we stumbled upon a bug.

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.

Wednesday, April 22, 2009

Lagrangian Relaxation with GAMS

> You had this example about Lagrangian Relaxation with GAMS?

Here you go: lagrange.pdf.

Monday, April 20, 2009

MS Chart for .NET 3.5

Free chart component for .NET:

Note: there is a download for the component itself and a separate download for the Visual Studio integration.

TSP Example

The graph in this model is somewhat misleading as the coordinates are in DD.MM format where DD are degrees and MM are minutes. A more realistic graph is created by using decimal degrees: DD + MM/60. Note: the optimal tour is calculated correctly as the distances were ok. Just the display of cities was slightly inaccurate. See below: