Wednesday, May 18, 2016

Finding all solutions of a system of linear inequalities III

An easier method to find all corner points for a system of linear inequalities is to encode the basis using binary variables and then use the Cplex solution pool. Here we used a loop where we added binary cuts to prevent an earlier found basis to be revisited again. A complete description of the problem can be found in that post also. With the Cplex solution pool we just use one solve and we don’t have to bother with specifying the cuts. In addition we expect to see much of a performance improvement. So how does performance stack up with the loop?

image

As we can see the performance of the solution pool is superior compared to the loop implementation. Note: it will give the same solution (although in a different order).

Model

The complete model listing is here:

*-----------------------------------
* 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'
));

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

* 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;

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

* solution pool solve
*-----------------------------------


$onecho > cplex.opt
SolnPoolAGap = 0.0
SolnPoolIntensity = 4
PopulateLim = 10000
SolnPoolPop = 2
solnpoolmerge solutions.gdx
$offecho

model1.optfile=1;
Solve Model1 using mip minimizing dummy;


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

* read solution data
*-----------------------------------

sets
  p
'solution points' /soln_model1_p1*soln_model1_p1000/
  c(p)
;
parameters
  v(p,k)
'collect values'
  bas(p,k,b)
'collect basis status'
;
execute_load "solutions.gdx",v=x,bas=bs,c=index;


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

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

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

set u(p) 'subset of index';
u(p) =
no
;

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

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

display
w;