I Introduction
Many realworld optimization problems can be summarized as optimizing a number of conflicting objectives simultaneously with a set of equality and/or inequality constraints, which are denoted as constrained multiobjective optimization problems (CMOPs). Without lose of generality, a CMOP considered in this paper can be defined as follows [1]:
(1) 
where is an
dimensional objective vector, and
. is an inequality constraint, and is the number of inequality constraints. is an equality constraint, and represents the number of equality constraints. is an dimensional decision vector.When solving CMOPs with both inequality and equality constraints, we usually convert the equality constraints into inequality ones by introducing an extremely small positive number . The detailed transformation is given as follows:
(2) 
To deal with a set of constraints in CMOPs, the overall constraint violation is a widely used approach, which summarizes them into a single scalar, as follows:
(3) 
Given a solution , if , is a feasible solution. All the feasible solutions constitute a feasible region , which is defined as . For any two solutions , is said to dominate if for each and , denoted as . If there is no other solution in dominating solution , then is called a Pareto optimal solution. All of the Pareto optimal solutions constitute a Pareto optimal set (). The mapping of the in the objective space is called a Pareto optimal front (), which is defined as .
A key issue in CMOEAs is to maintain a balance the objectives and constraints. In fact, most constrainthandling mechanisms in evolutionary computation are applied for this purpose.
For example, the penalty function approach adopts a penalty factor to maintain the balance between the objectives and constraints. It converts a CMOP into an unconstrained MOP by adding the overall constraint violation multiplied by a predefined penalty factor to each objective [2]. In the case of , it is called a death penalty approach [3], which means that infeasible solutions are totally unacceptable. If is a static value during the search process, it is called a static penalty approach [4]. If is changing during the search process, it is called a dynamic penalty approach [5]. In the case of is changing according to the information collected from the search process, it is called an adaptive penalty approach [6, 7, 8, 9].
In oder to avoid tuning the penalty factors, another type of constrainthandling methods is suggested, which compares the objectives and constraints separately. Representative examples include constraint dominance principle (CDP) [10], epsilon constrainthandling method (EC) [11], stochastic ranking approach (SR) [12], and so on.
In CDP [10], three basic rules are adopted to compare any two solutions. In the first rule, given two solutions , if is feasible and is infeasible, is better than . In the case of and are both infeasible, the one with a smaller constraint violation is better. In the last rule, and are both feasible, and the one dominating the other is better. CDP is a popular constrainthandling method, as it is simple and has no extra parameters. However, it is not suitable for solving CMOPs with very small and narrow feasible regions [13]. For many generations, most or even all solutions in the working population are infeasible when solving CMOPs with this property. In addition, the diversity of the working population can hardly be well maintained, because the selection of solutions is only based on the constraint violations according to the second rule of CDP.
In order to solve CMOPs with small and narrow feasible regions, the epsilon constrainthandling (EC) [11] approach is suggested. It is similar to CDP except for the relaxation of the constraints. In EC, the relaxation of the constraints is controlled by the epsilon level , which can help to maintain the diversity of the working population in the case when most solutions are infeasible. To be more specific, if the overall constraint violation of a solution is less than , this solution is deemed as feasible. The epsilon level is a critical parameter in EC. In the case of , EC is the same as CDP. Although EC can be used to solve CMOPs with small feasible regions, controlling the value of properly is nontrivial at all.
Both CDP [10] and EC [11] first compare the constrains, then compare the objectives. SR [12]
is different from CDP and EC in terms of the order of comparison. It adopts a probability parameter
to decide if the comparison is based on objectives or constraints. For any two solutions, if a random number is less than , the one with the nondominated objectives is better. The comparison is based on objectives. On the contrary, if the random number is greater than , the comparison is first based on the constraints, then on the objectives, which is the same as that of CDP. In the case of , SR is equivalent to CDP.In recent years, many work has been done in the field of manyobjective evolutionary algorithms (MaOEAs) [14], which gives us a new way to solve CMOPs. In order to balance the constraints and the objectives, some researchers adopt multiobjective evolutionary algorithms (MOEAs) or MaOEAs (in the case of the number of objectives is greater than three) to deal with constraints [15]. For an Mobjective CMOP, its constraints can be converted into one or extra objectives. Then the Mobjective CMOP is transformed into an  or objective unconstrained MOP, which can be solved by MOEAs or MaOEAs. Representative examples include Cai and Wang’s Method (CW) [16] and the infeasibility driven evolutionary algorithms (IDEA) [17].
To maintain a good balance between the constraints and the objectives, some researchers combine several constrainthandling mechanisms, which can be further divided into two categories, including adopting different constrainthandling mechanisms in different evolutionary stages or in different subproblems. For example, the adaptive tradeoff model (ATM) [18] uses two different constrainthandling mechanisms, including multiobjective approach and the adaptive penalty functions in different evolutionary stages. The ensemble of constrainthandling methods (ECHM) [19] uses three different constrainthandling techniques, including epsilon constrainthandling (EC) [11], selfadaptive penalty functions (SP) [9] and superiority of feasible solutions (SF) [20]. Three subpopulations are generated in ECHM, and each subpopulation uses a different constrainthandling method.
In this paper, we propose a biphasic CMOEA, namely push and pull search (PPS), to balance the objectives and the constraints. Unlike the abovementioned constrainthandling methods, the PPS divide the search process into two different stages. In the first stage, only the objectives are optimized, which means the working population is pushed toward the unconstrained PF without considering any constraints. Furthermore, the landscape of constraints in CMOPs can be estimated in the push stage, which can be applied to conduct the parameters setting of constrainthandling approaches applied in the pull stage. In the pull stage, an improved EC [11] constrainthandling approach is adopted to pull the working population to the constrained PF. In summary, it provides a new framework and has the following potential advantages.

It has the ability to get across large infeasible regions in front of the constrained PF. Since the constraints are ignored in the push stage, the infeasible regions in front of the real PF present no barriers for the working population.

It facilitates the parameters setting in the constrainthandling methods. Since the landscape of constraints has already been explored by the push process, a lots of information can be discovered and gathered to guide the search and parameter setting in the pull stage.
The rest of the paper is organized as follows. Section II introduces the general idea of PPS. Section III gives an instantiation of PPS in the framework of MOEA/D, called PPSMOEA/D. Section IV designs a set of experiments to compare the proposed PPSMOEA/D with five other CMOEAs, including MOEA/DIEpsilon, MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D. Finally, conclusions are made in section V.
Ii The General Framework of Push and Pull Search
Constraints define the infeasible regions in the decision space, which may also affect the PFs of CMOPs in the objective space. The influence of the infeasible regions on the PFs can be generally classified into four different situations. For each situation, the search behavior of PPS is illustrated by Fig.
14 respectively, which can be summarized as follows.
Infeasible regions block the way towards the PF, as illustrated by Fig. 1(a). In this circumstance, the unconstrained PF is the same as the constrained PF, and PPS has significant advantages compared with other CMOEAs. Since the constraints are ignored in the push stage of PPS, the infeasible regions have no effect on the searching of PPS. Fig. 1(a)(e) show the push process in different stages, where the working population crosses the infeasible regions in this case without any extra efforts. As the constrained PF is the same as the unconstrained PF, the true PF has already been approximated by the working population in the push process, the pull search has no affect on the working population, as shown in Fig. 1(f).

The unconstrained PF is covered by infeasible regions and becomes no more feasible. Every constrained Pareto optimal point lies on some constraint boundaries, as illustrated by Fig. 2(a). In this circumstance, PPS first approaches the unconstrained PF by using the push strategy as illustrated by Fig. 2(a)(c). After the working population approaches the unconstrained PF, the pull strategy is applied to pull the working population towards the real PF, as illustrated by Fig. 2(d)(f).

Infeasible regions make the original unconstrained PF partially feasible, as illustrated by Fig. 3(a). In this situation, some parts of the true PF have already been achieved during the push search as illustrated by Fig. 3(c). In the pull stage, infeasible solutions are pulled to the feasible and nondominated regions as illustrated by Fig. 3(d)(f). Finally, the whole real PF has been found by PPS.

Infeasible regions may reduce the dimensionality of the PF of a CMOP, as illustrated by Fig. 4(a). In the push stage, a set of nondominated solutions without considering constraints are achieved as illustrated by Fig. 4(a)(c) . Then, infeasible solutions in the working population are pulled toward the feasible and nondominated regions in the pull search process, as shown in Fig. 4(d)(f).
(a) (b) (c) 
(d) (e) (f) 
(a) (b) (c) 
(d) (e) (f) 
(a) (b) (c) 
(d) (e) (f) 
(a) (b) (c) 
(d) (e) (f) 
From the above analysis, it can be observed that PPS can deal with CMOPs with all situations of PFs affected by constraints.
The main steps of PPS includes the push and pull search processes. However, when to apply the push or pull search is very critical. A strategy to switch the search behavior is suggested as follows.
(4) 
where represents the max rate of change between the ideal and nadir points during the last th generation. is a userdefined parameter, here we set . The rates of change of the ideal and nadir points during the last th generation are defined in Eq. 5 and Eq. 6 respectively.
(5) 
(6) 
where are the ideal and nadir points in the th generation, and . are the ideal and nadir points in the th generation.
At the beginning of the search, is initialized to . At each generation, is updated according to Eq. 4. If is less or equal then the predefined threshold , the search behavior is switched to the pull search.
To summarize, PPS divides its search process into two different stages, including the push search and pull search. At the first stage, the push search without considering constraints is adopted to approximate the unconstrained PF. Once Eq. 4 is satisfied, the pull search is used to pull the infeasible solutions to the feasible and nondominated regions. In this situation, constraints are considered in the pull search process. PPS terminates when a predefined halting condition is met. In the following section, we will describe the instantiation of the push and pull strategy in a MOEA/D framework.
Iii An Instantiation of PPS in MOEA/D
In this section, the details of an instantiation of the push search method, as well as the method adopted in the pull search, and the proposed PPS in the framework of MOEA/D are described.
Iiia The push search
In the push search stage, an unconstrained MOEA/D is used to search for nondominated solutions without considering any constraints. When solving a CMOP by using MOEA/D, we can decompose the CMOP into a set of single constrained optimization subproblems, and optimize them simultaneously in a collaborative way. Each subproblem is associated with a decomposition function by using a weight vector . In the decompositionbased selection approach, an individual is selected into next generation based on the value of the decomposition function.
There are three popular decomposition approaches, including weighted sum [21], Tchebycheff [21] and boundary intersection approaches [22]. In this paper, we adopt the Tchebycheff decomposition method, with the detailed definition given as follows.
(7) 
where is a weight vector, and . is the ideal points as described in section II.
In the push search stage, a newly generated solution is kept to the next generation based on the value of as described in Algorithm 1.
IiiB The pull search
In this process, some infeasible solutions are pulled to the feasible and nondominated regions. To achieve this, a constrainthandling mechanism is adopted to punish the infeasible solutions in the pull search stage. An improved epsilon constrainthandling to deal with constraints is proposed with the detailed formulation given as follows.
(8) 
where is the ratio of feasible solutions in the th generation. is to control the speed of reducing the relaxation of constraints in the case of , and . is to control the searching preference between the feasible and infeasible regions, and . is to control the speed of reducing relaxation of constraints in the case of . is updated until the generation counter reaches the control generation . is set to the maximum overall constraint violation of the working population at the end of the push search.
The above proposed epsilon setting method in Eq. 8 is different from our previous work [23] and [24]. In this paper, is set to the maximum overall constraint violation of the working population at the end of the push search, while is set to the overall constraint violation of the top th individual in the initial population in [23] and [24]. In the case of , according to the epsilon setting in [23] and [24], which can not help to explore the infeasible regions and get across large infeasible regions. To avoid this situation, is increased if is greater than a predefined threshold in [23] and [24]. However, increasing properly is still an issue that needs to be further investigated. In this paper, does not need to be increased to get across large infeasible regions, because the working population has already gotten across large infeasible regions in the push stage.
In the pull stage, a newly generated solution is selected into the next generation based on the value of and the overall constraint violation as illustrated by Algorithm 2.
IiiC PPS Embedded in MOEA/D
Algorithm 3 outlines the psuecode of PPSMOEA/D, a PPS embedded in MOEA/D. A CMOP is decomposed into single objective subproblems, and these subproblems are initialized at line 1. The max rate of change of ideal and nadir points is initialized to , and the flag of search stage . Then, the algorithm repeats lines 327 until generations have been reached. The value of and the search strategy are set at lines 414 based on Eq. 4 and Eq. 8. and are initialized at line 7, and the function MaxViolation() calculates the maximum overall constraint violation in the current working population. Lines 1519 show the process of generating a new solution. The updating of the ideal points are illustrated by lines 2022. Lines 2333 show the updating of subproblems, and different search strategies are adopted according to the state of the searching. In the case of , the push search is adopted, otherwise, the pull search is used. The generation count is updated at line 35. The set of feasible and nondominated solutions is updated according to the nondominated ranking in NSGAII [10] at line 36.
Remark: The proposed PPS is a general framework for solving CMOPs. It can be instantiated in many different MOEAs, even though only PPSMOEA/D is realized in this paper. At each search stage, many information can be gathered to extract useful knowledge which can be used to conduct the searching in both search stages. In fact, knowledge discovery is a critical step in the PPS framework. In this paper, we only adopts some statistic information. For example, the maximum overall constraint violation at the end of the push search stage is adopted to set the value of . The ratio of feasible solutions at the pull search stage is used to control the value of
. In fact, many data mining methods and machine learning approaches can be integrated in the PPS framework for solving CMOPs more efficiently.
Iv Experimental study
Iva Experimental Settings
To evaluate the performance of the proposed PPS method, the other five CMOEAs, including MOEA/DIEpsilon, MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D, are tested on LIRCMOP114 [24]. The detailed parameters are listed as follows:

The mutation probability ( denotes the dimension of a decision vector). The distribution index in the polynomial mutation is set to 20.

DE parameters: , .

Population size: . Neighborhood size: .

Halting condition: each algorithm runs for 30 times independently, and stops until 300,000 function evaluations are reached.

Probability of selecting individuals from its neighborhood: .

The max number of solutions updated by a child: .

Parameters setting in PPSMOEA/D: , , , , .

Parameters setting in MOEA/DIEpsilon: , , , and .

Parameters setting in MOEA/DEpsilon: , , and .

Parameter setting in MOEA/DSR: .
IvB Performance Metric
To measure the performance of PPSMOEA/D, MOEA/DIEpsilon, MOEA/DEpsilon, MOEA/DCDP, MOEA/DSR and CMOEA/D, two popular metrics–the inverted generation distance (IGD) [25] and the hypervolume [26] are adopted.

Inverted Generational Distance (IGD):
The IGD metric reflects the performance regarding convergence and diversity simultaneously. The detailed definition is given as follows:
(9) 
where denotes a set of representative solutions in the real PF, is an approximate PF achieved by a CMOEA. denotes the number of objectives. For LIRCMOPs with two objectives, 1000 points are sampled uniformly from the true PF to construct . For LIRCMOPs with three objectives, 10000 points are sampled uniformly from the PF to constitute . It is worth noting that a smaller value of IGD may indicate better performance with regards to diversity and/or convergence.

Hypervolume ():
reflects the closeness of the set of nondominated solutions achieved by a CMOEA to the real PF. The larger means that the corresponding nondominated set is closer to the true PF.
(10) 
where is the Lebesgue measure, denotes the number of objectives, is a userdefined reference point in the objective space. For each LIRCMOP, the reference point is placed at 1.2 times the distance to the nadir point of the true PF. It is worth noting that a larger value of may indicate better performance regarding diversity and/or convergence.
IvC Discussion of Experiments
IvC1 A Comparison between PPSMOEA/D and MOEA/DIEpsilon
The difference between PPSMOEA/D and MOEA/DIEpsilon is that MOEA/DIEpsilon does not apply the push search without considering any constraints, and they adopt different epsilon setting methods. To illustrate the effect of the PPS framework, the performance regarding IGD and HV metrics of PPSMOEA/D and MOEA/DIEpsilon on LIRCMOP114 is compared.
The statistical results of the IGD values on LIRCMOP114 achieved by PPSMOEA/D and MOEA/DIEpsilon in 30 independent runs are listed in Table I. For LIRCMOP12, LIRCMOP58, LIRCMOP11 and LIRCMOP14, PPSMOEA/D performs significantly better than MOEA/DIEpsilon. For LIRCMOP9, MOEA/DIEpsilon is significantly better than PPSMOEA/D. For the rest of test instances, MOEA/DIEpsilon and PPSMOEA/D have no significant difference.
The statistical results of the HV values on LIRCMOP114 achieved by PPSMOEA/D and MOEA/DIEpsilon in 30 independent runs are listed in Table II. For LIRCMOP13, LIRCMOP58 and LIRCMOP14, PPSMOEA/D is significantly better than MOEA/DIEpsilon. There is no significantly difference between PPSMOEA/D and MOEA/DIEpsilon on LIRCMOP4 and LIRCMOP1013. For LIRCMOP9, PPS is significantly worse than MOEA/DIEpsilon.
Since the real PF of LIRCMOP9 locates on its unconstrained PF and has several disconnected parts, as illustrated by Fig. 5(a), PPSMOEA/D should performs well on this problems. However, the statistic results of both IGD and HV indicate that PPSMOEA/D is significantly worse than MOEA/DIEpsilon. One possible reason is that PPSMOEA/D has not converged to the whole unconstrained PF in the push search stage.
To verify this hypothesis, we adjust the value of in Eq. (5) and Eq. (6), which is a critical parameter to switch the search behavior. The mean values of IGD and HV values of MOEA/DIEpsilon and PPSMOEA/D with different values of are shown in Fig. 5(b) and Fig. 5(c), respectively. We can observe that, in cases of , PPSMOEA/D is better than MOEA/DIEpsilon.
(a) The landscape of LIRCMOP9 (b) IGD (c) HV 
To summarize, PPSMOEA/D is superior to MOEA/DIEpsilon on eight test problems, is not significantly different from MOEA/DIEpsilon on five test problems, is inferior to MOEA/DIEpsilon on only one test instance. It can be concluded that the PPS framework can help to enhance the performance of a CMOEA in terms of the IGD and HV metrics.
IvC2 A Comparison of PPSMOEA/D and the other four decompositionbased CMOEAs
The statistical results of the IGD values on LIRCMOP114 achieved by PPSMOEA/D and the other four decompositionbased CMOEAs in 30 independent runs are listed in Table I. According to the WilcoxonTest in this table, it is clear that PPSMOEA/D is significantly better than MOEA/DEpsilon, MOEA/DCDP, MOEA/DSR, CMOEA/D on all of the fourteen tested problems in terms of the IGD metric.
The statistical results of the HV values on LIRCMOP114 achieved by PPSMOEA/D and the other four decompositionbased CMOEAs in 30 independent runs are listed in Table II. It can be observed that PPSMOEA/D is significantly better than MOEA/DEpsilon, MOEA/DCDP, MOEA/DSR and CMOEA/D on all the test instances in terms of the HV metric. It is clear that the proposed PPSMOEA/D performs better than the other four decompositionbased CMOEAs on all of the test instances.
IvC3 The advantages of PPSMOEA/D
The box plots of IGD values for the six tested CMOEAs on LIRCMOP6, LIRCMOP7 and LIRCMOP11 are shown in Fig. 6. It is clear that the mean IGD value of PPSMOEA/D is the smallest among the six tested CMOEAs. To further illustrate the advantages of the proposed PPSMOEA/D, the populations achieved by each tested CMOEAs on LIRCMOP6, LIRCMOP7, and LIRCMOP11 during the 30 independent runs are plotted in Fig. 8, Fig. 9 and Fig. 10, respectively.
For LIRCMOP6, there are two large infeasible regions in front of the PF, and the unconstrained PF is the same as the constrained PF, as illustrated by Fig. 7(a). In Fig. 8, we can observe that PPSMOEA/D can get across the large infeasible regions at each running time, while the other five CMOEAs are sometimes trapped into local optimal.
For LIRCMOP7, there are three large infeasible regions, and the unconstrained PF is covered by infeasible regions and becomes no more feasible as illustrated by Fig. 7(b). In Fig. 9, we can observe that PPSMOEA/D can converge to the real PF at each running time, while the other five CMOEAs can not always converge to the real PF.
For LIRCMOP11, infeasible regions make the original unconstrained PF partially feasible, as illustrated by Fig. 7(c). The PF of LIRCMOP11 is disconnected and has seven Pareto optimal solutions. PPSMOEA/D can find the seven discrete Pareto optimal solutions at each running time, while the other five CMOEAs can not find all the seven Pareto optimal solutions sometimes, as shown in Fig. 10. In Fig. 10(b), it seems that MOEA/DIEpsilon has found all the Pareto optimal solutions. However, MOEA/DIEpsilon can not find all the Pareto optimal solutions more than fifteen out of 30 runs. For MOEA/DEpsilon, MOEA/DCDP, MOEA/DSR and CMOEA/D, some Pareto optimal solutions are not found during the thirty independent runs, as illustrated by 10(c)(f).
Test Instances  PPSMOEA/D  MOEA/DIEpsilon  MOEA/DEpsilon  MOEA/DCDP  MOEA/DSR  CMOEA/D  
LIRCMOP1  mean  6.413E03  7.972E03  5.745E02  1.107E01  1.811E02  1.256E01 
std  1.938E03  3.551E03  2.887E02  5.039E02  1.665E02  7.028E02  
LIRCMOP2  mean  4.673E03  5.234E03  5.395E02  1.434E01  9.634E03  1.397E01 
std  7.844E04  1.009E03  2.128E02  6.055E02  7.232E03  5.436E02  
LIRCMOP3  mean  8.6.0E03  1.131E02  8.806E02  2.608E01  1.775E01  2.805E01 
std  5.184E03  6.422E03  4.356E02  4.335E02  7.196E02  4.214E02  
LIRCMOP4  mean  4.677E03  4.848E03  6.506E02  2.527E01  1.949E01  2.590E01 
std  1.116E03  2.051E03  3.011E02  4.340E02  6.401E02  3.513E02  
LIRCMOP5  mean  1.837E03  2.132E03  1.150E+00  1.045E+00  1.040E+00  1.103E+00 
std  9.263E05  3.791E04  1.979E01  3.632E01  3.655E01  2.994E01  
LIRCMOP6  mean  2.490E03  2.334E01  1.269E+00  1.090E+00  9.427E01  1.307E+00 
std  3.399E04  5.062E01  2.947E01  5.198E01  5.899E01  2.080E01  
LIRCMOP7  mean  2.797E03  3.735E02  1.509E+00  1.460E+00  1.079E+00  1.564E+00 
std  9.854E05  5.414E02  5.094E01  6.082E01  7.581E01  4.235E01  
LIRCMOP8  mean  2.778E03  2.749E02  1.620E+00  1.375E+00  1.012E+00  1.579E+00 
std  7.558E05  5.917E02  3.054E01  6.149E01  7.244E01  3.706E01  
LIRCMOP9  mean  9.940E02  4.977E03  4.902E01  4.810E01  4.854E01  4.810E01 
std  1.519E01  1.367E02  4.221E02  5.244E02  4.775E02  5.243E02  
LIRCMOP10  mean  2.108E03  2.111E03  2.132E01  2.159E01  1.925E01  2.130E01 
std  7.754E05  7.111E05  5.315E02  6.810E02  6.812E02  4.634E02  
LIRCMOP11  mean  2.832E03  5.809E02  3.468E01  3.418E01  3.157E01  3.806E01 
std  1.359E03  5.793E02  9.285E02  9.216E02  7.491E02  8.949E02  
LIRCMOP12  mean  2.704E02  3.358E02  2.520E01  2.690E01  2.064E01  2.6.0E01 
std  5.002E02  5.178E02  8.980E02  9.058E02  5.613E02  9.628E02  
LIRCMOP13  mean  6.455E02  6.460E02  1.200E+00  1.210E+00  8.864E01  1.180E+00 
std  2.177E03  1.639E03  3.055E01  3.172E01  5.756E01  3.775E01  
LIRCMOP14  mean  6.419E02  6.540E02  1.021E+00  1.107E+00  1.027E+00  1.250E+00 
std  1.690E03  2.037E03  4.855E01  3.976E01  4.700E01  5.298E02  
WilcoxonTest (SDI)  –  851  1400  1400  1400  1400 
Wilcoxon’s rank sum test at a 0.05 significance level is performed between PPSMOEA/D and each of the other five CMOEAs. and denote that the performance of the corresponding algorithm is significantly worse than or better than that of PPSMOEA/D, respectively. Where ’SDI’ indicates PPSMOEA/D is superior to, not significantly different from or inferior to the corresponding compared CMOEAs.
Test Instances  PPSMOEA/D  MOEA/DIEpsilon  MOEA/DEpsilon  MOEA/DCDP  MOEA/DSR  CMOEA/D  
LIRCMOP1  mean  1.016E+00  1.014E+00  9.6.0E01  7.538E01  9.956E01  7.414E01 
std  1.580E03  2.433E03  3.280E02  8.947E02  2.915E02  1.222E01  
LIRCMOP2  mean  1.349E+00  1.348E+00  1.283E+00  1.061E+00  1.337E+00  1.070E+00 
std  1.010E03  1.318E03  2.880E02  1.078E01  1.466E02  9.099E02  
LIRCMOP3  mean  8.703E01  8.682E01  7.976E01  4.856E01  5.914E01  4.707E01 
std  2.650E03  3.919E03  3.929E02  4.305E02  1.072E01  4.089E02  
LIRCMOP4  mean  1.093E+00  1.093E+00  1.021E+00  7.349E01  8.147E01  7.308E01 
std  2.467E03  2.458E03  4.191E02  5.436E02  8.699E02  5.162E02  
LIRCMOP5  mean  1.462E+00  1.461E+00  4.297E02  1.630E01  1.817E01  9.721E02 
std  2.919E04  1.332E03  2.353E01  4.426E01  4.389E01  3.699E01  
LIRCMOP6  mean  1.129E+00  9.259E01  5.398E02  1.880E01  3.021E01  2.328E02 
std  1.771E04  4.235E01  2.214E01  3.873E01  4.616E01  1.275E01  
LIRCMOP7  mean  3.015E+00  2.863E+00  3.031E01  3.743E01  9.878E01  2.035E01 
std  2.663E03  1.958E01  9.070E01  9.577E01  1.271E+00  7.522E01  
LIRCMOP8  mean  3.017E+00  2.936E+00  1.060E01  5.173E01  1.099E+00  1.658E01 
std  1.139E03  1.855E01  5.488E01  1.054E+00  1.198E+00  6.114E01  
LIRCMOP9  mean  3.570E+00  3.707E+00  2.737E+00  2.770E+00  2.752E+00  2.770E+00 
std  2.242E01  1.878E02  1.483E01  1.843E01  1.640E01  1.841E01  
LIRCMOP10  mean  3.241E+00  3.241E+00  2.885E+00  2.879E+00  2.927E+00  2.888E+00 
std  3.077E04  2.482E04  1.023E01  1.364E01  1.347E01  9.765E02  
LIRCMOP11  mean  4.390E+00  4.232E+00  3.341E+00  3.353E+00  3.383E+00  3.243E+00 
std  2.217E04  1.840E01  2.573E01  2.568E01  2.899E01  2.551E01  
LIRCMOP12  mean  5.614E+00  6.093E+00  4.884E+00  4.827E+00  5.032E+00  4.890E+00 
std  1.6.0E01  1.578E01  3.168E01  3.284E01  1.755E01  3.453E01  
LIRCMOP13  mean  5.710E+00  5.711E+00  4.6.0E01  4.628E01  1.892E+00  6.293E01 
std  1.275E02  1.301E02  1.301E+00  1.423E+00  2.572E+00  1.712E+00  
LIRCMOP14  mean  6.193E+00  6.182E+00  1.334E+00  8.813E01  1.270E+00  1.795E01 
std  1.310E02  1.093E02  2.452E+00  1.969E+00  2.288E+00  2.598E01  
WilcoxonTest (SDI)  –  851  1400  1400  1400  1400 
Wilcoxon’s rank sum test at a 0.05 significance level is performed between PPSMOEA/D and each of the other five CMOEAs. and denote that the performance of the corresponding algorithm is significantly worse than or better than that of PPSMOEA/D, respectively. Where ’SDI’ indicates PPSMOEA/D is superior to, not significantly different from or inferior to the corresponding compared CMOEAs.
(a) LIRCMOP6 (b) LIRCMOP7 (c) LIRCMOP11 
(a) LIRCMOP6 (b) LIRCMOP7 (c) LIRCMOP11 
(a) PPSMOEA/D (b) MOEA/DIEpsilon (c) MOEA/DEpsilon 
(d) MOEA/DCDP (e) MOEA/DSR (f) CMOEA/D 
(a) PPSMOEA/D (b) MOEA/DIEpsilon (c) MOEA/DEpsilon 
(d) MOEA/DCDP (e) MOEA/DSR (f) CMOEA/D 
(a) PPSMOEA/D (b) MOEA/DIEpsilon (c) MOEA/DEpsilon 
(d) MOEA/DCDP (e) MOEA/DSR (f) CMOEA/D 
V Conclusion
This paper proposes a general PPS framework to deal with CMOPs. More specifically, the search process of PPS is divided into two stages, including the push and pull search processes. At the push stage, constraints are ignored, which can help PPS to get across infeasible regions in front of the unconstrained PFs. Moreover, the landscape affected by constraints can be estimated at the push stage, and this information, such as the ration of feasible solutions and the maximum overall constraint violation, can be applied to conduct the settings of parameters coming from the constrainthandling mechanisms in the pull stage. When the max rate of change between ideal and nadir points is less or equal than a predefined threshold, PPS is switched to the pull search process. The infeasible solutions achieved in the push stage are pulled to the feasible and nondominated area by adopting an improved epsilon constrainthandling technique. The value of epsilon level can be set properly according to the maximum overall constraint violation obtained at the end of the push search stage. The comprehensive experiments indicate that the proposed PPS achieves competitive or statistically significant better results than the other five CMOEAs on most of the benchmark problems.
It is also worthwhile to point out that there has been very little work regarding using information of landscape affected by constraints to solve CMOPs. In this context, the proposed PPS provides a viable framework. Obviously, a lot of work need to be done to improve the performance of PPS, such as, the augmented constrainthandling mechanisms in the pull stage, the enhanced strategies to switch the search behavior and the data mining methods and machine learning approaches integrated in the PPS framework. For another future work, the proposed PPS will be implemented in the nondominated framework, such as NSGAII, to further verify the effect of PPS. More other CMOPs and realworld constrained engineering optimization problems will also be used to test the performance of the PPS embedded in different MOEA frameworks.
Acknowledgment
This research work was supported by Guangdong Key Laboratory of Digital Signal and Image Processing, the National Natural Science Foundation of China under Grant (61175073, 61300159, 61332002, 51375287), Jiangsu Natural Science Foundation (BK20130808) and Science and Technology Planning Project of Guangdong Province, China (2013B011304002).
References
 [1] K. Deb, Multiobjective optimization using evolutionary algorithms. John Wiley & Sons, 2001, vol. 16.
 [2] C. A. C. Coello, “Theoretical and numerical constrainthandling techniques used with evolutionary algorithms: a survey of the state of the art,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 11–12, pp. 1245–1287, 2002.

[3]
T. BDack, F. Hoffmeister, and H. Schwefel, “A survey of evolution
strategies,” in
Proceedings of the 4th international conference on genetic algorithms
, 1991, pp. 2–9.  [4] A. Homaifar, C. X. Qi, and S. H. Lai, “Constrained optimization via genetic algorithms,” Simulation, vol. 62, no. 4, pp. 242–253, 1994.
 [5] J. A. Joines and C. R. Houck, “On the use of nonstationary penalty functions to solve nonlinear constrained optimization problems with ga’s,” in Evolutionary Computation, 1994. IEEE World Congress on Computational Intelligence., Proceedings of the First IEEE Conference on. IEEE, 1994, pp. 579–584.
 [6] J. C. Bean and A. ben HadjAlouane, A dual genetic algorithm for bounded integer programs, 1993.
 [7] D. W. Coit, A. E. Smith, and D. M. Tate, “Adaptive penalty methods for genetic optimization of constrained combinatorial problems,” INFORMS Journal on Computing, vol. 8, no. 2, pp. 173–182, 1996.
 [8] A. Ben HadjAlouane and J. C. Bean, “A genetic algorithm for the multiplechoice integer program,” Operations research, vol. 45, no. 1, pp. 92–101, 1997.
 [9] Y. G. Woldesenbet, G. G. Yen, and B. G. Tessema, “Constraint handling in multiobjective evolutionary optimization,” IEEE Transactions on Evolutionary Computation, vol. 13, no. 3, pp. 514–525, June 2009.
 [10] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGAII,” IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182–197, Apr 2002.

[11]
T. Takahama and S. Sakai, “Constrained optimization by
constrained particle swarm optimizer with
level control,” in Soft Computing as Transdisciplinary Science and Technology. Springer, 2005, pp. 1019–1029.  [12] T. P. Runarsson and X. Yao, “Stochastic ranking for constrained evolutionary optimization,” IEEE Transactions on evolutionary computation, vol. 4, no. 3, pp. 284–294, 2000.
 [13] Z. Fan, W. Li, X. Cai, K. Hu, H. Lin, and H. Li, “Anglebased constrained dominance principle in moea/d for constrained multiobjective optimization problems,” in IEEE Congress on Evolutionary Computation, 2016, pp. 460–467.
 [14] B. Li, J. Li, K. Tang, and X. Yao, “Manyobjective evolutionary algorithms: A survey,” ACM Computing Surveys (CSUR), vol. 48, no. 1, p. 13, 2015.
 [15] E. MezuraMontes and C. A. Coello Coello, “Constrainthandling in natureinspired numerical optimization: Past, present and future,” Swarm and Evolutionary Computation, vol. 1, no. 4, pp. 173–194, Dec. 2011.
 [16] Z. Cai and Y. Wang, “A multiobjective optimizationbased evolutionary algorithm for constrained optimization,” Evolutionary Computation, IEEE Transactions on, vol. 10, no. 6, pp. 658–675, 2006.
 [17] T. Ray, H. K. Singh, A. Isaacs, and W. Smith, “Infeasibility driven evolutionary algorithm for constrained optimization,” in Constrainthandling in evolutionary optimization. Springer, 2009, pp. 145–165.
 [18] Y. Wang, Z. Cai, Y. Zhou, and W. Zeng, “An adaptive tradeoff model for constrained evolutionary optimization,” IEEE Transactions on Evolutionary Computation, vol. 12, no. 1, pp. 80–92, 2008.
 [19] B. Y. Qu and P. N. Suganthan, “Constrained multiobjective optimization algorithm with an ensemble of constraint handling methods,” Engineering Optimization, vol. 43, no. 4, pp. 403–416, 2011.
 [20] K. Deb, “An efficient constraint handling method for genetic algorithms,” Computer methods in applied mechanics and engineering, vol. 186, no. 2, pp. 311–338, 2000.
 [21] K. Miettinen, Nonlinear Multiobjective Optimization. Springer Science & Business Media, 1999, vol. 12.
 [22] Q. Zhang and H. Li, “MOEA/D: A multiobjective evolutionary algorithm based on decomposition,” IEEE Transactions on evolutionary computation, 2007.
 [23] Z. Fan, W. Li, X. Cai, H. Li, H. Huang, and Z. Cai, “An improved epsilon constraint handling method embedded in MOEA/D for constrained multiobjective optimization problems,” in 2016 IEEE Symposium Series on Computational Intelligence, Dec 2016.
 [24] Z. Fan, W. Li, X. Cai, H. Huang, Y. Fang, Y. You, J. Mo, C. Wei, and E. D. Goodman, “An improved epsilon constrainthandling method in MOEA/D for cmops with large infeasible regions,” arXiv preprint arXiv:1707.08767, 2017.
 [25] P. A. Bosman and D. Thierens, “The balance between proximity and diversity in multiobjective evolutionary algorithms,” Evolutionary Computation, IEEE Transactions on, vol. 7, no. 2, pp. 174–188, 2003.
 [26] E. Zitzler and L. Thiele, “Multiobjective evolutionary algorithms: a comparative case study and the strength pareto approach,” IEEE Transactions on Evolutionary Computation, vol. 3, no. 4, pp. 257–271, 1999.
Comments
There are no comments yet.