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;

You really should get advise from your professor. He/she gets paid to help you with this.

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:



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

z can be continuous with z ∈ [0,1]. See also http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-binary-variables.html.

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.

See also [link].

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

GAMS/Gurobi link

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
/;