Comparing some modeling languages

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:

1. Is there support for non-linear programming 
2. Performance characteristics for very large models
3. Ability to do data manipulation efficiently
4. 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)
5. 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

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

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

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

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

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

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:

1. 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.
2. The initial solution was constructed by ordering the jobs by due date.
3. The Gantt chart for the machine can look like:
   
   Here the late jobs are colored red.
4. The model minimizes sum of lateness tardiness. Additional objectives could be number of late jobs and total makespan.
5. More difficult to model is this: scheduled downtime. If jobs can not be interrupted, this causes holes (idle time) in the schedule.
6. A different, slightly more complicated formulation would result from:
   x(i,j) = 1 if job i is scheduled immediately before job j 
              0 otherwise
7. 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).
8. Many papers suggest a meta-heuristic (genetic algorithm etc.) to solve this.
9. Of course a MIP model can still help to define the problem correctly and to solve small instances to optimality.
10. 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.
11. 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).
12. 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).