Saturday, January 2, 2016

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:

  1. Solve the MIP model
  2. If not solved to optimality : STOP(ERROR)
  3. If objective started to deteriorate: STOP(DONE)
  4. Add a cut by augmenting the set cut and the parameter α. 
  5. 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

 
Enumerate all optimal LP bases

$offtext


Sets
     k                      
/seattle, san-diego, new-york, chicago, topeka/
     i(k)  
canning plants   / seattle, san-diego /
     j(k)  
markets          / new-york, chicago, topeka /
;

Parameters

     a(i) 
capacity of plant i in cases
      
/    seattle     350
           
san-diego   600  /

     b(j) 
demand at market j in cases
      
/    new-york    325
           
chicago     300
           
topeka      275  / ;

Table d(i,j)  distance in thousands of miles

                 
new-york       chicago      topeka
   
seattle          2.5           1.7          1.8
   
san-diego        2.5           1.8          1.4  ;

Scalar f  freight in dollars per case per thousand miles  /90/
;

Parameter c(i,j)  transport cost in thousands of dollars per case
;

          c(i,j) = f * d(i,j) / 1000 ;
Scalars

     nb 
'number of basic variables'
     nnb
'number of non-basic variables'
  ;
nb =
card(i)+card(j);
nnb =
card(i)*card
(j);

parameter
Mx(i,j), Mk(k);
Mx(i,j) = min(a(i),b(j));
Mk(i) = a(i);
Mk(j) = b(j);

set
bs(*,*);
bs(i,j)=
yes
;
bs(
'-',k) = yes
;

set n /iter1*iter100/
;
set cn(n) 'used cuts'
;
parameter
alpha(n,*,*);
cn(n) =
no
;
alpha(n,bs) = 0;

Variables

     x(i,j) 
shipment quantities in cases
     z      
total transportation costs in thousands of dollars
     s(k)      
'slacks'
     beta(*,*) 
'basis'
;

Positive Variable x,s;
Binary Variable
beta;


Equations

     cost       
define objective function
     supply(i)  
observe supply limit at plant i
     demand(j)  
satisfy demand at market j
     basisx(i,j) 
'x=nb => x=0'
     basiss(k)   
's=nb => s=0'
     basis       
'basis accounting'
     cut(n)      
'0-1 cuts'
;

cost ..        z  =e= 
sum((i,j), c(i,j)*x(i,j)) ;
supply(i) ..  
sum
(j, x(i,j)) + s(i) =e=  a(i) ;
demand(j) ..  
sum
(i, x(i,j)) - s(j) =e=  b(j) ;

basis..       
sum
(bs,beta(bs)) =e= nb;
basisx(i,j)..  x(i,j) =l= Mx(i,j)*beta(i,j);
basiss(k)..    s(k) =l= Mk(k)*beta(
'-'
,k);


cut(cn)..     
sum
(bs, alpha(cn,bs)*beta(bs)) =l= nb-1;


Model transport /all/
;

scalar
zopt;
scalar continue /1/
;
parameter
reportv(*,*,*,*);
parameter
reportb(*,*,*,*);
option
reportv:2:3:1;
option
reportb:0:3:1;

option
optcr=0;
option
mip=cplex;
transport.solvelink=%Solvelink.LoadLibrary%;
transport.solprint=%Solprint.Quiet%;

loop
(n$continue,
  
Solve
transport using mip minimizing z ;
   zopt$(
ord
(n)=1) = z.l;
   continue$(transport.modelstat<>1
or
z.l > zopt + 0.0001) = 0;
  
if
(continue,
      alpha(n,bs) = round(beta.l(bs));
      cn(n) =
yes
;
      reportv(
'x'
,i,j,n) = x.l(i,j);
      reportv(
's','-'
,k,n) = s.l(k);
      reportv(
'z','-','-'
,n) = z.l;
      reportb(
'x'
,i,j,n) = beta.l(i,j);
      reportb(
's','-',k,n) = beta.l('-'
,k);
   );
);

Display
reportv, reportb;