Many real-life disciplines (e.g., optimal design , economics  and water management ) often involve optimization problems having multiple conflicting objectives, known as multi-objective optimization problems (MOPs). Rather than a global optimum that optimizes all objectives simultaneously, in multi-objective optimization, decision makers often look for a set of Pareto-optimal solutions which consist of the best trade-offs among conflicting objectives. The balance between convergence and diversity is the cornerstone of multi-objective optimization. In particular, the convergence means the approximation to the Pareto-optimal set should be as close as possible; while the diversity means the spread of the trade-off solutions should be as uniform as possible.
Evolutionary algorithm, which in principle can approximate the Pareto-optimal set in a single run, has been widely accepted as a major approach for multi-objective optimization . Over the past three decades, many efforts have been devoted in developing evolutionary multi-objective optimization (EMO) algorithms and have obtained recognized performance on problems with two or three objectives [5, 6, 7, 8, 9, 10, 11, 12]. With the development of science and technology, real-life optimization scenarios bring more significant challenges, e.g., complicated landscape, multi-modality and high dimensionality, to the algorithm design. As reported in [13, 14, 15], the performance of traditional EMO algorithms severely deteriorate with the increase of the number of objectives. Generally speaking, researchers owe the performance deterioration to three major reasons, i.e., the loss of selection pressure for Pareto domination 
, the difficulty of density estimation in a high-dimensional space and its anti-convergence phenomenon , and the exponentially increased computational complexity .
As a remedy to the loss of selection pressure for Pareto domination in a high-dimensional space, many modified dominance relationships have been developed to strengthen the comparability between solutions, e.g., -dominance , fuzzy Pareto-dominance , -optimality , preference order ranking , control of dominance area  and grid dominance . Very recently, a generalized Pareto-optimality was developed in  to expand the dominance area of Pareto domination, so that the percentage of non-dominated solutions in a population does not increase dramatically with the number of objectives. Different from , the expansion envelop for all solutions are kept the same in . Other than the Pareto dominance-based approaches, -optimality was proposed in  to help rank the solutions.
The loss of selection pressure can also be remedied by an effective diversity maintenance strategy. To relieve the anti-convergence phenomenon,  applied the maximum spread indicator developed in  to activate and deactivate the diversity promotion mechanism in NSGA-II . To facilitate the density estimation in a high-dimensional space,  proposed to replace the crowding distance used in NSGA-II by counting the number of associated solutions with regard to a predefined reference point. In particular, a solution is considered being associated with a reference point if it has the shortest perpendicular distance to this reference point. Instead of handling the convergence and the diversity separately,  proposed to co-evolve a set of target vectors, each of which represents a optimization goal. The fitness value of a solution is defined by aggregating the number of current goals it achieves, i.e., dominates. The population diversity is implicitly maintained by the goals widely spread in the objective space; while the selection pressure towards the convergence is strengthened by co-evolving the optimization goals with the population.
The exponentially increased computational costs come from two aspects: 1) the exponentially increased computational complexity for calculating the hypervolume, which significantly hinders the scalability of the indicator-based EMO algorithms; 2) the significantly increased computational costs for maintaining the non-domination levels  of a large number of solutions in a high-dimensional space. As for the former aspect, some improved methods, from the perspective of computational geometry [29, 30, 31] or Monte carlo sampling , have been proposed to enhance the efficiency of hypervolume calculation. As for the latter aspect, many efforts have been devoted to applying some advanced data structures to improve the efficiency of non-dominated sorting procedure [32, 33, 34]. It is worth noting that our recent study  showed that it can be more efficient to update the non-domination levels by leveraging the population structure than to sort the population from scratch in every iteration.
As reported in [36, 37, 38], decomposition-based EMO methods have become increasingly popular for solving problems with more than three objectives, which are often referred to as many-objective optimization problems (MaOPs). Since the decomposition-based EMO methods transform an MOP into several single-objective optimization subproblems, it does not suffer the loss of selection pressure of Pareto domination in a high-dimensional space. In addition, the update of the population relies on the comparison of the fitness values, thus the computational costs do not excessively increase with the dimensionality. As reported in , different subproblems, which focus on different regions in the objective space, have various difficulties. Some superior solutions can easily dominantly occupy several subproblems. This is obviously harmful to the population diversity and getting worse with the increase of the dimensionality. To overcome this issue,  built an interrelationship between subproblems and solutions, where each subproblem can only be updated by its related solutions. Based on the similar merit,  restricted a solution to only updating one of its closest subproblems. In , two metrics were proposed to measure the convergence and diversity separately. More specifically, the objective vector of a solution is at first projected onto its closest weight vector. Then, the distance between the projected point and the ideal point measures the solution’s convergence; while the distance between the projected point and the original objective vector measures the solution’s diversity. At the end, a diversity-first update scheme was developed according to these two metrics. Analogously,  developed a modified achievement scalarizing function as the convergence metric while an angle-based density estimation method was employed to measure the solution’s diversity.
Recently, there is a growing trend in leveraging the advantages of the decomposition- and Pareto-based methods within the same framework. For example,  suggested to use the Pareto domination to prescreen the population. The local density is estimated by the number of solutions associated with a pre-defined weight vector. In particular, a solution located in an isolated region has a higher chance to survive to the next iteration. Differently,  developed an angle-based method to estimate the local crowdedness around a weight vector. In addition to the density estimation, the weight vectors divide the objective space into different subspaces, which is finally used to estimate the local strength value  of each solution. Analogously, in , a non-dominated sorting is conducted upon all the subspaces, where solutions in different subspaces are considered non-dominated to each other.  developed a dual-population paradigm which co-evolves two populations simultaneously. These populations are maintained by different selection mechanisms respectively, while their interaction is implemented by a restricted mating selection mechanism.
Although the decomposition-based EMO methods have been extensively used for MaOPs, a recent comprehensive study  demonstrated that the performance of decomposition-based EMO methods strongly depends on the shape of the Pareto front (PF). This phenomenon can be attributed to two major reasons:
Most, if not all, decomposition-based EMO methods merely consider a single scalarizing function in fitness assignment. Since the contour line of a scalarizing function does not adapt to a particular problem’s characteristic, the flexibility is restricted.
As discussed in the previous paragraph, different regions of the PF have various difficulties. The balance between convergence and diversity of the search process can be easily biased by some dominantly superior solutions. The increasing dimensionality exacerbates this phenomenon.
Bearing the above considerations in mind, this paper develops a new decomposition-based method, called adversarial decomposition, for many-objective optimization. Generally speaking, it has the following three features:
It maintains two co-evolving populations simultaneously, each of which is maintained by a different scalarizing function. In particular, one population uses the boundary intersection-based scalarizing function, while the other one applies a modified achievement scalarizing function. In this regard, the search behaviors of these two populations become different, where one is convergence-oriented while the other is diversity-oriented.
In order to make these two populations complement each other, they use ideal and nadir points respectively as the reference point in their scalarizing functions. By doing this, the two populations are evolved following two sets of adversarial search directions.
During the mating selection process, two populations are at first stably matched to form a set of one-one solution pairs. In particular, solutions within the same pair concentrate on similar regions of the PF. Thereafter, each solution pair can at most contribute one mating parent for offspring generation. By doing this, we can expect an uniformly spread of the computational efforts over the entire PF.
The remainder of this paper is organized as follows. sec:preliminary provides some preliminaries useful to this paper. sec:proposal describes the technical details of our proposed method step by step. The empirical studies are presented and discussed in sec:setup and sec:experiments. Finally, sec:conclusion concludes the paper and provides some future directions.
This section first provides some basic definitions of multi-objective optimization. Afterwards, we briefly introduce the decomposition of multi-objective optimization.
2.1 Basic Definitions
The MOP considered in this paper can be mathematically defined as follows:
where is a -dimensional decision variable vector from the decision space and is a -dimensional objective vector in the objective space .
Given solutions , is said to dominate , denoted by , if and only if for all and .
A solution is called Pareto-optimal if and only if there is no other solution dominates it.
The Pareto-optimal set (PS) is defined as the set of all Pareto-optimal solutions. Their corresponding objective vectors form the Pareto-optimal front (PF).
The ideal point is , where for all . The nadir point is , where .
Under some mild conditions, the task of approximating the PF can be decomposed into several scalar optimization subproblems, each of which is formed by a weighted aggregation of all individual objectives . In the classic multi-objective optimization literature , there have been several established approaches for constructing scalarizing functions, among which the weighted Tchebycheff (TCH) and penalty-based boundary intersection (PBI)  are the most popular ones. More specifically, the TCH function is mathematically defined as:
where is a user specified weight vector, for all and . Note that is set to be a very small number, say , in case . The contour line of the TCH function is shown in fig:decomposition(a) where . From this figure, we can clearly see that the control area (i.e., the area that holds better solutions) of the TCH function is similar to the Pareto domination defined in def:pareto, e.g., solutions located in the gray shaded area (i.e., the control area of ) are better than . Note that the TCH function cannot discriminate the weakly dominated solution . For example, the TCH function values of and are the same, but according to def:pareto.
As for the PBI function, it is mathematically defined as:
As discussed in , and measure the convergence and diversity of with regard to , respectively. The balance between convergence and diversity is parameterized by , which also controls the contour line of the PBI function. In fig:decomposition(b), we present the contour lines of PBI functions with different settings.
3 Many-Objective Optimization Algorithm Based on Adversarial Decomposition
In this section, we introduce the technical details of our proposed many-objective optimization algorithm based on adversarial decomposition, denoted as MOEA/AD whose pseudo code is given in alg:moead-2p, step by step. At the beginning, we initialize a population of solutions via random sampling upon ; the ideal and nadir points; a set of weight vectors and build their neighborhood structure according to the method in . Afterwards, we assign directly to the two co-evolving populations, i.e., diversity population and convergence population . Note that and share the same weight vectors, each of which corresponds to a unique subproblem for and respectively. To facilitate the mating selection process, we initialize a matching array and a sentinel array , where indicates a solution in is paired with a solution in and indicates whether this pair of solutions work in similar regions of the PF. During the main while loop, the mating parents are selected from the solution pairs. The generated offspring solution is used to update and separately. After each generation, we divide the solutions in into different solution pairs for the next round’s mating selection process. The major components of MOEA/AD are explained in the following subsections.
3.1 Adversarial Decomposition
As discussed in sec:introduction, the flexibility of the decomposition-based method is restricted due to the use of a single scalarizing function. Bearing this consideration in mind, this paper develops an adversarial decomposition method. Its basic idea is to maintain two co-evolving and complementary populations simultaneously, each of which is maintained by a different scalarizing function.
More specifically, one population is maintained by the PBI function introduced in sec:decomposition, where we set as recommended in . The other population is maintained by a modified augmented achievement scalarizing function (MA-ASF) defined as follows:
where is an augmentation coefficient. Comparing with the TCH function in eq:TCH, the MA-ASF uses the nadir point to replace the ideal point and the absolute operator is removed to allow to be smaller than , where . Furthermore, the augmentation term in the MA-ASF helps avoid weakly Pareto-optimal solutions. As shown in fig:decomposition(c), the contour line of the MA-ASF is the same as that of the TCH function in case ; while the control area of the MA-ASF becomes wider when setting . In this case, the MA-ASF is able to discriminate the weakly dominated solution, e.g., the MA-ASF value of in fig:decomposition(c) is better than that of when . Here we use a sufficiently small as recommended in .
To deal with problems having different scales of objectives, we normalize the objective values before using the scalarizing function. By doing this, the PBI function becomes:
where and for all . The MA-ASF is re-written as:
In the following paragraphs, we will discuss the complementary effects achieved by the PBI function and MA-ASF.
As shown in fig:decomposition(b) and fig:decomposition(c), the control areas of the PBI function with is much narrower than that of the MA-ASF. In this case, the population is driven towards the corresponding weight vector by using the PBI function as the selection criteria. Moreover, comparing to the MA-ASF, the control areas shared by different weight vectors are smaller in the PBI function. Accordingly, various weight vectors have a larger chance to hold different elite solutions and we can expect an improvement on the population diversity. On the other hand, due to the narrower control area, the selection pressure, with regard to the convergence, of the PBI function is not as strong as the MA-ASF. In other words, some solutions, which can update the subproblem maintained by a MA-ASF, might not be able to update the subproblem maintained by a PBI function. For example, as shown in fig:decomposition(b), although , the PBI function value of is worse than that of with regard to . In this case, the PBI function has a high risk of compromising the population convergence.
The other reason, which results in the different behaviors of the PBI function and MA-ASF, is their adversarial search directions by using the ideal point and nadir point respectively. More specifically, the PBI function pushes solutions toward the ideal point as close as possible; while the MA-ASF pushes the solutions away from the nadir point as far as possible. Therefore, given the same set of weight vectors, the search regions of the PBI function and MA-ASF complement each other. For example, for a convex PF shown in fig:weights(a), solutions found by the PBI function using the ideal point concentrate on the central region of the PF; in contrast, those found by the MA-ASF using the nadir point fill the gap between the central region and the boundary. As for a concave PF shown in fig:weights(b), solutions found by the PBI function using the ideal point sparsely spread over the PF; while those found by the MA-ASF using the nadir point improve the population diversity.
In summary, by using the scalarizing functions introduced above, i.e., the PBI function and the MA-ASF, the adversarial decomposition method makes the two co-evolving populations become complementary, i.e., one is diversity-oriented (denoted as the diversity population ) and the other is convergence-oriented (denoted as the convergence population ). In addition, the search regions is also diversified so that the solutions are able to cover a wider range of the PF.
3.2 Population Update
After the generation of an offspring solution , it is used to update and separately. Note that the optimal solution for each subproblem is assumed to be along its corresponding reference line that connects the origin and the weight vector . Thus, to make as diversified as possible, we expect that different subproblems can hold different elite solutions. In this case, we restrict to only updating its closest subproblem (as shown in line 1 to line 5 of alg:update). As for , its major purpose is to propel the population to the PF. To accelerate the convergence progress, we allow a dominantly superior solution to take over more than one subproblem, say . In particular, we at first sort the priorities of different subproblems according to their distances to . It can update its closest subproblems in case has a better MA-ASF function value (as shown line 6 to line 15 of alg:update). It is worth noting that we reserve two additional terms, in line 13 of alg:update, to facilitate the mating selection procedure introduced in sec:mating. One is the degree of closeness of the updated solution to its corresponding subproblem, denoted as ; the other is the index of this solution’s closest subproblem, denoted as .
3.3 Mating Selection and Reproduction
The interaction between the two co-evolving populations is an essential step in algorithms that consider multiple populations [47, 51]. To take advantage of the complementary effects between and , the interaction is implemented as a restricted mating selection mechanism that chooses the mating parents according to their working regions. Generally speaking, it contains two consecutive steps: one is the pairing step that makes each solution in be paired with a solution in ; the other is the mating selection step that selects the appropriate parents for the offspring reproduction. We will illustrate them in detail as follows.
3.3.1 Pairing Step
To facilitate the latter mating selection step, the pairing step divides the two populations into different solution pairs, each of which contains two solutions from and respectively. This is achieved by finding a one-one stable matching between the solutions in and . As a result, solutions in the same pair are regarded to have a similar working regions of the PF.
To find a stable matching between solutions in and , we need to define their mutual preferences at first. Specifically, since each solution in is close to its corresponding subproblem, we define the preference of a solution in (denoted as ) to a solution in (denoted as ) as:
where is the weight vector of the subproblem occupied by . As for the preference of to , it is defined as:
Then, we sort the preference list of each solution in an ascending order and apply our recently developed two-level one-one stable matching method [52, 53] to find a stable matching. Note that the two-level stable matching method is able to match each agent with one of its most preferred partners. Since a solution of an -objective problem always locates within the local area between closest weight vectors, the length of the preference list is reduced to in the first-level matching process . As a result, the matched solutions in the first-level stable matching should work on the similar regions of the PF. For example, as shown in fig:match, the solution pairs formed in the first-level stable matching are surrounded by the red solid curves. From this figure, we can see that these matched solutions are close to each other and work on the similar regions. Therefore, we set the corresponding index of a sentinel array to 1, i.e., where and denote the corresponding subproblems have collaborative information. During the second-level matching process, the remaining solutions are matched based on the complete preference lists. Note that the matched solutions in the second-level stable matching are not guaranteed to work on the similar regions any longer. As shown in fig:match, the solution pair formed in second-level stable matching, surrounded by the red dotted curve, are away from each other. Thus, we set . The pseudo code of this pairing step is presented in alg:stm2l.
3.3.2 Mating Selection Step
The mating parents consist of two parts: one is the principal parent; the others are from its neighbors. Note that each solution pair is only allowed to contribute at most one mating parent to avoid wasting the computational resources upon the similar regions. Given a solution pair , the first step is to decide the population from which the principal parent is selected. This depends on the following three criteria.
The first criterion is the subproblem’s relative improvement. As for , it is defined as:
where and are respectively the current and previous solution of the th subproblem in , respectively. As for , it is defined as:
If and are in a tie, we need the following two secondary criteria for decision making.
As discussed in sec:ad, a solution found by the PBI function is not guaranteed to be non-dominated. If is dominated by a solution in , it is inappropriate to be chosen as a mating parent.
As discussed in sec:update, a dominantly superior solution can occupy more than one subproblem by considering the MA-ASF. In this case, such solution may occupy a subproblem away from its working region, which makes it inadequate to be a mating parent.
The pseudo code of the principal parent selection mechanism is given in alg:popSelect. If the subproblem’s relative improvement of is larger than that of , it means that the corresponding subproblem of has a higher potential for further improvement. And the principal parent should be selected from , i.e., . Otherwise, will be chosen as the principal parent (line 1 to line 4 of alg:popSelect). If and are in a tie, we need to check the domination status of and ’s closeness to the corresponding subproblem (line 6 to line 11 of alg:popSelect). By comparing the subproblems’ relative improvements, we can expect an efficient allocation of the computational resources to different regions of the PF. The two secondary criteria implicitly push the solutions towards the corresponding weight vectors thus improve the population diversity.
After the selection of the principal parent, the other mating parent are selected from the neighbors of the subproblem occupied by the principal parent. More specifically, if the principal parent is from , we store the solutions of its neighboring subproblems from both and into a temporary mating pool . Note that only those subproblems having collaborative information (i.e., ) are considered when choosing solutions from (line 5 and line 8 of alg:matingSelect). On the other hand, if the principal parent is from , only solutions from have the chance to be stored in . Note that we do not consider the solution that has the same closest subproblem as the principal parent (line 12 and line 14 of alg:matingSelect). Once is set up, the other mating parents are randomly chosen from it.
4 Experimental Setup
In this section, we describe the setup of our empirical studies, including the benchmark problems, performance indicator, peer algorithms and their parameter settings.
4.1 Benchmark Problems
Here we choose DTLZ1 to DTLZ4 , WFG1 to WFG9 , and their minus version proposed in , i.e., DTLZ1 to DTLZ4 and WFG1 to WFG9 to form the benchmark problems in our empirical studies. In particular, the number of objectives are set as . The number of decision variables of DTLZ and DTLZ problem instances  is set to , where for DTLZ1 and DTLZ1 and for the others. As for WFG and WFG problem instances , we set , where and . Note that the DTLZ and WFG benchmark problems have been widely used for benchmarking the performance of many-objective optimizers; while their minus version is proposed to investigate the resilience to the irregular PF shapes. All these benchmark problems are scalable to any number of objectives.
4.2 Performance Indicator
In our empirical studies, we choose the widely used Hypervolume (HV) indicator  to quantitatively evaluate the performance of a many-objective optimizer. Specifically, the HV indicator is calculated as:
where is the solution set, is a point dominated by all Pareto-optimal objective vectors and VOL indicates the Lebesgue measure. In our empirical studies, we set and the objective vectors are normalized to before calculating the HV. The larger the HV value is, the better the quality of is for approximating the true PF. Each algorithm is run 31 times independently and the Wilcoxon’s rank sum test at 5% significant level is performed to show whether the peer algorithm is significantly outperformed by our proposed MOEA/AD. Note that we choose the population that has the higher HV value as the output of our proposed MOEA/AD.
4.3 Peer Algorithms
Here we choose nine state-of-the-art many-objective optimizers to validate the competitiveness of our proposed MOEA/AD. These peer algorithms belong to different types, including two decomposition-based algorithms (MOEA/D  and Global WASF-GA ), two Pareto-based algorithms (PICEA-g  and VaEA ), two indicator-based algorithms (HypE  and KnEA ), two algorithms that integrates the decomposition- and Pareto-based selection together (NSGA-III  and -DEA ), and an improved two-archive-based algorithm (TwoArch2 ). Some further comments upon the peer algorithms are listed as follows.
Global WASF-GA is a decomposition-based algorithm that selects solutions to survive according to the rankings of solutions to each subproblem. Similar to MOEA/AD, it uses the ideal point and nadir point simultaneously in its scalarizing function. However, instead of maintaining two co-evolving populations, Global WASF-GA only has a single population, where half of it are maintained by the scalarizing function using the ideal point while the other are maintained by the nadir point.
PICEA-g co-evolves a set of target vectors sampled in the objective space, which can be regarded as a second population and is used to guide the environmental selection.
TwoArch2 maintains two archives via indicator-based selection and Pareto-based selection separately. In particular, an -norm-based diversity maintenance scheme is designed to maintain the diversity archive.
4.4 Parameter Settings
Weight vector: We employ the method developed in  to generate the weight vectors used in the MOEA/D variants. Note that we add an additional weight vector to remedy the missing the centroid on the simplex for the 8-, 10- and 15-objective cases.
Population size: We set the population size the same as the number of weight vectors. In particular, is set as 91, 210, 157, 276 and 136 for respectively.
Termination criteria: As suggested in , the termination criterion is set as a predefined number of generations, as shown in tab:evaluations.
Problem Instance =3 =5 =8 =10 =15 DTLZ1, DTLZ1 400 600 750 1,000 1,500 DTLZ2, DTLZ2 250 350 500 750 1,000 DTLZ3, DTLZ3 1,000 1,000 1,000 1,500 2,000 DTLZ4, DTLZ4 600 1,000 1,250 2,000 3,000 WFG, WFG 400 750 1,500 2,000 3,000 Table 1: Settings of the Number of Generations.
Neighborhood size: .
Probability to select in the neighborhood: .
The intrinsic parameters of the other peer algorithms are set according to the recommendations in their original papers.
5 Empirical Studies
In this section, we present and discuss the comparison results of our proposed MOEA/AD with the other state-of-the-art peer algorithms. The mean HV values are given in tab:DTLZ_2 to tab:DTLZ_v, where the best results are highlighted in boldface with a gray background.
5.1 Comparisons on DTLZ and WFG Problem Instances
According to Wilcoxon’s rank sum test, and indicates whether the corresponding algorithm is significantly worse or better than MOEA/AD respectively.
Generally speaking, MOEA/AD is the most competitive algorithm for the DTLZ problem instances. As shown in tab:DTLZ_2, it wins in 161 out of 180 comparisons. More specifically, for DTLZ1 and DTLZ3, MOEA/AD obtains the largest HV values on all comparisons except for the 15-objective DTLZ1 and the 3-objective DTLZ3 instances. As for DTLZ2, MOEA/D-PBI obtains the best HV value on the 3-objective case, while MOEA/AD takes the leading position when the number of objectives becomes large. For DTLZ4, the best algorithm varies with different number of objectives. Even though MOEA/AD loses in 4 out of 45 comparisons, the differences are very slight. In addition, as for two decomposition-based algorithms, MOEA/D-PBI obtains significantly worse HV values than MOEA/AD on 14 out 20 comparisons, while Global WASF-GA were significantly outperformed on all 20 comparisons. In particular, Global WASF-GA fails to approximate the entire PF on all DTLZ instances due to its coarse diversity maintenance scheme. As for the two recently proposed Pareto-based many-objective optimizers, the HV values obtained by PICEA-g are significantly worse than MOEA/AD on all problem instances. This can be explained by its randomly sampled target vectors which slow down the convergence speed. VaEA performs slightly better than PICEA-g but it is still outperformed by MOEA/AD on 18 out of 20 comparisons. As expected, the performance of two indicator-based algorithms are not satisfied. In particular, KnEA merely obtains the best HV values on the 15-objective DTLZ2 and DTLZ4 instances. -DEA and NSGA-III, which combine the decomposition- and Pareto-based selection methods together, achieve significantly better results than MOEA/AD in 2 and 3 comparisons respectively, where MOEA/AD beats them in 9 and 11 comparisons respectively. TwoArch2, which also maintains two co-evolving populations, is significantly outperformed by MOEA/AD on all DTLZ instances except for the 15-objective DTLZ1 and the 3-objective DTLZ3. Given these observations, we find that the genuine performance obtained by MOEA/AD does not merely come from the two co-evolving populations. The adversarial search directions and their collaborations help strike the balance between convergence and diversity.