Generating all non-dominated solutions in a multi-objective integer programming model

In http://ideas.repec.org/a/eee/ejores/v158y2004i1p46-55.html a simple MIP formulation is shown that enumerates all non-dominated solutions in a multi-objective integer programming model. As usual, it is instructive to see if we can reproduce the results. It was not difficult to verify the results by coding the algorithm in GAMS. Actually in GAMS it is quite easy to add the additional constraints using dynamic sets. I am sure this is much more compact (and easier to write) than the C code mentioned in the paper.

In this experiment we test the simple MOIP:

max x₁ − 2x₂
max -x₁+3x₂
x₁ − 2x₂ ≤ 0
x1,x2 ∈ {0,1,2}

$ontext

   This algorithm finds non-dominated solutions of an integer programming model.

   Erwin Kalvelagen, erwin@amsterdamoptimization.com, feb 2010.


   Reference:
      Sylva, Crema: "A method for finding the set of non-dominated vectors for multiple
      objective integer linear programs", EJOR 158 (2004), 46-55.

$offtext



set i 'variables' /i1,i2/
    k 'objectives' /k1,k2/
;

parameters
  f(k) 'minimal improvement' /k1 1, k2 1/
  lambda(k) 'initial weights' /k1 4, k2 3/
  M(k) 'big-M' /k1 4, k2 2/
  xup 'upperbounds' /i1 2, i2 2/
  ;

table C(k,i) 'cost coefficients'
       i1  i2
   k1   1  -2
   k2  -1   3
;


integer variables x(i);
x.up(i) = xup(i);

variables zw,z(k);

equations
   obj(k),e,wobj;

obj(k).. z(k) =e= sum(i, c(k,i)*x(i));

e.. x('i1')-2*x('i2') =L= 0;

wobj.. zw =e= sum(k, lambda(k)*z(k));

*
* this mip is just for checking (we don't use it in the algorithm)
*
option optcr=0;
model m1 /obj,wobj,e/;
solve m1 using mip maximizing zw;

parameter sols(*,*);
sols('MIP',i)=x.l(i);
sols('MIP',k)=z.l(k);
display sols;


*
* here we use the algorithm: find solutions that are not worse in all objectives
* than previously found solutions.
*

set
  tall /t1*t100/
  t(tall)  'dynamic set'
;



binary variable y(tall,k);
equations
  nondominance(tall,k)  'solution not dominated by earlier solutions'
  county(tall)          'count number of y-s'
;
parameter zmem(tall,k) 'memory';

nondominance(t,k).. z(k) =g= (zmem(t,k)+f(k))*y(t,k) - M(k)*(1-y(t,k));
county(t).. sum(k, y(t,k)) =g= 1;

model m2 /obj,wobj,e,nondominance,county/;


t(tall) = no;
zmem(tall,k) = 0;

option limrow=0;
option limcol=0;
m2.solprint=2;
m2.solvelink=2;

alias (iter,tall);
scalar ok /1/;
loop(iter$ok,
   solve m2 using mip maximizing zw;
   if (m2.modelstat = 1,

      t(iter) = yes;
      zmem(iter,k) = z.l(k);

      sols(iter,i)=x.l(i);
      sols(iter,k)=z.l(k);

   else
      ok = 0;
   );
);


display sols;

The reported five solutions correspond to the results reported in the paper:

----    114 PARAMETER sols               i1          i2          k1          k2 MIP       2.000       2.000      -2.000       4.000 t1        2.000       2.000      -2.000       4.000 t2        2.000       1.000                   1.000 t3        1.000       2.000      -3.000       5.000 t4        1.000       1.000      -1.000       2.000 t5                    2.000      -4.000       6.000

(The line with MIP results can be ignored as that is not part of the algorithm).

This method can work using a good MIP solver with moderate size problems.

Some notes:

• We can streamline this a bit for zero-one models. Especially in combination with cuts that prevent same solutions.
• In that case we just can do
  • z_k ≥ z^*_t,k − M_k (1 − y_t,k)
  • ∑_k y_t,k ≥ 1
  • + a cut of the form : ∑_{i ∈ P(t)} x_i − ∑_{i ∈ N(t)} x_i ≤ |P(t)| − 1
  This is somewhat more natural and intuitive.
• This would also relax the assumption in the paper that all cost-coefficients are integer valued (c_i ∈ Z)
• This approach seemed to work ok for some larger test data with more objectives.
• The M_k's can be found automatically by running some models so a bounding box is found.
• Although this method uses weights, they are not very important, as long as they are nonzero. I.e. a weight of all ones is perfectly appropriate for this algorithm.
• A modeling language that allows loops and multiple solves is in an advantage here. A system that allows only to model one instance of a model, will have troubles to implement these “mini-algorithms” in a clean way.
• A subsequent paper by the authors describes a way to generate not all, but well-dispersed non-dominated solutions. For larger problems where all non-dominated solutions are not possible to find quickly, this may be good approach.