GAMS code for solving a multi-objective problem using the improved version of the augmented e-constraint method can be found here:
https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_epscmmip.html
In this model, both objectives are used for maximization purpose. How do I convert this model for a bi-objective minimization integer programming problem?
AUGMECON2 code:
$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.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 '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.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' /;