Constrained Optimization Based on Ensemble Differential Evolution and Two-Level-Based Epsilon Method

Constrained optimization problems (COPs) are common in many fields, and the search algorithm and constraint handling technique play important roles in the constrained evolution algorithms. In this article, we propose a new optimization algorithm named CETDE based on ensemble differential evolution (DE) and a two-level epsilon-constrained method. In the ensemble DE variant, some promising parameters and mutation strategies constitute the candidate pool, and each element in the pool coexists throughout the search process and competes to generate new solutions. The two-level epsilon method is proposed by incorporating a generation and a population comparison level to retain more promising solutions without degrading the solution quality. Moreover, a diversity promotion scheme is developed to improve the population distribution when the search becomes trapped in a small region. The superior performance of CETDE is validated by comparison with some state-of-the-art COEAs over two sets of artificial benchmarks and five real-world problems. The competitive results show that CETDE is an effective method for solving COPs.


I. INTRODUCTION
Constrained optimization problems (COPs) are frequently encountered in many real-world design domains. Without loss of generality, a COP can be written as [1], [2]: 1 , · · · , x n ] T s.t. g j ( x) ≤ 0, (j = 1, 2, · · · , p) h j ( x) = 0, (j = p + 1, · · · , q) a i ≤ x i ≤ b i , (i = 1, 2, · · · , n) (1) where n is the number of the dimensions of the decision variable x, a = [a 1 , · · · , a n ] T and b = [b 1 , · · · , b n ] T are the lower and upper bounds, respectively, f ( x) is the objective function to be optimized, and g( x) and h( x) are the p inequality and q − p equality constraints, respectively. For COPs, The associate editor coordinating the review of this manuscript and approving it for publication was Juan Liu .
the constraint violation of solution x is defined as: where δ is a positive tolerance value for transforming the equality constraints into inequality constraints. We say that x is feasible if φ( x) = 0, otherwise, x is infeasible. Currently, evolutionary algorithms (EAs) are frequently adopted to solve optimization problems in various domains because of their advantages over traditional mathematical methods [3]- [8]. However, we should be aware that EAs were originally proposed for unconstrained global numerical optimization problems, and the main goal of constrained optimization algorithms is finding the global optimal solution in the feasible domain. Thus, to solve COPs effectively, the powerful EAs should be combined with suitable constraint handling techniques (CHTs) in a proper way [9]- [11].
To date, many CHTs have been designed for natureinspired search algorithms, and each has its own VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ merits [12], [13]. Based on the methods used to treat the constraints, we can group the existing CHTs into three different categories, i.e., penalty function (PF)-based methods, superiority of feasibility-based methods and multiobjective-based methods.
The main idea of the PF-based methods is to first transform the original problem into an unconstrained problem by adding a penalized constraint violation to the objective function so that any existing search algorithm can be adopted to solve the unconstrained problems. These methods are easy to use but have also been criticized because the overall performance depends strongly on the associated penalty factor. To avoid this limitation, advanced PF-based methods have been extensively studied. For example, Runarsson and Yao proposed stochastic ranking method which does not require a priori information about COPs and also does not use any additional penalty factors [14]. Gan et al adopted a dynamic exponential function-like penalty factor setting to convert COPs into unconstrained ones and then the simplex crossover-based EA was used as the search algorithm [15]. Wang and Cai proposed an adaptive trade-off model, which the feasible rate of the main population is utilized [9]. Tessema and Yen proposed an adaptive PF method to exploit infeasible solutions with low objective value and low constraint violation [16].
The main principle of the superiority of feasibility-based methods is to prefer the feasible solutions to the infeasible solution. These methods usually do not require any additional parameters, making them quite simple and easy to use [17]. However, these methods often display heavy selection pressure, which is likely to lead to premature convergence, particularly for COPs with complex constraints. To avoid this problem, some improvements have been proposed to alleviate the selection pressure. For example, Takahama and Sakai introduced new CHT named ε-constrained methods by using a relax comparison level in feasibility rule. This method treats some solution with low constrained violation as feasible ones [18]. Fan et al proposed a feasible ratio based epsilon level control method. This method is then incorporated into decomposition-based multiobjective EA to solve constrained multiobjective problems [19]. Zhang et al proposed a adaptive ε-method based on the knowledge of constraint violation in the main population. Meanwhile, a more simpler heuristic rule is designed control the comparison level instead of traditional exponential function-based methods [20]. To solve COPs effectively, Wang et al first adopted the common feasibility rule during the selection procedure and then applied an individual replacement technique to preserve some infeasible solutions with low objective values [21]. Xu et al proposed a cluster-replacement-based feasibility rule to alleviate the selection pressure behind the ordinary feasibility rule [22]. Wang et al proposed an individual-dependent feasibility rule to alleviate the preference from both the objective function and constrained violation [11].
The multiobjective-based methods usually introduce some additional objective function to form a new problem and then employ the concept in the multiobjective optimization community to solve it [23]- [25]. However, the application of these methods still faces some challenges because solving multiobjective optimization problems is still quite challenging in and of itself. For this kind of method, the key issue is how to use the multiobjective concept to balance the search between objective functions and constraints. For example, Cai and Wang tried to combine the original objective function and constrained violation together to form a bi-objective problem, i.e. [f ( x), G( x)] and then the domination concept in multiobjective optimization community is applied to tackle the new problem [23], [24]. Deb and Datta combined the elitist non-dominated sorting genetic algorithm with PF method and proposed a hybrid EA-cum-CHT in a manner complementary to each other [26]. In addition to using static bi-objective structures, some dynamic constrained multiobjective framework is also proposed. Zeng et al converted COP into an equivalent dynamic constrained multiobjective optimization problem with 3 objectives [27], i.e. original objective, constrained violation and niche-count objective. During the search process, the constraint boundary and niche size are dynamically reduced. The new framework is incorporated into Pareto ranking-, decomposition-and hypervolume indicator-based multiobjective evolutionary algorithms, showing some advantages. Later, Jiao et al proposed an improved dynamic constrained multiobjective framework by the Pareto dominance and a new feasible-ratio control mechanism [28]. One remarkable superiority of this new framework is allowing the search from both feasible and infeasible domains. Zeng et al proposed Pareto front transformation model to retain some promising dominated feasible vectors and nondominated infeasible ones with worse objective. Thus the search within both the feasible and infeasible region is allowed [25].
In addition to designing effective CHTs, study of powerful search algorithms is another important issue in the development of competitive COEAs. Many search algorithm are available to the EA community, such as the genetic algorithm (GA), particle swarm optimization (PSO), artificial bee colony algorithm (ABC), and teaching-learning-based optimization (TLBO). Among these types of algorithms, differential evolution (DE) is a simple yet effective search algorithm over a continuous space, and a large number of DE-based methods have been proposed to handle COPs [29], [30]. The simplest approach for the use of DE for COPs is to directly combine basic DE with various CHTs. However, the performance of DE strongly depends on the mutation and crossover strategies and the associated control parameters, i.e., the population size N , scale factor F and crossover rate CR. Thus, these strategies and parameters must be determined carefully and be appropriate for a specific COP. As mentioned above, DE was originally proposed for unconstrained problems. To enhance the search ability of DE when solving complex COPs, various kinds of approaches were proposed, such as the design of new mutation operator [10], [31]- [33], study of advanced parameter control strategies [10], [22], [34], [35], and introduction of additional diversity schemes [11], [21], [36].
In fact, there has been increasing interest in enhancing the search ability of DE by using multiple strategies and parameters. For these multi-strategy-based variants, each has its own features, For example, Wang et al proposed a composite DE variant by using three different offspring generation strategies and three control parameter settings according to some other researchers' studies [37]. One distinguishing feature of this method is three trial vectors are generated by three different trial vector generation strategies for each target vector in the main population and the best among the three will be used to update the corresponding parent vector during the selection process. Inspired by this framework, Wang et al proposed a constrained composite DE variant to solve COPs [38]. Xu et al proposed a cooperative ranking-based mutation strategy in DE for COPs [33]. Although the commonly used DE/rand/1 mutation strategy is adopted, it uses two different selection technique to determine parent vectors at two different situations. Sun et al proposed new DE variant by adopting two mutation strategies with different features [39]. One strategy is focusing on global exploration while the other on local exploitation. Among all these multi-strategy-based DE variants, ensemble of parameters and mutation strategies is a simple, yet powerful DE variant for global numerical optimization [40]. This approach attempts to provide a systematic framework that combines multiple mutation strategies with some parameter values. In this DE variant, some candidate pools are provided beforehand, and learning information is extracted for further processing. This scheme has been proved to be very effective when solving unconstrained global numerical optimization. How to extend it for COPs needs to be further studied.
The ε-constrained method proposed by Takahama and Sakai [18] seeks to relax the comparison in the feasibility rule and introduces a new comparison level, ε. If the constrained violation values of solutions are less than the ε level, they are considered to be feasible, and the comparison is based solely on the objective function values. Otherwise, the solutions are treated as infeasible, and the solution with the small constrained violation is deemed to be better. By introducing this relaxed ε level, the search algorithms will have more opportunities to explore the infeasible region around the feasible domain, which will be very helpful in the case where the feasible optimal solution lies at the boundary of constraints. However, the original ε-constrained method shares the same comparison level for all solutions in the main population. If the ε level is too large, the search will focus on minimizing the objective functions, degrading the exploitation around the feasible solution. By contrast, if the ε level is too small, more attention will be focused on decreasing the constraint violation.
Motivated by the above considerations, in this article, we propose a new constrained DE algorithm named CETDE by combining an ensemble DE variant with a new twolevel-based ε method to solve COPs. In CETDE, some well-studied mutation strategies and parameter values based on previous studies performed by other researchers are collected to form some candidate pools. If a combination of strategy and parameter can generate promising solutions in a previous search, it will have a greater likelihood to be used in the subsequent evolution process stage. Moreover, a diversity promotion scheme is used to improve the population diversity when the search is trapped into a narrow domain.
The remainder of this article is organized as follows. Section II briefly introduces the basic DE algorithm and the ε-constrained method. Section III presents our proposed CETDE algorithm in detail. Section IV presents the experimental results and performance comparisons. Finally, Section VI draws conclusions

II. DIFFERENTIAL EVOLUTION AND EPSILON-CONSTRAINED METHOD A. DIFFERENTIAL EVOLUTION
The traditional DE begins with randomly initializing a population with N n-dimensional vectors, i.e., 1 , · · · , x g i,n ] g is the ith target vector and g is the generation number. Then, DE employs mutation, crossover and selection operators sequentially and iteratively until a termination condition is met [29], [41].
The mutation operator generates a mutant vector for every target vector. The following three operators are used in this article: where indices r 1 , r 2 , r 3 , r 4 and r 5 are mutually exclusive integers randomly generated within [1,N ] that are also different from index i, x best is the current found best solution, F is the scale factor lying between 0 and 1, and K is a random number The crossover operator generates a trial vector u g i using the target vector x g i and mutant vector v g i . The most frequently used binomial crossover strategy is described by: where sn = rndint(1, n) is an arbitrary number in {1, 2, · · · , n} and u g i,j is the jth element of the ith trial vector. Parameter CR is called the crossover rate and is a real number between 0 and 1.
After generating u g i , DE uses a greedy selection strategy to form the target vector x g+1 i at generation g + 1: VOLUME 8, 2020

B. THE EPSILON-CONSTRAINED METHOD
The ε constrained method is a commonly used CHT. The main idea of this method is to treat some infeasible solutions as feasible solutions during the early stage based on a predefined comparison level ε controlled by Eq. (8). Meanwhile, the value of ε is gradually decreased until the generation number reaches T c .
where ε 0 is the initial comparison level, x θ is the top θ vector and θ = 0.2N , and cp ∈ [2, 10] is used to control the decreasing speed.
Let f 1 (f 2 ) and φ 1 (φ 2 ) be the objective value and constraint violation at solution x 1 ( x 2 ). Then, the ε-level comparison, denoted as ≤ ε , between x 1 and x 2 is given by: It can be observed from Eq. (9) that for ε = 0, ≤ 0 is equivalent to the feasibility rule [17]; by contrast, for ε = ∞, ≤ ∞ becomes the pure comparisons between the objective values.

III. PROPOSED CETDE A. BASIC IDEA OF CETDE
Many researchers have pointed out that the search algorithm and constraint handling technique (CHT) are highly important for solving COPs [9]- [11], [33], [42]. Thus, improvement in the search ability and development of effective CHTs are commonly considered when solving COPs. In our CETDE, DE with the ensemble of mutation strategies and parameters is adopted as the main search engine, and the proposed two-level ε-constrained method is used as the main CHT. Moreover, a diversity promotion scheme is used to improve the population diversity. The framework of CETDE is given in Algorithm 1. Next, we will describe these three key components of CETDE in detail.

B. ENSEMBLE DIFFERENTIAL EVOLUTION FOR COPs
The effectiveness of DE in solving a given problem strongly depends on the mutation strategy and its associated parameter values. While many mutation strategies are available, none of them can achieve high overall performance for all problems with different characteristics. Thus, DE with multiple mutation strategies and parameters is popular. To exploit the concept of multiple strategies, ensemble DE [40] is applied to design the search algorithm when solving the COPs. The main idea of ensemble DE is that three candidate pools, i.e., the mutation strategy pool (Mpool), scale factor pool (Fpool) and crossover rate pool (CRpool), are introduced. Generally, we expect that the elements presented in the each pool should have distinct characteristics and Algorithm 1 Pseudocode of CETDE Algorithm 1: Set the parameters: population size N , mutation strategy pool Mpool, scale factor pool Fpool, crossover rate pool CRpool, the maximize fitness evaluation number MaxFES. 2: Set the generation number g = 1 and randomly generate a population with N vectors uniformly in the Evaluate the objective function f ( u t i ) and constraint violation φ( u t i ). for i = 1 : N do 13: Get the ranking value r(i) for ith vector. 14: Calculate the ε g i comparison level for ith vector at gth generation using Algorithm 2. 15: Randomly select M g i , F g i and CR g i from the corresponding candidate pools or a combination from the stored successful combination pool. 22: end for 23: Perform the diversity promotion scheme using Algorithm 3. During the evolution stage, little learning information is available at the first stage, and much information can be used at the later stage. Thus, we randomly select the mutation strategy and parameters from the candidate pools for every vector at g = 0. Suppose that the mutation strategy, scale factor and crossover rate associated with vector x and CR g i may be not suitable, and they can be re-selected. In this case, we randomly select the mutation strategy and parameters from the corresponding candidate pools or from the successful combination pool with equal probability.

C. TWO-LEVEL EPSILON-CONSTRAINED METHOD
For the ordinary ε-constrained method, all of the vectors in the main population share the same comparison level at generation g. This configuration will introduce some more infeasible vectors to promote the diversity of main population, but it will be highly likely to reduce the convergence speed and degenerate the current feasible optimal solution. It will also most likely cause the search to become trapped in the infeasible region so that no feasible solutions will be obtained at the end of the search.
To overcome these drawbacks to some extent, we proposed a new two-level ε method. The main idea of this method is that two comparison levels, named the generation level and population level here, are introduced. The generation level functions as the ordinary ε method. For the population level, we hope that different vectors have different comparison levels. To achieve this goal, we first sort all the vector from the best to worst according to the following criteria: (1) Feasible vectors are better than infeasible ones; (2) Ranking feasible vectors according to their function values; (3) Ranking infeasible solutions based on their constraint violations. Based on the sorting values, we can obtain the ranking value for each vector x g i as r(i) = N − i, i = 1, 2, . . . , N . And then the combined level can be calculated using a calculation model [33]. Many models are available, and the simplest linear model is described by: where ε g is the generation level calculated by Eq. (8); ε g i is the combined level for ith vector at generation g; r(i) is the ranking value for the ith vector x g i in the main population. It can be observed from Eq. (10) that each vector has its own comparison level value and the superior ones have rigorous levels while the inferior ones own some relatively relaxed level. The main advantage of this idea is that the algorithm can retain much better solutions with respect to both objective function and constrained violation around the superior solutions. It will also make the search to focus on minimizing objective function values for the inferior vectors. Thus the combined two-level approach can improve the search algorithm's exploration and exploitation abilities.
To clearly illustrate the difference between the ordinary and two-level ε method, Fig. 1 provides some additional information for these two methods with N = 10 and T c = 20. From the perspective of generation level, both ordinary and two-level epsilon methods share the same parameter settings, as shown in Eq. (8) and Fig. 1(a). However, from the perspective of population level, these two methods are different. Supposing that the ith vector at generation g is ε g i , for the ordinary epsilon method, all vectors have equal ε values, i.e. ε g 1 = ε g 2 = · · · , = ε g N , at every generation, as illustrated in Fig. 1(b). By contrary, the two-level epsilon method allow each vector to have different ε values, i.e. ε g 1 = ε g 2 = · · · , = ε g N . More specifically, the best vector has lowest ε value and the worst has the largest value, as illustrated in Fig. 1(d). When both the population and generation levels are considered, we get the combined level for ordinary and two-level ε method, as illustrated in Fig. 1(c) and Fig. 1(e), respectively. The main procedure of the two-level ε method is given in Algorithm 2.

Algorithm 2 The Two-Level ε Method
Require: rank of the vector: r(i), the current generation number g, inital epsilon value ε 0 . Ensure: The epsilon value for vector x g i , i.e. ε g i . 1: Calculate the epsilon generation level value ε g at gth generation using Eq. (8). For some COPs with extremely complicated objective function or constraints, the feasible region becomes highly nonlinear and always exhibits a multimodal property. Thus, the search is highly likely to become trapped in the region of a local minimum. Introduction of new solutions into the population to enhance the population diversity is a common approach to alleviate this problem. Many mechanisms can be applied to achieve this goal, such as the mutation operator [9], [21], restart scheme [38], and solution replacement [22] techniques. To introduce new solutions, three issues must be considered, i.e., when to introduce new solutions, which solution should be introduced and how to introduce the solutions.  Intuitively, if all of the vectors are located in a very small search region, we can say that premature trapping occurs, and we can introduce some new solutions to enhance the population diversity. However, the measurement of the similarity of all of the vectors in the main population is still a tough problem. Here, a simple indicator, i.e., standard deviations of both constraint violations and objective values, is used to evaluate the similarity of a population. Based on the standard deviation values, some infeasible trial vectors with low objective functions are conditionally maintained, and then, a greedy selection strategy is applied to replace some target vectors in the population. The diversity scheme used in this article is given in Algorithm 3. In Algorithm 3, threshold η is applied to decides whether to use the diversity promotion scheme or not. If the η is too big, the diversity promotion scheme will be frequently applied and infeasible region with low objective value will be located. On the contrary, if η is too small, the diversity promotion scheme will be called rarely, which will lead to the failure of increasing population diversity. In fact, we have made some experiments with different values (i.e. 1 × 10 −4 , 1 × 10 −5 and 1 × 10 −6 ), on problems from IEEE CEC2006, and we did not notice any significant difference in results. Therefore, we arbitrarily peak 1 × 10 −4 in this article.

IV. EXPERIMENTAL STUDY A. TEST PROBLEMS AND EXPERIMENTAL SETTINGS
To evaluate the performance of CETDE, two well-known benchmark sets from IEEE CEC2006 and CEC2018 are selected as the test instances. The first set contains 24 problems that involve 5 types of objective functions (linear, nonlinear, polynomial, quadratic, and cubic), with different numbers of design variables (between 2 and 24) and 4 kinds of constraints (linear inequalities, linear equalities, nonlinear inequalities, and nonlinear equalities) with different numbers of constraints (between 1 and 38). The second set contains 28 problems with a wide variety of constraints and dimensions. The mathematical formulas and properties of these functions can be found in [1] and [43], respectively.
For all experiments, the positive tolerance value δ is set to 0.0001 to transform the equality constraints into inequality constraints. For the ε-constrained method, T c is set to 0.2T max , and cp is set to 6.

B. GENERAL PERFORMANCE OF CETDE ON IEEE CEC2006 PROBLEMS
In this section, the general performance of CETDE is tested over 24 problems collected in CEC2006. For this experiment, the population size N is set to 50 for all 24 test instances. First, we execute CETDE 25 times and record the optimal function error values, i.e., (f (x) − f (x * )), after 5 × 10 3 , 5 × 10 4 and 5 × 10 5 function evaluations, respectively. Then, we calculate some statistical indicators in terms of 'Best', 'Median', Worst, c, 'v', 'Mean' and standard deviation ('Std'), where c is a sequence of three numbers denoting the number of violated constraints for the median solution that are greater than 1.0, between 0.01 to 1.0, and between 0.0001 to 0.1, v is the mean value of all constraint violations for the median solution, and the numbers in the parenthesis following 'Best', 'Median' and 'Worst' are the numbers of constraints unable to satisfy the feasible condition for the best, median and worst solutions, respectively.
For this experiment, the optimal function error values after 5 × 10 5 FES over 25 runs are summarized in terms of 'Best', Median, Worst, Std, feasible rate, success rate and success performance values in Tables 5 and 6. As shown in Tables 5  and 6, the obtained best values for all 24 instances by CETDE on 25 independent runs are almost equal to the known optimal value f ( x * ), except G20 and G22. Since G20 and G22 are extremely hard and no feasible solutions have been found so far, we do not report the statistical values in these tables. It is observed that CETDE obtained a 100% feasible rate on all 22 instances. For the success rate, CETDE achieved 100% on 21 out of 22 instances. In particular, CETDE obtained a 96% success rate on G02.
The experimental results of the CETDE on 24 instances from IEEE CEC2006 are shown in Table 6, where the known optimal value and the achieved Best, Median, Worst, Mean and Std values according to 25 independent runs are presented. It is found from the results presented in this table that CETDE can find the optimal solutions in all 25 runs consistently for all of the instances, except G02. Although CETDE cannot consistently achieve the optimal value on G02, the obtained values in 25 runs are quite close to the known optimal value, as verified by the small Std value. Thus, we deduce that the proposed CETDE is an effective and robust method for COPs.  To provide additional information about the CETDE's search ability when solving instances from IEEE CEC2006, we show the convergence graphs of the function error values, i.e., lg(f ( x) − f ( x * )), versus fitness evaluations (FES) at the median run over 25 runs in Figs. 2-5. Since CETDE cannot find any feasible solutions on G20 and G22 over 25 runs, we do not show the convergence graphs of these two instances in the figures. It is observed from these figures that for most instances, the function error values decrease dramatically at the early stage and evolve slowly at the later stage. This is because the obtained values are quite close to the known optimal value. Based on these figures, we conclude that the proposed CETDE shows relatively fast convergence for most instances. Table 7 gives the mean and standard deviation of objective values with respect to CETDE and nine competing algorithms. Among all nine competitors, the first four are DE-based methods (i.e., DPDE [31], rank-iMDDE [36], eDE [44], and FRC-CEA [28]), the next four are non-DE-based methods (i.e., AIS [45], CMPSOWV [46], I-ABC [47] and ITLBO [48]), and the last one is a PSO-and DE-based hybrid method (i.e., DPD [49]). Since the best feasible optimal values are known, we compare the mean function values with the best known optima. If the obtained mean function error value is less than a predefined threshold, i.e., f ( x) − f ( x * ) ≤ 10 −4 , we denote it using symbol in Table 7. Because the optimal   feasible solutions of G20 and G22 are unknown, Table 7 does not contain the statistical information on these two instances. We note that the maximal function evaluation number for this experiment is set to 2.4×10 5 and that the Mean and Std values provided by all competing algorithms are taken directly from the corresponding references for fair comparison.

C. COMPARISON TO OTHER STATE-OF-THE-ART COEAs
As shown in Table 7, the hybrid DPD method ranks first on all ten algorithms because it obtained 20 successful runs out of 22. The proposed CETDE and other three competing algorithms (i.e., DPDE, eDE and FRC-CEA) rank second because they only obtained 19 successful runs out of 22. The non-DE-based algorithms (i.e., AIS, CMPSOWV, I-ABC and ITLBO) perform slightly worse than the DE-based approaches, obtaining 15, 16, 7 and 15 successful runs out of 22, respectively. Based on the results of this experiment, we conclude that CETDE either outperforms or performs similarly to most competitors on the 22 test instances from IEEE CEC2006.
Based on the Mean indicator summarized in Table 7, the average ranks obtained by applying the Friedman procedure are shown in Fig. 6. We can find that CETDE ranks fourth among all ten algorithms, which means that CETDE performs worse than DPDE, rank-iMDDE and FCR-CEA and better than the rest six COEAs on IEEE CEC2006 problems.

D. COMPARISONS WITH COEAs ON IEEE CEC2018 PROBLEMS
To further test the performance of CETDE when solving more challenging instances, in this section 28 instances selected from IEEE CEC2018 are used as test problems. Four different dimension levels, i.e. 10, 30, 50, 100, are applied in the experiment. Here we set the population size N to 50, 50, 80 and 100 when the dimension level is 10, 30, 50 and 100, respectively. Still 25 independent runs are carried out for every problem at each dimension level. Table 8 summarized the mean objective values, mean constraint violation values and feasibility rate for 28 problems over 25 runs.
Similarly, we choose three algorithms, i.e. MAg-ES [50], IUDE [51] and LSHADE [52], as the competitors because these algorithms achieved very competitive results in the IEEE CEC2018 competition. Because the optimal values of these problems are not provided in [43], we compare the  performance of all algorithms based on the mean objective values over 25 runs. We obtained the experimental results regarding to competitors from their corresponding original references. The multi-problem Friedman's test is carried out by the KEEL software [53] and the results are shown in Fig. 7. It can be found that CETDE performs worse than IUDE and better than MAg-ES at 4 different dimension levels. When compared with LSHADE, CETDE performs worse at D = 100 and performs better at other three dimension levels. Overall, CETDE achieves very competitive results when solving problems from IEEE CEC2018.   The statistical results for all problems are summarized in Table 9. Please note that the final results of the competing methods are directly taken from the corresponding references and the symbol '-' means the corresponding values are not available. It can be found from Table 9 that for most problems, CETDE obtained very similar results with other competing methods and CETDE is a potential method to deal with constrained real-world problems.

V. ADDITIONAL DISCUSSION
In our CETDE, three main components contribute to obtaining high performance when solving COPs, i.e., the multistrategy-based DE variant, two-level epsilon method VOLUME 8, 2020  and diversity enhancement scheme. To study whether these components are essential and whether they indeed can improve the search ability of DE, in this section, we provide some additional experiments to demonstrate the effectiveness and rationality of these components in CETDE.

A. EFFECTIVENESS OF THE ENSEMBLE DE METHOD
The proposed CETDE contains three different mutation operators when solving COPs, and here, we design three additional methods to test the necessity of using multiple strategies. To facilitate the experiments, we denote these three methods as CDE-1, CDE-2 and CDE-3 as follows: (1) CDE-1: DE using strategy ''DE/rand/2''; (2) CDE-2: DE using strategy ''DE/current-to-rand/1''; (3) CDE-3: DE using strategy ''DE/current-to-best/1''; We also run each method 25 times and compute the Mean, Std, feasible rate (FR) and success rate (SR) values after 5 × 10 5 FES. The experimental results are presented in Table 10. It should be noted that all of the methods obtained similar results on 11 instances (i.e., G01, G04, G05, G06, G08, G09, G11, G12, G15, G16 and G24) and could not find any feasible solutions on two instances (i.e., G20 and G22), and we do not show the results in Table 10 for these instances. An examination of the data presented in Table 10 shows that CETDE achieves better results than the other three methods in terms of most indicators, particularly for FR and SR. Based on this study, we conclude that constrained DE with multiple strategies performs better than DE with a single mutation strategy.
All algorithms are carried out 25 times, and the Mean, Std, FR and SR values after 5 × 10 5 FES are also summarized in Table 11. Here we do not present the results on 18 instances (i.e., G01, G04, G05, G07, G08, G09, G11, G12, G14, G15, G16, G18, G19, G20, G22, G23 and G24) because all methods achieve similar results. It is found from Table 11 that CETDE obtains slightly better/comparable results than those achieved by other three variants, showing that the two-level epsilon method is an effective CHT for solving COPs.

C. EFFECTIVENESS OF THE DIVERSITY PROMOTION SCHEME
In CETDE, a diversity scheme is proposed to enhance the population diversity. To experimentally validate this scheme's effectiveness, we implemented a new method by removing this scheme. We denote this method as CDE-6 and compare it with CETDE on 24 instances over 25 runs. Table 12 gives the Mean, Std, FR and SR values after 5 × 10 5 FES. We note that CDE-8 and CETDE obtained different results on 3 out of 24 instances. Thus, only these three instances, i.e., G02, G13 and G17, are summarized. It is found that CETDE performs better than CDE-6 in term of most indicators, particularly for the SR values. Thus, the proposed diversity scheme can improve the search ability of DE.

VI. CONCLUSION
The performance of a COEA strongly depends on the search algorithm and CHT. In this article, we proposed a new COEA, denoted as CETDE, to solve COPs by combining an ensemble DE variant and an improved two-level epsilon method. In CETDE, both the control parameters and mutation strategies with various characteristics are introduced into the DE algorithm. All the parameters and strategies coexist throughout the search process and compete to generate new solutions. We also introduced a generation and a population comparison level into the original epsilon method to retain more promising solutions without degrading the solution quality. Moreover, a diversity scheme is proposed to increase the population diversity when the search becomes trapped in a local domain.
The proposed CETDE method is compared with some state-of-the-art COEAs on two sets of artificial benchmarks and 5 mechanical engineering design problems. Extensive and systematic studies show that the proposed CETDE achieved better or comparable results on these selected problems. The effectiveness and benefits of the ensemble DE variant, two-level epsilon method and diversity scheme are experimentally validated by some additional studies. The experimental results demonstrate the effectiveness of these components for solving COPs.
In the future, we would like to further enhance the search ability of CETDE for more complex COPs. Another interesting research direction would be to extend our DE variant, two-level epsilon method and diversity scheme to constrained multiobjective optimization problems.
[58] V. V. de Melo and G. L. C. Carosio, ''Investigating multi-view differential evolution for solving constrained engineering design problems,'' Expert Syst. Appl., vol. 40 ZONGHAO ZHANG received the B.S. degree from Tongling University, in 2018. He is currently pursuing the M.S. degree with the Shanghai University of Engineering Science. His current research interests include intelligent computing algorithms and their application in path planning. VOLUME 8, 2020