An improved constrained multi-objective optimization evolutionary algorithm for carbon fibre drawing process

ABSTRACT In this paper, an improved ϵCMOEA/D-DE is proposed to improve the performance of CMOP algorithm and achieve parameter optimization in carbon fibre drawing process. In order to avoid overusing the infeasible solutions, two repair operators are introduced in the population evolution model. More specially, during the evolutionary process, when the constraint violation of infeasible solution exceeds a tolerance threshold, the proposed two repair operators are used to find a better solution to repair the infeasible solution. On the other hand, in order to enhance the convergence rate of the Differential Evolution (DE), a modified DE is proposed. Then, an ϵCMOEA/D-mDE-RO is proposed by incorporating these two improved strategies into the ϵCMOEA/D-DE. Subsequently, the performance of the proposed ϵCMOEA/D-mDE-RO is evaluated on the constrained test problem series and the experiment shows that the proposed algorithm outperforms the existing ϵCMOEA/D-DE. Finally, in order to further illustrate the application potential, the proposed algorithm is successfully applied in optimizing carbon fibre drawing process and the optimal draw ratio vector is obtained.


Introduction
The multi-objective optimization problem is common in the practical engineering area such as knapsack problems (Ponsich, Jaimes, & Coello, 2013), power system frequency control (Wang et al., 2017), carbon fibre drawing process (Chen, Ding, Jin, & Hao, 2013), and so on. It is well known that the evolutionary algorithm is an effective tool to solve the multi-objective optimization problem and, in the past decades, a variety of evolutionary algorithms have been proposed/developed such as genetic algorithm (Deb, Pratap, Agarwal, & Meyarivan, 2002), differential evolution (DE) algorithm (Bandyopadhyay & Mukherjee, 2015;Wang, Liao, Zhou, & Cai, 2014), extremal optimization algorithm (Li, Lu, Zeng, Wu, & Chen, 2016;Lu, Zhou, Zeng, & Du, 2018), PSO algorithm Zeng, Wang, Zhang, & Alsaadi, 2016;Zeng, Zhang, Liu, Liang, & Alsaadi, 2017), and immune algorithm (Khaleghi, Farsangi, Nezamabadi-Pour, & Lee, 2011). Note that the evolutionary algorithms mentioned above are only able to deal with the multi-objective optimization problems without constraints or with box-constraints. However, in most practical engineering applications, the objective functions to be optimized are accompanied with the CONTACT Bo Shen bo.shen@dhu.edu.cn, shenbodh@gmail.com complex equality or inequality constraints. Therefore, it is of practical significance to find a suitable constraint handling technology (CHT) so as to obtain an optimal solution to the constrained multi-objective problem (CMOP).
In recent years, a great deal of effort has been made on the CHTs. The developed CHTs mainly focus on the constrained single-objective optimization problem and the CHTs for the CMOP are relatively few. Among the existing CHTs for the CMOP, the penalty function is the most commonly used and the CMOP is usually solved by transforming it into an unconstrained multi-objective optimization problem. For example, in Riche, Knopf-Lenoir, and Haftka (1995), a genetic algorithm based segregated penalty function has been proposed and it has been shown that the performance is less sensitive to the penalty parameters. In Carlson and Shonkwiler (1998), an annealing penalty has been proposed by using the simulated annealing algorithm and a variable penalty function related to the temperature. Furthermore, in Hamida and Schoenauer (2002), an adaptive penalty function has been been proposed which is helpful for making the algorithm search in the regions around the boundary parts. Nevertheless, it should be pointed out that in the algorithms based on the above penalty function, the infeasible solutions are still overused (Jan, Khanum, & Tairan, 2016).
On the other hand, the ε-constraint theory based approach is capable of tuning parameters appropriately. For example, in Becerra and Coello (2016), an εconstraint theory based hybrid multi-objective optimization algorithm has been proposed by using the evolutionary single-objective optimizer and the diversity of Pareto non-dominated solutions has been enhanced. In Takahama and Sakai (2010), an ε-constraint DE has been proposed by combining the gradient-based mutation and the algorithm robustness has been improved. Recently, in Yang, Cai, and Fan (2014), the ε-constraint theory and decomposition technology have been combined and the εCMOEA/D-DE has been proposed to solve the CMOP. However, in the εCMOEA/D-DE, the infeasible solutions have still been overused in the evolutionary process which may reduce the searching efficiency and the convergence rate of algorithm.
As analysis above, for a CMOP, two important performance indices should be concerned, i.e. searching efficiency and convergence rate. In order to enhance the searching efficiency, an alternative method is to use the infeasible solutions reasonably. In this paper, two repair operators are proposed to avoid overusing the infeasible. Specifically, during the evolutionary process, when the constraint violation of infeasible solution exceeds a tolerance threshold, the proposed two repair operators are employed to find a better solution to repair the infeasible one. With respect to convergence rate, we propose a modified DE, which extends the traditional DE in εCMOEA/D-DE in order to improve the convergence rate. By combining the proposed repair operators and modified DE, an improved εCMOEA/D-DE, i.e. the εCMOEA/D-mDE-RO is proposed to enhance the algorithm performance for the COMPs.
On the other hand, the carbon fibre drawing process can be viewed as a typical CMOP. Generally, the drawing process consists of six drawing steps, i.e. spinneret drawing, air drawing, DMF coagulation drawing, hot water drawing, boiling water drawing and third-class drawing (Chen, Ding, & Hao, 2013). In order to realize the optimization in the characters of the carbon fibre such as density, strength and breaking elongation ratio, the multiple draw ratios must be designed simultaneously. Actually, for multi-objective optimization problem for the carbon fibre drawing process, there have been a number of advanced approaches available in the existing literatures, see e.g. (Chen et al., 2013;Kalantari, Dong, & Davies, 1999;Salmalian, Nariman-Zadeh, Gharababei, Haftchenari, & Varvani-Farahani, 2010). In order to test the performance of the εCMOEA/D-mDE-RO proposed in this paper, we try to apply the εCMOEA/D-mDE-RO into the multi-objective optimization problem for the carbon fibre drawing process.
Based on the above discussion, the main contributions of this paper can be summarized in the following two aspects.
(1) A new εCMOEA/D-mDE-RO is proposed where the improved two repair operators and the modified DE is capable of (a) tackling the problem of overusing infeasible solutions during the evolutionary process; (b) speeding up significantly the searching convergence; and (c) achieving a proper balance between the use of infeasible solutions and the optimal performance in the searching process.
(2) The proposed εCMOEA/D-mDE-RO is successfully employed in the optimization of the carbon fibre drawing process and a set of optimal drawing ratio vector has been obtained which provides a reference for the practical production of the carbon fibre.
The remainder of the paper is organized as follows. Section 2 provides the basic knowledge of COMP. Section 3 gives the description of εCMOEA/D-DE. Section 4 details the proposed εCMOEA/D-mDE-RO. In Section 5, the comparative experimental results are provided and discussed. Section 6 show the comparative experimental results for the carbon fibre drawing process by applying the proposed εCMOEA/D-mDE-RO. Finally, conclusions are presented in Section 7.

Preliminaries
The CMOP under consideration is described as follows x n ] T ∈ R n is an n dimensional decision vector and x lower and x upper are the lower and upper bounds of all decision variables x i (i = 1, 2, . . . , n). In solving the above MOPs, the Pareto theory is often used. Some notions associated with the Pareto theory are given as follows (Martinez & Coello, 2014;Miettinen, 1999).
Definition 2.1: ∀ x, y ∈ , x is said to dominate y or y is said to be dominated by x (denoted as x ≺ y) if and only if f i (x) ≤ f i (y) for all i = 1, 2, . . . , m and there is at least one i (i = 1, 2, . . . , m) such that f i (x) < f i (y). Definition 2.2: ∀ x * ∈ , x * is said to be a Pareto nondominated solution or Pareto optimal solution if and only if there is no other solution y ∈ such that y ≺ x * .

Definition 2.3:
The Pareto optimal set is expressed as PS = {x ∈ |x is Pareto non-dominated solution}, and the Pareto optimal front is defined by PF = {F(x)|x ∈ PS}.

The description of εCMOEA/D-DE
As mentioned above, the proposed algorithm is improved based on εCMOEA/D-DE (Martinez & Coello, 2014). Thus, in this subsection the εCMOEA/D-DE is described firstly. The εCMOEA/D-DE integrates classical MOEA/D-DE (Zhang & Li, 2007) and a an ε-constraint theory based CHT. The main idea of εCMOEA/D is to transform a CMOP into several scalar problems via the Tchebycheff approach and then to obtain the Pareto optimal solutions by using the ε-constraint theory. Meanwhile, the DE operator and the polynomial mutation (PLM) operator (Tang & Tseng, 2013) are used to obtain the corresponding PS and make the optimization objectives converge to the obtained PF. Some relevant theories and methods involved in the algorithm are outlined as follows.
The Tchebycheff method is adopted to transform a CMOP into several scalar problems. Specifically, by using the Tchebycheff approach, a CMOP is decomposed into N scalar optimization subproblems by the N well-distributed weight vectors, denoted as W = {w 1 , w 2 , . . . , w N }. The scalar function of the ith subproblem is defined by where w i ∈ W represents weight vector of the ith subproblem, w i j represents the weight w i corresponding to the jth dimensional objective, z * = (z * 1 , z * 2 , . . . , z * m ) T is the reference point and z * j = min{f j (x)|x ∈ R n } is the best value looked for so far for the objective f j .
When dealing with constraints in multi-objective optimization problems, the penalty function approach is widely adopted, which is defined as follows where s is a penalty scaling parameter and (x) is the constraint violation for the solution x. (x) is given by where β is a positive number, g i (x) and h j (x) represent inequality constraints and equality constraints, respectively. It is worth mentioning that the equality constraints h j (x) can be converted into the inequality constraints according to |h j (x) − ξ | ≤ 0, where ξ is a small threshold which can be selected similarly in Martinez and Coello (2014). If the constraint violation of decision vector x is greater than zero, the vector x is viewed as an infeasible solution. However, the algorithm based on penalty functions has some disadvantages such as difficulties in tuning parameters and overusing of the infeasible solutions. Therefore, we adopt ε-constraint theory to deal with the constraints. This approach depends on the constraint violation of each individual and the constraint violation degree is defined by the neighbourhood of this individual. To define the constraint violation degree, the definition of neighbourhood of each individual is firstly given by where α is an empirical parameter and the neighbourhood B(i) contains the indexes of M closest weight vectors of the ith subproblem. Then, the constraint violation degree ε (Martinez & Coello, 2014) can be defined as where τ is an adjustable threshold and i,max and i,min are the maximum and minimum constraint violation both taken from neighbouring x l , l ∈ E, and can be described by Based on the constraint violation of each individual, the ε-constraint theory is used to evaluate the individual quality. The evaluation function sequence (f , ) is composed of objective function f and constraint violation . Hence, for any ε ≥ 0, the ε-comparison '≤ ε ' between the sequence (f 1 , φ 1 ) and the sequence (f 2 , 2 ) can be defined by Note that, when ε = ∞, ε-comparison ≤ ε is simplified as the unconstrained multi-objective optimization problem where only the objective functions are compared.
The DE operator and the PLM operator are used to make the optimization objectives converge to PF. The DE operator (Storn & Price, 1997) is defined as where CR and F are two empirical parameters, x r 1 k is a current individual in the kth generation, and x r 2 k and x r 3 k are taken from the neighbourhood whose indexes are r 2 and r 3 , respectively.
The PLM operator (Tang & Tseng, 2013) is defined as: q is a positive number, r ∈ rand(0,1), , and x LB k and x UB k represent the lower and upper bounds of x k , respectively.

An novel εCMOEA/D-DE algorithm: εCMOEA/D-mDE-RO
In this section, a new εCMOEA/D-mDE-RO algorithm is proposed to further improve the search ability of the original εCMOEA/D-DE algorithm. The main novelty of the proposed εCMOEA/D-mDE-RO lies in the introduction of two repair operators and a modified DE operator into the population evolution model. More specifically, when the constraint violation of infeasible solution exceeds the tolerance threshold during the evolutionary process, the proposed two repair operators are used to find a better solution to repair the infeasible one. Furthermore, in order to overcome the disadvantage of the low convergence rate in traditional DE, a modified DE is proposed to improve the convergence rate. Compared to the original εCMOEA/D-DE algorithm, the newly introduced two repair operators and the modified DE operator in the population evolution model make it possible for us to (1) make better use of infeasible solutions during the evolutionary process; (2) pursue strong convergence capability; (3) keep an adequate balance between the convergence and the diversity. As such, the proposed εCMOEA/D-mDE-RO could not only exploit the infeasible solutions more reasonably than the εCMOEA/D-DE algorithm, but also has higher convergence rate.

The first repair operator
If the constraint violation of the ith individual x i (the individual a) satisfying i > γ × ε, we perform the repair operator on this individual. In the neighbourhood of x i , if the ratio of feasible solutions (RFS) satisfies RFS > 0.9, we consider this neighbourhood to be high-quality, and then we adopt the strategy that a solution y (the individual b) is selected via the minimum objective value g te from this neighbourhood, as shown in Figure 1(a). The RFS indicator refers to the number of feasible solutions found in the neighbouring R. This practice can be described as follows where E consists of the the subscript of each neighbourhood individual, R consists of each neighbourhood individual and y is the solution with minimum objective value g te in the neighbouring R.
If RFS ≤ 0.9, we select a set of solutions satisfying the constraint violation l < γ × ε, which is denoted by T. Then we use the objective function sorting to find a solution y (the individual b) via the minimum objective value g te , that is, the solution y is the one nearest to x i in T. Finally, let the solution y repairs a. This practice can be described as follows The principle of above approach is to make the algorithm concentrate to search in the infeasible region beside the boundary parts and then realize the purpose of using infeasible solutions efficiently and reasonably, which is depicted in Figure 1(b). However, if T couldn't be found, it means that the constraint violation of all solutions in neighbourhood E satisfies l > γ × ε or l = 0 and the number of infeasible solutions is not less than two (if M = 20). Therefore, we consider the neighbourhood E to be low-quality relatively, and then select a solution which is the farthest away from x i . In Figure 1, the solid line in upper picture is the feasible region (FR), the solid line in the lower picture is the infeasible region (IR). Figure 1(a) shows the repair process of the case of RFS > 0.9 where the individual b replaces the individual a. Figure 1(b) shows the replacement process satisfying preconditions that RFS < 0.9 where the neighbouring T can be found and then the individual b repairs the individual a. The region between solid line and dotted line is the infeasible region beside the boundary parts, and the infeasible solutions in this region are deemed the high-quality ones. By this approach, the infeasible solutions are used reasonably and the algorithm performance is hence improved.

The second repair operator
The main idea of second repair operator is to inhibit the excessive convergence in the process of population evolution. If constraint violation i of the infeasible solution x i exceeds γ × ε, an better solution y (the individual b) is taken from its neighbouring U to repair it. More concretely, the solution y that is the closest to the Pareto optimal solution in the previous generation is selected, and then repairs the inferior x i . The main steps are given as follows.
(1) If i > γ × ε, we select a set of solutions which satisfies the constraint violation < i , denoted by neighbouring U.
(2) Then, we use objective function sorting to find a solution y via the minimum objective value g te . (3) Finally, the solution y repairs the original solution x i .
The notations U and y are formulated as follows In the above equations, the weight of x i , w i remains to be constant which is different from the one in Equation (13) and Equation (15) where the weight are varying. The advantages of this repair operator are that the algorithm is inhibited to generate individuals within excessive constraint violations at the beginning of the algorithm evolutionary process. Additionally, the algorithm will accumulate the beneficial effect in the subsequent generations, and eventually improve the algorithm performance. The replacement process is depicted in Figure 2.

The modified DE operator
We propose a modified DE operator which has better search efficiency and convergence rate according to the design requirement of practical engineering problems. The modified DE operator and the two repair operators proposed above are combined in the proposed algorithm. The results on the CTP-series test instances show that the proposed DE operator can effectively speed up the convergence of the algorithm while improve the expense of diversity.
The main characteristic of the DE operator is mutation which is also seen as a variation or perturbation modelled by a random variable. The DE operator widely used in MOEA/D is described as follows.
where F is the scaling factor, x r 1 k is the decision variable vector, x r 2 k and x r 3 k are the donor vectors, r 2 and r 3 are the indexes that are randomly selected from the neighbourhood E and y k is a trial vector generated by the DE operator. It can be seen that, if rand (0, 1) < CR, the DE operator is performed on x r 1 k , see Figure 3(a). Considering the design requirement of the practical engineering problems, we propose a modified DE operator which has better search efficiency and convergence rate. The strategy can be outlined as follows. A base vector x best is selected as the best solution in current population. Then, by using DE operator generates the intermediate vector y * k , let y * k inherit any of decision variables from x r 1 k (the dimension n of decision variables satisfying n ≤ 10). Finally, the operator generates the trial vector y k(x 1 ) or y k(x 2 ) , as shown in Figure 3(b). Figure 3(a) represents the search space of a simple DE operator in 2-D and Figure 3(b) represents the search space of a modified DE operator. In Figure 3, y k(x 1 ) and y k(x 2 ) are the trial vectors which inherits x r 1 k and y * k is an intermediate vector.

The framework of εCMOEA/D-mDE-RO
The framework of εCMOEA/D-mDE-RO1 is given in Algorithm 1. Line 15 and Line 29 ∼ Line 33 are the modified parts where Line 15 describes the modified DE operator (mDE) and Line 29 ∼ Line 38 are the pseudo-code of the first repair operator (RO1) that we proposed. Line 1 ∼ Line 28 and Line 34 ∼ Line 41 are the εMOEA/D-mDE where if Line 15 performs the DE operator on (x r 1 , x r 2 and x r 3 ), the algorithmic framework can be simplified as the εMOEA/D-DE. Specifically, When the constraint violation of an individual satisfies > γ × ε (Line 29) and RFS satisfies RFS > 0.9 (Line 30), the neighbouring R (Line 31) is obtained in terms of the specified condition, and then a superior individual is selected from the neighbouring R (Line 32). If RFS ≤ 0.9 (Line 33), the neighbouring T (Line34) is obtained by the specified condition and then a superior individual with the minimum objective value g te is obtained from the neighbouring T (Line 35). At last, the superior individual repairs the original individual (Line 9).
The framework of εMOEA/D-mDE-RO2 can be described by replacing the pseudo-code part of Line 29 ∼ Line 38 in Algorithm 1 with Algorithm 2. When the constraint violation of an individual exceeds γ × ε (Line 1), the neighbouring U (Line 2) is obtained by the specified condition. Then, a superior individual with the minimum objective value g te is obtained from this neighbouring U (Line 3). Finally, the superior individual repairs the original one (Line 4).
Algorithm 1: The Framework of εMOEA/D-mDE-RO1 Input: G: the maximum iterations N: the number of subproblems in εMOEA/D-bDE-CHT2 W: a group of weight vectors w 1 , . . . , w N M: the number of weight vectors contained in the neighbourhood δ: a probability decides whether parent solutions are chosen from the neighbourhood B(i) or not n r : a parameter about replacement number in the neighbourhood τ : a parameter used in defining the constraint violation degree ε γ : tolerant threshold Output: P: the final population found by εMOEA/D-bDE-CHT2 1 : Initialize a random population P ← x 1 , . . . , x N 2 : the M closest weight vectors to w i 5 : end for 6 : z * ← (z * 1 , . . . , z * m ) T 7 : while the maximum iterations G is not met do 8 : for i ← 1 to N do 9 : if rand(0, 1) < δ then 10 : E ← B(i) 11 : else 12 : E ← {1, 2, . . . , N} 13 : end if 14 : Set r 1 ← i and select two indexes r 2 , r 3 ∈ E randomly 15 : y c ← Differential Evolution operator on (x r1 , x r2 , x r3 and x best ), x best is the best solution in current pollution 16 : y ← Polynomial mutation operator on y c 17 : Calculate F(y) 18 : for j ← 1 to m do 19 : if f j (y) < z j then 20 : z j ← f j (y) 21 : end if 22 : end for 23 : c ← 0 24 : for l ← 1 to M do 25 : if RFS > 0.9 then 31 : R = {x l , l ∈ E} 32 : y gd = {x l | min(g te (x l |w l , z * )), x l ∈ R}, where w, z * is same to the ith 33 : else 34 : y gd = {x l | min(g te (x l |w l , z * )), x l ∈ T} 36 : end if 37 : y ← y g d 38 : end if 39 : if (g te (y|w l , z * ), i (y)) ≤ ε (g te (y|w l , z * ), i (x l )) and c ≤ n r then 40 : x l ← y 41 : where w, z * is same to the ith 4 : y ← y g d 5 : end if

Experimental setting
In this subsection, we first introduce the experiment circumstance including the test instances with related parameters, the definition of performance index, the parameters used in the proposed algorithm and the operating circumstance.
Additionally, we use the popular performance metric Hypervolume (HV) (Liu et al., 2016) as the performance index. HV metric can demonstrate both the convergence and diversity of Pareto non-dominated solutions in a sense (see in Figure 4). In Figure 4, we assume that Q = A, B, C, D is the Pareto non-dominated solutions set and thus the HV is the ABCDW enclosed by the discontinuous boundary, where W is the reference point. If the reference point is (0, 0) and test instances are minimized, then the smaller value is the HV metric, the better convergence and diversity are the algorithm. This can be described by the following formula: where |Q| is the number of the Pareto non-dominated solutions, v i is the diagonal corners of hypercube constructed with the reference point W and solutions Q.

Parameters setting in εCMOEA/D-mDE-RO and operating circumstance
In this experiment, the parameters are given as follows: population size N = 200, the dimension of decision variable n = 10, the adjustable threshold τ of constraint violation degree ε = 0.2, the initial reference point z * = (+∞, . . . , +∞) T , the neighbourhood size M = 0.1N = 20, the maximum number of replacements in the neighbourhood n r = 0.01N = 2, the empirical parameter in Equation (5) α = 0.9, the empirical parameters in DE operator CR = 1.0 and F = 0.5, the empirical parameter and implemented probability of PLM q = 2 and p m = 1/n, and a user-defined maximum number of generations G = 50. We independently run the algorithm for 30 times for each test instance on an identical computer (Inter(R) Core(TM) i7−6700 CPU @ 3.40GHz, Matlab 9.0).

Experimental results and analysis
In this section, we implement three comparative experiments to evaluate the performance of the introduced εCMOEA/D-mDE-RO algorithm. The first experiment is to test the performance of the

Comparison εCMOEA/D-DE-RO and εCMOEA/D-DE
In order to test the sensibility of the tolerance value γ , we set the parameter γ = 2 and γ = 10. The simulation results are shown in Figure 5. It can be seen from Figure 5 that, equipped with the two repair operators as the CHT approach, the algorithm performance is significantly improved in CTP6 ∼ CTP8; and the second repair operator approach is better than the first repair operator approach.
The corresponding HV values are listed in Table 2. From Table 2, it is seen that the value with an unserscore is a better result obtained by εCMOEA/D-DE-RO1, the italic value is a better result obtained by εCMOEA/D-DE-RO2 and the value with overstriking is the best result in each CTP test instance. Table 2 indicates that εCMOEA/D-DE-RO2 (γ = 10) has the best algorithm performance in terms of HV values.

Comparison εCMOEA/D-mDE and εCMOEA/D-DE
By using the modified DE operator, the convergence is enhanced comparing with the traditional DE. The corresponding results are shown in Figure 6. It can be been from Figure 6 that the εCMOEA/D-mDE has higher search efficiency and convergence rate than εCMOEA/D-DE.

Comparison εCMOEA/D-mDE-RO and εCMOEA/D-DE
After the above experimental results and analysis, we incorporate the proposed two strategies into εMOEA/D (i.e. the εCMOEA/D-mDE-RO) and the simulation results of the εCMOEA/D-mDE-RO and the traditional εCMOEA/D-DE are given in Figures 7 and 8. From Figure 7, we can see that the proposed algorithm has better performance in all CTP-series test instances. Figure 8 shows a group of Pareto optimal fronts (PFs) obtained by εCMOEA/D-mDE-RO, when the number of generations satisfying G = 200 and from Figure 8, it can be seen that the negative effects of the traditional DE operator disappears when the number of generations is large. That is, the algorithm deeply converges to the global optimal solutions.

The application of εCMOEA/D-mDE-RO into optimization of carbon fibre drawing process
The carbon fibre drawing process is a typical CMOP. In the carbon fibre production process, precursor fabrication is a critical factor to ensure carbon fibre products to meet high-quality standards, and the drawing process is a crucial working procedure during precursor fabrication. Generally, the drawing process consists of six drawing steps, i.e. spinneret drawing ratio, air drawing ratio, DMF coagulation drawing ratio, hot water drawing ratio, boiling water drawing ratio and third-class drawing ratio (Chen et al., 2013;Wang, Ni, & Liu, 2000). In order to realize the purpose of optimizing the fundamental characters, e.g. density, strength and breaking elongation ratio, the multiple draw ratios should be selected reasonably. In this section, the proposed algorithms are employed to optimize the drawing ratios in carbon fibre drawing process, as an application example of the proposed algorithms.

Optimization problem
As analysis in Chen et al. (2013) and Wang et al. (2000), the optimization problem for the density and strength of x l ∈ [10, 18], l = 1, 2, . . . , 6 where x l is the drawing ratio of six drawing stages, x mul is the total drawing ratio, f 2LTH is the lower boundary and f 1 , f 2 are, respectively, the objective functions representing the density, strength of carbon fibre which are given by In this example, the weight of drawing process is taken as w = [0.3102, 0.1703, 0.0895, 0.1703, 0.1703, 0.0895] and the lower boundary is chosen as f 2LTH = 10 × 10 3 .

Experimental results
In this subsection, we solve the optimization problem and obtain the optimal parameters of drawing ratio and the Pareto fronts (PFs). Then, the obtained parameters are compared with the values given in Chen et al. (2013). Because in Chen et al. (2013) the results are obtained after iteration times G = 50, we set the same iteration times in all our experiments, and the algorithms are independently run for 30 times. Note that the optimization problem (22) is maximized. Therefore, the bigger average HV-metric values, the better algorithm performance. Firstly, we test the convergence on the carbon fibre drawing optimization problem, the results of all the proposed modified algorithms and εCMOEA/D-DE are shown in Figure 9.
It is seen from Figure 9 that all the modified algorithms have better convergence than the traditional εCMOEA/D-DE and the εMOEA/D-mDE-RO2 has the best performance in terms of efficiency and convergence. Secondly, in Figure 10, all PFs obtained by modified and εCMOEA/D-mDE-RO are presented. It can be seen from Figure 10 that, compared with εCMOEA/D-DE, the modified algorithms obtain better optimal solutions.
In Figure 10, points A, B, C are three optimal ones found by εCMOEA/D-mDE-RO1 where points A and C represent the best f 2 and f 1 , respectively and B is an intermediate . Furthermore, we compare them with the unconstrained optimal points obtained by SICSA, i.e. the two unconstrained optimal points G and H taken from (Chen et al., 2013) and (Chen et al., 2013), respectively. Comparison results are shown in Table 3. It can be seen from Table 3 that all points obtained by our proposed modified algorithms are better than the ones obtained by SICSA. Therefore, it is believed that the obtained drawing ratios can serve as a reference and a guidance for the practical production of carbon fibre, which also verifies the effectiveness of the proposed algorithms in this paper.  Table 3. Objective functions and drawing ratios obtained by different algorithms.

Results
Optimal Points x 1 x 2 x 3 x 4 x 5 x 6 x 7

Conclusion
In this paper, the improved εCMOEA/D-DE has been proposed to solve the CMOP. In order to avoid overusing the infeasible solutions in the evolutionary process, two repair operators have been proposed, and the DE operator has been modified to improve the search ability of convergence. Then the newly proposed repair operators and the modified DE operator have been incorporated into the εCMOEA/D-DE and an improved CMOP has been proposed. Then, three experiments have been implemented to demonstrate the performance of the improved CMOP. Moreover, the proposed algorithm has been successfully applied in optimizing carbon fibre drawing process and the optimal draw ratio vector has been obtained. Future works would include (1) how to further improve the diversity of the proposed εCMOEA/D-mDE-RO; and (2) how to apply the proposed algorithms to solve other maximized problems in carbon fibre and the more common industrial systems.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported in part by the National Natural