Monday, August 1, 2011

More on network models

I received some additional questions on solving some very large network models. Basically we used the formulation presented here: http://yetanothermathprogrammingconsultant.blogspot.com/2011/06/network-formulation.html. A slightly more optimized version is shown below:

$ontext

 
large sparse max flow network example


 
Performance

     
nodes  arcs    gams       solver
                    
execution  time

     
1k     10k     0.031      0.045
     
2k     20k     0.047      0.112
     
5k     50k     0.156      0.286
     
10k   100k     0.343      1.125

 
Erwin Kalvelagen,
 
Amsterdam Optimization,
 
erwin@amsterdamoptimization.com

$offtext

$set n 10000


option limrow=0, limcol=0, solprint=off;

sets

  i
  source(i)
  sink(i)
;

alias(i,j);

parameter capacity(i,j) 'capacity, but also defines network'
;

$gdxin network%n%
$loaddc i source sink capacity


abort$(card(source)<>1) "We want one source node"
;
abort$(card(sink)<>1) "We want one sink node"
;

set arcs(i,j) 'we have an arc if capacity<>0'
;
arcs(i,j)$capacity(i,j) =
yes
;

scalar narcs 'number of arcs'
;
narcs =
card
(arcs);
display
narcs;

parameter
rhs(i);
rhs(source) = -1;
rhs(sink) = 1;


variables

  x(i,j)
'flow along arcs'
  f     
'total flow'
;
positive variables x;
x.up(arcs) = capacity(arcs);


equations

   flowbal(i)
'flow balance' ;

flowbal(i)..
  
sum(arcs(j,i), x(j,i)) - sum
(arcs(i,j), x(i,j)) =e= f*rhs(i);

model m /flowbal/
;
solve
m maximizing f using lp;

display
f.l;

The performance of GAMS looks very good. GAMS execution time scales slightly better than the LP solver used to solve the problem:

image

GAMS almost scales linearly with the number of nonzero elements to be generated. In this case that means linearly with the number of nodes (because the number of nonzeroes is equal to 10*nodes in these models).

The main difference with the earlier model is that we replaced:

x.up(i,j) = capacity(i,j);

by

x.up(arcs) = capacity(arcs);

This actually saves some time and memory. In the first case we set many x.up’s to zero even though they are not used in the model. These variables with non-default bounds will need to be created and use up memory. In the second case, we only touch variables that are really used.

This small change makes a lot of difference for larger models. For n=5k we see GAMS execution time decrease from 3.775 seconds to 0.156 seconds. If things get large one needs to pay attention to detail!

AMPL Version

A direct translation to AMPL can look like:

set nodes;
set source within nodes;
set sink within nodes;
param capacity {(i,j) in (nodes cross nodes)}, default 0;

set arcs within (nodes cross nodes) := {i in nodes, j in nodes: capacity[i,j]}; # version 1 : slow calc of arcs
# set arcs := {i in nodes, j in nodes: capacity[i,j]};                                   # version 2 : slow calc of x
# set arcs within (nodes cross nodes);                                                    # version 3 : fast, read from data file

param rhs{i in nodes} := if i in source then -1 else if i in sink then +1;

var x {(i,j) in arcs} >= 0, <= capacity[i,j];
var f;

maximize obj:f;

subject to flowbal{i in nodes}:
   sum{(j,i) in arcs} x[j,i] - sum{(i,j) in arcs} x[i,j] = f*rhs[i];

This formulation is not handled very efficiently by AMPL: the timings for n=5k are

##genmod times:
##seq      seconds    cum. sec.    mem. inc.  name
## 77            0            0            0  derstage
## 81            0            0            0  sstatus
## 95            0            0            0  nodes
## 96            0            0            0  source
## 97            0            0            0  sink
## 98            0            0            0  capacity
## 99      10.4521      10.4521   1146012664  arcs
## 100            0      10.4521            0  rhs
## 101    0.0156001      10.4677        65544  x
## 103            0      10.4677            0  f
## 105            0      10.4677            0  obj
## 107    0.0624004      10.5301      4148296  flowbal

Especially the calculation of the set arcs is expensive.

Note that if we use an alternative for arcs:

set arcs := {i in nodes, j in nodes: capacity[i,j]};

the burden is moved to x:

##genmod times:
##seq      seconds    cum. sec.    mem. inc.  name
## 77            0            0            0  derstage
## 81            0            0            0  sstatus
## 95            0            0            0  nodes
## 96            0            0            0  source
## 97            0            0            0  sink
## 98            0            0            0  capacity
## 99            0            0            0  arcs
## 100            0            0            0  rhs
## 101      10.5301      10.5301   1147323392  x
## 103            0      10.5301            0  f
## 105            0      10.5301            0  obj
## 107    0.0780005      10.6081      4147928  flowbal

When we run a few different instances we see we do something causing non-linear scaling:

image

It is clear this is not a good approach. We need to do something about this set arcs. In AMPL we don’t really need to calculate this set, we can directly populate this from the data by replacing in the data file:

param: capacity :=
n1 n32 87
n1 n37 44
n1 n124 38

…

by

param: arcs : capacity :=
n1 n32 87
n1 n37 44
n1 n124 38
…

Now we get the set arcs basically for free, and we get linear scaling:

image

The total AMPL generation time vs. LP solution time is here:

image

This is very close to the GAMS timings given that this was on a different machine.

GLPK

The open source tool GLPK can handle many AMPL linear models. However, in some cases it is substantial slower. This model demonstrates this behavior. We run the model as:

C:\projects\tmp>timeit \Projects\glpk\glpk\bin\glpsol.exe -m network.mod -d network1000.dat
Reading model section from network.mod...
network.mod:24: warning: unexpected end of file; missing end statement inserted
24 lines were read
Reading data section from network1000.dat...
network1000.dat:10108: warning: unexpected end of file; missing end statement inserted
10108 lines were read
Generating obj...
Generating flowbal...
Model has been successfully generated
lpx_simplex: original LP has 1001 rows, 10001 columns, 20003 non-zeros
lpx_simplex: presolved LP has 1000 rows, 10001 columns, 20002 non-zeros
lpx_adv_basis: size of triangular part = 999
*     0:   objval =  0.000000000e+000   infeas =  0.000000000e+000 (1)
*   171:   objval =  4.780000000e+002   infeas =  0.000000000e+000 (1)
OPTIMAL SOLUTION FOUND
Time used:   0.1 secs
Memory used: 13.4M (14073536 bytes)

Version Number:   Windows NT 6.0 (Build 6002)
Exit Time:        0:32 am, Thursday, August 4 2011
Elapsed Time:     0:00:05.563
Process Time:     0:00:05.600
System Calls:     160954
Context Switches: 89933
Page Faults:      20453
Bytes Read:       248437
Bytes Written:    87315
Bytes Other:      130434

C:\projects\tmp>

Besides being much slower than AMPL, the performance of the model generation phase is showing more than linear scaling:

image

Model generation is here calculated as total elapsed time minus time used for the LP solver.

My conclusions:

  1. AMPL and GAMS are capable of generating simple network models faster than a good LP solver can solve them
  2. But a simple mistake and performance goes down the drain
  3. Direct translations between GAMS and AMPL models are not always a good idea.

Thursday, July 28, 2011

Scheduling of TV Advertisement: Theory vs Practice

I am working on a model for scheduling TV advertisements. Looking at the paper: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.292&rep=rep1&type=pdf this is easy to model (but difficult to solve). My model seems at least 10 times as complicated than what this academic study suggests. We needed more complexity to deal with:

  1. The number of spots to broadcast is not known in advance (depends on expected ratings of a break)
  2. There are complicated issues with spread over time
  3. There is a rather complex issue with the quality of the remaining unused space (this needs to be valued with respect to expected ratings and optimized)
  4. Different contracts have different requirements, and these requirements are not always easily modeled
  5. The model is huge (> 150K binary variables) for a monthly schedule
  6. The input data (from a database) is large and complex (this is sometimes called “data-intensive”)
  7. Numerous other things

Thursday, July 21, 2011

Pyomo

For a thoughtful comment on the design of Python based modeling system see: http://yetanothermathprogrammingconsultant.blogspot.com/2011/04/performance-of-python-based-modeling.html?showComment=1310240113570#c3299569615459362578.

My assessment was clearly not convincing. To summarize: essentially the modeling framework Pyomo is a scalar system. This compares somewhat to the difference between AMPL and GLPK (http://www.gnu.org/software/glpk/glpk.html). Both systems are written in C so we have taken out the language issue (i.e. Python is slow compared to C). The modeling system in GLPK (called Mathprog) is dense (except for the final generated matrix which is of course sparse). This actually work in many cases just fine: many sub-matrices used in the model are often small enough that dense handling is not really a problem. I have seen this in my practice: the majority of the models implemented in GLPK/Mathprog just work fine. Of course there are some models (the big sparse suckers) that bring GLPK/Mathprog on its knees while AMPL will happily generate the model very quickly. If I see this I generally suggest some reformulation or more often just give the advise: buy AMPL.

I suspect this can also happen with Pyomo. For some models a simple dense/scalar structure where we iterate over the carthesian product of indices will not perform. Of course with a system like Pyomo the modeler has more possibilities to implement better data-structures so we can skip zeros quickly and efficiently. But that means we have moved the burden from modeling system to the modeler.

Patent Application

One of my clients is applying for a patent on a model I helped them with. I am mentioned as co-inventor. That is mostly just an ego-booster I suspect, as this was work-for-hire. From what I understand the patent application is quite narrowly defined, so hopefully I am still able to use my bag of LP tricks for other clients.

Tuesday, July 5, 2011

Min sum or min num

In an advertisement model we may have to deal with underdelivery: the delivery of less impressions, visitors, GRPs than contractually agreed upon. So what to do. We can minimize the SUM of this quantity with the possible effect of many orders with small underdelivery numbers (spreading the pain). When minimizing the MAX this effect is even guaranteed. Or we can minimize the NUMber of orders with underdelivery. 

image

In practice, a good approach is to minimize a linear combination of both the SUM and the NUM objectives.

GAMS: GDXDUMP and Equations

There is a fairly deep issue when comparing GDX files if one of the GDX files has the following properties:

  1. It is produced by GDXDUMP
  2. It contains Equation records

The following illustrates the problem:

C:\projects\tmp>gamslib indus89
Copy ASCII: indus89.gms

C:\projects\tmp>gams indus89 lo=0 gdx=1

C:\projects\tmp>gdxdump 1.gdx output=2.gms
*  GDX dump of 1.gdx
*  Library in use : C:\PROGRA~1\GAMS23.7
*  Library version: GDX Library      BETA  8Jun11 23.7.0 WEX 25828.25830 WEI x86_64/MS Windows
*  File version   : GDX Library      BETA  8Jun11 23.7.0 WEX 25828.25830 WEI x86_64/MS Windows
*  Producer       : GAMS Base Module BETA  8Jun11 23.7.0 WEX 25828.25830 WEI x86_64/MS Windows
*  File format    :    7
*  Compression    :    0
*  Symbols        :  281
*  Unique Elements:  245

C:\projects\tmp>gams 2.gms lo=0 gdx=2

C:\projects\tmp>gdxdiff 1.gdx 2.gdx eps=1.0e-6
GDXDIFF          BETA  8Jun11 23.7.0 WEX 25828.25830 WEI x86_64/MS Windows
File1 : 1.gdx
File2 : 2.gdx
Summary of differences:
   bdraft   Data are different
   brepco   Data are different
  ccombal   Data are different
  consbal   Data are different
     cost   Data are different
   demnat   Data are different
divcnlsea   Data are different
   divsea   Data are different
   fodder   Data are different
   grnfdr   Data are different
   laborc   Data are different
     nbal   Data are different
  nwfpalc   Data are different
     objn   Data are different
    objnn   Data are different
     objz   Data are different
    objzn   Data are different
  protein   Data are different
   prseaw   Data are different
  qcombal   Data are different
  subirrc   Data are different
   tdraft   Data are different
watalcpro   Data are different
watalcsea   Data are different
  watalcz   Data are different
waterbaln   Data are different
Output: diffile.gdx
GDXDiff finished

C:\projects\tmp>gdxdump diffile.gdx Symbols
*  GDX dump of diffile.gdx
*  Library in use : C:\PROGRA~1\GAMS23.7
*  Library version: GDX Library      BETA  8Jun11 23.7.0 WEX 25828.25830 WEI x86_64/MS Windows
*  File version   : GDX Library      BETA  8Jun11 23.7.0 WEX 25828.25830 WEI x86_64/MS Windows
*  Producer       : GDXDIFF
*  File format    :    7
*  Compression    :    0
*  Symbols        :   26
*  Unique Elements:  249
   Symbol    Dim Type  Explanatory text
1 bdraft      4  Equ  Differences
2 brepco      3  Equ  Differences
3 ccombal     4  Equ  Differences
4 consbal     4  Equ  Differences
5 cost        3  Equ  Differences
6 demnat      3  Equ  Differences
7 divcnlsea   3  Equ  Differences
8 divsea      2  Equ  Differences
9 fodder      4  Equ  Differences
10 grnfdr      4  Equ  Differences
11 laborc      4  Equ  Differences
12 nbal        3  Equ  Differences
13 nwfpalc     2  Equ  Differences
14 objn        1  Equ  Differences
15 objnn       1  Equ  Differences
16 objz        1  Equ  Differences
17 objzn       1  Equ  Differences
18 protein     4  Equ  Differences
19 prseaw      3  Equ  Differences
20 qcombal     4  Equ  Differences
21 subirrc     4  Equ  Differences
22 tdraft      4  Equ  Differences
23 watalcpro   3  Equ  Differences
24 watalcsea   3  Equ  Differences
25 watalcz     4  Equ  Differences
26 waterbaln   4  Equ  Differences

C:\projects\tmp>

Apart from a small bug related to $ONEMPTY, it turns out that the way GDXDUMP exports equations records, it looses information about the equation type. This leads to the above results. Luckily in practice this will not be an issue, as equations are not often part of such an exercise.

Background.

During a course I taught about GAMS I was asked if it is possible to edit GDX files. The answer is: GDXDUMP will give you a GAMS representation of a GDX file which you can edit. Using the call:

GAMS gmsfile gdx=newgdxfile

you can make a new GDX file with the changes incorporated. The tool GDXDIFF will allow you to monitor the changes.

Although this is correct, I chose the INDUS89 example as illustration. That was unfortunate, as that failed because of the exotic Equation problem.

Excel issue

I was doing some reporting on the results of a large, complex MIP model. I made the unfortunate error of having some spurious blank cells. With the =SUM() function, skipping such cells is identical as considering them as zero, so there is no problem. With a =MIN() function however this is not the case. Even =MINA() will not help there:

image

It is sometimes good to work with different people on a project so that one of your collaborators can catch this before submitting results to the customer!

In this case the result was beneficial for us. With this, our estimated improvement in monthly revenue when using a MIP based model, increased from 250K euros to 300K euros.

Friday, July 1, 2011

How bumblebees tackle the traveling salesman problem

http://www.eurekalert.org/pub_releases/2011-06/qmuo-hbt062811.php

Sorting in GAMS

This question came up during a workshop I taught at a large power company:

I want to generate a number of load-duration and price-duration curves, but it takes too long.

Here is a small example demonstrating some formulations for sorting a 1d parameter:

option profile=1;
set i /i1*i10000/
;
parameter
p(i);
p(i) = uniform(1,100);

display
p;

execute_unload "p"
,p;
execute "gdxrank p.gdx p2.gdx"
;

parameter
pindex(i);
execute_load "p2.gdx"
,pindex=p;
display
pindex;


* slow

parameter p2(i);
alias
(i,j);
p2(i)=
sum(j$(ord
(i)=pindex(j)),p(j));
display
p2;

* almost as slow

parameter p3(i);
alias
(i,j);
loop((i,j)$(ord
(i)=pindex(j)),
   p3(i)=p(j);
);

display
p3;

* fast using trick:

parameter p4(i);
p4(i+(pindex(i)-
ord
(i))) = p(i);
display
p4;

The timings are:

Method Time (seconds)
SUM (assignment to p2) 28.407
LOOP (assignment to p3) 22.230
TRICK (assignment to p4) 0.265