A Constrained Many-Objective Optimization Evolutionary Algorithm With Enhanced Mating and Environmental Selections

Unlike the considerable research on solving many-objective optimization problems (MaOPs) with evolutionary algorithms (EAs), there has been much less research on constrained MaOPs (CMaOPs). Generally, to effectively solve CMaOPs, an algorithm needs to balance feasibility, convergence, and diversity simultaneously. It is essential for handling CMaOPs yet most of the existing research encounters difficulties. This article proposes a novel constrained many-objective optimization EA with enhanced mating and environmental selections, namely, CMME. It can be featured as: 1) two novel ranking strategies are proposed and used in the mating and environmental selections to enrich feasibility, diversity, and convergence; 2) a novel individual density estimation is designed, and the crowding distance is integrated to promote diversity; and 3) the $\theta $ -dominance is used to strengthen the selection pressure on promoting both the convergence and diversity. The synergy of these components can achieve the goal of balancing feasibility, convergence, and diversity for solving CMaOPs. The proposed CMME is extensively evaluated on 13 CMaOPs and 3 real-world applications. Experimental results demonstrate the superiority and competitiveness of CMME over nine related algorithms.

where m is the number of objective functions; x = (x 1 , . . . , x n ) T is an n-dimensional decision vector; n is the number of decision variables; x ∈ S, and S ⊆ R n denotes the search space. g j (x) and h j (x) are the jth inequality and equality constraints, respectively. When m > 3, (1) is a CMaOP [4].
In real-world applications, there are many problems that can be formulated as a CMaOP, such as multicast routing problems [5], software product line problems [6], interval optimization problems [7], and vehicle routing problems [8].
For a CMOP/CMaOP, the degree of constraint violation (CV j ) of x at the jth constraint is calculated as CV j (x) = max(0, g j (x)), j = 1, . . . , p max(0, |h j (x)| − η), j = p + 1, . . . , q where η is a very small numerical quantity to relax the equality constraints into inequality ones. The value 10 −4 is adopted in this work for the chosen test problems; however, other values could be used depending on the application domain. Finally, the overall constraint violation (CV) of x is calculated as The solution x is feasible if CV(x) = 0; otherwise, it is infeasible.
Due to the success of using evolutionary algorithms (EAs) for different optimization problems [9]- [12], the use of EAs for CMaOPs is possible and interesting. However, there are some challenges as follows.
1) Besides the optimization strategies for different objective functions in high-dimensional objective space, efficient constraint-handling techniques (CHTs), such as penalty function [13] and stochastic ranking [14], are required to minimize the CV to handle the constraints.
2) The infeasibility caused by constraints may result in different shapes of constrained Pareto fronts (CPFs), increasing the difficulty faced by constrained manyobjective optimization EAs (CMaOEAs). Therefore, 2168-2267 c 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information.
achieving a good balance among feasibility, convergence, and diversity during the optimization process is essential to handle different types of CPFs in CMaOPs [15].
3) The Pareto dominance pressure dramatically decreases when there is an increase of the objective space dimensionality [16], [17], and this problem becomes harder to solve with the presence of infeasibility caused by constraints. Therefore, the strengthened dominance is crucial to strengthen the selection pressure for CMaOPs. In the evolutionary computation community, research efforts on the use of EAs for solving the CMaOPs is lukewarm [18]. Although there are many existing techniques for CMOPs/CMaOPs, the efficiency of balancing convergence, diversity, and feasibility of these approaches is still insufficient, especially some difficulties or drawbacks in existing research. There is much room to improve the performance of EAs to effectively deal with the CMaOPs. Based on this consideration, in this article, we propose a novel CMaOEA with enhanced mating and environmental selections. The proposed method is referred to as CMME. Its main features include: 1) two new ranking strategies are proposed for feasibility, convergence, and diversity promotion; 2) a novel individual density estimation technique is designed and the crowding distance [19] is used to improve diversity; 3) θ -dominance [20] is applied to enhance the selection pressure on convergence and diversity; and 4) these three techniques are combined into the newly developed enhanced mating and environmental selections to achieve a good balance of feasibility, convergence, and diversity to effectively solve CMaOPs.
The main contributions are as follows. 1) In the enhanced mating selection, a "type 1 ranking" (R #1 ) strategy that prefers feasible or nondominated solutions is proposed for the promotion of feasibility, convergence, and diversity. The feasible solutions are preferred for promotion of feasibility, and the nondominated solutions, without considering any constraint, are preferred for promotion of convergence. Meanwhile, the designed individual density estimation strategy, which prefers solutions located in sparse areas, is applied for the promotion of diversity. 2) In the enhanced environmental selection, a "type 2 ranking" (R #2 ) strategy, which prefers feasible and nondominated solutions is proposed to promote convergence and feasibility. In addition, the crowding distance and θdominance are also used in the environmental selection to improve diversity and convergence. To extensively evaluate the performance of the proposed CMME algorithm, 13 CMaOPs with a variable number of objective functions, 3 three-objective CMOPs, and 3 realworld application CMOPs were chosen as the test suite. CMME was compared with nine related methods. The results indicate the superiority of CMME for solving different CMOPs/CMaOPs.
The remainder of this article is organized as follows. Section II introduces the background and related work. In Section III, the proposed CMME method is proposed in detail, followed by a presentation of the experiments and analysis in Section IV. Finally, Section V presents the conclusions and points out some possible future directions.

II. BACKGROUND AND RELATED WORK
This section first describes the background, including some representative CHTs and the θ -dominance. Then, related CMaOEAs are briefly reviewed.
A. Representatives of Constraint-Handling Techniques 1) Constrained Dominance Principle: Inspired by the feasibility-rule-based CHT [21], Deb et al. [19] developed the constrained dominance principle (CDP) for CMOPs. Specifically, solution x is said to constrained dominate y (denoted as x ≺ CDP y), if any of the following situations is satisfied.
1) x is feasible, while y is not.
2) Both x and y are feasible, and x ≺ y.
Here, x ≺ y indicates x dominates y based on the traditional Pareto dominance. Due to its simple definition and light implementation, CDP is widely used in various CMOEAs, such as NSGA-II [19] and C-MOEA/D [22]. Nevertheless, CDP prefers feasible solutions over infeasible ones [23], [24], which always results in premature convergence, especially when the feasible regions are discontinuous or faulted.
2) ε-Constrained Technique: To overcome the shortage of CDP, a relaxation factor ε is considered and ε-constrained technique has been designed [25]. For two solutions x and y, the ε-constrained dominance x ≺ ε y holds if one of the following conditions is satisfied.
Due to the existence of the relaxation factor ε, infeasible solutions are likely to be selected; thus, the exploration ability is enhanced. In recent years, the ε-constrained technique has been widely used [26]- [29]. However, although it performs well on CMOPs with large infeasible regions, the ε-constrained technique also leads to the instability of search preference. A too large ε favors the selection of infeasible solutions far from feasible regions, which degenerates the search process.
3) Penalty Function: Another possible way to alleviate the shortage of CDP is by using the penalty function [13]. Generally, the fitness F (x) of a solution x based on the penalty function is calculated as where β is the penalty factor. Ever since the proposal of the penalty function, it has received a lot of attention [30], [31]. However, the adaption of parameter β to fit different kinds of problems is still challenging [23]. 4) Two-Rankings-Based Fitness: Ma et al. [23] developed a two-rankings-based fitness (TRF). Two rankings R c (x) and R p (x) for a solution x are defined, where R c (x) and R p (x) are the rankings of x based on CDP and the traditional Pareto dominance, respectively. Note that, in R p (x), the constraints are not considered. Based on these two rankings, the weighted sum fitness of x is calculated as and where P f is the proportion of feasible solutions in the current population, which adaptively reflects the evolution status. It is used to achieve a balance between constraints and objective function values. But actually, if the population mainly consists of feasible solutions (α ≈ 1), then TRF has the same effect as CDP. In this situation, infeasible solutions are very likely to be rejected; thus, jumping out of the local feasible region becomes hard. Therefore, when dealing with some convergence-hard or diversity-hard CMOPs [4], TRF will fall into the same dilemma (trapped in local optimal) as CDP does. 5) Two-Archive Technique: Li et al. [10] introduced a constrained two-archive EA (C-TAEA) for CMOPs and CMaOPs, where the convergence, feasibility, and diversity are balanced through the two-archive technique (TAT). In TAT, two archives [i.e., convergence archive (CA) and diversity archive (DA)] co-evolve simultaneously. The archive CA is used to promote the convergence and feasibility, while DA is used to enhance diversity. Experimental results showed that the C-TAEA can handle CMOPs and CMaOPs with intersected and discontinuous feasible regions [10]. However, as pointed out in [24], its mating selection strategy usually chooses one parent individual from CA and the other from DA. The generated offspring will locate in the middle region between CA and DA; thus, neither possessing better convergence, feasibility than individuals in CA, nor being able to promote the diversity of DA. In other words, the generated offspring contain much useless information.
6) Two-Stage Technique: Some researchers have applied the two-stage strategy for constraint handling. ToP [32] deals with CMOPs by separating the search process into two phases. In the first phase, a reference vector is used to transform the original CMOP into a constrained single-objective optimization problem. Then in the second phase, the original CMOP is optimized. CMOEA-MS [33] adopts two adaptively switched stages. One solves the original CMOP, and the other ignores all constraints to solve the unconstrained MOP. Pull and push search (PPS) [34] includes push and pull stages. The push stage ignores all constraints to push the population near the PF, then the pull stage considers all constraints to pull the population near the CPF.
where i = 1, . . . , m; z * consists of min x∈S f i (x), and z nad consists of max x∈S f i (x). Let w be the reference vector of x; d 1 (x) be the projection distance, and d 2 (x) be the vertical distance. d 1 (x) and d 2 (x) can be formulated as Then, PBI is calculated as where θ > 0 is a predefined penalty factor, which is usually set to be 5.
2) M2M: To improve the performance of MOEA/D on problems with an extreme bias, Liu et al. [36] developed the M2M technique, which decomposes the objective space into a number of simple subregions. In M2M, each subregion has its own subpopulation. Based on M2M, the density of each subregion can be evaluated to reflect the crowding degree of each subpopulation. Fig. 1 shows the subregion division with M2M, where the objective space is divided into different subregions by the reference vectors.
Definition 1: Solution x is said to θ -dominate solution y, denoted as x ≺ θ y, if both x and y are in the same subregion, and PBI(x) < PBI(y).
In Fig. 2, the θ -dominance principle is illustrated. There are five reference vectors (w 1 , . . . , w 5 ) and seven individuals (x 1 , . . . , x 7 ). Based on M2M, x 1 , x 2 , and x 5 are associated with w 1 , w 2 , and w 4 , respectively; x 3 and x 4 are associated with w 3 ; and x 6 and x 7 are associated with w 5 . Then, based on the θ -dominance, x 1 , x 2 , x 3 , x 5 , and x 6 are assigned to the first front, while x 4 and x 7 are in the second front. Therefore, combining the advantages of M2M and PBI, the θ -dominance can provide the strengthened selection pressure on both convergence and diversity [20].

C. Existing CMaOEAs
This section briefly reviews some representative related CMaOEAs.
Combined with NSGA-II [19], Saxena et al. [18] presented an EA with ε-dominance, controlled infeasibility, and a nonlinear dimensionality reduction scheme for CMaOPs. Jain and Deb [22] embedded the CDP [19] into NSGA-III [37] and MOEA/D [35] to extend them for CMaOPs. In I-DBEA [38], an adaptive ε formulation was designed to handle the constraints. Cooperating with a preemptive distance comparison scheme and two independent distance measures, I-DBEA is capable of dealing with CMaOPs. Elarbi et al. [39], [40] introduced an isolated solution-based constrained Pareto dominance, which favors infeasible solutions associated with isolated subregions or with a smaller CV to exploit infeasible isolated solutions. Miyakawa et al. [41] proposed the CMOEA/D-DMA for CMaOPs, where the directed mating and archives of infeasible solutions are used. Li et al. [10] introduced the C-TAEA with TAT to balance convergence, feasibility, and diversity to deal with CMaOPs. In CCMODE [42], a co-evolutionary framework with (m + 1) subpopulations uses m subpopulations for constrained single-objective optimization and an archive subpopulation for constrained m-objective optimization. Zhou et al. [43] designed a tri-goal evolution framework, where three indicators for convergence, diversity, and feasibility are presented for CMaOPs. Fan et al. [29] used a PPS framework with a two-stage strategy, which handles the constrained and unconstrained MaOPs separately to accelerate the approach to the CPF. Wang and Xu [44] developed an angle-based EA with infeasibility information, in which different kinds of infeasible solutions are utilized in the environmental and mating selections. Jiao et al. [26] transformed CMaOP into a dynamic CMaOP with the -CDP for handling constraints and optimizing objectives simultaneously.
The features and drawbacks of these existing CMaOEAs are shown in Table I. III. OUR APPROACH: CMME In this section, the proposed method is elucidated. Sections III-A-III-E present the details of CMME. Then, Section III-F elaborates the significance and differences between our method and existing methods. Finally, Section III-G analyzes the computational complexity.  [45], much less research has been done on CMaOPs. There is still much room to improve the performance of EAs for dealing with CMaOPs.
Balancing feasibility, convergence, and diversity is the main challenge for handling CMaOPs. To effectively solve CMaOPs, one possible way is to design new enhanced mating and environmental selections [46] which can achieve the balance of feasibility, convergence, and diversity. Based on this consideration, we propose a novel CMaOEA with enhanced mating and environmental selections to achieve a good balance of feasibility, convergence, and diversity.
2) Expanding θ -Dominance to CMaOEA: Experimental results show that θ -dominance is effective on handling MaOPs. With the help of reference vectors and the strengthened dominance relation, θ -DEA, which applies θ -dominance, performs well on dealing with MaOPs. By optimizing PBI in each subregion, θ -dominance can make a good tradeoff between convergence and diversity. However, it cannot solve CMaOPs without a CHT to handle constraints. Therefore, expanding the effective θ -dominance to CMaOEA is worth trying.
3) Multiranking for Balancing Convergence, Diversity, and Feasibility: As mentioned in Section II-A and Table I, the existing CHTs, as well as the existing CMaOEAs, still have some shortages that must be overcome. Inspired by the idea of TRF and TAT, we use a multiranking strategy to achieve a balance of different goals. Further, Pareto dominance has a preference for well-converged solutions, while the CV value represents the feasibility performance. Since the core issue of CMaOEA is achieving a balance between convergence, diversity, and feasibility, designing different rankings for different goals is promising. Also, the application of different rankings in mating and environmental selection is promising.

B. Individual Density Estimation
Based on M2M, the objective space can be divided into different subregions. To promote the diversity of the solution set, the individual density estimation is proposed as shown in Algorithm 1. First, the objective function values of each individual in the solution set S are normalized by using (6) (line 1). Also, the N subregion density index is initialized, which is the same as the size of the reference vector set (line 2). Then, each individual is associated with its nearest reference vector based on the cosine similarity (line 4). The individual density is estimated by calculating the number of  6: c j = c j + 1; 7: end for 8: for i = 1 to |S| do 9: individuals in this individual's located subregion. Note that the constraints are not considered when estimating the individual density.
To better understand the principle of Algorithm 1, Fig. 3 gives an example for the individual density estimation according to M2M. In Fig. 3, the objective space is divided into five subregions based on the five reference vectors (i.e., w 1 , . . . , w 5 ). Based on the cosine similarity between the individual and its nearest reference vector, each individual is associated with the corresponding reference vector. Then, the density of the individual is represented by the number of individuals in its located subregion. In this way, d(

C. Enhanced Mating Selection
To effectively solve CMaOPs, besides convergence and diversity, feasibility is very important. In the mating selection, one of the main goals is to detect new feasible areas. Therefore, it is essential to select suitable parents to form the mating pool M for offspring generation. Based on this consideration, we propose the R #1 strategy to assign the rank of an individual to measure its suitability.
1) Type 1 Ranking (R #1 ): The R #1 rank of an individual x is defined as follows.
Definition 2: R 1 (x) = 0, if x is feasible or it is a nondominated individual; otherwise, R 1 (x) is equal to the number of individuals that dominate x according to CDP [19].
According to Definition 2, all feasible individuals and all nondominated individuals are assigned the best ranks. They have more chances to be selected into M. The pseudocode of ranking calculation is shown in Algorithm 2. First, the Algorithm 2 RankingCalculation(S, t r ) Input: S (solution set), t r (type of ranking) Output: R (ranking set) 1 9: s = 0; 10: for j ← 1 to |S| do 11: if x j ≺ CDP x i then 12: s = s + 1; 13: end if 14: end for 15: R(x i ) ← s; 16: end if 17: constraint violation and dominance front number of S are calculated and stored in CV and F (lines 1 and 2), respectively. Then, the rank of each individual in S is calculated according to a specific ranking strategy. t r is a flag to determine which type of ranking strategy is used. As depicted in Fig. 4, although x 1 and x 2 are infeasible solutions, they are ranked 0 because they are in the first dominance front. Their inclusion in M could guide the search toward undeveloped feasible regions. x 3 is not in the first front; however, it is also ranked 0 because of its feasibility. The inclusion of x 3 will promote the exploitation of developed feasible regions.
2) Mating Selection With R #1 : Combined with R #1 , the enhanced mating selection is given in Algorithm 3. It is based on the tournament selection. First, two randomly selected individuals x and y are compared based on their ranks and the better one is added into M. If both of them have the same ranks, the one with smaller density is added into M. If R(x) = R(y) and d(x) = d(y), then either x or y is randomly chosen for M. The above procedure continues until |M| ≥ N. Randomly pick x or y to M; 15: end if 16: end if 17: end while 18 The benefit from using the R #1 strategy and individual density estimation is that the individuals in M could guide the search to explore new feasible regions and exploit the existing feasible regions with diversity promotion. Therefore, the mating selection is able to make a good tradeoff between exploration and exploitation on feasibility. Note that if the population only contains feasible solutions, which usually appears in the late stage of evolution, the mating selection will promote convergence and diversity.

D. Enhanced Environmental Selection 1) Type 2 Ranking (R #2 ):
In environmental selection, a strengthened selection pressure is required to balance feasibility, convergence, and diversity to obtain the next parent population. Therefore, we propose the R #2 strategy to assign the rank of an individual x.
Definition 3: R 2 (x) = 0, if x is a feasible and nondominated individual; otherwise, R 2 (x) is equal to the number of individuals that dominate x according to CDP.
From Definition 3, the individuals are ranked 0, if they are feasible and in the first front. As depicted in Fig. 5, only the red points are assigned the best ranks. Therefore, the R #2 strategy while |S| < N do 14: 16: end while 17: while |S| > N do 18: D ← IndividualDensityEstimation(S, W) as shown in Algorithm 1; ⇐
2) Environmental Selection With R #2 : Based on the R #2 strategy, the enhanced environmental selection is developed in Algorithm 4, where "⇐" indicates the differences between the proposed environmental selection and the CA update strategy in C-TAEA [10]. First, the crowding distance of each individual in C is calculated by the technique presented in NSGA-II [19]. All crowding distance values are saved into D cr (line 2). Then, the distance of each individual is combined with its rank as the supplementary diversity (line 3). Subsequently, all individuals in C with ranks less than 1 are selected into the candidate set S c (lines [4][5][6][7][8]. Based on the size of S c , there are three situations as follows. 1) |S c | = N: In this case, S c can be used directly as the population for the next generation. 2) |S c | > N: This situation means that most of the individuals in C are feasible with good distribution. In this case, the selection strategy should be designed to promote convergence and diversity. 3) |S c | < N: On the contrary, this situation reflects that either the number of feasible solutions with good distribution is not enough or the infeasible solutions are of the majority. Thus, the selection strategy should be inclined to the feasibility. For |S c | > N, some poor performed feasible individuals are discarded from S c to obtain the next population P with balanced convergence and diversity.
1) First, S c undergoes the fast nondominated sorting with the θ -dominance to obtain the front F (line 12).  (10) where z * is the ideal point; w k is the reference vector corresponding to the most crowded subregion, and f j (x) represents the jth objective value of x. For |S c | < N, the CA update strategy proposed in C-TAEA [10] is used to explore the feasible regions.
1) First, the infeasible individuals in C are separated and saved in SI in line 31. 2) Then, the individuals in SI are evaluated by a new biobjective optimization problem formulated as The front FI of SI is obtained based on the fast nondominated sorting with the MOP shown in (11) (line 32). M ←MatingSelection(P, N, R 1 , D); 7: O ←OffspringGeneration(M); 1) The R #2 strategy is able to enhance both feasibility and convergence. 2) Both the crowding distance and individual density estimation can promote diversity.

Algorithm 5 Framework of CMME
3) The use of the θ -dominance is capable of enhancing the selection pressure on both the convergence and diversity. 4) The CA update strategy at the situation |S c | < N in [10] can pursue the feasibility.

E. Framework of CMME
Combining the proposed procedures, the framework of the proposed CMME method is shown in Algorithm 5. Initially, the reference vector set W with size N is randomly generated as suggested in [47], followed by the population P's initialization. At each generation, the following procedures are repeated until the termination condition is met. 1) In line 4, the density of each individual in P is estimated as shown in Algorithm 1. 2) In lines 5 and 6, first the individuals in P are ranked based on the type 1 ranking (R #1 ) in Algorithm 2. Then, based on the ranking R 1 and density D, the mating pool M is obtained based on Algorithm 3. 3) In line 7, the offspring population O is generated through the genetic operators. Then, in line 8, P and O are combined to form the combined set C. 4) In line 9, C is ranked based on the type 2 ranking (R #2 ). And in line 10, the environmental selection procedure is performed to generate the new population P in Algorithm 4.

F. Remarks
1) CMME includes two novel ranking mechanisms. The designed environmental and mating selection strategies are based on these two ranking strategies. Also, the individual density, crowding distance, and θ -dominance play essential roles. All these components work together to achieve the balance of convergence, diversity, and feasibility to effectively handle CMaOPs. 2) The proposed two ranking strategies are different from the TRF technique [23]. In TRF, the two rankings are used in the fitness calculation as formulated in (5). However, in CMME, the R #1 strategy is used in the enhanced mating selection which aims to explore the undeveloped feasible regions and to exploit the existing feasible regions. The R #2 strategy is tailored for the enhanced environmental selection which aims to promote both convergence and feasibility.
3) The R #1 strategy prefers to select nondominated or feasible solutions as mating parents, which can guide the generation of offspring to nondominated regions and feasible regions, thus improving the convergence and feasibility of population. Besides, choosing nondominated solutions as mating parents can help with jumping out of local feasible regions as shown in Fig. 4. Therefore, the performance on convergence [ Fig. 4

4) The proposed environmental selection in Algorithm 4
is similar to the CA update strategy in C-TAEA [10]. However, there are some important differences marked as ⇐ in Algorithm 4: a) the CA update strategy only considers the feasibility, while in our approach, the feasibility, convergence, and diversity are all taken into consideration. Based on R #2 ranking, the promotion of convergence and feasibility is ensured. Meanwhile, the crowding distance helps with diversity promotion; b) the traditional Pareto dominance for the CA update strategy is replaced by the θ -dominance in the nondominated sorting to strengthen the selection pressure for the CMaOPs, thus improving convergence and diversity; and c) individual density estimation is used in the enhanced environmental and mating selections, which has an impact on searching the sparse subregions, thus improving diversity. Therefore, the synergy of these techniques under the framework of CMME can achieve the balance among convergence, diversity, and feasibility to handle CMaOPs. 5) In the enhanced mating selection, the traditional Pareto dominance is used; while in the enhanced environmental selection, the θ -dominance is used. The reason is that in the mating selection, feasibility is more important; however, the θ -dominance may discard some feasible solutions due to its strengthened selection pressure as mentioned in Section II-B3. On the other hand, in the environmental selection, the convergence and diversity at the situation |S c | > N are more important, and thus, the θ -dominance is applied.

G. Computational Complexity of CMME
The computational complexity of CMME is composed of the following.
1) The complexity of generating reference vector set W and initial population P are O (mN) and O(mnN), respectively.
2) The complexity of individual density  O(mN 2 ). Therefore, the overall complexity of CMME is O(mN 2 ).

IV. EXPERIMENTAL RESULTS AND ANALYSIS
In this section, the performance of CMME is extensively evaluated by comparing with nine related methods on the C-DTLZ [22], DC-DTLZ [10], DAS-CMOP [4], and MW [15] test functions and three real-world CMOPs [48].
3) Genetic Operators and Parameter Settings: For CMME, the simulated binary crossover (SBX) [50] and polynomial mutation (PM) [19] were used as the genetic operators with the following parameter settings.
1) Crossover probability p c = 1 and distribution index η c = 20. 2) Mutation probability p m = 1/n and distribution index η m = 20. The population size N and maximal number of function evaluations E max for the ten methods on C-DTLZ, DC-DTLZ, and MW benchmark problems with more than three objectives are listed in Table S-I in the supplementary file. For the threeobjective DAS-CMOP and MW instances, E max was set to 500. It is guaranteed that all algorithms have the possibility to converge to CPF, and E max of all algorithms is consistent for a fair comparison. The parameter settings for the nine methods  II  AVERAGE RANKINGS OF HV AND IGD BY THE FRIEDMAN  TEST FOR DIFFERENT CMAOEAS used for comparison were the same as suggested in their corresponding original references. All the parameter settings in the following text remain the same unless a change is mentioned. 4) Performance Indicators: Two widely used indicators, inverted generational distance (IGD) [51] and hypervolume (HV) [52], are chosen to compare the performance of different algorithms. The detailed descriptions of them are given in the supplementary file.
Following the approach in [24], 10 000 uniformly distributed points were sampled on the true PF for the IGD calculation. For the HV calculation, the objective values were first normalized, and in the normalized objective space, (1.1, 1.1, . . . , 1.1) was adopted as the reference point.

B. Comparison With Other Algorithms
In this section, CMME is compared with the aforementioned nine algorithms. Each algorithm was executed over 30 independent runs on each test case. The mean and standard deviations of HV and IGD indicators were recorded. The KEEL software was used to calculate the statistical results [53]. The Wilcoxon test with a significance level of 0.05 and the Friedman test with the Bonferroni correction at a significance level 0.05 were adopted to perform statistical analysis on the experimental results. In addition, according to the Wilcoxon test, we use "+," "−," and "=" to show that the compared algorithm is significantly better than, significantly worse than, or statistically similar to CMME, respectively.
The detailed HV and IGD results are, respectively, reported in Tables S-II-S-V in the supplementary file, where "NaN" indicates that the algorithm cannot find feasible solutions on this problem, and "null" means that the obtained feasible solution set is far away from the true PF. In addition, the multiple-problem analysis is also performed to further verify the performance of CMME. 1 The average rankings of all the algorithms by the Friedman test are shown in Table II. The summarized results by the Wilcoxon test on the HV and IGD metrics are given in Tables S-VII and S-VIII in the  supplementary file, respectively. 1) Description About the Results: In terms of the HV indicator on C-DTLZ and DC-DTLZ, the results in Table S-II  in the supplementary file, Table II, and Table S-VII in the supplementary file indicate that: 1) CMME obtained the best HV values on 14 cases, followed by PPS (11), C-TAEA (9), CCMO (6), NSGA-III (3), and DCNSGA-III (2). For NSGA-II-ToR, ToP, I-DBEA, and TiGE-2, there was no best HV value on any case; 2) as shown in the last column of Table S-II in  the supplementary file, CMME significantly outperformed NSGA-II-ToR, ToP, CCMO, NSGA-III,  I-DEBA, PPS, C-TAEA, TiGE-2, and DCNSGA-III on 38,37,26,22,38,29,19,36, and 24 cases, respectively. However, CMME only lost on 2, 4, 8, 4, 3, 10, 8, 3, and 4 cases, respectively; 3) CMME has the best ranking among the ten algorithms as shown in Table II, followed by C-TAEA; 4) among all the compared CMaOEAs, CCMO performed well on most three-objective problems since it is tailored for CMOPs. PPS performed well on C1-DTLZ3 and DC-DTLZ3 instances, of which the PF and CPF are consistent and the ignorance of constraints at the push stage accelerates the convergence speed to the PF as well as the CPF in dealing with multimodal instances; 5) according to the multiple-problem analysis by the Wilcoxon test shown in Table S-VII in the supplementary file, it is clear that except for C-TAEA, CMME yielded significantly better results than the other algorithms on the whole. With respect to IGD, from the results in Tables S-II, S-III, and S-VIII in the supplementary file, we can observe that: 1) on 17 cases, CMME obtained the best IGD values, followed by CCMO (10), C-TAEA (9), NSGA-III (6), PPS (5), and TiGE-2 (4). We find that CCMO performed well on the three-objective functions; however, its performance deteriorated dramatically for the higher objective problems. CMME was more robust for the test CMaOPs at a variable number of objective functions; 2) in addition, CMME obtained significantly better results on the majority of the test cases than NSGA-II-ToR (47) (27); 3) CMME had the best IGD ranking among the ten algorithms; 4) CMME significantly outperformed the other nine algorithms based on the Wilcoxon test shown in Table S-VIII in the supplementary file. For DC3-DTLZ, the current E max does not seem to be enough for most CMaOEAs to converge to the CPF; therefore, the results of HV and IGD contain many NaN and null values. Next, the final solution sets on DC3-DTLZ are plotted to further study the efficiency of CMME.
To obtain a more intuitive observation, we select three-objective C1-DTLZ1, DC1-DTLZ3, DC2-DTLZ3, and DC3-DTLZ3; five-objective C2-DTLZ2 and DC2-DTLZ3; eight-objective DC1-DTLZ1; 10-objective DC1-DTLZ3 and DC2-DTLZ1, and 15-objective DC3-DTLZ3 as representatives. The obtained feasible PFs of the ten algorithms with the median HV value among 30 runs are plotted in Figs. S-1-S-10 in the supplementary file. From these figures, we can see that: 1) on three-objective C1-DTLZ1, CMME, CCMO, and DCNSGA-III obtained the best PFs with the most promising convergence and diversity. The obtained PFs of other seven algorithms are poor on either the convergence or diversity. On three-objective DC1-DTLZ3, only CMME finally converged to the CPF. On three-objective DC2-DTLZ3 and DC3-DTLZ3, CMME was stuck in the outer layer; 2) on five-objective C2-DTLZ2, NSGA-III yielded the best PF, CMME was slightly worse than NSGA-III on diversity; 3) on eight-objective DC1-DTLZ1, C-TAEA achieved the best results while CMME and DCNSGA-III performed slightly worse on diversity. On 10-objective DC1-DTLZ3, although all the algorithms failed to converge to the CPF, CMME obtained the best convergence in dealing with multimodal situation and discontinuous feasible regions; 4) on 10-objective DC2-DTLZ1, CMME performed the best on both convergence and diversity, while on fiveobjective DC3-DTLZ3, although all the CMaOEAs did not converge to the CPF, CMME and C-TAEA obtained the best convergence and diversity; 5) on 15-objective DC3-DTLZ3, CMME also obtained the best convergence. In terms of DAS-CMOP and MW shown in Tables S-IV (HV) and S-V (IGD) in the supplementary file, CMME outperformed other algorithms on almost all of the three instances on DAS-CMOP instances, and it can also obtain the best overall results on MW instances. Fig. S-11 in the supplementary file showed the feasible and nondominated solutions with the median HV value among 30 runs on 3-objective DAS-CMOP and MW instances obtained by CMME. It can be seen that CMME finally approached the CPF for the selected cases, except for DAS-CMOP9.
In addition, DAS-CMOP8 and three-objective MW8 were chosen to compare the performance of different methods. The final solution sets obtained by CMME and other CMaOEAs on DAS-CMOP8 and MW8 were presented in Figs. S-12 and S-13 in the supplementary file, respectively. As shown in these figures, CMME obtained better results than most of the other algorithms, except for CCMO.
2) Analysis and Summaries of the Results: According to the above results, we can summarize that: 1) CMME performed well on C1-DTLZ1, C2-DTLZ2, C3-DTLZ4, DC1-DTLZ, DAS-CMOP7, and DAS-CMOP8; 2) CMME was competitive with other algorithms on DC2-DTLZ1 and DC3-DTLZ; 3) CMME performed poorly on C1-DTLZ3, DC2-DTLZ3, and DAS-CMOP9; 4) CMME obtained the best overall results on the MW problems. On most of the MW instances, CMME outperformed other methods; 5) CMME performed well on MW4 and MW8 but performed poorly on MW14. For C1-DTLZ1, C2-DTLZ2, and C3-DTLZ4 with relatively simple constraints, the CDP is enough to handle their constraints. In addition, θ -dominance applied in CMME is efficient for dealing with regularly shaped CPFs. Therefore, CMME is suitable for these instances. For DC1-DTLZ1 and DC1-DTLZ3, the R #1 ranking increases the algorithm's ability to explore discontinuous feasible regions as shown in Fig. 4. Therefore, even when the CPFs were irregular, CMME obtained better results. For MW4, the PF is the same as the CPF, and there are no outer layer infeasible regions, so CDP is sufficient for handling constraints. The θ -dominance is also suitable for its regular CPF shape, so CMME obtained a well-distributed solution set. As for DAS-CMOP7 with convergence hardness, R #2 ranking and the enhanced environmental selection strengthened the selection pressure so that CMME obtained a better result. As for DAS-CMOP8 with feasibility hardness, the CDP's prioritization of feasibility ensures the feasibility promotion. Together with the convergence ability provided by R #2 ranking, CMME also obtained a good result.
DC2-DTLZ and DC3-DTLZ have inner outer layer feasible regions, and with the help of R #1 ranking and the designed mating strategy, CMME could across the infeasible regions and reach the inner layers. However, within the given times of function evaluations, none of these algorithms could deal with the multimodal DC2-DTLZ3 and DC3-DTLZ3 situations as shown in Figs. S-3 and S-4 in the supplementary file. As shown in Figs. S-9 and S-10 in the supplementary file, CMME reached closer to the CPF than other algorithms when m increased. Therefore, the trend is obvious that CMME suites many-objective situations. For MW8 and MW14, CMME only obtained competitive results due to the irregular CPF feature.
For C1-DTLZ3 and DC2-DTLZ3, CMME was worse than PPS. For these multimodal instances, the designed push stage of PPS can efficiently push the population near the PF by ignoring constraints. Therefore, although PPS could not obtain good diversity, it had better convergence. As for DAS-CMOP9, all these approaches could not obtain good results since DAS-CMOP9 has both feasibility and diversity hardness. The poor performance on multimodal and diversity hardness instances indicates the difficulty of CMME on diversity promotion.
According to the results of C-TAEA and CMME depicted in these figures, especially on three-objective instances, we can see that θ -dominance applied in CMME significantly promotes diversity performance of the final results. Although the proposed enhanced environmental selection is based on the CA update framework of C-TAEA, the performance of CMME was much better than C-TAEA due to the help of two rankings (i.e., R #1 and R #2 ) and the mating environmental selections.
Therefore, from our results and analysis, we conclude that benefited from the enhanced mating and environmental selections, the CMME algorithm achieved a good balance of feasibility, convergence, and diversity on the CMaOPs with variable numbers of objective functions. It provided better or highly competitive results compared with nine state-of-the-art algorithms. However, CMME lacks diversity promotion ability, which results in poor performance on multimodal and diversity hardness instances.

C. Influence of Different Components in CMME
This section empirically discusses the influence of different components in CMME. There are four variants to be compared with the proposed CMME method are as follows.  III  AVERAGE RANKINGS OF HV AND IGD BY THE FRIEDMAN TEST FOR  DIFFERENT CMME VARIANTS   TABLE IV  RESULTS OBTAINED BY THE WILCOXON TEST, WHERE CMME IS  COMPARED WITH ITS VARIANTS   1) CMME-1: The traditional Pareto dominance replaces the θ -dominance used in the environmental selection as shown in line 12 of Algorithm 4. 2) CMME-2: The TRF substitutes for the R #1 in the mating selection to evaluate a solution. 3) CMME-3: In the R #2 strategy, both the feasibility and nondominance are considered. In C-TAEA, only the feasibility is used. Therefore, in CMME-3, the nondomination is removed to verify the effectiveness of the R #2 strategy. 4) CMME-4: In the mating and environmental selections, the individual density and crowding distance are used for diversity promotion. In CMME-4, both of them are removed to show their influence on the performance of CMME. Table S-VI in the supplementary file reports the detailed HV and IGD results for the five CMME variants. The multipleproblem statistical results by the Friedman and Wilcoxon tests are given in Tables III and IV, respectively. From the results, it can be seen that which are described as follows.
1) CMME Versus CMME-1: CMME obtained significantly better results on 17 and 24 cases in terms of HV and IGD, respectively. CMME lost on nine and nine cases on HV and IGD. In the rest of the cases, both of them obtain similar results. The results in Table IV show that CMME obtained slightly better HV performance, however, it significantly outperformed CMME-1 on the IGD metric. This indicates that the strengthened selection pressure by the θ -dominance in the environmental selection can improve the performance of CMME on CMaOPs. 2) CMME Versus CMME-2 and CMME Versus CMME-3: CMME yielded significantly better results in the majority of the test cases compared with CMME-2 and CMME-3. It also consistently outperformed both CMME-2 and CMME-3 on HV and IGD according to the Wilcoxon test. The results mean that the proposed two ranking strategies play essential roles in CMME to effectively solve CMaOPs.
3) CMME Versus CMME-4: In the majority of the cases, both obtained similar results. CMME had slightly better HV results, but performed significantly better than CMME-4 on IGD. This phenomenon proves the effectiveness of the individual density used in the mating selection and the crowding distance used in the environmental selection. 4) The results in Table III show that CMME consistently obtained the best average rankings on both HV and IGD. This means that the synergy of different components can improve the performance of CMME for CMaOPs.

D. Experiments on Real-World Problems
In this section of experiments, the efficacy of CMME was verified through experiments on real-world applications. Three instances were chosen from [48] as CMaOPs from the real-world applications: 1) gear box design problem [54]; 2) multiproduct batch plant problem [55]; and 3) heat exchanger network design problem [56].
The parameter settings are shown in Table S-IX in the supplementary file. The results of HV are presented in Table S-X in the supplementary file, from which we can see that generally CMME outperformed the other algorithms on these three instances except for NSGA-III.
For three-objective Gear Box Design, CMME is only worse than NSGA-III. For three-objective Multiproduct Batch Plant and Heat Exchanger Network Design, CMME outperformed all the other algorithms.
To further analyze the results, we depict the final solution sets of CMME and the algorithms in comparison on multiproduct batch plant (denoted as MpBP in the figures) in Fig.  S-14 in the supplementary file. From Fig. S-14 in the supplementary file, we can see that CMME performed well on MpBP. The convergence and diversity performance of CMME were significantly better than the results of other methods.
In summary, CMME performed well on these threeobjective real-world application problems.

V. CONCLUSION AND FUTURE WORK
In this article, we have proposed a CMaOEA with enhanced mating and environmental selections, referred to as CMME, for CMaOPs. Specifically, two ranking strategies were designed for evaluating the convergence and feasibility of individuals. In addition, an individual density estimation was developed to improve diversity. Also, the θ -dominance has been embedded into the environmental selection to strengthen the selection pressure for better handling of the many-objective situations. We found that the synergy of different components in CMME achieves a good balance of feasibility, convergence, and diversity. The proposed CMME algorithm was extensively evaluated through 13 benchmark CMaOPs with a variable number of objective functions, three benchmark CMOPs with three objective functions and three real-world applications. CMME was compared with nine related algorithms. The results show that CMME is able to provide better or highly competitive performance on the whole. However, there remain some limitations to be studied in the future. For example, the performance of CMME on the multimodal test and high-dimensional instances, such as C-DTLZ3 and DC-DTLZ3, needs to be further improved. Moreover, the diversity performance of CMME also needs to be further improved. On CMaOPs with irregular CPF, such as MW14, the performance of CMME is poor due to the evenly distributed reference vector set. In addition, combining the proposed two ranking strategies with other enhanced dominance relations may develop other improved CMaOEAs for handling CMaOPs effectively.
The MATLAB source code of CMME can be obtained from the authors upon request.