Optimal solving of the pump scheduling problem by using a Harmony Search optimization algorithm

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.

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 2: HM is compiled with many randomly generated solution vectors as the HMS.

  

  (1)
• Step 3: A new harmony vector, x′= (x′₁, x′₂, … , x′_N) is generated, according to three rules: Memory Considerations, Pitch Adjustments and Randomization. The value of the first decision variable (x′_i) can be chosen from any value in the specified HM range (x_i¹ ~ x_i^HMS) and a probability of choosing a totally random value for the decision variable x′_i is also considered:

  

  (2)

  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.

  Values of the other decision variables (x′₂, … , x′_N) 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.

• Step 4: From the calculation of the OF, if the new harmony vector x′= (x′₁, x′₂, … , x′_N) 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.

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

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_i) to the active mode (during t′_i). As a consequence, the sum of all time intervals for each pump has to be smaller than or equal to the scheduling period T.

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

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.

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.

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.

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

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

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.