A modified indicator-based evolutionary algorithm (mIBEA)

Multi-objective evolutionary algorithms (MOEAs) based on the concept of Pareto-dominance have been successfully applied to many real-world optimisation problems. Recently, research interest has shifted towards indicator-based methods to guide the search process towards a good set of trade-off solutions. One commonly used approach of this nature is the indicator-based evolutionary algorithm (IBEA). In this study, we highlight the solution distribution issues within IBEA and propose a modification of the original approach by embedding an additional Pareto-dominance based component for selection. The improved performance of the proposed modified IBEA (mIBEA) is empirically demonstrated on the well-known DTLZ set of benchmark functions. Our results show that mIBEA achieves comparable or better hypervolume indicator values and epsilon approximation values in the vast majority of our cases (13 out of 14 under the same default settings) on DTLZ1-7. The modification also results in an over 8-fold speed-up for larger populations.


I. INTRODUCTION
In multi-objective optimisation, where multiple objectives are optimised simultaneously, the goal is to find a set of P areto-optimal solutions known as the Pareto front (PF).The PF consists of a set of solutions that are not dominated by each other, which are termed as non-dominated solutions, representing the trade-off that exists between different objectives.This dominance relation, also known as Pareto dominance relation, (≺) is defined between solutions x 1 and x 2 .We say w.l.o.g., in a minimisation problem that x 1 dominates x 2 (x 1 ≺ x 2 ) if and only if f i (x 1 ) ≤ f i (x 2 ) for all k objective functions (i ∈ {1, . . .k}), and f i (x 1 ) < f i (x 2 ) for at least one objective function.One of the difficulties in multi-objective search is to find a set of solutions to minimise the distance to the true Pareto front (PF) while maintaining the diversity of the solution set in the objective space.
Multi-objective evolutionary algorithms (MOEAs) are widely used to solve various multi-objective optimisation problems [16] and are considered to be general and robust search mechanisms [1].Evolutionary algorithms (EAs) are a class of stochastic optimisation methods that mimic the process of evolution in nature.Examples of EAs, such as genetic algorithms, evolutionary programming, and evolution strategies [1] operate on a set of solutions using the basic principles of natural evolution: selection, reproduction by means of recombination and mutation.The algorithmic difference between single-objective EAs and MOEAs is that additionally in an MOEA, the multiple objectives of a solution must be transformed into a single fitness value to facilitate the comparison of individual solutions [6].Thus, MOEAs often vary in the method of this transformation, considering the balance of convergence and diversification during the search.Some MOEAs, such as NSGA-II [7], SPEA2 [18], and AGE [3,13] incorporate Pareto-based ranking of the individuals and an additional density measurement (crowding distance in NSGA-II, k-th nearest neighbour in SPEA2,dominance in AGE) in the objective space.However, between two non-dominated solutions, purely Pareto-based MOEAs are not able to ascertain which solutions have better potential for convergence.IBEA [17] was one of the earliest indicator-based MOEAs proposed in the literature.Originally, it came in two variants: one using -dominance for guidance, denoted IBEA , and another based on hypervolume, denoted IBEA HD ('HD' stands for hypervolume difference) which will be referred to as IBEA from this point onward in this study.IBEA associates a fitness value with each solution based on the selected indicator (hypervolume or ), attempting to guide the search towards the true PF.IBEA was shown to achieve significantly better performance on various benchmark functions than NSGA-II and SPEA2 [17], however the distribution of the solutions found by IBEA has rarely been reported or discussed in detail.In this paper, we propose a modified variant of IBEA HD , termed mIBEA, which adds a Pareto-based element to this indicator-based method, analysing the distribution of nondominated solutions found.For further information about MOEAs and indicator-based MOEAs in particular, we refer the interested reader to [5,6,16].
The remainder of the paper is structured as follows.We first describe the original IBEA in Section II-A, then present observations of the solution distributions observed using existing MOEAs in Section II-B.The proposed mIBEA is introduced in Section III.Experimental results comparing IBEA and mIBEA are presented in Section IV.We draw conclusions and provide 978-1-5090-4601-0/17/$31.00 c 2017 IEEE ideas for future work in Section V.

II. INDICATOR-BASED EVOLUTIONARY ALGORITHM
Since the focus in the paper is on the hypervolume variant of IBEA, i.e.IBEA HD , we first give a detailed description of the original IBEA.We then provide visualisations and observations for the non-dominated solution sets found using IBEA and two other existing MOEAs, one Pareto dominance-based (NSGA-II [7]) and one indicator-based (SMS-EMOA [2]).

A. Description of IBEA
The core idea of IBEA HD is to employ a binary hypervolume indicator in the selection process, when determining which solutions survive to the next generation.The binary hypervolume indicator assigns a real-valued number to two solution sets with respect to a reference point.The formula of I HD (A, B) is defined as space that is dominated by population B, but not by A, shown in Equation ( 1) [17].The pseudocode of the original IBEA is given in Alg. 1. IBEA first randomly generates an initial population in Step1, then the following steps loop until the stopping criterion is satisfied.The objective values are scaled and the fitness is assigned to each individual in Step2 and Step3.Step4 performs environmental selection, iteratively removing the worst individual in the population P based on indicator value until µ individuals remain (this step will do nothing in the first iteration of the algorithm as P = µ).Upon removal of each solution, the indicator values of the remaining solutions must be updated.This step continues until the number of solutions in P does not exceed µ.The standard mating selection (Step6) and variation (Step7) steps are performed to generate new individuals and add them to the population P .

B. Solution distribution issues in existing MOEAs
It is common in the MOEA literature to use various quantitative indicators, such as hypervolume, -approximation and inverted generational distance (IGD) [19], to compare the performance of different MOEAs.However, the distribution of non-dominated solutions over the resultant solution fronts is relatively rarely reported, at least in part due to the lack of quantitative indicators to measure it in the objective space.
One method for measuring the spread of solutions over the objective space we will use here is Generalised Spread (GS) [15].GS computes the average Euclidean distance of any two consecutive solutions [7] (or nearest neighbours in [15]) in the population, and considers this average in the context of the minimal achieved distances to the extreme points of the true Algorithm 1: (Adaptive) IBEA HD as described in [4,17] Inputs: µ: population size; N : number of total solution evaluations; ρ: objective values scaling factor; κ: indicator value (hypervolume difference) scaling factor Output: A: Pareto set approximation Step1: Initialisation: Randomly generate the initial population P and set up evaluation counter m.Step2: Scale objective values: 1) find the lower (b i = min x∈P f i (x)) and upper 3) scale each objective to the interval [0, 1]; 1) calculate indicator values I(x 1 , x 2 ) using the scaled objective values f i instead of original objective values f i , and determine the maximum absolute indicator value Step4: Environmental Selection: iterate the following three sub-steps until the size of population P ≤ µ. 1) choose an individual x ∈ P with the smallest fitness value, i.e., F (x ) ≤ F (x) for all x ∈ P 2) remove x from P 3) update fitness values of all the remaining individuals x ∈ P as F (x) = F (x) + exp −I(x ,x)/c * κ Step5: Termination criteria: if m ≥ N , return the non-dominated solutions in P as A; otherwise, continue.Step6: Mating Selection: apply binary tournament selection operator to P to select two parents from P .Step7: Variation: perform crossover operator on the two parents to generate two offspring, and use an mutation operator on one of the offspring.The resulting offspring is then added to P .Increment iteration counter (m = m + 1).Go to Step 2.
PF.The smaller the value of GS, the better spread the resultant front has as the variance is reduced and/or points closer to the extreme points have been found.The ideal spread is achieved when GS equals to 0 as this means that the extreme points are found and the solutions are evenly spread out.In this study, for the purpose of analysis, we also use a new metric -border fraction (bf) to provide a certain quantitative measurement of the solution distribution of different PFs.We define border fraction as the fraction of non-dominated points lying in a pre-defined border or area of the resultant PF.More precisely, assume the lower bound of each objective function is 0, the border threshold for each objective is denoted as θ i , then the border is the volume defined by (0 where k is number of objectives in the multi-objective problem.Given ndi non-dominated solutions in the population, the border fraction bf value is the proportion of non-dominated individuals in the population which lie within in this border.Note that bf is not a general metric to measure the solution distribution of a Pareto optimal (or approx.)set.
As a motivating example, we apply three well-known MOEAs, NSGA-II, SMS-EMOA and IBEA, to a sample of three DTLZ problems [8], using the default settings for each method.Figure 1 shows the non-dominated solutions found in 100 independent runs for each MOEA on each problem, except for SMS-EMOA which was run only 10 times due to its high running time.The maximum number of evaluations for each run was 100,000.The bf and ndi values are also provided on the top right corner of each plot.Firstly, we can see that the solutions from NSGA-II can cover the whole fronts of all the problems DTLZ1-3.Secondly, SMS-EMOA can reach almost full coverage of the PF for DTLZ1, however, on DTLZ2 and 3 the PFs have clear gaps between the border and central area.Thirdly, similar to SMS-EMOA, IBEA also shows clear gaps (empty space) on the PF of DTLZ2.Also, solutions by the original IBEA are heavily concentrated in the corner points on DTLZ1 and DTLZ3.

NSGA-II SMS-EMOA IBEA DTLZ1
The GS measurement of the solution fronts obtained from NSGA-II, SMS-EMOA and IBEA (as illustrated in Figure 1) is shown in Table I.The best GS on each benchmark function is highlighted in bold.Not surprisingly, IBEA performs poorly w.r.t.GS among all the three algorithms on DTLZ1-3.However, interestingly, although SMS-MOEA shows gaps on the resultant PFs of DTLZ2 and 3, the GS indicator still shows that it performs best on both problems.This suggests that under certain situations (e.g., PFs of DTLZ2 and 3 from SMS-EMOA) GS cannot capture the distribution of the PF, i.e., the spread calculated is not equivalent to the distribution.
In the following, we use bf and ndi to capture the bias within the distribution towards extreme points on the solution fronts.The border threshold is set as 0.03 for DTLZ1 and 0.1 for both To the best of our knowledge, even though these uneven solution distributions have been plotted in articles before, e.g.[2,12], no further details are discussed in those papers.Although IBEA is often reported as having exceptional indicator value performance [9,17], the actual distribution of solutions that it obtains is surprisingly uneven.Additionally, mediocre performance [12,14] is reported on certain benchmark functions, possibly as a result of this uneven distribution.In the following section, we introduce a modified version of IBEA that attempts to overcome this issue.

III. PROPOSED MODIFIED IBEA (mIBEA)
The behaviour of IBEA observed in the previous section suggests that the search is trapped in certain regions of the search space when solving DTLZ1 and DTLZ3, and suffers from a lack of diversity in solutions.In other words, it seems that after a certain time point, IBEA is not able to produce solutions with high hypervolume contributions.A poor spread of non-dominated solutions when using IBEA was observed previously by Li et al. [11] in the context of a real world application, the multi-objective wind farm layout optimisation problem.The proposed method in that paper illustrated that the hybridisation of Pareto dominance-based MOEAs and IBEA could improve the overall performance, in terms of both convergence and diversification within the search.It is this observation that the modified IBEA, mIBEA, we introduce here is based on.
mIBEA excludes dominated solutions at each generation after the new population is created.To achieve this, we use the dominance-based sorting method of NSGA-II, however, other approaches for achieving this can be used as well.
The difference between mIBEA and IBEA is to add this one step before the scaling Step2 in IBEA (as described in Section II-A).The effect is that the scaling of the objective scores is no longer affected by dominated solutions that are far away from the best non-dominated solutions in the population, see Step2.1 of Algorithm 1 where minimum and maximum values across the entire population are determined.The pseudocode of mIBEA is given in Algorithm 2.

A. Experimental Design
The following experiments were performed on the full DTLZ problem set (DTLZ1-7) [8].Our analysis focuses on the three-objective DTLZ1 and 3, as qualitatively good coverage of the true Pareto front can be achieved relatively easily by different algorithms.Besides, these two problems exhibit different properties: DTLZ1 has a planar front in 3D, DTLZ3 has 1/8th of a sphere front in 3D.DTLZ2 and 4 share the same PF as DTLZ3 but differ in the overall domain of the objective values.For example, randomly drawn solutions in the initial population, have objective values of 100-150 in DTLZ1 and 500-1500 in DTLZ3 while the range is significantly smaller e.g., with 2 to 5 in DTLZ2.It is DTLZ1 and 3 where we are expecting the most benefit of the removal of dominated solutions, whereas the insertion of our processing step should be neutral on the other functions.
The settings for the 3D DTLZ problems follow the recommended settings from [8].We are using the algorithms and test functions as implemented in the jMetal framework [10].For a fair comparison, all algorithms utilise the same operators with the same corresponding parameters.Specifically, a binary tournament selection operator is used for selection.Simulated Binary Crossover (SBX) and polynomial mutation with probabilities of 0.9 and (1/number of parameters) respectively, are used to create new offspring.The distribution parameters of the crossover and mutation operators are η c = 20.0 and η m = 20.0,respectively.A population size of 100 is used with the total number of solution evaluations limited to 100,000.Other algorithmic parameters are set as the default settings from the original IBEA [17].100 runs of each algorithm are performed.The version of IBEA that we are using for comparison is the one provided by Brockhoff [4], which has fixed the hypervolume calculation bug in the original IBEA.
In the following, we will assess the performance of mIBEA qualitatively by inspecting the plots, quantitatively using performance indicators and by measuring the runtime.We analyse experimentally the effect that our patch has on IBEA.We do this on DTLZ1-7 and by varying the scaling parameter ρ and the population size µ.

B. Visualisation of Resultant Solution Distributions
The resultant populations under the default setting of mIBEA, with ρ = 2.0, µ = 100, are illustrated in Figure 2. Compared to the original IBEA presented in Figure 1, the coverage obtained by mIBEA on DTLZ1 and DTLZ3 is a clear improvement.More specifically, the 9890 non-dominated individuals (ndi) found by mIBEA for DTLZ1 spread cross the whole front surface, while the corresponding ndi from IBEA is only 616 and only distributed in the extreme points, resulting in much higher bf (0.9935) than mIBEA (0.5297).Similarly, mIBEA improves the original IBEA by generating many more non-dominated solutions (ndi = 5292) and solutions between the extreme points (bf = 0.7897) and the central area of the front are found, while the original IBEA only obtains 182 non-dominated solutions which only lie in the corner points.In order to provide more quantitative insights into the distribution of solutions found by each algorithm, Figure 3 compares the proportion of non-dominated solutions found in different regions on each DTLZ problem of the PF.For the three highlighted areas the proportion of the total number of non-dominated solutions contained in the area is provided, with the corresponding proportion of the actual PF surface in 3D the area represents given in the bracket.The boundaries from outermost triangle (blue area) to the innermost triangle (green area) of DTLZ1 are 0.03 and 0.1, for DTLZ3 are 0.1 and 0.2.
For example, in Figure 3b, 0.53(0.33)indicates that the fraction of ndi observed from mIBEA in the blue triangle area of DTLZ1 is 0.53, however the proportion of the total surface this area corresponds to is 0.33.In an ideal solution distribution, the fraction of ndi should be very close to the surface proportion, if one is interested in this coverage notion.
Figure 3 also shows that mIBEA could find new trade-off solutions between extreme points on both DTLZ problems that the original IBEA is not able to.In the case of mIBEA, 53% and 79% of non-dominated solutions are found in the edge regions of DTLZ1 and DTLZ3 respectively.Interestingly, a decreased fraction of ndi from the outside to the inside of the PF of DTLZ1 is observed, with a 'tiny' fraction (3%) of solutions lie in the middle area (green) area of DTLZ3.In the previous subsection, the objective values scaling factor ρ is fixed to 2.0.As the hypervolume contribution of each solution is calculated based on the scaled objective values (as described in Step2 and Step3 of IBEA, see Section II-A), here we examine the effect of modifying the value of ρ.Within IBEA and mIBEA, ρ directly affects the ranking of solutions, therefore influencing the preservation of elite solutions.
A series of ρ values within the range of 1.0 to 1024 are tested.From a set of preliminary experiments, we found that when ρ ∈ [1.0, 1.1], the PFs on DTLZ functions show no obvious improvement when using either the original IBEA or mIBEA.However when ρ > 1.1, some improvement was shown by IBEA and significant improvement by mIBEA.However, when ρ ≥ 2, the PFs remain relatively unchanged by both algorithms.Thus, we choose ρ = 1.0, 1.15, 2.0, 64.0 to plot the PFs from IBEA and mIBEA on DTLZ1 and 3 in Figure 4.
As we can see, the scaling parameter ρ in IBEA has little to no impact on the distribution of the solutions on DTLZ1 and DTLZ3, as the solution sets always collapse to the extreme points of the objective space no matter what the ρ value is.However, ρ can affect the performance of mIBEA.
We also compare the performance of IBEA and mIBEA w.r.t. the hypervolume and + indicators across the entire DTLZ family in Table III.The higher hypervolume or lower + , the better performance of the algorithm has.The reference point of hypervolume is set as 0.5 d for DTLZ1, where d is objective index, 1.0 d for the rest DTLZ problems, DTLZ7 where the reference point is (1,1,6) due to the third objective.Results that are better with statistical significance (Wilcoxon rank-sum test, 5% significance level) are highlighted in bold.It shows that mIBEA not only improves the coverage of the resultant fronts, but also the convergence, especially on DTLZ1, DTLZ3, and DTLZ7 under all the chosen ρ values, as well as most ρ values of DTLZ6, while performing similarly on DTLZ2, 4 and 5.Only in one of the 14 cases (with standard ρ = 2), mIBEA is statistically worse, while it is comparable or statistically better in 13 of 14 cases.
In summary, the observations of PF changes from original IBEA and mIBEA suggest the following: • The original IBEA performs poorly on benchmark functions with large objective spaces and comparatively small fronts, i.e.DTLZ1 and DTLZ3 and the change of objective value scaling factor ρ(≥ 1) does not affect (either positively or negatively) the performance; • Our mIBEA improved IBEA overall on DTLZ1-7 with wider PF coverage and better convergence w.r.t.hypervolume and + .Moreover, the performance comparison between NSGA-II and mIBEA using default settings (µ = 100, ρ = 2) (shown in Table II) suggests that mIBEA still does not beat NSGA-II on most DTLZ problems.Perhaps it is because the still existing gaps (Fig. 4) and non-uniform solution distribution (Fig. 3) on the mIBEA fronts.Further investigations are to be conducted.We also investigate the performance of both the original IBEA and mIBEA using larger population sizes under different ρ values.Figure 5 shows that when population size goes up to 1000, population size has little affect on PF coverage for DTLZ1 by the original IBEA.It even deteriorates for DTLZ3 on which the maximum of each objective of DTLZ3 is 1.0 on the true PF.However, the proposed mIBEA clearly covers greater parts of the PF surface on both DTLZ1 and DTLZ3 given various ρ values.

E. Running time of IBEA and mIBEA
As a positive side-effect, the focus on non-dominated solutions within mIBEA results in significant speed-ups of the runtime.Table IV shows that IBEA on average gives an over 8-fold speed-up compared to IBEA using population size of µ = 1000 (the ρ is set as default value 2.0).The reason is that Step2 in mIBEA (Section III) removes dominated solutions in every generation.Therefore, fewer solutions are persevered in the population, reducing the calculation time for indicator values i.e., hypervolume difference in the original IBEA.As the indicator value calculation is pairwise (Step3 in Section II-A), n 2 indicator value calculations are required, where n is the current population size.

V. DISCUSSION AND FUTURE WORK
The proposed mIBEA introduces the ranking from Pareto dominance-based multi-objective evolutionary algorithm into the indicator-based algorithm during the elite solution preservation process.The empirical results from mIBEA over the entire DTLZ benchmark functions show that mIBEA significantly improves the original IBEA on the coverage of resultant Pareto fronts, as well as hypervolume and + .Also, over 8fold speed-ups are obtained when using a larger population size.Also note that the proposed mIBEA does not introduce any additional parameter to the original IBEA.However, some gaps are still observable from original IBEA and mIBEA.Similar distributions can be observed from the hypervolume-based SMS-EMOA, although these are better than IBEA.We have performed an initial experiment on the SMS-EMOA's scaling parameter offset.However, no clear changes have been observed under different offset values.It suggests that the gaps from original IBEA, mIBEA and SMS-EMOA may be a consequence of using the hypervolume as an indicator in optimisation when faces are encountered that are (near-)parallel to at least one objective's axis.
In the future, we will explore ways to systematically examine the reasons for the gaps and improve mIBEA by redirecting its attention away from extreme solutions, without introducing arbitrary decisions, magic numbers or explicitly defined weights or preferences in the objective space.
) where I H (A) denotes the hypervolume formed by the solution set A. Correspondingly, I H (A + B) means the hypervolume of the union of solution set A and B. I HD (A, B) is negative if all solutions in B are dominated by solutions in A. Note that I H (A) = I H (B).

Fig. 2 :
Fig. 2: Visualisation of the non-dominated solutions found by mIBEA on DTLZ1 and DTLZ3 over 100 independent runs.

ACKNOWLEDGEMENTFig. 5 :
Fig. 5: PFs of DTLZ1 and DTLZ3 generated by IBEA using population size 1000 with different ρ values.100 runs are performed for each experiment.The threshold of border fraction (bf) for DTLZ1 is 0.03 and DTLZ3 is 0.1.

TABLE II :
Performance Comparison of NSGA-II and mIBEA