Saturday, December 29, 2012

Option optcr

Building and Solving Mathematical Programming Models in Engineering and Science (Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts)

Some of the MIP models in http://www.amazon.com/gp/product/0471150436 (these models are written in GAMS) seem to be missing OPTION OPTCR=0, so some non-optimal solutions are presented. By default GAMS uses a 10% relative gap criterion.

E.g.:

image

But I see:

image

Next month teaching a GAMS course in Germany. I’ll be sure to mention this (and how to spot this in the listing file).

GAMS is quite unique in having a default with nonzero gap optimality criteria. Other systems typically like to solve to optimality unless instructed otherwise.

PS. See comment for another system with a nonzero gap tolerance.

Cplex 12.5 vs 12.4 Out of memory handling

With Cplex 12.4:

image

Besides Cplex crashing here, we also see poor error recovery and exception handling by the calling problem. It is good to see that with Cplex 12.5 we see much better behavior:

image

This is with 32 bit software. We could use a 64-bit version instead to allow larger memories, but instead I opted for reformulating the model and making it smaller so I don’t leave this work to the presolver.

Thursday, December 27, 2012

Better element names

Sometimes I see GAMS code like:

SETS
  K
index of periods of time  /1*4/
  J
index of generators /1*3/

In general I don’t like plain numbers as set elements. I prefer something like:

SETS
  K
index of periods of time  /t1*t4/
  J
index of generators
/gen1*gen3/

This often improves the readability of the output. This is especially the case in the GDX viewer where we can swap rows and columns. Here is a some sample output:

Numerical

Id+Numerical

SETS
  K
index of periods of time  /1*4/
  J
index of generators /1*3/

SETS
  K
index of periods of time  /t1*t4/
  J
index of generators
/gen1*gen3/

----    110 VARIABLE p.L  output power of generator j at period k

 

            2           3           4

 

1     150.000     350.000     260.000

2                 100.000

3                  50.000     140.000

 

----    110 VARIABLE p.L  output power of generator j at period k

 

              t2          t3          t4

 

gen1     150.000     350.000     260.000

gen2                 100.000

gen3                  50.000     140.000

clip_image001

clip_image002

For this small example this is not terribly convincing, but for larger, high-dimensional data it surely helps in understanding output quickly.

Wednesday, December 26, 2012

Lazy Constraints and GAMS

> Can I use lazy constraints in my GAMS model?

Sorry, the answer is “no”. Lazy constraints are supported by the Cplex and Gurobi solvers. Other modeling systems (AMPL, AIMMS) support lazy constraints, but GAMS does not. There have been some attempts to implement this in GAMS but they were not useable in practice (using scalar models as intermediary). This was more for testing I think. 

Here is a discussion on lazy constraints: http://orinanobworld.blogspot.com/2012/08/user-cuts-versus-lazy-constraints.html

Long SQL queries

Over the holidays I ran a small script against SQL Server. Basically spatial queries. I ran the script in pieces to make sure I did not have to rerun the whole thing if something went wrong. That was a good idea: some of the queries took 6 or 7 hours or more. With memories getting larger, may be tasks like this are better done in-core than using a database. 

Tuesday, December 25, 2012

No error message

I have posted regularly about poor error messages from software. In general I find that the importance of formulating meaningful  error messages is underestimated by developers, leading to some level of frustration by users. Of course, issuing no error message at all is even easier for a developer (and even more difficult for users to deal with).

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

     MODEL   LANDfeas            OBJECTIVE  totallow
     TYPE    LP                  DIRECTION  MINIMIZE
     SOLVER  MOSEK               FROM LINE  165233

**** SOLVER STATUS     4 Terminated By Solver     
**** MODEL STATUS      14 No Solution Returned    
**** OBJECTIVE VALUE                0.0000

RESOURCE USAGE, LIMIT         20.311    900000.000
ITERATION COUNT, LIMIT         0        900000

    Copyright (C)   MOSEK ApS, Fruebjergvej 3, Box 16
                    DK-2100 Copenhagen, Denmark
                   
http://www.mosek.com
GAMS/MOSEK Extended license detected

No solution returned

I received a listing file with this fragment. I have no idea how I can help the user as I have no idea what went wrong. Now I need to get the model and data files and see if I can reproduce the problem. We see some strange trade-offs at work here: the developers saved some time by not providing code to print a readable message, and as a result (multiple) users have to spent extra time and effort in supporting and maintaining an application.

Sunday, December 23, 2012

MySQL and union

For a GAMS model I need to read data from a MySQL database. Some of the data can be easily imported into a GAMS gdx file using an SQL UNION construct:

SELECT id,xxx FROM tableA
UNION
SELECT `obj’,xxx FROM tableB

Now the interesting part is that ids are INTs. So basically the first column looks like

1 xxx
2 xxx
3 xxx
obj xxx

On my MySQL installation, the first column gets a correct type of VARCHAR(11) as the integer type was INT(10) – the extra char is for the sign I think (although we have only positive ids). However when the client was testing this, it went terribly wrong. Turned out that their MySQL version created a type of VARBINARY(11), a BLOB like type. This created havoc downstream. Luckily I was using $LOADDC in GAMS to read the GDX file produced by SQL2GMS, so we caught a domain error.

No doubt this was bad SQL produced by me. I tried something similar with SQL Server, and SQL Server complains directly about not being able to put a string into integer type. The fix: using a cast as in:

SELECT cast(id as char(10)),xxx FROM tableA
UNION
SELECT `obj’,xxx FROM tableB

A little bit troublesome that different versions of MySQL behave differently on this.

Monday, December 10, 2012

Fun with SQL Server and Quotes

I am trying to help with some SQL and we were really having fun today! A complicated dynamic SQL statement involving cursors was failing with the illuminating error: syntax error near ‘a’. When working with SQL Server it is not easy to parameterize for table names, hence the dynamic SQL instead of normal static SQL. Well, the problem turned out to be a quote in a name in the database. This was causing the generated SQL statement to be incorrect.

The best would be to fix all SQL queries, but that was not feasible right now. Instead we tried to replace all quotes by two quotes. (This was a staging table anyway, so no problem in polluting this table: it is not used by others). Sometimes you have to be pragmatic and take a short cut. The static SQL to do: replace a single quote by two single quotes looks already funky but the dynamic SQL version is really horrific:

image

Sunday, December 9, 2012

Low cost parallel computing

There are a few interesting projects going on:

Tuesday, December 4, 2012

History of LP/MIP

http://www.gurobi.com/pdf/a-brief-history-of-optimization.pdf. I am always fascinated by these stories… Every time when looking back I am impressed by the size and complexity of models we can solve these days on our laptops compared to what we could solve on mainframes and supercomputers a few decades ago (and with much more convenience).

PS. You also may want to look at the whole issue at: http://www.zib.de/groetschel/publications/OptimizationStories.pdf. More interesting historic details about optimization.

Monday, December 3, 2012

GAMS question

A question was posed at https://groups.google.com/forum/?fromgroups=#!topic/gamsworld/0qYEcVd9x8I:

set    j  /n1*n5/,
       k  /n5/;

variable pi(j);

scalar ab;
ab=sum (k, sum(j$(ord(j)=ord(k)), pi.L(j) )  );

The user want to pick the 5-th element. This implementation clearly does not work. Some valid answers were given:

Ab= sum(j$sameas(j,“n5“), pi.L(j) )  ;
or
ab=sum (k, sum(j$sameas(j,k), pi.L(j) )  );

However I would say this is not really good modeling practice. The best approach IMHO would be:

set   j     /n1*n5/,
      k(j)  /n5/;

variable pi(j);

scalar ab;
ab=sum(k, pi.L(k));

There are two advantages of this approach:

  • The set k is now domain checked against j; no elements outside j can be specified here without alarm bells going off. Domain checking is very important when developing large, complex models.
  • We can substantially simplify the assignment to ab and increase readability. The expression conveys immediately what it is supposed to do.

Setting up proper subsets will help you down the road.

Parallel LP

> Can I run LP solver in parallel on multiple cores? What about all running on GPUs?

All major LP solvers include a barrier/interior point code that can run on multiple cores. This includes Cplex, Gurobi and Mosek. Do not expect perfect scaling but some speed-up is possible. Sparse simplex algorithms are not ideally suited for parallelization. Some recent results were shown by John Forrest. See: ISMP 2012 talk.

About GPUs:

image

(From the same talk).

There are some other algorithms that can be implemented on GPUs. Many meta-heuristics have been ported to these devices. Here are some links:

Friday, November 30, 2012

Don’t always trust what a solver is saying

This is supposed to be a fairly difficult MIQP. However, when I run it I get quickly: INFEASIBLE:

Tried aggregator 2 times.
MIQP Presolve eliminated 1 rows and 1 columns.
MIQP Presolve modified 248 coefficients.
Aggregator did 52 substitutions.
Reduced MIQP has 125 rows, 548 columns, and 3696 nonzeros.
Reduced MIQP has 496 binaries, 0 generals, 0 SOSs, and 0 indicators.
Reduced MIQP objective Q matrix has 52 nonzeros.
Presolve time =    0.00 sec.
Probing time =    0.00 sec.
Clique table members: 654.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: none, using 1 thread.
Tried aggregator 1 time.
No QP presolve or aggregator reductions.
Presolve time =    0.00 sec.
Parallel mode: none, using 1 thread for barrier
Number of nonzeros in lower triangle of A*A' = 2360
Using Approximate Minimum Degree ordering
Total time for automatic ordering = 0.00 sec.
Summary statistics for Cholesky factor:
  Rows in Factor            = 125
  Integer space required    = 241
  Total non-zeros in factor = 3007
  Total FP ops to factor    = 90075
Itn      Primal Obj        Dual Obj  Prim Inf Upper Inf  Dual Inf         
   0  1.4172319e+007 -1.4172911e+007 2.45e+003 5.04e+002 2.77e+007
   1  1.2037355e+007  1.0009350e+007 1.78e+002 3.66e+001 2.01e+006
   2  1.1903244e+007  1.1658232e+007 2.15e+001 4.42e+000 2.43e+005
   3  1.1888608e+007  1.1839300e+007 4.33e+000 8.90e-001 4.89e+004
   4  1.1885153e+007  1.1882071e+007 2.71e-001 5.56e-002 3.05e+003
   5  1.1884924e+007  1.1884909e+007 1.33e-003 2.73e-004 1.50e+001
Barrier time =    0.05 sec.

Total time on 1 threads =    0.05 sec.
Root relaxation solution time =    0.05 sec.

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

      0     0    infeasible                                          5        

Root node processing (before b&c):
  Real time             =    0.05
Sequential b&c:
  Real time             =    0.00
                          -------
Total (root+branch&cut) =    0.05 sec.
MIQP status(119): integer infeasible or unbounded

Problem is integer infeasible.

**** SOLVER STATUS     1 Normal Completion        
**** MODEL STATUS      10 Integer Infeasible      
**** OBJECTIVE VALUE               NA

OK, but that sounds not really correct. I found a good integer solution with a linearized objective (same constraints). So lets fix the integer variables and see what happens:

* solve MIP
solve gms1 minimizing z using mip;
* fix integer variables

x.fx(i,t) = x.l(i,t);
* and feed it into MIQP solver
solve gms2 minimizing z using miqcp;

Now I see:

Tried aggregator 1 time.
MIQP Presolve eliminated 178 rows and 601 columns.
All rows and columns eliminated.
Presolve time =    0.00 sec.
MIQP status(101): integer optimal solution
Fixing integer variables, and solving final QP..
Tried aggregator 1 time.
QP Presolve eliminated 178 rows and 601 columns.
All rows and columns eliminated.
Presolve time =    0.00 sec.

Total time on 1 threads =    0.00 sec.
QP status(1): optimal

Proven optimal solution.

MIP Solution:     14817439.000000    (0 iterations, 0 nodes)
Final Solve:      14817439.000000    (0 iterations)

Hmm, so not integer infeasible after all. Some more experiments yielded:

  • Using MIPSTART with the MIP solution allowed the MIQP solver to proceed
  • Scaling the objective cause the MIQP to proceed also. The scaling looks like:
    objective2.. z =e= sum(t,sqr(reserves(t)))/1e6;

Don’t always trust the messages from the solvers…

Sunday, November 25, 2012

GAMS qp interface not very coherent

It is possible to write quadratically constrained problems (and QPs in particular) in GAMS and solve them efficiently with solvers like Cplex, Gurobi. But the rules are somewhat ad-hoc and incoherent.

  1. GAMS will not extract a QCP when using the notation x**2 (for a rather obscure reason; this should probably be dealt with in a smarter way)
  2. GAMS will extract a QCP with the notation x*x or sqr(x) or power(x,2)
  3. But GAMS will again not allow power(x,1) or power(x,0) so we can’t write the following:

sets
  g
'generators' /g1,g2/
  d
'degree'   /d0,d1,d2/
;

table cfuelcost(g,d)  'coefficients for fuel cost as quadratic function of output.'
     
d2        d1       d0
g1    0.00889   10.333   200
g2    0.00741   10.833   240
;


parameter v(d);
v(d) =
ord
(d)-1;

variables

   p(g)     
'output'
   fuelcost 
'total fuel cost'
;
positive variable p;

equations

   efuel   
'fuel consumption constraint'
;

* power(x,2) is allowed
* power(x,1) or power(x,0) is not recognized?????
*efuel.. fuelcost =e= sum(g, cfuelcost(g,'d2')*power(p(g),2) + cfuelcost(g,'d1')*p(g) + cfuelcost(g,'d0'));
efuel.. fuelcost =e=
sum((g,d), cfuelcost(g,
d)*power(p(g),v(d)));

model m /all/
;
option
qcp=cplex;
solve
m using qcp minimizing fuelcost;

In forbidden cases you’ll get the poorly worded error message:

*** Could not load data from file: Detected 1 general nonlinear rows in model

Friday, November 23, 2012

Smoothing

In a scheduling model for maintenance of power generators the objective is to smooth the reserves. The authors use a quadratic objective to model and the results over a year look like:

image   image

MIQP’s still solve significantly slower than MIP models, so I was looking for a linear approximation. The first thing that comes to mind is to minimize the max. This gives:

imageimage

The problem here is we get a few weeks with really low reserves. An alternative is to minimize the difference between the max and the min. This looks like:

imageimage

This actually looks quite similar to the first picture, and we have achieved our goal: make the model linear.

Wednesday, November 21, 2012

1 vs 4 threads

Running a small MIP using 1 thread:

Optimize a model with 177 rows, 549 columns and 3800 nonzeros
Presolve time: 0.02s
Presolved: 177 rows, 549 columns, 3944 nonzeros
Variable types: 1 continuous, 548 integer (496 binary)
Found heuristic solution: objective 949.0000000

Root relaxation: objective 4.990000e+02, 947 iterations, 0.03 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  499.00000    0   63  949.00000  499.00000  47.4%     -    0s
     0     0  499.00000    0   93  949.00000  499.00000  47.4%     -    0s
     0     0  499.00000    0   63  949.00000  499.00000  47.4%     -    0s
     0     2  499.00000    0   63  949.00000  499.00000  47.4%     -    0s
H  135   117                     910.0000000  499.00000  45.2%  47.3    0s
*  204    96              44     901.0000000  499.00000  44.6%  47.1    0s
  1963   797  891.00000   34   53  901.00000  578.12129  35.8%  46.1    5s
  6122  2135  822.02658   39   78  901.00000  720.84013  20.0%  34.5   10s
 12502  3801     cutoff   44       901.00000  796.00000  11.7%  27.4   15s
 21404  5575  891.00000   51   42  901.00000  862.00000  4.33%  21.7   20s
 34494  6391  891.00000   42   53  901.00000  891.00000  1.11%  17.1   25s
 52100  4430  891.00000   67   60  901.00000  891.00000  1.11%  12.8   31s
 68445  3439     cutoff   72       901.00000  891.00000  1.11%  10.7   35s
 94101  2702 infeasible   56       901.00000  891.00000  1.11%   8.9   40s
115791  2055     cutoff   55       901.00000  891.00000  1.11%   8.0   45s
141811  1165  891.00000   40   81  901.00000  891.00000  1.11%   7.2   50s
168722   165  891.00000   44   59  901.00000  891.00000  1.11%   6.7   55s

Explored 177443 nodes (1164346 simplex iterations) in 56.54 seconds
Thread count was 1 (of 8 available processors)

Optimal solution found (tolerance 0.00e+00)
Best objective 9.010000000000e+02, best bound 9.010000000000e+02, gap 0.0%
MIP status(2): Model was solved to optimality (subject to tolerances).

Now solve the same model with the same solver on the same machine using 4 threads:

Optimize a model with 177 rows, 549 columns and 3800 nonzeros
Presolve time: 0.02s
Presolved: 177 rows, 549 columns, 3944 nonzeros
Variable types: 1 continuous, 548 integer (496 binary)
Found heuristic solution: objective 949.0000000

Root relaxation: objective 4.990000e+02, 947 iterations, 0.03 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  499.00000    0   63  949.00000  499.00000  47.4%     -    0s
     0     0  499.00000    0   93  949.00000  499.00000  47.4%     -    0s
     0     0  499.00000    0   63  949.00000  499.00000  47.4%     -    0s
     0     3  499.00000    0   63  949.00000  499.00000  47.4%     -    0s
*  141    97              50     901.0000000  499.00000  44.6%  50.9    0s
  6876  3503     cutoff   64       901.00000  558.96236  38.0%  30.4    5s
 31149  7878 infeasible   58       901.00000  891.00000  1.11%  18.1   10s
 72874  4957     cutoff   63       901.00000  891.00000  1.11%  12.2   18s
 87864  4665  891.00000   55    8  901.00000  891.00000  1.11%  11.3   20s

. . . .
97009933 49020  897.00000   80   47  901.00000  897.00000  0.44%   5.8 9085s
97054184 47405 infeasible   86       901.00000  897.00000  0.44%   5.8 9090s
97103347 45433 infeasible   95       901.00000  897.00000  0.44%   5.8 9095s
97174397    76 infeasible   97       901.00000  900.91810  0.01%   5.8 9100s

Explored 97174474 nodes (560614904 simplex iterations) in 9100.06 seconds
Thread count was 4 (of 8 available processors)

Optimal solution found (tolerance 0.00e+00)
Best objective 9.010000000000e+02, best bound 9.010000000000e+02, gap 0.0%
MIP status(2): Model was solved to optimality (subject to tolerances).

Solution times go from 1 minute to 2.5 hours! Extrapolating this to 8 cores gives me the shivers….

Data errors

Even small tables can contain errors. An example data set in a published paper about scheduling maintenance of power generators looks like:

image

In GAMS I transcribed this as:

table data(i,*)
          
capacity   allowed-from  allowed-to  outage-weeks

*             MW        week           week    number of weeks
 
unit1      555          1             26           7
 
unit2      555         27             52           5
 
unit3      180          1             26           2
 
unit4      180          1             26           1
 
unit5      640         27             52           5
 
unit6      640          1             26           3
 
unit7      640          1             26           3
 
unit8      555         27             52           6
 
unit9      276          1             26          10
 
unit10     140          1             26           4
 
unit11      90          1             26           1
 
unit12      76         27             52           3
 
unit13      76          1             26           2
 
unit14      94          1             26           4
 
unit15      39          1             26           2
 
unit16     188          1             26           2
 
unit17      58         27             52           1
 
unit18      48         27             52           2
 
unit19     137         27             52           1
 
unit20     469         27             52           4
 
unit21      52          1             26           3
;

table manpower(i,t) 'Manpower required for each week'
       
week1 week2 week3 week4 week5 week6 week7 week8 week9 week10
unit1    10    10     5     5     5     5     3
unit2    10    10    10     5     5
unit3    15    15
unit4    20
unit5    10    10    10    10    10
unit6    15    15    15
unit7    15    15    15
unit8    10    10    10     5     5     5
unit9     3     2     2     2     2     2     2     2     2     3
unit10   10    10     5     5
unit11   20
unit12   10     1     5    15
unit13   15    15
unit14   10    10    10    10
unit15   15    15
unit16   15    15
unit17   20
unit18   15    15
unit19   15
unit20   10    10    10    10
unit21   10    10    10

;

Just to be sure I added some checks:

set error(i,t) 'check for mismatch between manpower and outage-weeks';
error(i,t)$(manpower(i,t)=0
and ord(t)<=data(i,'outage-weeks')) = yes
;
error(i,t)$(manpower(i,t)<>0
and ord(t)>data(i,'outage-weeks')) = yes
;
abort$card(error) "manpower mismatch"
,error;

Low and behold, we see:

----     93 manpower mismatch

----     93 SET error  check for mismatch between manpower and outage-weeks

             week4

unit12         YES

**** Exec Error at line 93: Execution halted: abort$1 'manpower mismatch'

Indeed we have for this unit a problem in the manpower allocation.