Development and Sensitivity Analysis of an Improved Harmony Search Algorithm with a Multiple Memory Structure for Large-Scale Optimization Problems in Water Distribution Networks

: The continuous supply of drinking water for human life is essential to ensure the sustainability of cities, society, and the environment. At a time when water scarcity is worsening due to climate change, the construction of an optimized water supply infrastructure is necessary. In this study, an improved version of the Harmony Search Algorithm (HSA), named the Maisonette-type Harmony Search Algorithm (MTHSA), was developed. Unlike the HSA, the MTHSA has a two-floor structure, which increases the optimizing efficiency by employing multiple explorations on the first floor and additional exploitations of excellent solutions. Parallel explorations enhance the ability in terms of exploration (global search), which is the tendency to uniformly explore the entire search space. Additional exploitations among excellent solutions also enhance the ability of local searches (effective exploitation), which is the intensive exploration of solutions that seem to have high possibilities. Following the development of the improved algorithm, it was applied to water distribution networks in order to verify its efficiency, and the numerical results were analyzed. Through the considered applications, the improved algorithm is shown to be highly efficient when applied to large-scale optimization problems with large numbers of decision variables, as shown in comparison with the considered optimizers.


Introduction
Following the introduction of the Genetic Algorithm (GA) [1], one of the most widely known meta-heuristic algorithms, Tabu Search (TS) [2], Simulated Annealing (SA) [3], Ant Colony Optimization (ACO) [4], Particle Swarm Optimization (PSO) [5], and Differential Evolution (DE) [6] were developed, mimicking natural phenomena, human behaviors, and physical phenomena.Subsequently, the Harmony Search Algorithm (HSA) was developed in 2001 [7,8], and many meta-heuristic algorithms followed, such as the Honeybee Algorithm (HBA) [9], the Firefly Algorithm (FA) [10], Cuckoo Search (CS) [11], the Bat Algorithm (BA) [12], the Water Cycle Algorithm (WCA) [13], and the Mine Blast Algorithm (MBA) [14]).Individual algorithms have been continuously enhanced through improvements, and recently they have actively been improved through combinations of different meta-heuristic algorithms or combinations with algorithms with different mathematical bases.Although many meta-heuristic algorithms have been developed over the last proposed a least-cost design method for WDNs under multiple loading conditions, and El-Ghandour and Elbeltagi [46] applied meta-heuristic algorithms, including the GA, PSO, and ACO, to three WDNs (Two Loop, New York, and El-Mostakbal networks) and compared the performance of these algorithms.Lee et al. [47] proposed a meta-heuristic algorithm motivated by a vision correction procedure, and they applied this to the Balerma network design problem.
However, among the WDNs that were applied in previous studies, those other than the Balerma network had small numbers of pipes (numbers of decision variables).In addition, the Balerma network took the form of a typical agricultural irrigation network.Going forward, diverse types of water distribution networks should additionally be considered in order to compare the efficiency levels of algorithms.Yoo et al. [48] implemented pipe diameter optimization of Saemangeum networks, considering diverse supply types (branched, looped, and pumped type), and compared and analyzed the optimization results by type with existing designs.
In the present study, a Maisonette-type Harmony Search Algorithm (MTHSA) was developed, based on the Harmony Search Algorithm (HSA).Using this MTHSA, optimal pipe diameter design problems, that is, finding low-cost pipe diameter combinations within limited hydraulic conditions, were solved.In particular, the MTHSA was applied to large-scale problems, and the effects of this application were analyzed, with the results presented.The first benchmark problem selected in this study was the Balerma network minimum-cost design, which is a representative example of an existing engineering optimization problem and is continuously being studied by various researchers.In this study, a sensitivity analysis was performed to evaluate the performance of the parameter combination of the newly developed optimization algorithm for the optimal design of the Balerma network, a representative optimization problem.In addition, in order to verify whether the developed optimization algorithm can be universally applied to pipe network optimization problems, the Saemangeum networks were additionally selected from among the existing optimization benchmark problems that have different characteristics from the Balerma network.

Problem Descriptions 2.1. Formulation of Optimal Pipe Design
In general, the minimization of pipe construction costs is considered an objective function, and the pressure and velocity requirements are used as a constraint.Equations ( 1)-(3) show the objective function with penalty functions.
In the above, C(D i ) denotes costs according to pipe diameters per unit length.Li, Di, hj, hmin/max, Pj, vj, vmin/max, Vj, np, and nn denote the pipe length, pipe diameter, nodal pressure at node j, minimum/maximum nodal pressure requirements, penalty function for checking the pressure constraint, link velocity at pipe j, minimum/maximum pipe velocity requirements, penalty function for checking the velocity constraint, number of pipes, and number of nodes, respectively.α, β, γ, and δ are constants in the penalty functions.
In addition, in order to obtain the nodal pressure, hydraulic simulations of the water distribution network should be performed.In this study, the proposed optimization method was coupled with a hydraulic simulator EPANET [49] for optimization purposes.This is a free-access and widely used water distribution network solver, which employs gradients to determine the mass and energy conservation equations [42].A pipe network design problem was solved by choosing the least-cost pipe diameter among commercial candidate pipes satisfying hydraulic flow demands, energy conservation, and head restriction.The mass conservation of node, energy conservation of loop, minimum/maximum pressure at node, minimum/maximum velocity in pipe, and pipe diameter size availability constraints can be defined by the following Equations ( 4)- (11): Here, Qin, Qout, and Qe denote flow into the node, flow out of the node, and the external outflow rate (demand) at the node, respectively.
Here, ∆Hk is the head loss in pipe k and NL is the total number of closed loops in the system.Equation ( 5) means that the sum of head losses considering the flow direction of individual pipes within the loop is "0".Head loss always occurs depending on the flow of water in individual pipes within the loop, but the loss depending on the flow direction appears as a positive or negative value.The head loss for each pipe is the head difference between inter-connected nodes, which can be calculated using head loss equations such as the Hazen-Williams, Darcy-Weisbach, and Manning equations, as shown in Equations ( 6)-( 8).
Here, Q is the volumetric flow rate, m 3 /s (cubic meters per second), L is length, D is diameter, and C is the pipe roughness coefficient in the Hazen-Williams equation.
Here, Hk is the head loss, f is the friction factor, L is length, D is diameter, V is velocity, and g is the gravity coefficient in the Darcy-Weisbach equation.
Here, V is the cross-sectional average velocity, n is the Gauckler-Manning coefficient, R is the hydraulic radius, and L is length in Manning equations.hmin ≤ hj ≤ hmax j = 1, 2, 3, . .., nn Here, hj, hmin, hmax, and nn denote the pressure head at node j, the minimum and maximum required pressure heads at node j, and the number of nodes in the network, respectively. vmin Here, vj, vmin, vmax, and np denote the velocity at pipe j, the minimum and maximum required velocities at node j, and the number of pipes in the network, respectively.
Here, Di is the diameter of pipe I, and {D} denotes the set of commercially available pipe diameters.In fact, all design variables (i.e., pipe diameter) should be chosen from the suggested commercial sizes (discrete optimization problem).

Studied Networks 2.2.1. Balerma Network
The Balerma network has four reservoirs, eight loops, 454 pipes, and 443 demand nodes, as shown in Figure 1.It was firstly presented by Reca and Martinez [50].The minimum nodal pressure requirement is 20 m at each node, and the pipe velocity requirements are not considered for the Balerma network.The Balerma network has four reservoirs, eight loops, 454 pipes, and 443 demand nodes, as shown in Figure 1.It was firstly presented by Reca and Martinez [50].The minimum nodal pressure requirement is 20 m at each node, and the pipe velocity requirements are not considered for the Balerma network.The discrete sizes of commercial diameters, accompanied by their respective cost per unit length, can be found in Table 1.Table 1 shows that each pipe can be selected from the 10 candidate diameters, resulting in a total of 10,454 potential network configurations.[48] as an irrigation water supply network in the Republic of Korea.While originally planned as a looped network, the Saemangeum network was also considered as branched and pumped types in this study.The branched network (Plan 1) consists of 335 pipes with a combined length of 37.88 km, shown in Figure 2a.On the other hand, the looped network (Plan 2) has 356 pipes; we converted some pipe configurations and generated 10 more loops from the branched network, as shown in Figure 2b.Furthermore, the pumped irrigation network The discrete sizes of commercial diameters, accompanied by their respective cost per unit length, can be found in Table 1.Table 1 shows that each pipe can be selected from the 10 candidate diameters, resulting in a total of 10,454 potential network configurations.[48] as an irrigation water supply network in the Republic of Korea.While originally planned as a looped network, the Saemangeum network was also considered as branched and pumped types in this study.The branched network (Plan 1) consists of 335 pipes with a combined length of 37.88 km, shown in Figure 2a.On the other hand, the looped network (Plan 2) has 356 pipes; we converted some pipe configurations and generated 10 more loops from the branched network, as shown in Figure 2b.Furthermore, the pumped irrigation network is made up of 345 pipes that are supplied through four pumping stations within the network (Plan 3, Figure 2c).The length of the pipes is 41.39 km, and based on four pumping stations, the whole network is divided into four sections.
Sustainability 2024, 16, x FOR PEER REVIEW 6 is made up of 345 pipes that are supplied through four pumping stations within the work (Plan 3, Figure 2c).The length of the pipes is 41.39 km, and based on four pum stations, the whole network is divided into four sections.The data regarding the pipe cost per unit length for the pipes of different diam used in this study are shown in Table 2.In total, 18 types of commercial pipes were sidered in the optimization process.The data regarding the pipe cost per unit length for the pipes of different diameters used in this study are shown in Table 2.In total, 18 types of commercial pipes were considered in the optimization process.These were obtained from the 'Water Facilities Construction Cost Estimation Report', from K-water [51] in the Republic of Korea, which provides cost data on the construction costs for different types of pipes [48].In addition to nodal pressure, flux conditions of the pipe line were also considered as hydraulic conditions for the Saemangeum networks.The reason why these water pressure and water flux standards exist in pipe network optimization problems is because water supply stability can be problematic if the water pressure is low or too high, or if the water flux is slow or very fast.Therefore, the purpose of optimal water pipe network design is to determine the pipe diameter that supplies water at the minimum cost within the appropriate water pressure and water flux range.In the case of nodal pressure, the minimum value and the maximum value are set as 10 m and 35 m, respectively.In the case of flux conditions of the pipe line, values should be maintained at between 0.01 m/s and 2.5 m/s.

Description of the Proposed Algorithm
A maisonette is defined as an apartment that is usually part of a larger building with two levels, often having its own entrance from the outside.The multi-floor structure has the advantage of making efficient use of space, because two floors can be utilized, although they are independent residential units.That is, the first floor and the second floor have the characteristic of being independent from each other, while still being mutually linked.Accordingly, the MTHSA in the present study, based on this structure, was developed as an optimization technology.Figure 3 shows the structure and components of the MTHSA.That is, multiple rooms can be arranged on the first floor, and independent spaces different from those on the first floor can be separately organized on the second floor.Therefore, in the case of the MTHSA, multiple Sub Harmony Memories (SHMs) are organized on the first floor.On the second floor, separate spaces are formed.Here, excellent solutions created in individual rooms on the first floor, and additional HSAs that can be implemented through the solutions, can be separately stored.That is, unlike existing HSAs, this structure increases optimization efficiency through multiple (parallel) searches on the first floor and additional searches among excellent solutions on the second floor.In the MTHSA, the Harmony Memory is divided into separated Sub Harmony Memories, and an independent search is performed in each Sub Harmony Memory, thereby improving the local search performance (exploitation).In addition, based on the best solutions of Sub Harmony Memories on the first floor, the solution memory of the second floor is formed, and further operations are performed, thereby improving the global search performance (exploration) in the MTHSA.
on the first floor and additional searches among excellent solutions on the second floor.In the MTHSA, the Harmony Memory is divided into separated Sub Harmony Memories and an independent search is performed in each Sub Harmony Memory, thereby improv ing the local search performance (exploitation).In addition, based on the best solutions o Sub Harmony Memories on the first floor, the solution memory of the second floor i formed, and further operations are performed, thereby improving the global search per formance (exploration) in the MTHSA.

Model Parameters
In the present study, five parameters were used.HMS, HMCR, PAR, and BW wer used as with HSAs, and Sub Harmony Memory Size (SHMS) was additionally employed HMS means the total number of Harmony Memories used to solve optimization problems The HMCR is a user parameter that shows the probability of randomly generating new solutions, or selecting new solutions from existing Harmony Memories.If the value of th HMCR is 0.7, this means that the probability of generating new solutions from existing HMs is 70%, and the probability of random selection of new solutions is 30%.The PAR i the probability of adjusting the value of each decision variable of new solutions within th designated range (BW).If the value of PAR is 0.1, the probability of adjusting the value o each decision variable is 10%.The designated range used in this case was determined by the bandwidth.In the case of water distribution network problems, because discrete var iables for selecting one of the available sizes of commercial pipes were used, the BW i defined as increasing or reducing pipe diameters by one level.SHMS means the size o rooms arranged on the first floor of the MTHSA.For instance, if the value of HMS is 6 and the value of SHMS is 20, three HMs of size 20 will be formed on the first floor.

Procedure of the Proposed Algorithm
A pseudo code that shows the process of implementing the MTHSA is shown in Al gorithm 1.First, the values of the HMS, SHMS, HMCR, PAR, and BW parameters of th MTHSA are set.
According to the set sizes of HMS and SHMS, the number of rooms on the first floo and the number of excellent solutions to be organized on the second floor are determined That is, if the values of HMS and SHMS are 60 and 20, respectively, the number o rooms on the first floor and the excellent solutions to be organized on the second floor wil both be three.The next stage is to perform parallel searches on the first floor.Using th three Harmony Memory sets organized on the first floor, HSAs are implemented

Model Parameters
In the present study, five parameters were used.HMS, HMCR, PAR, and BW were used as with HSAs, and Sub Harmony Memory Size (SHMS) was additionally employed.HMS means the total number of Harmony Memories used to solve optimization problems.The HMCR is a user parameter that shows the probability of randomly generating new solutions, or selecting new solutions from existing Harmony Memories.If the value of the HMCR is 0.7, this means that the probability of generating new solutions from existing HMs is 70%, and the probability of random selection of new solutions is 30%.The PAR is the probability of adjusting the value of each decision variable of new solutions within the designated range (BW).If the value of PAR is 0.1, the probability of adjusting the value of each decision variable is 10%.The designated range used in this case was determined by the bandwidth.In the case of water distribution network problems, because discrete variables for selecting one of the available sizes of commercial pipes were used, the BW is defined as increasing or reducing pipe diameters by one level.SHMS means the size of rooms arranged on the first floor of the MTHSA.For instance, if the value of HMS is 60 and the value of SHMS is 20, three HMs of size 20 will be formed on the first floor.

Procedure of the Proposed Algorithm
A pseudo code that shows the process of implementing the MTHSA is shown in Algorithm 1. First, the values of the HMS, SHMS, HMCR, PAR, and BW parameters of the MTHSA are set.
According to the set sizes of HMS and SHMS, the number of rooms on the first floor and the number of excellent solutions to be organized on the second floor are determined.
That is, if the values of HMS and SHMS are 60 and 20, respectively, the number of rooms on the first floor and the excellent solutions to be organized on the second floor will both be three.The next stage is to perform parallel searches on the first floor.Using the three Harmony Memory sets organized on the first floor, HSAs are implemented separately.In this case, the HMCR and PAR are used to create new solutions.The new solutions and the objective function of the worst solution in each set are then compared with each other, and if the values of the new solutions are better, the location of the Harmony Memory will be changed.The next stage is to perform searches among excellent solutions on the second floor.The best solutions in individual sets obtained on the first floor become the components of the initial Harmony Memories on the second floor.The HSA is implemented using the Harmony Memories on the second floor, organized as such.The process of optimization in the MTHSA is iterated until the terminating condition (e.g., the number of iterations) is satisfied.In general, whereas one new solution is generated at each iteration in the HSA, the number of new solutions generated at each iteration varies with the values of HMS and SHMS in the case of the MTHSA.That is, when the values of HMS and SHMS are 60 and 20, respectively, the HSA is implemented three times on the first floor and one time on the second floor in each iteration, so that a total of four new solutions are generated.

Sensitivity Analysis of Parameters
Before verifying the applicability of the MTHSA, an analysis of the sensitivity of each user parameter was carried out, and the results are presented here.The sensitivity analysis was applied to the Balerma network.In order to have a fair comparison, 30 runs were implemented for each parameter combination.The Number of Function Evaluations (NFE) of all simulations was set to 45,400.Here, the NFEs are general calculation termination conditions for evaluating the performance of each technique in the field of optimization under the same conditions.
In order to analyze the sensitivity of each parameter, among various parameter combinations, a combination of 60 HMSs, 20 SHMSs, a PAR of 0.02, and an HMCR (MTHSA (1)) of 0.90, which produces relatively good results, was set as a reference combination.The results of the sensitivity analyses are presented in Tables 3-6.Tables 3 and 4 show the sensitivity analysis results for HMS and SHMS.
It can be observed that the smaller HMS values and the larger SHMS values show better results.One thing that must be noted is that a ratio between HMS and SHMS (HMS/SHMS) that is closer to one shows better results.Since the maximum NFEs have been determined, the larger the value of the ratio between HMS and SHMS, the more times function evaluations are implemented per time of iteration.In other words, those types of the MTHSA that can be advantageous in terms of the number of implementations of function evaluation, by performing searches among excellent solutions on the second floor while reducing the number of solution searches on the first floor, show relatively higher optimization efficiency.However, it is judged that the ratio between HMS and SHMS (HMS/SHMS) that shows the highest efficiency may vary with the value of the Maximum NFEs.Therefore, in the case of the MTHSA, determining appropriate ratios between HMS and SHMS is judged to be an important criterion.Tables 5 and 6 show the results of HMCR and PAR sensitivity analysis.It can be observed that relatively larger HMCR values and relatively smaller PAR values show better results.In the case of HMCR, the larger the value, the higher the probability of generating new solutions in the existing Harmony Memory.Therefore, applying HMCR with a large value means generating a new optimal solution on the second floor by using the information of the best solutions of the first floor with a high probability, and therefore the optimal solution can be found faster.It can be seen that in the case of PAR, because the optimization problem involves a large number of decision variables (i.e., 454 decision variables), large PAR values act as an element that obstructs the convergence of solutions.In particular, it can be seen that deviations in optimal solutions determined by the value of PAR appear much larger compared to other parameters.Therefore, in the case of large-scale problems with large numbers of decision variables and alternative solutions, in order to obtain reasonable optimal results, using relatively lower PAR values is judged to be desirable.

Results of Balerma Network
Based on the results of the sensitivity analyses conducted Section 4.1, the results for the optimal design of the Balerma network were organized and analyzed, as presented in the following.Table 7 shows the combinations of parameters and the NFE used to compare and analyze the results.MTHSA (I) and MTHSA (II) denote those parameter combinations that showed excellent results among the sensitivity analysis results.MTHSA (III) shows a case with the same parameters as those used in MTHSA (I), but a different value for the NFEs.
Table 8 presents those results from the existing literature that are different from the results of the MTHSA.The optimal results of the GA, SA, MSATS [36], the HAS (1), HAS (II), PSHS [41], the genetic heritage evolution by stochastic transmission (GHEST) [52], the MBA, IMBA (I), IMBA (II) [43], the decomposition and multistage approach for optimizing the design of WDSs using the DE (Decompose-DE) [53], the GANOME [50], the HD-DDS-1, HD-DDS-2 [54], and other optimization techniques from the literature are presented.Table 8 presents the results, organized in ascending order of the Maximum NFE values used in individual pieces of the literature.First, in the case of Best Cost, Decompose-DE and NLP-DE2 show the least cost, at a value of 1.923 (106 EUR).It can be observed that the MTHSA (III) developed in the present study shows almost the same value, as it yields a value of 1.932 (106 EUR).In the case of the Decompose-DE method, before applying DE, pre-processing is implemented based on the engineering characteristics of water distribution networks.The pre-processing process used here should be preceded by gaining a sufficient understanding of water distribution networks.Therefore, the MTHSA has an advantage in that it shows sufficiently good results even without pre-processing based on the engineering knowledge of water distribution networks.Unlike the MTHSA, the NLP-DE2 is a combination of non-linear programming and DE, and it confirms that the MTHSA is efficient, as good optimal solutions were obtained when only meta-heuristic techniques were used in the case of the MTHSA.
On reviewing the results, when the value for the Maximum NFEs was fixed at 45,400, it can be seen that the proposed method and the IMBA show excellent results.For the case of the Best Cost, IMBA (I) showed that it can design solutions at a cost lower by 0.021 (106 EUR) than for MTHSA (II).However, the Average and Worst Costs of the proposed algorithm were shown to be lower by approximately 0.030 (106 EUR) and 0.063 (106 EUR), respectively, compared to IMBA (I).This means that the proposed algorithm achieved better optimal solutions on average compared to the IMBA.
When reviewing the results of HSA (I) and HSA (II), which are based on the HSA and PSHS, it can be seen that when a Maximum NFE value of 45,400 was applied, the MTHSA achieved much better results.Whereas a value of 2.018 (106 EUR) could be obtained by applying an NFE of 10,000,000 in the case of HSA (II), a value of 1.932 (106 EUR) was obtained by choosing 1,000,000 as the NFE in the case of MTHSA (III).
Figure 4 shows the results (nodal pressure) of a numerical analysis of WDNs, when the optimal pipe diameter combination from MTHSA (I) is chosen.The node that shows the closest result to the minimum required pressure of 20 m is no.55 (20.01 m), and it can be seen that this value satisfies the minimum required pressure.In particular, in the case of WDNs, an even pressure distribution is very important in terms of leakage management.Figure 4b shows the cumulative distribution of the pressure of each node, appearing as a result of the optimal design, from which an even distribution without any particular unequally distributed pressure section can be identified.
the closest result to the minimum required pressure of 20 m is no.55 (20.01 m), and it can be seen that this value satisfies the minimum required pressure.In particular, in the case of WDNs, an even pressure distribution is very important in terms of leakage management.Figure 4b shows the cumulative distribution of the pressure of each node, appearing as a result of the optimal design, from which an even distribution without any particular unequally distributed pressure section can be identified.

Results of Saemangeum Network
Yoo et al. [48] applied the HSA to the Saemangeum network for the first time, in order to present optimal design results.In the present study, the results of the application of the proposed MTHSA algorithm were compared with the results of Yoo et al. [48] and the results of that comparison are presented.The parameters used in the proposed algorithm were an HMS of 60, SHMS of 20, PAR of 0.03, and HMCR of 0.90.Note that the parameters were determined through a sensitivity analysis for the Saemangeum network.The Maximum NFE was set to 15.5 (106), in consideration of the scales and characteristics of the problems, although Yoo et al. [48] did not present their values.
Table 9 presents the results that Yoo et al. [48] obtained by applying the HSA, and the results of the MTHSA.According to the application, the results for the MTHSA were of a higher quality for all three design plans.
That is, the MTHSA is shown to be capable of designing optimal pipe diameters that satisfy the hydraulic conditions, at costs that are lower by approximately 1.77% to 2.79%.In the case of a real water distribution network design, which requires large construction costs, a small percentage of cost reduction can reduce the construction costs by a large amount.In particular, it is observed that the MTHSA can provide economic saving effects amounting to approximately 300 million KRW in the case of the Looped Type (Plan 2), in which relatively more diverse water supply paths can appear.
Table 10 present the arithmetic statistical values obtained when 30 independent runs were implemented for each plan.It can be observed that the MTHSA achieved even values for all plans, without any large deviation.

Results of Saemangeum Network
Yoo et al. [48] applied the HSA to the Saemangeum network for the first time, in order to present optimal design results.In the present study, the results of the application of the proposed MTHSA algorithm were compared with the results of Yoo et al. [48] and the results of that comparison are presented.The parameters used in the proposed algorithm were an HMS of 60, SHMS of 20, PAR of 0.03, and HMCR of 0.90.Note that the parameters were determined through a sensitivity analysis for the Saemangeum network.The Maximum NFE was set to 15.5 (106), in consideration of the scales and characteristics of the problems, although Yoo et al. [48] did not present their values.
Table 9 presents the results that Yoo et al. [48] obtained by applying the HSA, and the results of the MTHSA.According to the application, the results for the MTHSA were of a higher quality for all three design plans.That is, the MTHSA is shown to be capable of designing optimal pipe diameters that satisfy the hydraulic conditions, at costs that are lower by approximately 1.77% to 2.79%.In the case of a real water distribution network design, which requires large construction costs, a small percentage of cost reduction can reduce the construction costs by a large amount.In particular, it is observed that the MTHSA can provide economic saving effects amounting to approximately 300 million KRW in the case of the Looped Type (Plan 2), in which relatively more diverse water supply paths can appear.
Table 10 present the arithmetic statistical values obtained when 30 independent runs were implemented for each plan.It can be observed that the MTHSA achieved even values for all plans, without any large deviation.
In particular, as presented in the results of Yoo et al. [48], it can be seen that in the case of Plan 2 (Looped Type), in which more water supply paths can appear, relatively larger deviations in the minimum and maximum costs appeared compared to Plans 1 and 3.However, the value of the maximum cost of the MTHSA was shown to be smaller than or similar to that of the HSA, and this result can be interpreted as showing the optimization efficiency of the MTHSA.

Conclusions
In the present study, an improved version of the HSA, reinforced further in terms of diversification and intensification, was proposed.The MTHSA configures a two-storied structure, and it increases the optimization efficiency by implementing multiple searches on each floor and additional searches among the better solutions.In order to verify the proposed algorithm, it was applied to pipe diameter optimization problems for WDNs, which constituted representative discrete optimization problems.In particular, the algorithm was applied to large-scale WDNs, with at least 300 decision variables to verify in actual problems.First, through analyses of the sensitivity of the MTHSA, PAR was identified as the most sensitive parameter.However, it was observed that optimization efficiency could be increased by adjusting the ratio between HMS and SHMS (HMS/SHMS).In general, in the case of small NFEs, HMS/SHMS values closer to one showed a faster and more efficient convergence.However, when HMS was set to three times the value of SHMS, efficient optimization results could be obtained, on average.In the results, for the Balerma network optimal design problem, the developed MTHSA algorithm was found to reduce the minimum cost by about 19.85% compared to the existing HSA under the same termination conditions (45,400 NFEs).In addition, when the optimal solution search calculation was continued, it was confirmed that the MTHSA showed better results even under the condition of NFEs equivalent to 10% of the existing HSA.It could be identified that in the optimal design of pipe diameters, the MTHSA was advantageous over existing algorithms, as it does not involve a pre-processing process that requires expertise, as well as not using complicated mathematical technology.In particular, when compared to the existing literature in which only meta-heuristic algorithms were used, the MTHSA showed considerably more efficient optimization results.In the Saemangeum network, it was also confirmed that the newly developed algorithm showed a superior performance compared to the existing HSA.The results showed that design costs were reduced by 1.77 to 2.29% in all three types of Saemangeum networks.Additionally, with the MTHSA, the optimization calculation structure can be easily changed by adding only one parameter (SHMS).
However, the MTHSA developed here should be applied to diverse optimization problems (e.g., numerical benchmark problems, real-world problems), in order to further verify the quality of the algorithm.Single-objective problems were considered for the optimal cost design in this study; however, the MTHSA can also applied to WDNs with multi-objective features, which are more complex than single-objective problems.Therefore, in future studies, the MTHSA will be applied to multi-objective design problems for WDNs, and also other engineering problems with various problem characteristics.In addition, if technologies to reduce or automatically adjust the number of parameters using selfadaptation techniques considering problem characteristics are developed based on the results of the sensitivity analysis, the algorithm will become more practical.

Figure 3 .
Figure 3. Schematic structure and components of the MTHSA.

Figure 3 .
Figure 3. Schematic structure and components of the MTHSA.
(a) Contour plot for pressure (b) Frequency plot for pressure

Figure 4 .
Figure 4. Nodal pressure distribution of optimal results in MTHSA (I).

Figure 4 .
Figure 4. Nodal pressure distribution of optimal results in MTHSA (I).

Table 1 .
Candidate pipe diameters and costs for the Balerma network.

Table 1 .
Candidate pipe diameters and costs for the Balerma network.

Table 2 .
Commercial pipe diameters and their costs per length for the Saemangeum networks.
Max number of iterations) Generate new harmonies by accepting best harmonies in first floor Adjust pitch to get new harmonies (solutions) if (rand < HMCR), choose an existing harmony randomly in each of sub harmony memories else if (rand < PAR), adjust the pitch randomly within limits else generate new harmonies via randomization end if Accept the new harmonics (solutions) if better Gathering best harmonics in each of sub harmony memories Generate new harmony by accepting best harmony in second floor Adjust pitch to get new harmony (solution)

Table 3 .
Effects of the MTHSA parameter (HMS) for the Balerma network.

Table 4 .
Effects of the MTHSA parameter (SHMS) for the Balerma network.

Table 5 .
Effects of the HMCR parameter used in the MTHSA for the Balerma network.

Table 6 .
Effects of the PAR user parameter used in the MTHSA for the Balerma network.

Table 7 .
Initial parameters used for the Balerma network.

Table 8 .
Numerical optimization results for the Balerma network.

Table 9 .
Optimal design results and cost comparison with the previous study for the Saemangeum networks.

Table 10 .
Numerical optimization results for the Saemangeum networks.