## Abstract

Pumps are installed in water distribution networks (WDNs) to ensure adequate service levels in the case of poor water pressure (e.g. because of low elevation of reservoirs or high head losses within the WDN). In such cases optimal pump scheduling is often required for the opportunity of significant energy saving. Optimizing the pump operation also allows a reduction in damage and maintenance times. Among the approaches available in the literature to solve the problem, meta-heuristic algorithms ensure reduced computational times, although they are not able to guarantee the optimal solution can be found. In this paper, a modified Harmony Search Multi-Objective optimization algorithm is developed to solve the pump scheduling problem in WDNs. The hydraulic solver EPANET 2.0 is coupled with the algorithm to assess the feasibility of the achieved solutions. Hydraulic constraints are introduced and penalties are set in case of violation of the set constraints to reduce the space of feasible solutions. Results show the high performances of the proposed approach for pumping optimization, guaranteeing optimal (or near optimal) solutions with short computational times.

## INTRODUCTION

Development of optimization procedures for water distribution networks (WDNs) is a relevant topic in the field of the water systems management because of both technical and economic impacts. The aim is to improve the management of the systems by achieving the maximum benefit with minimum cost; thus effective design and management of systems is required.

Over the last decades, several applications were developed as optimization procedures for hydraulic systems. For complex systems, a high number of variables in optimization algorithms dramatically increases the number of functions to be evaluated, thus requiring huge computational times. The application of meta-heuristic approaches allows good results to be reached by improving the computational efficiency (Creaco & Pezzinga 2015). Several models are available in the literature to solve the pump scheduling problem (e.g. Beckwith & Wong 1995; Mackle *et al.* 1995; Savic *et al.* 1997; Lopez-Ibanez *et al.* 2008), also for pump-and-treat remediation of a contaminated aquifer (Liu *et al.* 2000) or irrigation systems (Moradi-Jalal *et al.* 2004).

In the field of natural-based meta-heuristic models, Geem *et al.* (2001) implemented an optimization procedure, named Harmony Search (HS), able to conceptualize the musical process of searching for a perfect state of harmony, in analogy with the musical harmony, soundly played by a jazz improvisation. The harmony in music was analytically represented by a solution vector and the musician's improvisation was modelled through local and global schemes of optimization techniques.

The HS optimization algorithm has been applied to several applications in different operative fields De Paola *et al.* (2016a, 2016b, 2017a, 2017b). Among these, Geem (2005) applied the HS algorithm to the pump switching problem in a water supply system, considering Objective Function (OF) as the minimization of the annual energy costs. He compared the reached results with further meta-heuristic approaches, such as Genetic Algorithm and branch and bound methods, achieving better results in terms of minimal energy, mean energy consumptions and computational times.

More recently, Kougias & Theodossiou (2013) extended the HS to optimize the pump scheduling in WDNs using the Polyphonic-HSA variant. De Paola *et al.* (2016b) coupled the Harmony Search Multi-Objective (HSMO) algorithm and EPANET 2.0 (Rossman 2000) to solve the pump scheduling problem. The ‘classical’ formulation of the HS proposed by Geem (2005) was used. The model was tested on the Anytown WDN (Pasha & Lansey 2009) as modified by Savic *et al.* (2011).

Lopez-Ibanez (2009) proposed a genetic algorithm to solve the pump scheduling problem, in which the maximum number of switches for each pump was set. An approach based on the Relative Time Trigger Control (RTTC) was proposed to improve the efficiency of the algorithm, although solutions with warnings resulting from the hydraulic simulations were not discarded, thus reducing the number of feasible solutions.

In the present work, a modified HSMO algorithm was implemented for Pump Scheduling optimization. Hydraulic simulations and physical constraints were performed by using the hydraulic solver EPANET 2.0. The maximum number of switches for each pump was set and unfeasible solutions violating hydraulic and physical constraints were *a-priori* discarded in order to generate only reliable solutions. Unlike the approach proposed by De Paola *et al.* (2016b), the decision variables were set considering relative-time intervals (Lopez-Ibanez 2009).

The model was tested on both the modified Anytown WDN and the network proposed by Van Zyl *et al.* (2004). Results were also compared with the GaNetXL model (Savic *et al.* 2011) in the first test case, with the results obtained by Van Zyl *et al.* (2004) and Lopez-Ibanez (2009) in the second test case.

## COMPUTATIONAL METHODOLOGY

The HS algorithm reproduces the research of a ‘perfect state’ of harmony, in such a way similar to what happens when a musical jazz improvisation is played. The fantastic musical harmony, generated by an aesthetic estimation, is similar to the algorithm research of the global optimum, evaluated through the calculation of the OF. Each set of values, attributed to the decision variables, corresponds to the set of pitches played by musical instruments. Hence, the perfected harmony is reached only by performing a sufficient number of practices, corresponding to iterations for the algorithm. The algorithm is numerically composed of five steps (Geem *et al.* 2001):

*Step 1*: Initialization of optimization problem and HS algorithm parameters;*Step 2*: Initialization of Harmony Matrix (HM);*Step 3*: Improvisation of a new harmony from HM;*Step 4*: Updating HM;*Step 5*: Checking the stopping criterion.*Step 1*: The parameters required for optimization are specified. They are the Harmony Memory Considering Rate (HMCR), the Pitch Adjusting Rate (PAR), the Harmony Memory Size (HMS: number of solution vectors) and the termination criterion, which corresponds to the total number of improvisations (i.e. the number of iterations). HMCR and PAR are set in Step 3, allowing the improvement of the solution vectors. As Multi-Objective algorithm, a Pareto optimal set is made up of solutions that are not dominated by any other solution in the population. The graphical representation in the space of all solutions, which belongs to the Pareto optimal set, determines the Pareto front. In the present work, an archive set named Harmony Matrix (HM) is introduced to record the non-dominated solutions. If one solution is found to be dominant on one or more solutions during the process, it will replace them into the HM.*Step 3*: A new harmony vector,*x′**=*(*x′*, … ,_{1}, x′_{2}*x′*) is generated, according to three rules: Memory Considerations, Pitch Adjustments and Randomization. The value of the first decision variable (_{N}*x′*) can be chosen from any value in the specified HM range (_{i}*x*) and a probability of choosing a totally random value for the decision variable_{i}^{1}~ x_{i}^{HMS}*x′*is also considered: being HMCR (between 0 and 1) the probability of selecting one value from the historical values stored in the HM and the complement (1−HMCR) the probability of choosing a random feasible value, not limited to those stored in the HM._{i}Values of the other decision variables (

*x′*, … ,_{2}*x′*) can be chosen in the same way. Following the Memory Considering Operation, the Pitch Adjusting Operation starts. In the classical formulation of the HS, it operates only for values chosen from the HM. Such an operator applies the PAR to set the rate of moving to neighbouring values for the initially chosen value from the HM._{N}*Step 4*: From the calculation of the OF, if the new harmony vector*x′**=*(*x′*, … ,_{1}, x′_{2}*x′*) has a better result than the worst harmony in the HM, it is included into the HM and the existing worst harmony is consequently excluded from it._{N}*Step 5*: Computation is stopped when the termination criterion (i.e. the number of iterations) is satisfied. If not, Steps 3 and 4 are repeated.

In the proposed approach (Figure 1), the Pitch Adjusting Operation has been included into the algorithm, because the candidate values are chosen using the relative-time procedure proposed by Lopez-Ibanez (2009).

*C*is the total energy cost for a 24-hour period,

_{T}*TN*is the total number of daily pump switches,

_{SW}*t*is the hourly time step,

*C*(

*t*) is the unit cost during time step

*t*,

*E*(

*t*) is the energy used in the time step

*t*(depending on the pump flow

*Q*(

_{p}*t*) and water level

*H*(

*t*) in tanks) and

*N*(

_{SW}*t*) is the number of pump switches during the time

*t*. In the model, the pump switch is intended as the turning on or off of a pump which is, respectively, turned off or on during the previous time step. To avoid the frequent turning on and off caused by a high number of switches, the algorithm tends to reduce the number of starts for each pump

*TS*and, consequently, the total starts in the system

^{P}*TS*, is evaluated as follows: where

*N*is the total number of pumps in the system.

^{p}*TS*, and the following constraint is

^{p}*a-priori*introduced to limit the number of starts of each pump to a fixed threshold

*ST*: A further constraint concerns the water level into the tank, which is set to range, for each time step

*t*, between a minimum storage head

*H*and a maximum storage head

_{min}*H*: The water level depends on both the tank level at the previous time step

_{max}*H*(

*t*–1) and the flow rate

*Q*(

_{p}*t*) at the current time step

*t*entering or leaving the tank.

## REPRESENTATION OF SCHEDULES

The most important change proposed in the paper is the use of relative-time intervals (Lopez-Ibanez 2009) instead of a binary approach to assess the pump state. Each pair of decision variables represents the time during which a pump is inactive and active, respectively. According to this definition, a pair of decision variables (*t _{i},t′_{i}*) represents a single pump start, because it specifies the transition from the inactive status (during

*t*) to the active mode (during

_{i}*t′*). As a consequence, the sum of all time intervals for each pump has to be smaller than or equal to the scheduling period

_{i}*T*.

*T*, schedules when the pump is not active at the end of the scheduling period can be represented as: When relative times are applied to the time-controlled triggers representation, the most significant implementation issue is the range of decision variables. In this case, zero-length intervals are allowed for any decision variable, and thus the range of feasible values is [0,

*T*]. Such a formulation enables the representation of schedules with pump starts smaller than or equal to

*ST*(

*TS*

^{p}*≤*

*ST*): where

*s*is the sequence of decision variables representing the schedule of pump

^{p}*p*and

*K*is an integer number between 0 and

*T*. As an example, in Figure 2 a time-controlled trigger representation is plotted for

*ST*= 3.

The three pairs of decision variables [(0,2), (8,3), (8,3)] indicate that the pump turns on at time 0 for 2 hours, then is turned off for 8 hours, turned on for 3 hours, turned off for other 8 hours and finally turned on for 3 hours. The sum of all the time intervals is 24 h, according to Equation (10).

^{24}= 16,777,216 reliable solutions. However, if the number of pump starts is restricted to three (

*TS*≤ 3), the feasible searching space, with respect to this constraint, is limited to 290,998 solutions, corresponding to the 1.74% of the total searching space (Lopez-Ibanez 2009). The number of solutions with

^{p}*TS*not higher than

^{p}*ST*can be calculated as:

Moreover, instead of adding penalty mechanisms to discard solutions with unfeasible tank levels and supplied demands as in Kougias & Theodossiou (2013), the range of feasible solutions was reduced by excluding all the unfeasible solutions (see Equations (7)–(9)) from the process.

## ANYTOWN WDN CASE STUDY: RESULTS AND DISCUSSION

To assess the effectiveness of the proposed model, numerical simulations were performed for the Anytown WDN (Pasha & Lansey 2009), as modified by Savic *et al.* (2011). The WDN (Figure 3) is supplied by a reservoir (ID 20) and four pumps (IDs P1 to P4) operating in parallel.

The characteristic curves of the four pumps are plotted in Figure 4(a) and 4(b).

A two-part daily tariff was considered for energy costs, for all pumps, with a unit price of 0.0244 £/kWh for time steps between 0:00 and 7:00 and 0.1194 £/kWh between 7:00 and 24:00.

The tank (ID 21) has an elevation of 65.53 m a.s.l. and diameter of 12.20 m. Minimum and maximum tank heads were set to 67.67 m a.s.l. and 76.20 m a.s.l., respectively. The solver EPANET 2.0 was used for hydraulic simulations. If the constraints given by Equations (7) and (8) were not satisfied, the algorithm added 9,999 to the total cost, so as to discard unfeasible solutions. If the constraint given by Equation (9) was not satisfied (i.e. the solution returned flow discharges running the pump outside its operating range), the algorithm went back to Step 1 to modify the starting solution. Simulations were performed considering a time interval of 1 hour for both time steps and pattern.

Node demand was modelled considering the daily pattern given in Figure 5, according to Savic *et al.* (2011). The pattern has demand factors ranging between 0.40 and 1.20.

According to Equation (6), four groups of simulations were performed by setting a maximum number of switches per pump equal to 3, 4, 5 and 6, respectively. As termination criterion, a maximum number of iterations equal to 0.5% of the solutions having *TS ^{p}* ≤

*ST*was set for each scenario, as summarized in Table 1. Simulations were performed using a workstation with a 32-bit operating system, 4GB RAM and Intel CORE 2 Quad Q8200 2.33 GHz processor. In Table 1, the maximum number of iterations, the lowest OF value, the total number of switches, the daily operating hours, the iteration number in which the lowest value of the OF is achieved and the duration of simulations are summarized for each scenario.

Max no. of starts per pump . | No. of iterations (–) . | Total daily energy costs (£/day) . | Total no. of switches (–) . | Daily operating time (h) . | No. of iterations (–) . | Time duration (hh:mm:ss) . |
---|---|---|---|---|---|---|

3 | 1,455 | 1,190.64 | 8 | 41 | 1,449 | 00:15:35 |

4 | 8,810 | 933.92 | 10 | 30 | 8,801 | 00:40:40 |

5 | 28,422 | 895.16 | 10 | 32 | 26,917 | 01:31:09 |

6 | 55,464 | 871.25 | 14 | 34 | 51,961 | 03:20:40 |

Max no. of starts per pump . | No. of iterations (–) . | Total daily energy costs (£/day) . | Total no. of switches (–) . | Daily operating time (h) . | No. of iterations (–) . | Time duration (hh:mm:ss) . |
---|---|---|---|---|---|---|

3 | 1,455 | 1,190.64 | 8 | 41 | 1,449 | 00:15:35 |

4 | 8,810 | 933.92 | 10 | 30 | 8,801 | 00:40:40 |

5 | 28,422 | 895.16 | 10 | 32 | 26,917 | 01:31:09 |

6 | 55,464 | 871.25 | 14 | 34 | 51,961 | 03:20:40 |

Table 1 points out the capability of the proposed approach, regardless of the maximum number of starts per pump. The lowest energy cost per day was achieved for the simulation with *ST* = 6. As expected, the worst solution was achieved for *ST* = 3 and the best for *ST* = 6, with a difference of 319.39 £/day, corresponding to 36.6% of the energy cost. Scenarios with *ST* = 4 and *ST* = 5 returned intermediate solutions, with costs 7.2 and 2.7% higher than *ST* = 6, respectively.

The scenario with *ST**=* 6 returned the highest total number of switches (14), whereas for *ST* = 3 only eight switches were calculated. Both *ST* = 4 and *ST* = 5 returned 10 total switches. For all scenarios, the best solution was reached almost at the end of the iterations. The computational effort varied exponentially at increasing *ST*, ranging between around 16 min for *ST* = 3 and 201 min for *ST* = 6.

For example, in Figure 6 the pumps running at any time step are plotted for the best solution achieved with *ST* = 6. As mentioned above, the model returned a total number of switches *ST _{t}* = 14.

All scenarios show a fairly similar trend in achieving the best solution. Figure 7 shows the cost rate (i.e. the ratio between the cost at a generic iteration and the minimum cost) at varying the iteration rate (i.e. the ratio between the generic iteration and the maximum iteration) for all the analyzed scenarios. Nevertheless, scenarios with a greater number of starts achieve the best solutions in a smaller (relative) time with respect to the scenarios with a lower number of starts.

Variation of the water level in Tank 21 during the 24-h simulations is given in Figure 8 for each scenario. Constraints (7) and (8) were always satisfied, since the water level at the end of the simulations was higher than the initial level.

Results were also analysed through the Empirical Attainment Function (EAF), which returns the boundary separating points that were known to be attainable (i.e. dominated in Pareto sense) in at least a fraction (quantile) of the runs from those that were not attainable. The EAF inferred from the four scenarios is plotted in Figure 9. The median EAF represents the curve where the fraction of attainable points is equal to 50%. These plots may be useful for exploring the performance of stochastic global search algorithms for two-objective optimization problems, helping to identify the algorithmic behaviours in a graphical way.

Results were compared with those returned by the GaNetXL model. For the same number of total switches, the proposed model always returns better solutions than GaNetXL. Values are given in Table 2 for the total number of switches *ST _{t}* ranging between 4 and 14.

Total no. of switches (–) . | Total daily energy costs (£/day) . | Relative errors (%) . | |
---|---|---|---|

Proposed approach . | GaNetXL . | ||

14 | 871.25 | 1,078.07 | –19.2 |

12 | 905.00 | 1,100.26 | –17.7 |

10 | 895.16 | 1,181.75 | –24.3 |

8 | 931.23 | 1,259.15 | –26.0 |

6 | 910.20 | 1,417.87 | –35.8 |

4 | 1,078.62 | 1,667.10 | –35.3 |

Total no. of switches (–) . | Total daily energy costs (£/day) . | Relative errors (%) . | |
---|---|---|---|

Proposed approach . | GaNetXL . | ||

14 | 871.25 | 1,078.07 | –19.2 |

12 | 905.00 | 1,100.26 | –17.7 |

10 | 895.16 | 1,181.75 | –24.3 |

8 | 931.23 | 1,259.15 | –26.0 |

6 | 910.20 | 1,417.87 | –35.8 |

4 | 1,078.62 | 1,667.10 | –35.3 |

## RESULTS AND DISCUSSION

A second case study was also developed, by considering the WDN proposed by Van Zyl *et al.* (2004) (Figure 10).

The WDN is supplied by a reservoir (ID 100) with a constant head of 20 m a.s.l. and two tanks (Tank A and Tank B) at elevations of 80 and 85 m a.s.l., respectively. Three pumps (Pump 1A, 2B, 3B) were installed: Pump 1A and 2B have the same characteristic curves and are supplied by the reservoir, whereas Pump 3B is a booster pump installed at higher elevation (100 m a.s.l.). The characteristic curves of the pumps are plotted in Figure 11.

Two demand nodes were considered (IDs 5, 6), with the 24-h demand pattern plotted in Figure 12 (Van Zyl *et al.* 2004). The energy tariff pattern was the same considered for the Anytown WDN.

Results were also compared with the model of Lopez-Ibanez (2009). Similarly to this one, *ST* = 4 was considered and 20,000 iterations were set as termination criterion. Simulations returned the best solution with eight total switches and a daily cost of 323.5 £/day.

Comparison with methodologies of both Van Zyl *et al.* (2004) and Lopez-Ibanez (2009) showed the effectiveness of the proposed approach (Table 3). The minimum cost is lower than the solution found by Van Zyl *et al.* (2004) and only slightly greater than the solution of Lopez-Ibanez (2009). Lopez-Ibanez (2009) found the same number of switches, whereas 10 switches resulted from the approach of Van Zyl *et al.* (2004).

Methodology . | Total daily energy cost (£/day) . | Total number of switches (–) . |
---|---|---|

Van Zyl et al. (2004) | 344.2 | 10 |

Lopez-Ibanez (2009) | 322.5 | 8 |

Proposed approach | 323.5 | 8 |

Methodology . | Total daily energy cost (£/day) . | Total number of switches (–) . |
---|---|---|

Van Zyl et al. (2004) | 344.2 | 10 |

Lopez-Ibanez (2009) | 322.5 | 8 |

Proposed approach | 323.5 | 8 |

The operation of the pumps is plotted in Figure 13, showing the pumps running at any time step for the best solution. As summarized in Table 3, such solution is obtained with a total number of switches equal to eight.

Finally, the evolution of the OF with the number of iterations is plotted in Figure 14. As for the Anytown WDN, because of the large number of starts (*ST* = 4), the OF rapidly decreases and achieves a value very close to the minimum cost in a small relative number of iterations. Also in this case, such a value has been achieved for an iteration rate around 30%.

## CONCLUSIONS

In this paper a HSMO algorithm was developed as an effective tool for pump scheduling optimization in WDNs. The algorithm improved pump scheduling using a RTTC instead of the classic binary representation. Moreover, unlike models available in the literature, hydraulic and physical constraints (water tank levels, pump performances, schedule representation) allowed the exclusion of *a priori* unfeasible solutions, thus reducing the space of solutions. An OF to minimize both the total daily energy costs and the daily number of switches for each pump was considered. Hydraulic and physical constraints were set to guarantee the operation of the network.

The effectiveness of the model was tested on two benchmark WDNs available in the literature. The comparison of results showed the capability of the proposed algorithm to achieve good solutions in low computational times. In both cases, the proposed model achieved better results (or very similar), i.e. with lower cost, in short computational times. Analysis also showed that for a greater number of starts, the algorithm achieves solutions very close to the minimum cost in a small (relative) number of iterations.

Application to more complex systems will demonstrate the suitability of the proposed approach for a wider set of operative conditions.

Proc. World Environ. Water Res. Congress