A Survey on Evolutionary Constrained Multiobjective Optimization

Handling constrained multiobjective optimization problems (CMOPs) is extremely challenging, since multiple conflicting objectives subject to various constraints require to be simultaneously optimized. To deal with CMOPs, numerous constrained multiobjective evolutionary algorithms (CMOEAs) have been proposed in recent years, and they have achieved promising performance. However, there has been few literature on the systematic review of the related studies currently. This article provides a comprehensive survey for evolutionary constrained multiobjective optimization. We first review a large number of CMOEAs through categorization and analyze their advantages and drawbacks in each category. Then, we summarize the benchmark test problems and investigate the performance of different constraint handling techniques (CHTs) and different algorithms, followed by some emerging and representative applications of CMOEAs. Finally, we discuss some new challenges and point out some directions of the future research in the field of evolutionary constrained multiobjective optimization.

dispatch (CEED) problem [2], the urban bus scheduling problem [3], and energy saving optimization problem [4]. This kind of problems is denoted as constrained multiobjective optimization problems (CMOPs). CMOPs bring more challenges than their unconstrained counterparts, due to the coexistence of multiple objectives and constraints. For CMOPs, because of the conflict among multiple objectives, there exist no single solution but a set of solutions that indicates the best possible tradeoffs among the objectives. The set of solutions forms the so-called Pareto-optimal set (PS), and its image in objective space is called the Pareto front (PF). Meanwhile, because the presence of the constraints can change the landscape of the search region and cause infeasible regions, the feasibility of the solutions cannot be neglected. Consequently, the final goal of solving CMOPs is to obtain a set of feasible solutions to approach the PF with good convergence and diversity.
In the past two decades, multiobjective evolutionary algorithms (MOEAs) have become a very popular alternative for addressing the multiobjective optimization problems [5], [6]. The most important reason is that MOEAs belong to population-based stochastic metaheuristics that utilize the principles of natural selection to guide a set of solutions toward the PF. As a result, MOEAs are capable of providing a set of solutions for a multiobjective optimization problem in a single run. Up to now, a large number of MOEAs have been proposed, and they can be roughly divided into the methods based on dominance [7]- [9], the methods based on decomposition [10]- [12], and the methods based on indicator [13]- [15]. These MOEAs have shown competitive performance in solving unconstrained multiobjective optimization problems. However, they can not be directly applied to CMOPs due to the absence of a constraint handling technique (CHT). In fact, CHT can be seen as a selection mechanism to cope with constraints. In various practical application scenarios, CMOPs are more common than the unconstrained multiobjective optimization problems. Hence, a large amount of effort has been devoted to designing CHTs and specific mechanisms, promoting the emergence and development of various constrained MOEAs (CMOEAs) [16]. Meanwhile, to investigate the performance of different CMOEAs, a number of benchmark test suites have also been designed in recent years [17], [18].
Despite numerous studies that have been conducted on constrained multiobjective optimization in the community of evolutionary computation, to the best of our knowledge, there has been no comprehensive survey of them presented so far. Hence, this article aims at providing a systematical survey on evolutionary constrained multiobjective optimization regarding the following aspects.
1) On the basis of the categorization shown in Fig. 1, existing CMOEAs are divided into seven classes, and each of them is introduced. Moreover, the advantages and limitations of each class are discussed. It is expected that this introduction of CMOEAs might help researchers to deeply and thoroughly understand the advanced techniques in solving CMOPs. Furthermore, the performances of different CHTs and algorithms are both investigated.
2) The benchmark test problems for evaluating the performance of CMOEAs are introduced and analyzed, and these problems are classified according to their respective characteristics. Hopefully, this might facilitate the performance comparison among different CMOEAs and the development of new CMOEAs.
3) The applications of evolutionary constrained multiobjective optimization in some typical fields are introduced. This not only provides a practical perspective for researchers but can also inspire practitioners to tackle CMOPs in other similar areas by means of CMOEAs. 4) Although the existing CMOEAs have exhibited competitive performance in solving many general CMOPs, as more complicated features are continually emerging in CMOPs, a number of challenges still need to be considered, such as the CMOPs with dynamic, multimodal, computationally expensive, or large-scale characteristics. Therefore, this article gives some topics of promising future research to promote the development of CMOEAs on more types of challenging CMOPs. The remainder of this article is structured as follows. The basic concepts on CMOPs are introduced in Section II. The classification and introduction of the existing CMOEAs are described in Section III. Section IV introduces the benchmark test problems in detail, followed by the performance comparison of different CHTs and algorithms in Section V. The specific practical application scenarios are introduced in Section VI. Some directions of potential future research are outlined in Section VII. Finally, the conclusions are drawn in Section VIII.

II. BACKGROUND
Generally, a CMOP can be described as where S is decision space and x is the decision vector with D-dimensions in S. F( x) is the objective vector, which contains m objectives that need to be solved. g j ( x) ≤ 0 represents the jth inequality constraint. h j ( x) = 0 is the (j − l)th equality constraint. l and (k − l) denote the number of inequality and equality constraints, respectively. Note that the decision vectors can be continuous or/and discrete and, thus, this article focuses on both continuous and discrete problems. For a decision vector x, its constraint violation degree on the jth constraint is calculated as where δ is a positive tolerance value to relax the equality constraints. By adding the constraint violation degree of x on each constraint, its total constraint violation degree can be obtained According to the total constraint violation CV, the feasibility of one solution can be identified. If the CV of one solution is 0, it is called feasible solution; otherwise, it is an infeasible solution.
For two feasible solutions x 1 and x 2 , x 1 is said to dominate for at least one j ∈ {1, . . . , m}. If a solution is not dominated by any solution, it is called the Pareto-optimal solution. The set composed of all Pareto-optimal solutions is PS, and the mapping of PS in objective space is PF.
The presence of constraints brings more difficulties in obtaining the PF of a CMOP than its unconstrained counterpart [i.e., ignoring the constraints in (1)]. The constraints can divide the feasible space into narrow disconnected regions, make a large proportion of the search space infeasible, and/or make the PF of unconstrained MOP partially or completely infeasible. For the sake of distinction, constrained PF (CPF) and unconstrained PF (UPF) are commonly used to denote the PFs for constrained and unconstrained multiobjective optimization problems, respectively [18]- [21]. Obviously, obtaining the UPF is the goal of unconstrained multiobjective optimization, while covering the CPF with good convergence and distribution is desired in CMOPs. To be specific, a good convergence means that the solutions can approach the PF accurately, and a good distribution stands for the solutions featuring good spread and evenness in the objective space.
Specifically, the following main challenges will be encountered when tackling the CMOPs [23].
1) Difficulty of Feasibility: As shown in Fig. 2(a), constraints may cause a larger proportion of the search space to become an infeasible region. As a result, a very small feasible region is hard to be located in.
2) Difficulty of Convergence: As shown in Fig. 2(b), the infeasible region may become a barrier on the way to the CPF. An algorithm is easy to be trapped in the feasible outside region, thus finding the true CPF is not easy. 3) Difficulty of Diversity: As seen in Fig. 2(c), the disconnected feasible region will make the true CPF consists of multiple disconnected segments, thus covering all of them is a challenging task. In view of the aforementioned challenges, achieving the balance among convergence, diversity, and feasibility is important and critical in solving CMOPs. This research conception has been widely adopted by researchers and a variety of CMOEAs has been tailored in the past years. The details of these CMOEAs will be introduced in the following section.

III. SUMMARY OF EXISTING CMOEAS
In this section, the currently available CMOEAs are classified and introduced. It is worth noting that CHTs are indispensable when combining MOEAs to form CMOEAs. Furthermore, for most of the existing CMOEAs, their efforts usually focus on the CHTs. Thus, as shown in Fig. 1, we use a taxonomy referring to the CHTs instead of the type of MOEAs. The existing CMOEAs can be divided into seven categories: 1) methods based on the penalty function [24], [25]; 2) methods based on the separation of objectives and constraints [26], [27]; 3) multiobjective methods (MOs) [28], [29]; 4) methods of transforming CMOPs into other problems [19], [30]; 5) hybrid methods [31], [32]; 6) methods of altering the reproduction operators [33], [34]; and 7) other methods.

A. Methods Based on Penalty Function
These methods are mainly based on the constraint violation degree to construct the penalty term. By adding the penalty term into objective functions, the constrained optimization problem is transformed into an unconstrained one.
How to set the penalty coefficient is the most critical issue in the penalty function methods, playing a decisive role in terms of efficiency of the algorithm. According to the different setting methods, the penalty method can be divided into the static method, dynamic method, and adaptive method. Obviously, the static method means that the penalty coefficient will not change during the evolution process, but in this way, the preference between objectives and constraints will remain unchanged in the early and late stages of evolution, which is not conducive to the balance between the objectives and constraints. This kind of methods effectively solves simple problems, but it cannot achieve the ideal effect when solving complex problems. Dynamic method is a method where the penalty coefficient changes with the increase of evolutionary generations or other indicators. Through the change of penalty parameters, the weight of objectives and constraints is different in different stages of evolution, so as to achieve the balance between them. However, the difficulty of this method is to design the appropriate change rules, and different change rules need to be set on different problems. The adaptive method preserves the information of the population in the evolution process, and then feeds the obtained information back to the population to adjust the penalty coefficient. Thus, compared with the above two methods, the adaptive strategy can obtain better performance for more complicated problems, since the feedback information can provide guidance of the potential evolution.
Jadaan et al. [35] combined a new nondominated ranked genetic algorithm (NRGA) with a nonparameter penalty method to continuously search for the better PS. Zhang et al. [36] designed a multihive colony artificial bee colony method to deal with CMOPs, by employing a general coevolution model. In this method, the adaptive penalty function is employed to deal with the constraints, which penalizes infeasible individuals by the number of feasible individuals in the ant colony. By using the proportion of current feasible solutions, Woldesenbet et al. [37] employed the distance measure in the objective space to adaptively guide the population evolution. Jan and Zhang [38] employed a threshold to control the degree of punishment for infeasible solutions. Based on this work, Jan et al. [39] put forward an adaptive and dynamic version for the threshold, and then embedded it into the MOEA based on decomposition (MOEA/D) framework to deal with CMOPs. Jiao et al. [40] suggested a feasible-guiding strategy, which utilizes the degree of constraint violation to modify the objective functions to get the new fitness. Furthermore, in order to achieve the balance between the objectives and constraints, the penalty coefficient is adjusted by the feasibility proportion in the current population. Maldonado and Zapotecas-Martínez [41] proposed a dynamic penalty function method, in which the parameters gradually change with generation number and are embedded into the MOEA/D framework to solve CMOPs. Vaz et al. [42] designed a threestep penalty, in which three different penalty coefficients were used at different stages of evolution. Inspired by deep learning, Fan et al. [43] proposed a learning-guided parameter setting method, which can adaptively generate penalty parameters. This penalty function combines constraint violation degree, objective function values, and current generation to design an exponential decay model, which is embedded into the push-and-pull search (PPS) [19] framework to solve CMOPs.
The following algorithms have similar mechanisms to the penalty function methods, so they are also classified in this category. Ma et al. [44] designed a CHT based on two rankings, named ToR, in which each individual is sorted in two ways: one based on the objective function and the other based on the constrained dominance principle (CDP) [45]. Then, according to the ratio of feasible solutions in the current population, the two rankings are adaptively weighted to obtain a new fitness function. However, its parameters change from 0.5, which means that the preference on constraints is always equal to or greater than the objectives. In addition, because only the proportion of feasible solutions is considered. If the proportion is equal to 0, the weight of objectives and constraints is always equal, which will make the algorithm fail to find feasible solutions on some complex problems. To alleviate these issues, Yu et al. [46] proposed a dynamic selection preference-assisted mechanism to transform the preference between objectives and constraints. To be specific, the preference for the objectives is gradually reduced from 1 to 0 according to the designed cosine function curve [47], [48]. Ming et al. [20] also carried out two kinds of rankings for each individual, one based on convergence and the other based on diversity. Finally, according to the generation information, these two kinds of rankings are weighted by the rule of the sigmoid function. In order to guide the population in entering the feasible region from different directions in the early stage of evolution, and approach the Pareto-optimal solutions in the later stage, Ma and Wang [25] proposed a shift-based penalty method, named ShiP, in which the infeasible solutions are first shifted according to the distribution of its adjacent feasible solutions, and the degree of migration is adaptively controlled by the proportion of feasible solutions. Then, the shifted infeasible solutions are punished considering their constraint violations. García et al. [49] used the penalty function method to select a group of points that are close to the feasible space and have small objective functions. Moreover, the method is embedded in the indicator-based framework to solve the equality CMOP.

B. Methods Based on the Separation of Objectives and Constraints
These methods compare objectives and constraints separately, mainly including the CDP [45], ε constrained method [50], [51], and stochastic ranking (SR) [52].
1) CDP: CDP was proposed by Deb et al. [45], which is most commonly used due to its simplicity and ease of implementation. CDP compares paired individuals A and B using the following criteria.
1) When both individuals A and B are feasible solutions, if A Pareto dominates B, A is selected to enter the next generation. 2) When A is the feasible solution and B is the infeasible solution, A is selected. 3) When both A and B are infeasible solutions, the individual with the smallest constraint violation degree is selected. The operation of CDP is relatively simple, but it prefers feasible solutions, which will cause the population to fall into some local feasible regions when CMOPs have discrete feasible regions or infeasible barriers. In order to remedy this shortcoming of CDP, Jimenez et al. [53] employed the niche technology to increase the diversity of the population. Saha and Ray [54] proposed an equality constraint repair method based on probability points, which combines CDP and clustering methods to repair infeasible solutions. With the aim of using the effective information of infeasible solutions, Wei and Wang [55] applied CDP to the infeasible elitist-based particle swarm algorithm. In the early stage, individuals with good objective values will be retained, and individuals with small objective values and constraint violations will be saved in the later stage. Based on the NSGA-II and MOEA/D framework, Jain and Deb [56] proposed NSGA-III and C-MOEA/D to constrained multiobjective optimization. Fan et al. [57] designed a new angle-based CHT, named ACDP, which integrates angle information into the CDP, and further uses the information of infeasible solutions. Later, the ACDP was embedded into the MOEA/D framework to solve CMOPs [58]. In addition, Ning et al. [59] proposed an improved hybrid multiobjective optimization algorithm to solve CMOPs, in which each solution is given a constrained nondominated rank according to its constraint violation degree and Pareto rank. Wang and Xu [27] used the ACDP to solve the constrained many-objective optimization problem, in which individuals with good diversity will be first selected to cope with the challenge caused by the large infeasible obstacles.
2) ε Constrained Method: This method was proposed by Takahama and Sakai [50], which uses a parameter ε to relax constraints. ε is gradually reduced, and when the constraint violation degree of an individual is less than ε, it is considered as a feasible solution. Obviously, when ε is reduced to 0, the ε constrained method is the same as CDP. When comparing individuals A and B, the ε constrained method uses the following criteria.

ε, and CV(A) < CV(B), then
A is selected. Saxena et al. [60] integrated the ε constrained method into the NSGA-II framework to control the infeasibility of the population and enhance the convergence of the population. Zapotecas-Martínez et al. [61] proposed a multiobjective particle swarm optimization algorithm (MOPSO) based on decomposition, and the ε constrained method is integrated into it to solve CMOPs. Yang et al. [62] also integrated the ε constrained method and adaptive operation selection into the multiobjective framework based on decomposition. Becerra et al. [63] used the ε constrained method to obtain several points on (or very close to) the PF, and then the method based on the rough set was employed to extend these solutions to cover the whole PF. Martinez and Coello Coello [64] showed a method based on the ε constrained method, in which the relevant information of the neighborhood in MOEA/D is used to obtain the promising solutions in the allowable feasible region. The ε constrained method was improved and embedded into the MOEA/D framework by Fan et al. [65], in which the proportion of feasible solutions in the current population is employed to dynamically adjust the ε parameter level. Yang et al. [66] designed a multiobjective differential evolution algorithm (MODE-SaE) based on the improved ε constrained method, in which the ε level can be adaptively adjusted based on the maximum and minimum constraint violation values of infeasible individuals. Zapotecas-Martínez and Ponsich [67] established the CMOP as a biobjective problem by considering the constraint violation degree and scalar function as objectives, and integrated the ε constrained method into MOEA/D to solve the constructed problem. Yang et al. [68] proposed a dynamic constraint handling mechanism, which divides the search process into two modes: 1) unconstrained search and 2) constrained search. In the constrained search mode, an improved ε constrained method is used to improve population diversity. To ensure the diversity of PF solutions, Wang et al. [69] combined the ε constrained method with a niche strategy to solve CMOPs. To prevent the population from falling into the local feasible region or the infeasible region, Zhu et al. [70] proposed a technique to help the population escape from the stagnation status, and then an improved ε constrained method is used to search for CPF.

3) SR:
In this method, a probability parameter pf is introduced. When two individuals are compared, the probability with pf only compares their objective function value, and the probability with (1-pf ) compares their constraint violation degree. This method is able to use the information of the objective function to some extent. Geng et al. [71] used SR to balance objectives and constraints when solving CMOPs. Jan and Khanum [72] embedded the modified SR into the MOEA/D framework to solve CMOPs. Ying et al. [73] proposed an adaptive stochastic ranking mechanism, which dynamically adjusts the probability parameter according to the current evolution stage and the difference of individuals' violation degree. Liu et al. [74] studied CMOEAs based on indicators by combining the indicator-based MOEA with CDP, the ε constrained method, and SR, respectively. Gu et al. [75] proposed an evolutionary algorithm based on the surrogate, in which an improved SR strategy based on fitness mechanism and adaptive probability operator was proposed. This strategy considers the convergence and diversity to improve the quality of candidate solutions. Fig. 3 shows the methodologies of CDP, the ε constrained method, and SR. x represents an individual and its objective value and constraint violation degree are f ( x) and CV( x), respectively. When the other individuals are compared with x, if they are located in the blue shadow, it means that they are better than x. Obviously, when CDP is used to select individuals, it will give priority to feasible solutions. In this way, the convergence rate of the population is very fast, but it will also cause the population to easily fall into the local optimum. While both the ε constrained method and SR can make up for the defect of CDP to a certain extent, since they consider the information of infeasible solutions. However, the parameter settings in the ε constrained method and SR are also challenging because the parameter may be problem-dependent.

C. Multiobjective Methods
In the MOs, the constraints are regarded as one additional objective or multiple additional objectives, the CMOP is transformed into an unconstrained counterpart, and then the MOEAs can be employed to solve the converted problem. The objective functions after transformation are as follows: In [76], the constraints were transformed into two new objectives: one is based on the penalty function, the other is equal to the number of violation constraints. Based on the constraint violation measure, Ray et al. [77] added a new objective, and designed an infeasibility-driven evolutionary algorithm to promote the population to approach the constraint boundary from the infeasible region. Isaacs et al. [78] regarded constraints as a new objective and redefined the constraint minimization problem of the original m objectives as an unconstrained minimization problem with m + 1 objectives. Long [79] constructed a novel CHT for solving CMOPs, in which the convergence, diversity, and feasibility of the obtained solutions are used as three new objectives for the multiobjective subproblem. Taking the constraint violation degree as a new objective, Peng et al. [29] designed a new CHT based on directed weights to solve CMOPs, in which two types of weights, respectively, distributed in feasible and infeasible regions are employed to guide the search toward the promising regions. Zhou et al. [80] proposed a tri-goal evolutionary framework for solving constrained many-objective problems. This framework transforms constraints into feasibility indicators, and then is combined with convergence indicator and diversity indicator to form three new objectives.

D. Methods of Transforming CMOPs Into Other Problems
In order to more efficiently solve CMOPs, many researchers convert CMOPs to other problems, such as converting a CMOP to collaborative optimization problems or two-stage optimization problems. After converting, some promising operators like coevolution can help the population to better explore the search space, discover some new and potential information, and obtain the complete CPF in the end.

1) Converting CMOPs Into the Collaborative Optimization Problems Based on Multipopulation:
Chafekar et al. [81] transformed a CMOP into multiple single-objective optimization problems, in which several genetic algorithms are implemented, and each of them optimizes an objective and then exchanges objective information. Wang et al. [82] proposed a cooperative differential evolution (DE) framework also in the use of m subpopulations, each of which optimizes an objective with constraints, i.e., single constrained optimization. In addition, an archive population is used to preserve the obtained constrained nondominated solutions to approach the CPF. Liu and Wang [83] proposed a CMOEA based on decomposition and temporary archive, which decomposes a CMOP into several subproblems, and each of the subproblems has its own subpopulation and temporary archive. Then, each subproblem is optimized by a coevolution strategy. Later, they proposed a constraint handing scheme based on boundary search and archive [84], in which a CMOP is decomposed into several subproblems, and each subproblem has its own archive. In addition, the boundary search strategy is constructed to improve the efficiency of the algorithm. In order to maintain the distribution of solutions, Yang et al. [85] decomposed CMOP into multiple subproblems by partitioning the objective space and used multiple CHTs to solve the optimization problems. Liu et al. [86] changed a CMOP into a two-population optimization problem, where one population only optimizes the constraints, and the other population focuses on optimizing the objectives. Meanwhile, the interaction of information and the transfer of knowledge exists between these two populations. Tian et al. [87] proposed a coevolutionary framework (CCMO) for solving CMOPs, in which one population is using to solve the original CMOP, that is, to search for CPF, and the other population ignores constraints to find UPF. In addition, these two populations assist each other to solve CMOPs. Li et al. [88] designed a dual-archive evolutionary algorithm (C-TAEA), one is a convergence-oriented archive (CA) that aims at pushing the population along the PF3 and the other is a diversity-oriented archive (DA) used to explore the undeveloped areas of CA and maintain population diversity. Wang et al. [89] designed a cooperative MOEA by employing two populations: 1) propulsive population and 2) normal population, Concentrating on convergence, the propulsive population considers no constraints in the early stage, while in the late stage, only constraints are considered. The normal population searches the whole CPF and prioritizes the feasibility and the diversity of the population. A bidirectional coevolution algorithm was designed by Liu et al., in which the main population and an archive population are employed simultaneously [90]. Specifically, the main population keeps the feasibility and moves from the feasible side to the CPF, while the archive population uses angle information to maintain population diversity and approximate the CPF from the infeasible side.
2) Converting CMOPs Into Two-Stage Optimization Problems: Santana-Quintero et al. [91] combined MOEA with the local search method based on rough set theory to solve CMOPs. In the first stage, MOEA was used to approach the PF. In the second stage, the fuzzy set theory was used to improve the diversity and the convergence of PF solutions. Fan et al. [19] proposed a PPS framework. The goal of the push stage is to cross the infeasible region to the UPF. In the pull stage, the improved ε constrained method is employed to search the CPF. Based on the PPS framework, a CMOP was decomposed into a group of simple subproblems [92], with each subproblem corresponding to a subpopulation, and the PPS framework is applied to each subpopulation to solve the related subproblem. Based on the diversity distance measure, Wang et al. [93] proposed an objective space modified mechanism, which enables the promising infeasible solutions to be more effective in finding the optimal solutions. In addition, the PPS search framework is used to adjust the position of the PF to prevent the population from falling into the local optimum, and the time complexity can be reduced. Garcia-Garcia et al. [94] utilized the characteristics of cellular genetic algorithms (CGAs) and combined PPS technology to solve CMOPs. Tian et al. [95] designed a two-stage evolutionary algorithm, named CMOEA-MS, in which one stage can help the population reach the feasible region, and the other stage can make the population spread along the feasible boundary. In addition, based on the status of the population, the algorithm can adaptively switch between these two stages. For CMOPs with constraints in both decision and objective spaces, Liu and Wang [96] suggested a two-stage optimization idea. The first stage is a single objective problem, which aims to search for the promising feasible region; the second stage requires searching for the final PF and reaching the CPF. Xiang et al. [21] proposed a two-stage algorithm, named CIC-MOEA/D, the goal of the first stage is to find UPF, so only the objectives are considered; moreover, the constraints are gradually emphasized with the goal of approaching the CPF well in the second stage. In [97], the whole evolutionary process was divided into two stages. The first stage aims to keep the balance between convergence and diversity, and the second stage is devoted to maintaining feasibility and diversity and, thus, covering a well-distributed PF. Yu and Lu [98] proposed a corner point algorithm based on DE by including two stages. The first stage is to find the corner points, and the second stage is to search the real CPF. Zhang et al. [30] used the framework of artificial bee colony to divide the optimization process into two stages. In their first stage, fast nondominated sorting is employed in promoting the population to reach the PF. In the second stage, the Tchebycheff method is employed to improve population diversity. Ming et al. [99] proposed a two-stage evolutionary algorithm, in which the first stage is used to find UPF and store the feasible solutions obtained, and the second stage concentrates on exploring CPF.

E. Hybrid Methods
Although EA can find promising areas from global perspectives, its subtle search in local regions may be inferior to mathematical programming (MP) [100], [101], especially facing equality constraints that usually reduce the dimensions of search space. Many constrained optimization algorithms relax equality constraints into inequality ones, which fails to ensure the strict feasibility of the obtained solutions. Some researchers have combined EA with MP to deal with CMOPs. Kelner et al. [31] proposed a hybrid optimization technology by combining genetic algorithm with local search strategy based on interior-point method. Datta et al. [102] proposed a new multiobjective optimization technology, which combines an interior-point method with the nonconvex radial boundary intersection. The radial boundary intersection decomposes the multiobjective optimization problem into several subproblems, which find the solutions of the nearest reference point radially outward along equidistant lines. Morovati and Pourkarimi [103] used the Zoutendijk method to solve CMOPs, where the algorithm considers all objectives and constraints, and a convex quadratic subproblem is proposed to generate a convenient improved feasible direction. Schütze et al. [104] proposed a method to calculate the search direction by using neighborhood information. In this method, with regards to a given point in the population, the greedy search direction of given data is calculated in the use of the neighborhood solution of this point. Lara et al. [32] proposed the construction method of subspace-based movements in the search process, in which the multiobjective stochastic local search is used to guide individuals. Uribe et al. [105] hybridized MOEAs with a special local search mechanism, proposed a descent directions calculation method for the biobjective optimization problem, in which the constraints information is considered in local search. Cuate et al. [106] combined MOEA with continuation-like technology to obtain a fast and reliable numerical solver. Hernández et al. [107] proposed a method of mixing the hypervolume Newton method and evolutionary strategy to obtain a fast and reliable algorithm to deal with CMOPs. Schütze et al. [108] combined the gradient subspace approximation with EA to solve CMOPs, which allows the descent direction to be calculated in a best-fitting manner from the given neighborhood information.

F. Methods of Altering the Reproduction Operators
This type of methods focuses on the design of efficient and specific operators in the process of reproduction. Qian et al. [109] combined an adaptive DE with the ε constrained method to solve CMOPs, in which the generation strategies of trial vector and DE parameters are gradually and adaptively adjusted according to the knowledge learned in the previous search process. Yu et al. [110] designed a new mutation mechanism to deal with infeasible and feasible solutions, which can produce a well-distributed CPF. Qu and Suganthan [111] proposed a diversity promotion mechanism to prevent the population from falling into the local optimum. Xu et al. [33] designed a new DE variant with an infeasibleguiding mutation operator to solve CMOPs, in which the good infeasible solutions are employed to guide the population into the promising region. Wang et al. [112] proposed an adaptive DE with Pareto dominance, in which the parameters are adjusted adaptively. Ramesh et al. [34] proposed an improved generalized DE for CMOPs by replacing the commonly used binomial/exponential crossover with the simulated binary crossover. Liu et al. [113] designed a new strategy based on group-sorting to maintain population diversity, and the crossover operator was designed for each subpopulation to improve the exploration ability. In order to solve CMOPs with a large number of objectives or variable dimensions, He et al. [114] emphasized the role of offspring generation when reproducing promising feasible or valuable infeasible offspring solutions. Miyakawa et al. [115] proposed a direct matching strategy, which can assign a unique search direction to individuals according to their positions in the objective space. Afterwards, they further proposed a method to control the selection area of promising infeasible solutions by determining the dominant region of the solutions [116].

G. Other Methods
Herein, some methods that do not belong to any of the above categories are discussed. For some problems, it may be more beneficial to allow the population to perform nongreedy or uphill movement, such as the simulated annealing (SA) [117]. Singh et al. [118] integrated the surrogate assisted method into SA to solve CMOPs. The optimization method based on the surrogate model can reduce the number of computationally expensive simulations. Singh et al. [119] combined the Kriging model with improved multiobjective probability and feasibility probability to solve CMOPs. Datta and Regis [120] proposed a surrogate-assisted evolution strategy to solve CMOPs with expensive black-box inequality constraints. Zhang and Qian [121] employed an artificial immune system model for solving nonlinear constrained multiobjective problems in a time-varying environment. Based on the super rectangle search, Wei and Wang et al. [122] proposed a PSO variant, in which the optimal solution of the next time step can be predicted. Considering that the single CHT cannot perform better than other CHTs in every problem, Qu and Suganthan [123] hybridized multiple CHTs to solve CMOPs. Zapotecas-Martínez and Coello Coello [124] proposed an algorithm for continuous box-CMOPs, in which the nonlinear simplex search scheme is used to obtain the PS. Sadollah et al. [125] used the water cycle algorithm to solve CMOPs, in which one archive is employed to save the obtained nondominated solutions. Based on the biological immune system, Qian et al. [126] designed an immune optimization method for solving CMOPs. Yuan et al. [16] introduced an indicator-based CHT for CMOPs to guide the population to uniformly explore the promising areas. Li et al. [127] designed a pioneer selection strategy, in which the pioneer individuals give no consideration to constraints, and the proportion of the pioneer individuals is gradually reduced to obtain a well-distributed CPF.

H. Discussion
The relationship between CHTs and MOEAs is discussed here. In general, one CMOEA is composed of CHT and MOEA. To be specific, the CHT is used to tackle the constraints and ensure the feasibility of solutions, while MOEA is taken as the main evolution power to approximate the PF with a set of well-converged and well-distributed feasible solutions. Naturally, designing CMOEAs by integrating MOEAs (i.e., dominance-based [9], decompositionbased [10], and indicator-based algorithms [15]) with CHTs is a practicable and simple way to cope with CMOPs, since all of them have been widely studied individually. Based on the current studies, it is observed that CHTs are often combined with dominance-based and decompositionbased MOEAs [44], [58], since they have better flexibility and generalization than the indicator-based MOEAs. Especially, for decomposition-based MOEAs, they transform a multiobjective optimization problem into multiple singleobjective optimization problems. Therefore, many CHTs applied to constrained single-objective optimization can be directly embedded. For the indicator-based MOEAs using quality indicators to define selection mechanisms, they have represented a promising way to deal with many-objective optimization problems due to the increase of selection pressure [15]. Thus, combining indicator-based MOEAs with CHTs is a viable alternative for solving constrained manyobjective optimization problems.
In addition, based on the above reviews on existing CMOEAs, the advantages and limitations of each category are summarized in Table I. Although the existing algorithms have exhibited excellent performance in solving CMOPs, some research gaps need to be considered. 1) CMOPs can be solved on the principle of multitask optimization [128]- [130]. Based on the knowledge transfer, evolutionary multitask (EMT) has shown exceedingly promising performance in solving multiple different but related tasks simultaneously. Naturally, the optimization of objectives and the satisfaction of constraints can be seen as two related tasks, thus some knowledge transfer mechanisms in EMT can be readily adopted to solve the CMOPs. 2) In order to assist the solving of CMOPs, some auxiliary tasks can be created by dynamically considering some ones of all constraints. Moreover, for different auxiliary tasks, different CHTs can be employed to exert the complementary effects.
3) The fitness landscape of CMOPs can be analyzed. Based on the characteristics of fitness landscape, more suitable CHTs and evolution operators can be recommended to a specific CMOP. 4) The idea of hyperheuristic algorithm [131], [132] is also beneficial to the utilization of CHTs. The hyperheuristic algorithm has the self-learning ability, which uses the feedback information of heuristic algorithms to improve its performance. Thus, some hyperheuristic strategies are expected to be designed to efficiently use the advantages of different CHTs in the evolution process.
SRN, TNK, and OSY are the first three test problems proposed in the literature, and are often employed for the performance comparison of algorithms. The constraints used in SRN eliminate the partial unconstrained PS, that is, its CPF is part of UPF. The feasible objective space of TNK is the same as the feasible decision variable space. It has a disconnected CPF on the boundary of nonlinear constraints. OSY contains six constraints, and its Pareto-optimal region consists of five segments, each segment located at the intersection of some constraints. The main limitations of these three test problems are as follows: 1) low dimension; 2) most objective functions and constraints are not completely nonlinear, so it is easy to find good solutions; and 3) the complexity of different constraints introduced by them in constrained optimization is not adjustable.
Based on the main defects of the above three problems, the CTP test suite was designed. In the CTP test function set, the difficulty of constraint functions can be adjusted. The constraint functions used mainly provide two types of difficulties: 1) the difficulty of approaching CPF and 2) the difficulty of the whole search space. The CTP test suite contains seven test problems, CTP1-CTP7. The landscape of the region near CPF in test problem CTP1 is complicated, because each constraint is an implicit nonlinear function of decision variables. The presence of constraints makes the search area that is near to the PF infeasible. In addition, by using more constraints, multimodal, and deceptive functions, the complexity of the test problem can be further increased. The test problems CTP2-CTP7 make it difficult to optimize the whole search space. The constraints of CTP2-CTP7 contain six parameters, and the constraint test problems with different difficulties can be obtained by adjusting these six parameters. Although the CTP problem set has a certain improvement compared with the previous three, there are still some defects, i.e., the dimension is not high enough, the feasible region is large, and all problems are restricted to two objectives. In addition, although the number of constraints of CTP1 is adjustable, the shape of its CPF remains unchanged with the increase of constraints. The shapes of CPFs for CTP2-CTP7 are adjustable, but these problems have only one constraint and usually have large feasible regions.
Based on the framework of CTP test suite, the CF test suite is designed. The CF test set has a nonlinear PF in decision space, its CPF has unconnected geometry, and the CF problems also introduce complex variable linkages, so the algorithm finds it difficult to converge when solving the CF problems. However, the CF problems have some defects. For example, the difficulty is not adjustable, the feasible regions are large, and the number of objectives is not extensible.
The NCTP test set is developed to address the defects of CTP problems, and has the following characteristics: 1) the distance function adopts the Rosenbrock function, so the convergence difficulty is increased; 2) high-dimensional decision space is considered; and 3) in order to explicitly reduce the feasibility ratio, an additional constraint is added. In the NCTP test suite, the CPFs of NCTP1-NCTP6 are discontinuous, the CPFs of NCTP7-NCTP12 are composed of continuous and discontinuous parts, and the CPFs of NCTP13-NCTP18 are continuous.
The C-DTLZ test suite is designed on the basis of the unconstrained DTLZ test problems [143]. The C-DTLZ test set mainly contains three kinds of constraints. An infeasible barrier can be provided by introducing the first kind of constraints, while the UPF is still feasible. In addition, several isolated feasible regions along the UPF are defined by the second kind of constraints, and the CPF is part of UPF. The third type of constraints makes the whole UPF infeasible, and CPF is formed by the boundary of the feasible region. Since the C-DTLZ test suite is scalable with multiple difficulties, such as multimodality, irregularity, and bias, it has been widely used to evaluate the performance of CMOEAs. However, the position variables on the PSs of the C-DTLZ test suite lack internal independence, which is also its main limitation.
Fan et al. [140] suggested eight constrained optimization problems, CMOP1-CMOP8. In these problems, a total of three types of constraint functions are defined. The first type of constraints can determine the shape of the CPF, the second type can determine the feasible proportion of the search space, and the third type can simultaneously focus on the CPF shape and the proportion of feasible solutions in the entire search space.
The LIR-CMOP test suite contains 14 CMOPs with small feasible regions, and there exists a complex relationship between location and distance variables. The general characteristic of LIR-CMOP test problems is that their true CPFs are blocked by a large number of infeasible regions, making it difficult to be discovered during the evolution process. Its constraint function is composed of a controllable shape function and distance function. Specifically, the shape function is used to make the shapes of CPFs convex and concave, and the distance function can adjust the convergence difficulty.
The MW test suite contains 14 test problems with various characteristics. Some problems have narrower parts in their continuous/discontinuous feasible regions, and the CPFs of some problems contain only a few isolated Pareto-optimal solutions. Hence, the CPFs of MW test problems are relatively difficult to find. In the construction of the MW test set, a global control process and a local adjustment process are introduced, the former can control the size of the feasible region and the latter is able to adjust the complexity of the boundary in the feasible region. The purpose is to generate CPF with different geometric constraints. In addition, in the MW test set, because the distance function can be extended to any number of decision variables, all test problems are scalable in terms of the number of decision variables. Regarding the number of objectives, only test problems MW4, MW8, and MW14 are scalable.
The previous constraint test sets have no consideration to the constraints both in decision and objective spaces. In light of this case, the DOC test set is designed including nine problems. The characteristics of the DOC test set are as follows: 1) it contains objective constraints and decision constraints simultaneously; 2) it contains both inequality constraints and equality constraints; 3) its CPFs have a variety of properties, for example: continuous, disconnected, convex, concave, linear, mixed, degenerate, and multimodal; and 4) the feasible region in the decision space also exhibits various properties, such as nonlinearity, minima, and multimodality. The disadvantage of DOC is that the dimensions are fixed and not adjustable. DAS-CMOP test suite has adjustable difficulty and expandable objective by employing three primary types of difficulties to characterize the constraints. The level of difficulty in DAS-CMOP is defined by a triplet with each of its parameters specifying each primary difficulty-type level, respectively. In addition, by combining the three primary constraint types with different difficulty triplets, a variety of constraints can be provided for CMOPs.
Zhou et al. [17] designed 16 scalable constrained test problems (CF1-CF16) by taking the variable dependencies into account. The constraints introduced by these problems can be divided into convergence difficulty and diversity difficulty. The former refers to the introduction of infeasible obstacles when approaching the CPF, and the position and the distance variables are related to each other. The diversity difficulty constraints will limit the feasible optimal regions, with the result that different shapes of CPFs can be obtained. Among these problems, CF1 and CF2 are the simplest, because they both have regular CPFs, while the CPFs of the remaining problems are irregular, causing difficulties in solving them.
Picard and Schiffmann [141] used the design of electromechanical actuators to develop realistic benchmark CMOPs. Twenty problems are derived with four different constraint levels and up to five objectives. These problems are representative of a variety of mechanical design applications with discontinuities in the objective space.
On the basis of the unconstrained test set ZDT [144], the constrained test suite CZDT was proposed in [54]. It contains five equality constrained test problems, and each test problem has two objectives. In addition, the only feasible solution sets of these problems are the global Pareto-optimal fronts.
Based on the DTLZ test problems [143], eight equality constrained optimization problems were proposed in the Eq-DTLZ test suite [142]. These test problems are adjustable in terms of the number of decision variables, objectives, and equality constraints. According to the position of the CPF and UPF of each test function in the objective space, the test problems can be divided into the following four types [18].
1) Type I: The UPF is also feasible, as shown in Fig. 4(a), and it is completely overlapped with CPF. 2) Type II: The UPF is partially feasible, as shown in Fig. 4(b), CPF is a part of UPF. 3) Type III: As shown in Fig. 4(c), the UPF is partly feasible and partially overlapped with the CPF. 4) Type IV: As seen in Fig. 4(d), the UPF is located in the infeasible region; thus, UPF and CPF are completely separated. Obviously, for Type I and Type II problems, the objectives can be given more preference since the CPF is part of the UPF. That is, finding the UPF is extremely beneficial to obtaining CPF. For Type III problems, the balance between objectives and constraints should be considered reasonably. While for the problems that belong to Type IV, it will be better to emphasize the importance of constraints since there is no intersection between UPF and CPF. Hence, different types of problems require different configurations of preference in objectives and constraints.
According to the above classification method, various test sets are classified in Table II. In addition, the characteristics of these test suites are presented in Table III. V. PERFORMANCE ANALYSIS In this section, the performance of different CHTs and CMOEAs on different types of problems are investigated. MW and LIR-CMOP test suites are employed because they are widely used and can be classified into four types of problems, as shown in Table II. All experiments are carried out on the PlatEMO platform [145], and the parameter settings for experiments are given in the supplementary material to retain space.

A. Comparison of Different CHTs
Five popular CHTs are chosen for performance comparison. They are self-adaptive penalty method (SP) [37], CDP [45], ε constrained method [109], SR [71], and MO [77]. For a fair comparison, these CHTs are embedded into the widely used NSGA-II framework to form CMOEAs. The detailed experimental results are shown in Table S  percentages of the best and second-best IGD results of all CHTs on each type of problems.
From Table IV, for Type I CMOPs, SP outperforms other CHTs since it performs either best or second-best on 83.3% of this type of problems. For this kind of problems, the UPF completely overlap with CPF, so the information of objectives and constraints on the UPF is very desirable. Since SP prefers individuals with better objectives and lower constraint violation, it encourages the population to overcome the infeasible obstacle to reach the CPF. During the evolution, the degree of punishment on infeasible solutions is increasing; thus, more efforts can be focused on approaching the CPF. Hence, the SP method has an advantage on solving Type I CMOPs.
For Type II CMOPs, the ε constrained method and SR achieve the better performance than other CHTs. ε constrained method takes some infeasible solutions that satisfy the ε-level as pseudo-feasible solutions. By doing this, these pseudofeasible solutions can be evolved toward the UPF. The CPF of Type II CMOPs is part of the UPFs, thus searching for the UPF is very beneficial for finding the CPF. As the ε value gradually decreases, the population is capable of gradually approaching the CPF. In addition, SR has a certain probability to evolve the population only considering objectives, thus some infeasible solutions will be utilized to cover the UPF. This property is beneficial for searching for CPF in this type of CMOPs.
For Type III and Type IV problems, MO obviously outperforms other CHTs. This method regards the constraint violation as an additional objective when comparing infeasible solutions, and some infeasible solutions close to the boundary of the feasible region will be preserved. Consequently, the population can simultaneously approach the UPF and the boundary of the feasible region. Since the CPF is part of the boundary of feasible region in Type III and Type IV CMOPs, MO has the ability to approximate the CPF from both the feasible and infeasible sides. Fig. S-1 of the supplementary material shows the convergence curves of IGD with different CHTs on MW2, MW6, MW7, and LIR-CMOP7. It is worth noting that these four problems belong to Types I-IV, respectively. We can see that different CHTs are suitable to different problems, which are consistent with the above observations. In addition, Fig. 5 shows the rankings of all CHTs on MW and LIR-CMOP test sets obtained by Friedman's test [146]. Overall, the  " INDICATES THAT THE TEST SET HAS THIS FEATURE; "×" MEANS THE TEST SET DOES  NOT HAVE THIS FEATURE; AND "⊗" INDICATES THAT SOME FUNCTIONS IN THE TEST SET HAVE THIS FEATURE   TABLE IV  PERCENTAGE OF THE BEST  performance of MO achieves the first ranking, followed by the ε constrained method.

B. Comparison of Different CMOEAs
In order to analyze the performance of CMOEAs on different types of problems, five representative CMOEAs are selected, which are BiCo [90], C-TAEA [88], CCMO [87], MOEA/D-DAE [70], and PPS-MOEA/D [19]. These CMOEAs have achieved exceedingly competitive performance on CMOPs and are widely used in comparison experiments. The detailed results are shown in Table S-2 of the supplementary material. Table V summarizes percentages of the best and second-best IGD results obtained by different algorithms for each type of problems.
From the results in Table V and Table S-2, BiCo exhibits competitive performance both on Type I and Type IV CMOPs. The reasons are as follows. BiCo includes one main population and one archive population. The main population uses the CDP to approach the CPF from the feasible region side. While in the archive population, the nondominated infeasible solutions are found and an angle-based selection scheme is developed to guide the population toward the CPF also from the infeasible side. In addition, for Type IV problems, the solutions between UPF and CPF are very effective to find CPF. Therefore, BiCo is capable of utilizing these solutions stored in the archive population to obtain good performance.
C-TAEA performs the second-best performance on Type II problems, this superiority mainly derives from MW1, MW5, and MW6. The CPFs of these MW problems are discontinuous, and are part of the UPFs. These problems require more diversity of population and information of objective function to find more CPF fragments. The CA in C-TAEA is used to find CPF, and the DA ignores the constraints, and more objective function information is considered. In addition, DA is also used to find areas that have not been searched by CA, which increases the diversity of population and makes more CPF fragments be achieved.
CCMO performs better than other algorithms on the first three types of problems, especially Type II and Type III problems. This is mainly because the UPFs and CPFs of these problems have a relatively high degree of overlap. CCMO uses a dual-population mechanism, in which population 1 focuses on approaching the CPF by both considering the objectives and constraints, while population 2 aims at searching the UPF by completely ignoring the constraints. Meanwhile, these two populations exchange information in the environment selection. By doing this, population 2 can help population 1 maintain population diversity and pass through the infeasible barrier, which encourages the search for CPF on the first three types of problems. While the superiority of CCMO degrades on Type IV problems with the completely detached UPF and MOEA/D-DAE mainly performs well on some LIR-CMOP problems, because its detect-and-escape strategy once detects the population falls into the local feasible or infeasible regions, the constraints will be relaxed by using the improved ε constrained method, to make the population escape from the local regions, and then gradually move closer to the CPF.
PPS-MOEA/D surpasses other competitors on Type IV problems. The evolution process of PPS-MOEA/D involves push and pull stages. The former prefers crossing the infeasible regions in front of the UPF, so as to prevent the population from falling into the local feasible regions. The latter focuses on driving the solutions obtained by the push stage toward the CPF. By doing this, the individuals will be gradually pulled back from the UPF to CPF. Fig. S-2 in the supplementary material shows the convergence curves of IGD with different CMOEAs on MW2, MW6, MW7, and LIR-CMOP7. In general, the convergence speed of CCMO and BiCo has a narrow advantage since their second population can help the main population quickly converge to the CPF. In addition, Fig. 6 presents the rankings of CMOEAs on MW and LIR-CMOP test suites obtained by Friedman's test [146]. Obviously, CCMO performs the best, and the performance rankings of other algorithms are comparable.

VI. APPLICATIONS
Many practical problems can be regarded as CMOPs, and the CMOEAs have achieved great success in a variety of fields. In this section, the applications of evolutionary constrained multiobjective optimization in some representative areas are introduced. The detailed information of the real-world problems solved by CMOEAs is summarized in Table VI.
Fan et al. [58] combined ACDP with the MOEA/D framework to solve the I-beam optimization problem. Qian et al. [109] used SADE-εDE to solve the welded beam design problem, speed reducer design problem, and disc brake design problem. Peng et al. [29] also applied their algorithm to these problems. Xu et al. [33] used IMDE to address four common engineering design problems, i.e., the design of welding beam, disk brake design problem, speed reducer design problem, and the spring design problem. Sadollah et al. [125] solved six practical engineering problems by using the MOWCA algorithm. Yang et al. [66] used the MODE-SaE algorithm to deal with the steady-state bi-source compressed-air pipeline optimization problem. Tian et al. [95] used the CMOEA-MS algorithm to solve the car side impact problem, vibration platform design problem, and water resource problem. Fan et al. [19] proposed the PPS framework and achieved good results in solving the optimization problem of robot gripper. Maminov and Posypkin [174] constructed the Pareto boundary of the CMOP to solve the robot design optimization problem. Xiang et al. [21] applied the CIC-MOEA/D algorithm to the search-based software engineering problem. Singh et al. [156] used CMOEA based on the surrogate model (ECMO) to solve antenna design problems. Zhang et al. [154] proposed a MOPSO and then applied it to solve the problem of multimodel switching controller based on the switched system model. The turning of composite materials needs to consider the feed speed, spindle speed, and other parameters, which can be regarded as a CMOP. Gadagi and Adake [157] solved this problem based on the genetic algorithm and particle swarm optimization algorithm. For the optimal design problem of reinforced concrete structures, Afshari et al. [158] tested a variety of multiobjective optimization algorithms, compared their performance, and analyzed the effects of different algorithms on this problem.

B. Scheduling Optimization Problems
The problems of dispatching optimization have emerged unremittingly in recent years, including the urban bus scheduling problem [3], CEED problem [159], the optimal reactive power dispatch (ORPD) problem [162], and so on.
The urban bus scheduling problem is a common practical problem, Ma and Wang [25] employed ShiP to solve this problem and achieved good results. In the economic emission scheduling problem of the power system, it is necessary to minimize the cost and emission while satisfying various constraints. Chen et al. [160] designed a constrained multiobjective population extremal optimization algorithm (CMOPEO-EED) to solve this kind of problems. El-Shorbagy and Mousa [161] proposed an algorithm, which combined the dominance principle, a clustering method, and a constraint repair method to solve the economic emission scheduling problem. Yu et al. [110] designed a new DE mutation strategy for solving the CEED problem. In order to solve the problem of reactive power optimal dispatch in power system operation, Mohseni-Bonab et al. [162] used the ε constrained method and a fuzzy satisfaction method to select the optimal compromise solutions.

C. Path Planning Problems
The optimization of the path problem can be used in all aspects of scientific and technological life, such as vehicle routing problems with time windows [163], robot trajectory planning [164], constrained multiobjective solid traveling salesman problems [165], and unmanned aerial vehicle (UAV) planning [166].
Vehicle routing problem with time windows is a widely studied combinatorial optimization problem with a complex decision space and strict constraints. Tian et al. [87] used the CCMO framework to solve the vehicle routing problems with time windows. In order to solve the trajectory optimization problem of a motor-driven parallel robot, the constrained multiobjective genetic algorithm (MOGA) was used by Chen and Pham [164]. The constrained multiobjective entity traveling salesman problems in rough, fuzzy rough, and random rough environments is a particularly challenging problem, Maity et al. [165] proposed a rough MOGA (R-MOGA) to solve this problem. UAV has been developing rapidly in recent years. Ramirez-Atencia and Camacho [166] combined a MOEA with a constraint satisfaction problem model (NSGA-II-CSP) to solve the UAV task planning problem.

D. Resource Optimization Problems
In order to make better use of resources or avoid the waste of resources, the optimization of resources is important. At present, it mainly includes water resources optimization [167], water distribution network optimization [168], water resource management problem [169], groundwater remediation design [170], wind farm layout optimization [171], thermocline filling layer thermal energy storage [172], and test resource allocation optimization [173].
Li et al. [88] proposed the C-TAEA algorithm and applied it to optimize the water distribution network. Wang and Xu [27] proposed a method, named C-AnEA, to solve the five-objective water resource management problem. Ouyang et al. [170] proposed an adaptive CMOEA based on the multialgorithm (AMALGAM), then employed it to solve the multiobjective optimization design of groundwater remediation in nonaqueous liquid contaminated sites. Sorkhabi et al. [171] combined NSGA-II with a penalty function method and constrained programming (CHCP) to solve the constrained multiobjective wind farm layout optimization problems. Marti et al. [172] used CMOEA to optimize the efficiency and the material cost of the thermocline packed bed energy storage system with air as heat transfer fluid. Chen and Zhou [175] transformed a complex CMOP into several simple subtasks, and then applied this thought to the constrained multiobjective portfolio problem. Su et al. [173] first proposed a multiobjective testing resource allocation problem model with predetermined reliability (ECHTs), which allocated the limited test time to each module, and optimized the system reliability, test cost, and test time effectively. Yu and Lu [98] suggested a corner point-based algorithm to cope with the problem of resource arrangement in emergency management that is built as a CMOP.

VII. FUTURE DIRECTIONS
Despite the numerous CMOPs that have been tackled by evolutionary algorithms, this is a research area with several potential future research topics. Some of them are discussed as follows.
1) Handling Dynamical CMOPs: In various real-world scenarios, the dynamic environment is involved in many CMOPs. That is, the objectives and/or constraints of CMOPs may change over time, such as the optimal control problems [176], the optimization of fluid catalytic cracking operation [177], and the operational indices optimization of beneficiation process [178]. These problems are called dynamic CMOPs (DCMOPs), and they are more challenging due to the simultaneous presence of objectives, constraints, and dynamism. For a DCMOP, its CPF will be changed due to the dynamic environment. As a result, the CPFs in different environments may differ in terms of position, shape, and geometry. The goal of solving DCMOPs is to rapidly obtain all CPFs in different environments. Hence, the main challenge in solving such problems is that the optimization algorithm needs to track changes in the environment and adaptively modify the search strategy when changes are detected. However, only a few efforts have been focused on developing the evolutionary algorithms to solve DCMOPs [179], [180]. Therefore, more efficient algorithms are expected to be developed for DCMOPs. 2) Handling Multimodal CMOPs: Multimodal optimization is a kind of problems, in which multiple Paretooptimal solutions correspond to the same objective, and all Pareto-optimal solutions are desired [181], [182]. Multimodal CMOPs exist in many practical applications, such as the speed reducer model [183] and the pressure vessel model [184]. Solving this type of problems requires a strong diversity to find multiple Paretooptimal solutions. Despite many MOEAs being customized to solve multimodal multiobjective optimization problems [185], it is difficult to address multimodal CMOPs, since constraints complicate the decision space.
To efficiently solve multimodal CMOPs, some attempts should be focused on the diversity preservation mechanisms in the decision space to obtain multiple feasible Pareto-optimal solutions simultaneously. 3) Handling Computationally Expensive CMOPs: Some practical optimization problems require large computing resources and an amount of time for function evaluation [186], [187], such as computational electromagnetics [188] and computational fluid dynamics [189]. They can consume several minutes to many hours for conducting one evaluation of a candidate solution. Such problems are computationally expensive. In order to solve these problems, surrogates, the computationally cheap models, are used to reduce the function evaluations. However, the existing surrogate model-based algorithms have not been well expanded to the computationally expensive CMOPs, which are often faced in real-world applications [120]. Therefore, more effective surrogate model-based methods should be developed and combined with the existing CMOEAs to solve computationally expensive CMOPs. 4) Handling Large-Scale CMOPs: Despite the existing CMOEAs achieving promising performance in solving general CMOPs with small scale decision variables and objectives, as the decision variables or objectives are scaled up, their effectiveness may dramatically deteriorate since the curse of dimensionality [190]. Large-scale CMOPs are also common, such as the multiobjective vehicle routing problem [191], time-varying ratio error estimation problem [192], and configuring software optimization [193]. To address these problems, the mechanisms [194] customized for large-scale problems can be embedded into the existing CMOEAs framework. 5) Handling Robust CMOPs: The objectives in many practical problems are easily affected by uncertain environmental factors [195], [196]. This type of problem is called the robust optimization problem [197], [198]. In robust CMOPs, both constraints and objectives may be affected by external factors. When solving such problems, it is necessary to obtain the solutions that can quickly return to the initial position after being disturbed by influencing factors. This kind of robust solutions is more acceptable than the Pareto-optimal solutions. However, the robust mechanism is more difficult to design due to the complicated uncertain environment. Therefore, it is expected to design reasonable robustness mechanisms and combine them with constraint processing technology to solve robust CMOPs. 6) Problem-Type-Guided Evolution: As mentioned in Section IV, there are different types in the relationship between UPF and CPF. If the relationship can be learned and then used to guide the subsequent evolution operators, the searchability will be enhanced efficiently. To be specific, for different types of problems, UPF has different guiding effects on CPF. For example, when the problem belongs to Type I or Type II, the information on the UPF can be directly used by the CPF. In this case, the UPF information can help the population find CPF through knowledge transfer [199]. Hence, developing the strategy to obtain the type of problems will be beneficial to the solving of CMOPs. 7) Algorithm Recommendation System: According to the principle of no free lunch [200], no algorithm can achieve good results on all types of problems. It is expected that each algorithm is allocated to its own suited problems. Hence, various CMOEAs can be combined with machine learning [201], deep learning [202], and other techniques to extract the characteristics of different problems, and then recommend the matched algorithms according to the characteristics of the CMOPs [203]. 8) Landscape-Based Evolution Operators: Different fitness landscapes and constraint landscapes [204], [205] exist in different problems. Landscapes have different characteristics and can reflect the internal properties of problems, so the required evolution operators are also different. Detecting the fitness landscape of a problem, the algorithm applied to it can be selected according to the landscape, so the landscape information can be considered and used to guide the design of evolution operators in CMOEAs [206].

VIII. CONCLUSION
This article reviewed the research of evolutionary constrained multiobjective optimization, covering the basic concepts, existing algorithms, benchmark test functions, applications, and future research directions. First, the relevant theoretical background and some concepts have been introduced. Second, this article has reviewed the existing CMOEAs and divided them into seven categories according to their CHTs and internal mechanisms. Third, some existing test function suites have been introduced and classified in detail, and the advantages and disadvantages of each test suite are summarized. Fourth, some popular applications of evolutionary constrained multiobjective optimization have been introduced and summarized. Fifth, the performance of different CHTs and CMOEAs on different types of problems was investigated. Finally, some future research directions have been discussed. It is hoped that the work presented in this article will help researchers to become familiar with evolutionary constrained multiobjective optimization, and promote the future development of this research direction.