## Thursday, April 30, 2009

### If-then-else in AMPL/Cplex

I´m making my final project working with AMPL and nobody in my university can help me. […] The question is the following: I must take a decision over a variable, DECISION[n]. This variable can be a number between 0 to 64. I want to obtain another variable, RESULT[n], that if the first variable is Zero, the second one is zero too. In the other cases, must be 1. In C language it could be: if DECISION[n] =0 then RESULT[n]=0 (else RESULT[n]=1). i´m working with CPLEX.

E.g. using Cplex indicator constraints:

EQ{n in N}: RESULT[n]=0 ==> DECISION[n]=0 else DECISION[n]>=0.0001;

## 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.

## 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:

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:

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:

In this posting we developed an equation

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:

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:

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:

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

Here you go: lagrange.pdf.

## Monday, April 20, 2009

### MS Chart for .NET 3.5

Free chart component for .NET:

### 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:

## Friday, April 17, 2009

### How to linearize

> I have this expression x1⋅x2⋅x3, where all variables are binary. How can I linearize this?

The nonlinear equation z=x1⋅x2⋅x3 with xi ∈ {0,1} can be linearized as:

z ≤ x1
z ≤ x2
z ≤ x3
z ≥ x1+x2+x3−2

## Thursday, April 16, 2009

### GAMS Solaris 64 bit problem

> When I run GAMS on a Solaris 64 bit SPARC system I encounter the following problem:
>
> ld.so.1: gmsco3ux.out: fatal: libfui.so.1: open failed: No such file or directory

This means you are missing Fortran run time support. A client of mine had the same problem when I sold him a GAMS system for Solaris 64-bit SPARC. In this case you try to run the CONOPT solver, which is written in Fortran. It is dynamically linked and therefore needs some run time libraries. If you don't have a Fortran compiler installed on your system you may get this error. Software suppliers are allowed to provide these files. The GAMS web site actually does provide a download for the 32-bit runtime here. Unfortunately for 64 bit no such download is provided. You can contact GAMS support to request the 64 bit Fortran run time libraries.

## Monday, April 13, 2009

### Speeding Up Gams

The following calculations took more than 4 hours:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =
sum(phzone\$CountryLGSxTRCrops(faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr), 1);

scalar max;
max = smax((faoqcountry,faocrop,grid,irrg_key,lgsxtr), count(faoqcountry,faocrop,grid,irrg_key, lgsxtr));
display max;

The parameter CountryLGSxTRCrops is very large (but sparse): 2.5 million elements. The reason for this fragment to be so slow is two-fold:

• In the assignment to parameter count we are doing things out of order. GAMS prefers to loop over sets at the end of the index list, but here we loop over an index in the middle (phzone). Sometimes GAMS is able to reorder things automatically to do things faster, but in this case apparently not.
• smax is executed dense, as zeroes can be significant. In GAMS zero and “does not exist” is the same. This can be exploited when sparse processing is used. Here GAMS is conservative and reverts to dense processing as the implicit zeros may be important in calculating the correct SMAX value.

This a direct rewrite into an explicit loop that executes in 30 seconds:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
scalar mx /0/;
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr)=0;
loop((faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr)\$CountryLGSxTRCrops(faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr),
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =  count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) + 1;
mx = max(mx,count(faoqcountry,faocrop,grid,irrg_key,lgsxtr));
);

Explicit loops are expensive in GAMS but still much cheaper compared to dense processing. Even faster is:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
parameter ReorderedCrops(faoqcountry,faocrop,grid,irrg_key,lgsxtr,phzone);
option ReorderedCrops < CountryLGSxTRCrops;

count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =
sum(phzone\$reorderedCrops(faoqcountry,faocrop,grid,irrg_key,lgsxtr,phzone), 1);

scalar max;
max = smax((faoqcountry,faocrop,grid,irrg_key,lgsxtr)\$count(faoqcountry,faocrop,grid,irrg_key,lgsxtr),
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr));
display max;

Here we use the (undocumented)  < option to reorder the big parameter such that the phzone index is last. This is the representation that GAMS likes most when looping over phzone. Now the calculation of count is very fast. Next we force the calculation of max to be performed over the nonzero elements in parameter count only. This formulation takes less than 5 seconds. We use some more memory by duplicating CountryLGSxTRCrops.

## Sunday, April 12, 2009

### SQL2GMS tricks

Using SQL2GMS to read a large set of CSV files has the big advantage to be able to repair things on the fly. We had a few NULL values in one of the CSV files and I fixed this by skipping these records:

Q39=select 'faotqcntry_'&intFAOLVSTQ_CNTRY,'faolivestp_'&intFAOLVST2_CODE,[dblFAOLVSTPRICE_\$S/MT] from LivestockPrices.csv where not isnull([dblFAOLVSTPRICE_\$S/MT])

Then I got the request to use a value of 99999 for those records. This can be easily done by:

Q39=select 'faotqcntry_'&intFAOLVSTQ_CNTRY,'faolivestp_'&intFAOLVST2_CODE,iif(isnull([dblFAOLVSTPRICE_\$S/MT]),99999,[dblFAOLVSTPRICE_\$S/MT]) from LivestockPrices.csv

Note: we generate more meaningful set elements than 1, 2, 3,… by prefixing them with an ID string so they become faotqcntry_1, faotqcntry_2, faotqcntry_3,… This helps when debugging GAMS models.

## Wednesday, April 8, 2009

### Unjust benchmark: TSP MIP/Gurobi vs GA/Octave

I was asked to have a look at the Genetic Algorithm from this site: http://www.mathworks.com/matlabcentral/fileexchange/13680. Just to get an idea I benchmark here the TSP model eil51 from TSPLIB against GAMS/Gurobi. I used Octave instead of Matlab. The solution TSP tour is quite different.

 MIP/Gurobi GA/Octave Source eil51.gms tsp_ga.m Language GAMS Matlab Solver Gurobi Octave Options threads 4 cuts 2 heuristics 0.2 defaults Objective 426.00 (Proven optimal) 437.23 Time (seconds) 214.255 1939.2

Warning: to a large extent this is not a fair comparison.
TSP Tour from GAMS/Gurobi:

TSP Tour from TSP_GA:

## Monday, April 6, 2009

### Binary QP

I have been trying to solve an Binary Integer Program using AMPL/CPLEX
with the following objective function

minimize Total_Cost:
sum{i in flights,k in gates} In[i,k] * Pnty * (Y[i,k]-X[i,k])^2 + sum
{i in flights, k in gates} ct*FTG[i,k]*(Y[i,k]-X[i,k])

where X[i,k] is the only decision variable.

I tried to solve it, but it gave me the following error:

CPLEX 11.2.0: 5 diagonal QP coefficients of the wrong
Variable Diagonal
_svar[8] -5
_svar[4] -2.5
_svar[3] -24.5
Diagonal QP Hessian has elements of the wrong sign.

When I tried to solve the problem using Minos (trial version) solver
for miniature problem, it solves it but the integrality of some
variables are ignored.

I am trying to linearize the above objective function and trying to
solve it through the CPLEX, but I am not successful yet. Can you
suggest for any solution or approach to solve it.

Assuming Y is a parameter and X is a binary variable (and all other identifiers are parameters), you can replace

(Y[i,k]-X[i,k])^2 = Y[i,k]-2*Y[i,k]*X[i,k]+X[i,k]^2 = Y[i,k]-2*Y[i,k]*X[i,k]+X[i,k]

which is linear in X. Note that x2=x when x is a binary variable.

### 24 cores

This [link] is to practice your French. Microsoft Solver Foundation on a machine with 24 cores to schedule the TechDays 2009 (Paris) conference.

## Sunday, April 5, 2009

During a demo of GAMS/Gurobi we found a minor error. It shows NOPT's for equations in optimal MIP models. No show-stopper: you should just ignore these messages for now. This is with model magic.gms from the model library:

S O L V E      S U M M A R Y

MODEL   william             OBJECTIVE  cost
TYPE    MIP                 DIRECTION  MINIMIZE
SOLVER  GUROBI              FROM LINE  81

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

---- EQU maxu  maximum generation level (1000mw)

LOWER          LEVEL          UPPER         MARGINAL

type-1.12pm-6am        -INF          -13.8000          .              .
type-1.6am-9am         -INF           -8.0000          .              .
type-1.9am-3pm         -INF          -13.0000          .              .
type-1.3pm-6pm         -INF           -2.7500          .              .
type-1.6pm-12pm        -INF          -12.7500          .              .
type-2.12pm-6am        -INF           -0.4500          .              .
type-2.6am-9am         -INF             .              .         -2100.0000  NOPT
type-2.9am-3pm         -INF             .              .         -4200.0000  NOPT
type-2.3pm-6pm         -INF             .              .         -2100.0000  NOPT
type-2.6pm-12pm        -INF             .              .         -4200.0000  NOPT
type-3.12pm-6am        -INF             .              .              .
type-3.6am-9am         -INF             .              .              .
type-3.9am-3pm         -INF             .              .              .
type-3.3pm-6pm         -INF           -5.0000          .              .
type-3.6pm-12pm        -INF             .              .              .

**** REPORT SUMMARY :        4     NONOPT ( NOPT)
Note: these NOPTS are used to mark rows and columns with the wrong sign for the marginals in LP's. In this case the signs for these row marginals (i.e. duals) look just fine. In GAMS, after a MIP is solved, an LP is formed by fixing the integer variables. The marginals reported in the listing file is for this fixed LP problem. In general, NOPT flags should never appear in a model that is declared optimal: they are reserved for models that are intermediate non-optimal, a status that can happen e.g. when hitting a time or iteration limit before the model was optimal. If we analyze this a little bit further, we note that for an optimal minimization problem all marginals corresponding to non-basic rows (and columns) at upper bound should be negative or EPS. This is the case for these rows, and the NOPT flag is therefore inaccurate. Note that a row or column can be basic while at bound (this is also known as degeneracy). We see a few of those degenerate rows also (they are reported correctly).

Update: this is fixed.

## Friday, April 3, 2009

### glpk/ampl to GAMS conversion

> I don’t have glpk2gams, how do I convert an ampl or glpk model to gams?

1. Generate an MPS file from AMPL or GLPSOL:
1. Use the command line flag –wmps in glpsol
2. In AMPL you can use
3. ampl: option auxfiles rc;
ampl: write mxxx;

Here xxx can be replaced by a better name

2. Use mps2gms to translate this MPS file into a GAMS file. Note: mps2gms is part of the GAMS distribution so you already have it.

## Thursday, April 2, 2009

### GAMS: set up a one-to-one mapping

In a model I was working on yesterday, I needed something like:

set i 'jobs' /start1*start10, job1*job100, end1*end10/;
set r 'resources' /r1*r10/;
set startmap(i,r) /
start1.r1
start2.r2
start3.r3
start4.r4
start5.r5
start6.r6
start7.r7
start8.r8
start9.r9
start10.r10
/;
set endmap(i,r) /
end1.r1
end2.r2
end3.r3
end4.r4
end5.r5
end6.r6
end7.r7
end8.r8
end9.r9
end10.r10
/;

This “diagonal mapping” can be written more compactly as using some new GAMS syntax:

set i 'jobs' /start1*start10, job1*job100, end1*end10/;
set r 'resources' /r1*r10/;
set startmap(i,r) /
start1*start10:r1*r10
/;
set endmap(i,r) /
end1*end10:r1*r10
/;