An Ensemble Framework of Evolutionary Algorithm for Constrained Multi-Objective Optimization

: In the real-world, symmetry or asymmetry widely exists in various problems. Some of them can be formulated as constrained multi-objective optimization problems (CMOPs). During the past few years, handling CMOPs by evolutionary algorithms has become more popular. Lots of constrained multi-objective optimization evolutionary algorithms (CMOEAs) have been proposed. Whereas diﬀerent CMOEAs may be more suitable for diﬀerent CMOPs, it is diﬃcult to choose the best one for a CMOP at hand. In this paper, we propose an ensemble framework of CMOEAs that aims to achieve better versatility on handling diverse CMOPs. In the proposed framework, the hypervolume indicator is used to evaluate the performance of CMOEAs, and a decreasing mechanism is devised to delete the poorly performed CMOEAs and to gradually determine the most suitable CMOEA. A new CMOEA, namely ECMOEA, is developed based on the framework and three state-of-the-art CMOEAs. Experimental results on ﬁve benchmarks with totally 52 instances demonstrate the eﬀectiveness of our approach. In addition, the superiority of ECMOEA is veriﬁed through comparisons to seven state-of-the-art CMOEAs. Moreover, the eﬀectiveness of ECMOEA on the real-world problems is also evaluated for eight instances.


Introduction
Constrained multi-objective optimization problems (CMOPs) widely exist in real-world applications and scientific research fields, such as the multi-objective caching optimization [1], vehicle routing problems [2], web service location allocation problems [3], 3D indoor redeployment in IoT collection networks [4], asymmetric UAV trajectory planning [5], multiobjective testing resource allocation problems [6], and ring road bus lines and fare design [7]. Generally, CMOPs can be referred to as the optimization problems with constraints and two or more conflicting objectives [8], which can be mathematically formulated as: where m denotes the number of objective functions, x = (x 1 , . . . , x n ) T represents an ndimensional decision vector, n denotes the number of decision variables, x ∈ S and S ⊆ R n are the search space. g i (x) and h j (x) are the i-th inequality and j-th equality constraints, respectively. q is the number of constraints. The constraint violation (cv j (x)) of x at the j-th constraint can be formulated as: cv j (x) = max(0, g j (x)), j = 1, · · · , p max(0, |h j (x)| − σ), j = p + 1, · · · , q ,

CMOEAs Based on Constrained Dominance Principle
The first category adopts constrained dominance principle (CDP) [12] as CHT. Assuming two solutions x and y, x is said to constrained dominate y (denoted as x ≺ c y), if one of the following conditions holds: • CV(x) < CV(y); • CV(x) = CV(y) and x ≺ y. where x ≺ y means that x dominates y with both of the following conditions are satisfied: • ∀i ∈ {1, · · · , m}, f i (x) ≤ f i (y); • ∃i ∈ {1, · · · , m}, f i (x) < f i (y). CDP is the most widely used CHT due to its simplicity on logic and implementation. Wang et al. [23] proposed a two phase based framework termed ToP in which a transformed constrained single objective optimization problem is optimized in the first phase and the original CMOP is optimized in the second phase. Yao et al. [20] proposed a two archive framework in which the convergence archive aims to promote convergence and feasibility while the diversity archive aims to promote diversity. Liu et al. [24] proposed a bidirectional coevolution framework in which one population approximates constrained Pareto front (CPF) from the feasible regions and the other from the infeasible regions. Tian et al. [17] proposed a weak coevolutionary framework in which one population solves the original CMOP while the other solves the constraints-ignored helper problem, and the two coevolved populations generate offspring on their own. Latter, Tian et al. [25] proposed an adaptively switched two stage framework in which one stage handles the original CMOP while the other ignores constraints. Wang et al. [18] proposed a two ranking based fitness evaluation mechanism in which CDP and Pareto dominance are both considered.

CMOEAs Based on ε-Constrained
The second category adopts constraint relaxation technique i.e., ε-constrained technique [26] as CHT. With the relaxation ε, some infeasible solutions might survive to the next generation. Therefore, it enhances the ability of the algorithm to jump out of local feasible regions. Specifically, for two solutions x and y, x ε-constrained dominates y (denoted as x ≺ ε y), if one of the following conditions holds: • CV(x) < ε & CV(y) < ε and x ≺ y; • CV(x) = CV(y) = ε and x ≺ y; • CV(x) > ε & CV(y) > ε and CV(x) < CV(y).
ε-constrained technique is also widely used in CMOP community since its relaxation of constraints can alleviate the shortages of CDP. Fan et al. [21] developed an update strategy of ε for the CMOPs with large infeasible regions. Then, they further did some research works on this topic in PPS [22] and PPS-M2M [27] by introducing new update strategies of ε. Zhang et al. [16] developed the detect-and-escape framework in which a detect strategy is designed to tell whether the search is trapped in local optimal and the escape strategy is designed to escape from the local optimal. Zeng et al. [28] proposed to transform a CMOP to a dynamic MOP and use an ε-constrained technique to handle constraints.

CMOEAs Based on Penalty Function
The third category adopts penalty function [29] as CHT. The penalty function combines the CV and objective function values to evaluate the fitness of a solution. Generally, the fitness of solution x, F (x), based on penalty function is calculated as where β is the penalty factor. Jiao et al. [30] developed a feasible-guiding strategy to preserve some valuable infeasible solutions. In c-DPEA [31], Ming et al.designed a self-adaptive penalty function to preserve competitive infeasible solutions in one of its co-evolutionary populations. In ShiP [32],

Evaluating the Performance of a CMOEA by HV
Hypervolume (HV) [36] is a performance indicator used for evaluating the performance of an algorithm [37]. HV calculates the volume or hypervolume of the area surrounded by the obtained solution set and the reference point. Due to the independence of true CPF, HV is widely used in the CMOEA community for the performance evaluation when comparing a proposed approach to existing one, especially for real-world application problems whose true CPF is mostly unknown. Therefore, through the comparison of HV value, we can easily tell which CMOEA is better and which is worse. So in this paper, we propose to use HV as the selection criteria in designing the ensemble framework.

Ensemble of CMOEAs
As mentioned and analyzed above, one single CMOEA is not likely to be omnipotent for all kinds of CMOPs. Given the fact that there have been already many CMOEAs that can suit CMOPs with different features, developing an ensemble framework which can adaptively select the most effective CMOEA for a given CMOP is valuable. To be more specific, if an ensemble framework can integrate different CMOEAS such as CCMO and PPS, and tell which one can solve the current CMOP well, then, the ensemble obtains better versatility. In the next section, our proposed ensemble framework, as well as the generated ECMOEA, will be detailed.

The Proposed Ensemble Framework
The proposed ensemble framework of CMOEAs can be expressed as the flowchart shown in Figure 2. After initialization, the three embedded CMOEAs (denoted as CMOEA1, CMOEA2 and CMOEA3) are evolved simultaneously. Before g < 0.5 * G max (i.e., the first half of evolutionary process), all the three CMOEAs are evolved. When g = 0.5 * G max , the worst CMOEA is detected by comparing the HV values of the three CMOEAs' populations. Then, the CMOEA with worst performance (i.e., HV value) is deleted. The evolution of the rest two CMOEAs (assume they are CMOEA1 and CMOEA2) are continued until g < 0.8 * G max . When g = 0.8 * G max , the worse CMOEA is detected through the comparison between their populations' HV values. Then, the one with smaller HV value is deleted. Afterwards, the last survived CMOEA is evolved to the end of the evolutionary process. Finally, the population of the last CMOEA is output as the final solution set.
Generally, the search process is divided into three stages. In the first stage (g < 0.5 * G max ), all the CMOEAs embedded are evolved. In the second stage (0.5 * G max ≤ g < 0.8 * G max ), one CMOEA is deleted and the remained two are continued. In the third stage (0.8 * G max ≤ g < G)max), the last CMOEA is evolved to the end.

The Proposed ECMOEA
Based on the proposed ensemble framework as related in the above part, a new CMOEA termed ECMOEA is generated. In this paper, the generated ECMOEA is embedded with ICMA, CCMO and PPS. Specifically, these three state-of-the-art CMOEAs work as the CMOEA1, CMOEA2 and CMOEA3 in the framework. The pseudo code is shown in Algorithm 1. The details are as follows. Similarly, these three CMOEAs are denoted as CMOEA1, CMOEA2 and CMOEA3. Besides, the deleted CMEOAs are assumed to be CMOEA3 and CMOEA2 without losing generality.

Algorithm 1 The Framework of ECMOEA
Input: N, termination condition (g max ) Output: P 1: P 1 , P 2 , P 3 Initialize the populations and parameters for all CMOEAs embedded; 2: Initialize the determined mark t 1 = t 2 = 0, initialize the generation index g = 0; 3: while g ≤ g max do 4: if g < 0.5 × g max then 5: else if g < 0.8 × g max then 9: if t 1 == 0 then 10: Calculate the HV values of all populations P 1 , P 2 , P 3 ;

11:
Delete the CMOEA with minimum HV value (assume it is CMOEA3); 12: Set t 1 = 1; First of all, the populations and parameters (include the reference vector set if needed) are initialized (line 1). Furthermore, the determined mark t 1 , t 2 and the generation index g are initialized (line 2). Then, until the termination condition g ≥ G max met, the following steps are performed.

•
If 0 ≥ g < 0.5 * G max , lines 5-7 are performed, which is the first stage as mentioned above. • If 0.5 * G max ≥ g < 0.8 * G max , lines 9-15 are performed, which is the second stage as mentioned above. The determined mark t 1 is used in this stage to determine whether the worst CMOEA is deleted. • If 0.8 * G max ≥ g < G max , lines 18-26 are performed, which is the third stage as mentioned above. The determined mark t 2 is used in this stage to determine whether the worse CMOEA is deleted.
With ICMA, CCMO and PPS embedded, an instantiation is presented in Algorithm 2 to better understand the principle and procedure of the proposed framework, as well as the generated ECMOEA. In this example, without losing generality, we assume that ICMA is the most suitable CMOEA for the CMOP, while PPS is the most unsuitable CMOEA. Then the algorithm goes as follows: First, three populations and all the parameters needed are initialized (lines 1-2). Then, in the first stage (lines 4-7), all the three CMOEAs are co-evolved to evolve P 1 , P 2 and P 3 , respectively. When g = 0.5 × g max , the HV values of P 1 , P 2 and P 3 are calculated. Since PPS is the most unsuitable CMOEA, is obtained the minimum HV value. Therefore, at the first decision time, PPS is deleted (lines 9-12).Then, since t 1 is set 1, during 0.5 × g max < g < 0.8 × g max , ICMA and CCMO are co-evolved (lines [14][15]. When g = 0.8 × g max , it comes to the second decision time (lines [18][19][20][21]. Since ICMA is the most suitable CMOEA for the CMOP to be solved, the HV value of P 2 , obtained by CCMO, is less than the HV value of P 1 . Therefore, at this time, CCMO is deleted and in the rest evolutionary process, only ICMA is evolved (lines 22-23). Finally, P 1 , the population obtained by the most suitable CMOEA (i.e., ICMA) is regarded as the final solution set.

Computational Complexity
As mentioned above in this Section, ECMOEA mainly includes the embedded CMEOAs and the HV calculation. The computational complexity of the embedded CMOEAs are the same to the original CMOEAs. Since the calculation of HV values only performed at g = 0.5 * G max and g = 0.8 * G max , the time consumption can be ignored. Therefore, the computational complexity of the proposed ensemble framework is determined by the embedded CMOEAs. In this paper, since ICMA, CCMO and PPS are embedded, whose computational complexities are O(mN 2 ), O(mN 3 ) and O(mN 2 ), the generated ECMOEA's computational complexity is O(mN 3 ).

Experimental Results and Analysis
In this section, experimental studies are conducted on PlatEMO [38] matlab platform in order to evaluate the effectiveness of the proposed ensemble framework and the generated ECMOEA. Section 4.1 details the experimental settings, then, Section 4.2 demonstrates the effectiveness of the proposed ensemble framework. Section 4.3 presents the comparison results of ECMOEA to other methods and the analysis. Next, Section 4.4 evaluates the effectiveness of ECMOEA on real-world CMOPs. Finally, the effectiveness of ECMOEA on constrained many-objective optimization problems (CMaOPs) are tested in Section 4.5.

Performance Indicators
The performance of a CMOEA is evaluated by two performance indicators: inverted generational distance (IGD) [10] and hypervolume (HV) [36]. IGD is calculated by Equation (5) when S is the final solution set, W is the reference vector set sampled on the true PF and d(w i , x j ) is the Euclidean distance in the objective space between the solution x j in S and the vector w i in W: Suppose λ is the Lebesgue measure, v i is the hypervolume shaped by the solution x i and the reference point, then HV is calculated by Equation (6): A smaller IGD value means better performance of the algorithm. In contrast, the greater the HV value, the better the performance of the algorithm.
According to the approaches in [41], 10,000 uniformly distributed points were sampled on the true CPF for IGD calculation. As 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 for the HV calculation. Each algorithm was run independently 30 times on each test case. The mean and standard deviations of each metric (i.e., HV and IGD) were recorded. Then the Wilcoxon rank sum test [42] with a significance level of 0.05 and Friedman test [43] with Bonferroni correction at a significance level 0.05 were adopted to perform statistical analysis on the experimental results. In all the tables showing statistical results, we use "+", "−", and "=" to show that the result obtained by another CMOEA is significantly better, significantly worse or statistically similar to that obtained by ECMOEA.

Genetic Operators and Parameter Settings
PPS and MOEA/D-DAE adopt the DE operator since they are embedded with MOEA/D [13]. ToP also adopts DE as genetic operator. For these three methods, the parameter settings are α = 0.95, τ = 0.1, cp = 2, and l = 20. For ECMOEA and the other CMOEAs, Polynomial Mutation (PM) [12] and Simulated Binary Crossover (SBX) [44] are applied in the GA evolutionary operator with parameters set as follows.
All other parameters for the MOEAs are the same as in their original literatures, which are the default settings in PlatEMO [38]. These parameter settings did not change in our experiments.
The population size N is set to 91 for all problems in the experiments. The maximum function evaluation times is set to 100,000 for C-DTLZ, DC-DTLZ, FCP and MW and 150,000 for DAS-CMOP and LIR-CMOP to ensure the convergence of these CMOP instances.

Demonstration of the Effectiveness of the Proposed Ensemble Framework
First of all, we designed the experiment to demonstrate the effectiveness of the proposed ensemble framework. Specifically, three state-of-the-art CMOEAS, ICMA, CCMO and PPS are embedded into the proposed framework, and a new CMOEA termed ECMOEA is generated as related above. Then, ECMOEA is compared with ICMA, CCMO and PPS in this part. Thus, the effectiveness of the proposed framework can be verified through comparison to its original implants.
The statistical results of ECMOEA and CCMO, ICMA, PPS on DAS-CMOP, DTLZ, FCP, LIR-CMOP and MW in terms of HV and IGD are reported in Tables A1-A10. Assuming >, < and = mean the former CMOEA performed better than, worse than and equal to the latter one. From these tables we can see that: From these results, we can see that the most suitable CMOEA for DAS-CMOP, DTLZ, FCP, LIR-CMOP and MW are ICMA, CCMO, ICMA, ICMA (and PPS) and ICMA, respectively. Under the proposed framework, the generated ECMOEA finally obtained better results than those CMOEAs except for the most suitable one. This is because the evaluations are aparted to those worse CMOEAs, thus, ECMOEA can not concentrate on evolving the most suitable CMOEA. However, the results indicate that the proposed framework can finally select the suitable CMOEA.

Comparison to Other CMOEAs
Then, in this part, ECMOEA is compared with seven state-of-the-art CMOEAs: CMOEA-MS, C-TAEA, DCNSGA-III, MOEA/D-DAE, NSGA-II-ToR, TiGE-2 and ToP. The statistical results of ECMOEA and CMOEAs in comparison on DAS-CMOP, DTLZ, FCP, LIR-CMOP and MW in terms of HV and IGD are reported in Tables A11-A20, respectively. From the tables, the proposed ECMOEA obtained the best overall performance on DAS-CMOP, FCP, LIR-CMOP and MW benchmarks. As for DTLZ, ECMOEA was slightly worse than MOEA/D-DAE, but was better than other methods. Especially, ECMOEA performed well on FCP and LIE-CMOP benchmarks.
To better illustrate the results, the feasible and non-dominated solution sets obtained by ECMOEA on all DAS-CMOP, DTLZ, FCP, LIR-CMOP and MW instances with the median IGD value among 30 runs are depicted in Figures A1-A5. From these figures we can see ECMOEA can finally approximate CPF and obtain good distribution in face of the different features of these benchmarks. When dealing with DAS-CMOP with adjustable feasibility, convergence and diversity, the ensemble framework can select ICMA and CCMO to adapt to the requirements. When dealing with DTLZ with regular and irregular CPF, discontinuous feasible regions and inner-outer layer feasible regions, the weak coevoluion based CCMO can be selected to adapt to these features. When dealing with FCP, ICMA can be selected whose indicator-based CHT can preserve some infeasible solutions to better deal with FCP instances. Similarly, when dealing with LIR-CMOP, ICMA and PPS, which is tailored for LIR-CMOP, can be selected. When dealing with MW, which contains different geometric relations between PF and CPF, CCMO can be selected to adapt to the CMOPs whose CPF is close to PF, and ICMA or PPS can be selected when CPF is far from PF.
The KEEL [45] software was used to perform the non-parameter statistical analysis on HV and IGD results. The Friedman test and the Wilcoxon test were conducted. The results of Friedman test are reported in Table 1, where the grey-blodface indicates the best ranking obtained by an algorithm and the blodface means the second ranking of an algorithm. As we can see, ECMOEA obtained the first ranking, while ICMA and CCMO obtained the second and the third, respectively. The results of Wilcoxon test are reported in Tables 2 and 3, respectively. As can be seen, ECMOEA was significantly better than other all methods. Therefore, ECMOEA outperformed other methods statistically. Upper diagonal of level significance α = 0.9,Lower diagonal level of significance α = 0.95. (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11)

ECMOEA on Real-World CMOPs
In this part, experiments of ECMOEA and CMOEAs in comparison on eight mechanical and chemical applications from RWMOP [46] are performed. The detailed parameter settings and information about these real-world CMOPs are shown in Table 4. The original references can be found in RWMOP literature. Since the true CPF of these real-world CMOPs are unknown, only HV is used for the comparison of ECMOEA and other methods. The statistical results of HV obtained by ECMOEA and CMOEAs in comparison on these eight instances are reported in Table A21. It is clear that ECMOEA obtained the best overall results on all the eight instances. ECMOEA outperformed most other methods when dealing with various real-world problems with various features. This is because ECMOEA can adapt to different problems with the help of the ensemble framework.
Two bar truss design (TBT) problem is chosen for further studies. The mathematical formulations of objective functions and constraints are as follows: Minimize : x 3 x 1 subject to : with bounds : The feasible and non-dominated solution sets obtained by ECMOEA and other methods on TBT with the median IGD value among 30 runs are presented in Figure A6. From the results we can find that ECMOEA finally approximated CPF and obtained good distribution.

ECMOEA on CMaOPs
In this part, the effectiveness of ECMOEA on CMaOPs is tested through comparisons on C-DTLZ and DC-DTLZ test suites. For m = 5, 8, n = 9, 12 for C1-DTLZ1 and DC1-DTLZ1, n = 14, 17 for other instances, the population size is N = 126, 156 and the maximum evaluation is E max = 100,000, 140,000, respectively. The experimental settings are the same to which introduced in Section 4.1.
The statistical results on HV and IGD are listed in Tables A22 and A23, respectively. Generally, ECMOEA obtained slightly better overall performance than ICMA, CCMO and PPS. The reasons can be that, since three CMOEAs (i.e., ICMA, CCMO and PPS) are all coevolved during the evolutionary process, the function evaluations are divided. Therefore, in dealing with some CMOPs that the embedded CMOEA is good at (e.g., ICMA and CCMO on C-DTLZ), even ECMOEA finally find the most suitable CMOEA, less evaluations are used. On the contrary, ECMOEA can outperform ICMA and CCMO on some DC-DTLZ instances, because ECMOEA finally found that PPS is more suitable for these DC-DTLZ instances as DC-DTLZ instances are of more complex infeasible regions.
However, the other algorithms, especially the constrained many-objective optimization evolutionary algorithms (CMaOEAs) in comparison, obtained better performance than ECMOEA.

Conclusions
In this paper, an ensemble framework of CMOEA is proposed, into which any existing CMOEA or new CMOEA can be easily embedded. The framework possesses the ability to tell which CMOEA is the best for the given CMOP, and thus, choosing the most suitable CMOEA as the solver. Based on this ensemble framework and three state-of-the-art CMOEAs, a new CMOEA termed ECMOEA is proposed. The proposed ECMOEA can gradually delete the poorly performed CMOEAs and finally choose the most suitable one.
Experimental results on five benchmarks of totally 52 instances have verified the effectiveness of the proposed framework, and demonstrated the superiority or at least competitiveness of ECMOEA compared to seven state-of-the-art CMOEAs. Furthermore, the effectiveness of ECMOEA on real-world applications have been testified through experiments on eight real-world application CMOPs.
In the future, the following attempts are worth trying: • Embedding other CMOEAs or mathematical methods into the framework to further enhance the versatility of ECMOEA; • Designing an adaptive strategy instead of the current 0.8E max , 0.5E max strategy can save the evaluations for the most suitable CMOEA; • Adopting advanced techniques in Machine learning field such as the k-nearest neighbors algorithm (KNN) [47] to improve the efficiency of the ensemble framework is also worth trying; • Embedding effective CMaOEAs in the ensemble framework to enhance the performance on handling CMaOPs is also expected.

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

Abbreviations
The following abbreviations are used in this manuscript: