Integrated Multiobjective Optimization and Simulation Model Applied to Drinking Water Treatment Placement in the Context of Existing Infrastructure

Drinking water systems provide safe drinking water, and consist of treatment and distribution infrastructure to transport clean water throughout their service areas (Blake 1956; Melosi 2011). The infrastructure necessary to ensure that water is sufficiently treated and transported and its quality is maintained until consumption vary based on the existing regional conditions and the objectives of the decision makers and stakeholders. Most urban areas in the United States have existing public drinking water supply systems; thus decisions typically involve modification of or additions to existing systems (Mays 2000). To properly identify the engineering solutions that best meet stakeholder objectives and assess the effect of these solutions on existing infrastructure performance, models and other decision tools are required (Loucks 2002).

Selection of an appropriate model framework depends on the needs of the specific stakeholders involved, the existing infrastructure to be integrated, and the type of decision to be made. Some common considerations of stakeholder groups (e.g., utility operators, customers, and city officials) include water quality, budget, redundancy, and source variability (Loucks 2002). A wide selection of models has been demonstrated on water system engineering decisions, providing a suite of potential models. Therefore, for any decision context, the model framework used should be informed not only by the engineering problem considered, but also by the specific decision makers involved. For example, simulation models are commonly used by utility operators to diagnose system issues, and a familiar and proven simulation model already in use by decision makers can be leveraged in a larger model framework. Furthermore, the level of interaction that a decision maker prefers will lend itself to different model selections. Highly interactive models allow a decision maker’s tacit knowledge to be embedded into the model (Greco et al. 2017). This can be more efficient and can produce alternative solutions preferred by the decision maker, but it can also lead to bias in the selection method and can be overly time consuming for decision makers (Greco et al. 2017). Furthermore, in situations with multiple stakeholders and goals, conflicts are common. Balancing the trade-offs among the solutions under consideration can be difficult. In this context, models, particularly multiobjective optimization models, can be a critical tool to assess these trade-offs and facilitate the selection of an appropriate solution (Revelle et al. 2004).

A variety of models have been demonstrated on many relevant water system engineering problems in the research literature. The most common approach includes a combination of simulation and optimization models applied through an integrated model framework (Goulter 1992; Hart and Murray 2010; Maier et al. 2014). In recent years, this optimization–simulation model framework has been extended to make integration of expert decision maker knowledge easier (Greco et al. 2017). However, within each model component, selections and assumptions are required. A more complex model is not better in all contexts, and model assumptions and level of detail should be tailored to the problem and the decision maker. For example, simulation models vary in the level of detail they capture, representing every pipe or just the primary ones, and in the selection of the water quality constituent kinetic reaction equations (Rossman 2000). Optimization model formulations can vary from linear representations of simplified water flow dynamics to highly nonlinear representations of water quality reactions or energy kinetics (Alperovits and Shamir 1977; Maier et al. 2014).

Although there is no shortage of optimization formulations or demonstrations of case studies applying these models to complex drinking water systems, the uptake of these models by drinking water managers has been slow (Goulter 1992; Hart and Murray 2010; Maier et al. 2014). In recent years the dominant modeling approach has included a fully integrated simulation model or surrogate simulation model (i.e., a numerical model fitted to predict the performance of the original simulation model) with a nonlinear multiobjective optimization model, typically solved using evolutionary algorithmic approaches (Behzadian et al. 2012; Bi and Dandy 2014; Broad et al. 2005; Nicklow et al. 2010). Although the value of this approach to find near optimal and creative solutions to specific research problems has been demonstrated, these models require selection of many modeling parameters and require many assumptions in the solution scheme. The implications of these algorithmic assumptions or how to select the best parameter values for a specific network are not always well understood (Maier et al. 2014). Depending on the system, the decision makers, and how well the decision space is understood, the uncertainty in these assumptions may be discouraging use of these models.

Schwetschenau et al. (2019) described a deterministic sequential optimization–simulation model framework that can be applied to the engineering problem of identifying the best location for a treatment plant or new supply tie-in within an existing distribution system based on a simple water quality objective. In that framework, a physicochemical simulation model predicts water quality and flow behavior within an existing distribution system for selection of feasible treatment solutions. A single-objective optimization model assesses these results to identify the treatment solution that best meets the identified objective. As described by Schwetschenau et al., compared with other approaches reported in the literature, the relative simplicity of the methodology has intuitive appeal; it also overcomes the artificial limits on the decision space that are a product of the computational burden of simulation models. However, the price is that approximate hydraulics are included in the optimization model, which require iteration to assure realistic solutions. Results from an application to a sample network indicate that the approach is viable.

The present work extends the optimization component of the framework to a multiobjective optimization model to allow for multiple conflicting objectives to be considered when identifying a preferred (rather than optimal) treatment configuration. The model expansion is described and then illustrated with a case study.

Over time, many objectives have been used by decision makers to select drinking water treatment alternatives. Providing a reliable supply of microbe-free water has been a primary objective since the initial design of water treatment systems (Blake 1956). A free chlorine residual reacts with pathogens, neutralizing them within the drinking water system. However, natural organic material and biofilm growth also decrease the residual concentration overtime (Haas 1990, 1999; LeChevallier 1999). Residual levels must be maintained such that biofilms and other residual reactants do not prematurely reduce a system’s protection against pathogenic intrusion events. Reducing contaminants formed within distribution systems, such as disinfection byproducts (Boccelli et al. 1998; USEPA 2002) and lead (Del Toral et al. 2013; Lacey et al. 1985; Levallois et al. 2014; Potash et al. 2015), also have been important goals for decades. As additional contaminants of concern have been identified, water age (i.e., the amount of time water spends in the distribution system moving from the point of treatment to point of consumption) has emerged as a general proxy for water quality (AWWA 2002; Machell and Boxall 2014). Recent contamination threats to source water supply have shifted focus to supply redundancy as a primary objective for water utilities. For example, Lake Erie summer algal blooms threatened Toledo’s water supply in 2011 and in subsequent summers (Michalak et al. 2013), and a chemical spill compromised Charleston’s water supply in January 2014 (Gabriel 2014). Source quality and quantity redundancy is also critical for regions that may experience climate change–induced source water changes (Georgakakos et al. 2014; Li et al. 2014). These new concerns do not replace the goals of consistent supply, disease prevention, and control of distribution system contaminants. Rather, they are added requirements, and they often conflict with existing goals. Individually, each objective might lead to a different engineering solution. Conflicting objectives can be difficult to balance or compare without optimization models designed to evaluate myriad solutions against competing goals.

Multiobjective programming (MOP) techniques have been applied to water systems for nearly as long as optimization models have been applied to water systems (Cohon 2003; Shamir 1983). Given the number of objectives of interest, it is not surprising that multiobjective optimization models have become more popular in recent years. Multiobjective applications have included cost versus energy efficiency trade-offs in pumping schedules (Ostfeld et al. 2014) and parallel water system design (Kang and Lansey 2012), cost versus water quality for booster station placement (Ayvaz and Kentel 2015; Prasad et al. 2004), and leak detection and control (Creaco et al. 2016; De Paola et al. 2017). Many of the preceding applications solve MOP using evolutionary algorithm techniques, and more recently many of these have been extended to include varying levels of decision maker interaction (Greco et al. 2017). Multiobjective models that have demonstrated the use of mixed-integer or linear optimization formulations, or the use of deterministic optimization approaches integrated with simulation models, are limited. This work extends the single-objective formulation presented by Schwetschenau et al. (2019) into a multiobjective framework.

One of the major benefits of multiobjective optimization techniques is the ability to compare the trade-offs among objectives and to evaluate a set of engineering alternatives. To this end, the solution method selected in this work had to identify nondominated solutions to estimate the Pareto frontier and the trade-offs. The noninferior set estimation (NISE) method (Cohon 2003; Solanki et al. 1993) was selected. NISE is a set generating technique that uses a geometrical approach to find a good approximation of the set of nondominated (i.e., noninferior) solutions. The NISE method represents a way to approximate the Pareto frontier and is related to the polyhedron surface approximation (Greco et al. 2017). NISE was selected over the more recent surface approximation methods due to its simplicity and the visualizations it provides. NISE systematically explores the objective space and quantifies the associated trade-offs among objectives. This process improves decision-maker understanding of the effect of engineering solutions across the objective space, leading to better-informed decisions.

The present work expands the model framework designed by Schwetschenau et al. (2019) to integrate multiple objectives into the mixed-integer linear program (MILP) facility-assignment optimization model component. The model framework has five steps (Fig. 1): (0) identify and describe an existing distribution system; (1) simulate system hydraulic and water quality behavior; (2) select potential system configurations through optimization; (3) validate the system hydraulic and water quality behavior through simulation; and (4) select preferred network treatment configuration. Fig. 1 was modified from Schwetschenau et al. (2019) to reflect changes made to Step 2. In Step 2 the single-objective optimization from Schwetschenau et al. (2019) is replaced with a multiobjective mixed-integer linear program optimization model (MOMILP). Similar to Schwetschenau et al. (2019), the MOMILP identifies the location and capacity of treatment that best meets the defined objectives. The Appendix describes this model framework in detail.

In the multi-objective optimization model, the decision variables include five binary variables, one for each of the five potential treatment nodes, defining whether that treatment node is selected to supply the system or not; and continuous variables for capacities of the selected supply locations. For each demand node in the system, the model selects a supply location by identifying the location that provides the best water quality to that node. The greater the total demand of all nodes assigned to a specific potential plant location, the greater is the capacity assigned to that plant location. Constraints, including node demands and capacity sufficient to meet demands, ensure that physical system requirements are met, including satisfaction of water demand at each node and chlorine concentrations restricted to the range specified by the Safe Drinking Water Act (USEPA 2003) ($0.02–4.0\mathrm{mg}/\mathrm{L}$). The lower bound [which is considered to be the minimum detection level (Rice et al. 2012)] is enforced in the model by prohibiting a demand node to be assigned to a specific plant location if the residual disinfectant concentration to that demand node would be below the lower bound. For cases in which a single plant cannot maintain the minimum residual at all demand nodes, additional treatment locations are identified, and single-plant solutions are considered infeasible (i.e., not a viable engineering alternative based on the defined constraints).

The detailed optimization model formulation selected for demonstration is summarized in Fig. 2. Objectives were selected based on a review of the literature and include cost and two measures of water quality, one indicative of public health and one representing the overall water quality throughout the distribution system. Specific formulations of each objective can vary and should be tailored to the goals of a specific system. Specific to this application, the three objectives include minimizing cost, maximizing chlorine residual concentration, and minimizing population-weighted water age.

The cost objective is defined as the sum of capital and operating costs based on typical unit-process costs (Clark 1982). Capital and operating costs are modeled as the sum of fixed and flow-based cost components with costs differentiated by supply source type (groundwater or surface water). The cost parameters used were defined by Schwetschenau et al. (2019).

The chlorine residual objective is defined as the sum of the chlorine residual concentrations across all system demand nodes; this leads to a value that is different from mean chlorine residual, which might be reported by a utility. The networkwide sum of chlorine residuals is maximized, because network protection from contaminants is improved with a higher chlorine residual. This objective is designed to promote improved protection against pathogenic intrusion in the event of pipe breaks or pressure transients, beyond what is provided by the regulatory minimum detectable, or approximately $0.02\mathrm{mg}/\mathrm{L}$ (Rice et al. 2012; USEPA 2006).

The water age objective is defined as the sum of the water age at each node, weighted by the demand at that node to ensure the best water quality for the most users. Ostfeld et al. (2008) used a similar objective to minimize the population affected by contaminated water consumption. This objective serves as a proxy for minimizing the exposure to trihalomethane (THM) species of disinfection by-product (DBP)s and other constituents that form over time in the distribution system. When used in conjunction with the chlorine residual objective described previously, the trade-off between increasing disinfectant residual to prevent microbial illness events and reducing water age to minimize exposure to carcinogenic THM DBP species can be evaluated.

Water age, and chlorine residual are computed at each node in the system based on the preliminary simulations conducted in Step 1. A different parameter value is assigned to each node for each treatment plant configuration simulated, and these are referred to as parameter sets. Parameter sets are populated using the same methods described by Schwetschenau et al. (2019). Briefly, parameter sets enable the prediction of system water quality by assigning demand nodes to specific treatment plants and determining the resulting capacity required for each treatment plant location in the system. Assignments are made based on maximizing the defined water quality objectives. This method does not account for system flow behaviors; therefore, hydraulic simulations in Step 3 are used to verify the system performance as defined by water age and chlorine residual concentration. To improve accuracy, the model is solved parametrically, and each iteration is associated with a different number of treatment plants utilized.

The NISE algorithm, a precursor to current polyhedral surface approximation methods (Greco et al. 2017), is applied to each pair of objectives in turn, generating paired trade-off curves or frontiers (Solanki et al. 1993). Results from the parametric iteration of the optimization model are aggregated for each objective pair to generate a single Pareto frontier (or nondominated set). This pairwise treatment of objectives guarantees that no portion of the feasible region will be overlooked (Cohon 2003). In integer programming, the feasible region is nonconvex and the resulting optimization is NP-hard (Greco et al. 2017). For such a problem, the NISE method can only find solutions on the convex hull. As a result, there may be nondominated solutions that the method cannot identify (Cohon 2003). Although the NISE method can miss select integer solutions, this is not a weakness of the method per se. Rather, it is a result of the nonconvexity of the feasible region, which poses a problem for all surface approximation methods, including heuristic approaches. The benefits of the NISE method include that the full breadth of the feasible region is evaluated systematically, and it is easy to visualize. These advantages outweigh the downside that it cannot guarantee that every integer solution is identified. Furthermore, in large problems for which an approximation of the nondominated set is sought, missing a few nondominated solutions is considered to be an acceptable trade-off for the other benefits of the method.

An advantage of a generating technique such as the NISE method is that decision makers are able to see the choices before having to set priorities among the objectives. The added insight into the trade-offs among the objectives, and the associated range of alternatives, lay the foundation for better-informed decisions (Reed et al. 2013). The present analysis evaluated three trade-off pairs among the three objectives considered: total cost, systemwide sum of chlorine residual, and demand-weighted cumulative water age. The formulated optimization model was solved with a branch and bound algorithm in GAMS version 25.1 (Rosenthal 1988).

As major infrastructure changes are made to existing drinking water systems, models are useful to address the effects of treatment and distribution choices on water quality. Modeling approaches are needed that provide a diverse set of optimal or near-optimal alternatives to decision makers and allow quantitative comparison across potential solutions. Furthermore, these models need to be capable of integrating multiple, potentially conflicting stakeholder objectives. Multiobjective generating techniques, such as NISE, produce trade-off curves and several nondominated alternatives which allow decision makers to evaluate their choices without having to define their preferences beforehand. This is the basis for better-informed decisions. The modeling framework developed by Schwetschenau et al. (2019) and expanded here presents an alternative framework that leverages a simpler form of optimization model (deterministic MILP), uses a native physicochemical simulation model that is familiar to most decision makers, and combines them with a deterministic sequential integration approach that is simple and allows for visualization of tradeoffs at each stage. The framework is designed to improve decision maker use of these models and to be flexible, allowing modification of submodels as needed.

The multiobjective formulation of this model framework was evaluated on the Net3 demonstration network. Three objectives were considered: cost, and two water quality objectives. The existence of a trade-off curve between the two water quality objectives identifies the need to evaluate multiple water quality objectives when evaluating treatment configurations. Of the five possible treatment locations considered, no more than three were required to provide optimal drinking water quality. Five different treatment configurations, using between one and three treatment locations, with different capacity allocations, defined the Pareto frontier surface.

The integrated optimization and physicochemical simulation model framework initially presented by Schwetschenau et al. (2019) is described to provide context for the multiobjective extension presented in this work. The model framework has five steps: (0) identify and describe an existing distribution system; (1) simulate system hydraulic and water quality behavior; (2) select potential system configurations through optimization; (3) validate system hydraulic and water quality behavior through simulation; and (4) select preferred network treatment configuration.

Step 0 includes identifying and specifying the physical parameters of the distribution system to be evaluated. Decision makers are consulted at this stage to determine the specific sites at which treatment can feasibly be located and to identify the objectives of primary concern to the decision makers.

Step 1 simulates the specified distribution system in EPANET 2.0 (Rossman 2000), a hydrodynamic water quality model. Each simulation requires the system’s treatment configuration to be specified (i.e., location and capacity of treatment water supplies added to the distribution system). Simulations generated at this stage are used to specify an input parameter set used in the optimization model. Simulations should be representative of the range of water quality results that could be expected from the network based on the range of possible treatment configurations. For example, the location at which supply capacity is placed in the system affects the water quality at different demand points in the network. Although the preliminary configurations simulated could be selected based on expert knowledge, here they are randomly selected, an uninformed set. The resulting water quality predicted at each node for each simulation run is used to generate a parameter set relating the predicted water quality at each node to the location of treatment in the system. This parameter set is intended to capture the nonlinear hydraulic and water quality behavior of the network and used in the MILP optimization.

Step 2 solves the MOMILP iteratively through application of the NISE method. The NISE method is a geometric approach to defining the Pareto frontier, quantifying the trade-offs among objectives and identifying the associated nondominated solutions. In Step 2, the NISE method is applied parametrically. For each parametric iteration, the total number of treatment locations is held constant and the Pareto Frontier is defined for variations in capacity allocation among the specified number of treatment locations. The Pareto frontier generated from each unique parametric iteration is aggregated, and the overall Pareto frontier can be determined from comparison of the objective values for each solution. The nondominated solutions from each parametric case are then validated with the simulation model in Step 3.

Step 3 simulates the water quality for each solution identified by the MOMILP in Step 2. This resimulation step acts as a validation of the optimization step. The simulated water quality at each node is used to calculate the value of each objective function. These validated objective function values define the validated Pareto frontier. Comparison of the Pareto frontier generated in Step 2 with the validated frontier generated in Step 3 is a measure of the accuracy of the model framework. An acceptable level of error is relative, and is based on the specific intended use of the model and the expectations of the decision maker. If the error at this stage is considered unacceptably high, then an informed set of preliminary simulations should be used in Step 1 and the model framework should be rerun.

If the error calculated in Step 3 is considered acceptable, the validated Pareto frontier and the associated treatment alternatives can be discussed with the decision maker in Step 4. Step 4 integrates any nonnumeric objectives or constraints important to the decision maker into the model framework by allowing the decision maker to select the recommended alternative from the set of model generated alternatives.

The level of integration between the physicochemical simulation and optimization stages in these types of models must be decided by the decision maker. Although some simulation–optimization approaches use a fully embedded simulation model, iterative solutions also are common for a variety of optimization formulations. Razavi et al. (2012) discussed this decision in detail as it relates to surrogate models, and 11 of the 32 studies reviewed in detail utilized a sequential framework (i.e., iterative). This work used a native physicochemical simulation model and a sequential integration approach to link the models. This decision was based on the desire to keep the model flexible and applicable for modeling larger systems and a sequential framework gives the user more control over the computational intensity of the model and allows for the decision maker to interact with results at each stage.

Step 2 includes a multiobjective optimization formulation that was not considered in the previous work (Schwetschenau et al. 2019). The multiobjective optimization formulation used in the demonstration case is shown in Fig. 2. Four cost parameters are used in the cost objective function and are defined by supply source type. Cost function coefficients are derived using the approach developed by the USEPA (Clark 1982; Gumerman et al. 1979) and applied by Boccelli et al. (2007) and shown in Schwetschenau et al. (2019).

Some or all data, models, or code generated or used during the study are available from the corresponding author by request. The specific Net3 network specification files are publicly available and included with the EPANET Programmers toolkit (USEPA 2014). Minor modifications were made to the input files; the specific files used 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 selected to run the optimization software, and the optimization formulation is provided in the Fig. 2. Specific GAMS scripts can be provided upon request. The methods by which the model components are integrated are described in the Appendix and Schwetschenau et al. (2019).

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.