Multi-Objective Optimization (moo)

Table of Contents

Introduction

In Multi-Objective Optimization (MOO) there is more than one objective function and there is no single optimal solution that simultaneously optimizes all the objective functions. In MOO the concept of optimality is replaced by Pareto efficiency or optimality. Pareto efficient (or optimal, nondominated, etc.) solutions are solutions that cannot be improved in one objective without deteriorating in at least one of the other objectives. The mathematical definition of an efficient solution is the following:

Without loss of generality assume that all objective functions \(f_k\) with \(k=1,...,p\) are for minimization. A feasible solution \(x\) of a MOO problem is (strongly) efficient if there is no feasible solution \(x'\) such as \(f_k(x') \leq; f_k(x)\) for \(k=1,...,p\) with at least one strict inequality. A feasible solution \(x\) is weakly efficient if there is no feasible solution \(x'\) such as \(f_k(x') < f_k(x)\) for \(k=1,...,p\).

In case of a non-convex feasible objective and decision space (e.g. MIP), the set of efficient solutions can be further partitioned into supported and non-supported efficient solutions. Supported efficient solutions are optimal solutions of the weighted sum single objective problem for any non-negative weights (corner solutions of the feasible space). See the following graph that illustrates the difference between supported and non-supported efficient solutions:

Supported and non-supported efficient solutions in the objective space Z. Graph taken from Li, Chendong & Park, Chulwoo & Pattipati, Krishna & Kleinman, David. (2011). Distributed algorithms for biobjective assignment problems. Proceedings of the IEEE Conference on Decision and Control. 5893-5898. 10.1109/CDC.2011.6161434.

A simple MOO method is lexicographic optimization. Lexicographic optimization presumes that the decision-maker has lexicographic preferences, ranking the efficient solutions to a lexicographic order of the objective functions. A lexicographic minimization problem can be written as:

\(lex\) \(min\) \(f_1(x),f_2(x),...,f_p(x)\)
s.t.
\(x \in S\)

where \(x\) is the vector of decision variables, \(S\) is the feasible region, and \(f_1(x),f_2(x),...,f_p(x)\) are the \(p\) objective functions ordered from the most to the least important. A lexicographic minimization problem with \(p\) objective functions can be solved using a sequence of \(p\) single-objective optimization problems as follows:

For \(i=1,...,p\) do

  • Solve single-objective optimization problem:

    \(min\) \(z_i\)
    s.t.
    \(f_k(x)=z_k\) \(\forall k=1,...,i-1\)
    \(x \in S\)

    where \(z_k\) are the objective function variables

  • Fix the value of objective function variable \(z_i\) for the following iterations

End for

The MOO methods provided by this LibInclude file use lexicographic optimization to calculate a payoff table to obtain e.g. ranges for the objectives. The \(p \times p\) payoff table contains the results of \(p\) lexicographic optimizations with orders \( \lbrace f_1(x),f_2(x),...,f_p(x) \rbrace , \lbrace f_2(x),f_3(x),...,f_p(x),f_1(x) \rbrace , ..., \lbrace f_p(x),f_1(x),...,f_{p-1}(x) \rbrace \).

This LibInclude file provides a collection of methods to generate Pareto efficient solutions for MOO models. All implemented methods are so-called a-posteriori (or generation) methods for generating Pareto efficient solutions. After the Pareto efficient set or a representative subset of the Pareto efficient set has been generated, it is returned to the user, who then analyses the solutions and selects the most preferred among the alternatives.

Usage

LibInclude file moo.gms provides methods for multi-objective optimization in GAMS and can be used as follows:

$libinclude moo method modelName modelType objectiveSet direction objectiveVar points paretoPoints [optList]

Where:

Argument Description
method Specify the name of the method to use. Available methods are EpsConstraint, RandomWeighting, and Sandwiching. See Methods for details.
modelName Specify the name of the multi-objective optimization model that should be solved.
modelType Specify the corresponding model type, e.g. LP.
objectiveSet Specify the set of objectives.
direction One-dimensional parameter that contains the optimization direction for each objectiveSet element - min (-1) or max (+1).
objectiveVar Specify the objective function variable. Compared to the traditional scalar objective function variable with a single objective optimization problem in GAMS, this objective function variable is one-dimensional.
points Specify the set of pareto points. The number of members in the set defines the maximum number of returned pareto points. Therefore, make sure to define the set large enough.
paretoPoints Specify the two-dimensional (points,objectiveSet) parameter that will be used to save the objective values of pareto points.
optList Optional list of -optkey=optval option pairs, e.g. -savepoint=1 -optfile_main=99. See available options below.
Method Option Description
All: debug Set to 1 to activate debugging. This will create additional information being written to the log as well as creating files to support debugging, e.g. if savepoint is activated a file savepoint_iteration_map.csv is written to the savepoint_dir that contains the mapping of savepoint files to the corresponding iteration. Default: 0.
optfile_init Solver option file for initial solves. In particular, solves to calculate the payoff table. Default: 0.
optfile_main Solver option file for main solves, i.e. the main part of the methods. Default: 0.
plot Set to 1 to activate plots for multi objective problems with two or three objective functions, e.g. a plot of the pareto points. Python package matplotlib needs to be installed. Default: 0.
plot_dir Set the directory where plots should be saved. Per default the GAMS working directory will be used.
pre Prefix for names of GAMS symbols used in moo.gms. Allows to make GAMS symbol names unique and avoids problems with redefined symbols. The default is MOO_<MOO___COUNTER>, where MOO___COUNTER is a (global) compile-time variable that counts the MOO runs in the GAMS program, so each run has unique GAMS symbol names. This option allows to change the default.
savepoint If set to 1, a point format GDX file that contains information on the solution point will be saved for each pareto point. Default is 0, where no solution points will be saved.
savepoint_dir Set the directory where save point files should be saved. Per default the GAMS working directory will be used.
savepoint_filename Set the (base) name of the save point files. The file name gets extended by the points set to generate unique file names as follows: <savepoint_filename>_<points>.gdx. Per default this option is set to the name of the method. If an empty string is provided the file name matches exactly with the points set: <points>.gdx.
solver Default solver for all model types that the solver is capable to process.
EpsConstraint: gridpoints Specify the number of grid points per objective function used as a constraint (default: 10). As an example, consider a model with three objective functions, the Epsilon Constraint method will optimize one objective function using the other two as contraints. With 10 grid points (default) the maximum number of runs will therefore be 10x10. In general, increasing the number of grid points increases the density of the approximation.
RandomWeighting: iterations Terminating condition in terms of iterations. Default: inf.
seed Set the seed for the pseudo random number generator that generates the weights. Set to None to set the seed to a random number. Default: 3141.
time Terminating condition in terms of time in seconds. Default: 3600.
Sandwiching: gap Terminating condition in terms of a gap (default: 0.01). The gap is calculated by the current maximum distance between the inner and outer approximation divided by the initial maximum distance.
time Terminating condition in terms of time in seconds. Default: 3600.


Methods

This section gives an introduction and details on the implemented methods for MOO. The implemented methods are so-called a-posteriori (or generation) methods for generating Pareto efficient solutions. After the Pareto efficient set or a representative subset of the Pareto efficient set has been generated, it is returned to the user, who then analyses the solutions and selects the most preferred among the alternatives.

Augmented Epsilon Constraint [EpsConstraint]

The Epsilon Constraint method optimizes one of the objective functions using the other objective functions as constraints. By parametrical variation in the Right-Hand Side (RHS) of the constrained objective functions the efficient solutions of the problem are obtained. The implemented Augmented Epsilon Constraint method is based on:

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

The Augmented Epsilon Constraint method is a variation of the Epsilon Constraint method that produces only efficient solutions (no weakly efficient solutions) and also avoids redundant iterations as it can perform early exit from the relevant loops (that lead to infeasible solutions) and exploits information from slack variables in every iteration.

Compared to the Random Weighting method the runs are exploited more efficiently and almost every run produces a different efficient solution, thus obtaining a better representation of the efficient set. Allowed model types: LP, MIP as well as corresponding Non Linear cases. The method is able to find the exact Pareto set for integer programming problems by appropriately tuning the number of grid points. The Augmented Epsilon Constraint method is able to produce supported and unsupported efficient solutions.

Implementation details:

  • The Augmented Epsilon Constraint method optimizes one of the objective functions using the other objective functions as constraints, incorporating the constraint part of the model as follows:

    \( max (f_1(x) + eps \times (s_2/r_2+ 10^{-1} \times s_3/r_3+...+10^{-(p-2)} \times s_p/r_p)) \)

    s.t.
    \(f_2(x) - s_2 = e_2\)
    \(f_3(x) - s_3 = e_3\)
    ...
    \(f_p(x) - s_p = e_p\)
    \(x \in S\) and \(s_k \in \mathbb{R}^+\)

    where \(x\) is the vector of decision variables, \(S\) is the feasible region, and \(f_1(x),f_2(x),...,f_p(x)\) are the \(p\) objective functions. \(e_k\) are the parameters for the RHS for the specific iteration drawn from the grid points of the objective functions \(k\) ( \(k=2,3,...,p\)). \(s_k\) are slack or surplus variables, \(r_k\) is the range of the \(k\)-th objective function (calculated from the payoff table), and \(eps \in [10^{-6}, 10^{-3}]\).

  • The method is implemented as follows: From the payoff table we obtain the range of the \(p-1\) objective functions that are used as constraints. The range of the \(k\)-th objective function is divided into \(q_k\) equal intervals using \((q_k-1)\) intermediate equidistant grid points. Thus, we have in total \((q_k+1)\) grid points that are used to vary parametrically the RHS ( \(e_k\)) of the \(k\)-th objective. The total number of runs becomes \((q_2+1)\times (q_3+1)\times ...\times (q_p+1)\).

    The discretization step for objective function \(k\):

    \(step_k=r_k/q_k\)

    The RHS of the corresponding constraint in the \(i\)-th iteration in the specific objective function will be:

    \(e_{ki} = fmin_k + i_k\times step_k\)

    where \(fmin_k\) is the minimum obtained from the payoff table and \(i_k\) the counter for the specific objective function.

  • At each iteration we check the surplus variable that corresponds to the innermost objective function (in this case \(p=2\)) and calculate the bypass coefficient as:

    \(b = floor(s_2/step_2)\)

    Where \(floor()\) is a function that always rounds down and returns the largest integer less than or equal to a given number. When the surplus variable \(s_2\) is larger than the \(step_2\) it is implied that in the next iteration the same solution will be obtained with the surplus variable being \(s_2-step_2\). Therefore, the bypass coefficient actually indicates how many consecutive iterations can be bypassed.

  • Early exit from the loop when the problem becomes infeasible for some combination of \(e_k\). The bounding strategy for each one of the objective functions starts from the more relaxed formulation (lower bound for maximization objective function or upper bound for a minimization) and move to the most strict. If we arrive to an infeasible solution there is no need to perform the remaining runs of the loop, as the problem will become even stricter.

Random Weighting [RandomWeighting]

With the Random Weighting method, the weighted sum of the objective functions is optimized. By randomly varying the weights different pareto points are obtained.

Allowed model types: All model types with an objective variable, in particular LP, MIP, NLP, and MINLP. Note that the Random Weighting method may spend a lot of runs that are redundant in the sense that there can be a lot of combinations of weights that result in the same efficient point. Also note that the Random Weighting method cannot produce unsupported efficient solutions.

Implementation details:

  • For Random Weighting the scaling of the objective functions has a strong influence on the obtained results. Therefore, the objective functions are scaled between 0 and 1 using the objective function ranges obtained from the payoff table.

Sandwiching

Sandwiching algorithms are used to approximate the pareto front by sandwiching the nondominated set between an inner and outer approximation. Sandwiching algorithms iteratively improve both the Inner and Outer approximation of the Pareto Set (IPS and OPS) to minimize the distance between the approximations. The implemented Sandwiching method is based on:

Koenen, M., Balvert, M., and Fleuren, HA. (2023). A Renewed Take on Weighted Sum in Sandwich Algorithms: Modification of the Criterion Space. (Center Discussion Paper; Vol. 2023-012). CentER, Center for Economic Research.

An advantage of the Sandwiching method is that the distance between the IPS and OPS provides an upper bound on the approximation error which can also be used as a termination criterion for the algorithm. Please be aware that for more than two objectives the distance using the definition of the Koenen et al. (2023) can sometimes increase over the iterations. Allowed model type: LP.

Implementation details:

  • The Sandwiching algorithm of Koenen et al. (2023) is based on a weighted sum approach, i.e. the weighted sum of the objective functions is optimized. By varying the weights at each iteration the efficient points are obtained. The algorithm selects weights derived from the IPS that have the potential to reduce the approximation error the most.
  • The IPS is a convex hull encompassed in the objective space formed by a set of efficient points and the OPS is formed by supporting halfspaces of those efficient points.
  • The basic procedure of the algorithm can be described as follows:
    1. Initialize the IPS and OPS using the pseudo nadir point and so called anchor points. The pseudo nadir point and anchor points are obtained from the payoff table and also define the initial bounding box of the objective space.
    2. Expand the IPS based on its distance to the OPS. For each relevant hyperplane of the convex hull of the IPS determine the distance to the OPS and find the hyperplane for which the current distance from the IPS to the OPS is maximal. The normal of this hyperplane will be used as weights in the current iteration. If a new efficient point is found, update the IPS and OPS, otherwise set the distance of the current hyperplane to zero. The procedure is repeated until a stopping criterion is met (see options gap or time).
  • The calculation of the distance between IPS and OPS is based on SUB (see Koenen et al. (2023)).
  • The scaling of the objective functions has a strong influence on the obtained results. Therefore, the objective functions are scaled between 0 and 1 using the objective function ranges obtained from the payoff table.

Examples

Example 1: Solve scalable multi-objective knapsack model

The following example demonstrates how to solve a multi-objective knapsack model (MIP) using the moo LibInclude:

$if not set NBITM $set NBITM 50
$if not set NBDIM $set NBDIM  2
$if not set NBOBJ $set NBOBJ  2

$if not set GRIDPOINTS $set GRIDPOINTS 20

Set
   i 'items'                / i1*i%NBITM% /
   j 'weight dimensions'    / j1*j%NBDIM% /   
   k 'value dimensions'     / k1*k%NBOBJ% /
   p 'pareto points'        / point1*point1000 /;

Parameter
   a(i,j)           'weights of item i'
   c(i,k)           'values of item i'
   b(j)             'knapsack capacity for weight j'
   dir(k)           'direction of the objective functions 1 for max and -1 for min' / #k 1/   
   pareto_obj(p,k)  'objective values of the pareto points' 
;
a(i,j) = UniformInt(1,100);
c(i,k) = UniformInt(1,100);
b(j)   = UniformInt(1,100) * %NBITM%/4; 

Variable
   Z(k) 'objective variables'
   X(i) 'decision variables';

Binary Variable X;

Equation
   objfun(k) 'objective functions'
   con(j)    'capacity constraints';

objfun(k).. sum(i, c(i,k)*X(i)) =e= Z(k);

con(j)..    sum(i, a(i,j)*X(i)) =l= b(j);

Model example / all /;

$onEcho > cplex.opt
threads 1
$offEcho

$libInclude moo EpsConstraint example MIP k dir Z p pareto_obj -gridpoints=%GRIDPOINTS% -solver=cplex -optfile_init=1 -optfile_main=1

display pareto_obj;

The example is easily scalable by allowing to set the number of items, the number of weight dimensions and the number of objectives through double dash parameters. This example is also part of the GAMS Data Utilities Library, see model [moo01] for reference. The model is solved using EpsConstraint method with 20 gridpoints. The objective function values of the pareto points found are saved in parameter pareto_obj:

----    618 PARAMETER pareto_obj  objective values of the pareto points

                 k1          k2

point1      953.000    1076.000
point2      956.000    1062.000
point3      972.000    1052.000
point4     1005.000    1037.000
point5     1030.000    1024.000
point6     1049.000    1004.000
point7     1063.000     990.000
point8     1088.000     955.000
point9     1091.000     946.000
point10    1105.000     932.000
point11    1121.000     876.000
point12    1130.000     827.000
point13    1137.000     783.000

If plot is set to 1 and Python package matplotlib is installed, a plot of the pareto points is generated:

Pareto points for knapsack problem with two objectives

Setting --NBOBJ=3 with model [moo01] allows to solve a knapsack problem with three objectives. These are the resulting pareto points:

Pareto points for knapsack problem with three objectives

Example 2: Solve multi-objective power generation model

The following example demonstrates how to solve a multi-objective power generation model (LP) using the moo LibInclude:

$inlineCom [ ]

$if not set NBOBJ $set NBOBJ  2
$if not set METHOD $set METHOD Sandwiching

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 /
   points  'pareto points' / point1*point1000 /;

$set min -1
$set max +1

Parameter dir(k) 'direction of the objective functions 1 for max and -1 for min'
                 / cost %min%, CO2emission %min%, endogenous %max% /
          pareto_obj(points,k) 'objective values of the pareto points';

Set pheader / capacity, cost, CO2emission /;

Table pdata(pheader,p)
                         Lignite    Oil    Gas    RES
   capacity [GWh]          61000  25000  42000  20000
   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';

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

Set kk(k) 'active objective functions';
kk(k) = yes;
$if %NBOBJ%==2 kk('endogenous') = no;

$onEcho > cplex.opt
threads 1
$offEcho

$libInclude moo %METHOD% example LP kk dir z points pareto_obj -iterations=20 -gridpoints=5 -savepoint=1 -savepoint_filename= -savepoint_dir=savepoints -solver=cplex -optfile_init=1 -optfile_main=1
execute 'gdxmerge savepoints%system.DirSep%*.gdx > %system.NullFile%';

display pareto_obj;

This example is also part of the GAMS Data Utilities Library, see model [moo02] for reference. The model is solved using the Sandwiching method with a gap of 0.01 (default). The objective function values of the pareto points found are saved in parameter pareto_obj:

----   1164 PARAMETER pareto_obj  objective values of the pareto points

              cost  CO2emissi~

point1 2112000.000   85824.000
point2 3180000.000   50580.000
point3 4470000.000   20340.000

In addition, savepoint is activated and thus, for each pareto point a point format GDX file that contains information on the solution point will be written. Using gdxmerge all savepoint GDX files can be merged into a single GDX file.

If plot is set to 1 and Python package matplotlib is installed, a plot of the pareto front is generated:

Pareto front for power generation problem with two objectives

The plot shows the pareto efficient extreme points as well as the lines connecting these points (the IPS). An animation of the algorithm shows how the IPS and OPS are updated at each iteration:

Sandwiching animation for power generation problem with two objectives

Starting with the initial points from the payoff table, the IPS and the OPS are constructed (iteration 0). In iteration 1, a new point is found and the IPS and OPS are updated. In the last two iterations no new point is found but the OPS is updated so the upper bound on the approximation error (distance between the IPS and OPS) decreases and the gap of 0.01 is reached.

Setting --NBOBJ=3 with model [moo02] allows to solve a power generation problem with three objectives. This is the resulting pareto front:

Pareto front for power generation problem with three objectives

Example 3: Solve Bensolvehedron model

The following example demonstrates how to solve an instance of Bensolvehedron (LP) using the moo LibInclude:

* number of objectives
$if not set p $set p 2
* complexity of the polyhedron
$if not set m $set m 3

$eval p2m %p%+2*%m%
$eval n (%p2m%)**%p%
$eval n2 2*%n%

Set p               'objectives' / p1*p%p% /
    m               'complexity' / m1*m%m% /
    n2              'variables 2x' / n1*n%n2% /
    n(n2)           'variables' / n1*n%n% /
    cel             'elements in C' / c1*c%p2m% /
    comb(n,p,cel)   'all combinations of elements in C' / system.powerSetRight /
    points          'pareto points' / point1*point1000 /;
    
Parameter C(p,n)                'objective matrix'
          C_val(cel)            'value of elements in C' 
          A(n2,n)               'constraint matrix'
          b(n2)                 'constraint vector'
          dir(p)                'direction of the objective functions 1 for max and -1 for min' / #p -1 /   
          pareto_obj(points,p)  'objective values of the pareto points';

A(n,n) = 1;
A(n2,n)$(ord(n2)=ord(n)+card(n)) = -1;
b(n) = 1;
C_val(cel) = ord(cel)-((%p2m%-1)/2)-1;
loop(comb(n,p,cel),
  C(p,n) = C_val(cel);
);

Variable Z(p)           'objective variables';
Positive Variable X(n)  'decision variables';

Equations objectives(p)
          constraints(n2);
          
objectives(p).. Z(p) =e= sum(n, C(p,n)*X(n));
constraints(n2).. sum(n,A(n2,n)*X(n)) =l= b(n2);

Model bensolvehedron / all /;

$onEcho > cplex.opt
threads 1
$offEcho

$libinclude moo Sandwiching bensolvehedron LP p dir Z points pareto_obj -solver=cplex -optfile_init=1 -optfile_main=1

display pareto_obj;

Bensolvehedron is a class of multi-objective optimization problems introduced on besolve.org. An instance of Bensolvehedron is generated based on the number of objectives p and the complexity of the polyhedron m.

This example is also part of the GAMS Data Utilities Library, see model [moo03] for reference. The model is solved using the Sandwiching method with a gap of 0.01 (default). The objective function values of the pareto points found are saved in parameter pareto_obj:

----   1266 PARAMETER pareto_obj  objective values of the pareto points

                 p1          p2

point1      -64.000         EPS
point2      -63.000      -7.000
point3      -61.000     -15.000
point4      -58.000     -22.000
point5      -55.000     -27.000
point6      -50.000     -34.000
point7      -34.000     -50.000
point8      -27.000     -55.000
point9      -22.000     -58.000
point10     -15.000     -61.000
point11      -7.000     -63.000
point12         EPS     -64.000

If plot is set to 1 and Python package matplotlib is installed, a plot of the pareto front is generated:

Pareto front for Bensolvehedron problem with two objectives

The plot shows the pareto efficient extreme points as well as the lines connecting these points (the IPS). An animation of the algorithm shows how the IPS and OPS are updated at each iteration:

Sandwiching animation for Bensolvehedron problem with two objectives

Setting --p=3 with model [moo03] allows to solve the Bensolvehedron problem with three objectives. With a gap of 0.005 the Sandwiching method generates the following pareto front:

Pareto front for Bensolvehedron problem with three objectives