Description
The Augmented Epsilon Constraint Method version 2 (AUGMECON2) The method is applied to a Multi-Objective Integer Programming problem (specifically a Multi-Objective Multi-Dimensional Knapsack Problem) with 50 binary variables X, 2 objective functions and 2 constraints The AUGMECON2 can be used to generate the exact Pareto set (all the Pareto optimal solutions) if the step size (i.e.the interval between the grid points of the objective functions that are used as constraints) is appropriately chosen. For problems with integer objective function coeffcients the step size should be at most equal to unity. The exact Pareto set of the specific problem consists of 35 Pareto optimal solutions. The solution time is approximately 4 secs on a 2012 machine. The gridpoints are set to 491 = the second objective function range. The output file 2kp50_augmecon2_results.txt contains the payoff table, the gridpoints and the Pareto optimal solutions. The indication "jump" is used to flag when one or more grid points are skipped. The model is separable. The first part of the model (till line 129) is the problem description and the second part (from line 130) is the implementation of the method Additional information can be found at: http://www.gams.com/modlib/adddocs/epscmmip.pdf
Small Model of Type : MIP
Category : GAMS Model library
Main file : epscmmip.gms
$title Improved Version of eps-Constraint Method for Multiobjective Optimization (EPSCMMIP,SEQ=384)
$onText
The Augmented Epsilon Constraint Method version 2 (AUGMECON2)
The method is applied to a Multi-Objective Integer Programming problem
(specifically a Multi-Objective Multi-Dimensional Knapsack Problem)
with 50 binary variables X, 2 objective functions and 2 constraints
The AUGMECON2 can be used to generate the exact Pareto set (all the Pareto
optimal solutions) if the step size (i.e.the interval between the grid points
of the objective functions that are used as constraints) is appropriately
chosen. For problems with integer objective function coeffcients the step size
should be at most equal to unity.
The exact Pareto set of the specific problem consists of 35 Pareto optimal
solutions. The solution time is approximately 4 secs on a 2012 machine. The
gridpoints are set to 491 = the second objective function range. The output file
2kp50_augmecon2_results.txt contains the payoff table, the gridpoints and the
Pareto optimal solutions. The indication "jump" is used to flag when one or
more grid points are skipped.
The model is separable. The first part of the model (till line 129) is the
problem description and the second part (from line 130) is the implementation
of the method
Additional information can be found at:
http://www.gams.com/modlib/adddocs/epscmmip.pdf
Mavrotas, G, Effective implementation of the eps-constraint method in
Multi-Objective Mathematical Programming problems.
Applied Mathematics and Computation 213, 2 (2009), 455-465.
Mavrotas, G, and Florios, K, An improved version of the augmented
eps-constraint method (AUGMECON2) for finding the exact Pareto set in
Multi-Objective Integer Programming problems.
Applied Mathematics and Computation 219, 18 (2013), 9652-9669
Keywords: mixed integer linear programming, multiobjective optimization,
eps-constraint method, exact Pareto set, mathematics
$offText
$eolCom //
$sTitle Example Model Definitions
Set
I 'constraints' / i1* i2 /
J 'decision variables' / j1*j50 /
K 'objective functions' / k1* k2 /;
Parameter
dir(k) 'direction of the objective functions 1 for max and -1 for min' / k1 1, k2 1 /
b(I) 'RHS of the constraints' / i1 1445, i2 1502.5 /;
Table c(K,J) 'matrix of objective function coefficients C'
j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16 j17
k1 21 69 26 92 77 30 96 80 60 61 52 92 19 10 63 34 100
k2 24 92 53 25 10 31 83 34 64 69 95 40 59 87 13 94 53
+ j18 j19 j20 j21 j22 j23 j24 j25 j26 j27 j28 j29 j30 j31 j32 j33 j34
k1 60 11 12 37 100 74 17 60 69 49 69 49 59 17 21 74 85
k2 52 61 53 78 34 89 32 28 56 52 40 41 59 35 96 72 55
+ j35 j36 j37 j38 j39 j40 j41 j42 j43 j44 j45 j46 j47 j48 j49 j50
k1 83 41 29 63 56 38 66 92 25 84 89 21 46 94 96 92
k2 100 44 90 66 59 22 72 25 36 16 56 91 61 56 66 53 ;
Table a(I,J) 'matrix of technological coefficients A'
j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16 j17
i1 84 49 68 20 97 74 60 30 13 95 19 41 17 95 73 12 66
i2 19 96 93 64 72 91 32 96 44 76 69 82 51 38 52 22 83
+ j18 j19 j20 j21 j22 j23 j24 j25 j26 j27 j28 j29 j30 j31 j32 j33 j34
i1 55 75 20 56 80 59 66 25 70 95 96 62 74 31 59 21 85
i2 27 70 56 29 89 86 48 13 95 66 94 16 44 67 90 48 29
+ j35 j36 j37 j38 j39 j40 j41 j42 j43 j44 j45 j46 j47 j48 j49 j50
i1 45 97 23 53 51 95 58 68 62 45 83 82 47 15 52 72
i2 90 54 77 28 100 86 51 62 40 54 21 55 50 62 51 77 ;
Variable
Z(K) 'objective function variables'
X(J) 'decision variables';
Binary Variable X;
Equation
objfun(K) 'objective functions'
con(I) 'constraints';
objfun(K).. sum(J, c(K,J)*X(J)) =e= Z(K);
con(I).. sum(J, a(I,J)*X(J)) =l= b(I);
Model example / all /;
$sTitle eps-constraint Method
Set
k1(k) 'the first element of k'
km1(k) 'all but the first elements of k'
kk(k) 'active objective function in constraint allobj';
k1(k)$(ord(k) = 1) = yes;
km1(k) = yes;
km1(k1) = no;
Parameter
rhs(k) 'right hand side of the constrained obj functions in eps-constraint'
maxobj(k) 'maximum value from the payoff table'
minobj(k) 'minimum value from the payoff table'
numk(k) 'ordinal value of k starting with 1';
Scalar
iter 'total number of iterations'
infeas 'total number of infeasibilities'
elapsed_time 'elapsed time for payoff and e-sonstraint'
start 'start time'
finish 'finish time';
Variable
a_objval 'auxiliary variable for the objective function'
obj 'auxiliary variable during the construction of the payoff table'
sl(k) 'slack or surplus variables for the eps-constraints';
Positive Variable sl;
Equation
con_obj(k) 'constrained objective functions'
augm_obj 'augmented objective function to avoid weakly efficient solutions'
allobj 'all the objective functions in one expression';
con_obj(km1).. z(km1) - dir(km1)*sl(km1) =e= rhs(km1);
* We optimize the first objective function and put the others as constraints
* the second term is for avoiding weakly efficient points
augm_obj..
a_objval =e= sum(k1,dir(k1)*z(k1))
+ 1e-3*sum(km1,power(10,-(numk(km1) - 1))*sl(km1)/(maxobj(km1) - minobj(km1)));
allobj.. sum(kk, dir(kk)*z(kk)) =e= obj;
Model
mod_payoff / example, allobj /
mod_epsmethod / example, con_obj, augm_obj /;
Parameter payoff(k,k) 'payoff tables entries';
Alias (k,kp);
option optCr = 0, limRow = 0, limCol = 0, solPrint = off, solveLink = %solveLink.loadLibrary%;
* Generate payoff table applying lexicographic optimization
loop(kp,
kk(kp) = yes;
repeat
solve mod_payoff using mip maximizing obj;
payoff(kp,kk) = z.l(kk);
z.fx(kk) = z.l(kk); // freeze the value of the last objective optimized
kk(k++1) = kk(k); // cycle through the objective functions
until kk(kp);
kk(kp) = no;
* release the fixed values of the objective functions for the new iteration
z.up(k) = inf;
z.lo(k) = -inf;
);
if(mod_payoff.modelStat <> %modelStat.optimal% and
mod_payoff.modelStat <> %modelStat.integerSolution%,
abort 'no optimal solution for mod_payoff';);
File fx / 2kp50_augmecon2_results.txt /;
put fx ' PAYOFF TABLE'/;
loop(kp,
loop(k, put payoff(kp,k):12:2;);
put /;
);
minobj(k) = smin(kp,payoff(kp,k));
maxobj(k) = smax(kp,payoff(kp,k));
* gridpoints are calculated as the range (difference between max and min) of
* the 2nd objective function from the payoff table
$if not set gridpoints $set gridpoints 491
Set
g 'grid points' / g0*g%gridpoints% /
grid(k,g) 'grid';
Parameter
gridrhs(k,g) 'RHS of eps-constraint at grid point'
maxg(k) 'maximum point in grid for objective'
posg(k) 'grid position of objective'
firstOffMax 'some counters'
lastZero 'some counters'
* numk(k) 'ordinal value of k starting with 1'
numg(g) 'ordinal value of g starting with 0'
step(k) 'step of grid points in objective functions'
jump(k) 'jumps in the grid points traversing';
lastZero = 1;
loop(km1,
numk(km1) = lastZero;
lastZero = lastZero + 1;
);
numg(g) = ord(g) - 1;
grid(km1,g) = yes; // Here we could define different grid intervals for different objectives
maxg(km1) = smax(grid(km1,g), numg(g));
step(km1) = (maxobj(km1) - minobj(km1))/maxg(km1);
gridrhs(grid(km1,g))$(dir(km1) = -1) = maxobj(km1) - numg(g)/maxg(km1)*(maxobj(km1) - minobj(km1));
gridrhs(grid(km1,g))$(dir(km1) = 1) = minobj(km1) + numg(g)/maxg(km1)*(maxobj(km1) - minobj(km1));
put / ' Grid points' /;
loop(g,
loop(km1, put gridrhs(km1,g):12:2;);
put /;
);
put / 'Efficient solutions' /;
* Walk the grid points and take shortcuts if the model becomes infeasible or
* if the calculated slack variables are greater than the step size
posg(km1) = 0;
iter = 0;
infeas = 0;
start = jnow;
repeat
rhs(km1) = sum(grid(km1,g)$(numg(g) = posg(km1)), gridrhs(km1,g));
solve mod_epsmethod maximizing a_objval using mip;
iter = iter + 1;
if(mod_epsmethod.modelStat<>%modelStat.optimal% and
mod_epsmethod.modelStat<>%modelStat.integerSolution%,
infeas = infeas + 1; // not optimal is in this case infeasible
put iter:5:0, ' infeasible' /;
lastZero = 0;
loop(km1$(posg(km1) > 0 and lastZero = 0), lastZero = numk(km1));
posg(km1)$(numk(km1) <= lastZero) = maxg(km1); // skip all solves for more demanding values of rhs(km1)
else
put iter:5:0;
loop(k, put z.l(k):12:2;);
jump(km1) = 1;
* find the first off max (obj function that hasn't reached the final grid point).
* If this obj.fun is k then assign jump for the 1..k-th objective functions
* The jump is calculated for the innermost objective function (km=1)
jump(km1)$(numk(km1) = 1) = 1 + floor(round(sl.L(km1),6)/step(km1));
loop(km1$(jump(km1) > 1), put ' jump';);
put /;
);
* Proceed forward in the grid
firstOffMax = 0;
loop(km1$(posg(km1) < maxg(km1) and firstOffMax = 0),
posg(km1) = min((posg(km1) + jump(km1)),maxg(km1));
firstOffMax = numk(km1);
);
posg(km1)$(numk(km1) < firstOffMax) = 0;
abort$(iter > 1000) 'more than 1000 iterations, something seems to go wrong';
until sum(km1$(abs(posg(km1) - maxg(km1))<1e-6),1) = card(km1) and firstOffMax = 0;
finish = jnow;
elapsed_time = (finish - start)*60*60*24;
put /;
put 'Infeasibilities = ', infeas:5:0 /;
put 'Elapsed time: ',elapsed_time:10:2, ' seconds' /;