Yet Another Math Programming Consultant

Comparing some modeling languages

2012-01-23T09:13:00.001-05:00

Not an easy task, as can be seen here:

http://lists.gnu.org/archive/html/help-glpk/2012-01/msg00033.html

To really appreciate the features of a modeling system requires quite some investment. When comparing more systems, also often a “least common denominator” model is used. Some issues are more specialized, such as:

Is there support for non-linear programming
Performance characteristics for very large models
Ability to do data manipulation efficiently
Ability to write tailored algorithms if the model is not amenable to standard solution methods (see e.g. http://yetanothermathprogrammingconsultant.blogspot.com/2012/01/dinkelbachs-algorithm.html)
Ability to debug large models and view large data sets

Some of these points are more difficult to evaluate and to weigh in a comparison.

PS: see comment below for some other interesting capabilities.

network

2012-01-21T16:18:00.001-05:00

I have a question: In GAMS, is there something similar to "structure" in C ?
I need to define a network optimization problem and I need to define: A node as the structure, incoming arcs as one attribute, and outgoing arcs as another attribute.

GAMS does not have that and most likely you don’t “need” that anyway. A compact way to describe a network is shown here:

http://yetanothermathprogrammingconsultant.blogspot.com/2011/06/network-formulation.html

Some small improvements over this formulation for even larger problems are documented here: http://yetanothermathprogrammingconsultant.blogspot.com/2011/08/more-on-network-models.html.

Dinkelbach’s Algorithm

2012-01-16T23:46:00.001-05:00

The paper:

Dinkelbach’s Algorithm as an Efficient Method for Solving a Class of
MINLP Models for Large-Scale Cyclic Scheduling Problems [link]

by Prof. Grossmann e.a. investigates some mixed-integer linear fractional programming problems. It compares some standard MINLP codes against an implementation of Dinkelbach’s algorithm and concludes that this algorithm performs very good. The algorithm consists of a solving a series of MIP problems with each time a different weighted objective function. This method is so simple it can be coded in GAMS in just a few lines.

If the problem is something like:

then we need to solve a series of MIP problems:

for different values of λ.

I wanted to see if I could use this on a (very) large mixed-integer linear fractional programming problem. Of course first try it on some very small examples. The first test is a continuous problem:

Example 1: Linear fractional programming problem

$ontext

   References:

     Fengqi You, Pedro M. Castro, Ignacio E. Grossmann1
     Dinkelbach’s Algorithm as an Efficient Method for Solving a Class of
     MINLP Models for Large-Scale Cyclic Scheduling Problems

     Said Tantawy
     AN ITERATIVE METHOD FOR SOLVING LINEAR FRACTION
     PROGRAMMING (LFP) PROBLEM WITH SENSITIVITY ANALYSIS

$offtext

set
   i 'products' /a1,a2/
;

parameters
   uprofit(i) 'unit profit' /
      a1 4
      a2 2
   /
   ucost(i) 'unit cost' /
      a1 1
      a2 2
   /
   fixedcost /5/
   fixedprofit /10/
   matuse(i) 'raw material usage' /
      a1 1
      a2 3
   /
   rawavail 'raw material available' /30/
;

variables
   x(i) 'production'
   profit
   cost
   z    'objective variable'
;
positive variables x,profit,cost;

equations
   eprofit 'calculate profit'
   ecost   'calculate cost'
   eraw    'raw material usage'
   eprodcon 'production constraint'
   eratio   'minlp objective'
;

eprofit.. profit =e= fixedprofit + sum(i, uprofit(i)*x(i));

ecost..   cost =e= fixedcost + sum(i, ucost(i)*x(i));

eraw..    sum(i, matuse(i)*x(i)) =l= rawavail;

eprodcon.. x('a1') + 5 =g= 2*x('a2');

eratio..   z =e= profit/cost;
cost.lo = 1;

*-------------------------------------------------------------
* solve as nlp
*-------------------------------------------------------------

model m1 /all/;
m1.solprint = 2;
solve m1 maximizing z using nlp;

display "------------------------------------",
        "NLP Solver",
        "------------------------------------",
        z.l,x.l;

*-------------------------------------------------------------
* Dinkelbach's algorithm
*-------------------------------------------------------------

scalars
q 'optimal objective at end of algorithm' /0/
continue /1/
tol /0.1/
iterations
;

equation linobj;
linobj .. z =e= profit - q*cost;

model m2 /eprofit,ecost,eraw,eprodcon,linobj/;
m2.solprint = 2;

set iter /iter1*iter5/;

loop(iter$continue,
   solve m2 maximizing z using lp;
   if (z.l < tol,
      continue = 0;
      iterations = ord(iter);
      display "------------------------------------",
              "Dinkelbach algorithm",
              "------------------------------------",
              q,iterations,x.l;
   else
      q = profit.l/cost.l;
   );

);

This continuous problem is very easy:

----     73 ------------------------------------
            NLP Solver
            ------------------------------------
            VARIABLE z.L                   =        3.714 objective variable

---- 73 VARIABLE x.L production

a1 30.000

----    102 ------------------------------------
            Dinkelbach algorithm
            ------------------------------------
            PARAMETER q                    =        3.714 optimal objective at end of algorithm
            PARAMETER iterations           =            2

---- 102 VARIABLE x.L production

a1 30.000

The next small example is an integer problem:

Example 2: an integer fractional programming problem

$ontext

   References:

     Fengqi You, Pedro M. Castro, Ignacio E. Grossmann1
     Dinkelbach’s Algorithm as an Efficient Method for Solving a Class of
     MINLP Models for Large-Scale Cyclic Scheduling Problems

     Ildiko Zsigmond
     Mixed Integer linear Fractional Programming By A Branch-and-Bound
     Technique 1985

$offtext

variables
   x1,x2
   num 'numerator'
   denom 'denominator'
   z    'objective variable'
;
integer variables x1,x2;
x1.up = 5;
x2.up = 4;

equations
   enum    'numerator'
   edenom   'denominator'
   e1,e2,e3
   eratio
;

enum.. num =e= 2*x1+x2-2;

edenom.. denom =e= x1-x2+1;

e1.. -5*x1 + 4*x2 =l= 0;
e2.. -x1+x2 =l= 0.5;
e3.. 2*x1+x2 =l= 11;

eratio..   z =e= num/denom;
denom.lo = .1;

option optcr=0;

model m1 /enum,edenom,e1,e2,e3,eratio/;
m1.solprint = 2;
solve m1 maximizing z using minlp;

display "------------------------------------",
        "MINLP Solver",
        "------------------------------------",
        z.l,x1.l,x2.l;

*-------------------------------------------------------------
* Dinkelbach's algorithm
*-------------------------------------------------------------

scalars
q 'optimal objective at end of algorithm' /0/
continue /1/
tol /0.1/
iterations
;

equation linobj;
linobj .. z =e= num - q*denom;

model m2 /enum,edenom,e1,e2,e3,linobj/;
m2.solprint = 2;

set iter /iter1*iter5/;

display "------------------------------------";
loop(iter$continue,
   solve m2 maximizing z using mip;
   if (z.l < tol,
      continue = 0;
      iterations = ord(iter);
      display "------------------------------------",
              "Dinkelbach algorithm",
              "------------------------------------",
              q,iterations,x1.l,x2.l;
   else
      q = num.l/denom.l;
   );

);

The listing file shows:

----     51 ------------------------------------
            MINLP Solver
            ------------------------------------
            VARIABLE z.L                   =        7.000 objective variable
            VARIABLE x1.L                  =        3.000
            VARIABLE x2.L                  =        3.000
----     82 ------------------------------------
            Dinkelbach algorithm
            ------------------------------------
            PARAMETER q                    =        7.000 optimal objective at end of algorithm
            PARAMETER iterations           =        4.000
            VARIABLE x1.L                  =        3.000
            VARIABLE x2.L                  =        3.000

For our very large problem (approx. 1 million rows) this method was even more successful. All NLP solvers I have access to had troubles with the sheer size of the problem, even though the NLP relaxations are linearly constrained. But the MIP problems generated inside Dinkelbach’s algorithm turned out to be large but easy to solve. In addition the algorithm converged in about 5 major iterations.

Reviewing NLP model

2012-01-10T11:09:00.001-05:00

In a fairly complex model, we ended up trying a nonlinear objective:

where the rest of the model was linear (after applying some reformulations). As the model contains both integer and binary variables, this now has become an MINLP.

For several reasons this is actually not a very good expression:

this expression generates a lot of non-linear variables and non-linear non-zero elements in the matrix
as some z(i)’s can be zero it is difficult to protect this against division-by-zero (although the way we ran it had a linear model in front of this model, so our starting point was reasonable and no division-by-zero occurred)

A better formulation would be:

Here we only have two non-linear variables: the rest of the variables is now appearing linearly in the constraints. In addition, although we cannot bound each z(i) away from zero, we could introduce a nonzero lower bound on w.

Although we are adding extra constraints and variables to the MINLP model, this actually makes the model easier to solve and more robust.

PS1. The model is now becoming somewhat complex as we try to deal with different types of decisions in the same model. We call this an “integrated model” rather than just messy!

PS2. In some cases fractions can be reformulated using a “fractional programming” technique. In practice for larger models this often turns out to be a difficult reformulation to implement.

Question: it is suggested in a comment below to use (2d?) bisection. Is that really a good idea? Any references for such an approach applied in a similar situation?

Savings from applying MIP models

2011-12-25T06:23:00.001-05:00

In the first chapter of “Applied Integer Programming: Modeling and Solution”, by Der-San Chen; Robert G. Batson; Yu Dang (2010), there is a table with successful mip applications taken from articles from Interfaces. See http://media.wiley.com/product_data/excerpt/67/04703730/0470373067.pdf.

Some numbers seem low: $0.25 million savings (EK:annual?) from a least-cost crew scheduling model at American Airlines in 1981, while others seem rather high: a coal transportation study by the World Bank and Chinese State Planning Commission (1995): $6.4 billion in savings (over what period is not specified).

Bad NLP modeling

2011-12-22T16:39:00.001-05:00

In http://groups.google.com/group/knitro/browse_thread/thread/a1785c5c6f64b62e/ we see the following constructs being used in an NLP model:

c6.. y1*(y1-1) =e= 0;
c7.. y2*(y2-1) =e= 0;
c8.. y3*(y3-1) =e= 0;

This is not really a good idea. These constraints essentially say: y1,y2 and y3 can assume the values 0 or 1, i.e. they are binary variables.

If you really want to use binary variables in an NLP, form a proper MINLP model. If this idea of using non-convex constraints to model binary variables was really so hot, then we would no longer need a MIP solver and just use an NLP solver.

It looks a bit that this poster is attempting to run before (s)he can walk. It would be much better if you have some knowledge about NLP modeling before before actually working on “stochastic global optimization”. May be the advisor or teacher should play a more constructive role here.

Notes:

in general NLP solvers don’t “understand” the functions you give them. They merely sample points, function values and derivatives as they go along. One exception is the global solver Baron: this solver really tries to work with and understand the functional form of the NL functions.
although in general it is better to model binary decisions as binary variables in an MINLP (or MIP if the rest is linear), there are some attempts to use nonlinear programming techniques to solve integer models:

W. Murray and K. M. Ng, “An Algorithm for Nonlinear Optimization Problems with Binary Variables,” Computational Optimization and Applications, Vol. 47, No. 2, 2008, pp. 257-288
Roohollah Aliakbari Shandiz, Nezam Mahdavi-Amiri, “An Exact Penalty Approach for Mixed Integer Nonlinear Programming Problems”, American Journal of Operations Research, 2011, 1, 185-189

Sequence dependent setup times: difficult MIPs

2011-12-19T19:36:00.001-05:00

A simple single machine scheduling model with sequence dependent setup times is a good example where MIP solvers are not very good yet.

The simplest formulation for the 15 job problem we show below, is probably where we use the following definition:

x(i,j) = 1 if job i is scheduled before job j
0 otherwise

Obviously we don’t need the diagonal: x(i,i). Also if we know that x(i,j) = 1-x(j,i), so we don’t need both the lower- and upper-triangular part. A model now can look like:

$ontext

    Single Machine with Sequence Dependent Setup times.

    Data from: https://www.msu.edu/~rubin/files/research.html

$offtext

sets
   i0    'all jobs' /job0*job15/
   i(i0) 'jobs to schedule' /job1*job15/
;

alias(i,j), (i0,j0);

table data(*,*)
         job1 job2 job3 job4 job5 job6 job7 job8 job9 job10 job11 job12 job13 job14 job15
proctime   88   95   98 102 100   97 100   96 102    97   108   107   102   100    96
duedate 1162 1194 1195 1196 1229 1236 1240 1301 1323 1353 1368 1402 1421 1432 1433
    job0   12    7   15   19    2   14    7   19    9     3    20     0    17    13     4
    job1    0   16   19   10    2   13   17   20    3     9     6     3    19     9    13
    job2    1    0   10   19    3    7    9   18   20    14    10    11    17     5     7
    job3    8   17    0    7   19    6   10    3    3     6    14     4    12    20     8
    job4   18   13   10    0   18    1   14    4   16     3     6    13     2     2    15
    job5    2   14   14    3    0   10   12    6   20    12     1     7    16     6     3
    job6    4    2   10    1   16    0   20   13   17    10    15    19    12    18     6
    job7    5    5   19   18    5    2    0    8    5    15    17     0     0    11     2
    job8   13   18    4   15   19   13    3    0   18    14    12    16    15    16     5
    job9   19   17    9   11   17   20   20    3    0     2    10     4     0     8    13
   job10   13    6   20   19   15   13   12    7    5     0    11     5     5     8     7
   job11 11    9    3   10   11    0    9   19   18     4     0    14     3    19     3
   job12   14   17   18   10    0   13    0   20   19     8     4     0     1    18    15
   job13    7   20    1   10   12    0    7   13    5    13     5    11     0    20    4
   job14    7   14   17   13   19   10   11   10   15     9     5     2    10     0     5
   job15    0    5    6   18   12    8   20    0    1    17     7    18    12    18     0
;

parameters
proctime(i)
setuptime(i0,j)
duedate(i)
;
proctime(i) = data('proctime',i);
setuptime(i0,j) = data(i0,j);
duedate(i) = data('duedate',i);

scalar M;
~~M =~~ ~~sum((i0,j),setuptime(i0,j))+sum(i,proctime(i));~~
M = sum(j, smax(i0,setuptime(i0,j)) + proctime(j));

variables
   x(i,j)        'j later than i'
   starttime(i) 'start of job'
   endtime(i)    'end of job'
   makespan      'total makespan'
   z             'objective variable'
   late(i)       '~~lateness~~ tardiness of job'
;
binary variables x;
positive variables late,starttime,endtime;

equations
   calcstart1    'bound on starttime'
   calcstart2    'bound on starttime'
   calcstart3    'bound on starttime'
   ontime        'compare endtime with duedate'
   emakespan     'bound on makespan'
   calcendtime   'end=start+proc'
   obj           'objective'
;

set lt(i,j);
lt(i,j)$(ord(i)<ord(j))=yes;

calcstart1(lt(i,j)).. starttime(j) =g= endtime(i) + setuptime(i,j) - M*(1-x(i,j));
calcstart2(lt(j,i)).. starttime(j) =g= endtime(i) + setuptime(i,j) - M*x(j,i);
calcstart3(i).. starttime(i) =g= setuptime('job0',i);

calcendtime(i).. endtime(i) =e= starttime(i) + proctime(i);

ontime(i).. endtime(i) =l= duedate(i) + late(i);

emakespan(i).. makespan =g= endtime(i);

obj.. z =e= 0*makespan + 1*sum(i,late(i));

* initial solution: use mipstart to load it
x.l(lt(i,j)) = 1;

option optcr=0, reslim=3600;
model m1 /all/;
solve m1 minimizing z using mip;

display x.l,starttime.l,endtime.l;

We can find the optimal solution of 90 in about an hour or so (running on one thread):

However, as we can see we are far from proving this is indeed the optimal solution: the best bound has not moved at all.

Notes:

Constraint calcstart3 can be written as a bound. Of course these days, a MIP solver will do that for you in the pre-solve phase.
The initial solution was constructed by ordering the jobs by due date.
The Gantt chart for the machine can look like:

Here the late jobs are colored red.
The model minimizes sum of ~~lateness~~ tardiness. Additional objectives could be number of late jobs and total makespan.
More difficult to model is this: scheduled downtime. If jobs can not be interrupted, this causes holes (idle time) in the schedule.
A different, slightly more complicated formulation would result from:
x(i,j) = 1 if job i is scheduled immediately before job j
0 otherwise
IBM Ilog’s CP solver using OPL can model sequence dependent setup times directly (see: http://publib.boulder.ibm.com/infocenter/cosinfoc/v12r3/topic/ilog.odms.cpo.help/Content/Optimization/Documentation/Optimization_Studio/_pubskel/usrcpoptimizer2793.html).
Many papers suggest a meta-heuristic (genetic algorithm etc.) to solve this.
Of course a MIP model can still help to define the problem correctly and to solve small instances to optimality.
Paper for data file is: “SCHEDULING IN A SEQUENCE DEPENDENT SETUP ENVIRONMENT WITH GENETIC SEARCH”, Paul A. Rubin, Gary L. Ragatz, Computers & Operations Research, 22, no 1, jan 1995, pp 85-99.
The paper notes that according to a survey, 70% of scheduling problems may be subject to sequence dependent setup times (this number looks a little bit high to me).
A generalization of setup times that are sequence dependent is to consider them history dependent, i.e. dependent on a number of preceding jobs. See for instance http://onlinelibrary.wiley.com/doi/10.1002/nav.21472/abstract.

PS. After the comment from Paul Rubin, I updated the big-M calculation. This reduces the size of M from 3819 to 1775. Although this is a better formulation, an example run with Cplex, 1 thread actually shows somewhat poorer performance of the new M. In the picture below we have M1 as the old M=3819 and M2 the results for M=1775.

Sometimes with a MIP we see that a better formulation can have worse performance on certain instances.

PS2. See comment below. Officially we have tardiness = max(0,lateness). I hardly use the term tardiness: “normal people” tend to use the more intuitive earliness and lateness (both nonnegative).

Optimization for Dummies

2011-12-19T16:07:00.001-05:00

Gurobi came out with a nice little booklet explaining optimization in basic non-geek terms. It may be useful to show your boss what you want him (or her) to spend some money on. If you want a copy, contact the Gurobi people (www.gurobi.com).

Actually I have purchased a few of these “for dummies” books – there must be a ton of them by now –, and in general they are quite well written.

Logs in NLP

2011-12-14T23:55:00.001-05:00

In http://www.or-exchange.com/questions/4214/entropy-maximization-log0-ampl-tricks someone seems to claim that:

y=log(x),
x.lo = 0.001

is not good modeling practice. I don’t agree: the better NLP solvers don’t evaluate nonlinear functions outside their bounds. And from my more than 20 year experience in nonlinear modeling: I have used this construct all the time.

A large percentage of NLPs that fail can be fixed by looking into specifying:

better bounds
better starting point
and better scaling

The comment below talks about IPOPT. No, even IPOPT will not evaluate nl funcs outside their bounds (come on, it is an interior point code!!). See http://drops.dagstuhl.de/volltexte/2009/2089/pdf/09061.WaechterAndreas.Paper.2089.pdf . Slight complication is bound_relax_factor, but in practice that is not a problem. Maybe the poster is confusing infeasibility and interior points: a point can be inside the bounds but still be infeasible.

Apart from this remark, indeed IPOPT is an excellent solver. It is a fantastic complement to the well-known active-set solvers like CONOPT, MINOS and SNOPT. If you have large models with many superbasics: these solvers tend to have problems with that, while an interior point algorithm does not really care.

Student posts (2)

2011-12-14T23:47:00.001-05:00

One more for my collection (http://yetanothermathprogrammingconsultant.blogspot.com/2010/09/student-mails.html):

Bonjour,

Je viens de voir votre blog sur la programmation linéaire. En fait je débute sur GLPK et je voulais
vous poser la question: comment exprimer une condition dans une contrainte?
J'ai utilisé if mais ça marche pas.

Merci d'avance.
Cordialement,

Actually I understand enough French to understand this, but nevertheless it is somewhat indicative of how little effort is being employed before sending emails around asking for help.

Analytics in Retail

2011-12-14T23:41:00.001-05:00

The last few day I was working with some people on some reasonably complicated optimization problems in the retail industry. One of the great resources on using analytics in that environment is:

http://www.amazon.com/New-Science-Retailing-Transforming-Performance/dp/1422110575/

On my way back from Chicago there was bad weather and planes breaking down, so I had quite some time to read this. Highly recommended to get a feeling about some issues in this industry that can be tackled by some math.

During the work we found again that it really helps to try to develop some models very early in the project. Even with invented data, it allows you to really focus on the problem. As a result you will find the weak spots in your knowledge about the problem and will generate often very pointed questions. The exercise of inventing data is also a very good tool to get a better understanding.

Modeling is more difficult than many people think…

2011-12-01T07:53:00.001-05:00

From a Bachelors thesis:

In this simple transportation model I see at least two problems.

The capacity constraint should probably look like:

forall (o in DCs)
ctCapacity:
sum(p in Products, d in Stores)
Trans[p][o][d] <= Capacity[o];

The care needed to write correct models is often underestimated.

Matrix multiplication

2011-12-01T00:23:00.001-05:00

Theory: http://www.scottaaronson.com/blog/?p=839. This really seems to excite the theorists.
In practice: http://yetanothermathprogrammingconsultant.blogspot.com/2010/06/sparse-matmul-revisited.html

MPS files: this does not look right

2011-11-30T16:46:00.001-05:00

Writing an MPS for a QP model from GAMS using the Cplex option writemps I saw this:

* ENCODING=ISO-8859-1
NAME          gamsmodel
ROWS
N obj
E obj
E budget
G retcon
COLUMNS
    z         obj                             1
    z         obj                             1
    x(GAB)    budget                          1
    x(GAB)    retcon                     0.0025
    x(GAP)    budget                          1
    x(GAP)    retcon                     0.1375

…

This does not look right. Probably a problem in the GAMS link where a new free row is added to the problem with the same name as an existing row.

This does not compute

2011-11-23T11:05:00.001-05:00

GLPK doesn't accept an constraint in which a variable is devided by an variable, i would like to know if there is some solution to linearize this constraint. knowning that var2 is binary and var 1 is reel in case var1 / var 2.

This means var2=1 and we can replace the expression var1/var2 by var1.

Recently read

2011-11-23T05:57:00.001-05:00

During a transatlantic flight I had the opportunity to read this book on parallel programming for .Net. It describes some of the higher level constructs available in .Net that make it easier to implement parallel algorithms. Of course the big issue is that you need to (re)design your algorithms to make parallel operations effective or even possible. Some algorithms are just not very amenable to be parallelized (the Simplex method for LP being a good example). The book contains enough warnings to make you think twice before starting to work on parallel implementations.

MPS Files (2)

2011-11-21T16:41:00.001-05:00

As follow up on http://yetanothermathprogrammingconsultant.blogspot.com/2011/10/mps-files.html it is noted that Microsoft Solver Foundation also reads some MPS files incorrectly. E.g. consider the files:

C:\Users\erwin\Documents>type mip4mps.oml
Model[
Decisions[Integers[0, Infinity], x],
Constraints[x <= 5],
Goals[Maximize[x]]
]

C:\Users\erwin\Documents>type output.mps
NAME          MODEL
ROWS
L constrai
N goal_8ef
COLUMNS
    INTMARK   'MARKER'                 'INTORG'
    x         constrai   1             goal_8ef -1
    INTMARK   'MARKER'                 'INTEND'
RHS
    RHS       constrai   5
BOUNDS
ENDATA

C:\Users\erwin\Documents>

The MPS file was generated from the OML file. If we inspect the MPS file we see no bounds on x. This means the variable x should be binary.

If we solve the MPS file with Solver Foundation we see:

C:\Users\erwin\Documents>"\Program Files (x86)\Microsoft Solver Foundation\3.0.2.10889\Bin\MSFCli.exe" +verbose 3 output.mps
===== Processing .\output.mps =====
Solution Quality: Optimal
===Solver Foundation Service Report===
Date: 11/21/2011 10:27:00 PM
Version: Microsoft Solver Foundation 3.0.2.10889 Standard Edition
Model Name: DefaultModel
Capabilities Applied: MILP
Solve Time (ms): 367
Total Time (ms): 489
Solve Completion Status: Optimal
Solver Selected: Microsoft.SolverFoundation.Solvers.SimplexSolver
Algorithm: Primal
Arithmetic: Hybrid
Variables: 1 -> 1 + 2
Rows: 2 -> 2
Nonzeros: 2
Eliminated Slack Variables: 0
Pricing (exact): SteepestEdge
Pricing (double): SteepestEdge
Basis: Slack
Pivot Count: 1
Phase 1 Pivots: 0 + 0
Phase 2 Pivots: 1 + 0
Factorings: 4 + 1
Degenerate Pivots: 0 (0.00 %)
Branches: 0
===Solution Details===
Goals:
goal_8ef: -5

Decisions:
x: 5
===== Finished .\output.mps: 00:00:00.5501834 =====

This is incorrect. The reported value of x=5 indicates x is not read as a binary variable. If we solve the MPS file with CBC we see:

C:\Users\erwin\Documents>..\Downloads\cbc output.mps -solve -solution stdout
Welcome to the CBC MILP Solver
Version: 2.7.5
Build Date: Nov 10 2011
Revision Number: 1759

command line - ..\Downloads\cbc output.mps -solve -solution stdout (default strategy 1)
At line 1 NAME MODEL
At line 2 ROWS
At line 5 COLUMNS
At line 9 RHS
At line 11 BOUNDS
At line 12 ENDATA
Problem MODEL has 1 rows, 1 columns and 1 elements
Coin0008I MODEL read with 0 errors
Continuous objective value is -1 - 0.00 seconds
Cgl0004I processed model has 0 rows, 0 columns (0 integer) and 0 elements
Cbc3007W No integer variables - nothing to do
Cuts at root node changed objective from -1 to -1.79769e+308
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                -1.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Optimal - objective value -1.00000000
0 x 1 -1
Total time (CPU seconds): 0.02 (Wallclock seconds): 0.02

C:\Users\erwin\Documents>

The output is a little bit confusing, but indeed x assumes the value 1. This confirms x is binary.

Rick Rosenthal: Ten Keys to Success in Optimization Modeling

2011-11-16T21:10:00.001-05:00

I recently talked to a an ex-student of Rick Rosenthal. The following presentation contains some concise descriptions of his ideas and beliefs on modeling:

http://faculty.nps.edu/rrosent/docs/rosenthalAtlanta.ppt

HPC applications in the Amazon cloud

2011-11-16T21:03:00.001-05:00

From http://aws.amazon.com/hpc-applications/:

“For example, a 1064 instance (17024 cores) cluster of cc2.8xlarge instances was able to achieve 240.09 TeraFLOPS for the High Performance Linpack benchmark, placing the cluster at #42 in the November 2011 Top500 list.”

Sulum Optimization

2011-11-14T17:56:00.001-05:00

Bo Jensen at http://www.sulumoptimization.com announced the availability of a new commercial LP solver. Some first impressions:

LP (+network) only for now, MIP will come later.
Aggressively priced compared to the big ones (Cplex, Gurobi). I can see this can be a very attractive offering in certain cases.
But it has a nasty license manager (somehow I always have troubles with these). The other commercial guys have similar drawbacks.
I think if you have a single license you can solve one LP at the time (a result of the previous point).
Low level API. In many cases, to be efficient you will need to use some more advanced modeling tools.

The product seems to be well-positioned between the Open Source offerings such as Lp-Solve, GLPK and CBC and the heavy guns (Cplex, Gurobi, Xpress).

Automatic Reformulations

2011-11-09T07:04:00.001-05:00

More and more we see that modeling systems apply automatic reformulations.

For instance I see that Cplex’s Concert accepts

model.add (IloMin (x, y) == 0)

for x≥0,y≥0. Essentially this is “x=0 or y=0”. This construct seems to get translated into

a+b≥1
a=1 ↔ x = 0
b=1 ↔ y = 0
x,y≥0, a,b∈{0,1}

Here the indicator constraints of Cplex are used, so no upper-bounds are needed. For other MIP solvers one would likely use something like:

x ≤ b * x.up
y ≤ (1-b) * y.up
x,y≥0, b∈{0,1}

Similarly in the language CMPL we see that x*y is accepted and translated if x or y is a binary or integer variable (see http://www.coliop.org/download/cmpl/CMPL-1-5-2.pdf section 7.3).

Solver Options

2011-11-08T23:26:00.001-05:00

In general I don’t like to use solver options to improve performance. I would hope that a solver would do a good job with default settings. Hopefully it will figure out from the model what settings it should use, an if not, will detect slow performance and does something about it. In addition I often feel that my own solver settings are not that stable. A next version of the solver may well show slower performance with my own settings, and a next version of the model may have similar issues. Also there is this statement I have heard someone saying: “each solver option is a situation where the developer did not know what to do”. Sometimes I notice users are very impressed by the many options a solver offers (may be also intimidated by the sometimes obscure but very technically advanced sounding descriptions). Here I almost argue the opposite. As a result, I often spent more time on model reformulations than on solver options.

But sometimes the differences are so dramatic, I have to break my own rules:

C:\projects\lsi\src\Development\TMS\Source\TestResults\erwin_ERWIN-THINK 2011-11-07 14_32_27\Out>type model.oml
Model[
Decisions[
Integers[0, 58], This model is a little bit ugly as it was generated automatically by SaveModel.
Foreach[ However it shows it is pretty small.
{iter1, 59},
x[iter1]
]
],
Constraints[
prereq0 -> x[13] < x[0],
prereq1 -> x[20] < x[1],
prereq2 -> x[17] < x[2],
prereq3 -> x[7] < x[3],
prereq4 -> x[19] < x[4],
prereq5 -> x[26] < x[5],
prereq6 -> x[8] < x[6],
prereq7 -> x[21] < x[7],
prereq8 -> x[4] < x[8],
prereq9 -> x[12] < x[10],
prereq10 -> x[0] < x[11],
prereq11 -> x[18] < x[12],
prereq12 -> x[3] < x[13],
prereq13 -> x[9] < x[14],
prereq14 -> x[27] < x[15],
prereq15 -> x[5] < x[16],
prereq16 -> x[15] < x[17],
prereq17 -> x[6] < x[18],
prereq18 -> x[16] < x[19],
prereq19 -> x[22] < x[20],
prereq20 -> x[23] < x[21],
prereq21 -> x[14] < x[22],
prereq22 -> x[25] < x[23],
prereq23 -> x[10] < x[24],
prereq24 -> x[24] < x[25],
prereq25 -> x[2] < x[26],
prereq26 -> x[1] < x[27],
prereq27 -> x[29] < x[28],
prereq28 -> x[31] < x[29],
prereq29 -> x[32] < x[30],
prereq30 -> x[33] < x[31],
prereq31 -> x[28] < x[32],
prereq32 -> x[36] < x[34],
prereq33 -> x[34] < x[35],
prereq34 -> x[38] < x[36],
prereq35 -> x[35] < x[37],
prereq36 -> x[42] < x[39],
prereq37 -> x[39] < x[40],
prereq38 -> x[44] < x[40],
prereq39 -> x[45] < x[40],
prereq40 -> x[49] < x[40],
prereq41 -> x[53] < x[40],
prereq42 -> x[54] < x[40],
prereq43 -> x[56] < x[40],
prereq44 -> x[56] < x[41],
prereq45 -> x[56] < x[42],
prereq46 -> x[56] < x[43],
prereq47 -> x[50] < x[44],
prereq48 -> x[52] < x[45],
prereq49 -> x[57] < x[45],
prereq50 -> x[56] < x[46],
prereq51 -> x[56] < x[47],
prereq52 -> x[39] < x[48],
prereq53 -> x[44] < x[48],
prereq54 -> x[45] < x[48],
prereq55 -> x[49] < x[48],
prereq56 -> x[53] < x[48],
prereq57 -> x[54] < x[48],
prereq58 -> x[56] < x[48],
prereq59 -> x[43] < x[49],
prereq60 -> x[46] < x[49],
prereq61 -> x[56] < x[50],
prereq62 -> x[56] < x[52],
prereq63 -> x[41] < x[53],
prereq64 -> x[47] < x[54],
prereq65 -> x[40] < x[55],
prereq66 -> x[48] < x[55],
prereq67 -> x[51] < x[56],
prereq68 -> x[56] < x[57],
alldiff -> Unequal[Foreach[
{iter2, 59},
x[iter2]
]]
]
]

C:\projects\lsi\src\Development\TMS\Source\TestResults\erwin_ERWIN-THINK 2011-11-07 14_32_27\Out>"\Program Files (x86)\Microsoft Solver Foundation\3.0.2.10889\Bin\MSFCli.exe" +verbose 0 model.oml
===== Processing .\model.oml =====
Solution Quality: Feasible
===Solver Foundation Service Report===
Date: 11/8/2011 11:09:55 PM
Version: Microsoft Solver Foundation 3.0.2.10889 Express Edition
Model Name: DefaultModel
Capabilities Applied: CP
Solve Time (ms): 184798 Timings with default settings
Total Time (ms): 185031
Solve Completion Status: Feasible
Solver Selected: Microsoft.SolverFoundation.Solvers.ConstraintSystem
===Solution Details===
Goals:
===== Finished .\model.oml: 00:03:05.0714386 =====

C:\projects\lsi\src\Development\TMS\Source\TestResults\erwin_ERWIN-THINK 2011-11-07 14_32_27\Out>"\Program Files (x86)\Microsoft Solver Foundation\3.0.2.10889\Bin\MSFCli.exe" +variable conflictdriven +verbose 0 model.oml
===== Processing .\model.oml =====
Solution Quality: Feasible
===Solver Foundation Service Report===
Date: 11/8/2011 11:06:38 PM
Version: Microsoft Solver Foundation 3.0.2.10889 Express Edition
Model Name: DefaultModel
Capabilities Applied: CP
Solve Time (ms): 203 Timings with option
Total Time (ms): 433
Solve Completion Status: Feasible
Solver Selected: Microsoft.SolverFoundation.Solvers.ConstraintSystem
===Solution Details===
Goals:
===== Finished .\model.oml: 00:00:00.4739553 =====

Note: it is not surprising that Cplex added a tuning option to suggest some good settings. A related effort is: http://www.cs.ubc.ca/labs/beta/Projects/MIP-Config/.

Finding feasible schedules for a column generation algorithm

2011-11-07T18:57:00.001-05:00

For a complex scheduling application, I need to generate a large number of feasible schedules for a class. Basically:

find x[1]…x[n] in {1..n}
under some precedence constraints x[i] < x[j]
all-different(x[1]…x[n])

MIP models are not very good for this (see http://yetanothermathprogrammingconsultant.blogspot.com/2011/10/sorting-by-mip.html). But a constraint programming model should be very easy to formulate and may perform quite good.

With Microsoft Solver Foundation I had to use some nondefault settings (VariableSelection=ConfictDriven) to get good performance. Genererating multiple solutions is easily performed with Solution.GetNext(). I don’t think it is easy to generate a subset of solutions that has a certain spread. The only thing that comes to mind is doing many calls to Solution.GetNext() and just using each K-th solution. As I don’t know in advance how many solutions there are, I am thinking about increasing K as we go.

For some classes I have simply a very large number of possible schedules (or differently put: not enough precedence constraints to reduce this number). To find a reasonable subset we could augment the above system of equations with a random objective. Now we need to restate the objective each cycle, so we no longer can use Solution.GetNext().

From what I see most modeling systems have available constraint programming support (e.g. Ilog OPL, MS Solver Foundation) or have something announced (AMPL, AIMMS).

Update: link from comment below: www.nicta.com.au/pub?doc=4759.

Slacks (with penalties) should not be integer variables

2011-10-31T03:17:00.001-04:00

In the paper http://www.cs.nott.ac.uk/~jxm/timetabling/gor2007-paper.pdf a number of integer models for timetabling are discussed. The author also published the models used in this paper http://www.cs.nott.ac.uk/~jxm/timetabling/gor2007-empirical.zip (this should be very much applauded). Of course that means that I can provide some criticism.

In the scheduling models some soft constraints are modeled by introducing slacks, in this case CourseMinDayViolations:

subto ctCountWorkingDays:
forall <c> in Courses:
sum <d> in Days:
CourseSchedule[c,d] >= CourseHasMindays[c] - CourseMinDayViolations[c];

This variable can then be minimized in the objective:

minimize value:
    (sum <p> in Periods: sum <r> in Rooms: sum <c> in Courses with CourseHasStudents[c] > RoomHasCapacity[r]:
      (Taught[p,r,c] * (CourseHasStudents[c] - RoomHasCapacity[r]))
    )
+ 2 * (sum <cu> in Curricula: sum <d> in Days: sum <sc> in DaySingletonChecks: SingletonChecks[cu,d,sc])
+ 5 * (sum <c> in Courses: CourseMinDayViolations[c])
;

That looks all good. But I don’t agree with the declaration:

var CourseMinDayViolations[Courses] integer >= 0 <= nbDays;

A slack variable like this should be declared as a positive continuous variable. It will be integer automatically. With this declaration there is the possibility the solver will branch on this variable.

My simplest argument would be: a model with fewer integer variables is likely to solve faster.

But sometimes solvers actually add integer variables! See:

Tried aggregator 1 time.
MIP Presolve eliminated 1634 rows and 5 columns.
MIP Presolve modified 10 coefficients.
Reduced MIP has 2401 rows, 11512 columns, and 92222 nonzeros.
Reduced MIP has 11028 binaries, 0 generals, 0 SOSs, and 0 indicators.
Probing time = 0.02 sec.
Tried aggregator 1 time.
MIP Presolve eliminated 10 rows and 0 columns.
Reduced MIP has 2391 rows, 11512 columns, and 91732 nonzeros.
Reduced MIP has 11032 binaries, 480 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.20 sec.

(This model had originally only binary and continuous variables).

Still not proven optimal…

2011-10-26T23:14:00.001-04:00

This is a very small MIP:

Starting Gurobi...
Optimize a model with 24 rows, 54 columns and 132 nonzeros
Presolve removed 6 rows and 6 columns
Presolve time: 0.11s
Presolved: 18 rows, 48 columns, 114 nonzeros
Variable types: 12 continuous, 36 integer (0 binary)

so I expected to run this very quickly on a fast 8 core Xeon machine. The time limit of 9999 seconds should be no problem even with the OPTCR (rel. gap tolerance) of 0.0. But I see:

835689342 7572026 infeasible   63         0.00062    0.00000   100%   1.0 9980s
836077602 7576642    0.00000   65   14    0.00062    0.00000   100%   1.0 9985s
836454289 7582365    0.00000   56   15    0.00062    0.00000   100%   1.0 9990s
836850333 7587882    0.00000   45   15    0.00062    0.00000   100%   1.0 9995s

Cutting planes:
Gomory: 1
MIR: 2

Explored 837155963 nodes (841647794 simplex iterations) in 9999.00 seconds
Thread count was 8 (of 8 available processors)

Time limit reached
Best objective 6.179540598832e-04, best bound 0.000000000000e+00, gap 100.0000%
MIP status(9): Optimization terminated due to time limit.

This is not a Gurobi problem: Cplex is just as bad. Apparently these branch and bound methods are just not suited very well for this problem. I have solved some huge problems with hundreds of thousands of binary variables, but here we see some evidence, not all problems are that easy.

Recently I saw someone writing:

gurobi> m.setParam('Threads',80)
You may all drool now.

Guess we need something like that to solve my tiny 36 integer variable problem….

Notes:

The OPTCR (gap tolerance) setting is not an issue as it turns out the best bound remains at 0.0, so a gap of 100% is reported.
Gurobi is solving the nodes pretty efficiently. Approx. 1 simplex iteration per node, and > 80k nodes per second.

Update

Some Cplex runs solve this fast (see comments below), but on some machines I see the same problem:

Cplex 12.3.0.0

Reading data...
Starting Cplex...
Tried aggregator 2 times.
MIP Presolve eliminated 1 rows and 1 columns.
MIP Presolve modified 36 coefficients.
Aggregator did 6 substitutions.
Reduced MIP has 18 rows, 48 columns, and 114 nonzeros.
Reduced MIP has 0 binaries, 36 generals, 0 SOSs, and 0 indicators.
…
Elapsed real time = 9880.28 sec. (tree size = 1251.36 MB, solutions = 70)
Nodefile size = 1123.85 MB (527.18 MB after compression)
478309344 6481362        0.0000    15        0.0012        0.0000 5.26e+008 100.00%
479049990 6483695        0.0000    15        0.0012        0.0000 5.27e+008 100.00%
479749224 6499596        0.0000    15        0.0012        0.0000 5.27e+008 100.00%
480524806 6502892    infeasible              0.0012        0.0000 5.28e+008 100.00%
481284148 6500253        0.0000    16        0.0012        0.0000 5.29e+008 100.00%
482022525 6501343        0.0000    16        0.0012        0.0000 5.30e+008 100.00%
482764367 6511470        0.0000    16        0.0012        0.0000 5.31e+008 100.00%

Mixed integer rounding cuts applied: 2
Zero-half cuts applied: 2
Gomory fractional cuts applied: 2

Root node processing (before b&c):
Real time             =    0.03
Parallel b&c, 4 threads:
Real time             = 9998.98
Sync time (average)   = 420.45
Wait time (average)   =   45.69
                          -------
Total (root+branch&cut) = 9999.01 sec.
MIP status(107): time limit exceeded
Fixing integer variables, and solving final LP...
Parallel mode: deterministic, using up to 2 threads for concurrent optimization.
Tried aggregator 1 time.
LP Presolve eliminated 25 rows and 55 columns.
All rows and columns eliminated.
Fixed MIP status(1): optimal

Resource limit exceeded.

MIP Solution: 0.001230 (530937464 iterations, 482992718 nodes)
Final Solve: 0.001230 (0 iterations)

Best possible:           0.000000
Absolute gap:            0.001230
Relative gap:            1.000000

This is a worse solution than Gurobi found. But on a different machine using 8 threads:

Cplex 12.3.0.0

Mixed integer rounding cuts applied: 3
Gomory fractional cuts applied: 1

Root node processing (before b&c):
Real time             =    0.12
Parallel b&c, 8 threads:
Real time             =    4.10
Sync time (average)   =    0.39
Wait time (average)   =    0.02
                          -------
Total (root+branch&cut) =    4.23 sec.
MIP status(101): integer optimal solution
Fixing integer variables, and solving final LP...
Parallel mode: deterministic, using up to 2 threads for concurrent optimization.
Tried aggregator 1 time.
LP Presolve eliminated 25 rows and 55 columns.
All rows and columns eliminated.
Fixed MIP status(1): optimal

Proven optimal solution.

MIP Solution: 0.000616 (175023 iterations, 203294 nodes)
Final Solve: 0.000616 (0 iterations)

Best possible:           0.000616
Absolute gap:            0.000000
Relative gap:            0.000000

So we have some very different behavior here: optimal in just 4 seconds.

Special case of integer cuts

2011-10-20T21:12:00.001-04:00

In http://yetanothermathprogrammingconsultant.blogspot.com/2011/10/integer-cuts.html we derive a constraint that forbids an integer configuration y^*.

If the model contains a constraint like:

then we can simplify the integer cut to:

as we only need to keep track of the variables flipping from 1 → 0. This is because any variable moving from 0 → 1 will cause another variable to move from 1 → 0.

This simpler constraint has fewer nonzeroes which is important if we add a lot of such cuts, as is the case in one of my models.

MPS files

2011-10-18T05:35:00.001-04:00

Yesterday we found out that Microsoft Solver Foundation is writing incorrect MPS files when the model has general integer variables without an upper bound. See: http://social.msdn.microsoft.com/Forums/en-US/solverfoundation/thread/7669081f-18ac-445b-850e-434bff80da38.

Right after I submitted an error report, the following e-mail arrived:

I trying to convert one problem on gmpl format to free mps format, both in attached, and the solve it.
Unfortunately, when solve the mps file I get an answer different from the answer that I get when solve using the gmpl file.

I was using the following commands:
$ glpsol --math exemp01.mod
$ glpsol --math exemp01.mod --check --wfreemps exemp01.mps
$ glpsol --freemps exemp01.mps

See: http://lists.gnu.org/archive/html/help-glpk/2011-10/msg00049.html

Lo and behold, this is exactly the same problem!

The problem is shown in the following MPS file (generated by Cplex, comments added by me):

NAME          gamsmodel
ROWS
N obj
E e1
L e2
COLUMNS
    MARK0000 'MARKER'                 'INTORG'
    x         obj                            -1
    x         e2                              1
    MARK0001 'MARKER'                 'INTEND'
    z         obj                            -1
    z         e1                              1
RHS
    rhs       e2                            200
BOUNDS
* column x is a general integer variable
* it needs an LI card to make sure it is not
* interpreted as a binary variable
LI bnd       x         0
FR bnd       z
ENDATA

This LI card is missing from the GLPSOL and Solver Foundation generated files. The rule is: if no bounds are specified, an integer variable has default bounds 0 and 1. Obviously this will change the meaning of the model.

Note: In GAMS by default all integer variables have a finite upper bound, so no problem there. It is possible (with some effort) to create a model with +INF bounds on integer variables (using PF4=0 command line option). It turns out that with such a model, the tool CONVERT (this tool can write MPS files), produces an MPS file with the same error. Admittedly somewhat contrived, but not correct anyway.

Sorting by MIP

2011-10-16T20:31:00.001-04:00

Here is a toy model that illustrates that certain models are poorly suited for solving by MIP solvers.

Sorting data is a rather elementary operation. Surprisingly, GAMS has nothing built in for this. There is an external tool for sorting (gdxrank), but in some environments (e.g. NEOS) this is not allowed. There are some slower GAMS coding tricks that can help. Of course we can also write a simple MIP model that performs sorting:

$ontext

   mip model for sorting

$offtext

set i /i1*i100/;

parameter p(i) "original data";
p(i) = uniform(0,1);

alias (i,j);
binary variable x(i,j) "permutation matrix";
variable y(i) "sorted version of p";
variable z;

equations
dummy    "dummy objective"
assign1 "assignment constraint"
assign2 "assignment constraint"
ydef     "y is permuted version of x"
sort     "require y to be sorted"
;

dummy.. z =e= 0;

assign1(i).. sum(j, x(i,j)) =e= 1;
assign2(j).. sum(i, x(i,j)) =e= 1;
ydef(i)..    y(i) =e= sum(j, x(i,j)*p(j));
sort(i+1).. y(i+1) =g= y(i);

model m /all/;
solve m minimizing z using mip;

display
"-----------------------------------------------------------------------",
"Results",
"-----------------------------------------------------------------------",
   p,y.l;

This is quite slow however. n=100 is no problem, but n>200 will require significant time (>1000 seconds). Constraints assign1 and assign2 form an assignment problem, which is easy to solve. The simple extension will cause relaxation based MIP solvers to slow down enormously. Obvious question: are there better formulations? I would guess constraint programming solvers would do better on this.

For some instances Gurobi is doing very good, probably because of its excellent heuristics. But this is not the case for all instances: in some cases it is just as slow as other solvers.

The above model contains two noteworthy modeling tricks:

a dummy objective so that we stop after the first integer solution is found
the equation sort is indexed over i+1 so that the last one is not generated: this equation has card(i)−1 members.

MIQP: optimality tolerances vs objective scaling

2011-10-16T01:22:00.001-04:00

I was working on a MIQP model in finance: the optimal allocation of contracts to client accounts so that each client observes an average price that is as close as possible to the average contract price. The model was implemented in Solver Foundation (for the client) and GAMS (for prototyping).

With one data set we observed Gurobi was giving slightly sub-optimal solutions. I.e. the obj was 1.396560E-5 instead of the number found by other solvers: 0. Most likely this is a tolerance problem. I tried the optimality tolerance (from 1e-6 to 1e-8 and 1e-9) but that did not help. That option is just for the LP/QP subproblems. What actually did the trick was multiplying the objective by a 1000:

Goals[Minimize["Objective"->
1000*Sum[{a,accounts}, (PriceAlloc[a]-AvgPrice)^2]
]],

Apparently this is not scaled away (somewhat to my surprise), allowing us to convey the solver not to stop too early.

I also checked: reset the relative gap tolerance to 0 from its default of 1.0e-4. In the GAMS and Solver Foundation version of the model this did not help.

Update:

On some new data sets we still see somewhat unstable behavior: different solutions with different objs by just reordering data.
Gurobi solves KKT conditions, so may need to tighten feasibility tolerance. This helps somewhat, but not completely. Question: I suppose this introduces some danger of declaring models infeasible.
I get better stability when formulated as a MIP (using an absolute value formulation).

Solving a large sparse system of equations with an LP solver

2011-10-04T12:40:00.001-04:00

When we want to solve a large, sparse system of linear equations

using an LP solver we expect this to be reasonably fast. After all, solving systems of equations is what an LP solver does all the time. So we introduce a zero objective and solve. With Gurobi a sample problem gives:

Starting Gurobi...
Optimize a model with 9801 rows, 9801 columns and 48609 nonzeros
Presolve removed 4599 rows and 4599 columns
Presolve time: 0.06s
Presolved: 5202 rows, 5202 columns, 47242 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s
    3958      handle free variables                          5s
    4975      handle free variables                         11s
    5199    0.0000000e+00   0.000000e+00   0.000000e+00     25s
    5199    0.0000000e+00   0.000000e+00   0.000000e+00     25s

Solved in 5199 iterations and 25.19 seconds

Actually we can much improve on this. By specifying an advanced basis. We know the equations (slacks) will be non-basic in the solution and all the variables will be basic. If we set up an advanced basis for this we see for the same problem:

Starting Gurobi...
Optimize a model with 9801 rows, 9801 columns and 48609 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 0.0000000e+00 0.000000e+00 0.000000e+00 0s

Solved in 0 iterations and 0.45 seconds

Indeed no simplex iterations were needed: just an invert.

Update

In response to a comment about Octave: that is still faster:

octave-3.2.4.exe:3> tmp=load("-ascii","mat.txt");
octave-3.2.4.exe:4> b=load("-ascii","rhs.txt");
octave-3.2.4.exe:5> a=sparse(tmp(:,1),tmp(:,2),tmp(:,3));
octave-3.2.4.exe:6> nnz(a)
ans = 48609
octave-3.2.4.exe:7> tic();x=a\b;toc()
Elapsed time is 0.057 seconds.
octave-3.2.4.exe:8>

Notes:

Not completely fair of course: the LP solver does still a lot of extra things (scaling, pricing etc.).
Note that I used a sparse matrix in Octave. Not sure what method they use to solve this system.
This was a structured test matrix: see: http://www.public.iastate.edu/~akmitra/aero361/design_web/Laplace.pdf. For more unstructured matrices Gurobi may do much better relative to Octave.

MIP caveats

2011-10-03T13:41:00.001-04:00

In http://www.or-exchange.com/questions/3316/metrics-rules-of-thumb-tool-support-to-assess-mip-models a number of caveats when doing MIP modeling are listed:

All 0 objective: although it should be easier to find any solution than one specific, in certain situations it might pay off to add a fake objective function to guide the search. Anything reasonably stable should work and preferable all variables would have different objective coefficients. This sometimes works, sometimes not. I guess it also has to do with the fact that solvers usually detect this and do something about it.
Unstable objective: Very big and very small coefficients in the objective function, especially values below the dual tolerance, can cause issues with dual feasibility and very unexpected results.
One or few objective coeffcients: One can observe that MIP solvers have trouble solving instances with only one (or a few) non-zeros in the objective function. At least part of the problems with these models is dual degeneracy.

Point 1: I am not sure about this one. With Obj=0 we really want to stop as soon as a feasible integer solution is found. With an obj<>0 this is no longer the case (or am I misreading this; may be additionally the stopping criterion “stop after first feasible solution” is implied). I don’t have many MIPs without an objective (and first integer solution found is not always very good), so my experience is limited. Somehow I feel this would need more empirical evidence.

Point 2: agreed. E.g. with long planning horizons (like in forestry models) some of the reduced cost can become very small. Playing with LP optimality tolerances can help in such cases (I usually don’t want to touch those, but this is an exception).

Point 3: I don’t understand that completely. I have many MIPs with very few cost coefficients (often just one: e.g. max profit). I have seen cases where small “random” obj coefficients can help to break symmetry and force solutions to be different. But that is different than is described here.

I.e. in general these are not the problem areas I would quickly consider if I have a difficult MIP. Either this user has access to different solvers than I do or he is solving different types of models. Or both.

Anyway, it is interesting to see I have such a different experience about problem areas in MIP modeling. I also think that many “older” tricks (such as I mentioned under point 3, another one is using priorities) are more and more becoming obsolete due to smarter and more sophisticated MIP solvers. I have a few “old” models with some tuning tricks, that really nowadays solve faster without all the tricks.

Ateji defunct: OptimJ now free

2011-10-02T18:25:00.001-04:00

Unfortunately, the company Ateji has ceased to operate (http://www.ateji.com/blog/ateji-is-closed/, http://www.ateji.com/). They had an interesting modeling framework based on Java.

Sometimes I divide modeling tools up as follows:

Matrix generators: write LP or MPS file from application. This is done since the very early times. The main difference is that Fortran is no longer the obvious choice for programming language to use.
Specialized modeling languages: AMPL, GAMS, …
Low level solver APIs. This is essentially like writing a matrix generator, minus the part where you write a file. In many cases the solver API can help a little bit by providing routines like addrow() and addcol().
Hybrids, somewhere between 2 and 3, where modeling constructs are added to traditional programming languages. E.g. Cplex Concert, and OptimJ.

I suspect that OptimJ is mainly attractive for the intersection of the set of Java programmers and the set of modelers of optimization models. Besides that this subset may be somewhat small (many modelers are not proficient Java programmers), there is also a lot of competition targeting this particular group of users.

I believe that some of the more complex data issues in practical modeling are often better (i.e. more efficiently) dealt with in a specialized language than in a traditional programming language. In statistics this has largely been accepted. Much statistical work is done inside specialized “language” environments such as SAS, SPSS and Stata. In optimization this is not (yet) the case: the lower level APIs still have a large following. This while data is often more irregular in practical optimization than in statistics. Of course much statistical computing is done “off-line” so the integration issues that plague modeling languages are just not there.

Integer cuts

2011-10-02T17:57:00.001-04:00

I received a few questions about how to form constraints that exclude a known 0-1 integer solution. Assume our solution we want to exclude is called y^*then we can write:

This can be linearized as follows:

The final result is also sometimes written as the following constraint:

Note that y is a binary decision variable.

Time schedules

2011-09-20T04:30:00.001-04:00

Here is a good example how a web based rostering system can look like: http://www.timeforge.com. The videos at http://www.timeforge.com/manual/video-manual/ are an easy way to see how this all works.

They are in fact using MIPs in some of their algorithms!

Cplex MIP and time limits

2011-09-19T19:37:00.001-04:00

Looks like Cplex tries to do something smart just before it hits a time limit: see the logs in http://www.ibm.com/developerworks/forums/thread.jspa?messageID=14684289&tstart=0#14684289. In there you see that the gap is reduced just before the time limit kicks in, almost independent of the value of the time limit itself.

Above: run with 100 seconds time limit and below: 1000 seconds time limit.

You would almost think it would be beneficial for Cplex to have this end-of-time behavior every now and then.

Note: these runs are on one thread (thus hopefully deterministic – may be apart from the exact iteration number where the time limit kicks in). The statement “The answer to this is an implementational detail: when you use time limits, the solving process becomes non-deterministic close to the end of the time limit.” does not completely satisfy me. I would like to have a better understanding about what is happening here.

GPU computing on Amazon cloud

2011-09-19T05:49:00.001-04:00

May be this can help bringing parallel computation on GPUs to business applications:

Dear Amazon Web Services Customer,

We are excited to announce that starting today you can run Microsoft Windows Server 2008 R2 and Microsoft Windows SQL Server Standard 2008 R2 on Cluster Compute and Cluster GPU instances for Amazon EC2. Cluster Compute Instances provide customers with high CPU capability within a high bandwidth, low latency network for IO intensive computing. Cluster GPU instances allow customers to additionally take advantage of NVidia Tesla GPUs for general purpose GPU (GPGPU) computing using the CUDA or OpenCL programming models. A number of customers with database workloads and applications in areas such as media processing, rendering, and computational finance, have expressed interest in running Microsoft Windows Server on the most powerful Amazon EC2 instances. Starting today, customers can run applications on Microsoft Windows Server across all Amazon EC2 instance types.
Cluster Compute and Cluster GPU instances are currently available in the US-East (N. Virginia) region. For more information, please visit http://aws.amazon.com/ec2 andhttp://aws.amazon.com/windows.

Sincerely,
The Amazon EC2 Team

Big-M and Scaling

2011-09-13T09:29:00.001-04:00

I was working on a client MIP model that showed strange behavior: with Cplex and Mosek optimal values were found while Gurobi and Xpress reported the model to be infeasible. There were a few Big-M values in the model, in equations like:

Some of the big M values were very large e.g. 99999999. The actual value of this M should be equal to the upper bound of x. Unfortunately x can assume very large values, something like $100,000,000. This meant we could only reduce the big M values by a small amount. With some trial and error we found some big M values that worked. Of course this was just a band-aid solution.

In the end we chose to scale the model by changing units for x. Instead of $ we now use millions of $. As a result the new big M values can just be M=100. As the model was rather complex we needed to spend half a day on this exercise. After this work, the model solved with all solvers without a problem.

Open source MIP solvers

2011-09-06T17:29:00.001-04:00

Hi Erwin,
I have a 2 questions, could you please point me into a good direction?

1. Which open source optimization solver do you recommend? (One that can withstand a large and complex problem)
I've been researching online and I've come across with lpsolve, coin-or lp (CLP), and GLPK.

2. I've been having a really tough time to get GLPK to work. I'm new to this field of open source software but I've found out that I need an interface from where I should call the GLPK library. My question is, which application (or interface) should I use?

Thanks in advance!

CBC (https://projects.coin-or.org/Cbc) seems to be the fastest open source MIP solver (see http://plato.asu.edu/ftp/milpc.html). The academic code SCIP is also quite good but is not open source (it requires a license for any non-academic application). Of course the commercial codes are often (much) faster.

The GLPK API is simpler to use than CBC however.

In general I would say the low level solver API are often not a good idea to use for larger, more complex models. For those a modeling language is often more appropriate.

You may want to consider making your application solver independent. Then you can pick the solver most suited for the occasion. Again this is most easily done using a modeling framework.

Implementing Parsers

2011-09-06T00:26:00.001-04:00

During my transatlantic flight I was reading

Although it was a good read I expected a little bit more info about implementing funky languages and some advice about that. Actually this book is largely about parsing (parsing using the ANTRL parser generator to be more precise) and largely about traditional programming languages.

I usually write simple recursive-descent parsers directly (i.e. without the aid of a parser generator) as typically I only parse very simple things. in addition I feel a little more sure-footed this way as I have a little bit more control to handle an occasional dirty thing. So this is a little bit outside my direct interest. However, I learned a bit about memoizing parsers, a term I was not yet familiar with.

More on network models (2)

2011-08-22T17:10:00.001-04:00

Re: http://yetanothermathprogrammingconsultant.blogspot.com/2011/08/more-on-network-models.html.

What if you go bigger? Here some results with two of the leading commercial LP solvers (using default settings) compared to GAMS:

LP solution times are a little bit more unpredictable (in this case this is positive: somewhat quicker than thought by extrapolating from first picture in more-on-network-models.html), but GAMS execution time (including generation time) is still nicely linear.

Shipping problem

2011-08-22T10:36:00.001-04:00

I am starting to work on a problem that looks like a shipping problem I worked on earlier: garments (shirts) are sent from China to stores in the US either directly in a number of packages containing standard contents or through a US warehouse. The cheapest is to send the standard boxes from China but as there are just a few standard configurations, stores will need to be supplied with additional shirts from the US warehouse. Shirts come in different sizes.

The basic model looked like:

All quantities with bars are parameters (constants). Even though this is a simple model it was not completely trivial to get good solutions:

The model as presented is an MINLP. We linearized to make it a MIP.
The data set was rather large because of a large number of stores.
We could preprocess the demand data as some stores have the same demand; this reduced the size by about a third. We added weights to the objective to represent how many times a demand pattern was used.
In practice over-supplying a store was expensive, so we can fix OverShipped(Store,Size)=0.
It was not practical to solve this model in one swoop. We found good solutions by first finding a good content for a single standard box, and then (after fixing the content of this first box) finding a good proposal for a second standard box. Adding more standard boxes was not reducing the objective by much.
The set sizes can also include other attributes such as color etc.

Bin Packing

2011-08-15T13:55:00.001-04:00

For pure bin packing problems some really good heuristics exist. For the problem illustrated in http://www.developerfusion.com/article/5540/bin-packing/ one of the best solutions is:

So if we solve to optimality, do we gain a lot? From http://en.wikipedia.org/wiki/Bin_packing_problem is looks doubtful. Indeed when I solve:

$ontext

   Bin packing example

$offtext

sets
   b 'bins' /bin1*bin19/
   i 'items' /item1*item45/
;

scalar binsize /80/;
parameter itemsize(i) /
item1   26, item2   57, item3   18, item4    8, item5   45
item6   16, item7   22, item8   29, item9    5, item10 11
item11   8, item12 27, item13 54, item14 13, item15 17
item16 21, item17 63, item18 14, item19 16, item20 45
item21   6, item22 32, item23 57, item24 24, item25 18
item26 27, item27 54, item28 35, item29 12, item30 43
item31 36, item32 72, item33 14, item34 28, item35   3
item36 11, item37 46, item38 27, item39 42, item40 59
item41 26, item42 41, item43 15, item44 41, item45 68
/;

variables
   x(i,b) 'assignment'
   y(b)    'whether bin is used'
   z       'count used bins'
;
binary variable x,y;

* relax y: automatically integer
y.prior(b) = INF;

equations
   onebin(i)
   cap(b)
   binusd
   obj
   symmetry(b)
;

onebin(i)..    sum(b, x(i,b)) =e= 1;
cap(b)..       sum(i, x(i,b)*itemsize(i)) =l= binsize;
binusd(i,b).. y(b) =g= x(i,b);
obj..          z =e= sum(b, y(b));
symmetry(b+1).. y(b) =g= y(b+1);

model m1 /onebin,cap,binusd,obj/;
model m2 /onebin,cap,binusd,obj,symmetry/;
option optcr=0;
solve m1 minimizing z using mip;

we get an optimal solution with 17 bins.

Of course as discussed in http://yetanothermathprogrammingconsultant.blogspot.com/2011/07/scheduling-of-tv-advertisement-theory.html my real models are almost always much more complicated than shown in these stylized bin packing examples. In those cases a MIP-based solution algorithm can be even more worthwhile than a (simple) heuristic as there is more room for improvement. In the real case we dealt with a complicated multi-objective model where our model seems to outperform the heuristic in almost all objectives. This is the case even after the heuristic was improved using some fine-tuning after looking at our solutions.

More on network models

2011-08-01T16:52:00.001-04:00

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:

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

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:

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

…

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:

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

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:

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

My conclusions:

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

Scheduling of TV Advertisement: Theory vs Practice

2011-07-28T00:06:00.001-04:00

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:

The number of spots to broadcast is not known in advance (depends on expected ratings of a break)
There are complicated issues with spread over time
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)
Different contracts have different requirements, and these requirements are not always easily modeled
The model is huge (> 150K binary variables) for a monthly schedule
The input data (from a database) is large and complex (this is sometimes called “data-intensive”)
Numerous other things

Pyomo

2011-07-21T13:39:00.001-04:00

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

2011-07-21T12:48:00.001-04:00

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.

Dashboards in Powerpoint

2011-07-06T10:36:00.001-04:00

I have built several “dashboards” in Excel (e.g. http://yetanothermathprogrammingconsultant.blogspot.com/2010/01/excel-reporting.html), but here is a different approach: use Powerpoint.

http://office.microsoft.com/en-us/support/interactive-dashboards-using-powerpoint-and-excel-HA102298884.aspx?redir=0

Min sum or min num

2011-07-05T18:49:00.001-04:00

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.

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

GAMS: GDXDUMP and Equations

2011-07-05T17:17:00.001-04:00

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

It is produced by GDXDUMP
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

2011-07-05T15:59:00.001-04:00

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:

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.

How bumblebees tackle the traveling salesman problem

2011-07-01T00:07:00.001-04:00

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

Sorting in GAMS

2011-07-01T00:06:00.001-04:00

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

Cplex polishing

2011-06-30T07:45:00.001-04:00

For some difficult large models, the Polishing facility in Cplex can sometimes do extremely well. This is an example:

Elapsed real time = 874.90 sec. (tree size = 11.31 MB, solutions = 2)
    233   235     -158.9737 3551      621.0895     -159.3086   258914 125.65%
    234   236     -158.9448 3578      621.0895     -159.3086   258964 125.65%
    266   268     -158.9498 3578      621.0895     -159.3086   265158 125.65%
    306   308     -158.9338 3568      621.0895     -159.3086   271534 125.65%
    310   312     -158.9315 3511      621.0895     -159.3086   272354 125.65%
    312   314     -158.9489 3581      621.0895     -159.3086   272367 125.65%
    322   324     -158.9310 3574      621.0895     -159.3086   273646 125.65%
    346   348     -158.8857 3561      621.0895     -159.3086   277217 125.65%
    354   356     -158.8368 3532      621.0895     -159.3086   278697 125.65%
    370   372     -158.8741 3567      621.0895     -159.3086   281105 125.65%
Elapsed real time = 995.15 sec. (tree size = 49.71 MB, solutions = 2)

Starting condition for polishing reached (PolishAfterTime).
Starting solution polishing.

*   390+ 390                          360.3876     -159.3086   285698 144.20%
*   390+ 390                          231.7287     -159.3086   285698 168.75%
*   390+ 390                          100.2251     -159.3086   285698 258.95%
    390   392     -158.8957 3503      100.2251     -159.3086   285698 258.95%
    391   393     -158.9880 3601      100.2251     -159.3086   286729 258.95%
    392   394     -158.9935 3613      100.2251     -159.3086   286790 258.95%
    393   395     -158.9735 3574      100.2251     -159.3086   287216 258.95%
    394   396     -158.9763 3595      100.2251     -159.3086   287247 258.95%
    395   397     -158.9980 3601      100.2251     -159.3086   287328 258.95%
    396   398     -158.9806 3600      100.2251     -159.3086   287599 258.95%
    397   399     -159.0641 3615      100.2251     -159.3086   296751 258.95%

GUB cover cuts applied: 854
Clique cuts applied: 14
Cover cuts applied: 1540
Implied bound cuts applied: 20
Flow cuts applied: 922
Mixed integer rounding cuts applied: 67
Zero-half cuts applied: 1169
Gomory fractional cuts applied: 236

Root node processing (before b&c):
Real time             = 127.44
Parallel b&c, 8 threads:
Real time             = 7070.70
Sync time (average)   = 168.52
Wait time (average)   =    0.00
                          -------
Total (root+branch&cut) = 7198.14 sec.
MIP status(107): time limit exceeded

Although the relative gap looks bad, this is in fact a very good solution. The reason for the gap problem is that we are relatively close to zero in our objective (compared to solutions found by heuristics) and the fact that we are crossing zero in the gap calculation does not help.

Advanced basis with Cplex

2011-06-18T02:58:00.001-04:00

Usually Cplex behaves very predictable with LP’s that have an advanced basis to start from. But now I have something that puzzles me. Basically I do a solve from an optimal solution. This should normally run in just a few iterations. But with the default setting I see:

First model

Iteration log . . .
Iteration:     1   Dual objective     =             1.005320
Iteration:   581   Dual objective     =           296.751221
Iteration: 1247   Dual objective     =           511.623938
Iteration: 1681   Dual objective     =           514.072257
LP status(1): optimal

Restart

Iteration log . . .
Iteration:     1    Objective     =           515.104554
Perturbation started.
Iteration:   402    Objective     =           515.104554
Iteration:   806    Objective     =           515.104088
Iteration: 1107    Objective     =           515.103794
Iteration: 1406    Objective     =           515.103557
Iteration: 1684    Objective     =           515.103314
Iteration: 1939    Objective     =           515.103216
Iteration: 2186    Objective     =           515.103124
Iteration: 2434    Objective     =           515.103037
Iteration: 2672    Objective     =           515.102875
Iteration: 2951    Objective     =           515.102790
Iteration: 3210    Objective     =           515.102730
Iteration: 3472    Objective     =           515.102670
Iteration: 3709    Objective     =           515.102630
Iteration: 3943    Objective     =           515.102592
Iteration: 4199    Objective     =           515.102525
Iteration: 4445    Objective     =           515.102493
Iteration: 4692    Objective     =           515.102468
Iteration: 4933    Objective     =           515.102448
Iteration: 5181    Objective     =           515.102428
Iteration: 5431    Objective     =           515.102411
Iteration: 5689    Objective     =           515.102395
Iteration: 5962    Objective     =           515.102379
Iteration: 6224    Objective     =           515.102362
Iteration: 6487    Objective     =           515.102353
Iteration: 6765    Objective     =           515.102342
Iteration: 7041    Objective     =           515.102335
Removing perturbation.
LP status(1): optimal

The restart from an optimal basis takes more iterations than solving from scratch. This is with a default setting of LPMETHOD=0 (automatic). When I use LPMETHOD=1 (Primal Simplex) or LPMETHOD=2 (Dual Simplex) I see:

Restart with LPMETHOD=1

Primal crossover.
Primal: Fixed no variables.
Dual: Fixing 993 variables.
      992 DMoves: Infeasibility 6.76704726e+000 Objective 5.15104554e+002
      394 DMoves: Infeasibility 1.33823843e+000 Objective 5.15104554e+002
        0 DMoves: Infeasibility 7.17019911e-001 Objective 5.15104554e+002
Dual: Pushed 0, exchanged 993.

Iteration log . . .
Iteration: 1 Objective = 515.104554
LP status(1): optimal

Restart with LPMETHOD=2

Dual crossover.
Dual: Fixing 994 variables.
      993 DMoves: Infeasibility 2.07167616e-013 Objective 5.15104554e+002
      414 DMoves: Infeasibility 7.97140132e-014 Objective 5.15104554e+002
        0 DMoves: Infeasibility 0.00000000e+000 Objective 5.15104554e+002
Dual: Pushed 1, exchanged 993.
Primal: Fixed no variables.
LP status(1): optimal

These seem indeed to take fewer iterations to solve the restarted model. I always assumed LPMETHOD=0 would just choose primal or dual simplex, but apparently it also does other things, and in this case that is not helping.

Update

I did some more experiments. With LPMETHOD=0 I can now also achieve a restart if I set all variable levels to zero. I.e.

solve m minimizing z using lp;
solve m minimizing z using lp;

gives a bad restart, while:

solve m minimizing z using lp;
x.l(i,j) = 0;
z.l=0;
solve m minimizing z using lp;

gives me zero iterations for the restart. Something is not completely right. As Bo suggested in his comments, the reason may be related to Cplex trying to do a presolve on an advanced basis. (I have heard about some ancient solvers that tried a basis preserving presolve, may be more limited than what current presolvers do). The difference between LPMETHOD=0 and both LPMETHOD=1 and LPMETHOD=2 is still a mystery to me.

Update 2

The behavior can be partially explained by looking at the Cplex documentation describing the Advind parameter:

For continuous models solved with simplex, setting 1 (one) will use the currently loaded basis. If a basis is available only for the original, unpresolved model, or if CPLEX has a start vector rather than a simplex basis, then the simplex algorithm will proceed on the unpresolved model. With setting 2, CPLEX will first perform presolve on the model and on the basis or start vector, and then proceed with optimization on the presolved problem.

I was using the model from a GAMS environment. There the default is 2. The documentation for GAMS/Cplex is somewhat more limited than the above paragraph:

Advind (Integer)
Use an Advanced Basis. GAMS/Cplex will automatically use an advanced basis from a previous solve
statement. The GAMS Bratio option can be used to specify when not to use an advanced basis. The Cplex
option advind can be used to ignore a basis passed on by GAMS (it overrides Bratio).
(default = determined by GAMS Bratio)
0 Do not use advanced basis
1 Use advanced basis if available
2 Crash an advanced basis if available (use basis with presolve)

Looks like 1 is better for a very good basis, 2 may be better for models with a basis that is not very good (in that case presolving will pay off even when destroying some basis information).

Network formulation

2011-06-17T01:00:00.001-04:00

I was looking at a model containing underlying traffic network. Such a network can be quite elegantly formulated in GAMS, but apparently many users still have problems with this. The mathematical formulation is:

The GAMS formulation can look like:

$ontext
max flow network example

Data from example in
     Mitsuo Gen, Runwei Cheng, Lin Lin
     Network Models and Optimization: Multiobjective Genetic Algorithm Approach
     Springer, 2008

Erwin Kalvelagen,
Amsterdam Optimization,
May 2008

$offtext

sets
i 'nodes' /node1*node11/
source(i) /node1/
sink(i) /node11/
;
alias(i,j);

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

parameter capacity(i,j) /
   node1.node2 60,   node1.node3 60,   node1.node4 60,   node2.node3 30
   node2.node5 40,   node2.node6 30,   node3.node4 30,   node3.node6 50
   node3.node7 30,   node4.node7 40,   node5.node8 60,   node6.node5 20
   node6.node8 30,   node6.node9 40,   node6.node10 30, node7.node6 20
   node7.node10 40, node8.node9 30,   node8.node11 60, node9.node10 30
   node9.node11 50, node10.node11 50 /;

set arcs(i,j);
arcs(i,j)$capacity(i,j) = yes;
display arcs;

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

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

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;

MIP Gap not decreasing

2011-06-14T12:46:00.001-04:00

Here is a small fragment of a Cplex log:

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

   3186 2520     -666.6217 4096      956.6330     -667.2010 1313338 169.74%
   3226 2560     -666.6205 4097      956.6330     -667.2010 1323797 169.74%
   3266 2600     -666.6201 4095      956.6330     -667.2010 1335602 169.74%
Elapsed real time = 2801.61 sec. (tree size = 77.54 MB, solutions = 2)
* 3324+ 2656                         -125.5775     -667.2010 1363079 431.31%
   3334 2668     -666.5811 4052     -125.5775     -667.2010 1370748 431.31%
   3380 2714     -666.5799 4017     -125.5775     -667.2010 1388391 431.31%
   3422 2756     -666.5791 4011     -125.5775     -667.2010 1403440 431.31%

Usually we expect the gap to decrease while the MIP solver is working its way. Here we see a different behavior. Because we go through zero the gap% is actually increasing here.

I believe the gap calculation used by Cplex is:

Obviously this is sensitive w.r.t. to the sign of BestInteger. Of course in practice we don’t see this behavior very often (in many models BestInteger and BestNode have the same sign).

Writing solutions back into a database

2011-06-07T14:13:00.001-04:00

The issue of writing optimization results back into a database is often a little bit more complex than reading input data (see http://yetanothermathprogrammingconsultant.blogspot.com/2011/04/live-databases.html for some thoughts on reading data).

First with writing more things can go wrong: you need additional permissions, you can overwrite data, and conceptual there are some issues. Here are some questions and possibilities that come to mind:

Do we have a say in how the tables look like or do we need to follow some existing data model?
Update or Insert. Do we use an SQL INSERT or SQL UPDATE to populate solution tables?
Some databases have SQL extensions that implement an UPSERT. E.g. MySQL has INSERT INTO … ON DUPLICATE KEY UPDATE, or REPLACE INTO, Oracle and SQL Server have MERGE.
It may be better to erase the solution tables first.
Or always INSERT with a new key indicating the run id. A new solution gets a new run id.
For very large data some RDBMS systems have a BULK INSERT facility (SQL Server: BCP or BULK INSERT, SQL*Loader for Oracle).
We can write to temporary tables and deal with the problem of merging results into the real tables further down the road.
Just write a well-designed temp database (e.g. in Access) and let the IT people deal with getting the data into their systems.

In general I prefer 8, let others take responsibility for plugging things into the database. Note that this step can be expensive. A few weeks ago I was discussing a complex model where we needed to get the results into the RDBMS. To achieve a reasonable turnaround time, I was given a time limit of about an hour to come up with solutions for the model, while the system guy allocated 2 hours for himself to put the solution back into the database (much time is spent on checks).

Remote Desktop adventure

2011-06-05T22:50:00.001-04:00

Suddenly my remote desktop window went black. Even after reconnecting, I got again just a black screen. Took me a little while to figure out that with Ctrl-Alt-End I could send a Ctrl-Alt-Del to the remote session, so I could log out. That solved the problem: after logging on again everything was fine.

GAMS/Gurobi and Ctrl-C

2011-05-30T00:47:00.001-04:00

I am seeing some strange results when I interrupt a GAMS/Gurobi run with Ctrl-c (in the IDE using the Interrupt button). The GAMS link should give in that case the best integer solution found so far, basically identical to an iteration or time limit. However the solution seems wrong: many of the binary variables x(i,j) in a large model are off (they are reported as zero instead of one). I suspect this is a problem in the GAMS link (i.e. not in Gurobi itself). From what I can see the link it passing back levels of integer variables incorrectly in these circumstances (I have seen this before: integer variables should have basis status SUPERBASIC to prevent only a basis status being sent back).

Btw, for this (huge) model the Gurobi heuristics do a fantastic job finding good solutions:

     0     0 179791.544    0 6744 1511067.07 179791.544 88.1%     -   91s
H    0     0                    247836.12285 179791.544 27.5%     -   93s
H    0     0                    232489.21120 179791.544 22.7%     -   96s
     0     0 181272.276    0 7366 232489.211 181272.276 22.0%     - 397s
H    0     0                    226736.43968 181272.276 20.1%     - 399s
H    0     0                    218612.71830 181272.276 17.1%     - 404s

Just need to remember not to use the Interrupt button.

PS. Another minor issue with the link:

MIP Solution:        7259.007800    (-2147483648 iterations, 0 nodes)
Best possible:        6916.876514
Absolute gap:          342.131286
Relative gap:            0.047132

it has troubles with the iteration count (which actually is 120570, a number that should not overflow).

Optimal Spread (2)

2011-05-19T09:08:00.001-04:00

Is turns out we also have to distribute x(t)’s such that we follow a non-uniform parameter. In the previous section http://yetanothermathprogrammingconsultant.blogspot.com/2011/05/optimal-spread.html we used a function v(t) = 1 with cumulative function w(t) = t.

In this case we consider a generic parameter v(t), with its cumulative cousin w(t)=w(t-1)+v(t). To trace this function as closely as possible we can use:

E.g. when we use v(t)=t, then we see:

---- 77 PARAMETER result

v x

i1        1.000
i2        2.000
i3        3.000
i4        4.000
i5        5.000       1.000
i6        6.000
i7        7.000
i8        8.000
i9        9.000       1.000
i10      10.000
i11      11.000
i12      12.000       1.000
i13      13.000
i14      14.000       1.000
i15      15.000
i16      16.000       1.000
i17      17.000
i18      18.000       1.000
i19      19.000
i20      20.000       1.000

Optimal spread

2011-05-14T18:55:00.001-04:00

In a multi-objective MIP model I am working I have a binary variable x(t). This variable is mostly zero but now and then it is one. One of the goals is to achieve an even spread of the ones. I.e.

The first row has the ones clustered in the middle. The second one looks pretty good.

The following idea is based on how Q-Q plots are used to find distributions in statistical data (http://en.wikipedia.org/wiki/Q-Q_plot).

First we form a running sum:

The “predicted’ running sum for an evenly distributed x(t) can be expressed as tn/T where n is the total number of ones, t is the current index and T is the last t:

Now we try to keep the values s(t) as close as possible to the predicted values tn/T. I.e.:

Note that n is a variable in my model: I don’t know in advance how many x(t)’s are turned on. Indeed when I run this (with n=3) I see:

---- 42 VARIABLE x.L

t3 1.000, t7 1.000, t10 1.000

Here are the results when we compare the two solutions (click to enlarge):

Indeed this seems to confirm this approach could work, and we can add d(max)-d(min) as an extra objective to minimize.

Excel: Trace Dependents

2011-05-04T16:22:00.001-04:00

I often use the “Trace Dependents” tool in Excel to quickly see which formulas use a given cell. Here is an example where I was a bit overwhelmed:

Yet Another Math Programming Consultant

Comparing some modeling languages

network

Dinkelbach’s Algorithm

Example 1: Linear fractional programming problem

Example 2: an integer fractional programming problem

Reviewing NLP model

Savings from applying MIP models

Bad NLP modeling

Sequence dependent setup times: difficult MIPs

Optimization for Dummies

Logs in NLP

Student posts (2)

Analytics in Retail

Modeling is more difficult than many people think…

Matrix multiplication

MPS files: this does not look right

This does not compute

Recently read

MPS Files (2)

Rick Rosenthal: Ten Keys to Success in Optimization Modeling

HPC applications in the Amazon cloud

Sulum Optimization

Automatic Reformulations

Solver Options

Finding feasible schedules for a column generation algorithm

Slacks (with penalties) should not be integer variables

Still not proven optimal…

Update

Special case of integer cuts

MPS files

Sorting by MIP

MIQP: optimality tolerances vs objective scaling

Solving a large sparse system of equations with an LP solver

Update

MIP caveats

Ateji defunct: OptimJ now free

Integer cuts

Most popular programming languages

Time schedules

Cplex MIP and time limits

GPU computing on Amazon cloud

Big-M and Scaling

Open source MIP solvers

Implementing Parsers

More on network models (2)

Shipping problem

Bin Packing

More on network models

AMPL Version

GLPK

Scheduling of TV Advertisement: Theory vs Practice

Pyomo

Patent Application

Dashboards in Powerpoint

Min sum or min num

GAMS: GDXDUMP and Equations

Background.

Excel issue

How bumblebees tackle the traveling salesman problem

Sorting in GAMS

Cplex polishing

Advanced basis with Cplex

Update

Update 2

Network formulation

MIP Gap not decreasing

Writing solutions back into a database

Remote Desktop adventure

GAMS/Gurobi and Ctrl-C

Optimal Spread (2)

Optimal spread

Excel: Trace Dependents