Multiobjective Optimization on a Budget

Indicator-based evolutionary algorithms are among the most powerful multi-objective algorithms, in particular when using hypervolume (HV) contribution as indicator. They are not really suitable for many-objective problems, as the computational cost for computing HV contributions becomes prohibitive as the number of objectives increases. This working group was looking at ways to make HV-indicator based evolutionary algorithms computationally feasible also for many-objective problems. Since steady-state indicator-based EAs in every iteration only generate one new solution and remove one solution from the population, we saw the following opportunities:

1. 1.

   Identifying the solution with the least HV contribution may be easier than computing the HV contribution for all solutions.

2. 2.

   The problems to be solved in each iteration are very similar, with only one solution replaced by another. A lot of the computational results may be transferred from one iteration to the next.

3. 3.

   Evolutionary algorithms are stochastic algorithms and inherently tolerant to noise. Thus it is not clear whether it is actually necessary to always correctly identify the solution with the minimal HV contribution. Perhaps an approximation with tuneable precision would be sufficient and not jeopardise the optimisation behaviour of the MOEA.

4. 4.

   HV calculations, in particular Monte Carlo approximations, could benefit from the use of GPUs.

Hypervolume contribution of a point $s$ defined by equation (4) could alternatively be calculated as the difference of hypervolume of $\{s\}$ and hypercuboid-bounded hypervolume of $S\setminus \{s\}$ within $[r_{\ast },s]$, i.e.:

where $ℒ([r_{\ast },s])$ is the hypervolume of hypercuboid $[r_{\ast },s]$. In practice, the use of equation (5) allows for a faster calculation of hypervolume contribution than equation (4), since hypervolume is calculated just once (the time of calculation of $ℒ([r_{\ast },s])$ is negligible) and many points in $S$ may become dominated after projection onto $[r_{\ast },s]$.

There is a lot of literature on how to efficiently compute the HV, and also some literature on how to compute HV contributions. In fact, it is too much to list here, so the interested reader is instead referred to the surveys in [6, 11].

[8] proposed the Quick Extreme Hypervolume Contribution (QEHC) algorithm which may be used to efficiently find a point (solution) with either the minimum or the maximum contribution for any number of objectives. Since calculation of hypervolume contribution boils down to the calculation of hypervolume (see (5)), QEHC uses the algorithm introduced in [7] to calculate concurrently HV contribution of each solution with quickly converging guaranteed lower and upper bounds for hypervolume contribution. It then uses these bounds to stop calculation of HV contributions for solutions that may not give the unique minimum or maximum contribution.

[3] and [9] proposed greedy lazy approaches for hypervolume subset selection problem with either incremental or decremental approach. These algorithms in each iteration select the solution with the maximum ([3]) or minimum ([9]) hypervolume contribution and utilize the fact that hypervolume is a non-decreasing submodular function:

Since these algorithms only add ([3]) or remove ([9]) a solution in each iteration, they may use contributions calculated in a previous iteration as the upper or lower bounds in a subsequent iteration. Note, however, that these bounds cease to be valid if solutions are both added and removed form a set of solutions.

[4] extends [3] and exploits submodular properties of the HV indicator to reduce the number of HV contribution calculations when selecting a subset from a large number of Pareto-optimal solutions. [10] proposes a local search method for selecting the subset of size $k$ with maximum HV from a larger set of Pareto optimal solutions. Among other things, they show that the number of solutions whose HV contribution is affected by removal of one solution grows very quickly with increasing number of objectives (Fig. 4 in [10]).

[2] propose a fast approximation algorithm to determine the solution with the smallest HV contribution. For given $ϵ,\delta >0$ it identifies, with probability at least $(1-\delta )$, a solution with contribution at most $(1+ϵ)$ times the true minimal HV contribution. It is shown to work on very large problem instances with thousands of solutions and hundred dimensions.

[12] proposes a method to efficiently approximate a solution’s HV contribution using line segments.

[5] develop a neighborhood structure among local nadir points (also referred to as local upper bounds) in order to compute the entire nondominated set of a discrete multi-objective optimization problem in an efficient way. The neighborhood structure is updated with every new nondominated point. An advantage of this neighborhood structure is that once one local nadir point is known that has to be updated due to the insertion of a new point, one can easily navigate through the list of local nadir points to find all those that have to be updated in this iteration, too.

We concluded that the following combination of algorithms might be promising.

As a baseline, we could use the algorithm from [8] to quickly identify the solution with the minimum HV contribution. As a result, we obtain bounds of the HV contribution of each solution in the population. After removing the solution with the minimum HV contribution and adding a new solution, we would have to update the information to quickly identify the next solution with minimum HV contribution. This would be done in two steps:

1. 1.

   First, we would check whether there are any dominated solutions. If there is at least one dominated solution, we can remove one of them at random as all of their HV contributions are equal to zero, and thus minimal. If several dominated solutions are found, they can be removed iteratively without the need to recompute any HV contributions.

2. 2.

   If there are no dominated solutions, we can use the algorithm in [5] to identify the solutions whose HV contribution may have changed. Only those solutions need to have their bounds re-set, while the other solutions may keep their bounds from the previous iteration. With this update, the algorithm from [8] can be used again to quickly identify the solution with minimal HV contribution.

We also noted that it is easy to adapt the algorithm in [8] to work with approximations, rather than running it until the solution with the guaranteed smallest HV contribution remains. One could stop the algorithm earlier, while some intervals still overlap, with the tolerated overlap controlling the approximation error.

The number of solutions that need to be updated in one iteration of Step 2 may vary from one iteration to another, but it can generally not be large in all iterations.

Additional ideas discussed at the working group include

• $\mathrm{■}$

  The paper by [2] proposes an efficient algorithm to identify the solution with minimal HV contribution and approximation guarantee. It would be worthwhile to explore whether it is possible to speed this up by transferring computational results from one iteration to the next.

• $\mathrm{■}$

  For algorithms that add or remove more than one solution, [4] discusses possibilities to speed up computations.

• $\mathrm{■}$

  For ray-based approximations: when increasing the budget to increase accuracy, is it better to use more rays or more points along each ray?

• $\mathrm{■}$

  Can we cleverly employ GPUs for MC estimation approaches for HV? Which way to cut? Should each sample be evaluated in a core containing the reference set, or each reference set member on a core compared to a set of samples?

• $\mathrm{■}$

  How accurate do we need the estimation to be for effective use in (expensive) optimisation algorithms – i.e., what budgets do we need to employ. Is this fixed throughout the run or should it vary? Can we reproduce the observations from [1] for noisy single-objective problems: Accuracy is important at the start, not really in the middle, very important at the end?

• $\mathrm{■}$

  What is our budget for approximation given known computation time for fitness evaluation, i.e., should we rather spend more time on accurate HV computations, or work with crude approximations and instead do more iterations?

• $\mathrm{■}$

  If we only add points to the Pareto front, or only remove points, then we can make use of the fact that HV contributions can only decrease or increase, respectively. But only adding points will mean that we need to have an unconstrained size of the Pareto front.

This report was conducted as part of the “surrogates” subgroup of the Dagstuhl Seminar 23361 Multiobjective Optimization on a Budget, and the work benefited from wider discussions with Jonathan Fieldsend, Hisao Ishibuchi, Pascal Kerschke, Joshua Knowles, Boris Naujoks, Alma Rahat, Vanessa Volz, Hao Wang, Kaifeng Yang.

Multiobjective optimization refers to finding optimal solutions under the presence of multiple, often conflicting objective functions. Often, this is to be done within a predefined budget arising from resource limitations. The typical process consists of solving a multiobjective optimization problem specified by problem owners (POs) followed by decision makers (DMs) interacting with the solutions and studying the trade-offs between objective functions, and considering constraints (e.g., manufacturability of a mechanical apparatus), and ancillary information to identify solutions that meet their requirements. This is a rich field of study with many publications focusing on solving specific aspects of the whole process, see, for example, [1] for a survey on the various interactive multiobjective optimization methods, and [2] for a recent overview and survey of how to assess performance in this context. Nonetheless, there is a lack of an overall framework and taxonomy that allow discussion and contributions to be directed within a well-thought-out structure. As pointed out, e.g., in [3], the whole process with the problem formulation is rarely discussed in publications. In this report, we aim to propose an initial sketch for such a frame of reference that can be used as a tool by researchers to facilitate discussions and identify the scope of contributions within this structure. Next, we briefly discuss the proposed frame of reference.

Multiobjective optimization and related decision-making is an iterative process, and there are many elements to it. Firstly, we have the DMs or POs, who are attempting to solve a decision-making problem with many objective functions, supported by an (or a team of) analyst(s), or in other words an expert in optimization and decision-making methodologies. Together they are the elements of the human side of the process. Naturally, the other component of this process is then the computational side consisting of a mathematical or computational model devised through a specification, a suite of solvers appropriates for addressing the multiobjective optimization problem (MOP), and a module for information and knowledge extraction (e.g., for visualising the trade-offs) with an interactive user interface.

There are different forms of interactions between the elements. On the human side, the interactions are between the DM and the analyst. DMs use requests, to inform the analyst of what would be required of the model. Requests are related to the objective functions, constraints and decision variables. They can be reductionist (e.g., locate the best subset), expansive (e.g., include a new aspect), or a form of amendments (e.g., refinement or reformulations). Furthermore, the requests tend to be triggered by some recognition or need for change. The instigator in this case can be either the DM or the analyst. When the DM instigates, apart from the initial model specification, it may be because of new information becoming available, for example, a recognition of a new element of the problem that was not defined initially, a change in preference information relative to objective functions, new knowledge about the problem becoming available, or discrepancy with the real world requirements. On the other hand, an analyst may be able to instigate a change request after observing the behaviour of the DM. The latter may be automated through anomaly detection techniques.

At the boundary of the human-computer ecosystem, the analyst uses the requests to (re-)formulate model specifications and configure the appropriate solvers. On the other hand, the DMs interact with the computational side through a user interface where they can query the (estimated) Pareto front, and iteratively identify interesting (regions) of solutions. Moreover, they may instigate another run of the solver to investigate particular regions of interest through preference elicitation.

These core elements and interactions can be used to construct an frame of reference for knowledge exchange (FKE) that captures the solution process. An illustration of the IDF is shown in Figure 6. We now discuss the framework through an example.

Different interactions may be associated with different types of budgets, but they are primarily due to resource limitations. The DM must identify the ultimate solution within any such budget. Below, we present a few examples of budgets in this context.

Wall clock time.

The overall time before the DM must finalise their decision may be well-defined.

Model related budget.

Many real world problems would require substantial time for each function evaluation if it requires numerical simulations. For example, a computational fluid dynamic simulation of a draft tube may take a thousand seconds [6]. So, the budget may be about a few hundred function evaluations, and thus impact the interactions and their nature.

Solver budget.

The solver may itself be expensive. For example, entropy search for multiobjective optimization requires numerically approximating an acquisition function that can discriminate between solutions, but optimizing such acquisition functions to identify the next best solution may be exorbitant [7]. Thus, there are practical limitations on how much time we can spend before evaluating the next solution.

Proprietery software.

Many professional simulation software are proprietary, and therefore, there may be limits to how many licences are available to a DM.

Preparation budget.

Interactions between the DM and analysts may take an insignificant amount of time for discussions. In addition, for the analyst to evaluate and prepare models/solvers with DM guidance may take some time.

In this report, we briefly proposed and discussed a framework for iterative discovery of final solutions by a DM while interacting with multiobjective optimization methods supported by analysts. Future work involves expanding the framework and validating it with multiple real examples, with the possibility of incorporating multiple DMs working independently.

Optimization problems have often to be solved on a limited budget, e.g. due to time restriction or restrictions on the number of function evaluations. Multiobjective optimization problems are in general more costly to solve compared to single-objective optimization problems. Specifically, the number of objective functions and the type of objective functions directly influence the effort needed to solve the problem. Hence, in case of budget constraints for solving the problem, an important first step is to reduce the complexity of the multiobjective model as far as possible. In our working group, we have discussed and analyzed different ways to simplify a multiobjective optimization problem using reduction approaches.

We considered both reduction of structural complexity (e.g. by low-order polynomial approximation of objective functions) and reduction of problem size, in particular with respect to objective functions (e.g. by scalarization and aggregation of different objective functions). We refer to Figure 7 for an overview of different possible reduction scopes.

As mentioned above, reducing the complexity of a multiobjective problem is relevant in the context of optimization on a budget, with the idea of making better use of the limited number of evaluations available. The price of the reduction in general as well as measures for the quality of some reduction in particular are other interesting topics that deserve to be explored.

In this report, we present the main results achieved during the seminar:

• $\mathrm{■}$

  Reducing the number of objective functions by scalarizing some of the convex objectives – refer to Section 4.4.2

• $\mathrm{■}$

  Sufficient conditions for reducing the number of objectives locally using gradients – refer to Section 4.4.3

For completeness, we list other interesting discussion topics that we did not have time to delve into due to the limited time resources available:

• $\mathrm{■}$

  Reducing the complexity of the problem by (local) low-order polynomial approximation of objective functions (linearization or quadratic approximation)

• $\mathrm{■}$

  Best balance of high- and low-fidelity models during optimization

• $\mathrm{■}$

  Possible use of surrogate models to check local convexity

We aimed at the examination of the relation between the efficient set, i.e. the set of efficient solutions of a tri-objective problem and a family of bi-objective problems, obtained by combining two of the functions using weighted sums and by varying the weights. To be more specific, we examine the problem

with $f_i:{ℝ}^n\to ℝ$, $i\in [3]:=\{1,2,3\}$ continuous functions and a feasible set $S\subseteq {ℝ}^n$. The set of efficient solutions of (MOP) is denoted by $ℰ$ and the set of weakly efficient solutions by ${ℰ}_w$.

With this problem, we associate the family of problems

where $k\in [3]$ and $w\in [0,1]$. The set $\{i,j\}:=[3]\setminus \{k\}$ refers to the remaining indices, i.e. the indices other than $k$. The set of efficient solutions of (BOP($k,w$)) is denoted by $ℰ(k,w)$ and the set of weakly efficient solutions by ${ℰ}_w(k,w)$. In the literature, the effects of adding or deleting an objective function have been studied, see [2]. Moreover, for a fixed choice of $k$ and $w$, (BOP($k,w$)) and its relation to (MOP) has been addressed among others in [1]. Here, we examine the family of problems (BOP($k,w$)).

The following result relates the set of optimal solutions of the problems mention above to each other:

We start with the first inclusion. W.l.o.g. let $k=1$ and then $i=2$, $j=3$. Let $w\in (0,1)$ and $\overline{x}\in ℰ(1,w)$. Assume $\overline{x}\notin ℰ$. Then, there exists $\tilde{x}\in S$ with $f_1(\tilde{x})\le f_1(\overline{x})$, $f_2(\tilde{x})\le f_2(\overline{x})$, and $f_3(\tilde{x})\le f_3(\overline{x})$, with strict inequality for at lest one of the inequalities. Thus, it holds $w{f}_2(\tilde{x})+(1-w)f_3(\tilde{x})\le w{f}_2(\overline{x})+(1-w)f_3(\overline{x})$.

If , then this is a contradiction to $\overline{x}\in ℰ(1,w)$.

Otherwise, or holds true. We have in contradiction to $\overline{x}\in ℰ(1,w)$.

For the second inclusion, let $\overline{x}\in ℰ$. Then, due to convexity of $f(S)+{ℝ}_{+}^3$, there exists a vector $v\in {ℝ}_{+}^3\setminus \{0\}$ such that $\overline{x}$ is a minimal solution of

Set $\alpha :=v_2+v_3\ge 0$. If $\alpha =0$, then $v_2=v_3=0$, $v_1>0$ and $\overline{x}$ is a minimal solution of

and, thus, $\overline{x}\in {ℰ}_{\omega }(k,w)$. Otherwise, $\alpha >0$ and $\overline{x}$ is also a minimal solution of

and for $w:=\frac{v_2}{\alpha }\in [0,1]$, the point $\overline{x}$ is a minimal solution of a weighted sum ob the objectives $f_1(x)$ and $w{f}_2(x)+(1-w)f_3(x)$ with the two weights $v_1/\alpha \ge 0$ and $1$. Thus, $\overline{x}\in {ℰ}_{\omega }(k,w)$. ∎

However, the biobjective problems (BOP($k,w$)) are in general not capable of covering the full complexity of the three-objective problem (MOP) as the following counterexample within the next proposition shows.

For the following counterexample, we consider the points in outcome space. Let four feasible points in outcome space be given by $P_1={(5,0,0)}^{\top }$, $P_2={(0,5,0)}^{\top }$, $P_3={(0,0,5)}^{\top }$, $P_4={(4,4,4)}^{\top }$, see Figure 8. Then, $P_4$ is dominated for all problems of the form (BOP($k,w$)) (i.e., for all $k\in \{1,2,3\}$):

For $k=1$, the four points are mapped to $\left(\begin{array}{c}5\\ 0\end{array}\right)$, $\left(\begin{array}{c}0\\ 5w\end{array}\right)$, $\left(\begin{array}{c}0\\ 5\left(1-w\right)\end{array}\right)$, $\left(\begin{array}{c}4\\ 4\end{array}\right)$. Then, $\left(\begin{array}{c}4\\ 4\end{array}\right)$ is dominated for all $w\in [0,1]$ by $\left(\begin{array}{c}0\\ 5w\end{array}\right)$ or by $\left(\begin{array}{c}0\\ 5\left(1-w\right)\end{array}\right)$. The cases, $k=2,3$ yield the same result due to symmetry. ∎

The proof of Proposition 1 shows that unsupported efficient solutions may not be obtained as optimal solutions of one of the associated biobjective subproblems (BOP($k,w$)). Thus, additional convexity assumptions seem to be necessary. However, as the following proposition shows, it is not sufficient if only two of the three objectives are convex.

Consider the tri-objective problem (MOP) with $n=2$, $f_1(x)={(1-{(x_1-1)}^2-x_2)}^2$, $f_2(x)=x_1+\epsilon x_2$, $f_3(x)=2-x_1+\epsilon x_2$ for some $\epsilon >0$. Moreover, let $S=\{{(x_1,x_2)}^{\top }:x_1,x_2\ge 0\}$, i.e. we study

Then, $\overline{x}={(1,1)}^{\top }$, $\widehat{x}={(0,0)}^{\top }$ and $x^{\prime }={(2,0)}^{\top }$ are efficient solutions, since the non-negative objective function $f_1$ equals zero for all of them (i.e., $f_1(\overline{x})=f_1(\widehat{x})=f_1(x^{\prime })=0$). The corresponding vectors in the outcome space are

For $\overline{x},\widehat{x},x^{\prime }$, we determine the corresponding outcome vectors for the bi-objective optimization problem (BOP(1,$w$))

One can easily verify, that $\overline{x}$ is a dominated solution of (BOP(1,$w$)) for all values $w\in [0,1]$, since $\overline{x}$ is a dominated by $\widehat{x}$ if $w>\frac{1}{2}(1-\epsilon )$ or by $x^{\prime }$ if . ∎

Note that similar results can be easily shown for $k=2,3$, since the weighted sum in the second objective of (BOP($k,w$)) involves a potentially non-convex function.

Let us consider an algorithm for solving the multiobjective optimization problem (MOP) and assume that this algorithm relies on iteratively computing descent directions of the individual objective functions in a local optimization procedure. Then, if the gradient of one individual objective function is locally a convex combination of the others, this objective function does not have to be considered in the optimization process. More precisely, the following statement holds true.

We use the definition that a direction $d$ is a descent direction for a continuously differentiable function $g:{ℝ}^n\to ℝ$ in $\overline{x}$ if .

Using that $d$ is a descent direction for $f_1$ and $f_2$ we immediately get

∎

Some decisions within current (evolutionary) multiobjective optimization algorithms are solely based on rankings among solutions. The environmental selection step is one classical example. In the case where these decisions are costly, for example because the objective functions playing a role in the decision are expensive to evaluate, surrogate models for this decision process can be learned to replace the costly exact decision process by a cheaper, hopefully adequate enough process. This is in contrast with most existing surrogate-assisted (multiobjective) algorithms in which the objective functions are modelled by surrogates and the decision process is run on the surrogate-evaluated solutions.

In the subgroup on rank-based surrogates, we explored the possibilities of how rank-based surrogates can be used when rank-based decisions have to be made within an algorithm. As a first simple example, we propose to use a support vector machine to predict the environmental selection decisions in the NSGA-II algorithm [2]. In the following, we discuss the main ideas behind this rank-surrogate based NSGA-II and provide its pseudocode. Implementation, testing, and numerical benchmarking of the proposed algorithm remains a task for future work.

Following the ideas of the lq-CMA-ES algorithm [3] for single-objective surrogate-based optimization, we discuss a very basic multiobjective algorithm with a rank-based surrogate. The main working principle of the algorithm is to learn a rank-based surrogate model for the decision which solutions to keep and which to abandon at each step (environmental selection). The pseudocode in Algorithm 2 uses as baseline algorithm the well-known NSGA-II [2] and a Support Vector Machine (SVM, [1]) as the surrogate model but other multiobjective algorithms and surrogate models could be used instead as well.

In order to save as many expensive function evaluations as possible, we evaluate, in each iteration, only a small proportion $p_{\lambda }$ of the $\lambda$ newly sampled solutions. Based on the last $\lambda$ evaluated solutions in an archive $𝒜$ and the corresponding environmental selection decision on them, a surrogate model is learned to predict, for each solution $x_i$ of a new set of $\lambda$ solutions $x_1,\mathrm{\dots },x_{\lambda }$ whether the solution $x_i$ is to be kept (value 1) or not (value 0) for the next iteration. If the new model is predicting the environmental selection decisions well enough compared to the previous model (in terms of the Kendall-tau rank correlation coefficient, [4], line 15), we save the remaining function evaluations of the iteration and continue with the original algorithm, here NSGA-II. Only if the predictions between the old and the new surrogate model differ too much, i.e., if the Kendall rank correlation coefficient is smaller than a given threshold $T_{\tau }$, the next $\lceil p_{\lambda }\lambda \rceil$ solutions of the current iteration are evaluated successively on the true objective functions, and a new model is learned until the rank correlation coefficient between the old and the new model is larger than the target $T_{\tau }$ (or until all $\lambda$ solutions in the iteration are evaluated).

Note that in the pseudocode of Algorithm 2, the classifier $C_t$ returns $\mu$ solutions out of a set of $\mu +\lambda$ solutions (i.e., the ones that are selected for survival). To achieve this with SVMs, we return the $\mu$ solutions with the largest predicted value among all $\mu +\lambda$ solutions.

Rank-based decisions in randomized algorithms such as the environmental selection of evolutionary algorithms might profit from rank-based surrogates that are only trained on and can only provide rankings. If the algorithm is invariant to monotonous transformations of the objective functions, for example, it will keep this property with a rank-based surrogate. Our proposal of using a support vector machine to predict the environmental selection within NSGA-II is the first step towards this goal.

It remains to be shown that such a simple rank-based surrogate actually works in practice. Potentially, we need to feed more information to the surrogate model than just the decision of whether a solution is kept or not. Pairwise rankings between solutions or a total ranking on the input solution set can be imagined here.

Decision-making problems under conflicting objectives can be modeled and solved using multiobjective optimization. In general, such optimization problems with two or more objectives do not lead to a single optimal solution; instead, they lead to a whole set of so-called efficient solutions and corresponding nondominated points in the objective space. The set of all corresponding nondominated points is often called the Pareto front. Since each point cannot be improved in one or more objectives without worsening in at least one of the other objectives, all of these points are incomparable with respect to the optimization problem itself. Thus, usually in practice a decision maker must decide which solution is the best; such a decision should be made within the application context.

Decision-makers may bring inherent limitations to the optimization process, for example, in their availability and responsiveness. Furthermore, the set of nondominated points of a multiobjective optimization problem may become very large and, therefore, difficult for a decision maker to assort. In the following work, we present ideas that assist the decision-making process by creating a concise representation of the Pareto front that contains only a pre-determined number of nondominated points. That is, we assume the decision maker knows in advance a “budget” of how many nondominated points to consider when making a decision. The task is then to find a “good” representation of the Pareto front subject to the budget, so that the representation includes at most the given number of points. We consider a representation to be “good” if it is optimal with respect to some performance indicator of interest; further, the representation should not require advance computation of the complete Pareto front (or advance computation of a representation or approximation containing far more than the given number of points). The ability to compute such representations may improve decision-making by presenting a few options rather than all options, and may provide improved storage performance when multiple Pareto fronts must be computed as sub-problems of some larger multiobjective optimization (e.g., stochastic multistage) problem.

The remainder of this report is organized as follows. In the rest of the introduction, we provide a short overview on basic definitions and notation of multiobjective optimization and representation (Section 4.6.1.1), state the representation problem on a budget in a formal way (Section 4.6.1.2), and present previous work and existing literature (Section 4.6.1.3). In Section 4.6.2, we concentrate on the R2-indicator as quality measure for the representation. We provide an overview of the R2-indicator (Section 4.6.2.1), present a problem definition where the whole representation is computed “all-at-once” (Section 4.6.2.2), and present a mixed-integer reformulation of the problem (Section 4.6.2.3). The special case of multiobjective (mixed-integer) linear problems is discussed in Section 4.6.2.4. Computational results are included in Section 4.6.3, and Section 4.6.4 contains concluding remarks and future research.

The selection of an appropriate performance indicator is crucial when computing representations of the Pareto front. In this context, the R2-indicator has recently received attention since it is relatively easy to evaluate, and since it yields representations that also perform well with respect to the hypervolume indicator [7]. We provide a formal definition below and refer to [12] and [4] for a more detailed introduction to the concept.

We have implemented the MILP model (24)-(28) for a biobjective binary knapsack problem in AMPL and solved it with Gurobi 10.0.2. Figure 9 shows R2-optimal representations $ℛ$ for an instance of a biobjective binary knapsack problem with $30$ items together with the complete Pareto front ${𝒴}_{\text{Par}}$. We chose $T\in \{3,4,5\}$ and considered a weight set $𝒲$ given by

with $K=11$ for all shown representations. The reference point $u$ is set to the ideal point, shifted by a multiplicative factor of $1.01$.

In this work, we provide a closed form mixed integer programming formulation for the computation of (globally) optimal R2-representations for general multiobjective optimization problems. The complexity of the formulation depends on the underlying multiobjective optimization problem. For example, if the underlying problem is a multiobjective linear or mixed integer linear programming problem, then our formulation is a mixed integer linear programming problem that can be solved by available solvers.

Future research should discuss multiobjective problems with a more complex structure or with an increasing number of objective functions. In this case, the development of efficient solution heuristics, including iterative and greedy approaches, could be complemented by further improvements of the presented mixed integer programming formulations.

Given that, as it was mentioned above, the R2-indicator becomes equivalent to hypervolume under appropriate settings [7, 22], our formulation could probably also be adapted to finding an approximately optimal hypervolume representation. The representation would be only approximately optimal since we use a finite number of weight vectors. This would require adaptation of our model to a slightly different, inverse version of the Chebyshev function.

In this report, we focus on the all-at-once approach, which may be difficult for available solvers. If the model is too difficult for available solvers, this approach could also be easily adapted to locating the solutions one-by-one in a greedy manner. We would just need to treat values of already selected solutions as fixed parameters and optimize location(s) of just a single (or several) new solution(s). Since R2 most likely shares with the hypervolume the property of being a non-decreasing submodular function, the greedy approach would probably give some approximation guarantee. Another approach that could be used if the model is too difficult for available solvers would be to solve it with some single-objective metaheuristics or a hybrid approach combining metaheuristics with solvers called for some smaller subproblems. Furthermore, even if the proposed model could not be solved for some practical problems in acceptable time, it could still be applied to generate benchmarks for benchmarking heuristic methods aiming at finding a given number of efficient solutions.

This work has benefited substantially from Dagstuhl Seminar 23361: Multiobjective Optimization on a Budget, and from wider discussions with Matthias Ehrgott within the working group on representations. S. R. Hunter gratefully acknowledges support from the United States Air Force Office of Scientific Research under award number FA9550-23-1-0488.