Toward Sustainable Urban Drainage Infrastructure Planning: A Combined Multiobjective Optimization and Multicriteria Decision-Making Platform

First, stakeholders or decision makers introduce all types of urban drainage systems and technologies that they wish to be considered in the design process in addition to all technical and practical considerations. Any screening tool for GBI type or sanitation technology selection, location identification, and narrowing of dimensioning variables without detailed optimization can be employed in this step. Primary screening restricts the search space and enhances optimization efficiency (Zhang and Chui 2018). Examples of screening tools can be found in Eaton (2018), Inamdar et al. (2013), Johnson and Sample (2017), Kuller et al. (2019), Martin-Mikle et al. (2015), Mitchell (2005), and Spuhler et al. (2018).

Next, all combinations of the proposed systems and technologies are generated systematically to perform the mathematical optimization. To do that, Step 1a outlines a base graph for the pipe network which includes all drainage feasibilities concerning street alignments, topology, barriers, watercourses, and existing sewers in the area under design. In Steps 1b and 1c, candidate locations of outlets; potential GBI type, size, and location; and, if necessary, treatment solutions such as on-site facilities, constructed wetlands, or traditional wastewater treatment plants are determined. These steps provide the potential decision variables $\mathbf{d}$ of the optimization problem [Eqs. (1) and (2)]. Step 1d defines all relevant performance indicators, such as reliability, resilience and sustainability, for MOOs. To prepare for the simulation-optimization step (Step 2), it is crucial to define simple performance indicators that are inexpensive to evaluate. Otherwise, the optimization in Step 2 may take a long time to converge and come close to acceptable solutions. For example, when evaluating hydraulic reliability or functional resilience, one or more design storms (which can also be considered for a more robust design) or multiple extreme events derived from historical rain data should be used instead of recorded time series. The section “Results and Discussion” and the studies mentioned in the previous section can aid selection of appropriate indicators for a specific problem.

In Step 1e, a numerical UDS model (EPA SWMM in this study) is constructed to evaluate the predefined performance indicators. Here the model parameters (e.g., impervious area, manning roughness, and soil characteristics) are calibrated or estimated, and design conditions and physical and technical constraints (e.g., design loadings, maximum velocity, minimum slope, maximum buried depth, and water quality standards) are defined. Based on the resulting data, a base model is constructed from which different design schemes are generated during the simulation-optimization step (Step 2).

Mathematically, the multiobjective optimization of HGBGIs can be formulated as follows:where ${\mathbf{d}}_{\mathrm{opt}}$ = optimal choice for decision variables that define the UDS; and $\mathbf{d}$ contains elements of at least three subproblems [Eq. (2)]:

To perform a systematic optimization, it is vital to systematically generate various HGBGI schemes that satisfy all technical and physical constraints. Three different modules are employed within the proposed framework to do this: (1) the hanging gardens algorithm to generate feasible layouts with an arbitrary degree of (de)centralization; (2) an adaptive algorithm to hydraulically design the generated layouts; and (3) an algorithm to define GBI type, size, and location. Fig. 2 schematically shows the simulation-optimization procedure proposed in this study and demonstrates the connection between different algorithms.

The hanging gardens algorithm, recently introduced by Bakhshipour et al. (2019a), generates a sewer layout based on graph theory, generating all possible sewer layouts and exploring different DCs. The 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. Then it adds other arbitrary outlets from the candidates to the generated layout considering the desired DC. Using graph theory, it assigns parts of the layout to different outlets and generates a decentralized layout. The algorithm needs $2\times (NL+N_{PO})$ 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. Full details on the hanging gardens algorithm are available in Bakhshipour et al. (2019a).

For each generated layout, CGI specifications such as pipe diameters, pump stations, and invert elevation are designed in a way that satisfies all hydraulic and technical constraints. Technical constraints such as telescopic pattern, minimum cover depth, maximum excavation depth, and minimum and maximum slope, are satisfied via the adaptive approach introduced in Haghighi and Bakhshipour (2015a, 2012). The hydraulic constraints of maximum velocity and no flooding are handled by a penalty function during optimization. If $NP$ is the number of pipes in the UDS, this algorithm needs $3NP$ decision variables to design the system hydraulically: one variable per pipe to assign the size of a pipe, one to assign the slope, and one to determine whether there is a pump station upstream of a pipe. Finally, distributed measures, here GBIs, are added to the designed UDS using an adaptive algorithm to construct a hybrid green-blue-gray UDS, as explained in the following paragraph.

Generating a hybrid alternative requires three extra decision variables for each subcatchment ($3NS$ decisions, where $NS$ is the number of subcatchments). For each subcatchment, a binary variable decides whether it is equipped with any GBI. The second variable determines the type of GBI(s) from a list of feasible candidates for each subcatchment using Eq. (4). The third variable defines the size of the GBI(s) considering the minimum and maximum feasible size of each measure using Eq. (5)where $gbi\text{\_}{type}_i$ = uniform random variable generated by the optimization engine; ${NGBI}_i$ = number of GBI candidates in subcatchment $i$; and $GBI\text{\_}{Type}_i$ = decoded GBI type in that subcatchment. As an example, suppose the GBI candidate list in one subcatchment is as follows, and ${gbi}_i$ randomly generated by the optimization engine is 0.67:

Using Eq. (15), ${GBI}_i=round(1+(7-1)\times 0.67)=5$, which means that $Rainbarrel+Greenroof$ is selected for that catchment.

$\mathbf{GBI}\mathbf{\text{\_}}{\mathbf{Size}}_{\mathbf{i}}$ is handled as follows:where $gbi\text{\_}{size}_i$ = uniform random number generated by the optimization engine; $GBI\text{\_}{Size}_{\min ,i}$ and $GBI\text{\_}{Size}_{\max ,i}$ = minimum and maximum permissible sizes of $GBI\text{\_}{Type}_i$, respectively; and $GBI\text{\_}{Size}_i$ = decoded size of $GBI\text{\_}{Type}_i$ in that subcatchment.

Using the algorithms, any arbitrarily generated set of decision variables $\mathbf{d}$ constructs a feasible HGBGI scheme. These algorithms, together with performance indicators and a numerical model for evaluating system performance (calculating values of indicators), provide essential tools and materials to perform MOOs. The MOO engine generates a set of Pareto-optimal solutions using different sets of decision variables that satisfy the defined constraints. There are several genetic algorithms (GAs) based MOO engines (e.g., NSGA-II, BORG, GALAXY) in the literature that can be employed for the current optimization (BORG in our study).

The MOO in Step 2 identifies the solutions that dominate other solutions and constructs the Pareto front. As there is no clear preference between the solutions on the Pareto front according to the optimization objectives from Step 2, a final decision step is required to help decision makers select an appropriate option. This step can account for aspects not covered by the MOO in Step 2.

The MOO compares different solutions by evaluating only a limited number of objective functions. However, some criteria cannot be considered for practical reasons such as computational power limitations. For example, applying the global resilience analysis introduced by Mugume et al. (2015b) to evaluate the structural resilience of UDSs requires generating and evaluating many pipe failure scenarios for each solution. This increases the computational burden exponentially. Second (and again due to the computational burden), modeling flooding of specific locations in the urban area has to be done separately using a 2D model based on an accurate detailed digital terrain model (DTM). This is feasible only for a few selected solutions. Third, there are other aspects that cannot be represented mathematically and criteria that cannot be quantified adequately (e.g., risk of human casualties or public nonacceptance).

Multicriteria decision analysis (MCDA) techniques [e.g., analytic hierarchy process (AHP) and analytic network process (ANP), technique for order of preference by similarity to ideal solution (TOPSIS)] provide the opportunity to include numerous ranges of indicators and thoroughly analyze the limited number of solutions that exceed this step (Malczewski and Rinner 2015; Wang et al. 2017). Hence, in this step of the proposed framework, a full range of desired indicators is first defined. Then, the performance of all selected solutions is evaluated. Finally, the solutions are ranked employing an MCDA. No specific MCDA technique is prescribed here for this purpose, as the choice must be suitable for the problem. For our application, we use TOPSIS in the next step.

We considered two simulation-optimization scenarios for the case study. In the first scenario, the proposed framework is applied to the four fixed layouts with different DCs recommended by Bakhshipour et al. (2019b). Therefore, the problem is simplified to the optimal sizing of the sewers and the location of the GBIs. In the second scenario, the layout and DC variables are also considered. This scenario provides the material to explore the effect of layout configuration on the different performance criteria.

In both scenarios, minimizing LCC is used as an optimization objective because LCC is the most determinative parameter for the stakeholders in the area and is different from other performance criteria, thus presenting an independent problem dimension. The second optimization objective is maximizing total sustainability ($Su{s}_{\mathrm{Total}}$), here defined as the geometric mean of reliability, resilience, and environmental sustainability. The reasons are as follows: (1) considering each performance indicator as a separate objective function increases the computational effort exponentially due to the large scale of the test case (more than 1,000 decision variables); and (2) as discussed earlier, these performance indicators have a pyramidal structure with an unknown relationship that can be obtained through the proposed formulation, as explained in the following paragraphs

As mentioned before, sustainability has three aspects: public, economic, and environmental. The economic aspect of sustainability is regarded here as a separate objective function (LCC). In addition, increasing reliability and resilience in Eq. (13) results in decreasing urban flood probability, which automatically enhances public sustainability. Therefore, we can guarantee that all aspects of sustainability are acknowledged in our proposed MOO formulation. In this study, for convenience of interpretation all introduced indicators are a real number between zero and one; zero indicates the lowest performance of each indicator; one, the highest performance. Generally, the geometric mean is used when considering several criteria that cannot compensate for each other; that is, one aspect rated zero sets the overall performance to zero.

By using the geometric mean, total sustainability is zero when one or more of the indicators are equal to zero, and one if and only if all indicators are equal to one. Therefore, the suggested pyramidal structure of the performance indicators can be obtained. The remaining technical and social aspects and indices are treated in the decision-making step. It is worth mentioning that the choice of performance indicators and MOO problem formulation depend on the problem at hand. We cannot prescribe any general formulation here. However, any MOO formulation can be handled using our proposed framework.

For the second scenario, maximizing DC is also used as an extra objective function. This means that, for a certain amount of LCC and total sustainability, solutions with greater DC (lower number of parts or outlets) are selected. This allows the optimization engine to fully explore the feasible space and determine the influence of DC on the performance indices. Therefore, the general optimization problem for the test case is reformulated and simplified as follows:

Table 1 summarizes the number, type, and function of decision variable vector $\mathbf{d}$ for the test case.

In this study, the Borg multiobjective evolutionary algorithm (Hadka and Reed 2013) is employed for the MOO based on its successful application to water resource problems (Eckart et al. 2018). More detail about this algorithm is given in Appendix.

As described earlier, some solutions from the Pareto front obtained in the previous step are selected in this step for more comprehensive assessments and then for the final decision. The SWMM simulations for the selected solutions calculate all additional indicators and rank the solutions. In addition to the simple indicators used in the simulation-optimization step, 15 additional indicators (Table 2) are used here for final decision making. These indicators are evaluated under continuous simulation, which uses six months of rainfall data from October 2018 to March 2019, a period recognized for its extreme events in recent years that caused many problems in the area, such as flooded streets, infectious diseases, traffic jams, and public protests. The total rainfall depth during this period was 268.4 mm, which was 45% more than the 30-year average. The maximum precipitation in 24 h during this period was 45.1 mm.

The proposed indicators are divided into four categories:

Here TOPSIS is used to aid in prioritizing solutions, and the entropy method is adopted to assign indicator weights (Li et al. 2014). TOPSIS, first proposed by Hwang and Yoon (1981) and later developed by Yoon (1987) and Hwang et al. (1993), ranks alternatives based on relative similarity to ideal solutions (Roszkowska 2011). It is a practical technique with an intuitive and clear logic that represents the rationale of human choice. It can be straightforwardly applied and enjoys high computational efficiency (Hung and Chen 2009).

As the determination of a specific weight for each index in TOPSIS is usually subjective, the entropy method was utilized to calculate the weights and reduce the subjectivity (Li et al. 2014). Entropy is a term in information theory introduced by Shannon (1948). The calculation steps of the so-called entropy method for weight determination are as follows, supposing a decision matrix $x_{ij}$ with $m$ alternatives (rows) and $n$ indicators (columns) (Wang et al. 2017):

As a first analysis, we look at the influence of DC on the Pareto front. For the fixed layout scenarios, DC has two main effects. First, all solutions with a certain DC are dominated by solutions with a lower DC, so the lower DC values are beneficial. This is more obvious for the most centralized alternative (e.g., compare $\mathrm{DC}=100\%$ and $\mathrm{DC}=55\%$). However, the more decentralized scenarios ($\mathrm{DC}=33\%$ and $\mathrm{DC}=0\%$) are close to each other, indicating that more structural decentralization does not cause a recognizable additional benefit. The same trend is observed for the second scenario, where the layouts are not fixed. The distance between Pareto fronts is high for more centralized solutions and decreases as DC decreases. Second, the range of solutions discovered by the MOO engine is highly dependent on the DC. As can be observed in Fig. 4, the broadest range of solutions belongs to the layout with $\mathrm{DC}=66\%$ (total sustainability ranging from 40 to 95), followed by $\mathrm{DC}=33\%$, $\mathrm{DC}=0\%$, and $\mathrm{DC}=100\%$. As a result, a very centralized or very decentralized layout might restrict the feasible search space, so more investigation is needed to determine whether there is a meaningful and insightful relation between the feasible search space or number of solutions on the Pareto fronts and system flexibility.

Our next analysis examines the effect of layout configuration flexibility on the quality of solutions for the optimum HGBGI. As seen in Fig. 4, for a fixed DC all optimal solutions of the second scenario (layout optimized) dominate solutions of the first scenario (fixed layout). The solutions of the second scenario with $\mathrm{DC}=88\%$ (two outlets) practically overlap with the solutions of the first scenarios with $\mathrm{DC}=66\%$ (four outlets). Particularly for $\mathrm{DC}=100\%$, solutions from the second scenario are observed to explore a broader range of objectives than the solutions from the first scenario. Only for $\mathrm{DC}=33\%$ are the results of the two scenarios close. This similarity might be explained by the well-designed layout of the first scenario for this DC. Overall, we conclude that joint global optimization, as is possible with our formulation, is advisable. Here, we omit the fixed layout cases (Scenario 1) from further analysis.

From the Pareto front of the second scenario, seven solutions are selected to proceed to the final decision-making stage. Five solutions with different DC and LCC lower than 300,000 million rials and a total sustainability greater than 70%—also the solution with the lowest LCC and the greatest total sustainability—are chosen from the Pareto fronts. Table 3 presents the indicator values, entropy weights, TOPSIS scores, and final solution rankings. The results show that the solution with $\mathrm{DC}=33\%$ has the greatest score (Rank 1) among all selected solutions. Fig. 5 shows this design and its specifications. Other selected solutions ($\mathrm{DC}=100\%$, 88%, 55%, 0%) are shown in Figs. 6–9. The TOPSIS scores for solutions with Ranks 1–5 are relatively close. This suggests that most of the solutions on the Pareto front can provide satisfactory quality. It must be noted that all of these conclusions are based on a specific case study and no general conclusion can therefore be made. However, the proposed framework can be applied to different case studies with completely different situations without restriction.

As mentioned before, the proposed framework enables stakeholder involvement in the complex decision-making processes. In the problem definition stage, these stakeholders can suggest their desired technologies, such as different types of GBIs, treatment facilities, or conventional systems. In the simulation-optimization stage, all possible combinations of these technologies are systematically generated and evaluated to provide a Pareto front of nondominated solutions. Finally, in the final decision-making stage, stakeholders participate in selecting and ranking solutions. They can do so either by determining the weights of the indicators or by ranking the solutions based on their preference as an extra indicator. This approach might reduce the unfavorable procedural outcomes, resistance, and conflicts that stakeholders often cause when they feel undervalued and unheard (Pigmans et al. 2019). However, this framework restricts the domain of their contribution to only feasible, previously optimized, and plausible solutions. We also calculate solution ranks without considering the stakeholders’ opinions (without including ${\mathrm{Stk}}_{\mathrm{Rank}}$ in Table 3 to calculate TOPSIS scores). The only change in comparison with the ranking in Table 3 is that the order of the first and second solutions is reversed.

Finally, the role of structural resilience in the layout configuration of the final design is investigated. Structural resilience is often neglected in the design of UDSs. However, it can play a significant role in system performance during extreme events. For instance, the fluvial flood that occurred in the study area in the winter of 2020 dramatically increased the groundwater level in the riversides, resulting in the choking and blockage of several main pipes of the existing centralized wastewater collection network. As a consequence, several parts of the network, even far away from the riversides, were out of service, causing wastewater to spill from manholes and overflow in the streets and resulting in serious disturbances and health problems for citizens.

In Fig. 5, the structural resilience of each outlet [using Eq. (8)] is presented. These values can be interpreted as the percentage of total area not affected by any failure in that outlet. As an example, the minimum structural resilience among all outlets is 79%, which belongs to Outlet 4. This means that only 21% percent of the total area is connected to that outlet. The critical pipes, which are the first pipe in each part of the system with a structural resilience less than 90%, are also shown in Fig. 5. The pipes downstream of the critical pipes have less structural resilience. However, increasing the system’s redundancy in weak points can resourcefully increase network resilience by for example, adding some loops or pipes that divert the flow direction to other parts of the network in emergency conditions. In this way, the effects of single pipe failures can be significantly restricted. Of the 530 pipes in the case study, 4 are recognized as critical, as seen in Fig. 5. Therefore, only by introducing 4 additional elements can a minimum (90% here) structural resilience be achieved for all pipes in the system. Interestingly, the structural resilience of Outlet 9 is 84%; however, none of its upstream pipes has a structural resilience of less than 90%. The reason is that, in this part of the network, the stormwater is smartly collected from three different main collectors, as depicted in Fig. 5.