Exploring Optimal Pump Scheduling in Water Distribution Networks with Branch and Bound Methods

Water utilities can achieve significant savings in operating costs by optimising pump scheduling to improve efficiency and shift electricity consumption to low-tariff periods. Due to the complexity of the optimal scheduling problem, heuristic methods that cannot guarantee global optimality are often applied. This paper investigates formulations of the pump scheduling problem solved using a branch and bound method. Piecewise linear component approximations outperform non-linear approximations within application driven accuracy bounds and demand uncertainties. It is shown that the reduction of symmetry through the grouping of pumps significantly reduces the computational effort, whereas loops in the network have the opposite effect. The computational effort of including convex, non-linear pump operating, and maintenance cost functions is investigated. Using case studies, it is shown that linear and fixed-cost functions can be used to find schedules which, when simulated in a full hydraulic simulation, have performances that are within the solver optimality gap and the uncertainty of demand forecasts.


Introduction
In the US up to 4 % of all energy demand is consumed by water distribution and treatment works (Pasha and Lansey 2014). In particular 70 % of the life cycle cost of a pump system can be attributed to its electricity consumption (Nault and Papa 2015). Optimal pump scheduling has been shown to reduce the energy cost of a water distribution system (WDS) by 10-20 %, this by shifting consumption to time periods with lower electricity cost or improving operational efficiency (Boulos et al. 2001).
Pump scheduling for fixed energy tariffs has been studied extensively, detailed reviews can be found in D' Ambrosio et al. (2015) and Singh (2014). However, as energy grids are changing with the introduction of renewables, different operating strategies may be favourable compared to optimising for a fixed energy tariff. This paper outlines in detail a selection of optimisation approaches for such analyses. Then the performance of a range of mathematical optimisation methods suitable for different analyses is explored. Following an extensive review of published literature on optimising pump schedules, a mixed integer problem formulation solved with a branch and bound method is selected as optimisation method most suited to achieve the objectives outlined above. The performance of different model approximations for the objective function, which describes the operating cost, and the constraints, describing the hydraulic model are investigated. The approximations investigated include quadratic equations for modelling pumps and pipes, as well as convex and piecewise linear approximations. The improvements in computational speed for different approximations of system components are assessed with respect to the resulting loss in accuracy of the hydraulic solution. Results of a preliminary analysis of this method by the same authors (Menke et al. 2015) are further explored on a larger benchmark network. Here the changes in computational performance are investigated when taking into account maintenance constraints on pumps, grouping pumps to reduce the symmetry of the model and changing the problem size by varying the number of time steps and the size of the network. This analysis is performed using a benchmark network and skeletonised network based on an operational network model supplying roughly a million customers in the UK.

Literature Review
Optimal pump schedules for a water distribution system (WDS) are derived using either heuristic methods or mathematical optimisation. In the former methods, heuristics are used to generate schedules and evolutionary methods guide the search for optimality. Mathematical optimisation solves the problem by using information of the objective function to guide the search and a set of constraints to limit the search to feasible states of the real system.

Review of Heuristic Methods
The complexity and non-linearity of the equations describing the hydraulic system has led to the application of heuristic methods, which solve the hydraulic problem separately from the optimisation problem, in a dedicated hydraulic simulation. Genetic algorithms (GA) have been applied to optimise pump scheduling (Mackle et al. 1995), and have been combined with improvements for local searches (Van Zyl et al. 2004;Reis et al. 2006;Siew et al. 2016). The large number of solutions generated by a GA lends itself to the application of multi-objective optimisation (Savic et al. 1997;Siew et al. 2014). To reduce the computational effort of each iteration the hydraulic simulation can be replaced with an artificial neural network (Jamieson et al. 2007). Warm starts for GAs can improve the computational effort and have been applied to pump scheduling (Pasha and Lansey 2014). Simulated annealing, a heuristic method inspired by the cooling processes of metals, has also applied to pump scheduling (Goldman and Mays 1999;Pedamallu and Ozdamar 2008;Teegavarapu and Simonovic 2002;Samora et al. 2016). Ant colony optimisation, a nature inspired optimisation technique, has also been applied optimise pump scheduling (López-Ibáñez et al. 2008;Afshar et al. 2015). Other evolutionary methods applied to pump schedule optimisation in WDS include harmony search optimisation (Kougias and Theodossiou 2013) and particle swarm optimisation (Wegley et al. 2000;Ostadrahimi et al. 2012). Heuristic methods can provide optimised pump schedules for hydraulic systems, but cannot provide bounds on the global optimality of the solution.

Review of Mathematical Optimisation Approaches
Mathematical optimisation approaches that have been applied to pump scheduling include dynamic programming (Dreizin 1970), model predictive control (Sun et al. 2015), iterative methods (Price and Ostfeld 2013), benders decomposition (Cai et al. 2001) and branch and bound methods (Costa et al. 2016).
Dynamic Programming (DP) has been applied to the optimisation of pump schedules but for larger systems it suffers from the "curse of dimensionality" (Powell 2005). To overcome this, the WDS can be decomposed into sub-networks or pump switching decisions can be aggregated (Joalland and Cohen 1980;Zessler and Shamir 1989).
Model Predictive Control (MPC) is an iterative control method applied in receding horizon fashion, widely used in the control of plant operations (Camacho and Alba 2013). MPC has been applied to the optimisation of WDS operations (Ocampo-Martinez et al. 2013;Fiorelli et al. 2013). Although it does not guarantee global optimality of operations, it has yielded good control strategies in practice (Nikolaou 2001), acting against uncertainties in demand at each schedule update. It has also been applied as a control technique to mitigate maximum demand charges when operating a WDS (van Staden et al. 2011). Other iterative methods, decompose the scheduling problem into sub-problems that can be solved by a range of mathematical optimisation methods. This includes the application of linear programming (Bunn and Helms 1999;Price and Ostfeld 2013), or graph search methods (Price and Ostfeld 2016).
The pump scheduling problem can be posed as a mixed integer non-linear problem (MINLP) and solved using a branch and bound algorithm (Gleixner et al. 2012;Burgschweiger et al. 2008). For a special network structure the components can be approximated with convex relaxations, yielding a more tractable convex MINLP (Bonvin et al. 2016).
Other methods to improve the tractability are a Lagrangian decomposition of the problem (Ghaddar et al. 2014), or generalized benders decomposition (Verleye and Aghezzaf 2016). The computational effort can be reduced, however these decompositions can lead to a loss of the global guarantee of the optimality.  (2004). It is comprised of two storage tanks and a demand node supplied by two pump stations with a total of three pumps

Problem Formulation
An example of a WDS, used to illustrate pump scheduling, is shown in Fig. 1. The system transports water from a source reservoir, J1, to demand nodes, for example J5 in this network. The required pressure head is provided by pumps such as main 1 or booster and stored in elevated tanks such as J6. Pipes labelled P1, P2 . . . P8 connect the components.
The optimal pump scheduling problem can be formulated as the mathematical optimisation problem: Minimise : Pumping cost subject to : Energy balance, where the energy and mass balance are the governing physical constraints for the WDS. As outlined the optimisation method presented here aims to enable analysis of energy usage strategies of WDS and this requires the guarantee of a bound on the global optimality gap. A convex branch and bound problem can provide such a guarantee (Garfinkel and Nemhauser 1972), a general outline of a branch and bound algorithm applied to pump scheduling is given by Algorithm 1.

Algorithm 1 Branch and bound
Given a MILP of the form: Minimise: 1 2 subject to: where 1 , 2 , 1 , 2 , define the problem, is set of continuous variable and is a set of binary variables. The outline below gives a general procedure, however to reduce computational effort, solvers may use a different order of steps or additional steps, to find solutions and determine branching directions (IBM 2009;Achterberg 2009).
Step 0: Select a binary variable and partition the problem into two branches and setting the selected variable to its extreme values.
Step 1: Set an upper bound to or find a feasible solution for (2) through a heuristic and set . Set a lower bound by evaluating (2) with a relaxed integrality requirement (by letting [0 1]). Set .
Step 2: Evaluate with each branch by relaxing the integrality requirement for the other binary variables and solve the relaxed problem.
Step 3: If a solution from Step 1 is a feasible solution of (2) (respects the integrality constraint) set it to .
Step 4: Set the infimum (greatest lower bound) of the solutions to Step 1 as if it is greater than the current .
Step 5: Discard any branches where of the branch is greater than of another branch.
Step 6: If the optimality gap is lower than the desired level, or the computational time limit has elapsed then stop else go back to Step 1.

The solution forms the lower bound for each branch
Finding more feasible solutions through heuristics can improve the convergence performance.

Objective Function
For a fixed speed pump, the pump's ON-OFF state can be described by an integer variable. The state of a pump, i p at time step j ∈ [0, N] it is given by T i p ,j ∈ {0, 1}. The cost of operating a fixed speed pump can be expressed as: where P i p ,j is the cost of energy in having pump i p ON at time j . This approximation assumes a fixed power consumption. While the power consumption may be relativity constant in the operating range of the pump, the power consumption of a pump can vary with flow. To account for this the power consumption can be approximated: where Eq. 4 is a linear approximation with coefficients c e i p ,1 and m e i p ,1 and a quadratic approximation is given by: with coefficients a e i p ,2 , b e i p ,2 and c e i p ,2 . For fixed speed pumps, switching the pump state causes fatigue loading on the water distribution system and the pump it self. To account for these, pump switches are often charged with a proxy cost (Lansey and Awumah 1994;Savic et al. 1997).
A penalty function that penalises ON-to-OFF and OFF-to-ON switches equally is added to the objective function.
where P s is a switching penalty and the equality in Eq. 6 holds because T ∈ {0, 1}.
In summary, the pump cost are considered either as fixed during the time interval the pump operates as described by Eq. 3, or proportional to the flow rate with either a linear dependence on the flow through the pump as given by Eq. 4 or a quadratic dependence on the flow rate as described by Eq. 5. In addition, maintenance cost can be accounted for by adding costs for switching as in Eq. 6.

Pump Approximations
When a pump is ON, the characteristic curve of the pump describes the relationship between the head difference across the pump and the flow rate. When the pump is OFF, the flow through the pump is zero. For pump i p that connects nodes J 1 and J 2 with flow q i p , the pressure difference across the pump at a time step can be represented with a polynomial function: This results in quadratic equality constraints that make problem (1) a non-convex mixedinteger non-linear program (MINLP). However, convex relaxation can be formed by neglecting the concave constraints. These constraints are implemented using the big-M method detailed in Section 3.3 and Gleixner et al. (2012). A set of linear constraints describing a convex set can also be considered to approximate the characteristic curve. The constraints are enforced if T i p = 1, and are described by: where m i p ,1 . . . m i p ,N con and c i p . . . c i p ,N con are the linear coefficients and N con is the number of these constraints and Δh ub is an upper bound on the head difference across the pump. These constraints are also implemented using the big-M method.
To further simplify the pump description, a formulation prescribing a fixed head difference is implemented assuming that pumps operate near their best efficiency point (BEP): where h f ixed is the head difference at the BEP and is enforced if the pump is on, q min and q max are limits close to the BEP flow rate. Such a formulation was used by Gleixner et al. (2012).
Often multiple pumps are installed in parallel or in series at pumping stations, such as pumps main 1 and main 2 in the Van Zyl network shown in Fig. 1. These pump configurations can be modelled as individual pumps or as a single pump represented with a composite pump performance curve (Larock and Jeppson 2000).
For pump station i p with a set of pumps with a total of k p composite pump curves, the different composite curves of the pump station can, for a given time step, be approximated by quadratic equality constraints as given by: If a hierarchy between the pumps can be enforced, the number of pump combinations can be reduced to the number of pumps in the station, i.e. k p = n p . The ON-OFF states of the pumps are linked explicitly through (Gleixner et al. 2012): Representation (12) can group pumps that are not immediately connected. Both pump grouping approaches were applied to the Van Zyl network and other test networks.
In summary the characteristic curve of the pumps in the network are approximated using either quadratic equality constraints as given by Eq. 7, a set of linear constraints forming a convex set as given by Eq. 8 or by simply enforcing a head difference and flow limits as shown in Eq. 9. The pumps in pump stations are either left as individual pumps or grouped by modelling their composite curves as given by Eqs. 10 and 11 or by explicitly enforcing a hierarchy through Eq. 12.

Pipe Approximations
Frictional head loss of pipe flows is often modelled using the Hazen-Williams or Darcy-Weisbach formulae. These non-linear empirical equations can be approximated to reduce the computational effort and formulate a convex problem. Since the Hazen-Williams and Darcy-Weisbach equations both have a power of ≈ 2 with respect to the flow rate, a quadratic polynomial approximation can provide a close fit. When considering the example from Fig. 1, a quadratic constraint for a given time step for pipe P 2 , connecting nodes J 3 & J 4 is given as: For the MIP solver, Eq. 13 is implemented using the big-M constraints: where, h j is the head at node j and a P 2 , b P 2 & c P 2 are the fitted coefficients for the pipe, c P 2 represents the elevation difference. λ P 2 is a binary switching parameter with λ P 2 = 1 ⇒ q P 2 ≥ 0 and λ P 2 = 0 ⇒ q P 2 ≤ 0 and M is a sufficiently large constant. This results in quadratic equality constraints and renders problem (1) a non-convex MINLP. However, as for Eq. 7, by neglecting the concave constraints (lines 2 and 3 in Eq. 14), a convex relaxation of the problem can be found. In order to replace quadratic constraints and generate a convex problem, piecewise linear approximations are used to provide an approximation of the head loss in a pipe. For the pipe P 2 , the piece-wise linear constraints are implemented using big-M constraints: where λ P 2, (1,2...N piece ) ∈ {0, 1} is a switch for each corresponding linear section. A formulation with N pieces introduces N pieces of integer variables per pipe at each time step. Consequently, approximations with more pieces that more closely follow the head loss curve may scale badly for larger network models.
In summary, the pipes in the system have head loss modelled with either quadratic equality constraints as shown in Eq. 14 or a set of piecewise linear constraints as given by Eq. 16. By increasing the number of linear pieces given by N con , it is possible to more closely approximate the head loss curve; in this paper either N con = 3 or N con = 7 are used.

Mass balance at network nodes
Since steady-state approximations of the hydraulic conditions are used, the mass flow is equal to the volume flow. For a network node joining links P 1 , P 2 , · · · P n , the mass flow must balance at each time step j as: where d j is the demand at the node at time step j .A minimum hydraulic head can also be enforced at nodes if desired. For a tank, such as J 4 in Fig. 1, with flows q P 2 and q P 3 the mass balance for time steps j = 1 . . . N − 1 is given by: where the horizontal surface area of the tank is given by A J 4 and the tank is assumed to have vertical walls. The constraints on the final time step level of the tank are described by: where δ V defines an acceptable limit of volumetric difference between initial and final fill levels. By defining the final state in terms of the initial state, they aren't specified by the input data as this would limit the feasible search space and would lead to a sub-optimal final solution. This is a common technique to set the initial and final tank levels (Price and Ostfeld 2013).

Problems Considered and Investigation Procedure
The component approximations and the associated objective functions and constraint equations are used to construct a number of MIP problems as summarised in Table 1. By  (14), (17) For each optimisation problem defined the constraints apply to all relevant components: (7, 8 or 9) apply to all pumps, (16 or 14)apply to all pipes and (17) to all nodes and (18) to all tanks comparing the computational effort required and the level of optimality of the solutions among the different problem formulations, the individual effects of each modelling decision can be assessed. The investigation presented in this paper is done in three stages. At first, the effect of symmetry is investigated by considering network formulations with and without grouping constraints on pump stations. Secondly, the optimal schedules found are simulated using a fully nonlinear hydraulic simulation. This simulation provides data on the "true operational cost" C and hydraulic conditions of the WDS when operated using the schedules from optimisation. This data is then compared to the computed operating cost C and hydraulic conditions modelled in the optimisation. The mean error is given by: where X f is the set of feasible solutions computed by the MIP solver using one of the formulations and |X f | is the size of the set. Finally, the errors in the estimation of the hydraulic operation by the optimisation and the computational effort of the optimisation can be weighted up to compare trade-offs and guide the selection of suitable methods for further analysis. The hydraulic solution computed in the optimisation (u) is compared to a high-accuracy hydraulic simulation solution (u ) produced via a null space method  using the pump control profiles (x) from the optimisation as inputs. The mean error in hydraulic analysis is given by: where E hyd is the mean error in flow rates for all flows in a given network and over all time steps.

Networks Analysed
When analysing the performance of pump scheduling optimisation methods, it is important to consider different networks since the suitability of certain methods and algorithms vary with the network topology, controlled components and operational conditions (Jolly et al. 2014). The Van Zyl network shown in Fig. 1 is a test network with one loop. The performance of the method presented here on a larger network is analysed using the network shown in Fig. 2a. To analyse the impact of loops in a WDS on the computational effort, three artificial networks with zero (Fig. 2b), one (Fig. 2c) and two loops (Fig. 2d) are investigated. (a)

Results and Discussion
The computational effort in solving MIPs cannot always be predicted accurately; however, the results show that certain aspects of the problem formulation should be considered carefully to improve performance.

Problem Symmetry
In the Branch and Bound algorithm, grouping reduces the number of possible search tree branches with equal objective function values. This is known as symmetry and it is considered good practice to remove it from the model (Sherali and Smith 2001). Grouping limits the combinatoric possibilities for pump settings from: for the ungrouped case to: for the grouped pumps, where n H is the number of pump stations and k p,j the number of pump combinations in station j . Removing symmetry has a very pronounced effect on the computational performance of the solution. The Van Zyl Network with twelve time steps and approximation number 7 is solved in 373 s when the pumps are not grouped, but in less than 4 s when the pumps in the main pump station are grouped. Similar performance gains were observed for different networks and time step sizes. As removing symmetry does not affect the optimality of the solution, grouping is applied to all simulations.

Optimisation Constraints
For a given network, the size of the MIP problem (i.e. the number of variables and constraints) is dictated by the number of time steps and component approximations chosen. To investigate the influence of the overall size of the problem on the computational effort the analysis is performed on a range of simulation time steps. These are doubled three times from 6 to 48 representing four hour to half hour time-steps for the daily pump schedules. The observed increase in computational effort as number of time steps are increased is exponential for all approximation methods and networks. However, the number of links and nodes of the network graph, does not have such a distinct effect. The number of integer variables and the number of possible configurations that are explored in the branch and bound algorithm increase, but as linked components have linked states, the the size of the network in terms of pipes and nodes has less influence on the computational effort. The computational effort measured in CPU time taken to get to a solution with a 5 % bound on gap to global optimality, or in the optimality bound after 600 s, are shown in Fig. 3 along with the operating cost of the schedules found using a full hydraulic simulation. The formulations 1a, 2 and 3 are computationally the most expensive. These formulations are intractable for all but the smallest number of time steps and network sizes. The piecewise linear pipe approximations are significantly faster to evaluate than the quadratic equality constraints. This brings a reduction in computational effort that is amplified through the use of a mixed integer solver tailored for linearly constrained MILP/MIQP problems. However,  Table 1 as the differences in computational effort between approximations 1a, 2 and 3 shows, the problems with only some piecewise approximations are faster to evaluate with the same solver. WA linearly constrained, quadratic cost problem is formed by formulations 4-6, 7a and 7d. A linearly constrained, linear cost problem is formed by formulations 7b, 7c and 8. The number of pieces in the piecewise approximations Eq. 16 (i.e. how closely the approximations fit the nonlinear pipe headloss equations) has a significant impact on the computational effort, much more so than the number of piece-wise linear constraints that define the convex sets for the pump models Eq. 8. This is expected because the number of pipes np in a WDS is larger than the number of pumps n p , even by many orders of magnitude for nonskeletonised models, and so impacts the number of integer variables for the MIP through a larger proportional constant (i.e. np × N pieces for pipes versus n p × N con for pumps).
As the constraints represent the hydraulic model of the network, the error in network states from the schedule optimisation compared to a hydraulic simulation of the schedules is governed by the choice of component approximation. The non-linear approximations 1-3 provide a close fit to the hydraulic properties of the network components. The small geometric mean errors of 1.74 %, 3.22 %, 3.99 % for approximations 1a ,2 and 3, respectively, demonstrate this as well as the small spread in errors observed with the largest errors in flow rate of 3.2 %, 6.0 % and 6.2 %, respectively. For the piecewise linear approximations the mean error ranges from 2.8 to 4 % with approximation 4 having the smallest error and approximation 7a the largest. The maximum error observed is 7 %. This indicates that the number of pieces in the pipe piecewise linear approximations dictates the size of the error in the hydraulic approximation. The error in demand estimate for 15 minute time steps is approximately 5-10 % (Herrera et al. 2010;Donkor et al. 2014). These results lead to the conclusion that a piecewise approximation with only a few pieces will often have sufficient detail to provide a schedule with a hydraulic error of size similar to the error in the demand forecast while significantly reducing the computational effort.
The simple fixed head representation for the pump, approximation 8, gave results with the smallest solve time. The results however have the largest hydraulic error, which may lead to schedules that are not repeatable in a real network. For example some flow rates through the pumps returned by the MIP solver for the 48 time steps on the Van Zyl network were higher than those achievable in hydraulic simulation, leading to an infeasible solution.

Objective Function
The objective function formulations are investigated by comparing formulations 1a and 1b as well as 7a to 7d. The accuracy is given by (20), comparing the objective value of the optimisation to the operating cost found by a full hydraulic simulation of the schedule. The approximation 7b has the largest average error in estimated objective function E OF = 7.1 %. Together with approximation 7c it has the lowest computational effort, reaching 10 % optimality bound after 600s compared to 13.7 % for approximation 7a and 55.5 % for approximation 7d. However, the problem formulated with approximation 7d resulted in E OF = 2.8 %, only approximation 1b achieved a similarly low error. Approximations 7a and 7c both have an error of ∼ 4 %. This shows that the linear objective functions are slightly faster to solve, while the quadratic objective function can negatively impact the solution speed, especially when a more detailed representation of the power consumption is required. With all constraints kept the same, the results do not conclusively relate the computational effort with different objective functions. When compared to minimizing the linear energy cost function f 1 (·) together with the quadratic switching penalties f 4 (·), minimizing f 1 (·) only results in a fractional increase in the computational effort or a marginal increase in the optimality gap of computed objective function when a 5 % gap cannot achieved within the set maximum solve time of 600s. However, as the error caused by approximating the power consumption curve of the pump is dependent on the pumps in the network and the shape of their power curves, these can vary significantly for different pumps. If the quadratic coefficient of f 3 (·) is negative, the objective function is no longer positive semi-definite or convex and may make the MIP only solvable to local optimality. The quadratic term for the switching constraints however, not considered in approximation 7b does not have a negative effect on performance. For the network models used the errors are small, the simplest pump cost approximation with a consideration for pump switching is used for all further analysis as outlined in Table 1. Bunn and Reynolds (2009) showed that in practice savings have been made through optimising the pump operations near the best efficiency point (BEP) while other applications such as Gleixner et al. (2012) have used less detailed representations with good success.

Loops in the WDS
Loops in the network have a significant impact on the computational effort in solving the scheduling problem. Despite representing an only marginally larger problem than the loop 0 network, the schedules for loop 1 and loop 2 networks are computed considerably slower with all formulations. For a simulation with 24 time steps the computational effort for the problem described by approximation 7 increases from 1 s for loop 0 to 9 s to 43 s for loop 2 . Similar changes in computational effort are observed for other time step sizes or approximations. The negative computational performance impact of loops could also be a contributing factor, as to why the Van Zyl network is only solved marginally faster than considerably larger Networks as highlighted by Fig. 3. Through network simplification methods, it is possible to remove loops from network models and replacing them with equivalent pipes. Such preliminary work could yield significant performance improvements for the optimisation methods presented here (Simpson et al. 2012).

Further Work
The motivation for this work was to enable further research on energy usage strategies and the need for a bound on the global optimality of solutions for comparisons of energy usage strategies. An example of such comparisons can evaluate the potential of WDS to provide reserve energy to a grid (Menke et al. 2016a). Further modifications could enable analysis of the scheduling problem as a two stage stochastic problem. Energy recovery from WDS with pumps as turbines or large pressure drops in networks also offer potential for optimisation and comparison of operation strategies (Kougias et al. 2014;Carravetta et al. 2012). While energy consumption of WDS are of concern, measuring the performance of a WDS may not be so straight forward and a range of performance indicators (Kanakoudis et al. 2011) could be used to optimise the operations for different priorities, possibly with a multi-objective optimisation approach (Pecci et al. 2016).
The problem formulation presented here can be expanded, as it only addresses fixed speed pumps, which make up the majority of installed pumps for utilities. Mathematical formulations and optimization for variable speed pumps, which are more popular in new pump installations, can be added to provide a wider scope of analysis (Menke et al. 2016b). Any loops in the network significantly reduces the computational performance of the formulations presented here. However, since loops are a common feature of networks, there are two possible approaches for further work. One alternative is to design tailored optimization solvers that take advantage of structure in networks with loops, for example through further constraints or relaxations. Another approach can consider methods for the reduction of networks to remove loops while generating functionally equivalent networks that are easier to evaluate for the MIP solvers. An example of such a reduction is described in Deuerlein (2008), and the performance of optimal pump schedules obtained from such a system could be the focus of further work.

Conclusions
In practice, for a piecewise linear approximation of pipe head loss formulae, the minimum number of pieces that generate an acceptable level of error in the hydraulics should be used. If the pumps' power consumption characteristics have a significant gradient with respect to flow rate, linear or quadratic approximations of the power consumption may be required. In optimizing schedules over 24 hours of operation, the computed operational cost was shown to improve by 5-10 % when using smaller time steps although the computational time for the MIP increases exponentially. Finally, the network model supplied to the optimisation should be simplified as much as possible, removing loops and symmetry where possible.
Using the method presented here optimised pump schedules with bounded optimality gaps can be generated. These can be used for further analysis of the energy consumption of the WDS. The application of such analysis has already been demonstrated and will provide a basis for further analysis of the operation of WDS in times of large changes in the energy supply.