$title eps-Constraint Method for Multiobjective Optimization (EPSCM,SEQ=319) $onText The eps-Constraint Method This is a GAMS implementation of the augmented eps-constraint method for generating the efficient (Pareto optimal, nondominated) solutions in multiobjective problems. The eps-constraint method optimizes one of the objective functions using the remaining objective functions as constraints, varying their right hand side. The generated optimal solutions proved to be efficient solutions of the multiobjective problem under certain conditions. The eps-constraint method consists of two phases: 1. Creation of the payoff table 2. Use the ranges from the payoff table in order to apply the method The augmented eps-constraint uses lexicographic optimization in the construction of the payoff table (in order to secure the Pareto optimality of the individual optima) and a slightly modified objective function in order to ensure the production of Pareto optimal (and not weakly Pareto optimal) solutions. In addition, it performs early exit from infeasible loops improving the performance of the algorithm in multi-objective problems with several objective functions. The algorithm can work also with MIP models. Actually the advantages of the eps-constraint method over the weighting method are more apparent for MIP problems where the non supported Pareto optimal solutions can be produced. A simplified power generation problem is used to illustrate the method: Four types of power generation units are considered in a region, namely, lignite fired, oil fired, natural gas fired and units exploiting renewable energy sources (RES) which are mostly small hydro and wind. The power generation characteristics of these units are shown in table pdata. The yearly demand is 64000 GWh and is characterized by a load duration curve which can be divided into three type of loads: base load (60%), medium load (30%) and peak load (10%). The lignite fired units can be used only for base and middle load, the oil fired units for middle and peak load, the RES units for base and peak load and the natural gas fired units for all type of loads. The endogenous sources are lignite and RES. We consider three objective functions: the minimization of production cost, the minimization of CO2 emissions and the minimization of external dependence (i.e. oil and natural gas) by maximizing the endogenous energy sources. The task is to generate the payoff table and the Pareto optimal (efficient, non-dominated) solutions of the problem. Additional information can be found at: /modlib/adddocs/epscm.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. Keywords: linear programming, eps-constraint method, multiobjective optimization $offText $log --- Using Python library %sysEnv.GMSPYTHONLIB% $inlineCom [ ] $eolCom // $sTitle Example Model Definitions Set p 'power generation units' / Lignite, Oil, Gas, RES / i 'load areas' / base, middle, peak / pi(p,i) 'availability of unit for load types' / Lignite.(base,middle) Oil.(middle,peak), Gas.set.i RES.(base, peak) / es(p) 'endogenous sources' / Lignite, RES / k 'objective functions' / cost, CO2emission, endogenous /; $set min -1 $set max +1 Parameter dir(k) 'direction of the objective functions' / cost %min%, CO2emission %min%, endogenous %max% /; Set pheader / capacity, cost, CO2emission /; Table pdata(pheader,p) Lignite Oil Gas RES capacity [GWh] 31000 15000 22000 10000 cost [$/MWh] 30 75 60 90 CO2emission [t/MWh] 1.44 0.72 0.45 0; Parameter ad 'annual demand in GWh' / 64000 / df(i) 'demand fraction for load type' / base 0.6, middle 0.3, peak 0.1 / demand(i) 'demand for load type in GWh'; demand(i) = ad*df(i); Variable z(k) 'objective function variables'; Positive Variable x(p,i) 'production level of unit in load area in GWh'; Equation objcost 'objective for minimizing cost in K$' objco2 'objective for minimizing CO2 emissions in Kt' objes 'objective for maximizing endogenous sources in GWh' defcap(p) 'capacity constraint' defdem(i) 'demand satisfaction'; * Objective functions objcost.. sum(pi(p,i), pdata('cost',p)*x(pi)) =e= z('cost'); objco2.. sum(pi(p,i), pdata('CO2emission',p)*x(pi)) =e= z('CO2emission'); objes.. sum(pi(es,i), x(pi)) =e= z('endogenous'); defcap(p).. sum(pi(p,i), x(pi)) =l= pdata('capacity',p); defdem(i).. sum(pi(p,i), x(pi)) =g= demand(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'; k1(k)$(ord(k) = 1) = yes; km1(k) = yes; km1(k1) = no; Set kk(k) 'active objective function in constraint allobj'; 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'; Variable a_objval 'auxiliary variable for the objective function' obj 'auxiliary variable during the construction of the payoff table'; Positive Variable sl(k) 'slack or surplus variables for the eps-constraints'; 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.. sum(k1,dir(k1)*z(k1)) + 1e-3*sum(km1,sl(km1)/(maxobj(km1) - minobj(km1))) =e= a_objval; allobj.. sum(kk, dir(kk)*z(kk)) =e= obj; Model mod_payoff / example, allobj / mod_epsmethod / example, con_obj, augm_obj /; option limRow = 0, limCol = 0, solPrint = off, solveLink = %solveLink.callModule%; Parameter payoff(k,k) 'payoff tables entries'; Alias (k,kp); * Generate payoff table applying lexicographic optimization loop(kp, kk(kp) = yes; repeat solve mod_payoff using lp 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.feasibleSolution%, abort 'no feasible solution for mod_payoff'); display payoff; minobj(k) = smin(kp,payoff(kp,k)); maxobj(k) = smax(kp,payoff(kp,k)); $if not set gridpoints $set gridpoints 10 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 'counter' lastZero 'counter' numk(k) 'ordinal value of k starting with 1' numg(g) 'ordinal value of g starting with 0'; 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)); gridrhs(grid(km1,g))$(%min% = dir(km1)) = maxobj(km1) - numg(g)/maxg(km1)*(maxobj(km1) - minobj(km1)); gridrhs(grid(km1,g))$(%max% = dir(km1)) = minobj(km1) + numg(g)/maxg(km1)*(maxobj(km1) - minobj(km1)); display gridrhs; * Walk the grid points and take shortcuts if the model becomes infeasible posg(km1) = 0; $ifThen set SAVESOL * Poor man's strings in GAMS file fSTRING; fSTRING.nw=0; $macro STRINGDEF(sym) singleton set sym / sym / $macro STRING(sym) sym.te(sym) $macro STRINGASSIGN(sym,text) put_utility fSTRING 'assignText' / 'sym' / text $macro STRINGAPPEND(sym,text) put_utility fSTRING 'assignText' / 'sym' / sym.te(sym) text STRINGDEF(fname); $endIf parameter zl(k); embeddedCode Python: sol = set() pauseEmbeddedCode repeat rhs(km1) = sum(grid(km1,g)$(numg(g) = posg(km1)), gridrhs(km1,g)); solve mod_epsmethod maximizing a_objval using lp; if(mod_epsmethod.modelStat <> %modelStat.optimal%, // not optimal is in this case 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 zl(k) = z.l(k); continueembeddedCode: sol.add(tuple(gams.get('zl'))) pauseEmbeddedCode $ ifThen set SAVESOL STRINGASSIGN(fname,'sol'); loop(k, STRINGAPPEND(fname,'_' z.l(k))); STRINGAPPEND(fname,'.gdx'); put_utility 'gdxout' / STRING(fname); execute_unload x.l; $ endIf ); * Proceed forward in the grid firstOffMax = 0; loop(km1$(posg(km1) < maxg(km1) and firstOffMax = 0), posg(km1) = posg(km1) + 1; firstOffMax = numk(km1); ); posg(km1)$(numk(km1) < firstOffMax) = 0; until sum(km1$(posg(km1) = maxg(km1)),1) = card(km1) and firstOffMax = 0; Set s 'solutions' / s1*s50 /; Parameter solutions(s,k) 'unique solutions'; continueembeddedCode: solutions = [] for r in zip(gams.get('s'),sorted(sol)): for rp in r[1]: solutions.append((r[0],*rp)) gams.set('solutions', solutions) endEmbeddedCode solutions display solutions; $ifThen set SAVESOL set ss(s) 'active solutions'; option ss