MAGNets: Model Reduction and Aggregation of Water Networks

MAGNets offers the function reduce_model that allows the user to submit a water network model and customize the model by choosing: (1) a specific operating point (around which the model will be linearized); (2) a list of nodes that should remain in the reduced model; and (3) a maximum nodal degree value for the nodes removed from the model. Fig. 2 demonstrates the reduce_model function with the following inputs: (1) the inp_file is a .INP file containing the water network model object; (2) the op_pt is the operating point, which represents a multiple of the reporting time step. If the user does not provide an operating point, it will be set by default to 0, i.e., the start time of the simulation. Users are encouraged to check if the reporting time step of the .INP file aligns with the simulation period they are interested in investigating and then adjusting accordingly; (3) nodes_to_keep is a list of nodes to keep in the model. If the user does not provide a list of nodes to be retained in the model, every node that is not a special element or that is not connected to a special element will be removed from the model (as will be described below); and (4) max_nodal_degree, which dictates the maximum nodal degree of nodes being removed from the model. The nodal degree of a node is equal to the number of pipes incident to the node. For example, if max_nodal_degree = 1, all nodes that are not in nodes_to_keep and that are connected to only one other node by one pipe will be iteratively removed from the model. The reduce_model function returns two outputs: (1) a .INP of the reduced network model that will be created in the destination folder of the original model, and (2) a water network model object representing the reduced network model (shown as wn2 in Fig. 2) that the user can access using the USEPA WNTR Python package (Klise et al. 2017) without having to import the new .INP file.

Prior to executing the main model reduction algorithm, MAGNets performs several preprocessing steps. First, the .INP file provided by the user is converted into a water network model object that is compatible with WNTR. Basic characteristics of the network are extracted and stored and are used as the basis to determine which nodes will be removed from the model. Nodes are classified as special, non-removable, or removable, and including a node in subsequent calculations is dependent of this initial categorization. Special nodes include reservoirs, tanks, source nodes, and junctions that have been assigned control rules. Special links include pumps, valves, and links that have been assigned control rules. Non-removable nodes include nodes that are directly connected to special nodes or special links. If the user provides a list of nodes to keep in the model, all nodes in that list are assigned non-removable status. All nodes in the model that are not assigned special or non-removable status are categorized as removable and will be removed during the model reduction process. Second, before the aggregation and reduction process is executed, the model is trimmed by removing branches and parallel pipes in order to enhance the efficiency of the reduction process. Every dead-end node connected by a single pipe is removed from the model and the demand attributed to that node is allocated to the upstream node. Several iterations are performed until all dead-end nodes and branched pipes are removed from the model. After branched pipes are removed from the model, all parallel pipes are identified and replaced by a single equivalent pipe (Boulos et al. 2006).

The main steps of the model reduction procedure are: (1) non-linear equations describing the hydraulics of the system are constructed and linearized around the operating point; (2) the linearized system of equations is reduced through the method of Gaussian elimination, resulting in the redistribution of nodal demands, the creation of new fictitious links, and the deletion of unnecessary nodes from the model; and (3) the updated linearized equations are converted back into non-linear form to derive the characteristics of the reduced water network model, which is subsequently returned to the user. The reduction process is briefly described below and the full details of the process can be found in Ulanicki et al. (1996) and Alzamora et al. (2014).

Step 1: Compute linear conductances

After the model has been trimmed, all non-special pipes are assigned conductance values that retain the diameter, length, and roughness of the pipes. The head loss equation for a single pipe with diameter $D$, length $L$, and roughness coefficient $C$ can be formulated aswhere $Q_{ij}$ is the flow rate in the pipe connecting nodes $i$ and $j$, $K_{ij}=\alpha L_{ij}/C_{ij}^{e_1}{D}_{ij}^{e_2}$, $e_1=1.852$, and $e_2=4.87$, and $\alpha$ is the Hazen-Williams unit constant (Boulos et al. 2006). To reduce rounding errors, the value of $\alpha$ is calculated as the mean value of the unit coefficients for each pipe in the network based on the values extracted from simulation results using the equation ${\alpha }_{ij}=\mathrm{\Delta }{h}_{ij}^0{C}_{ij}^{e_1}{D}_{ij}^{e_2}/L_{ij}{Q}_{ij}^{e_1}$, where $\mathrm{\Delta }{h}_{ij}^0=h_i^0-h_j^0$ is the head loss at the operating point. Rearranging Eq. (1), we getwhere $g_{ij}=K_{ij}^{-e_3}$ and $e_3=1/e_1$.

Linearizing Eq. (2) around an operating point, we getwhere ${\tilde{g}}_{ij}=g_{ij}\mathrm{\Delta }{h}_{ij}^{0^{e_3-1}}$ is linearized pipe conductance. Additionally, for each node $k$, the total linearized conductance is computed as the sum of all linear conductances of the incident pipes aswhere pipe $i$ belongs to the set of all pipes incident to node $k$.

Step 2: Reduce the linear model

After all conductance values are linearized, the list of nodes to be removed is sorted in ascending order based on the node degree. The node with the lowest node degree is deemed the “removal node” $k$ and all nodes connected to the removal node are labeled “adjacent nodes”. For each node $k$ that is being removed, two actions occur: (1) Nodal demand of node $k$ is allocated to its adjacent nodes proportional to the linear conductance, and the new demand of each adjacent node $i$ is updated aswhere $d_i^0$ and $d_k^0$ are the base demands at nodes $i$ and $k$, respectively, and $d_i^{\prime }$ is the new adjusted base demand at node $i$. (2) The linear conductances of pipes connecting all the adjacent nodes will be updated as follows:

If two adjacent nodes were not previously connected, then ${\tilde{g}}_{ij}=0$ and a new fictitious pipe is created. The code optimizes the order in which the nodes are removed from the model to minimize the number of fictitious pipes that are added to the model. After all adjacent nodes and connected links are updated, node $k$ and all the incident pipes are removed from the model, and the list of nodes and pipes remaining in the model is updated and again sorted in ascending order based on node degree, with the “removal node” for the next iteration being the new node with the lowest node degree. The process repeats in an iterative manner until all nodes in the list of nodes to be removed are removed from the model.

Step 3: Retrieve properties of the reduced non-linear model

At the final stage, the non-linear properties of the reduced model, i.e., aggregated demands and pipe properties, are retrieved from the linearized conductance as follows:

We assume that all new fictitious pipes have roughness coefficient $C=100$ and length $L$ equals to the average length of all pipes in the original model, which can be substituted into Eq. (7) to calculate the new diameter of each link. Notably, the new pipe diameters, which are calculated according to the pipe conductances, are not expected to be reasonable values on par with the original model. The final outcome is a new hydraulic model with a compatible .INP file that can be used to simulate network hydraulics in EPANET or WNTR.

The performance of the package was tested using 12 networks of varying complexity obtained from Rossman et al. (2020), Hernandez et al. (2016), and Ostfeld and Salomons (2008). Table 1 summarizes the tested networks, number of nodes in the original and reduced models, running times, and the maximum error. The MAGNets repository (Thomas and Sela 2022) contains all the tested networks and several example codes that demonstrate the functionality of the package. The networks were modified prior to reduction to ensure that all nodes to be removed had the same demand pattern. The error was computed as the maximum relative error (%) during the entire duration of the simulation among all remaining nodes in the reduced model, where the relative error was computed as the difference in nodal heads in the reduced and original models divided by the nodal head in the original model. Results demonstrate that the reduction process running times for all networks were under 3.3 s with the exception of BWSN2, which took approximately 4 min. Additionally, the maximum observed errors in all networks were under 5.5%. As an example, Fig. 3 shows the full and reduced network topology for KY2, and Fig. 4 shows the flow rate and head at a selected pump and tank. The codes to reproduce Table 1 and Figs. 3 and 4 are available on GitHub Repository https://github.com/meghnathomas/MAGNets/tree/master/publications.

For most networks, the deviation of heads in the reduced model from those of the original model can be attributed to slight temporal misalignment in control rules activation causing changes in tank filling times. An example of this can be demonstrated by running the MAGNets reduce_model function on network KY2 and comparing the heads of tank T-2 in the original and reduced models. The control rules operating on the tanks in the network dictate that each tank will be opened and closed at specific pressure values. As the network has been linearized around a specific operating point, slight changes in heads at other times along the simulation duration are to be expected. However, these slight changes then affect the time step in which the tank achieves the head that will result in it opening or closing. Fig. 4 demonstrates how the choice of operating point can affect the reduced model accuracy: if operating point 0 is chosen for network KY2 [as is shown in Figs. 4(a and c)], the hydraulic performance of the reduced model begins to deviate from the original model as the tank head begins to drop earlier than in the original model, thus causing a temporal shift in the head values at tank T-2. This difference in tank head is mirrored in the flow rate of Pump-1. Errors in pressure heads along the rest of the time duration remain within reasonable limits. Upon testing and calculating the error for all possible operating points, it could be determined that operating point 13 produced very similar simulation results for the original and reduced models, as can be seen in Figs. 4(b and d).

Errors in heads and flows between original and reduced network models can be minimized by testing a wide range of operating points and using different reduced models for specific operating conditions. It is also possible to reduce errors by leaving more nodes in the reduced model; however, this requires careful consideration. Users may also employ the max_nodal_degree input to control the level of reduction, e.g., trim only branches and merge pipes connected in series. Another strategy to only partially reduce the model is to halt the reduction process after a desired fraction of nodes in the original model have been removed. Figs. S1 and S2 show the median error, reduction time, and number of pipes in the reduced models for KY2 when controlling the maximum degree of the nodes removed (i.e., 1, 2, 3, and no restriction) and controlling the number of removed nodes (25%, 50%, and 75%), respectively.

We investigated the impact of the order in which the nodes are removed on the accuracy and running times of model reduction. Upon determining the list of nodes that can be removed from the model, the list was reordered in three ways: randomly, statically, and dynamically, each time removing 25%, 50%, 75%, and 100% of removable nodes. The results of the random ordering are the mean results of ten different randomized ordered lists. The statically ordered list was determined by arranging the list of nodes to be removed in ascending order of node degree at the beginning of the reduction process. In the dynamically ordered list, the nodes are removed in ascending order by the node degree, which is constantly updated each time a node is removed from the network.

We defined three performance criteria: (1) the accuracy of the reduced model, which was computed as the median relative error between the pressure heads of nodes remaining in the reduced network and the equivalent nodes in the original network; (2) model reduction time; and (3) number of links remaining in the reduced model. We show results for three networks of increasing size, NET3 (92 nodes), KY1 (791 nodes), and BWSN2 (12,523 nodes). Fig. 5 shows the three performance criteria under different removal sequences and node list partitions (25%, 50%, and 75%) for the test networks. Each row represents a performance criterion and each column represents a network. For example, the figure on the top left shows the median percentage error between the original NET3 model and reduced models of NET3 when 25%, 50%, and 75% of the nodes are removed in static, random, and dynamic order. The figure in the middle left demonstrates the reduction time in seconds, and the bottom left figure shows the number of pipes left in the reduced model for NET3 under the same criteria. Overall, the dynamically ordered lists outperformed the other orders in all three metrics of interest, resulting in lower median error, significantly fewer pipes retained in the reduced model, and a lower reduction time, with the exception of greater reduction time for the small network NET3. Results for NET3 differ from those of KY1 and BWSN2 most likely because the network was considerably smaller in size and the time required for reordering the nodes in each iteration exceeded that of removing the elements; however, the reduction times for all ordering methods and list sections were very similar. Surprisingly, the statically ordered list performed the worst in all three performance criteria. Statically ordered reduction created more pipes in the model for all three networks, thus significantly affecting running times because the number of elements in the model increased with each iteration. Statically ordered removal also generated high median errors, on par with or surpassing that of the random ordered removal. Additionally, the results revealed that error does not necessarily correlate with the size of the network or with the number of nodes left in the model. The statically ordered removal of 75% of nodes for BWSN2 had a running time of more than 72 hours, so those results were not included in the analysis.

We can conclude that dynamically ordering the nodes to be removed best serves the purpose of maximizing model accuracy and minimizing reduction time and model complexity. Randomly ordered lists would be the next best alternative with lower reduction times and model complexities compared to statically ordered lists while maintaining a comparable level of accuracy. Following these results, dynamic ordering of nodes was implemented in MAGNets.

In this section, we compare the proposed approach with the WNTR skeletonization function. WNTR is an open-source Python package that offers a skeletonization functionality, which allows the user to remove branches and merge pipes in series and in parallel (Klise et al. 2017). WNTR skeletonization function requires the user to provide a network .INP file and a threshold for maximum pipe diameter for pipe removal, where pipes with diameters less than the threshold will be removed from the model. When pipes are removed/merged, the demands of the incident nodes that are removed are assigned to the closest neighboring node by appending the base demand value and demand pattern of the removed node as a new demand category for the receiving neighbor node. Then, the properties of the new equivalent pipe are updated based on the characteristics (length, diameter, roughness) of the removed pipes. However, WNTR does not allow the user to control which nodes to keep in the model or remove nodes with nodal degree greater than 2. In contrast, MAGNets reallocates the base demands of nodes that are removed from the model to its adjacent nodes proportionally to the conductances of the incident pipes, enables the user to specify which nodes should remain in the model, and allows for the removal of nodes with a nodal degree greater than 2.

We demonstrate the comparison of the performance of WNTR and MAGNets using network KY2 by removing nodes with degrees 1 and 2. For WNTR, the largest pipe diameter in the network was supplied as the diameter threshold. For MAGNets, an operating point of 13 hours was provided as an input. The original KY2 model contained 811 nodes (see Table 1 and Fig. 3). Removing only branches from the KY2 model resulted in a model with 593 nodes (a 27% reduction) and yielded identical reduced models with WNTR and MAGNets, in which the absolute relative error between nodal heads in the reduced and original models was lower than ${10}^{-4}\%$ for all nodes. Removing additional nodes with degree 2 resulted in a reduced model with 459 nodes (a 43% reduction) for both WNTR and MAGNets. Although the size of the reduced models is the same, model properties are different since different methods are applied in the reduction process. Fig. S3 shows the absolute relative error for the reduced models resulting in removing nodes with degrees 1 and 2 using MAGNets (top left) and WNTR (top right). The shaded areas represent the 75%, 95%, and 100% percentiles, and the red line represents the mean absolute error. We observe that both reduced models exhibit a very good match with the original models in terms of the errors, where MAGNets slightly outperforms WNTR with a maximum error of 0.07% compared to 0.14%. However, as demonstrated earlier, MAGNets is capable of reducing the model even further. Using the reduce_model function without limiting the node degree of the removed nodes, the model size is reduced to only five nodes (a 99% reduction) (Fig. 3) and the largest absolute relative error for the heads of the remaining nodes is 0.56% (Table 1). Fig. S3 (bottom) shows the distribution of error for all nodes throughout the simulation duration, where shaded areas represent the 75%, 95%, and 100% percentiles, and the red line represents the mean absolute error.

A set of scenarios for the C-Town network was evaluated to test the capabilities of MAGNets’ reduce_model function when it is confronted with different demand patterns and when the user wishes to reduce specific sections of the network. The C-Town network model contains five different District Metered Areas (DMAs), each with a unique demand pattern (Ostfeld et al. 2012). In order to adapt the reduction strategy to C-Town, a list of boundary nodes that lay on the borders of DMAs was provided as the nodes_to_keep input for each call of the reduce_model function. A stage-wise reduction of the C-Town model was performed by first removing all nodes with nodal degree 1 from the model and then removing nodes with nodal degree 2 from DMA1, DMA3, and DMA5. The number of nodes in the final reduced model decreased from 388 to 176 (a 55% reduction). Fig. S4 shows the network layout of the full and reduced models. The time series of the relative error of the heads at each node reveals that even at this high level of reduction the reduced model performs very similarly to the original model (see Fig. S5). The relative errors are negligible for the majority of the extended period simulation with the exception of a few large spikes which can be attributed to time step delays in control rules activation. For example, Fig. S6 shows that spike 1 and spike 2 in Fig. S5 can be explained by delays in the operational status of pump PU2 based on the water level of tank T1. The maximum difference between heads of nodes in the original and reduced models ranges from $-3m$ to $12m$, while the mean difference ranges from $-0.02m$ to $0.08m$. Fig. S7 shows the distribution of errors and demonstrates that the mean absolute relative percentage error is 0.06% or lower for all nodes remaining in the reduced model. Finally, Fig. S8 shows the spatial distribution of the errors, demonstrating that nodes with similar error are located close together in the reduced C-Town model, and the nodes with the largest errors after reduction lie in DMA1 or are adjacent to pumps. These examples show that, if applied carefully, MAGNets can also be used to significantly reduce models with varying demand patterns. Users are advised to attempt a range of reduction strategies in order to maximize the usefulness and accuracy of the reduced models.