An Improved Epsilon Constraint-handling Method in MOEA/D for CMOPs with Large Infeasible Regions
Zhun Fan Wenji Li Xinye Cai Han Huang Yi Fang Yugen You Jiajie Mo Caimin Wei Erik Goodman
cs.NE













An Improved ε-constrained Method in MOEA/D for CMOPs with Large Infeasible Regions        













[1]Department of Electronic Engineering, Shantou University, Guangdong, 515063, China

[2]College of Computer Science and Technology, Nanjing University of Aeronautics and Astronautics, Jiangsu, 210016, China

[3]School of Software Engineering, South China University of Technology, Guangdong, 515063, China

[4]Department of Mathematics, Shantou University, Guangdong, 515063, China

[5]BEACON Center for the Study of Evolution in Action, Michigan State University. East Lansing, Michigan, USA.








An Improved Epsilon Constraint-handling Method in MOEA/D for CMOPs with Large Infeasible Regions
    Zhun Fan  [1]         Wenji Li  [1]         Xinye Cai [2]         Han Huang [3]         Yi Fang   [1]         Yugen You [1]         Jiajie Mo [1]         Caimin Wei [4]        Erik Goodman  [5]

    Received: date / Accepted: date
======================================================================================================================================================================================================




This paper proposes an improved epsilon constraint-handling mechanism, and combines it with a decomposition-based multi-objective evolutionary algorithm (MOEA/D) to solve constrained multi-objective optimization problems (CMOPs). The proposed constrained multi-objective evolutionary algorithm (CMOEA) is named MOEA/D-IEpsilon. 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/D-IEpsilon, 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 LIR-CMOPs. Then the fourteen benchmarks, including LIR-CMOP1-14, are used to test MOEA/D-IEpsilon and four other decomposition-based CMOEAs, including MOEA/D-Epsilon, MOEA/D-SR, MOEA/D-CDP and C-MOEA/D. The experimental results indicate that MOEA/D-IEpsilon is significantly better than the other four CMOEAs on all of the test instances, which shows that MOEA/D-IEpsilon is more suitable for solving CMOPs with large infeasible regions. Furthermore, a real-world problem, namely the robot gripper optimization problem, is used to test the five CMOEAs. The experimental results demonstrate that MOEA/D-IEpsilon also outperforms the other four CMOEAs on this problem.  









§ INTRODUCTION


Real-world 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 (<cit.>): 


    𝐅(𝐱) = (f_1(𝐱),…,f_m(𝐱)) ^ T
        g_i(𝐱) ≥ 0, i = 1,…,q 
        h_j(𝐱) = 0, j= 1,…,p 
       𝐱∈ℝ^n

where F(𝐱) = (f_1(𝐱),f_2(𝐱), … ,f_m(𝐱)) ^ T ∈ℝ ^m is an m-dimensional objective vector, g_i(𝐱) ≥ 0 is an inequality constraint, and h_j(𝐱)=0 is an equality constraint. 𝐱∈ℝ^n is an n-dimensional decision vector. The feasible region S is defined as the set {𝐱 | g_i(𝐱) ≥ 0, i = 1,…,q  and h_j(𝐱) = 0, j= 1,…,p}.

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: 

    ϕ(𝐱) = ∑_i=1^q |min(g_i(𝐱),0)| + ∑_j = 1^p |h_j(𝐱)|

If ϕ(𝐱) = 0, 𝐱 is feasible; otherwise, it is infeasible. Any solution in set S is feasible, and for any two solutions 𝐱^1 ∈ S and 𝐱^2 ∈ S, 𝐱^1 is said to dominate 𝐱^2 if f_i(𝐱^1) ≤ f_i(𝐱^2) for each i ∈{1,...,m} and f_j(𝐱^1) < f_j(𝐱^2) for at least one j ∈{1,...,m}, denoted as 𝐱^1 ≼𝐱^2. For a solution 𝐱^* ∈ S, if there is no other solution in S dominating 𝐱^*, then 𝐱^* is called a Pareto optimal solution. A set including all of the Pareto optimal solutions is called a Pareto optimal set (PS). Mapping the PS into the objective space obtains a set of objective vectors, which is called a Pareto optimal front (PF), and PF = {F(𝐱)| 𝐱∈ PS}. 

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 constraint-handling mechanism, and the other is the multi-objective evolutionary algorithm (MOEA).

In terms of constraint-handling, many methods have been proposed in evolutionary optimization (<cit.>). They can be roughly divided into penalty function methods, special representations and operators, repair methods, separation of objectives and constraints and hybrid methods (<cit.>). The penalty function method is widely used due to its simplicity in the constraint handling (<cit.>). 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 constraint-handling techniques have had a relatively high impact in evolutionary optimization, including feasibility rules, stochastic ranking, ε-constrained method, novel penalty functions, novel special operators, multi-objective concepts and ensemble of constraint-handling techniques (<cit.>). 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 non-dominated-based methods, and representative examples include NSGA-II (<cit.>), PAES-II (<cit.>), SPEA-II (<cit.>), NSGA-III (<cit.>) and so on. The second type is decomposition-based approaches, and typical 
examples include MOEA/D (<cit.>), MOEA/D-DE (<cit.>), EAG-MOEA/D (<cit.>), MOEA/D-M2M (<cit.>), MOEA/D-SAS (<cit.>) and so on. Currently, MOEA/D is a popular algorithm to solve unconstrained multi-objective 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 indicator-based methods. This type of MOEAs selects solutions based on the improvement of a performance metric. Representative methods include IBEA (<cit.>), SMS-EMOA (<cit.>), HypE (<cit.>), FV-MOEA (<cit.>) and so on. 

There are two commonly used test suites of CMOPs, including CTP (<cit.>) and CF test instances (<cit.>). For CTP1-CTP5 and CF1-CF10, the feasible regions are relatively large, and a CMOEA can approximate their PFs without encountering any infeasible obstacles during the entire evolutionary process. Thus, CTP1-5 and CF1-10 are not good test problems to evaluate the performance of constraint-handling mechanisms. For the remaining test problems CTP6-8, 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 constraint-handling techniques. When solving CTP (<cit.>) and CF (<cit.>) test instances, the constraint dominance principle (CDP) (<cit.>) 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 LIR-CMOP1-14). 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 real-world 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 <ref>. 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 (<cit.>), 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 <ref> introduces related work on MOEA/D and the existing CMOEAs based on MOEA/D. Section <ref> illustrates the improved epsilon constraint-handling method as here embedded in MOEA/D. Section <ref> designs a set of new CMOPs (LIR-CMOPs) with large infeasible regions. Section <ref> describes a comprehensive set of experiments to compare the proposed CMOEA (MOEA/D-IEpsilon) with four other CMOEAs, including MOEA/D-Epsilon, MOEA/D-SR, MOEA/D-CDP and C-MOEA/D. In Section <ref>, a robot gripper optimization problem is used to test MOEA/D-IEpsilon and the other four CMOEAs. Finally, Section <ref> presents the conclusions.



§ RELATED WORK




 §.§ MOEA/D



MOEA/D (<cit.>) 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 λ^i. In MOEA/D, a set of N uniformly spread weight vectors λ^1,…,λ^N are adopted to formulate N subproblems. The weight vectors λ^i satisfy ∑_k = 1 ^ mλ_k^i = 1 and λ_k^i ≥ 0 for each k ∈{1,…,m}. In terms of decomposition methods, there are three commonly used approaches, including weighted sum (<cit.>), Tchebycheff (<cit.>) and boundary intersection approaches (<cit.>).

In the weighted sum approach, each subproblem is defined by summing each objective weighted by a different weight. The j-th subproblem with the weighted sum decomposition method is defined as follows:

    g^te(𝐱|λ) = ∑_i = 1 ^mλ_i^jf_i(𝐱)
          𝐱∈S

For a minimizing MOP, in the case of a convex PF, the weighted sum approach can work well. However, if the PF is non-convex, only a part of PF can be found by this approach.

In the Tchebycheff decomposition method, the j-th subproblem is defined as follows: 

    g^te(𝐱|λ,z^*) = max_1 ≤ i ≤ m{λ_i^j |f_i(𝐱) - z_i^*| }
          𝐱∈S

where z^* = (z_1^*,…,z_m^*) is the ideal point, and z_i^* = min{f_i(𝐱 | 𝐱∈ S}. 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 d_1 and d_2 are defined to evaluate the convergence and diversity respectively. The j-th subproblem is defined as follows: 

    g^pbi(𝐱|λ^j,z^*) = d_1 + θ d_2
          𝐱∈S
           d_1 = ‖ (F(𝐱) - z^*)^Tλ^j ‖/‖λ^j ‖
           d_2 = ‖ (F(𝐱) - z^* ) - d_1 λ^j/‖λ^j ‖) ‖

The boundary intersection method is able to solve MOPs with any shape of PFs. However, the penalty factor θ must be set in advance.



 §.§ Decomposition-based CMOEAs



In decomposition-based 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 C-MOEA/D (<cit.>), MOEA/D-Epsilon (<cit.>), MOEA/D-CDP (<cit.>) and MOEA/D-SR (<cit.>).

C-MOEA/D (<cit.>) embeds an epsilon constraint-handling approach into MOEA/D, and the epsilon value is set adaptively. To be more specific, the epsilon level is set to CV_mean * FR. CV_mean denotes the mean value of the overall constraint violation in the current population, and FR (Number of feasible solutions/Population size) denotes the feasible ratio of solutions in the current population. For two solutions, if their overall constraint violations are both less than CV_mean * FR 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/D-Epsilon (<cit.>) also adopts the epsilon method to handle constraints. Unlike C-MOEA/D, the epsilon value in MOEA/D-Epsilon is set dynamically with the increase of generation counter K. The detailed setting of the epsilon value can be found in (<cit.>).

MOEA/D-CDP (<cit.>) adopts CDP (<cit.>) 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/D-SR (<cit.>) embeds the stochastic ranking method (SR) (<cit.>) in MOEA/D to deal with constraints. A parameter p_f ∈ [0,1] is set to balance the selection between the objectives and the constraints in MOEA/D-SR. For two solutions, if a random number is less than p_f, the one with the better aggregation value is selected into the next generation. If the random number is greater than p_f, the solutions selection is similar to that of MOEA/D-CDP. In the case of p_f = 0, MOEA/D-SR is equivalent to MOEA/D-CDP. 

In summary, C-MOEA/D and MOEA/D-Epsilon both adopt the epsilon constraint-handling 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 C-MOEA/D, ε is always less or equal than CV_mean, and in MOEA/D-Epsilon, ε is always decreasing during the evolutionary process. In MOEA/D-CDP, 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/D-SR applies a parameter p_f to balance the searching between the feasible and infeasible regions. In order to get across large infeasible regions, p_f should be set dynamically. However, p_f is a static parameter in MOEA/D-SR. To overcome the shortcomings of the four decomposition-based CMOEAs discussed above, an improved epsilon constraint-handling method embedded in MOEA/D is proposed.



§ 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.



 §.§ Epsilon Level Comparison

In the epsilon constraint handling approach (<cit.>), the relaxation of constraints is controlled by the epsilon level ε. For two solutions 𝐱^1 and 𝐱^2, their overall constraint violations are ϕ^1 and ϕ^2. Then, for any ε satisfying ε≥ 0, the epsilon level comparison ≼_ε is defined as follows:


    (𝐱^1,ϕ^1) ≼_ε  (𝐱^2,ϕ^2)⇔   𝐱^1 ≼𝐱^2, if ϕ^1, ϕ^2 ≤ε
    𝐱^1 ≼𝐱^2, if ϕ^1 = ϕ^2 
    ϕ^1 < ϕ^2, otherwise

In Eq. (<ref>), the epsilon comparison approach is equivalent to CDP (<cit.>) when ε = 0. In the case of ε = ∞, it does not consider any constraints. In other words, the comparison between any two solutions is based on their non-dominated ranks on objectives when ε = ∞.




 §.§ Epsilon Level Setting

In the epsilon constraint-handling method, the setting of ε is quite critical. In (<cit.>), an epsilon level setting method is suggested as follows:


    ε(k) = ε(0)(1 - k/T_c)^cp, 0 < k < T_c, ε(0) = ϕ(𝐱^θ)  
    
    0, k ≥ T_c

where 𝐱^θ is the top θ-th individual of the initial population sorted by overall constraint violations in a descending order. cp is to control the speed of reducing relaxation of constraints. ε(k) is updated until the generation counter k reaches the control generation T_c. When k ≥ T_c, ε(k) = 0. The recommended parameter ranges in (<cit.>) are listed as follows: θ = (0.05 * N), cp ∈ [2,10] and T_c ∈ [0.1T_max, 0.8T_max]. N denotes the population size, and T_max represents the maximum evolutionary generation.



 §.§ Improved Epsilon Level Setting

The setting of ε(k) in Eq.(<ref>) 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:



    ε(k) = 
    rule 1:ϕ(𝐱^θ), if  k = 0
    
    rule 2:(1 - τ)ε(k-1), if  r_k < α and  k < T_c
    
    rule 3:(1 + τ)ϕ_max, if  r_k ≥α and  k < T_c 
    
    rule 4:0, if  k ≥ T_c

where ϕ_k(𝐱^θ) is the overall constraint violation of the top θ-th individual in the initial population, r_k is the ratio of feasible solutions in the k-th generation. τ ranges between 0 and 1, 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 α∈[0,1]. ϕ_max is the maximum overall constraint violation found so far.

The ε(0) setting method in Eq. <ref> is sometimes the same as that in Eq. <ref>. If ε(0) = 0, ε(k) in Eq. <ref> is identically equal to zero, which tends to hinder a CMOEA's exploration of the infeasible regions. However, ε(k) in Eq. <ref> is not identically equal to zero when ε(0) = 0 according to the third rule of the proposed epsilon setting approach.

In the case k > 0, three rules are adopted to control the value of ε in Eq. <ref>. Rule 2 is adopted to strengthen the searching in the feasible regions. Rule 3 is used to strengthen the exploration in the infeasible regions. The last rule 4 is same as in the CDP (<cit.>) constraint-handling method.

Two parameters k and r_k are applied to choose the right control rule for ε(k). If k < T_c and r_k < α, rule 2 for setting ε(k) is adopted. In this circumstance, ε(k) is set to (1 - τ)ε(k-1), which has an exponential decreasing rate. It has a faster descent rate than the epsilon setting in Eq. (<ref>), which can help to enhance the searching in the feasible regions more effectively. If k < T_c and r_k ≥α, rule 3 for setting ε(k) 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 rule 3, ε(k) = (1 + τ)ϕ_max, 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 r_k is less than α, rule 2 is adopted to enhance the exploration in the feasible regions. Otherwise, rule 3 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. (<ref>), the proposed method in Eq. (<ref>) has the ability to increase ε(k) during the evolutionary process, which can help to solve CMOPs with large infeasible regions.   

In the case of k ≥ T_c, rule 4 is applied. In this situation, ε(k) = 0, and the epsilon constraint-handling method exerts the highest selection pressure toward the feasible regions.




 §.§ Embedding the improved epsilon method in MOEA/D


The proposed MOEA/D-IEpsilon integrates the improved epsilon constraint-handling method in Eq. <ref> into the framework of MOEA/D. In MOEA/D-IEpsilon, 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. (<ref>).

For a given weight vector λ, there exists an optimal solution of Eq. (<ref>), and this optimal solution is also a Pareto optimal solution of Eq. (<ref>). Therefore, we can achieve different Pareto optimal solutions of Eq. (<ref>) by setting different weight vectors.




The psuecode of MOEA/D-IEpsilon is listed in Algorithm 1. It is almost the same as that of MOEA/D, except for the method of subproblem updating. Lines 1-6 initialize a number of parameters in MOEA/D-IEpsilon. First, a CMOP is decomposed into N subproblems which are associated with λ^1,…, λ^N. Then the population P, the initial epsilon value ε(0), the ideal point z^* and the neighbor indexes B(i) are initialized.

Lines 11-22 generate a set of new solutions and update the ideal point z^*. To be more specific, a set of solutions which may be updated by a newly generated solution 𝐲^j is selected (lines 11-17). In line 18, the differential evolution (DE) crossover is adopted to generate a new solution 𝐲^j. The polynomial mutation operator is executed to mutate 𝐲^j in line 19. The ideal point z^* is updated (lines 20-22).

Lines 23-30 implement the updating process of subproblems. In line 26, the subproblems are updated based on the improved epsilon constraint-handling approach, and the detailed procedures are listed in Algorithm 2. Finally, a set of non-dominated solutions (NS) is selected based on the non-dominated sort in line 33.




In Algorithm 2, there are three basic rules to update a subproblem. For two solutions 𝐱^j and 𝐲^j, if their overall constraint violations are less than or equal to ε(k), and 𝐲^j has a smaller aggregation value (the value of the decomposition function) than that of 𝐱^j, then 𝐱^j is replaced by 𝐲^j (lines 3-7). If 𝐱^j and 𝐲^j have the same overall constraint violation, and 𝐲^j has a smaller aggregation value than that of 𝐱^j, then 𝐱^j is replaced by 𝐲^j (lines 8-12). Otherwise, if 𝐲^j has a smaller overall constraint violation than that of 𝐱^j, then 𝐱^j is replaced by 𝐲^j (lines 13-14). When the subproblem is updated, the function UpdateSubproblems(𝐱^j,𝐲^j,ε(k)) returns true, otherwise, it returns false. 



§ TEST INSTANCES



To evaluate the performance of the proposed MOEA/D-IEpsilon, a set of new CMOPs with large infeasible regions (named LIR-CMOPs) is designed according to our previous work (<cit.>). 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 (<cit.>).

The shape functions are applied to set the shape of the PFs. In the LIR-CMOP 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 LIR-CMOP5-14, the distance functions are multiplied by a scale factor, which increases difficulty of convergence. The detailed definitions of LIR-CMOPs are listed in the Appendix.  

In this test suite, four test problems, including LIR-CMOP1-4, have large infeasible regions. Fig. <ref>(a)-(d) plot the feasible regions of LIR-CMOP1-4, 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.  

LIR-CMOP5 and LIR-CMOP6 have convex and concave PFs, respectively, as shown in Fig. <ref>(e)-(f) , and their PFs are the same as those of their unconstrained counterparts. The PFs of LIR-CMOP5 and LIR-CMOP6 can be achieved by a MOEA without any constraint-handling mechanisms.

In order to expand the test scope, LIR-CMOP7 and LIR-CMOP8 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 constraint-handling methods cannot find the real PFs for LIR-CMOP7 and LIR-CMOP8, which are shown in Fig. <ref>(g)-(h).

LIR-CMOP9-12 have two different types of constraints. The first type creates large infeasible regions as shown in the black ellipses in Fig. <ref>(i)-(l). The second type creates difficulty in the entire objective space, as it divides the PFs of LIR-CMOP9-12 into a number of disconnected segments. For LIR-CMOP9-10, their PFs are a part of their unconstrained PFs, and for LIR-CMOP11-12, their PFs are situated on their constraint boundaries.

In the LIR-CMOP test suite, CMOPs with three objectives are also designed. Two CMOPs, including LIR-CMOP13 and LIR-CMOP14, have three objectives as shown in Fig. <ref> (a)-(b) . The PF of LIR-CMOP13 is the same as that of its unconstrained counterpart. The PF of LIR-CMOP14 is located on the boundaries of its constraints.









§ EXPERIMENTAL STUDY





 §.§ Experimental Settings


To evaluate the performance of the proposed MOEA/D-IEpsilon, four other CMOEAs (MOEA/D-Epsilon, MOEA/D-SR, MOEA/D-CDP and C-MOEA/D), with differential evolution (DE) crossover, are tested on LIR-CMOP1-14. The detailed parameters of these five CMOEAs are listed as follows:



  * Mutation probability Pm = 1/n (n is the number of decision variables) and its distribution index is set to 20. CR = 1.0, f = 0.5.

  * Population size: N = 300. Neighborhood size: T = 30.

  * 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: δ = 0.9.

  * The maximal number of solutions replaced by a child: nr = 2.

  * Parameter setting in MOEA/D-IEpsilon: T_c = 800, α = 0.95, τ = 0.1 and θ = 0.05 N.

  * Parameter setting in MOEA/D-Epsilon: T_c = 800, cp = 2, and θ = 0.05 N. 

  * Parameter setting in MOEA/D-SR: S_r = 0.05.




 §.§ Performance Metric

To measure the performance of MOEA/D-IEpsilon, C-MOEA/D, MOEA/D-CDP, MOEA/D-SR and MOEA/D-Epsilon, two commonly used metrics–the inverted generation distance (IGD) (<cit.>) and the hypervolume (<cit.>) 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:

    IGD(P^*,A) = ∑_y^* ∈ P^*d(y^*,A)/| P^* |
    
    
    d(y^*,A) = min_y ∈ A{√(∑_i = 1 ^m (y^*_i - y_i)^2)}

where P^* is a set of representative solutions in the ideal PF, A is an approximate PF achieved by a CMOEA. The value of IGD denotes the distance between P^* and A. For CMOPs with two objectives, 1000 points are sampled uniformly from the true PF to construct P^*. (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 P^*. It is worth noting that a smaller value of IGD represents better performance with regards to both diversity and convergence.



  * Hypervolume (HV):

HV reflects the closeness of the non-dominated set achieved by a CMOEA to the real PF. The larger HV means that the corresponding non-dominated set is closer to the true PF. 


    HV(S)=VOL(⋃_x∈ S [f_1(x),z_1^r]× ...[f_m(x),z_m^r] )

where VOL(·) is the Lebesgue measure, 𝐳^r=(z_1^r,...,z_m^r)^T is a reference point in the objective space. For a LIR-CMOP, 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 HV represents better performance regarding both diversity and convergence. 



 







 §.§ Discussion of Experiments



  §.§.§ Performance comparison on LIR-CMOP test suite

The statistical results of the IGD values on LIR-CMOP1-14 achieved by five CMOEAs in 30 independent runs are listed in Table <ref>. As discussed in Section <ref>, LIR-CMOP1-4 have large infeasible regions in the entire search space. For these four test instances, MOEA/D-IEpsilon is significantly better than the other four tested CMOEAs in term of the IGD metric. Fig. <ref>(a)-(b) shows the final populations achieved by each CMOEA with the best IGD values during the 30 runs on LIR-CMOP1 and LIR-CMOP4. It is clear that MOEA/D-IEpsilon has the best performance regarding diversity among the five CMOEAs under test.

LIR-CMOP5-12 have large infeasible regions, as discussed in Section <ref>. It can be observed that MOEA/D-IEpsilon is significantly better than the other four tested CMOEAs on NCMOP5-12. The final populations achieved by each CMOEA on LIR-CMOP9 and LIR-CMOP11 with the best IGD values are plotted in Fig. <ref>(c)-(d). For LIR-CMOP9, MOEA/D-Epsilon, MOEA/D-SR, MOEA/D-CDP and C-MOEA/D only achieve a part of the real PF. However, MOEA/D-IEpsilon can obtain the whole real PF. Thus, MOEA/D-IEpsilon performs better than the other four CMOEAs in terms of diversity. For LIR-CMOP11, the proposed method MOEA/D-IEpsilon can achieve the whole PF. However, the other four CMOEAs do not converge to the whole PF. Thus, MOEA/D-IEpsilon has better convergence performance than the other four CMOEAs. For three-objective test instances (LIR-CMOP13 and LIR-CMOP14), MOEA/D-IEpsilon is also significantly better than the other four CMOEAs. 

Table <ref> shows the results of the HV values of LIR-CMOP1-14 achieved by five CMOEAs in 30 independent runs. It is clear that MOEA/D-IEpilon is significantly better than the other four CMOEAs on all of the fourteen test instances in terms of the HV metric. 



  §.§.§ Analysis of Experimental Results
 
From the above performance comparison on the fourteen test instances LIR-CMOP1-14, it is clear that MOEA/D-IEpsilon has better diversity and convergence performance than the other four decomposition-based 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 constraint-handling method can deal with the large infeasible regions very well using its automatic adjustment of the epsilon level.














§ ROBOT GRIPPER OPTIMIZATION



To verify the capability of MOEA/D-IEpsilon to solve real world optimization problems, a robot gripper optimization problem with two conflicting objectives and eight constraints is explored.



 §.§ Definition of the robot gripper optimization



The robot gripper optimization problem is defined in (<cit.>). 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. <ref>.








The robot gripper optimization problem considered in this paper is defined as follows:

    f_1(𝐱) = P/min_zF_k(𝐱,z)
       f_2(x)=a+b+c+e+l 
        c_1(𝐱) =Y_min-y(𝐱,Z_max) ≥ 0 
        c_2(𝐱) = y(𝐱,Z_max) ≥ 0
        c_3(𝐱) = y(𝐱,0)-Y_max≥ 0
        c_4(𝐱) = Y_G-y(𝐱,0) ≥ 0
        c_5(𝐱) = (a+b)^2-l^2-e^2 ≥ 0
        c_6(𝐱) = (l-Z_max)^2+(a-e)^2 ≥ b^2
        c_7(𝐱) = l-Z_max≥ 0
        c_8(𝐱) = min F_k(𝐱,z)-F_G ≥ 0

where 𝐱=[a, b, c, e, l, f, δ]^T has seven decision variables, and each variable is shown in Fig. <ref>. The range of each decision variable is as follows: 10mm≤ a≤150mm, 10mm≤ b≤150mm, 100mm≤ c≤200mm, 0mm≤ e≤ 50mm, 10mm≤ f≤ 150mm, 100mm≤ l≤ 300mm and 1.0≤δ≤ 3.14. Two rules are applied to fix the value of f, and they are defined as follows:



    Rule1: if (a<4b and c<a+b) then f=2e+10


    Rule2: if (a<4b and c>a+b) then f=e+50





According to the geometric analysis, the gripping force F_k in Fig. <ref> can be defined as follows:

    F_k=Pbsin(α+β)/2ccosα.


The displacement of the gripper end is defined as follows:


    y(𝐱,z)=2[e+f+c+sin(β+δ)].

where g=√((l-z)^2+e^2)+ϕ, α=arccos(a^2+g^2-b^2/2ag), β=arccos(b^2+g^2-a^2/2bg)-ϕ,ϕ=arctane/l-z and z denotes a dynamic displacement of the gripper actuator in the range of 0 to 100 mm.

The first objective f_1(x) represents a force transmission ratio between the actuating force P and the minimum gripping force min F_k(𝐱,z). We prefer to transform more actuating force into the gripper force. Thus, this objective should be minimized.

The second objective f_2(x) is the sum of all elements of the robot gripper. It is relevant to the weight of the robot gripper, and minimizing f_2(x) 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/D-IEpsilon. In Fig. <ref>, we can observe that the robot gripper optimization problem has large infeasible regions (RFS=0.1396), which can be solved well by the proposed method MOEA/D-IEpsilon according to our previous analysis. To verify this hypothesis, MOEA/D-IEpislon and the other four decomposition-based CMOEAs are tested on the robot gripper optimization problems.









 §.§ Experimental study





  §.§.§ Experimental settings

To solve the robot gripper optimization problem and evaluate the performance of the proposed MOEA/D-IEpsilon, five decomposition-based CMOEAs, including MOEA/D-IEpsilon, MOEA/D-Epsilon, MOEA/D-SR, MOEA/D-CDP and C-MOEA/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 <ref> 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 (<cit.>) to measure the performance of the five tested CMOEAs. In the robot gripper optimization case, the reference point z^r = [5,800]^T.




  §.§.§ Analysis of experiments


Table <ref> shows the statistical results of HV values of MOEA/D-IEpsilon and the other four CMOEAs on the robot gripper optimization problem. It is clear that MOEA/D-IEpsilon is significantly better than the other four CMOEAs. To further demonstrate the superiority of the proposed method MOEA/D-IEpsilon, the non-dominated solutions achieved by each CMOEA during the 30 independent runs are plotted in Fig. <ref>(a)-(e). The box plot of HV values of the five CMOEAs is shown in Fig. <ref>(f). From Fig. <ref>, we see that MOEA/D-IEpsion has better performance than the other four CMOEAs. 










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 non-dominated solutions achieved by MOEA/D-IEpsilon as shown in Fig. <ref>. The configurations of the robot gripper mechanism at each point are also plotted in Fig. <ref>. 

To measure the minimum gripping force min_z F_k(𝐱,z), 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 10^13 N/m.

Table <ref> shows the simulation results of the minimum gripping force min_z F_k(𝐱,z) with three different configurations of the robot gripper. The relative errors between the theoretical gripping forces and the simulated gripping forces are less than 0.1 %. Thus, we can conclude that the optimization results of the robot gripper optimization problem achieved by MOEA/D-IEpsilon are correct. 






§ CONCLUSION


This paper proposes an improved epsilon constraint-handling method embedded in the framework of MOEA/D. A new CMOEA named MOEA/D-IEpsilon has been proposed. The comprehensive experimental results indicate that MOEA/D-IEpsilon has the ability to cross the large infeasible regions. Compared with the other four decomposition-based CMOEAs including MOEA/D-Epsilon, MOEA/D-SR, MOEA/D-CDP and C-MOEA/D, MOEA/D-IEpsilon has following advantages:


  
  * The performance of MOEA/D-IEpsilon is not sensitive to the initial epsilon value.
  
  * MOEA/D-IEpsilon has the ability to explore the feasible and infeasible regions simultaneously during the evolutionary process.
  
  * MOEA/D-IEpsilon 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/D-IEpsilon is suitable for solving CMOPs with large infeasible regions.


In terms of CMOPs, a new set of CMOPs named LIR-CMOP1-14 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/D-IEpsion is significantly better than the other four CMOEAs on this test suite. Thus, we hypothesize that MOEA/D-IEpsilon is better than the other four CMOEAs in solving CMOPs with large infeasible regions, in general. To demonstrate the capacity of MOEA/D-IEpsilon 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/D-IEpsilon outperformed the other four CMOEAs.

Proposed further work includes studying new constraint-handling 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 constraint-handling methods in different evolutionary stages. 



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 High-Level University Project “Green Technologies” for Marine Industries.













§ APPENDIX

   
In this section, the detailed definitions of LIR-CMOP1-14 are listed in Table <ref>.






§ COMPLIANCE WITH ETHICAL STANDARDS


Conflict of Interest The authors declare that they have no conflict of interest.

Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.




























spbasic


