The reliability and ease of applying metaheuristic methods in solving large and complex equation systems make it interesting to be applied as an alternative solution to solving problems in various fields. This article proves the effectiveness of an optimization model based on the metaheuristic method for the analysis of hydraulic parameters of drinking water distribution pipes. The metaheuristic methods explored are Differential Evolution (DE) algorithm, Particle Swam Optimization (PSO) algorithm and CODEQ algorithm. The effectiveness of the three methods is measured relative by comparing the results of the analysis of the three models with the results from Newton Raphson method and Monte Carlo simulation method. The analysis shows that the optimization model based on the DE, PSO and CODEQ algorithms is very effective for solving problems on a simple network that has 6 pipe elements and 5 service nodes. The results obtained have a level of accuracy as good as Newton Raphson method. In the case of complex networks that have 32 pipe elements and 21 service nodes, there is an indication of performance degradation which is indicated by a decrease in fitness value. In this case, Newton Raphson method still shows its consistency. The optimization model based on the metaheuristic method is still far more effective than the Monte Carlo simulation method, although it is not as effective as Newton Raphson method. The Monte Carlo simulation method is not recommended for hydraulic pipe network analysis, even for simple networks.

Keywords

de-pso-codeq Algorithm, Hydraulic Analysis, Network Pipe

Share and Cite:

1. Introduction

Analysis of the distribution network hydraulic parameters is an important part of the design and maintenance of drinking water supply systems. In the case of a branched network, the hydraulic analysis can be done more simply because even though the hydraulic equation system is non-linear but it is still possible to be transformed into a linear function. Thus the solution can be solved analytically. But in a closed network system (loop) solving the equation system by an analytical way is not possible. In the case of a complex network, the equation system formed is non-linear and complex so that it requires special methods to solve it.

Until now there have been many methods for hydraulic analysis of drinking water distribution pipelines proposed by world researchers. Method development starts from the graphical method, the use of physical analogies to the use of mathematical models. Development of hydraulic analysis methods for drinking water distribution pipelines in general include; the Hardy Cross method (Cross, 1936; Hoag and Weinberg, 1957), Simultaneous Node method (Martin and Peters, 1963; Shamir and Howard, 1968), Simultaneous Loop method (Epp and Fowler, 1970; Jeppson, 1976), Simultaneous Pipe petode namely the Linear Method (Wood and Charles, 1972), the Simultaneous Network method or the Global Gradient Method (Todini and Pilati, 1987). These methods are developed from the basis of physical analogies and mathematical modeling [1]. Efforts to improve the performance of hydraulic analysis methods were proposed by several researchers, including; nodal analysis models of looped water distribution networks [2], and hydraulic analysis of water supply networks using a Modified Hardy Cross method [3].

There are currently many application packages offered by manufacturers for pipeline hydraulic analysis. A commercial application package called EPANET has been developed not only for hydraulic simulations for pipelines, but includes water quality analysis [4] [5]. EPANET has been successfully applied for pipe hydraulic analysis at “Teiul Doamnei” in Bucharest which has 250 pipe elements and 212 nodes to serve a population of 40,000 [6]. During its development, EPANET became one of the most popular applications because of its reliability and ease of modification to solve certain problems. EPANET can be applied to solve problems in pipelines that have discharge and pressure gaps at each service node [7]. ELGTnet, a commercial application package proposed by Gupta and Prasad can complete a pipeline hydraulic analysis with a very short iteration time [8]. The WaterCAD Simulator was successfully applied to the water distribution network in the Sakwa Region of Nigeria [9]. Loop 4.0 and Permata Air V8i were successfully applied for hydraulic analysis and optimization of pipe diameters in rural networks in zone 1 Nava shihora state of Gujarat, India [10].

Along with the development of the computing world, the metaheuristic method is considered quite reliable because of the ease of its implementation and its ability to find “good” and fast solutions, especially to solve systems of equations that have high dimensions. Metaheuristics is a method for finding solutions that integrate interactions between local search procedures and higher strategies to create processes that can get out of the local optima point and search in the solution space to find global solutions. The shuffled complex evolution (SCE) algorithm can be used as an alternative solution for pipeline hydraulic analysis. Under special conditions, the SCE Algorithm can solve problems that cannot be solved by the Global Gradient Algorithm (GGA) [11]. Application of Differential Evolution (DE) Algorithms for hydraulic analysis on complex pipe network that have 12 (twelve) pipe elements can show excellent performance. The absolute head at each node has similarities with the results obtained from the Newton Raphson method [12] [13] [14]. Genetic Algorithms (GAs), Harmony Search (HS) and Simulated Annealing (SA) have excellent performance for optimizing pipe diameters in drinking water distribution pipelines [15] [16].

This article presents a new model as an alternative solution to solving the problem of hydraulic parameter analysis in complex drinking water distribution networks by utilizing the advantages of metaheuristic methods. In this model, the hydraulic equation system in a pipeline is brought into the form of a non-linear matrix equation system, and then solved using the optimization concept based on the Metaheuristic method. The optimization process uses 3 methods, namely; DE algorithm, PSO algorithm and CODEQ algorithm. As the optimized variable is the head value at each node in the pipe network, as an objective function is to minimize the difference between the discharge outflow at each service node and the targeted outflow or required flow, and as a constraint function is a system of hydraulic equations that applies to the pipe network. The iteration process to find the optimal conditions is done using the M-FILE program code from the MATLAB software. To measure the effectiveness of the model developed from three methods, the results of the analysis obtained are compared with the results of the analysis of the Newton Raphson method and conservative simulations using the Monte Carlo simulation method. The model reliability test uses two hypothetical data sets, namely; 1) simple network data sets and 2) complex network data sets. In addition to knowing the level of effectiveness of metaheuristic methods in solving hydraulic analysis problems in networks, through the results of the analysis of two data sets that have different characters, it is expected to identify the strengths and weaknesses of the models developed.

2. Method and Materials

2.1. Hydraulic Equation System

2.1.1. The Basic Principle of Flow on the Network Pipe

Flow analysis in network pipe must meet the following basic principles of continuity and energy conservation; 1) the flow in a pipe must meet the laws of head loss for a single pipe flow, 2) the amount of flow that enters the network system must be the same as the flow that leaves it, and 3) the flow that enters a node must be the same as the number of flows leave it.

2.1.2. Head Loss in a Pipe Element

The pipe element with the notation i connects vertices j and k as shown in Figure 1. If the head at node j and k is Z_j and Z_k, then the head loss of element i is:

$h_f=Z_k-Z_j=\alpha \cdot L_i\cdot \frac{v_i^2}{2g}$ (1)

where,

h_f—head loss (m);

Z_k—absolute head at node k (m);

Z_j—absolute head at node j (m);

L_i—length of pipe element i (m);

v_i—velocity in pipe element i (m/s);

α—roughness coefficient of pipe element i.

If v_i = velocity in the pipe element i, A_i = cross-sectional area of the pipe i, and Q_i = discharge in element i, then this applies:

$Q_i=A_i\cdot v_i$ (2)

Head in the pipe element is needed to overcome the friction resistance that occurs when flowing. The relationship between high energy loss and discharge in pipe flow, according to Hazen-William is stated by [17]:

$Q=A\cdot v$ (3)

$Q=0.2785\cdot C_{HW}\cdot D^{2.63}\cdot S^{0.54}$ (4)

$Q=\frac{0.2785\left(C_{HW}\cdot D^{2.63}\right)\cdot h{f}^{0.54}}{L^{0.54}}$ (5)

$S=\frac{hf}{L}$ (6)

with,

Q—discharge (m³/s);

D—pipe diameter (cm);

S—slope of the energy grade line;

hf—head loss due to limit friction (m);

L—length of the pipe element (m);

CHW—Coefficient of pipe wall roughness by Hazen William. CHW values can be estimated using Table 1.

Furthermore if,

$k=0.2785\cdot \left(C_{HW}\cdot D^{2.63}/L^{0.54}\right)$, dan $Z=hf$ (7)

then Equation (7) can be simplified into:

$Q^i=k^i\cdot {\left(Z^i\right)}^{0.54}$ (8)

or,

$Q^i=\frac{k^i}{{\left(Z^i\right)}^{0.46}}\cdot \left(Z^i\right)$ (9)

Source: [18].

and,

$Q^i=K{t}^i\cdot \left(Z^i\right)$ (10)

$K{t}^i=\frac{k^i}{{\left(Z^i\right)}^{0.46}}$ (11)

Equation (10) shows that the Ktⁱ value is not constant, but is also influenced by the head loss (Z) so that the analysis will provide a system of non-linear equations. By analogy, the above equation system can also be derived using the basic Darcy-Weisbach equation or other similar equations.

2.1.3. Flow Equation System for a Pipe Element

By using Equation (11) it can be calculated that the discharge of element i is a function of the head in that element. Discharge in an element marked positive (+) if it leaves node k and is marked negative (−) if it goes to node j which corresponds to the head at node k greater than node j. In an incoming debit node is positive (+) and negative (−) when leaving the node.

$Q_{kj}=k^i\cdot {\left(Z_k-Z_j\right)}^{0.54}$ (12)

$Q_{kj}=\frac{k^i}{{\left(Z_k-Z_j\right)}^{0.46}}\cdot \left(Z_k-Z_j\right)$ (13)

then,

$Q_k^i=k{t}^i\cdot d{Z}_i=k{t}^i\left(Z_k-Z_j\right)$ (14)

$Q_j^i=k{t}^i\cdot d{Z}_j=k{t}^i\left(Z_j-Z_k\right)$ (15)

with,

$Q_k^i$ —discharge in node k element i;

$Q_j^i$ —discharge in node j element i.

If the equation system is transformed into a matrix, then it is obtained:

$k{t}^i\cdot \left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\times \left[\begin{array}{c}{Z}_k\\ Z_j\end{array}\right]=\left[\begin{array}{c}{Q}_k^i\\ Q_j^i\end{array}\right]$ (16)

or,

$Q^i=K^{i\ast }\cdot Z^i$ (17)

with,

$Q^i=\left[\begin{array}{c}{Q}_k^i\\ Q_j^i\end{array}\right]$ is a discharge vector element;

$K^{i\ast }=k^i\cdot \left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$ is a characteristic matrix of element I;

$Z^i=\left[\begin{array}{c}{Z}_k^i\\ Z_j^i\end{array}\right]$ is a difference vector of absolute pressure at nodes in element i.

2.1.4. Hydraulic Equation System in a Simple Network

A simple network model is shown in Figure 2. In the figure, Q_n contributes to each node. This Q_n will be positive (+) if the flow enters and negative (−) if the flow leaves the node (there is real use). Based on the description, the summary discharge from the elements connected to the node value must be equal to the number of discharges that exit the node, thus for each node applies:

${\displaystyle {\sum }_i^N{Q}_j^i}=F$ (18)

where the numbers 1, 2, 3, 4 indicate the index node and the letters a, b, c, d indicate the index of the pipe element. Furthermore, by means of a matrix then:

$F_i={\displaystyle \sum Q_k}+{\displaystyle \sum Q{n}_i}=0$ (19)

with, F_i = balance discharge at node i dan ${\displaystyle \sum Q{n}_i}$ = increase or decrease in flow at node i. If the equation is changed in a matrix symbol then it applies:

$F=K\cdot \left(Z\right)\cdot Z$ (20)

The characteristic matrix includes the supporting factors for each element. These supporters are presented in a table according to the symbols that are linked. To facilitate implementation, it is necessary to make a table of the relationship of elements with their nodes. According to the simple model above, we can present the relationship of elements with their nodes as shown in Table 2 and Table 3.

2.2. Hydraulic Equation System Solution Using Newton Raphson Method

This method is a numerical method based approach. The solution to Newton Raphson’s equation for the pipeline hydraulic analysis graphically is shown in Figure 3.

The calculation steps for the Newton Raphson method are systematically described as follows:

1) Set the initial absolute head at each node (Z_initial), [Z_i];

2) Calculate and arrange the gradient line i, or matrix characteristic of network pipe [Kt_i];

3) Calculate external discharge at F_i service points, with the equation $\left[F_i\right]=\left[K{t}_i\right]\ast \left[Z_i\right]$ ;

4) Calculate discharge balance at each node $\left[\Delta F_i\right]=\left[F\right]-\left[F_i\right]$ ;

5) Calculate [ΔZ_i] by solving the equation $\left[K{t}_i\right]\ast \left[\Delta z_i\right]=\left[\Delta F_i\right]$ ;

6) Calculate [Z_i+1] with the equation $\left[Z_{i+1}\right]=\left[Z_i\right]+\left[\Delta Z_i\right]$ ;

7) Repeat step 2) to step 6) until the desired tolerance is found;

8) Set [Z_i+1] at the end of the iteration as the absolute head at the node and is the variable sought;

9) Calculate the hydraulic parameters of the network pipe according to the value of Z_i produced in step 8), including; discharge and absolute head at the nodes, as well as the discharge, velocity and head loss in each pipe element.

2.3. Hydraulic Quation System Solution Using Metaheuristic Method

2.3.1. Optimization Process

In the implementation of metaheuristic methods, the solution of the network pipe hydraulic equation system can be approached through the optimization process. In the metaheuristic method, the objective function is stated as a fitness function. The objective function of the optimization process is the water balance at each service node or can be expressed as the minimization of the deviation of the service nodes from the hydraulic analysis with the targeted node discharge, i.e. according to the data. Thus the objective function or fitness function of the optimization process can be stated as:

$\text{fitness}=F={\text{RMSE}}_{\text{minimum}}={\displaystyle {\sum }_{i=1}^n\frac{\sqrt{{\left(Q_i\left(Z\right)-Q{o}_i\right)}^2}}{n}}$ (21)

where:

F—water balance at service nodes (l/s);

RMSE—root mean square error (l/s);

Z—absolute head at node, (m);

Q_i(Z)—the outflow discharge at node i (l/s) from which is the function Z, and is calculated based on Equation (17);

Qo_i—discharge targeted node i (l/s);

i—node index;

n—number of node.

The optimization process using the metaheuristic method to find the height of head at nodes (Z_i) which produces a minimum RMSE value schematically is shown in Figure 4.

2.3.2. Particle Swam Optimization (PSO) Algorithm

The PSO algorithm as an optimization tool provides a population-based search procedure where each individual called a particle changes their position with respect to time. In the PSO system, each particle flies around a multi dimensional search space and adjusts its position based on personal experience and the experience of the particles next to it. In the case of minimization of a function, the analysis steps are described as follows [19] [20] [21]:

1) Assuming that the size of the group or number of particles is N. To reduce the number of evaluation functions needed to find a solution, it is better if the size of N is not too big but not too small so that there are many possible positions towards the best or optimal solution.

2) Generating an initial population of x with ranges of x^(B) and x^(A) randomly so that we get $x_1,x_2,\cdots ,x_N$. Where x is an optimized variable. In the case of network pipe hydraulic analysis, the optimized variable is the head at each node (Z_i), so x^(B) = the lower limit of the value of Z_i and x^(A) = the upper limit of the value of Z_i. Particle j and its velocity in iteration i are denoted as $x_j^{\left(i\right)}$ and $v_j^{\left(i\right)}$, so these initial particles are denoted: $x_1\left(0\right),x_2\left(0\right),\cdots ,x_N\left(0\right)$.

The vector, $v_j\left(0\right),\left(j=1,2,\cdots ,N\right)$ is the particle or coordinate vector of the particle. Then an evaluation of the objective function value for each particle is stated by:

$f\left(x_1\left(0\right)\right),f\left(x_2\left(0\right)\right),\cdots ,f\left(x_N(\; 0\; )\right)$

The objective function is calculated using Equation (21).

3) Calculate the velocity of all particles. All particles move to the optimal point with a certain velocity. Initially, all the velocity of the particle is assumed to be zero. Set iteration i = 1.

4) In the iteration i, 2 important parameters are found for each particle j, namely:

a) The best value so far from $x_j^{\left(i\right)}$ (particle coordinates j in iteration i) and expressed as P_best,j, with the lowest objective function value, $f\left[x_j^{\left(i\right)}\right]$, which a particle j found in the previous iteration. The best value for all particles $x_j^{\left(i\right)}$ found up to the iteration i, G_best, with the smallest objective function (fitness function) value among all particles for all previous iterations, $f\left[x_j^{\left(i\right)}\right]$.

b) Calculate the velocity of particle j on the iteration i with the following formula:

$v_j\left(i\right)=v_j\left(i-1\right)+c_1{r}_1\left[P_{best,j}-x_j\left(i-1\right)+c_2{r}_2\left[G_{best}-x_j\left(i-1\right)\right]\right],j=1,2,\cdots ,N$ (22)

where c₁ and c₂ are learning rates for individual abilities (cognitive) and social influences (herd), r₁ and r₂ random numbers are uniformly distributed in the 0 - 1 interval. So the parameters c₁ and c₂ indicate the weight of the memory (position) of a particle against the memory of the group (swarm). Values of c₁ and c₂ are usually 2, so the multiplication of c₁r₁ and c₂r₂ ensures that the particles will approach the target by about half the difference.

c) Calculate the position or coordinates of particle j in the iteration i by:

$x_j\left(i\right)=x_j\left(i-1\right)+v_j\left(i\right),j=1,2,\cdots ,N$ (23)

Evaluate the value of the objective function for each particle and is expressed as:

$f\left[x_1\left(i\right)\right],f\left[x_2\left(i\right)\right],\cdots ,f\left[x_N\left(i\right)\right]$ (24)

The objective function is calculated using Equation (21).

5) Check the convergence of the solutions obtained. If all particles go to the same value, it is called convergent. If it hasn’t converged then step 4 is repeated by updating i = i + 1, by calculating the new values of P_best,j and G_best. This iteration process continues until all the particles go to the same solution point. Usually stopping criteria will be found, for example, the number of differences between the current solution and the previous solution is very small.

6) The optimum Z_i is determined based on the best fitness value obtained in the last generation. Hydraulic parameters of the network pipe include; the absolute discharge and energy of the node, and the discharge, flow velocity and compressive height loss in each pipe element are calculated based on the optimum Z_i value.

2.3.3. Differential Evolution (DE) Algorithm

Differential Evolution (DE) algorithm includes stochastic and population-based search methods. The DE algorithm has similarities with other evolutionary algorithms (EA), but is different in terms of distance and direction information from the population that is now used to guide the process of finding a better solution. DE was developed by Reiner Storn and Kenneth Price in 1996. The analysis in the DE Algorithm contains 4 (four) components, namely 1) initialization, 2) mutation, 3) crossover and 4) selection [19] [22] [23]. The relationship of the four components is shown in Figure 5.

1) Initialization

In the case of network pipe hydraulic analysis, the optimized vector variable is the absolute head at each node (Z_j), and j is the node index. Before initializing the vector variable that is sought, it is necessary to determine the lower limit (lb_j) and the upper limit (ub_j) of the optimized variables. The lower limit and the upper limit will be used as a first step to generate the value of the variable being sought. For generation of initial values of the 0^th generation variable, variables to j and vector i can be represented by the following notation.

$x_{j,i,0}=l{b}_j+ran{d}_j\left(0,1\right)\left(u{b}_j-l{b}_j\right)$ (25)

Random numbers are generated by the rand function, where the resulting numbers lie between (0, 1). The index j shows the variable j. In the case of minimization of functions with 2 variables, then j will be worth 1 and 2. Determination of the upper and lower limits is very dependent on the problem being solved.

2) Mutation

After initialization, DE will mutate and combine the initial population to produce a population with the size of an N vector experiment. In DE mutation is done by adding the difference of two vectors to the third vector by:

$v_{i,g}=x_{r_0,g}+F\left(x_{r_1,g}-x_{r_2,g}\right)$ (26)

It appears that the difference between the two vectors chosen at random needs to be scaled before being added to the third vector, $x_{r_0,g}$. The scale factor $F\in \left(0,1\right)$ has a positive real value to control the population growth rate. The base vector index r₀ can be determined in various ways, in general, a different random method is used with the index for the target vector, i. Besides being different from each other and different from the index for the base vector and target vector, the difference vector index r₁ and r₂ are also selected once per mutant.

3) Crossover

At this stage DE crosses each vector x_i,g, with the mutant vector v_i,g, to form the crossing vector, u_i,g with the formula.

$u_{i,g}=u_{j,i,g}=\{\begin{array}{l}{v}_{j,i,g}\to \text{jika}\left(rand\left(0,1\right)\le Cratauj=j_{rand}\right)\\ x_{j,i,g}\to \text{jika}\left(rand\left(0,1\right)>Cratauj\ne j_{rand}\right)\end{array}$ (27)

4) Selection

If the trial vector u_i,g has a value of the objective function that is smaller than the target vector’s objective function x_i,g, then u_i,g will replace the position of x_i,g in the population in the next generation. If the opposite occurs, the target vector will remain in its position in the population.

5) Presentation of results

- The iteration process will stop at the specified stopping criteria, i.e. the maximum number of generations given;

- Optimum Z_i is the x_i value determined based on the best fitness value obtained;

- Calculate pipeline hydraulic parameters at optimum Z_i conditions, including; the absolute discharge and energy of the node, as well as the discharge, flow velocity and compressive height loss in each pipe element.

2.3.4. CODEQ Algorithm

CODEQ is an algorithm proposed by Omran (2008). This algorithm is a synthesis from chaotic search, opposition-based learning, differential evolution and quantum mechanism. CODEQ algorithm is a metaheuristic method that involves population as a solution for continuous problems. Omran (2008) describes the steps of the CODEQ algorithm as follows [24].

1) Optimized variables are expressed as x_i. In the case of a hydraulic network analysis, x_i is the pressure height at each node (Z_i).

2) Initialize population from random vector s in the solution search area.

3) Generating a trial vector v_i(t) for each vector x_i(t) in the iteration t by mutation, according to the formula:

$v_i\left(t\right)=x_i\left(t\right)+\left[x_{i1}\left(t\right)-x_{i2}(t)\right]\ln \left(\frac{l}{u}\right)$ (28)

U is obtained randomly U(0, 1) and $i_1\ne i_2\ne i$.

4) If the fitness function of the trial vector v_i(t) is better than vector x_i(t), then the vector v_i(t) replaces the position x_i(t), if not then vice versa.

5) Generating new vector w(t) on each t iteration with the formula:

$w\left(t\right)=LB+UB-R\cdot x_b\left(t\right)$, jika rand ≤ 0.5 (29)

or,

$w\left(t\right)=x_g\left(t\right)+\left|x_{i1}\left(t\right)-x_{i2}\left(t\right)\right|\cdot \left(2c\left(t\right)-1\right)$, jika rand > 0.5 (30)

r is obtained randomly U(0, 1), LB and UB are the lower and upper limits of the problem, x_b(t) is the worst (least fit) vector in the iteration t, x_g(t) is the best (fitness) vector in the iteration i, x_i1(t) and x_i2(t) are randomly selected vectors with $i_1\ne i_2\ne i$ and c(t) are chaotic variables obtained from the formula:

$c\left(t\right)=\frac{c\left(t-1\right)}{p}jikac\left(t-1\right)\in \left(0,p\right)$ (31)

or,

$c\left(t\right)=\frac{1-c\left(t-1\right)}{1-p}\text{jika}c\left(t-1\right)\in \left(p,1\right)$ (32)

where c(0) and p are obtained randomly at intervals (0, 1).

6) If the fitness function of the new vector w(t) is better than the worst vector x_b(t), then the vector w(t) replaces the position of x_b(t), if not then vice versa. Fitness function is calculated using Equation (21).

7) Steps 2-5 are repeated until the stopping criteria is met.

8) Determine the optimum Z_i is x_b(t) obtained from the best fitness value.

9) Calculate the hydraulic parameters of the pipeline at the optimum Z_i conditions, including; the discharge and absolute head at node, as well as the discharge, velocity and head loss in each pipe element.

2.4. Hydraulic Equation System Solution Using Monte Carlo Method

Monte Carlo simulations have been successfully applied in various fields. The application of the Monte Carlo simulation method to solve the pipeline hydraulic equation system is analyzed through the following five stages:

1) Formulate an equations system of performance for the optimization model to be simulated. In this case, the model performance function is the same as the fitness function in the metaheuristic method as stated in Equation (21), namely the value of the root mean square error (RMSE) which is an indicator of the deviation between the calculated discharge and the targeted node discharge.

2) Generating random numbers with normal, uniform, triangular, beta or other probabilistic distributions.

3) Calculate the appropriate random variable for each optimum variable sought based on the desired number and sample space. As the random variable in this case is the absolute head at the service nodes (Z_i).

4) Evaluate the performance of the model by using input random variable values generated from step 3), including.

- calculation of discharge at nodes as a function of absolute head at nodes $\left[Q_i=f\left(Z_i\right)\right]$ using Equation (17).

- calculation of RMSE values according to equation (21).

5) Repeat steps 3) and 4) as many samples as desired.

6) Determine the optimum Z_i value that occurs at minimum RMSE.

7) Calculate the hydraulic parameters of the network pipe at the optimum Z_i conditions, including; the discharge and absolute head at the nodes, discharge, velocity and head loss in each pipe element.

2.5. Data for Model Testing

2.5.1. Data Characteristic on the Netwok-1

The effectiveness of the hydraulic analysis model developed in this study was tested using a hypothetical data set. There are 5 (five) models tested, namely the hydraulic analysis model based on Newton Raphson method, DE Algorithm, PSO Algorithm, CODEQ Algorithm and Monte Carlo simulation method. The five models are tested on simple network pipe and complex network pipe. A simple network pipe is called network-1, schematically shown in Figure 6. Network-1 consists of 5 service nodes and 6 pipe elements. Water in network pipe comes from reservoirs that have a higher water level than node elevations, so the flow system in the network pipe is only influenced by gravity. Reservoir water level at network-1 is at elevation +10.00 m and service nodes are at elevation +0.00 m. Total targeted outflow discharges at all nodes are 8.00 l/s, consisting of node 3, node 4, node 5 respectively 4.00, 2.00 and 2.00 l/s. The characteristics of nodes and pipe elements in network-1 are presented in Table 4 and Table 5 in detail.

2.5.2. Data Characteristic on the Netwok-2

Complex network data is called network-2, schematically shown in Figure 7. Networks consist of 21 nodes and 32 pipe elements. The water level in the reservoir is +50.00 m, node 1 is +20.00 m and the other nodes are +10.00 m. The total targeted discharge is 165.00 l/s distributed at 30 service nodes as shown in Figure 7. Pipe elements 1 and 2 have a diameter 30 cm and other elements are uniform with diameter 20 cm. Pipe element 1 has length 200 m and 31 other elements have the same length, which is 500 m. The characteristics of the nodes and pipe elements in the networks are detailed in Table 6 and Table 7.

3. Result and Discussion

3.1. Hydraulic Analysis of Simple Network (Network-1)

A comparison of the hydraulic analysis results on network-1 using the 5 methods developed in this research is summarized in Table 8 and Table 9. Table 8 column [2] to column [7] presents a comparison of discharge at node from the

results of analysis with discharge targeted, according to the data shown in Table 5. Hydraulic analysis based on Newton Raphson method, DE Algorithm, PSO Algorithm and CODEQ Algorithm produces a discharge value at the node equal to the targeted discharge value. This shows that the four methods succeeded in completing the hydraulic equation system very well. The discharge at node from the results of a hydraulic analysis based on the Monte Carlo Simulation method shows different results and tends to have significant deviations. Table 9 presents a comparison of the absolute head at nodes, discharge, velocity and head loss in each pipe element produced by the five methods developed. The table shows that the analysis of the Newton Raphson method, the DE Algorithm, the PSO Algorithm and the CODEQ Algorithm tend to produce the same value, but the results of the analysis from the Monte Carlo simulation method show quite significant different results. This indicates that the Monte Carlo simulation method is not relevant to be applied to solve hydraulic equation systems in network pipe even in simple network pipe, and three other metaheuristic methods can demonstrate their effectiveness in solving the problem of simple network pipe hydraulic analysis.

Newton Raphson method for hydraulic analysis in network-1 uses the input of initial head at nodes (Z_i) as shown in Table 5, the maximum number of iterations is set 100 times or stops on stopping criteria (dZ) ≤ 10⁻⁶. dZ is the difference between the absolute number of Z values in the iteration i (Z_i) with the absolute number of Z values in the iteration i + 1 (Z_i+1). Running the program with the input parameter produces a value of dZ close to “0”, which means that the discharge at the service nodes is the same as the targeted discharge. In this regard, it can be said that the estimated Z value based on the results of Newton Raphson’s analysis approaches the exact value. Table 8 column [2] and column [3] present a comparison of the targeted discharge at nodes and the approach discharge at nodes. Table 9 columns [2] [7] and [12] present the hydraulic parameters from the Newton Raphson analysis results. Furthermore, the hydraulic parameters from the Newton Raphson method analysis are used to measure the performance of other methods in solving the problem of hydraulic analysis on network-1.

Hydraulic analysis on network-1 uses an optimization method based on the DE, PSO and CODEQ algorithms involving the same parameter values, namely; the number of individuals in the population (N) = 520, the maximum number of generations (iteration) = 150, the lower boundary of the head at node (lb) = −30.00 m and the upper boundary of the head at node (ub) = 30.00 m. The development of fitness values from generation to generation is shown in Figure 8. Although starting from the best fitness values are different, convergent conditions can be achieved in the 150^th generation. The best fitness values or minimum RMSE that can be achieved by DE, PSO and CODEQ Algorithms are 0.003 l/sec, 0.000 and 0.004 l/sec respectively. Quantitatively it appears that the PSO algorithm is more effective than the other 2 methods, although the difference is not significant. The value of the hydraulic parameters at each node and each pipe element resulting from 3 metaheuristic methods have the same values as the results of the Newton Raphson analysis, as shown in Table 8 and Table 9. This shows that the two methods have the same level of effectiveness in terms of accuracy.

Hydraulic analysis on network-1 using the Monte Carlo simulation method involving parameters, including; population size (nn) = 700,000, number of samples (N) = 500,000, the lower boundary value of head at nodes (lb) = −30.00 m and the upper boundary of head at nodes (ub) = 30.00 m. Figure 9 shows the RMSE values for each sample, after being sorted from the largest value to the smallest value. The best RMSE value is 0.440 l/s. The outflow discharge at each node in converging condition is shown in Table 8 column [7]. Discharge at each node have a significant difference to the targeted discharge. Other hydraulic parameters include; absolute head at each node, as well as discharge, velocity and head loss on each pipe element also show significant differences when compared with the results of the analysis of the Newton Raphson method and the metaheuristic method, as shown in Table 8 and Table 9. These conditions indicate that the Monte simulation method Carlo is not relevant to be applied for hydraulic analysis of pipe network, even in simple networks. This is due to the high sensitivity of the head value at the node to the outflow discharge at each node. The application of Monte Carlo simulation methods with large sample inputs may show better results, but it becomes inefficient in terms of time to achieve convergent conditions.

There are 2 factors that cause the Monte Carlo simulation method is not relevant to be applied for the analysis of hydraulic parameters in the network, namely; 1) variable of head at nodes (Z_i) and variable of outflow discharge at nodes (Q_i) in the network have a relationship with a very high level of sensitivity, so that a small change in the value of the variable Z_i in one node will have a major effect on the value of Q_i in all node, and 2) the value of the Z_i variable at each node is continuous, if the value must be determined conventionally through a trial and error approach it will be difficult to find a convergent condition, because of the unlimited population. In the Newton Raphson method, the Z_i value is found iteratively but systematically through a mathematical approach, so that the greater the number of iterations, the higher the accuracy of the resulting Z_i value. In this case, the Metaheuristic method is used to find the combination of Z values at each node under optimum conditions, to produce the minimum RMSE value. The optimum variable search process in the metaheuristic method combines the local and global search process using certain algorithms according to the method applied, so that the resulting Z_i value will improve from generation to generation, as shown in Figure 8. The greater the number of generations (iteration) involved, systematically will produce a better Z_i value.

The Monte Carlo method has a slight difference from the trial and error method. The effectiveness of the calculation process on the Monte Carlo simulation method is done by limiting the population by taking a number of samples. The sample value is chosen randomly by generating random variables with normal distribution or other distributions according to the specified lower boundary and upper boundary values. The number of samples and the width of the range between upper boundary and lower boundary will greatly determine the effectiveness and accuracy of the resulting Z_i variable. Various combinations of Z_i values according to the number of samples are used to calculate the Q_i value at each node to produce an RMSE value. Because the Z_i value is determined randomly, the resulting RMSE value is also random. The process of finding the Z_i value is not systematic, the Z_i value of the n^th sample does not guarantee that it can improve the Z_i value of the n + 1 sample. In Figure 9, the RMSE values have been sorted from the largest value to the smallest value. Monte Carlo simulation methods may be able to show effective work if the lower and upper limits are given precisely, and the number of samples is given in large quantities. However this has become inefficient in terms of iteration time.

3.2. Hydraulic Analysis of the Complex Network (Network-2)

The results of analysis network-2 from the 5 methods developed indicate that Newton Raphson method is still more effective than 4 other methods. Figure 10 shows the comparison of the mean deviation of discharge at nodes from the analysis results with the targeted discharge according to the data shown in Table 7 in the line [3]. In terms of the accuracy of the discharge at node from Newton Raphson method gives a very small deviation value, which means it can show very good performance. The Monte Carlo method produces the greatest deviation which means it shows very poor performance. The analysis results from method base on DE algorithm and CODEQ algorithm show an equivalent deviation value and are slightly better than the PSO Algorithm. Figure 11 shows the mean deviation between the absolute head at the node, and the discharge, velocity and head loss in the pipe element from the results of the metaheuristic method analysis of the same parameters from the Newton Raphson method. The graph shows that the DE, PSO and Codeq algorithm-based methods produce an equivalent deviation but the results of the analysis of the Monte Carlo method show much greater results.

The analysis of network-2 using Newton Raphson method is done by inputting the initial absolute head at the nodes (Z) as shown in Table 7, the maximum number of iterations is set 100 times or at the stopping criteria dZ = 10⁻⁶. From the results of the analysis along 100 iterations, convergent conditions can be achieved. Table 8 column [2] and column [3] show that the discharge node from the analysis results is close to the targeted discharge as shown in Table 7 in the column [3]. The mean of deviation of discharge at node is obtained 0.186 l/s with a standard deviation of 0.37 l/s indicating that the deviation is very small and occurs evenly at all nodes. This indicates that the Newton Raphson Method is very relevant and effective in solving flow analysis problems in networks. The results of the analysis of other hydraulic parameters are presented in Table 10 columns [8] and Table 11 columns [2] [7] and [12]. The results of the analysis from Newton Raphson method are then used as a benchmark to measure the performance of other methods in the analysis of network-2.

The application of DE algorithm on network-2 analysis uses the parameter input the number of individuals in the population (N) = 2200, the maximum number of generations (iteration) = 2000, the lower boundary (lb) = −30.00 m and the upper boundary (ub) = 30.00 m. The analysis result gives the best fitness value 1.7215, the progress of fitness value from generation to generation is shown in Figure 12. The discharge at node from the analysis results is shown in Table 10 column [4]. There is a difference in the value of the discharge at node from the analysis with the targeted discharge, namely with an average deviation of 6.02 l/s and a standard deviation of 5.22 l/s. This value is greater than the value generated by Newton Raphson method. This indicates that the DE algorithm is quite effective in solving the problem of network-2 hydraulic analysis, but its performance is not as good as the results of the Newton Raphson method. The absolute head at nodes, discharge, velocity, and head loss in each pipe element also shows different results when compared to the results of Newton Raphson method.

The mean deviations of the four hydraulic parameters were 6.06 m, 34.2 l/s, 1.09 m/s, 6.94 m, respectively. The comparison of the hydraulic parameter values from the results of the analysis based on the DE Algorithm with the results of the analysis of other methods is shown in Table 11 and Figures 13-17.

The hydraulic analysis on the network-2 using PSO Algorithm involves parameters of the number of individuals in the population (N) = 1000, the maximum number of generations (iteration) = 2000 or at control stopping 10⁻⁸, the lower boundary (lb) = −30.00 m and the upper boundary (ub ) = 30.00 m. In the 2000 iteration the best fitness value was 4.0575. This value is greater than the fitness value generated from the analysis based on the DE Algorithm. The progress of fitness values from generation to generation is shown in Figure 12 and the discharge at nodes from the analysis are presented in Table 10 column [5]. The mean deviation of discharge at node from the analysis with the targeted discharge at nodes is 15.79 l/s with a standard deviation of 10.07 l/s. This shows that in terms of the accuracy of the discharge at nodes the PSO algorithm can work quite well, but it is not better than Newton Raphson method and the DE algorithm. The absolute head at nodes, discharge, velocity and head loss in each pipe element also shows different results when compared with the results of the analysis of Newton Raphson method. The mean deviations of the four hydraulic parameters were 4.29 m, 23.16 l/s, 0.7 m/s, 3.87 m respectively. This value is slightly better when compared to the hydraulic parameters generated by method based on the DE algorithm. The comparison of the hydraulic parameter values from the PSO Algorithm analysis with the analysis from other methods is shown in Table 11 and Figures 13-17.

The application of the CODEQ algorithm for hydraulic analysis on the network-2 uses parameters of the number of individuals in the population (N) = 1000,

the maximum number of generations (iteration) = 2000, the lower boundary (lb) = −30.00 m and the upper boundary (ub) = 30.00 m. Analysis until the 2000^th iteration produced the best fitness value (RMSE) of 1.5728 l/s, and the progress of fitness values from generation to generation is shown in Figure 12. Discharge at nodes from the analysis as shown in Table 10 column [6] shows the difference between the discharge at nodes from the analysis results with the targeted discharge (according to the data), ie with an average deviation of 5.87 l/s and a standard deviation of 4.28 l/s. This shows that in terms of the accuracy of the node discharge the CODEQ algorithm has better performance than the DE and PSO algorithms, but it is not more effective than Newton Raphson method. The absolute head at nodes, discharge, velocity, and head loss in each pipe element also shows different results when compared with the analysis results from Newton Raphson method. The mean deviation of the four hydraulic parameters was 3.13 m, 26.00 l/s, 0.82 m/s, 3.21 m, respectively. The average value is slightly smaller when compared with the hydraulic parameters produced by the DE algorithm and the PSO algorithm. The comparison of hydraulic parameter values from the CODEQ algorithm analysis results with the analysis results from other methods is shown in Table 11 and Figures 13-17.

The hydraulic analysis of network-2 using the Monte Carlo Simulation method involves parameters of the amount of data (N) = 800,000, the number of samples (nn) = 400,000, the lower boundary of Z_i (lb) = −30.00 m and the upper boundary of Z_i (ub) = 30.00 m. The analysis produces the best performance indicators or produces a minimum RMSE value of 23.8981 l/s obtained from the 74,567 sample. The graph of the RMSE values analyzed from each sample and sorted from the maximum to the minimum value is shown in Figure 18. The analysis results of the discharge at nodes are shown in Table 10 column [7]. There is a very significant difference between the analyzed discharge at nodes and the targeted discharge, namely the average deviation of 76.64 l/s and the standard deviation of 80.17 l/s. The magnitude of the deviation indicates that the Monte Carlo simulation method has a poor performance, even worse than 4 other methods. Absolute head at nodes, discharge, velocity and head loss in each pipe elements from the analysis results of the Monte Carlo simulation method

with the analysis results from Newton Raphson method also have a very large deviation, respectively 8.85 m, 54.16 l/s, 1.71 m/s, 10.97 m. Comparison of hydraulic parameter values from the network hydraulic analysis with the Monte Carlo simulation method with the analysis results from other methods is shown in Figures 13-17. In these graphs, it appears that the curve of the five hydraulic parameters produced by the Monte Carlo simulation method does not show a good trend when compared to the other four analysis methods. This further reinforces the notion that the Monte Carlo simulation method is not relevant to be applied for flow analysis in network pipe.

4. Conclusions

Newton Raphson method can show very consistent performance in solving the problem of hydraulic parameter analysis in drinking water distribution networks. In terms of accuracy and speed towards convergent conditions, this method is reliable. The optimization model based on DE, PSO and CODEQ algorithms for hydraulic parameter analysis in pipelines has an equivalent level of performance, but not as effective as Newton Raphson method. In the case of simple pipe network, three optimization models based on metaheuristic methods are reliable, but in the case of complex pipe network the three methods are inefficient in terms of accuracy and speed of reaching convergent conditions. So, the effectiveness of applying the metaheuristic method is largely determined by the complexity of the pipe network analyzed.

The Monte Carlo simulation method shows very poor performance, even in the case of a simple network. This shows that the Monte Carlo Simulation method is not relevant for the analysis of hydraulic parameters in pipe network for drinking water distribution. There are at least 2 factors that cause the metaheuristic method and the Monte Carlo simulation method to be less effective in solving hydraulic equation systems in complex pipe network, namely; 1) the head at node variable (Z_i) as input and the outflow discharge variable at node (Q_i) as output in the pipe network have a relationship with a very high level of sensitivity, so that a small change in the value of the variable Z_i at one node will have a large effect on the value of Q_i in all nodes, and 2) the Z_i variable in each node is a continuous variable. On the other hand, the metaheuristic method and the Monte Carlo simulation method are search methods based on sampling resulting from random variables with certain value constraints, while the optimized variables are continuous and have a high level of sensitivity. This contrast difference makes it difficult for the Metaheuristic and Monte Carlo simulation methods to find accurate solutions, especially in complex pipe network. The use of multiple precision variables, narrowing the search space by providing appropriate lower and upper bound values, and increasing the number of iterations can improve the performance of both methods.

Acknowledgements

The author would like to thank the Directorate of Research and Community Service at the University of Muhammadiyah Malang, hopefully this small work will contribute positively to the development of science and technology.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ormsbee, L.E. (2006) The History of Water Distribution Network Analysis: The Computer Age. 8th Annual Water Distribution Systems Analysis Symposium, Cincinnati, 27-30 August 2006, 1-6.
[2]	Sarbu, I. (2011) Nodal Analysis Models of Looped Water Distribution Networks. ARPN Journal of Engineering and Applied Sciences, 6, 115-125. http://www.arpnjournals.com
[3]	Moosavian and Jaefarzadeh (2014) Hydraulic Analysis of Water Supply Networks Using a Modified Hardy Cross Method. IJE Transactions C: Aspects, 27, 1331-1338. https://doi.org/10.5829/idosi.ije.2014.27.09c.02
[4]	Halagalimath, S.S., Vijaykumar, H. and Patil, N.S. (2016) Hydraulic Modeling of Water Supply Network Using EPANET. International Research Journal of Engineering and Technology, 3, 1022-1027. http://www.irjet.net
[5]	Al-Hashim, Nassrin and Al-Mansori (2011) Analysis and Evaluation of the Potable Water Network and Water Quality in Al-Dinawiya City. Diyala Journal of Engineering Sciences, 5, 191-204.
[6]	Georgescu, A.M., Perju, S., Georgescu, S.C. and Anton, A. (2014) Numerical Model of a District Water Distribution System in Bucharest. Procedia Engineering, 70, 707-714. http://www.sciencedirect.com https://doi.org/10.1016/j.proeng.2014.02.077
[7]	Abdy Sayyed, M.A.H., Gupta, R. and Tanyimboh, T.T. (2014) Modelling Pressure Deficient Water Distribution Networks in EPANET. Procedia Engineering, 89, 626-631. http://www.sciencedirect.com https://doi.org/10.1016/j.proeng.2014.11.487
[8]	Ayad, A., Awad, H. and Yassin, A. (2013) Developed Hydraulic Simulation Model for Water Pipeline Networks. Alexandria Engineering Journal, 52, 43-49. http://www.elsevier.com/locate/aej http://www.sciencedirect.com https://doi.org/10.1016/j.aej.2012.11.005
[9]	Izinyon, O.C. and Anyata, B.U. (2011) Water Distribution Network Modelling of a Small Community Using Watercad Simulator. Global Journal of Engineering Research, 10, 35-47.
[10]	Mehta, V.N. and Joshi, G.S. (2019) Design and Analysis of Rural Water Supply System Using Loop 4.0 and Water Gems V8i for Nava Shihora Zone I. International Journal of Engineering and Advanced Technology, 9, 2258-2266. https://doi.org/10.1155/2014/979706
[11]	Moosavian, N. and Jaefarzadeh, M.R. (2014) Hydraulic Analysis of Water Distribution Network Using Shuffled Complex Evolution. Journal of Fluids, 2014, Article ID: 979706. https://doi.org/10.1155/2014/979706
[12]	Sulianto (2015) Linear Programming to Search Optimum Diameter Pipe in Network Pipe Open in Water Supply System. Media Teknik Sipil, 13, 91-98. http://ejournal.umm.ac.id/index.php/jmts/article/view/1216
[13]	Sulianto (2015) The Search Optimum Diameter on Open Network Pipe System Using GA. Proceeding National Conference of Civil Engineering XI 2015, Surabaya, 28 January 2015, 191-204. https://fdocumente.com/document/prosiding-semnas-xi-2015.html
[14]	Sulianto (2015) Flow Analysis on Pipe Distribution Network Using Differential Evolution. Proceeding the 1st Young Scientist International Conference of Water Resources Development and Environmental Protection, Malang, 5-7 June 2015, B54-B62.
[15]	Saldarriaga, J., Páez, D., León, N., López, L. and Cuero, P. (2014) Historical Development of Power Use Methods for WDS Design and Their Evolution towards Optimization Metaheuristics, Numerical Model of a District Water Distribution System in Budapest. Procedia Engineering, 70, 707-714. http://www.sciencedirect.com https://doi.org/10.1016/j.proeng.2014.02.162
[16]	Reca, J., Gil, C., Martínez, J. and Banos, R. (2008) Application of Several Metaheuristic Techniques to the Optimization of Real Looped Water Distribution Networks. Water Resources Management, 22, 1367-1379. https://www.researchgate.net/publication/248070757 https://doi.org/10.1007/s11269-007-9230-8
[17]	Adeleke and Olawale (2013) Computer Analysis of Flow in the Pipe Network. Transnational Journal of Science and Technology, 3, 45-71.
[18]	Jolly, M.D., Lothes, A.D., Bryson, L.S. and Ormsbee, L. (2014) Research Data Base of Water Distribution System Models. Journal of Water Resources Planning and Management, 140, 410-416. https://doi.org/10.1061/(ASCE)WR.1943-5452.0000352
[19]	Santosa, B. and Willy, P. (2011) Metoda Metaheuristik Konsep dan Implementasi, Indonesian. Guna Widya Press, Surabaya.
[20]	Santos, C.A.G., Macedo, P.K.D., Freire, M., Mishra, S.K. and Junior, A.S. (2011) Application of a Particle Swarm Optimization to the Tank Model. Proceedings of Symposium H03, Melbourne, July 2011, IAHS Publ. 347.
[21]	Mtolera, I., et al. (2014) Optimization of Tree Pipe Networks Layout and Size Using Particle Swarm Optimization. WSEAS Transactions on Computers, 13, 219-230.
[22]	Karabo, G.D. and Okdem, S. (2004) A Simple and Global Optimization Algorithm for Engineering Problems: Differential Evolution Algorithm. Turkish Journal of Electrical Engineering and Computer Science, 12, 53-60.
[23]	Moosavian, N. and Jaefarzadeh, M.R. (2013) Pressure-Driven Demand and Leakage Simulation for Pipe Networks Using Differential Evolution. World Journal of Engineering and Technology, 1, 49-58. https://doi.org/10.4236/wjet.2013.13008
[24]	Omran, M.G.H. and Salman, A. (2009) Constrained Optimization Using CODEQ. Journal of Chaos, Solutions and Fractals, 42, 662-668. http://www.elsevier.com/locate/chaos https://doi.org/10.1016/j.chaos.2009.01.039