Speeding up the water distribution network design optimization using the ΔQ method

Q_ij⁽⁰⁾

initial flow distribution (m³/s)

ΔQ

flow correction (m³/s)

R_ij

pipe flow resistance

n

flow exponent (-)

L_ij

pipe length (m)

D_ij

pipe inside diameter (m)

Λ

friction factor (-)

g

gravity acceleration (≅9.81 m s⁻²)

C

Hazen–Williams roughness coefficient (-)

ΔH_ij

pipe headloss (m)

ΔH_res

head difference between reservoirs (m)

v_k

flow velocity (m/s)

p_j

pressure head (m)

f

fitness function ($),(€)

I

investment in the network ($),(€)

I_p

pressure head condition penalty ($),(€)

I_v

flow velocity condition penaly ($),(€)

C_p

specific value of pressure head penalty function ($),(€)

C_v

specific value of velocity penalty function (€)

t

computation time (s)

In the last decades, modern engineering practice has been moving towards a higher degree of computer and software usage. As computer sciences experience fast progress, engineering practice is trying to catch up and utilize stronger processors and new software capabilities that are introduced. One of the engineering fields which benefits from this progress is hydraulics and water system analysis, especially the numerical modelling of fluid flow and optimization of the design and performance of water systems. In this paper the focus is on the optimization of water distribution networks (WDN), specifically the long computation time needed for multiple runs of the hydraulic solver and the possibility to reduce it.

In order to accelerate the computation or optimization process, different approaches are possible: the use of more powerful computers (faster processors), parallel computing systems with adapted hydraulic solvers, hybrid algorithms, more efficient hydraulic methods or any combination of the previously mentioned options. The first approach is strictly hardware-driven, in which the focus is on the IT industry and its capability to develop hardware with more raw power. The second approach has been addressed by numerous authors, computer and communication scientists, showing a variety of success in the results (Di et al. 2009; Artina et al. 2012; Marques et al. 2012; Smith et al. 2013). The SWMM numerical model for hydrodynamic rainfall–runoff and urban drainage simulations has been parallelized with minimal changes to the original code by Burger et al. (2013), where the achieved speedup was from six up to 10 times. On the other hand, several attempts to parallelize hydraulic solvers inside the EPANET2 have not been that successful (Alonso et al. 2000; Crous et al. 2012; Burger 2014) although it seemed a straightforward task. The idea was to tackle the matrix operations performed in each iteration step of the solver with several available packages that support parallel linear algebra (von zur Gathen 1993; Van de Gejin 1997; Agullo et al. 2009). In the work of Diao et al. (2013), significant code changes were made in order to implement a domain decomposition as an idea for EPANET2's acceleration. The resulting speedup was encouraging, reaching the value of eight times on one tested network but the issue of severe and complicated code changes makes this approach hard to implement. Zecchin et al. (2012) showed that iterative solvers preconditioned with the algebraic multigrid method (AMG) are faster than the current EPANET2 solvers only in cases of large artificial networks, while for typical EPANET2 problems speedup was not obtained. Different optimization algorithms have been used so far for the WDN design optimization, such as differential evolutionary algorithms (Suribabu 2010) and modified genetic algorithms (Montesinos et al. 1999; Neelakantan & Suribabu 2005). Hybrid algorithms, in general, as the tools for multi-objective design, have outperformed the non sorting genetic algorithm (NSGA II), mostly in terms of diversity of obtained Pareto fronts (Wang et al. 2015a), while lowest computational burden has been achieved with the low-level hybrid algorithm proposed by Craeco & Franchini (2012). Since the focus of the paper is not on the optimization method, the standard genetic algorithm (GA) (Savic & Walters 1997), efficient in handling a single objective optimization problem (Holland 1992), is used in combination with more efficient hydraulic methods.

Currently, when a water supply network is to be optimized, usually EPANET2 is used as the solver (Rossman 2007). EPANET2 is a reliable, free software package, available as a standalone EXE or DLL version, which can be easily integrated with most optimization programs. When used for single runs, EPANET2 is quite efficient in terms of the computation time. Large networks can be solved in a few minutes at maximum. However, a research scenario such as the optimization of a network design, involves multiple runs of a hydraulic model, requiring significant computation time for a single run of the optimization algorithm (Diao & Rauch 2013).

To solve the continuity equation EPANET2 uses a so-called global gradient algorithm (GGA) originated by Todini & Pilati (1987). The hydraulic basis of this method, and all other node methods, originates from Hardy Cross and the method he originally called the method of balancing flows (Cross 1936). At that time, Hardy Cross also introduced the ΔQ method or the method of balancing heads, for hydraulic calculation of a network (Cross 1936). The ΔQ method was derived from Cross's original moment distribution method, which he used for structural analysis. The method was in common use in the 1950s and 1960s until the arrival of computers and the node methods took over the throne. Both methods introduced by Hardy Cross originate from ‘paper and pen’ solvers time, and thus are classified as the local approaches for a network solution (Todini & Rossman 2013). The most significant difference is that in the node-based method, the number of non-linear equations is equal to the number of nodes while in the ΔQ method it is equal to the number of loops in the network. The number of loops is usually much smaller than the number of nodes. This already implies a modest reduction of computation time if it is presumed that a similar amount of resources is needed for both types of equations. For example, the large BWSN2 benchmark network (Ostfeld et al. 2008) with 12,527 nodes has only 2,308 loops. The authors' experience is that even aggregated, real size water distribution networks usually have about five time less loops than nodes.

Craeco & Franchini (2014) presented a comparison of the Newton–Raphson global algorithm and the ΔQ method (loop flows algorithm), in which the latter had shown a slightly smaller computation burden on fictitious networks. In the same paper, the ΔQ method is recommended for network design optimization algorithms.

In this paper, the possibility of using the ΔQ method in hydraulic calculations of the network inside the optimization algorithm's evaluation function was explored. An iterative solver of the non-linear equations based on the ΔQ method was programmed in C ++ as a DLL and called from the optimization algorithm to make it comparable with the EPANET solver. The base point for the analysis was the upgraded ΔQ method in which the hydraulics was solved using the ΔQ method in which the shared flow corrections are updated in the equations as the computation progresses. Furthermore, the following three modified variations were derived: the fixed ΔQ method, the variable ΔQ method and the fixed iteration ΔQ method. With the first two variations, an attempt was made to bypass the computational burden of the iterative solver. The idea was to promote intelligent, hydraulically based solutions. The first results presented by the authors (Stanić et al. 2012; Ivetić et al. 2014), on the well-known problem of the New York tunnels (NYT), were promising both in computational time and in reaching the suboptimal solution. A step forward in terms of the network complexity was made, and the intermediate network of Fossolo (FOS) (Bragalli et al. 2012) was tested. All of the obtained results are presented in this paper.

The paper is structured as follows: in the ‘Methods’ section, the basics of the ΔQ method and minimal loop detection algorithm are recapitulated, followed by a brief description of the used optimization algorithm and variations of ΔQ method implementation. The section is concluded with a brief formulation of the two case networks and the performance indicators used for the comparison with the EPANET2-based algorithm. The obtained improvements of the performance indicators are presented and deliberated in the ‘Results and discussion’ section. Finally, conclusions are drawn and directions for further analysis are formed.

In order to employ the ΔQ method for hydraulic calculations, an analysis of network graph topology has to be performed in the pre-processing stage. It is necessary to detect all cycles, or loops in the network which are later split and thus the network structure is changed (a tree-structured graph is obtained). Tree-structured or branched networks are quite easy to handle in a water distribution analysis and unknown flows and nodal heads can be obtained in a double sweep algorithm (Stanić et al. 1998).

The ΔQ method was originally presented by Hardy Cross in 1936. It was proposed as one of the two possible methods for flow analysis in networks of conduits. In his original work, the ΔQ method was called the ‘Method of balancing heads’ (Cross 1936). The method is based on the fact that in every closed loop (circuit) of the water supply network the sum of total head loss is equal to zero. This is derived from the conservation of energy equation for a closed loop.

The number of equations that need to be solved corresponds to the number of loops in the network. In the work of Todini & Rossman (2013), this method is called the loop flows algorithm and it solves equations for all the loops in the network simultaneously. In this research, equations were solved in succession as the speed of the convergence can be improved if the flow correction for one loop, calculated in the current iteration, is used to calculate the flow correction for another loop with which it shares a common pipe. This way the ΔQ solver is upgraded. Take Figure 1(c), for example: If the flow correction for the first loop is calculated in the i-th iteration ⁠, when the calculation for the flow correction in the second loop is conducted, flow in their common pipe 2–3 can be expressed as ⁠.

Aside from the ‘ordinary loops’, the term ‘pseudo loop’ is introduced (Larock et al. 2000). This is a loop that is formed between two reservoirs/tanks with defined heads. A number of such loops is one less than the number of reservoirs in the network. In that case Equation (3) becomes:

Prior to conducting the calculations using the ΔQ method, network loops need to be detected. Initial loops' detection is done based on the results of the BFS propagation algorithm previously mentioned. The number of loops corresponds to the number of unused links during the BFS propagation, which are not a part of the spanning tree so they must be closing the circuit and creating the loops. The initial loops are not likely to be geometrically minimal. The geometrical minimal loop is defined as the one that cannot be presented as a union of any other loops. Detecting these loops is not an easy task. Some algorithms are based on using the outer or the ‘back edges’ of the network (Jha 2007) but these have to be predefined. It is clear that in the case of thousands of pipes this would be a very demanding job. Craeco & Franchini (2014) presented an algorithm which utilizes the Dijkstra algorithm to search for the shortest path (from the topological viewpoint, meaning that all graph links have the same weight of one) between two nodes.

In this research, another approach based on the graph theory algorithms is presented. A network loop is considered to be minimal if it shares common pipes with a minimum number of other loops. Also, the number of these shared pipes should be minimal. The authors refer to these loops as the topologically minimal loops and in the following sections they are referred to as the ‘minimal basis loops’. This algorithm does not guarantee that identified minimal basis loops will be the absolutely topological minimal or geometrically minimal, which after all is not necessary for the ΔQ method to be employed. However, minimal basis loops defined in this manner will enable obtainment of the simplest form of nonlinear equations to be solved. The algorithm has been tested on numerous examples and benchmark problems used in urban water modelling and optimization literature and, in most cases, it found minimal basis loops that are both topologically and geometrically minimal. The presented algorithm for the minimal basis loops' identification takes three steps. First, the BFS propagation algorithm is run to obtain the initial spanning tree (ST) and the initial set of loops (InitSet). The second step is the transformation of the ST in order to get the set of loops with the smallest number of pipes (STSet). Finally, decomposition is done if needed and a set of final minimal basis loops (FSet) is extracted from the graph. The algorithm steps are explained in the following section.

GAs are employed as an optimization method and are efficient in terms of running time and finding suboptimal solutions (Holland 1992). GAs are called inside an optimization algorithm, after the pre-processing stage, with the task to generate coded solutions which are to be tested inside an evaluation function. The EPALG software, developed at the Faculty of Civil Engineering at the University of Belgrade, was chosen for the GA. For every solution examined, a new value of the fitness function is computed. Based on this value, examined solutions are ranked, where the solution with the lowest value of the fitness function is ranked as suboptimal.

In this paper, some of the EPALG settings were kept fixed throughout all of the optimization algorithm runs to have comparable results between different methods. These settings regard the GAs' ability to converge to the best suboptimal solution and affect the way mutation, crossover, selection and replacement are done inside GAs themselves. The settings used were:

• Mutation: Type: Reinit number of bits affected: 1

• Crossover: Type: Two point Probability: 0.8

• Selection: Type: Tournament Number of solutions competing: 2

• Replacement: Type: Uniform Portion of population being replaced: 1

In the pre-processing stage after the FBN is derived, unknown flows were calculated for the new network using a back sweep algorithm. These flows satisfy the continuity equation in the network nodes, but are not an exact solution for the pipe flows, which means that they can be used as a starting assumption of the flows for the ΔQ method. This is done only once, so later calls of the evaluation function will just include an iterative solution for the flow corrections, in the manner explained in ‘The ΔQ method’ section. When the convergence condition is met, final flows can be calculated and pressure head distribution is obtained using a forward sweep algorithm. This idea was named the upgraded ΔQ method (Method A), and it became the starting point of the investigation.

Several variations of optimization algorithm were derived and examined in order to reduce computation time. Tests were made with the version of the optimization algorithm in which values of the flow corrections were kept fixed, equal to the initial values obtained in the pre-processing stage – the fixed ΔQ method (Method B). There was a major applicability issue with this method in the more complex and new networks. In complex networks, with a large number of loops with diameters inside a single loop varying significantly, the starting assumptions for flow corrections, if they can be computed, can have a significant difference from their final, exact values. For the new networks for which there are no predefined pipe diameters, starting assumptions for flow corrections cannot be calculated in the pre-processing stage.

In order to overcome these shortcomings, another approach was tested, the variable ΔQ method (Method C). The initial values of flow corrections, as calculated during the pre-processing stage, were assumed to be new unknown variables, subject to optimization together with other unknowns (e.g., pipe diameters). In addition, one more penalty function has to be added in the fitness function calculation, which will guide the optimization of the flow corrections. Mutual for these two approaches (Methods B and C) is the fact that in each evaluation function call, no time-consuming iterative hydraulic computation has to be done. Furthermore, based on the results of previous tests (Ivetić et al. 2014) with fixed ΔQ method (B), the idea to run only a few iterations inside ΔQ method solver in order to obtain some ‘near to exact’ flow corrections (fixed iteration ΔQ method (Method D) was explored.

In order to make a comparison with the EPANET2 solver, an unbiased, iterative solver for the ΔQ method was made in the form of a DLL file which will be integrated inside the evaluation function of the optimization algorithm.

The key points for the hydraulic solver implementation inside the optimization algorithm are the pre-processing stage and the evaluation function. The ways that these two parts of the algorithm are formed will be the focal points for all of the presented methods' description. Here the ΔQ method iterative solver for Equation (1) is programmed as a DLL file. Equation (1) is solved using the flow corrections obtained through an iterative computation of Equation (3). In every pass through the evaluation function, a solver is run externally and the results of the hydraulic calculation (pressure and flow distribution) are used to compute a value of the fitness function. The pre-processing stage, which is mutual for all alternatives, and the evaluation function computation are described below.

• Pre-processing stage (preceded by a call of GA):

1. The BFS algorithm is run, a directed ST is obtained, the FBN is formed and the initial flow distribution is determined.

2. Minimal basis loop detection as explained in the previous section (Figure 4).

3. If the starting disposition includes prebuilt pipes (pipes with known diameter e.g., NYT), a network flow distribution is computed and the resulting flow corrections ΔQ can be stored as starting values for further GA calls. It is presumed that flow corrections calculated this way will not differ much from their final values which will be calculated for every call of the evaluation function, thus, this should speed up the convergence of Equation (3), if it is possible.

• Evaluation function calculation:

1. For every tested alternative, the ‘exact’ values of the flow corrections are computed by calling the iterative solver (DLL file) from the optimization algorithm. In order to improve the convergence and thus to reduce the computational burden, flow corrections from current iteration are used as explained in ‘The ΔQ method’ section.

2. Using exact flow corrections, the pressure head distribution for the network is calculated in only one pass.

3. The fitness function is calculated.

The fixed ΔQ method is developed and implemented in a similar manner to Method A. It was differentiated from Method A by omitting the iterative calculation of Equation (1), which means that there is no need to call the DLL solver inside the evaluation function.

• The evaluation function calculation:

1. For every tested alternative, the previously determined values of flow corrections ΔQ are used. This variation of ΔQ method implementation is feasible only for networks with pre-existing pipes (pipes with known diameters) due to the fact that the initial values of the flow corrections must be computed in the pre-processing stage.

2. Using the fixed flow corrections, the pressure head distribution for the fictitious branch network is calculated in only one pass.

3. The fitness function is calculated.

In Method A, each evaluation function's calculation involves the iterative calculation of the ‘exact’ flow corrections. In the variable ΔQ method, it is assumed that flow corrections are unknown variables, whose values are, together with pipes' diameters values, optimized.

• The evaluation function calculation:

1. For every tested alternative (having pipes' diameters and flow corrections as variables to optimize) only the pressure head distribution for the FBN is calculated in one pass.

2. The fitness function is calculated.

3. An additional penalty function is added to the value of the fitness function, to compensate for the fact that the used flow corrections are not ‘exact’. If the pressure drop between the neighbouring nodes where the loop is split is positive in the direction of water flow, the value of the penalty function is zero. This solution is feasible, meaning that the pressure reducing valve could be installed in the place where the pressure drop occurs. However, if the pressure drop between neighbouring nodes where the loop is split is negative, the solution would require pumping between these nodes, and it has to be penalized. The same penalty value as for pressure head deficiency is used.

To reduce the search space for flow corrections, the value is entered as multiplication of the flow corrections computed in the pre-processing stage ΔQ⁰, M × ΔQ⁰. The range of possible values for multiplication factors can vary, meaning that the used search space is problem dependent and should be handled with care.

Method C does not solve the network for the ‘exact’ flows, instead it produces some sort of suboptimal network flows which further implies that the accuracy of the hydraulic results is degraded compared to Method A.

This approach uses the DLL iterative solver with a predefined number of iterations.

• The evaluation function calculation:

1. For every tested alternative, near to exact values of flow corrections are computed by a predefined number of iterations in a solver. As in the previous case, convergence is improved through a current iteration flow correction update.

2. Using near to exact flow corrections, the pressure head distribution for the network is calculated in only one pass.

3. The fitness function is calculated.

The reference optimization algorithm, used for comparison purposes, in terms of both suboptimal solutions and needed computation time, is based on the EPANET2 hydraulic solver. The solver is called as a DLL file in every pass through the evaluation function, and obtained results are used for the fitness function value computation. In this case the pre-processing stage is not needed, while the evaluation function processing has the following form.

• The evaluation function calculation:

1. For every tested alternative, an EPANET2 DLL file is called for hydraulic simulation of the network, in the same way as the DLL for the ΔQ method is called in optimization algorithms (A) and (D). The results, needed for the fitness function calculation, are extracted from the solver and stored.

2. The fitness function is calculated.

This is an example of a gravitational WDN made out of n_n = 20 nodes with n_r = 1 source node or reservoir. The nodes and the reservoir are interconnected with n_p = 21 large pipes forming n_l = 2 loops. The current disposition of the system cannot satisfy the minimal nodal pressure head values of 20 m of water column. Therefore, it is necessary to reconstruct the network in order to meet the given condition in nodes. Changing the diameters of existing pipes is not possible due to the problem of excessive water demand shortfall so the remaining options are either: (a) to duplicate the existing tunnel with some of the diameters offered; or (b) to do nothing. The number of diameters in the catalogue for new pipes is 15. Together with the do nothing option, this makes 16 possible solutions for every pipe in the network, forming a search space of size 1.93 × 10²⁵. The fitness function f for this example is made up of two parts; the first one is investment in the water network I (Equation (6)) and the second is the penalty I_p for failing to meet the minimal nodal pressure head values (p^min = 20 m H₂O).

In the above equation, C_k, D_k and L_k are cost of the new pipe per meter, diameter and length of the pipe, respectively, j is the number of a node, p_j is the pressure head value in the node j while O is the function whose value is above zero if the p_j < p^min. The specific value of the penalty function is C_p = 15,000,000 $/m.

With the NYT reconstruction issue being classified as a medium sized optimization problem (Wang et al. 2015b), a step forward was to test an intermediate size optimization problem. The FOS network, another benchmark example from the literature, was chosen. The starting network topography was defined from the EPANET2's inp file (http://emps.exeter.ac.uk/media/universityofexeter/emps/research/cws/downloads/data/3-epanet/FOS.inp).

The pressure head penalty function has an expanded form compared to Equation (6), with a maximal pressure head condition introduced, where p_j^max is the maximal pressure head for j-th node. In the velocity penalty function v_k is the flow velocity in k-th pipe, v_k^max is the maximal velocity for k-th pipe. The specific values of penalty functions for pressure and velocity are C_p = 15,000,000 €/m and C_v = 50,000,000 €/m, respectively.

The results of all presented ΔQ methods (A, B, C and D) are compared to the results obtained using the EPANET2 hydraulic solver (R). The first tests of the ΔQ method implementation, including all variations (A, B, C and D) were run on the NYT network. These different approaches share the same modified FBN (Figure 5(b)). There are two loops in the network: the larger loop 1 contains a reservoir along with the nodes numbered from 2 to 15, while loop 2 has nodes 11, 9, 16 and 20. The location of the split is arbitrary and inconsequential. Loop 1 was split in the proximity of node 6 where a new node 6′ is introduced as a start node for the downstream pipe. Loop 2 was split close to node 20, with the new node 20′ being generated. Flow corrections for loops 1 and 2, ΔQ₁ and ΔQ₂, respectively, are introduced as demands in the nodes 6 and 20, and in the nodes 6′ and 20′ as negative demands or inflows. A comparison of obtained results, for 1,000 generations and population of 100, is given in Table 2.

The computer used for testing was an Intel i7-2630QM CPU with 6 GB of RAM memory. The obtained optimal solution using the EPANET2 DLL, for the NYT problem, is the same as one taken from the literature f_opt = 38.6 × 10⁶ $ (Dandy et al. 1996), proving it as a valid reference for this investigation.

The computed performance indicators for the NYT network (Table 2) showed significant reduction of the computation time for the optimization algorithm in all of the examined variations (A, B, C and D). In terms of the speedup factor, best results were derived with the variation B, with the value of 78.0. This was expected since no iterative computation was performed inside the evaluation function. Approach C had a slightly smaller speedup factor of 70.9, probably due to the fact that an extra penalty function had to be computed for every examined alternative. It is to be expected to have a further reduction of this factor in real applications, since we have knowledge from the previous tests that the correct values of ΔQ₁ and ΔQ₂ deviate not more than 25% from the first iteration (Ivetić et al. 2014), and we could make the search space rather narrow. This is due to the fact that it is a case of an existing, simple network with just two loops. It is obvious that this approach needs to be adjusted in order to allow the GA to perform a successful search for the suboptimal values of flow corrections, in the rather broad search space. Methods A and D had shown a similar computational burden, with the speedup factors of 21.1 and 21.2, respectively. The similar computational burden is probably due to the fact that for the NYT example, a prefixed number of iterations (10) for variation D was sufficient to compute the ‘exact’ or at least near to the ‘exact’ values of flow corrections. Furthermore, approach A produced the best-known solution, while approach D was the second best with almost the same solution as the previous. Methods B and C produced slightly worse solutions. The suboptimal solution degradation occurs mostly due to the hydraulic results' inaccuracy, which resulted in higher penalty function value. The acceleration in algorithms using ΔQ methods is owed to the fact that less equations need to be solved, as well as the way the network is processed. The pre-processing stage, called just once, performs a large proportion of the necessary analysis, therefore later calls for the evaluation function take much less computation time than in reference algorithm.

In the case of the FOS network, tests with four different GA evaluation function settings using two variations A and D were done. Methods B and C could not be used in the presented manner since none of the pipe diameters was known at the beginning. From Figure 6, it is clear that in all tests, the values of the fitness function were similar for reference algorithm R and the two examined variations of the ΔQ method implementation A and D. For this example, none of the used algorithms reached the best known solution, and it can be concluded that in terms of the suboptimal solution, algorithms A, D and R do not outperform each other. In case of A and D, finding these solutions took much less computation time (Figure 7). The speedup factor for algorithm A ranged from 39.3 to 62.6 and for algorithm D from 57.6 to 105.6. Between approaches A and D, roughly the same value of the fitness function was achieved, with algorithm D taking significantly less time. It is appropriate to point out the benefits of the current iteration flow correction update. For a single test run of evaluation function in approach A for the FOS network, it took 23 iterations to reach the desired accuracy for the pipe flows of 10⁻⁶ m³ s⁻¹, as opposed to 39 iterations for solution without this modification (about 40% reduction in number of iterations). Similar numbers were obtained for different test runs but this percentage clearly depends on the loops; configuration and pipe diameters.

During the optimization of the WDN, hydraulic computation of a network inside an evaluation function consumes most of the computation time. In this paper, use of the ΔQ method for hydraulic computation inside an optimization algorithm is presented through several different approaches. Each of these approaches requires a pre-processing stage, in which the loops in the original water network are detected through minimal basis loop detection algorithm, and network split into the FBN. If a prebuilt network exists, as in the NYT case, initial values of the flow corrections are assumed through hydraulic computation of this branch network.

The upgraded ΔQ method (A) approach computes the correct values of the flow corrections inside each evaluation function, the fixed ΔQ method (B) omits the iterative computation of the flow corrections by using the initial values throughout the optimization run, while the variable ΔQ method (C) includes the values of the flow corrections, through multiplicative factors, as additional variables for optimization. Finally, the fixed iteration ΔQ method (D) is actually the approach A used with a fixed number of iterations of Equation (3). Only methods A and D use the hydraulic solver programmed as a DLL, as well as the reference algorithm (R) which uses EPANET2 in a DLL form. The first tests of the presented variations were done on the NYT problem. In all of the cases, significant computation time reduction was achieved, primarily due to the fact that the ΔQ method solves fewer equations than the GGA integrated in EPANET2. These type of results are expected for real size networks in which the number of loops is rarely higher than 20% of the number of nodes (e.g., BWSN2) (Ostfeld et al. 2008). This implies 80% less equations to be solved with the ΔQ method.

On top of this, in the variations B and C, the network hydraulics is solved only once in the pre-processing stage which makes them even faster. In terms of the quality of the suboptimal solutions obtained, only algorithm A managed to compute a global optimum. The others, B, C and D, have shown a slight degradation in this performance indicator, which is caused by the hydraulic inaccuracy.

Further testing was undertaken on the benchmark example of the FOS network, where approaches A and D were compared with the reference EPANET2 based algorithm R. The tests, with four different GA configurations, have shown a major speedup, with the speedup factor reaching the value of over 100. Since method D had a fixed number of iterations it was the fastest method tested. In terms of finding the suboptimal solutions, both A and D found similar solutions as algorithm R.

The obtained results have shown that the ΔQ method is remarkably faster than EPANET2 when used in medium and intermediate optimization problems. Further investigation will be undertaken on real size water distribution networks. Apart from just changing the hydraulic solver, this approach can be utilized combined with other ways to reduce computation time, such as parallelization, network decomposition, etc.

The authors express their gratitude to the Serbian Ministry of Education and Science for support through the project TR37010: ‘Rain water drainage systems as part of the urban and transport infrastructure’.