Adaptive Truncation technique for Constrained Multi-Objective Optimization

Abstract

The performance of evolutionary algorithms can be seriously weakened when constraints limit the feasible region of the search space. In this paper we present a constrained multi-objective optimization algorithm based on adaptive ε-truncation (ε-T-CMOA) to further improve distribution and convergence of the obtained solutions. First of all, as a novel constraint handling technique, ε-truncation technique keeps an effective balance between feasible solutions and infeasible solutions by permitting some excellent infeasible solutions with good objective value and low constraint violation to take part in the evolution, so diversity is improved, and convergence is also coordinated. Next, an exponential variation is introduced after differential mutation and crossover to boost the local exploitation ability. At last, the improved crowding density method only selects some Pareto solutions and near solutions to join in calculation, thus it can evaluate the distribution more accurately. The comparative results with other state-of-the-art algorithms show that ε-T-CMOA is more diverse than the other algorithms and it gains better in terms of convergence in some extent.

1. Introduction

 In general, constrained multi-objective optimization problems (CMOP) can be defined as follows.

\(\begin{aligned} &\operatorname{Min} \quad F(X)=\left[f_{1}(X), f_{2}(X), \ldots, f_{m}(X)\right]\\ &\text { s.t. }\left\{\begin{array}{ll} \boldsymbol{X} \in \boldsymbol{\Omega}, & l_{k} \leq x_{k} \leq u_{k} \\ g_{i}(\boldsymbol{X}) \leq 0, & i=1,2, \ldots, p \\ h_{j}(\boldsymbol{X})=0, & j=1,2, \ldots, q \end{array}\right. \end{aligned}\)       (1)

 where X=(x₁,x₂,…,xn) is a n-dimensional decision vector, and Ω is the decision space. l_k and u_k are the lower and upper bound of xk for k=1,…n. The objective vector F(X): Ω→R^m consists of m real-valued objective functions and R^m is the objective space. gi(X)≤0, i=1,…p is the i^th inequality constraint, and hj(X)=0, j=1,…q is the j^th equality constraint.

 The constraint violation is a scalar value of the constraint violation function which can be formulated as:

\(G(\boldsymbol{X})=\sum_{i=1}^{p} \max \left(0, g_{i}(\boldsymbol{X})\right)+\sum_{j=1}^{q} \max \left(0,\left|h_{j}(\boldsymbol{X})\right|\right)\)       (2)

 Therefore, a solution X is called the feasible solution when it satisfies all constraint conditions. Otherwise, X is called the infeasible solution. All the feasible solutions constitute the feasible region S and S∈Ω.

 Definition 1 (Pareto domination). Given two decision vectors X₁, X₂∈S, X1 is said to dominate X₂, denoted by X₁ X₂, if and only if f_i(X₁)≤fi(X₂) for every i∈{1,2,...,m} and f_j(X₁)<f_j(X₂) for at least one index j∈{1,2,...,m}.

 Definition 2 (Pareto solution). A decision vector X*∈S is Pareto optimal to (1) if there is no solution X∈S such that X> X*.

 Definition 3 (Pareto set, PS). Pareto set is defined as: PS={X∈Ω|X is Pareto optimal}.

 Definition 4 (Pareto front, PF). Pareto front is defined as: PF={F(X)∈R^m|X∈PS}.

 Compared with unconstrained multi-objective optimization problems (MOPs), less work have been done for solving CMOPs. In fact, almost all CMOPs arising from the real-world applications have multiple objectives and exist some inequality or enquality constraints. Yet it is difficult to solve CMOPs, because finding feasible solutions may require additional computational resources. The main challenges in CMOPs are distinct limits on constraints, objectives, the interrelationship between objectives and constraints. Due to the extensive applications and great challenges, CMOPs have gradually caused concern.

 The research of CMOPs mainly refers to the evolutionary algorithms (EAs) and the constraint handling techniques [1]. The EAs mainly include genetic algorithm, evolutionary programming and evolutionary strategy [2]. There were some other swarm intelligence EAs like ant colony algorithm [3], particle swarm [4], cultural algorithm [5], differential evolution [6], artificial bee colony [7], gravitational search algorithm [8], biogeography based optimization [9] and teaching-learning-based optimization [10]. In contrast, there’s not a lot of research about the following EAs yet very promising. Deb [11] proposed an improved NSGA-II as NSGA-III, which replaced the crowding distance in NSGA-II with a clustering operation by series of evenly distributed reference lines. Hereafter, a fitness function similar to penalty-based boundary intersection (PBI) function was developed to strength selection pressure for NSGA-III [12]. Zhang suggested a multi-objective evolutionary algorithm based on decomposition (MOEA/D) [13], which decomposed a MOPs into a group of sub-problems with single-objective and then simultaneously evolved these sub-problems. So far, studies based on MOEA/D have been gradually received attention. For instance, it has combined with swarm intelligence [14], local search strategy [15], and adaptive mechanism [16].

 Constraint handling techniques, as the other important branch in CMOPs, have developed from penalty function method [17], [18], stochastic ranking [19], [20], feasibility rule [21], [22] to multi-objective optimization [23], [24], ε-constraint [25], dual population storage [26], and so on. The penalty function methods had the advantage of simple operation [27], yet it was difficult to set up suitable penalty coefficient, which was too small might lead to local optimum or too large might slow down evolution. To overcome above shortcomings of penalty function methods, stochastic ranking [28] was proposed to balance the objective functions and the constraints. But key parameter pf played a great impact on the results of different problems. Feasibility rule [29] emphasized that the feasible solutions were always better than the infeasible solutions. The popularity of this simple constraint handling method lied on its ability to be coupled to various algorithms, without introducing additional parameters. However it was likely to fall into premature convergence on highly constrained problems. Multi-objective optimization method [30] transformed the constraints into one or more new objectives, so the CMOPs was changed into the unconstrained MOPs. Some research results [31] had shown that multi-objective optimization method was not the best choice for handling CMOPs. The ε-constraint [32] coordinated the relationship between feasible solutions and infeasible solutions, so it could ensure diversity of the evolution population. The feasible solutions and the infeasible solutions were respectively stored by two populations in the dual population storage [33], which could avoid the direct comparison between them. Therefore, the diversity would be strengthened by infeasible solutions.

 The core of CMOP algorithms is how to coordinate diversity and convergence. We need to make the balance between the objectives and the constraints, the coordination between feasible solution and infeasible solution. To this end, a constrained multi-objective optimization algorithm based on adaptive ε truncation is proposed to achieve effective balance between diversity and convergence. Our main contributions have three aspects. Firstly, ε-truncation operation is motivated by strength and weakness of two advanced constraint handling techniques (ε-constraint and dual population storage). The ε-constraint emphasizes one-to-one replacement. However when substituted individual is located in sparse region, the selection strategy lacks ability of diversity maintenance. So ε-constraint emphasizes more convergence than diversity. Dual population storage implicitly maintains the diversity by getting infeasible solutions involved in the evolution. But whether the infeasible solutions are conserved in the infeasible solution set is completely determined by their constraint violation. In such case, the presence of infeasible solution with poor objective values will lead to reducing search efficiency. We aim to coordinate the relationship between feasible solutions and infeasible solutions, through preferentially selecting the feasible Pareto optimum solutions and infeasible solutions which own the lower constraint violation and the better objective value. Thus, both diversity and convergence are considered in the proposed ε-truncation technique. Secondly, the improved crowding density estimation only selects a part of Pareto individuals and the near individuals to participate in calculation. It not only reduces amount of computation, but also eliminates harmful effect of distant individuals and weak individuals. So it can assess the diversity more accurately. Finally, exponential variation is introduced after crossover and mutation operation. Some studies have revealed that crossover and mutation operation in differential evolution algorithm can keep good diversity. And exponential variation is characterized by local search. Hence, integrating them will be able to balance the global exploration and the local exploitation.

 The remainder of this paper is organized as follows. Section 2 reviews the previous work in CMOPs. Section 3 describes the proposed ε-T-CMOA in detail. Section 4 presents the test problems, performance metric, parameter settings, parameter analysis, comparative experimental results and discussions, innovative point analysis. Eventually, conclusions and future work are drawn in Section 5.

 

2. Related Work

 Despite a number of latest achievements, the research on CMOPs is far from being completely explored. Especially for CMOPs with severely constrained, multi-modal and high-dimensional, the existing CMOPs algorithms are still not powerful enough, and the requirement for more effective algorithms is pressing. In this section we mainly review related literature based on EAs and constrained handling techniques.

 An adaptive penalty function method was investigated in literature [34], where different fitness functions were constructed for the feasible solution and the infeasible solution according to their corresponding objective values and constraint violation. To some extent, the exploited infeasible solutions contributed to enlarging exploration scope. A constrained evolutionary multi-objective algorithm based on dual population storage was proposed in literature [26]. The feasible solutions were retained in feasible solution set. Meanwhile, the infeasible solutions with the smallest constraint violation were preferentially selected into infeasible set. Similarly, literature [35] suggested a constrained multi-objective algorithm with particle swarm. The Pareto optimal solutions in sparse region were favorably selected to enter feasible solution set, and the infeasible solutions with both the better objective value and the lager crowing distance were choose to enter infeasible solution set. Thanks to the participation of the excellent infeasible solutions and the feasible Pareto solutions, diversity and convergence have been effectively balanced. A novel constraint handling technique called optimal sequence was presented in literature [36], where objective function, diversity metric and constraint violation were aggregated to form a new fitness function, so elite selection, diversity and feasibility had got coordinated. Although it saved the amount of computation, convergence was not accurate. Literature [37] developed a collaborative evolutionary algorithm based on dual population storage. The feasible solutions evolved toward optimality while the infeasible solutions evolved toward feasibility. It realized information sharing, subsequently the diversity was reinforced.

 Up to present, MOEA/D has been successfully applied to MOPs. However, only little research solves CMOPs by means of MOEA/D. This may be owing to the trouble in coordinating diversity with convergence when fitness landscape is severely constrained. Literature [38] had combined ε-constraint with MOEA/D for the first time. But it missed a part of Pareto solutions on the problems that their PF consisted of a set of discrete points. Similar ideas could be found in literature [23], where CMOPs was decomposed into several sub-problems, and each sub-problem corresponded a sub population. Then, all populations were optimized collaboratively. Literature [39] integrated MOEA/D with stochastic ranking and feasibility rule. Yet the issue that multiple sub problems corresponded to the same optimal individual would seriously affect diversity.

 

3. Proposed Algorithm

 Double population storage and ε-constraint are two kinds of advanced constraint handling techniques. The double population storage only uses constraint violation to update the infeasible solution set, which leads to preserving the individuals of poor objective value. Thereby, evolutionary efficiency and convergence speed are hindered. The ε-constraint coordinates feasible solutions and infeasible solutions by adjusting parameter ε. But its one-to-one replacement mechanism may result in loss of a part of Pareto solutions. Because if replaced individuals are located in sparse areas, the search will lose these areas. In this section, a new constrained handling technique called adaptive ε-truncation is suggested to better balance diversity and convergence.

 

3.1 Adaptive ε-Truncation

 Definition 5 (ε-truncation set). Given a population P⊂Ω, ε-truncation set (P≤ε) is defined as: P_≤ε={X∈P|G(X)≤ε}, and P>ε is defined as: P_>ε={X∈P|G(X)>ε}.

 Definition 6 (ε-truncation operation). Let P⊂Ω be a population at generation t, dividing P into P_≤ε and P_>ε is called ε-truncation operation.

\(\begin{aligned} &\varepsilon(t)=\left\{\begin{array}{ll} \varepsilon(0) \times\left(1-2 t / G_{\max }\right)^{2}, & t \leq 0.5 G_{\max } \\ 0, \quad e l s e \end{array}\right.\\ &\varepsilon(0)=0.4 \times \frac{1}{N} \sum_{i=1}^{N} G\left(X_{i}\right) \end{aligned}\)       (3)

where G_max represents the maximum generation, ε(t) is the truncation level in generation t and ε(0) is initial truncation level, N is the population scale.

 Property 1. When ε=0, P_≤ε is a feasible solution set, and P_>ε is an infeasible solution set.

 Proof. As ε=0, through definition 5 and definition 6, P_≤ε indicates the set in which constraint violation equals to zero. Thus, P_≤ε is a feasible solution set, similarly P_>ε represents an infeasible solution set in which constraint violation is greater than zero.

 Property 2. The double population storage is a special case of adaptive ε-truncation. When ε=0, the properties of adaptive ε-truncation is same as double population storage.

 Proof. Suppose ε=0, according to properties 1, P_>ε is a feasible solution set and P_>ε is an infeasible solution set. Therefore, the adaptive ε-truncation can be seen as double population storage.

 The comparison between adaptive ε-truncation and double population storage are discussed as follow. Similarities: 1) Both of them store individuals in two populations, and diversity can be enhanced by allowing infeasible solutions taking part in the evolution; 2) They can avoid direct comparison of the feasible and infeasible solutions. Actually it’s difficult to compare feasible solutions with excellent infeasible solutions in a quantitative way. Differences: 1) Dual population storage is the special case of adaptive ε-truncation when ε is equal to zero. 2) The adaptive ε-truncation retains a small number of excellent infeasible solutions in the evolution population, however dual population storage conserves infeasible solutions by infeasible solution set. 3) In the early, ε is relatively large in adaptive ε-truncation. So the infeasible solutions with the better objective value can be exploited to expand the exploration range, and the diversity will be enhanced. Meanwhile, feasible Pareto solutions are reserved to ensure convergence. In the later, ε gradually tends to zero with evolutionary generation increasing, there will be more Pareto optimal solutions in the population. Therefore, much more attention is paid to convergence than diversity. We can see that ε-truncation strategy is self-adapting and owns the strong ability of balancing diversity and convergence.

 The comparison between ε-truncation and ε-constraint is described as follow. Similarities: Both of them can effectively coordinate feasible solutions and infeasible solutions by permitting some infeasible solutions to participate in the evolution, and then they can improve the diversity; Differences: 1) It’s convenient for adaptive ε-truncation to integrate other tactics such as niche and crowding density estimation, which help to improve distribution of the evolution population. 2) The ε-constraint is based on one-to-one comparison rule. Although it accelerates the convergence rate, it causes to lose some individuals located in sparse region, and hampers the population diversity. Nevertheless, adaptive ε-truncation bases on set selection operation which can maintain good diversity by introducing improved crowding density estimation.

 

3.2 Evolutionary Process

 The key to evolutionary process lies in balance of exploration and exploitation. However, generally these two are conflicting. Therefore, it not only needs to strengthen the exploration and improve the diversity, but also enhance the local exploitation for enhancing the search efficiency. To better coordinate exploration and exploitation, we adopt random selection of basis vectors and difference vectors to increase disturbance, so as to improve diversity. Besides, after the differential mutation and crossover [40], we introduce exponential variation [41] to enhance exploitation. The mutation operator and crossover operator are respectively shown as formula (4) and (5).

\(\boldsymbol{V}_{i}=\boldsymbol{X}_{r 1}+F_{1} \times\left(\boldsymbol{X}_{r 2}-\boldsymbol{X}_{r 3}\right)+F_{2} \times\left(\boldsymbol{X}_{r 4}-\boldsymbol{X}_{r 5}\right)\)       (4)

where F is a scaling factor, \(\boldsymbol{X}_{r 1}, \boldsymbol{X}_{r 2} \boldsymbol{X}_{r 3}, \boldsymbol{X}_{r 4} \text { and } \boldsymbol{X}_{r}\) are distinct vectors and r1, r2, r3, r4 and r5 are diverse integers chosen from the set {1,2,…N}.\(V_{i}=\left(V_{i, 1}, V_{i, 2}, \ldots, V_{i, n}\right)\) is a mutation vector.

\(\boldsymbol{U}_{i, j}=\left\{\begin{array}{l} \boldsymbol{V}_{i, j}, \quad { rand }(j) \leq C R \\ \boldsymbol{X}_{i, j}, \text { otherwise } \end{array}\right.\)       (5)

where CR is a single parameter within the interval (0, 1), which controls the fraction of vector components inherited from the mutation vector.\(V_{i}=\left(V_{i, 1}, V_{i, 2}, \ldots, V_{i, n}\right)\) is a crossover vector. V_i,j is the j^th dimension of V_i , and X_i,j is the j^th dimension of X_i.

 Some studies have shown that differential mutation and crossover possess good ability of diversity maintenance. However, when dealing with complex CMOPs, differential evolution strategy has the problems of low convergence accuracy and slow convergence rate. The reason is that exploitation ability of differential evolution is insufficient. For the purpose of improving local exploitation ability, exponential variation is used after differential mutation and crossover, as shown in the formula (6).

\(X_{i, j}^{*}=\left\{\begin{array}{l} \boldsymbol{U}_{i, j}+\left[1-(2-2 \times {rand})^{\frac{1}{\eta+1}}\right] \times\left(u_{j}-l_{j}\right), \quad {rand} 0<0.5 \\ \boldsymbol{U}_{i, j}+\left[(2 \times {rand})^{\frac{1}{\eta+1}}-1\right] \times\left(u_{j}-l_{j}\right), \text { otherwise } \end{array}\right.\)       (6)

Nevertheless, the jth dimension of   may be out of range [lj,uj], so a repairing operation is required, as shown in the formula (7).

\(\boldsymbol{X}_{i, j}^{\prime}=\left\{\begin{array}{ll} l_{j}+r a n d 0 \times\left(u_{j}-l_{j}\right), & \text { if } \quad \boldsymbol{X}_{i, j}^{*}u_{j} \\ \boldsymbol{X}_{i, j}^{*}, \text { otherwise } \end{array}\right.\)       (7)

 

3.3 Elitist Selection Based ε-Truncation

 The N individuals are generated through evolutionary operation to form an offspring population (Ot). Then the offspring population and the parent population (Pt) are merged to produce a combined population (C_t) which comprises 2N individuals. Subsequently adaptive ε-truncation divides C_t into two populations (P_≤ε and P_>ε).

\(\underbrace{\left\{{X}_{1}, \ldots, X_{N}\right\}}_{P_t}+\underbrace{\left\{X_{1}^{\prime}, \ldots, X_{N}^{\prime}\right\}}_{O_t}=\underbrace{\left\{X_{1}, \ldots, X_{N}, X_{1}^{\prime}, \ldots, X_{N}^{\prime}\right\}}_{C_t}\)       (8)

\(\left\{\boldsymbol{X}_{1}, \ldots, \boldsymbol{X}_{N}, \boldsymbol{X}_{1}^{\prime}, \ldots, \boldsymbol{X}_{N}^{\prime}\right\}=\underbrace{\left\{\boldsymbol{X}_{1}, \ldots, \boldsymbol{X}_{N 1}\right\}}_{G(X_i, 1) \leq \varepsilon(t)}+\underbrace{\left\{\boldsymbol{X}_{N 1+1}, \ldots, \boldsymbol{X}_{2 N}\right\}}_{G(X_i )>\varepsilon(t)}\)       (9)

where N1 is the scale of the ε-truncation set (P_≤ε).

 So an important question to ask is ‘how to select N individuals as the next generation population (P _t+1)?’. It consists of three kinds of cases:

 Case one: When |P_≤ε|>N, we need to choose N individuals from P_≤ε as a next generation population.

 The individuals with the higher Pareto level and the greater crowding density should be chosen preferentially. First of all, the fast non-dominated sorting [42] is used to classify P_≤ε into F1,F2,…,Fk as shown in Fig. 1.

Fig. 1. Pareto level graph

 Next, the several levels F1,F2,...,Fs (1≤s≤k) are identified to make the total number (N_s) of individuals greater than or equal to N. If N_s is exactly equal to N, F₁,F₂,...,F_s are regarded as the next generation population. Otherwise, it needs to calculate the crowding density of all the individuals in F_s by formula (10). Afterward, N-(N_s-|F_s|) individuals that possess the most largest crowding density in Fs are combined with F₁,F₂,…,F_s-1 to become the next generation population. It is noteworthy that, the crowding density of all border solutions (block label in Fig. 1) in Fs is defined as infinity so as to ensure them to enter the next generation, and thereby the exploration range is broadened.

\(d_{i}=\frac{N-1}{\frac{1}{d_{i, 1}}+\ldots+\frac{1}{d_{i, j}}+\ldots+\frac{1}{d_{i, N-1}}}\)       (10)

where N is the scale of the evolution population, di,j is on behalf of Euclidean distance between individual X_i and X_j in objective space.

 Because the formula (10) has considered Euclidean distance between the individual X_i and all the other individuals in the population, the amount of calculation is heavy, and distant individuals or weak Pareto level individuals will produce negative influence. Therefore, it leads to the calculation deviation and reduces the diversity.

For this purpose, an improved crowding density formula is proposed as formula (11).

\(d_{i}=\frac{T 1}{\frac{1}{d_{i, 1}}+\ldots+\frac{1}{d_{i, j}}+\ldots+\frac{1}{d_{i, T 1}}}+\frac{T2}{\frac{1}{d_{i, 1}}+\ldots+\frac{1}{d_{i, j^{\prime}}}+\ldots+\frac{1}{d_{i, T 2}}}\)       (11)

where T1=T2=|F_s|^0.5/2, and d_i,j represents the j^th nearest Euclidean distance between X_i and the all individuals in the population, and di,j’ indicates the j^th nearest Euclidean distance between X_i and these individuals whose Pareto level is not inferior to X_i.

 It is clear that T1 and T2 are much smaller than N in formula (11), thereby computing cost is effectively decreased compared to formula (10). Furthermore, selecting the nearest T1 individuals will decrease influence of the distant individuals, and choosing the T2 individuals will eliminate effect of the individuals with weak Pareto level. Hence, the distribution of the individuals can be accurately reflected via the formula (11), and diversity will be kept well subsequently.

 Case two: When |P_≤ε|=N, P_≤ε is used as the next generation population.

 Case three: When |P_≤ε |<N and N-|P_≤ε| are regarded as the next generation population. This case indicates that the feasible region is likely to be relatively small, so only a small number of feasible solutions are generated in the evolution. Thus, it should put more emphasis on the search for the feasible solutions. Then the retention of infeasible solutions with the lower constraint violation can urge the evolution to close to the feasible region in a certain extent.

 

3.4 Process of the Proposed Algorithm

 For ease of understanding, this section provides specific steps of the ε-T-CMOA, and the flowchart is as shown in Algorithm 1.

Algorithm 1: Framework of the proposed ε-T-CMOPA

Input    Gmax: the maximum evolution generation    N: the size of population    F: the scaling factor    CR: the crossover probability    η: the exponential factor    pm: the variation factor Output    P: the final population found by ε-T-CMOPA 1: P0 ← Initialize Population 2: t=0 3: While t≤Gmax 4:    Ot ← Create offspring population by Evolution process 5:    Ct ← Pt∪Ot 6:    P≤ε and P>ε ← Execute ε-truncation 7:    Pt+1 ← Elitist selection based ε-truncation 8: End

 Step 1: Parameters Initialization, including the maximum generation (G_max), the size of population (N), the scaling factor (F), the crossover probability (CR), the exponential factor (η), and variation factor (p_m);

 Step 2: Generate the initial population P₀ that includes N individuals, the j_th dimension xj of each individual Xi=(x₁,x₂,...,x_n) is produced by \(x_{j}=l_{j}+r a n d 0 \times u_{j}, j=1,2, \ldots, n, i=1,2, \ldots, N, r a n d 0\)is a random number in the interval [0,1];

 Step 3: Execute the evolution operation, and generate the offspring population (O_t);

 Step 4: Combine the parent population (P_t) with Ot, and calculate the constraint violation and objective function value of the combined population (C_t);

 Step 5: Utilize adaptive ε-truncation to divide C_t into two sets (P_≤ε and P_>ε);

 Step 6: Implement elitist selection based ε-truncation and select N individuals as the next generation population (P_t+1);

 Step 7: Judge whether the G_max is satisfied? If so, output the Pareto optimal individual in the population; if not, go to step 3.

 

4. Experimental Classification Results and Analysis

4.1 Test Problems

 The PC Configuration are system: Windows XP_SP3; RAM: 2G; CPU: G620; CPU 2.60 GHz; Computer Language: MATLAB 2010. As a basis for the comparisons, some well-known test problems, CTP2-CTP7 [43] and CF1-CF5 [44], are involved.

 In CTP2-CTP7 series, parameter a controls distance from feasible region boundary to the PF; Parameter b controls the number of discrete region of the PF; Parameter c controls uniformity degree of the PF, if c=1, the PF is evenly distributed, else if c<1, the PF is close to the direction in f₁, else the PF moves toward to f₂; Parameter d controls length of the PF; Parameter e controls position of the feasible boundary; Parameter θ controls slope of the PF.

 In CF1-CF5 series, the number of decision variable is 10, and the number of objectives is 2. Thecorresponding PF lies in fi∈[0,1]. The test problems relate to the PF that has discrete points (CF1), convex (CF2), concave (CF3), continuous (CF4, CF5), disjointed (CF1, CF2, CF3), differently scaled (CF2, CF3, CF4, CF5).

 

4.2 Performance Metric

 In order to evaluate the performance of the related algorithms, the performance metrics are required. However, none of the existing performance metrics can authoritatively appraise both the diversity and the convergence of the obtained non-dominated solutions. Therefore, one must exercise various metrics in the performance evaluation of algorithms. In this paper, we have chosen the generational distance (GD) [34] which is widely used as an indicator for assessing the convergence, and the spacing (SP) [34] which reflects the diversity on CTP-series test problems. Simultaneously, the inverse generational distance (IGD) [44] which indicates a combined measure of both diversity and convergence are employed on the CF-series test problems.

\(\mathrm{GD}=\frac{1}{N}\left(\sum_{i=1}^{n} d_{i}^{p}\right)^{1 / p}\)       (12)

where N is the number of the obtained solution set \(\left(P_{\text {konw }}\right), p=2, d_{i}, i=1,2, \ldots\) indicates the nearest Euclidean distance between Xi in P_konw and the all solutions in PF. It is needless to say that the smaller GD denotes better convergence.

\(\mathrm{SP}=\sqrt{\frac{1}{N} \sum_{i=1}^{N}\left(\bar{d}-d_{i}^{\prime}\right)^{2}}\)       (13)

where N is the number of the obtained solution set, \(d_{i}^{\prime}=\operatorname{Min}_{j}\left(\sum_{k=1}^{m}\left|f_{k}^{i}(x)-f_{k}^{j}(x)\right|\right), i, j=1,2, \ldots, N_{\mathrm{PF}} \wedge i \neq j, \bar{d}=\sum_{i=1}^{N} d_{i}\) Obviously, the SP is smaller, the distribution is better.

 Supposing that P is a set of uniformly distributed points in the PF, and S is the obtained solution set, IGD metric is defined as follows:

\(\mathrm{IGD}=\frac{1}{N_{\mathrm{PF}}}\left(\sum_{i=1}^{n} d\left(v_{i}, S\right)\right)\)       (14)

where N_PF is the number of the points in P, and d(v_i,S) is the minimum Euclidean distance between v_i and the solutions in S.

 Since PF of CTP3 and CTP4 are composed of multiple discrete points, and PF of CTP5 consists of a continuous region and a plurality of discrete points, so it is inappropriate for using the SP to measure distribution. Thus, CTP2-CTP7 series test problems are divided into two groups in this paper: CTP2, CTP6 and CTP7 are taken as group 1, CTP3, CTP4 and CTP5 are seen as group 2. The number of discrete points and the GD are respectively employed to measure distribution and convergence in group 1, and similarly the SP and the GD in group 2. In CF1-CF5 series test problems, we will use IGD statistics for comparing results.

 

4.3 Experimental Settings

 All the experimental settings in the proposed algorithm are listed as the following.

 1) Number of runs: Each algorithm is run 30 times independently for each test instance.

 2) Termination criterion: The termination condition of an algorithm is specified in the form of the maximum number of generations (G_max). We use the same G_max for 3000 in all problems.

 3) Population size: The setting of population size N is 100 in all problems.

 4) Parameters for evolution process: F=0.5, CR=0.9, p_m=1/n, η=20.

 5) Significance test: To test differences in statistical significance, the Conover-Inman procedure is applied with 5% significance level in the Kruskal-Wallis test [45] for all pairwise comparisons.

 

4.4 Parameter Analysis

 The key parameters F and CR have a direct impact on evolution. Thus, this section discusses how F and CR affect the results of SP and GD on the test problem CTP2. All experiments run independently for 10 times and contrast experiment results have shown in Fig. 2- Fig. 5.

 F controls search scope, so it modulates the exploration ability of population. The lager F value introduces better diversity, but it will hamper the convergence. In contrast, although the smaller F value speeds up the convergence, it weakens the diversity and the algorithm may fall into local optima. When F ranges from 0.4 to 0.6 in Fig. 2, the distribution of obtained solution set is the best. Beside, as F equals to 0.5, ε-T-CMOA provides with the optimal stability. In Fig. 3, when F is between 0.2 and 0.5, the convergence of the obtained non-domination solutions looks better. However, the performance deteriorates with increasing F.

 CR controls the weight of variation individual in the crossover individual. As CR becomes smaller, the proportion of variation individual becomes fewer. Correspondingly diversity will be decreased. However, bigger CR slows down the convergence. In Fig. 4, the distribution of non-domination solutions get better when CR is from 0.6 to 1.0. However as CR becomes smaller, the distribution tends poorer. In Fig. 5, when CR is at interval [0.5, 1.0], the convergence performs well. Therefore, the bigger CR contributes to improve both diversity and convergence.    

Fig. 2. Different F for SP when CR=0.9 on CTP2

Fig. 3. Different F for GD when CR=0.9 on CTP2

Fig. 4. Different CR for SP when F=0.5 on CTP2    

Fig. 5. Different CR for GD when F=0.5 on CTP2

 

4.5 Comparative Experiment on CTP2-CTP7

 The mean and standard deviation of the SP, the number of discrete points and GD indicator value over 30 independent runs are used in our experiment. In addition, to test the differences in statistical significance, the performance score [46] is adopted to reveal how many other algorithms significantly outperform the corresponding algorithm on the considered test problem, and the smaller score denotes the better algorithm.

 To verify ε-T-CMOA, the following six state-of-the-art algorithms are considered as the peer algorithms: 1) NSGA-II with constrained multi-objective optimization algorithm [42], is denoted as NSGA-II; 2) BB-MOPSO [35]. It makes use of dual populations storage as the constraint handling technique and PSO as the evolutionary process; 3) Literature [26]. The DE and dual population storage are combined to solve CMOP; 4) Literature [34]. The SBX crossover and mutation are applied to generate new solutions, and adaptive penalty function method is used to handle the constraints. 5) Literature [47]. The new objective values are served to select individuals to the next generation. 6) Literature [48]. A hybrid genetic algorithm is suggested based on the boundary simulation method and trie-tree data structure.

 Table 1 shows both the mean and the standard deviation (St.de) of SP and GD. It can be observed that, the mean of the SP of ε-T-CMOA are better than the other six kinds of algorithms on the test problems CTP2 and CTP6. They indicate that ε-T-CMOA keeps the better population diversity in the evolutionary process, and gains the well-distributed solutions. For test instance CTP7, ε-T-CMOA displays the optimal performance in the mean of GD, which discloses that the obtained solution set of ε-T-CMOA is nearest to PF. However, Literature [47] performs best in terms of mean of SP on CTP7, and it is slightly better than ε-T-CMOA. Literature [34] algorithm achieves the best mean of GD on the test problems CTP2, and Literature [48] gets the smallest mean of GD on the test problems CTP6. Finally, ε-T-CMOA has found the smallest St.de of the SP and GD in the group 1 that indicates the optimal stability of our algorithm.

Table 1. Performance metrics on the first group of test problems

 Table 2 presents the mean and the St.de of the number of discrete points and GD in the group 2. It is clear that literature [26], BB-MOPSO and ε-T-CMOA have consistently found all 14 discrete points in all 30 independent runs on the test problem CTP3, yet GD mean of ε-T-CMOA are much better than these algorithms. Literature [34] gets the good mean of GD on the test problem CTP3, but the mean of discrete points is just 9.90. Because of loss of diversity maintenance in Literature [34] algorithm, some search area has not been fully explored, as a result, it misses a portion of Pareto solutions in evolution. Literature [47] and Literature [48] obtain the same excellent GD as our algorithm on CTP3, while their SP are relatively weak. On the test problem CTP4, both the number of discrete points and GD of ε-T-CMOA are superior to other algorithms. For the test problem CTP5, ε-T-CMOA has found all 15 discrete points in 30 independent runs, which is significantly better than the other algorithms. Literature [47] acquires the best convergence precision and Literature [48] gets second, but their number of discrete points are poor. It reveals that the diversity maintenance in Literature [47] and Literature [48] is not enough.

 As can be seen from Table 1 and Table 2, ε-T-CMOA performs the best GD on CTP3, CTP4 and CTP7. These results empirically illustrate that the final obtained non-dominated solutions by ε-T-CMOA for these test problems approximate the PF very well. Meanwhile, ε-T-CMOA offers the best distribution on CTP2, CTP3, CTP4, CTP5 and CTP6 except CTP7. These results demonstrate that ε-T-CMOA has the good ability of maintaining diversity. Since CTP3, CTP4 and CTP7 consist of discrete solutions or disconnected regions, in order to approach their PF much diversity of population is required right from the beginning of the algorithmic runs. For this purpose, ε-T-CMOA allows excellent infeasible solution to join in the evolution, so search scope is enlarged. Besides, the feasible Pareto solutions provide the power to move forward the PF. Therefore, the balance between diversity and convergence could be one of the reasons for better performance in ε-T-CMOA.

Table 2. Performance metrics on the second group of test problems

Fig. 6. Average score over CTP2-CTP7 test problems

Fig. 7. Ranking in the overall performance score for seven algorithms

 In order to better visualize the performance score in Table 1 and Table 2, Fig. 6 displays the performance over CTP2-CTP7. For the all test problems, ε-T-CMOA obtains good score. Finally, to quantify the overall performance of each algorithm on the all test problems, Fig. 7 shows the average performance score of all six test problems for seven algorithms. According to Fig. 7, we conclude that ε-T-CMOA could get the optimal results.

 Fig. 8-Fig. 13 have shown the Pareto solutions obtained by ε-T-CMOA on the test problems CTP2-CTP7. Among them, the gray area represents the feasible region, and the asterisk (electronic version in blue) represents Pareto optimal solutions.

 From Fig. 8 we can see that the PF of the test problem CTP2 consists of 13 disconnected areas. The non-dominated solutions found by ε-T-CMOA converge well to the PF and gain the good distribution. The PF of the test problem CTP3 comprises 14 discrete points as shown in Fig. 9. All the algorithms have found the optimal solutions in all 30 independent runs except NSGA-II, Literature [47] and Literature [48], which indicates a need to strength the diversity for these algorithms. In Fig. 10, the PF of the test problem CTP4 is composed by 14 discrete points, which are located in boundary of feasible region with longer tip compared to the test problem CTP3. Thus, it is more difficult to be solved. The ε-T-CMOA has found all the discrete points in 29 independent runs. The Pf of the test problem CTP5 includes a continuous region and 15 discrete points (see Fig. 11). The ε-T-CMOA consistently obtains all the discrete points in 30 independent runs. Fig. 12 shows that the feasible region of the test problem CTP6 involves a series of discrete strip-shaped region and the PF of CTP6 locates at the lowest position of the feasible region. Hence, the search must pass through different feasible region and eventually approach to PF, so it is easy to fall into local convergence. However, all algorithms have converged to the PF with good distribution. As can be seen from Fig. 13, the PF of the test problem CTP7 is formed by five regional components and a discrete point. And all the algorithms have made good distribution and convergence.

Fig. 8. Obtained solutions by ε-T-CMOA on CTP2

Fig. 9. Obtained solutions by ε-T-CMOA on CTP3

Fig. 10. Obtained solutions by ε-T-CMOA on CTP4 

Fig. 11. Obtained solutions by ε-T-CMOA on CTP5

Fig. 12. Obtained solutions by ε-T-CMOA on CTP6

Fig. 13. Obtained solutions by ε-T-CMOA on CTP7

 

4.6 Comparative Experiment on CF1-CF5

 In this section, we have carried out further experiments between ε-T-CMOA and two good performers (LIULI [49] and MTS [50] on CEC 2009) on CF1-CF5. Table 3 compares the IGD statistics (all 30 independent runs) acquired from LIULI, MTS and ε-T-CMOA.

Table 3. The comparison of IGD for three algorithms

 From Table 3, we can see that ε-T-CMOA has found substantially better mean and St.de of IGD than LIULI on the test problems CF1, CF3 and CF5. However, LIULI accomplishes better performance on the test problems CF2 and CF4. The ε-T-CMOA obtains the better results in terms of minimum and mean of IGD values than MTS on the test problems CF1-CF3, and the case is vice versa on the test problems CF4 and CF5. LIULI and MTS are essentially two advanced algorithms based on multi sub-populations. They evolve these sub-problems simultaneously by uniform weight vectors methods, which can maintain the diversity of population and approach the evenly distributed Pareto optimal solutions. Therefore, it may be the internal causes that LIULI and MTS gains slightly better results on continuous functions CF4 and CF5. However, the obtained St.de values by ε-T-CMOA on CF2, CF4 and CF5 have clear advantage, which indicates the optimal stability of ε-T-CMOA.

 Fig. 14-Fig. 18 present the non-dominated solutions obtained by ε-T-CMOA on the test problems CF1–CF5. These solutions are selected from the final population of the run with the best IGD among the 30 independent runs. It is obvious that ε-T-CMOA has gotten good convergence on the test problems CF1 and CF2 from Fig. 14 and Fig. 15.

Fig. 14. Obtained solutions by ε-T-CMOA for CF1

Fig. 15. Obtained solutions by ε-T-CMOA for CF2

Fig. 16. Obtained solutions by ε-T-CMOA for CF3    

Fig. 17. Obtained solutions by ε-T-CMOA for CF4

Fig. 18. Obtained solutions by ε-T-CMOA for CF5

 In Fig. 16, the PF of the test problem CF3 is discontinuous and concave, so it could be harder than all other test problems to solve. Three algorithms do not completely converge to the PF. Thus, these algorithms need to further strengthen the exploration and exploitation capability. In Fig. 17 and Fig. 18, ε-T-CMOA acts relatively poorly similar to LIULI and MTS on test problems CF4 and CF5, because non-dominated solutions are not distributed evenly along the entire PF.

 In summary, compared with the test problems CTP2-CTP7, much more efforts should be paid to the test problems CF1–CF5.

 

4.7 Validity verification of improved crowding density estimation

 To verify the effectiveness of the improved crowding density estimation, the original crowding density estimation (Equation 10) and improved crowding density estimation (Equation 11) are compared based the framework of ε-T-CMOA on the test problems CTP2-CTP7. They are abbreviated as original and improvement respectively. The contrastive experimental results are presented in Table 4-Table 6.

 Table 4 clearly shows that the SP of improvement is better than original on the test functions CTP2, CTP6 and CTP7. Especially on the test problem CTP7, distribution has been improved significantly. Thus, the improved crowding density estimation really performs well on diversity maintenance. Table 5 displays that the number of discrete points is almost the same between original and improvement. Both of them have found all the discrete points in 30 times independent runs. In Table 6, all the mean GD values obtained by improvement are smaller than original, particularly on the test problems CTP2 and CTP7. Finally the St.de values found by Equation 11 are generally better, which demonstrates enhancement of the algorithm stability. So the improved crowing density estimation plays an important role in improving diversity in ε-T-CMOA. Besides, the good diversity ensure eventually to come close to PF.

Table 4. The comparison of SP between original and improvement

Table 5. The number of discrete points between original and improvement

Table 6. The comparison of GD between original and improvement

 

4.8 Validity verification of exponential variation

 To verify the effectiveness of the exponential variation, we carry out the comparative experiments between removal of exponential variation of ε-T-CMOA and ε-T-CMOA on the test problems CTP2-CTP7. They are abbreviated as original and improvement respectively. The contrastive results are illustrated in the Table 7-9.

 From Table 7, we can find that the SP values of improvement perform better than original on the test functions CTP2, CTP6 and CTP7. And the distribution on CTP7 has been well increased. In Table 8, the number of discrete points gained by both original and improvement is almost the same, which demonstrates no significant interaction of exponential variation on CTP3, CTP4 and CTP5. Table 9 shows that all the mean GD values obtained by improvement get better, except on CTP7. Also the results indicate that the exponential variation as a local search strategy can strength the exploitation ability, which ensures the better convergence.

Table 7. The comparison of SP between original and improvement

Table 8. The number of discrete points between original and improvement

Table 9. The comparison of GD between original and improvement

5. Conclusion and future work

 In this paper, we have suggested a constrained multi-objective optimization algorithm, called ε-T-CMOA. From the perspective of the complementary advantages of dual population storage and ε-constraint, a novel constraint handling technique, ε-truncation, is presented to coordinate diversity and convergence by adaptively exploiting feasible Pareto solutions and infeasible solutions with both the lower constraint violation and the better objective value. Besides, the improved crowding density can evaluate the distribution accurately, and cut down the computation load via selecting a part of Pareto solutions and near solutions to participate in the calculation. Moreover, for achieving a better compromise between global exploration and local exploitation, exponential variation is introduced following crossover operation and mutation operation.

 To validate the competitiveness, we have conducted a comprehensive experimental comparison with eight state-of-the-art algorithms that subject to different kinds of technologies. Many benchmark test problems (CTP2-CTP7 series and CF1-CF5 series) are selected to challenge different capabilities of the algorithms. In CTP-series test problems, the proposed ε-T-CMOA has been able to successfully find a well-converged and well-diversified set of solutions over multiple independent runs. However, on some CF-series test problems, like LIULI and MTS, ε-T-CMOA encounters with an increasingly difficult task of maintaining diversity and converging to the PF. Although different algorithms have exhibited their working on different problems, we have observed that none of these algorithms is capable of treating all the test problems, which implies that a careful choice of algorithms is still required at present when we deal with a complicated CMOP.

 In the future, it will be attractive to extend ε-T-CMOA to handle constrained many-objective problems by incorporating advanced EAs. Moreover, we will apply ε-T-CMOA to real-world problems for further confirming its effectiveness.