To optimize the design of a water distribution network (WDN), a large number of possible solutions need to be examined; hence computation efficiency is an important issue. To accelerate the computation, one can use more powerful computers, parallel computing systems with adapted hydraulic solvers, hybrid algorithms, more efficient hydraulic methods or any combination of these techniques. This paper explores the possibility to speed up optimization using variations of the *ΔQ* method to solve the network hydraulics. First, the *ΔQ* method was used inside the evaluation function where each tested alternative was hydraulically solved and ranked. Then, the convergence criterion was relived in order to reduce the computation time. Although the accuracy of the obtained hydraulic results was reduced, these were feasible and interesting solutions. Another modification was tested, where the *ΔQ* method was used just once to solve the hydraulics of the initial network, and the unknown flow corrections were added to the list of other unknown variables subject to optimization. Two case networks were used for testing and were compared to the results obtained using EPANET2. The obtained results have shown that the use of the *ΔQ* method in hydraulic computations can significantly accelerate the optimization of WDN.

## NOMENCLATURE

*Q*_{ij}^{(0)}initial flow distribution (m

^{3}/s)- Δ
*Q* flow correction (m

^{3}/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

^{−2})*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)

## INTRODUCTION

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.

## METHODS

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

*ΔQ*

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.

*Q*is introduced. The Hazen–Williams (HW) equation is used to calculate pressure head loss in a pipe. The expanded form of the sum of the head losses for the loop (Figure 1(a)), written in the adopted clockwise direction for the flow correction Δ

*Q,*becomes:where

*R*is the HW pipe flow resistance characteristic, is the initial pipe flow (m

_{ij}^{3}s

^{−1}) and

*n*= 1.85 is the flow exponent (–). The pipe flow resistance characteristic is calculated as:

*L*is the pipe length (m),

_{ij}*D*is the pipe inside diameter (m) and

_{ij}*C*is the HW roughness coefficient (–). Equation (1) is a nonlinear equation that needs to be solved for the flow correction Δ

*Q*. For this purpose, the Newton (also known as Newton–Raphson) iterative method (Hoffman 2001) was employed, which converges quadratically if the initial approximation is sufficiently close to the solution. The general form of the solution is:

*i*is the iteration number,

*ij*is the

*ij*-th pipe in the loop and Δ

*Q*is the

_{p}*p*-th flow correction that corrects the initial flow in the pipe

*ij*(there can be more than one, if the pipe is common for two loops (Figure 1(c)). The

*sign*equals one (1) if the direction of the introduced correction Δ

*Q*is the same as the direction of the initial flow and minus one (−1) if otherwise. The iterative calculation is done until the desired precision is reached.

_{p}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:

### Minimal basis loops' detection

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.

*InitSet*with such configuration that there are two pipes (p1 and p4) in which flow needs to be corrected with both flow corrections Δ

*Q*and Δ

_{1}*Q*, since both of these pipes are common to the two loops. The total number of pipes in the loops for this configuration is 4 + 3 = 7. In order to get a better configuration, unused pipes will be swapped with adjacent ones that are used during the propagation. For loop1, pipe p3 can be swapped either with pipe p2 or p4. In both cases the ST configuration would be such that the number of links in the loops would stay the same, 4 + 3 = 7. Now the second loop, loop2, is tried. Its swapping pipe is p5 and it can be swapped with p1 or p4. If pipe p5 is swapped with p4, a new configuration of the ST gives the total number of links in loops as 3 + 3 = 6 (Figure 2(b)). The same configuration would be obtained if p5 was swapped with p1. So one of these configurations is chosen (no matter which one) as it is better than the previous one and it is marked as the new best configuration. The procedure is repeated until a better configuration cannot be reached. The last configuration is remembered as the

_{2}*STSet*. For the considered simple example, this configuration is at the same time the one with the minimal basis loops. This means that in the last step of the algorithm, which is the extraction of the final set of minimal basis loops,

*STSet*will be identified as

*FSet.*This, however, does not have to be the case, so in the following paragraph the last step of the algorithm (extracting the

*FSet*) will be explained in detail and then illustrated on a more complex example.

*FSet*) is based on the sorting of the

*STSet*set of loops by their length (number of pipes) and creation of the overlapping matrix which shows the number of common (shared) pipes between any two loops. The overlapping matrix

*A*is a square matrix with the dimensions

*NL × NL*, where

*NL*is the number of the loops. The diagonal elements of the matrix A are

*Aii =*0 since it represents an overlapping of loop with itself and

*Aij*is equal to the number of shared pipes between the loops

*i*and

*j*. The decomposition considers two loops at the time and it starts from the loops with the highest number of shared pipes. The loop with more pipes, or with more shared pipes, is considered to be the large one (

*Lloop*) and the other one is the small one (

*Sloop*). Once the candidates for the decomposition have been detected, they are combined to create a new loop (

*Nloop*) which could possibly replace one of them in the

*STSet*of loops if it is a ‘more minimal basis’. First, it is checked whether the

*Sloop*is a minimal one. It is defined as minimal if it shares one pipe with the others, or if it does not share any pipes, in which case it is called a ‘hanging loop’ (

*Hloop*). If it satisfies this condition it is moved to the final set of loops (

*FSet*). If

*Nloop*has a smaller number of pipes than

*Lloop*it replaces it in the

*STSet*, which now has one loop less since the

*Sloop*is transferred to the

*FSet*. If it has more pipes than the

*Lloop*, it still can be a candidate to stay in the

*STSet*if it shares smaller (or the same) number of links with the other loops from the

*FSet*and

*STSet*. This comes from the definition of a minimal basis loop in which another criterion, besides the minimum number of pipes in a loop, is the minimum number of shared pipes. After this, the current

*STSet*is searched for the loops (

*Hloops*). If found, they are also moved from the

*STSet*to the

*FSet*. The described steps of the algorithm are repeated until the number of loops in the

*FSet*corresponds to the number of loops

*NL*. A pseudo code for the algorithm is presented in Figure 3.

*InitSet*which is shown in Figure 4(a). This set consists of four loops with a total of 16 pipes and nine of them are shared between the loops. The second step of the algorithm, described earlier, will provide an improved set of the loops –

*STSet*(Figure 4(b)). This set has a smaller number of pipes (14) and smaller number of pipes that are shared (4). Now the previously described third step is employed to extract the final set of loops. The first matrix A is created. Loop1 (p3, p8, p9, p10, p4) is chosen for the

*Lloop*since it is the longest one. It shares pipes with all other loops but with loop3 (p3, p7, p4) it shares most of them, two to be precise (p4 and p3). Thus, loop3 is chosen as the

*Sloop*. Combining these two loops the

*Nloop*is obtained (p8, p9, p10, p7).

*Sloop*is minimal as it shares only one pipe with other loops, p3 with loop4 and p4 with loop2, and is transferred from the

*STSet*to the

*FSet*.

*Nloop*has four pipes which is smaller than the five pipes of the

*Lloop*, so

*Nloop*replaces

*Lloop*in the

*STSet*. Now

*STSet*is modified and it has three loops that share no links (Figure 4(c)). All of them are

*HLoops*and are transferred to the

*FSet*, which now has four loops and the algorithm is finished.

*FSet*has the loops with a total of 13 pipes and three (p3, p4, p7) of them are shared between the loops.

### Optimization method

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

### Implemented optimization algorithms

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 upgraded ΔQ method inside an optimization algorithm (Method A)

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):

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

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

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:

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.

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

The fitness function is calculated.

#### The fixed ΔQ method inside an optimization algorithm (Method B)

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:

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.Using the fixed flow corrections, the pressure head distribution for the fictitious branch network is calculated in only one pass.

The fitness function is calculated.

#### The variable ΔQ method inside an optimization algorithm (Method C)

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:

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.

The fitness function is calculated.

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 ^{0}*,

*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.

^{0}*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*.

#### The fixed iteration ΔQ method inside an optimization algorithm (Method D)

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

The evaluation function calculation:

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.

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

The fitness function is calculated.

#### The reference method: EPANET2 DLL inside an optimization algorithm (Method R)

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:

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.The fitness function is calculated.

*f*; 2) computation time

*t*; and 3) speedup factor, expressed as the ratio of the reference optimization algorithm (

*R*) computation time and examined algorithm (

*A, B, C*or

*D*) computation time.In the case of algorithm

*A*, the convergence criterion was 10

^{−6}m

^{3}s

^{−1}in two successive iterations. For

*Method D*, the number of iterations was fixed to 10 if the convergence criterion is not met.

### Case networks

#### NYT network

*ΔQ*method application was extracted from the literature (Dandy

*et al.*1996). The starting network for optimization is presented in Figure 5(a).

This is an example of a gravitational WDN made out of *n _{n}* = 20 nodes with

*n*= 1 source node or reservoir. The nodes and the reservoir are interconnected with

_{r}*n*= 21 large pipes forming

_{p}*n*= 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

_{l}^{25}. 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*for failing to meet the minimal nodal pressure head values (

_{p}*p*= 20 m H

^{min}_{2}O).

In the above equation, *C _{k}*,

*D*and

_{k}*L*are cost of the new pipe per meter, diameter and length of the pipe, respectively,

_{k}*j*is the number of a node,

*p*is the pressure head value in the node

_{j}*j*while

*O*is the function whose value is above zero if the

*p*. The specific value of the penalty function is

_{j}< p^{min}*C*= 15,000,000 $/m.

_{p}#### FOS network

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).

*n*= 36 nodes with

_{n}*n*= 1 source node or reservoir, interconnected with

_{r}*n*= 58 pipes forming

_{p}*n*= 22 loops. The minimum pressure head of all the demand nodes is 40 m, while the maximum pressure head of each node is given in Table 1. Besides pressure head condition, flow velocity condition needs to be satisfied, where in each pipe it has to be less than 1 m/s. For the network pipes, 22 different diameters are available in catalogue. This means that the search space has the size of 7.25 × 10

_{l}^{77}. The fitness function is made of three parts: investment in the distribution network

*I*, penalty

*I*for unfulfilled pressure head conditions and

_{p}*I*penalty for failing to meet maximum velocity condition.

_{v}N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} |
---|---|---|---|---|---|---|---|---|---|---|---|

1 | 55.85 | 7 | 53.10 | 13 | 59.10 | 19 | 58.10 | 25 | 56.6 | 31 | 56.60 |

2 | 56.60 | 8 | 54.50 | 14 | 58.40 | 20 | 58.17 | 26 | 57.6 | 32 | 56.80 |

3 | 57.65 | 9 | 55.00 | 15 | 57.50 | 21 | 58.20 | 27 | 57.1 | 33 | 56.40 |

4 | 58.50 | 10 | 56.83 | 16 | 56.70 | 22 | 57.10 | 28 | 55.35 | 34 | 56.30 |

5 | 59.76 | 11 | 57.30 | 17 | 55.50 | 23 | 56.80 | 29 | 56.5 | 35 | 55.57 |

6 | 55.60 | 12 | 58.36 | 18 | 56.90 | 24 | 53.50 | 30 | 56.9 | 36 | 55.10 |

N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} | N
. _{i} | P (m)
. _{max} |
---|---|---|---|---|---|---|---|---|---|---|---|

1 | 55.85 | 7 | 53.10 | 13 | 59.10 | 19 | 58.10 | 25 | 56.6 | 31 | 56.60 |

2 | 56.60 | 8 | 54.50 | 14 | 58.40 | 20 | 58.17 | 26 | 57.6 | 32 | 56.80 |

3 | 57.65 | 9 | 55.00 | 15 | 57.50 | 21 | 58.20 | 27 | 57.1 | 33 | 56.40 |

4 | 58.50 | 10 | 56.83 | 16 | 56.70 | 22 | 57.10 | 28 | 55.35 | 34 | 56.30 |

5 | 59.76 | 11 | 57.30 | 17 | 55.50 | 23 | 56.80 | 29 | 56.5 | 35 | 55.57 |

6 | 55.60 | 12 | 58.36 | 18 | 56.90 | 24 | 53.50 | 30 | 56.9 | 36 | 55.10 |

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*is the flow velocity in

_{k}*k-*th pipe,

*v*is the maximal velocity for

_{k}^{max}*k-*th pipe. The specific values of penalty functions for pressure and velocity are

*C*= 15,000,000 €/m and

_{p}*C*= 50,000,000 €/m, respectively.

_{v}## RESULTS AND DISCUSSION

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 _{1}* 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.

_{2}GA optimization algorithm . | Fitness function f [10^{6} $]
. | CPU time t [s]
. | Speedup factor [-] . |
---|---|---|---|

EPANET2 DLL (R) | 38.6 | 390 | / |

Upgraded ΔQ (A) | 38.6 | 18.5 | 21.1 |

Fixed ΔQ (B) | 40.2 | 5 | 78.0 |

Variable ΔQ (C) | 39.8 | 5.5 | 70.9 |

Fixed iteration ΔQ (D) | 39.0 | 18.4 | 21.2 |

GA optimization algorithm . | Fitness function f [10^{6} $]
. | CPU time t [s]
. | Speedup factor [-] . |
---|---|---|---|

EPANET2 DLL (R) | 38.6 | 390 | / |

Upgraded ΔQ (A) | 38.6 | 18.5 | 21.1 |

Fixed ΔQ (B) | 40.2 | 5 | 78.0 |

Variable ΔQ (C) | 39.8 | 5.5 | 70.9 |

Fixed iteration ΔQ (D) | 39.0 | 18.4 | 21.2 |

*InitSet*had a total of 175 pipes, the

*STSet*143 and the

*FSet*101 pipes. Optimization algorithms

*A*and

*D*were tested for four different GA settings (Table 3). Comparisons of computed performance indicators are shown in Figures 6 and 7, respectively.

No. test . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|

GA generations | 1,000 | 1,000 | 2,000 | 2,000 |

GA population | 100 | 200 | 100 | 200 |

No. test . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|

GA generations | 1,000 | 1,000 | 2,000 | 2,000 |

GA population | 100 | 200 | 100 | 200 |

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

^{6}$ (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 _{1}* and

*ΔQ*deviate not more than 25% from the first iteration (Ivetić

_{2}*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^{−6} m^{3} s^{−1}, 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.

## CONCLUSIONS

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.

## ACKNOWLEDGEMENTS

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’.