Optimization of water distribution networks using simulated annealing (SA) with variable neighbourhood search (VNS): a case study in the city of Malard

Urban water consumption is rising globally, and water distribution network systems (WDNs) represent a significant portion of city infrastructure costs. The pipe distribution network of a water supply can constitute a major part of the project's capital cost. Numerous arrangements of pipe sizes within the network may exist to satisfy field requirements. Each arrangement results in a different project total cost. Therefore, attempts to obtain optimal solutions have involved hydraulic analysis of the network and the application of optimization techniques. The optimal solution aims to minimize the network's cost. Efficient WDN design is crucial for meeting demand and maintaining pressure while minimizing costs.

Several optimization models have been reported for the optimal design of water distribution systems. Two different optimization techniques are applied to solve these models: the deterministic optimization techniques (including linear, dynamic, and non-linear programming) and the stochastic optimization techniques such as Evolutionary algorithms (EAs). A deterministic approach for solving WDN typically involves using methods based on mathematical modeling and algorithmic techniques where the outcomes are entirely predictable and not influenced by random variables (El-Ghandour & Elbeltagi 2018). In contrast, a probabilistic approach for solving WDN incorporates uncertainty and randomness into the decision-making process. Unlike deterministic methods, probabilistic approaches recognize that certain parameters (e.g. water demand, pipe failures, pressure variations) are subject to variability, which is modeled and addressed through statistical or stochastic techniques.

EAs are fundamentally based on exploratory search and natural processes or even artificial intelligence. EAs, inspired by natural processes, have been widely applied for this purpose. Among these, genetic algorithms (GAs) were pioneering methods and have evolved into techniques such as NSGA-II for multi-objective optimization (El-Ghandour & Elbeltagi 2018).

Other notable algorithms include differential evolution (DE), shuffled complex evolution (SCE), particle swarm optimization (PSO), ant-colony optimization (ACO), harmony search (HS), and simulated annealing (SA). Each algorithm has unique strengths in addressing various optimization problems within WDNs (El-Ghandour & Elbeltagi 2018).

The latest research in the field of WDN optimization has been conducted by Sampaio Caradot Guilbert & Parez (2024). That study aims to employ a pipe deterioration model based on Markov chains, with transition matrices estimated from survival curves for different pipe cohorts. The proposed approach seeks to determine the appropriate levels of investment (CAPEX) and operational expenses (OPEX) levels in the coming decades. It was tested with real-world data from a sewer network in Sofia, Bulgaria, and the results show that it provides efficient long-term rehabilitation plans.

Despite these advancements, integrating variable neighborhood search (VNS) with other algorithms has shown promise but remains underexplored in WDN optimization. VNS is a framework for building heuristics based on systematic changes in neighborhoods, both in the descent phase, to find a local minimum, and in the perturbation phase, to emerge from the corresponding valley. A review of the literature revealed that no solution method exists that combines a heuristic algorithm with the VNS framework. This is because the VNS framework significantly improves the solutions (Hansen Mladenović & Moreno Pérez 2010). It is worth noting that the combination of SA and VNS for optimizing other problems has been used previously (Jandaghi Divsalar & Emami 2021), demonstrating improved solution quality and solving speed. However, this combination has not yet been applied to water network optimization. Furthermore, a review of the literature revealed that none of the presented heuristic algorithms have addressed the issue of parameter tuning.

As mentioned earlier, several heuristic approaches have been developed in the literature to optimize WDNs. Among these approaches, algorithms that combine variable neighborhood search (VNS) with different metaheuristic algorithms have garnered significant attention due to their promising results. However, VNS has not been applied to WDN optimization. Therefore, this study introduces a new approach that combines SA with variable neighborhood search (VNS), termed SA-VNS, incorporating an efficient local search for WDN optimization. Additionally, parameter tuning is employed to guide the solution process toward the optimal solution and reduce convergence time. The performance of this algorithm was compared with similar algorithms across different networks. Furthermore, the model's ability to address a real case study was investigated. The coupled hydraulic modeling and optimization procedure was used to assess the design for network extension. The results indicate that, as expected, improvements were observed in the cost function, convergence rate, and overall algorithm performance.

The algorithm comprises a local search and a shaking phase, as indicated in Figure 2, both of which are detailed in the next sections. The SA-VNS algorithm is coded in C + +, and the EPANET C ++ toolkit is linked for evaluating hydraulic and steady-state constraints.

The initial solution was generated using a random structure. First, random diameters for the pipes were selected, and then the hydraulic constraints were evaluated. If the hydraulic constraints were satisfied, these random diameters were used as the starting solution for the SA-VNS algorithm.

In the shaking phase, a percentage (β) of the total number of pipes was selected, and a new random diameter was assigned to each pipe, ensuring that the hydraulic constraints were satisfied. The shaking percentage is a parameter of the algorithm, with a value between 0 and 1, and will be discussed in the next section.

Several experiments were conducted on the proposed algorithms. First, this section details the test instances and their creation process. Benchmarks from the literature were also performed. Next, the Taguchi design of experiments (DOE) method was used to tune the parameters of the proposed algorithm. The experimental results were then analysed to compare the performance of the SA-VNS approach on the benchmark instances. All experiments were conducted on a PC with a 2.4 GHz Core 2 Duo CPU, 6 GB RAM, and Windows 10. The SA-VNS algorithm was implemented in Visual C ++ using Visual Studio 2010, and the steady-state constraints were evaluated using the EPANET toolkit for C + +. Additionally, statistical analysis was performed with Minitab 18 for the Taguchi DOE.

Several parameters were employed in the architecture of the proposed SA-VNS algorithm that must be initialized before running the main experiments, as discussed in Section 3. Fine-tuning these parameters may have a significant impact on the algorithm's performance in terms of both the objective value and the runtime (number of evaluated functions). The parameters of the SA-VNS algorithm include the initial temperature (T_max), the temperature decreasing coefficient (α), the shaking percentage (β), the maximum number of iterations, and the percentage of sorted pipes where the diameter decreases by 3 units (γ). In this section, the Taguchi design of the experiment method is used to design the necessary experiments to select the optimal value for each of these parameters among various predetermined levels. The Taguchi technique works by first performing some tests based on these specified parameter values, then using the outcomes of all experiments as the parameters. After that, this method employs a loss function to find the optimal combination of parameters. The approach applied to analyse the tests in this research is the signal-to-noise (S/N) ratio method. This approach expresses the amount of dispersion around a particular value, which means how results vary across different experiments. In Tables 1 and 2, the tuning parameters and the levels considered in the tests for small networks (fewer than 30 nodes) and large networks (more than 30 nodes) are presented. For the SA_VNS algorithm, five values for each parameter were considered.

In large networks, it was observed that SA-VNS was sensitive to shaking, with β set to 0.2 being appropriate for such networks. On the other hand, a higher number of iterations suggests that more local search can lead to better solutions. For both small and large networks, tuning the γ parameter to 0.6 indicates that lowering the diameter of all pipes does not always result in a better solution. Therefore, this parameter, which is innovatively designed in the SA-VNS algorithm, is important.

The performance of an optimization algorithm depends on three categories: effectiveness (i.e. the measure of how close the algorithm gets to the global optimum), efficiency (i.e. the computational effort required to obtain the solution), and reliability (i.e. the ability to achieve the same optimal solution by adjusting the control parameter settings). To compare the performance of different metaheuristic algorithms, many runs and some statistical metrics of the results are required.

The proposed definition of η_algorithm in Equation (7) is based on the individual influences of N_gen and N_obj−eval. The equation satisfies the upper extreme ideal condition which gives the ideal performance indicator (η_algorithm = 100%). When N_gen = N_obj−eval = 1, the ideal number of generations matches the ideal minimal number of evaluations. Beside the satisfied condition, the mathematical formulation of Equation (7) takes consideration for large networks, the value of N_gen is very large and the common logarithm with base 10 simplifies the correlation with each of N_gen and N_obj−eval and decreases the floating-point errors. The values 0.99 and 0.01 (i.e. their sum equals 1) used in the second and third terms of Equation (7) balance the importance of N_gen over N_obj−eval by making an algorithm with low N_gen and high N_obj−eval yields higher values than an algorithm with high N_gen and low N_obj−eval. Here all studied benchmarks have been performed 15 times and based on Equation (7), the standard deviation and performance indicator are calculated as Tables 16–18.

The performance of the SA-VNS algorithm for the two-loop network shows that both the average and maximum costs have improved compared to previous studies. Additionally, the SA-VNS has the lowest standard deviation among those reported in the literature. Furthermore, the performance indicator of the SA-VNS algorithm was superior to 11 previous studies, with only one study showing slightly different results. These findings demonstrate that the SA-VNS algorithm delivers high-quality results in terms of both the difference between optimal solutions and efficiency for the studied network.

The performance results of the SA-VNS algorithm for the New York City network indicate that the average cost was lower than that of all other algorithms, and the standard deviation was also better than in previous studies. Additionally, the performance indicator of the SA-VNS algorithm surpassed that of other algorithms.

Finally, the performance analysis of the SA-VNS algorithm for the El-Mostakbal City network, tested by three researchers (as mentioned in Section 4.3), shows that the SA-VNS achieved better results in terms of average and maximum costs, standard deviation, and performance indicator, outperforming all other methods. It is worth noting that the improved results compared to previous studies are due to two main factors: first, parameter tuning was not performed in previous investigations, and second, this novel algorithm incorporates a new and efficient local search mechanism.

In this paper, the SA-VNS algorithm is applied to a real WDN. This section first describes the city properties and network specifications, then explains the demand estimation method, and finally presents the network optimization and the obtained results.

Tank and reservoir characteristics are detailed in Tables 19 and 20. The reservoirs are equipped with four types of pumps, designated as A, B, C, and D in the tables. The number of pumps in each reservoir is also listed in Table 19. The water network consists of two types of pipelines:

• Distribution Pipelines: 5,285 pieces with a total length of 257 km.

• Transmission Pipelines: 1,007 pieces with a total length of 52 km, as shown in Figure 18.

In some areas, where demand is high, multiple pipes are used in parallel between two nodes. To simplify the model, parallel pipes between two nodes are treated as a single pipe. It is important to note that head losses due to elbows, tees, valves, expanders, reducers, etc., are not considered in this analysis. The pipelines are made from various materials: polyethylene, ductile iron, asbestos cement, and galvanized steel. The Hazen–Williams coefficients used for these materials are as follows:

• Polyethylene: 140, Asbestos cement: 120, Ductile iron: 100, Galvanized steel: 100.

To estimate the demand at each node, the Thiessen polygon method, introduced by Thiessen (1911) was implemented. This approach can be summarized as:

• Water consumption: based on statistical data, the average water consumption per person in Malard City is 220 L/day.

• Population density: Malard City is categorized into three density levels, each with different consumption rates: low density: 60,000 people/km²; middle density: 30,000 people/km²; high density: 15,000 people/km².

• Drawing the polygons: the Thiessen polygons are constructed using the perpendicular bisectors of triangles formed by connecting network nodes. In total, 74 polygons related to 74 nodes are obtained. It is important to note that only relevant areas with water consumption were included, and irrelevant areas were excluded.

• Calculating node demand: the water consumption at each node is determined using the following formula: