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:

<a href="http://www.gams.com/modlib/adddocs/epscmmip.pdf">http://www.gams.com/modlib/adddocs/epscmmip.pdf</a>

References

  • Mavrotas, G, Effective implementation of the ε-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.

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
$offtext

$eolcom //
$STitle Example model definitions

Sets
   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
;

Variables
   Z(K)      objective function variables
   X(J)      decision variables
Binary Variables X;

Equations
   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
Variables
   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 Variables sl
Equations
   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 / ;
Model 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.Integer Solution%,
   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, 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.Integer Solution%,
    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 reach 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(sl.L(km1)/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$(posg(km1)=maxg(km1)),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' / ;