Optimal scheduling of variable speed pumps with mixed integer linear programming

Pump scheduling plays a crucial role in optimising water distribution network (WDN) operation and has been an important subject of research over the last decades. The purpose of pump scheduling is to find the sequence of pump ON/OFF statuses and pump speeds, in case of variable speed pumps (VSPs), that minimise the pumping costs. This is achieved via (a) pumping during time periods with lower electricity tariffs and utilising water storage, and (b) routing the flows through the network such that hydraulic losses are minimised. Optimal pump schedules depend on the electricity tariff profile, the demand profile, and the hydraulic characteristic of network components. Since up to 70% of total operating costs of WDNs are attributed to electricity consumption for pumping (Menke et al. 2015), pump scheduling is an economically important task. Past research and industrial case studies have shown that optimisation of pump schedules can lead to up to 10–20% reduction in pumping costs, i.e. up to 7–14% of the total operating costs of WDNs.

Nevertheless, finding optimal pump schedules is a difficult dynamic mixed integer nonlinear and non-convex optimisation problem (Bonvin et al. 2017). Nonlinearity and non-convexity arise from the nonlinearities imposed by the network equations that act as constraints. Integer variables used for selecting statuses of pumps and other active network components such as valves, as well as, in case of piece-wise approximations and relaxations, selection of active (linear) domains, make the problem mixed integer. Dynamics are introduced by storage tanks and require the optimisation problem to consider variables and constraints from all time steps, what substantially increases the problem size. Consequently, pump scheduling problems are inherently difficult to solve and, depending on the size of the model and on the optimisation method, may get stuck at local minima, fail to find a global optimum within the allocated time or at all, or fail to find a feasible solution altogether.

Optimisation methods, including those applied to pump scheduling, can be broadly classified into two categories: (a) deterministic mathematical programming, such as mixed integer nonlinear programming and mixed integer linear programming (MILP) and (b) stochastic evolutionary searches, such as genetic algorithms, particle swarm optimisation, and Ant Colony Optimisation (ACO). The former uses full or partial information about the model to guide the optimisation process, while the latter explores the objective space without reliance on the information about the model. Consequently, mathematical programming tends to be significantly faster at the cost of being more difficult to formulate, while evolutionary algorithms (EAs) tend to be slow but relatively easy to set up. On top of this, they are able to work in conjunction with an input–output model of any mathematical form. Slower convergence speeds are leveraged, to some degree, by some genetic algorithms (GAs)’ massive parallelisation capabilities, as demonstrated by Reed & Hadka (2014). In addition, many GAs support multiple objectives, while multi-objective extensions within mathematical programming frameworks are far less common. Last but not least, EAs, including GAs, are very good at exploring large decision spaces. Unfortunately, they cannot guarantee that the global optimal solution has been found nor provide bounds on global optimality (Menke et al. 2016) which is in contrast to (convex) mathematical programming which can provide such guarantees.

The literature on the subject of pump scheduling is extensive and voluminous. For a more in-depth study, the readers are referred to the most recent review papers on the topic by Mala-Jetmarova et al. (2017) and Wu et al. (2018). Here, we shall only mention a handful of publications on pump scheduling that pertain to MILP. Applications of MILP have recently become more popular, most likely as a result of the improvements in speed, reliability, scalability, and stability of MILP solvers, such as CPLEX (Cplex, IBM ILOG 2009), GUROBI (Gurobi Optimization, LLC 2023), MOSEK (MOSEK 2023), and SCIP (Bestuzheva et al. 2021).

The originally non-convex mixed integer nonlinear programming (MINLP) pump scheduling problems need to be numerically simplified to improve the speed and the robustness of finding globally optimal solutions. The literature distinguishes between (a) model simplification/reduction that is applied to the model before optimisation problem formulation and (b) simplification of the optimisation problem using various mathematical techniques from the operational research field. Model network simplifications were described by Anderson & Al-Jamal (1995), Deuerlein (2008), and Alzamora et al. (2014) and in the context of real-time pump scheduling, by Shamir & Salomons (2008). The most popular methods for problem simplification are (a) various relaxations and approximations of constraints, objective(s), and relaxations of decision variables types such that the original problem is turned into a convex nonlinear or linear problem, (b) problem decomposition into several easier to solve problems, and (c) relaxation of the optimality criterion e.g. Gleixner et al. (2012).

The issue of non-convexity was addressed by Fooladivanda & Taylor (2018) and Singh & Kekatos (2019) who turned the non-convex MINLP formulation into a convex MINLP formulation via different second-order cone relaxations. In addition, Fooladivanda & Taylor (2018) incorporated VSP and pressure reducing valves (PRVs) in their work. Bonvin et al. (2017) addressed the non-convexity in a less formal and more heuristic way that required restricted formulations that only supported network topologies without loops.

Problem decomposition into short- and long-term optimisations was investigated by Pulido-Calvo & Gutierrez-Estrada (2011). Similarly, Ulanicki et al. (2007) proposed a time decomposition method via a solution of a relaxed continuous problem to find optimum reservoir trajectories, followed by a solution of a mixed integer pump scheduling problem that tracks these trajectories. Lagrangian decomposition and Benders decomposition were successfully applied to MINLP problems by Ghaddar et al. (2015) and Naoum-Sawaya et al. (2015), respectively. Alternatively, decomposition can also be approached by dividing the network into smaller isolated subsystems.

Various mixed integer pump scheduling problem formulations involving approximations and relaxations were recently proposed by a number of authors. Vieira et al. (2020) used an iterative approach with a feedback loop from EPANET to iteratively limit the error from MILP with component relaxations. In a conceptually similar way, Liu et al. (2020) created a MILP formulation with component relaxations which are tightened by the solutions of a series of EPANET simulations. The authors claimed that their formulation with pipe characteristic relaxations over-performs one with piece-linear approximations due to reduction of the number of binary variables. Salomons & Housh (2020) tested different thresholds put on the numbers of binary variables in order to align the computational times of MILP pump scheduling schemes with the requirements of real-time pump scheduling applications. Bonvin et al. (2021) developed a relaxation of non-convex constraints of the original non-convex MINLP problem using polyhedral outer approximations and solved the relaxed convex problem with the branch and bound method for convex MINLPs. Their method also supports VSP. Most recently, Tasseff et al. (2022) developed tight polyhedral relaxations of the original MINLP, derived novel cuts using duality theory, added novel optimisation-based bound tightening and cut generation procedures, and implemented their method in an open-source Julia package (Tasseff et al. 2019). The authors addressed two main deficiencies of current state-of-the-art MILP solvers: the slow improvement of dual bounds and the difficulty in generating feasible primal solutions. The authors considered fixed-speed pumps only.

Although many advancements in pump scheduling using mixed integer programming (MIP) have been introduced recently, the literature on the subject is still fragmented with papers addressing some important aspects of the methodology while omitting others. In addition, treatment of crucial network components such as VSP is under-represented. Meanwhile, the findings reported by many authors suggest that the current state-of-the-art MILP solvers are able to find optimal pump schedules for networks with around 100 nodes Liu et al. (2020), i.e. of sufficient complexity to make them practical. This suggests that an automated method for solving pump scheduling problems with MIP could be of practical value to the community. The research is moving in the right direction, and recently, a comprehensive MIP framework for optimisation of water network operation problems has been published by Thomas & Sela (2024); however, although extensible, it does not yet feature VSP.

In this article, we describe our methodology for automatic network conversion into a MILP problem and subsequent solution using one of the available solvers. We provide a mathematical description of a MILP formulation of a WDN which offers a formal but gentle introduction to the subject. The formulation includes VSP, which is novel. The linearised VSP model is derived from the equations published in the thematic literature (Ulanicki et al. 2008) which were shown to integrate well with EPANET Janus & Ulanicki (2020). The validity of the approach and its reliability and robustness are tested on a small network with two parallel VSP and one tank. The methodology is accompanied with an open-source software package written in OCTAVE/MATLAB and Python for solving pump scheduling problems with MILP solvers (Janus & Ulanicki 2023) accessible on https://github.com/tomjanus/milp-scheduling. The software can be used to replicate the results of this study or incorporated into new studies on pump scheduling.

Water network is modelled with standard network equations describing: (1) headlosses in pipes, (2) flow continuity equations in nodes, including tanks, (3) pump characteristics describing pump head and pump power consumption vs. flow, and (4) head vs. volume relationships in tanks. More information about modelling of WDNs can be found in the study by Strafaci et al. (2007). The equations are briefly listed below since their understanding is needed for the understanding of the subsequent model approximation steps.

When approximating nonlinear functions with linear functions, the choice of linearisation technique is often left to the user. The possible choices are: (1) tangent line approximation using first, e.g. two, terms of Taylor’s expansion, (2) linearisation by substitution via introduction of additional variables and transformations, (3) piece-wise linearisation, and (4) convex hull approximation. The choice of method should be based on the specific characteristics of the nonlinear function and the context of the problem. Each method has its own limitations and applicability and, in the context of pump scheduling, will affect the solution accuracy and the size of the formulation. In the following sections, the choice of linearisation techniques was made by the authors, but the methodology is not limited to those choices and the readers are encouraged to try different techniques to fine-tune their problem formulations.

Equation (22) ‘binary linearises’ the flow with respect to segment selection. When is zero, is forced to become zero. Otherwise, is bound between and ⁠. Equation (23) relates the pipe flow variable q to the auxiliary flow variable ww. Equation (24) ascertains that only one segment in the piece-linear equation is active at a time. Finally, Equation (25) models linear pipe headloss.

The top surface of the polyhedron is described with equations of four Euclidean planes ⁠. For each plane ⁠, the point must lie within the triangular domain defined by the plane’s projection onto the space – see Figure 3 (right). The equations of each of the four planes are derived from the three vertices of the polyhedron that are contained within it. The domain for each plane is defined by three inequality constraints derived from the equations of the lines containing the three sides of its projection – see Section 5.2 of the Appendix. Direction of each inequality is indicated in Figure 3 with a short blue arrow. The projections of points ⁠, ⁠, ⁠, ⁠, and onto the space are denoted as ⁠, ⁠, ⁠, ⁠, and ⁠, respectively.

The formulation is composed of (1) objective function describing the total pumping cost, (2) a set of equality constraints modelling the linear and linearised network component equations plus additional auxiliary relationships, e.g. for setting bounds of the segment and segment selection variables in piece-linear approximations, (3) a set of inequality constraints describing the binary-linearised component equations plus additional constraints, e.g. symmetry breaking, (4) a set of lower bound and upper bound, also known as box constraints on selected decision variables, and (5) indices of binary decision variables.

The study involved two analyses. In the first analysis, pump schedules were optimised for a 24-h time horizon with a default set of inputs and parameters: tank elevation m, final tank level m = initial tank level ⁠, average demand L/s, and tank diameter m. The aim of the analysis is to study the behaviour of the MILP solver on our problem formulation and to verify the correctness of results. In the second analysis, a batch of pump schedule optimisation was performed for 81 combinations of network parameters and inputs. The ranges of parameters were as follows: tank elevations m, tank level differences ⁠, average demands L/s, and tank diameters m. The goal was to test the reliability and the robustness of the method over a range of operating points and to measure the calculation times required by the CPLEX solver to find optimal solutions. In both studies, the solver was set to terminate upon achieving the MIP optimality gap of 0.05. MIP optimality gap is a difference between the optimal objective value of a MIP or MILP problem and that of its linear programming relaxation. It is used to to adjust early termination criterion by setting a threshold on the minimum accepted quality of solution. For more details, please refer to Cplex, IBM ILOG (2009) or a documentation of any state-of-the-art mixed integer solver.

Flows in the selected network elements and heads in the selected network nodes are significantly altered by the new pump schedule as demonstrated in the middle and the right subplots vs. the left subplot in Figured 7 and 8, respectively. As the operation switched from the initial schedule with one pump working at a constant speed of 0.8 to the schedule in which the tank’s storage capacity is utilised in order to reduce pumping during high tariff periods, the flows in the elements between the pumps and the tank exhibit more variability. Consequently, the flows are higher, in absolute values, during the times when the tank is filling and emptying.

It is interesting to observe the discrepancies in heads and flows between MILP, which uses linear approximations of the network components, and the simulator which uses a complete (nonlinear) network model. We can see an offset between head at the demand node returned by the optimiser (middle) and the simulator (right), and a relatively higher pump outlet heads returned by the optimiser compared to the simulator. These differences stem from the inaccuracies introduced by the linear and piece-linear approximations These approximations can be tightened via introduction of a larger number of piece-linear segments or by iterative adjustment of the locations of break points. Both approaches normally come at a cost of increasing running times. Similarity in tank levels is preserved to a greater degree than similarity in the heads of non-storage nodes. This is due to the fact that tank levels change in response to changes in flows, which are affected by approximations to a lesser degree than pressures are. While flows are forced inputs in demand driven simulations, heads are the outputs and thus are determined by the model and hence subject to its inaccuracies.

We ran 81 pump schedule optimisations at different sets of inputs and produced 80 optimal solutions within the MIP optimality gap of 0.05. One solution was infeasible due to conflicting constraints imposed by the inputs. All calculations were performed on a laptop equipped with Intel(R) Core(TM) i5-6300U CPU, 2.40 GHz processor and 16 GB of RAM, and IBM(R) ILOG(R) CPLEX(R) v. 22.1.1.0 MILP solver. Seventy-five percent of calculations returned (globally) optimal solution within of the median time of 1.0 s, while none of the calculations took longer than 1.5 s. A total of 96.6% of solutions returned results in which the mean average error between the reservoir level time-series from the simulator and the optimiser was between 0.1 and 0.3 m – a rather robust outcome.

The results demonstrate that explicit implementation of VSPs in pump scheduling with MILP via piece-linear approximation of the pump characteristic is a viable approach, and it can be introduced into the currently developed pump scheduling optimisation frameworks such as MILPNet of Thomas & Sela (2024) or WaterModels (Tasseff et al. 2019). Moreover, this approach is consistent with past thematic literature (Ulanicki et al. 2008) and uses the pump model that was introduced into EPANET dev-battery branch (Janus & Ulanicki 2020). Hence, it can be included in frameworks which automatically convert EPANET models to mathematical programming formulations, e.g. Thomas & Sela (2024). Contrary to Thomas & Sela (2024), we construct the optimisation problem in Equation (1) manually without the use of an AML. This lower level implementation may lead to faster conversion from nonlinear network models to their mixed integer linear representations, which could be beneficial in real-time applications. A disadvantage of this approach stems from the additional time that needs to be spent on extending the software to support new solvers and network components. It is also not suited for writing extensible optimisation frameworks where new features could be included as add-ons and plug-ins. Treatment of VSPs in this article is different to other published approaches, e.g. Bonvin et al. (2021). While we explicitly linearise the pump characteristic, Bonvin et al. (2021) consider VSPs as nonlinear components and requires solving a nonlinear programme on top of a (relaxed) MILP formulation. Such a solution might indeed scale better to larger networks, as mentioned by the authors, albeit at a cost of increased complexity, which could be a barrier to adoption by end users. Furthermore, this two-layer approach could be more difficult to implement in multi-threaded applications, while multi-threading is already supported off the shelf by some of the MILP solvers such as CPLEX, GUROBI, and MOSEK and could be implemented at ease. With many improvements and automated optimisations added to MILP solvers, formulations that use linear and piece-linear approximations of network components, e.g. Shamir & Salomons (2008), and Thomas & Sela (2024) and this work, or convex relaxations, e.g. Vieira et al. (2020), Liu et al. (2020), and Tasseff et al. (2022) sound practically viable while being conceptually simple. The limitations imposed by network size can be partially rectified through network reduction, e.g. using a recently developed package by Thomas & Sela (2023).

In light of the extensive recent literature on MILP approaches in pump scheduling, we opted to concentrate on incorporating VSPs, an area traditionally given less focus, instead of contributing yet another option for whole-system optimisation. Minimising the number of other network components was a deliberate choice to enhance clarity. Instead, we strived to offer a detailed and transparent mathematical description of the fundamental parts of the network. Hence, this article also serves as an introduction to pump scheduling with MILP for water engineers and scientists unfamiliar with the method. It is complemented by closely matched source code, which aides in understanding the manuscript. The code facilitates automatic translation of network models to MILP formulations and provides variable mapping between the linear programme and the network model. This allows MILP solvers to be integrated into automated in-loop simulation–optimisation solutions.

The study demonstrates that MILP with linear and piece-linear approximations of the objective and of the model components is a valid method for finding globally optimal pump schedules in networks containing VSPs. Although the method was tested on a very small network, the average calculation time of approximately 1 s is a promising result, indicating that the same approach can be adapted to solving more complex networks. Similar studies showed that this indeed is the case. The method proved to be robust and able to arrive at global optimal solution in a similar time for a range of operating points. The article is accompanied by an open-source toolbox written in OCTAVE/MATLAB and Python (Janus & Ulanicki 2023) which can be used as an additional source of information supporting the manuscript, for further experimentation with the method and for reproducing the results of this study.

While the existing literature extensively covers technical aspects of formulating and solving pump scheduling problems with MILP, future efforts should prioritise merging these findings into user-friendly, flexible open-source tools for practitioners. Recent developments, such as MILPNet by Thomas & Sela (2024) in Python, and WaterModels by Tasseff et al. (2019) in Julia, alongside our manuscript’s MATLAB code, exemplify this trend. Such tools can facilitate further experimentation, with potential directions including expanding component libraries to include multitude of parameterised approximation and relaxation approaches for each component, enhancing speed via decomposition techniques, e.g. Lagrangian or Benders decomposition, and leveraging parallel execution capabilities of solvers like CPLEX, GUROBI, and MOSEK. These tools could support offline as well as real-time pump optimisation. Another application could be integrated design-operation optimisation with a hierarchical structure in which the inner pump schedule optimisation loop is solved using MILP while higher level decisions, e.g. long-term policies, design options, are calculated with EAs.

All relevant data are available from an online repositoryor repositories (please ensure the DOI/URL has been provided as a submission item).

The authors declare there is no conflict.