A Mixed‐Integer Linear Programming Framework for Optimization of Water Network Operations Problems

Water distribution systems (WDSs) are critical infrastructure used to convey water from sources to consumers. The mathematical framework governing the distribution of flows and heads in extended period simulations of WDSs lends itself to application in a wide range of optimization problems. Applying the classical mixed integer linear programming (MILP) approach to model WDSs hydraulics within an optimization framework can contribute to higher solution accuracy with lower computational effort. However, adapting WDSs models to conform to a MILP formulation has proven challenging because of the intrinsic non‐linearity of system hydraulics and the complexity associated with modeling hydraulic devices that influence the state of the WDS. This paper introduces MILPNet, an adjustable framework for WDSs that can be used to build and solve an extensive array of MILP optimization problems. MILPNet includes constraints that represent the mass balance and energy conservation equations, hydraulic devices, control rules, and status checks. To conform to MILP structure, MILPNet employs piece‐wise linear approximation and integer programming. MILPNet was implemented and tested using Gurobi Python API. Modeling accuracy was shown to be comparable to EPANET, a public domain software for hydraulic modeling, and sensitivity analyses were conducted to examine the impacts of the modeling assumptions on the performance of MILPNet. Additionally, application of the framework was demonstrated using pump scheduling optimization examples in single and rolling horizon scenarios. Our results show that MILPNet can facilitate the construction and solution of optimization problems for a range of applications in WDSs operations.


Introduction
Water distribution systems (WDSs) can be large, complex systems containing thousands of components.The simplest WDSs contain sources of water (such as reservoirs), consumers (or junctions), and pipes to convey water from the sources to the consumers by gravity.The junction pressure heads and pipe flow rates of such gravity-fed systems can be computed by solving the mass balance at the nodes and energy conservation along pipes, provided that junction demands and pipe characteristics are known in advance (Boulos et al., 2006).However, large urban WDSs can be vast and usually require the introduction of additional control elements to ensure that all demands and pressures are met: tanks can allow for water storage in times of low demand or high supply, pumps boost pressure head, and valves can limit pressure head and flow.Tanks, pumps, and valves add flexibility to the operation of a system, but at the same time also introduce new levels of complexity.
Modeling WDSs to predict and control system behavior is a challenging task.First, the energy conservation equations that govern frictional head loss in pipes, which depends on the flow rate, diameter, length, and roughness, are non-linear (Boulos et al., 2006).The equations relating pump head gain to flow rate are also non-linear.These non-linearities greatly complicate the feasible space for solutions and make it challenging to determine the hydraulic state of the system.Second, pump and valve behavior is typically determined by their operating states.For example, pumps can assume either open or closed status, which must be modeled in different ways.If the pump is open, head is added to the downstream node according to the pump curve, and if the pump is closed, the upstream and downstream nodes are disconnected.Valves can inherit three different statuses: active, open, and closed.These discrete states, combined with the fact that flow through pipes can often be in either direction, make the problem non-convex with continuous and discrete variables (D'Ambrosio et al., 2015;Cassiolato et al., 2021).Third, boundary conditions impact the operation of several hydraulic elements, for example, tanks have finite storage capacities, and pumps and valves do not allow back flow.Finally, the operation of control elements in real WDSs is highly influenced by economic and supply constraints.For example, pump operation is typically desirable only at times of low electricity prices while minimizing the number of pump switches to reduce maintenance cost and extend asset life (Housh & Salomons, 2019).The complexity associated with operating Abstract Water distribution systems (WDSs) are critical infrastructure used to convey water from sources to consumers.The mathematical framework governing the distribution of flows and heads in extended period simulations of WDSs lends itself to application in a wide range of optimization problems.Applying the classical mixed integer linear programming (MILP) approach to model WDSs hydraulics within an optimization framework can contribute to higher solution accuracy with lower computational effort.However, adapting WDSs models to conform to a MILP formulation has proven challenging because of the intrinsic non-linearity of system hydraulics and the complexity associated with modeling hydraulic devices that influence the state of the WDS.This paper introduces MILPNet, an adjustable framework for WDSs that can be used to build and solve an extensive array of MILP optimization problems.MILPNet includes constraints that represent the mass balance and energy conservation equations, hydraulic devices, control rules, and status checks.To conform to MILP structure, MILPNet employs piece-wise linear approximation and integer programming.MILPNet was implemented and tested using Gurobi Python API.Modeling accuracy was shown to be comparable to EPANET, a public domain software for hydraulic modeling, and sensitivity analyses were conducted to examine the impacts of the modeling assumptions on the performance of MILPNet.Additionally, application of the framework was demonstrated using pump scheduling optimization examples in single and rolling horizon scenarios.Our results show that MILPNet can facilitate the construction and solution of optimization problems for a range of applications in WDSs operations.
• A mixed-integer linear programming framework (MILPNet) for formulating and solving water distribution system optimization problems is presented • MILPNet models system dynamics, hydraulic devices, control rules, and status checks and is flexible to adding more devices and conditions • The optimization model can be generated from a .INP file and case-specific objectives and constraints can be specified via Python interface

Supporting Information:
Supporting Information may be found in the online version of this article. 10.1029/2023WR034526 2 of 23 WDSs thus leads to their application in a wide range of optimization problems.Common optimization problems involve pump scheduling (Housh & Salomons, 2019), pressure management (Price et al., 2022), state estimation (Vrachimis et al., 2019), and real-time control (Preis et al., 2011).A review of common optimization problems in the context of WDSs can be found in Mala-Jetmarova et al. (2017); Bello et al. (2019); Creaco et al. (2019).This paper introduces MILPNet-an adjustable optimization framework based on mixed integer linear programming for water network operations problems.
In the context of previous studies, the approaches for solving the various WDS optimization problems can be broadly classified into two main categories: meta-heuristic and mathematical programming methods.Meta-heuristic methods rely on external solvers to evaluate the governing equations of WDSs dynamics and performance of the proposed solutions (Chopard & Tomassini, 2018).These methods decouple WDS simulation and optimization and develop solutions based on exploratory search, but are often heavily time-and computation-intensive, and can only guarantee local solutions.On the other hand, mathematical programming methods allow to explicitly model the governing hydraulic equations within the optimization problem as constraints, and represent the performance as objective functions (Coelho & Andrade-Campos, 2014;Collins et al., 1978;Sela Perelman & Amin, 2015).Comprehensive and efficiently designed mathematical programming models can provide optimal or near-optimal solutions even for large networks without relying on a multitude of hydraulic simulations, thus typically requiring shorter running times and lower computational effort (Bragalli et al., 2012;Wang et al., 2020).Therefore, due to advantages with regards to accuracy and computational complexity, we examine only mathematical programming techniques in this paper.
The non-linear system of equations that represents WDS hydraulics, combined with the discrete states and controls of hydraulic elements, are most accurately captured in a mixed-integer non-linear programming (MINLP) formulation.A wide variety of optimization problems can be modeled as MINLP and a range of global and local solvers (commercial and open-source) have been developed to solve general MINLP problems (Kronqvist et al., 2019).However, despite extensive progress, these problems are generally very difficult to solve.In the context of optimization in WDSs, there are two dominant paradigms to address the challenge of MINLP problems: (a) create approximations and assumptions to reduce model complexity for the application of modern solvers, such as CPLEX and Gurobi (CPLEX, 2009;Gurobi Optimization, LLC, 2022), and (b) develop tailored algorithms or in-house solvers to obtain accurate solutions.
Many studies attempt to solve a variety of WDS optimization problems by relying on approximations or network simplifications to fit linear programming (LP) or mixed-integer linear programming (MILP) problems that can be solved using off-the-shelf solvers.For example, Pasha and Lansey (2009) generated an LP approximation for the pump energy balance equations for simplified versions of WDSs to optimize pump scheduling.Liu et al. (2020) used a pre-processing method to refine their MILP formulation for networks containing pumps, tanks, and check valves using piece-wise linearization for demand response optimization.Menke et al. (2016) also used piecewise linearization to construct a problem with MILP constraints and a quadratic objective function to optimize pump scheduling to provide demand response services.Vieira et al. (2020) explored an approach to optimize WDS power consumption costs with linear envelope relaxations and include constraints to model tank storage limitations.Other studies applied alternative approximation and relaxation methods to solve non-linear WDS optimization problems using off-the-shelf solvers.Eck and Mevissen (2012), Gleixner et al. (2012), andBonvin et al. (2017) used quadratic approximations to represent the pipe head loss equations to construct optimization problems that could be solved with modern solvers.Cassiolato et al. (2022) applied convex hull relaxations and disjunctive programming techniques to solve an MINLP problem to minimize WDS operation and installation costs using the global solvers BARON and SBB (Bussieck & Drud, 2001;Tawarmalani & Sahinidis, 2005).
Other researchers have developed tailor-made algorithms that can solve MINLP WDS optimization formulations faster and more accurately than off-the-shelf solvers by using new solution approaches or by utilizing relaxation, reformulation, and penalization techniques.Eiger et al. (1994) and Sherali et al. (1998) approximated the non-linear head loss equations with polyhedral outer approximations, the bounds for which were derived from branch-and-bound procedures.Outer approximations have also been employed by Shi and You (2016) and Tasseff et al. (2022) to optimize pump scheduling.Pecci et al. (2019) previously developed a novel branch-and-bound method to solve the MINLP valve placement problem for WDSs, and this work was further developed by Ulusoy et al. (2020) to enhance system resilience.Branch-and-bound techniques were also adopted by Bonvin et al. (2021) within a two-step algorithm incorporating integer linear relaxation to optimize pump scheduling. 10.1029/2023WR034526 3 of 23 Pecci et al. (2021) developed an algorithm combining polyhedral relaxation, bound tightening, and a rounding heuristic to optimize an MINLP problem describing the placement of pressure-reducing valves (PRVs) and chlorine booster stations in WDSs.Ghaddar et al. (2017) utilized semidefinite programming relaxations to solve the MINLP valve setting problem to reduce network leakage using a new quadratic energy conservation approximation.Dai and Li (2014) reformulated the WDS MINLP formulation optimizing PRV placement into a mathematical program with complementarity constraints that is penalized and solved as a series of NLP problems.While these tailored algorithms potentially provide quicker and more accurate solutions compared to approximated models solved using off-the-shelf solvers, these in-house tailored algorithms and solvers are often restricted to the research group that developed the procedures, and replication and generalization of the approaches to other case studies can prove challenging.
It follows that although a wide range of mathematical programming approaches have been applied to WDS optimization problems, the resulting methods have had to trade off between accessibility and accuracy.Basic LP and MILP formulations of WDS hydraulics do not guarantee accurate results due to approximation and incomplete modeling of all facets of the WDS.On the other hand, although tailored algorithms preserve non-linearities and may outperform off-the-shelf solvers, they require complex, heavy duty optimization procedures and are not easily transferable to other problems and case studies.More importantly, most works that rely on mathematical programming to model WDS optimization overlook control devices (such as PRVs that have multiple states) and status checks (such as the closing of the tank inlet link when the tank is full or empty).
This paper aims to address the gap between accessibility and accuracy by proposing MILPNet, a general optimization framework for optimization problems related to WDSs operations.MILPNet automatically reformulates WDS hydraulics, hydraulic device states, and controls into a MILP structure which can be integrated into a wide range of optimization problems.Leveraging the MILP formulation enables using modern off-the-shelf solvers that, combined with hardware speedups, can solve large problems efficiently and avoid the need to develop highly specialized solvers.MILPNet can model: (a) system hydraulics (i.e., mass and energy conservation); (b) common hydraulic elements, such as tanks and pumps; (c) control elements, such as gate valves (GVs) and pressure reducing valves (PRVs) that can have multiple states with different characteristics; and (d) status checks and control rules common to many hydraulic systems, such as pumps turning on/off when tanks are empty/full.We demonstrate the application of our approach using Python with the Gurobi solver.Users can provide a hydraulic model using a standardized input (e.g., EPANET) or create the model directly via the programming interface, and augment the optimization model by specifying objectives and additional constraints for case-specific applications.A reproducible, open-source set of scripts that can be generalized to other network models and replicated using the Gurobi Python API can be found in the MILPNet Github repository: https://github.com/meghnathomas/MILPNet.
This paper is structured as follows.In Section 2, we show a detailed formulation of all the hydraulic equations, devices, and controls using MILP.To demonstrate that MILPNet accurately models WDSs hydraulics and controls, in Section 3.1, we first compare MILPNet with EPANET simulations, and then in Section 3.2 explore sensitivity of solution accuracy and running times to problem size and approximation tightness.In Section 3.3, we further validate MILPNet ability to emulate WDS dynamics using eight different benchmark networks.In Section 3.4, we demonstrate the application of MILPNet to optimal pumping scheduling problems.Finally, Section 4 outlines future research directions and concludes this paper.

Methodology
In this section, we describe in detail the formulation of the MILPNet optimization problem starting with the decision variables, followed by a description of the problem constraints based on system hydraulics.Next, we describe a set of constraints to check (and, if necessary, modify) the status of control elements at each time step in accordance with conditions and control rules placed upon them.Finally, implementation aspects are specified, for example, setting up the model and specifying solver parameters.

Decision Variables
For an optimization problem that is formulated over a specified horizon the decision variables in MILPNet include continuous and binary variables at each time step t = 1, 2, …, T. Continuous variables represent the flow in each link (i.e., pipes, pumps, and valves) and pressure head at each node (i.e., junction, tank) in the network.Binary variables include: (a) status variables representing the operational status of the hydraulic devices that can be controlled (e.g., pipes, pumps, PRVs, GVs), and (b) auxiliary variables that inform change in operational status.Table 1 summarizes the core decision variables of the MILPNet optimization framework.

Modeling System Hydraulics
The first set of constraints includes the fundamental mass and energy principles that govern the flows through links (i.e., pipes, pumps, and valves) and pressure in nodes (i.e., junctions, tanks, and reservoirs) in WDSs (Boulos et al., 2006).For each node i adjacent to node j = 1, …, N, the mass balance equation can be formulated as where     is the flow rate in pipe l connecting nodes i and j, and     is the demand at node i at time step t.Pipes in a WDS typically impart a loss of energy (i.e., head loss) to the water they convey due to friction.For each pipe l, the head loss  Δ   at time step t is formulated as: where     and     are the heads at nodes i and j, respectively.
MILPNet employs the Hazen-Williams friction model to derive head loss as a function of pipe properties and flow rate, which is formulated as: where pipe l connecting nodes i and j has diameter , e 1 = 1.852, and e 2 = 4.87, and α is the unit constant (Boulos et al., 2006).Since flow can be bidirectional, the inclusion of the head loss equations results in an optimization problem that is non-linear and non-convex.A range of linearization methods have been used in the literature to represent the non-linear relationship between head loss and flow rate, for example, piece-wise linearization (Geißler et al., 2011;Menke et al., 2015), using the Jacobian (Martinez Alzamora et al., 2014;Ulanicki et al., 1996), and Taylor series expansion (Xu & Goulter, 1998).MILPNet approximates the above head loss equation by employing piece-wise linearization, that is, the relationship between head loss and flow is approximated using a series of consecutive linear segments.Figure 1 shows piece-wise linearization of the pipe head loss equation with 1, 3, and 5 segments.Details regarding the linearization procedure can be found in Bradley et al. (1977).

Modeling Hydraulic Devices
In this section, we describe the modeling of the additional hydraulic devices that are common to many WDS networks: storage tanks, pumps, GVs, and PRVs.The table shown in Figure 2 summarizes the modeling conditions and constraint formulations for each hydraulic device.

Tanks
Tanks allow storage within the system and can fill or empty to satisfy system dynamics or follow operating rules.
If a system includes storage tanks, then additional constraints are defined, which account for changes in tank water level between consecutive time steps.The mass balance for tanks is formulated as where   −1  and     are state variables representing the tank head at time steps t − 1 and t respectively, Δt is the time step length,   −1  is the flow rate of water into the tank through link tl at time t − 1, and A tk is the cross-sectional area of the tank.

Pumps
Pumps are modeled as links that add pressure head to the system.If a pump p connects nodes i and j and has open status at time step t, head gain  Δ   is awarded to the downstream node j according to a head-flow curve of the form: where     is the flow rate through the pump, and A and B are constants unique to each pump (Boulos et al., 2006).When the pump has closed status, the start and end nodes are disconnected from each other, that is, there is zero flow through the pump, and the start and end node heads, H i and H j , have no functional relationship and can take on values independent of each other.In order to model the two distinct states of the pump, we introduce new continuous auxiliary variables     and status variables     , modify Equation 2 as follows and add the following big-M constraints: where M 1 is an arbitrarily large number.Equation 7 models the following conditions: at each time step, the pump status is determined by the binary variable  6is allowed to assume any value and disconnect nodes i and j.In order to maintain the linearity of the problem, the quadratic head gain function is approximated 10.1029/2023WR034526 6 of 23 using a piece-wise linear formulation similarly to the pipe head loss equations.Pumps do not allow back flow, and therefore     is lower-bounded to restrict the flow and allow only non-negative values.The logic, as well as implicit and explicit constraint formulations, are shown in Figure 2. Figure S1 in the Supporting Information S1 shows piece-wise linearization of the pump head loss equation with 1, 3, and 5 segments.These constraints were implemented in MILPNet using Gurobi indicator constraints (Gurobi Optimization, LLC, 2022).

Gate Valves
Gate valves (GVs) are links that can be assigned open status to allow flow under regular head loss conditions or closed status to disconnect the start and end nodes.The MILP formulation is achieved similarly to pumps: we introduce new auxiliary and status variables, modify Equation 2, and add big-M constraints to enable link closure.For GV g connecting nodes i and j, we add a continuous variable     to Equation 2 to derive the following: ) .These relationships are explicitly modeled for each time step in Equation 11, shown in Figure 2, where M 2 is an arbitrarily large number.

Control Valves: Pressure Reducing Valves
PRVs are designed to limit the head at the end node of the valve based on a specific pressure setting value.Each PRV connecting nodes i and j can have either active, open, or closed status at each time step depending on the upstream, downstream, and setting heads of the valve, that is, H i , H j , and H set , respectively.The common PRV operating rules are as follows (Rossman et al., 2020).If the PRV upstream node head, H i , is greater than the head setting value H set , then the PRV status is set to active and additional head loss is induced such that H j is equal to H set .If H i is lower than H set , then the PRV status is set to open and H i and H j are set to be equal (MILPNet assumes no minor head losses in the valves).Reverse flow is not allowed, hence if H j is greater than H i , then the PRV status is set to closed and the start and end nodes are disconnected.In order to determine the status of the PRV at each time step t, we introduce new continuous variables that determine PRV status within a set of big-M constraints.
The modified energy balance equation governing a PRV v is modeled as follows: where  Δ   is the head loss along the PRV from start node i to end node j,     is an auxiliary variable that enables to disconnect the two nodes, and     is a non-negative auxiliary variable that generates additional head loss if the PRV is in active status.The table shown in Figure 2 lists the logic and implicit and explicit constraint formulation (see Equation 12), where ɛ is an arbitrarily small number and M 3 is an arbitrarily large number.Other common control valves, for example, pressure sustaining valves and flow control valves, can be modeled similarly and added to the MILPNet framework.

Modeling Status Checks and Control Rules
To maintain proper functioning and engineering constraints, modeling WDSs involves including implicit and explicit rules that determine the states of common control devices (e.g., tanks, pumps, and valves) based on the system state (e.g., pressures and flows).For example, links leading into a tank (henceforth referred to as tank-links) are typically set to close at time step t + 1 if the tank is completely full or empty in time step t, or a set of control rules can force a pump to turn on (or off) when the tank level reaches specific low (or high) values (Creaco et al., 2019).These conditions are important to accurately represent WDSs, but are often overlooked in many optimization problems.In this section, we review the implicit and explicit constraint formulations for status checks and control rules.
MILPNet models the following status checks and control rules: (a) implicit event-based status check to prevent back flow in pumps and control valves; (b) implicit event-based status check and change of a tank-link if the tank is full or empty; (c) explicit event-based control rules changing the status of a link depending on the water level in a tank; and (d) explicit time-based control rules changing the status of a link at specific time steps.Status checks preventing back flow for pumps and PRVs are implicitly modeled in Equations 10 and 12 since flow rate through these devices can only assume non-negative values.The other types of status check and control rule conditions are summarized in the table shown in Figure 5 and are described below.

Tank-Link Status Check
MILPNet prevents flow into a tank through a tank-link (Q tl ) by closing the link if the water level in the tank (H tk ) reaches the maximum allowable level and prevents flow out of the tank if the water level reaches the minimum allowable level . This status check is modeled for each time step by adding continuous variable u tl to the tank-link head loss constraint in Equation 2 similarly to Equations 6 and 8, introducing status variable y tl to inform tank-link status, and introducing auxiliary variables x 1 and x 2 to determine tank level violation.The status check condition is modeled explicitly in Figure 5, Equation 13. Figure 3 schematically illustrates tanklink conditions modeled in Equation 13by showing the time-series of H tk and y tl and the corresponding values of x 1 and x 2 for each H tk level.The values of the x variables at each time step correspond to the tank head regions: 3 illustrates that when the water level reaches the maximum or minimum levels, that is,    1 = 1 or    2 = 0 , y tl is set to 1, indicating that the tank-link will close, and when only    2 = 1 , then y tl is set to 0, indicating that the tank-link is open.Note that in the time step following a tank-link closure, the tank-link may re-open if the hydraulic solutions indicates that water is released from the tank (if it was previously full) or filled into the tank (if it was previously empty).

Event-Based Control Rules
The MILPNet framework allows modeling event-based control rules that inform link status changes according to tank levels.These rules are common in WDS operation; for example, it is common to turn on or off a pump based on a water level in a tank (Creaco et al., 2019).These rules are typically structured as pairs of instructions such as below, where the status of a control-link cl (such as a pump, valve, or pipe) is set to closed status after the tank level H tk meets or exceeds an upper control level  min  that represent tank capacity limits.Explicit event-based control rules differ from implicit tank-link status checks in the following ways: (a) control rules can relate the status of any link in the system to a tank level, and not just the tank-link; and (b) tank-link status changes are typically short-lived, and are usually reversed in the subsequent time step, that is, the tank-link usually closes only at the time step that the tank is completely full or empty, and will typically open in the next time step to allow the tank to release/draw water from the system.On the other hand, control-link status is maintained until a rule in the control rule pair is violated.
Event-based control rules are modeled for each time step by adding continuous variable u cl to the control-link head loss constraint in Equation 2, introducing status variable y cl to inform control-link status, and introducing status variables z 1 , z 2 , and z 3 to determine control rule violation.Equation 14 in Figure 5 describes how a combination of the control-link status at the previous time step and z status variables inform the control-link status .Figure 4 illustrates the control-link conditions by showing the time-series of tank head H tk and control-link status y cl and the corresponding values of z 1 , z 2 , and z 3 for each H tk level.We can clearly see from Figure 4 that the control-link status has closed status when the head upper control level is reached and the status changes to open status only when the lower control level is reached.

Time-Based Controls
MILPNet models time-based control rules as link status changes at specific time steps in the optimization horizon.This can be expressed for a control-link cl through rule statements of the form: The MILPNet framework includes constraints that set the status of a control-link (y cl ) to open/closed at the specified time step and maintains this status until another status change is enforced by a subsequent control rule.The constraint formulation of these time-based control rules are explicitly modeled in Equation 15 shown in Figure 5.

Setting up the Optimization Problem
Finally, all elements described previously can be combined to formulate an optimization problem modeled for a WDS.The constraints ensure that the hydraulics of a system containing junctions, pipes, tanks, reservoirs, pumps, PRVS, and GVs is modeled and governed by hydraulic equations, status checks, and control rules for each time step t = 1, 2, …, T of the operating horizon.The heads of reservoirs at every time step, initial tank heads, initial pump and valve statuses, and time-varying nodal demands are considered known inputs to the optimization problem.Given some user-defined objective function f, the general MILP optimization can be summarized as follows:  Eqs. ( 14), (15) ()

Solver Implementation
In this paper, all optimization problems are solved using the Gurobi API for Python (Gurobi Optimization, LLC, 2022).Gurobi employs the branch-and-bound algorithm, which first solves the optimization problem with relaxed integer variable restrictions and then searches for an optimal solution by re-introducing branches of paired integrality inequality constraints to solve mixed-integer programming (MIP) problems (Lawler & Wood, 1966).
In order to speed up and improve problem solutions, Gurobi also utilizes techniques such as problem reductions, cutting planes, various heuristics, and parallel solution approaches.To facilitate usability, MILPNet was developed for users that are familiar with EPANET.Users can define the WDS model using a .INP file in EPANET that is received by MILPNet and parsed to obtain the majority of model inputs and parameters (e.g., pump curves, tank minimum and maximum levels, control rules, time-varying demands and reservoir heads, etc.).Users can specify additional inputs, such as the number of segments used for piece-wise linear approximation, or the pipes that should be modeled as GVs (i.e., can have open or closed status) directly in the Python interface.MILPNet utilizes built-in functions offered by Gurobi to model piece-wise linear approximations, logic constraints (e.g., AND, OR, XOR), and indicator constraints (e.g., given linear constraint A and binary variable x, if x = 1, then A is true) (Gurobi Optimization, LLC, 2022).Big-M constraints are used extensively in this framework to model conditional statements.The M value for each set of constraints is tightly chosen to be just greater than the greatest bound of the associated continuous variable and is automatically generated by the MILPNet code.Similarly, since Gurobi does not model strict inequality constraints, MILPNet adds small ɛ value to strict inequality constraints (i.e., >, <), which assumes different values in flow rate and head constraints.If the solver is unable to find a feasible solution that aligns with the rigid model inputs, then MILPNet allows the Gurobi API to find a feasible solution by employing its built-in constraint relaxation technique.The stopping criterion for the optimization problem solution is the MIP optimality gap, which is the difference between the upper bound (the objective value of the incumbent solution) and the lower bound (the minimum of all current solutions).The default MIP optimality gap value used by Gurobi (0.01%) is used for all optimization problems solved in this paper, including problems that require relaxation.
We used the Python package WNTR (Klise et al., 2017) to parse the .INP file of the water network model to extract the network topology and run extended period simulations (EPS) with the EPANET engine to generate inputs and bounds for the problem.We used Gurobi version 9.5.1, on a machine with Intel(R) Core(TM) i7-7700 CPU @ 3.60 GHz and 32 GB RAM through a Python API (Python version 3.7.4).All MILPNet scripts can be found in the Github repository: https://github.com/meghnathomas/MILPNet.

Results
In this section, we first validate that MILPNet accurately models WDS hydraulics and controls by comparing the results of the feasibility problem with EPANET, test the sensitivity of solution accuracy and running times to problem size and approximation tightness, and test performance using eight different benchmark networks (Sections 3.1-3.3).Then, in Section 3.4 we demonstrate the application of MILPNet to optimal pumping scheduling problem using two WDSs under different energy cost structures.Table S1 in Supporting Information S1 summarizes the objective and constraints for each optimization problem in the results section.

MILPNet Capabilities
The MILPNet framework is first used to demonstrate modeling functionality and accuracy.We use network NET2 as an example and validate the optimization results against results obtained from EPANET (Rossman et al., 2020).NET2 is adopted from Rossman et al. (1994), and contains 35 junctions, 40 pipes, and one tank.One junction is assigned a negative demand flow and serves as a source.Figure 6 shows the layout of the system and additional details can be found in Rossman et al. (1994).In this section, all optimization problems are performed for an operating horizon of 24 hr and time step Δt of 1 hr.The non-linear energy conservation equations for pipes are modeled with three piece-wise linear segments (n p = 3) and the pump equations are modeled with five piece-wise linear segments (n pu = 5).
To get solutions only to the hydraulic problem, the optimization problem is solved as a feasibility problem where the objective function is set to f = 0. We demonstrate different conditions that can be modeled with MILPNet, including: Case (A): basic system hydraulics, Case (B): pumps, Case (C): control rules, and Case (D): PRVs.
It is important to note that we expect slight differences in results obtained with MILPNet and standard hydraulic solvers when modeling WDSs hydraulics and controls.First, the accuracy of the results depends on the number of segments in the piece-wise linearization.Second, in hydraulic solvers that perform EPS, such as EPANET (Rossman et al., 2020), the state of the system is determined using the forward solver approach, where the results are calculated based on steady-state equations that are solved sequentially for each hydraulic time step.In this approach, the set of equations does not include the state equations for the tank heads, which serve as boundary conditions and ultimately determine the flows and pressures in the system, that is, the state of the tanks is updated after the flows are determined.In MILPNet, all flows, and heads in the systems, and heads in the tanks are determined simultaneously for the entire operating horizon.Third, the MILPNet formulation performs operational changes (such as tank status change) at predefined time steps t, unlike hydraulic solvers (such as EPANET) which allow status changes of hydraulic devices within time steps.Fourth, rounding errors and coefficients, such as the unit constant (α) and diameter coefficient (e 2 ) in the Hazen-Williams equation (Equation 3), may cause small differences between EPANET and MILPNet results.Finally, to enable solution convergence EPANET allows PRVs to be closed even when the head upstream is greater than the head downstream, while MILPNet strictly enforces the rules modeled in Equation 12, that is, if the head upstream is greater than the head downstream the PRV has to be either open or active depending on the settings.Next, we show rigorous comparison with EPANET.

Case (A) Basic system hydraulics
We first examine the network model NET2 to evaluate the accuracy of MILPNet in a system with a single storage tank.This optimization problem includes constraints that model mass balance at nodes, energy balance over pipes, and tank-link status checks.We solve the optimization problem for the network model and compare our results for pipe flow rates and junction and tank heads with those obtained using EPANET.Figure 7 displays the tank level (left), the pressure head at the source junction (center), and the flow rate in Pipe-21 located in a loop away from the source (right) based on MILPNet (solid) and EPANET (dashed) simulation results.The close match between the pressure heads and flow rates illustrates the accuracy of the results obtained using MILPNet while satisfying the minimum and maximum allowed storage levels in the tank.The slight difference between source pressure head and pipe flow rate results can be attributed to the number of piece-wise segments, n p . Figure S2 in Supporting Information S1 shows MILPNet results for different values of n p .The maximum difference between MILPNet and EPANET does not exceed 10 −5 m for tank level, 0.62 m for the source node pressure head, and 0.6 m 3 /hr for pipe flow rate.

Case (B) Pumps
In this scenario, we replace the negative demand node in NET2 with an equivalent reservoir-pump system in order to demonstrate how MILPNet models pump behavior.In EPANET, the initial status of a pump will be maintained until an explicit control rule or implicit hydraulic condition (e.g., head exceeds shutoff head) forces a change in pump status (Rossman et al., 2020).This condition is an implicit assumption that EPANET makes in the simulation approach.Here, we solve an optimization problem in which the initial open pump status is 10.1029/2023WR034526 13 of 23 maintained throughout the simulation duration in order to mimic EPANET behavior.Figure 8 compares EPANET simulation results for tank level and pump flow rate with MILPNet results under the default EPANET modeling assumptions.We observe that MILPNet closely matches EPANET results and the pump remains open throughout the operating horizon while the tank fills up to its maximum level.Note that the optimization problem needed relaxation in order to obtain a feasible solution.We observe that the relaxation of the problem has a negligible effect on solution accuracy.The maximum difference between MILPNet and EPANET results does not exceed 0.15 m for tank level (below 0.8% of the average tank level) and 3.1 m 3 /hr for pump flow rate (below 2% of the average flow), which is below the precision level of measurement instruments (AWWA, 2018;Walski, 2021).Text S1 in the SI describes an alternate optimization problem formulation in which the pump is allowed to switch between open or closed states at each time step while still satisfying system hydraulics.Figure S3 in Supporting Information S1 shows the results for tank head and pump flow rate based on this modified modeling approach.

Case (C) Control rules
Next, we demonstrate MILPNet's ability to model event-based and time-based control rules.We begin with a typical event-based control rule scenario.The NET2 network (with the uncontrolled pump-reservoir system) was assigned the following event-based control rules linking pump status to tank level: The top row in Figure 9 shows the EPANET and MILPNet results for tank level and flow rate when these control rules are introduced to the model.We observe that both solvers close the pump once the tank water level reaches   The bottom row in Figure 9 shows the EPANET and MILPNet results for tank head and flow rate when the time-based control rules are introduced to the model.We observe that the MILPNet results closely match those of EPANET; the pump is open and the tank fills up at the beginning of the horizon until time T1; the pump closes at time T1 following the control rule, and subsequently, the tank begins to drain; the pump re-opens at time T2, at which point the tank begins to fill up again; and finally, the pump closes again at time T3 and the tank drains again.Figure 9 exemplifies the difference between event-and time-driven controls.The maximum difference between MILPNet and EPANET results does not exceed 0.42 m for tank level (below 3% of the average tank level) and 3.1 m 3 /hr for pump flow rate (below 2% of the average flow).

Case (D) PRVs (active, open, and closed)
We demonstrate the MILPNet model in the presence of PRVs with different states by adding a PRV to NET2 along Pipe-7 and modifying the start and end node elevations of the PRV in order to induce different PRV states in the NET2 model.Figure S4 in Supporting Information S1 shows, from left to right, the flow rate through the PRV, and the head at the start and end nodes when the PRV achieves (a) active, (b) open, and (c) a mix of active and closed states over the course of the operating horizon as a result of the hydraulic conditions in the system based on the rules previously mentioned in Equation 12. Figure S4 Supporting Information S1 demonstrates that MILPNet closely matches EPANET when the PRV achieves active or open status throughout the simulation duration.In the third scenario, both MILPNet and EPANET results demonstrate that the PRV is initially closed and then switches to active status, but there is a difference in the timing of the switching.As mentioned before, this discrepancy is attributed to differences in modeling assumptions by EPANET and MILPNet, number of segments, and solution horizon.The maximum difference between MILPNet and EPANET results does not exceed 0.5 m for node pressure heads (below 1% of the average pressure head) and 4.5 m 3 /hr for valve flow rate (below 5% of the average flow).

Sensitivity Analysis
In this section, we examine the effects of the (a) number of segments for piece-wise linearization, and (b) size of the WDS network model on the size, accuracy, and running time of the optimization problem.

Piece-Wise Linearization Approximation Accuracy
We test the impact of the number of segments in the head loss and pump curve approximation on MILPNet model accuracy and running times.Here, accuracy is quantified as the mean absolute percentage deviation between EPANET and MILPNet pressure head results.The MILPNet formulation is applied to NET2 (with a controlled pump and a reservoir in place of the negative demand node), and we vary between 1, 3, 5, 7, and 9 segments in the linear approximations of the pipe head loss and pump curve equations.This experiment is repeated 10 times for an operating horizon of 10 hr, and Figure 10 depicts the average running time (left) and accuracy of the MILPNet model (right) for each combination of number of pipe head loss linearization segments, n p , and number of pump curve linearization segments, n pu .We observe the following: (a) the running time is most sensitive to n p .There are several pipes in the model, and so, the linearization of each head loss equation for each pipe at each time step significantly increases the total running time; (b) there is no straightforward correlation between n pu and running time.For example, choosing n pu = 3 seems to be a better option than choosing n pu = 9, and surprisingly, n pu = 1, regardless of choice of n p .This performance could be attributed to the fact that the 1-segment approximation provides a less accurate representation of the pump curve and flow rate than the 3-segment approximation, and more computational effort is needed to find the optimal solution.On the other hand, the 9-segment approximation adds many constraints to the optimization problem and, therefore, also increases the running time; and (c) the pump approximation appears to contribute more to the accuracy of the overall simulation compared to the pipe approximation.The pump dominates most of the hydraulics in the system, while the pipes mostly have local effects; hence, including a more precise pump approximation in the optimization model is more important than improving the pipe approximation.Notably, there is little improvement achieved with respect to accuracy or running time when increasing n pu from 5 to 9.

Size of the Water Network
Here, we test the effect of the size of the water network on the number of variables, running time, and accuracy of the MILPNet formulation.Recall that the number of junctions in the network affects the number of mass balance equations in the optimization problem, and the number of pipes affects the number of head loss equations.However, all other conditions (e.g., status checks and control rules) remain the same.Since in large water systems the number of control elements and sources typically remains small, the following analysis primarily assesses the impact of the number of junctions and pipes in WDSs.We used the Python package MAGNets (Thomas & Sela, 2023a) to reduce the network model of NET2 to different degrees.Here, a reduced network model contains fewer junctions but is hydraulically equivalent to the full-size network.The reduction methodology is based on work by Ulanicki et al. (1996).We compare the full-size NET2 network (with controlled pump and reservoir) with reduced network models from which 25%, 50%, 75%, and 100% of the non-essential junctions were removed, using an operating horizon of 24 hr and n p = 3 segments and n pu = 5 segments for the piece-wise linear approximation.
Figure 11 (left) shows the number of status variables (representing element status) and continuous variables (representing nodal heads and link flow rates) present in the optimization model for the full-sized and reduced NET2 networks.We observe that (a) the number of status variables remains constant for each network.Since the status variables for all networks correspond to tank-link and pump status checks, and since the tank and pump remain unchanged in each reduced model, we can expect the associated status variables to remain unchanged despite changes in network size; and (b) the number of continuous variables decreases as the number of junctions and pipes in the network decreases.This decrease is expected because the continuous variables are linearly associated with the number of elements in the network.Figure 11 (right) displays the impacts of network size on MILPNet running times and model accuracy.We observe that as the number of junctions in the network decreases, the running time also decreases, likely due to the reduction of constraints in the optimization problem.The MILPNet model running time decreased from 121 s for the original network (with 34 junctions) by 97% to 4 s for the fully reduced network (containing 2 junctions).Importantly, the size of the network does not seem to impact the accuracy of the MILPNet results, as the deviation from nodal pressure heads in the full model simulated using EPANET remains approximately constant.This finding could be significant for optimization models that can benefit from simplified network model size and complexity without compromising accuracy.

Application to Benchmark Networks
We tested the MILPNet framework by solving the feasibility problem for 8 benchmark WDSs varying in complexity with respect to system size, hydraulic devices, and degrees of control.Table 2 summarizes the system characteristics, controls and devices, optimization properties, and original source for each system.Note that the original names of the systems were altered for illustration purposes.Each model was constructed with n p = 3 segments, n pu = 5 segments, and a 10 hr operating horizon.The results in Table 2 show that MILPNet is able to quickly solve the optimization problems for various system hydraulics.The running times range from less than 0.2 s for B-NET (a small network with no hydraulic devices or control rules) to under 2 min for D-NET (a network with changing PRV statuses).We observe that the inclusion of multiple control elements and control rules has a greater impact on the  running time compared to network size in terms of the number of junctions and pipes.For example, applying MILPNet to F-NET is 7 times faster compared to G-NET, despite the fact that F-NET has about 6 times as many junctions as G-NET.This significant difference in running times is due to the large number of hydraulic devices and control rules in the G-NET system, and consequently, the large number of binary variables included in the optimization problem.Note that control elements and control rules constrict the feasible space of the solution, thereby often prompting the need for model relaxation in order to obtain feasible solutions.MILPNet employed relaxation methods offered by Gurobi to obtain feasible and accurate solutions for the D-NET, G-NET, and H-NET models.Solution accuracy was evaluated in all cases similar to that demonstrated in Section 3.2, and resulted in a similar high performance.

Application to Optimal Pump Scheduling
In this section, we demonstrate the application of MILPNet to optimal pump scheduling problems.The optimization problems involve determining the operational status and flow rates supplied by the pumps to minimize the energy costs while subject to system hydraulics.We investigate the impact of varying cost structure patterns, penalizing number of pump switches, and selecting time step size on the problem solution.To test the effectiveness of our approach, we apply the scheduling problem to two networks: a small network, A-NET (Rossman et al., 2020), with one pump, a single source, and a storage tank, and a larger network, I-NET (Clark et al., 1995), with two pumps, two sources, and three storage tanks, over a single 24-hr operating horizon.Additionally, we provide an exploratory example demonstrating the application of MILPNet to achieve model predictive control.This is demonstrated by solving a week-long rolling horizon optimization problem for the A-NET network.

Single Horizon Optimization for a Small Network
We first test the performance on the A-NET system (Rossman et al., 2020) under three cost structures, c, to observe how the solution varies under (a) a time-varying cost structure, c 1 , which represents a typical energy tariff (here, as an example, we adopt the City of Austin energy tariff (City of Austin, 2021)); (b) a fixed cost structure, c 2 , calculated as the average daily cost; and (c) c 3 , a hypothetical time-varying cost structure to test model response.Figure S5 in Supporting Information S1 depicts the network topology and Figure 12 shows the cost structure patterns (red dashed lines).Our objective is to minimize the following: , the status at times t and t + 1 are the same, else the two statuses are different).The optimization problem is solved for each cost structure with n p = 3 segments, n pu = 5 segments, and an operating horizon of 24 hr.The solid blue line in Figure 12 shows the optimization results for pump flow rate under cost structure c 1 (top), c 2 (middle), and c 3 (bottom) with 1-hr time steps.We observe that all three optimization problems aim to minimize the pumping costs and number of switches, and the resultant pumping schedule is affected by the energy cost structure.In all three cases, the MILPNet optimization solution responds as expected by operating the pump when energy cost is low, and not operating the pump when energy cost is high, while satisfying the hydraulic constraints.The pump also rarely switches status due to the switching penalty (see Equation 16). Figure S6 in Supporting Information S1 shows the pump flow rates for the default pump setting (i.e., no optimization) and for the optimized scheduling problem under the three cost structures considering only the energy cost and without the switching penalty.Table 3 summarizes the pumping cost with and without applying the switching penalty, the number of switches, and the running times of each optimization problem.Results demonstrate that the number of pump switches is significantly higher when pump switching penalty is overlooked.
In order to investigate how the choice of time step affects the results, we resolve the 24-hr optimization problem (with switching penalties) using 30-min, 1-hr, and 2-hr time steps.Table 4 summarizes the pumping cost, number of switches, and running time for optimization problems with the three different time steps.Additionally, Figure 12 shows the optimization results for pump flow rate under cost structure c 1 (top), c 2 (middle), and c 3 (bottom) with 30-min (dashed line) and 2-hr (dotted line) time steps.We observe that increasing the time step from 1 to 2 hr significantly reduces running time, while decreasing the time step to 30 min considerably increases running times.The inverse relationship between time step size and running time can be attributed to the size of the resultant optimization problem.In a fixed-duration optimization problem, a larger time step results in fewer time steps and, therefore, a smaller optimization problem.Conversely, a smaller time step provides a more detailed solution at the expense of a larger optimization problem.Reducing the time step size leads to either an unchanged or decreased pumping cost, with the number of switches remaining constant.This phenomenon can be attributed to the fact that a more granular optimization problem allows for a broader array of possible feasible solutions and, thus the objective function value will either improve or remain constant.The choice of time step applied in the optimization problem should reflect the typical control time step for operational decisions.A manageable optimization problem size can be achieved by solving shorter duration problems (e.g., 12 hr instead of 24 hr) if smaller time step is desired.

Single Horizon Optimization for a Medium-Sized Network
The pump scheduling optimization problem is next applied to a larger network, I-NET (Clark et al., 1995), with 92 junctions, 116 pipes, two pumps, and three tanks.Figure S7 in Supporting Information S1 depicts the network topology.The objective is again to minimize the cost of operating both pumps where each pump is subject to cost structures c 4 , c 5 , and c 6 and penalize the number of switches (Equation 16).The cost structures follow the same pattern as c 1 , c 2 , and c 3 , but have different magnitudes and are applied to 2-hr time steps.The optimization problem is solved for each cost structure with n p = 1 segment, n pu = 5 segments, and an operating horizon of 24 hr with 2 hr time steps.Figure 13 shows the flow rates for Pump-1 and Pump-2 under the three cost structures and Table 3 lists the cost and number of switches with and without the switching penalty for each case.Results demonstrate that, similar to A-NET, all three optimized pumping schedules are affected by the energy cost structure applied to the problem, and that pumping costs and the number of pump switches are minimized.Figure S8 in Supporting Information S1 shows the optimization results for pump flow rates under all three cost structures without a switching penalty, and Table 3 summarizes the cost and number of switches with and without the switching penalty for each case.Similar to A-NET, we observe that even though the costs of operating the pumps are lower, the number of pump switches are significantly higher when pump switching is not penalized.These illustrative examples demonstrates how the proposed framework can be used to formulate and solve optimization problems by adding custom objective and constraint functions.

Rolling Horizon Optimization for a Small Network
The optimal pumping scheduling.There are many avenues to extend the current research for additional applications, such as the inclusion of more advanced approximations of the cost function, variable speed pumps, control settings, and real-time control.Here we demonstrate a straightforward application of MILPNet to facilitate the solution of a dynamic optimization problem where inputs may vary over time.For example, consumer demands are inherently uncertain, and hence, operational strategies that rely on prior assumed demands may not lead to optimal operation when the actual demands are revealed.Model predictive control (MPC) is a common control methodology that informs and periodically updates future operational strategies based on newer data in near-real time in a rolling horizon fashion (Castelletti et al., 2023;Salomons & Housh, 2020).We applied MILPNet to optimize the pumping schedule for network A-NET for a 7-day period under the daily cost structure c 1 using the MPC-based approach.Here, we assume that at each hour t a forecasting model generates a prediction for the demands for the following 24-hr period.The pump scheduling optimization problem is then solved for the next 24 hr (t + 24) given the forecasted demands and the known tank level and pump status from the previous timestep.The resulting pump controls for the succeeding time step t + 1 are implemented, and then the optimization problem is re-solved for the next 24 hr period with new inputs.This process repeats in a rolling horizon fashion as new data becomes available.Figure 14 shows the resulting tank levels (top) and pump status and cost (bottom) for a 7-day operating horizon.We observe that the pump operating strategy varies in each day due to the changing demand forecasts, but the pump is always turned off at times of high cost.Figure S9 in Supporting Information S1 further details the average demand forecast across the network (top) and the changes in the optimized pumping schedule as the demands are updated (bottom).The rolling horizon approach allows generating an operational strategy for the immediate future while taking into account the changing costs, uncertainty, and operational constraints in the next 24 hr.This method can be further adapted for larger networks by maintaining a tight approximation for the constraints in the near-sighted horizon and relaxing some of the approximations in the far-sighted horizon (Salomons & Housh, 2020).

Conclusions
This work introduces MILPNet, a mixed-integer linear programming framework that can be applied to a range of WDS optimization problems.MILPNet is shown to model WDS hydraulics, incorporate hydraulic control devices and rules, and achieve a level of accuracy comparable to widely used hydraulic solvers like EPANET.To facilitate usability, MILPNet was developed for users that are familiar with EPANET (Rossman et al., 2020), such that users can define the WDS model using a .INP file and specify additional inputs directly in the Python interface.The MILPNet framework is generalizable to a range of water networks and can be used to solve optimization problems using off-the-shelf solvers, thus making WDS optimization more accessible to modelers while accounting for comprehensive representation of system dynamics.Further development of the MILPNet model could include the modeling of additional hydraulic devices, such as variable speed pumps, pressure sustaining valves, and flow control valves.
In addition to expanding the scope of MILPNet and extending to other applications, future work is needed to improve upon the current framework.First, due to the differences between modeling approaches offered by hydraulic solvers and MILPNet (as outlined in Section 3.1), differences between simulation results are inevitable.These differences can be addressed by shortening simulation duration or using smaller time steps.Second, MILPNet employs a number of binary variables for each time step in order to model piece-wise linearization of head loss equations, hydraulic devices, status checks, and control rules.The addition of these binary variables considerably slows down model running times, especially for larger networks.A promising method to reduce running times is by incorporating reduced-size network models into the optimization procedure as our results show that MILPNet solution accuracy appears to be consistent for reduced and full-sized models (see Figure 11).Other potential approaches include increasing the MIP optimality gap criterion, reducing the number of approximation segments (n p and n pu ), and varying the big-M values.Future work may also investigate the integration of methods such as Bender's decomposition or Lagrangian decomposition to improve the efficiency of the MILP formulation (Ghaddar et al., 2015;Verleye & Aghezzaf, 2016).Finally, the MILPNet framework has only been tested using the commercial solver Gurobi (which provides a free academic license), and further research could explore the implementation using open-source solvers such as HiGHS (Huangfu & Hall, 2018), SYMPHONY (Ralphs & Güzelsoy, 2005), and GLPK (Makhorin, 2020).MILPNet includes Python codes and .INP network files to encourage users to further develop, improve, and extend the optimization framework.MILPNet is under continuous maintenance, improvement, and development.With the constant need for creating new methods for operating urban WDSs efficiently, this work contributes toward overcoming some of the existing challenges with modeling flexibility and accessibility.

Figure 2 .
Figure 2. Constraints to model hydraulic devices at each time step.

Figure 3 .
Figure 3. Tank-link status check: a time series of tank head H tk (top) and status variable y tl (bottom).

Figure 4 .
Figure 4. Event-based control rules: a time series of tank head H tk (top) and status variable y cl (bottom).

Figure 5 .
Figure 5. Constraints to model status checks and control rules.
Pump if Tank level above 20.4 m Pump if Tank level below 16.5 m
level and re-open the pump when the tank empties to the lower control level.A small delay can be observed between the EPANET and MILPNet results due to the discrepancy in solution approach; EPANET allows status changes within simulation time steps, whereas MILPNet only enforces status changes at discrete points within the simulation.Next, the following time-based controls were added to the network:

Figure 10 .
Figure 10.Heat maps depicting the impact of the number of segments in pipe and pump piece-wise linearization on: problem running time (left), and mean absolute nodal pressure head deviation (right).

Figure 11 .
Figure 11.Effects of NET2 system size on number of optimization problem variables (left) and problem running time and accuracy (right).

Figure 13 .
Figure 13.Pump flow rate for I-NET under pumping energy cost structure c 4 (top), c 5 (middle), and c 6 (bottom).

Figure 14 .
Figure 14.Week-long tank level (top) informed by model predictive control pumping schedule (bottom, blue solid line) for A-NET under pumping energy cost structure c 1 (bottom, red dotted line).
ul and is only set back to open status if H tk drops below the lower control level H ll :

Table 2
Summary of MILPNet Application to Benchmark Networks

Table 3
Case Study Optimization Problem Results for A-NET and I-NET Using MILPNet results of this work demonstrate the current capabilities of MILPNet framework leveraging MILP for operational problems in WDSs, such as the