Integrating Structural Resilience in the Design of Urban Drainage Networks in Flat Areas Using a Simplified Multi-Objective Optimization Framework

Structural resilience describes urban drainage systems’ (UDSs) ability to minimize the frequency and magnitude of failure due to common structural issues such as pipe clogging and cracking or pump failure. Structural resilience is often neglected in the design of UDSs. The current literature supports structural decentralization as a way to introduce structural resilience into UDSs. Although there are promising methods in the literature for generating and optimizing decentralized separate stormwater collection systems, incorporating hydraulic simulations in unsteady flow, these approaches sometimes require high computational effort, especially for flat areas. This may hamper their integration into ordinary commercially designed UDS software due to their predominantly scientific purposes. As a response, this paper introduces simplified cost and structural resilience indices that can be used as heuristic parameters for optimizing the UDS layout. These indices only use graph connectivity information, which is computationally much less expensive than hydraulic simulation. The use of simplified objective functions significantly simplifies the feasible search space and reduces blind searches by optimization. To demonstrate the application and advantages of the proposed model, a real case study in the southwest city of Ahvaz, Iran was explored. The proposed framework was proven to be promising for reducing the computational effort and for delivering realistic cost-wise and resilient UDSs.


Introduction
Stormwater drainage is an essential service in safeguarding public health and safety [1][2][3]. Flooding not only causes an inconvenience to local inhabitants but also can cause significant economic damage and even mortality. Due to climate change, more extreme weather events are anticipated in the future. This, together with continuous urbanization, the increase in impervious surfaces, and population growth, can debilitate a city during a major storm event [4,5]. Therefore, urban drainage systems (UDSs) should be resilient to extreme loading conditions to lessen the span and extent of failure [4][5][6].
Two types of resilience pertain to UDSs: structural resilience and functional resilience. Structural resilience is the UDS's ability to minimize the consequences of structural failures, such as clogging and pipe cracking [7]. A UDS's functional resilience is its ability to maintain service or to return quickly to service during extreme loading conditions. In other words, the system can handle loads for which it was not designed [8].
Concurrently, the need to construct new UDSs in developing countries is increasing [9][10][11]. Nonetheless, planning an efficient UDS can be tricky as various cities do not have sufficient resources for stormwater management improvements and may not have the capital to address all regulatory requirements. As municipalities often work under budgetary limitations, they must allocate their resources wisely and make decisions to find the optimal tradeoff between the service objectives and the life cycle costs (LCC) [12,13].
Conventional approaches to design UDSs in "non-flat" areas include some fundamental stages. Firstly, an outlet is identified and its contributing area (catchment) is marked on a topographical map. Next, a preliminary horizontal alignment (layout) is produced. As pipe networks are supposed to collect wastewater and stormwater gravitationally, the designer can rely on the topography of the area and can follow natural ground slopes in the direction of the determined outlet. This approach might lead to a near-optimal layout depending on the designer's experience and the steepness of the area [14]. In flat areas, the problem is much more challenging to solve. There is no considerable variation in topography elevations; subsequently, the designer cannot see and trace prominent natural ground slopes to a distinguished outlet. In such areas, there are numerous possibilities for connectivity of the sewers and the outlet(s) location. In contrast with steep areas, engineering experience and judgment are not adequate to design an optimal sewer layout in flat areas. In such areas, the number of possible layouts exponentially increases with the number of sewers [15,16]. The lack of natural slopes in the network highly increases the sensitivity of design specifications and the construction and operational costs for layout configuration [14][15][16][17]. Therefore, utilizing optimization methods to generate optimal layouts is needed [3,16,18].
A limited number of sewer layout generation and optimization methods can be found in the literature, for example, in [15,17,[19][20][21][22][23]. Some of these works use graph theory to represent the layout of UDSs. Manholes, lift pump stations, storage tanks, and outlets of the system are depicted as nodes (vertices) in the graph, and sewers, under-pressure pumps, and weirs are represented as the edges.
Most of the methods presented so far for UDS layout generation can only generate and optimize centralized layouts with one main artery, which conveys the flow to a single outlet. However, centralized UDSs are costly to construct, maintain, and operate and have limited resilience and adaptability to future challenges [9,24,25]. Along with economic deterrents, there are concerns about the resilience of centralized UDSs or their ability to provide long-term sustainability [26][27][28][29][30]. Recently, the idea of centralized urban drainage networks has increasingly been questioned. The latest investigations suggest a transition from centralized to decentralized or hybrid schemes.
The current literature supports structural decentralization as a way to introduce structural resilience into an UDS [24]. Structural decentralization strategies reduce the degree of centralization (DC) by reducing the number of elements that are linked and interconnected by releasing stormwater by various outlets to different water bodies locations (see Section 2.2) [10]. Structural resilience is often neglected in the design of UDSs [31][32][33]. However, it can play a significant role in system performance during extreme events.
After layout design, the size of the pipes and their slopes as well as the size and location of pumping and storage facilities are calculated based on the estimated flows from the contributing area in a way that satisfies all technical and practical design constraints and criteria [15,34]. The so-called rational method is the most common approach to sizing the sewers of UDSs. The sewers are sized such that they convey the design discharge in the steady uniform flow [35,36]. Open-channel unsteady flow should be considered in the simulation of sized sewers to prevent possible shortcomings of some conventional approaches in finding optimal UDSs [35].
Sometimes, a trial-and-error design procedure is used to size the system's components. The size of components is suggested and verified for compliance with technical and hydraulic constraints. In the case of a constraint violation, a new design is proposed, and this procedure is iterated until a feasible design is obtained. This approach only provides a limited number of feasible designs that are not necessarily optimal, especially in the case of large systems in flat areas [16,18]. Besides not providing optimal solutions, trial-and-error approaches are time-consuming and require considerable human effort. Additionally, hydraulic reliability-based approaches, which assure adequate hydraulic capacity for the conveyance of design, do not equip the system for future challenges or guarantee ongoing reliability and resilience of services as the system ages and incurs wear and tear [8]. Therefore, recent academic literature has tried to utilize modern optimization methods (1) to deliver optimal solutions, (2) to make the design procedure systematic, and (3) to integrate different performance indicators (e.g., resilience and sustainability) into the design procedure [37,38].
These approaches mainly call for simulation-optimization procedures in which (1) a mathematical model is employed to mimic the real sewer network and (2) an optimization solver is linked to the mathematical model to aid decision-making by automatically generating and evaluating a large number of scenarios [36]. Although there is a massive number of academic works on the optimal design of UDSs over many years, the concepts and techniques developed have not yet found their way into routine practice [1]. One reason is that the aforementioned approaches are complicated and time-consuming, making them unfavorable for integration into conventional UDS design software like SWMM (Storm Water Management Model). Conjointly, most available optimization frameworks and tools for planning UDSs consider only the second sub-problem (hydraulic design) for a fixed (predesigned) centralized layout. However, as mentioned before, layout design and DC play significant roles in the ultimate design costs and structural resilience.
Recently, Bakhshipour et al., (2019) [24] introduced a UDS layout generator algorithm, namely the hanging gardens algorithm (HGA), to generate and optimize (de)centralized urban drainage systems for both flat and steep terrains. The proposed method was applied against a real case study: a section of the city of Ahvaz, Iran. Although their results showed that the proposed model can generate realistic layouts and can find near-optimal solutions, they only considered construction costs as the objective function. Moreover, their proposed optimization framework required about 100,000 hydraulic simulations to converge to the optimal solution for that test case.
In a follow-up study, Bakhshipour et al. (2021) [39] included other performance indicators such as resilience and sustainability in their design. To do that, they introduced a multicriteria decision-making model for sustainable planning of UDSs considering different centralized or decentralized strategies. Although their proposed framework is promising in handling many decisions, objective functions, and indicators for designing sustainable UDSs, they did not manage to solve the required high computation effort problem. For the same case study, their framework required about 800,000 hydraulic simulations to deliver optimal solutions. This study aims to improve the algorithms and frameworks introduced in the above two works to significantly lessen the burden of computations and to increase the efficiency of separate stormwater collection networks' simulation-optimization. To do that, this study introduces a simplified cost function and structural resilience indices that can be used as heuristic parameters for finding the optimal DC and layout configuration. These indices only use graph connectivity information, which is computationally less expensive than hydraulic simulation, to estimate the construction cost and structural resilience of each generated layout. The use of a simplified objective function allows for so-called prescreening of layout designs before the hydraulic optimization step. This leads to a significant reduction in the number of solutions to be computed.

Problem Formulation and the Proposed Method
Mathematically, the optimization of separate stormwater collection networks can be formulated as follows [39]: where d opt is the optimal choice for decision variable d that define the UDS. Decision variable d contains the elements of two sub-problems (Equation (2)): DC and layout parameters that determine the connectivity between pipes in the system in each part, and the distribution of the system as a whole when multiple outlet candidates are available.
where N SO is the number of selected outlets from a list of candidates and N PO is the number of possible candidate outlets.  [24,39] developed two simulation-optimization frameworks to solve this highly complex optimization problem. They tried to simultaneously find the optimal combination of layout, DC, and system hydraulic specifications among numerous candidates. In large networks, as the number of possible layouts exponentially increases with the number of sewers in the system [16], the problem becomes more intractable to solve and the simulation-optimization approach becomes computationally more and more expensive. For instance, the layout generator hanging gardens algorithm works on a random basis. As a consequence, it generates many irrational but feasible, in terms of graph theory, layouts. All the irrational layouts were evaluated using the hydraulic simulation software. This not only increased needless computational effort but also prevented the optimization solver from efficiently exploring the search space by forcing it to explore unfeasible parts of the search space.
By first designing the layout, the number of solutions taken forward for hydraulic optimization can be reduced. This is beneficial when model complexity and computational demand are high. A prescreening of the layout designs can be done regarding one or more chosen performance metrics such as structural resilience and/or construction costs. To do that, an appropriate and simple objective function is needed to evaluate and compare different generated layouts during optimization, based only on the sewer network characteristics-connectivity between edges (or sewers) and nodes (or manholes).
Hence, this study's first contribution is to develop, evaluate, and choose appropriate simplified objective functions. To do that, first, several simplified cost indices were introduced (see Section 2.2). Then, each index was used as the objective function to conduct a single objective optimization (SOO) to optimize the layout configuration. In other words, when designing the layout, the only objective function was the particular cost index. A second SOO was then performed to design the optimized layouts hydraulically, and each fully designed system was evaluated. After evaluation, one superior cost index was chosen.
The second contribution is the integration of structural resilience into the design of UDSs by proposing a fast multi-objective optimization (MOO) framework based on the simplified objective functions developed in the previous step (see Section 3.2). In this step, the MOO phase was conducted. The layout generation problem then had the following properties: (1) the chosen superior simple cost index and (2) the structural resilience index. As both objectives were based on the graph properties, this MOO was done very quickly in comparison to the works that had tried to solve the layout and hydraulic design sub-problems simultaneously. A postprocessing step was performed to select potentially interesting layouts from the resulting Pareto front of non-dominating solutions. Next, an SOO was done to optimize the selected layout(s) hydraulically.
The proposed framework employed several algorithms, including an algorithm to generate (de)centralized layouts (hanging gardens algorithm, Section 2.2), an algorithm to systematically size the sewers (adaptive algorithm, Section 2.3), and an SOO and a MOO engine (simulated annealing and Borg Multiobjective Evolutionary Algorithm (MOEA). Figure 1 demonstrates the relationship between these algorithms (proposed framework). In the next section, each algorithm is described in short. More detail about each method can be found in the given references. In this study, Simulated Annealing (SA) was employed for SOO (Section 3.2) [40], and the Borg Multi-Objective Evolutionary Algorithm [41] was employed for the MOO (Section 3.3) based on its successful application for water resource problems [42].

Layout Generator (Hanging Gardens Algorithm)
The layout of an UDS can mathematically be presented as a graph with specific properties. Graph theory is a branch of discrete mathematics that studies the principles of the mathematical expression of graphs. The hanging gardens algorithm, recently introduced by [24], is a sewer layout generator algorithm based on graph theory that can generate all possible sewer layouts and can explore different DCs. The hanging gardens algorithm starts with nominating several outlet candidates in the area under design and generating a centralized layout with an arbitrary outlet from the candidate list.
It uses loop-by-loop cutting to generate a feasible centralized layout. Then, other arbitrary outlets from the candidates are added to the generated layout considering the desired DC. The hanging gardens algorithm uses a graph-theory-based approach to assign parts of the layout to different outlets and to generate a decentralized layout. The algorithm needs 2 × (NL + N PO − 1) decision variables to generate one feasible layout. NL is the number of loops in the base graph, and N PO is the number of possible outlets. Figure 2 summarizes this procedure. For more in-depth information, please see [17] and [24].

Sizing the Network's Components (Adaptive Algorithm)
The next step for the sewer network with a given layout generated by the hanging gardens algorithm is to size the sewers and pumps. Sewer diameters and slopes as well as the number and location of pumping stations must be designed in a way that satisfies all the hydraulic and technical constraints. To satisfy technical constraints like the telescopic pattern, minimum cover depth, maximum excavation depth, and minimum and maximum slope, the adaptive approach introduced in [14,34] was employed. This algorithm can automatically handle all constraints in designing sewage collection networks. Designing a sewage collection network is based on steady-state flow, but as discussed in the introduction, designing stormwater collection networks must be done using full 1D Saint Venant equations for more accurate results. When unsteady flow equations are used to simulate flow in the pipes, flow-related constraints such as keeping flow velocity between the minimum and maximum bounds or maintaining the proportional water depth under the specified maximum value are not satisfied. In this study, the flow-related constraints were handled using penalty functions during optimization, and the following constraints were automatically satisfied using the adaptive algorithm described in the next paragraph: 1.
For each manhole, place upstream pipes at elevations higher than the downstream one; 2.
choose sewer diameters from the commercial list; 3.
maintain the minimum buried depth to prevent damages from the traffic loads and other surface activities; and 4.
for each manhole, assign the outlet pipe's diameter equal to or greater than the upstream inlet pipes' (telescopic pattern).
If NP is the number of pipes in the UDS, the adaptive algorithm needs 3NP decision variables to design the system hydraulically: one per pipe to assign the size of a pipe [D], one to assign the slope of it [S], and one that determines whether there is a lift pump station upstream of a pipe [P]. Herein, a vector having 3NP members in the range of (0, 1) named as the normal design vector is introduced to the model. This vector indirectly represents pipe diameters (the first NP members) and slopes (the second NP members) and directly gives the pump indicators (the third NP members). For more in-depth information, please see [34]. Here, it is explained how the adaptive algorithm determines the pipe sizes and adaptively satisfies the telescopic pattern constraint. Suppose that d is a vector generated during optimization by the optimization engine (e.g., SA or BORG). The adaptive algorithm uses this vector (d) to assign the diameter of pipes using the following equation [34]: where D max is the largest commercially available size and D min is determined concerning the telescopic pattern [20]. The size of each pipe must be equal to or larger than that of its upstream pipes, which mathematically means the following: where [DU] contains the pipe diameters connected to the upstream end of the pipe at hand. The hanging gardens algorithm calculates [DU] by simply using the layout connectivity matrix [24]. The layout generator algorithm and the adaptive algorithm can provide the essential materials for optimization. They need in total 2(NL + N PO − 1) + 3NP decision variables (generated by the optimization engine) to produce one structurally feasible (not necessarily hydraulically) UDS. The only remaining step is to define the objective function to evaluate different generated designs during the optimization. The next section introduces three simplified cost functions and one simple indicator to evaluate structural resilience.

Simplified Cost Functions (Proposed Indices) 2.4.1. Cost of UDSs
The cost of a UDS is dependent upon pipe diameters, slopes, and the specifications of pumping facilities. The pumps were disregarded in this study (for designing separate storm collection networks) as all layouts that required pumps were penalized. The pipe diameters and slopes are only known after the hydraulic design step but are dependent upon layout characteristics known during the layout design step. Therefore, assumptions about the cost of a UDS can be made based upon its layout characteristics. The proposed cost indices were developed based on hydraulic assumptions and sewer layout characteristics, such as pipe length and contributing area connected to each sewer. All assumptions were made to develop simplified indices further and to reduce the computational complexity of the UDS design optimization problem. Next, three simple cost indices are introduced: the length area index, the area diameter index, and the elevation rank index.

Elevation Rank Index (ERI)
Although it is difficult to determine the actual cost or resilience of a UDS layout without completing the hydraulic specifications, one method of determining relative cost and resilience is through heuristic techniques, which act to indirectly reach a desired goal (i.e., minimized construction cost and maximized resilience) through correlated relations of layout parameters with desired objectives. Villiers et al. (2017) [43] developed and proposed a list of 13 heuristic influence factors which can guide the optimization procedures towards UDS designs with these desired characteristics. Although decentralized layouts (layouts with multiple outlets) were not considered, it was concluded that three proposed factors, namely distance rank, elevation rank, and betweenness centrality, have a promising correlation with construction cost. They act by shaping different layout parameters to promote a reduction of the average diameter size of the system. Here, the elevation rank index (ERI) was adopted as the first simplified cost metric.
To calculate the elevation rank, all nodes or manholes were reordered from the lowest elevation to the highest elevation and then ranked accordingly, with the node with the lowest elevation ranked as 1 and so on. Elevation Rank Difference (ERD) is the difference between the ranks of connecting manholes, i.e., each pair of the upstream and downstream manhole, and is formulated in Equation (6). This will result in an NP of values. Then, the standard deviation of the ERDs is computed. The standard deviation of the ERDs is one component of the elevation rank index.
where UER i is the elevation rank of the upstream manhole for pipe i and DER i is the elevation rank of the downstream manhole for pipe i. The standard deviation of ERDs has been proven to have a positive correlation with system construction cost. In other words, when the standard deviation of ERDs decreases, so does the construction cost [40].
Initially, elevation rank was not developed considering decentralized sewer system layouts. To adapt the index to decentralized UDS layouts, two additional values, Std[Nn] and r, (Equation (7)) were incorporated into the ERI. Std[Nn] is the standard deviation of the number of nodes (manholes) in each tree. This will inherently favor layouts with smaller DC, as layouts with fewer trees will have higher variances. r is the ratio of the number of nodes in the largest tree to the total number of nodes. This ratio promotes more evenly distributed trees with respect to the number of nodes per tree.
The fully proposed ER index is as follows: ERI can be used as a simple metric that would minimize the elevation rank of a sewer layout design and inherently the construction costs.

Length Area Index
The length area index (LAI) relates the relative system cost to the length and contributing area (the cumulative summation of the impervious area upstream, which drains to a pipe) of each pipe. This index was inspired by the index that Walters and Smith (1995) [44] proposed. Their simple objective function reflected the layout cost in terms of the length and concave function of the flow rate of each link. Here, as calculating the flow rate requires hydraulic simulation, it was proposed to use the cumulative impervious area connected to each pipe as a surrogate for the real flow rate. The length of pipes (L i ) can be retrieved from the graph connectivity matrix (see [24]). The accumulated connected impervious area (A cum,i ) can be implicitly obtained using the layout configuration. The final index value is the sum of the length of each pipe by the square root of the contributing area (Equation (9)).

Area Diameter Index (ADI)
The area diameter index (ADI) was formed from hydraulic assumptions based on the total impervious area connected to a single pipe. As mentioned previously, the cost of a pipe depends predominantly on its diameter, and diameter is dependent on runoff volume draining into the sewer. Consequently, larger contributing areas lead to larger pipe design diameters. Therefore, an estimation of the cost of each pipe can be obtained by relating the diameter of a single pipe directly to its contributing area. The contributing area can be easily calculated using the layout configuration and is usually already known when designing a sewer network. To overcome the inability to calculate the exact diameter during the layout design stage, a relation between the contributing area and design diameter was assumed.
To do that, for each available commercial pipe size, the rational method and a predefined design storm were used to calculate the maximum impervious area that can be connected to a pipe with a certain diameter and fixed slope. In such areas, satisfying the maximum excavation depth constraint and avoiding too many lift pump stations are very problematic tasks. Therefore, the minimum allowable slope is usually considered the optimal slope in flat areas [14,24]. By fixing the slope and using the rational method, a relation between the contributing area and design diameter can be obtained. Alternatively, the ADI can be calculated using a simple trial-and-error approach in SWMM. To do that, it is enough to construct a simple model with a sub-catchment and a pipe connected to it. Then, for each available commercial pipe with a fixed slope, the maximum area that can be drained by that can be determined. In this approach, it is possible to consider different types of design storms (e.g., Block Rain and Euler-type II) for a more robust design. The maximum area connected to each pipe size is the minimum calculated impervious area for all different design storms. It is also possible to perform a long-term simulation using historical precipitation data. After calculating the ADI using different types of design storms, each diameter and its corresponding impervious area as calculated are tested using the long-term simulation. If the capacity of the calculated diameter(s) is not enough, the maximum impervious area must be updated. The single events from historical data that cause flood problems in this step are extracted. For hydraulic optimization, any desirable type of design storm and several extreme events can be used. As our proposed approach (see Section 3.2) can significantly reduce the number of hydraulic simulations, this is not computationally problematic. However, doing that can increase the robustness of the final design meaningfully.
As explained before, it is easy to retrieve the length of pipes (L i ) and the accumulated connected impervious area (A cum,i ) from the graph connectivity matrix. Upon having the length, diameter, and slope of each pipe, it is possible to estimate the construction cost of a generated layout using a cost function. An example of this method is given in Section 3.2.

Structural Resilience Index (SRI)
The maximum structural resilience is obtained when clogging in a sewer has the least effect on its upstream lines [45]. The structural resilience can be considered during the layout design problem, as will be explained in the following paragraphs. On the other hand, evaluation of the UDS's functional resilience can only be done using hydraulic simulations (after designing the layout and sizing the sewers).
To quantify the clogging consequences in sewer networks, Haghighi and E. Bakhshipour (2015) [45] proposed a simple criterion that uses layout design properties. This criterion can indirectly evaluate how well a sewer system can react to sewage flow blockage. The principle adopted is that "Averagely, less population affected by a pipe clogging the sewer network would be more reliable". If it is assumed that the probability of simultaneous sewer clogging in the network is extremely low, one-at-a-time clogging is adapted here to evaluate the resilience of the system's performance. The resilience of each individual sewer (Res str,i ) is mathematically measured by the following index.
where A cum,i is the area connected to pipe i and A total is the total area. For the entire network, the averaged resilience index (Res structrual ) is measured as follows: Res structrual can be very insensitive to large networks as there are a large number of upstream pipes with low discharges and high resilience. This restricts Equation (11) from comparing alternative layouts. To reduce this effect, for large sewer networks, using only sewers with a resilience index less than the threshold 90% was suggested [45]. Therefore, Res str,i became Res str,i<90% . This removed the upstream branches from consideration. As a result, this caused the layout optimization to be more sensitive to designing a layout with one main collector, which conveys high sewage discharges. To account for the effect of DC on the resilience of the layout, Bakhshipour et. al., (2020) [39] proposed Equation (12)  (%) 100% i f NP Res str,i<90% = 0 (12) in which NP Res str,i>90% is the number of pipes with structural resiliency more than 90% and NP Res str,i<90% is the number of pipes with structural resiliency less than 90%. SRI is equal to zero when all sewers are connected to more than 10% percent of the total area (NP Res str,i>90% = 0) and 100% if each outfall is connected to up to 10% percent of the total area (NP Res str,i<90% = 0).

Case Study
To show how the proposed framework works, a real case study introduced in [46] was considered. The case study is a section of the city of Ahvaz, which is in the southwest of Iran. With a semi-dessert climate, entirely flat topography, and high groundwater level, Ahvaz is highly urbanized and has no stormwater management system (it has only an aged sewage collection system). Due to this, Ahvaz encounters urban flooding, which causes public inconvenience, and economic and environmental destruction. Due to the aforementioned issues and the unaffordable high initial investment costs, it is not feasible for the city to construct a centralized stormwater management system. Consequently, a decentralized separate stormwater collection system is the favored approach for this area.
The section of Ahvaz considered here is 500 ha, and the total area of the city is 18,500 ha. The base graph (Figure 3) has 181 sub-catchments (loops in the base graph), 530 pipes (total length of 75 km), and 10 candidate outlets [46]. The design constraints, including minimum slopes for different diameters, are presented in Table 1. As the minimum required slope for smaller diameters is larger than that for larger pipes, the optimal solution is implicitly related to both parameters. That means that the model sometimes prefers using larger diameters in some branches to avoid deep excavation and the lowering of all other parts of the system.   [24]. For the hydraulic simulations to design the sewers (second sub-problem), the dynamic wave approach in SWMM was used with a 5-year design storm (6-hour duration with a total depth of 30.2 mm) for a corresponding hydraulic simulation. Surcharging was allowed to use the full capacity of the pipe networks, but flooding was not acceptable. The maximum allowable velocity in the pipes was 4 m/s, and the maximum allowable excavation depth was fixed at 5 m. Layouts that would need lift stations or designs that cannot satisfy the abovementioned hydraulic constraints were automatically omitted using a penalty function.

Maximum
For evaluating the functional resilience of the selected solutions in the last part of this study, design storms with return periods of 10, 20, 25, and 50 years were used (with respective total depths of 38.3, 46.7, 49.5, and 57.5, as shown in Figure 4). Besides the storm designs, the selected designs were evaluated using six months of recorded rainfall data from October 2018 to March 2019. This period was selected because of several extreme events that caused a lot of trouble in the area. Total precipitation during this period was 268.4 mm, which was 45% more than the 30-year average precipitation. The maximum precipitation in 24 hours during this period was 45.1 mm, which was close to the 10 years design storm but happened in 24 hours instead of 6 hours ( Figure 5).

Analyzing the Cost Indices
To evaluate each cost index, a single objective optimization problem coupled with the simulated annealing algorithm was formulated. The problem was structured as two separate optimization problems: the first for layout design and the second for hydraulic design. The layout design was optimized considering a simple cost index, and then, the resulting optimized layout configuration was progressed to hydraulic optimization. The life cycle cost (LCC) was the objective function of the hydraulic optimization, which included the capital and operation and maintenance costs for a 30-year service life. The cost function used in this study can be found in [24]. The HGA was applied to optimize the layout, and then, the resulting layout was hydraulically designed using the adaptive algorithm and SWMM hydraulic simulation. This was done for the LAI, ADI, and ERI. After the hydraulic optimization was completed, each fully designed UDS was evaluated for average and maximum diameter and buried depth. Table 2 presents a summary of these parameters for each optimized design. After evaluating the LCC, it can be said that the ADI was the superior simple cost index for this special case study. For further analysis, the LAI and ERI were skipped and only ADI and SRI were used. Table 3 presents the calculated AD parameters for the test case. To calculate the ADI for the test case, the 5-year design storm, the recorded precipitation (Figures 4 and 5), and the approach described in Section 2.4.4 were used. For all commercial pipe sizes, the design storm was decisive in our case. Therefore, no single extreme event from the recorded time-series was taken here for hydraulic optimization.

Introducing the Fast MOO Framework
First, a MOO for the layout design was done with the objective functions of the ADI and SRI. As both ADI and SRI work only based on the graph connectivity matrix (inexpensive calculations), this step can be done quickly. Depending on how large the area under design is and how fast the computer, is it can take from minutes to hours.
In the next step, selected layouts from the resulting Pareto front were hydraulically optimized using the adaptive algorithm. The layout design optimization step was coupled with BORG MOEA, and the hydraulic optimization step was coupled with the simulated annealing algorithm (SA).
For hydraulic optimization, the results of the AD index (from Table 4) were used as the initial solution for the SA. This means that the maximum diameter (D max ) was assigned for each pipe as the diameter suggested by the ADI, as presented in Equation (13). For example, if the total impervious area connected to a particular pipe in the network is 60 ha, the D max from Table 3 where D ADI is the size of each pipe in the network based on ADI and D min is determined concerning the telescopic pattern and smallest commercially available size, as explained in Equations (4) and (5). D commercial is the commercially available sizes for this case study. Doing this significantly reduces the search space and accelerates the optimization procedure. As the ADI assigns each pipe's diameter based on the rational method and the minimum allowable slopes, it usually tends to overestimate the size of pipes. Therefore, by starting the optimization from the design that the ADI gives (set d for all pipes in Equation (13)), the optimization is simplified from assigning the size of all pipes from a completely random space (that has many infeasible solutions) to reducing the size of pipes from a feasible near-optimal solution. In our test case, D max and D min for many upstream pipes were equal or there were only one or two alternatives for each of them. For example, considering the available diameters in Table 1, if D max that the ADI dictates for a particular pipe is equal to 0.4, the only available sizes for this pipe during optimization are 0.25, 0.35, and 0.4. Without including this, the SA must choose randomly from the whole list of available sizes (10 alternatives in Equation (13)). Considering the telescopic pattern, if the optimization engine assigns a larger diameter (the probability of doing that in this example is 70%), all pipes downstream will be overdesigned. This approach, not reducing the search space, significantly increases the number of iterations that the optimization engine needs for convergence.

MOO Results Analysis
The result of the MOO for the case study was a Pareto front containing 33 nondominated layout alternatives, as can be seen in Figure 6.
Three promising layouts from the Pareto front were chosen to be hydraulically designed. The following three layout alternatives were chosen based upon their SRI and ADI value (LCC): (1) the least costly layout (design 1), (2) the most resilient layout (design 3), and (3) the median-valued alternative (design 2), almost directly in between the first two selected layouts (design 1 and design 3). The purple squares present the final cost of the selected layouts after hydraulic optimization with SA. The hydraulically optimized designs for designs 1-3 can be seen in Figures 7-9. Figure 10 shows the optimal design suggested in   [24], where only the construction cost is considered as the objective function (design 4). In the following paragraphs, these designs are analyzed and compared in depth.

Analyzing the Computational Efficiency
The LCC of the first design is identical to the LCC of design 4. Recall that design 4 was obtained from Bakhshipour et. al., (2019) [24], where they formulated a single objective optimization (only for cost) using a simultaneously layout-hydraulic optimization method. In theory, the simultaneous single optimization might result in more promising solutions, but it did not here. The reason is that the predefined maximum number of iterations (100,000 in that study) might not be enough to reach the global optimization.
Optimizing the first scenario hydraulically using the proposed algorithm in the current framework only required approximately 1500 hydraulic simulations. Additionally, this design guaranteed a much higher structural resilience (85.3% instead of 79.9%). This proved that the proposed framework can integrate structural resilience into the UDS design procedure and can significantly reduce the computational burden. The number of required hydraulic simulations for the design 2 and 3 hydraulic optimizations were about 2000 and 2500, respectively. It can be said that, on average, the proposed framework reduced the computational effort by 98% for the hydraulic simulations in comparison with the framework proposed in Bakhshipour et. al (2019) [24] for this specific test case. The proposed framework in Bakhshipour et. al. (2020) [39] can also consider structural resilience, but for the same case study, it needed 800,000 hydraulic simulations until convergence. However, as they integrated other objectives like reliability, functional resilience, and sustainability and green-gray-blue infrastructures as decisions instead of pipe-only in this study, the results cannot be compared directly.
By looking carefully at the Pareto front, it can be interpreted that, although the hydraulic optimization step can reduce LCC, this reduction is not substantial. For the three selected layouts in this study, the average reduction was less than 5.0%, indicating that the ADI can generate near-optimal pipe diameters for designing stormwater collection systems in flat areas. Therefore, in such areas, the DC and layout configurations mostly determine the cost and structural resilience of the system. Notwithstanding, currently in the literature, most efforts in UDS optimization are focused on sizing the sewers. It can be concluded that it is possible to skip the second optimization (hydraulic design) for the preliminary design aims.

Analyzing the Structural Resilience
The Pareto front in Figure 6 shows the relationship between LCC and structural resilience that can be described using a polynomial equation. This trend remained identical for the selected layouts after hydraulic optimization (purple squares). To gain 1.2% structural resilience, from design 1 to design 2, 3.5% more capital investment is needed, and to gain only 0.5% more structural resilience, from design 2 to design 3, 4.5% more capital investment (8% in total) is required.
In Figures 7-10, the structural resilience of critical pipes and outfalls is demonstrated. These values implicitly indicate the fraction of the total area not affected by any structural failure. Design 1 has a total of 12 pipes, with structural resilience of less than 90%. However, there are 21 critical elements in design 4 with very similar LCC. Comparing the structural resilience of designs 1 and 4 (85.3 and 79.9 in order) proves that the proposed index is a reliable indicator for assessing and optimizing the structural resilience of UDSs. The critical elements in designs 2 and 3 are very identical. Interestingly, in design 3, there is no critical non-outfall element. This shows that the total area is well distributed between different outlets and in each part of the system. The minimum structural resilience in both designs 2 and 3 is 84%, which is 10% more than that for design 4. Design 2 has only 1 non-outfall critical element very close to one outfall. Therefore, spending 4.5% more capital to go from design 2 to design 3 might not be reasonable.
Recognizing the critical points in the system can be crucial. By strengthening the system in these points, one can make sure no single structural failure can affect a predefined fraction of the serving area (10% percent in this study). Strengthening the weak points can be achieved by introducing extra elements to the system. For example, by connecting a manhole near a critical point using an extra sewer to a manhole in an adjacent part (tree).
Not considering the structural resilience in the design of the existing sewage collection network caused trouble for citizens during the fluvial flood that happened in winter 2020. The fluvial flood dramatically increased the groundwater level in the riversides. Although the existing levees successfully protected the city from a destructive flood, the infiltration of river water and groundwater into the pipes resulted in choking and blockage of several main collectors of the existing centralized sewage collection network (structural failure) in parts near the riversides. Consequently, several parts of the network, even far away from the riversides, were out of service. This led to spilling from the manholes and overflowing wastewater into the streets, resulting in serious disturbances and health problems for the citizens [29]. Therefore, by integrating the structural resilience in the design of new UDSs using the proposed framework, the adverse consequences of such extreme events can be lessened.

Analyzing the Functional Resilience
In the end, the relation between structural resilience and LCC with the functional resilience of UDSs was investigated. The functional resilience of each scenario was assessed with the hydraulic performance indicator (HPI) using equation 14 adopted from [47]. For this purpose, the design storms with the return periods of 10, 20, 25, and 50 years were used. The HPI is defined by the following [47]: (14) where V f looding is the total water that overflows the manholes and V runo f f is the total runoff volume. Table 5 summarizes this analysis. After an initial inspection of the functional resilience for all designs during each design storm, designs 1 to 3 seemed to exhibit similar values. Although when comparing both designs with similar DC (designs 2 and 3) it was noticed that design 2 has slightly higher functional resilience for each design storm despite its lower structural resilience, in general, no distinct relationship can be identified between structural resilience and functional resilience. However, it can be seen that design 1 (DC = 22%) and design 4 (DC = 33%) have lower functional resilience than designs 2 and 3 (DC = 11%). It can be concluded that functional resilience is more a function of DC and storage volume in the system than a function of structural resilience. Also, the performance of each design was evaluated using the recorded time-series represented in Figure 5. The results showed that all design scenarios can handle this challenging time period without any flooding.

Summary and Conclusions
This study introduced a simplified multi-objective optimization framework for integrating structural resilience in the layout design of urban drainage systems (UDSs). First, three simplified cost indices that only use simple layout graph connectivity information were introduced. Next, each index was analyzed and evaluated using a single optimization method. Then, one superior simplified cost index (Area Diameter Index (ADI) in this study) was chosen. To integrate structural resilience in the design of UDSs, a multi-objective optimization of the layout design considering structural resilience and the ADI index was conducted. As both ADI and the structural resilience index work only based on the graph connectivity matrix (inexpensive calculations), this step was not computationally expensive. Finally, some selected layouts from the resulting Pareto front were hydraulically optimized using a novel computationally efficient approach based on the ADI. The key conclusions resulting from this research are as follows:

•
The proposed framework can significantly reduce the computational effort needed for optimizing UDSs in flat areas without a noticeable sacrifice in the quality of solutions. Doing that will increase the potential of the proposed frameworks and algorithms to be incorporated into commercial UDS design software to deliver more sustainable and less expensive designs.

•
As the number of required hydraulic simulations for optimizing sewers' sizes was significantly reduced (98% in our test case), it is possible to consider different types of design storms and historical precipitation data within the proposed framework. That leads, apparently, to more robust designs.

•
The proposed indicator of structural resilience can reliably evaluate the structural resilience of different UDSs. Furthermore, the proposed framework can integrate structural resilience into the UDS design procedure.

•
In flat areas, the layout configuration and the degree of centralization are the most challenging and decisive problems for optimizing the UDSs, and sizing the sewers can be handled with simple optimization methods, as proposed in this study.

•
There is no apparent relation between functional and structural resilience in UDSs. Therefore, to build these different types of resilience into the system, completely different strategies must be taken.
Finally, based on this study's results and a review of the existing literature in the field, we suggest the following Table 6 as a guideline to choose an appropriate approach for different UDS design optimization problems. Although the focus of the present study was on the design of separate stormwater collection networks, layout selection is very relevant to the design of combined and separate sewage collection systems. Other crucial criteria must be considered for sewage collection layout design. These include sulfide control, air and water pollution, and requirements for wastewater treatment, with all of its high costs of construction and operation associated. Sulfide control may be particularly important in large systems and/or hot climates. For more information, we highly encourage the readers to refer to [48,49]. It worth mentioning that the proposed model for sewage collection systems design in [15] included suitable pumping systems or lift stations and that the costs associated with WWTPs is deterministic, efficient, and very fast, at least when the hydraulic simulations in unsteady flow are not considered or coupled with the design hydraulic model in steady uniform flow.