Wednesday, May 18, 2016

Finding all solutions of a system of linear inequalities II

In this post (a followup to this one), I want to find all corner points related to these inequalities:
\[\boxed{\begin{align} &3 \le 5x_1+3x_2+3x_3+5x_4 \le 5 \\ &0 \le x_2+x_3 \le 1 \\ &0 \le x_1+x_4 \le 1 \\ &1 \le 3x_3+5x_4 \le 3 \\ &0 \le 5x_1+3x_2 \le 2 \\ &0 \le x_1,x_2,x_3,x_4 \le 1 \end{align}}\]
The approach I follow is based on this post. That is, we use binary variables to explicitly encode the LP  basis, and then solve a series of MIP models where we add cuts to the model to exclude the currently found basis. There are two main differences with the model shown  in the original post:

  • The current problem has no objective. This is not a major hurdle: we can add a dummy objective. The algorithm becomes actually slightly easier as we don't have to check if the objective starts to deteriorate.
  • The current problem has bounded variables. In the original problem we assumed we have positive variables without an upper bound. That meant we could use a single binary variable to encode the basis status: a variable is basic or non-basic. In the current case we could do the same thing and implement the bounds as explicit constraints. However, instead I choose for a slightly more sophisticated approach: we allow the basis status to assume three different values: basic, non-basic at lower bound and non-basic at upper bound. This is exactly how Simplex based linear programming solvers deal with bounded variables.

Positive variables/slacks

Let’s recapitulate what we did in the case of positive variables. When we have positive variables (or slacks) \(x_k\ge 0\), we can use the following basis encoding:
\[\beta_k = \begin{cases} 1 \> \text{if \(x_k\) is basic}\\
                                   0 \> \text{if \(x_k\) is non-basic} \end{cases}\]
If a variable is non-basic it should be equal to zero. This is implemented as
\[0 \le x_k \le M \beta_k\]
In addition we know that we have \(m\) (number of LP rows) basics:
\[\sum_k \beta_k = m\]

Bounded variables/slacks

To handle bounded variables (and slacks): \(\ell_k \le x_k \le u_k\), we introduce an index \(b=\{\text{B,NBLO,NBUP}\}\), where B means basic (variable is between its bounds), NBLO means nonbasic at lower bound (i.e. \(x_k = \ell_k\)) and NBUP means nonbasic upper (\(x_k=u_k\)). Now we can write:
\[\beta_{k,b} = \begin{cases} 1 \> \text{if \(x_k\) has basis status \(b\)}\\
                                                 0 \> \text{otherwise}\end{cases}\]
Of course each variable can assume only one (and exactly one) basis status:
\[\sum_b \beta_{k,b} = 1 \> \forall k\]
If a variable is non-basic, it should be at bound:
\[\begin{align} & \ell_k \le x_k \le \ell_k + (u_k-\ell_k)(1-\beta_{k,\text{NBLO}})\\
                     & u_k \ge x_k \ge u_k - (u_k-\ell_k)(1-\beta_{k,\text{NBUP}}) \end{align}\]
Again, we have \(m\) basics:
\[\sum_k \beta_{k,\text{B}} = m\]

Solution

After running the model we have a large number of feasible bases, 767 to be precise. Some partial results are listed below:

image

The top part are the values of the columns and rows. The bottom part shows the basis status. We see we get some duplicate values (they have a different basis, but the same values; this is called degeneracy in linear programming).

We can easily add some code to the model to find the unique values and filter out the duplicates:

image

There are only 20 unique corner point solutions.

We had to run 768 models for that (the last model was infeasible so we stopped), and every model became larger by one constraint. That also means models are becoming more expensive to solve. This is illustrated in this graph: both the number of Simplex Iterations needed to solve each model (orange line) and the solution times (blue line) show an increase over time.

image

Model

Below the complete model is listed:

*-----------------------------------
* original model equations
*-----------------------------------

set
  k   
'rows+columns' /1*5,r1*r5/
  j(k)
'columns'      /1*5/
  i(k)
'rows'         /r1*r5/
  bnd 
'bounds'       /lo,up/
;


variables
   x(j) 
'original decision variables'
   r(i) 
'row levels'
;
equations
  e1,e2,e3,e4,e5
;

e1.. r(
'r1') =e= 5*x('1')+3*x('2')+3*x('3')+5*x('4');
e2.. r(
'r2') =e= x('2')+x('3'
);
e3.. r(
'r3') =e= x('1')+x('4'
);
e4.. r(
'r4') =e= 3*x('3')+5*x('4'
);
e5.. r(
'r5') =e= 5*x('1')+3*x('2'
);

table
rowbound(i,bnd)
   
lo up

r1   3  5
r2      1
r3      1
r4   1  3
r5      2
;

r.lo(i) = rowbound(i,
'lo');
r.up(i) = rowbound(i,
'up'
);
x.lo(j) = 0;
x.up(j) = 1;

*-----------------------------------

* basis encoding using
* binary variables
*-----------------------------------

scalar m 'number of LP rows';
m =
card
(i);

set b 'basis status'

     
/NBLO  'nonbasic at lower bound'
      
NBUP  'nonbasic at upper bound'
      
B     'basic'/
;
binary variables
   bs(k,b)
'basis status of each row and column'
;

equations
  onebs(k)  
'only basis status is current'
  countb    
'number of basics = m'
  vlower(j) 
'variable is nb at lower'
  rlower(i) 
'row is nb at lower'
  vupper(j) 
'variable is nb at upper'
  rupper(i) 
'row is nb at upper'
;

onebs(k)..  
sum(b, bs(k,b)) =e= 1;
countb..    
sum(k, bs(k,'B'
)) =e= m;
vlower(j)..  x(j) =l= x.lo(j) + (x.up(j)-x.lo(j))*(1-bs(j,
'NBLO'
));
rlower(i)..  r(i) =l= r.lo(i) + (r.up(i)-r.lo(i))*(1-bs(i,
'NBLO'
));
vupper(j)..  x(j) =g= x.up(j) - (x.up(j)-x.lo(j))*(1-bs(j,
'NBUP'
));
rupper(i)..  r(i) =g= r.up(i) - (r.up(i)-r.lo(i))*(1-bs(i,
'NBUP'
));


*-----------------------------------

* cuts to prevent revisiting
* same basis
*-----------------------------------

set
   iter
'iterations' /iter1*iter1000/
   c(iter)
'used cuts'
;
parameter cc(iter,k,b) 'cut-coefficients';
scalar
nn;
nn =
card
(k);

equation
cuts(iter);
cuts(c)..
sum
((k,b),cc(c,k,b)*bs(k,b)) =l= nn-1;

c(iter) =
no
;
cc(c,k,b) = 0;

*-----------------------------------

* dummy objective
*-----------------------------------

equation edummy;
variable
dummy;
edummy.. dummy =e= 0;

*-----------------------------------

* model statement
*-----------------------------------

model Model1 /all/;
Model1.solvelink=%Solvelink.LoadLibrary%;
Model1.solprint=%Solprint.Quiet%;

option
optcr=0;
option
mip=cplex;

*-----------------------------------

* solve loop
*-----------------------------------

parameters
  v(k,iter)
'collect values'
  bas(k,b,iter)
'collect basis status'
;

scalar continue /1/;
loop
(iter$continue,
  
Solve
Model1 using mip minimizing dummy;
   continue$(Model1.modelstat<>%ModelStat.Optimal%) = 0;
  
if
(continue,
      cc(iter,k,b) = round(bs.l(k,b));
      c(iter) =
yes
;
      v(j,iter) = x.l(j);
      v(i,iter) = r.l(i);
      bas(k,b,iter) = bs.l(k,b);
   );
);

display
v, bas;

*-----------------------------------

* remove duplicate values
*-----------------------------------

parameters
  w(k,iter)
'unique values'
  found
/0/
;

set u(iter) 'subset of c';
u(iter) =
no
;

loop
(c,
 
if(card
(u)=0,
    u(c) =
yes
;
    w(k,c) = v(k,c);
 
else

     found = 0;
    
loop(u$(not found),
       found$(
sum
(k,abs(w(k,u)-v(k,c)))<1.0e-5) = 1;
     );
    
if(not
found,
        u(c) =
yes
;
        w(k,c) = v(k,c);
     );
  );
);

display
w;