I would guess there are problems where primal simplex is doing well. Apparently not enough of those.
Tuesday, January 12, 2016
Matlab primal simplex solver
I would guess there are problems where primal simplex is doing well. Apparently not enough of those.
Simulating a color-bar in Bokeh
Bokeh does not (yet) have built-in support for color-bars. Below we simulate this by adding a figure below the scatter plot.
Interesting to see how this is on the to-do list for a while: https://github.com/bokeh/bokeh/issues/1441. Especially intriguing is the statement that a constraint solver is used to place things on the canvas.
Tuesday, January 5, 2016
Book buyers price index
We have the well-known CPI index, but may be we also need a separate index for people who buy books.
Of course this also applies to journal papers. In this post I was encouraged to have a look at this paper:
$37.95 for 11 pages that is about $3.50 per page for a downloadable PDF file.
GAMS to Python (2)
A simple, convenient way to move model results from a GAMS model into Python for further postprocessing is via a SQLite database. This looks like a roundabout way but in practice it works quite well.
To export GAMS data to a SQLite database first save all data to a GDX file and then convert this to SQLite. This can be automated using two lines of GAMS code:
execute_unload 'results.gdx';
execute 'gdx2sqlite -i results.gdx -o results.db -fast';
There are some advantages to save just everything in the database:
- All the names (tables, columns) are the same as in the GAMS model.
- Everything is saved. If you need additional data, it is already there. No need to rerun the GAMS model.
- These steps are quite fast.
- A SQLite database is just a standalone file. You can email it to colleagues.
On the Python side it is quite easy to read data from an SQLite database:
The location and size of the circles are taken from the GAMS result set. Inspecting data in a tabular form is also very easy:
For exploratory work the ipython notebook interface is very convenient.
See also: http://yetanothermathprogrammingconsultant.blogspot.com/2014/09/gams-to-python.html
Sunday, January 3, 2016
A variant of the Lights Out game
Suppose we have a grid like this, with some binary values turned on:
This is our initial state s0. The goal is to turn all lights off, but there is a catch. When flipping a light switch all lamps in the same row and in the same column also flip. E.g. when we turn off light (1,1) we get:
Which lights do we need to switch to make it dark? We want to minimize the number of needed switches. See this post.
Here is a MIP model:
The idea is to add up the switches for each position (i,j) and then say: we want he final state to be even (i.e. 0,2,4,…). A simple way to require that a variable is even is to write it as 2*y where y is an integer variable. Interestingly, we are not interested in the value of y.
Note that we subtract s(i,j) because we double counted s(i,j) when adding up the switches in the same column and the same row.
The result is:
---- 39 VARIABLE s.L switch state
c1 c2 c3 c4
r1 1 1
r2 1 1 1
r3 1 1
r4 1 1
In the MIP model we apply these switches all at once, but if we do them one by one we would see:
Now, a real difficult question is: what is the initial state s0 that requires the most switches?
Update
Hint: just adding s0 as variables and maximizing is not what I am after here (that would be easy but somewhat uninteresting). We want to find the initial configuration that requires the largest optimal switching plan. I.e. we end up with a bilevel mip. That is a notorious difficult thing.
More info on mixed-integer bilevel programming can be found in J.F.Bard, “Practical Bilevel Optimization: Algorithms and Applications”, 1999. Now on sale at Amazon for $459!!
Constraint Solving and Planning with Picat
Interesting book about an interesting language. I once had a copy of Turbo Prolog; Picat looks and feels similar in many respects.
Saturday, January 2, 2016
Finding all optimal LP solutions using the Cplex Solution Pool
In http://yetanothermathprogrammingconsultant.blogspot.com/2016/01/finding-all-optimal-lp-solutions.html a MIP formulation was presented to enumerate all optimal LP bases. A loop construct was used to find new basic feasible solutions. To top things off a cutting plane technique was implemented to prevent rediscovering solutions we already knew about.
Here we show a GAMS model that uses the Cplex solution pool. This is a somewhat esoteric option in Cplex. For a discussion of this option see: http://orinanobworld.blogspot.com/2013/01/finding-all-mip-optima-cplex-solution.html. Using this technique will simplify the process as we now only have a single solve. The solution looks very similar. The only thing that is different is the order of the solutions.
Here we show a GAMS model that uses the Cplex solution pool. This is a somewhat esoteric option in Cplex. For a discussion of this option see: http://orinanobworld.blogspot.com/2013/01/finding-all-mip-optima-cplex-solution.html. Using this technique will simplify the process as we now only have a single solve. The solution looks very similar. The only thing that is different is the order of the solutions.
 |
| Solutions found with Cplex solution pool technology |
The table has the same layout as the one in the earlier post. The first part shows the values. First we have the values for the decision variables xi,j followed by the slacks si and sj. The last row is the objective z. The second part reports the basis status. This table has first the decision variables xi,j followed by the slacks si and sj. Remember that a 1 indicates ‘basic’.
GAMS Model
| $ontext |
Note: the solutions will be written to the GDX file solutions.gdx.
Performance
The Cplex solution pool shows some amazing performance. To test this we duplicate the above model tenfold by adding an index to all variables and equations. Essentially we create a staircase model:
In our version we have 10 of those blocks and they are all identical. This gives us a model with many, many optimal bases. We set a limit of 10k solutions (PopulateLim). To find all these 10x models takes just 23 seconds (total turnaround time of the whole GAMS job).
Reference
Emilie Danna, Mary Fenelon, Zonghao Gu, Roland Wunderling, “Generating Multiple Solutions for Mixed Integer Programming Problems," Proceedings IPCO '07 Proceedings of the 12th international conference on Integer Programming and Combinatorial Optimization, Pages 280 - 294, Springer, 2007
Finding all optimal LP solutions
I am regularly confronted with the question: My LP model has multiple optimal solutions. Can I retrieve all optimal LP solutions? Probably the correct answer on this question is: no. If there are really multiple solutions there are really infinitely many of them. A better (?) question is: Can I get all optimal basic feasible solutions (bfs)? (i.e. only the corner points).
Note that many modelers think this is a wrong question. The argument goes like this: If you have a preference for one solution above another, you should give the model an incentive to find the better solution (e.g. by adding more detail to the objective). If you did not specify such an incentive, you really say that you don’t care. Another argument that is often used is that the number of optimal bases may be large. For an interesting discussion see (4,5,6).
Enumerating all optimal bases is not so easy either and I don't know a simple algorithm to do this. A somewhat complicated, low-level algorithm is found in (1). Another idea is from (2): encode the basis using binary variables and use a MIP formulation to enumerate them. This last approach is especially appealling as we can use a high level modeling language to implement the algorithm. In the paper (2) a GAMS formulation is used. Here we will develop a very similar implementation.
Lets have a look at the transportation model from the GAMS model library. This is actually a model from (3) but the data is slightly different. It has two different solutions and different LP solvers may find either of these solutions. The basic model looks like:
First some accounting. For a general LP problem with n (structural) variables and m constraints, m rows or columns will be basic, and n rows or columns will be non-basic. If we ignore the objective, our transportation network model has I+J constraints and I*J variables (here I and J indicate the number of elements in sets i and j). This means in any Simplex basis solution we have NB=I+J basics and NNB=I*J nonbasics.
In this example all variables are non-negative: x ≥ 0. This means all non-basic variables are zero. All basic variables can assume any value ≥ 0.
To enumerate all optimal bases for this model, the first thing we do is: introduce explicit slack variables s. After adding these slacks, the equations look like:
The next step is to encode the basis. We introduce some binary variables β that indicate whether a variable is basic (one) or non-basic (zero). If non-basic we know the variable must be zero. (Note: we ignore here the more general case where a non-basic variable can be at lower or upper bound). We also know the number of β's that should be one.
Finally we need to add some dynamic cuts to the model to forbid earlier discovered basic feasible solutions. This cut looks like:
Here α are coefficients that we set in our algorithm. Any time a variable is basic in an optimal solution, its corresponding coefficient is set to one.
The general way to cut off an existing 0-1 solution is derived here: http://yetanothermathprogrammingconsultant.blogspot.com/2011/10/integer-cuts.html. However in this case we can use a special simplified version. This is shown here: http://yetanothermathprogrammingconsultant.blogspot.com/2011/10/special-case-of-integer-cuts.html.
The model is done. Now we need a small loop to drive this thing. Basically the algorithm looks like:
- Solve the MIP model
- If not solved to optimality : STOP(ERROR)
- If objective started to deteriorate: STOP(DONE)
- Add a cut by augmenting the set cut and the parameter α.
- Go to step 1.
Below are the optimal solution values and the optimal bases.
The reportv table shows the values. First we have the values for the decision variables xi,j followed by the slacks si and sj. The last row is the objective z. The basis table reportb has first the decision variables xi,j followed by the slacks si and sj. Remember that a 1 indicates ‘basic’. Indeed all columns have 5 basic entries as the problem has 5 equations. We see there are eight different optimal bases but they correspond to only two different solutions. This illustrates that the same solution can be encoded by different bases. The z row in the values report confirms these are all optimal solutions with the same objective.
Note that GAMS by default returns a basis for MIP models. They do this as follows. After a MIP was solved, fix the integer variables to their current values and resolve the model as an LP. Return the duals of this LP. This GAMS basis is not used in this discussion. That means that this algorithm can be used with any MIP solver and can be directly translated into different environments (e.g. R, Matlab, Python etc)
A version of this model using the Cplex solution pool is here: http://yetanothermathprogrammingconsultant.blogspot.com/2016/01/finding-all-optimal-lp-solutions-using.html.
References
(1) Ralph E. Steuer, Multiple Criteria Optimization: Theory, Computation, and Application, Wiley, 1986.
(2) Sangbum Lee, Chan Phalakornkule, Michael M. Domach, Ignacio E Grossmann, Recursive MILP model for finding all the alternate optima in LP models for metabolic networks, Computers & Chemical Engineering, Volume 24, Issues 2–7, 15 July 2000, Pages 711-716,
(3) George B. Dantzig, Linear Programming and Extensions, Princeton University Press, 1963.
(4) Quirino Paris, "Multiple Optimal Solutions in Linear Programming Models." Amer. J. Agr. Econ. 63 (1981):724-727
(5) Bruce A. McCarl and Carl H. Nelson, “Multiple Optimal Solutions in Linear Programming Models: Comment,” Amer. J. Agr. Econ. 65 (1983):181-183
(6) Quirino Paris, "Multiple Optimal Solutions in Linear Programming Models: Reply,” Amer. J. Agr. Econ. 65 (1983):184-186
Appendix: the complete GAMS model
| $ontext |
Friday, January 1, 2016
Internet Explorer WebGL bug
I was doing some experiments with a 2d scatter plot with more than 5k points (see http://yetanothermathprogrammingconsultant.blogspot.com/2015/12/visualization-of-large-multi-criteria_18.html for more information). The Bokeh visualization tool allows to speed things up by using WebGL in the browser. This will allow some graphics related operations to be executed on the GPU (Graphics Processing Unit). For some operations this makes quite some difference in speed. However, this showed a bug in MS Internet Explorer (or MS Edge which I am using). The correct plot after selecting some points is on the left. MS IE draws all points however with a size of one pixel if WebGL is used. Chrome does the right thing (whether WebGL is used or not). MS IE also works correctly if the plot is created without WebGL support.
The Bokeh docs seem to confirm this: "Making a selections of markers on Internet Explorer will reduce the size of the markers to 1 pixel (looks like a bug in IE)."
|
|
Subscribe to:
Posts (Atom)