A Dual-Population-Based NSGA-III for Constrained Many-Objective Optimization

The main challenge for constrained many-objective optimization problems (CMaOPs) is how to achieve a balance between feasible and infeasible solutions. Most of the existing constrained many-objective evolutionary algorithms (CMaOEAs) are feasibility-driven, neglecting the maintenance of population convergence and diversity when dealing with conflicting objectives and constraints. This might lead to the population being stuck at some locally optimal or locally feasible regions. To alleviate the above challenges, we proposed a dual-population-based NSGA-III, named DP-NSGA-III, where the two populations exchange information through the offspring. The main population based on the NSGA-III solves CMaOPs and the auxiliary populations with different environment selection ignore the constraints. In addition, we designed an ε-constraint handling method in combination with NSGA-III, aiming to exploit the excellent infeasible solutions in the main population. The proposed DP-NSGA-III is compared with four state-of-the-art CMaOEAs on a series of benchmark problems. The experimental results show that the proposed evolutionary algorithm is highly competitive in solving CMaOPs.


Introduction
Constrained multi-objective optimization problems (CMOPs) widely exist in the real world, such as the design of the water distribution network [1] and the vehicle routing problem with time windows [2]. In this paper, a dual-population-based NSGA-III with an auxiliary population that ignores constraints is proposed. The offspring of its auxiliary population can guide the evolution of the main population.
In the last decade, several new methods have been proposed to solve unconstrained multi/many-objective optimization problems, such as a new competitive mechanism [3], adaptive operator selection [4], and adaptive weight vector generation strategy [5]. In recent years, CMOPs have attracted increasing research interest in the field of evolutionary computation. Constrained multi-objective evolutionary algorithms (CMOEAs), such as PPS [6], ToP [7], and c-DPEA [8], have achieved good performance on CMOPs. However, CMaOPs involve both complex constraints and highdimensional objective spaces, so it is inefficient to solve CMaOPs directly with existing CMOEAs. Specifically, the difficulties in solving CMaOPs are as follows: first, even in the unconstrained case, balancing convergence and diversity in a high-dimensional objective space is difficult to achieve [9]. Although CMOEAs can handle constraints, they cannot efficiently select excellent solutions in a high-dimensional objective space. Meanwhile, the algorithm should balance the search between feasible and infeasible regions. Particularly, some CMaOPs have complex constraints. It is difficult for individuals to cross multiple disjoint infeasible regions to approach the constraint Pareto front (CPF), which leads to poor convergence and diversity of the population [10]. In addition, most evolutionary algorithms just simply push solutions toward the feasible boundary, and it is difficult to balance the search between feasible and infeasible regions [11]. This may cause the population to fall into locally optimal or locally feasible areas, especially when the feasible region is discrete, narrow, or far from the unconstrained Pareto front (UPF) [12]. To alleviate the above problems, this article proposes a dual-population-based constrained many-objective evolutionary algorithm, named DP-NSGA-III. It includes two collaborative and complementary populations. The main characteristics of DP-NSGA-III are as follows: 1.
To take advantage of the infeasible solution with excellent objectives, we introduce an auxiliary population, and its environmental selection [13] is different from the main population, using p-norm to select individuals. Two populations cooperate by exchanging offspring. Each population can evolve more efficiently using information from the other population. Specifically, Population1 is the main population, which finds the CPF by the ε-constraint handling method. Population2 with less convergence pressure is the auxiliary population, which will eventually converge to UPF due to the neglected constraints.

2.
From the point of view of the problem, neither population solves the original CMaOPs.
Since the ε-constraint boundary replaces the original constraint, Population1 solves a problem with looser constraints than the original problem and is equivalent to the original problem at the end of evolution. Population2 solves the unconstrained problem corresponding to the original problem. This strategy indirectly solves the original CMaOPs. 3.
To make the main population of DP-NSGA-III infeasibility-driven, we design a monotonically decreasing ε-constraint handling function that takes no additional parameters. That is, the ε-bounded boundary is larger in the early period and close to 0 in the late period. Therefore, this function has a greater impact in the early evolutionary period. The function aims to involve partially infeasible solutions in the evolution, balancing convergence and diversity by using excellent infeasible solutions.
The rest of this article is organized as follows: Section 2 briefly introduces the existing CMOEAs or CMaOEAs. The details of the proposed DP-NSGA-III are described in Section 3. Sections 4 and 5 present the experimental setup and experimental results, respectively. Finally, the conclusions and future work are given in Section 6.
The constraint-dominance principle(CDP) [15] was one of the earlier well-known constraint handling methods. NSGA-II [15] employs the CDP to compare pairs of solutions. For two infeasible solutions, the one with less constraint violation is better; the two feasible solutions are still compared by the Pareto dominance. Both C-NSGA-III [20] and C-MOEA/D [20] are well-known representative algorithms that embed this technique. Fan et al. [21] proposed an angle-based CDP, which considers the ratio of feasible solutions and the angular relationship between individuals. Among the above algorithms, NSGA-III is designed for many-objective optimization problems, so C-NSGA-III can solve some CMaOPs.
In recent years, some evolutionary algorithms based on multiple stages have been designed to solve CMOPs. PPS [6] divides the search process into two stages: Push and Pull. The Push phase does not consider constraints, which helps to cross large infeasible areas and push the population to UPF. In the Pull phase, an improved ε-constraint handling method is used to pull the population back to the CPF. Liu et al. [7] proposed a phased framework ToP.
In ToP, the first stage uses a weighted sum to transform the CMOP into an unconstrained single objective optimization problem to find a promising feasible region. The second stage uses the existing CMOEA to find the CPF of the original problem. CMOEA-MS [22] is also a two-stage evolutionary algorithm, which adjusts the fitness evaluation strategies during the evolutionary process to adaptively balance objective optimization and constraint satisfaction.
In addition, dual-population-based or coevolutionary technique has also been applied to CMOEAs or CMaOEAs. c-DPEA [8] is a typical dual-population-based evolutionary algorithm. Population1 uses an adaptive penalty function named saPF to handle infeasible solutions, and Population2 is oriented to feasibility. At the same time, an adaptive fitness function called bCAD is designed to better balance the convergence and diversity in the evolution process. Although c-DPEA can handle CMOPs well, it is inefficient in solving CMaOPs. C-TAEA [11] is an evolutionary algorithm with two archives. It maintains two collaborative archives simultaneously. The convergence-oriented archive (CA) aims to optimize constraints and objectives to ensure feasibility and convergence, while the diversity-oriented archive (DA) mainly tends to explore the areas that have not been exploited by the CA without considering feasibility. It is worth noting that C-TAEA is one of the few evolutionary algorithms designed for high-dimensional constrained objective spaces.

Proposed Algorithm: DP-NSGA-III
This section presents the proposed dual-population-based many-objective evolutionary algorithm, named DP-NSGA-III. The procedure of DP-NSGA-III is shown in Figure 1. DP-NAGS-III is a cooperative evolutionary algorithm that evolves two populations, namely, Population1 and Population2. The populations in DP-NSGA-III are evolved separately, only sharing all the offspring in each generation, where Population1 uses the ε-constraint handling method and Population2 ignores the constraint. Population2 can provide solutions for Population1 that violate the constraint but have excellent objective values to participate in evolution. The final output population of the algorithm is Population1. The ε-constraint handling method and the detailed evolution process of each population will be described below.

Constraint Handling
A CMOP can be mathematically defined as: where x is a decision vector in decision space Ω. F(x) is the objective vector that consists of M conflicting objective functions in objective space. g j (x) is the j-th inequality constraint.
A solution x is called a feasible solution in case g(x) ≤ 0; otherwise, it is an infeasible solution. In particular, when M > 3, the problem is called a constrained many-objective optimization problem (CMaOP), and the objective space is said to be high-dimensional.
In this article, the constraint violation is defined as the normalized violation of x on all constraints: where P denotes the initial population, and G i (x) represents the degree of the constraint violation of x on the i-th constraint: To avoid a population excessively weighted towards feasibility and to be able to involve infeasible solutions with excellent objective values in evolution, a parameter-free ε-constraint boundary function to deal with constraints is proposed, with the detailed formulation given as follows: where cv 0 (x) is the constraint violation of population initialization, and t is the number of iterations. The constraint boundary is a monotonically decreasing function with an initial value approximating cv 0 (x). As the population iterates, the constraint boundary value approaches 0. Therefore, at the end of evolution, this ε-constraint handling method is consistent with Deb's CDP, i.e., feasible solutions are preferred over infeasible solutions. Based on the above designed ε(t), the evolutionary algorithm solves the problem as shown below: The problem is equivalent to the original CMaOPs with relaxed constraints in the early stage. In the late stage, since g(x) → 0, it is equivalent to the original problem.
Based on the proposed ε-constraint handling method, the proposed evolutionary algorithm will solve Equation (1) indirectly. The main population solves the problem represented by Equation (5), and the auxiliary population solves the unconstrained problem corresponding to Equation (1).

Mating Selection and Offspring Reproduction
Population1 uses the ε-constraint handling method to deal with constraints. To make the offspring satisfy the ε-constraint boundary to the greatest extent, Population1 prefers an ε-feasible solution to an ε-infeasible solution, or a solution with a small constraint violation to a solution with a large constraint violation as the parent and puts it into the mating pool1. Algorithm 1 describes the above tournament selection procedure in detail.
Population2 does not consider constraints, so the nondominated solution is preferred as the parent and puts it into the mating pool2. If the two solutions do not dominate each other, then the selection is made by survival score [13].
Once the mating pool is constructed, the offspring of the population of size N is created using simulated binary crossover (SBX) [23] and polynomial mutation (PM) [24] operators. Finally, coevolution is achieved by merging Population1 (or Population2) with offspring 1 and offspring 2.

Environmental Selection of Population1
The environment selection of Population1 is based on NSGA-III combined with the ε-constraint handling method. This stage will select N outstanding individuals from Population1, offspring1, and offspring2 as new Population1.
The simplest case is that the number of ε-feasible solutions N f is less than N. Then, the new population consists of all ε-feasible solutions and N − N f solutions with the smallest constraint violation.
Otherwise, N f is more than N. The reference-point-based nondominated sorting with two steps is performed. First, different nondominated levels are classified by using the nondominated sorting based on M + 1 objectives: Note that the constraint violation is treated as a new objective in this process. Then, Second, reference-point-based elite selection is used to maintain population diversity. N − |F 1 ∪ F 2 ∪ · · · ∪ F k−1 | solutions are selected in F k based on the association between solutions and the predefined reference points (see Algorithm 2).
It is worth noting that, since the environmental selection of Population1 tends to ε-feasible solutions, the excellent infeasible solutions in offsping2 also have the possibility of entering the next generation. That is, the ε-constraint handling method provides the possibility of retaining the excellent offspring in coevolution.

Algorithm 2:
Reference-point-based nondominated sorting input : ε-feasible set S, reference points Z output : Population1 P t+1 1 Performing nondominated sorting on S based on Association Operation: 6 Associate each solution of P t+1 with the nearest reference point. The number of solutions associated with each reference point r is called niche count (r).

Environmental Selection of Population2
Population2 does not consider the constraint and eventually aims to converge to the UPF. Its environment selection [13] evaluates the solutions by p-norm, which can well distinguish the solutions in the high-dimensional objective space. This stage will select N outstanding individuals from Population2, offspring1, and offspring2 as new Population2. First, the nondominated levels are still divided by nondominated sorting. The following steps will be completely different from Population1. In each generation, Population2 is assigned fitness according to a fitness function called survival score. The survival score is composed of the ratio of the diversity score and the convergence score.
The distance of the solution to the real PF reflects the convergence. However, the real PF is unknown, so the convergence is usually expressed as the Euclidean distance from the solution to the ideal point Z min = (0, 0), i.e., the 2-norm of the solution. In Figure 2a, B is located at f 1 + f 2 = 1, C is located at f 1 + f 2 = 1.2, and A is located at the intersection of the red dashed line and the blue line. According to the Euclidean distance, since points A and B lie on the same circle, 1 = A 2 = B 2 > C 2 . However, the green line shows that A and C have the same distance to the real PF (assuming we know it), while B is the farthest. Therefore, in this case, 2-norm does not correctly represent the distance from the solution to the real PF. However, 1-norm can derive the correct distance relationship because A and C lie on 1-norm's contour f 1 1 + f 1 2 = 1.2 and B is on the outside of the contour, i.e., Therefore, a suitable p-norm is chosen to measure the convergence of the solution. Since the real PF is unknown, assume that the shape of the first front F 1 is the same as real PF and use the p-norm of F 1 to replace the p-norm of the real PF. A method to find the p-norm of F 1 is given below.
After normalizing the populations, the solution x = (x 1 , · · · , x D ) in the F 1 can be regarded as on the unitary hypersurface, satisfying the unique contour equation (x The function is second order derivable and can be solved by Halley's method [25] as follows: The p-norm of F 1 can be solved by substituting any solution in F 1 into Equation (6). Figure 2b was used to verify the correctness of Equation (6), which is the case of p = 2. D cos 7π 16 , sin 7π 16 , E cos π 4 , sin π 4 and F cos π 8 , sin π 8 are the normalized solutions in F 1 . Let the initial value p 0 = 1, then substitute D, E or F into Equation (6) to obtain p = p 2 ≈ 2 in only two iterations. p ≈ 2 is consistent with the actual situation, since the three solutions lie on the contour of the 2-norm.
The computational complexity of Halley's method is O(MT), where T is the number of iterations, generally around 3, much smaller than the size of the population N. Therefore, there is no impact on the time complexity of the algorithm.

Survival Score
Once the p-norm of the first nondominated front F 1 is calculated using Equation (6), we can measure both the convergence and diversity of F 1 [13] accordingly: The convergence score of an individual is the Minkowski distance of the individual to the ideal point Z min = (0, 0). Since the ideal point is the origin, the convergence is the p-norm of the individual. The diversity score of an individual is computed as the minimum Minkowski distance with the other solutions in the front F 1 . The survival score of each solution x ∈ F 1 consists of the ratio of diversity to convergence as follows: The procedure for calculating the survival score of the solution in F 1 is shown in Algorithm 3. Note that, if the solution x / ∈ F 1 , then its survival score is the reciprocal of the convergence score.

Discussion
Population2 has a different population distribution than Population1 because it ignores the constraint. On the one hand, Population2 can provide offspring with good diversity. In Figure 3a, Population1 is trapped in a local optimum, and its Offspring1 have difficulty finding other feasible areas. Population2 is widely distributed around the UPF due to ignoring constraints. Through the coevolution of offspring combinations, the better-diversified Offspring2 can lead Population1 to explore other feasible regions. On the other hand, Population2 can provide offspring with good convergence. In Figure 3b, Population1 falls into a local optimum away from the CPF and Population2, which does not consider constraints and explores towards the UPF. The Offspring2 near the CPF can pull Population1 near the CPF.

Time Complexity of DP-NSGA-III
The proposed ε-constraint handling method has no effect on the time complexity of NSGA-III, so the time complexity of Population1 is still O(MN 2 ) [9,26].
The time complexity of Population2 is influenced by nondominated sorting and the calculation of survival score. The time complexity of nondominated sorting is the same as Population1 with O(MN 2 ) [26]. The computational complexity of calculation of survival score of

Test Instances
In this article, C-DTLZ [20], DC-DTLZ [11], and MW [27] suites have been selected for examining the performance of DP-NSGA-III. The C-DTLZ suite is based on the DTLZ suite, and four constraints are added to the objective space, while the DC-DTLZ suite introduces three types of constraints in the decision space. The DC-DTLZ suite involves many local optima in infeasible regions, which causes the algorithm may easily fall into some infeasible local optima and cannot enter the feasible region. The objective number of the C-DTLZ suite and the DC-DTLZ suite is scalable. The MW suite was recently proposed, and it contains a variety of features, such as a small feasibility area, multiple complex nonlinear constraints, a scalable number of objectives, and high-dimensional decision vectors. Since some of the test functions in the MW suite have a fixed objective number, in this section, we only test MW4, MW8, and MW14 with variable objective space dimensions.
According to the relationship between UPF and CPF, the above test problems can be divided into three categories [20]: Type-I: The CPF is the same as the UPF. MW4, MW14, C1DTLZ1, C1DTLZ3, DC2DTLZ1, and DC2DTLZ3 are Type-I problems.
Type-II: Infeasible area makes the UPF partly feasible and the CPF is a part of the UPF. MW8, C2DTLZ2, DC1DTLZ1, and DC1DTLZ3 are Type-II problems.
Type-III: The UPF is completely infeasible and the CPF is located on the boundary of the feasible region. C3DTLZ4, DC3DTLZ1, and DC3DTLZ3 are Type-III problems.

Performance Metrics
To measure the performance of DP-NSGA-III and other peer CMaOEAs, the following two commonly used metrics are adopted. Furthermore, the Wilcoxon rank sum test with a significance level of 0.05 is adopted to perform statistical analysis.
(1) Inverted Generational Distance (IGD) [28]: Given P * as a set of points uniformly sampled along the PF and P as the set of solutions obtained from an evolutionary algorithm. The IGD value of P is calculated as where dist(x, P) is the Euclidean distance between x and its nearest neighbor in P.
(2) Hypervolume (HV) [29]: Let z r = z r 1 , · · · , z r m T be the worst point dominated by all the Pareto optimal objective vectors. The HV of P is defined as the volume of the objective space dominated by solutions in P and bounded by z r

Algorithms for Comparison
Four state-of-the-art and representative CMaOEAs were chosen for performance comparison: C-NSGA-III [20] (Abbreviated as NSGA-III), C-TAEA [11], DCNSGA-III [30], and MOEA/D-2WA [31]. NSGA-III and C-TAEA have been briefly introduced in Section 2. Thus, we briefly describe others' basic ideas as follows: DCNSGA-III: It is an improved NSGA-III that uses dynamic multi-objective techniques to solve CMaOPs. It is more suitable for processing CMaOPs than NSGA-III due to the upgrade of the constraint processing method. MOEA/D-2WA: It is a modified MOEA/D using two-type weight adjustments. During the search, the number of infeasible weights is dynamically lowered to help infeasible solutions with greater convergence transcend the infeasible areas, as well as to help infeasible solutions with better diversity identify many feasible sub regions.

Parameter Settings
(1) Algorithms and reproduction operators: To make a fair comparison and ensure the performance of each algorithm, the parameters of all the compared algorithm and their reproduction operator parameters, such as the simulated binary crossover and polynomial mutation, are set as suggested in their original papers. All the compared algorithms were implemented in PlatEMO [32]. (2) Population size and the number of function evaluations: Table 1 shows the number of reference points and population size of each algorithm for different numbers of objectives M [11,27]. The number of function evaluations is the termination condition of all algorithms [11,27]. The details are shown in Table 2.

Experimental Results
This section focuses on analyzing the performance of the proposed DP-NSGA-III with the peer CMaOEAs on three test suites. Note that the following population distribution plots or Parallel coordinate plots (PCPs) are taken for the median IGD of 30 independent experiments.

Performance on C-DTLZ Suite
Wilcoxon's test results on C-DTLZ suite are summarized in Table 3. According to the IGD values, the results show that DP-NSGA-III is significantly better than its competitors on at least 7/16 and at most 14/16 test problems. The results indicate that DP-NSGA-III performs much better than its rivals on at least 8/16 and at most 12/16 test problems, based on the HV values. Table 3. Wilcoxon's test results on C-DTLZ suite 1 .
Tables A1 and A2 show the mean value and standard deviation of the IGD and HV values obtained from the above-mentioned five algorithms over 30 runs on the C-DTLZ suite with 3, 5, 10, and 15 objectives. The C-DTLZ suite is moderately difficult, and the proposed algorithm achieves good results on all but C3DTLZ4.
For C1DTLZ1, only a narrow section near PF is feasible. Both the DP-NSGA-III and comparison algorithms converge on this problem. Figure 4 shows that differences in performance metric are mainly influenced by the uniformity of population distribution. C-TAEA has a messy population distribution and therefore has the worst performance metric.
The C1DTLZ3 problem has a highly multimodal landscape and infeasible region barriers that pose a challenge to the convergence of CMaOEA. Populations have difficulty crossing infeasible barriers and are easily trapped in localized feasible areas. The IGD of the proposed DP-NSGA-III is better than other peer algorithms by more than an order of magnitude in all four dimensions. The global optimal value of C1DTLZ3 is between 0 and 1. Figure 5 indicates that DP-NSGA-III converges to the global optimum except for a very small number of solutions. The other algorithms only converge to the boundary of the local feasible region and fall into the local optimum with a large number of solutions with objective values greater than 1. For DP-NAGA-III, the auxiliary populations that ignore the constraints help the main population to cross a huge infeasibility region. The CPF of C2DTLZ2 is part of the UPF, and the areas are not connected. This problem mainly tests the diversity maintenance ability of the algorithm. Most evolutionary algorithms found the CPF, and DP-NSGA-III performed relatively well.
The CPF of C3DLTZ4 is far from the UPF and is the boundary of the feasible area. In this issue, DP-NSGA-III was able to find CPF, but the population diversity was poor. Figure 6 shows the most uniform population distribution obtained by MOEA/D-2WA with two-type weight.

Performance on DC-DTLZ Suite
Wilcoxon's test results on DC-DTLZ suite are summarized in Table 4. According to the IGD values, the results show that DP-NSGA-III is significantly better than its competitors on at least 14/22 and at most 20/22 test problems. The results indicate that DP-NSGA-III performs much better than its rivals on at least 12/22 and at most 21/22 test problems, based on the HV values.
Tables A3 and A4 show the mean value and standard deviation of the IGD and HV values obtained from the above-mentioned five algorithms over 30 runs on the DC-DTLZ suite with 3, 5, 10, and 15 objectives. DC1DTLZ1 and DC1DTLZ3 are the Type-II problems. The CPF consists of multiple disjoint feasible areas. The obtained population is easily trapped in the local optimal region, and it is difficult to find the complete PF. Performance metrics show that DP-NSGA-III has a greater advantage over DC1DTLZ3 than DC1DTLZ1. Figure 7 illustrates that the population of DP-NSGA-III is evenly distributed on two separated PFs of DC1DTLZ1, while MOEA/D-2WA fails to find the PF above. DC2DTLZ1 and DC2DTLZ3 belong to the Type-I problems. Although it has the same CPF as UPF, the infeasibility region is huge and only the region near PF is feasible. Moreover, there exist many local optima in infeasible regions, making it much harder for solutions to escape. DP-NSGA-III achieved the best results in all four dimensions of DC2DTLZ1. Figure 8 demonstrates that DP-NSGA-III and MOEA/D-2WA perform best on DC2DTLZ3; in particular, the convergence and diversity of DP-NSGA-III are improved over NSGA-III and DCNSGA-III. NSGA-III and DCNSGA-III have objective values greater than 1 in multiple dimensions, indicating that there are solutions that fail to converge to the global optimum. Ignoring the few solutions with objective values greater than 1, the polyline in parallel coordinates for C-TAEA was not uniformly distributed compared to DP-NSGA-III, indicating an uneven distribution of the population.    DC3DTLZ1 and DC3DTLZ3 can be generated as Type-III problems, which are a mixture of the first two kinds of test instances. As the number of objectives increases, the feasible region becomes smaller and smaller. Note that no feasible solutions are found for all algorithms on the DC3DTLZ1 problems when M ≥ 10. Figure 9 shows that, except for DP-NSGA-III, all the algorithms fail to jump out of the infeasible region on DC3DTLZ3. DP-NSGA-III has the best convergence and distribution.

Performance on MW Suite
Wilcoxon's test results on MW suite are summarized in Table 5. According to the IGD values, the results show that DP-NSGA-III is significantly better than its competitors on at least 6/12 and at most 10/12 test problems. The results indicate that DP-NSGA-III performs much better than its rivals on at least 6/12 and at most 10/12 test problems, based on the HV values.
Tables A5 and A6 show the mean value and standard deviation of the IGD and HV values obtained from the above-mentioned five algorithms over 30 runs on the MW suite with 3, 5, 10, and 15 objectives.
MW4, MW8, and MW14 are the problems with the adjustable number of objectives in this suite. From Figure 10, we can see that DP-NSGA-III and MOEA/D-2WA perform comparably and better than other algorithms on the MW8 problem. NSGA-III and DCNSGA-III do not cover areas with a value of 1 in the first few dimensions. Different objectives of the MW14 have different ranges of value domains, with a maximum value of 5 when M = 15. Figure 11 shows that MOEA/D-2WA, which relies entirely on the reference vector, has degraded performance in this problem. When M ≤ 14, the population obtained by MOEA/D-2WA fails to find areas with an objective value of 1, while DP-NSGA-III still performs well in higher dimensions.

Stability of DP-NSGA-III
To evaluate the stability of each algorithm, box plots are used to represent the statistical results of IGD metrics for 30 experiments of the algorithms on each problem (M = 15). For convenience, each algorithm is abbreviated with its initials. Each subplot corresponds to the case of M = 15 in Tables A1, A3, and A5. Figure 12 shows that, while the other algorithms have more outliers and larger deviations, DP-NSGA-III has only a few outliers and small deviations on individual problems. DP-NSGA-III has the smallest inter-quartile range and the best stability for most problems. In general, DP-NSGA-III has obvious advantages in median, inter-quartile range and outliers, reflecting that DP-NSGA-III can better balance convergence and diversity on 15-dimensional CMaOPs and has better algorithm stability.

Conclusions
This paper proposed a dual-population-based constrained many-objective evolutionary algorithm, named DP-NSGA-III. Compared with existing CMaOEAs, DP-NSGA-III coevolves by exchanging the offspring of two populations. Two populations use different environmental selection to select superior individuals. Among them, the Auxiliary population deals with CMaOP that ignores constraints, and the main population focuses on solving the original CMaOP. The main population can make full use of the information of infeasible solutions in the auxiliary population to solve the constraint problem more effectively. From the problem perspective, the designed dual population solves the original CMaOPs indirectly. The main population solves the problem of loose constraints and the auxiliary population solves the problem of ignoring constraints. In addition, to improve the ability of the main population to handle constraints, this paper designed a parameter-free ε-constraint handling method. In the main populations driven by infeasible solutions, ε-feasible solutions can provide the impetus for evolution.
The performance of DP-NSGA-III has been studied for a series of suites with complex constraints with up to 15 objectives. Thanks to the proposed constraint handling methods that reduce the difficulty of handling constraints, the experimental results fully demonstrate the competitiveness of DP-NSGA-III in CMaOPs compared to the four state-of-the-art CMaOEAs, especially in DC-DTLZ suite with highly complex constraints. However, if the CPF is not a regular surface, or if the UPF is separated from the CPF, DP-NSGA-III performs poorly. On C3DTLZ4 with both features, DP-NSGA-III converged to CPF, but the populations were not as well distributed as NSGA-III. The solutions provided by the auxiliary population have a side effect instead. In the future, it is necessary to further explore the underlying mechanisms of dual-population and design better dual-population strategies. In terms of environment selection, more efficient scalar methods need to be designed to replace Pareto dominance relationships that make it difficult to distinguish individuals in a high-dimensional objective space.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:  Table A1. IGD values obtained by DP-NSGA-III and other peer algorithms in the C-DTLZ suite 1 .  1 "+", "−", or "=" indicates that the corresponding algorithm is significantly better than, worse than, or comparable to DP-NSGA-III. The best average value among these algorithms on each test problem is highlighted in bold. Table A2. HV values obtained by DP-NSGA-III and other peer algorithms in the C-DTLZ suite 1 . +/−/= 4/11/1 3/11/2 2/12/2 6/8/2 1 "+", "−", or "=" indicates that the corresponding algorithm is significantly better than, worse than, or comparable to DP-NSGA-III. The best average value among these algorithms on each test problem is highlighted in bold. Table A3. IGD values obtained by DP-NSGA-III and other peer algorithms in the DC-DTLZ suite 1 . 1 "+", "−", or "=" indicates that the corresponding algorithm is significantly better than, worse than, or comparable to DP-NSGA-III. The best average value among these algorithms on each test problem is highlighted in bold. Table A4. HV values obtained by DP-NSGA-III and other peer algorithms in the DC-DTLZ suite 1 .  Table A6. HV values obtained by DP-NSGA-III and other peer algorithms in the C-DTLZ suite 1 . Problem  M  NSGA-III  C-TAEA  DCNSGA-III  MOEA/D-2WA  DP-NSGA-III   MW4