Hybrid particle swarm optimization and tabu search for the design of large-scale water distribution networks

Water distribution network (WDN) optimization has received special attention from various technicians and researchers, mainly due to its high costs of implementation, operation and maintenance. However, the low computational efficiency of most developed algorithms makes them difficult to apply in large-scale WDN design problems. This article presents a hybrid particle swarm optimization and tabu search (H-PSOTS) algorithm for WDN design. Incorporating tabu search (TS) as a local improvement procedure enables the H-PSOTS algorithm to avoid local optima and show satisfactory performance. Pure particle swarm optimization (PSO) and H-PSOTS algorithms were applied to three benchmark networks proposed in the literature: the Balerma irrigation network, the ZJ network and the Rural network. The hybrid methodology obtained good results when seeking an optimal solution and revealed high computational performance, making it a new option for the optimal design of real water distribution networks


INTRODUCTION
Water distribution networks (WDNs) are hydraulic systems composed of reservoirs, pumps, pipes, valves, sensors and other accessories designed to transport drinking water with sufficient flow and pressure to meet consumers' water needs continuously and appropriately. WDNs are an important part of the water supply system, mainly due to their high cost of implementation, operation and maintenance.
WDN designers seek to assign pipe diameters to minimize investment costs, including those related to implementing and purchasing materials and equipment. Traditionally, engineers design pipe networks using trial and error guided by experience. However, the state variables of the problem are linked by the mass and energy conservation laws of physics, which make the design of large systems complex and exhaustive.
The term "optimization", present in several engineering studies and projects, refers to the search for the best solution to a problem. WDN design is classified as a large combinatorial discrete nonlinear nondeterministic polynomial-hard (NP-hard) optimization problem (Moosavian & Lence, 2020). Such problems are nonlinear and non-convex, and global optimization techniques cannot usually solve them. Therefore, the development of new models makes a substantial contribution to solving this type of problem (Cassiolato et al., 2020).
The evaluation of complex real-world problems using large search spaces is currently the main difficulty addressed in the specialized literature. The characterization of the size of a WDN is dictated by its number of nodes and meshes, due to the growing number of independent and variable equations to be determined during its dimensioning, in addition to the iteration between them, which makes the problem more complex and difficult to solve . In this context, this work presents a hybrid particle swarm optimization and tabu search (H-PSOTS) algorithm used to design large-scale WDNs. In recent years, the hybridization of these techniques has been applied in several areas (Chentoufi & Ellaia, 2018;Tang & Lee, 2018;Ahmadian et al., 2019;Alharkan et al., 2020;Ebadi & Jafari Navimipour, 2019;Toshev, 2019), but a literature search revealed no such applications to the proposed problem type.
Particle swarm optimization (PSO) is used widely in WDN design problems due to its ease of implementation and the good results reported in the literature. At each iteration of the PSO process, numerous particles explore the search space to find the best results for the objective function. During this process, the particles can resume the positions visited and revisit the positions already reached by other. In large problems, this iterative process results in high computational cost and execution time. Therefore, a tabu search (TS) algorithm was inserted into the PSO process to restrict repetitive movements. The TS algorithm gained popularity in the scientific community as an auxiliary adaptive procedure, as it improves the exploration of local search spaces and facilitates greater solution diversification in the search of the global optimum.

METHOD
The H-PSOTS algorithm was implemented using Matlab software. According to Chandramouli (2019), Matlab is best suited for all mathematical operations and can link external libraries. The EPANET input (.inp) file is first read to create the network model and obtain all network parameters (e.g. layout, node demands, pipe lengths). Matlab was combined with the functions of the EPANET toolkit to simulate the various solutions indicated by the optimization algorithm. On Windows systems, the library is precompiled into a dynamic link library file called EPANET2.dll. All analyses were performed on a personal computer with an Intel Core i7 1.88-GHz processor and 8 GB of RAM.

Optimal design of WDN
In this article, WDN design was formulated as a least-cost optimization problem with a selection of pipe diameters as the decision variables. In general, WDN optimization was defined as follows: Minimize: Total cost of the WDN Subject to: 1. Conservation of mass: Inflows and outflows must balance at each node; 2. Conservation of energy: Head loss around a closed loop must be equal zero, and head loss along a path between two reservoirs must be equal to the reservoirs' elevation difference; 3. Minimum pressure requirements: Minimum pressure must be provided at network locations for a given set of demands; 4. Acceptable pipe sizes: Diameters must be selected from an admissible set.
The objective function (Equation 1) was assumed mathematically to be a cost function of pipe diameters and lengths plus the penalty. Whenever the pressure at a node violated the minimum pressure requirement, a penalty was added to the total cost. This penalty made solutions that did not meet the minimum pressure requirement infeasible or infactible. where ( ) i f D is the total cost of the network, ( ) i f D is the cost in pipe diameter per unit length, i L is the pipe length, m and n are the number of pipes and nodes of the network, and j PP is a penalty function intended to ensure that pressure constraints are satisfied.
During the optimization process, EPANET received the set of new diameters as input and returned the pressures on the network nodes to the algorithm. The head loss equation depended on the problem evaluated; available choices included the Darcy-Weisbach, Hazen-Williams and Chezy-Manning equations. The original information from the case studies was maintained so the results could be compared with the works published in the literature. In addition, the optimization models included integer variables, and no approximations were used to find the diameter values in the discrete set of commercially available diameters.

Particle Swarm optimization
The PSO algorithm is widely used because of its simplicity and attractive search efficiency. Kennedy & Eberhart (1995) developed the parallel evolutionary computation technique based on the social behavior metaphor to solve unconstrained optimization problems (Yang et al., 2019). In the PSO algorithm, a set of particles, called a population or swarm, is used to "research" for the best solution to the optimization problem. The search space is dictated by the objective function and the restrictions imposed, and each particle represents a unique set of solutions to the problem. The particle movement cognitive system allows the best solutions to be identified more quickly, decreasing the computational cost of the exploration process. However, a larger swarm size, dictated by the complexity of the problem, results in increased execution time.
The particles' movement is associated with Gbest, the best position reached within the search space for the swarm particles, and Pbest, the best individual position reached within the search space for each particle. In successive iterations, the particle's new velocity is calculated using its previous velocity (weighted by an inertia factor) and the distance from its current position to its Pbest position and the Gbest position (weighted by social and cognitive factors), as shown in Equation 2. To manage Gbest or Pbest's influence on particle motion, Kennedy & Eberhart (1995) adopted two positive constants, called acceleration coefficients c1 and c2, which are directly linked to "social and cognitive factors". In this work, following the recommendation of these authors, value 2 was adopted for the coefficients. The current velocity vector is added to the previous particle position, t i X , giving rise to the new position, t 1 Where: i v and i X are the velocity and position of the i th particle, w is inertia weight that introduces the weighting of the current velocity on the particle in the next generation, i Pbest represents the position with the best fitness found so far by the i th particle, i Gbest is the position with the best fitness found so far by all the particles in the population, c 1 and c 2 are positive constants, 1 r and 2 r are random numbers between 0 and 1, t is the number of the iteration, and w is the inertia weight that controls the impact of the velocity history on the current velocity to influence the trade-off between global and local experiences. The inertia factor (Equation 4) is related to the portion of the velocity vector that is connected directly to the particles' movement. The balance between the global and local exploration of the particles is related to this factor's value. Eberhart & Shi (2000) proposed that the factor varies between an initial value of 0.9 and a final value of 0.4 in the iterations and decreases gradually. In addition, they noted that high factor values are more effective at the beginning of the search, exploration, because they accelerate convergence and spread the particles over a broader search space. Smaller values favor local exploration by improving the refinement of the search for the final solution.
Where: ini w and fin w are the initial and final inertia factors; and N is the total number of iterations. In order to determine the values of initial and final inertia factors, a sensitivity analysis was performed on the first network evaluated, Balerma irrigation network.
Iteration tolerance (IT) is one of the stopping criteria attributed to the PSO algorithm to prevent simulations from entering an infinite loop. This factor interrupts the optimization process when the number of iterations with the same result exceeds a certain proportion of the total number of iterations that have already occurred. In this work, an IT of 30% was adopted in all case studies. Additionally, it is highlighted that this work adopted a maximum number of iterations, in each simulation, equal to 1000. Glover (1989) proposed the concept of the tabu search (TS) heuristic. The basic principle of the TS is to pursue a local search by allowing non-improving moves whenever it encounters a local optimum. Five key elements (i.e. search space, neighborhood structure, tabu list, aspiration criteria and termination criteria) must be designed when developing a TS (Yang et al., 2017). A fundamental element of TS is using flexible memory functions (tabu list) to forbid transitions, called "moves", from the current solution to other previously visited candidate solutions. The list is updated with each move of the particles and remains in memory for a number of iterations. The longer the list, the more restrictive the search process and the longer the computer runtime. Another factor that improves the method's quality is restricting moves that guide the particles to new search space neighborhoods, avoiding premature convergence in a local optimum.

H-PSOTS hybrid algorithm
Inserting the TS algorithm into the PSO algorithm's search process increased particle movement efficiency because the TS improves local exploration and decreases premature algorithm convergence and computational cost. The convergence speed of TS depends on the initial solution, and the parallelism of the PSO population helps the TS find promising search space regions very quickly. Incorporating TS into a PSO algorithm as a local improvement procedure enables the algorithm to maintain the population's diversity (Shen et al., 2008). During the present study's exploration process, the algorithm reduced the repetition of executed movements, allowing a greater diversification of solutions in the search for the global optimum.
The initial position of the particles in the PSO process comprised the tabu list's movement restrictions. As the position of the particles was updated and before the calculation process in the hydraulic simulator and the objective function, the algorithm verified whether a particle of the swarm had already reached the position. If it was a new position, the particle did not move, the swarm update process continued and the list received that position as an increment. However, if the particle's position was already on the tabu list, it was updated. This process was repeated until the whole swarm was updated when its particles visited new positions. When every particle in a swarm has been updated, the list can be updated equally with the new positions, depending on whether the particles were able to keep the previous ones. The number of times or positions stored in the list was determined based on a sensitivity analysis carried out in one of the case studies. To guarantee algorithm convergence regardless of the situation, the Gbest position was not part of the list.
In addition to the tabu list size, the aspiration criterion was defined. This criterion removed the tabu list status from particles with high values in the final optimization phase (close to the maximum number of iterations or iteration tolerance). Thus, the particles had more freedom to seek solutions, preventing the algorithm from stagnating in a local optimum. This criterion directed the particles to the best identified cost, returning solutions that could be useful in swarm diversification. Another criterion adopted for the list was related to the particles that presented the same position as Gbest, facilitating better particle positions and the convergence of the algorithm. In all case studies, the criteria for stopping the search were an iteration tolerance of 30% and a maximum of 1000 iterations. The general flowchart of the H-PSOTS algorithm is shown in Figure 1.
Because evolutionary algorithms are not deterministic, a sensitivity analysis was performed to define the number of particles in the PSO algorithm and the tabu list size. The other input parameters of the PSO algorithm were adopted based on the literature and remained as a standard configuration (default) for all case studies. The adoption of specific parameters for each case study would present better results, but a single algorithm configuration was chosen to determine a typical configuration suitable for any network.

RESULTS AND DISCUSSION
For comparison, the PSO and H-PSOTS algorithms were applied to three complex WDNs referenced in the literature. The first network studied was the large irrigation system located in Balerma (Reca & Martínez, 2006), in the province of Almeria (Andalusia, Spain). The second, known as the ZJ network (Zheng et al., 2011), was a provincial water supply network in China. The last case study, the Rural network, was also a complex irrigation system, initially presented by Marchi et al. (2014). Network layouts are shown in Figure 2.

Balerma irrigation network
The Balerma irrigation network (Reca & Martínez, 2006) is an adaptation of the WDN for the irrigation system in the Sol-Poniente district, located in Balerma, province of Almeria, Spain. This irrigation network includes 454 PVC pipes, 443 nodal demands, eight loops and four reservoirs. The diameter set consisted of 10 diameter types ranging from 113 to 581.8 mm: D(mm) = {113. 0, 126.6, 144.6, 162.8, 180.8, 226.2, 285.0, 361.8, 452.2, 581.8}. Table 1 presents the costs of the specified diameters.  The minimum pressure required at each node was 20 meters, and pipe friction head loss was computed for all pipes using the Darcy-Weisbach equation with a roughness coefficient of 0.0025 mm.
The algorithm was applied to this network with two objectives: (a) to determine the parameters: inertia factors, number of particles, and the tabu list size, which defined the number of restrictive positions that the particles occupied in each iteration, and (b) to verify the proposed algorithm's performance. The complexity of the WDN design problem is related to the number of sections in the network and the diameters available for the design. Therefore, this network has a search space with 10 454 possible solutions (viable and non-viable).
To determine the inertia factor (Equation 4), the PSO algorithm was applied in 50 simulations for different ini w and fin w values, keeping the number of particles (100 particles) and the social and cognitive factors (equal to two) constant. The ini w and fin w values ranged from 0.3 to 1.0. Figure 3 shows the results of the sensitivity analysis, where the center of each circle indicates the average values of runtime and cost obtained in 50 simulations, while the diameter represents the standard deviation of costs. Some simulations with values less than 0.3 and greater than 1.0 were performed; however, it was decided not to include them in the text, as they obtained premature convergence in local minimums or high simulation runtime. Thus, this work adopted 0.4 and 0.9 for the ini w and fin w values, respectively. With the inertia factor defined, the next step was to determine the number of particles used during the optimization process. For this, several simulations were performed with the PSO algorithm by varying the number of particles from 10 to 200. The best values of the simulations' results are shown in Figure 4. The best solution, € 2.10 million, remained stable from 160 particles. Therefore, this value was adopted for the H-PSOTS algorithm because the increase in particles would result in a higher computational cost and a consequent increase in execution time. The relationship between the number of particles and the number of variables (pipe diameters) was 35%. This relationship was maintained in the other networks studied, even though increasing the number of particles can yield better results.
To analyze and determine the hybrid algorithm's tabu list size, simulations were performed in which the list retained up to four previous positions for each particle until the complete swarm update. The best values of the analysis for each the list size increment are shown in Figure 5. It is noteworthy that, regardless of the tabu list size, each new position was added to the list during the particle update process. In addition, at the end of the swarm update, the newly generated positions replaced were added to the old ones, depending on the choice adopted.
Increasing the tabu list worsened the algorithm's solution execution time, demonstrating that a simple list is the best option. The time increase resulted from a growth in the number of checks performed on the list during each particle's updates. In terms of cost, large lists limited the swarm's diversification, which is, in many cases, essential for the improvement of results in the process's advanced stages. During advanced stages of simulation, this diversification helped to decrease the chances of premature convergence in local minima.   After all parameters were determined, the H-PSOTS algorithm was applied, and its results compared with those of the PSO algorithm and the literature (Zheng et al., 2011;Sheikholeslami et al., 2016;Bi et al., 2020) (Table 2). The H-PSOTS algorithm obtained a minimum cost of € 1.998 million, close to the best-known minimum cost of € 1.923 million found by Zheng et al. (2011) using NLP-DE2. Although the result obtained by hybrid algorithm might not represent the best cost in the literature, it determined a good solution with a low simulation runtime in relation to the PSO algorithm. Figure 6 shows the evolution of the costs of the best solutions found by algorithms PSO and H-PSOTS during the optimization   processes, while Figure 7 shows the graphical representations of the results obtained (costs and simulation runtimes) in the 100 simulations performed (50 for each algorithm). The H-PSOTS algorithm determined the optimum solution with a simulation runtime of 19 minutes, while the PSO algorithm required 43 minutes. The H-PSOTS algorithm resulted in greater exploration of the search space in a shorter execution time because calculation processes were avoided when moving the particles during each iteration. The algorithm's runtime offers the designer more flexibility to scale real systems, and it offers a set of viable solutions with similar prices and reliability. Figure 8 shows the diameters of the optimized network and the nodal pressures obtained by the H-PSOTS algorithm.

ZJ network
The ZJ network, originally described by Zheng et al. (2011), is a real network in an eastern province of China. It comprises 50 loops, 164 pipes, 113 nodes and a single reservoir. The minimum pressure required at each node was 22 meters, which can be considered the lower pressure limit of the WDN's nodes. Pipe friction head loss was computed using the Hazen-Williams equation, whose coefficients were ω = 10.667, α = 4.871 and β = 1.852; the dimensionless roughness coefficient C was 130 for all pipes. Fourteen potential pipe diameters between 150 and 1000 mm were available for the network design (Table 3); therefore, the search space for this optimization problem consisted of 14 164 = 9.226 × 10 187 possible solutions.
The number of particles adopted in this case study was 58, which corresponds to 35% of the number of pipes in the network. This follows the relationship determined in the analysis of the Balerma network. The results of the design using the PSO and H-PSOST algorithms and the works in the literature are shown in Table 4.   The value obtained using the H-PSOTS algorithm was $ 7.70 million, which corresponds to a reduction of approximately 7.5% compared to the PSO algorithm. The current best-known solution for the ZJ network reported in the literature is $ 7.082 million, as reported by Zheng et al. (2011), which was preceded by Bi et al. (2020), with a solution of $ 7.431 million. Figure 9 shows the evolution of the costs of the best solutions found by algorithms PSO and H-PSOTS during the optimization processes, while Figure 10 shows graphical representations of the results obtained (costs and runtimes) in the 100 simulations performed. For the best results obtained, the simulation runtimes were 20 and 24 minutes for the H-PSOTS and PSO algorithms, respectively. As expected, the H-PSOTS algorithm presented a more efficient search, which resulted in a reduction in execution time of about 20% compared to the PSO algorithm. This decrease occurred because of the reduced number of mathematical operations and iterative processing between Matlab and EPANET.   Optimal WDN design presents a good configuration of diameters and nodal pressures, with gradual variations in diameters. In other words, it always presents the flow direction from the largest diameter to the smallest. Figure 11 shows the network diameters and nodal pressures obtained using the H-PSOTS algorithm.

Rural network
The Rural network is a real WDN based on an irrigation system and consists of 379 junction nodes, 98 loops, and two reservoirs linked by 476 PVC pipes. Because it contains numerous loops, the Rural network is a good representative example of a highly complex WDN. The network was initially optimized by Marchi et al. (2014) using the genetic algorithm, PSO and differential evolution algorithms. The minimum pressure required at each node is zero, and the Darcy-Weisbach equation was adopted to calculate head loss using a roughness of 0.012 mm for all pipes.
The available set of diameters used to design the network comprised 15 possible pipe diameters (Table 5). The large number of possible diameters increased the search spaces, which presented 15 476 (6.61 × 10 559 ) possible solutions, making it a good challenge for the algorithm proposed in this research.
The swarm consisted of 167 particles in the simulations with the PSO and H-PSOTS algorithms, and the results are shown in Table 6. Although the PSO algorithm had several iterations with values close to those of the H-PSOTS algorithm, its solution and execution time were worse. Figure 12 shows graphical representations of the results obtained (costs and simulation runtimes) in the 100 simulations performed, while Figure 13 shows the evolution of the costs of the best solutions found by algorithms PSO and H-PSOTS during the optimization processes.
The H-PSOTS algorithm reached values close to the current best-known solution ($ 30.99 million) in 29 minutes, with an execution time that was 19% shorter than that of the PSO algorithm (36 minutes). The algorithms' evolutionary behavior profile is similar, with small evolutions at each iteration. Figure 14 shows the pipe diameters and nodal pressures of the network optimized using the H-PSOTS algorithm.
The methodology proposed in this study is aimed at improving the design of WDNs based on the execution time and the quality of the results. This technique, which involves two high performance algorithms in applications of this type, creates a set of solutions that meet the water demands and minimum pressures at a minimum cost and in an acceptable execution time.