{'fraction_non_alphanumeric': 0.03160642570281125, 'fraction_numerical': 0.008493975903614458, 'mean_word_length': 4.867915635678096, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Characterization of Constrained Continuous Multiobjective Optimization Problems: A Feature Space Perspective

Despite the increasing interest in constrained multiobjective optimization in recent years, constrained multiobjective optimization problems (CMOPs) are still unsatisfactory understood and characterized. For this reason, the selection of appropriate CMOPs for benchmarking is difficult and lacks a formal background. We address this issue by extending landscape analysis to constrained multiobjective optimization. By employing four exploratory landscape analysis techniques, we propose 29 landscape features (of which 19 are novel) to characterize CMOPs. These landscape features are then used to compare eight frequently used artificial test suites against a recently proposed suite consisting of real-world problems based on physical models. The experimental results reveal that the artificial test problems fail to adequately represent some realistic characteristics, such as strong negative correlation between the objectives and constraints. Moreover, our findings show that all the studied artificial test suites have advantages and limitations, and that no"perfect"suite exists. Benchmark designers can use the obtained results to select or generate appropriate CMOP instances based on the characteristics they want to explore.


Introduction
Many real-world continuous optimization problems involve multiple conflicting objectives that need to be optimized simultaneously and constraints that need to be respected [1]. Such optimization problems are called constrained multiobjective optimization problems (CMOPs) and have recently gained much interest in the optimization community. Indeed, many novel approaches for constraint handling and new benchmark suites of CMOPs have been proposed in the last three years [2].
For this reason, never before has there been such a strong need for problem analysis that would inform the selection of CMOPs for benchmarking multiobjective optimization algorithms (MOA). Currently, the benchmark suites are still insufficiently understood, and their characteristics remain unrevealed. In addition to characteristics describing the number of variables, objectives and constraints, type of objective and constraint functions, and geometric properties of Pareto front shapes, there are only few and limited techniques proposed to explore the search space of CMOPs, especially from the field of exploratory landscape analysis (ELA) [3,4,5].
There are only limited studies on ELA in the multiobjective combinatorial context [6,7,8], an initial study in continuous multiobjective optimization [9], and some first attempts to visualize biobjective continuous problem landscapes with up to three variables [10,4,11,12]. The situation is even worse in constrained optimization, with a single preliminary study on the characterization of constrained continuous single-objective optimization problems [5] and an attempt to extend this work to multiobjective optimization [3]. The latter work studied five ELA features, including the feasibility ratio, disjointedness of feasible regions, and relationships between the objectives and constraints. Although the study considers ELA techniques for multiple objectives and constraints simultaneously, many other crucial aspects still need to be considered. For example, recently proposed constraint handling techniques assume that CMOP violation landscapes consist of several "suboptimal" subregions-what we will call violation multimodality-that severely aggravate algorithm performance [13]. While the multimodality for unconstrained multiobjective optimization has been extensively studied in [14], the violation multimodality has not yet been studied from the ELA perspective [2].
Due to the lack of appropriate techniques to characterize CMOPs, preparing a sound and well-designed experimental setup for constrained multiobjective optimization is a challenging task. A poorly designed benchmark might lead to inadequate conclusions about CMOP landscapes and MOA performance. Two steps are required to resolve this issue: 1) providing new ELA techniques specialized for CMOPs and 2) using them to reveal the relationship between CMOP characteristics and MOA performance [15]. In this paper, we address step 1) by extending various ELA techniques to constrained multiobjective optimization, and provide a thorough analysis of the available test suites of CMOPs. In contrast, step 2) is currently under investigation and not a part of this paper due to limited space.
The contributions of this work are as follows: • We generalize the concept of a fitness landscape to CMOPs and introduce two mathematical definitions to formally describe violation multimodality in constrained multiobjective optimization.
• We use various ELA techniques including space-filling design, information content, random and adaptive walks to derive 29 ELA features, of which 19 are novel. These features are then used to characterize violation multimodality and smoothness and measure the conflicting nature between the objectives and constraints. In contrast to the related work, we also discuss the sensitivity and scalability of the proposed techniques.

Theoretical background
This section provides the theoretical background for this study. After the definitions of the constrained multiobjective optimization problem and constraint violation, we present the concept of the fitness and violation landscape. Finally, we discuss the topological property of connectedness.

Constrained multiobjective optimization problems
A CMOP is, without loss of generality, formulated as: In continuous optimization, the equality constraints are usually reformulated into inequality constraints as follows: where η > 0 is a user-defined tolerance value to relax equality constraints (in our study it is set to 10 −4 ). Using this definition, all constraints can be treated as inequality constraints. A solution x is said to be feasible if it satisfies all the constraints g i (x) ≤ 0 for i = 1, . . . , I + J.
One of the most important concepts in constrained optimization is the notion of the constraint violation. For a single constraint g i it is defined as v i (x) = max(0, g i (x)). For all constraints together it is combined as into the overall constraint violation. A solution x is feasible iff its overall constraint violation equals zero (v(x) = 0). Note that other definitions for overall constraint violation exist, and their use would impact the analysis performed in this study. However, the proposed definition for the overall constraint violation is by far the most commonly adopted in constrained optimization [2], and as such, it represents the most appropriate choice. A feasible solution x ∈ S dominates another solution y ∈ S iff f m (x) ≤ f m (y) for all 1 ≤ m ≤ M and f m (x) < f m (y) for at least one 1 ≤ m ≤ M . Additionally, a solution x * ∈ S is Pareto optimal if there is no solution x ∈ S such that it dominates x * .
The set of all feasible solutions is called the feasible region and is denoted by F = {x ∈ S | v(x) = 0}. All nondominated feasible solutions represent a Pareto-optimal set, S o . The image of the Pareto-optimal set in the objective space is the Pareto front and is denoted here by

Fitness and violation landscapes
The concept of a fitness landscape was formally defined in [17] by identifying the following three key elements: The distance metric d is used to quantify the similarity between solutions in the search space and to define neighborhoods. In continuous optimization, the Euclidean distance is often used. Close solutions are those whose Euclidean distance is within a predefined threshold, δ. Formally, the solution y ∈ S is part of the neighborhood of the solution x ∈ S, denoted by N (x, δ), iff d(x, y) ≤ δ.
In [5], an adaptation of the fitness landscape was proposed for constrained optimization problems by replacing the objective function with the overall constraint violation function. The resulting idea is known as the violation landscape.

Connectedness
A set X ∈ S of a topological space is said to be disconnected if it can be represented as the union of two or more disjoint non-empty open subsets from X. Otherwise, X is said to be connected.
The maximal connected subsets (with respect to the inclusion order) of a non-empty space are called connected components. The components of any topological space X form a partition of X-they are disjoint, non-empty, and their union is the entire space.

Theoretical advances
In this section, we overview theoretical advances in landscape analysis for CMOPs. After defining constrained multiobjective problem landscapes, we introduce feasible components and basins of attraction in the violation landscape. Finally, we conclude with three examples depicting these novel concepts.

Constrained multiobjective problem landscapes
Following the definitions of fitness and violation landscapes from Section 2.2, we generalize these concepts to CMOPs by including multiple objective functions and the overall constraint violation function simultaneously. A constrained multiobjective problem landscape is defined with the following four elements: • v : S → R is the overall constraint violation function, The notion of a neighborhood can be applied to CMOP landscapes without any modification. For example, a solution x * ∈ S is a local Pareto-optimal solution if it is feasible and there exists a δ > 0 such that no feasible solution from the neighborhood N (x * , δ) dominates x * . The set of all local Pareto-optimal solutions will be denoted by S l . It is obvious that S o ⊆ S l .
In [13], two topological structures concerning violation multimodality were identified that might cause premature algorithm stagnation: feasible subregions and areas that locally violate the constraints the least. The feasible region can consist of multiple disconnected subregions. The more such subregions there are, the harder it is for the algorithm to find a good approximation for the whole Pareto-optimal set. Moreover, the search space of a CMOP can have many areas that locally violate the constraints the least and are infeasible. Again, many such areas might aggravate the performance of an algorithm since it can get stuck in the infeasible region. Although these concepts are essential, they are not formally defined, which can lead to ambiguous interpretations. To overcome this situation, we provide a rigorous mathematical formulation and, in this way, facilitate their further consideration.

Feasible components
A connected component F ⊆ F of the feasible region is called a feasible component. Note that all feasible components form a partition of the feasible region. With feasible components, we provide a formal definition of feasible subregions. A feasible component can be treated as "suboptimal" when it contains no Pareto-optimal solutions. Obviously, if a problem landscape contains several "suboptimal" feasible components, it might be harder for an optimizer to find good approximations for all Pareto-optimal solutions.

Basins of attraction in the violation landscape
Basins of attraction in the fitness landscape are well-known objects widely used to characterize unconstrained single-objective optimization problems. In this paper, we reuse this concept in the violation landscape. However, we need to introduce some crucial notions first.
Similarly to the local Pareto-optimal solution, we can define a local "optimal" solution in the violation landscape-a solution that locally violates the constraints the least. Formally, The set of all local minimum-violation solutions will be denoted by F l and of all minimum-violation solutions by F m . It is obvious that F m ⊆ F l and if there exist feasible solutions in S, then F m = F .
A connected component M ⊆ F l is called a local minimum-violation component and it provides a formal definition for areas that locally violate the constraints the least. All the local minimum-violation components form a partition of F l . Now, we can finally define a basin of attraction. Let us consider an abstract local search procedure as a mapping from the search space to the set of local minimum-violation solutions, µ : S → F l , such that µ(x) = x for all x ∈ F l . Then, a basin of attraction of a local minimumviolation component M and local search µ is a subset of S in which µ converges towards a solution from M, i.e., B(M) = {x ∈ S | µ(x) ∈ M}. In this case, the local minimum-violation component M is said to be an attractor of the basin B(M). If there is only one basin in S, then the corresponding violation landscape is said to be unimodal. Otherwise, it is multimodal.
Each feasible component is a subset of one of the local minimum-violation components. The reverse is, in general, not true. However, if a local minimum-violation component contains feasible solutions, then it contains at least one feasible component as well. Consequently, one can discover feasible components by first identifying local minimum-violation components that contain feasible solutions. Figure 1 illustrates problem landscapes and violation landscapes for problems C2-DTLZ2 [18], MW6 [1], DAS-CMOP1 [19] and MW7 [1] with two variables and two objectives. All figures were obtained by evaluating a grid of 501 × 501 solutions from the search space. The first column shows the dominance ratio for each grid solution expressed as a proportion of grid solutions that dominate it [10]. Those solutions that are not dominated by any grid solution represent approximations of the Pareto-optimal solutions and are shown in black in the plots. The second column depicts the violation landscape in terms of constraint violation values with feasible components presented in white. The third column shows the problem landscape by visualizing the dominance ratio of the feasible regions in blue hues and the infeasible regions in pink. Again, black denotes approximations of Pareto-optimal solutions.

Examples
The Pareto-optimal set of C2-DTLZ2 without considering constraints is the horizontal line L : x 2 = 0.5 (Figure 1a). Similarly holds for MW6 (Figure 1d), while the Pareto-optimal sets of DAS-CMOP1 and MW7 without constraints are curves in the search space (Figures 1g and 1j).
As we can see from Figure 1b, C2-DTLZ2 has three feasible components that are also the only three local minimum-violation components. The corresponding basins are separated with two vertical lines L 1 : x 1 = 0.25 and L 2 : x 1 = 0.75. All the feasible components and basins intersect with the Pareto-optimal set ( Figure 1c). In other words, each basin attracts the search towards nondominated feasible solutions, making this problem easy for constraint handling.
On the other hand, MW6 consists of 35 feasible components and 72 basins ( Figure 1e). Between two adjacent horizontal components, there is an additional basin that contains no feasible solutions. It is also interesting that not all the feasible components contain Paretooptimal solutions (Figure 1f). This makes the MW6 problem harder to handle, especially for separation-based constraint handling techniques that separately handle feasible and infeasible solutions and strictly favor feasible solutions over the infeasible ones. The Pareto-optimal set is a subset of the Pareto-optimal set corresponding to MW6 without considering constraints.
The DAS-CMOP1 problem has three feasible components: a small one in the upper left corner, a medium-sized one in the upper right corner, and a large one in the middle. In contrast to C2-DTLZ2 and MW6, the feasible components have irregular shapes. The feasible components are also the only local minimum-violation components. The Pareto-optimal set is a subset of the Pareto-optimal set corresponding to DAS-CMOP1 without considering constraints.
In contrast, the Pareto-optimal set of the MW7 problem consists of two disconnected parts of the Pareto-optimal set belonging to the MW7 problem without constraints and parts of the boundary between the feasible and infeasible regions ( Figure 1l). This problem has one feasible component and also only one basin of attraction ( Figure 1k).
To identify and count the basins of attraction, we use the truncated Newton method [20] as a local search procedure µ. This method approximates the Newton's direction by approximately solving the Newton equation using an iterative technique. In this study, the conjugate gradient technique is used as the iterative solver. Note that any other deterministic local search could be used instead and possibly result in a different number of basins. Nevertheless, the truncated Newton method proved to be efficient and reliable for identifying basins of attraction.  Figure 1: Plots of the dominance ratio (the first column), the violation landscape (the second column), and the problem landscape (the third column) of three benchmark problems. The approximation of the Pareto-optimal set is depicted in black in dominance ratio and problem landscape plots.

Proposed ELA features
The ELA features presented in this study can be, based on the methodology used to derive them, categorized into four groups: space-filling design, information content, random walk and adaptive walk features. They are all summarized in Table 1.
Size of the "optimal" basin New a For the first time applied in multiobjective optimization in our work. b For the first time applied on the violation landscape in our work.

Space-filling design features
The space-filling design approach employed in this paper can be described with the following three steps: 1. Generate an initial sample X S ⊆ S of solutions following a selected space-filling design and derive feasible solutions.
2. Cluster feasible solutions to obtain approximations for feasible components.
3. Find nondominated solutions among the obtained feasible solutions.
Following this approach, various features can be derived to characterize feasible components. The first feature of this group is the number of feasible components, N F , and it is estimated as the total number of the obtained clusters. Next, the proportion of solutions from X S that form a particular cluster is used to measure the size of the corresponding feasible component. Three features were derived with this respect: the smallest feasible component, F min , the median feasible component, F med , and the largest feasible component, F max .
The following two features considered within this group express the "optimality" of feasible components. The first feature reflects the "optimality" of the largest feasible component, O(F max ), and it is expressed as the proportion of nondominated feasible solutions in the largest cluster. The second feature is the size of the "optimal" feasible component, F opt . It is approximated as the proportion of all solutions from the initial sample that form a cluster with the maximum number of nondominated feasible solutions.
The proportion of feasible solutions in the initial sample is used to estimate the feasibility ratio, ρ F , of the corresponding CMOP.
In addition, two features express the correlations between the objectives and constraints. They are the minimum and maximum correlations among objectives and the overall constraint violation and are denoted by ρ min and ρ max , respectively.
Finally, the proportion of Pareto-optimal solutions on the feasible region boundary is approximated as the proportion of nondominated solutions on the boundary of one of the obtained clusters and is denoted by ρ ∂So .

Information content features
The information content approach was originally used to analyze a sequence of fitness values obtained by sorting a sample of solutions and quantifies the "information content" of a fitness landscape (smoothness, ruggedness, and neutrality). We extend this approach to the violation landscape. Instead of analyzing fitness sequences, we analyze sequences of overall constraint violation values.
The approach follows the description provided in [21] and it is here presented only briefly. The interested reader can find a detailed explanation in the original paper. First, a sample of solutions X I = {x (1) , . . . , x (n) } is generated following a space-filling design. The generated solutions from X I are sorted in a sequence as follows: The first solution is selected at random, while each subsequent solution is selected as the nearest solution to the current one. Then, a new sequence, φ(λ) = {φ (1) (λ), . . . , φ (n−1) (λ)}, is generated from the original one following the rule where . In more detail, the symbol indicates a decrease of the overall constraint violation, → a neutral area within the threshold λ, and an increase of the overall constraint violation. In particular, two consecutive symbols compose a block, ab, where a, b ∈ { , →, }, which represents an object in the violation landscape, i.e., a slope, peak, etc. For example, indicates a constant increase of the overall constraint violation, while indicates a bottom. In particular, two identical consecutive symbols represent a smooth area in the violation landscape, while two different consecutive symbols indicate a rugged area.
Finally, the information content of X I is defined as the Shannon entropy of the distribution for the blocks ab: where p ab is the probability of finding one of possible blocks ab, where a = b, within the sequence φ(λ). The logarithm base is six as this is the number of different blocks where a = b. The main intuition behind information content is that a large value of H(λ) is an indication of perpetual changes in the violation landscape and therefore its ruggedness. Another important measure is the partial information content defined as M (λ) = |φ (λ)| /(n − 1), where φ (λ) is a sequence derived from φ(λ) by removing all → and repeated symbols.
From H(λ) and M (λ), various features can be derived representing smoothness of the violation landscape. In this work, we study three features often used in the literature [15]: the maximum information content H max = max λ {H(λ)}, the settling sensitivity ε S = log 10 (min λ {λ | H(λ) < 0.05}), and the initial partial information

Random walk features
A random walk through the problem search space is a well-known and often used ELA technique in single-objective optimization. A walk is an ordered sequence of solutions X R = {x (1) , . . . , x (n) } such that x i ∈ N (x (i−1) , δ) for all i ∈ {2, . . . , n}. During a random walk, there is no particular criterion to select the neighboring solution at each step, i.e., a random neighbor is selected [8].
For constrained single-objective optimization, a ratio of feasible boundary crossings, ρ ∂F , was introduced in [5]. This feature quantifies the number of boundary crossings from feasible to infeasible space and vice versa, as encountered by a random walk through the search space. Given a sequence X R generated from a random walk, a binary sequence is derived B = {b (1) , . . . , b (n) } such that b i = 0 iff x i is a feasible solution. The feature ρ ∂F is defined as the proportion of steps in the random walk that cross the boundary of the feasible region: Since the ratio of feasible boundary crossings does not consider the objective values, it is a feature characterizing violation landscapes. Therefore, it can be used in multiobjective optimization without any modifications [3].
To provide a more robust measure, [5] defines ρ ∂F as the average value over N independent random walks. However, this aggregate value provides no local information at all. To avoid this shortcoming, we also consider the maximum and minimum values of ρ ∂F observed in N random walks. The resulting features are: the minimum ratio of feasible boundary crossings (ρ ∂F ) min , the median ratio of feasible boundary crossings (ρ ∂F ) med , and the maximum ratio of feasible boundary crossings (ρ ∂F ) max . A value of (ρ ∂F ) min close to zero indicates that there exists at least one larger feasible (or infeasible) region. Similarly, a large value of (ρ ∂F ) max indicates that there is an area in the violation landscape with a heavily dissected feasible component.

Adaptive walk features
In contrast to random walks, an improving neighbor is selected at each step of an adaptive walk, as a local search would do. In unconstrained single-objective optimization, adaptive walks are often used to derive various features to characterize basins of attraction of the fitness landscape. In this paper, we extend these approaches to the violation landscape to quantitatively describe its basins of attraction.
The approach considered in this study is similar to the technique described in [22]. However, various novel features specialized for CMOPs are introduced here for the first time. The approach consists of four steps: 1. Generate an initial sample X A ⊆ S of solutions following a selected space-filling design.
2. Perform a local search for each solution x ∈ X A to obtain local minimum-violation solutions.
3. Cluster local minimum-violation solutions to obtain approximations for local minimumviolation components and basins of attraction.
4. Find nondominated feasible solutions among the obtained local minimum-violation solutions.
Using the above approach, we can derive various features to characterize basins of attraction of the violation landscape. The resulting total number of the found clusters provides an approximation for the number of basins, N B , the first feature of this group. The proportion of all solutions from the initial sample that converge to a specific cluster is used to measure the size of the corresponding basin. The resulting features are: the smallest basin, B min , the median basin, Finally, the last two features considered in this group indicate the "optimality" of basins of attraction. The first feature is the "optimality" of the largest basin, O(B max ), expressed as the number of nondominated feasible solutions in the largest cluster. The second feature is the size of the "optimal" basin, B opt , expressed as the proportion of solutions that converge to the cluster with the maximum number of nondominated feasible solutions.

Experimental analysis
In this section, we first introduce the test suites chosen for the experiments. Next, we discuss the experimental setup, including the specific techniques and parameter settings used to derive the ELA features. The derived features are then used to evaluate the existing test suites. Additionally, we present the sensitivity analysis of the parameters used for feature extraction, and finally, we comment on the scalability of the applied space-filling design techniques. Hydro-static thrust bearing design 4 2 7
All the mentioned test suites consist of artificially designed CMOPs and are not derived from real-world applications. To overcome this weakness, a novel suite named RCM was proposed in [16]. The RCM suite collects 50 real-world optimization problems based on physical models, including problems from mechanical design, chemical engineering, power electronics, etc. The problem instances come in various dimensions and numbers of objectives. After excluding all the problems containing discrete variables and those with more than five variables, 11 RCM problems remained for our experimental analysis and are presented in Table 2. As we can see, ten of the selected RCM problems have two objectives while the reaming one has five objectives. The number of constraints varies from one to eight, and all the constraints are inequalities.
Note that an additional real-world multiobjective optimization problem suite including eight CMOPs was proposed in [28]. Nevertheless, all relevant CMOPs from this suite are also included in the RCM suite, and it was therefore omitted from this study.
The basic characteristics of the test suites used in this study are summarized in Table 3.

Experimental setup
The proposed features are demonstrated on the test suites listed in Section 5.1. In particular, bi-objective C-DTLZ and DC-DTLZ problems were considered with the default number of constraints. Additionally, a difficulty triplet of (0.5, 0.5, 0.5) was used for the DAS-CMOP suite as this is by far the most frequently used difficulty triplet in the literature. Three dimensions of the search space (2,3,5) were considered to measure the scalability of various techniques used to derive the ELA features. For each dimension, initial sample sizes were decided based on some initial experiments and were selected as the minimum number of solutions needed for the feature values to converge. The details are presented in Table 4. Note that because the real-world problems from the RCM test suite are not scalable, the parameter values used for artificial problems with five variables were used also for the four-dimensional RCM problems.
To produce the initial samples used by all the proposed techniques to derive ELA features, we used the space-filling design based on Latin hypercube sampling.
We used the density-based spatial clustering of applications with noise (DBSCAN) [29] as the clustering algorithm for the identification of feasible components and basins, as well as to discover nondominated solutions on the boundary of the feasible region. The following DBSCAN configuration was used in the experiments: distance metric was set to the Euclidean distance, the number of samples in a neighborhood for a point to be considered a core point was set to five. At the same time, the maximum distance between two solutions for one to be considered as in the neighborhood of the other, ε, was defined based on initial experiments for each dimension separately (see Table 4).
To express the relationship between the objectives and constraints (ρ min and ρ max ), the Spearman's rank correlation coefficient was used.
Following the recommendation from [21], the initial sample size X I for the information content features was set to 1 000D.
Simple random walks [30] were employed for calculating features representing ratios of feasible boundary crossings. The number of steps X R was set to 10 000, each with maximum step size, δ, of 1% of the range of the domain for each test problem [5]. Thirty independent runs were conducted to obtain the values for (ρ ∂F ) min , (ρ ∂F ) med and (ρ ∂F ) max .
Finally, we used the truncated Newton method [20] as a local search procedure used in adaptive walks.

Implementation details
All CMOPs and techniques to derive the ELA features were implemented in the Python programming language [31]. We used the pymoo [32] implementation for CTP, DAS-CMOP and MW, while the rest of the suites was reimplemented from scratch. Next, pyDOE [33] was used for Latin hypercube sampling, scikit-learn [34] for DBSCAN, pflacco [35] for the calculation of information content features, and SciPy [36] for the truncated Newton method and to calculate Spearman's rank correlation coefficients. The rest of the functionalities, including the simple

Figures 2 and 3 show violin plots of distributions for seven selected features
. The first column represents the distribution for the set of all the considered CMOPs, while the rest of the columns correspond to each suite separately. Each black dot represents one problem instance. Its y-axis depicts the feature value, while the x-axis has no specific meaning and is used for better visualization only. The violin plot (colored area) approximates the probability density function for the feature distribution. For example, Figure 2b shows there are more problem instances in the CF suite (third column, violin plot in orange) with ρ min ≈ 0 than those with ρ min ≈ −0.75. In addition, there are no problem instances with ρ min ≈ −1. The violin plot in light blue behind each suite corresponds to the distribution of RCM in Figure 2, and in light gray to the distribution of the set of all CMOPs in Figure 3. The idea is to show the coverage of RCM characteristics by the artificial test suites (Figure 2) and to expose the differences between a single test suite and the rest of the suites (Figure 3). In addition, Figure 4 shows two heat maps depicting the coverage of RCM characteristics by the artificial test suites, and the coverage of characteristics present in the set of all the considered CMOPs. The coverage is expressed as one minus the mean value of all distances between the feature values from the RCM suite (or the set of all CMOPs) and the nearest feature value from a selected single suite. Note that feature values need to be normalized to [0, 1] before calculating the coverage values. A coverage value near one indicates excellent coverage. This means that for each feature value from RCM (or the set of all CMOPs), there exists at least one problem from the selected single suite whose feature value is close to that one from RCM (or the set of all CMOPs). In Figure 4, a lighter color indicates that a feature value is well covered by a test suite, and vice versa. For example, in Figure 4a, N F (first row) as represented by RCM is well covered by all the artificial suites, while F max (fourth row) is poorly covered, especially by DC-DTLZ and MW.

Test suites evaluation
According to [37], a good benchmark suite should "include the difficulties that are typical of real world instances of the problem class under investigation". To evaluate whether considered artificial test suites contain the characteristics observed in real-world CMOPs, we use 11 continuous and low-dimensional problems from the RCM suite and the coverage metric proposed in Section 5.4. It is to be noted that the low number of considered real-world problems limits the generality of the conclusions that can be drawn from these results. For example, we cannot claim that some problem characteristics that are missing from these 11 problems are in fact unrealistic. Nevertheless, we can identify which characteristics of these real-world problems are not present in the artificial problems.  The results show that several features characterizing feasible components are not adequately covered by the studied artificial test suites (Figure 4a, features from N F to F opt ). Especially the recently proposed artificial test suites NCTP, DC-DTLZ, LIR-CMOP, and MW insufficiently represent problem characteristics of the studied RCM problems. For example, many of the considered RCM problems have only one feasible component, while recent artificial test problems have many feasible components. As a consequence, the RCM problems result in feasible components of any size, while, for example, the MW suite contains problems with rather small feasible components only. This shortcoming of the artificial test suites can be observed nicely for F opt (Figure 2a).
Similarly, the existing artificial test suites fail to satisfactorily represent the conflicts between the objectives and the overall constraint violation. As we can see in Figure 2b, strong negative correlations, ρ min ∈ [−1.00, −0.75], observed in the studied RCM problems, are found in the C-DTLZ suite only. We can argue that a general multiobjective optimization algorithm would most likely struggle to maintain good convergence and diversity when the objectives and the overall constraint violation are strongly negatively correlated. However, our findings show that the studied artificial test suites do not cover this particular characteristic well.
The results also indicate that the feasibility ratio and information content features are poorly covered by the artificial test suites except for CFs and LIR-CMOPs (Figure 4a, features ρ F , H max , ε s and M 0 ). While the RCM problems take a full range of feature values, the artificial test problems have often extremely small or large feasible ratios and extremely smooth or very rugged violation landscapes (Figure 2c).
On the other hand, the features related to basins of attraction are mostly covered well by the studied artificial test suites (Figure 4a, features from N B to B opt ). The suites that are not diverse enough considering the basin-related features are CF, NCTP, and DAS-CMOP. Additionally, the ratios of feasible boundary crossings are well covered by the considered artificial test suites, while ρ ∂So is very well covered by all the artificial test suites except for DAS-CMOP.
Finally, the studied artificial test suites have many basin-related characteristics that cannot be observed in the selected RCM problems. While the considered RCM problems often contain only one basin where all Pareto-optimal solutions are located, many recently proposed artificial test problems have several basins, and the Pareto-optimal solutions are spread among multiple basins (Figures 2d, 3a and 3c). This is in line with the recent development in constrained multiobjective optimization where the primary source of complexity is severe violation multimodality (e.g., Figure 1f) [2]. However, our study suggests that other relevant aspects, such as the sizes of feasible components and basins of attraction, negatively correlated objectives and the overall constraint violation, diverse ruggedness and feasibility ratios, are frequently insufficiently addressed by the literature.
An additional analysis is devoted to measuring whether the considered test suites contain a wide variety of problems with different characteristics. This is done in the same way as the comparison with RCM, except that in this case the test suites are contrasted against the set of all CMOPs. The result are shown in Figure 4b.
Interestingly, the results are not drastically different from the previous ones, indicating that RCM contains problems with various characteristics. For example, the coverage of most of the space-filling design and information content features is very similar (Figure 4b, features from N F to M 0 ). This is also true for seven features concerning basins of attraction B med , B max , (B F ) min , (B F ) med , (B F ) max , ∪B F and B opt , except for the NCTP suite which has significantly worse coverage when contrasted against the set of all CMOPs.
Within the space-filling design, the only true exception is O(F max ). The coverage for this feature is significantly worse when comparing the suites against the set of all CMOPs. This happens because RCM does not cover O(F max ) particularly well (Figure 3a), which has already been explained before. As we can see, this aspect is even more drastically manifested for O(B max ) ( Figure 3c).
Finally, the coverage of random walk features is slightly worse when the suites are compared against the set of all CMOPs (Figure 4b, features (ρ ∂F ) min , (ρ ∂F ) med and (ρ ∂F ) max ). This is true because RCM contains only problems with small values of these particular features. However, in general, this set of features is well covered. Almost all the problems have only small ratios of feasible boundary crossings except for few problems from CTP, CF and DAS-CMOP (Figure 3b). The same observation holds for the space-filling design feature N F , and for the following adaptive walk features N B , B min , v(B) med , v(B) max and v(B max ).
In summary, we can see that all the studied test suites have some advantages and some limitations. For this reason, we recommend to include CMOPs from various suites to construct a comprehensive and sound benchmark.

Sensitivity analysis
The crucial step in deriving space-filling design and adaptive walk features is the accurate identification of feasible components and basins of attraction, respectively. Once the feasible components and basins are identified, all other features can be unequivocally derived. On the other hand, if the clustering step incorrectly determines the feasible components and basins of attraction, the rest of the features is also affected. The identification of both feasible components and basins of attraction strongly depends on two parameters: initial sample size (|X S | and |X A |) and the maximum distance between two solutions ε used by the DBSCAN algorithm.
The goal is to find the smallest sample size and the smallest ε that would still yield reliable results. Small sample sizes are preferred since they require fewer function evaluations. Additionally, a small ε is necessary to not erroneously merge two or more disjoint feasible components (or basins of attraction). On the other hand, a too small ε can result in the separation of feasible components (or basins) into two or more components (or basins). For this reason, a careful selection of these parameters is necessary.
Space-filling design features are more sensitive to parameter setting than adaptive walk features. All infeasible solutions in the space-filling design approach are discarded and used only for enumeration. On the other hand, almost all initial solutions for adaptive walk features produce local-minimum violation solutions, thus facilitating the clustering procedure. For this reason, we investigate the sensitivity of space-filling design only. The presented analysis is very similar for adaptive walk features, except that smaller initial sample sizes are required in this case.
The goal of the presented analysis is to measure the parameter sensitivity for the identification of feasible components. More specifically, we study the extraction of the number of feasible components N F with respect to 1) different parameter values and 2) repetition of the experiments, since a stochastic approach is used to obtain the initial sample. To address 1), we test the space-filling design approach using all combinations of the following parameter values: |X S | ∈ {10 000, 25 000, 50 000, 100 000, 250 000, 500 000} (7) and ε ∈ {0.01, 0.02, . . . , 0.14, 0.15}.
For 2), we repeat all the experiments 30 times, each time with a different initial sample. The analysis is performed for D = 2, 3, as the correctness of the derived components can be visually checked, while this is not possible for D = 5. The sensitivity analysis is performed on two CMOPs: C2-DTLZ2 and DAS-CMOP1. On the one hand, C2-DTLZ2 represents a simple problem concerning feature extraction. It has three well-shaped feasible components that are spaced apart (Figure 1b). On the other hand, the feasible components of DAS-CMOP1 have irregular shapes, and the small one in the upper left corner is located close to the large one in the middle (Figure 1h). This is precisely what makes the identification of components hard for DBSCAN. The 2-D DAS-CMOP1 problem with two variables has three components, while the DAS-CMOP1 with three variables has five components.
The results are shown in the form of heat maps in Figure 5, where the median logarithmic difference between the exact number of feasible components and the derived number of clusters is depicted. The values are normalized with the number of feasible components and logarithmic to obtain better visualization. Smaller values in lighter colors indicate better performance, and vice versa. Specifically, white color indicates a perfect match.
As we can see from Figure 5a, for C2-DTLZ2 (D = 2) the exact number of feasible components can be obtained for all |X S | ≥ 25 000 and all ε ≥ 0.02, while for DAS-CMOP1 (D = 2) already the combination |X S | = 10 000 and ε = 0.02 is sufficient to identify the feasible components ( Figure 5c). Interestingly, for DAS-CMOP1 (D = 2) we can see that for ε ≥ 0.12 one component is recognized incorrectly. What happens is that the small component in the top left corner is merged with the large one due to the oversized ε ( Figure 6). Nevertheless, a combination of |X S | = 25 000 and ε = 0.02 proved to be the most reliable for D = 2, resulting in a perfect match for all 30 repetitions on C2-DTLZ2 and DAS-CMOP1. Moreover, Figure 6 shows that DBSCAN is really effective also when dealing with CMOPs with irregularly-shaped feasible components.
For D = 3, Figures 5b and 5d clearly show we need larger |X S | and ε to obtain the exact number of components. For both CMOPs a combination of |X S | = 100 000 and ε = 0.04 proved to be the most effective, correctly matching all the components in all 30 repetitions. Additionally, we can see that for DAS-CMOP1 two or more components are incorrectly merged for ε ≥ 0.06.  An interesting observation is that the technique seems to work better for smaller initial sample sizes when using ε ≤ 0.02. The reason behind this behavior is that DBSCAN finds no cluster at all in such cases and the number of identified feasible components is thus zero.
In conclusion, for feasible components that are spaced apart, the space-filling design and similarly adaptive walk approaches are highly reliable and insensitive to small changes in parameter values. Furthermore, the performance is not affected by the shape of the feasible components. On the other hand, for problems with feasible components that are placed close together, the techniques are more sensitive to the selection of parameter values.
Finally, a discussion on the parameter sensitivity for information content and random walks can be found in [21,5]. For this reason, we just briefly summarize the key findings. In [21], a comprehensive sensitivity analysis was performed to measure the robustness of the information content approach. Among others, the authors addressed the effect of varying the size of X I . The experimental results showed that already a sample size of 100D proved to be robust for deriving information content features. Similarly, in [5], the authors demonstrated that the random walk approach was reliable for deriving ratio of feasible boundary crossings. This was shown on a wide variety of violation landscapes by obtaining small standard deviations while repeating the experiments 30 times.

A note on scalability
For problems with up to five variables the derived ELA feature values closely coincided with actual values. In contrast, it was virtually impossible to investigate high-dimensional CMOPs since the size of the initial sample used by the space-filling design techniques grows exponentially. Indeed, a more detailed experimental analysis showed that to obtain meaningful results, we would need at least 2 000 000 solutions for problems with seven variables and 15 000 000 solutions for problems with ten variables. This is not unexpected since it is a well-known consequence of the curse of dimensionality.
In addition, the feasible components of the artificial test problems often have the shape of an ellipsoid (Figures 1b and 1e). As the dimension of the search space increases, a feasible component (ellipsoid) becomes an insignificant volume relative to that of the decision space (rectangular area). For this reason, it is very hard to obtain relevant results using space-filling design techniques for such large-scale problems.

Conclusions
In this work, we have presented a comprehensive investigation of CMOPs using landscape analysis. Firstly, we have extended the concepts of fitness and violation landscapes to constrained multiobjective optimization and proposed a mathematical formulation to describe multimodality in the violation landscape. Secondly, four ELA techniques, namely space-filling design, information content, random walks, and adaptive walks, have been adapted for CMOPs. Next, we have used the resulting features to analyze eight artificial test suites and compare them against the RCM suite, which consists of real-world problems based on physical models. In particular, the artificial test suites were assessed with respect to the representativeness of various RCM characteristics. Additionally, by contrasting the test suites against a set of all the considered CMOPs, we have also assessed whether the test suites contain problems with various characteristics. Finally, we have discussed the advantages and limitations of the existing artificial test suites as well as the sensitivity and limitations of our methodology. The experimental results show that the proposed ELA techniques are adequate to characterize constrained multiobjective optimization problem landscapes with up to five decision variables.
Our findings indicate that the artificial test problems fail to satisfactorily represent some real-world problem characteristics, e.g., strong negative correlation between the objectives and the overall constraint violation. Moreover, the predominant source of complexity in novel artificial test problems is very often violation multimodality. This can lead to biased interpretations and conclusions since many other relevant characteristics are missed. All this clearly suggests that further advances in landscape analysis for CMOPs are required to overcome these drawbacks.
This study is an important step towards enhancing the understanding of test suites developed for constrained multiobjective optimization. It has revealed that all the studied artificial test suites have certain advantages and limitations and showed that no "perfect" test suite exists.
However, we are aware that the presented work has two shortcomings. Firstly, the proposed approach to characterize violation multimodality is not sufficiently scalable, and thus only low-dimensional problems were addressed in this study. Secondly, the correlations between the proposed features and multiobjective optimization algorithm performances were not evaluated. Consequently, the effect of characteristics on constraint handling techniques could not be adequately assessed. To overcome these limitations, we will experiment with ELA techniques to characterize large-scale optimization problems. We will also comprehensively measure the algorithm performance while solving test CMOPs. Special attention will be devoted to understanding the behavior of constraint handling techniques subject to various problem characteristics.