Preference-based cone contraction algorithms for interactive evolutionary multiple objective optimization

Journal Pre-proof Preference-based cone contraction algorithms for interactive evolutionary multiple objective optimization Miłosz Kadziński, Michał K. Tomczyk, Roman Słowiński PII:

S2210-6502(18)30878-2

DOI:

https://doi.org/10.1016/j.swevo.2019.100602

Reference:

SWEVO 100602

To appear in:

Swarm and Evolutionary Computation BASE DATA

Received Date: 17 October 2018 Revised Date:

25 July 2019

Accepted Date: 23 October 2019

Please cite this article as: Mił. Kadziński, Michał.K. Tomczyk, R. Słowiński, Preference-based cone contraction algorithms for interactive evolutionary multiple objective optimization, Swarm and Evolutionary Computation BASE DATA (2019), doi: https://doi.org/10.1016/j.swevo.2019.100602. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Preference-based cone contraction algorithms for interactive evolutionary multiple objective optimization Milosz Kadzi´ nskia , Michal K. Tomczyka,∗, Roman Slowi´ nskia,b a Institute

of Computing Science, Poznan University of Technology, Piotrowo 2, 60-965 Pozna´ n, Poland Research Institute, Polish Academy of Sciences, Newelska 6, 01-447 Warsaw, Poland

b Systems

Abstract We introduce a family of interactive evolutionary algorithms for Multiple Objective Optimization (MOO). In the phase of preference elicitation, a Decision Maker (DM) is asked to compare some pairs of solutions from the current population. Such holistic information is represented by a set of compatible instances of achievement scalarizing or quasi-convex preference functions, which contribute to the construction of preference cones in the objective space. These cones are systematically contracted during the evolutionary search, because an incremental specification of the DM’s pairwise comparisons is progressively reducing the DM’s most preferred region in the objective space. An inclusion of evolved solutions in this region is used along with the dominance relation to revise an elitism principle of the employed optimization algorithm. In this way, the evolutionary search is focused on a subset of Pareto optimal solutions that are particularly interesting to the DM. We investigate, moreover, how the convergence is influenced by the use of some pre-defined and newly proposed self-adjusting (dynamic) interaction patterns. We also propose a new way for visualizing the progress of an evolutionary search. It supports understanding the differences between effects of a selection pressure imposed by various optimization algorithms. Keywords: Multiple objective optimization, Interactive evolutionary hybrids, Cone contraction, Pairwise comparisons, Interaction patterns

1. Introduction Evolutionary Multiple Objective Optimization Algorithms (EMOAs) are heuristic-based approaches which mimic the process of natural evolution in a population of solutions with the aim of solving an optimization problem involving multiple conflicting objectives [57, 73]. Such algorithms have been primarily designed for finding a set of well-distributed non-dominated solutions approximating an entire Pareto-front (PF). In this context, EMOAs have proven their usefulness in deriving a variety of efficient solutions to real-world MOO problems, such as supply chain design [36], pick-up and delivery [24], electoral redistricting [66], environmental/economic dispatch [54], or joint maintenance and production planning [55]. The focus of Interactive Multiple Objective Optimization (IMO) algorithms is quite different from EMOAs, because they involve DM’s preferences in the search process which aims at identifying a single most preferred solution. This is attained by means of alternating stages of decision making and optimization. In the course of an interaction with a DM, in the decision making phase, (s)he is asked to provide preference information which gives an indication for finding – in the subsequent optimization stage – the solutions which better correspond to his/her value system. These solutions are often called options or alternatives, as they compete for the best compromise choice. The major advantage of the interactive optimization algorithms derives from offering to the DM a possibility for systematically learning about the problem and its feasible solutions on the way to finding the most preferred ∗ Corresponding author: Institute of Computing Science, Pozna´ n University of Technology, Piotrowo 2, 60-965 Pozna´ n, Poland. Ph. +48-61 665 2950, Fax +48-61 8771525. Email addresses: [email protected] (Milosz Kadzi´ nski), [email protected] (Michal K. Tomczyk), [email protected] (Roman Slowi´ nski)

Preprint submitted to Swarm and Evolutionary Computation

October 25, 2019

option. In fact, in such an interactive process, both the DM and the machine learn, because the machine is adjusting the preference model to incoming preference information, which ensures convergence of the algorithm [3]. The incorporation of DM’s preferences into EMOAs can be conducted in different ways. On the one hand, these preferences can be specified a priori, and a suitably adjusted evolutionary algorithm can use them from the beginning of an optimization process to direct the search. The example forms of such preference information include a reference point [21, 50], extreme trade-offs [8], or desirability functions [68]. However, specification of preference information before the optimization may often be troublesome due to a limited DM’s knowledge on the solution space. On the other hand, a posteriori provided preference information can be associated with any EMOA approximating an entire PF. In this stream, the representative works are based on the dominance paradigm (e.g., NSGA-II [15] or SPEA2 [75]), indicator (e.g., IBEA [74] or HypE [1]), or decomposition (e.g., NSGA-III [13], MOEA/D [72]). In all above cases, the most preferred solution can be identified by progressively reducing, with the use of interactively provided preferences, a subset of efficient options returned by EMOA. Nonetheless, the quality of such an a posteriori optimization vastly depends on the standard of a finite set of solutions contained in the exploited population. In case its size is limited or the number of objectives is high, EMOAs may fail to successfully approximate an entire PF, which, in turn, negatively affects the final choice [7, 35, 36]. To address these problems, EMOAs have been revised to incorporate the interactively provided preference information already during the evolutionary search [5] (for a detailed survey on the interactive EMOAs, see [52, 71]). A progressive specification of preferences allows the DM to learn about the problem and its feasible solutions, making him/her convinced about the psychological convergence of the optimization process. The goal of using such interactive evolutionary hybrids is neither to generate a good approximation of an entire PF nor to find a single most preferred solution, but rather to identify the region(s) on PF containing the DM’s most favourable option(s). This is attained by using the DM’s preferences to bias the evolutionary search to a set of solutions that are particularly interesting for him/her. In this way, the DM is not required to analyse a large set of mostly irrelevant options, which gives him/her a better position to make a right decision. The information required by the preference-based EMOAs may take either direct or indirect form. The former corresponds to the values of parameters of an assumed preference model. The example forms of such a direct preference information include objective weights [17], comparison thresholds [14], aspiration or reservation levels [56, 70], parameters of a desirability function [68], a reference point [17, 48, 53, 61, 67], a target region [46], a size of territory preventing crowding [38], scores of solutions [44, 62], or fitness intervals assigned to the solutions contained in a small population [60]. However, a direct specification of such parameter values – even if some imprecision or uncertainty is admitted – is considered to be cognitively demanding for the DM. For this reason, the algorithms incorporating indirect preference information in form of decision examples are prevailing. In particular, in the course of the evolutionary search, the DM may be expected to compare pairs of solutions from the current population [7, 31, 51], order a small subset of solutions [16, 28], select the most and the least preferred options from such a subset [23, 25, 26], assign solutions to one of the pre-defined classes [27]. Such holistic statements reduce the cognitive load on the part of the DM, allowing to infer the compatible values of model parameters via preference disaggregation [30]. When using indirect preference information, there may exist many parameter values for which the preference model reconstructs the decision examples (for a detailed survey on preference modelling in MOO, see [69]). In other words, there may exist many instances of the preference model that would be compatible with the preference information. The way of dealing with such an indeterminacy of the model is essential for biasing the evolutionary search. Some algorithms [2, 7, 16, 51] select a single compatible model instance using some protocol, and use its indications to impose a selection pressure guiding the algorithm to the DM’s most preferred solution. This idea has been employed in view of using different preference models such as, e.g., additive [7, 51] or polynomial value functions [16], or a highly non-linear complex function [2]. However, since each compatible instance may rank the solutions that were not directly judged by the DM in a different way, the results of such an optimization algorithm vastly depend on the protocol of selecting a represen2

tative model instance. Consequently, a recent trend in the interactive evolutionary hybrids consists in accounting for all compatible preference model instances, hence incorporating the prudence principle in the optimization [4, 6]. Such a robustness concern can be implemented in two different though interrelated ways. On the one hand, all compatible preference model instances can be exploited to derive robust results (e.g., the necessary preference relation [6, 28], fronts of potential optimality [7, 63], or stochastic acceptability indices [64]) that are subsequently used to promote the solutions that proved their superiority over other available options with high certainty. This ensures a greater variability of solutions constructed by the algorithm in view of incompleteness of the DM’s preference information. On the other hand, the model instances compatible with such indirect preferences can be used to indicate the most desirable regions or to eliminate the irrelevant parts of the objective space. A popular way of implementing the latter approach consists in modelling the DM’s information by means of preference cones and guiding the optimization process with the cone contraction procedure. Preference cones were first introduced in [41] for dealing with IMO. The underlying assumption was that the DM’s preferences could be represented with some quasi-convex function, and a convex cone was generated based on the DM’s pairwise comparisons so that to span over all solutions judged worse than the solution indicated as a less preferred one in the pair compared by the DM. In this way, the set of admissible solutions can be reduced by excluding subsets of feasible options which are certainly not preferred over the solutions already identified as acceptable ones. This idea was subsequently extended in a few algorithms. In particular, [23] incorporated the preference cones within the interactive evolutionary hybrid to rank the solutions not directly compared by the DM. Furthermore, [37] adapted the preference cones for dealing with mixed-integer optimization. In addition, [59] used the preference cones as a tool for eliminating non-preferred solutions and additionally boosting the search conducted by the evolutionary algorithm by approximating the extreme gradient directions. Moreover, [39] showed how to generate more powerful cones, assuming that the DM’s decision policy can be represented with an Lα -norm, and incorporated the underlying technique for reducing the solution space to find the DM’s most preferred option. Yet different interpretation of preference cones for dealing with interactive optimization was introduced in [33]. In the proposed cone contraction algorithm, the DM’s pairwise comparisons are represented with a set of all compatible Achievement Scalarizing Functions (ASFs). The directions of the isoquants of these functions form a cone in the objective space, which is contracted with each DM’s judgement. In the course of an optimization, the solutions contained in the cone are considered as potential best options, whereas the solutions outside the preference cone are eliminated from further consideration. In what follows, we discuss the main contributions of the paper. They refer to proposing novel preference-based EMOAs using preference cones, their comprehensive evaluation on a set of benchmark problems and experimental comparison with the state-of-the-art interactive evolutionary hybrids, verifying an impact of the newly proposed interaction patterns and different evolutionary bases on the algorithms’ performance, and introducing new ways of visualizing the progress of an evolutionary search. Firstly, we propose a novel interactive EMOAs based on preference cones. During the evolutionary search, the DM is asked to compare some pairs of non-dominated solutions from the current population. Such holistic information is used to construct preference cones with the aim of converging a population of solutions towards a region of the PF being highly preferred to the DM. For this purpose, we use a procedure incorporating either all compatible ASFs [33] or distance-based Lα -norms [39]. The cones defined using these two preference models are employed to indicate, respectively, the desirable and non-relevant regions in the objective space. Nonetheless, they contract such a space in their unique ways. Since none of these procedures has been so far implemented within an EMOA, we are the first to explore the usefulness of thus constructed preference cones in the evolutionary search. Let us emphasize that although many recently proposed preference-based EMOAs are interactive, they fail to incorporate jointly the desired properties of conducting robustness analysis and not requiring the DM to directly express his/her preferences in form of precise values of model parameters. For example, I-MOEA/D-PLVF and I-NSGA-III-PLVF [44] are interactive and robust. However, the underlying robust optimization procedures involve a method-specific parametrization. Moreover, these approaches expect the DM to score a few solutions in a pop3

ulation, which is perceived to be more cognitively demanding. The additional parametrization is also required by, e.g., InDM2 [50], expecting the DM to interactively specify a reference point. In turn, ECC-MRW [35] and NEMO0 [7] are based on user-friendly pairwise comparisons. However, they employ only one representative preference model instance to direct the evolutionary search, hence failing to address the robustness concerns and being more susceptible to misguiding the search to a less preferred region on the Pareto front. This issue has been tackled by iMOEA/D [26], which expects the DM to indicate his/her best option in a small set of solutions and finds an entire region of the PF being highly preferred by the DM, instead of using just a single preference preference. Yet, similarly to I-MOEA/D-PLVF and I-NSGA-III-PLVF, iMOEA/D requires an additional parametrization to perform such a robust search. The methods proposed in this paper are most similar to NEMO-II and its counterpart for group decision optimization problems – NEMO-GROUP [34]. These methods are interactive, perform a robust search, and do not require any additional parametrization on the part of DM. The main algorithmic differences come from using a different preference model (i.e., Lα -norm rather than AVF) as well as guiding the evolutionary search by means of uniquely conducted robustness analysis by means of preference cones. The proposed algorithms are tested within an extensive study involving the DTLZ [18] and WFG [29] benchmark optimization problems with the number of objectives ranging from 2 to 5. The experimental study investigates which implementation of preference cones is more advantageous in view of finding a region of the PF being highly relevant to the DM. Furthermore, the best performing among the newly proposed interactive evolutionary hybrids is compared against the state-of-the-art optimization algorithms which incorporate the same form of indirect preference information [7, 35, 51] but differ in terms of the incorporated preference model and its exploitation. This comparison is oriented towards investigating an impact that both robustness preoccupation and consistency of the underlying preference model with the DM’s decision policy have on the algorithms’ performance. Let us observe that the majority of studies concerning interactive evolutionary hybrids assume that the questioning of the DM is performed at regular intervals. However, the quality of populations evolved by the preference-based EMOAs may vastly depend on when such interactions are conducted. This paper contributes to the literature by proposing novel interaction schemes, which distribute the questioning of the DM in different ways. Specifically, the interaction with the DM may be delayed until the population is already well-converged and well-distributed or the DM may be asked to compare more pairs of solutions within a single preference elicitation step. Furthermore, we propose a self-adjusting (dynamic) interaction pattern, which automatically determines when the DM should be questioned, based on the analysis of properties of the evolved population. Moreover, we investigate whether the use of a particular evolutionary base has a significant impact on the performance of the preference-based EMOAs. In the basic variant, the search conducted by the proposed algorithm is guided by considering both the inclusion of evolved solutions in the most preferred regions as well as the dominance relation, hence revising the front-based, generational scheme originally suggested in NSGA-II [15]. However, inspired by IEMO/I [65], we also propose an indicator-based counterpart of the interactive evolutionary hybrid using preference cones. This allows us to compare front- and indicator-based schemes incorporated into preference-based cone contraction algorithms in terms of the quality of constructed populations and computation time. Finally, we introduce a new way of visualizing the progress of an evolutionary search. It aggregates several independent runs of the algorithm, revealing when different regions in the objective space are discovered during the search. Analysis of such an expected spatial distribution of solutions allows for a better understanding of the selection pressure imposed by various optimization methods. The method is also adapted to represent different statistics which are of interest for depicting the course of optimization, e.g., crowding of the solutions in a particular region.

4

2. Search space reduction methods Let us consider a Multiple Objective Optimization problem defined as follows: M inimize {f1 (x), f2 (x), . . . , fM (x)} , subject to: x ∈ X,

(1)

where M is the number of objectives, x is a decision vector belonging to a non-empty feasible region X, and fi is an objective function used to evaluate x in terms of the i-th objective. A potential solution to such a problem is denoted by s = [f1 (x) = s1 , f2 (x) = s2 , . . . , fM (x) = sM ], and a set of all such solutions – Ω. A pair of solutions is related by the dominance relation ∆ (sa ∆sb ) iff ∀i∈{1,...,M } sai ≤ sbi ∧ ∃i∈{1,...,M } sai < sbi . Solution sa ∈ Ω is called Pareto-optimal if ¬∃sb ∈Ω,sb 6=sa sb ∆sa , and a set of all Pareto-optimal solutions is called a Pareto Front. In the remainder of this section, we recall a pair of cone contraction methods that will be used to interactively select the most preferred solution among a finite set of solutions. This can be attained by means of processing the pairwise comparisons systematically provided by the DM. Firstly, we recall a method, called CON EASF , which generates a preference cone by means of all ASFs compatible with such comparisons [33]. Secondly, we remind an approach, called CON ELα , which eliminates the solutions that are not recommended by any compatible instance of the Lα -norm, hence being located in the undesired region of the objective space [39]. When presenting the methods, we will assume that the set of R ordered pairwise comparisons provided by the DM at the current stage of interaction is denoted by H = {(sa  sb )1 , (sa  sb )2 , . . . , (sa  sb )R }, where  denotes the preference relation and the elements in H are ordered from the most to the least recently specified. Although these methods allow the DM to provide an arbitrary reference point z, in this paper we will assume that it is set to → − the ideal point z = 0 to focus solely on investigating the impact of the DM’s pairwise comparisons. 2.1. Interactive cone contraction method using Achievement Scalarizing Function The CON EASF method employs the DM’s pairwise comparisons to systematically contract a preference cone, indicating the most preferred region in the objective space [33]. The solutions contained in such a cone pass to the next iteration, hence being considered as the potential best options, whereas the solutions located outside the cone are eliminated. The algorithm is terminated once a single most preferred solution or a small subset of satisfying solutions is left. To explain the idea of constructing a preference cone, let us consider Figure 1a. It involves a set of solutions a {s , . . . , sg } and a single pairwise comparison se  sg provided by the DM. A preference cone formed by such a comparison is distinguished with a grey colour. It is composed of the directions of all isoquants of the ASFs originating in the reference point z = [0, 0], and giving advantage to se over sg , i.e., compatible with se  sg .

s

b

sc

0.6

0.8 0.6

f2

sd

a boundary isoquant t uan of on isoq ASF i t c ire ary ble 0.2 a d ound pati a b a com of

0.4

0.0 0.0

1.0

Solutions: remaining excluded

s

e

0.4

sf

0.4

0.6

0.8

Solutions: remaining excluded

sb sc sd

se sf

0.2

sg

eliminated region

0.2

sa eliminated region

0.8

sa

f2

1.0

0.0 0.0

1.0

f1

sg eliminated region 0.2

0.4

0.6

0.8

1.0

f1

(a) H = {(se  sg )1 }

(b) H = {(sc  sb )1 , (se  sg )2 }

Figure 1: CON EASF applied to a two-objective problem involving 7 example solutions sa = [0.1, 0.9], sb = [0.2, 0.8], sc = [0.4, 0.7], sd = [0.5, 0.5], se = [0.6, 0.4], sf = [0.7, 0.3], and sg = [0.8, 0.1], and two pairwise comparisons (a) se  sg and (b) sc  sb .

5

A boundary direction of an isoquant is depicted in Figure 1a. When it comes to the reduction of the set of potentially best options, solutions sf and sg , being exterior to the cone, are excluded from further consideration. With more comparisons, each one concerning a pair of solutions interior to the preference cone, the favourable space of the objective space is reduced. For example, when accounting for another example pairwise comparison sc  sb (see Figure 1b), the preference cone is contracted, and hence solutions sa and sb are eliminated. In what follows, we discuss how the algorithm verifies whether a particular solution is exterior or interior to the preference cone. Formally, CON EASF assumes that the DM’s preferences can be modelled with an ASF defined in the following way: M X d(s, w, z) = maxi=1,...,M {wi (si − zi )} + ρ (si − zi ) , (2) i=1

P where w is a normalized weight vector ( i=1,...,M wi = 1) being unknown to the method, z is a reference point, and ρ is an augmentation multiplier arbitrarily set to a small positive value. The lesser is the distance d(s, w, z) from s to z, the more preferred is s. Thus, when considering an example pairwise comparison sa  sb , it needs to imply: → − → − d(sa , w, 0 ) < d(sb , w, 0 ).

(3)

The above constraint can be reformulated as follows: sa  sb ⇒

∃

∀

i∈{1,...,M } j∈{1,...,M }

wj saj + γa,b < wi sbi ,

(4)


PM where γa,b = ρ i=1 sai − sbi . We call an ASF incorporating weight vector w = [w1 , w2 , . . . , wM ] compatible with the DM’s pairwise comparisons contained in H if:

∀

∃

∀

(sa sb )∈H i∈{1,...,M } j∈{1,...,M }

⇒


∀ w1 sa1 + γa,b < w1 sb1 ∧ (s s )∈H a

b

wj saj + γa,b < wi sbi ⇒

(5)



 w2 sa2 + γa,b < w1 sb1 ∧ . . . ∧ wM saM + γa,b < w1 sb1 ∨




 w1 sa1 + γa,b < w2 sb2 ∧ w2 sa2 + γa,b < w2 sb2 ∧ . . . ∧ wM saM + γa,b < w2 sb2 ∨ . . . ∨ 


 w1 sa1 + γa,b < wM sbM ∧ w2 sa2 + γa,b < wM sbM ∧ . . . ∧ wM saM + γa,b < wM sbM .



(6)

A set of directions of the isoquants of all ASFs reproducing the DM’s pairwise comparisons constitutes a preference cone. Example 1. Let us consider a pairwise comparison se = [0.6, 0.4]  sg = [0.8, 0.1] depicted in Figure 1. Then, following Eq. (6): se  sg =⇒ [(0.6w1 + γe,g < 0.8w1 ) ∧ (0.4w2 + γe,g < 0.8w1 )] ∨ [(0.6w1 + γe,g < 0.1w2 ) ∧ (0.4w2 + γe,g < 0.1w2 )] .

(7)

When assuming γe,g = 0 or negligibly small, a set of weight vectors which, when incorporated in ASF, are compatible with se  sg is defined as follows: w ∈ {[w1 , w2 ] : w1 > 0.5w2 ∧ w1 + w2 = 1}. To determine whether solution s is interior or exterior to the preference cone, one needs to determine a direction s ws = [w1s , . . . , wis , . . . , wM ] of the isoquant passing through s and z in the following way [33]: 1/si wis = PM , for i = 1, . . . , M . j=1 1/sj

(8)

If ASF parametrized with ws reproduces all pairwise comparisons in H (i.e., if Eq. (3) or, equivalently, Eq. (6) is 6

satisfied), s is interior to the preference cone. Otherwise, s is eliminated, because a direction of an isoquant passing through s compares the solutions differently than judged by the DM. Example 2. Let us continue Example 1 and focus on solution sd = [0.5, 0.5] depicted in Figure 1. Following Eq. (8), d d d d d ws = [w1s = 0.5, w2s = 0.5]. Note that constraint (7) is satisfied because w1s > 0.5 · w2s . Moreover, assuming d → d → d − − ρ = 0, d(se = [0.6, 0.4], ws , 0 ) = 0.3 < d(sg = [0.8, 0.1], ws , 0 ) = 0.4. Thus, an ASF parametrized with ws compares se and sg in the same way as the DM, and sd is interior to the preference cone. On the contrary, solution f f sf = [0.7, 0.3] is eliminated, because w1s = 0.3 ≤ 0.5 · w2s = 0.5 · 0.7. 2.2. Interactive cone contraction using the Lα -norm The CON ELα method assumes that the DM’s preferences can be modelled with an Lα function, which is defined as follows:  1/α α  P for α ∈ [1, ∞), i=1,...,M (wi |si − zi |) Lα (s, w, z) = (9) max {w |s − z |} for α = ∞, i=1,...,M i i i where α is a compensation parameter that can be either set a priori or adjusted during the method’s execution. In particular, L1 corresponds to a weighted sum with w being a vector of weights representing the trade-offs between criteria. On the other extreme, L∞ corresponds to a Chebyshev function. With greater α, it is more difficult to compensate a loss on one objective with a gain on another. To illustrate a reduction of the search space conducted by CON ELα by means of the DM’s pairwise comparisons, let us consider Figure 2a. When considering a pairwise comparison sc  se , a cone is spanned over all solutions → − that cannot be judged better than se for any compatible Lα function. Obviously, sc  se =⇒ Lα (sc , w, 0 ) < → − Lα (se , w, 0 ). Let us emphasize that unlike for CON EASF , a cone constructed by CON ELα indicates an undesired region in the objective space, and hence the solutions contained in it should be eliminated. To illustrate the impact that different values of α may have on the search space reduction, let us refer to Figure 2a. It demonstrates the cones implied by sc  se for α ∈ {1, 2, ∞}. For α = 1, solution sh is eliminated as there does not exist any compatible instance of the L1 -norm favouring sh over the non-preferred solution se . For α = 2, solutions sg and sh are rejected. Finally, for α = ∞, solutions sf , sg , and sh are eliminated on the base of the DM’s preferences. Hence, with greater α, more solutions are eliminated. Nonetheless, it is important to set α in a way that guarantees a high degree of consistency between Lα and the DM’s indirect judgements [39]. When assuming α = ∞ and considering an additional pairwise comparison sc  sb (see Figure 2b), solutions sa and sb are additionally excluded as they are located in the undesired region in the objective space.

s sd

f2

0.6

se

Solutions: excluded (α = 1) 0.2 excluded (α = 2) excluded (α = ∞) remaining

0.2

0.4

sf α=∞

0.4

0.0 0.0

α=1

0.6

sb sc sd

0.6

f2

c

0.8

sa

se

0.4

sg

0.8

sf

s

h

0.2 Solutions: 0.0 0.0

1.0

eliminated region

sb

eliminated region

0.8

1.0

sa

α=2

1.0

sg

sh

excluded (α = ∞) remaining

0.2

0.4

0.6

0.8

1.0

f1

f1

(b) H = {(sc  sb )1 , (sc  se )2 }

(a) H = {(sc  se )1 }

Figure 2: CDEMOLα applied to a two-objective problem involving 7 example solutions sa = [0.1, 0.9], sb = [0.2, 0.8], sc = [0.4, 0.7], sd = [0.5, 0.5], se = [0.6, 0.4], sf = [0.62, 0.25], sg = [0.7, 0.3], and sh = [0.8, 0.1]. The DM specifies (a) sc  se and (b) sc  sb .

7

In what follows, we remind a pair of constraint sets that should be considered to verify if solution s is interior to the cone [39]. In case α < ∞, one should verify if Eq. (10) is feasible; otherwise, one should check if Eq. (11) is feasible. The feasibility of the respective constraint set indicates that s is located inside the cone and should be eliminated. In turn, the infeasibility speaks in favour of s, identifying it as one of the potential best options at the current stage of interaction.

for all (sa  sb )j ∈ H :

∀ µj i∈{1,...,M }

h

sbi


α

 α i  α
α − sai ≤ si − sbi ,

(10)

µj ≥ 0;

for all (sa  sb )j ∈ H :

∀

b i∈{1,...,M }:sa i
sbi ≤ si .

(11)

3. Preference-based evolutionary cone contraction algorithms In this section, we introduce the basic variant of preference-based cone contraction algorithm for interactive evolutionary MOO. This algorithm asks the DM for pairwise comparisons of solutions from the current population, uses such holistic statements for constructing a preference cone, and drives the evolutionary search by means of both the indications of a preference cone and a dominance relation to the DM’s most preferred region in the objective space. In Section 3.1, we briefly remind NSGA-II [15], whose main idea of non-dominated fronts has been adopted with a suitable modification in the basic variant of the proposed algorithm. We chose the NSGA-II sorting scheme for the selection of parents based on the concept of non-dominated fronts because it is the most natural classification of solutions in a multiple objective optimization problem. When there is no preference information, this sorting scheme, together with the crowding distance criterion, make a pressure towards converging a well-spread population of solutions to the PF. When some preference information appears, the dominance relation is enriched, but the concept of non-dominated fronts still works with the only difference that now the pressure is towards a concentration of the population in the most preferred region of the PF. Adopting this scheme, we direct our primary focus on the decision aiding perspective rather than on the integration of more advanced evolutionary schemes. In this way, we can investigate which implementation of preference cones and which prioritization of sorting criteria are more advantageous in view of biasing an entire population to the region of PF being the most relevant for the DM. Moreover, implementation of a generational framework allows us to compare the proposed methods with the stateof-the-art algorithms, which incorporate the same type of preference information and a similar front-based scheme. Note that such an evolutionary strategy has been prevailing in the existing interactive evolutionary hybrids (see, e.g., [7], [35]). In Section 3.2, we discuss four interaction patterns, which decide upon when the questioning of the DM should be performed. In Section 3.3, we introduce a preference function, used by the proposed method for quantifying the consistency of solutions with the DM’s preference information. This function refers to the inclusion of solutions in the preference cone, which is contracted after each interaction with the DM. Finally, different variants of a novel interactive evolutionary algorithm incorporating a selected interaction pattern and a preference function are presented in Section 3.4. 3.1. The NSGA-II algorithm NSGA-II is a MOO algorithm evolving a population of N solutions Pg = [s1 , s2 , . . . , sN ], which is initially generated randomly or with some problem-specific method. In each generation g, NSGA-II alternates the following steps: (1) selection, (2) reproduction, (3) sorting, and (4) survival of the fittest, to approximate an entire PF. In the selection step, NSGA-II constructs a mating pool containing well-ranked solutions of population Pg . This can be performed by means of, e.g., random, proportional, or tournament selection. Then, NSGA-II applies genetic operators 8

on the selected solutions in order to construct an offspring population Pgo of size N . In this regard, mutation and crossover operators can be used for ensuring, respectively, exploration and exploitation in the evolutionary search. These can be, e.g., a polynomial mutation [12] or a simulated binary crossover [11]. Further, the algorithm sorts a combined population P = Pg ∪ Pgo , and selects the best N solutions, which are passing to next generation. The entire process is repeated for a pre-defined number of generations GEN . A sorting step performed by NSGA-II involves a pair of criteria oriented towards converging a well-spread population to the PF. Firstly, the algorithm divides a population P into the non-dominated fronts F. The first front F1 contains all non-dominated solutions in P , the second front F2 is composed of all non-dominated solutions in P \ F1 , etc. The lower is a number of the non-dominated front in which a solution is contained, the more preferred such a solution is. A secondary sorting criterion applied by NSGA-II is a crowding distance, being a measure of density of solutions surrounding a particular solution in the population. Within each non-dominated front, solutions with a greater crowding distance are preferred. 3.2. Interaction patterns In the course of an evolutionary search, the proposed optimization algorithm asks the DM to compare pairs of non-dominated solutions. Performing the preference elicitation in each generation would be unrealistic due to a high cognitive effort on the part of DM [49]. For this reason, we assume a limited budget IN T of the DM calls. Since IN T << GEN , a relevant research question consists in determining when the questions to the DM should be asked so as to avoid redundant interactions and maximize an information gain from the DM’s holistic judgements. Although in what follows we discuss four interaction patterns (denoted by START-ESISTR , ESIPCS , PP-ESIPOS , and RDS), when the algorithm is run, it should be parametrized to use only one of them. Let P AT T ERN (g, Pg , Pgo ) ∈ {0, 1, . . . , IN T } denote a function indicating a number of pairs the DM is expected to compare in the g-th generation, and EI – an elicitation interval corresponding to a number of generations the algorithm waits until the next interaction with the DM is performed. Specifically, we consider the following three pre-defined interaction patterns: • START-ESISTR : the DM is asked to compare ST R pairs of solutions in the first generation, and the remaining IN T − ST R interactions concerning a single pair of solutions are equally distributed over the optimization process:

P AT T ERN (g, Pg , Pgo )

= ST ART -ESI

ST R

(g, Pgp , Pgo )

=

   ST R  1    0

for g = 1, for g = 1 + (q + 1) · EI,

(12)

otherwise,

where EI = bGEN/(IN T − ST R + 1)e and q = 0, . . . , IN T − ST R − 1. This pattern is motivated by offering the algorithm more arguments to focus the population already at the early stage of optimization. • ESIPCS : the interactions with the DM are equally distributed over GEN generations, and (s)he is asked to compare P CS pairs of solutions in each interaction:  P CS P AT T ERN (g, Pg , Pgo ) = ESI P CS (g, Pg , Pgo ) = 0

for g = 1 + q · EI,

(13)

otherwise,

where EI = bGEN/(IN T /P CS)e and q = 0, . . . , IN T /P CS − 1. Specifically, application of ESI with P CS = 1 consists in asking the DM to compare a single pair of solutions in evenly distributed interactions. It is a default preference elicitation strategy to be used with the proposed interactive evolutionary algorithm. Note that the greater P CS, the less frequent, though more informative, is a given interaction with the DM who is asked to compare more pairs of solutions at the same time. 9

The main motivation for using ESI is to ensure a proper balance between interacting with the DM and performing the optimization. On the one hand, when the interval between subsequent interactions is too large, the algorithm may be prone to stagnation. On the other hand, when such an interval is too short, the method may not have enough time to properly converge the population, using fully the DM’s pairwise comparison(s). To avoid both problems in the context of eliciting a pre-defined number of pairwise comparisons, the ESI pattern imposes even distribution of the interactions with the DM. • PP-ESIPOS : the first interaction with the DM is postponed until the P OS-th generation, and the IN T interactions involving a single pairwise comparison are evenly distributed over the remaining GEN − P OS generations:  1 for g = P OS + 1 + q · EI, P AT T ERN (g, Pg , Pgo ) = P P -ESI P OS (g, Pgp , Pgo ) = 0 otherwise,

(14)

where EI = b(GEN −P OS)/IN T e and q = 0, . . . , IN T −1. This pattern is motivated by letting the algorithm construct a well-distributed population due to not involving any interaction at an early stage of the evolutionary search. Apart from the above three pre-defined (static) interaction patterns, we introduce a novel dynamic pattern, denoted by RDS. It determines when the DM’s questioning should be performed based on the progress attained in the evolutionary search. Specifically, we refer to a ratio of solutions in Pg being dominated by at least one offspring solution in Pgo , i.e.: ηg (Pg , Pgo ) =

|{s ∈ Pg : ∃so ∈Pgo so ∆s}| |Pg |

.

(15)

A small value of ηg (Pg , Pgo ) can be interpreted as a stagnation of the algorithm since offspring solutions do not dominate their parents. To identify the stagnation, one can verify if ηg (Pg , Pgo ) is smaller than some threshold th ∈ [0, 1], i.e., if ηg (Pg , Pgo ) < th. 1.0

1.0 1.0

µg (Pg , Pgo ) = 8/10

µg (Pg , Pgo ) o= 2/10

µg (Pg , Pg ) = 8/10

0.8

0.6

0.6 0.6 dominated parent solutions (Pg ) 0.4 0.4 nondominated parent solutions (Pg ) 0.2 0.2 offspring solutions (Pgo ) Pareto front Pareto front Pareto front 0.0 0.00.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

0.4 0.2 Pareto front 0.0 0.0

0.2

0.4

0.6

0.8

1.0

0.8

f2

f2

f2

0.8

f1

f1 (a) Interaction should be postponed

dominated dominated parent solutions(P(P parent solutions g )g ) nondominated nondominated parent solutions parent solutions(P(P g )g ) o o offspring solutions g )g ) offspring solutions(P(P Pareto front Pareto front

f1

(b) Interaction should be performed

Figure 3: Example parent and offspring populations and their corresponding ηg (Pg , Pgo ) values.

In Figure 3, we depict two example populations consisting of parent and offspring solutions (each sub-population of size 10), and their corresponding ηg (Pg , Pgo ) values. When it comes to Figure 3a, 8 parent solutions are dominated by at least one offspring solution. Hence ηg (Pg , Pgo ) = 8/10. A great value of ηg (Pg , Pgo ) indicates that the pace of the evolutionary search is high as the offspring population yields an improvement and thus many parent solutions would be removed when implementing the survival of the fittest. In this case, there is no need to supply the method with another DM’s pairwise comparison. In turn, Figure 3b illustrates a stagnation as only 2 out of 10 parent solutions are dominated by at least one offspring solution. In this case, it may be worth to interact with the DM, because 10

his/her another pairwise comparison may impose an additional evolutionary pressure. Such an interaction would be performed in case the employed threshold th would be greater than 0.2. Note that, since EMOAs are susceptible to some random fluctuations, we compute an average η over τ last generations, i.e.: ηgτ (Pg , Pgo )

=

g X

!, ηz (Pz , Pzo )

τ.

(16)

z=g−τ +1

In this paper, we consider τ = b1% · GEN e. Then, we define the following pattern:  1 if η τ (P , P o ) < th and |H| < IN T , g g g P AT T ERN (g, Pg , Pgo ) = RDS th (g, Pg , Pgo ) = 0 otherwise.

(17)

Hence, if ηgτ (Pgp , Pgo ) < th, the DM is asked to compare a pair of solutions unless a comprehensive number of already conducted comparisons is equal to IN T . Note, however, that in case ηg (Pg , Pgo )τ < th, it is also highly probable τ o ) < th. Consequently, to avoid an excessive questioning of the DM, the preference elicitation is that ηg+1 (Pg+1 , Pg+1 not performed for the next τ generations once ηg (Pg , Pgo )τ < th was triggered. 3.3. Ranking solutions with interactive cone contraction In this section, we introduce a preference function that will be incorporated in the proposed interactive evolutionary algorithm for quantifying the quality of each solution with respect to DM’s preferences revealed by progressively provided pairwise comparisons. This function is inspired by the mechanism implemented in the interactive cone contraction method based on ASF [33]. The underlying idea consisted in removing unsupported solutions after each incrementally provided pairwise comparison. This allowed to maintain during the optimization process only these solutions that could be considered as potentially most preferred to the DM, being interior to the preference cone implied by all preference information pieces provided so far. In this paper, we adapt this idea to the context of EMO and different preference models. Specifically, instead of removing the solutions which are inconsistent with some of the DM’s pairwise comparisons, we define a preference function that assigns fitness to solutions based on their compatibility with the incrementally supplied preference information pieces. Intuitively, the more pairwise comparisons provided from the beginning of the optimization process support a given solution, the more preferred it should be. Let us first define an auxiliary boolean function CON EM ODEL (s, sa  sb ), where M ODEL ∈ {ASF, Lα }. This function is instantiated with true if solution s is supported by the underlying cone, being consistent with pairwise comparison sa  sb . Specifically, CON EASF (s, sa  sb ) = true if s is contained in the preference cone, whereas CON ELα (s, sa  sb ) = true when s is exterior to the cone, hence not being eliminated. Finally, we define the preference function as follows: RAN KM ODEL (s, H) = argminj∈{1,...,R}

∀

(sa sb )k ∈H,k∈{j,...,R}

CON EM ODEL s, sa  sb



k

= true.

(18)

Overall, the quality of solution s is increased if it is located in the preference cone delimited by more oldest pairwise comparisons in H. Hence, lower values of RAN KM ODEL (s, H) are preferred. In particular, if all pieces of preference information in H supported s, then RAN KM ODEL (s, H) = 1. If solution s dropped out of the preference cone after the most recent pairwise comparison (sa  sb )1 , then RAN KM ODEL (s, H) = 2. On the other extreme, if s was supported only by the least recent pairwise comparison in H, then RAN KM ODEL (s, H) = |H|, whereas in case no decision example in H supported s, the result is RAN KM ODEL (s, H) = |H| + 1. The solutions with lower values of RAN KM ODEL (s, H) can be judged as more consistent with the DM’s pairwise comparisons. By promoting them during the evolutionary search, we aim at biasing the latter toward the DM’s highly preferred region of the PF.

11

RAN K ASF = 3

0.8

RAN K ASF = 2

1.0

RAN K ASF = 1

0.8

sd sc sf

0.4

sd sc

0.6

f2

f2

0.6

RAN KL∞ = 3

sa

0.4

sf

se 0.2

RAN K ASF = 4

0.0 0.0

0.2

0.4

Pareto front 0.6 0.8

RAN KL∞ = 2

RAN KL∞ = 4

1.0

sa se

sb

sb

0.2 RAN KL∞ = 1 0.0 0.0 0.2 0.4

1.0

f1

Pareto front 0.6 0.8

1.0

f1

(a) RAN KASF (s, H)

(b) RAN KL∞ (s, H)

Figure 4: Objective space partitioned by RAN KM ODEL (s, H) for H = {(se  sf )1 , (sc  sd )2 , (sa  sb )3 }. Table 1: The primary and secondary sorting criteria used by CDEMO and DCEMO (M ODEL ∈ {ASF, Lα }). Algorithm DCEMOM ODEL & CDEMOM ODEL if |H| = 0 CDEMOM ODEL if |H| > 0 DCEMOM ODEL if |H| > 0

primary-sort non-dominated fronts RAN KM ODEL non-dominated fronts

secondary-sort crowding-distance non-dominated fronts RAN KM ODEL

Example 3. Let us consider three ordered pairwise comparisons H = {(se  sf )1 , (sc  sd )2 , (sa  sb )3 } depicted in Figure 4. Let us remind that (se  sf )1 is the most recent pairwise comparison while (sa  sb )3 is the least recent one. In Figures 4a and 4b, we illustrate how functions RAN KASF and RAN KL∞ , respectively, divide the objective space in terms of the qualities of solutions measured with respect to preference information included in H. The darker the shade of grey is, the more preferred is the solution contained in the respective area. Additionally, each region is annotated with the corresponding fitness value. For example, in Figures 4a and 4b – RAN KM ODEL (sb , H) = 4, RAN KM ODEL (sd , H) = 3, RAN KM ODEL (sc , H) = 2, and RAN KM ODEL (se , H) = 1 for M ODEL ∈ {ASF, L∞ }. However, when it comes to sf , in Figure 4a – RAN KASF (sf , H) = 3, whereas in Figure 4b – RAN KL∞ (sf , H) = 2. 3.4. Preference-based cone contraction algorithm for interactive evolutionary multiple objective optimization In this section, we discuss the basic variant of the preference-based evolutionary cone contraction algorithm for MOO. Its sketch is presented as Algorithm 1. The solutions in the initial population Pg of size N are generated randomly (line 3). In the main loop (line 4), an offspring population Pgo is constructed (line 5) and combined with the parent population Pg (line 6) to form population P . Then, the algorithm verifies if the DM should be questioned. For this reason, P AT T ERN is used to determine the number P C of pairwise comparisons expected from the DM (line 7). If P C > 0, the method selects P C pairs of distinct non-dominated solutions from P (line 10), and asks the DM to indicate which one (s)he prefers in each pair (line 11). The answer is incorporated into the history H of preference elicitation (line 12). Having elicited the DM’s preference information, the algorithm selects N best solutions from P (lines 13-19). This process involves a pair of sorting criteria (lines 13-14). We distinguish two sorting schemes, called DCEMO and CDEMO, which are summarized in Table 1. In particular, DCEMO applies the non-dominated sorting as the primary criterion, while employing the cone-based preference function as the secondary criterion to rank solutions in the same non-dominated front. On the contrary, the order of sorting criteria used by CDEMO is inverse, hence giving the priority to the quality of solutions defined by the preference function rather than the dominance relation. Let us emphasize that until no pairwise comparison is elicited (H = ∅), the employed sorting criteria are the same as in NSGA-II (see Table 1), meaning that the search is not biased to any particular region. For the sake of brevity, we will refer to the two variants of the proposed algorithm by CDEMOM ODEL and DCEMOM ODEL 12

Algorithm 1 A general scheme of the proposed interactive evolutionary cone contraction algorithm. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14:

g=1 H=∅ Create initial population Pg repeat Pgo = make-new-pop(Pg ) P = Pg ∪ Pgo P C = P AT T ERN (g, Pg , Pgo ) if P C > 0 then for n = 1, . . . , P C do (sa , sb ) = randomly select a pair of non-dominated solutions from P (sa ? sb ) = ask the DM to indicate which solution (s)he prefers (? ⇔  or ? ⇔ ≺) H = (sa ? sb ) ∪ H P sort = primary-sort(P) for Pisort ∈ P sort do secondary-sort(Psort ) i

15: Pg+1 = ∅ and i = 1 16: while |Pg+1 | + |Pisort ∈ P sort | ≤ N do 17: Pg+1 = Pg+1 ∪ Pisort 18: i=i+1 19: Pg+1 = Pg+1 ∪ Pisort [1 : (N − |Pg+1 |)] 20: g =g+1 21: until g > GEN

f2

s

0.4

sc

e

s

sa

0.4

0.6

s ∈ F1 s ∈ F2 s ∈ F3

0.6

f

0.4

sd

b

RAN K ASF = 4 0.2

0.8

s

0.2 0.0 0.0

RAN K ASF = 1

sg

sc

RAN KL∞ = 2

sf

sd

sb

0.2

s ∈ F1 s ∈ F2 s ∈ F3

0.8

RAN KL∞ = 3

se

f2

sg

0.8 0.6

1.0

RAN K ASF = 2

RAN K ASF = 3

sa

RAN KL∞ = 1 0.0 0.0 0.2 0.4

1.0

RAN KL∞ = 4

1.0

0.6

0.8

1.0

f1

f1 (a) M ODEL = ASF

(b) M ODEL = L∞

Figure 5: Objective space partitioned by CDEMO and DCEMO incorporating RAN KM ODEL (s, H) for H defined in Example 3.

with M ODEL ∈ {ASF, Lα } indicating an applied CON E method. When selecting N best solutions which pass to the next generation (lines 16-19), the potential ties on the two sorting criteria are broken randomly. The main loop is repeated GEN times (line 21). In e-Appendix 1 (supplementary material available on-line), we discuss the computational complexities of different sorting procedures (see Table 1). Example 4. Let us illustrate the differences between the sorting procedures applied by the proposed variants of the cone contraction algorithm. In Figure 5, we depict a population composed of 7 solutions. The three non-dominated fronts are as follows: F1 = {sa , sb , sc }, F1 = {sd , se }, and F3 = {sf , sg }. The darker is the shade of grey, the more preferred is the solution contained in the respective area, in terms of RAN KM ODEL (s, H) for H defined in Example 3, M ODEL = ASF (Figure 5a) and M ODEL = Lα (Figure 5b). In Table 2, we report the pre-orders imposed by four different variants of the proposed method. In particular, although both CDEMOASF and DCEMOASF rank sb first, the second most preferred solution for CDEMOASF is sd , whereas DCEMOASF favours more sc . This is due to giving a priority to solutions contained in, respectively, the first non-dominated front or the preference cone consistent with all pairwise comparisons. 13

Table 2: Pre-orders of solutions imposed by different variants of CDEMO and DCEMO in Example 4. Algorithm CDEMOASF CDEMOLα DCEMOASF DCEMOLα

Pre-order (sb )  (sd )  (sf , sg )  (se )  (sc )  (sa ) (sa , sb , sc )  (sd )  (se )  (sf )  (sg ) (sb )  (sc )  (sa )  (sd )  (se )  (sf , sg ) (sa , sb , sc )  (sd )  (se )  (sf )  (sg )

4. Visualization of the Progress of Evolutionary Multiple Objective Optimization Visualization plays an essential role in understanding an optimization process and performance of a particular method. Specifically, in studies incorporating a posteriori methods, one illustrates a final population to prove that the PF has been approximated well (see, e.g., [13], [43], [72]). When it comes to the interactive methods, one visualizes populations obtained at different stages, which allows to capture a systematic focus on some region in the objective space, being more relevant to the DM according to the progressively specified preferences (see, e.g., [4], [34], [36]). However, both visualization concepts have some limitations in the context of studying the performance of EMOAs. Firstly, they reflect only a single optimization run. When using the evolutionary algorithms, which are susceptible to some random fluctuations, visualization of a single optimization may not be representative for the general performance of the method. Secondly, an illustration of just few selected optimization stages neglects a vast majority of solutions constructed throughout the evolutionary search. We propose a novel visualization method, called Visualization of the Progress of Evolutionary Multiple Objective Optimization (ViPEMO), tackling the two aforementioned limitations. In particular, it aggregates several independent runs of an EMOA, hence illustrating an average method’s performance (similarly to the attainment surfaces introduced in [22]). Moreover, it accounts for all solutions generated throughout the evolutionary search. As an outcome, the method provides a plot revealing when different regions in the objective space were discovered during the search. To distinguish between various generations, the method uses different shades or colours. An outline of the proposed method is presented as Algorithm 2. Although we focus on a two-objective case, the idea can be generalized to account for more objectives. ViPEMO incorporates bounds L and R (L < R) to constrain an objective space which is visualized in the plot. Then, it divides the domain of a particular objective into λ ∈ N+ equal segments of size Λ (line 1). Subsequently, the method computes a vector D of ordered segments (line 2), and uses it to divide the plot area into λ2 distinct subregions O (line 3). Solution s is contained in subregion Ow,h (s ∈ Ow,h ) if s1 ∈ [L1 + (w − 1)Λ1 , L1 + wΛ) ∧ s2 ∈ [L2 + (h − 1)Λ2 , L2 + hΛ2 ).

Algorithm 2 An algorithm for generating a expected spatial distribution of generated solutions.

1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

Input: EMOA; L, R – bounds of an objective domain; λ ∈ N+ – a number of buckets for each objective; T – a number of runs of the EMOA. Output: a plot Λi = (Ri − Li ) /λ for i = 1, 2 Split the objective domain: Di = [Di1 , Di2 , . . . , Diλ ] for i = 1, 2, where Dij = [Li + (j − 1)Λi , Li + jΛi ) 
Initialize the subregion set: O = D1 × D2 = D1w , D2h = Ow,h : D1w ∈ D1 ∧ D2h ∈ D2 for t = 1, . . . , T do t Initialize the generation matrix: Cw,h = ∞, for w = 1, . . . , λ, h = 1, . . . , λ Initialize EMOA for g = 1, . . . , GEN do Update population Pt maintained by EMOA for s ∈ Pt do for Ow,h ∈ O : s ∈ Ow,h do t t Cw,h = attribute-value(Cw,h , g)

12: Set the plot area according to L and R 13: for (w, h) ∈ {1, . . . , λ} × {1, . . . , λ} do t t 14: C 0 = {Cw,h : t ∈ {1, . . . , T }, Cw,h 6= ∞} 15: if |C 0 | 6= 0 then 16: Set a colour of the corresponding subregion according to some statistic of C 0

14

ViPEMO derives the results from T independent runs of the EMOA (line 4). Let us denote by C t a matrix t representing the results obtained in the t-th run. Its element Cw,h stores some computed statistic. In the default t implementation, Cw,h is set to g if the EMOA constructed a solution contained in subregion Ow,h after g generations. t t Initially, Cw,h is set to ∞ (line 5), and hence Cw,h = ∞ if subregion Ow,h was never discovered by the algorithm. t In each run, population P is evolved for GEN generations (lines 7-8). After each generation, ViPEMO identifies subregions Ow,h being occupied by different solutions (lines 9-10). Then, some pre-selected attribute value is comt puted and stored in Cw,h . In this paper, we tuned ViPEMO to identify when – for the first time during the search – different regions in the objective space were discovered by an EMOA. Having set a plot area according to the objective bounds L and R (line 13), ViPEMO iterates over all distinct regions Ow,h to aggregate results from the T independent runs (line 14). Thereafter, for each (w, h), the method t accounts for all Cw,h not instantiated with ∞ (line 15), hence identifying the runs in which region Ow,h contained t some solution. If there are no such Cw,h for (w, h) (line 16), Ow,h is left uncoloured. Otherwise, the method t computes some pre-selected statistic of Cw,h (line 17). In this paper, we compute a mean value of stored Cw,h . This value decides upon the colour of Ow,h according to some pre-defined colouring pattern. In e-Appendix 2, we discuss different extensions of ViPEMO along with some ways for enhancing comprehensibility of the generated plots. Since EMOAs are heuristics, it is not guaranteed that the newly generated offspring will be of high quality, even if parent solutions were considered to be satisfactory. In fact, due to being susceptible to some random fluctuations, EMOAs can sometimes generate poor solutions, which is typical for all heuristic-based methods. Such solutions can be considered as outliers. Obviously, due to their relatively lower quality, they would be quickly removed from the population as only the fittest solutions can survive from one generation to another. If these solutions were taken into account by ViPEMO, a clarity of the final plot would be negatively affected. To prevent this, one may employ some rules for eliminating these solutions from being visualized. In this study, we employ two such rules. Firstly, we eliminate solutions which survived for only one generation. The latter can be treated as a sign of their unsatisfactory quality as they did not compare positively to other solutions even for a very short period. In this regard, ViPEMO aims at illustrating the solutions whose quality proved to be satisfactory. This is in line with the survival of the fittest rule, which is a key concept of the evolution. Obviously, the exclusion period can be extended to more than one generation. In such a case, even more solutions would not be reflected in the final plot. Secondly, we colour only these regions in the objective space that contained some solution in at least 5% · T of the algorithm’s runs. The outliers that are sometimes generated by EMOAs can be located in different regions in the objective space, being far form the region currently exploited by the population. When running the algorithm many times, it is expected that such outliers would rarely occupy the same regions. Hence, we assumed that in case a region would not contain at least one solution for at least 5% · T of runs, we would leave it uncoloured, hence neglecting some potential outlier solutions. Clearly, the threshold can be set to a value greater or lesser than 5% · T , which would lead to neglecting, respectively, more or less solutions from the final plot. 5. Experimental evaluation In this section, we report the results of an extensive experimental evaluation of the proposed algorithm. In what follows, let us first discuss an experimental setting. Benchmark problems. We evaluated the proposed method on DTLZ1-7 [18] and WFG1-9 [29] benchmarks. The number k of position-related variables and the number l of distance-related variables were set according to [18], [29], and [45] (see Table 3). The problems involved from M = 2 to 5 objectives. Consideration of up to 5 objectives does not pose a cognitive burden to the DM who is asked to holistically compare the solutions pairwise [20]. A similar limitation was assumed in other related studies (see, e.g., [4], [47]). Evolutionary mechanisms. We used a tournament selection of size 5 to construct a mating pool. To generate an offspring, we used a polynomial mutation with a probability of 1/(k + l) and a simulated binary crossover with 15

a probability of 1.0. The distribution index of both operators was set to 10.0. The population size N was made dependent on the number of objectives M , and was set to N = 30(M − 1). The maximal numbers of generations GEN were set according to Table 3. Preference elicitation. We assumed a budget of DM calls IN T equal to 10 interactions. Similar limit was assumed in, e.g., [23], [47], and [58]. Obviously, the greater the number of interactions would be, the more information on the DM’s value system would be available to the algorithm. However, this number should be limited so as to avoid fatigue and motivational biases [49]. Whenever different is not explicitly stated, we use ESI 1 interaction pattern to decide upon when the DM should be questioned. Hence, EI = GEN/10 (see Table 3). To simulate the DM’s answers within the experiments, we constructed an artificial DM incorporating a predefined preference function. By default, we used for this purpose the L∞ function with a pre-defined weight vector wDM . This choice was motivated by a high consistency of L∞ with ASF for ρ ≈ 0. Thus, whichever cone contraction technique was incorporated by the proposed method, it was compatible with the DM’s decision policy (note that a different scenario not ensuring consistency between the two models is accounted for in Section 5.3). To simulate the interaction with the DM, we randomly selected a pair of non-dominated solutions sa and sb from the current population and used L∞ with weights wDM to indicate the more preferred solution. Evaluation strategy. For each test problem, we ran the proposed method 100 times. Each run involved an artificial DM using a different objective weight vector wDM . These vectors had been randomly pre-generated before the experimental evaluation. When interacting with different artificial DMs in various runs, the method was expected to converge to different regions of the PF. In this regard, such an experimental strategy exhibits a general performance of the method. Furthermore, the credibility of results is enhanced as the objective weight vectors are not selected arbitrarily. A similar evaluation strategy was used in [36], [47], and [65]. Quality indicators. For each run involving a unique DM, we reported the following quality indicators [36]: • Best Relative Score Difference (BRSD): a relative score difference between the optimal solution sopt of the artificial DM and the best constructed solution s ∈ P , according to the DM’s preference model: BRSD(P, w

DM

) = minsj ∈P

n . o − − − opt L∞ (sj , wDM , → 0 ) − L∞ (s , wDM , → 0) L∞ (sopt , wDM , → 0) ,

(19)

• Average Relative Score Difference (ARSD): an average relative score difference between the optimal solution sopt of the artificial DM and each solution s ∈ P , according to the DM’s preference model: ARSD(P, w

DM

)=

X  .  . − − − opt L∞ (sj , wDM , → 0 ) − L∞ (s , wDM , → 0) L∞ (sopt , wDM , → 0) |P |,

(20)

sj ∈P

Table 3: Characteristics of the benchmark optimization problems used in the experimental verification of the proposed method: the shape of PF, the number of objectives M ∈ {2, 3, 4, 5}, the number of position-related variables k, the number of distance-related variables l, the size N of population, the maximal number GEN of generations, the elicitation interval EI (* - the shape of PF is concave for M = 2 and 3, but according to [29], it is not known for M > 3; # – since we used a standard 15 decimal digits of precision, for WFG1 we set α of b poly function to 0.2). Test problem DTLZ1 DTLZ2 DTLZ3 DTLZ4 DTLZ5 DTLZ6 DTLZ7 WFG1# WFG2 WFG3 WFG4-9

Shape of PF Flat Concave Concave Concave * * Mixed, disconnected Mixed Convex, disconnected Flat Concave

M M M M M M

k −1 −1 −1 −1 −1 −1

M −1

2(M − 1)

2(M − 1)

2(M − 1) 2(M − 1)

16

l M +5−1 M + 10 − 1 M + 10 − 1 M + 10 − 1 M + 10 − 1 M + 10 − 1

GEN 1000 200 1000 200 200 1000 200

20

20

1500

150

M + 20 − 1

EI 100 20 100 20 20 100

20

500

50

20 20

500 500

50 50

where sopt denotes a Pareto optimal solution that attains the best score according to the DM’s preference function (i.e., Lα parametrized with wDM ). To generate such a solution for each benchmark problem, we made this problem easy to solve by setting the numbers of position- and distance-related variables to their minimal values (e.g., k = M −1 and l = 1 for DTLZ2). Then, we run 100 times a preference-based variant of NSGA-II with a crowding distance replaced with the pre-defined DM’s preference function, and assumed sopt is equivalent to the most advantageous solution identified across all runs. For each problem, we report absolute mean values of BRSD and ARSD for all accounted algorithms as well as the average ranks R they attained in terms of these two measures. Note that R = 1 means that a specific algorithm identified either the best single solution or the best average solution in all 100 independent runs. When presenting these results in a tabular form, we distinguish with a bold font the algorithm which attained the best rank as well as all other methods whose ranks did not differ significantly according to the Wilcoxon signed-rank test for 5% significance level. Furthermore, since the mean values of the quality indicators are very small, for clarity of presentation we multiplied them by 10ρ and reported the underlying ρ in the respective tables. In e-Appendix 1, we discuss the results of additional experiments concerning computation times of the proposed algorithms. In particular, we confirm that they scale well with the number of objectives, being applicable in a truly interactive manner, i.e., without forcing the DM to wait excessively long for the termination of subsequent optimization stages. 5.1. Experimental comparison of different variants of the interactive cone contraction evolutionary algorithm In this section, we report the results of an experimental study involving four variants of the proposed interactive EMOA. When it comes to CDEMO and DCEMO incorporating a cone contraction based on Lα , we set α = ∞. In Figure 6, we illustrate the average BRSD and ARSD attained throughout the evolutionary search for DTLZ1 and DTLZ2 with M = 3 and 5 objectives. Note that we use a logarithmic scale for the y-axis. The illustrations reveal that CDEMOASF and DCEMOASF outperformed their counterparts incorporating the cone contraction procedure 102

102

relative score difference

10

relative score difference

relative score difference

102

1

100 10−1 10−2 10−3

0

250

500

750

1 1010 1

1000

00

generation

0

10 102

relative score difference relative score difference

relative score difference

200 200

500 500

generation generation

750 750

1000 1000

(b) DTLZ1; M = 5

(a) DTLZ1; M = 3 10

CDEMO BRVD ASF CDEMO BRVD ASF CDEMO ARVD CDEMO ARVD ASF ASF DCEMO BRVD DCEMO BRVD ASF ASF DCEMO ARVD DCEMO ARVD ASF ASF CDEMO BRVD CDEMO BRVD L∞L∞ CDEMO ARVD CDEMO ARVD L∞L∞ DCEMO BRVD L∞L∞ DCEMO BRVD DCEMO L∞LARVD DCEMO ∞ ARVD

CDEMOASF BRVD CDEMOASF ARVD DCEMOASF BRVD 0 0 1010 DCEMOASF ARVD CDEMOL∞ BRVD CDEMOL∞ ARVD DCEMOL∞ BRVD −1 −1DCEMO 10 L∞ ARVD 10

10−1

0

101

CDEMO BRVD CDEMO ASFBRVD ASF CDEMO ARVD CDEMO ASFARVD ASF DCEMO BRVD DCEMO ASFBRVD ASF DCEMO ARVD DCEMO ASFARVD ASF CDEMO BRVD CDEMO BRVD L∞ L∞ CDEMO ARVD CDEMO ARVD L∞ L∞ DCEMO DCEMO BRVD L∞ L∞BRVD DCEMO L∞ DCEMO ARVD L∞ARVD

10−1CDEMOASF BRVD

CDEMOASF ARVD DCEMOASF BRVD 100 DCEMOASF ARVD CDEMOL∞ BRVD CDEMOL∞ ARVD 10−2DCEMOL∞ BRVD 10−1 DCEMOL∞ ARVD

10−2 0

50

100

generation

(c) DTLZ2; M = 3

150

200

00

200 50

100 500

generation generation

150 750

200 1000

(d) DTLZ2; M = 5

Figure 6: Convergence plots for average BRSD and ARSD for different algorithms applied to DTLZ1 and DTLZ2 with M = 3 and 5 objectives.

17

DTLZ1

Table 4: Average BRSD (first row) and ARSD (second row) for the populations generated in the last iteration by four algorithms for the DTLZ and WFG test problems with M = 2, . . . , 5 objectives. Average ranks R attained by the algorithms according to either BRSD or ARSD (ρ = 2).

M 2 3 4

DTLZ2

5 2 3 4

DTLZ3

5 2 3 4

DTLZ4

5 2 3 4

DTLZ5

5 2 3 4

DTLZ6

5 2 3 4

DTLZ7

5 2 3 4

WFG1

5 2 3 4

WFG2

5 2 3 4

WFG3

5 2 3 4

WFG4

5 2 3 4 5

CDEMOASF StD Mean 0.08 0.09 0.14 0.17 0.57 0.32 0.80 0.83 3.04 1.64 10.38 5.79 9.13 6.22 44.61 28.70 0.03 0.02 1.03 0.44 2.37 1.08 3.34 2.82 2.31 1.18 3.22 4.27 1.71 1.12 3.71 5.03 2.73 1.81 2.91 2.40 8.54 2.60 6.07 10.52 18.30 30.47 57.88 149.19 102.93 88.94 763.11 508.03 0.47 0.17 5.59 2.21 3.00 1.23 3.85 3.05 1.93 1.15 3.34 4.01 2.83 1.42 4.60 5.53 0.04 0.02 1.03 0.48 0.96 0.26 2.01 0.88 6.29 4.54 8.96 7.24 8.47 6.66 9.86 10.79 7.44 4.18 7.97 3.98 8.55 7.22 8.41 8.45 52.65 33.88 73.90 46.05 20.89 21.84 33.26 28.60 0.89 0.33 1.01 0.55 11.87 6.54 13.70 9.77 13.30 10.66 18.49 17.05 28.01 19.45 34.50 27.07 0.40 0.37 0.48 0.57 1.45 1.49 1.92 2.23 1.87 2.68 2.19 3.57 2.07 3.09 2.45 3.93 2.98 4.40 3.11 4.41 2.12 3.81 2.51 3.74 2.38 1.62 2.32 2.25 1.83 1.26 1.88 1.92 0.37 0.66 0.46 0.84 0.88 0.83 1.86 1.36 1.06 1.22 1.48 2.43 1.02 1.86 1.61 3.18 0.15 0.15 1.89 1.09 1.15 1.94 3.89 3.86 2.43 2.40 3.38 5.34 2.28 3.53 3.36 6.85

R 1.51 1.33 1.41 1.44 1.58 1.62 1.45 1.46 1.60 1.74 1.87 1.49 1.98 1.52 2.06 1.59 1.75 1.50 1.40 1.38 1.29 1.25 1.50 1.29 2.07 2.04 1.89 1.48 1.90 1.47 1.95 1.50 1.69 1.59 2.05 1.72 1.85 1.48 2.03 1.48 2.08 1.62 1.83 1.45 1.77 1.49 1.70 1.45 1.42 1.37 1.92 1.69 1.62 1.63 1.57 1.75 1.57 1.50 1.50 1.50 1.48 1.53 1.45 1.45 2.78 1.97 1.83 1.53 1.79 1.54 1.62 1.46 1.99 1.68 1.74 1.59 1.48 1.50 1.48 1.44 1.86 1.84 1.55 1.49 1.55 1.51 1.60 1.43

DCEMOASF StD Mean 0.11 0.13 0.19 0.23 0.86 0.34 3.90 1.38 1.21 0.89 8.37 3.82 16.81 10.02 90.50 45.24 0.08 0.02 1.23 0.61 2.34 1.05 3.71 3.12 2.38 1.32 3.48 4.31 1.75 1.16 2.98 4.24 1.89 2.03 2.67 2.96 17.26 6.97 46.39 338.25 42.50 88.95 130.59 191.98 116.50 96.29 995.57 732.27 0.17 0.05 0.81 0.56 2.85 1.38 3.92 3.56 2.01 1.16 3.50 4.40 2.30 1.46 4.01 5.21 0.04 0.02 1.20 0.65 7.17 2.14 10.60 3.84 10.42 5.03 13.19 9.74 5.69 4.77 10.59 11.62 7.44 4.19 8.25 3.86 10.52 7.46 12.60 11.89 83.75 36.70 95.30 47.00 24.40 23.65 41.43 34.47 4.78 1.48 7.71 3.90 15.15 9.79 15.66 12.35 14.26 10.03 18.21 15.32 27.87 19.26 32.53 24.28 0.39 0.28 0.57 0.72 1.30 1.48 1.83 2.32 1.79 2.58 2.05 3.40 1.87 3.11 2.11 3.87 2.49 3.94 3.37 4.20 2.46 4.04 3.03 4.38 2.78 1.81 2.80 2.48 2.71 1.77 2.70 2.39 0.52 0.64 0.67 0.92 0.86 0.76 1.52 1.17 1.15 1.32 1.50 2.58 1.08 1.93 1.47 3.34 0.13 0.16 2.01 1.38 1.25 2.24 4.01 3.76 2.66 2.57 3.54 5.51 2.36 3.61 2.99 7.07

R 1.98 1.78 1.62 1.57 1.42 1.38 1.57 1.54 1.98 2.00 1.85 1.59 2.00 1.48 2.09 1.42 1.85 1.66 1.60 1.63 1.72 1.75 1.55 1.71 1.70 1.94 2.05 1.66 1.91 1.56 2.11 1.53 1.88 2.12 2.04 2.00 1.78 1.69 1.68 1.52 2.07 1.80 1.62 1.60 1.76 1.64 1.78 1.55 2.42 2.41 2.04 1.89 1.55 1.54 1.44 1.52 1.60 1.64 1.50 1.50 1.53 1.47 1.55 1.55 2.42 2.36 1.88 1.79 1.86 1.68 1.79 1.65 1.96 1.80 1.56 1.43 1.62 1.50 1.58 1.57 1.83 2.18 1.62 1.51 1.74 1.49 1.74 1.60

18

CDEMOLα StD Mean 0.81 0.62 3.06 2.98 18.93 18.35 95.66 109.38 44.20 49.14 321.85 604.71 90.05 95.61 633.48 1733.19 0.23 0.14 2.56 1.89 3.42 2.06 9.76 16.65 2.08 1.69 6.04 20.10 1.06 1.21 4.38 20.21 34.89 17.56 55.02 36.96 335.41 411.03 1658.82 840.43 552.93 465.70 5286.51 2303.74 736.92 790.53 10077.27 2892.83 2.48 0.56 3.46 2.49 2.34 1.48 8.78 15.28 2.06 1.64 6.45 19.81 1.17 1.19 4.56 21.06 0.26 0.14 2.26 1.99 0.46 0.13 4.09 3.07 10.52 10.13 15.83 20.98 8.83 9.79 18.65 27.49 8.65 5.22 16.30 11.85 11.19 20.22 27.84 70.67 80.84 88.86 155.61 204.70 50.39 28.97 228.67 101.69 11.21 2.76 11.29 4.72 11.21 9.60 25.45 25.13 30.31 27.87 43.57 45.57 90.79 64.74 126.21 90.99 1.34 0.99 2.04 2.96 3.83 5.50 5.95 10.78 4.45 7.22 6.12 11.58 4.52 7.14 5.97 10.67 2.08 3.26 3.71 3.73 2.93 3.92 7.09 4.59 2.70 2.89 3.45 7.53 1.74 2.59 2.71 6.80 0.66 1.01 1.52 2.35 1.97 2.10 8.77 5.03 2.38 3.82 4.92 12.84 2.02 4.31 3.17 11.60 0.79 0.45 3.02 2.26 3.87 3.90 8.48 15.80 4.01 5.29 6.44 16.77 3.82 6.07 7.81 16.09

R 3.45 3.41 3.46 3.47 3.49 3.53 3.55 3.54 3.13 3.01 3.21 3.51 3.12 3.53 2.94 3.55 3.02 3.26 3.54 3.45 3.38 3.49 3.46 3.50 2.92 2.97 2.94 3.47 3.07 3.53 2.90 3.44 3.21 3.17 2.98 3.24 3.16 3.32 3.23 3.50 2.96 3.31 3.29 3.57 3.31 3.49 3.28 3.54 2.79 2.99 2.79 3.17 3.45 3.43 3.72 3.48 3.49 3.53 3.44 3.44 3.61 3.58 3.49 3.49 2.49 2.79 3.23 3.41 3.19 3.50 3.33 3.48 2.82 3.37 3.30 3.47 3.55 3.64 3.42 3.47 3.23 3.01 3.45 3.51 3.29 3.47 3.26 3.49

DCEMOLα StD Mean 0.57 0.48 2.03 2.52 20.59 19.85 126.89 111.03 41.17 48.20 299.20 574.33 79.10 89.47 634.98 1649.28 0.86 0.27 4.90 3.06 1.53 1.23 8.76 15.01 1.30 1.29 5.82 20.49 0.97 1.24 4.57 20.38 53.98 21.43 73.48 47.01 362.21 413.45 1704.37 1021.86 643.60 585.22 5374.66 2013.31 660.91 762.81 10070.96 2576.34 18.29 5.15 18.21 7.21 2.22 1.50 10.37 14.99 1.70 1.51 6.25 19.23 1.53 1.35 5.17 21.60 0.34 0.17 3.01 2.10 1.30 0.31 4.83 3.24 11.06 11.14 15.95 21.79 8.84 9.88 18.10 26.80 8.09 4.95 14.83 9.79 11.50 19.22 25.44 65.75 88.03 84.47 140.39 188.45 52.57 32.77 217.98 87.03 16.51 8.57 16.42 10.80 28.97 18.56 35.89 30.84 34.49 28.63 45.79 45.54 78.95 55.92 108.09 80.64 1.28 0.77 1.82 2.95 3.87 5.74 5.55 10.99 4.56 7.12 5.96 11.30 4.48 7.13 5.80 10.70 2.75 4.57 4.62 4.83 2.51 3.62 6.62 4.70 1.94 2.27 2.70 7.11 2.23 2.56 3.05 6.94 0.78 1.18 1.55 2.34 1.89 1.99 8.62 4.50 1.74 3.21 3.23 11.06 1.83 4.30 3.24 11.48 0.69 0.44 3.68 2.83 3.80 4.28 9.47 15.87 3.78 5.48 7.31 17.69 3.40 6.11 6.72 16.00

R 3.06 3.48 3.51 3.52 3.51 3.47 3.43 3.46 3.29 3.25 3.07 3.41 2.90 3.47 2.91 3.44 3.38 3.58 3.46 3.54 3.61 3.51 3.49 3.50 3.31 3.05 3.12 3.39 3.12 3.44 3.04 3.53 3.22 3.12 2.93 3.04 3.21 3.51 3.06 3.50 2.89 3.27 3.26 3.38 3.16 3.38 3.24 3.46 3.37 3.23 3.25 3.25 3.38 3.40 3.27 3.25 3.34 3.33 3.56 3.56 3.38 3.42 3.51 3.51 2.31 2.88 3.06 3.27 3.16 3.28 3.26 3.41 3.23 3.15 3.40 3.51 3.35 3.36 3.52 3.52 3.08 2.97 3.38 3.49 3.42 3.53 3.40 3.48

WFG5

Table 4: (continued)

M 2 3 4

WFG6

5 2 3 4

WFG7

5 2 3 4

WFG8

5 2 3 4

WFG9

5 2 3 4

Average ranks (R)

5 2 3 4 5 All

CDEMOASF Mean R StD 1.34 1.85 1.79 1.89 2.32 1.87 1.68 2.29 1.68 4.17 4.01 1.55 2.98 2.79 1.71 5.99 3.51 1.53 3.99 2.39 1.50 7.27 3.08 1.52 1.27 1.96 2.21 2.17 3.49 1.77 2.09 3.08 1.78 4.59 4.45 1.63 3.19 3.10 1.68 6.43 4.28 1.53 4.43 2.35 1.60 7.79 3.03 1.53 0.09 0.15 1.82 0.77 1.27 1.76 2.43 3.91 1.82 4.58 4.68 1.49 2.44 2.50 1.62 5.27 3.29 1.51 3.49 2.05 1.48 6.84 2.92 1.49 1.26 0.88 1.68 1.88 1.74 1.64 2.26 2.38 1.67 4.15 3.51 1.48 3.30 3.10 1.61 5.66 3.48 1.41 3.74 2.16 1.62 7.09 3.07 1.55 0.70 1.25 1.61 1.37 2.28 1.57 1.78 2.76 1.68 4.46 4.64 1.52 2.50 2.23 1.49 5.78 3.15 1.42 3.79 2.14 1.63 7.90 3.17 1.52 Mean StD 1.84 0.33 1.67 0.20 1.73 0.19 1.53 0.09 1.65 0.18 1.50 0.09 1.64 0.20 1.49 0.10 1.71 0.24 1.55 0.15

DCEMOASF Mean R StD 1.32 1.79 1.62 1.76 2.11 1.74 1.95 3.16 1.74 4.48 4.41 1.55 2.76 2.73 1.48 5.51 3.69 1.47 3.93 2.32 1.59 7.17 2.86 1.48 0.88 0.55 2.84 1.64 1.35 2.08 1.99 3.44 1.71 3.96 4.10 1.47 3.13 2.74 1.64 6.32 3.98 1.47 4.32 2.44 1.53 7.62 3.03 1.47 0.19 0.83 1.61 1.10 2.28 1.84 1.96 3.16 1.72 4.40 4.39 1.58 2.36 2.36 1.64 5.16 2.98 1.50 3.74 2.27 1.58 6.98 2.96 1.51 1.36 0.95 2.16 2.04 1.57 2.06 2.66 3.45 1.69 4.67 4.33 1.55 3.52 3.02 1.72 6.77 4.16 1.60 3.66 2.11 1.48 6.85 3.02 1.45 0.98 1.94 2.11 2.06 3.79 1.97 1.74 2.75 1.66 4.47 4.51 1.48 2.58 2.35 1.69 6.36 3.79 1.58 3.87 2.11 1.79 7.68 3.02 1.48 Mean StD 2.00 0.34 1.96 0.23 1.74 0.18 1.61 0.16 1.69 0.16 1.55 0.10 1.68 0.20 1.53 0.07 1.78 0.26 1.66 0.23

CDEMOLα Mean R StD 1.50 1.99 3.28 3.45 3.38 3.27 3.28 3.32 3.15 13.39 6.28 3.31 6.51 3.53 3.51 21.78 7.17 3.56 7.17 3.07 3.34 19.57 5.52 3.49 0.85 0.62 1.98 3.11 3.39 3.03 3.19 3.21 3.15 14.76 7.66 3.25 5.93 3.48 3.37 21.39 6.54 3.53 7.87 3.44 3.40 21.89 6.14 3.54 0.90 2.33 3.41 4.22 5.86 3.42 4.32 4.55 3.02 15.21 9.19 3.39 5.44 3.15 3.49 19.34 6.13 3.46 6.82 2.99 3.49 20.36 6.51 3.51 1.58 1.15 2.88 3.74 3.43 3.10 4.78 4.23 3.30 15.90 8.68 3.55 6.59 3.75 3.35 19.72 6.35 3.55 7.05 3.07 3.46 19.75 5.55 3.63 1.61 3.15 3.01 3.98 5.31 3.11 4.45 4.14 3.32 19.20 8.55 3.66 5.18 3.36 3.33 20.53 6.83 3.48 6.28 2.77 3.28 20.56 5.90 3.47 Mean StD 3.00 0.38 3.17 0.20 3.22 0.21 3.43 0.13 3.35 0.16 3.51 0.07 3.35 0.21 3.51 0.04 3.23 0.29 3.40 0.19

DCEMOLα Mean R StD 1.45 1.93 3.31 2.69 2.58 3.12 3.93 4.09 3.43 15.70 8.34 3.59 6.36 3.57 3.30 20.21 6.90 3.44 7.28 3.25 3.57 20.08 6.06 3.51 1.36 2.61 2.97 3.76 4.54 3.12 4.30 4.37 3.36 18.27 9.29 3.65 5.81 3.27 3.31 20.81 6.29 3.47 7.89 3.25 3.47 22.01 6.48 3.46 0.49 1.23 3.16 3.15 4.49 2.98 4.92 5.13 3.44 16.38 9.33 3.54 5.16 3.16 3.25 19.38 6.07 3.53 6.83 3.17 3.45 20.21 6.81 3.49 1.69 1.46 3.28 4.20 4.06 3.20 4.78 4.13 3.34 15.82 8.58 3.42 6.41 3.78 3.32 19.48 6.97 3.44 6.77 2.66 3.44 19.27 5.77 3.37 1.83 2.44 3.27 5.28 5.79 3.35 4.64 4.41 3.34 17.29 8.26 3.34 5.36 3.20 3.49 21.74 7.12 3.52 6.29 2.69 3.30 20.79 6.25 3.53 Mean StD 3.15 0.27 3.19 0.19 3.31 0.18 3.43 0.15 3.30 0.17 3.45 0.07 3.34 0.19 3.46 0.07 3.28 0.21 3.38 0.17

based on L∞ . This suggests that CDEMOL∞ and DCEMOL∞ cannot impose as sufficient evolutionary pressure as CDEMOASF and DCEMOASF , and thus their convergence is slower. The differences in performances are more evident for DTLZ1 (which is more challenging when compared to DTLZ2 [18]) for M = 3 and 5 objectives, where the average qualities of solutions contained in the final populations constructed by CDEMOASF and DCEMOASF were better than the qualities of the best solutions identified by CDEMOL∞ and DCEMOL∞ . Let us conclude that the convergence plots do not reveal any difference between CDEMO and DCEMO incorporating the same cone contraction technique. In Table 4, we provide comprehensive results concerning 16 test problems for M = 2, 3, 4, and 5 objectives. They refer to the quality of solutions contained in the population constructed in the final generation. Apart from reporting the outcomes for each problem individually, we computed the ranks averaged across all benchmarks considered jointly. The results confirm the superiority of CDEMOASF and DCEMOASF over the remaining algorithms. Specifically, CDEMOASF attained the average ranks of 1.71 (BRSD) and 1.55 (ARSD) for all considered problem instances, while the respective ranks for DCEMOASF are 1.78 (BRSD) and 1.66 (ARSD). These two variants performed similarly for the vast majority of benchmarks with the proviso that CDEMOASF performed slightly better than DCEMOASF when few objectives were involved (see, e.g., DTLZ1-2, DTLZ7, WFG6, WFG8-9 with M = 2, and DTLZ3 with M = 2, 3, 4). These problems are easier to solve, diminishing the role of non-dominance sorting. As a result, prioritizing the solutions highly preferred to the DM over non-dominated ones may be advantageous. This observation is supported by the average ranks attained by these methods according to BRSD and ARSD, for problems 19

with M = 2 objectives. In particular, CDEMOASF attained the average ranks of 1.84 (BRSD) and 1.67 (ARSD), while DCEMOASF attained slightly worse ranks of 2.00 (BRSD) and 1.96 (ARSD). For problems with M = 4 or 5 objectives, the dominance relation is scarce, increasing the role of the CON E method in differentiating the solutions. This, in turn, implies that CDEMO and DCEMO sort the solutions in a similar way, which explains the lack of any significant difference in the performance of these methods. The average ranks attained by DCEMOL∞ and CDEMOL∞ are vastly worse than for their counterparts using ASF. The proposed method incorporating cone contraction based on L∞ does not impose a sufficient evolutionary pressure when tackling difficult problems involving many objectives. To support this claim, let us compare ARSDs attained by CDEMOASF and CDEMOL∞ when applied to DTLZ1. For M = 2, these methods attained mean ARSD of 0.0017 and 0.0298, respectively. The difference in performance of both variants of CDEMO becomes more evident with the increase in the number of objectives. For example, CDEMOASF attained mean ARSD of 0.2879 for M = 5, whereas the respective score for CDEMOL∞ was equal to 17.3319. 5.2. Visualization of solutions constructed by the proposed method In this section, we provide examples of populations constructed by the considered algorithms, hence supporting the understanding of numerical results. The visualizations have been derived from applying ViPEMO to DTLZ2 with the objective bounds L0 = L1 = 0.0 and R0 = R1 = 1.5, and a partition of the objective space into λ2 = 10000 equal sub-areas. To emphasize the differences in the evolutionary search performed by various methods, we slowed down its pace by setting the number of distance-related variables l to 50. The algorithms were run T = 100 times for the two objective weight vectors wDM = [0.5, 0.5] and wDM = [0.7, 0.3] of the simulated DM. The generated expected spatial distributions of the evolutionary optimization are presented in Figure 7. We used four distinct shades of grey to differentiate between generations. The area containing an optimal Pareto optimal solution is distinguished with a black circle, whereas the PF is marked with a dashed lined. In e-Appendix 2, we present the plots using richer colouring patterns and illustrating other statistics that provide additional insights on the algorithm’s performance. Let us compare the expected spatial distributions for CDEMOASF and DCEMOASF . The solutions constructed by these algorithms progressively converged towards the optimal solution (see Figures 7a, 7b, 7e and 7f). The distributions support the claim that when dealing with few objectives DCEMOASF may attain slightly worse results than CDEMOASF . Indeed, the population constructed by the latter algorithm was more focused due to promoting in the search highly preferred solutions over the non-dominated ones. This implied stronger evolutionary pressure and faster convergence. The distributions reveal that DCEMOLα and CDEMOLα did not manage to focus the search around the DM’s most preferred region in the objective space so well as DCEMOASF and CDEMOASF (see Figures 7c, 7d, 7g, and 7h). To explain this, we investigated the search properties of the underlying cone contraction methods. For this purpose, we assumed a fixed set of 3 pairs of solutions to be compared by the DM with wDM = [0.5, 0.5]. Specifically, the DM judged that s1a = [1.20, 0.35]  s1b = [1.35, 0.20] after 50-th generation, s2a = [0.85, 1.10]  s2b = [0.70, 1.25] after 100-th generation, and s3a = [0.85, 0.85]  s3b = [0.95, 0.75] after 150-th generation. In Figure 8, we illustrate the expected spatial distributions for CDEMOASF and CDEMOL∞ interacting with thus simulated DM in 100 runs. We marked the boundaries of systematically contracted preference cones with a solid black line and the DM’s most preferred region with the black circle. Let us emphasize that both CON EASF and CON EL∞ were originally developed for exploiting a finite set of solutions. However, in evolutionary MOO, a population (i.e., a set of solutions) changes from one generation to another. In this perspective, the discriminatory power of CON EL∞ strongly depends on both the pairs of solutions compared by the DM as well as his/her answers. For illustrative purpose, let us consider Figure 8b. It reveals that only 3rd pairwise comparison differentiated between Pareto-optimal solutions, whereas the previous two DM’s judgements had no discriminatory power after the population had already converged towards the PF. In turn, CON EASF used each pairwise comparison to focus the search on the part of PF aligning with the DM’s preferences (see Figure 8a). Thus, the advantage of CDEMOASF and DCEMOASF over their counterparts incorporating CON EL∞ derives from 20

1.2

f2

f2

0.9

0.0

0.3

0.6

0.9

1.2

1.5

100 0.6

50

0.3

0

0.0 0.0

0.6

1.2

0.9

1.2

1.5

150

0.9 100 0.6

50

0.3

200

f2

f2 0.6

50

0.3 0.0 0.0

0

0.3

0.6

200

1.2

0.9

f2

f2 0.6

100 0.6

50

1.5

150

50

0.3

0

0.0 0.0

0.3

0.6

f1

1.2

1.5

0

(f) DCEMOASF (wDM = [0.7, 0.3]) 200

1.5 1.2

1.2 0.9

f2

f2 0.6

0.9

1.2

1.5

100 0.6

50

0.3

150

50

0.3 0.0 0.0

0

generation

100

generation

0.9

200

1.5

150

0.6

0.9

f1

(e) CDEMOASF (wDM = [0.7, 0.3])

0.3

generation

100

1.2

0

200

1.2

generation

0.9

0.9

1.5

1.5

150

0.3

1.2

(d) DCEMOL∞ (wDM = [0.5, 0.5])

(c) CDEMOL∞ (wDM = [0.5, 0.5]) 1.5

0.9

f1

f1

0.0 0.0

0

generation

100

0.6

1.5

1.2

generation

0.9

0.3

1.2

1.5

150

0.0 0.0

0.9

(b) DCEMOASF (wDM = [0.5, 0.5]) 200

1.5

0.6

0.3

f1

(a) CDEMOASF (wDM = [0.5, 0.5])

0.3

50

0.3

f1

0.0 0.0

150

generation

0.6

generation

100

0.0

1.2

150

0.9

200

1.5

200

1.5

0.3

0.6

0.9

1.2

1.5

0

f1

f1

(h) DCEMOL∞ (wDM = [0.7, 0.3])

(g) CDEMOL∞ (wDM = [0.7, 0.3])

Figure 7: Expected spatial distributions of evolutionary multiple objective optimization for different algorithms applied to DTLZ2 (M = 2, l = 50).

21

200

1.5

2-nd PC

1.2

2-n d

f2

f2

PC

1.2

1.5

1.0 1.0 0.0

0

0.0

f1

0.3

0.6

(a) CDEMOASF

1.5 1.5

0.0 0.0 1.5 1.5 1.0

0.5 0.5

0.0 0.0 1.5 1.5 1.0

0.0 0.0 1.5 1.5 1.0

1.0 1.0 1.01.5 1.5 f2 0.5 f2 0.5 0.0 0.0 0.0 0.00.5 0.5 f1 f1

1.5 1.5 1.5 1.5

1.5 1.5 1.5 1.5

1.5 1.5

1.5 1.5

ParetoPareto front front solution solution

1.0 1.0

f3

f3

f3

f3

0.5 0.5

0.0 0.0 0.0 0.0 1.5 1.0 1.5 1.0 1.5 1.0 1.5 1.0 1.5 1.01.5 1.5 0.5 1.0 1.0 1.01.5 0.5 0.5 f f2 0.5 0.0 0.00.5 f 0.0 0.00.5 0.5 f2 2 0.0 0.0 f2 0.5 0.0 0.0 f1 f1 1 f1

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.5 1.0 1.5 1.5 1.5 1.0 1.5 1.5 1.0 1.0 1.5 1.5 1.0 1.01.5 1.0 1.01.5 1.5 1.5 0.5 0.5 0.5 0.5 1.5 1.0 1.0 1.5 f2 1.0 f2 0.5 1.0 0.5 0.0 0.0 0.0 0.0 f f2 0.5 1.0 f 1.0 0.5 0.0 0.0 0.0 0.5 0.5 f f11.0 f11.01.5 1.5 0.0 2 0.5 1 0.5 0.5 0.5 f2 f2 0.0 0.0 f2 0.5 f2 0.5 0.0 0.0 0.0 0.0 0.0 0.0 f1 f11 f1 f1

(f) g = 50; w2DM

Pareto front

(d) g = 200; w1DM

ParetoPareto front front solution solution

0.5 0.5

0.5 0.5 0.5 0.5

(e) g = 20; w2DM

ParetoPareto front front ParetoPareto front front solution solution solution solution

0.5 0.5 0.5 0.5

1.0 1.0

f3f3

f3f3

f3f3

f3f3

0.5 0.5 0.5 0.5

1.5 1.5

1.0 1.0 1.0 1.0

(c) g = 100; w1DM

ParetoPareto front front ParetoPareto front front solution solution solution solution

1.0 1.0 1.0 1.0

1.0

0.0 0.0 0.0 0.0 0.0 0.01.5 0.0 0.01.5 1.5 1.0 1.0 1.5 1.5 1.0 1.5 1.0 1.01.5 1.5 1.5 0.5 1.0 1.01.5 0.51.0 0.5 0.5 1.5 1.5 1.5 1.5 1.5 1.0 1.0 f0.5 f2 0.5 0.5 0.0 0.0 0.00.5 0.5 ff2 1.0 ff2 1.0 0.0 2 0.5 1.0 0.0 0.0 0.0 0.0 0.5 1 0.5 ff1 f1.0 ff11.0 f11.0 1.5 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.5 0.5 2 f2 2 f2 1 f1 1 f1

(b) g = 50; w1DM

ParetoPareto front front ParetoPareto front front solution solution solution solution

1.0 1.0 1.0 1.0

1.0 1.0 1.0 1.0

1.5 1.5 1.5 1.5

0.5 0.5 0.5 0.5

1.0 1.0 1.01.5 1.5 f2 0.5 f2 0.5 0.0 0.0 0.0 0.00.5 0.5 f1 f1

(a) g = 20; w1DM

ParetoPareto front front ParetoPareto front front solution solution solution solution

1.0

1.0 0.5 0.5 0.5 f2 50) f2 0.5 0.0 0.0 0.0 2, l = with the DM’s f1 f1 objective

f3f3

0.5 0.5

0

f3 f3

f3 f3

f3 f3

f3f3

1.0 1.0

1.5 1.5 1.5 1.5

f3 f3

1.0 1.0

ParetoPareto front front solution solution

1.5

(b) CDEMOLα

0.0 0.0 1.5 1.5 1.0

Figure 8: Expected spatial distributions for CDEMOASF and weight vector wDM = [0.5, 0.5], P C denotes a pairwise comparison provided by the DM. ParetoPareto front front solution solution

1.2

50

0.5 0.5

1.5 1.5 1.0 0.5 0.5 0.51.0 1.0 f2 0.5 0.0 0.0to CDEMOL∞f2 applied 0.0DTLZ2 f1 f1 (M =

1.5 1.5

ParetoPareto front front solution solution

1.0 1.0

0.9

f1

0.5 0.5

1.5 1.5

1-st PC

0.9

front front Pareto 0.3Pareto solution solution

f3 f3

0.6

1.5 1.5

3-rd PC

50

C 1-st P

100 0.6

generation

PC

150

0.9

f3 f3

0.3

100

generation

rd 3-

0.3 0.0 0.0

1.2

150

0.9 0.6

200

1.5

(g) g = 100; w2DM

(h) g = 200; w2DM

Pareto front

Pareto front Pareto front Figure 9: solution Solutions constructed by CDEMO solution ASF after 20-th, 50-th, 100-th, and 200-th generation for DTLZ2 (M = 3) with the DM’s solution solution 1.5 1.5 1.5 1.5 [1/3, 1/3, 1/3] (a, b, c, d) and w2DM = [0.5, 0.3, 0.2] (e, f, g, h). objective weight vector set to w1DM =

1.0 1.0

1.0 1.0

f3

f3

f3

f3

the0.5ability to better discriminate 0.5 0.5 between the solutions, in particular, when the algorithms already converged close 0.5 to the PF. 0.0 0.0 0.0 0.0 1.5 1.0 1.5 1.0 1.5above 1.5 1.0 two-objective 1.5 1.5 1.5 problems. 1.5 All illustrations concerned In Figure 9, we visualize populations generated by 1.0 1.0 0.5 1.0 0.5 0.5 0.5 0.51.0 1.0 f2 f2 0.5 f2 0.5 0.0 0.0 0.0 0.0 f2 0.5 0.0 0.0 0.0 0.0 f1 f1 f1 f1 CDEMOASF applied to DTLZ2 with M = 3 objectives. Specifically, we illustrated solutions generated after 20th, 50-th, 100-th, and 200-th generations. Furthermore, the method was run for two artificial DMs using wDM = [1/3, 1/3, 1/3] and wDM = [0.5, 0.3, 0.2]. The provided illustrations reveal how the search space was systematically reduced, focusing the population on the PF’s region being relevant to the DM. 5.3. Comparison with selected interactive evolutionary algorithms for multiple objective optimization In this section, we compare CDEMOASF , i.e., the best performing variant of the proposed method, with the selected state-of-the-art interactive EMOAs. Even though all accounted methods are based on holistic judgements in form of pairwise comparisons, they differ in terms of the assumed preference model and accounting for the robustness concerns. In particular, ECC-MRW [35] incorporates ASF as the preference model and the non-dominated fronts as the primary sorting criterion (hence, being similar to DCEMO). However, it ranks the solutions contained in the same non-dominance front using a single arbitrarily selected compatible instance of ASF, thus failing to account 22

Table 5: The primary and secondary sorting criteria, robustness concern, and preference model incorporated by different EMOAs guided by interactively provided DM’s pairwise comparisons. Algorithm CDEMOASF ECC-MRW [35] NEMO-II [7] NEMO-0 [7] P&K [51]

primary-sort

secondary-sort

RAN KASF non-dominated fronts non-dominated fronts the most discriminative ASF fronts of potentially optimal solutions crowding-distance non-dominated fronts the most discriminative AVF the most discriminative ALVF

Robustness analysis yes no yes no no

Preference model ASF ASF AVF AVF ALVF

Table 6: The average ARSD and respective ranks R attained by different methods for the populations generated in the last generation for the DTLZ2(C) and WFG1 test problems with M = 2, 3, 4, and 5 objectives. The DM’s decision policy was modelled with L1 and L∞ (ρ = 3 and 2 for, respectively, the DM modelled with L1 and L∞ ). CDEMOASF Mean R StD DTLZ2C

WFG1

M 2 3 4 5 2 3 4 5

Average R

DTLZ2

WFG1 Average R

2 3 4 5 2 3 4 5

1.64 13.58 48.69 95.98 3.58 54.60 77.65 86.81 Mean 3.45

3.92 20.89 50.40 91.33 3.82 24.95 28.02 29.35 StD 1.00

3.22 2.82 2.83 2.77 2.14 4.60 4.62 4.56

0.44 2.82 4.27 5.03 0.57 2.23 3.57 3.96 Mean 1.93

1.03 3.34 3.22 3.71 0.48 1.92 2.19 2.34 StD 0.51

1.35 1.86 2.10 2.27 1.04 1.94 2.27 2.61

ECC-MRW [35] NEMO-II [7] NEMO-0 [7] Mean R Mean R Mean StD StD StD Simulated Decision Maker’s model = L1 10.00 0.03 3.96 29.73 4.21 0.09 1.25 10.99 73.72 3.23 39.16 98.54 4.02 6.43 1.52 62.69 148.34 152.27 21.48 19.17 2.01 101.85 130.80 4.06 191.47 161.59 56.95 37.86 2.06 144.59 165.61 3.67 7.32 2.40 4.16 4.63 4.40 1.24 1.75 2.04 30.03 4.24 16.98 56.97 3.27 7.40 1.55 28.26 28.29 8.62 37.82 33.49 2.90 8.64 1.92 56.70 31.92 19.71 14.17 2.50 46.42 47.99 2.72 100.06 Mean Mean Mean StD StD StD 3.66 1.82 2.99 0.63 0.39 0.26 Simulated Decision Maker’s model = L∞ 6.41 32.07 8.56 8.70 2.57 14.04 4.45 11.36 6.90 29.88 8.99 8.93 2.47 9.17 4.60 10.43 8.24 56.47 5.74 9.29 2.75 15.38 4.98 8.09 8.10 69.65 6.36 9.81 2.50 14.22 5.00 7.05 5.61 27.60 5.17 6.50 2.86 8.12 4.68 7.23 5.71 20.73 6.08 9.07 2.59 3.86 4.78 8.83 5.27 15.94 7.08 7.31 2.45 3.17 4.69 9.65 3.95 12.95 5.55 4.82 2.25 2.78 4.68 5.67 Mean Mean Mean StD StD StD 2.56 4.73 2.63 0.19 0.19 0.18

P&K [51] StD

R

Mean

R

3.32 3.17 3.08 3.00 3.22 2.68 2.82 2.62

2.60 51.74 79.79 172.42 4.55 29.63 34.15 34.07 Mean 3.09

5.80 95.98 92.53 163.69 2.33 65.10 62.16 58.04 StD 0.36

3.00 3.47 3.02 3.50 3.49 2.90 2.74 2.60

2.74 2.70 2.40 2.30 2.66 2.64 2.80 2.79

26.74 17.64 9.97 13.23 17.22 9.78 7.43 7.38 Mean 3.15

28.41 21.51 13.47 14.33 20.25 13.79 9.99 11.62 StD 0.47

3.89 3.37 2.77 2.93 3.76 3.05 2.79 2.67

for the robustness preoccupation. Furthermore, we consider a pair of algorithms using an AVF [40] with piece-wise linear marginal value functions as a preference model, NEMO-0 and NEMO-II [7]. The former is a counterpart of ECC-MRW using a single compatible instance of AVF to rank the solutions which are incomparable in terms of dominance, whereas the latter replaces the non-dominated fronts with the fronts of potential optimality with respect to any compatible instance of AVF. In this regard, NEMO-II exploits a set of compatible instances of AVF to promote the solutions which are ranked first for at least one compatible model instance. Lastly, we consider a method proposed in [51], denoted by P &K. This algorithm assumes that the DM’s preferences are modelled with an Additive Linear Value Function (ALVF) and selects – similarly to NEMO-0 and ECC-MRW – a single representative instance of the model to rank solutions in the population. However, unlike the remaining methods, it uses only one- rather than two-tier ranking for selection. For the summary of properties of these algorithms, see Table 5. We evaluated the performance of the algorithms on DTLZ2 and WFG1 with M = 2, . . . , 5 objectives. For each problem, we considered the artificial DMs using either L1 or L∞ . When it comes to the parametrization of NEMO, we set the number of characteristic points of the marginal value functions to 2 or 5 when dealing with the DM’s decision policy simulated with L1 and L∞ , respectively. In the former case, this guarantees a full consistency between the preference model and the DM’s pairwise comparisons. In the latter case, the selection of 5 characteristic points was motivated by the results presented in [32] so that to ensure a balance between the model’s expressiveness and the robustness of provided results. Since the PF of DTLZ2 is concave, we used its convex counterpart [19], denoted by DTLZ2C, for the DM modelled with L1 In Table 6, we report ARSD for the populations constructed in the final generation by all four methods. The results indicate that CDEMOASF was superior to the remaining methods for all benchmarks involving the DM modelled with L∞ . Specifically, CDEMOASF attained the best average rank of 1.93, followed by ECC-MRW with R = 2.56. However, in case the simulated DM compared pairs of solutions in line with some AVF, both methods incorporating ASF as a preference model were outperformed by NEMO-II, NEMO-0, and P&K. In particular, NEMO-II, NEMO-0, and P&K attained the average ranks of 1.88, 2.99, and 3.09, respectively. The above results indicate that the performance of a preference-based EMOA can be vastly improved in two ways. 23

1.2

1.2

ECC-MRW g=50

1.2

0.9

f2

0.6

ECC-MRW g=200 CDEMO ASF g=20 ECC-MRW 0.6 CDEMO ASF g=50 g=200 ECC-MRW g=100

f2

ECC-MRW g=200

CDEMOASF g=20

CDEMOASF g=100 0.6 ECC-MRW g=100 CDEMOASF g=200 g=20 ECC-MRW g=50 0.3 ECC-MRW g=100 ECC-MRW g=200

ECC-MRW g=100

0.3

Optimum

0.9

Optimum f2

Optimum

0.9

0.3 ECC-MRW

0.0

0.0 0.0

0.3

0.6

f1

(a) wDM = [0.5, 0.5]

0.9

1.2

0.0 0.0 0.0

0.3 0.3

0.6

0.9

1.2

f0.6 1

0.9

1.2

CDEMOASFg=20 g=50 CDEMO ASF CDEMO ASF g=100 CDEMO ASF g=50 CDEMO ASF g=200 CDEMO ASF g=100 ECC-MRW g=20 CDEMO ASF g=200 ECC-MRWg=20 g=50 ECC-MRW ECC-MRWg=50 g=100 ECC-MRW ECC-MRWg=100 g=200 ECC-MRW ECC-MRW g=200

f1

(b) wDM = [0.7, 0.3]

Figure 10: Populations constructed by CDEMO and ECC-MRW [35] after 20-th, 50-th, 100-th, and 200-th generations for DTLZ2 with M = 2 objectives. The DM’s decision policy was simulated using L1 with wDM .

Firstly, a model employed by the method should align with the DM’s value system. As indicated in Table 6, the methods whose model was consistent with the DM’s preference function were ranked more advantageous. Secondly, accounting for the robustness preoccupation improves the quality of constructed solutions. The latter was confirmed with CDEMOASF being superior to ECC-MRW as well as NEMO-II outperforming NEMO-0 and P&K. The sole exception was the advantage of NEMO-0 and P&K over NEMO-II for the DM incorporating L∞ . This derives from using – by NEMO-II – a flexible AVF with 5 characteristic points, which may indicate many solutions as potentially optimal and thus deteriorate the pace of the evolutionary search. To illustrate a difference in the performances of methods either incorporating or not accounting for the robustness concern, we applied CDEMOASF and ECC-MRW to DTLZ2 with M = 2 objectives, and assumed a DM incorporating L∞ with wDM = [0.5, 0.5] and wDM = [0.7, 0.3]. In Figure 10, we illustrate the solutions constructed by both methods after 20-th, 50-th, 100-th, and 200-th generations. One may observe that CDEMOASF was systematically focusing the solutions in the relevant region of the PF. In turn, ECC-MRW was concentrating the solutions in different regions during the evolutionary search, and finally did not manage to reach close enough to the optimum. To explain this, let us remind that ECC-MRW uses a single ASF instance to differentiate between the solutions which are incomparable in terms of dominance. Such an instance is selected arbitrarily and can simply change with each new pairwise comparison provided by the DM, thus orienting the evolutionary search towards different, potentially less relevant Pareto-optimal solutions. Yet, ECC-MRW constructs a more crowded population when compared to CDEMOASF . This is especially true for problems involving many objectives, where the reduction of the search space imposed by CDEMOASF is less efficient. It explains why ECC-MRW could attain as satisfactory results as CDEMOASF according to ARSD for WFG1 with M = 5 and the DM incorporating L∞ (see Table 6). 5.4. Interaction patterns In this section, we evaluate different interaction patterns for deciding when the DM should be questioned to provide pairwise comparisons of solutions from the current population. In particular, we compared the performance of CDEMOASF employing different (both static and dynamic) interaction patterns P AT T ERN = {ST ART -ESI ST R , ESI P CS , P P -ESI P OS , RDS th } with CDEMOASF incorporating a default pattern ESI 1 with equally distributed questions starting from the first generation. We evaluated all algorithms on DTLZ2 and WFG1 involving from M = 2 to 5 objectives, and reported average ARSD for (a) the populations generated in the last generation in a tabular form, and (b) the populations constructed throughout the evolutionary search in the form of convergence plots. Furthermore, we illustrate how the number of collected pairwise comparisons |H| changed throughout the optimization. For brevity, we refer to a particular variant of CDEMOASF by recalling solely a name of the employed interaction pattern. 24

Table 7: Average ARSD for CDEMO using ST ART -ESI ST R parametrized with different values of ST R, applied to DTLZ2 and WFG1 with M = 2, . . . , 5 (ρ = 2).

ST R

Mean

M =2 StD

R

ESI 1 (EI = 20)

(1) 2 3 4 5 6 7 8 9 10

0.44 0.66 0.72 1.68 1.04 2.52 2.88 3.73 5.99 9.66

1.03 1.51 1.38 3.51 1.83 4.21 3.98 4.98 6.83 9.50

2.67 3.31 3.84 4.52 4.56 5.60 6.39 7.08 8.03 9.00

(1) 2 3 4 5 6 7 8 9 10

0.57 0.67 0.83 1.00 1.47 1.69 2.48 3.77 6.15 10.29

0.48 0.80 0.72 0.94 1.88 1.55 2.15 3.48 5.28 8.56

2.90 2.77 3.61 3.78 4.53 5.42 6.65 7.38 8.40 9.56

ESI 1 (EI = 150)

average relative score difference

ST ART -ESI ST R

10−1

0

500

1000

generation

M =3 Mean StD DTLZ2 2.82 3.34 3.35 3.63 3.51 4.84 3.80 4.53 3.67 3.79 4.08 3.96 4.93 4.69 5.48 4.93 7.18 6.00 7.20 5.14 WFG1 2.23 1.92 2.53 2.07 2.63 2.07 3.10 2.47 3.47 2.61 4.28 3.00 4.96 3.98 5.72 4.20 7.03 5.02 8.91 5.59

R

Mean

M =4 StD

R

Mean

M =5 StD

R

3.36 4.28 4.35 4.59 4.81 5.25 6.00 6.65 7.54 8.17

4.27 4.72 4.98 4.59 5.06 5.01 5.45 6.79 6.08 8.22

3.22 4.14 4.31 3.70 3.83 3.76 3.90 4.67 3.99 5.24

4.26 4.64 4.76 4.46 4.99 5.19 5.45 6.72 6.62 7.91

5.03 5.52 5.49 5.41 5.39 6.74 6.06 6.86 7.96 9.21

3.71 3.58 3.59 3.93 3.89 4.02 3.86 3.68 4.13 5.08

4.35 4.67 4.61 4.72 4.64 5.95 5.22 6.13 7.12 7.59

2.77 3.02 3.28 3.98 4.73 5.72 6.61 7.11 8.30 9.48

3.57 3.38 3.45 3.97 4.44 4.66 5.11 6.07 6.19 7.05

2.19 2.05 2.04 2.52 2.75 3.00 3.27 3.91 4.22 4.49

3.63 3.27 3.48 4.35 4.85 5.44 6.22 7.26 7.71 8.79

3.96 4.02 4.12 4.22 4.47 4.79 5.14 5.46 5.75 6.07

2.34 2.21 2.42 2.31 2.77 2.86 3.10 3.64 3.47 3.71

3.82 3.90 4.03 4.43 4.79 5.40 6.44 6.61 7.43 8.15

10

8

10−1

10−2

history size |H|

ST ART -ESI ST R

average relative score difference

Pattern

ESI 1 (EI = 150) 2

-ESI ESI 1ST (EIART = 150) 3 ST-ESI ART(ST -ESI ST ART R = 2) ST ART R = 43) ST-ESI ART(ST -ESI ST ART R = 54) ST-ESI ART(ST -ESI ST ART -ESI (ST R = 65) ST ART -ESI ST ART -ESI (ST R = 76) ST-ESI ART(ST -ESI ST ART R = 87) ST-ESI ART(ST -ESI ST ART R = 8) ST ART R = 99) ST-ESI ART(ST -ESI ST ART R = 10 10) ST-ESI ART(ST -ESI

ESI 1 (EI = 150) ST ART -ESI 2 ST ART -ESI 3 4 4 ST ART -ESI 5 ST ART -ESI ST ART -ESI 6 ST ART -ESI 7 2 ST ART -ESI 8 ST ART -ESI 9 ST ART -ESI 10 6

0

0

(a) Average Relative Score Difference

1500

0

500

500

1000 generation

1000

1500

1500

generation (b) History size |H|

Figure 11: Average ARSD and history size (the number of so far collected pairwise comparisons) for CDEMOASF with different interaction patterns, applied to WFG1 with M = 3 objectives.

Focus on the ST ART -ESI ST R interaction pattern. Let us first focus on the comparison of results attained with the interaction pattern ST ART -ESI ST R , which elicits ST R pairs of solutions in the first generation, whereas distributing the remaining 10 − ST R interactions concerning a single pair of solutions equally over the optimization process (see Section 3.2). In Table 7, we provide the results attained with ST ART -ESI ST R for ST R = 2, . . . , 10 (note that ST ART -ESI 1 is equivalent to ESI 1 ). The results indicate that application of ESI 1 was superior to the remaining patterns. Moreover, comparable results were attained with ST ART -ESI ST R for relatively small values of ST R. The less objectives were involved, the more evident the advantage of ESI 1 over other patterns was. In particular, when considering DTLZ2 with M = 2, . . . , 5 objectives there was no significant difference in the performance of ESI 1 and ST ART -ESI ST R for the following ST R values: for M = 2 → ST R = 2, M = 3 → ST R ∈ {2, 3}, M = 4 → ST R ∈ {2, 3, 4}, and M = 5 → ST R ∈ {2, 3, 4, 5}. This was due to imposing a stronger evolutionary pressure by ST ART -ESI ST R with ST R > 1 in the beginning of optimization, which helped converging faster to the PF when dealing with more objectives. For illustrative purpose, let us consider Figure 11 reflecting the performance of different patterns applied to WFG1 with M = 3. One may observe that all variants of ST ART -ESI ST R were performing better than ESI 1 until around 400-th generation. This was due to a stronger evolutionary pressure imposed by multiple preference elicitations performed in the first generation (see the history size |H| depicted in Figure 11b). Yet, let us remind that the reference solutions for such pairwise elicitation questions were selected from the same population, which diminished 25

Pattern

P CS

Mean

M =2 StD

R

ESI P CS

1 (EI = 20) 2 (EI = 40) 5 (EI = 100)

0.44 1.44 5.23

1.03 2.87 6.40

1.27 1.96 2.77

ESI P CS

1 (EI = 150) 2 (EI = 300) 5 (EI = 750)

0.57 0.95 3.56

0.48 0.91 3.18

1.35 1.77 2.88

500

Mean

M =4 StD

R

Mean

M =5 StD

R

1.58 2.00 2.42

4.27 5.16 6.28

3.22 3.61 4.69

1.67 1.96 2.37

5.03 4.94 6.69

3.71 3.39 5.01

1.91 1.81 2.28

1.45 1.83 2.72

3.57 3.79 4.79

2.19 2.36 2.75

1.63 1.79 2.58

3.96 3.94 4.58

2.34 2.25 2.53

1.82 1.77 2.41

6

2

0

1000

1500

ESI 1 (EI = 150) ESI 2 (EI = 300) ESI 5 (EI = 750)

4

0

0

R

8

10−2 10−2

M =3 Mean StD DTLZ2 2.82 3.34 3.65 4.82 4.72 4.17 WFG1 2.23 1.92 2.71 1.84 5.06 3.33

10

10−1

history size |H|

10−1

average relative score difference

average relative score difference

Table 8: Average ARSD for CDEMO using ESI P CS parametrized with different values of P CS, applied to DTLZ2 and WFG1 with M = 2, . . . , 5 (ρ = 2).

ESI 1 (EI = 150) ESI 2 (EI 500= 300) ESI 5 (EI = 750)

1500

generation

0

500

1000

ESI 1 (EI = 150) ESI 2 (EI = 300) ESI 5 (EI = 750)

1500

generation

generation

(a) Average Relative Score Difference

1000

(b) History size |H|

Figure 12: Average ARSD and history size (the number of so far collected pairwise comparisons) for CDEMOASF with different interaction patterns, applied to WFG1 with M = 3 objectives.

the pace of preference cone contraction and its impact on the algorithm’s performance. Thus, with greater ST R, the average quality of solutions constructed in the final generation was deteriorated. Specifically, ST ART -ESI ST R with ST R = 10 attained the worst average ranks ranging from 7.59 (DTLZ2 with M = 5) to 9.56 (WFG1 with M = 2). Note that this variant can be perceived as an a priori method, since all 10 pairwise comparisons are provided in the first generation, which prevents further interaction with the DM. In turn, ESI 1 ensured a better balance between the population’s convergence and the delay between subsequent interactions. This allowed to impose a more efficient search space reduction, and consequently to generate a more crowded population. Specifically, ESI 1 attained the best average ranks ranging from 2.67 (DTLZ2 with M = 2) to 4.35 (DTLZ2 with M = 5). Focus on the ESI P CS interaction pattern. Let us now consider the patterns requiring the DM to provide P CS pairwise comparisons in each interaction. In Table 8, we provide the results attained with ESI P CS for P CS = 1 (i.e., a default interaction pattern), 2, and 5. Similarly to ST ART -ESI ST R , ESI 1 outperformed the remaining patterns for the vast majority of benchmark problems. In particular, ESI 1 attained the average ranks ranging from 1.27 (DTLZ2 with M = 2) to 1.91 (DTLZ2 with M = 5). The exceptions in this regard can be observed for DTLZ2 with M = 5 and WFG1 with M = 4 and M = 5. For these problems, there were no significant differences between ESI 1 and ESI 2 . On the one hand, smaller P CS values imply higher frequency of interactions, which results in the systematic convergence of the population by means of progressively renewed pressure. On the other hand, higher P CS makes a single interaction with the DM more informative as (s)he is expected to compare more pairs of solutions. For example, ESI 1 asked the DM to compare a single pair of solutions every 150 generations, while ESI 5 selected 5 pairs to be compared in the first and 750-th generations (see history size |H| in Figure 12b). However, as previously noted, the search space reduction may be inefficient when multiple pairwise comparisons are provided in the same generation. This is confirmed with the analysis of, e.g., average ARSD throughout the evolutionary search for ESI 1 and ESI 5 applied to WFG1 with M = 3 (see Figure 12a). Although ESI 5 collected 5 pairwise comparisons in the first generation, it was outperformed by ESI 1 already after 300 generations. It means that ESI 1 required only 2 interactions concerning a single pair of solutions to attain similar results to ESI 5 . However, as indicated in Table 8, 26

Table 9: Average ARSD for CDEMO using P P -ESI P OS parametrized with different values of P OS, applied to DTLZ2 and WFG1 with M = 2, . . . , 5 (ρ = 2).

Pattern

P OS

Mean

M =2 StD

R

ESI 1 (EI = 20)

(0) 10 20 30 40 50 60 70 80 90 100

0.44 0.69 0.58 0.45 0.38 0.49 0.63 0.56 0.47 0.43 0.43

1.03 1.73 1.29 0.94 0.75 1.12 1.22 1.70 0.87 0.82 0.75

5.57 6.13 6.57 5.80 5.62 5.36 6.80 6.22 6.15 5.82 5.96

(0) 75 150 225 300 375 450 525 600 675 750

0.57 0.55 0.54 0.59 0.57 0.62 0.72 0.71 0.71 0.85 0.93

0.48 0.49 0.44 0.66 0.36 0.41 0.54 0.49 0.49 0.52 0.61

5.11 4.47 4.71 5.07 5.44 5.72 6.30 6.37 6.83 7.95 8.03

P P -ESI P OS

ESI 1 (EI = 150)

P P -ESI P OS

M =3 Mean StD DTLZ2 2.82 3.34 3.23 3.64 3.02 2.98 3.12 3.51 3.05 3.69 3.23 3.59 3.45 3.89 3.07 3.74 3.04 3.97 3.34 3.71 2.96 3.07 WFG1 2.23 1.92 2.20 1.63 2.30 1.82 2.25 1.66 2.20 1.74 2.28 1.63 2.17 1.40 2.34 1.65 2.28 1.60 2.47 1.52 2.54 1.79

Mean

M =4 StD

R

Mean

M =5 StD

R

5.29 6.04 6.23 6.04 5.98 6.22 6.52 5.88 5.59 6.49 5.72

4.27 4.31 4.67 4.84 4.19 4.91 4.67 4.34 4.75 4.63 4.54

3.22 3.39 3.80 4.18 3.34 3.58 3.72 3.61 3.69 3.56 3.87

5.70 5.65 6.14 5.89 5.60 6.65 6.19 5.83 6.19 6.24 5.92

5.03 4.58 5.60 4.66 5.16 5.69 4.88 5.97 6.02 5.65 6.73

3.71 3.21 4.35 3.69 3.68 4.23 3.67 4.33 5.09 4.53 6.26

5.67 5.45 6.09 5.17 6.25 6.25 5.41 6.54 6.16 6.42 6.59

5.41 5.57 6.00 5.76 5.55 5.77 5.87 6.41 6.14 6.72 6.80

3.57 3.51 3.22 3.26 3.11 3.31 3.15 3.05 3.23 3.29 3.04

2.19 2.12 2.19 1.96 1.95 1.81 1.94 1.81 1.85 2.00 1.76

6.97 6.82 5.61 5.91 5.64 6.49 5.60 5.45 6.01 5.92 5.58

3.96 3.75 3.65 3.32 3.44 3.43 3.23 3.13 3.13 2.97 3.24

2.34 2.34 2.06 1.51 1.71 1.77 1.59 1.57 1.69 1.53 2.07

7.79 6.82 6.64 6.12 6.35 6.25 5.71 5.23 5.05 4.72 5.32

2 × 10−1

2 × 10−1

10

8

history size |H|

average relative score difference

average relative score difference

R

10−1

6 × 10−2 4 × 10−2 3 × 10−2

1 ESI (EI = 150) 10−1 P6P -ESI 75 P P -ESI 150 P P -ESI 225 300 −2 6P ×4P10-ESI P P -ESI 375 P P -ESI 450 P2P -ESI 525 −2 600 4P ×P10-ESI P P -ESI 675 750 −2 3P ×P10-ESI

ESI 1ESI (EI 1=(EI 150)= 150) 75 P P -ESI P P 150 -ESI 75 P P -ESI P P -ESI 150 P P -ESI 225 225 P P 300 -ESI P P -ESI P P 375 -ESI 300 P P -ESI 450 P P -ESI 375 P P -ESI P P 525 -ESI 450 P P -ESI P P 600 -ESI 525 P P -ESI 675 P P -ESI P P -ESI 600 750 P P -ESI P P -ESI 675

0

0

500

1000

generation

1500

0

0

500

1000

generation 500

P P -ESI 750

1500

1000

1500

generation

(b) History size H

(a) Average relative score difference

Figure 13: Average ARSD and history size (the number of so far collected pairwise comparisons) for CDEMOASF with different interaction patterns, applied to WFG1 with M = 5 objectives.

less frequent but more informative interactions may yield satisfactory results when many objectives are involved. It suggests that there is a trade-off between the frequency and intensity of interactions if the underlying problem is more challenging to solve. Focus on the P P -ESI P OS interaction pattern. Let us now examine a pattern postponing the preference elicitation until P OS-th generation and distributing all (in our case – 10) interactions with the DM over the remaining GEN − P OS generations. We used P OS ∈ {5%·GEN, 10%·GEN, . . . , 50%·GEN }, thus making the value of P OS dependent on the total number of generations GEN , which is different for DTLZ2 (GEN = 200) and WFG1 (GEN = 1500). We provide the results of the experimental comparison in Table 9 (note that P P -ESI 0 is equivalent to ESI 1 ). The results indicate no significant difference between application of P P -ESI P OS with P OS ∈ {0, 10, . . . , 100} to DTLZ2 for various numbers of objectives M . Let us remind that DTLZ2 is relative easy to solve, having no, e.g., local optima or flat regions. As a result, (a) new non-dominated solutions can be generated relatively fast and (b) it is easy to generate a variety of solutions representing different trade-offs between objectives. Consequently, the population can quickly converge to a region on the PF indicated by the preference cone. Thus, there is a great flexibility in adjusting P P -ESI when solving DTLZ2. In turn, WFG1 has a biased fitness landscape [29]. Specifically, the closer the algorithm gets to the PF, the more difficult it is to generate an offspring yielding an improvement. In this regard, it may be advantageous to postpone 27

the preference elicitation until the population is already converged to the PF so as to impose a stronger evolutionary pressure at the stage when the exploration of the objective space gets more challenging. The results provided in Table 9 confirm this claim. Although ESI 1 performed satisfactorily for WFG1 with M = 2 and 3 objectives, it was significantly outperformed by P P -ESI P OS with P OS >> 0 when 4 or 5 objectives were involved. For example, P P -ESI 675 and ESI 1 attained the average ranks equal to, respectively, 4.72 and 7.79 for WFG1 with M = 5. In Figure 13a, we illustrate the performance for different patterns when applied to WFG1 with M = 5. Let us focus on P P -ESI P OS with P OS >> 0 (e.g., P OS = 675 or 750). These methods evolved the population by means of NSGA-II until P OS-th generation was reached. As evident in Figure 13a, such an exploration was efficient until around 500-th generation, but in the subsequent generations there was no significant improvement in ARSD. Then, 10 preference elicitations were equally distributed over the remaining generations (see history size |H| in Figure 13b). Consequently, the interactions were more frequent when the population was closer to the PF, which, in turn, provided a stronger evolutionary pressure at the stage when exploration of the objective space became more challenging. It helped to improve the final results. In particular, all variants of P P -ESI P OS outperformed ESI 1 after 1000 generations. Intriguingly, ESI 1 had no initial advantage over the remaining methods, despite collecting a single pairwise comparison already in the first generation. This is due to a population being focused in a relatively small region of the objective space at the initial stage of the optimization when approaching WFG1 [29]. This, in turn, makes a variety of trade-off solutions limited, diminishing the role of the evolutionary pressure imposed by the specification of DM’s preferences. Focus on the RDS th interaction pattern. Let us now examine the results attained by CDEMOASF using dynamic interaction pattern RDS th . In this study, we considered th ∈ {0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5}. Note that the greater th is, the more likely is the preference elicitation to be triggered. Furthermore, let us remind that we limited the total number of interactions to 10. The results of an experimental comparison are provided in Table 10. When it comes to DTLZ2, the results show that RDS th parametrized with th > 0.20 and th > 0.10 performed similarly to ESI 1 for all benchmarks with M = 2, 5 or M = 3, 4, respectively. For the remaining variants of the method, the lower th was incorporated, the worse was the average rank R according to ARSD. In particular, RDS 0.01 attained the worst average ranks ranging from 10.49 to 10.67 for different values of M . Low values of th impose rare preference elicitation. It is thus possible that the DM was questioned less than the maximal allowed number (10) in the course of the evolutionary search, which may explain a relatively worse performance of the method. Furthermore, RDS 0.01 performed poorly on WFG1, attaining the average ranks ranging from 7.59 to 9.97 for different M . However, in contrast to DTLZ2, RDS th with higher th values did not manage to attain satisfactory results. In particular, RDS 0.50 attained the average ranks between 6.51 and 8.54. This, in turn, was caused by performing the preference elicitation too frequently. As the number of preference elicitations was limited to 10, most interactions with the DM were simply performed at the early stage of optimization. However, as already discussed, it is reasonable to postpone the preference elicitation when solving WFG1. This explains why RDS th parametrized with intermediate thresholds (e.g., th = 0.03) attained the best results. Importantly, RDS th with th ∈ {0.02, 0.03, 0.04} performed better than ESI 1 for WFG1 with M = 4 and 5 objectives. On the one hand, such thresholds were restrictive enough to avoid redundant pairwise elicitation questions. On the other hand, they guaranteed that all 10 interactions were performed during the optimization. In Figure 14, we illustrate the performance on WFG1 with M = 5 as well as the average number of collected pairwise comparisons when using the dynamic interaction pattern with different thresholds th. Let us start with discussing the performance of RDS 0.01 . This variant performed, on average, 8.96 interactions in the course of the evolutionary search. Due to fewer interactions, the method did not manage to focus the population as satisfactorily as other variants did (see ARSD for different methods in Figure 14a). Furthermore, RDS th incorporating th = 0.02, 0.03, and 0.04 proved to perform better than ESI 1 (see Table 10). These variants managed to question the DM 10 times. Yet, the preference elicitation was less frequent at the initial stage of the optimization when compared to ESI 1 . For instance, RDS 0.02 performed the same number of interactions as ESI 1 after around 250 generations. This suggests that RDS th adapted to the fitness landscape and postponed the preference elicitation when it was easier to 28

Table 10: Average ARSD for CDEMO using RDS th parametrized with different values of th, applied to DTLZ2 and WFG1 with M = 2, . . . , 5 (ρ = 2).

Pattern

th

Mean

M =2 StD

R

ESI 1 (EI = 20)

− 0.01 0.02 0.03 0.04 0.05 0.10 0.20 0.30 0.40 0.50

0.44 30.93 23.43 13.60 9.35 5.42 0.91 0.55 0.48 0.55 0.56

1.03 18.82 19.01 13.27 10.01 6.26 2.77 1.47 1.05 1.78 1.21

3.19 10.49 9.90 8.68 8.19 7.13 3.98 3.55 3.72 3.57 3.60

− 0.01 0.02 0.03 0.04 0.05 0.10 0.20 0.30 0.40 0.50

0.57 1.26 0.68 0.61 0.64 0.67 0.90 1.06 1.25 1.39 1.47

0.48 1.16 0.49 0.75 0.47 0.55 1.01 1.05 1.33 1.41 1.37

4.05 7.59 5.00 3.66 4.34 4.57 5.84 6.90 7.75 7.76 8.54

RDS th

ESI 1 (EI = 150)

RDS th

average relative difference historyscore size |H|

average relative score difference

Mean

M =4 StD

R

Mean

M =5 StD

R

3.38 10.49 10.03 8.74 7.79 6.50 3.63 3.99 4.06 3.65 3.74

4.27 23.73 21.06 13.81 9.61 6.31 4.53 4.54 4.10 4.12 4.41

3.22 6.83 6.55 6.32 5.15 4.20 3.72 4.95 3.18 3.58 3.75

4.05 10.61 10.02 8.83 7.34 5.71 4.06 3.85 3.85 3.71 3.97

5.03 22.63 19.19 13.32 8.84 5.62 5.43 4.83 5.04 4.80 5.17

3.71 5.02 4.97 5.31 4.98 3.74 4.10 3.57 3.34 3.85 3.67

4.23 10.67 9.94 8.59 6.78 4.86 4.62 4.03 4.27 3.61 4.40

4.84 9.97 4.93 4.62 4.46 4.73 5.28 5.75 6.56 7.10 7.76

3.57 4.67 3.54 3.26 3.36 3.74 3.59 4.15 4.28 4.10 4.04

2.19 2.96 2.23 1.99 1.94 2.57 2.39 2.73 2.73 2.82 2.56

5.55 8.44 5.22 4.51 4.59 5.45 5.37 6.60 7.09 6.67 6.51

3.96 4.53 3.55 3.58 3.53 3.92 3.94 3.85 4.19 4.06 4.36

2.34 2.39 2.10 1.88 1.91 2.41 2.32 1.99 2.39 2.12 2.55

5.99 8.09 4.57 4.68 4.78 5.69 5.97 5.83 6.92 6.39 7.09

8

ESI 1 (EI = 150) ESI 1 (EI = 150) RDS 0.01 0.01 RDS 0.02 RDS 0.02 RDS 0.03 RDS 0.03 RDS 0.04 RDS 0.04 0.05 RDS RDS 0.10 0.05 RDS RDS 0.20 0.10 RDS RDS 0.30 0.20 RDS RDS 0.40 0.30 RDS RDS 0.50 0.40 RDS RDS RDS 0.50

ESI101 −1 (EI = 150) 6 0.01 RDS RDS 0.02 RDS 0.03 0.04 RDS 4 0.05 −2 RDS 6 × 100.10 RDS 0.20 RDS 2 0.30 RDS 0.40 RDS 4 × 10−2 RDS 0.50

−1

6 × 10−2 4 × 10−2 3 × 10−2

R

210× 10−1

2 × 10−1

10

M =3 StD DTLZ2 2.82 3.34 27.14 10.67 23.71 11.97 15.64 9.48 10.90 7.94 7.03 5.94 3.01 3.36 3.46 3.96 3.45 4.25 2.70 3.02 3.06 3.35 WFG1 2.23 1.92 5.31 2.99 2.27 1.81 2.13 1.74 2.20 2.06 2.25 1.97 2.62 2.53 2.76 2.36 3.27 3.22 3.27 2.54 3.74 3.18

Mean

0

500

1000

generation

1500

0 3 ×010−2

500

0

500 generation

1000 1000

1500 1500

generation

(a) Average relative score difference

(b) History size |H|

Figure 14: Average ARSD and history size (the number of so far collected pairwise comparisons) for CDEMOASF with different interaction patterns, applied to WFG1 with M = 5 objectives.

generate offspring yielding an improvement. In turn, the preference elicitation was more frequent thereafter as the problem got more challenging, and, e.g., RDS 0.02 questioned the DM for the 10-th time around 900-th generation. As we previously discussed, postponing the preference elicitation may improve the results when the fitness landscape is particularly challenging near the PF. Finally, Figure 14b indicates that the greater th was, the sooner all 10 pieces of preferences were collected by RDS th . However, RDS th with th ≥ 0.05 performed worse then RDS th incorporating th = 0.02, 0.03, and 0.04. Furthermore, RDS th with th ≥ 0.30 was outperformed even by ESI 1 (see Table 10). This was due to performing preference elicitation when it was relatively easy to construct better offspring, instead of supporting the search when the problem got more challenging. 5.5. Impact of the evolutionary base on the performance of the proposed method So far we have focussed on studying the impact of different preference-based sorting procedures incorporated within a relatively simple generational front-based evolutionary framework modifying NSGA-II [15]. In this section, we investigate whether the use of a particular evolutionary base has a significant influence on the performance of preference-based EMOAs. For this purpose, inspired by IEMO/I [65], we propose an indicator-based counterpart of the interactive evolutionary hybrid using preference cones. Specifically, we replaced the non-dominated sorting – originally used in CDEMO – with a hypervolume measure [76], obtaining in this way a new variant of the algorithm, called CDEMO-HV. Moreover, instead of using a generational scheme, CDEMO-HV works in a steady-state envi29

Table 11: Characteristics of CDEMOASF and CDEMO-HVASF . Algorithm CDEMOASF CDEMO-HVASF

Evolutionary base front-based/generational indicator-based/steady-state

primary-sort RAN KM ODEL RAN KM ODEL

secondary-sort non-dominated fronts hypervolume

Table 12: Average ARSD attained by the populations constructed by CDEMOASF and CDEMO-HVASF applied to DTZL2 and WFG1 with M = 2, 3, 4, 5 (ρ = 2). NSGA-II [15] StD Mean

Mean

3.57 4.00 3.99 4.00

19.32 19.76 383.85 1559.63

9.76 7.11 220.51 1444.32

35.88 32.61 57.88 69.87

15.20 8.79 13.61 14.51

3.63 3.95 4.00 4.00

35.55 27.54 23.05 19.60

16.10 9.48 5.23 3.64

31.93 21.81 15.97 12.97

8.34 4.56 3.33 2.87

3.98 4.00 4.00 4.00

26.93 15.45 10.36 7.31

8.21 3.33 1.97 1.44

27.91 19.78 17.58 14.72

12.02 7.39 3.05 3.35

3.64 3.70 4.00 4.00

27.32 15.83 100 10.51 8.22

15.55 3.32 2.46 2.38

10−1

−2

10−1

10−2 0

0

500

1000

1500

average relative score difference

5.07 365.24 1126.90 1332.14

100

10

HypE [1] StD

15.17 659.72 3095.35 8768.99

average relative score difference

average relative score difference

M 2 3 4 5 M 2 3 4 5 M 2 3 4 5 M 2 3 4 5

R

CDEMOASF StD R Mean DTLZ1 3.33 0.14 0.17 2.99 0.80 0.83 2.62 10.38 5.79 2.80 44.61 28.70 DTLZ2 3.37 1.03 0.44 3.05 3.34 2.82 3.00 3.22 4.27 2.99 3.71 5.03 WFG1 3.02 0.48 0.57 3.00 1.92 2.23 3.00 2.19 3.57 2.96 2.45 3.93 WFG2 3.35 4.41 3.11 3.30 3.74 2.51 3.00 2.32 2.25 2.98 1.88 1.92

CDEMO-HVASF StD Mean R

R 1.60 1.45 1.00 1.00

9.38 1.42 480.44 471.22

56.54 4.98 340.52 458.37

1.50 1.56 2.39 2.20

1.53 1.44 1.52 1.38

0.50 2.59 3.71 5.40

1.27 2.60 2.45 2.88

1.47 1.56 1.48 1.63

1.69 1.78 1.85 1.86

0.33 0.91 1.48 1.96

0.24 0.59 0.87 1.13

1.31 1.22 1.15 1.18

1.32 1.35 1.44 1.64

5.84 4.14 3.28 1.30

8.42 4.64 3.25 0.81

1.69 1.65 1.56 1.38

100

10−1

NSGA-II HypE CDEMOASF NSGA-II HypECDEMO-HVASF

NSGA-II HypE 10 CDEMOASF CDEMO-HVASF −2

500

0

1000

500 generation

1000

1500

CDEMOASF CDEMO-HVASF

1500

generation

generation

(a) M = 3

(b) M = 5

Figure 15: Average ARSD attained by different methods applied to WFG1 with (a) M = 3 and (b) M = 5 objectives.

ronment, updating only one solution in a single iteration. Lastly, since CDEMOASF was confirmed to outperform CDEMOL1 , we parametrized CDEMOM ODEL and CDEMO-HVM ODEL with M ODEL = ASF . The differences between CDEMOASF and CDEMO-HVASF are summarized in Table 11. Apart from comparing the two interactive evolutionary hybrids, we account for the performance of NSGA-II and HypE [1], which can be seen as the non-interactive counterparts of, respectively, CDEMOASF and CDEMO-HVASF . This allows us to demonstrate how the results are influenced not only by the use of either front- or indicator-based scheme, but also by the employment of interactive and a posteriori methods. We evaluated the algorithms on DTLZ12 and WFG1-2 with M = 2, . . . , 5 objectives. Since CDEMO-HVASF and HypE run in a steady-state environment, to make their comparison with the generational algorithms CDEMOASF and NSGA-II fair, we suitably increased the number of iterations for the indicator-based methods. Specifically, CDEMO-HVASF and HypE were run for GEN × N iterations. In this way, both types of algorithms were allowed to generate the same number of solutions throughout the evolutionary search. However, for clarity and consistency of discussion on the experimental results, we refer to N iterations performed by the steady-state algorithms as to a single generation. In Table 12, we provide average ARSD for the final populations constructed by the four accounted algorithms. The results indicate that the interactive methods outperformed the a posteriori EMOAs, For example, NSGA-II attained the average rank around 4.00 for all problems, whereas HypE performed slightly better, though still being significantly worse than the preference-based algorithms. These results were expected as both NSGA-II and HypE approximated an entire PF and thus a vast majority of solutions that they constructed were not relevant for the 30

DM. In turn, the interactive methods focused the search around the DM’s most preferred option. In Figure 15, we illustrate the average ARSD attained by all methods when applied to WFG1 with M = 3 and M = 5 objectives. These plots reveal that NSGA-II and HypE could not improve their results after 500th generation. In turn, the interactive methods, being driven by the progressively supplied DM’s pairwise comparisons, were able to bias the search to the DM’s relevant region in the objective space. When it comes to comparing CDEMOASF and CDEMO-HVASF , there was no clear winner. The relative performances of these algorithms vastly depended on the problem’s characteristics. For example, in case of DTLZ1, both CDEMO methods performed similarly when up to 3 objectives were involved. However, for M = 4 and M = 5 – CDEMOASF was superior to CDEMO-HVASF , attaining average ranks of 1.00. Let us remind that DTLZ1 involves a multi-modal distance function which has to be minimized during the optimization [18]. In this perspective, the advantage of CDEMOASF over CDEMO-HVASF can be attributed to a less restrictive search imposed by the generational evolutionary scheme, which could allow to escape the local optima and hence to get closer to the PF. When it comes to DTLZ2, both interactive methods performed equally well. This result was expected as DTLZ2 is the easiest benchmark, and thus CDEMOASF could attain as satisfactory results as its indicator-based counterpart. In turn, CDEMO-HVASF was superior to CDEMOASF when applied to WFG1. This problem is particularly challenging due to incorporating a bias [29]. Specifically, the closer a population to the PF, the more difficult it is to generate an offspring yielding an improvement. In this perspective, CDEMO-HVASF could reach the PF faster than CDEMOASF , offering solutions being more relevant to the DM. In particular, for M = 5 – the average ranks for CDEMO-HVASF and CDEMOASF were equal to, respectively, 1.18 and 1.84. The advantage of incorporating an indicator- rather than front-based evolutionary scheme is also visible in Figure 15. As far as WFG2 is concerned, CDEMOASF performed slightly better than CDEMO-HVASF when M = 2 or M = 3 objectives were involved. However, for M = 4 there was no significant difference between these two methods, while for M = 5 – CDEMO-HVASF outperformed CDEMOASF (the average ranks equal to, respectively, 1.38 and 1.64). Let us remind that WFG2 uses a multi-modal function to model the shape of the PF. In this regard, it is challenging to explore different regions in the PF. CDEMO-HVASF incorporating hypervolume to construct a more diversified population should have a greater chance than CDEMOASF for reaching any region in the PF. Hence, with the increase in the number of objectives, CDEMO-HVASF attains a competitive advantage over CDEMOASF , which is confirmed by the experimental results. Although CDEMO-HVASF could perform significantly better than CDEMOASF for some specific optimization problems, its computation time was much higher due to employing the hypervolume measure. In e-Appendix 1, we discuss comprehensive computation times supporting this claim. For instance, for W F G1 and M = 5, CDEMOHVASF was around 180 times slower than CDEMOASF . 6. Conclusions In this paper, we proposed a novel interactive evolutionary algorithm for multiple objective optimization. Its basic variant revises NSGA-II to incorporate the progressively provided DM’s pairwise comparisons. Such indirect preference information is represented with a systematically contracted preference cone constructed by means of a set of either achievement scalarizing or Lα preference functions. The inclusion of the solutions in the preference cone is used along with the non-dominated fronts to bias the evolutionary search towards the DM’s most favourable region in the objective space. We distinguished a pair of sorting schemes, called DCEMO and CDEMO, differing in terms of the applied order of sorting criteria, i.e., non-dominated fronts and quality of solutions defined by the cone-based preference function. Moreover, the proposed algorithm was coupled with some static or novel dynamic interaction patterns that decide upon when the DM should be questioned for additional pairwise comparisons. We also introduced a novel visualization technique, called ViPEMO, which illustrates a typical performance of the interactive EMOA, clearly indicating when and which regions in the objective space were discovered by the optimization algorithm. The algorithm’s performance was verified within a thorough experimental study. Firstly, we compared different variants of the proposed method, revealing that it is more beneficial to apply preference cones based on ASF 31

rather than Lα -norm. The explanation of numerical outcomes was supported with the graphical analysis, revealing significant differences in the imposed evolutionary pressure. Moreover, the comparison of DCEMO and CDEMO indicated that in case of dealing with few objectives or problems being easy to solve, it is more beneficial to promote solutions contained in the systematically updated preference cones. However, with the increase in the number of objectives, when more solutions become incomparable in terms of dominance, the role of non-dominated sorting is diminished. As a result, both DCEMO and CDEMO perform similarly. Secondly, the best performing variant of the proposed algorithm CDEMOASF , was compared with some stateof-the-art interactive evolutionary hybrids: ECC-MRW, NEMO-II, NEMO-0, and P&K, when simulating different DM’s answering policies. These methods differed in terms of the incorporated preference model and accounting for the robustness concerns. The results revealed some useful remarks for designing interactive EMOAs, indicating that better solutions could be obtained when the model incorporated by the method and the DM’s preference function aligned and the search was guided by means of a set of preference model instances rather than a single arbitrarily selected one. Thirdly, we assessed the performance of CDEMOASF incorporating different interaction patterns. We observed that the results improved with the interactions equally distributed over the optimization run. Such a strategy ensured a systematic renewal of the evolutionary pressure, while giving the algorithm enough time to evolve a variety of solutions in the sub-region indicated by the contracted preference cone. We also demonstrated that it was beneficial to delay the preference elicitation when the fitness landscape was getting more complex with the algorithm’s progression towards the PF. This was confirmed by the application of both static (P P -ESI P OS ) and dynamic (RDS th ) interaction patterns, e.g., to WFG1 with at least 4 objectives. The former pattern postponed the first interaction with the DM due to some pre-defined rule, whereas the latter one delayed the DM’s questioning in case it was possible to generate an offspring yielding an improvement without his/her impact. Let us emphasize that all variants of CDEMOASF performed well when supplied with a very scarce preference information, i.e., 10 pairwise comparisons over all generations. This means that the cone contraction methods make a good use of few preference information pieces for directing the search to the DM’s most preferred of the objective space. Fourthly, although our main aim was to assess the benefits of incorporating the cone contraction procedures into front-based generational EMO algorithms and to compare the novel methods with the state-of-the-art algorithms using the same evolutionary framework, we also investigated to what extent the employed evolutionary base affects the performance of preference-based EMOAs. For this reason, we revised the indicator-based IEMO/I algorithm to incorporate preference cones. Thus obtained method, called CDEMO-HV, was compared against its generational front-based counterpart CDEMO. The experimental evaluation confirmed the advantage of CDEMO-HV over CDEMO in terms of discovering a population better corresponding to the DM’s preferences for more challenging optimization problems which involved a bias or a greater number of objectives. Nonetheless, such a competitive advantage was attained at the cost of increased computational complexity. We envisage the following future developments extending the current study. First, we will avoid an a priori selection of a preference model through maintaining a pool of such models and learning – during an optimization run – which model is more suitable for representing the DM’s preferences. We will also develop some active learning algorithms [9, 10, 64] for selecting pairs of solutions to be compared by the DM in view of maximizing the information gain from his/her answers and minimizing the number of interactions needed for constructing a satisfactory solution. Moreover, the dynamic interaction patterns can be enriched to account for the characteristics of solutions with respect to their inclusion in the contracted preference cone rather than non-dominance relation. Furthermore, we will extend the study concerning an impact of the evolutionary base on the performance of preference-based EMOAs. In this regard, we will propose and test a variety of indicator- and decomposition-based variants of the cone contraction algorithms for interactive evolutionary MOO. What is more, we will adapt the proposed algorithm for dealing with group optimization problems in order to handle potentially conflicting preferences of multiple DMs [34]. Finally, while we evaluated the proposed preference-based cone contraction algorithms on the DTLZ and WFG benchmark problems, in future studies we plan to consider a more challenging MaOP test suite [42]. 32

Acknowledgements The work of Milosz Kadzi´ nski was supported by the Polish Ministry of Science and Higher Education, grant no. 0296/IP2/2016/74. Michal K. Tomczyk acknowledges support from the Polish National Science Centre, grant no. DEC-2016/23/N/ST6/03795. The authors thank three anonymous referees and the editor for their remarks which helped us to significantly improve the previous versions of the paper. [1] Bader, J., Zitzler, E., 2011. HypE: An algorithm for fast hypervolume-based many-objective optimization. Evolutionary Computation 19 (1), 45–76. [2] Battiti, R., Passerini, A., 2010. Brain-computer evolutionary multiobjective optimization: a genetic algorithm adapting to the decision maker. IEEE Transactions on Evolutionary Computation 14, 671–687. [3] Belton, V., Branke, J., Eskelinen, P., Greco, S., Molina, J., Ruiz, F., Slowi´ nski, R., 2008. Interactive multiobjective optimization from a learning perspective. In: Branke, J., Deb, K., Miettinen, K., Slowi´ nski, R. (Eds.), Multiobjective Optimization: Interactive and Evolutionary Approaches. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 405–433. [4] Branke, J., Corrente, S., Greco, S., Slowi´ nski, R., Zielniewicz, P., 2016. Using Choquet integral as preference model in interactive evolutionary multiobjective optimization. European Journal of Operational Research 250 (3), 884–901. [5] Branke, J., Deb, K., Miettinen, K., Slowi´ nski, R. (Eds.), 2008. Multiobjective Optimization: Interactive and Evolutionary Approaches. Springer-Verlag Berlin Heidelberg, Germany. [6] Branke, J., Greco, S., Slowi´ nski, R., Zielniewicz, P., 2010. Interactive evolutionary multiobjective optimization driven by robust ordinal regression. Bulletin of the Polish Academy of Sciences Technical Sciences 58 (3), 347–358. [7] Branke, J., Greco, S., Slowi´ nski, R., Zielniewicz, P., 2015. Learning value functions in interactive evolutionary multiobjective optimization. IEEE Transactions on Evolutionary Computation 19 (1), 88–102. [8] Branke, J., Kauler, T., Schmeck, H., 2001. Guidance in evolutionary multi-objective optimization. Advances in Engineering Software 32, 499–507. [9] Ciomek, K., Kadzi´ nski, M., Tervonen, T., 2017. Heuristics for prioritizing pair-wise elicitation questions with additive multi-attribute value models. Omega 71, 27–45. [10] Ciomek, K., Kadzi´ nski, M., Tervonen, T., 2017. Heuristics for selecting pair-wise elicitation questions in multiple criteria choice problems. European Journal of Operational Research 262 (2), 693–707. [11] Deb, K., Agrawal, R. B., 1994. Simulated binary crossover for continuous search space. Tech. Rep. IITK/ME/SMD-94027, Indian Institute of Technology, Kanpur, India. [12] Deb, K., Goyal, M., 1996. A combined genetic adaptive search (GeneAS) for engineering design. Computer Science and Informatics 26, 30–45. [13] Deb, K., Jain, H., 2014. An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: Solving problems with box constraints. IEEE Transactions on Evolutionary Computation 18 (4), 577–601. [14] Deb, K., Kumar, A., 2007. Light beam search based multi-objective optimization using evolutionary algorithms. In: IEEE Congress on Evolutionary Computation. pp. 2125–2132. [15] Deb, K., Pratap, A., Agrawal, S., Meyarivan, T., 2000. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6 (2), 182–197. [16] Deb, K., Sinha, A., Korhonen, P. J., Wallenius, J., 2010. An interactive evolutionary multiobjective optimization method based on progressively approximated value functions. IEEE Transactions on Evolutionary Computation 14 (5), 723–739. [17] Deb, K., Sundar, J., Rao, N. U. B., Chaudhuri, S., 2006. Reference point based multi-objective optimization using evolutionary algorithms. International Journal of Computational Intelligence Research. 2 (3), 273–286. [18] Deb, K., Thiele, L., Laumanns, M., Zitzler, E., 2005. Scalable test problems for evolutionary multiobjective optimization. In: Evolutionary Multiobjective Optimization: Theoretical Advances and Applications. Springer London, London, pp. 105–145. [19] Deb, K., Zope, P., Jain, A., 2003. Distributed computing of Pareto-optimal solutions with evolutionary algorithms. In: Evolutionary Multi-Criterion Optimization: Second International Conference, EMO 2003, Faro, Portugal, April 8–11, 2003. Proceedings. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 534–549. [20] DeShazo, J. R., Fermo, G., 2002. Designing choice sets for stated preference methods: The effects of complexity on choice consistency. Journal of Environmental Economics and Management 44 (1), 123–143. [21] Fonseca, C. M., Fleming, P., 1993. Genetic algorithms for multiobjective optimization: Formulationdiscussion and generalization. In: 5th International Conference on Genetic Algorithms. pp. 416–423. [22] Fonseca, C. M., Fleming, P. J., 1996. On the performance assessment and comparison of stochastic multiobjective optimizers. In: Voigt, H.-M., Ebeling, W., Rechenberg, I., Schwefel, H.-P. (Eds.), Parallel Problem Solving from Nature — PPSN IV. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 584–593. [23] Fowler, J. W., Gel, E. S., K¨ oksalan, M., Korhonen, P., Marquis, J. L., Wallenius, J., 2010. Interactive evolutionary multi-objective optimization for quasi-concave preference functions. European Journal of Operational Research 206 (2), 417–425. [24] Garc´ıa-N´ ajera, A., L´ opez-Jaimes, A., 2018. An investigation into many-objective optimization on combinatorial problems: Analyzing the pickup and delivery problem. Swarm and Evolutionary Computation 38, 218 – 230. [25] Gong, D., Sun, J., Ji, X., 2013. Evolutionary algorithms with preference polyhedron for interval multi-objective optimization problems. Information Sciences 233, 141–161.

33

[26] Gong, M., Liu, F., Zhang, W., Jiao, L., Zhang, Q., 2011. Interactive MOEA/D for multi-objective decision making. In: Proceedings of the 13th Annual Conference on Genetic and Evolutionary Computation. GECCO ’11. pp. 721–728. [27] Greco, S., Matarazzo, B., Slowi´ nski, R., 2010. Dominance-based rough set approach to interactive evolutionary multiobjective optimization. In: Greco, S., Marques Pereira, R. A., Squillante, M., Yager, R. R., Kacprzyk, J. (Eds.), Preferences and Decisions: Models and Applications. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 225–260. [28] Greenwood, G. W., Hu, X. S., D’Ambrosio, J. G., 1996. Fitness functions for multiple objective optimization problems: Combining preferences with Pareto rankings. In: Foundations of genetic algorithms. pp. 437–455. [29] Huband, S., Hingston, P., Barone, L., While, L., 2006. A review of multiobjective test problems and a scalable test problem toolkit. IEEE Transactions on Evolutionary Computation 10 (5), 477–506. [30] Jacquet-Lagr` eze, E., Siskos, Y., 2001. Preference disaggregation: 20 years of MCDA experience. European Journal of Operational Research 130 (2), 233–245. [31] Jaszkiewicz, A., 2007. Interactive multiobjective optimization with the Pareto memetic algorithm. Foundation of Computing and Decision Sciences 32 (1), 15–32. [32] Kadzi´ nski, M., Ghaderi, M., Wasikowski, J., Agell, N., 2017. Expressiveness and robustness measures for the evaluation of an additive value function in multiple criteria preference disaggregation methods: An experimental analysis. Computers & Operations Research 87, 146–164. [33] Kadzi´ nski, M., Slowi´ nski, R., 2012. Interactive robust cone contraction method for multiple objective optimization problems. International Journal of Information Technology & Decision Making 11 (02), 327–357. [34] Kadzi´ nski, M., Tomczyk, M., 2017. Interactive evolutionary multiple objective optimization for group decision incorporating valuebased preference disaggregation methods. Group Decision and Negotiation 26 (4), 693–728. [35] Kadzi´ nski, M., Tomczyk, M., Slowi´ nski, R., 2018. Interactive cone contraction for evolutionary multiple objective optimization. In: Gaw¸eda, A. E., Kacprzyk, J., Rutkowski, L., Yen, G. G. (Eds.), Advances in Data Analysis with Computational Intelligence ˙ Methods: Dedicated to Professor Jacek Zurada. Springer International Publishing, Cham, pp. 293–309. [36] Kadzi´ nski, M., Tervonen, T., Tomczyk, M. K., Dekker, R., 2017. Evaluation of multi-objective optimization approaches for solving green supply chain design problems. Omega 68, 168–184. [37] Kallio, M., Halme, M., 2013. Cone contraction and reference point methods for multi-criteria mixed integer optimization. European Journal of Operational Research 229 (3), 645–653. [38] Karahan, I., K¨ oksalan, M., 2010. A territory defining multiobjective evolutionary algorithms and preference incorporation. IEEE Transactions on Evolutionary Computation 14 (4), 636–664. [39] Karakaya, G., K¨ oksalan, M., Ahipa¸sao˘ glu, S., 2018. Interactive algorithms for a broad underlying family of preference functions. European Journal of Operational Research 265 (1), 248–262. [40] Keeney, R., Raiffa, H., 1976. Decisions with Multiple Objectives: Preferences and Value Trade-offs. Cambridge University Press, Cambridge. [41] Korhonen, P., Wallenius, J., Zionts, S., 1984. Solving the discrete multiple criteria problem using convex cones. Management Science 30 (11), 1336–1345. [42] Li, H., Deb, K., Zhang, Q., Suganthan, P., Chen, L., 2019. Comparison between MOEA/D and NSGA-III on a set of novel many and multi-objective benchmark problems with challenging difficulties. Swarm and Evolutionary Computation 46, 104–117. [43] Li, H., Zhang, Q., 2009. Multiobjective optimization problems with complicated Pareto sets, MOEA/D and NSGA-II. IEEE Transactions on Evolutionary Computation 13 (2), 284–302. [44] Li, K., Chen, R., Savi´ c, D., Yao, X., 2019. Interactive decomposition multiobjective optimization via progressively learned value functions. IEEE Transactions on Fuzzy Systems 27 (5), 849–860. [45] Li, K., Deb, K., Zhang, Q., Kwong, S., 2015. An evolutionary many-objective optimization algorithm based on dominance and decomposition. IEEE Transactions on Evolutionary Computation 19 (5), 694–716. [46] Li, L., Wang, Y., Trautmann, H., Jing, N., Emmerich, M., 2018. Multiobjective evolutionary algorithms based on target region preferences. Swarm and Evolutionary Computation 40, 196–215. [47] Marquis, J., Gel, E. S., Fowler, J. W., K¨ oksalan, M., Korhonen, P., Wallenius, J., 2015. Impact of number of interactions, different interaction patterns, and human inconsistencies on some hybrid evolutionary multiobjective optimization algorithms. Decision Sciences 46 (5), 981–1006. [48] Mohammadi, A., Omidvar, M. N., Li, X., Deb, K., 2014. Integrating user preferences and decomposition methods for many-objective optimization. In: IEEE Congress on Evolutionary Computation. pp. 421–428. [49] Montibeller, G., Winterfeldt, D., 2015. Cognitive and motivational biases in decision and risk analysis. Risk Analysis 35 (7), 1230– 1251. [50] Nebro, A. J., Ruiz, A. B., Barba-Gonz´ alez, C., Garc´ıa-Nieto, J., Luque, M., Aldana-Montes, J. F., 2018. InDM2: Interactive dynamic multi-objective decision making using evolutionary algorithms. Swarm and Evolutionary Computation 40, 184–195. [51] Phelps, S., K¨ oksalan, M., 2003. An interactive evolutionary metaheuristic for multiobjective combinatorial optimization. Management Science 49 (12), 1726–1738. [52] Purshouse, R. C., Deb, K., Mansor, M. M., Mostaghim, S., Wang, R., 2014. A review of hybrid evolutionary multiple criteria decision making methods. In: 2014 IEEE Congress on Evolutionary Computation (CEC). pp. 1147–1154. [53] Qi, Y., Li, X., Yu, J., Miao, Q., 2018. User-preference based decomposition in MOEA/D without using an ideal point. Swarm and Evolutionary Computation.

34

[54] Qu, B., Zhu, Y., Jiao, Y., Wu, M., Suganthan, P., Liang, J., 2018. A survey on multi-objective evolutionary algorithms for the solution of the environmental/economic dispatch problems. Swarm and Evolutionary Computation 38, 1–11. [55] Rahmati, S. H. A., Ahmadi, A., Karimi, B., 2018. Multi-objective evolutionary simulation based optimization mechanism for a novel stochastic reliability centered maintenance problem. Swarm and Evolutionary Computation 40, 255–271. [56] Saborido, R., Ruiz, A. B., Luque, M., Miettinen, K., 2019. IRA-EMO: Interactive method using reservation and aspiration levels for evolutionary multiobjective optimization. In: Deb, K., Goodman, E., Coello Coello, C. A., Klamroth, K., Miettinen, K., Mostaghim, S., Reed, P. (Eds.), Evolutionary Multi-Criterion Optimization. Springer International Publishing, Cham, pp. 618–630. [57] Ser, J., Osaba, E., Molina, D., Yang, X., Salcedo-Sanz, S., Camacho, D., Das, S., Suganthan, P. N., Coello, C., Herrera, F., 2019. Bio-inspired computation: Where we stand and what’s next. Swarm and Evolutionary Computation 48, 220–250. [58] Sinha, A., Korhonen, P., Wallenius, J., Deb, K., 2014. An interactive evolutionary multi-objective optimization algorithm with a limited number of decision maker calls. European Journal of Operational Research 233 (3), 674–688. [59] Sinha, A., Malo, P., Kallio, M., 2018. Convex preference cone-based approach for many objective optimization problems. Computers & Operations Research 95, 1–11. [60] Sun, X., Gong, D., Jin, Y., Chen, S., 2013. A new surrogate-assisted interactive genetic algorithm with weighted semisupervised learning. IEEE Transactions on Cybernetics 43 (2), 685–698. [61] Thiele, L., Miettinen, K., Korhonen, P. J., Molina, J., 2009. A preference-based evolutionary algorithm for multi-objective optimization. Evolutionary Computation 17 (3), 411–436. [62] Todd, D. S., Sen, P., 1999. Directed multiple objective search of design spaces using genetic algorithms and neural networks. In: Genetic and Evolutionary Computation Conference. Morgan Kaufmann, pp. 1738–1743. [63] Tomczyk, M., Kadzi´ nski, M., 2019. Decomposition-based interactive evolutionary algorithm for multiple objective optimization. IEEE Transactions on Evolutionary Computation (accepted for publication), dOI:10.1109/TEVC.2019.2915767. [64] Tomczyk, M., Kadzi´ nski, M., 2019. EMOSOR: Evolutionary multiple objective optimization guided by interactive stochastic ordinal regression. Computers & Operations Research 108, 134–154. [65] Tomczyk, M. K., Kadzi´ nski, M., 2019. Robust indicator-based algorithm for interactive evolutionary multiple objective optimization. In: Proceedings of the Genetic and Evolutionary Computation Conference. GECCO ’19. ACM, New York, NY, USA, pp. 629–637. [66] Vanneschi, L., Henriques, R., Castelli, M., 2017. Multi-objective genetic algorithm with variable neighbourhood search for the electoral redistricting problem. Swarm and Evolutionary Computation 36, 37–51. [67] Vesikar, Y., Deb, K., Blank, J., 2018. Reference point based nsga-iii for preferred solutions. In: 2018 IEEE Symposium Series on Computational Intelligence (SSCI). pp. 1587–1594. [68] Wagner, T., Trautmann, H., 2010. Integration of preferences in hypervolume-based multiobjective evolutionary algorithms by means of desirability functions. IEEE Transactions on Evolutionary Computation 14 (5), 688–701. [69] Wang, H., Olhofer, M., Jin, Y., 2017. A mini-review on preference modeling and articulation in multi-objective optimization: current status and challenges. Complex & Intelligent Systems 3 (4), 233–245. [70] Wang, R., Purshouse, R. C., Fleming, P. J., 2013. Preference-inspired coevolutionary algorithms for many-objective optimization. IEEE Transactions on Evolutionary Computation 17 (4), 474–494. [71] Xin, B., Chen, L., Chen, J., Ishibuchi, H., Hirota, K., Liu, B., 2018. Interactive multiobjective optimization: A review of the state-of-the-art. IEEE Access 6, 41256–41279. [72] Zhang, Q., Li, H., 2007. MOEA/D: A multiobjective evolutionary algorithm based on decomposition. IEEE Transactions on Evolutionary Computation 11 (6), 712–731. [73] Zhou, A., Qu, Y. Y., Li, H., Zhao, S. Z., Suganthan, P. N., Zhang, Q., 2011. Multiobjective evolutionary algorithms: A survey of the state of the art. Swarm and Evolutionary Computation 1 (1), 32–49. [74] Zitzler, E., K¨ unzli, S., 2004. Indicator-based selection in multiobjective search. In: Parallel Problem Solving from Nature - PPSN VIII: 8th International Conference, Birmingham, UK, September 18-22, 2004. Proceedings. pp. 832–842. [75] Zitzler, E., Laumanns, M., Thiele, L., 2002. SPEA2: Improving the Strength Pareto Evolutionary Algorithm. In: Evolutionary Methods Design, Optimization, and Control with Applications to Industrial Problems (EUROGEN). pp. 95–100. [76] Zitzler, E., Thiele, L., 1999. Multiobjective evolutionary algorithms: a comparative case study and the Strength Pareto Approach. IEEE Transactions on Evolutionary Computation 3 (4), 257–271.

35

Research highlights

• We propose interactive evolutionary algorithms for multiple objective optimization • Interactive cone contraction is used to bias the evolutionary search • The proposed methods outperform selected state-of-the-art algorithms • We introduce a novel approach for visualizing a progress in the evolutionary search • Novel interaction patterns are used to control the questioning of a decision maker