An Improved Epsilon Constrainthandling Method in MOEA/D for CMOPs with Large Infeasible Regions
Abstract
This paper proposes an improved epsilon constrainthandling mechanism, and combines it with a decompositionbased multiobjective evolutionary algorithm (MOEA/D) to solve constrained multiobjective optimization problems (CMOPs). The proposed constrained multiobjective evolutionary algorithm (CMOEA) is named MOEA/DIEpsilon. It adjusts the epsilon level dynamically according to the ratio of feasible to total solutions (RFS) in the current population. In order to evaluate the performance of MOEA/DIEpsilon, a new set of CMOPs with two and three objectives is designed, having large infeasible regions (relative to the feasible regions), and they are called LIRCMOPs. Then the fourteen benchmarks, including LIRCMOP114, are used to test MOEA/DIEpsilon and four other decompositionbased CMOEAs, including MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D. The experimental results indicate that MOEA/DIEpsilon is significantly better than the other four CMOEAs on all of the test instances, which shows that MOEA/DIEpsilon is more suitable for solving CMOPs with large infeasible regions. Furthermore, a realworld problem, namely the robot gripper optimization problem, is used to test the five CMOEAs. The experimental results demonstrate that MOEA/DIEpsilon also outperforms the other four CMOEAs on this problem.
Keywords:
Constrained Multiobjective Evolutionary Algorithms Epsilon Constrainthandling Constrained Multiobjective Optimization Robot Gripper Optimization∎
1 Introduction
Realworld optimization problems usually involve the simultaneous optimization of multiple conflicting objectives with a number of constraints. Without loss of generality, a CMOP considered in this paper is defined as follows (Deb (2001)):
(1) 
where is an dimensional objective vector, is an inequality constraint, and is an equality constraint. is an dimensional decision vector. The feasible region is defined as the set and .
In CMOPs, there are usually more than one constraint. The overall constraint violation is a widely used approach to deal with constraint violations, as it summarizes them into a single scalar, as follows:
(2) 
If , is feasible; otherwise, it is infeasible. Any solution in set is feasible, and for any two solutions and , is said to dominate if for each and for at least one , denoted as . For a solution , if there is no other solution in dominating , then is called a Pareto optimal solution. A set including all of the Pareto optimal solutions is called a Pareto optimal set (). Mapping the into the objective space obtains a set of objective vectors, which is called a Pareto optimal front (), and .
CMOEAs aim to find a representative set of Pareto optimal solutions. They have to tackle the multiple conflicting objectives with a number of constraints simultaneously, and to maintain a good balance between convergence and diversity of the achieved solutions. In CMOEAs, there are two basic components: one is the constrainthandling mechanism, and the other is the multiobjective evolutionary algorithm (MOEA).
In terms of constrainthandling, many methods have been proposed in evolutionary optimization (Cai et al (2013); Hu et al (2013)). They can be roughly divided into penalty function methods, special representations and operators, repair methods, separation of objectives and constraints and hybrid methods (Coello (2002)). The penalty function method is widely used due to its simplicity in the constraint handling (Runarsson and Yao (2005)). However, the ideal penalty factors cannot be known in advance for an arbitrary CMOP, and tuning the penalty factors can be a very tedious task.
In recent years, a number of other constrainthandling techniques have had a relatively high impact in evolutionary optimization, including feasibility rules, stochastic ranking, constrained method, novel penalty functions, novel special operators, multiobjective concepts and ensemble of constrainthandling techniques (MezuraMontes and Coello Coello (2011)). However, most of them aim to solve constrained scalar optimization problems when they are first proposed.
MOEAs can be classified into three different types according to their selection approaches. The first type is nondominatedbased methods, and representative examples include NSGAII (Deb et al (2002)), PAESII (Corne et al (2001)), SPEAII (Zitzler et al (2001)), NSGAIII (Deb and Jain (2014)) and so on. The second type is decompositionbased approaches, and typical examples include MOEA/D (Zhang and Li (2007)), MOEA/DDE (Li and Zhang (2009)), EAGMOEA/D (Cai et al (2015)), MOEA/DM2M (Liu et al (2014)), MOEA/DSAS (Cai et al (2016)) and so on. Currently, MOEA/D is a popular algorithm to solve unconstrained multiobjective optimization problems (MOPs). MOEA/D decomposes a MOP into many scalar optimization subproblems, and optimizes them simultaneously in a collaborative way. The last type is indicatorbased methods. This type of MOEAs selects solutions based on the improvement of a performance metric. Representative methods include IBEA (Zitzler and Künzli (2004)), SMSEMOA (Beume et al (2007)), HypE (Bader and Zitzler (2011)), FVMOEA (Jiang et al (2015)) and so on.
There are two commonly used test suites of CMOPs, including CTP (Deb (2001)) and CF test instances (Zhang et al (2008)). For CTP1CTP5 and CF1CF10, the feasible regions are relatively large, and a CMOEA can approximate their PFs without encountering any infeasible obstacles during the entire evolutionary process. Thus, CTP15 and CF110 are not good test problems to evaluate the performance of constrainthandling mechanisms. For the remaining test problems CTP68, the feasible regions are relatively large, and the population of a CMOEA can reach these regions with high probability. Thus, CTP and CF test suites can not effectively measure the performance of constrainthandling techniques. When solving CTP (Deb (2001)) and CF (Zhang et al (2008)) test instances, the constraint dominance principle (CDP) (Deb et al (2002)) is good enough to handle the constraints.
To overcome the shortcomings of the CTP and CF test suites discussed above, we propose a set of new CMOPs (named LIRCMOP114). Each of them has a number of large infeasible regions, and the feasible regions are relatively small. The population of a CMOEA cannot easily discover these small feasible regions, which brings new challenges to the existing CMOEAs. In fact, many realworld optimization problems also have this characteristic. For example, the robot gripper optimization problem considered in this paper has large infeasible regions as illustrated in Section 6. Thus, it has important significance in practice to design specific mechanisms for solving CMOPs with large infeasible regions.
In this paper, we propose an improved constrained version of MOEA/D to deal with CMOPs. Compared with the original constrained method (Takahama and Sakai (2006)), the proposed method can keep a good balance in the search between the feasible and infeasible regions. It uses the information of the feasible ratio of the population to dynamically balance the exploration between the feasible regions and infeasible regions.
The remainder of the paper is organized as follows. Section 2 introduces related work on MOEA/D and the existing CMOEAs based on MOEA/D. Section 3 illustrates the improved epsilon constrainthandling method as here embedded in MOEA/D. Section 4 designs a set of new CMOPs (LIRCMOPs) with large infeasible regions. Section 5 describes a comprehensive set of experiments to compare the proposed CMOEA (MOEA/DIEpsilon) with four other CMOEAs, including MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D. In Section 6, a robot gripper optimization problem is used to test MOEA/DIEpsilon and the other four CMOEAs. Finally, Section 7 presents the conclusions.
2 Related work
2.1 Moea/d
MOEA/D (Zhang and Li (2007)) decomposes a MOP into a number of scalar optimization subproblems and optimizes them simultaneously in a collaborative way. Each subproblem is defined by a decomposition function with a weight vector . In MOEA/D, a set of uniformly spread weight vectors are adopted to formulate subproblems. The weight vectors satisfy and for each . In terms of decomposition methods, there are three commonly used approaches, including weighted sum (Miettinen (1999)), Tchebycheff (Miettinen (1999)) and boundary intersection approaches (Zhang and Li (2007)).
In the weighted sum approach, each subproblem is defined by summing each objective weighted by a different weight. The th subproblem with the weighted sum decomposition method is defined as follows:
minimize  (3)  
subject to 
For a minimizing MOP, in the case of a convex PF, the weighted sum approach can work well. However, if the PF is nonconvex, only a part of PF can be found by this approach.
In the Tchebycheff decomposition method, the th subproblem is defined as follows:
minimize  (4)  
subject to 
where is the ideal point, and . The Tchebycheff method is a widely used decomposition approach. It can approximate both concave and convex parts of PFs.
In the boundary intersection approach, two distances and are defined to evaluate the convergence and diversity respectively. The th subproblem is defined as follows:
minimize  
subject to  
where  
The boundary intersection method is able to solve MOPs with any shape of PFs. However, the penalty factor must be set in advance.
2.2 Decompositionbased CMOEAs
In decompositionbased CMOEAs, a CMOP is decomposed into a set of constrained scalar optimization subproblems, and these subproblems are solved in a collaborative way simultaneously. Representative methods include CMOEA/D (Asafuddoula et al (2012)), MOEA/DEpsilon (Yang et al (2014)), MOEA/DCDP (Jan and Khanum (2013)) and MOEA/DSR (Jan and Khanum (2013)).
CMOEA/D (Asafuddoula et al (2012)) embeds an epsilon constrainthandling approach into MOEA/D, and the epsilon value is set adaptively. To be more specific, the epsilon level is set to . denotes the mean value of the overall constraint violation in the current population, and () denotes the feasible ratio of solutions in the current population. For two solutions, if their overall constraint violations are both less than or their overall constraint violations are equal, the one with the better aggregation value is selected. Otherwise, the one with the smaller overall constraint violation is selected.
MOEA/DEpsilon (Yang et al (2014)) also adopts the epsilon method to handle constraints. Unlike CMOEA/D, the epsilon value in MOEA/DEpsilon is set dynamically with the increase of generation counter . The detailed setting of the epsilon value can be found in (Takahama and Sakai (2006)).
MOEA/DCDP (Jan and Khanum (2013)) adopts CDP (Deb et al (2002)) to deal with constraints in the framework of MOEA/D. There are three basic rules to select solutions. For two feasible solutions, the one with the better aggregation value is selected. For two infeasible solutions, the one with the smaller overall constraint violation is selected. For a feasible and an infeasible solution, the feasible one is selected.
MOEA/DSR (Jan and Khanum (2013)) embeds the stochastic ranking method (SR) (Runarsson and Yao (2000)) in MOEA/D to deal with constraints. A parameter is set to balance the selection between the objectives and the constraints in MOEA/DSR. For two solutions, if a random number is less than , the one with the better aggregation value is selected into the next generation. If the random number is greater than , the solutions selection is similar to that of MOEA/DCDP. In the case of , MOEA/DSR is equivalent to MOEA/DCDP.
In summary, CMOEA/D and MOEA/DEpsilon both adopt the epsilon constrainthandling approach to solve CMOPs. To get across large infeasible regions, should be increased at sometimes, and be greater than the maximum overall constraint violation in the current population. However, in CMOEA/D, is always less or equal than , and in MOEA/DEpsilon, is always decreasing during the evolutionary process. In MOEA/DCDP, feasible solutions are always better than infeasible solutions. Thus, the infeasible solutions which can help to get across large infeasible regions are difficult to survive. MOEA/DSR applies a parameter to balance the searching between the feasible and infeasible regions. In order to get across large infeasible regions, should be set dynamically. However, is a static parameter in MOEA/DSR. To overcome the shortcomings of the four decompositionbased CMOEAs discussed above, an improved epsilon constrainthandling method embedded in MOEA/D is proposed.
3 The Proposed method
In this section, the concept of epsilon level comparison, the original epsilon level setting method and the improved epsilon level setting approach are described.
3.1 Epsilon Level Comparison
In the epsilon constraint handling approach (Takahama and Sakai (2006)), the relaxation of constraints is controlled by the epsilon level . For two solutions and , their overall constraint violations are and . Then, for any satisfying , the epsilon level comparison is defined as follows:
3.2 Epsilon Level Setting
In the epsilon constrainthandling method, the setting of is quite critical. In (Takahama and Sakai (2006)), an epsilon level setting method is suggested as follows:
(7) 
where is the top th individual of the initial population sorted by overall constraint violations in a descending order. is to control the speed of reducing relaxation of constraints. is updated until the generation counter reaches the control generation . When , . The recommended parameter ranges in (Takahama and Sakai (2006)) are listed as follows: , and . denotes the population size, and represents the maximum evolutionary generation.
3.3 Improved Epsilon Level Setting
The setting of in Eq.(7) is always decreasing during the evolutionary process, which may not be suitable to solve CMOPs with large infeasible regions. To overcome this problem, an improved epsilon setting approach is suggested as follows:
(8) 
where is the overall constraint violation of the top th individual in the initial population, is the ratio of feasible solutions in the th generation. ranges between and , and has two functions. One is to control the speed of reducing the relaxation of constraints, and the other is to control the scale factor multiplied by the maximum overall constraint violation. is to control the searching preference between the feasible and infeasible regions, and . is the maximum overall constraint violation found so far.
The setting method in Eq. 8 is sometimes the same as that in Eq. 7. If , in Eq. 7 is identically equal to zero, which tends to hinder a CMOEA’s exploration of the infeasible regions. However, in Eq. 8 is not identically equal to zero when according to the third rule of the proposed epsilon setting approach.
In the case , three rules are adopted to control the value of in Eq. 8. is adopted to strengthen the searching in the feasible regions. is used to strengthen the exploration in the infeasible regions. The last is same as in the CDP (Deb et al (2002)) constrainthandling method.
Two parameters and are applied to choose the right control rule for . If and , for setting is adopted. In this circumstance, is set to , which has an exponential decreasing rate. It has a faster descent rate than the epsilon setting in Eq. (7), which can help to enhance the searching in the feasible regions more effectively. If and , for setting is applied. In this situation, most solutions are feasible. Thus, strengthening the exploration in the infeasible regions may help a CMOEA to get across a number of large infeasible regions. In , , which strengthens the exploration in the infeasible regions. Thus, the improved epsilon method has the balanced ability to explore the feasible and infeasible regions simultaneously.
is a critical parameter to balance the searching between the feasible and infeasible regions. If the RFS is less than , is adopted to enhance the exploration in the feasible regions. Otherwise, is applied to enhance the exploration in the infeasible regions. Thus, the proposed epsilon constraint method can keep a good balance of exploration between the feasible and infeasible regions. It utilizes the RFS to dynamically balance the exploration between the feasible regions and infeasible regions.
Compared with the setting in Eq. (7), the proposed method in Eq. (8) has the ability to increase during the evolutionary process, which can help to solve CMOPs with large infeasible regions.
In the case of , is applied. In this situation, , and the epsilon constrainthandling method exerts the highest selection pressure toward the feasible regions.
3.4 Embedding the improved epsilon method in MOEA/D
The proposed MOEA/DIEpsilon integrates the improved epsilon constrainthandling method in Eq. 8 into the framework of MOEA/D. In MOEA/DIEpsilon, a CMOP is decomposed into a number of constrained scalar subproblems, and these subproblems are optimized simultaneously in a collaborative way. In our experimental studies, the Tchebycheff approach is adopted, and its detailed definition is listed in Eq. (4).
For a given weight vector , there exists an optimal solution of Eq. (4), and this optimal solution is also a Pareto optimal solution of Eq. (1). Therefore, we can achieve different Pareto optimal solutions of Eq. (1) by setting different weight vectors.
The psuecode of MOEA/DIEpsilon is listed in Algorithm 1. It is almost the same as that of MOEA/D, except for the method of subproblem updating. Lines 16 initialize a number of parameters in MOEA/DIEpsilon. First, a CMOP is decomposed into subproblems which are associated with . Then the population , the initial epsilon value , the ideal point and the neighbor indexes are initialized.
Lines 1122 generate a set of new solutions and update the ideal point . To be more specific, a set of solutions which may be updated by a newly generated solution is selected (lines 1117). In line 18, the differential evolution (DE) crossover is adopted to generate a new solution . The polynomial mutation operator is executed to mutate in line 19. The ideal point is updated (lines 2022).
Lines 2330 implement the updating process of subproblems. In line 26, the subproblems are updated based on the improved epsilon constrainthandling approach, and the detailed procedures are listed in Algorithm 2. Finally, a set of nondominated solutions () is selected based on the nondominated sort in line 33.
In Algorithm 2, there are three basic rules to update a subproblem. For two solutions and , if their overall constraint violations are less than or equal to , and has a smaller aggregation value (the value of the decomposition function) than that of , then is replaced by (lines 37). If and have the same overall constraint violation, and has a smaller aggregation value than that of , then is replaced by (lines 812). Otherwise, if has a smaller overall constraint violation than that of , then is replaced by (lines 1314). When the subproblem is updated, the function returns , otherwise, it returns .
4 Test instances
To evaluate the performance of the proposed MOEA/DIEpsilon, a set of new CMOPs with large infeasible regions (named LIRCMOPs) is designed according to our previous work (Fan et al (2016)). In terms of constraint functions, all of them have large infeasible regions. In term of objective functions, there are two components: shape functions and distance functions (Huband et al (2006)).
The shape functions are applied to set the shape of the PFs. In the LIRCMOP test suite, two types of shape functions, including both convex and concave shapes, are designed. Distance functions are adopted that test the convergence performance of a CMOEA. In LIRCMOP514, the distance functions are multiplied by a scale factor, which increases difficulty of convergence. The detailed definitions of LIRCMOPs are listed in the Appendix.
In this test suite, four test problems, including LIRCMOP14, have large infeasible regions. Fig. 1(a)(d) plot the feasible regions of LIRCMOP14, respectively. It can be seen that the feasible regions of these test instances are very small. In other words, there are a number of large infeasible regions.
LIRCMOP5 and LIRCMOP6 have convex and concave PFs, respectively, as shown in Fig. 1(e)(f) , and their PFs are the same as those of their unconstrained counterparts. The PFs of LIRCMOP5 and LIRCMOP6 can be achieved by a MOEA without any constrainthandling mechanisms.
In order to expand the test scope, LIRCMOP7 and LIRCMOP8 are designed. For these two test instance, their unconstrained PFs are located in the infeasible regions, and their PFs are situated on their constraint boundaries. Thus, a MOEA without constrainthandling methods cannot find the real PFs for LIRCMOP7 and LIRCMOP8, which are shown in Fig. 1(g)(h).
LIRCMOP912 have two different types of constraints. The first type creates large infeasible regions as shown in the black ellipses in Fig. 1(i)(l). The second type creates difficulty in the entire objective space, as it divides the PFs of LIRCMOP912 into a number of disconnected segments. For LIRCMOP910, their PFs are a part of their unconstrained PFs, and for LIRCMOP1112, their PFs are situated on their constraint boundaries.
In the LIRCMOP test suite, CMOPs with three objectives are also designed. Two CMOPs, including LIRCMOP13 and LIRCMOP14, have three objectives as shown in Fig. 2 (a)(b) . The PF of LIRCMOP13 is the same as that of its unconstrained counterpart. The PF of LIRCMOP14 is located on the boundaries of its constraints.
(a) LIRCMOP1 (b) LIRCMOP2 (c) LIRCMOP3 
(d) LIRCMOP4 (e) LIRCMOP5 (f) LIRCMOP6 
(g) LIRCMOP7 (h) LIRCMOP8 (i) LIRCMOP9 
(j) LIRCMOP10 (k) LIRCMOP11 (l) LIRCMOP12 
(a) LIRCMOP13 (b) LIRCMOP14 
5 Experimental study
5.1 Experimental Settings
To evaluate the performance of the proposed MOEA/DIEpsilon, four other CMOEAs (MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D), with differential evolution (DE) crossover, are tested on LIRCMOP114. The detailed parameters of these five CMOEAs are listed as follows:

Mutation probability ( is the number of decision variables) and its distribution index is set to 20. , .

Population size: . Neighborhood size: .

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

Probability of selecting individuals in the neighborhood: .

The maximal number of solutions replaced by a child: .

Parameter setting in MOEA/DIEpsilon: , , and .

Parameter setting in MOEA/DEpsilon: , , and .

Parameter setting in MOEA/DSR: .
5.2 Performance Metric
To measure the performance of MOEA/DIEpsilon, CMOEA/D, MOEA/DCDP, MOEA/DSR and MOEA/DEpsilon, two commonly used metrics–the inverted generation distance (IGD) (Bosman and Thierens (2003)) and the hypervolume (Zitzler and Thiele (1999)) are adopted. The definition of IGD is shown next.

Inverted Generational Distance (IGD):
The IGD metric reflects the performance regarding convergence and diversity simultaneously. The detailed definition is as follows:
(9) 
where is a set of representative solutions in the ideal PF, is an approximate PF achieved by a CMOEA. The value of IGD denotes the distance between and . For CMOPs with two objectives, 1000 points are sampled uniformly from the true PF to construct . (Note that this measure cannot be used if the true Pareto front is unknown, so it is used primarily for benchmarking purposes.) For CMOPs with three objectives, 10000 points are sampled uniformly from the PF to constitute . It is worth noting that a smaller value of IGD represents better performance with regards to both diversity and convergence.

Hypervolume ():
reflects the closeness of the nondominated set 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, is a reference point in the objective space. For a 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 represents better performance regarding both diversity and convergence.
5.3 Discussion of Experiments
5.3.1 Performance comparison on LIRCMOP test suite
The statistical results of the IGD values on LIRCMOP114 achieved by five CMOEAs in 30 independent runs are listed in Table 1. As discussed in Section 4, LIRCMOP14 have large infeasible regions in the entire search space. For these four test instances, MOEA/DIEpsilon is significantly better than the other four tested CMOEAs in term of the IGD metric. Fig. 3(a)(b) shows the final populations achieved by each CMOEA with the best IGD values during the 30 runs on LIRCMOP1 and LIRCMOP4. It is clear that MOEA/DIEpsilon has the best performance regarding diversity among the five CMOEAs under test.
LIRCMOP512 have large infeasible regions, as discussed in Section 4. It can be observed that MOEA/DIEpsilon is significantly better than the other four tested CMOEAs on NCMOP512. The final populations achieved by each CMOEA on LIRCMOP9 and LIRCMOP11 with the best IGD values are plotted in Fig. 3(c)(d). For LIRCMOP9, MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D only achieve a part of the real PF. However, MOEA/DIEpsilon can obtain the whole real PF. Thus, MOEA/DIEpsilon performs better than the other four CMOEAs in terms of diversity. For LIRCMOP11, the proposed method MOEA/DIEpsilon can achieve the whole PF. However, the other four CMOEAs do not converge to the whole PF. Thus, MOEA/DIEpsilon has better convergence performance than the other four CMOEAs. For threeobjective test instances (LIRCMOP13 and LIRCMOP14), MOEA/DIEpsilon is also significantly better than the other four CMOEAs.
Table 2 shows the results of the HV values of LIRCMOP114 achieved by five CMOEAs in 30 independent runs. It is clear that MOEA/DIEpilon is significantly better than the other four CMOEAs on all of the fourteen test instances in terms of the metric.
5.3.2 Analysis of Experimental Results
From the above performance comparison on the fourteen test instances LIRCMOP114, it is clear that MOEA/DIEpsilon has better diversity and convergence performance than the other four decompositionbased CMOEAs on these fourteen test instances. A common feature of these test instances is that each of them has a number of large infeasible regions, which demonstrates that the proposed epsilon constrainthandling method can deal with the large infeasible regions very well using its automatic adjustment of the epsilon level.
Test Instances  MOEA/DIEpsilon  MOEA/DEpsilon  MOEA/DSR  MOEA/DCDP  CMOEA/D  

LIRCMOP1  mean  7.213E03  7.432E02  1.719E02  1.163E01  1.290E01 
std  2.425E03  3.538E02  1.554E02  7.265E02  8.055E02  
LIRCMOP2  mean  5.461E03  6.407E02  9.274E03  1.244E01  1.627E01 
std  1.520E03  3.869E02  9.723E03  5.492E02  5.819E02  
LIRCMOP3  mean  1.117E02  9.570E02  1.792E01  2.460E01  2.751E01 
std  6.856E03  4.529E02  7.306E02  4.444E02  3.895E02  
LIRCMOP4  mean  4.859E03  6.141E02  2.034E01  2.486E01  2.631E01 
std  1.591E03  4.127E02  6.038E02  3.858E02  3.331E02  
LIRCMOP5  mean  2.107E03  9.455E01  1.041E+00  9.827E01  8.637E01 
std  2.616E04  4.705E01  3.833E01  4.140E01  5.071E01  
LIRCMOP6  mean  2.058E01  1.177E+00  8.699E01  1.224E+00  1.277E+00 
std  4.587E01  4.376E01  5.992E01  3.726E01  2.587E01  
LIRCMOP7  mean  4.598E02  1.475E+00  1.074E+00  1.402E+00  1.511E+00 
std  6.855E02  5.309E01  7.606E01  6.226E01  5.032E01  
LIRCMOP8  mean  3.445E02  1.522E+00  1.253E+00  1.361E+00  1.575E+00 
std  6.002E02  4.716E01  6.597E01  5.888E01  3.849E01  
LIRCMOP9  mean  1.290E02  4.902E01  4.883E01  4.994E01  4.902E01 
std  3.300E02  4.221E02  4.130E02  2.526E02  4.221E02  
LIRCMOP10  mean  2.143E03  2.202E01  1.898E01  2.042E01  2.114E01 
std  1.261E04  3.589E02  6.277E02  6.573E02  5.641E02  
LIRCMOP11  mean  4.713E02  3.809E01  2.911E01  3.221E01  3.321E01 
std  5.410E02  1.131E01  3.525E02  7.723E02  7.109E02  
LIRCMOP12  mean  4.711E02  2.574E01  2.045E01  2.289E01  2.472E01 
std  5.662E02  8.768E02  6.771E02  7.823E02  8.883E02  
LIRCMOP13  mean  6.447E02  1.239E+00  1.059E+00  1.190E+00  1.215E+00 
std  1.844E03  2.555E01  5.033E01  3.290E01  3.140E01  
LIRCMOP14  mean  6.502E02  1.172E+00  9.005E01  1.204E+00  1.054E+00 
std  1.635E03  3.043E01  5.455E01  2.244E01  4.515E01 
Wilcoxonâs rank sum test at a 0.05 significance level is performed between MOEA/DIEpsilon and each of the other four CMOEAs. and denote that the performance of the corresponding algorithm is significantly worse than or better than that of MOEA/DIEpsilon, respectively. The best mean is highlighted in boldface.
Test Instances  MOEA/DIEpsilon  MOEA/DEpsilon  MOEA/DSR  MOEA/DCDP  CMOEA/D  

LIRCMOP1  mean  1.015E+00  9.413E01  9.840E01  7.499E01  7.344E01 
std  1.490E03  3.751E02  4.630E02  1.202E01  1.269E01  
LIRCMOP2  mean  1.348E+00  1.267E+00  1.337E+00  1.093E+00  1.033E+00 
std  1.717E03  5.526E02  2.252E02  1.016E01  9.522E02  
LIRCMOP3  mean  8.686E01  7.964E01  5.892E01  5.034E01  4.715E01 
std  3.373E03  3.618E02  1.105E01  5.141E02  3.786E02  
LIRCMOP4  mean  1.093E+00  1.024E+00  8.048E01  7.397E01  7.203E01 
std  1.910E03  5.903E02  8.956E02  5.264E02  4.480E02  
LIRCMOP5  mean  1.461E+00  2.833E01  1.773E01  2.428E01  3.870E01 
std  9.488E04  5.766E01  4.619E01  5.031E01  6.151E01  
LIRCMOP6  mean  9.412E01  1.255E01  3.341E01  8.582E02  3.750E02 
std  3.848E01  3.325E01  4.458E01  2.707E01  1.446E01  
LIRCMOP7  mean  2.847E+00  3.516E01  9.943E01  4.811E01  2.933E01 
std  2.205E01  9.304E01  1.268E+00  1.083E+00  8.776E01  
LIRCMOP8  mean  2.905E+00  2.690E01  7.043E01  5.223E01  1.788E01 
std  2.103E01  8.100E01  1.094E+00  9.669E01  6.669E01  
LIRCMOP9  mean  3.692E+00  2.737E+00  2.733E+00  2.705E+00  2.737E+00 
std  6.318E02  1.484E01  1.284E01  8.883E02  1.483E01  
LIRCMOP10  mean  3.241E+00  2.874E+00  2.929E+00  2.899E+00  2.886E+00 
std  3.537E04  7.851E02  1.064E01  1.207E01  1.126E01  
LIRCMOP11  mean  4.263E+00  3.218E+00  3.479E+00  3.406E+00  3.386E+00 
std  1.685E01  3.542E01  1.252E01  2.135E01  1.831E01  
LIRCMOP12  mean  5.552E+00  4.858E+00  5.059E+00  4.972E+00  4.902E+00 
std  1.730E01  3.280E01  2.103E01  2.596E01  3.233E01  
LIRCMOP13  mean  5.710E+00  3.097E01  1.184E+00  5.320E01  4.642E01 
std  1.084E02  1.048E+00  2.250E+00  1.442E+00  1.426E+00  
LIRCMOP14  mean  6.184E+00  5.617E01  1.912E+00  4.032E01  1.162E+00 
std  1.053E02  1.540E+00  2.705E+00  1.127E+00  2.281E+00 
Wilcoxonâs rank sum test at a 0.05 significance level is performed between MOEA/DIEpsilon and each of the other four CMOEAs. and denotes that the performance of the corresponding algorithm is significantly worse than or better than that of MOEA/DIEpsilon, respectively. The best mean is highlighted in boldface.
(a) LIRCMOP1 (b) LIRCMOP4 (c) LIRCMOP9 (d) LIRCMOP11 
6 Robot Gripper Optimization
To verify the capability of MOEA/DIEpsilon to solve real world optimization problems, a robot gripper optimization problem with two conflicting objectives and eight constraints is explored.
6.1 Definition of the robot gripper optimization
The robot gripper optimization problem is defined in (Saravanan et al (2009); Datta and Deb (2011)). Five objectives are formulated in these papers. The robot gripper optimization problem considered in this paper has two conflicting objectives and eight constraints. The geometrical structure of the gripper is plotted in Fig. 4.
The robot gripper optimization problem considered in this paper is defined as follows:
(11) 
where has seven decision variables, and each variable is shown in Fig. 4. The range of each decision variable is as follows: , , , , , and . Two rules are applied to fix the value of , and they are defined as follows:
According to the geometric analysis, the gripping force in Fig. 4 can be defined as follows:
(12) 
The displacement of the gripper end is defined as follows:
(13) 
where , , , and denotes a dynamic displacement of the gripper actuator in the range of 0 to 100 mm.
The first objective represents a force transmission ratio between the actuating force and the minimum gripping force . We prefer to transform more actuating force into the gripper force. Thus, this objective should be minimized.
The second objective is the sum of all elements of the robot gripper. It is relevant to the weight of the robot gripper, and minimizing can lead to a lightweight design.
To study the distribution of solutions in the objective space for the robot gripper optimization problem, 3,000,000 solutions are generated, where 1,500,000 solutions are generated randomly, and the other 1,500,000 solutions are generated by MOEA/DIEpsilon. In Fig. 5, we can observe that the robot gripper optimization problem has large infeasible regions (), which can be solved well by the proposed method MOEA/DIEpsilon according to our previous analysis. To verify this hypothesis, MOEA/DIEpislon and the other four decompositionbased CMOEAs are tested on the robot gripper optimization problems.
Symbol  Meaning of parameter  Value 

Minimal dimension of  50mm  
object to be gripped  
Maximal range of the  150mm  
gripper ends displacement  
Maximal dimension of  100mm  
object to be gripped  
Maximal displacement of  100mm  
the gripper actuator  
Actuating force of the  100N  
gripper  
The lower bound of  50N  
gripping force 
6.2 Experimental study
6.2.1 Experimental settings
To solve the robot gripper optimization problem and evaluate the performance of the proposed MOEA/DIEpsilon, five decompositionbased CMOEAs, including MOEA/DIEpsilon, MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D with the differential evolution (DE) crossover, are tested on the robot gripper optimization problem. The detailed parameters of these five CMOEAs are the same as listed in Section 5.1 except for the number of function evaluations. In the case of the robot gripper optimization problem, each CMOEA stops when 600,000 function evaluations are reached. As the ideal PF of the gripper optimization problem is not known in advance, we use only the hypervolume metric (Zitzler and Thiele (1999)) to measure the performance of the five tested CMOEAs. In the robot gripper optimization case, the reference point .
6.2.2 Analysis of experiments
Table 4 shows the statistical results of values of MOEA/DIEpsilon and the other four CMOEAs on the robot gripper optimization problem. It is clear that MOEA/DIEpsilon is significantly better than the other four CMOEAs. To further demonstrate the superiority of the proposed method MOEA/DIEpsilon, the nondominated solutions achieved by each CMOEA during the 30 independent runs are plotted in Fig. 6(a)(e). The box plot of values of the five CMOEAs is shown in Fig. 6(f). From Fig. 6, we see that MOEA/DIEpsion has better performance than the other four CMOEAs.
Test Instances  MOEA/DIEpsilon  MOEA/DEpsilon  MOEA/DSR  MOEA/DCDP  cMOEA/D 

mean  1.897E+03  1.891E+03  1.889E+03  1.869E+03  1.865E+03 
std  3.510E+00  7.151E+00  9.839E+00  8.124E+00  9.048E+00 
Wilcoxonâs rank sum test at a 0.05 significance level is performed between MOEA/DIEpsilon and each of the other four CMOEAs. and denote that the performance of the corresponding algorithm is significantly worse than or better than that of MOEA/DIEpsilon, respectively. The best mean is highlighted in boldface.
(a) MOEA/DIEpsilon (b) MOEA/DEpsilon (c) MOEA/DSR 
(d) MOEA/DCDP (e) CMOEA/D (f) The box plots of each CMOEA 
In order to verify the correctness of the optimization results of the robot gripper optimization problem, three representative individuals (A, B and C) are selected from the nondominated solutions achieved by MOEA/DIEpsilon as shown in Fig. 7. The configurations of the robot gripper mechanism at each point are also plotted in Fig. 7.
To measure the minimum gripping force , a spring with a large stiffness coefficient is set vertically at the end of the robot gripper during the simulation process. The spring force is regarded as the gripping force when the robot gripper is balanced by the spring. The simulation tool is ADAMS 2013, and the stiffness coefficient of the spring is N/m.
Table 5 shows the simulation results of the minimum gripping force with three different configurations of the robot gripper. The relative errors between the theoretical gripping forces and the simulated gripping forces are less than . Thus, we can conclude that the optimization results of the robot gripper optimization problem achieved by MOEA/DIEpsilon are correct.
Sampled point  The theoretical gripping force (N)  The simulated result (N)  Relative error 

A  50.0000  50.0002  0.0004% 
B  142.3168  142.4582  0.0994% 
C  92.5285  92.5877  0.0639% 
7 Conclusion
This paper proposes an improved epsilon constrainthandling method embedded in the framework of MOEA/D. A new CMOEA named MOEA/DIEpsilon has been proposed. The comprehensive experimental results indicate that MOEA/DIEpsilon has the ability to cross the large infeasible regions. Compared with the other four decompositionbased CMOEAs including MOEA/DEpsilon, MOEA/DSR, MOEA/DCDP and CMOEA/D, MOEA/DIEpsilon has following advantages:

The performance of MOEA/DIEpsilon is not sensitive to the initial epsilon value.

MOEA/DIEpsilon has the ability to explore the feasible and infeasible regions simultaneously during the evolutionary process.

MOEA/DIEpsilon utilizes the feasible ratio of the current population to dynamically balance the exploration between the feasible regions and infeasible regions. It keeps a good balance of the searching between infeasible and feasible regions.

MOEA/DIEpsilon is suitable for solving CMOPs with large infeasible regions.
In terms of CMOPs, a new set of CMOPs named LIRCMOP114 was designed and presented in this paper. A common feature of these test instances is that they have large infeasible regions. The experimental results show that MOEA/DIEpsion is significantly better than the other four CMOEAs on this test suite. Thus, we hypothesize that MOEA/DIEpsilon is better than the other four CMOEAs in solving CMOPs with large infeasible regions, in general. To demonstrate the capacity of MOEA/DIEpsilon to solve real engineering problems, a robot gripper optimization problem with two conflicting objectives and eight constraints was used as a test problem. The experimental results also demonstrated that MOEA/DIEpsilon outperformed the other four CMOEAs.
Proposed further work includes studying new constrainthandling mechanisms to solve CMOPs with different types of difficulty. One possible way is to collect more information about the working population, and utilize such information to guide a CMOEA to select appropriate constrainthandling methods in different evolutionary stages.
Acknowledgements.
This work was supported in part by the National Natural Science Foundation of China (NSFC) under grant 61300159, 61473241 and 61332002, by the Natural Science Foundation of Jiangsu Province of China under grant BK20130808, by the Project of Internation as well as Hongkongï¼Macao&Taiwan Science and Technology Cooperation Innovation Platform in Universities in Guangdong Province under grant 2015KGJH2014, by China Postdoctoral Science Foundation under grant 2015M571751, by the Science and Technology Planning Project of Guangdong Province of China under grant 2013B011304002, by Educational Commission of Guangdong Province of China under grant 2015KGJHZ014, by the Fundamental Research Funds for the Central Universities of China under grant NZ2013306, and by the Guangdong HighLevel University Project “Green Technologies” for Marine Industries.8 Appendix
In this section, the detailed definitions of LIRCMOP114 are listed in Table 6.
Problem  Objectives  Constraints 

LIRCMOP1  
LIRd2  
LIRCMOP3  
LIRCMOP4  
LIRCMOP5  
LIRCMOP6  
LIRCMOP7 