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