Friday, February 26, 2010

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

In 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 x1 − 2x2
max -x1+3x2
x1 − 2x2 ≤ 0
x1,x2 ∈ {0,1,2}


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

Erwin Kalvelagen,, feb 2010.

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


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

'minimal improvement' /k1 1, k2 1/
'initial weights' /k1 4, k2 3/
'big-M' /k1 4, k2 2/
'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);


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(*,*);
display sols;

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

'dynamic set'

binary variable y(tall,k);
'solution not dominated by earlier solutions'
'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));
sum(k, y(t,k)) =g= 1;

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

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

option limrow=0;
option limcol=0;

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

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


      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
    • zk ≥ z*t,k − Mk (1 − yt,k)
    • k yt,k ≥ 1
    • + a cut of the form : ∑i ∈ P(t) xi − ∑i ∈ N(t) xi ≤ |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 (ci ∈ Z)
  • This approach seemed to work ok for some larger test data with more objectives.
  • The Mk'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.