Birmingham Stochastic Ranking Algorithm for Many-Objective Optimization Based on Multiple Indicators

—Traditional multiobjective evolutionary algorithms face a great challenge when dealing with many objectives. This is due to a high proportion of nondominated solutions in the population and low selection pressure toward the Pareto front. In order to tackle this issue, a series of indicator-based algorithms have been proposed to guide the search process toward the Pareto front. However, a single indicator might be biased and lead the population to converge to a subregion of the Pareto front. In this paper, a multi-indicator-based algorithm is proposed for many-objective optimization problems. The proposed algorithm, namely stochastic ranking-based multi-indicator Algorithm (SRA), adopts the stochastic ranking technique to balance the search biases of different indicators. Empirical studies on a large number (39 in total) of problem instances from two well-deﬁned benchmark sets with 5, 10, and 15 objectives demonstrate that SRA performs well in terms of inverted generational distance and hypervolume metrics when compared with state-of-the-art algorithms. Empirical studies also reveal that, in the case a problem requires the algorithm to have strong convergence ability, the performance of SRA can be further improved by incorporating a direction-based archive to store well-converged solutions and maintain diversity.

[Pareto front (PF)] of different objectives.MOPs with more than three objectives are often called many-objective optimization problems (MaOPs) [1].MaOPs appear in various real-world applications such as car controller optimization [2], software engineering [3], and water supply portfolio planning [4].
In recent decades, many multiobjective evolutionary algorithms (MOEAs) have been proposed.However, traditional dominance-based MOEAs, such as NSGA-II [5] and SPEA2 [6], have been shown to be inefficient when dealing with MaOPs [7], [8].With an increasing number of objectives in the problem, the proportion of nondominated solutions is quite large and the traditional Pareto dominance loses its efficiency to push the population toward the PF.This is referred to as dominance resistance phenomenon [9].When the dominance-based (primary) selection criterion fails to differentiate the nondominated solutions, the diversity-based (secondary) criterion plays a vital role during environmental selection.Thus, the final population may spreads all over the objective space but fails to converge to the PF.
In order to overcome this obstacle, researchers have proposed various many-objecitve evolutionary algorithms (MaOEAs) in the literature.Based on the key ideas used, these methods can be categorized into the six classes [10].
1) Relaxed dominance-based algorithms try to alleviate the inefficiency of dominance by enlarging the dominated area of a solution.A series of approaches have been proposed, e.g., -dominance [11], controlling dominance area of solution [12], and L-dominance [13].Under these relaxed definitions, a solution has a higher chance to be dominated by other solutions and thus the selection pressure toward the PF is increased.One representative relaxed dominance-based algorithm is the GrEA [14] which uses grid-based convergence and diversity measurements to compare nondominated solutions.A difficult issue of these methods is to determine the extent of relaxation of the new dominance definitions for different problems and dynamically tuning method has been studied [15].2) Diversity-based methods try to improve the performance through more advanced diversity maintaining strategy.In general, by considering convergence to a certain amount, these methods aim to reduce the detrimental impact of diversity maintenance on the selection pressure.For example, diversity management mechanism DM1 [16] deactivates the diversity promotion mechanism once the This work is licensed under a Creative Commons Attribution 3.0 License.For more information, see http://creativecommons.org/licenses/by/3.0/population is excessively diverse.Shift-based density estimation (SDE) strategy [17] shifts the position of solutions according to their convergence information during the density estimation and it shows competitive performance when integrated into SPEA2 [17], [18].However, one might wonder what is the cost of these strategies with the gain of selection pressure (for example, SPEA2+SDE tends to concentrate on the central area of PF when tackling DTLZ1 problems [19]).
3) The aggregation-based methods use a series of scalarizing functions to decompose the MaOPs into a set of single objective subproblems.Various scalarizing functions and weight vector generation methods have been studied in [20] and [21].MOEA/D [20] is a popular aggregation-based method.The main advantage of MOEA/D over other aggregation-based methods lies in that it incorporates the neighborhood of subproblems to improve the efficiency of both generating new solutions and selecting solutions for the next generation.However, it maintains relatively poor diversity for high-dimensional problems and more advanced solutionvector mapping method has been proposed [22].Besides, although the contour lines of different scalarizing functions have been studied [23], [24], more in-depth studies on the scalarizing functions for different problems are still needed.4) Indicator-based methods take advantage of indicator values to guide the search process when optimizing an MaOP.Various indicators have been used to design MaOEAs, such as hypervolume [25], generational distance (GD) [26], and inverted generational distance (IGD) [27], and so on.The final solution set depends mainly on the characteristics of the indicators incorporated.However, hypervolume-based algorithms (evolutionary multi-objective optimization algorithms based S metric selection [28] and HypE [29]) are usually more time-consuming compared with other algorithms.Evaluating GD and IGD values needs a reference set to serve as the PF.Since the true PF is usually unknown a priori, maintaining the reference set is also a nontrivial task.5) Preference-based algorithms focus on the region of interest according to the user's preference information.In order to select solutions, a series of preference models have been studied in the literature such as goal specification, preference polyhedron, objective weighting, and so on [30], yet how to choose the appropriate preference model may be a problem-dependent task.An interesting algorithm, namely preference-inspired coevolutionary algorithms, model preference information as a set of solutions which coevolve along with the population [31], [32].6) More recently, hybrid MaOEAs that combine two or more of the abovementioned techniques have also been proposed.Two representative algorithms are NSGA-III [21] and Two_Arch2 [33].NSGA-III is based on Pareto and aggregation where the Pareto dominance-based nondomination sorting is used to drive the population toward the PF and a set of reference directions are used to maintain the diversity of the population.1 Two_Arch2 maintains two archives, i.e., the convergence archive (CA) and the diversity archive (DA), to aim at convergence and diversity, respectively.The additive indicator I + is used to update CA while Pareto dominance and a L p -norm based nearest neighbor distance is used to update DA.However, Two_Arch2 might fail to preserve the extreme points of the PF, while NSGA-III may struggle to converge to the PF on multimodal problems, as can be seen from our experimental results.In order to tackle an MaOP, the indicator-based algorithms seem to be quite straightforward since the final solution set is evaluated according to the indicators.However, a single indicator might bias the search toward a certain subregion.For example, indicator-based evolutionary algorithm (IBEA), an algorithm based on the I + indicator, struggles to maintain the diversity for solution set when dealing with many-objective problems [33].This phenomenon indicates that the I + indicator prefers convergence to diversity.Other indicators (e.g., the crowding distance [5], I SDE [17], etc.) may prefer diverse solutions instead.Since indicators may have different biases, which might complement each other, using multiple indicators rather than a single one for environmental selection may result in an even better algorithm.Motivated by this, a multiindicator algorithm is proposed in this paper.In particular, the key question to develop such an algorithm is how to carry out environmental selection based on multiple indicators that are inconsistent with each other.We employ the stochastic ranking technique, which was originally designed to balance the fitness and constraint violation in constrained optimization, to address this difficulty and show that the resultant algorithm, namely stochastic ranking-based multi-indicator algorithm (SRA), can achieve competitive performance on a large variety of test problems.
The main contributions of this paper are summarized as follows.
1) We introduce a new technique for balancing the influence of different indicators in a many-objective optimization algorithm.The balance is achieved through a ranking procedure based on the stochastic bubblesort.With a tunable parameter, it enables an algorithm designer to make an appropriate tradeoff between different indicators.2 2) In order to deal with problems that require the algorithm to have strong convergence ability, a directionbased archive (DBA) is incorporated into SRA to store well-converged solutions and maintain diversity.3) Comprehensive experimental studies covering both DTLZ and WFG problems reveal that different benchmark sets favor different algorithms and show the necessity of using both test suites to evaluate the performance of many-objective evolutionary algorithms.
The rest of this paper is organized as follows.Section II introduces the preliminary background of many-objective optimization.Section III presents the details of the SRA.Section IV illustrates SRA with archive.Section V is devoted to experimental setup.Section VI discusses the empirical results.Section VII concludes this paper and indicates some future directions.

A. MultiObjective Optimization Problem
Generally, an MOP can be stated as follows [34]: where x = (x 1 , x 2 , . . ., x n ) is the decision vector, is the (nonempty) decision space, F : → is the objective function vector, is the objective space.Here, we refer to MOPs as problems whose number of objectives m is larger than 1 and MaOPs as problems whose m is larger than 3.

C. Pareto Optimal Solution
A solution x * ∈ f is said to be Pareto optimal if no other solution x ∈ f can dominate it.

D. Pareto Set
The solution set consists of all the Pareto optimal solutions is called the Pareto Set: PS = {x ∈ f |∀y ∈ f , y ⊀ x}.

E. Pareto Front
The corresponding objective vector set of the Pareto Set is called the PF.

F. Approximation Set
An approximation set A ⊂ is defined as a solution set where no member of A dominates or is equal to any other member in A [36].
The goal of optimizing an MaOP is to obtain an approximation set A considering the following two subgoals [37].
1) Convergence: A are as close as possible to the PF.
2) Diversity: A spreads as diversely as possible.

G. Quality Indicator
An k-ary quality indicator I is a function I : k → R, which assigns each vector (A 1 , A 2 , . . ., A k ) of k approximation sets a real value I(A 1 , A 2 , . . ., A k ) [36].Quality indicators quantify the goodness of a solution set in terms of convergence, diversity, or both.

A. Overview
In the proposed SRA, a stochastic ranking procedure is employed to carry out environmental selection based on multiple indicators.The framework of SRA is described in Algorithm 1. First, N solutions are randomly created as the initial population.At each generation, randomly picked parent individuals are used to create offspring.After fitness evaluation, the offspring population is merged with the parent population.Then the indicator values of the merged population are computed.After that, the stochastic ranking based procedure is implemented for environmental selection.When the iterative optimization is finished, the nondominated solutions in the final population is returned as the output.

B. Indicators
Computationally, Algorithm 1 can accommodate any computable indicators (i.e., the I 1 and I 2 in line 11).However, since we expect to acquire additional benefits by involving two indicators, intuitively they should show different biases, e.g., one favors convergence and the other prefers diversity.Hence, the indicators I + [38] and I SDE [17] are chosen for SRA since: 1) these two indicators have been shown to be effective in terms of convergence or diversity and 2) they do not involve the nontrivial task of setting appropriate reference set/point.
The additive indicator I + and the corresponding I 1 (x) for comparing solutions are defined as where P is the population3 that includes x and y.
The I SDE and the corresponding I 2 (x) for comparing solutions are defined as (4) where Here, y precedes x means that the original index of the position of y in the population P is smaller than x.It should be noted that ( 5) is different from the original definition in [17], which computes ( 4) for all pairs of x and y.The reason is that adopting multiple indicators introducing additional computational costs, i.e., two indicators (instead of one) need to be computed.The computational overhead can be reduced by using ( 5) instead of the original definition.Meanwhile, this modification did not deteriorate the performance of SRA according to our preliminary experimental studies.

C. Stochastic Ranking-Based Environmental Selection
Suppose N individuals need to be selected from 2N individuals.Intuitively, the environmental selection can be viewed as comprising two steps.That is, the 2N individuals are first sorted according to a criterion, and then the best N are selected.The sorting process is trivial when only one indicator is employed, since it only takes the indicator value of each individual as the input.In case of multiple indicators, however, sorting becomes much more complicated, because different indicators might assign different ranks to the same individual.In fact, this is likely to happen especially when the indicators favor different aspect of multiobjective optimization, i.e., convergence and diversity.The sorting result in such a case should represent good balance between two indicators that produce inconsistent ranks to the same population.This scenario is similar to that commonly encountered in constrained (single objective) optimization, in which one may get quite different results by sorting the population according to fitness and constraint violation.Therefore, SRA utilizes the stochastic ranking technique [39], an approach that has been shown to be effective for the sorting problem in constrained optimization, to address the sorting problem in multiple indicator MaOEAs.
To be specific, stochastic ranking is applied once the indicator values of all individuals are obtained (line 13 in Algorithm 1).As shown in Algorithm 2, stochastic ranking is a stochastic bubble-sort algorithm.It ranks the individuals by sweeping the whole population (of size 2N) for N times. 4During each sweep, all adjacent individuals are compared according to values of a randomly chosen indicator.The sweep stops if there is no change in the rank ordering.The tradeoff of different indicators is controlled by the parameter p c ∈ (0, 1).After the ranking procedure terminates, the top N individuals will be selected. 4

D. Computational Complexity Analysis
Given an MaOP with m objectives and a population size of N, the time complexity of each generation of SRA is as follows: first, the time complexity of creating a new population (lines 5-8 in Algorithm 1) is O(N).Second, the fitness evaluation and population merging steps (lines 9 and 10 in Algorithm 1) need O(mN) and O(N), respectively.Third, the time complexity of computing I 1 (u i ) and I 2 (u i ) for all u i ∈ U t is O(mN 2 ).Fourth, the stochastic ranking-based environmental selection procedure (line 13 in Algorithm 1) takes O(N 2 ) time complexity.In summary, the computational complexity is O(mN 2 ), which is at the same level as the complexity of NSGA-III and Two_Arch2 [33].

IV. SRA WITH ARCHIVE
The SRA presented in Section III randomly chooses two individuals from a population to generate an offspring.Meanwhile, some other schemes, such as maintaining an external archive have shown their advantages in a number of existing MaOEAs [33], [40].Thus, a variant of SRA, namely SRA with Archive (SRA2) is further developed by incorporating a DBA into SRA.Here, the search directions are defined by a set of weight vectors, 5 which are set according to the method in [21].Given a set of weight vectors, SRA2 first generates 2N individuals, N of which are randomly assigned to each weight vector and are treated as the initial archive.The other N individuals are adopted as the initial population.In each generation, each member in the archive is recombined with an individual randomly picked from the current population to generate an offspring.Then, the archive is updated after the environmental selection step using the new population.
The pseudocode of the archive update procedure is shown in Algorithm 3.For weight vector-based methods, two key aspects are the association step and the replacement step.In SRA2, an individual is associated to the weight vector with the minimum perpendicular distance according to [21] and [22].An illustration of the association step is shown in Fig. 1.After that, an individual can replace at most n r archive members corresponding to the neighbors of the associated weight vector, as long as it obtains better fitness value than the archive members along the directions [20].The penalty-based boundary intersection (PBI) [20] fitness function, defined as a weighted sum of the perpendicular distance and the projection length, is adopted in SRA2 to compare an individual with the archive member.
In terms of computational complexity, the main additional cost induced by the use of archive lies in the archive updating phase.For each of the N solutions, the complexity of associating the solution to a weight vector and of updating the archive members is O(Nm) and O(Tm), respectively.Since T N, the complexity for the archive updating phase is O(mN 2 ), In other words, SRA and SRA2 have the same computational complexity, while the latter is more time-consuming than the former in practice.

A. Test Problems
In order to evaluate the performance of our algorithm, we have tested it on two well-defined benchmark problem sets: DTLZ [42] and WFG [43] test suites. 6Specifically, DTLZ1-DTLZ4 and WFG1-WFG9 with 5, 10, and 15 objectives are used for empirical studies.The parameter settings and the characteristics of the problems are listed in Table I. 6 For DTLZ1, we use the normalized form as in [20] and the PF is the hyper-plane satisfying that m i=1 f i = 1.

B. Performance Metrics
In the experiment, we use two indicators, the hypervolume indicator (HV) [25] and IGD [27].Both of them are widely used in [21] and [33].
IGD [27] is a distance-based metric, which is defined as where A = {a 1 , a 2 , . . ., a |A| } is an approximation set of the PF, PF = {p 1 , p 2 , . . ., p |PF | } is a subset of the true PF.The distance between a solution and a solution set is defined as the minimum distance between the solution and all solution set members.Here p is set to 1. Since the PF is known a priori for both DTLZ1-4 and WFG1-9 problems, we sample 500 000 points on the Pareto surface as PF according to [33].The IGD metric can measure the convergence and diversity of an solution set simultaneously.The smaller the IGD value is, the better the solution set is.
The HV measures the volume of solutions that is dominated by the approximation set [25], [28].Given a reference point z † and an approximation set A, HV(z † , A) is defined as where L is the Lebesgue measure, z † is the reference point in the objective space.It measures both convergence and diversity of an approximation set.HV is the only indicator which is Pareto-compliant [29].A higher HV value indicates a larger volume is dominated by the approximation set, thus corresponding to better performance.In our experiments, the output population is first normalized using the 1.1 × z nadir .After that, the hypervolume is computed using (1.0, . . ., 1.0) T as the reference point.In order to estimate the true HV value, we use a method from [44] according to [45].

C. State-of-the-Art Algorithms
In order to verify the performance of the proposed algorithm, the following state-of-the-art many-objective evolutionary algorithms (implemented in the jMetal framwork [46] on a 32-core 2.60 GHz Intel Xeon CPU with 64 GB RAM) are considered.
1) HypE [29]: HypE is an indicator-based many-objective algorithm.It uses Monte Carlo sampling to approximate the exact Hypervolume values.This enables the user to make a tradeoff between accuracy and computation time.Moreover, solutions are compared with the expected hypervolume loss for environmental selection as well as mating selection, aiming to maximizing the hypervolume of the solution set.2) AGE-II [47]: AGE-II is an approximation indicatorbased algorithm [48].It adopts -dominance to store all the -nondominated solutions into an archive. 7During environmental selection, the solution set that approximates the archive best is kept for the next generation.3) NSGA-III [21]: NSGA-III use a set of weight vectors to indicate different search directions.In order to maintain diversity, it associate solutions to weight vectors and compares niche counts of different weight vectors.
The algorithm has been shown to be highly competitive according to [21] and [33].4) Two_Arch2 [33]: Two_Arch2 is a newly proposed algorithm that uses two archives to take care of convergence and diversity, respectively.The CA is updated according to the additive indicator I + , while the DA is updated according to the L p -norm-based nearest neighbor distance.[20]: MOEA/D is a weight vector-based algorithm that decomposes an MaOP into a set of single-objective subproblems.The convergence of the algorithm is improved by the neighborhood definition and the diversity is maintained by the search directions represented by the weight vectors.7) IBEA [38]: IBEA is based on the I + indicator.It has been shown that it may converge quickly but leads to a solution set with poor diversity, which indicates that it might cause biased search on some problems.

D. Parameter Settings
The general and algorithm-specific parameter settings are summarized as follows.
1) Population Size: The population size of representative algorithms is shown in Table II.For the other algorithms, the population size is the same with that of NSGA-III.2) Reproduction Operators: The simulated binary crossover operator and polynomial mutation are used for reproducing offspring solutions [49].The mutation probability is set to 0.1, the crossover probability is set to 1.0.The mutation distribution index η m and the crossover distribution index η c are set to = 15.0 [33].
3) The neighborhood size is set to 20 and the maximum replacement number is set to 2 for MOEA/D and SRA2.
4) The in AGE-II is set to 0.1 [45].5) Termination Criterion: All algorithms are allowed for a maximum of 90 000 fitness evaluations for all the problem instances.6) Number of Runs: For all the problems instances, all the algorithms are repeated 20 times independently.7) Statistical Test: In order to test the difference of algorithms, the Wilcoxon rank sum test [50] (0.05 significance level) is applied for analysis.8) The value of parameter p c [p c lies in (0, 1)] needs to be set for SRA and SRA2.This parameter is an inherent parameter of stochastic ranking and controls the balance between indicators.The optimal value for p c is problem-dependent and may vary as the search progresses.To make a fair comparison, p c was not fine-tuned for each test instance.Instead, preliminary experiment was conducted on DTLZ1 with 5, 10, and 15 objectives.The performance of SRA and SRA2 with p c setting to 0.3 to 0.7 with step size 0.1 were compared.It was found that p c ∈ [0.4,0.6] works well for both algorithms.Hence, p c was set to a random number generated in this range at each generation of SRA and SRA2.

A. Bias of Indicators
In this experiment, we choose the ten objective DTLZ1 problem for investigation [33].In order to show the search bias of indicators, we implement three variants of stochastic ranking with the parameter p c set to 1, 0.5, and 0 (corresponding to I + only, I + and I SDE , I SDE only scenarios).The direction-based archive is not included in the three variants in order to keep the algorithm behavior away from the archive's impact.Fig. 2 shows a typical solution set in parallel coordinates obtained by the three algorithm variants.As can be seen from the figure, I + leads the population to convergence to a single solution in the objective space; I SDE results in a diverse population yet it is unable to push the population to converge to the PF; the combination of the two indicators with a probability of 0.5 results in a well-converged and well-spreading population.

B. Effectiveness of Archive
In order to show the effectiveness of the archive, we compared SRA and SRA2 (SRA with archive) on DTLZ1-4 problems with 5, 10, and 15 objectives.The IGD and HV results are shown in Table III.From the table, we can see that the overall performance of SRA2 is better than SRA in terms of both IGD and hypervolume, since the right column corresponding to SRA2 contains more gray items than the left one.For multimodal problems DTLZ1 and DTLZ3, SRA2 performs better than SRA in terms of IGD.The major credit goes to the direction based archive, which helps to tackle the multimodal issue and improve the convergence of the algorithm.For unimodal problems DTLZ2 and DTLZ4, the contradiction between different indicators seems to be less than that on multimodal problems.SRA2 is better than SRA in terms of hypervolume, which indicates that the archive helps to maintain diversity for the population.

C. Comparison With State-of the-Art Algorithms
The performance comparisons of the peer algorithms on DTLZ and WFG problems in terms of IGD and HV are shown in Tables IV-VII, respectively.Table VIII presents the pair-wise comparison results according to Wilcoxon rank sum test.Based on the statistical test results, the average performance score (APS) [29] is further calculated to reveal the overall performance of the algorithms.Given a certain problem instance, the performance score of an algorithm is defined as the number of algorithms that beat the algorithm on the specific instance significantly.The smaller the performance score is, the better the algorithm performs on that instance.Here, the APS over all test instances for all these algorithms in terms of IGD and HV is summarized in Fig. 3.
From the results of Wilcoxon rank sum test and APS, it can be observed that both SRA and SRA2 achieve competitive performance since they are both among the top three algorithms in terms of the overall IGD (rank 2nd and 3rd, respectively) as well as the overall HV (rank 1st and 2nd, respectively).Moreover, one should note that the best algorithms in terms of the overall IGD and HV are two different algorithms.These results demonstrate that the proposed SRA and SRA2 are able to achieve good performance in terms of different performance metrics through adopting the stochastic ranking strategy with multiple indicators.
The performance of algorithms not only vary over performance metrics, but also vary over problem classes.From the statistical test table and Fig. 3, one can see that our algorithms generally perform better on WFG problems than on DTLZ problems.This is because that the DTLZ test suite requires an algorithm to have strong convergence abilities but does not require strong diversification properties [51].In SRA and SRA2, on the other hand, both indicators (emphasizing convergence and diversity, respectively) are given the same importance.We can also observe that SRA2 achieves better performance than SRA on the DTLZ test suite, while this is not the case for WFG problems.The direction-based archive in SRA2 seems to be too greedy for WFG problems in terms of convergence over diversity.This also verifies the reason why SRA does not perform as well on DTLZ problems to some extent.1) DTLZ Test Suite: A closer look of the results in Tables IV and V can tell us that other algorithms achieve mixed performance on DTLZ problems.To be specific, according to Table IV, AGE-II and SPEA2+SDE are the best two algorithms in terms of IGD, yet the order of these two algorithms may vary with different objective numbers.Similar results are obtained for MOEA/D and NSGA-III in terms of HV on Table V.Thus, we show the parallel coordinates of the solution sets on 10-objective problems with the best IGD value among all the 20 runs in Figs.4-7 8 and discuss the results on DTLZ problems in detail.
For the DTLZ1 problem, SPEA2+SDE and NSGA-III achieve the best IGD and hypervolume performance, respectively.From Fig. 4, we can see that HypE, AGE-II, and IBEA lead the population to sub-regions of the PF in different forms, which demonstrates the search bias of indicators to some    DTLZ3 includes a great number of local optima, and imposes a great challenge for the algorithms to push the population toward the PF.For the DTLZ3 problem, the best IGD performance is achieved by MOEA/D and AGE-II while the best hypervolume value is achieved by SRA.From Fig. 6, one can see that HypE and IBEA still have issues with maintaining diversity.Both NSGA-III and Two_Arch2 have trouble converging within the given number of fitness evaluations.
One reason for the poor convergence of NSGA-III is that the Pareto dominance-based nondomination sorting loses its efficiency in high dimensional problems.The remaining algorithms show good convergence and diversity on this test problem.
For the DTLZ4 problem, the main challenge, maintaining diversity of the population, is created by a nonuniform density of solutions.AGE-II and MOEA/D perform best in terms of IGD and HV, respectively.Since this is a unimodal problem, it does not cause much trouble for the convergence of the tested algorithms.From Fig. 7, we can see that HypE performs poorly in terms of maintaining diversity.Two_Arch2 has trouble to obtain the boundary solutions on the PF.  2) WFG Test Suite: The results on WFG problems look more clear, so relatively less discussion is devoted to the WFG test suite.Our SRA and SRA2 performs well on WFG problems in terms of both IGD and HV.In fact, SRA ranks second in terms of overall IGD and first in terms of the overall hypervolume.The other algorithm which also rank 1st in terms of HV, i.e., HypE, obtains the second worst IGD.Table VIII also shows similar results of SRA and HypE, where the pairwise win-tie-loss counts of HypE against SRA in terms of IGD and HV are 0-0-27 and 13-4-10, respectively.This is another evidence that shows different indicators may have different biases.Not surprisingly, HypE performs well in terms of hypervolume since it is guided by the HV.Yet this comes with a loss in performance in terms of the IGD metric.
In order to reveal more details of the population of the algorithms, the results on the ten objective WFG4 problem is shown by parallel coordinates in Fig. 8. From the figure, one can see that the solution sets of HypE and IBEA are rather poorly distributed, while the proposed SRA provides a more evenly distributed population over different objective dimensions.This is because that emphasizing a single indicator may cause performance deterioration in terms of other indicators.

VII. CONCLUSION
In this paper, a multi-indicator-based algorithm is proposed for many-objective optimization problems.The proposed algorithm, namely SRA, employs the stochastic ranking technique to balance the search biases of different indicators.Moreover, in the case a problem requires the algorithm to have strong convergence ability, the performance of SRA can be further improved by incorporating an direction-based archive, which preserves well-converged solutions and maintains diversity.Experimental results on 39 test instances with 5, 10, and 15 objectives from two benchmark problem suites indicate that the proposed algorithm is competitive compared with stateof-the-art algorithms in terms of both IGD and hypervolume.As a side note, we could like to encourage all researchers to carry out experimental studies on both DTLZ and WFG problems because they have very different characteristics.Using only one of them could potentially bias conclusions drawn from the experiments.
Although the overall performance of SRA is very promising, more studies need to be carried out in the future.First, a simple parameter setting scheme is adopted in the current version of SRA.Yet the effect of the parameters on SRA's performance should be thoroughly studied and some adaptive parameter strategies can be proposed accordingly.Second, only one pair of indicators is studied in this paper.It should be interesting to examine the behavior of SRA with other indicators.Third, SRA is outperformed by MOEA/D on the DTLZ test suite.The main reason for this needs to be studied.Fourth, the performance of SRA on real-world applications and combinatorial problems with more objectives needs to be tested.

Fig. 1 .
Fig. 1.Associate solutions to weight vectors (denoted as V1, . . ., V5).All the solutions (denoted as S1, . . ., S5) that survive the stochastic ranking-based multi-indicator environmental selection in the previous step are associated with the weight vector with the minimum perpendicular distance.

Fig. 2 .
Fig. 2. Solution set obtained by three variants of stochastic ranking with the parameter p c set to 0 (left), 0.5 (middle), and 1 (right), respectively, (corresponding to I + only, I + and I SDE , I SDE only scenarios) on 10-objective DTLZ1 problem.

Fig. 3 .
Fig. 3. APS for all the compared algorithms on the instances from (a) DTLZ (left), (b) WFG (middle), and (c) both test suites (right).The smaller the APS is, the better the algorithm performs.

Fig. 4 .
Fig. 4. Solution set corresponding to the best IGD values for all algorithms on DTLZ1 problems with ten objectives.The x-axis is the objective number, while the y-axis is the objective value.

Fig. 5 .
Fig. 5. Solution set corresponding to the best IGD values for all algorithms on DTLZ2 problems with ten objectives.The x-axis is the objective number, while the y-axis is the objective value.

Fig. 6 .
Fig.6.Solution set corresponding to the best IGD values for all algorithms on DTLZ3 problems with ten objectives.The x-axis is the objective number, while the y-axis is the objective value.

Fig. 7 .
Fig. 7. Solution set corresponding to the best IGD values for all algorithms on DTLZ4 problems with ten objectives.The x-axis is the objective number, while the y-axis is the objective value.

Fig. 8 .
Fig.8.Solution set corresponding to the best IGD values for all algorithms on WFG4 problems with ten objectives.The x-axis is the objective number, while the y-axis is the objective value.

Algorithm 1 :
Main Loop of SRA input : an MaOP, population size N output: an approximation set A out 1 Randomly generate the initial population P 0 2 Evaluate all the individuals in P0 3 Set Q 0 ← Ø, t ← 0Create an offspring q i with randomly chosen p 1 ,p 2 ∈ P t 4 while t < MaxGen do 5 for each i ∈ N do 6

7
Set Q t ← Q t ∪ {q i } Obtain the combined population: U t ← P t ∪ Q t 11Compute the indicator values I 1 (u i ) and I 2 (u i ) for all u i ∈ U t Return the non-dominated solutions of A t+1 as A out

11 end 12 end 13 end 14 if no swap done then 15 break; 16 end 17 end 18
Copy top |U t |/2 solutions of U t to P t+1

TABLE III PERFORMANCE
COMPARISON OF SRA AND SRA2 IN TERMS OF THE AVERAGE IGD AND HV VALUES ON DTLZ PROBLEMS.THE BETTER AVERAGE VALUES OF THE TWO ALGORITHMS FOR EACH INSTANCE ARE HIGHLIGHTED IN GRAY

TABLE IV PERFORMANCE
COMPARISON OF SRA AND SRA2 WITH STATE-OF-THE-ART ALGORITHMS IN TERMS OF THE AVERAGE IGD VALUES ON DTLZ PROBLEMS.THE BEST AND SECOND BEST AVERAGE IGD VALUES AMONG ALL THE ALGORITHMS FOR EACH INSTANCE ARE HIGHLIGHTED IN GRAY AND LIGHT GRAY

TABLE V PERFORMANCE
COMPARISON OF SRA AND SRA2 WITH STATE-OF-THE-ART ALGORITHMS IN TERMS OF THE AVERAGE HV VALUES ON DTLZ PROBLEMS.THE BEST AND SECOND BEST AVERAGE HV VALUES AMONG ALL THE ALGORITHMS FOR EACH INSTANCE ARE HIGHLIGHTED IN GRAY AND LIGHT GRAY TABLE VI PERFORMANCE COMPARISON OF SRA AND SRA2 WITH STATE-OF-THE-ART ALGORITHMS IN TERMS OF THE AVERAGE IGD VALUES ON WFG PROBLEMS.THE BEST AND SECOND BEST AVERAGE IGD VALUES AMONG ALL THE ALGORITHMS FOR EACH INSTANCE ARE HIGHLIGHTED IN GRAY AND LIGHT GRAY

TABLE VII PERFORMANCE
COMPARISON OF SRA AND SRA2 WITH STATE-OF-THE-ART ALGORITHMS IN TERMS OF THE AVERAGE HV VALUES ON WFG PROBLEMS.THE BEST AND SECOND BEST AVERAGE HV VALUES AMONG ALL THE ALGORITHMS FOR EACH INSTANCE ARE HIGHLIGHTED IN GRAY AND LIGHT GRAY

TABLE VIII SUMMARY
OF THE WILCOXON TEST IN TERMS OF PERFORMANCE METRICS BETWEEN DIFFERENT ALGORITHMS ON DTLZ AND WFG PROBLEMS.THE PAIRWISE WIN-TIE-LOSS COUNTS OF ROWS AGAINST COLUMNS ARE SHOWN IN THE TABLE