A constrained multi-objective optimization algorithm using an efficient global diversity strategy

When solving constrained multi-objective optimization problems (CMOPs), multiple conflicting objectives and multiple constraints need to be considered simultaneously, which are challenging to handle. Although some recent constrained multi-objective evolutionary algorithms (CMOEAs) have been developed to solve CMOPs and have worked well on most CMOPs. However, for CMOPs with small feasible regions and complex constraints, the performance of most algorithms needs to be further improved, especially when the feasible region is composed of multiple disjoint parts or the search space is narrow. To address this issue, an efficient global diversity CMOEA (EGDCMO) is proposed in this paper to solve CMOPs, where a certain number of infeasible solutions with well-distributed feature are maintained in the evolutionary process. To this end, a set of weight vectors are used to specify several subregions in the objective space, and infeasible solutions are selected from each subregion. Furthermore, a new fitness function is used in this proposed algorithm to evaluate infeasible solutions, which can balance the importance of constraints and objectives. In addition, the infeasible solutions are ranked higher than the feasible solutions to focus on the search in the undeveloped areas for better diversity. After the comparison tests on three benchmark cases and an actual engineering application, EGDCMO has more impressive performance compared with other constrained evolutionary multi-objective algorithms.


Introduction
Constrained multi-objective optimization problems (CMOPs) widely exist in many practical applications [1,2], which need to optimize multiple conflicting objectives and meet one or more constraints simultaneously. The maximum problem can usually be transformed into the minimum problem. Without loss of generality, a CMOP can be mathematically formulated as follows. where X (x 1 , x 2 , . . . , x n ) is a candidate solution consisting of n decision variables, ⊆ R n is the decision space, F : → R m is the objective vector consisting of m conflicting objectives,g i (X ) ≤ 0 is the ith inequality constraint, q is the number of inequality constraints, h j (X ) 0 is the jth equality constraint and h is the number of equality constraints. The constraint violation of a solution x is calculated by A solution x is feasible if CV(x) 0, otherwise, x is called an infeasible solution. Given two feasible solutions x 1 , x 2 ∈ , x 1 is said to dominate x 2 (denoted as x 1 ≺ x 2 ), if the following conditions are met, f i (x 1 ) ≤ f i (x 2 ) for every i ∈ {1, . . . , m}, f i (x 1 ) ≤ f i (x 2 ) for at least one j ∈ {1, . . . , m}. A solution x * is called Pareto optimal if there is no solution in that Pareto dominates it. The set of all Pareto optimal solutions is called the Pareto set (PS), and PF {F(x)|x ∈ PS}, namely the Pareto set in the objective space is called the Pareto front (PF).
Evolutionary algorithms (EAs) do not have any mechanisms to handle constraints while most real-world design optimization problems often encounter many constraints. A large number of constraint handling techniques have been proposed in the field of evolutionary computing [3]. The existing constraint-handling techniques (CHTs) in constrained multi-objective optimization can be divided into six groups, including penalty functions [4], constrained dominance principle (CDP) [5], ε constrained method [5], stochastic ranking [7], multi-objective optimization-based methods [8], and hybrid methods [9].
In recent years, some constrained multi-objective evolutionary algorithms (CMOEAs) have also been developed to solve CMOPs. Fan et al. [10] proposed a biphasic CMOEA, namely push and pull search (PPS), which divides the search process into two different stages: push and pull search stages. In the push stage, a multi-objective evolutionary algorithm is used to explore the search space without considering any constraints. In the pull stage, the constraints are considered to help the population pull back to the feasible and non-dominated regions. Liu et al. [11] proposed a two-stage framework called ToP. In ToP, the first phase, the original CMOP is transformed into a constrained single-objective optimization problem. In the second phase, all the constraints and objectives are considered to obtain the Pareto optimal solutions. Zhu et al. [12] developed a decomposition-based constrained MOEA called MOEA/D-DAE. In MOEA/DDAE, a detect-and-escape strategy was proposed, which uses the feasible ratio and the change rate of overall constraint violation to detect stagnation, and adjusts the weight of the constraint violation for guiding the search to escape from stagnation states. Ma et al. [13] proposed a multi-stage CMOEA called MSCMO. In MSCMO, in the early stages, only a small number of constraints are considered, which can make the population efficiently converge to the potential feasible region with good diversity. As the algorithm moves to the later stages, more constraints are considered to search for the optimal solutions based on the solutions obtained in the previous stages.
Although constrained multi-objective optimization has been developed for two decades, there exist some limitations in solving problems with small feasible regions and complex nonlinear constraints, which may lead to the population being unable to jump out of the local optimal or local feasible region, and may get trapped into a narrow feasible region, especially when feasible region consists of several disjoint parts and narrow in the search space. Some recent studies do suggest that reasonably utilizing the information on objectives is beneficial to balance the importance of constraints and objectives in the evolutionary process [14,15], and the information of infeasible solutions can help to resolve the dilemma between exploration versus exploitation [16]. Following this idea, this paper introduces an evolutionary algorithm that uses an efficient global diversity strategy to maintain a certain number of infeasible solutions with welldistributed feature to solve CMOPs. The main contributions of this work are summarized as follows: 1. An efficient global diversity CMOEA (EGDCMO) is proposed to solve CMOPs with small feasible regions and complex constraints. In this algorithm, an efficient global diversity strategy is proposed to select infeasible solutions for the next generation, which helps to maintain infeasible solutions with well-distributed feature at the early stage of evolution. In this way, the population can effectively converge to the potentially feasible region with global diversity.

A new fitness function is used in the proposed EGDCMO
to evaluate infeasible solutions, which can balance the importance of constraints and objectives. 3. The proposed EGDCMO is tested on benchmark CMOPs and a real-world application to verify its effectiveness. According to the experimental results, the proposed EGDCMO outperforms some state-of-the-art CMOEAs on both the benchmark CMOPs and the engineering problem.
The remainder of this paper is organized as follows. In "Related works", we first introduce the existing evolutionary approaches developed for CMOPs. "The proposed algorithm" explains the proposed algorithm in detail. The experiments of the benchmark CMOPs and a real-world application are carried out in "Experimental study". "Conclusions" makes a summary.

Related works
In this section, we briefly review the existing multi-objective evolutionary algorithms (MOEAs) to solve CMOPs. The existing MOEAs that maintain a certain proportion of infeasible solutions in the population during evolution are also reviewed.

Existing MOEAs
Generally speaking, existing MOEAs can be divided into the following six categories.
The first category is the penalty function method-based MOEAs [4]. The penalty function method was first used to deal with single objective constrained optimization problems. Its basic principle is to add a penalty function to the original objective function to get an augmented objective function and convert a constrained optimization problem into an unconstrained one [17]. Because of its simple form, it has been widely used in EAs. However, penalty functions require careful tuning of penalty factors to properly punish the infeasible solution and penalty factors always depend on the problem [18]. Researchers have developed some types of penalty functions for EAs, such as static penalty [19], dynamic penalty [20], adaptive penalty [21], and so on. In addition, the penalty function is also used in MOEAs. Woldesenbet et al. [22] proposed a constraint handling technique for multi-objective evolutionary algorithms based on an adaptive penalty function and a distance measure. Jan et al. [23] introduced a penalty function in the MOEA/D-DE updating scheme, which penalizes infeasible solutions based on an adaptive threshold.
The second category is dominance-based approaches. Dominance-based approaches are mainly driven by feasibility information where feasible solutions are always better than infeasible solutions to survive to the next generation. CDP is obtained by extending the feasibility rule, which was proposed by Deb et al. [5] in NSGA-II, and it is one of the most popular and effective CHTs. In this method, given two solutions, x 1 is said to dominate x 2 , if one of the following conditions is met.
Oyama et al. [24] developed a modified dominance relation according to which solutions that violate a fewer number of constraints are preferred. Ma et al. [14] proposed a new fitness function with two rankings, its first ranking is the solution's ranking based on CDP, and the other is the solution's ranking based on Pareto dominance.
The third category is the ε constrained method. The idea was originally proposed by Takahama and Sakai [5] in 2006. In this method, for a solution whose constraint violation degree is smaller than ε, the solution is regarded as feasible, otherwise, the solution is regarded as infeasible. The key to the ε constraint method also requires careful tuning of parameters to properly control the value of ε. Yang et al. [25] used the epsilon method to handle constraints and applied it to multi-objective evolutionary algorithm based on decomposition (MOEA/D). Fan et al. [26] proposed an improved epsilon constraint-handling mechanism and combined it with a decomposition-based multi-objective evolutionary algorithm, which adjusts the value of ε dynamically according to the ratio of feasible to total solutions in the current population.
The fourth category is the stochastic sorting (SR) approach. The idea of SR to deal with constraints is to introduce a random probability P f to select the solutions with good fitness or low CV. In other words, given any pair of adjacent solutions, if both solutions are feasible, then the probability of comparing them (to determine which one is more suitable) according to the objective function is 1, otherwise, the probability of comparing solutions according to CV is P f . Geng et al. [27] combined stochastic ranking with an infeasible elitist-based MOEA. Zhang et al. [28] proposed the dynamic stochastic selection, in which the probability parameter P f in the stochastic ranking method could change dynamically during the evolution, and applied it to the framework of multimember differential evolution.
The fifth category is the multi-objective optimization approach. The main idea of multi-objective constraint handling is to treat the constraints as objectives, that is, a CMOP is transformed into an unconstrained multi-objective problem, and then the transformed problem is solved by the multi-objective optimization technique. Long et al. [29] designed a multi-objective optimization-based new constraint handling technique, by taking closeness, diversity, and feasibility as three objectives in a multi-objective subproblem. Zhou et al. [30] proposed a tri-goal evolution framework and designed two indicators for convergence and diversity, converting the constraints into the third indicator for feasibility to solve constrained many-objective problems.
The last category is the hybrid method approach. The hybrid method is the combination of two or more different CHTs to deal with constraints. Deb and Datta [31] proposed a CHT method that combines multi-objective optimization and penalty function methods to deal with constrained optimization problems. Qu et al. [9] employed an ensemble of three different constraint handling methods, including adaptive penalty function, the superiority of feasible solution, and the ε constrained method. In this method, the population is divided into three subpopulations, each of which adopts different CHTs. Wang and Cai [32] proposed an adaptive tradeoff model (ATM), which divides the evolutionary process of the population into three stages and adopts different strategies according to the characteristic of each stage.

Existing MOEAs with maintaining infeasible solutions
The information on infeasible solutions helps to explore undeveloped regions [12,33], and to solve many challenging problems, including problems consisting of multiple disjoint feasible regions, or a huge infeasible barrier, and so on. Most evolutionary optimization algorithms focus on finding the optimal feasible solutions of the problem. Therefore, in the process of evolution, feasible solutions are always considered to be better than infeasible solutions.
The method of preserving the infeasible solution is adopted for the first time in solving constrained singleobjective optimization problems. Mezura-Montes and Coello [34] suggested a simple multimembered evolutionary strategy (SMEs) that allows the best infeasible solution to remain in the population for the next generation. Ray et al. [35] proposed an infeasibility-driven evolutionary algorithm to solve CMOPs which maintains a small percentage of infeasible solutions close to the constraint boundaries during its course of evolution. Isaacs et al. [36] proposed a constraint handling method that maintains infeasible solutions in the population by reformulating the original m objective constrained minimization problem into an unconstrained m + 1 objective minimization problem, and an additional objective function is the number of constraint violations.
In recent years, some coevolutionary algorithms are proposed to solve CMOPs. In coevolutionary constraint handling techniques, there is usually an archive to save infeasible solutions. Tian et al. [33] proposed an evolutionary algorithm evolving two populations to solve CMOPs, in which a population evolves without considering constraints to maintain diversity, and the other population to solve the original CMOP to maintain feasibility as well as convergence. Li et al. [16] proposed a two-archive EA (C-TAEA) for constrained multi-objective optimization. In C-TAEA, the convergenceoriented archive (CA) is evolved to optimize the constraints and objectives and aims to push the population toward the Pareto front. The diversity-oriented archive (DA) is evolved to optimize only the objectives and aims to maintain the population diversity by exploring areas underexploited by the CA, including those infeasible regions. Sorkhabi et al. [37] proposed an efficient approach for constraint handling in multi-objective particle swarm optimization, which divides particles population into two non-overlapping populations, one is feasible population, the other is infeasible population, where feasible particles are evolved toward Pareto optimality, and infeasible particles are evolved toward feasible particles. Wang et al. [38] proposed a cooperative differential evolution framework (CCMODE) for constrained multi-objective optimization, including M subpopulations for constrained single-objective optimization and an archive population for constrained M-objective optimization.
Very recently, some CMOEAs using valuable information existing in objective functions are proposed, which a certain number of infeasible solutions are preserved in the population to enhance the convergence pressure at the early stage. Yu et al. [15] proposed a novel CMOEA named DSPCMDE, in which some valuable infeasible solutions are maintained, and the balance in objective functions and constraints is adjusted dynamically by considering the desired searching emphasis in different evolutionary stages. Tian et al. [39] proposed a two-stage CMOEA to solve CMOPs with different feasible regions. In one stage, some infeasible solutions are maintained to help the population cross infeasible regions, and the fitness evaluation strategies are adjusted during the evolutionary process to adaptively balance objective optimization and constraint satisfaction. Peng et al. [40] proposed a novel constraint-handling technique based on directed weights to deal with CMOPs, which utilizes the useful information contained in infeasible individuals to guide the search to the promising region.

The proposed algorithm
In this section, the proposed EGDCMO will be explained in detail. The general flow chart of this proposed EGDCMO is given in Fig. 1. The proposed EGDCMO starts with the random initialization of a population P with size N, and then the algorithm evolves from generation to generation until the maximal number of function evaluations (maxNFE) is reached. In each generation, the algorithm repeats the main following operations: (1) reproduction, (2) external archive update, (3) environmental selection of feasible solutions, and (4) an efficient global diversity strategy is used to select infeasible solutions for the next generation. It is worth noting that, EGDCMO maintains a certain number of infeasible solutions with diversity, which provides as much diversified information as possible. The specific demonstrations are summarized in the following contents.

Procedure of EGDCMO
Algorithm1 provides the pseudo code of the proposed EGD-CMO. In Line 1, the proposed EGDCMO starts with the random initialization of a population P with size N. The main loop of the proposed EGDCMO is shown in Lines 4-33. Especially, for each generation, the reproduction is shown in Lines 5-6, which uses the binary tournament selection and the polynomial mutation to generate offspring. The external archive update is shown in Line 7, where the non-dominated feasible solutions are saved in the external archive. In Line 8, the population is divided into an infeasible set and a feasible set. Line 9 shows the number of infeasible solutions that need to be maintained. If the infeasible set has no more than N inf solutions, all solutions from the infeasible set are selected. If the feasible set has no more than N f solutions, all the feasible solutions are selected and the rest are filled with infeasible solutions, which are formulated as follows.
where N 1 and N 2 are the number of infeasible solutions and feasible solutions respectively, α is the proportion of infeasible solutions to maintain, and N is the population size. We set β 0.8, and why we set it to 0.8 will be discussed in "Analysis of the proposed EGDCMO". The uniformly distributed weight vectors are generated as shown in Line 10, and the number of evenly distributed weight vectors is equal to N inf . Then, the fitness and the truncation method are used to select N f solutions from S f in the environment selection shown in Line 11. The determination of whether the value of the flag is 1 is shown in Lines 12-19. Specifically, when the average change of feasible population (S f ) within 150 gener-ations on each objective is less than the given threshold value of 0.01, the value of the flag is 1. The flag is judged again after 150 generations. For the average value of S f on the ith objective, we used Formula (4) to normalize the objective value of each feasible solution in S f . Then, we used Formula (5) to calculate the average value of feasible solutions in S f on each objective [13], and |S f | is the number of solutions in feasible population S f : In Lines 20-21, the N inf infeasible solutions are selected from S inf . Afterward, the fitness value of each solution in S inf is calculated. In Lines 22-23, feasible and infeasible solutions are combined. To add selection pressure to generate better infeasible solutions [36,41], the infeasible solutions are ranked higher than the feasible solutions to focus on the search in the undeveloped areas for better diversity. When the evolution has entered the final stage, that is, NFE 0.8*maxNFE, and NFE is the number of function evaluations, the population is replaced by the solutions in the external archive obtained in the previous generation and make α 0. The loop executes until the maximal number of function evaluations is reached, and then returns S as the final result, as shown in line 32.

A new fitness function to evaluate infeasible solutions
An example ranking procedure using CDP is presented in Fig. 2. For the problem with strict constraints, all solutions in the population may infeasible in the initial stage. Then sort the solutions according to CV, that is, the smaller the CV, the higher the ranking, and select the best N solutions. Therefore, solutions with smaller CV can survive to the next generation, and make the population evolve to the feasible region promptly. Due to the lack of mechanisms to maintain population diversity, the population cannot approach the feasible region from different directions and may get trapped in a narrow feasible region. It is worth mentioning that retaining Fig. 2 Ranking procedure using CDP some "good" infeasible solutions in each generation helps to solve the problem where the optimal solutions lie on constraint boundaries.
To achieve the tradeoff between objective functions and constraints, and choose "good" infeasible solutions for the next generation, a natural way is to use the information of objective functions. Fortunately, we can use the method in SPEA2 [42] to calculate the fitness of infeasible solutions without considering any constraints, and the ranking of infeasible solution x i can be obtained (S F i ) through Lines 1 of Algorithm 2, a smaller value of S F i represents the better performance of x i . Afterward, the other ranking of infeasible solution x i can be obtained (S NSCV i ), in which the number of satisfied constraints (NS) and CV are used as two objectives (minimize CV, and maximize NS) for sorting based on Pareto dominance, as shown in Line 3-10 of Algorithm 2. Subsequently, to balance the importance of objectives and constraints, we use a new fitness function to evaluate infeasible solutions.
Note that, P f is the proportion of feasible solutions. When handling infeasible solutions, the information of objectives is utilized and NS is also considered, and the larger the NS value is, the more likely the infeasible solution is to be located at the constrained boundary, especially when the problem has many constraints. Compared with CDP, which only considers minimizing CV, the new fitness function can evaluate the infeasible solution more comprehensively.
It is necessary to point out that although the fitness function to evaluate infeasible solutions in this paper has some similarities to that of the ToR in [14]. The main differences are as follows. In ToR, the fitness function of each solution is defined as the weighted sum of two rankings: one is the solution's ranking based on CDP and the other is the solution's ranking based on Pareto dominance. While in our method, the fitness function of infeasible solutions is defined as the weighted sum of two rankings: one is the solution's ranking based on the fitness of infeasible solutions obtained by SPEA2 without considering any constraints and the other is a non-dominated sorting of these two aspects (minimize CV, maximize NS). In our method, when handling infeasible solutions, the information of objectives is utilized and NS is also considered, and the larger the NS value is, the more likely the infeasible solution is to be located at the constrained boundary [43], especially when the problem has many constraints. Compared with CDP, which only considers minimizing CV, the new fitness function can evaluate the infeasible solution more comprehensively.
In ToR the search is more biased toward constraints than objectives as the proportion of feasible solutions increases. In our method, when the proportion of feasible solutions increases, α 1 is greater than α 2 . which means that the search is more biased toward objectives than constraints. That is to say, with the increase of α 1 , less and less information about constraints is considered, which helps explore the undeveloped areas for better diversity. Considering that more information on objective functions is helpful to the exploration both in the feasible and infeasible regions, and avoid stagnation in some local regions [15]. The fitness of currently infeasible solutions is related to the proportion of feasible solutions, and when the proportion of feasible solutions increases, the selection of infeasible solutions tends to utilize more objective function information.

An efficient global diversity strategy to select infeasible solutions for the next generation
The key component in EGDCMO is how to select N inf infeasible solutions from S inf . Before explaining the efficient global diversity strategy to select infeasible solutions, we first introduce the idea from MOEA/D-M2M [44] to divide the objective space into several subregions, each of which is represented by a unique weight vector generated by the canonical simplex-lattice design method [45]. Specifically, W {w 1 , w 2 ,…,w Ninf } is a set of uniformly distributed weight vectors, and a subregion i is defined as where f (x), w is the acute angle between f (x) and w. After the setup of subregions, solutions are assigned to them as follows. First, objective values of all solutions at the current generation are translated, which is calculated as where f * i and f nad i are the estimated ideal and nadir points, which are the minimum and maximum values of objective values of f i at the current generation, respectively. Then, each solution X of the population is associated with a unique subregion, and the index of the subregion is determined as The efficient global diversity strategy to select infeasible solutions for the next generation has two characteristics: (1) in the process of evolution, the way of selecting infeasible solutions according to the value of the flag and (2) at the early stage of evolution, the diversity of the population will be maintained when most solutions are infeasible. The pseudo code of selecting infeasible solutions is shown in Algorithm 3. Specifically, the value of the flag is judged according to the method introduced in "Procedure of EGDCMO". If flag 1, that is, the change of the average value of the feasible population (S f ) on each objective is less than the given threshold value of 0.01 within 150 generations, which means the feasible population obtained in the process of evolution changes little in the objective value, and the feasible population may fall into a narrow feasible region. Then, the fitness and the truncation method are used to select N inf solutions from S inf without considering the constraints (lines 3 of Algorithm 3), which helps explore the undeveloped areas for better diversity and jump over the infeasible band to evolve toward the global PF. At the early stage of evolution, flag 0, the population may only contain infeasible solutions. Then, we separately associate each solution in S inf with its corresponding subregion according to the method mentioned above. The next step is to determine which solutions should be selected, we iteratively investigate each subregion and decide the survival of solutions in S inf to the next generation. In particular, in each iteration, only one solution can survive in each subregion. In other words, for i W , i ∈ {1, . . . , N inf }, the solutions in S inf associated with i W , denoted as X i , their fitness values will be calculated by Formula (6), and the solution with the smallest fitness value will be chosen to survive to the next generation (lines 6-18 of Algorithm 3). Here, the best solution x b is identified as This iterative investigation continues until the infeasible population S inf is filled.

Analysis of the proposed EGDCMO
To analyze the effectiveness and understand the mechanisms to maintain population diversity of the proposed EGDCMO, we presented a hypothetical population with ten candidate solutions for the MW11 [46] problem, and the best five solutions will be selected. As shown in Fig. 3a and Table 1, when all solutions are infeasible, if CDP is adopted, A, B, C, D, and F are selected, since they have a smaller CV. Solutions that perform well in the object space are eliminated, this will make the population unable to approach the feasible region from different directions and may get trapped into a narrow feasible region. From Fig. 3b, the objective space is divided into several subregions, note that, the number of subregions depends on the number of infeasible solutions that need to be maintained, and each solution is associated with a subregion. Afterward, we traverse each subregion and select a solution from each subregion until the quantity requirements are met.
If more than one solution has been associated with i C , the solution with the smallest fitness value will be chosen. For example, three infeasible solutions A, B, and C in 1 C , their fitness values are calculated by the formula (6), and solution A with the smallest fitness value is selected. It can be found that since formula (6) introduces the information of objectives, the selected solution A is located near the constraint PF. Figure 3 and Table 2 present an example to illustrate this issue. Finally, A, F, H, I, and J are selected based on the proposed EGDCMO, where H, I, and J have the opportunity Fig. 3 Hypothetical distributions of ten candidate solutions A-J on MW11, the black dots in a and b are the best five solutions selected by CDP and EGDCMO, respectively Table 1 The best five solutions selected by CDP for the next generation on MW11 are highlighted in gray Table 2 The best five solutions selected by EGDCMO for the next generation on MW11 are highlighted in gray to approach other feasible regions and maintain population diversity. Figure 4 depicts the hypothetical distributions of ten candidate solutions on 2-objective LIR-COMP10, which has a band of the infeasible region, posing challenges to MOEAs in terms of convergence. As shown in Fig. 4a and Table 3, if CDP is adopted, F, G, H, I, and J are selected. Since CDP prefers feasible solutions and ignores some infeasible solutions that perform well in the objective space, it has trouble in jumping over the infeasible band. The proposed EGDCMO uses an efficient global diversity strategy to maintain N inf infeasible solutions, where N inf is calculated by the formula (3), and α is set to 0.5. From Fig. 4b, the objective space is divided into several subregions, and the solution with the smallest fitness value in each subregion will be selected. For example, infeasible solutions B, and C in 2 C , if only the constraint value is considered to determine the survival of the solution, C with the smallest CV value is selected instead of solution B with better convergence. We use the method introduced in "A new fitness function to evaluate infeasible solutions" to evaluate infeasible solutions in each subregion, which helps balance the importance of constraints and objectives. In this way, solution B with better convergence is selected (Table 4).
Next, to show the superiority of the proposed EGDCMO dealing with CMOPs with complex constraints. The simulated binary crossover [47] and polynomial mutation [48] are adopted for generating offspring. The probability of simulated binary crossover is set to 1, and the probability of polynomial mutation is set to 1/D (D denotes the number of decision variables). The proposed EGDCMO is compared to NSGA-II for detailed discussion. The populations in the early, middle, and last of NSGA-II and EGDCMO are plotted in Fig. 5. The MW11 problem has a very small feasible region (the size of feasible region < 0.1‰). As shown in the first column of Fig. 5, the solutions of NSGA-II are infeasible in the Fig. 4 Hypothetical distributions of ten candidate solutions A-J on LIR-CMOP10, the black dots in a and b are the best five solutions selected by CDP and EGDCMO, respectively Table 3 The best five solutions selected by CDP on LIR-CMOP10 for the next generation are highlighted in gray Table 4 The best five solutions selected by EGDCMO on LIR-CMOP10 for the next generation are highlighted in gray early. Due to the lack of mechanisms to maintain population diversity, the population only converges to a single feasible region. In the middle generations, some "good" infeasible solutions appear. Since CDP prefers feasible solutions, these infeasible solutions cannot survive to the next generation, failing to explore new areas of the feasible region until the algorithm stops. When EGDCMO handles infeasible solutions, the objective space is divided into several subregions, and the information on objectives is utilized, which achieves the tradeoff between objectives and constraints in each subregion. Thus, the population can approach the feasible region from different directions. Some "good" infeasible solutions that balance the importance of constraints and objectives are retained, which is helpful to search the Pareto optimal solutions located on the boundary of the feasible region and explore the undeveloped areas for better diversity. Figure 6 shows the populations in the early, middle, and last generations of NSGA-II and EGDCMO on LIR-CMOP10 [26], which has a particularly complex constrained landscape and a huge infeasible barrier. As shown in the first column of Fig. 6, when the number of feasible solutions is greater than the population size, the infeasible solutions with better convergence are eliminated since CDP prefers feasible solutions, which is difficult to generate offspring that can jump out of the infeasible barrier. Finally, the solutions obtained by NSGA-II are completely stuck at the boundary of the infeasible barrier. EGDCMO explicitly maintains N inf infeasible solutions which with better convergence, and N inf is calculated by the formula (3). From the second column in Fig. 6, we can find that many infeasible solutions are located in the infeasible barrier, which helps to generate offspring To analyze the influence of parameters on the performance of the proposed EGDCMO, β is set to 0.1, 0.3, 0.5, 0.7, 0.8, and 0.9, and we performed a comparative study on the problems with a highly multimodal landscape and a band of the infeasible region, that is DC2-DTLZ3 and DC3-DTLZ3 problems. Figure 7 plots Pareto fronts with median IGD value among 30 runs obtained by EGDCMO on DC3-DTLZ3 corresponding to different values of β, it can be observed that when β is set to 0.8 and 0.9, EGDCMO can jump over the infeasible band and find feasible solutions with the optimal objective value on DC3-DTLZ3. It can be concluded that the larger the β value, the better EGDCMO can fully explore in underdeveloped areas and jump out of local optimal. Figure 8 plots the changes in the inverted generational distance (IGD) values of the proposed EGDCMO on DC2-DTLZ3 and DC3-DTLZ3 corresponding to different values of β, it can be observed that when β is set to 0.8 and 0.9, EGD-CMO can obtain good results. Considering that the last state is to strengthen the convergence of the population, β is set to 0.8 in this paper to ensure that the EGDCMO has sufficient convergence.
The same comparative study is used to analyze the effect of α (i.e., the proportion of infeasible solutions to maintain) on EGDCMO, α is set to 0.1, 0.3, 0.5, 0.7, and 0.9. The Pareto  Fig. 9. From the figure, it is obvious that when α is set to 0.5, 0.7, and 0.9, EGDCMO can find the optimal feasible region and has better convergence performance. According to the convergence profiles of IGD values shown in Fig. 10, as can be seen from the figure, EGDCMO has poor convergence performance on DC2-DTLZ3 when α is set too small or too large, and EGDCMO converges faster when α is set to 0.5. Considering that some offspring solutions are generated around infeasible solutions, it is likely to generate some solutions that can jump out of the infeasible region, which is helpful to explore feasible and infeasible regions. At the same time, to ensure that there are enough feasible solutions to explore the feasible region and consider the tradeoff between diversity and convergence, α is set to 0.5 in this paper.

Experimental study
In this section, to verify the performance of the proposed EGDCMO, the proposed EGDCMO is compared with six other CMOEAs, including NSGA-II, ToP, MOEA/D-DAE, C-TAEA, PPS, and DSPCMDE on the constrained DTLZ test suite, the MW test suite and the LIR-CMOP test suite. All our experiments are executed on the evolutionary multiobjective optimization platform PlatEMO [49].

Parameter settings
1. Problems: The MW [46] problems, constrained DTLZ [50] problems and LIR-CMOP [26] problems are the benchmark problems adopted in this paper. The number of decision variables D and the number of objectivesM of each benchmark problem are set as follows. For the constrained DTLZ, D 7 and M 3 for C1-DTLZ1, the remaining problems. The number of function evaluations for constrained DTLZ, MW, and LIR-CMOP is set to 100,000, 100,000, and 300,000, respectively. 2. Compared algorithms: Four other CMOEAs, including NSGA-II [5], ToP [11], MOEA/D-DAE [12], C-TAEA [16], PPS [10], and DSPCMDE [15] are adopted in this paper to verify the performance of the proposed EGD-CMO. The parameters of all the compared CMOEAs are set as suggested in their original papers. For ToP, it is embedded in the NSGA-II-CDP, where the first phase ended when the feasibility proportion P f is larger than 1/3 or the difference δ is less than 0.2. In addition, the population size is set to 100 on problems with two objectives and 105 on problems with three objectives of all the comparison algorithms and the proposed algorithm in our experiment. 3. Genetic operators: NSGA-II, MOEA/D-DAE, C-TAEA, and EGDCMO adopt the simulated binary crossover and the polynomial mutation to generate offspring, DSPCMDE uses two popular DE [15] operators to generate offspring, whereas PPS and ToP use differen- Performance metrics: Due to the true pareto fronts (PFs) of benchmark problems are known, the inverted generational distance (IGD) [52] metric is employed to assess the performance of all the CMOEAs, and approximately 10,000 reference points are sampled on the PF of the problem to cal-culate IGD. Hypervolume (HV) [53] metric is also used in our experiment. A smaller IGD or a larger HV value indicates a better approximation to the true PF. Table 5 shows the mean value and standard deviation of the IGD values obtained by NSGA-II, ToP, MOEA/D-DAE, C-TAEA, PPS, DSPCMDE, and EGDCMO on the DTLZ and MW problems for 30 independent runs. The Wilcoxon rank sum test at a 0.05 level is presented, where "+", "−" and "≈" denote that the result is significantly better, significantly worse, and significantly similar to that obtained by EGDCMO, respectively. As shown in Table 5, the proposed EGDCMO obtains 18 best results in all 26 problems, which is followed by C-TAEA and MOEA/D-DAE achieving five and three best results respectively, whereas both NSGAII, ToP, PPS, and DSPCMDE achieve no best results. It also can be found that the results of EGDCMO are significantly better than that of NSGA-II, ToP, MOEA/D-DAE, C-TAEA, PPS, and DSPCMDE on 26,26,21,20,26, and 26 problems, respectively. Table 6 presents the comparison results of EGDCMO and its compared algorithms in terms of HV. According to the statistical results shown in Table   Table 5 Statistics of the IGD metric obtained by the proposed EGDCMO and all the compared CMOEAs on the DTLZ and MW problems
From the results regarding HV shown in Table 8 of all the CMOEAs on LIR-CMOP problems, we can see that EGD-CMO are significantly better than that of NSGA-II, ToP, MOEA/D-DAE, C-TAEA, PPS, and DSPCMDE on 14, 12, 7, 13, 5, and 8 test problems, respectively. Figure 15 plots Pareto fronts with median IGD value among 30 runs obtained by the seven CMOEAs on LIR-CMOP9. Although NSGA-II and ToP can cross the infeasible barrier, they cannot jump out of the local optimum. The EGDCMO and all the other four CMOEAs can cross the infeasible barrier, and evolve toward the global PF in the end.

Experimental results on an actual engineering application
Finally, to verify the performance of EGDCMO in solving practical applications, we performed experiments on a multiobjective engineering problem, i.e., structure optimization of a blended-wing-body underwater glider (BWBUG) [54]. When the BWBUG is lifted from the water, stress concentration may arise in the skeleton structure because of the vertical downwards force on the two wings. In the meantime, the frame structure of the glider is equivalent to the cantilever beam structure, and three parts should be considered in the analysis, which involves the gravity of the skeleton, equipment, and buoyancy material. For the equipment, given the specific tasks and functions of this BWBUG, the concentrated force is set to 500 N and the distributed force is set to 4500 N. On the other hand, for the buoyancy material and other accessories, the distributed force is set to 1000 N.
For BWBUG, the larger the internal volume, the more energy, and other equipment it can carry, and at the same time, the weight of the skeleton will increase. In practical use, we hope that the mass of the skeleton will be as light as possible when the stress constraints are met. Therefore, the skeleton mass should be minimized and the internal volume  should be maximized in the optimization of BWBUG. In this structure optimization, the internal layout of BWBUG can be changed, that is, the internal layout will affect the shape. Internal devices are not allowed to exceed the envelope surface of the shape, which is called layout constraint. In order to obtain minimizing skeleton mass and maximize the internal volume, and meanwhile satisfy the stress and layout constraints, the specific design parameters and optimization formula are summarized as follows. Figure 16 shows the shape and layout of the BWBUG parameter diagram, and the shape of the BWBUGs can be divided into plane shape and wing shape. The class and Shape Transformation (CST) method [55]. Reference [56] is used for wing parameterization. There are 20 design variables, including 5 section CST variables, 11 plane shape design variables, and 6 internal layout design variables. Parameter definitions are listed in Table 9. To be specific,D 1 and L are fixed, that is, the total length and width of this BWBUG are 1000 mm and 3000 mm, respectively. The channel of connection buoyancy-regulating cabin and battery cabins are distributed on both sides of the buoyancy-regulating cabin, and the radius and length are set to 20 mm. Figure 17 shows the skeleton full parameter diagram of BWBUG. The skeleton is generated based on the shape, and the overall frame structure of the glider is assembled from the internal frame structure. The internal skeleton could fall into two parts, i.e., fuselage and wing. The generated fuselage is combined with the wing beam and wing rib structure to get the internal structure of the fuselage and wing skeleton.
Subsequently, the skeleton is defined by 5 geometric parameters in Fig. 17. To be specific, t denotes the decision variable, and others relate to geometrical parameters of the shape. Finally, the minimum mass of the skeleton and the maximum internal volume are taken as the two objectives and the stress and layout as two constraints. The expression is as follows: where m sk is the mass of the skeleton structure, V is the internal volume, InterV is the interference volume, σ max is the maximal equivalent stress and σ s is the strength of the aluminum alloy. Considering the safety factor σ s is set to 232 Mpa. According to the initial shape of BWBUG, NACA0020 is selected as the basic airfoil of each section in the improved CST parameterization method. Each airfoil design has the where A Low1 denotes the lower boundaries, and A Up1 represents the upper boundaries of the design scopes, respectively. In addition, the ranges of plane shape variables, layout variables, and skeleton variables are shown in (13), (14), and Skeleton S low [5] S up [10] .
In the optimization process, finite element analysis is used to simulate this actual case. In the present study, the mate- Fig. 19 Solutions obtained by EGDCMO on structure design rial is aluminum alloy (material density: 2700 kg/m 3 ; elastic modulus: 71,000 Mpa; Poisson's ratio: 0.33). Since the wing and fuselage are separated in modeling, the contact between the wing and fuselage is performed in finite element analysis. In terms of IGD, the top three algorithms for obtaining the best results in the benchmark cases are EGDCMO, C-TAEA, and MOEA/D-DAE, we use them to solve this engineering problem. Both EGDCMO, C-TAEA, and MOEA/D-DAE used 3000 FEs to attain the Pareto set, and the population size is set to 30. For a fair comparison, EGDCMO, C-TAEA, and MOEA/D-DAE used the same initial population. Furthermore, the results are summarized in Fig. 18, where the blue dots, black plus signs, and green squares are Pareto front of EGDCMO, C-TAEA, and MOEA/D-DAE, respectively. The red dots are the Pareto frontier obtained from the union of EGDCMO, C-TAEA, and MOEA/D-DAE's Pareto sets. Relatively, 30 Pareto solutions are obtained by EGDCMO, and EGDCMO generated more non-dominated points with larger volume values. 20 Pareto solutions are obtained by C-TAEA, and C-TAEA generated more non-dominated points with larger volume values, while MOEA/D-DAE obtained no Pareto solutions. Moreover, 4 representative points located in the final Pareto front are selected, and their optimized geometric models are shown in Fig. 18 as well. Furthermore, Fig. 19 shows the detailed solutions obtained by EGDCMO, where the blue dots represent feasible designs, pink stars refer to the infeasible solutions, red squares are the infeasible solutions of DOE, and green squares are the infeasible solutions of DOE, and the red circles are Pareto front. Figure 20 shows the plane and section of 4 representative results. Specifically, section1, section2, section3, and section4 of each optimized glider have the same normalized section.
Correspondingly, the equivalent stress diagrams are provided in Fig. 21. In summary, EGDCMO cannot only deal

Conclusions
This paper has proposed an efficient global diversity CMOEA to solve CMOPs with small feasible regions and complex constraints. Specifically, an efficient global diversity strategy is proposed to maintain infeasible solutions well-distributed feature at the early stage of evolution. To this end, a set of weight vectors are used to specify several subregions in the objective space, we iteratively investigate each subregion. A new fitness function is suggested in the proposed algorithm to evaluate infeasible solutions and decide the survival of infeasible solutions in each subregion, which is capable of balancing the importance of constraints and objectives. In addition, the infeasible solutions are ranked higher than the feasible solutions to focus on the search in the undeveloped areas for better diversity.
In the experimental part, the proposed algorithm has been tested on three benchmark cases and has been compared with four existing CMOEAs. The results show that the proposed algorithm has got impressive results and has better overall performance than the compared MOEAs on benchmark CMOPs, especially when the CMOPs are composed of multiple disjoint parts in the feasible region. Moreover, the proposed algorithm has been tested on an actual engineering application, which has also shown better performance than existing MOEAs.