Improved selection strategy for multi-objective evolutionary algorithms with application to water distribution optimization problems

Multi-objective evolutionary algorithms (MOEAs) have been applied to water distribution system (WDS) optimization problems for over two decades. The selection strategy is a key component of an MOEA that determines the composition of a population, and thereby the evolutionary search process, which imitates natural selection by granting fitter individuals an increasing opportunity to reproduce. This paper proposes the convex hull contribution (CHC) selection strategy for generational MOEAs (CHC Gen ) as a novel selection strategy that is based on the CHC of solutions to the Pareto front in the objective space. Numerical experiments using a general MOEA framework demonstrate that the CHC Gen selection strategy is able to outperform existing popular selection strategies (e.g., crowding distance, hypervolume contribution, and hybrid replacement selection). Moreover, it is illustrated that the CHC Gen selection strategy is able to improve the performance of existing MOEAs such as NSGA-II and GALAXY. The conclusions are based on the results of six bi-objective WDS problems.


INTRODUCTION
Over the past decades, nature-inspired computing (such as evolutionary algorithms) has become increasingly popular in civil engineering optimization applications. Since the early pioneering work applying genetic algorithms to the optimal design of engineering structures (Adeli & Cheng, 1993;Adeli & Kumar, 1995) and pipe network design (Simpson et al., 1994), their usage in the field of civil engineering has grown to a broad array of applications, including the design of large steel structures (Park & Adeli, 1997), numerous water resources applications (Nicklow et al., 2010), structural optimization of free-form space frame roof structures (Kociecki & Adeli, 2014, health monitoring (Oh et al., 2017), expert systems (Wang et al., 2018), and sensor placement (Civera et al., 2021), to mention a few. In addition to the growth of application areas, there has also been a significant growth in the development of varied algorithm types. Survey papers for nature inspired algorithms abound (e.g., Siddique & Adeli, 2015a), but algorithms of note include: ant colony optimization (a population-based meta-heuristic [PBMH] based on the foraging behavior of ants, and their ability to find the shortest path to a food source, e.g., Maier et al., 2003); the multi-point spiral dynamics algorithm (Siddique & Adeli, 2014a); the water drop algorithm (a PBMH inspired by river dynamics, where flowing water particles redistribute virtual "soil" on a graph representation of the optimization problem, e.g., Siddique & Adeli, 2014b); the harmony search algorithm (a PBMH based on mimicking the process of musical improvization to achieve an optimal harmony, e.g., Siddique & Adeli, 2015b); the gravitational search algorithm (a PBMH based on Newtonian gravity and the laws of motion and a population of "mass particles" interacting as they traverse the search space, e.g., Siddique & Adeli, 2016); particle swarm optimization (PSO, a PBMH based on the behavior of swarming insects and the gradual diffusion of information through local neighbourhoods); spider monkey optimization (a PBMH based on the sub-group foraging behavior of spider monkeys, e.g., Akhand et al., 2020).
The popularity of the application of nature-inspired algorithms to civil engineering optimization problems is due to the complexity of these problems and the relative ease with which PBMHs can be applied. Civil engineering optimization problems are commonly characterized by non-linearity, multi-modality, high variable interdependencies, search space discontinuities, high variable dimensionality and multiple objectives, which can pose problems for traditional optimization strategies (Maier et al., 2019;Szemis et al., 2013). PBMHs are able to overcome many of these problems as they are inherently global search algorithms owing to their population-based strategies, meaning that they are suited to multi-modal and multi-objective problems and can scale with increasing problem dimension. Additionally, being effectively sample-based strategies, they do not require gradient information, meaning they are suited to non-linear problems with high variable interdependencies and discontinuities in the search space.
In the last 30 years, the use of multi-objective evolutionary algorithms (note that MOEAs a subset of PBMHs that evolve populations of solutions [e.g., PSO] as opposed to evolving a problem representation [e.g., the water drop algorithm]) have arguably become the preferred optimization approach in civil engineering water resource problems (e.g., Maier et al., 2014;Mala-Jetmarova et al., 2017. The MOEA process can be viewed as the selection and perturbation of an existing population of solutions by sequential application of specific operators (e.g., crossover and mutation), with the aim of generating new and improved solutions to iteratively refine the quality of solutions in the population set (Asadzadeh et al., 2014).
A key advantage of MOEAs is the identification of a set of solutions that provide near optimal trade-offs between competing objectives in a single optimization run . These solutions are called the approximate Pareto optimal solutions, and the set of these solutions is called the approximate Pareto front (as they generally provide an approximation to the true Pareto front). Typically, the shape of the approximate front (e.g., whether it is convex or non-convex) is unknown a priori and depends on a problem's characteristics. However, the shape of the approximate front is found to be convex for many problems in practice. For example, Asadzadeh et al. (2014) stated that for multi-objective hydrologic model calibration, the approximate front shape is usually convex (i.e., a knee-bend in the middle and long tails), as the front can be reasonably well approximated by a convex sequence of piecewise linear curves or hyperplanes. This statement is supported by numerical experiments involving hydrological model calibration (Fenicia et al., 2007b;Kollat et al., 2012;Lee et al., 2011). For water distribution system (WDS) problems, Wang et al. (2015) evaluated 12 WDS case studies, and within this work, the shapes of the approximate front were all found to be convex. Moreover, the shape of the 24 bi-objective hydraulic model calibration problems' approximate Pareto fronts, reviewed by Asadzadeh et al. (2014), confirms this pattern also.
As the actual Pareto front is not known for most reallife problems, the ability of different MOEAs to generate a high-quality front is often compared using metrics that correspond to different attributes of the front. For example, the average distance of the approximate front to the "ideal point" (i.e., the point in the objective space that dominates all other points) is compared to assess an MOEA's convergence ability, where the front closest to this point is preferred. Alternatively, the extent of the front is measured by evaluating the diversity of the set of solutions (e.g., the degree of spreading).
The ability of an MOEA to generate a high-quality front is affected by its search strategy, including the following factors: operators, which govern how offspring solutions are produced; a hyperheuristic, which manages the utilzation of each operator throughout the search; and a selection strategy, which determines which solutions are selected to be used by the operators to produce the offspring solutions and which offspring solutions survive to join the population. Based on how the selection strategy is utilized, there are two types of MOEAs, namely, generational and steady-state. A generational MOEA generates multiple offspring solutions within each generation, where generally the number of offspring solutions is equal to the size of the population. As a consequence of this process, the offspring solutions compete with each other and the existing population, resulting in the survival of the best solutions (Zapotecas-Martinez & Menchaca-Mendez, 2020). A steady-state MOEA involves the selection of two individuals from the population to generate a singular offspring within each generation, which replaces the worst-performing solution from the population. Within both of these algorithms, the selection strategy is particularly important, as it needs to be designed to drive the population to converge to increasingly fit regions of the search space while avoiding premature convergence to sub-optimal regions (Back, 1996;Hanne, 1999).
In the past 20 years, many selection strategies have been proposed and compared in the literature. For example, Emmerich et al. (2005) applied a "hypervolume contribution" (HVC 1 ) selection strategy (Knowles et al., 2003) to a steady-state multi-objective selection based on the dominated hypervolume. The HVC 1 caters to both the convergence and diversity of a Pareto front by measuring the contribution of each individual solution to the hypervolume (Zitzler & Thiele, 1999). In comparative experiments, HVC 1 outperformed the "non-dominated sorting genetic algorithm II" (NSGA-II; i.e., a generational algorithm) with a "crowding distance" (CD) selection strategy, which focuses on maintaining Pareto front diversity . Recently, Laszczyk and Myszkowski (2019) proposed a clone prevention method as an alternative, which focuses on increasing diversity in the solution space, as opposed to other methods that focus on diversity in the objective space.
However, these studies did not isolate the impact of the selection strategy from that of the other algorithm processes, which poses a difficulty in terms of being able to attribute any performance differences to the different selection strategies used. This limitation was addressed in part by Igel et al. (2007), who compared two selection strategies (HVC 1 and CD) by applying them within the same generational algorithm, the multi-objective covariance matrix adaptation evolution strategy. The results indicated that the algorithm with HVC 1 outperformed that with a CD selection strategy on a range of test functions.
Building on this work, Asadzadeh and Tolson (2013) compared the influence of four selection strategies using a steady-state MOEA termed the Pareto archived dynamically dimensioned search (PA-DDS). The selection strategies considered in this work were HVC 1 , an alternative HVC termed HVC 2 based on the work of Bader and Zitzler (2011), CD, and a purely random strategy (RND), which did not take account of any of the solutions' attributes. The mechanism of HVC 2 is akin to that of HVC 1 but has an additional parameter that defines the maximum number of solutions that should be considered in the calculation of HVC 2 (the interested reader should refer to Bader & Zitzler 2011). In the comparative study, it was shown that the HVC 1 selection strategy achieved the best overall performance over all test functions and water resources problems considered.
The above selection strategies are designed to be applied to a wide range of problems, regardless of their char-acteristics. However, it has been shown that algorithm performance can be improved if selection criteria are tailored to certain problem characteristics. Feng et al. (1997) and Cococcioni et al. (2007) showed that giving greater selection priority to solutions that are closer to the convex approximation of the Pareto fronts can improve the performance of an MOEA to problems with convex-shaped Pareto fronts. Given this finding, to tailor an algorithm to such convex-shaped problems, Asadzadeh et al. (2014) proposed a novel convex-hull contribution (CHC) selection strategy for PA-DDS (a steady-state algorithm). In this comparative study, the CHC selection strategy outperformed the HVC 1 selection strategy for both test functions and water resources problems and is therefore considered the current best-performing selection strategy for such problems. The approach proposed by Asadzadeh et al. (2014) gives a high selection priority to the non-dominated solutions that have greater CHC values. Typically, a solution's CHC is the difference in size (e.g., area or hypervolume) of the convex hull set between the approximate front with and without that solution. This selection strategy not only accounts for convergence and diversity in generating the approximate front, but its search behavior is well-suited and has been shown to outperform other search strategies for problems with Pareto fronts that have a convex shape, such as the WDS problem (Jahanpour et al., 2018). The reason for this is that more solutions in the knee region are improved during the search, as the solutions in this region normally have a greater CHC value, thereby having a greater chance to be selected for reproduction and to be retained in the population (Asadzadeh & Tolson, 2013;Jahanpour et al., 2018). In addition, better performance in the knee region is critical when MOEAs are used to support real-world decisionmaking, as these regions are of most interest as they provide a locally distinct compromise of each objective (Hadka & Reed, 2012;Mala-Jetmarova et al., 2018). While the use of the CHC selection strategy is only applicable to a sub-class of problems, most real-world optimization problems are likely to fall within this sub-class as discussed previously.
Although the use of the CHC selection strategy has been shown to improve algorithm performance for convex problems, this assessment has only been undertaken for steady-state MOEAs, as the current formulation of the CHC selection strategy cannot be directly applied to generational MOEAs. However, many studies have shown that generational MOEAs perform favorably in comparison to steady-state MOEAs. For example, Wang et al. (2015) compared the performance of 5 MOEAs (three generational and two steady-state MOEAs) for 12 WDS problems, where the generational algorithms typically achieved a better performance than the steady-state ϵ-MOEA, which is a fast MOEA for finding well-spread Pareto-optimal solutions (Deb et al., 2003). Also, in this work, the steady-state algorithm Borg (an auto-adaptive many-objective evolutionary computing framework;  exhibited effective performance on the largest case study of the 12 WDS problems, but its overall performance was found to be inferior to that of NSGA-II and AMALGAM (the latter, a genetically adaptive multi-objective multi-algorithm; Vrugt & Robinson, 2007), which are generational algorithms. Thereafter, Wang et al. (2017) proposed a new generational MOEA, called GALAXY (i.e., genetically adaptive leaping algorithm for approximation and diversity), that outperformed two steady-state MOEAs (ϵ-MOEA and Borg) on 5 WDS problems. In addition, Wang et al. (2020) tested a general MOEA framework with 12 operators on six WDS problems. The results show the proposed MOEA outperformed Borg and GALAXY. Consequently, given the favorable performance of generational MOEAs and that of the CHC selection strategy, there is value in developing a generic formulation of the convex hull selection criterion that can be applied to all types of MOEAs so that the benefit of the CHC selection strategy can be implemented within, and improve, the performance of generational MOEAs for problems with convex Pareto fronts.
The objectives of this study are therefore to (1) develop a new convex hull selection strategy formulation for generational MOEAs (CHC Gen ), (2) compare the performance of the CHC Gen selection strategy with that of existing selection strategies, and (3) test the potential of the proposed CHC Gen selection strategy to improve the performance of existing MOEAs. The above objectives are achieved by conducting an extensive numerical experimental program involving numerous WDS case studies as in the work of Wang et al. (2015Wang et al. ( , 2017, Zheng et al. (2016Zheng et al. ( , 2017, and Wang et al. (2020).
The remainder of this paper is structured as follows. The Methodology section formulates the CHC Gen selection strategy, and a general MOEA framework is introduced to provide a generic generational algorithm structure to investigate the impact of different selection strategies. Additionally, two existing MOEAs are outlined with modifications proposed by including the CHC Gen as their selection strategy to evaluate CHC Gen 's potential benefits for incorporation within other generational MOEAs. The structure of the computational experiments is also outlined. In the Results and Discussion section, the results of the numerical experiments are reported, which highlight the effectiveness of the proposed CHC Gen selection strategy, followed by the Conclusions section. F I G U R E 1 Generic multi-objective evolutionary optimization process (after Wang et al., 2020) 2 METHODOLOGY

Background
The general steps within an MOEA are shown in Figure 1 ( Wang et al., 2020). As can be seen, at the commencement of the optimization process, an initial set of solutions is randomly generated and forms the population x t at generation t with size N. Then, subject to the parent selection process, some solutions are selected as parent solutions that will have the opportunity to reproduce offspring y t . The reproduction process is facilitated by one or more operators (typically a recombination operator and a mutation operator). In addition, the degree to which an operator contributes to the search at each generation can be controlled with the aid of hyperheuristics, which are high-level automated search strategies for selecting the most appropriate lower-level operators (or heuristics, such as mutation and cross-over; Burke et al., 2013;Drake et al., 2019). Thereafter, x t and y t are combined to form the combined set c t . Replacement is carried out to select the next generation x t+1 . The above process is repeated until certain termination criteria are met, such as the execution of a fixed number of generations or until no better solutions are identified.
The purpose of the process of parent selection and replacement is to identify better solutions from the current solution set. A general process of parent selection and replacement is summarized as follows. For parent selection, in order to select better solutions from population set x, the primary criterion value for solution x i is given by l(x i |x), where l indicates the relative rank of x i with respect to the population x (e.g., non-dominance rank). Given this criteria value assignment, the population is organized into the sets 1 , ⋯, , based on their values, where if x i and F I G U R E 2 Replacement step for a generational multi-objective evolutionary algorithm (MOEA; after Deb et al., 2002) x j are in , then l(x i |x) = l(x j |x) and l(x i |x) < l(x k |x) for all ∈ for which G > K (i.e., each contains solutions of the same criterion value and indicates a better rank for lower K). The secondary selection criterion applies to each solution in each and is denoted as v(x| ), ∈ , which enables the ordering of solutions in as The sorting of the population x is then based on the overall rank of a solution with respect to the entire population, which is first based on l (i.e., the set that a solution is a member of) and secondarily on v. As indicated, in general, the primary selection criterion is associated with solution quality (e.g., non-dominance rank), and the secondary selection criterion is associated with the diversity of solutions (e.g., CD).
Replacement follows a similar approach to that used for parent selection as outlined in Figure 2. This stage involves the sorting of c t to select N solutions to form x t+1 . First, c t undergoes a non-dominance sorting process, which involves the allocation of solutions in c t into the ordered subsets 1 , 2 ⋯, where indicates the solutions from c t that form the ith order Pareto front (e.g., 1 are the Pareto optimal solutions for the entire set c t , 2 are the Pareto optimal solutions for the reduced set c t \ 1 , etc.). As |µ 1 | is typically smaller than N, all of the solutions in µ 1 are used to form the first |µ 1 | solutions in x t+1 . The remaining available N-|µ 1 | positions in x t+1 are filled by subsequent solutions in µ 2 . This process is repeated until the point K, for which |µ 1 |+ ⋅⋅⋅ +|µ K+1 | > N. At this point, the set µ K+1 is sorted according to the second selection criterion, where the top-ranked solutions are used to fill the remaining positions in x t+1 . This step represents the final step for a given generation.

F I G U R E 3
Convex hull of non-dominated points in a bi-objective space, where the axes correspond to the two objective values (after Asadzadeh et al., 2014). CHC, convex hull contribution

CHC selection strategy for generational MOEA
To outline the new CHC Gen selection strategy, the concepts of convex hull and convex hull contribution (CHC) are introduced, and the approach of calculating the CHC value for a solution, as outlined in Asadzadeh et al. (2014), is reviewed. Thereafter, the detailed procedure of CHC Gen for generational MOEAs is introduced.

Convex hull background
For an m-dimensional space, the convex hull of a set of points S∈ℝ m is the intersection of all convex sets containing S (Barber et al., 1996). An example of the convex hull (gray-filled area) of a given set of points (empty circles) set in a two-dimensional space is given in Figure 3. It can be seen that the convex hull area contains all of the given points. From Figure 3, a convex hull is bound by the segment line of two vertices (dot points) called facets (solid and dash lines). The CHC is described in the following. Each facet divides the space into two sides, one side contains the convex hull, and the other does not. The area of the region bounded by the facets (or area in a twodimensional space) is called the convex hull size. If any of the vertices are removed from a convex hull set, the size of the convex hull will change (Asadzadeh et al., 2014). For instance, if the vertex p is removed from the convex hull in Figure 3, the new facet (the dot-dash line) is the new boundary of the new convex hull, thereby reducing the size of the convex hull by an amount given by the sparse lined area in Figure 3. In m-dimensional space, the reduced volume (or area) that results from removing any vertices (e.g., that given by the lined area) is called the CHC and is given by: where S is a set of vertices, and V(S) is the convex hull size of the set S.
The key steps for CHC calculation for a non-dominated set are as follows (interested readers should refer to Asadzadeh et al., 2014, for more details). For the first step, the axes of a non-dominated set are normalized in the objective space to eliminate any potential bias resulting from varying scales of the objective functions. For the following explanation, consider that the two-dimensional space in Figure 3 is representative of an objective space, where the criterion for both objectives is a minimization (i.e., the lower-left side of the polygon is associated with the Pareto front). The solutions are divided into four groups as shown in Figure  To calculate the CHC for each solution, zero CHC values are assigned to the solutions in groups (i) and (ii) mentioned above. The calculation of CHC is only conducted for solutions in group (iii), where, for example, the CHC for point p is given by the sparse lined area in Figure 3. The CHC for the solutions in group (iv) is assigned as the CHC value of the closest solutions' CHC values from group (iii) (Asadzadeh et al., 2014). The original CHC selection strategy is designed for a steady-state MOEA (the aforementioned PA-DDS), for which the population set only contains non-dominated solutions. However, for a generational MOEA, the CHC selection strategy is not applicable for two reasons. First, a generational MOEA's population set may contain dominated solutions during the search , which are beneficial for solution diversity. In this case, dominated solutions may have greater CHC values than non-dominated solutions. For example, as outlined in Davoodi et al. (2011), the solution at the top facet is dominated but has a positive CHC value. Applying the CHC selection strategy directly within a generational MOEA would result in a higher selection probability to such dominated solutions, thereby potentially resulting in poor performance. Second, the portion of solutions with a non-zero CHC from a given set is usually small (Jahanpour et al., 2018). In other words, a great portion of solutions' CHC values are zero, and as a result, they become noncompetitive in a tournament selection process. This will result in a highly biased selection process and may not benefit an algorithm's performance. In summary, there is a need to propose a new CHC selection strategy for generational MOEAs.

2.2.2
Algorithm of the proposed CHC for a population-based MOEA The CHC Gen selection strategy for population-based MOEAs at each generation t is outlined as follows: 1: INPUT the combined set c t 2: CONSTRUCT the non-dominated sets 1 , 2 , … of c t 3: SET = 1, = 1, = ∅ 4: WHILE i ≤ r 5: SET z = 6: WHILE ≠ ∅ 7: CONSTRUCT the set The input to the algorithm is the combined set c t that comprised the parent and offspring populations (Line 1). This set is then organized into a sequence of nondominated sets 1 , 2 , ⋯, (Line 2). The loop variables are then initialized (Line 3) and the algorithm then loops through all non-dominated sets (Lines 4 to 13). The working variable z for the inner loop is defined as in Line 5, and then the inner loop builds the modified non-dominated sets from z (the u l variables, which for CHC Gen , replace the in Figure 2) until z is empty (Lines 6 to 11). The working set w is constructed as all elements of z for which the CHC value is positive (Line 7). This is to ensure that no solutions with zero or negative CHC values enter the modified non-dominated sets. The modified non-dominated set u l is then constructed from w as the ordered set, where the ordering of the set is based on the highest CHC values (i.e., the secondary criterion, ν; Line 8). The new u l is then concatenated to the end of the working ordered set u (Line 9). The set w is removed from y, the primary selection criterion l is updated (Line 10), and the inner loop is iterated (Line 11). When z is empty, the counter i is incremented (Line 12), the process is looped (Line 13), and the next non-dominance set is selected (Line 5). The outer loop is iterated for all non-dominated sets. At the completion of the outer loop, the working ordered set u comprised all elements of c t ordered first by the primary criterion and then by the secondary criterion. The new population +1 is then constructed as the first N elements of (Line 14), and then the new population is returned (Line 15).
In summary, the operation of CHC Gen is, starting from the combined set c t , to incrementally construct the F I G U R E 4 Overview of computational experiments. CD, crowding distance; GALAXY, genetically adaptive leaping algorithm for approximation and diversity; HR, hybrid replacement; HVC, hypervolume contribution; NSGA-II, non-dominated sorting genetic algorithm II; RND, random strategy; WDS, water distribution system. modified non-dominated sets u 1 , u 2 . . . , whose elements are ordered according to their respective CHC values. Of these modified non-dominated sets, the top N solutions will be retained and will comprise the next generation x t+1 .
The reason why it is necessary to use the modified non-dominated sets, as opposed to the standard nondominated sets, is because of the presence of solutions in the non-dominated set for which the CHC value is zero. That is, it is possible to have a non-dominated solution that possesses a zero CHC value (i.e., see the purely white point in Figure 3). As these points have a zero CHC value, the secondary selection criterion is unable to discriminate between the value of these points, and any ability to rank these solutions by the secondary criterion is lost. As such, it is necessary to incrementally remove the positive CHC values from the non-dominated sets (and store them in the modified non-dominated sets) in order to be able to re-valuate the CHC values and rank the solutions in the new modified non-dominated set by the secondary selection criterion. This approach enables all solutions within the non-dominated sets to be effectively utilized in the ranking process.

Numerical experiments
An overview of the computational experiments is shown in Figure 4. These were designed to compare the CHC Gen selection strategy with other existing selection strategies (Objective 2) and to investigate the effectiveness of applying the new CHC Gen selection strategy to existing MOEAs (Objective 3). In addition, it is noted that Objective 1 has been addressed in Section 2.2. The comparative study was undertaken through the application of a number of algorithms to six WDS case studies. A range of metrics was used to assess the algorithms' performance and the behavior of different selection strategies. Details of the numerical experiments are outlined below.
To adopt a systematic approach to compare the performance of different selection strategies (Objective 2, Figure 4), it is required to isolate the impact of the selection strategy on algorithm performance. As explained below, this is achieved through the adoption of a common algorithm framework to incorporate a given selection strategy while holding all other algorithm operations constant. The performance comparison is based on the selection strategy that results in better metric values, thereby indicating better performance. Moreover, to understand the behavior of each selection strategy, a selection metric is proposed in Section 2.6.2 to characterize the rank of solutions along the front. This is then applied in the analysis of the strategies' approximate fronts to characterize the distribution of solutions along the front.
The general MOEA framework (Wang et al., 2020) is a generational MOEA that has interchangeable components (e.g., operators, hyperheuristic, and selection strategy), thereby enabling the impact of different selection strategies to be assessed explicitly. As part of this framework, the operator set and hyperheuristic remain unchanged, and the strategies of parent selection and replacement are varied. Each constructed MOEA is named by the selection strategy utilized. For example, the general MEOA framework utilizing the CD selection strategy is denoted as Algorithm-CD. In addition, to assess the absolute performance of the general MOEA framework with the CHC Gen selection strategy, two popular generational MOEAs applied to water resources problems, NSGA-II and GALAXY, were modified by embedding the new CHC Gen selection strategy and comparing their performance with that of the original version of the algorithms.
In order to assess the influence of different selection strategies on MOEA performance, the CHC Gen and four existing selection strategies (which are hyper volume contribution [HVC], crowding distance [CD], hybrid replacement [HR], and random [RND]) were adopted within the general MOEA framework and compared with each other (i.e., in the largest solid block in Figure 4). These selection strategies were considered as they have been implemented in popular generational MOEAs (e.g., NSGA-II, GALAXY) and have shown effective performance for WDS problems (Wang et al., 2015(Wang et al., , 2017.
As introduced above, in order to test the utility of the new CHC Gen selection strategy to improve existing MOEAs, the CHC Gen selection strategy was adopted within the two existing MOEAs, which are denoted as NSGA-II-CHC Gen and GALAXY-CHC Gen . The two modified algorithms were compared with the original versions of these algorithms, NSGA-II and GALAXY. The general MOEA framework with the CHC Gen strategy was also included in the comparison to investigate if the relative influence of the selection strategy outweighs that of the hyperheuristic.
The above numerical experiments were conducted by using a bi-objective optimization problem (i.e., minimizing network cost and maximizing network resilience) for the six WDS problems used by Wang et al. (2020). These case studies include: the New York tunnel network (NYT), the Hanoi network (HAN), the Fossolo network (FOS), the Pescara network (PES), the Modena network (MOD), and the Balerma irrigation network (BIN). These case studies have been widely used to assess MOEA performance (Wang et al., 2015(Wang et al., , 2017Zheng et al., 2016). As shown in the central block of Figure 4, all optimization runs were duplicated 10 times with different starting positions in the decision domain.
For the "Result Assessment" block at the bottom of Figure 4, the results for the different objectives were compared by applying the one-way Kruskal-Wallis test (Kruskal & Wallis, 1952) with Dunn's D post-test (Dunn, 1964) to three end-of-run performance metrics. Moreover, a novel visualization metric, called the selection metric, is introduced and implemented to understand how different search strategies affect an algorithm's search.
The code of the general MOEA framework was adopted from Wang et al. (2020); the code of the selection strategies considered in this study was written as a MATLAB m-script; and the CHC Gen selection strategy, the convex hull's size (i.e., Lebesgue measure), vertices, facets, and so forth, can be acquired by implementing the "qhull" code, which is available at: http://www.qhull.org/ ( Barber et al., 1996). The best-known solutions were compared with reference Pareto fronts (the best-known Pareto fronts) for the six WDS case studies. Each best-known Pareto front is a result of working with merged solutions from several MOEAs with high computational budgets. The details of the fronts can be found in Wang et al. (2015) and Jahanpour et al. (2018). All simulations were run on the Phoenix high-performance computer (HPC) at the University of Adelaide, Australia. The Phoenix HPC is a heterogeneous hardware system that includes a mix of CPU-only and CPU/GPU-accelerated nodes. It has 260 nodes in total, which are equipped with 2x Intel Gold 6148, 40 cores @ 2.4 GHz, and 384GB memory for CPU nodes, with a maximum of 125 GB RAM per node.

Selection strategies considered
Details of the four selection strategies considered (see Section 2.3) are given in the following sub-sections.

Crowding Distance
The CD selection strategy is a popular selection strategy that has been widely applied in many MOEAs since it was proposed by Deb et al. (2002). It uses non-dominace ranks as (primary selection criterion) and CD as (secondary selection criterion). CD measures the diversity of solutions within the same non-dominated front in the objective space. For the non-dominated front, the CD values are set as infinity for extremal solutions and for the sum of side lengths of the segment lines that touch the neighboring points for the non-extremal solutions. Generally, a solution with a greater CD, which indicates a solution of greater uniqueness, is far from its neighbors in the objective space .

Hypervolume Contribution
The HVC selection strategy implements the hypervolume measure proposed by Zitzler and Thiele (1998). Solutions with a greater value of HVC are those that are normally located within the knee region of a Pareto front (Emmerick et al., 2005). The HVC selection strategy uses non-dominated rank as and HVC as . The HVC for a target solution within a non-dominated front ( ) is the difference between the hypervolume of all solutions in the front and the hypervolume of the all solutions in the front excluding the target solution. In addition, the HVC of the extremal solutions is usually dependent on reference point selection (Igel & Hansen, 2007). However, in this study, the HVC for the extremal solutions are equal to the HVC of the solutions next to each extremal solution, which is consistent with Asadzadeh et al. (2014).

Hybrid Replacement
The HR selection strategy was first introduced in GALAXY, which is a new hybrid MOEA that is tailored to WDS problems (Wang et al., 2017). The HR selection strategy contains two secondary selection criteria of the CD and the ϵ-dominance criterion (Deb et al., 2005) as . In the parent selection, the CD is activated for selecting parent solutions for reproducing offspring solutions. However, in the replacement step, the secondary selection criterion being used depends on the number of non-dominated solutions in the first front ( 1 ). If the number of elements of 1 is smaller than N, CD will be activated; otherwise, the ϵ-dominance selection is activated, forcing the ϵ-nondominated solutions to be included in x t+1 . For the later criterion, the objective domain that is bounded by the two solutions on the tail of 1 is divided into (N-2) by (N-2) grids. Each grid is denoted as an ϵ-box, where the dominant region within the ϵ-box is at the right bottom corner in this study (i.e., maximize resilience and minimize cost objective values). Solutions with the smallest Euclidean distance to the dominant corner are recognized as ϵ-non-dominated solutions and are included in x t+1 . The slots remaining in x t+1 are filled by the remaining ϵ-dominated solutions that are ordered by the distance to the corresponding dominant corner.

Random
The RND selection strategy was used as a control selection strategy in this study. It uses the non-dominated fronts as . The for each solution is randomly sampled in the range [0, 1] from a uniform distribution, removing any bias in the solution ordering.

Two-objective WDS optimization problem
In order to enable the performance of different optimization algorithms to be assessed in an objective, transparent, consistent and repeatable manner, Wang et al. (2015) collated a set of case studies for WDS optimization problems that could be used as performance assessment benchmarks. These are classified into different levels of complexity according to the size of the search space. In this paper, six WDS case studies were selected from that set.
The WDS optimization problem involves the selection of pipe diameters for a WDS network to optimize set criteria. This selection aims to achieve the design of the lowestcost network that satisfies performance constraints, such as the minimum pressure for each node in the network and constraints on the fluid velocity in each pipe. This type of problem is difficult to solve as it is NP-hard, nonconvex, high dimensional, multi-modal, and non-linearly constrained (Zecchin et al., 2012).
The objective functions used in this study are consistent with those used commonly for this problem type (Wang et al., 2015, Jahanpour et al., 2018, namely, maximiz-ing network resilience and minimizing network cost. The cost objective is given by where F c is the total network cost, which is determined by pipe diameter and pipe length ; a and b = specified cost function coefficient and exponent (i.e., the values refer to the settings in Wang et al., 2015); = the total number of pipes in the network. The network resilience objective is given by: where I n is the network resilience; m is the total number of demand nodes; , , and * are the demand, actual head, and minimum head required at each node j, respectively; N R is the total number of reservoirs; , are the actual discharge and actual head at the reservoir r; and is an indicator of diameter uniformity for pipes that are connected to node j and is defined by: where is the diameter of the pipe i connected to node j; and N pj is the total number of pipes that are connected to node j. Note that a larger value of represents a higher reliability of the network node since the diameter variations between these pipes are lower overall ( = 1 when all pipe diameters are identical; Prasad & Park, 2004).
In this study, a set of integer options, ranging from 1 to the number of commercially available diameter sizes, is used as the decision variable. The constraints of the WDS optimization problem in this study are the flow velocity in each pipe and pressure head at each node as specified by each case study. The EPANET 2.0 (Rossman, 2000) hydraulic simulation software is used to evaluate the flow rates and pressure heads for each pipe and node, respectively, to compute the constraints and resilience of the network.

2.5.2
WDS case studies and parameter setup As mentioned above, six WDS design problems commonly adopted in a wide range of MOEA studies (Jahanpour et al., 2018;Wang et al., 2015Wang et al., , 2017Zheng et al., 2016) were considered. These were chosen to provide a range of problem Abbreviations: MOEAs, multi-objective evolutionary algorithms; WDS, water distribution system. characteristics and sizes (see Wang et al., 2015). Table 1 provides details of the six WDS case studies, where the number of pipes varies from 21 to 454. The configurations of population size and computational budget (i.e., the number of function evaluations, N NFE ) were consistent with the configurations in Wang et al. (2017), as these were found to provide satisfactory outcomes for the case studies. The population size was set as 100 for the NYT, HAN, FOS, and PES problems, and for the large-scale WDS case studies, MOD and BIN, a population size of 200 was used.

2.6
Result assessment

End-of-run performance metrics
Three end-of-run performance metrics, "hypervolume" (I HV ; Zitzler & Thiele, 1999), "generational distance" (I GD ; Veldhuizen, 1999), and the "ϵ-indicator" (I ɛ+ ; Zitzler et al., 2003), were used to assess the relative performance of the algorithms. These metrics effectively capture both the convergence and diversity of an algorithm's approximate front. I HV is the ratio of the dominated volume of an approximate front, compared with a reference Pareto front representing both the convergence and diversity of solutions. I GD is the average distance between an approximate front and the reference Pareto front, in terms of evaluating convergence. I ϵ+ evaluates the minimum distance required to shift the approximate front to dominate the reference Pareto front, which measures the convergence and consistency of a solution set. As mentioned in Section 2.3, to yield a robust comparison among the algorithms' performance metrics, the one-way Kruskal-Wallis test (Kruskal & Wallis, 1952), with Dunn's D post-test (Dunn, 1964), were implemented to evaluate if a pair of algorithms' end of run data significantly differ from each other. This non-parametric analysis provides a statistical test of whether the means of two or more groups of data are equal (Ameca-Alducin et al., 2018;Hadka & Reed, 2012). If the difference is not statistically significant, the pairs' data are assigned as being equivalent.
Otherwise, the algorithm with the better metric median value is assigned as the better-performing algorithm. The pairwise statistical tests were conducted for all algorithm group pairs, where the number of times an algorithm was recorded as being the better performer was recorded and aggregated to yield the overall rank.

Selection metric
As mentioned in Section 2.3, to provide insight into how a selection strategy affects MOEA search, a novel visual metric is proposed, termed the selection metric. The metric is able to characterize the high-quality solutions along the front by combining the primary and secondary selection criteria values, as mentioned in Section 2.1, to give the solutions an overall rank. Solutions with a higher overall rank indicate a higher likelihood to be successful in the selection step. For example, for the CHC Gen selection strategy, a solution with a higher rank of the modified nondominated front (primary selection criterion) and a lower value of the CHC value (secondary selection criterion) has a higher overall rank than a solution with a lower rank of the non-dominated front with a higher value of the CHC value. The selection metric (Sel) is defined for solution x from population x as where norm ( | ) is the normalized primary selection criterion value for solution x from population x, given by: where ( | ) is the primary selection criterion value (i.e., the ranking of x in x); l max is the greatest primary selection criterion value; ( | ) is the secondary selection criterion value for solution x in sub-population y; and is the subpopulation of solutions of primary selection criterion rank i. The purpose of this metric is to rank solutions by taking F I G U R E 5 Ranks of the general MOEA frameworks with the different selection strategies for the six case studies considered. Algorithm-CHC Gen is the general MOEA using the generational convex hull contribution (CHC) selection strategy; Algorithm-HR is the general MOEA using the hybrid replacement (HR) selection strategy; Algorithm-HVC is the general MOEA using the hypervolume contribution (HVC) selection strategy; Algorithm-CD is the general MOEA using the crowding distance (CD) selection strategy; and Algorithm-RND is the general MOEA using the random selection strategy. BIN, Balerma irrigation network; FOS, Fossolo network; HAN, Hanoi network; MOD, Modena network; NYT, New York tunnel network; PES, Pescara network account of the influence of both primary and secondary selection criteria values. In this study, the primary selection criterion values are greater than one and have to be normalized to the range [0,1], which is the same as the secondary selection criterion value range.

Performance comparison
The rank of the three end-of-run metrics and the average ranks of the general MOEA framework algorithms with the different selection strategies for the six WDS problems are shown in Figure 5 as categorical surface plots, where shades from white to black indicate the best to the worst ranks, respectively. In Figure 5a,c, the case studies are presented in ascending order of complexity on the horizontal axis, and the general MOEA framework algorithms that are embedded with the different selection strategies are listed on the vertical axis. In Figure 5d, the first three columns are the average ranks across the six case studies for each algorithm for the associated metrics. In addi-tion, the fourth column is the overall average rank, which assesses the performance of each algorithm, across the three metrics' average ranks. As shown in Figure 5, it can be observed that Algorithm-CHC Gen generally performs better than the other algorithms using existing selection strategies with respect to the three end-of-runs metrics, as well as average ranks. In contrast, Algorithm-RND is the worst-ranked algorithm for all metrics and case studies (Figure 5a-c). Moreover, Algorithm-HR achieves the second-best average rank (Figure 5d). This implies the effectiveness of the recently developed HR selection strategy may have been an important contributing factor to the performance of GALAXY (Wang et al., 2017). In addition, Algorithm-HVC outperforms Algorithm-CD, which is consistent with the findings by Igel and Hansen (2007;Figure 5d). Given that the only component that is varied in the general MOEA framework algorithms is the selection strategy, these results show that applying different selection strategies affect MOEA performance and, furthermore, that the proposed CHC Gen selection strategy outperforms the existing selection strategies investigated for the WDS problems considered.
In addition, the benefit of adopting CHC Gen increases with the problem scale. For example, From Figure 5b,c, apart from Algorithm-RND, the algorithms using the CHC Gen , HVC, CD, and HR selection strategies typically achieve similar rank values in the small-scale problems (NYT and HAN). However, for large-scale problems, such as MOD and BIN, Algorithm-CHC Gen outperforms other algorithms most of the time. This finding suggests that adopting the CHC Gen selection strategy can benefit MOEA performance on problems with a high degree of complexity for this type of problem.

Discussion
To understand how different selection strategies' features affect an algorithm's search performance, the behavior of the algorithms with the four selection strategies (i.e., CHC Gen , HVC, HR, and CD) was assessed with the aid of the selection metric. Figure 6 shows the population of solutions after 20,000 N FE corresponding to a run applying Algorithm-RND to the BIN problem. The color spectrum indicates the assigned selection metric of the solutions for assessment based on the different selection strategies (the shading from blue to orange indicates the best to worst selection metric values, respectively). To aid in identifying the solutions with higher values, the top 20 ranked solutions, which have a higher probability to be selected as parent solutions for reproducing offspring solutions, are highlighted by red circles and denoted as preferred solutions. In

F I G U R E 7
Approximate front solutions of the BIN problem for the general MOEA frameworks with different selection strategies. The insert is a zoomed in view of the "knee region" of the Pareto fronts. this example, overall, the four population sets produced by conducting replacement as part of the four selection strategies are similar; however, the distribution of the selection metric for the example solution set is different, resulting in different Pareto fronts (Figure 7), as discussed below.
From Figure 6a, it can be seen that the CHC Gen selection strategy typically prefers solutions (blue color) on the convex hull, which are located in the following: the low cost and knee region (M€5 to M€8); the medium cost range (about M€10.5; knee region of the front); and the highcost region (about M€20). This results in a front (blue front in Figure 7) that has the greatest degree of diversity in resilience values (from around 0.45 to over 0.95) and is closer to the reference Pareto front (the gray front in Figure 7) within the low and mid-cost region in comparison to the other selection strategies considered. However, the preferred solutions are not uniformly distributed across the front but clustered into smaller distinct regions along the front. As shown in Figure 6a, the least preferred solutions (yellow and orange colors) of Algorithm-CHC Gen extend from the cost range of M€14 to M€18. The reason is that the convex solutions are clustered at both ends and the knee region of the approximate front. This behavior results in the gap of the final approximate Pareto front in the cost range between M€14 and M€18 (Figure 7). Overall, based on the success of Algorithm-CHC Gen , it is implied that the portion of the solutions within the convex hull in the population set is sufficient to direct the solutions toward the best-known Pareto front, which is consistent with the findings of Jahanpour et al. (2018).
For the HVC selection strategy, the distribution of the preferred solutions (Figure 6b) is much more even in the objective space than that obtained using CHC Gen . For example, the preferred solutions are observed in the cost range between M€8 and M€10; M€14 and around M€18 (Figure 6b). The above-preferred solutions are on the nonconvex portion of the front and are not preferred by the CHC Gen selection strategy. This feature resulted in a Pareto front generated by the HVC selection strategy that is more evenly distributed than that obtained using the CHC Gen selection strategy as shown in Figure 7. An interesting and unexpected result is that the more uniformly distributed selection preferences from HVC are not as effective as the more clustered selection preferences by CHC Gen .
The differences in the ranking of the solutions in accordance with the selection metric values produced by the CD and HR selection strategies are subtle as shown in Figure 6c,d. However, the difference is made clearer by considering more closely the distribution of solutions within a selected region of the front (see the zoom-in of the ɛ-boxes inserted in Figure 6c,d). It is observed that, for the non-dominated solution set, the HR selection strategy preserves more ɛ-non-dominated solutions than the CD selection strategy (i.e., the three orange color solutions with the lowest CD values are excluded by the CD selection strategy). Due to the inclusion of ɛ-non-dominance (Laumanns et al., 2002), the HR selection strategy preserves a good balance between convergence and diversity (Wang et al., 2017). As shown in Figure 7, Algorithm-HR outperforms Algorithm-CD in terms of the convergence and diversity of its approximate front in comparison to the F I G U R E 8 Ranks of non-dominated sorting genetic algorithm II (NSGA-II), NSGA-II-CHC, genetically adaptive leaping algorithm for approximation and diversity (GALAXY), GALAXY-CHC Gen and the general MOEA framework with CHC Gen selection strategy best-known front. However, despite the fact that the CD and HR selection strategies are able to generate fronts with a relatively even distribution of solutions along the entire front (i.e., no significant gaps in the front as with CHC Gen ), the two fronts resulting from the application of the CD and HR selection strategies did not outperform the one resulting from the application of the CHC Gen selection strategy. For completeness, the approximate Pareto front plots for all other problems (in addition to BIN as in Figure 7) are given in the Appendix Figures A1-A5, where similar trends are observed.
In summary, the selection strategy has a clear effect on MOEA performance. To achieve a better performance of an MOEA, a selection strategy requires not only good diversity along the front and convergence to the best-known front but also the ability to shape the solution set to approximate the shape of the Pareto front. Given the convexity of the best-known Pareto front, a key advantage of the CHC Gen selection strategy is that it is not only able to maintain a good diversity and convergence along the front, but it also focuses the search around the convex hull solutions, which enables the front to maintain its convexity as the search evolves.

Application of CHC Gen selection strategy to existing MOEAs
The algorithm ranks for NSGA-II and GALAXY, as well as their variants using the CHC Gen selection strategy instead of their native selection strategy, are shown in Figure 8. As can be seen, the use of the proposed CHC Gen selection strategy improves the performance of NSGA-II and GALAXY. This is indicated by the fact that NSGA-II-CHC Gen and GALAXY-CHC Gen outperform NSGA-II and GALAXY (Figure 8d), respectively. Moreover, the influence of the selection strategy does not outweigh the hyperheuristic search strategy. For example, it can be seen that GALAXY outperforms Algorithm-CHC Gen but not GALAXY-CHC Gen (Figure 8d). Thus, the influence of the hyperheuristic outweighs that of the selection strategy. In addition, a large size of the operator set is also a benefit to algorithm performance. This is evidenced by the higher average rank of Algorithm-CHC Gen in comparison to that of NSGA-II-CHC Gen (Figure 8d) and agrees with the findings of Wang et al. (2020). This is demonstrated by GALAXY-CHC Gen outperforming Algorithm-CHC Gen ; and Algorithm-CHC Gen outperforming NSGA-II-CHC Gen , respectively (Figure 8d).

CONCLUSION
A novel selection strategy called the generational CHC (CHC Gen ) is developed for generational MOEAs in this study. The CHC Gen selection strategy not only accounts for convergence and diversity in generating the approximate front, but its search behavior is well-suited to problems for which the Pareto front is convex. Moreover, the CHC Gen selection strategy improves the performance of existing MOEAs.
The general MOEA framework with the CHC Gen selection strategy was found to outperform the framework with other selection strategies in a numerical study involving six WDS problems. The CHC Gen selection strategy showed the overall best convergence, diversity, and consistency of the approximate fronts that were generated. From a detailed comparative analysis of a range of selection strategies applied to the same starting population, it was observed that a key advantage the CHC Gen selection strategy has is that it is not only able to maintain good diversity and convergence along the front, but focussing the search around the convex hull solutions enables the evolution of the front to maintain its convexity as the search evolves. Given the nature of the convex hull solutions within a non-dominated set that are closer to the "ideal point" and in distinct regions along the approximate front, this type of selection preference leads to an improved convergence and diversity of the search. Therefore, the CHC Gen selection strategy allows the algorithm to effectively explore the search space, thereby outperforming the existing selection strategies.
To further investigate the potential of the CHC Gen selection strategy to improve the performance of existing generational MOEAs, the current best generational MOEA for solving WDS problems (GALAXY) and the industry standard generational MOEA (NSGA-II) were modified to incorporate the CHC Gen selection strategy. The CHC Gen selection strategy was found to be able to boost the performance of these two algorithms, suggesting that the proposed selection strategy could benefit other existing MOEAs.
Future work should involve the applications of MOEAs with CHC Gen and other selection strategies (such as that proposed by Laszczyk & Myszkowski (2019)), to a wider range of civil engineering problems. This will enable the study of the impact of different selection strategies on an MOEA's performance in these other domains. In addition, added insight could be achieved by studying the impact of different selection strategies on the run-time evolution of an algorithm's population.

A C K N O W L E D G M E N T S
This work was supported by the PhD research program at the University of Adelaide. We would like to thank the anonymous reviewers, whose comments have helped to improve the quality of this paper.
Open access publishing facilitated by The University of Adelaide, as part of the Wiley -The University of Adelaide agreement via the Council of Australian University Librarians.

F I G U R E A 2
Approximate front solutions of the HAN problem for the general MOEA frameworks with different selection strategies. The insert is zoomed in view of the "knee region" of the Pareto fronts.

F I G U R E A 3
Approximate front solutions of the FOS problem for the general MOEA frameworks with different selection strategies. The insert is zoomed in view of the "knee region" of the Pareto fronts.

F I G U R E A 4
Approximate front solutions of the PES problem for the general MOEA frameworks with different selection strategies. The insert is zoomed in view of the :knee region" of the Pareto fronts.

F I G U R E A 5
Approximate front solutions of the MOD problem for the general MOEA frameworks with different selection strategies. The insert is zoomed in view of the "knee region" of the Pareto fronts.