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, 1998G 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.009GAMS/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 TERMINATEDSolving 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.