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

Optimization for Dummies

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

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)

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

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…

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

• 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

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

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

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)

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

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

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

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

Update: Unfortunately Sulum Optimization Tools are discontinued as per June 2015.

Automatic Reformulations

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

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

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

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…

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 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. … 129959  4046        0.0003    10        0.0031        0.0003   113627   91.54% *165794  5594      integral     0        0.0027        0.0005   146606   81.09% *170030  4859      integral     0        0.0019        0.0005   149634   72.63% *193953  3702      integral     0        0.0007        0.0006   171409   15.85% *195840   418      integral     0        0.0006        0.0006   172394    2.03% 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

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

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.