Optimal Resilience Enhancement of Water Distribution Systems

Abstract

Water distribution systems (WDSs) as critical infrastructures are subject to demand peaks due to daily consumption fluctuations, as well as long term changes in the demand pattern due to increased urbanization. Resilient design of water distribution systems is of high relevance to water suppliers. The challenging combinatorial problem of high-quality and, at the same time, low-cost water supply can be assisted by cost-benefit optimization to enhance the resilience of existing main line WDSs, as shown in this paper. A Mixed Integer Linear Problem, based on a graph-theoretical resilience index, is implemented considering WDS topology. Utilizing parallel infrastructures, specifically those of the urban transport network and the water distribution network, makes allowances for physical constraints, in order to adjust the existing WDS and to enhance resilience. Therefore, decision-makers can be assisted in choosing the optimal adjustment of WDS depending on their investment budget. Furthermore, it can be observed that, for a specific urban structure, there is a convergence of resilience enhancement with higher costs. This cost-benefit optimization is conducted for a real-world main line WDS, considering also the limitations of computational expenses.

1. Introduction

Fresh water is an essential resource, not only for human survival and well-being but also as crucial element in many sectors of the economy. Therefore, water distribution systems (WDSs) are critical infrastructures as they maintain multiple and vital societal functions, e.g., health, safety, security, economic and social well-being, characteristic of critical infrastructures as defined by the EU Council [1]. The increase in the world’s population combined with ongoing urbanization makes the urban WDS an essential infrastructure with a future-oriented focus [2], even if its components are aging [3]. Future water demand is predicted to rise by 30% between the years 2000 and 2050 [4]. Therefore, WDS adaptations leading to system resilience, i.e., the capability to withstand and recover from disturbances [5], has been of particular interest. Hereby, water suppliers are faced with the challenge of achieving high-quality and, at the same time, low-cost water supply. This leads to the necessity of the herein conducted cost-benefit optimization for the enhancement of the resilience of existing WDSs.

The method presented in this work can be identified as predictive asset management [3] in which the WDS is assessed by its key performance indicator, i.e., demand fulfillment, within the proposed resilience metric. Instead of performing a risk estimation and optimizing the system based on its robustness when facing a specific failure scenario, as typical for established asset management methods, resilience is the focus here. To analyze system behavior and to assist decision makers, the cost, considering the proposed preventive action, can be varied. Furthermore, by classifying it as proactive asset management, optimal infrastructure performance enhancement for a given cost can be identified. In turn, the operator has to choose the standard of function to be aimed at, considering related costs.

The resilience of different infrastructures to meet demand is the subject of current research [6]. Resilient infrastructures have a higher capability to meet an acceptable level of demand when faced with challenges [7] and are able to “efficiently reduce both, the magnitude and duration of the deviation from targeted system levels” [8].

To assess the resilience of infrastructures, networks can be mathematically described in planar graphs [9]. In this way, the WDS is translated into a mathematical graph, in which the pipes make up the edges connecting the nodes, which can be distinguished as either consumer or source nodes. By unfolding the edges, a planar graph can be achieved in which pipes are not allowed to cross under or over each other. In turn, intersections are always represented by a node connecting several edges, therefore allowing flow from one edge to continue towards any of the other adjoining edges, respecting the law of continuity.

Based on this mathematical description, different graph-theoretical resilience metrics exist. These metrics can be clustered into three groups. First are resilience metrics based on WDS topology. These metrics can assess resilience independent of a specific disruption scenario. Secondly, resilience metrics can focus on demand satisfaction for a specific disruptive event. This can be seen as contradictory to a resilience concept depending on realization, since a clear distinction with a robustness assessment, focusing on one specific disruption scenario, has to be ensured. Thirdly, there are resilience metrics which assess the recovery process after a disruptive event and therefore focus on the dynamic aspects of the infrastructure. Again, a clear distinction between this approach and robustness assessment is crucial.

There are different levels of sophistication in the first approach of assessing the resilience of a WDS based on its topology. There are simpler metrics, such as spectral and statistical metrics, for network assessment [10]. As an example, the meshed-ness coefficient gives the relative number of independent loops between the actual number and the maximum possible number, in a planar graph with the same number of nodes. This metric assesses the network’s redundancy instead of its resilience. Another network metric is central-point dominance, which assesses the average difference in centrality for each node of the network and therefore its robustness due to the high vulnerability of the central nodes. Resilient networks in turn are a result of both redundant and robust networks, as studied by Yazdani et al. [10,11]. Therefore, more sophisticated graph-theoretical metrics can assess the network by combining these two aspects. This work underlies the well-established graph-theoretical resilience index introduced by Herrera et al. considering both aspects [12]. The physical considerations underlying this approach, i.e., the hydraulic resistance of pipes, have been carefully derived [13].

When assessing the functional performance of a WDS by resilience metrics of the second type, specific disruption scenarios need to be taken into account. As stated before, it is of importance to differentiate from robustness analysis by considering only one disruption scenario. This is well accomplished by Todini, by defining a resilience index which is applicable independent of the disruption scenario [14]. Based on this resilience index, many extensions have been proposed [15]. Diao et al. propose resilience assessment of WDSs for different disruption scenarios [16]. Studying the change in water availability by, on the one hand, considering increased demand and, on the other, studying decreased precipitation, leads to the demand satisfaction resilience study conducted by Amarasinghe et al. [17].

Metrics using the third approach assess the dynamic behavior of a WDS exposed to a disruptive event by studying time-varying demand satisfaction, starting at the moment of the disruptive event. This is implemented by Hashimoto by considering the reciprocal time span during which water demand satisfaction is disturbed in order to assess the resilience of a WDS [18]. Since Japan is prone to disruptive events due to natural disasters, it has taken governmental and legal steps to prioritize the resilient design of its WDS, thereby, among others, aiming for a decrease in the time span of disturbed water demand in the case of natural hazards [19]. Likewise, Meng et al. study the recovery process after a disruptive event up to the point of completely restored functionality by various evaluation criteria instead of by one specific resilience measure [9].

The aim of this work is to not only assess the resilience of a present WDS but, furthermore, to enhance its resilience by adapting the WDS. As the addition of pipes to an existing WDS has been shown to be an effective measure to enhance WDS resilience, this approach is pursued [13]. Thereby, the physical feasibility of pipe addition within the existing urban structure is an important aspect. This can be ensured by taking into account the parallelism of infrastructures, i.e., the urban transportation network and the WDS [20]. To handle the challenging combinatorial problem of deciding which pipes lead to optimal resilience enhancement, a Mixed Integer Nonlinear Optimization Program is implemented. With the help of this optimization model, a cost-benefit analysis is conducted for a real-world WDS. The conducted cost-benefit analysis shows the convergence of the relative resilience enhancement due to the urban structure in which the WDS is embedded.

2. Materials and Methods

In this work a cost-benefit optimization of WDS resilience is implemented and conducted for a real-world main line WDS of a city district. Therefore, a Mixed Integer Linear Program is derived based on graph-theoretical resilience assessment considering the WDS topology introduced by Herrera et al. [12]. This resilience index considers both redundancy and robustness of the WDS. First, by considering all possible paths connecting demand and source nodes, the consequences of possible component failure or local contamination can be counteracted. Second, by ranking the paths by physical feasibility considering their hydraulic properties, the robustness of individual paths is accounted for. To do so, each path is assessed for its hydraulic resistance, computing the non-dimensional pressure loss of each edge making up the path, as derived by Lorenz et al. [13]. In the following, the method of assessing and optimizing the main line WDS resilience is described in detail. This method is applicable to any WDS while yielding the best results when considering a main line WDS. In this work it is applied to a real-world main line WDS of a German city district.

2.1. Resilience Assessment

To assess the resilience, the real-world WDS is presented as a planar graph in which the pipes are represented as edges connecting source and consumer nodes. The set of these nodes is given as $S$ and $C$, respectively, as well as the set of junction nodes $J$, while the unified set is given as $N$. The demand of a specific consumer node $c$ is given by $q_c$, while the overall demand of the WDS is given by $Q={\sum }_{c\in C}^{}{q}_c$. The set of existing edges is given in a binary matrix $T^0$ of size $\left|N\right|\times \left|N\right|$ in which an existing edge $T_{i,j}^0$ connecting node $i$ with node $j$ is represented by $1$ and all inexistent edges by $0$. Since the direction of flow in a pipe is not fixed but depends on the actual pressure gradient, all edges are undirected. Therefore, the matrix describing the WDS is symmetric. Furthermore, there are no self-loops within the network, meaning there are no pipes starting and ending at the same node, and therefore its diagonal is strictly $0$. Each edge has different characteristics, also presented in this matrix notation. Of interest for this work are the pipe diameters, lengths and roughness given in the matrices $D$, $L$, and $a$, respectively.

When considering the addition of pipes, the edges representing them must be noted in a similar way. To consider the physical feasibility of the addition of a pipe, the current urban structure is of great importance as pipes cannot be laid between all nodes. Different studies show, and legal standards command, that infrastructures are highly parallel [20,21]. In general, there is a high parallelism between supply infrastructures, i.e., the water, energy and information and communications technology, and the urban transportation network, i.e., streets and pathways. Therefore, the addition of pipes is restricted by the urban transportation network embedded in the present urban structure. Information on the urban transportation system is mostly open access data and can for example be extracted from OpenStreetMap [22]. With this data, the maximum possible edge matrix can be derived as $T^{\max }$ which is equally of size $\left|N\right|\times \left|N\right|$ and symmetrical. The edge characteristics of the added pipes can be written into the same matrices $D$, $L$, and $a$.

To assess the graph-theoretical resilience of the WDS, the overall resilience $I_{\mathrm{GT}}$ is given as the sum of each consumer node’s resilience, as defined by Herrera et al. [12]. To furthermore consider the demand and therefore the number of consumers depending on the demand fulfillment of a specific node, in this work it is proposed to weigh the resilience of each consumer node $c$ by its relative water demand $q_c/Q$.

$$I_{\mathrm{GT}}={\displaystyle \sum _{c\in C}}\frac{q_c}{Q}{\displaystyle \sum _{s\in S}}\frac{1}{K_{\max }}{\displaystyle \sum _{k\in K_{\max }}}\frac{1}{r_{s,c,k}}$$

(1)

This topology based resilience metric considers the hydraulic resistance $r_{s,c,k}$ of each $k$-shortest path feeding the consumer node $c$ from source $s$. As high hydraulic resistance paths have little impact on the resilience index, a maximum number of shortest paths to be considered, $K_{\max }$, can be determined for the critical transfer node of the complete WDS, following the approach introduced by Lorenz et al. [13]. The weight of the edges making up the shortest paths is, despite what the name might suggest, not its length but the hydraulic resistance of the pipe connecting the two nodes under consideration. Each path $k$ connecting source node $s$ and consumer node $c$ can be represented by a two dimensional binary matrix $A_{s,c,k}$ similar to the matrix $T$ describing the network. In this matrix $A_{s,c,k}$ sequence of nodes, part of the path starting at the source $s$ and ending at the consumer node $c$ is stored. Edges which are part of this path connecting two subsequently following nodes making up the path, have the value $1$ while all other edges are of value $0$. In contrast to matrix $T$ describing all pipes in a network, the matrix describing a specific path is not symmetrical, due to the direction of the path represented by the sequence of nodes connected by the edges making up the feeding path.

The physical considerations underlying the determination of the hydraulic resistance of each pipe is deduced by Lorenz et al. [13]. Applying an order of magnitude analysis, it was derived that in the case of a WDS turbulent flow in hydraulic rough pipes can be assumed. The hydraulic resistance $r\left(s,c,k\right)$ of a path made up of pipes $m\in M_{s,c,k}$ represents the dimensionless pressure loss along this path and can be determined as follows.

$$r\left(s,c,k\right)=\frac{\Delta p_{\mathrm{l}}}{\mathsf{\varrho }/2u_0^2}\approx {\displaystyle \sum _{m\in M_{s,c,k}}}f\left(m\right)\frac{L_m}{D_m}={\displaystyle \sum _{i\in N}\sum _{j\in N}}{A}_{s,c,k,i,j}f\left(i,j\right)\frac{L_{i,j}}{D_{i,j}}$$

(2)

Therefore, the pressure loss $\Delta p_{\mathrm{l}}$ along the path is related to the characteristic kinematic pressure given by the characteristic velocity $u_0$ as well as to the density of the fluid in the pipes $\mathsf{\varrho }$, herein water. In the case of turbulent flow in hydraulic rough pipes, the friction factor $f\left(m\right)$ is given by Prandtl-Karman’s law as $f\left(m\right)={\left(2\lg \left(D_m/a_m\right)+1.74\right)}^{-2}$ [23]. The calculation of hydraulic resistance can also be rewritten in matrix notation, making use of the representation of each shortest path for every consumer and source node in the five-dimensional matrix $A$.

2.2. Resilience Optimization

The resilience index introduced can be applied to assess the resilience of any WDS for which the network topology, as well as the pipe diameter, length, and roughness and the demand of all consumer nodes, is known. Therefore, present and further adaptations to the present WDS can be assessed. Adapting the WDS by adding pipes results in alternative feeding paths of different hydraulic resistance for different consumer nodes. These alternative feeding paths therefore impact the overall resilience differently, also depending on the weight of the consumer node being fed. As the problem of determining which pipe is best to add for maximum resilience enhancement rises to become a challenging combinatorial problem, it is useful to apply mathematical optimization in order to solve it.

The problem can be formulated as a maximization problem, as the resilience index of the WDS is to be maximized. The objective function is therefore the maximization of the resilience index.

$$I_{\mathrm{GT},\mathrm{opt}}=\max {\displaystyle \sum _{c\in C}}\frac{q_c}{Q}{\displaystyle \sum _{s\in S}}\frac{1}{K_{\max }}{\displaystyle \sum _{k\in K_{s,c}}}\frac{1}{r_{s,c,k}}$$

(3)

This objective function considers the reciprocal of the hydraulic resistance of each path feeding a consumer node. This leads to a nonlinear optimization problem with a concave solution space, which are in general computationally expensive to solve. The formulation to find the shortest path in a network by a discrete optimization problem can also be easily implemented [24]. Yet, in turn, finding the second, third and $k$-th shortest path is usually solved by heuristics, as for example by Yen’s K-shortest paths algorithm [25]. This makes the formulation of all shortest paths $A_{s,c,k}$ connecting a source and consumer node as an optimization problem infeasible for a dynamic graph, to the best of the authors knowledge. Identification or knowledge of all the shortest paths, in turn, is necessary to determine the hydraulic resistance as defined by Equation (2).

Taking into consideration the constraints on physical feasibility of a pipe addition to an existing WDS given by the parallelism of infrastructures, the possible addition of pipes is highly restricted compared to when considering a complete graph, i.e., a graph in which every pair of nodes is connected by an edge. Therefore, the calculation of all feeding paths in the maximum possible WDS, for each consumer node fed by any source, considering that all physically feasible pipes are added, is practicable for smaller networks. When calculating all possible feeding paths in a preprocessing step, the optimization problem can be rewritten as

$$I_{\mathrm{GT},\mathrm{opt}}=\max {\displaystyle \sum _{c\in C}}\frac{q_c}{Q}{\displaystyle \sum _{s\in S}}\frac{1}{K_{\max }}{\displaystyle \sum _{k\in K_{s,c}}}{x}_{s,c,k}{g}_{s,c,k},$$

(4)

where the hydraulic conductance $g\left(s,c,k\right)$, i.e., the reciprocal of the hydraulic resistance, is already known for each path. The binary decision variable $x\left(s,c,k\right)$ permits the choice of the best possible feeding paths as part of the adapted WDS, and thereby the pipes to be added to the existing WDS for a set maximum length added.

Based on the linearized resilience optimization, taking into consideration only the single shortest path of each source-consumer node combination, conducted in [13], the resulting Mixed Integer Linear Program considering the $K_{\max }$-shortest paths is given in Equations (5)–(11). This optimization program relies on the following variables, sets and parameters, given in

Table 1. Decision variables of the optimization program.

Variable	Domain	Description
$b_{i,j}$	$\left\{0,1\right\}$	Decision variable for choosing the pipes to be added to the existing water distribution system (WDS).
$x_{s,c,k}$	$\left\{0,1\right\}$	Decision variable for choosing the possible feeding paths for the adapted WDS.

,

Table 2. Sets of the optimization program.

Set	Description
$S$	Set of source nodes.
$C$	Set of consumer nodes.
$J$	Set of junction nodes.
$N$	Unified set of nodes—source, consumer or junction nodes, $N=S\cup C\cup J$.
$K_{s,c}$	Set of the number of maximum possible paths, $\forall s\in S,c\in C$.

and

Table 3. Parameters of the optimization program.

Parameter	Domain	Description
$T_{i,j}^0$	$\left\{0,1\right\}$	Edges of the existing WDS, $\forall i,j\in N$.
$T_{i,j}^{\max }$	$\left\{0,1\right\}$	Edges in the maximum possible WDS restricted by the urban structure, $\forall i,j\in N$.
$A_{s,c,k,i,j}$	$\left\{0,1\right\}$	All possible feeding paths in the maximum possible WDS $T_{i,j}^{\max }$, $\forall s\in S,c\in C,k\in K_{s,c},i,j\in N$.
$g_{s,c,k}$	${ℝ}_0^{+}$	Hydraulic conductance of each feeding path determined in the preprocessing, $\forall s\in S,c\in C,k\in K_{s,c}$.
$q_c$	${ℝ}_0^{+}$	Demand of each consumer node, $\forall c\in C$.
$K_{\max }$	$\mathbb{N}$	Maximum number of paths to be considered for the resilience assessment.
$L_{i,j}$	${ℝ}_0^{+}$	Length of each pipe in the maximum possible WDS $T_{i,j}^{\max }$, $\forall i,j\in N$.
$L_{\max }$	${ℝ}_0^{+}$	Maximum allowed length for the sum of all added pipes.

. It is implemented in Python using the Gurobi solver [26]. For the preprocessing of the WDS, the python packages WNTR [27] in combination with NetworkX [28] are applied.

As described in Table 1, two binary decision variables are introduced: on the one hand, choosing the pipes to be added to the existing WDS and therefore of similar form as the matrix describing the edges of the existing WDS, $b_{i,j}$; on the other hand, selecting the paths which are of the highest hydraulic conductance possible for the adapted WDS, $x_{s,c,k}$.

The unified node set $N$ contains the source node set $S$ as well as the consumer node set $C$ and the junction node set $J$ (see Table 2). Each possible path connecting source and consumer nodes is also given in the path set $K_{s,c}$.

In addition to variables and sets, the optimization program relies on multiple parameters together with the sets which define the exact instance for an optimization, i.e., the unique optimization problem for which a specific optimization result can be achieved (see Table 3). There are several parameters specific to the present WDS and its urban structure but also one parameter specific to the cost-benefit analysis conducted by this optimization program. This is the maximum allowed length for the sum of all pipes added $L_{\max }$ to determine the adapted WDS’s enhanced resilience under consideration of a specific investment budget, which can be varied as any positive real number including zero ${ℝ}_0^{+}$. For a set pipe diameter, the costs are assumed to be linearly dependent on the overall pipe length added, since they consist of material costs and construction costs, both approximately proportional to the length of the pipe. Therefore, the variation of this Parameter $L_{\max }$ allows the determination of a Pareto front for the WDS’s resilience under varying maximum overall pipe lengths added to the existing WDS. This is shown in detail in the Results, Section 3.

The remaining parameters of the optimization program are given by the existing WDS and the urban structure it is embedded in and are determined in the complex preprocessing of the optimization program. First, the main line WDS is determined from the existing overall WDS by making use of a cut-off criterion for the pipe diameter. This main line WDS is further processed by merging the demand of consumer nodes locally, i.e., $q_c$, so that the overall demand is still true to the original WDS but the number of consumer nodes is reduced and therefore the computational expenses can be reduced in the subsequent optimization. For this main line WDS the edges are noted in the matrix $T^0$. To determine the edges of the maximum possible WDS, OpenStreetMap data for the district containing the existing WDS are studied in combination with the existing WDS. By doing so all streets or pathways in which a pipe does not already lie are added as an edge to the existing WDS. This maximum possible WDS is represented by $T^{\max }$. For all possible additional edges, the important pipe characteristics, i.e., diameter, length and roughness, are chosen consistently to further diminish computational expenses. The length of the edges can be extracted from the OpenStreetMap data while diameter and roughness are chosen according to the typical pipe characteristics of newer pipes in the existing WDS. With this information, the hydraulic resistance of each edge in the maximum possible WDS can be determined.

Additionally, for this maximum possible WDS the critical transfer node and the maximum number of shortest paths to be considered in the resilience assessment of 1% accuracy $K_{\max }$ can be determined. This limit of 1% accuracy can be adapted in the preprocessing. Subsequently, the $K_{\max }$-shortest paths for the existing WDS are heuristically identified, making use of Yen’s K-shortest path algorithm. For the maximum possible WDS ideally all possible paths with less hydraulic resistance than the highest hydraulic resistance path should be determined, as all these paths can increase the WDS’s resilience. Since this is computationally very expensive, a cut-off of $b{K}_{\max }$-shortest paths for the maximum possible WDS has to be identified. This number varies with the size of the maximum possible WDS given by its consumer nodes and with the number of edges. The undertaken cut-off of the number of shortest paths considered for further optimization leads to limitation of the possible solution space but makes the problem feasible. Therefore, global optimization is only ensured if the feeding paths’ cut-off is worse than the already existing feeding paths. To further decrease computational expenses, this fact is made use of by cutting off the computation of any further feeding paths once the highest hydraulic resistance feeding path of the maximum possible WDS exceeds the highest hydraulic resistance feeding path of the existing WDS. To ensure that all $K_{\max }$-shortest paths of the existing WDS are considered in the optimization, the determined $2K_{\max }$-shortest paths for the maximum possible WDS are complemented by all computed paths of the existing WDS which are not already part of the $2K_{\max }$-shortest paths for the maximum possible WDS. The complemented $K_{s,c}$-shortest paths for the maximum possible WDS is then used to calculate the conductance $g_{s,c,k}$ for each path.

Following, the optimization program, Equations (5)–(11) can be set up and solved. The objective is the maximization of resilience when a fixed length of overall pipe addition is allowed (Objective 5). The constraint linking the choice of pipes to be added to the existing WDS and therefore the shortest possible paths feeding a consumer node is given by Constraint 6. This constraint states that all edges making up a possible feeding path have to be part of either the existing WDS or have to be added to the WDS. The length of the sum of all added pipes is further restricted by twice the maximum length of the sum of all added pipes, as the undirected nature of the WDS allows for the addition of an edge connecting node $i$ and $j$ and vice versa (Constraint 7). This interrelation is stated in Constraint 8. Thoughts on the physical feasibility of the addition of pipes is given by Constraint 9, as pipes may only be added if they are part of the maximum possible WDS but not already part of the existing WDS. To pre-limit the pipes under consideration to be added to the existing WDS for a certain maximum cost of $L_{\max }$, all pipes longer than $L_{\max }$ are directly eliminated from the set of pipes under consideration (Constraint 10). It is also ensured that only the best possible $K_{\max }$-shortest paths of the adapted WDS are considered for the resilience assessment of the adapted WDS (Constraint 11).

$$\max {\sum }_{c\in C}^{}\frac{q_c}{Q}{\sum }_{s\in S}^{}\frac{1}{K_{\max }}{\sum }_{k\in K_{s,c}}^{}{x}_{s,c,k}{g}_{s,c,k}$$

(5)

subject to

$$A_{s,c,k,i,j}{x}_{s,c,k}\le T_{i,j}^0+b_{i,j}\forall s\in S,c\in C,k\in K_{sc},i,j\in N$$

(6)

$${\displaystyle \sum _{i\in N}\sum _{j\in N}}{L}_{i,j}{b}_{i,j}\le 2L_{\max }$$

(7)

$$b_{i,j}=b_{j,i}\forall i,j\in N$$

(8)

$$b_{i,j}\le T_{i,j}^{\max }-T_{i,j}^0\forall i,j\in N$$

(9)

$$L_{i,j}{b}_{i,j}\le L_{\max }\forall i,j\in N$$

(10)

$${\sum }_{k\in K_{s,c}}^{}{x}_{s,c,k}\le K_{\max }\forall s\in S,c\in C$$

(11)

As stated earlier, the global optimization of the solution highly depends on the computational capacities available in preliminary examination of the solution space during preprocessing, when determining the shortest paths of the maximum possible WDS. Therefore, the distinction between consumer nodes and junction nodes is of importance, as only nodes with a demand, i.e., consumer nodes, are considered in the customized resilience index introduced in Equation (1). Junction nodes do not have an external consumer, but only allow for flow to be distributed to the neighboring edges.

This optimization model assists decision-makers in how to adapt an existing main line WDS for optimal resilience improvement for a set investment budget, taking into consideration the physical feasibility of the solution. When calculated for different instances of the same WDS, but with different investment budgets in mind, a Pareto front for the cost-benefit analysis can be determined. This also allows for an analysis of the maximum achievable resilience for the present urban structure.

3. Results

The previously presented optimization model is applied to a real-world main line WDS of a German city district. First, the existing WDS data provided by a local water supplier is prepared as described in the preprocessing paragraph of Section 2.2. The resulting main line WDS is reduced to 120 nodes and 127 pipes with diameters ranging between 100 mm and 240 mm. The roughness of each pipe depends on its material and age and ranges between 0.5 mm and 1.5 mm while the added or replaced pipes are mostly made of polyethylene and have a roughness of 0.5 mm. From this number of nodes, there are 80 consumer nodes and one source node. The remaining nodes are junction nodes, which remained in the WDS to represent its topology and allow for the branching and intersection of pipes considering the current urban structure. In further preprocessing and optimization, only 80 consumer nodes have to be considered leading to decreased computational expenses for the shortest path determination and the optimization of the objective function.

The maximum possible WDS is determined in a further preprocessing step in which the existing main line WDS is matched with the urban transportation network retrieved from OpenStreetMap using the OSMnx python package [29]. The maximum possible WDS is composed of 227 nodes and 274 pipes resulting in 147 possible additional pipes compared to the existing WDS. The existing main line WDS and maximum possible WDS are shown in Figure 1, in which existing pipes are represented by solid lines while possible additional pipes are represented by dashed lines.

Considering the characteristic pipe roughness and diameters of the existing WDS, it is assumed that all WDS adaptations are equally performed with polyethylene pipes of roughness 0.5 mm. The pipe diameter of additional pipes is chosen as 100 mm since additional pipes which are too large lead to a significant decrease in the flow velocity due to the laws of flow continuity. This is undesired since too little flow velocities brings the risk of biological build up in stagnating water. To show the physical feasibility of the additional pipes, a rough assessment of the maximum pressure losses compared to the available pressure is conducted, as presented in previous work [13]. Herein, losses of double the maximum spanning width of the maximum possible WDS, as well as the smallest pipe diameter of 100 mm, are considered. Thereby, the worst case scenario is covered. The order of magnitude analysis shows that this leads to roughly 3 bar pressure losses along this worst case scenario feeding path. The available pressure head of 53 m minimum height difference between the source and the highest consumer node still allows for a 3 bar feeding pressure at the consumer node, which is within the common feeding pressure range for fresh water facilities.

A further step in the preprocessing is the determination of the shortest feeding paths for all consumer nodes in the existing and the maximum possible WDS. As proposed by Herrera et al. the resilience index doesn’t change significantly for high resistance feeding paths [12]. Considering a 1% threshold for change in the resilience index for the critical transfer node in the maximum possible WDS, a critical number of feeding paths to be considered in the resilience assessment can be determined as $K_{\max }=24$, following the proposed method in [13]. Therefore, the 24 shortest feeding paths of the existing WDS and, due to the limited computational capacities, the 48 shortest feeding paths in the maximum possible WDS are determined. The unified set of these paths, considering the same paths once only, then allows for the optimization of the resilience of the existing WDS considering different maximum overall added pipe lengths for different optimization instances.

In this work to conduct a cost-benefit analysis, the instances are given by the variation of $L_{\max }=\{0,100,200,300,400,500,600,700,800,900,1000,1250,1500,1750,2000,2500,3000\}$ m. The first instance of $L_{\max }$ = 0 m determines the resilience of the existing WDS. Adaptations to the existing WDS depending on the respective instance are shown in Figure 2.

Herein, the existing main line WDS is marked in blue while the edges chosen by the optimization to be added to the WDS are marked in yellow. The presented adaptations to the existing main line WDS show that the addition of pipes does not follow a strict sequence regarding which pipes to add first, but the location of the pipes to be added varies with the overall length allowed to be added to the existing WDS. Two different types of pipes added to the existing WDS can be identified. On the one hand, pipes are added to close outside loops. On the other hand, pipes are added close to the first branching of the existing WDS after the source. The sequence regarding which kind of pipes will be added does not seem to follow a strict pattern but is individual to each instance. It can be observed that pipes added close to the first branching of the existing WDS do not seem to be discarded once the critical overall pipe length allowed is reached, so they can be added to the WDS. However, pipes closing outer loops in the existing WDS switch depending on the overall pipe length allowed to be added, and the already realized refinement of the area close to the first branching.

Regarding the resilience enhancement of the different instances, a Pareto front of the optimal relative resilience enhancement for the allowed overall pipe length added to the existing WDS can be mapped for the computed instances (Figure 3). This represents a cost-benefit analysis of the resilience enhancement for the currently existing WDS embedded in its urban structure, reflected upon by determining the maximum possible WDS. The relative resilience enhancement is determined as $I_{\mathrm{GT},\mathrm{rel}}\left(L_{\max }\right)=\left[I_{\mathrm{GT}}\left(L_{\max }\right)-I_{\mathrm{GT}}\left(0\right)\right]/I_{\mathrm{GT}}\left(0\right)$. The results of these instances show a steady resilience enhancement for longer overall pipe lengths added. For the present WDS and the urban structure it is embedded in, the close to linear resilience enhancement is limited to 130%. The maximum resilience enhancement seems to converge at 180%. Even though the number of pipes and the overall pipe length added to the existing WDS continues to increase, resilience stagnates. This seems to represent a behavior of diminishing returns and, even more so, of resilience saturation.

The implemented optimization model meets the challenging problem of which pipes to add for an optimal resilience enhancement of the existing main line WDS. The results give a clear indication of which pipes should be added considering the investment budget proportional to $L_{\max }$.

4. Discussion

The results of the cost-benefit optimization for an existing WDS leading towards enhanced resilience show that the challenging combinatorial problem cannot be solved intuitively. Interestingly, the resulting adapted main line WDS is not a result of the addition of more pipes to the previous optimum but depends on the overall pipe length allowed to be added to the existing WDS. The pipes chosen to be added can be completely different pipes, as can be noted for $L_{\max }$ = 200 m or $L_{\max }$ = 500 m in Figure 2a,b, respectively. In the case of less overall pipe length added, the closing of an outside loop seems to influence resilience more than adding an edge close to the first branching of the WDS. In turn, for a larger overall pipe length added, this outside loop is discarded and instead, branching close to the first branching is intensified, therefore leading to more low resistance alternative feeding paths. Similarly, this phenomenon can be observed for the different kinds of adapted WDS for $L_{\max }$ = 1000 m and $L_{\max }$ = 1750 m in Figure 2c,d, respectively. In general, pipes added close to the first branching of the existing WDS leading to a refinement and higher meshedness of this area are more likely to remain the same for the increasing instance of $L_{\max }$. However, pipes closing outer loops in the existing WDS switch depending on the overall pipe length allowed to be added and the already realized refinement of the area close to the first branching. Therefore, these pipes can be interpreted as a possibility for fine tuning of optimal resilience improvement, while pipes meshing the area of the first branching seem to have a significant influence on the overall resilience of the WDS. The latter observation is in accordance with the observations made by a brute force resilience adaptation study previously conducted by the authors [13].

The edges identified to add to the existing main line WDS by the optimization most likely already have pipes underlying the urban transportation structure, as most possible additional edges are along inhabited streets. Therefore, it can be useful to replace existing pipes with pipes of higher diameter. Thereby the resistance of the overall WDS is reduced and its feasibility is inherently given.

The conducted cost-benefit analysis for these instances, and the optimal relative resilience enhancement, show a clear convergence towards maximum possible resilience improvement, regardless of the allowed overall pipe length added to adapt the existing WDS. The results of pipe addition for different instances shows that, even for the highest overall allowed added length, not all possible pipes are added, but the stagnation of resilience improvement starts even for shorter overall allowed pipe lengths added. This indicates a limiting function of the urban structure for resilience enhancement. The discovered limit of resilience enhancement of WDS for the present urban structure leads to the further research question of how different urban structures compare considering their maximum possible resilience.

In conclusion, the developed optimization program and the conducted cost-benefit analysis of the addition of pipes to an existing main line WDS in order to improve its resilience has proven to be a potent tool. On the one hand, decision-makers such as water suppliers or governmental institutions are assisted when seeking to increase the resilience of the critical infrastructure for fresh water supply. The results show that the adaptation of the WDS towards optimal resilience enhancement for a given investment budget, proportional to the maximum overall added pipe length, is neither intuitive nor the result of a strict sequence of which pipes should be added first. The respective resilience enhancement aimed at by decision-makers can also be achieved by analyzing a coarse Pareto front and then deriving the respective maximum length added to the existing main line WDS. Furthermore, the observed limit of resilience improvement depending on the present urban structure shows that there exists a threshold of investment cost for which the resilience is further improved. On the other hand, the presented tool can additionally be applied for further urban structure analysis leading to maximum resilience of a WDS. This tool allows us to analyze as well as to assess the potential for resilience enhancement of a given WDS while focusing on physical feasibility within the present urban structure.