Good polling is difficult

From the Huffington Post, a surprisingly detailed exposé on non-response in polling on the US precedential race: http://www.huffingtonpost.com/2012/06/17/gallup-poll-race-barack-obama_n_1589937.html

Multiple solutions

Many questions are posed about finding a number or all optimal solutions of an LP / MIP problem. Of course this is not always a well-defined problem especially for the LP case. One could reinterpret the question as: give me a number of different optimal bases, but even that may give a large number of identical primal solutions (degeneracy). Probably what is meant is: give me different optimal bases where the primal solution is different.

For a MIP probably what is meant: give me optimal (or close to optimal) solutions that have different integer configurations. I.e. ignore any continuous variables that are possibly present in the model.

Modelers have been struggling with this for a long time. Here are some notable papers:

• Chemical Engineering Application:

  Recursive MILP model for finding all the alternate optima in LP models for metabolic networks
  by: S. Lee, C. Phalakornkule, M. M. Domach, I. E. Grossmann
  Computers and Chemical Engineering, Vol. 24 (2000),

• Agricultural Economics:
  
  Multiple Optimal Solutions in Linear Programming Models, by Quirino Paris, American J. of Agricultural Economics, vol 63 (1981). Also see the comments on this paper by McCarl e.a. and others and the rebuttals.

For a MIP some solvers have tools to help with this. Cplex has a solution pool facility (see the posts http://orinanobworld.blogspot.com/2012/05/k-best-solutions-in-amplcplex.html and http://orinanobworld.blogspot.com/2012/04/k-best-solutions.html). In addition there are interesting modeling tricks to forbid an existing integer solution (but that would require to repeatedly solve a model).

I often believe that asking for many solutions may be the wrong question. If you want to present a (large) number of solutions and let the user pick the best one it may be better to investigate how this works, and incorporate this in the model.

As an example, scheduling models often have (very) different solutions with almost the same objective. It may make sense to choose then a schedule that is close to the schedule calculated in a previous run. To keep solutions close to each other -- unless there is a big cost involved -- is sometimes called “persistency”. It often make sense to add a facility for this in scheduling models, as operators don’t like to look at schedules that change a lot without reason. I don’t see this aspect mentioned much in the literature

Solver Foundation Enterprise

I would have expected a more formal announcement than just this:

http://social.msdn.microsoft.com/Forums/en-US/solverfoundation/thread/aedcbb42-d3df-4180-bbaf-4db2e66acd9d

Looks to me that the best way forward is to open source the Solver Foundation code. Of course assuming there are no third-party components used directly in the Solver Foundation code base (third party solvers can be easily excluded).

Update: Gurobi has retracted the original statement. It still would be good to have some more visibility on the status of MSF.

Update 2: from http://social.msdn.microsoft.com/Forums/en-US/solverfoundation:

As users have pointed out, Microsoft has not been active on the Solver Foundation forums since Nate left.   We have been quiet while we have gone through restructuring and planning.  Some would say we have been too quiet.  We know we have very loyal and enthusiastic users who want to know the future of Solver Foundation.  So, here is a long overdue statement about our plans for Solver Foundation.

The current 3.1 release of MSF will be the last release as a standalone install.  We are working hard on integrating Microsoft Solver Foundation into a larger analytics framework that will help users build both prescriptive and predictive analytics.  We look forward to releasing this new product for your use as soon as we are able to do so.  This new product will provide a migration path for current Solver Foundation users and partners.  

We would like to continue to keep the current forum open to the community to discuss MSF until the release of the new product. However, Microsoft will be providing limited support of MSF in terms of monitoring the forums and providing bug fixes during this transition time.                                                                                                                               

We have been responding to email and will continue to do so.  If you have feedback on issues/bugs/improvements, we welcome your feedback via msfsupport@microsoft.com.  Please check back on the forum for future announcements with regards to the new analytics product.

Thanks for your support.

The Solver Foundation Team

I am not 100% what this means. Looks like it is going away from a somewhat API centered product largely aimed to programmers who also do want to do some optimization, towards something that has more enterprise umpf (“analytics framework”). I assume this means moving more in the direction of SAS and Fico, i.e. more emphasis on other things such as data base, statistics, data mining than on pure optimization.

Model improvements

Looking at the models used in the paper http://www.ime.usp.br/~egbirgin/publications/rob.pdf I noticed that a few standard “improvements” could be introduced. Here is an example:

First equations (13) and (14) have the same RHS. Hence it makes sense to introduce a new variable for this.

Secondly, when we look at the AMPL models provided by the authors, we see:

R1 {j in 1..n}:     T[j] >=   ( S[j,m] + sum{i in 1..n} ( x[i,j] * p[i,m] ) )             - sum {i in 1..n} ( x[i,j] * d[i] ); R2 {j in 1..n}:     E[j] >= - ( S[j,m] + sum{i in 1..n} ( x[i,j] * p[i,m] ) )             + sum {i in 1..n} ( x[i,j] * d[i] );

I.e. they substituted out C[j,m]. I would probably leave that as it was formulated in the paper.

Now, do these improvements really work? As is sometimes the case: the results on a single instance can be misleading. So here were my experiments:

		old version	new version
size	rows	335	500
	cols	405	570
	nz	5635	3865
gurobi	iterations	1060935	453173
	nodes	31220	5971
	time	60	38
cbc	iterations	4589173	4790425
	nodes	36154	49851
	time	395	478

We see that my version has fewer nonzero elements at the expense of more variables and equations. But not with every solver (or with every instance) this pays off. CBC actually becomes slower. Most likely this can be explained with the simple argument that the solver follows a different path. To really assert that a formulation is better one really needs to run a large number of models and see how they behave on average.

Note: the reformulated model (and the computational results) already included the suggestion in the comment from Paul. 

Gurobi version 5.0 is out

New Features and Enhancements

Some of the new additions and enhancements include:

• New QCP and MIQCP solvers - Allow you to solve quadratically constrained models at more than twice the speed of the leading competitor
• New MATLAB and R interfaces - Concise and tightly integrated for ease of use
• Support for Lazy Constraints - Uses a superior implementation that doesn't require you to turn off other key features
• Enhancements for handling infeasible and unbounded problems - Includes a barrier homogenous algorithm designed to solve infeasible and unbounded problems and a new FeasRelax procedure for finding solutions that minimize constraint violations
• Enhancements to overall performance - Includes faster presolve times, the ability to warm start the simplex algorithms from solution vectors, enhanced crossover numerical stability, and many more enhancements for speed and robustness

For more info see gurobi.com.

Binary OR

Is it possible to implement a kind of a boolean addition to binary variables.

For example:

Let X, Y and Z be binary variables.

X and Y ara inputs, Z is a result.

X + Y = Z

Such that

0 + 0 = 0

0 + 1 = 1

1 + 0 = 1

1 + 1 = 1

So this would be sometning like a logical .OR. addition?

This can be formulated as a set of inequalities:

The variable Z can be relaxed to a continuous variable between 0 and 1 (it will be integer automatically).

Longest flow shop (2)

Permutations in GAMS

This is a follow up to http://yetanothermathprogrammingconsultant.blogspot.com/2012/04/longest-flow-shop.html.

We consider again the small data set:

table proctime(m,j) 'processing times'             job1  job2  job3  job4  job5  job6  job7  job8  job9  job10   machine1   4     2     1     5     4     3     5     2     1     8   machine2   2     5     8     6     7     4     7     3     6     2   machine3   2     8     6     2     2     4     2     7     1     8 ;

This data set comes from:

The first approach to find the permutation with the longest schedule is to enumerate all permutations. GAMS has a built-in option to generate permutations of a set. It is not completely intuitive, in part because the way GAMS stores sets. All sets have fixed ordering with respect to a pool of set elements (the universe). To create a different ordering, an extra index is used. E.g. for 3 jobs we have:

$set n 3 $eval k fact(%n%) set    j           'jobs'                    /job1*job%n%/    p           'permutations'            /p1*p%k%/    allp(p,j,j) 'store all permutations'    m           'processing stage'        /machine1*machine3/ ; option allp>j; option allp:0:1:1; display allp;

The option > creates the permutations. Of course the syntax is not very readable: one would expect a function with a name somehow related to ‘permutation’. One of the major goals when using a modeling system is to produce readable models, and that is certainly not achieved with this syntax.

The result is:

----     13 SET allp  store all permutations INDEX 1 = p1             job1        job2        job3 job1         YES job2                     YES job3                                 YES INDEX 1 = p2             job1        job2        job3 job1         YES job2                                 YES job3                     YES INDEX 1 = p3             job1        job2        job3 job1                     YES job2         YES job3                                 YES INDEX 1 = p4             job1        job2        job3 job1                     YES job2                                 YES job3         YES INDEX 1 = p5             job1        job2        job3 job1                                 YES job2         YES job3                     YES INDEX 1 = p6             job1        job2        job3 job1                                 YES job2                     YES job3         YES

Although I don’t like the syntax, it is actually quite fast. If we set n to 10 and disable the display statement – it would take a long time to print all these permutations – we can generate all 10!=3,628,800 permutations in less than 2.5 seconds. This compares quite good to the following R session:

We can create and evaluate the schedules as follows:

loop(p,    putclose$(mod(ord(p),1000) = 0) "iteration " ord(p):0:0/;    completion(m,j) = 0;    loop((m,allp(p,jj,j)),       completion(m,jj) = max(completion(m,jj-1),completion(m-1,jj)) + proctime(m,j);    );    obj(p) = smax((m,j),completion(m,j)); );

Unfortunately this loop is very slow, it took almost 7 hours on my machine.

Results

The following picture shows the three schedules we discussed:

• optimal
• initial schedule: formed by sequence job1,job2,job3,…
• a longest schedule

A MIP Formulation

Of course with a little thought we can improve on this solution technique and write a MIP model:

* * assignment constraints * assign1(j).. sum(jj, position(j,jj)) =e= 1; assign2(j).. sum(jj, position(jj,j)) =e= 1; processing(m,j).. completion(m,j) =e= start(m,j) + sum(jj,position(jj,j)*proctime(m,jj)); * * start(m,j) = max(completion(m,j-1),completion(m-1,j)) * max1(m,j)$(ord(j)>1).. start(m,j) =g= completion(m,j-1); max2(m,j)$(ord(m)>1).. start(m,j) =g= completion(m-1,j); max3(m,j)..            start(m,j) =l= completion(m,j-1)+H*b(m,j); max4(m,j)..            start(m,j) =l= completion(m-1,j)+H*(1-b(m,j)); start.fx('machine1','job1')=0; cmaxdef.. cmax =e= completion('machine3','job10'); model flowshop /all/; option optcr=0; solve flowshop maximizing cmax using mip;

It is not 100% trivial as we need to make sure the max is attained exactly. We use the following set of inequalities:

In the case of minimization as formulated in http://yetanothermathprogrammingconsultant.blogspot.com/2012/04/longest-flow-shop.html we just could use bounds as the minimization would make sure that

start(m,j) = max(completion(m,j-1),completion(m-1,j))

would hold exactly (for the places where it matters – if there is slack in some machine this relation may not hold).

This model would look much simpler with a constraint programming formulation. There we could have used constructs like all-different and max to simplify things. Many modern modeling languages support constraint programming (unfortunately GAMS does not). Also note that AMPL using Cplex indicator variables allows this to be formulated quite elegantly (see http://yetanothermathprogrammingconsultant.blogspot.com/2009/07/formulation-ciminaibi.html). GAMS does not directly support indicator variables in the language, but instead some option files must be used.

The above MIP model solves in 0.2 seconds, clearly beating our brute force method.

Longest flow shop

I was trying to demonstrate some simple scheduling models. The permutation flow shop has an obvious first solution: run jobs in the order given in the problem data: job1,job2,job3,…. We compared that with the optimal solution:

The model I presented is fairly simple:

Someone came up with the question: “what would the picture look like for the worst permutation?”.

I think this is actually more difficult to model. We assume in several places that jobs are automatically pushed to the left and if we just start to maximize this won’t work anymore.

Of course just trying all solutions would be feasible in this case. We have 10 jobs so 10! = 3,628,800 permutations (just type 10! in google).

Model is infeasible or unbounded

My favorite return code from an LP/MIP solver is:

MIP status(4): Model was proven to be either infeasible or unbounded.

Come on guys, you really can do better than this! OK, what to do as a user? My simplest suggestion – using GAMS – is to add a bound on the objective variable and solve again:

z.lo=−1e10;
solve m minimizing z using mip;

Now at least we see:

MIP status(3): Model was proven to be infeasible.

If the model was actually unbounded you would have seen large numbers in the solution, making it easy to find out what was wrong.

Note: Please don’t tell me “dual infeasible” is a good return code. That is just another name for the same thing, and just as useless for a modeler. If you are working on a model you really want to know if the model is infeasible or unbounded. I understand it is an easy way out for the algorithm-developer to return something like “dual infeasible” (I have heard some users saying that is just laziness, and I understand that sentiment). The return code should really relate to the model and not what is the easiest for the algorithm-developer to implement.

Not everyone agrees with these opinions (see the comments).

The model was fairly small. If we compare the time both the user needed to email me the question and the model, for me to look at it and respond (say in total half an hour) and the time the solver would need to spend to get a better diagnostic (let’s say 0.1 seconds) then we can say a huge efficiency gain is possible here.

Parallel GAMS jobs (2)

In http://yetanothermathprogrammingconsultant.blogspot.com/2012/04/parallel-gams-jobs.html I described a simple approach I suggested to a client allowing to run multiple scenarios in parallel.

For a different client we needed to run a randomized algorithm that solves many small MIP models. They are so small that using multiple threads inside the MIP solver does not give much performance boost (much of the time is spent outside the pure Branch & Bound part – such as preprocessing etc.). However as the MIP problems are independent of each other we could generate all the necessary data in advance and then call the scenario solver (http://www.gams.com/modlib/adddocs/gusspaper.pdf). This will keep the generated problem in memory, and does in-core updates, so we don’t regenerate the model all the time.

When running the algorithm with n=100,000 MIP models we see the following performance. Note that besides the MIP models there is also a substantial piece of GAMS code that implements other parts of the algorithm.

Implementation	number of MIP models	solve time	rest of algorithm	total time
Traditional GAMS loop (call solver as DLL)	100,000	1068 sec	169 sec	1237 sec
Scenario Solver	100,000	293 sec	166 sec	459 sec

To get more performance I tried to run the scenario solver in parallel. That is not completely trivial as the solver has a number glitches (e.g. scratch files with fixed, hard coded names). I also run parts of the GAMS algorithm in parallel, but some parts had to be done in the master model after merging the results.

Implementation	number of MIP models	Worker threads	parallel sub-problem time	rest of algorithm (serial)	total time
Parallel + Scenario Solver	100,000	4	116 sec	67 sec	183 sec

The implementation does not win the beauty contest, but it could be developed quickly. For these larger problems you always have to watch out for performance bottlenecks. One wrong line in the GAMS model and the performance drops to hours of computation time. In addition the GAMS tools here are not very refined, but if implemented correctly quite effective. As an example consider how I merge the results of the sub jobs. Each sub job writes a results file (GDX file) into its own directory and the master model will read those and merge them together so we have a single solution set for further processing:

* * read and merge the GDX files * loop(job,    put_utility 'gdxin' / 'jobdir',(jobinfo(job,"jobno")):0:0,"\results.gdx"/;    execute_loadpoint solpoollevel,solpoolobj,pp; ) ;

This is not exactly readable and intuitive code (caused by lack of design when these features were added to the GAMS language in a hurry), but it is quite fast.