epscmmip.gms : Improved version of eps-Constraint Method for Multiobjective Optimization

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' /;