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.
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 |
No comments:
Post a Comment