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 x1 − 2x2
max -x1+3x2
x1 − 2x2 ≤ 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 |
(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 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.
