Integrated Simulation and Optimization Models for Treatment Plant Placement in Drinking Water Systems

Drinking water systems are critical to society. They protect residents from waterborne diseases (Blake 1956; Melosi 2011; Sedlak 2014) and encourage economic success of businesses by providing consistent water supplies to industries and supporting a healthy work force (Blake 1956; Melosi 2011; Tarr 1996). Drinking water systems include treatment infrastructure to remove contaminants from water supplies and water distribution system (WDS) infrastructure to transport treated water to households and businesses. Drinking water quality declines as water moves through a distribution network and is sensitive to changes in water flow, driven primarily by temporal and spatial variation in consumer demand (Babayan et al. 2005; Clark 2012).

Maintaining high quality and sufficient quantity of consumed water is a priority of any drinking water utility and often requires adaptation of the WDS over time. Within an existing WDS, new infrastructure is added over time to meet a variety of needs, from new demands to improving resiliency to loss of a supply. These changes can include the design and addition of new infrastructure (e.g., tanks and pumps), the design of replacement or modified infrastructure components (e.g., resizing storage or adding in-network treatment), or the adjustment of existing WDS operating protocols (e.g., modifying operating schedules for pumps and tanks). Any system modification affects water flow and quality in the network, sometimes in unanticipated ways (Loucks 2002). Therefore, it is important to consider the water quality implications of modifying distribution infrastructure. Previous work considered the placement of chlorine boosters stations (Ayvaz and Kentel 2015; Behzadian et al. 2012; Boccelli et al. 1998; Lansey et al. 2007; Ohar and Ostfeld 2014; Prasad et al. 2004; Tryby et al. 2002, 1999) and water quality sensors (Berry et al. 2006, 2009; Hart and Murray 2010) within distribution systems.

Resiliency and redundancy are key metrics for drinking water system design (Committee on the Beneficial Use of Graywater and Stormwater et al. 2016; Dandy and Engelhardt 2006; Gheisi et al. 2016; Herman et al. 2015; Loucks et al. 2005; Milly et al. 2008). The addition of multiple sources has been suggested as a way to improve resiliency (Barker et al. 2016; Committee on the Beneficial Use of Graywater and Stormwater et al. 2016; Kang and Lansey 2012; Newman et al. 2014; Okun 1997; USEPA 2012). Incorporation of new sources can involve aggregation and consolidation of neighboring systems, incorporation of new freshwater sources, such as groundwater (Vieira et al. 2011; Vieira and Cunha 2016), and integration of treated stormwater or wastewater sources (Barker et al. 2016; Committee on the Beneficial Use of Graywater and Stormwater et al. 2016; Kang and Lansey 2012, 2014). In each of these cases, the location and sizing of new sources are key design questions and have significant effects on the water quality.

Planning models are often used to evaluate such major infrastructure changes (Brown et al. 2011; Pahl-Wostl et al. 2011; Sedlak 2014). There has been extensive research in WDS planning models over the last 30 years using physically-based simulation models (Loucks et al. 2005; Mays 2000), numerically-based optimization models (Biegler and Grossmann 2004; Grossmann and Biegler 2004; Liebman 1976; Loucks et al. 2005; Revelle et al. 2004), and combinations of the two (Goulter 1992; Grayman 2006; Maier et al. 2014; Razavi et al. 2012a; Shamir 1983). The purpose of these models ranges from the selection of water sources to design and placement of treatment choices for new (Kang and Lansey 2012; Newman et al. 2014) and modified (Boccelli et al. 1998; Lansey et al. 2007; Tryby et al. 2002) water systems. Because centralized treatment and distributed delivery is the standard model for design and engineering of drinking water systems (Gleick 2003; Sharma et al. 2010), planning models usually evaluate a diverse range of treatment approaches, but simplify or exclude the distribution system effects associated with those alternatives. When distribution water quality effects are considered, only a small range of treatment solutions is considered (Tryby et al. 2002). Expanding water planning models to allow for integrated consideration of multiple treatment options (especially placement of treatment or integration of alternative supplies) while optimizing delivered water quality through the WDS will improve design of new systems as well as inform long-term planning for systems undergoing updates and expansions. The present work introduces a coupled simulation–optimization model structure designed to assess how the placement of new treatment infrastructure affects the potential for water quality challenges within the distribution system.

Simulation models attempt to predict consumed water quality under defined conditions for source water quality, treatment plant choices, and consumer demands (Clark 2012; Mays 2000; Rossman 2000). Optimization models seek to select treatment choices that ensure the highest water quality under constraints of source water quality, limits of treatment technology, and requirements of meeting regulations and continuous and variable consumer demands. When applied to infrastructure decisions, simulation models can be used to predict expected water quality changes for multiple alternatives, and optimization models can be used to narrow these potential choices to those that contribute to the system’s objectives.

Simulation models are most accurate for drinking water systems that are well defined, in which detailed physical system configuration and specific user demand patterns are known (Clark 2012; Rossman 2000). A system configuration includes pipe sizes and connections that compose the distribution system, as well as the location and capacity of treatment facilities. Because simulation models predict system performance and water quality, they can be used to determine the likelihood of meeting regulatory requirements before construction investments are made (Boccelli et al. 1998; USEPA 2003; Tryby et al. 2002; Xu et al. 2010). Simulation models are, however, limited due to their computational intensity (Broad et al. 2005; Isovitsch and VanBriesen 2008; Loucks et al. 2005; Razavi et al. 2012a) and the applicability of results to only the specific system configuration(s) considered (Clark 2012). When changes to the physical components of the systems are made, simulation model results must be revised (Boccelli et al. 1998; Rossman 2000). Even when considering many potential choices, engineering judgement is used to select a limited subset of possible configurations to make the modeling and analysis tractable. This can lead to overlooking effective configurations (Reed et al. 2013).

In contrast, optimization models compare multiple infrastructure configurations simultaneously and provide evaluation of each decision alternative with respect to specific objectives. An objective is selected based on the primary goal of the utility at the time. The supply of microbiologically safe water, which often has been represented by ensuring sufficient quantities of a disinfectant are added to the water (Boccelli et al. 1998; USEPA 2002), has been an objective since the inception of treatment (Blake 1956). Another water quality measure is water age, which serves as proxy for overall system water quality (USEPA 2003). Cost minimization is also a common objective, and can include both capital and operating costs (Farmani et al. 2005; Kang and Lansey 2012; Morgan and Goulter 1985; Padula et al. 2013; Shamir 1974). The objective function and decision variables are defined by the problem and decision context. For example, Boccelli et al. (1998, 2003) and Tryby et al. (2002) minimized the chlorine disinfectant dose per day (to meet a water quality objective of minimizing disinfection byproduct formation) through the selection of disinfectant injection locations and the dose required to meet regulatory residual concentration standards. To evaluate water quality sensor placement, Berry et al. (2006) and Ostfeld et al. (2008) considered the impact of a specified contamination event or set of events affecting water quality. Kang and Lansey (2012) and Newman et al. (2014) used optimization techniques for system-level planning to size reservoirs and new piping infrastructure based on capital and operating cost objectives.

The formulation of the optimization model is designed with the specific objectives and planning stage in mind and can be formulated in multiple ways. Optimization models can typically be described as either deterministic or heuristic. Deterministic optimization approaches solve the optimization model directly and identify an optimal solution. Deterministic approaches were common during the initial demonstrations of optimization models applied to drinking water infrastructure problems (Morgan and Goulter 1985; Shamir and Howard 1979; Walski et al. 1987), but their use has been limited in recent years. Boccelli et al. (1998) and Tryby et al. (2002, 1999) demonstrated how a mixed-integer linear program (MILP) model and deterministic solution approach could be applied to placement and operation of booster stations, to minimize the chlorine mass required. Berry et al. (2006, 2009) used a similar formulation to determine optimal placement of sensors and demonstrated both deterministic and heuristic approaches to minimize the effects of a contamination event. Berglund et al. (2017) used an iterative deterministic optimization model framework to predict leak detection.

Heuristic approaches, of which evolutionary algorithms are the most common, allow for nonlinearities (e.g., maximizing distribution system water quality) to be included directly in the optimization model (Bi and Dandy 2014; Broad et al. 2005; Maier et al. 2014). However, heuristic approaches are search techniques, which are designed to explore the decision space, and they do not guarantee an optimal solution. Furthermore, they require assumptions in the selection of the search parameters that are not always well understood (Maier et al. 2014). Although these approaches are used widely in research, their implementation has been limited in the engineering community (Goulter 1992; Hart and Murray 2010; Maier et al. 2014).

One benefit of MILP models and deterministic optimization solutions is their simplicity. In addition, these types of models are more straightforward to communicate to decision makers and any assumptions made in the model are incorporated in the formulation step, and therefore can be clearly communicated through discussion with the decision maker at this stage. With the proper framing and by pairing the simplified MILP optimization model with improved simulation models, integrated deterministic optimization with simulation models may present a useful alternative to drinking water system design, and the benefits of these optimization models are worth reconsidering.

However, system hydraulics are nonlinear and must be captured while maintaining a linear optimization formulation. The nonlinear optimization form including system hydraulics can be simplified to a linear form (Alperovits and Shamir 1977). The simplified linear optimization model can be used as a screening model (Cohon 2003; Loucks et al. 2005) in which the decision space is refined. The optimal alternatives found by the screening model can then be evaluated with a hydraulic simulation model to verify water quality performance. This approach links a simplified linear optimization model for screening to a more complex simulation model.

Facility siting models often take this approach and have been used to determine the placement and sizing of many types of infrastructure (e.g., fire stations and municipal waste sites) (Dresner and Hamacher 2002). This type of formulation, which often takes the form of a MILP, is efficient for placement problems and provides solutions that are guaranteed optimal under the defined constraints. Drinking water infrastructure components analyzed with facility siting models include reservoirs and treatment systems (Kang and Lansey 2012; Newman et al. 2014), chemical dosing boosters (Boccelli et al. 1998, 2003; Tryby et al. 2002), and water quality sensors (Berry et al. 2006).

There are two primary ways to link a simulation model with an optimization model, embedded and iterative. Fully embedding the simulation model requires the simulation model to be run every time the objective function is evaluated. This approach is more common with heuristic optimization models and can be exceptionally computationally intensive. Both embedded and iterative approaches are used in the peer-reviewed literature (Razavi et al. 2012b). In the iterative (e.g., sequential) approach there are two ways to link the optimization and simulation models: (1) constraining the possible set of decision alternatives through scenario analysis, which reduces the number of computationally intensive scenarios which must be simulated, and then those scenarios alone are evaluated in the optimization model (Barker et al. 2016; Herman et al. 2015; Kang and Lansey 2014); or (2) using a simplified representation of system-level water quality behavior in the simulation model for all possible decision alternatives to inform the optimization model (Berglund et al. 2017; Boccelli et al. 1998; Condon and Maxwell 2013; Tryby et al. 2002).

In the first technique, the simulation model is run for each defined scenario and the optimization model identifies a solution only from the predefined scenario set. For example, Newman et al. (2014) used scenario analysis to conduct an initial screening of the decision space for placement of new infrastructure on greenfield sites (sites with no prior infrastructure), and used optimization to choose among the screened scenarios. Kang and Lansey (2012) made similar greenfield assumptions to determine least-cost and low-emissions designs for dual water systems (in which treated storm or gray water treatment is placed parallel to a freshwater source) for drought-prone regions. Although the accuracy of the optimization result is not in question, this approach yields an overly constrained decision space because the feasible region is restricted to only the scenarios that have been simulated, which may result in the identification of suboptimal decision alternatives if the optimal alternative is not among the original scenarios.

In the second technique, the optimization model is structured to extrapolate water quality relationships from a subset of simulations and identifies solutions from a feasible region unconstrained by the simulated alternatives. When changes in the decision variables can be separated from the nonlinear system hydraulics, such as in the cases of chlorine booster placement (Boccelli et al. 1998; Tryby et al. 2002, 1999) and sensor placement (Berry et al. 2006; Isovitsch and VanBriesen 2008; Xu et al. 2008), this approach is straightforward. Because chlorine addition and sensors do not affect system hydraulics, mathematical simplifications (known as linear superposition) enables the integration of the simulated hydraulic flow behavior into the optimization model. Similar approaches are used when flow behavior does not induce changes in the water quality (Ayvaz and Kentel 2015; Helbling and VanBriesen 2009; Lansey et al. 2007; Propato and Uber 2004; Tryby et al. 2002). In each example, the optimization model can be run using the output from a single simulation model scenario. However, placement of source water supply and treatment infrastructure is expected to change the hydraulic behavior of the system and affects water quality. For these cases, multiple simulations must be run and combined to form the water quality data set used to inform the optimization model.

In the alternative approach presented here, the linear optimization screening model captures nonlinear system behavior by incorporating data sets as model parameters. The exogenous water quality data set is developed based on the simulation model output. From these simulations, the relationship between water quality and the treatment location variables is derived and the input data set is created. The number of simulations required affects the computational burden of the framework. Therefore, the number of simulations must be selected to balance computational challenges with how accurately the data set captures the water quality versus water supply location relationship. The combination of an efficient optimization model paired with an accurate hydraulic simulation model has the potential to address the analytical challenge that even a relatively small placement planning problem presents.

The present work provides a framework to evaluate the placement and sizing of flow-altering supply and treatment infrastructure within an existing distribution system without unnecessarily restricting the decision alternatives. Potential locations are restricted, but rather than use predetermined facility capacity configuration scenarios [as was done by Kang and Lansey (2012) and Newman et al. (2014)], which narrows the decision space and may cause inferior solutions to be considered, the selection of potential locations and their capacities are treated as decision variables in the MILP screening model proposed here (Loucks et al. 2005; Stedinger et al. 1983). To inform the optimization model, a physicochemical simulation model evaluates treatment decisions within an existing distribution system. The simulation model is run on a set of randomly selected configurations to predict water quality at specific locations within the network. The optimization model uses this output to interpolate water quality characteristics throughout the distribution system and select the optimal configurations of location and capacity. Interpolation allows for the decision space to remain unrestricted by the initial scenarios simulated. Following the optimization, the simulation model is used to verify optimization results. This structure allows the simulation model to be used in a limited capacity, in which accuracy is most needed, thus reducing the computational burden for the overall solution.

This integrated model framework is used to identify the location and select the capacity of treatment infrastructure within an existing distribution system so as to maximize the system-wide chlorine residual (as protection against microbial regrowth or intrusion). A schematic of the model framework is shown in Fig. 1.

The framework begins with a specification of the existing distribution system structure, including pipe sizes, connectivity, and existing water demands, and identification of where and what types of treatment infrastructure could be added (Fig. 1, Section 0). From this starting point, the framework proceeds through four steps: (1) simulate system hydraulic and water quality behavior and set the parameters required in the optimization model; (2) select potential system configurations through optimization; (3) validate the system hydraulic and water quality behavior through simulation; and (4) select network treatment configuration. Water quality predictions from Step 3 are compared with optimization model predictions from Step 2; if the difference in predicted water quality from each model is beyond the selected error tolerances, Steps 1–3 are repeated with an expanded set of preliminary simulations. The function and assumptions of each model and how they fit together within the model framework are described in the following sections.

Following the modeling framework described in sections “Background” and “Model Framework,” a small WDS commonly used in the literature, the EPA sample network Net3, was evaluated. The following sections describe each model step for this application in detail.

Drinking water systems are designed and operated to meet community water needs and protect public health. Identifying the best modifications and changes for existing drinking water infrastructure is difficult. Water quality relies not only on source-water treatment, but also on the ability to maintain high water quality throughout the distribution system. The modeling approach described in the present work evaluates a large decision space and identifies treatment or water supply configurations for existing systems and evaluates the performance of various configurations given the selected objective. Each stage of the framework provides additional information and insight into how the water quality performance of a drinking water system is related to treatment location and capacity decisions. The modeling framework presented here is proposed as an alternative to current evolutionary algorithmic approaches to integrated optimization and simulation models, which have not seen uptake by decision makers. The proposed approach integrates a simple deterministic mixed-integer linear optimization form with a traditional simulation model. The model framework allows for the decision maker to be included throughout the process, and all required assumptions are included in the formulation stages.

Demonstration of the method on Net3 illustrated the approach and provided an example of the information that can be drawn from the components of this model framework. For Net3, dividing treatment capacity across more than three plants resulted in a needless increase in infrastructure that did not improve water quality. The one-plant solution identified Site A as the optimal treatment location rather than a more centrally located site. The three-plant solution selected distributed capacity at the edges of the network. Assessing the effects of capacity allocation and system performance is also useful to water utility decision makers. Thus, the benefits of this modeling approach lie in how it elucidates the relationship among predicted system performance for the selected objective without prematurely narrowing the decision space. In a decision context, subsequent stages would reintroduce storage and evaluate how operations related to storage would have to be adjusted to avoid unnecessary water quality degradation.

Demonstration of the framework described here showed promise when tested on a simple network, EPA Net3. However, Net3 is a relatively small WDS; the framework needs to be demonstrated on a larger network to ensure it provides robust results for a range of WDS sizes. The approach described for selection of the preliminary simulation scenarios was sufficient for a small network, but it is possible that some networks may require additional preliminary simulations or an alternative approach to estimate the water quality variation associated with changes in treatment location across the network.

Water utilities rarely make decisions based on a single objective. Here, a chlorine residual objective was selected to prioritize microbiologically safe drinking water in selection of a treatment configuration. However, objectives change based on the distribution system and the specific goals a utility may be trying to achieve. In some instances, multiple conflicting objectives might be considered. The framework presented here is flexible and can be configured with a multiobjective optimization in Step 2. These modifications are included in Schwetschenau et al. (2019). Cost is typically a major driver for infrastructure planning decisions. In the case of Net3, due to economies of scale the single-plant solution would be expected to be more cost-efficient than the three-plant solution. Therefore, future work should include consideration of additional objectives, such as cost, that were not included in this version of the framework model.

The formulation of the MILP optimization model used for the demonstration case is shown in Fig. 2. The objective function is the sum of the minimum chlorine residual concentrations across all system demand nodes. Values of the parameter ${cl}_{i·j}$, the residual at node $j$ from a chlorine dose introduced at node $i$ are derived from the results of the simulations performed in Step 1.

This formulation, a capacitated facility location problem (MILP), has been the subject of many other studies and applied to other types of systems [an overview was discussed by Dresner and Hamacher (2002)]. The model decision variables include a binary variable ($y_i$) which is equal to 1 if the potential treatment facility, $i$, is selected, and 0 otherwise, and a continuous variable (${\mathit{cap}}_i$) representing the capacities of the selected treatment locations. The model also assigns a supply location to each demand node ($x_{i,j}$) by identifying the supply location that provides the best quality water to that node. The greater the number of demand nodes assigned to a specific potential supply location, the greater the capacity assigned to that location. The optimization model does not represent the hydraulic behavior of water through a network; thus Step 3 is intended to check that this simplification provides feasible results. Capacity constraints ensure that all demands are met and avoid results that propose very small plants. The capacity constraints set minimum and maximum feasibility bounds on the allowable plant capacity. There are fixed costs associated with any treatment plant, and the minimum capacity prevents an unrealistically small plant from being constructed. The maximum capacity is a nonbinding constraint in this case, but an upper bound may be necessary for some systems.

Fig. 3 describes the inputs and outputs from each step of the model framework and how they are used to generate the inputs for subsequent steps of the model.

Cost is a critical consideration when evaluating water treatment configurations, and although cost was not used as the primary objective here, it was computed after the fact to allow decision makers to used cost information in their selection of a configuration in Step 4.

For the configurations evaluated here cost can be estimated as follows:where $C_i^{\mathrm{cap}exF}$ and $C_i^{OMF}$are the fixed cost coefficients for plant $i$ capital and operating cost components; $C_i^{\mathrm{cap}ex}$ and $C_i^{OM}$ are the variable cost coefficients for each plant location $i$; $y_i$ is a binary variable, where 1 indicates that plant location $i$ is used; and ${\mathit{cap}}_i$ is the capacity allocated to location $i$. The budget parameter is the upper bound on cost set by the decision maker.

Cost parameters were estimated in Table 1.

Using this equation and the specified parameters (Table 1), the cost for each configuration presented in the paper was estimated (Table 2). It is not surprising that the costs increase with the number of plants in the configuration, but the cost of a specific configuration is useful information for a decision maker comparing the marginal improvement in water quality between configurations.

A database of water distribution systems (Jolly et al. 2014) was reviewed in the selection of a demonstration network for this model framework. Networks were considered based on several criteria. The first was a network for which an EPANET input file already existed and on which prior work was published. Second, a network must have a sufficient number of pipes and loops such that the solution would not be obvious or immediately apparent without a model. Therefore, commonly used networks such as the New York Tunnel System, Anytown, or the Hanoi water system networks, although common demonstration networks in the water distribution system design space (Andrade et al. 2016; Babayan et al. 2005; Basupi and Kapelan 2015; Fujiwara and Khang 1990; Walski et al. 1987), were not selected due to the limited number of pipes represented. Third, given the constraints associated with land acquisition and the assumption that any treatment facility or major pipe tie-in connection requires access, networks needed to have multiple locations and locations were only considered feasible if there was existing water distribution system infrastructure (i.e., a tank, treatment plant, or existing source). This limitation ensured that only feasible locations were considered for additional treatment infrastructure. In order to ensure a robust set of alternatives, only networks with a minimum of three feasible locations were considered. A set of candidate networks was selected based on these criteria, and then each was evaluated to ensure that pressure violations would not result from the use of any of the candidate feasible locations. For example, each network had to be able to meet all demands without error if all flow was routed through any of the candidate locations. Minor pipe-size adjustments were considered acceptable, but the addition of pipes and major sizing changes were not the focus of this work. After assessing the networks, the water quality variation across the network was visually compared as treatment was moved between candidate locations to ensure that differences were observed. Without these differences, the optimization could not be demonstrated. Several of the candidate networks were eliminated because of a lack of change in water quality based on treatment location. Net3 was the best network based on this evaluation.

Although Net3 is not currently the most popular system, there are many examples of its use to address these types of water systems questions in the past. Net3 was used in several water sensor placement studies (Berry et al. 2006; Yang and Boccelli 2014; Zechman and Ranjithan 2009). Other systems could have been be used to demonstrate the model framework, but for the sake of demonstration, Net3 had the necessary properties.

Some changes were made to Net3 to make it suitable for evaluation of the primary research question, demonstrating where treatment should be located within the existing distribution network based on the associated water quality. Changes include:

Some or all data, models, or code generated or used during the study are available from the corresponding author by request. Network data used here are available through the EPA and included with the EPANET Programmers toolkit (USEPA 2014). A few minor modifications were made to the default Net3 input files based on the needs of this work. The exact input files can be provided upon request. EPANET was run using a MATLAB version 2018A wrapper to allow for expanded coding capabilities (Hatchett and Uber 2015). GAMS was used as optimization software and the optimization formulation is provided in Appendix I. Specific GAMS scripts can be provided upon request. The methods by which the model components are integrated are described within the manuscript.

This work was performed with funding and support of the National Science Foundation Integrative Graduate Education Research Traineeship in Nanotechnology Environmental Effects and Policy fellowship program, Grant No. DGE0966227, the Pugh Fellowship, and a Dean’s Fellowship from the College of Engineering at Carnegie Mellon University. Access to the optimization software, GAMS, used in part of this analysis was provided by the Center for Climate and Energy Decision Making (SES-0949710) through a cooperative agreement between the National Science Foundation and Carnegie Mellon University.