# Mathematical Optimization Method of Low-Impact Development Layout in the Sponge City

Authors: Hong Men; Hao Lu; Wenjuan Jiang; Duo Xu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Aiming at the optimization layout of distributed low-impact development (LID) practices in the sponge city, a new mathematical method combining Stormwater Management Model (SWMM) and preference-inspired co-evolutionary algorithm using goal vectors (PICEA-g) was developed and was applied in the Ximen waterlogged area of Pingxiang City. Firstly, a block-scaled rainfall-runoff model was built in the study area by using SWMM. Then, an LIDs area optimization model was established by linking the SWMM and the PICEA-g based on the Matlab platform, which took the area ratios of various LIDs in each block as decision variables and took the total runoff, peak flow, suspended substance (SS) pollutant, and LIDs cost as objective functions. Thus, the problem of LIDs layout was turned into a mathematical optimization issue. So the cost-benefit optimal solutions with different emphases were found by using this algorithm, and the LIDs layout optimal scheme for this area was further analysed and verified by rainfall-runoff model. The results show that the total runoff reduction rates of the system reach a maximum of 21.8%, the peak flow reduction rates of the system are more than 10%, and the SS pollutant reduction rates are reduced by about 30% compared with before LIDs under the design storms of different return periods. The reduction rates of each runoff index are higher than the nondominated sorting genetic algorithm II (NSGA-II) method, and decision-makers can more effectively analyse the cost-benefit optimal solution from the Pareto solution sets. Therefore, the LIDs layout optimization method proposed in this paper has obvious advantages in solving similar many-objective optimization problems (MOOPs) in sponge city construction.


---

## Body

## 1. Introduction

In the past 40 years, China’s urbanization rate has increased by nearly 40% under the global trend of rapid urbanization and modernization [1]. The high degree of urbanization has changed many ecological green spaces into impermeable pavements, resulting in a weakening of urban natural ability to absorb rainwater. Once heavy rainfall exceeds the drainage capacity of the original rainwater pipe network, urban storm runoff cannot be discharged timely, which will cause road ponding and even urban waterlogging [2, 3]. In addition, the impact of human activities on the natural water cycle has generated a series of “urban water problems,” such as water pollution and water shortage, which are not conducive to the sustainable development of society and the economy [4, 5]. One of the effective methods to solve the above problems is the low-impact development proposed by the US Environmental Protection Agency (EPA), which focuses on the reasonable layout of distributed small-scale green infrastructure (such as green roof, permeable pavement, rain barrel, etc.), collects and processes rainwater from the source to achieve the purpose of controlling storm runoff [6, 7]. Based on the concept of low-impact development distributed treatment, “sponge cities” was first proposed by China in 2013, which co-ordinates the construction of sponge facilities from three aspects: source reduction, process control, and system governance [8]. It has been launched in 30 pilot cities in China, and the assessment standard for sponge city construction effect was officially implemented on August 1, 2019 [9].Researchers use storm management model tools, such as SWMM, InfoWorks CS, and MIKE URBAN, to simulate the effect of runoff control in a certain area of the city after applying some LIDs or LID combination practices and formulate rainwater system planning solutions [10–12]. In view of the distributed feature of sponge facilities, the rational layout of various LID practices has become a core issue. The general process and method of LIDs design in urban areas are shown in Figure 1. At first, the Urban Stormwater Treatment and Analysis Integration (SUSTAIN), a system to select best management practices (BMPs) developed by US EPA, is used to analyse the optimal combination of LIDs and BMPs [13]. Jia et al. [14] estimate the location and layout of LIDs BMPs taking SUSTAIN system as a planning support tool. The LIDs area ratios optimization are commonly further considered after determining the location of various kinds of LIDs. Bai et al. [15] set different LIDs layout ratios at a certain interval in different design scenarios. Environmental and economic benefits evaluation of the listed schemes is obtained by combining the analytical hierarchy process (AHP) and SWMM. Zhang et al. [16] set up four scenarios with different LIDs area ratios, construct the optimal objective function, and calculate the value of each scenario. The LIDs combination plan is obtained with the consideration of highest total runoff control rate per unit cost. Determining the ratios of LIDs layout area in a rainwater system is actually a multiobjective optimization problem (MOP). The intelligent optimization algorithms with powerful search ability are used to find the best LIDs layout scheme. Ghodsi et al. [17] propose a multievent single objective genetic algorithm, which considers the optimal solution of runoff volume and LIDs cost under various climate scenarios. According to their previous research, the weight of LIDs construction cost, runoff quantity and runoff quality could be set as 0.3, 0.35 and 0.35, respectively, when they are regarded as equally important objectives for optimization [18]. As more and more multiobjective intelligent evolution algorithms are proposed, various solving methods have been widely used in LIDs layout optimization and obtained some advances. Paola et al. [19] propose a decision Support System (DSS) based on harmony search algorithm to reduce both the flooded volumes and the expense of LIDs. Xu et al. [20] couple nondominated sorting genetic algorithm II (NSGA-II) to SWMM through Python and comprehensively consider the runoff index, pollutant load, and construction cost index to select a cost-effective solution as the final site-scale LIDs layout planning. Considering two objectives of the total cost of the LIDs and the annual runoff total control rate, Tao et al. [21] establish the multiobjective optimization mathematical model of the LIDs design in the west of Shanghai, China. The cost-benefit optimal plan of LIDs layout area is obtained by using NSGA-II.Figure 1
Process and method of LIDs design.The above researches of the LIDs optimization are generally aimed at two or three runoff evaluation objectives. However, according to the assessment standard for the effect of sponge city construction, assessment requirements have been formulated from various aspects such as runoff volume, peak runoff, and pollutant. It is necessary to consider many (more than three) objectives during LID practices design [22]. Therefore, the problem of reasonable layout of LIDs is transformed into a many-objective optimization problem (MOOP). One of the effective methods to solve the MOOP is to introduce the decision-makers’ preference information into the Pareto optimal solving process. Based on NSGA-II [23], this type of algorithm improves the Pareto domination rule by changing the fitness assignment to enhance the selection ability of the algorithm [24]. Among them, the representative algorithm is preference-inspired co-evolutionary algorithm using goal vectors (PICEA-g). Wang et al. [25] propose PICEA-g algorithm and prove that it is superior to some multiobjective evolutionary algorithms in solving MOOP. Li et al. [26] apply the PICEA-g to the optimization of combined cooling, heating, and power (CCHP) systems to optimize the system reliability, system cost, and environmental sustainability. According to the research of Ghaffari et al. [27], PICEA-g can also be applied to the evaluation of complex decision-making mode in project management.In addition, the study area of this paper is located in the old urban area of Pingxiang City, which has the problems of dense population and old pipe network facilities [28]. It is urgent to reduce the pressure of urban drainage pipe network and alleviate the social and economic problems caused by pipe network reconstruction through LIDs optimization [29]. Considering the above analysis, we establish an optimization model of LIDs layout by using the Matlab platform to combine SWMM with a many-objective optimization algorithm. The optimization model aims to reduce the total runoff, peak flow, SS pollution, and the cost of LIDs, and the PICEA-g algorithm is used to solve the optimal LIDs area ratio in each block.


## 2. Materials and Methods

Taking the Ximen waterlogged area of Pingxiang City as an example, a block-scaled many-objective optimization model of LIDs layout was established, where SWMM is used as a tool for calculating the optimization objective functions. Combining SWMM with a many-objective optimization algorithm, the layout ratio of different LIDs combinations in the urban area under the comprehensive consideration of multiple sponge city indicators can be optimized to achieve the all indicators as optimal as possible.


### 2.1. Research Area Modeling

#### 2.1.1. Case Study Area

Pingxiang City, Jiangxi Province, as one of the pilot cities of the sponge city in China, is known as the “Pearl of Ganxi” and is the center of economic development in Ganxi area. Pingxiang City has a long history and is one of the red revolutionary bases in China. It has gone through multiple stages of development. Especially in recent years, the process of urbanization has accelerated. But at the same time, high urbanization has also brought a series of water problems, such as urban waterlogging and nonpoint source pollution [30]. Since China proposed the construction of the sponge city in 2013, with the country’s attention and resource input, the public’s cooperation and recognition of the “sponge city” has increased gradually, which is more conducive to the construction and development of the sponge city.Pingxiang City is located in a subtropical monsoon humid climate zone. According to the rainfall data of Pingxiang Meteorological Bureau for many years (1955∼2009), the average annual rainfall is as high as 1596.7 mm. Most of the rainfall is concentrated in summer, and it is dry and less rainy when summer and autumn alternate. Recently, due to the impact of global climate change, extreme rainfall events have occurred from time to time [31]. The soil in Pingxiang City is mostly clay soil with poor permeability, which is not conducive to timely infiltration of rainwater. The Ximen waterlogged area studied in this paper is located in the southwest of the old city (As shown in Figure 2). The Pingshui River flows from north to south in the east and south, which accepts rainwater and sewage discharge in the city. The geomorphology of the study area is relatively complex, mostly hilly terrain. The ground elevation ranges from 89.12 to 106.45, showing high in the north and low in the south, and high in the east and low in the west. The total area of the study area is 94 ha, where there are densely distributed building units such as old communities, government units, schools, and commercial districts. The planning of green gardens is inadequate; especially the old communities have almost no green space.Figure 2
Location of Ximen waterlogged area of Pingxiang City.


#### 2.1.2. Rainfall-Runoff Model

Before setting up a SWMM model of the study area, it is necessary to collect the underlying surface, the pipe network, rainfall, and other data of the study area and import them into the SWMM. According to the percent of impervious area, the distribution of blocks, and the flow direction of the pipeline, the study area is divided into 59 subcatchments. In order to simplify the model running process, the drainage pipe network is reasonably generalized into 58 conduits and 4 outfalls. The generalization results are shown in Figure3. The setting method of hydrology, hydraulics, water quality parameters, and rain gauge of each module are, respectively, as follows.Figure 3
Generalization of the study area.When setting up a new model in SWMM, the project defaults could be set in advance. There are two preset options in this simulation, which are infiltration model and routing model. For the former, three main models of SWMM can simulate the infiltration process: Horton, Green-Ampt, and Curve number (CN). Among them, the principle of Horton model is to compare rainfall intensity with changing soil permeability. The core idea is to ensure that the actual infiltration rate is the smaller one of the rainfall intensity and the soil infiltration rate during the entire infiltration process [32]. This method requires only three input parameters, which are determined by the soil characteristics of study area and are easy to obtain. Therefore, this paper chose the Horton method which is widely used in the establishment of SWMM model to simulate the infiltration process. Parameters setting is shown in Table 1. For the other preset options, according to the different methods of solving the flow state in the conduits, there are three methods for simulating the flow of the conduits: Steady flow, Kinematic wave, and Dynamic wave. The simulation in this paper chooses the dynamic wave method that can calculate the complex flow state, such as pressure flow and backflow [33].Table 1
Parameters setting of model hydrological and hydraulic.
Parameters
Description
Values
Max. infil. rate
Maximum rate on the horton infiltration curve (mm/hr)
78.1
Min. infil. rate
Minimum rate on the horton infiltration curve (mm/hr)
3.82
Decay constant
Decay constant for the horton infiltration curve (1/hr)
2
N-imperv
Mannings N for impervious area
0.013
N-perv
Mannings N for pervious area
0.15
Dstore-imperv
Depth of depression storage on impervious area (mm)
3.15
Dstore-perv
Depth of depression storage on pervious area (mm)
5.08
%Zero-imperv
Percent of impervious area with no depression storage (%)
25
Conduit roughness
Manning’s roughness coefficient
0.015According to the method of parameters setting, the main hydrological and hydraulic parameters of the model can be divided into two categories.The first category is parameters that can be measured or calculated from actual data. (1) Area: by taking the 1: 200 satellite map of the study area as the background and combining the automatic calculation function of SWMM, the area of each subcatchment area was obtained. (2) Slope: from the collected elevation data, calculate the ratio of the elevation difference to the horizontal distance of the runoff paths to obtain the average surface slope of each subcatchment. (3) Percent of impervious area: according to the land use map (as shown in Figure4) in the study area, the land use types can be divided into buildings, roads, and green spaces. The runoff coefficients of them are specified as 0.75, 0.85, and 0.25, and they are used as weight to calculate the final percentage of impervious area in the subcatchment. (4) Properties of the drainage conduits and juctions. These parameters are determined by the information provided by the Pingxiang Surveying and Mapping Institute. The above parameters are presented in Tables 2–4.Figure 4
Land use map of the study area.Table 2
Subcatchment parameters.
Name
Area (ha)
%Imperv
Width (m)
%Slope
S1
0.4418
55.5
66.468
0.36
S2
0.1934
78.83
43.977
0.36
S3
0.4349
77.7
65.947
0.38
S4
1.6089
56
126.842
0.38
S5
1.61
55.7
126.886
0.38
S6
0.8104
79.5
90.022
1.5
S7
1
44
100
1
S8
1
55.7
100
1
S9
3.215
51.5
179.304
1
S10
2
46
141.421
0.35
S11
1.801
55.5
134.201
0.2
S12
0.81
58.5
90
0.2
S13
1.801
58.5
134.201
0.13
S14
0.5171
83.01
71.91
0.13
S15
1.79
49.5
133.791
0.13
S16
0.7
52
83.666
0.4
S17
2.2251
53.5
205.55
0.38
S18
0.9606
69.99
98.01
0.38
S19
2
55.5
141.421
0.2
S20
2.32
64
152.315
0.4
S21
1.4408
69.99
120.033
0.38
S22
0.68
56
82.462
0.4
S23
0.68
53.5
82.462
0.5
S24
0.2308
76.03
48.042
0.11
S25
1.5758
75.79
125.531
0.11
S26
0.2426
80.95
49.254
0.5
S27
0.2656
76.73
51.536
1.1
S28
2.5989
46
161.211
0.5
S29
2.58
41.5
160.624
0.11
S30
1.5821
53.5
125.782
0.38
S31
1.58
53.5
125.698
0.38
S32
2.3336
53.5
152.761
0.5
S33
0.89
54
94.34
0.5
S34
3.54
64
188.149
0.4
S35
1.1323
60.87
106.41
0.4
S36
0.358
71.5
59.833
5.04
S37
10.8248
75.87
329.01
0.1
S38
0.358
56
59.833
0.4
S39
5
48.5
223.607
0.1
S40
1.5779
53.5
125.614
0.5
S41
1.2224
80.83
110.562
0.11
S42
0.7574
82.52
87.029
1
S43
0.2
64.5
44.721
0.5
S44
0.1047
78.86
32.357
1.18
S45
0.4356
56.5
66
1.18
S46
0.1677
80.06
40.951
0.03
S47
2.335
58.5
115.542
0.5
S48
0.33
54.7
57.446
1
S49
1.2142
79.42
110.191
1.1
S50
1.335
53
115.542
1.1
S51
0.3358
71.49
57.948
0.23
S52
0.305
77.67
55.227
0.5
S53
0.45
54.7
67.082
0.5
S54
7.59
44.7
275.5
0.5
S55
0.4058
77.52
63.702
0.5
S56
1.45
54.7
120.416
0.23
S57
3.65
43.5
191.05
0.5
S58
2.33
44.7
115.326
0.5
S59
1.45
43.5
120.416
0.5Table 3
Juction parameters.
Name
Elevation (m)
J1
103.981
J2
102.89
J3
101.918
J4
98.102
J5
95.115
J6
93.913
J7
103.71
J8
95.115
J9
95.115
J10
98.102
J11
101.6
J12
100.001
J13
102.226
J14
101.6
J15
100.001
J16
98.889
J17
103.71
J18
102.226
J19
101.6
J20
104.505
J21
109.424
J22
108.338
J23
102
J24
105.356
J25
105.803
J26
108.719
J27
107
J28
98.102
J29
104.506
J30
102.226
J31
102.226
J32
101.001
J33
98.102
J34
95.115
J35
106.003
J36
106.003
J37
104.905
J38
103.71
J39
98.898
J40
95.115
J41
92.713
J42
90.443
J43
92.713
J44
94.5
J45
98.102
J46
95.115
J47
95.368
J48
91.908
J49
88.341
J50
87.5
J51
86.287
J52
85.215
J53
84.143
J54
86.287
J55
85.215
J56
88.341
J57
97.87
J58
107.126
PFK3
83.034
PFK2
89.190
PFK4
90.443
PFK1
97.300Table 4
Conduits parameters.
Name
From node
To node
Length (m)
pipe diameter (m)
1
J1
J2
120
0.6
2
J2
J3
120
0.6
3
J3
J4
120
0.6
4
J4
J5
120
1
5
J5
J6
180
1.2
7
J7
J13
120
1
11
J14
J15
150
1.2
12
J15
J16
150
2.4
13
J10
J9
90
1
14
J9
J6
180
1.2
15
J8
J6
130
1.2
16
J28
J8
90
1
17
J17
J18
120
1
18
J18
J19
140
1.2
19
J20
J17
120
1
21
J26
J27
120
2
22
J27
J24
150
1.2
23
J25
J24
150
1.2
24
J24
J23
180
2
25
J35
J37
120
0.8
26
J36
J37
120
1
27
J37
J38
120
1
28
J29
J17
120
1
29
J30
J19
120
1
30
J33
J34
120
1
31
J45
J46
120
1
32
J46
J44
120
1
33
J44
J42
120
1.2
34
J34
J43
150
1
35
J43
J42
150
1.2
36
J23
J39
150
2
38
J40
J41
100
1
39
J47
J48
180
2
40
J48
J49
180
2
41
J49
J50
180
2.4
42
J50
J51
180
2.4
43
J51
J52
180
2.4
44
J52
J53
200
2.4
45
J54
J55
120
2
46
J55
J53
120
2
47
J53
PFK3
180
2.4
48
J42
PFK2
150
1.2
49
J6
PFK4
180
2
50
J39
J57
150
2
51
J57
J47
150
2
52
J41
J42
120
1.2
53
J56
J50
180
2.4
54
J16
PFK1
160
2.4
56
J21
J22
120
0.8
57
J22
J58
120
0.8
58
J58
J35
120
0.8
59
J38
J31
160
1.2
60
J31
J32
180
1.2
C64
J13
J11
120
1
C65
J11
J12
120
1.2
C66
J12
J16
120
1.2
C69
J19
J15
150
1.2
C70
J32
J15
180
1.2The second category is parameters that refer to the SWMM User’s Manual or other References. (1) The values of N for impervious (pervious) area and depth of depression storage on impervious (pervious) area were determined according to different types of underlying surface by referring to the SWMM User’s Manual.(2) Percent of impervious area with no depression storage is generally taken as 25% [34]. (3) Most of the conduits are concrete and the roughness coefficients are set by referring to the user’s manual. The above approximate empirical parameters are also shown in Table 1.When it comes to water quality simulation, suspended substance (SS) is the main physical pollutant and most of the pollutants are positively correlated with SS, so it can represent the regional pollution situation [35]. In this simulation, the water quality of the study area was reflected by the total amount of SS. Referring to the relative articles, the concentration of pollutant in rainwater was set to 10 mg/L [36]. Because the existing pollutant monitoring data is incomplete, in this simulation, the initial pollutant buildup on subcatchment is set to zero. The parameters setting for simulating pollutant buildup and washoff in different types of land use are shown in Table 5. One thing to mention is that this model did not perform parameter verification due to the lack of actual monitoring data of pipe network flow and pollutant concentration changes under real rainfall events, and only an approximate model of the study area was constructed.Table 5
Parameters setting of model pollutant buildup and washoff
Buildup (saturation function)
Washoff (exponential function)
Maximum buildup (kg)
Rate constant
Saturation constant (day)
Coefficient
Index
Cleaning efficiency (%)
Building
140
0
4
0.007
1.8
0
Road
270
0
4
0.008
1.8
70
Green space
60
0
4
0.004
1.2
0In order to study the effect of LID practices on runoff control in different rainfall situations, especially short-duration heavy rainfall, design storms with different return periods were used to replace actual precipitation events [37]. The storm intensity formula in Pingxiang City is.(1)q=26191+0.78lgPt+100.79,where q is the average rainfall intensity (L/s·ha), t is the rainfall duration (min), and P is the rainfall return period (year).The commonly used methods of rainfall pattern decomposition include uniform method and nonuniform method [38]. The uniform method is relatively simple, but its design storm results are relatively small, which does not conform to the actual rainfall law [39]. The nonuniform method can be divided into Chicago rainfall pattern, Huff rainfall pattern, and triangular rainfall pattern according to the different design approaches of the rain peak [40]. This study chose the Chicago Law, which is widely applicable to short-duration rainfall and recommended by China’s “Code for Design of Outdoor Wastewater engineering.” By setting the rain peak 0.4 and the return periods 2-year, 10-year and 50-year, the design storms that last for 2 hours are generated and shown in Figure 5.Figure 5
Design rainfall process lines for different return periods.


#### 2.1.3. LIDs Layout Planning

Deploying LIDs control in the SWMM model requires two steps: first, set up several sets of LIDs design regardless of size, and then assign any required design combinations and design sizes to the selected subcatchment. SWMM5.1.13 has 8 kinds of LIDs to choose from. The LIDs selected for this simulation in the study area are vegetative swale (VS), permeable pavement (PP), green roof (GR), and rain barrel (RB). Among them, the vegetative swale is a shallow swale with herbs on the surface, which transmits rainwater runoff and has both infiltration and filtration functions [41]. It can trap some rainwater and particulate pollutants to reduce runoff discharge. Permeable pavements include permeable concrete and pavement brick paving in sidewalk, and colored permeable pavement in parking lots, which can quickly infiltrate direct rainfall or runoff formed by impervious pavements [42]. However, the rainwater that cannot infiltrate to permeable pavement in time overflows into the vegetative swale for infiltration and filtering and finally flows into the rainwater pipe network. As a low-impact development measure for buildings, green roof and rain barrel can be selected based on the age and type of building. Generally speaking, the older buildings are not suitable for setting green roof considering the load bearing capacity of roof [43]. The layout of rain barrels are more flexible and can be set on the side of buildings such as houses, factories, and schools, and rainwater stored in it could be used to irrigate green space or be used as domestic water for residents. The parameter settings of the above four types of LIDs in the model referred to the general LID facilities technical specifications in China [44] and was fine-tuned according to the actual deployment environment. The specific parameters setting of LID practices are shown in Table 6.Table 6
Parameters setting of LID practices.
LIDs type
Construction
Parameter option
Parameter value
Permeable pavement
Surface
Berm height (mm)
3
Surface roughness (Mannings N)
0.012
Surface slope (%)
1
Pavement
Thickness (mm)
120
Void ratio (voids/Solids)
0.2
Permeability (mm/hr)
2450
Storage
Thickness (mm)
230
Void ratio (voids/Solids)
0.75
Seepage rate (mm/hr)
300
Green roof
Surface
Berm height (mm)
6.35
Vegetation volume fraction
0.85
Surface roughness (Mannings N)
0.23
Surface slope (%)
1
Soil
Thickness (mm)
600
Porosity (volume fraction)
0.5
Conductivity (mm/hr)
12.7
Drainage mat
Thickness (mm)
200
Void fraction
0.75
Roughness (Mannings N)
0.014
Vegetative swale
Surface
Berm height (mm)
6
Vegetation volume fraction
0.85
Surface roughness (Mannings N)
0.24
Surface slope (%)
1
Swale side slope (run/rise)
7.33
Rain barrel
Storage
Barrel height (mm)
500With the help of siting tool in System of SUSTAIN, the layout locations of the four LIDs were evaluated and a total of 20 blocks were determined. Then, the LIDs combination schemes in each block were obtained as shown in Figure6. Then, it is necessary to determine the way of the LIDs deployment in each subcatchment of the model. There are two different ways: the first way is to add one or more LID practices, wherein each LID controls the runoff generated by the impervious zone (or permeate zone) at a specific ratio in the subcatchment. This way is suitable for large research area, and different combinations of LIDs can be deployed in subcatchments with multiple land use types [45]. The second way is to decompose the original subcatchments into several small catchment areas with separate LIDs control and underlying surface type. This deployment mode is also called LID series mode, that is, the outflow of one LID can be the inflow of another LID, which is suitable for small area or a situation where a detailed analysis of certain LID facility is required [34]. In this simulation, the second LIDs layout was selected. According to the LIDs combination in each subcatchment, the original subcatchments were decomposed into several small catchment areas with a single LID facilities or land use type. Figure 7 shows a schematic diagram of one LIDs deployment way. Among them, the permeable pavement can handle runoff from the hard surface and direct rainfall. Green roof (rain barrel) handles rainwater runoff from the building roof (rainwater stored in the rain barrel can be used to water green spaces). As a transmission processing unit, the vegetative swale receives the runoff from the original green space, permeable pavement, and green roof before the runoff discharged into the outlet. After formulating the combined deployment plan of the LID practices in each subcatchment, how to design the layout size of the LIDs in each subcatchment to achieve the best cost-effective is the issue to be discussed next.Figure 6
LIDs combination schemes in each block to be transformed.Figure 7
Schematic diagram of LIDs deployment.


### 2.2. Optimization Model

#### 2.2.1. PICEA-g Algorithm

Preference-inspired co-evolutionary algorithm (PICEA) was first proposed in 2013 to solve many-objective optimization problem (MOOP). This algorithm introduces several random preference sets and co-evolves with candidate solution sets, so that it can make candidate solutions approach the Pareto frontier while looking for decision-makers’ preference areas [46]. On this basis, Wang et al. [25] propose preference-inspired co-evolutionary algorithm using goal vectors (PICEA-g), which uses the goal vectors as preference sets. When solving MOOP, fitness allocation rules based on the original Pareto dominance relationship are added to enhance the ability of selecting individuals in the algorithm. It is verified by the author that the algorithm performs are better than the widely used NSGA-II [23], ε-MOEA [47], and other multiobjective evolutionary algorithms when solving many-objective (three or more) optimization problems. The core of the PICEA-g algorithm is to specify the fitness allocation rules for the goal vectors and the candidate solutions, so that the Pareto-dominated solutions can be retained as the parent population of the next generation [48]. The fitness calculation formulas of candidate solution Fs and goal vector Fg are [25](2)Fs=0+∑g∈G∪GCs≤g1ng,Fg=11+α,(3)α=1,ng=0,ng−12N−1,others,where G and Gc represent the parent and offspring populations of the goal vector g, respectively, ng is the number that satisfies the goal vector g, and N is the population size. If the candidate solution s does not satisfy any goal vector g, then the fitness value Fs of s is taken as 0.For the sponge city low-impact development and construction-planning field, it can be essentially attributed to the MOOP of complex systems. Combined with the superiority of the PICEA-g algorithm, this paper performed many-objective optimization of LID practices planning in order to find a more cost-effective layout solution efficiently.


#### 2.2.2. Decision Variables

The blocks in the study area can be divided into four types: administration (A), business (B), green space or square (G), residence (R), and the LIDs combinations deployed in each of the blocks to be transformed are different. Taking the layout ratios of LIDs in each subcatchment as the optimization decision variables, there are 58 decision variables in 20 blocks, and the specific form of decision variableD is shown in(4)D=Num,Type,where Num represents the serial number of the blocks, Type represents four LID practice types: PP, GR, VS, and RB. According to the requirements for the layout of LIDs in the “Technical Guide for Sponge City Construction (Trial)” and the investigation of research area [49], the maximum LIDs layout ratios on the corresponding underlying surface in different types of blocks are shown in Table 7.Table 7
Maximum LIDs layout ratios for different types of blocks.
Land classification
Maximum LIDs layout ratios
PP (%)
GR (%)
VS (%)
A
60
95
86
B
65
95
86
G
75
—
95
R
45
95
75


#### 2.2.3. Objective Functions

For assessing the effect of LID practices on the original urban drainage system, it is generally considered from several aspects: the effect of water volume control, the degree of water quality improvement, and the construction cost. In this study, four objective functions that can evaluate the cost-benefit of sponge city construction are selected: total runoff, which represents the sum of surface runoff for each subcatchment; peak flow, which indicates the maximum discharge flow at the system outfall; the SS pollutant is the total mass of SS discharged from the outlet; the total cost of the LID facilities includes construction costs but excludes other costs such as maintenance. The first three objective function values can be extracted from the report file of the SWMM, and the calculation formula for the total costCsum function of the LIDs can be expressed as(5)Csum=∑i=120∑j=13Sij×Cj+∑i=120Di4×C4,where Sij represents the layout ratio of a certain LIDs in a certain block, Cj represents LIDs (PP, GR or VS) construction cost per unit area, Di4 indicates the number of rain barrels in a certain block, and C4 represents the construction cost of RB per unit (the construction cost of the unit LID facility is shown in Table 8).Table 8
Construction cost of unit LID practices [50, 51].
LIDs type
PP (CNY/m2)
GR (CNY/m2)
VS (CNY/m2)
RB (CNY/each)
Construction cost
193
500
220
350


#### 2.2.4. Connection of SWMM and Optimization Model

Because the project file of the SWMM is stored in the form of texts as a configuration file, an input file, an output file, and a report file. We first use Matlab to call the text, update the values of the various decision variables, and related parameters generated by the compilation algorithm and call the dynamic link library files of SWMM. Then, by calling the calculation engine, the calculation result of the model which the user interests is read for the calculation of the objective function. In this way, the interactive connection between the SWMM model and the PICEA-g many-objective optimization model is completed, and the specific calculation process is shown in Figure8.Figure 8
Joint calculation process of SWMM and PICEA-g algorithm.


## 2.1. Research Area Modeling

### 2.1.1. Case Study Area

Pingxiang City, Jiangxi Province, as one of the pilot cities of the sponge city in China, is known as the “Pearl of Ganxi” and is the center of economic development in Ganxi area. Pingxiang City has a long history and is one of the red revolutionary bases in China. It has gone through multiple stages of development. Especially in recent years, the process of urbanization has accelerated. But at the same time, high urbanization has also brought a series of water problems, such as urban waterlogging and nonpoint source pollution [30]. Since China proposed the construction of the sponge city in 2013, with the country’s attention and resource input, the public’s cooperation and recognition of the “sponge city” has increased gradually, which is more conducive to the construction and development of the sponge city.Pingxiang City is located in a subtropical monsoon humid climate zone. According to the rainfall data of Pingxiang Meteorological Bureau for many years (1955∼2009), the average annual rainfall is as high as 1596.7 mm. Most of the rainfall is concentrated in summer, and it is dry and less rainy when summer and autumn alternate. Recently, due to the impact of global climate change, extreme rainfall events have occurred from time to time [31]. The soil in Pingxiang City is mostly clay soil with poor permeability, which is not conducive to timely infiltration of rainwater. The Ximen waterlogged area studied in this paper is located in the southwest of the old city (As shown in Figure 2). The Pingshui River flows from north to south in the east and south, which accepts rainwater and sewage discharge in the city. The geomorphology of the study area is relatively complex, mostly hilly terrain. The ground elevation ranges from 89.12 to 106.45, showing high in the north and low in the south, and high in the east and low in the west. The total area of the study area is 94 ha, where there are densely distributed building units such as old communities, government units, schools, and commercial districts. The planning of green gardens is inadequate; especially the old communities have almost no green space.Figure 2
Location of Ximen waterlogged area of Pingxiang City.


### 2.1.2. Rainfall-Runoff Model

Before setting up a SWMM model of the study area, it is necessary to collect the underlying surface, the pipe network, rainfall, and other data of the study area and import them into the SWMM. According to the percent of impervious area, the distribution of blocks, and the flow direction of the pipeline, the study area is divided into 59 subcatchments. In order to simplify the model running process, the drainage pipe network is reasonably generalized into 58 conduits and 4 outfalls. The generalization results are shown in Figure3. The setting method of hydrology, hydraulics, water quality parameters, and rain gauge of each module are, respectively, as follows.Figure 3
Generalization of the study area.When setting up a new model in SWMM, the project defaults could be set in advance. There are two preset options in this simulation, which are infiltration model and routing model. For the former, three main models of SWMM can simulate the infiltration process: Horton, Green-Ampt, and Curve number (CN). Among them, the principle of Horton model is to compare rainfall intensity with changing soil permeability. The core idea is to ensure that the actual infiltration rate is the smaller one of the rainfall intensity and the soil infiltration rate during the entire infiltration process [32]. This method requires only three input parameters, which are determined by the soil characteristics of study area and are easy to obtain. Therefore, this paper chose the Horton method which is widely used in the establishment of SWMM model to simulate the infiltration process. Parameters setting is shown in Table 1. For the other preset options, according to the different methods of solving the flow state in the conduits, there are three methods for simulating the flow of the conduits: Steady flow, Kinematic wave, and Dynamic wave. The simulation in this paper chooses the dynamic wave method that can calculate the complex flow state, such as pressure flow and backflow [33].Table 1
Parameters setting of model hydrological and hydraulic.
Parameters
Description
Values
Max. infil. rate
Maximum rate on the horton infiltration curve (mm/hr)
78.1
Min. infil. rate
Minimum rate on the horton infiltration curve (mm/hr)
3.82
Decay constant
Decay constant for the horton infiltration curve (1/hr)
2
N-imperv
Mannings N for impervious area
0.013
N-perv
Mannings N for pervious area
0.15
Dstore-imperv
Depth of depression storage on impervious area (mm)
3.15
Dstore-perv
Depth of depression storage on pervious area (mm)
5.08
%Zero-imperv
Percent of impervious area with no depression storage (%)
25
Conduit roughness
Manning’s roughness coefficient
0.015According to the method of parameters setting, the main hydrological and hydraulic parameters of the model can be divided into two categories.The first category is parameters that can be measured or calculated from actual data. (1) Area: by taking the 1: 200 satellite map of the study area as the background and combining the automatic calculation function of SWMM, the area of each subcatchment area was obtained. (2) Slope: from the collected elevation data, calculate the ratio of the elevation difference to the horizontal distance of the runoff paths to obtain the average surface slope of each subcatchment. (3) Percent of impervious area: according to the land use map (as shown in Figure4) in the study area, the land use types can be divided into buildings, roads, and green spaces. The runoff coefficients of them are specified as 0.75, 0.85, and 0.25, and they are used as weight to calculate the final percentage of impervious area in the subcatchment. (4) Properties of the drainage conduits and juctions. These parameters are determined by the information provided by the Pingxiang Surveying and Mapping Institute. The above parameters are presented in Tables 2–4.Figure 4
Land use map of the study area.Table 2
Subcatchment parameters.
Name
Area (ha)
%Imperv
Width (m)
%Slope
S1
0.4418
55.5
66.468
0.36
S2
0.1934
78.83
43.977
0.36
S3
0.4349
77.7
65.947
0.38
S4
1.6089
56
126.842
0.38
S5
1.61
55.7
126.886
0.38
S6
0.8104
79.5
90.022
1.5
S7
1
44
100
1
S8
1
55.7
100
1
S9
3.215
51.5
179.304
1
S10
2
46
141.421
0.35
S11
1.801
55.5
134.201
0.2
S12
0.81
58.5
90
0.2
S13
1.801
58.5
134.201
0.13
S14
0.5171
83.01
71.91
0.13
S15
1.79
49.5
133.791
0.13
S16
0.7
52
83.666
0.4
S17
2.2251
53.5
205.55
0.38
S18
0.9606
69.99
98.01
0.38
S19
2
55.5
141.421
0.2
S20
2.32
64
152.315
0.4
S21
1.4408
69.99
120.033
0.38
S22
0.68
56
82.462
0.4
S23
0.68
53.5
82.462
0.5
S24
0.2308
76.03
48.042
0.11
S25
1.5758
75.79
125.531
0.11
S26
0.2426
80.95
49.254
0.5
S27
0.2656
76.73
51.536
1.1
S28
2.5989
46
161.211
0.5
S29
2.58
41.5
160.624
0.11
S30
1.5821
53.5
125.782
0.38
S31
1.58
53.5
125.698
0.38
S32
2.3336
53.5
152.761
0.5
S33
0.89
54
94.34
0.5
S34
3.54
64
188.149
0.4
S35
1.1323
60.87
106.41
0.4
S36
0.358
71.5
59.833
5.04
S37
10.8248
75.87
329.01
0.1
S38
0.358
56
59.833
0.4
S39
5
48.5
223.607
0.1
S40
1.5779
53.5
125.614
0.5
S41
1.2224
80.83
110.562
0.11
S42
0.7574
82.52
87.029
1
S43
0.2
64.5
44.721
0.5
S44
0.1047
78.86
32.357
1.18
S45
0.4356
56.5
66
1.18
S46
0.1677
80.06
40.951
0.03
S47
2.335
58.5
115.542
0.5
S48
0.33
54.7
57.446
1
S49
1.2142
79.42
110.191
1.1
S50
1.335
53
115.542
1.1
S51
0.3358
71.49
57.948
0.23
S52
0.305
77.67
55.227
0.5
S53
0.45
54.7
67.082
0.5
S54
7.59
44.7
275.5
0.5
S55
0.4058
77.52
63.702
0.5
S56
1.45
54.7
120.416
0.23
S57
3.65
43.5
191.05
0.5
S58
2.33
44.7
115.326
0.5
S59
1.45
43.5
120.416
0.5Table 3
Juction parameters.
Name
Elevation (m)
J1
103.981
J2
102.89
J3
101.918
J4
98.102
J5
95.115
J6
93.913
J7
103.71
J8
95.115
J9
95.115
J10
98.102
J11
101.6
J12
100.001
J13
102.226
J14
101.6
J15
100.001
J16
98.889
J17
103.71
J18
102.226
J19
101.6
J20
104.505
J21
109.424
J22
108.338
J23
102
J24
105.356
J25
105.803
J26
108.719
J27
107
J28
98.102
J29
104.506
J30
102.226
J31
102.226
J32
101.001
J33
98.102
J34
95.115
J35
106.003
J36
106.003
J37
104.905
J38
103.71
J39
98.898
J40
95.115
J41
92.713
J42
90.443
J43
92.713
J44
94.5
J45
98.102
J46
95.115
J47
95.368
J48
91.908
J49
88.341
J50
87.5
J51
86.287
J52
85.215
J53
84.143
J54
86.287
J55
85.215
J56
88.341
J57
97.87
J58
107.126
PFK3
83.034
PFK2
89.190
PFK4
90.443
PFK1
97.300Table 4
Conduits parameters.
Name
From node
To node
Length (m)
pipe diameter (m)
1
J1
J2
120
0.6
2
J2
J3
120
0.6
3
J3
J4
120
0.6
4
J4
J5
120
1
5
J5
J6
180
1.2
7
J7
J13
120
1
11
J14
J15
150
1.2
12
J15
J16
150
2.4
13
J10
J9
90
1
14
J9
J6
180
1.2
15
J8
J6
130
1.2
16
J28
J8
90
1
17
J17
J18
120
1
18
J18
J19
140
1.2
19
J20
J17
120
1
21
J26
J27
120
2
22
J27
J24
150
1.2
23
J25
J24
150
1.2
24
J24
J23
180
2
25
J35
J37
120
0.8
26
J36
J37
120
1
27
J37
J38
120
1
28
J29
J17
120
1
29
J30
J19
120
1
30
J33
J34
120
1
31
J45
J46
120
1
32
J46
J44
120
1
33
J44
J42
120
1.2
34
J34
J43
150
1
35
J43
J42
150
1.2
36
J23
J39
150
2
38
J40
J41
100
1
39
J47
J48
180
2
40
J48
J49
180
2
41
J49
J50
180
2.4
42
J50
J51
180
2.4
43
J51
J52
180
2.4
44
J52
J53
200
2.4
45
J54
J55
120
2
46
J55
J53
120
2
47
J53
PFK3
180
2.4
48
J42
PFK2
150
1.2
49
J6
PFK4
180
2
50
J39
J57
150
2
51
J57
J47
150
2
52
J41
J42
120
1.2
53
J56
J50
180
2.4
54
J16
PFK1
160
2.4
56
J21
J22
120
0.8
57
J22
J58
120
0.8
58
J58
J35
120
0.8
59
J38
J31
160
1.2
60
J31
J32
180
1.2
C64
J13
J11
120
1
C65
J11
J12
120
1.2
C66
J12
J16
120
1.2
C69
J19
J15
150
1.2
C70
J32
J15
180
1.2The second category is parameters that refer to the SWMM User’s Manual or other References. (1) The values of N for impervious (pervious) area and depth of depression storage on impervious (pervious) area were determined according to different types of underlying surface by referring to the SWMM User’s Manual.(2) Percent of impervious area with no depression storage is generally taken as 25% [34]. (3) Most of the conduits are concrete and the roughness coefficients are set by referring to the user’s manual. The above approximate empirical parameters are also shown in Table 1.When it comes to water quality simulation, suspended substance (SS) is the main physical pollutant and most of the pollutants are positively correlated with SS, so it can represent the regional pollution situation [35]. In this simulation, the water quality of the study area was reflected by the total amount of SS. Referring to the relative articles, the concentration of pollutant in rainwater was set to 10 mg/L [36]. Because the existing pollutant monitoring data is incomplete, in this simulation, the initial pollutant buildup on subcatchment is set to zero. The parameters setting for simulating pollutant buildup and washoff in different types of land use are shown in Table 5. One thing to mention is that this model did not perform parameter verification due to the lack of actual monitoring data of pipe network flow and pollutant concentration changes under real rainfall events, and only an approximate model of the study area was constructed.Table 5
Parameters setting of model pollutant buildup and washoff
Buildup (saturation function)
Washoff (exponential function)
Maximum buildup (kg)
Rate constant
Saturation constant (day)
Coefficient
Index
Cleaning efficiency (%)
Building
140
0
4
0.007
1.8
0
Road
270
0
4
0.008
1.8
70
Green space
60
0
4
0.004
1.2
0In order to study the effect of LID practices on runoff control in different rainfall situations, especially short-duration heavy rainfall, design storms with different return periods were used to replace actual precipitation events [37]. The storm intensity formula in Pingxiang City is.(1)q=26191+0.78lgPt+100.79,where q is the average rainfall intensity (L/s·ha), t is the rainfall duration (min), and P is the rainfall return period (year).The commonly used methods of rainfall pattern decomposition include uniform method and nonuniform method [38]. The uniform method is relatively simple, but its design storm results are relatively small, which does not conform to the actual rainfall law [39]. The nonuniform method can be divided into Chicago rainfall pattern, Huff rainfall pattern, and triangular rainfall pattern according to the different design approaches of the rain peak [40]. This study chose the Chicago Law, which is widely applicable to short-duration rainfall and recommended by China’s “Code for Design of Outdoor Wastewater engineering.” By setting the rain peak 0.4 and the return periods 2-year, 10-year and 50-year, the design storms that last for 2 hours are generated and shown in Figure 5.Figure 5
Design rainfall process lines for different return periods.


### 2.1.3. LIDs Layout Planning

Deploying LIDs control in the SWMM model requires two steps: first, set up several sets of LIDs design regardless of size, and then assign any required design combinations and design sizes to the selected subcatchment. SWMM5.1.13 has 8 kinds of LIDs to choose from. The LIDs selected for this simulation in the study area are vegetative swale (VS), permeable pavement (PP), green roof (GR), and rain barrel (RB). Among them, the vegetative swale is a shallow swale with herbs on the surface, which transmits rainwater runoff and has both infiltration and filtration functions [41]. It can trap some rainwater and particulate pollutants to reduce runoff discharge. Permeable pavements include permeable concrete and pavement brick paving in sidewalk, and colored permeable pavement in parking lots, which can quickly infiltrate direct rainfall or runoff formed by impervious pavements [42]. However, the rainwater that cannot infiltrate to permeable pavement in time overflows into the vegetative swale for infiltration and filtering and finally flows into the rainwater pipe network. As a low-impact development measure for buildings, green roof and rain barrel can be selected based on the age and type of building. Generally speaking, the older buildings are not suitable for setting green roof considering the load bearing capacity of roof [43]. The layout of rain barrels are more flexible and can be set on the side of buildings such as houses, factories, and schools, and rainwater stored in it could be used to irrigate green space or be used as domestic water for residents. The parameter settings of the above four types of LIDs in the model referred to the general LID facilities technical specifications in China [44] and was fine-tuned according to the actual deployment environment. The specific parameters setting of LID practices are shown in Table 6.Table 6
Parameters setting of LID practices.
LIDs type
Construction
Parameter option
Parameter value
Permeable pavement
Surface
Berm height (mm)
3
Surface roughness (Mannings N)
0.012
Surface slope (%)
1
Pavement
Thickness (mm)
120
Void ratio (voids/Solids)
0.2
Permeability (mm/hr)
2450
Storage
Thickness (mm)
230
Void ratio (voids/Solids)
0.75
Seepage rate (mm/hr)
300
Green roof
Surface
Berm height (mm)
6.35
Vegetation volume fraction
0.85
Surface roughness (Mannings N)
0.23
Surface slope (%)
1
Soil
Thickness (mm)
600
Porosity (volume fraction)
0.5
Conductivity (mm/hr)
12.7
Drainage mat
Thickness (mm)
200
Void fraction
0.75
Roughness (Mannings N)
0.014
Vegetative swale
Surface
Berm height (mm)
6
Vegetation volume fraction
0.85
Surface roughness (Mannings N)
0.24
Surface slope (%)
1
Swale side slope (run/rise)
7.33
Rain barrel
Storage
Barrel height (mm)
500With the help of siting tool in System of SUSTAIN, the layout locations of the four LIDs were evaluated and a total of 20 blocks were determined. Then, the LIDs combination schemes in each block were obtained as shown in Figure6. Then, it is necessary to determine the way of the LIDs deployment in each subcatchment of the model. There are two different ways: the first way is to add one or more LID practices, wherein each LID controls the runoff generated by the impervious zone (or permeate zone) at a specific ratio in the subcatchment. This way is suitable for large research area, and different combinations of LIDs can be deployed in subcatchments with multiple land use types [45]. The second way is to decompose the original subcatchments into several small catchment areas with separate LIDs control and underlying surface type. This deployment mode is also called LID series mode, that is, the outflow of one LID can be the inflow of another LID, which is suitable for small area or a situation where a detailed analysis of certain LID facility is required [34]. In this simulation, the second LIDs layout was selected. According to the LIDs combination in each subcatchment, the original subcatchments were decomposed into several small catchment areas with a single LID facilities or land use type. Figure 7 shows a schematic diagram of one LIDs deployment way. Among them, the permeable pavement can handle runoff from the hard surface and direct rainfall. Green roof (rain barrel) handles rainwater runoff from the building roof (rainwater stored in the rain barrel can be used to water green spaces). As a transmission processing unit, the vegetative swale receives the runoff from the original green space, permeable pavement, and green roof before the runoff discharged into the outlet. After formulating the combined deployment plan of the LID practices in each subcatchment, how to design the layout size of the LIDs in each subcatchment to achieve the best cost-effective is the issue to be discussed next.Figure 6
LIDs combination schemes in each block to be transformed.Figure 7
Schematic diagram of LIDs deployment.


## 2.1.1. Case Study Area

Pingxiang City, Jiangxi Province, as one of the pilot cities of the sponge city in China, is known as the “Pearl of Ganxi” and is the center of economic development in Ganxi area. Pingxiang City has a long history and is one of the red revolutionary bases in China. It has gone through multiple stages of development. Especially in recent years, the process of urbanization has accelerated. But at the same time, high urbanization has also brought a series of water problems, such as urban waterlogging and nonpoint source pollution [30]. Since China proposed the construction of the sponge city in 2013, with the country’s attention and resource input, the public’s cooperation and recognition of the “sponge city” has increased gradually, which is more conducive to the construction and development of the sponge city.Pingxiang City is located in a subtropical monsoon humid climate zone. According to the rainfall data of Pingxiang Meteorological Bureau for many years (1955∼2009), the average annual rainfall is as high as 1596.7 mm. Most of the rainfall is concentrated in summer, and it is dry and less rainy when summer and autumn alternate. Recently, due to the impact of global climate change, extreme rainfall events have occurred from time to time [31]. The soil in Pingxiang City is mostly clay soil with poor permeability, which is not conducive to timely infiltration of rainwater. The Ximen waterlogged area studied in this paper is located in the southwest of the old city (As shown in Figure 2). The Pingshui River flows from north to south in the east and south, which accepts rainwater and sewage discharge in the city. The geomorphology of the study area is relatively complex, mostly hilly terrain. The ground elevation ranges from 89.12 to 106.45, showing high in the north and low in the south, and high in the east and low in the west. The total area of the study area is 94 ha, where there are densely distributed building units such as old communities, government units, schools, and commercial districts. The planning of green gardens is inadequate; especially the old communities have almost no green space.Figure 2
Location of Ximen waterlogged area of Pingxiang City.


## 2.1.2. Rainfall-Runoff Model

Before setting up a SWMM model of the study area, it is necessary to collect the underlying surface, the pipe network, rainfall, and other data of the study area and import them into the SWMM. According to the percent of impervious area, the distribution of blocks, and the flow direction of the pipeline, the study area is divided into 59 subcatchments. In order to simplify the model running process, the drainage pipe network is reasonably generalized into 58 conduits and 4 outfalls. The generalization results are shown in Figure3. The setting method of hydrology, hydraulics, water quality parameters, and rain gauge of each module are, respectively, as follows.Figure 3
Generalization of the study area.When setting up a new model in SWMM, the project defaults could be set in advance. There are two preset options in this simulation, which are infiltration model and routing model. For the former, three main models of SWMM can simulate the infiltration process: Horton, Green-Ampt, and Curve number (CN). Among them, the principle of Horton model is to compare rainfall intensity with changing soil permeability. The core idea is to ensure that the actual infiltration rate is the smaller one of the rainfall intensity and the soil infiltration rate during the entire infiltration process [32]. This method requires only three input parameters, which are determined by the soil characteristics of study area and are easy to obtain. Therefore, this paper chose the Horton method which is widely used in the establishment of SWMM model to simulate the infiltration process. Parameters setting is shown in Table 1. For the other preset options, according to the different methods of solving the flow state in the conduits, there are three methods for simulating the flow of the conduits: Steady flow, Kinematic wave, and Dynamic wave. The simulation in this paper chooses the dynamic wave method that can calculate the complex flow state, such as pressure flow and backflow [33].Table 1
Parameters setting of model hydrological and hydraulic.
Parameters
Description
Values
Max. infil. rate
Maximum rate on the horton infiltration curve (mm/hr)
78.1
Min. infil. rate
Minimum rate on the horton infiltration curve (mm/hr)
3.82
Decay constant
Decay constant for the horton infiltration curve (1/hr)
2
N-imperv
Mannings N for impervious area
0.013
N-perv
Mannings N for pervious area
0.15
Dstore-imperv
Depth of depression storage on impervious area (mm)
3.15
Dstore-perv
Depth of depression storage on pervious area (mm)
5.08
%Zero-imperv
Percent of impervious area with no depression storage (%)
25
Conduit roughness
Manning’s roughness coefficient
0.015According to the method of parameters setting, the main hydrological and hydraulic parameters of the model can be divided into two categories.The first category is parameters that can be measured or calculated from actual data. (1) Area: by taking the 1: 200 satellite map of the study area as the background and combining the automatic calculation function of SWMM, the area of each subcatchment area was obtained. (2) Slope: from the collected elevation data, calculate the ratio of the elevation difference to the horizontal distance of the runoff paths to obtain the average surface slope of each subcatchment. (3) Percent of impervious area: according to the land use map (as shown in Figure4) in the study area, the land use types can be divided into buildings, roads, and green spaces. The runoff coefficients of them are specified as 0.75, 0.85, and 0.25, and they are used as weight to calculate the final percentage of impervious area in the subcatchment. (4) Properties of the drainage conduits and juctions. These parameters are determined by the information provided by the Pingxiang Surveying and Mapping Institute. The above parameters are presented in Tables 2–4.Figure 4
Land use map of the study area.Table 2
Subcatchment parameters.
Name
Area (ha)
%Imperv
Width (m)
%Slope
S1
0.4418
55.5
66.468
0.36
S2
0.1934
78.83
43.977
0.36
S3
0.4349
77.7
65.947
0.38
S4
1.6089
56
126.842
0.38
S5
1.61
55.7
126.886
0.38
S6
0.8104
79.5
90.022
1.5
S7
1
44
100
1
S8
1
55.7
100
1
S9
3.215
51.5
179.304
1
S10
2
46
141.421
0.35
S11
1.801
55.5
134.201
0.2
S12
0.81
58.5
90
0.2
S13
1.801
58.5
134.201
0.13
S14
0.5171
83.01
71.91
0.13
S15
1.79
49.5
133.791
0.13
S16
0.7
52
83.666
0.4
S17
2.2251
53.5
205.55
0.38
S18
0.9606
69.99
98.01
0.38
S19
2
55.5
141.421
0.2
S20
2.32
64
152.315
0.4
S21
1.4408
69.99
120.033
0.38
S22
0.68
56
82.462
0.4
S23
0.68
53.5
82.462
0.5
S24
0.2308
76.03
48.042
0.11
S25
1.5758
75.79
125.531
0.11
S26
0.2426
80.95
49.254
0.5
S27
0.2656
76.73
51.536
1.1
S28
2.5989
46
161.211
0.5
S29
2.58
41.5
160.624
0.11
S30
1.5821
53.5
125.782
0.38
S31
1.58
53.5
125.698
0.38
S32
2.3336
53.5
152.761
0.5
S33
0.89
54
94.34
0.5
S34
3.54
64
188.149
0.4
S35
1.1323
60.87
106.41
0.4
S36
0.358
71.5
59.833
5.04
S37
10.8248
75.87
329.01
0.1
S38
0.358
56
59.833
0.4
S39
5
48.5
223.607
0.1
S40
1.5779
53.5
125.614
0.5
S41
1.2224
80.83
110.562
0.11
S42
0.7574
82.52
87.029
1
S43
0.2
64.5
44.721
0.5
S44
0.1047
78.86
32.357
1.18
S45
0.4356
56.5
66
1.18
S46
0.1677
80.06
40.951
0.03
S47
2.335
58.5
115.542
0.5
S48
0.33
54.7
57.446
1
S49
1.2142
79.42
110.191
1.1
S50
1.335
53
115.542
1.1
S51
0.3358
71.49
57.948
0.23
S52
0.305
77.67
55.227
0.5
S53
0.45
54.7
67.082
0.5
S54
7.59
44.7
275.5
0.5
S55
0.4058
77.52
63.702
0.5
S56
1.45
54.7
120.416
0.23
S57
3.65
43.5
191.05
0.5
S58
2.33
44.7
115.326
0.5
S59
1.45
43.5
120.416
0.5Table 3
Juction parameters.
Name
Elevation (m)
J1
103.981
J2
102.89
J3
101.918
J4
98.102
J5
95.115
J6
93.913
J7
103.71
J8
95.115
J9
95.115
J10
98.102
J11
101.6
J12
100.001
J13
102.226
J14
101.6
J15
100.001
J16
98.889
J17
103.71
J18
102.226
J19
101.6
J20
104.505
J21
109.424
J22
108.338
J23
102
J24
105.356
J25
105.803
J26
108.719
J27
107
J28
98.102
J29
104.506
J30
102.226
J31
102.226
J32
101.001
J33
98.102
J34
95.115
J35
106.003
J36
106.003
J37
104.905
J38
103.71
J39
98.898
J40
95.115
J41
92.713
J42
90.443
J43
92.713
J44
94.5
J45
98.102
J46
95.115
J47
95.368
J48
91.908
J49
88.341
J50
87.5
J51
86.287
J52
85.215
J53
84.143
J54
86.287
J55
85.215
J56
88.341
J57
97.87
J58
107.126
PFK3
83.034
PFK2
89.190
PFK4
90.443
PFK1
97.300Table 4
Conduits parameters.
Name
From node
To node
Length (m)
pipe diameter (m)
1
J1
J2
120
0.6
2
J2
J3
120
0.6
3
J3
J4
120
0.6
4
J4
J5
120
1
5
J5
J6
180
1.2
7
J7
J13
120
1
11
J14
J15
150
1.2
12
J15
J16
150
2.4
13
J10
J9
90
1
14
J9
J6
180
1.2
15
J8
J6
130
1.2
16
J28
J8
90
1
17
J17
J18
120
1
18
J18
J19
140
1.2
19
J20
J17
120
1
21
J26
J27
120
2
22
J27
J24
150
1.2
23
J25
J24
150
1.2
24
J24
J23
180
2
25
J35
J37
120
0.8
26
J36
J37
120
1
27
J37
J38
120
1
28
J29
J17
120
1
29
J30
J19
120
1
30
J33
J34
120
1
31
J45
J46
120
1
32
J46
J44
120
1
33
J44
J42
120
1.2
34
J34
J43
150
1
35
J43
J42
150
1.2
36
J23
J39
150
2
38
J40
J41
100
1
39
J47
J48
180
2
40
J48
J49
180
2
41
J49
J50
180
2.4
42
J50
J51
180
2.4
43
J51
J52
180
2.4
44
J52
J53
200
2.4
45
J54
J55
120
2
46
J55
J53
120
2
47
J53
PFK3
180
2.4
48
J42
PFK2
150
1.2
49
J6
PFK4
180
2
50
J39
J57
150
2
51
J57
J47
150
2
52
J41
J42
120
1.2
53
J56
J50
180
2.4
54
J16
PFK1
160
2.4
56
J21
J22
120
0.8
57
J22
J58
120
0.8
58
J58
J35
120
0.8
59
J38
J31
160
1.2
60
J31
J32
180
1.2
C64
J13
J11
120
1
C65
J11
J12
120
1.2
C66
J12
J16
120
1.2
C69
J19
J15
150
1.2
C70
J32
J15
180
1.2The second category is parameters that refer to the SWMM User’s Manual or other References. (1) The values of N for impervious (pervious) area and depth of depression storage on impervious (pervious) area were determined according to different types of underlying surface by referring to the SWMM User’s Manual.(2) Percent of impervious area with no depression storage is generally taken as 25% [34]. (3) Most of the conduits are concrete and the roughness coefficients are set by referring to the user’s manual. The above approximate empirical parameters are also shown in Table 1.When it comes to water quality simulation, suspended substance (SS) is the main physical pollutant and most of the pollutants are positively correlated with SS, so it can represent the regional pollution situation [35]. In this simulation, the water quality of the study area was reflected by the total amount of SS. Referring to the relative articles, the concentration of pollutant in rainwater was set to 10 mg/L [36]. Because the existing pollutant monitoring data is incomplete, in this simulation, the initial pollutant buildup on subcatchment is set to zero. The parameters setting for simulating pollutant buildup and washoff in different types of land use are shown in Table 5. One thing to mention is that this model did not perform parameter verification due to the lack of actual monitoring data of pipe network flow and pollutant concentration changes under real rainfall events, and only an approximate model of the study area was constructed.Table 5
Parameters setting of model pollutant buildup and washoff
Buildup (saturation function)
Washoff (exponential function)
Maximum buildup (kg)
Rate constant
Saturation constant (day)
Coefficient
Index
Cleaning efficiency (%)
Building
140
0
4
0.007
1.8
0
Road
270
0
4
0.008
1.8
70
Green space
60
0
4
0.004
1.2
0In order to study the effect of LID practices on runoff control in different rainfall situations, especially short-duration heavy rainfall, design storms with different return periods were used to replace actual precipitation events [37]. The storm intensity formula in Pingxiang City is.(1)q=26191+0.78lgPt+100.79,where q is the average rainfall intensity (L/s·ha), t is the rainfall duration (min), and P is the rainfall return period (year).The commonly used methods of rainfall pattern decomposition include uniform method and nonuniform method [38]. The uniform method is relatively simple, but its design storm results are relatively small, which does not conform to the actual rainfall law [39]. The nonuniform method can be divided into Chicago rainfall pattern, Huff rainfall pattern, and triangular rainfall pattern according to the different design approaches of the rain peak [40]. This study chose the Chicago Law, which is widely applicable to short-duration rainfall and recommended by China’s “Code for Design of Outdoor Wastewater engineering.” By setting the rain peak 0.4 and the return periods 2-year, 10-year and 50-year, the design storms that last for 2 hours are generated and shown in Figure 5.Figure 5
Design rainfall process lines for different return periods.


## 2.1.3. LIDs Layout Planning

Deploying LIDs control in the SWMM model requires two steps: first, set up several sets of LIDs design regardless of size, and then assign any required design combinations and design sizes to the selected subcatchment. SWMM5.1.13 has 8 kinds of LIDs to choose from. The LIDs selected for this simulation in the study area are vegetative swale (VS), permeable pavement (PP), green roof (GR), and rain barrel (RB). Among them, the vegetative swale is a shallow swale with herbs on the surface, which transmits rainwater runoff and has both infiltration and filtration functions [41]. It can trap some rainwater and particulate pollutants to reduce runoff discharge. Permeable pavements include permeable concrete and pavement brick paving in sidewalk, and colored permeable pavement in parking lots, which can quickly infiltrate direct rainfall or runoff formed by impervious pavements [42]. However, the rainwater that cannot infiltrate to permeable pavement in time overflows into the vegetative swale for infiltration and filtering and finally flows into the rainwater pipe network. As a low-impact development measure for buildings, green roof and rain barrel can be selected based on the age and type of building. Generally speaking, the older buildings are not suitable for setting green roof considering the load bearing capacity of roof [43]. The layout of rain barrels are more flexible and can be set on the side of buildings such as houses, factories, and schools, and rainwater stored in it could be used to irrigate green space or be used as domestic water for residents. The parameter settings of the above four types of LIDs in the model referred to the general LID facilities technical specifications in China [44] and was fine-tuned according to the actual deployment environment. The specific parameters setting of LID practices are shown in Table 6.Table 6
Parameters setting of LID practices.
LIDs type
Construction
Parameter option
Parameter value
Permeable pavement
Surface
Berm height (mm)
3
Surface roughness (Mannings N)
0.012
Surface slope (%)
1
Pavement
Thickness (mm)
120
Void ratio (voids/Solids)
0.2
Permeability (mm/hr)
2450
Storage
Thickness (mm)
230
Void ratio (voids/Solids)
0.75
Seepage rate (mm/hr)
300
Green roof
Surface
Berm height (mm)
6.35
Vegetation volume fraction
0.85
Surface roughness (Mannings N)
0.23
Surface slope (%)
1
Soil
Thickness (mm)
600
Porosity (volume fraction)
0.5
Conductivity (mm/hr)
12.7
Drainage mat
Thickness (mm)
200
Void fraction
0.75
Roughness (Mannings N)
0.014
Vegetative swale
Surface
Berm height (mm)
6
Vegetation volume fraction
0.85
Surface roughness (Mannings N)
0.24
Surface slope (%)
1
Swale side slope (run/rise)
7.33
Rain barrel
Storage
Barrel height (mm)
500With the help of siting tool in System of SUSTAIN, the layout locations of the four LIDs were evaluated and a total of 20 blocks were determined. Then, the LIDs combination schemes in each block were obtained as shown in Figure6. Then, it is necessary to determine the way of the LIDs deployment in each subcatchment of the model. There are two different ways: the first way is to add one or more LID practices, wherein each LID controls the runoff generated by the impervious zone (or permeate zone) at a specific ratio in the subcatchment. This way is suitable for large research area, and different combinations of LIDs can be deployed in subcatchments with multiple land use types [45]. The second way is to decompose the original subcatchments into several small catchment areas with separate LIDs control and underlying surface type. This deployment mode is also called LID series mode, that is, the outflow of one LID can be the inflow of another LID, which is suitable for small area or a situation where a detailed analysis of certain LID facility is required [34]. In this simulation, the second LIDs layout was selected. According to the LIDs combination in each subcatchment, the original subcatchments were decomposed into several small catchment areas with a single LID facilities or land use type. Figure 7 shows a schematic diagram of one LIDs deployment way. Among them, the permeable pavement can handle runoff from the hard surface and direct rainfall. Green roof (rain barrel) handles rainwater runoff from the building roof (rainwater stored in the rain barrel can be used to water green spaces). As a transmission processing unit, the vegetative swale receives the runoff from the original green space, permeable pavement, and green roof before the runoff discharged into the outlet. After formulating the combined deployment plan of the LID practices in each subcatchment, how to design the layout size of the LIDs in each subcatchment to achieve the best cost-effective is the issue to be discussed next.Figure 6
LIDs combination schemes in each block to be transformed.Figure 7
Schematic diagram of LIDs deployment.


## 2.2. Optimization Model

### 2.2.1. PICEA-g Algorithm

Preference-inspired co-evolutionary algorithm (PICEA) was first proposed in 2013 to solve many-objective optimization problem (MOOP). This algorithm introduces several random preference sets and co-evolves with candidate solution sets, so that it can make candidate solutions approach the Pareto frontier while looking for decision-makers’ preference areas [46]. On this basis, Wang et al. [25] propose preference-inspired co-evolutionary algorithm using goal vectors (PICEA-g), which uses the goal vectors as preference sets. When solving MOOP, fitness allocation rules based on the original Pareto dominance relationship are added to enhance the ability of selecting individuals in the algorithm. It is verified by the author that the algorithm performs are better than the widely used NSGA-II [23], ε-MOEA [47], and other multiobjective evolutionary algorithms when solving many-objective (three or more) optimization problems. The core of the PICEA-g algorithm is to specify the fitness allocation rules for the goal vectors and the candidate solutions, so that the Pareto-dominated solutions can be retained as the parent population of the next generation [48]. The fitness calculation formulas of candidate solution Fs and goal vector Fg are [25](2)Fs=0+∑g∈G∪GCs≤g1ng,Fg=11+α,(3)α=1,ng=0,ng−12N−1,others,where G and Gc represent the parent and offspring populations of the goal vector g, respectively, ng is the number that satisfies the goal vector g, and N is the population size. If the candidate solution s does not satisfy any goal vector g, then the fitness value Fs of s is taken as 0.For the sponge city low-impact development and construction-planning field, it can be essentially attributed to the MOOP of complex systems. Combined with the superiority of the PICEA-g algorithm, this paper performed many-objective optimization of LID practices planning in order to find a more cost-effective layout solution efficiently.


### 2.2.2. Decision Variables

The blocks in the study area can be divided into four types: administration (A), business (B), green space or square (G), residence (R), and the LIDs combinations deployed in each of the blocks to be transformed are different. Taking the layout ratios of LIDs in each subcatchment as the optimization decision variables, there are 58 decision variables in 20 blocks, and the specific form of decision variableD is shown in(4)D=Num,Type,where Num represents the serial number of the blocks, Type represents four LID practice types: PP, GR, VS, and RB. According to the requirements for the layout of LIDs in the “Technical Guide for Sponge City Construction (Trial)” and the investigation of research area [49], the maximum LIDs layout ratios on the corresponding underlying surface in different types of blocks are shown in Table 7.Table 7
Maximum LIDs layout ratios for different types of blocks.
Land classification
Maximum LIDs layout ratios
PP (%)
GR (%)
VS (%)
A
60
95
86
B
65
95
86
G
75
—
95
R
45
95
75


### 2.2.3. Objective Functions

For assessing the effect of LID practices on the original urban drainage system, it is generally considered from several aspects: the effect of water volume control, the degree of water quality improvement, and the construction cost. In this study, four objective functions that can evaluate the cost-benefit of sponge city construction are selected: total runoff, which represents the sum of surface runoff for each subcatchment; peak flow, which indicates the maximum discharge flow at the system outfall; the SS pollutant is the total mass of SS discharged from the outlet; the total cost of the LID facilities includes construction costs but excludes other costs such as maintenance. The first three objective function values can be extracted from the report file of the SWMM, and the calculation formula for the total costCsum function of the LIDs can be expressed as(5)Csum=∑i=120∑j=13Sij×Cj+∑i=120Di4×C4,where Sij represents the layout ratio of a certain LIDs in a certain block, Cj represents LIDs (PP, GR or VS) construction cost per unit area, Di4 indicates the number of rain barrels in a certain block, and C4 represents the construction cost of RB per unit (the construction cost of the unit LID facility is shown in Table 8).Table 8
Construction cost of unit LID practices [50, 51].
LIDs type
PP (CNY/m2)
GR (CNY/m2)
VS (CNY/m2)
RB (CNY/each)
Construction cost
193
500
220
350


### 2.2.4. Connection of SWMM and Optimization Model

Because the project file of the SWMM is stored in the form of texts as a configuration file, an input file, an output file, and a report file. We first use Matlab to call the text, update the values of the various decision variables, and related parameters generated by the compilation algorithm and call the dynamic link library files of SWMM. Then, by calling the calculation engine, the calculation result of the model which the user interests is read for the calculation of the objective function. In this way, the interactive connection between the SWMM model and the PICEA-g many-objective optimization model is completed, and the specific calculation process is shown in Figure8.Figure 8
Joint calculation process of SWMM and PICEA-g algorithm.


## 2.2.1. PICEA-g Algorithm

Preference-inspired co-evolutionary algorithm (PICEA) was first proposed in 2013 to solve many-objective optimization problem (MOOP). This algorithm introduces several random preference sets and co-evolves with candidate solution sets, so that it can make candidate solutions approach the Pareto frontier while looking for decision-makers’ preference areas [46]. On this basis, Wang et al. [25] propose preference-inspired co-evolutionary algorithm using goal vectors (PICEA-g), which uses the goal vectors as preference sets. When solving MOOP, fitness allocation rules based on the original Pareto dominance relationship are added to enhance the ability of selecting individuals in the algorithm. It is verified by the author that the algorithm performs are better than the widely used NSGA-II [23], ε-MOEA [47], and other multiobjective evolutionary algorithms when solving many-objective (three or more) optimization problems. The core of the PICEA-g algorithm is to specify the fitness allocation rules for the goal vectors and the candidate solutions, so that the Pareto-dominated solutions can be retained as the parent population of the next generation [48]. The fitness calculation formulas of candidate solution Fs and goal vector Fg are [25](2)Fs=0+∑g∈G∪GCs≤g1ng,Fg=11+α,(3)α=1,ng=0,ng−12N−1,others,where G and Gc represent the parent and offspring populations of the goal vector g, respectively, ng is the number that satisfies the goal vector g, and N is the population size. If the candidate solution s does not satisfy any goal vector g, then the fitness value Fs of s is taken as 0.For the sponge city low-impact development and construction-planning field, it can be essentially attributed to the MOOP of complex systems. Combined with the superiority of the PICEA-g algorithm, this paper performed many-objective optimization of LID practices planning in order to find a more cost-effective layout solution efficiently.


## 2.2.2. Decision Variables

The blocks in the study area can be divided into four types: administration (A), business (B), green space or square (G), residence (R), and the LIDs combinations deployed in each of the blocks to be transformed are different. Taking the layout ratios of LIDs in each subcatchment as the optimization decision variables, there are 58 decision variables in 20 blocks, and the specific form of decision variableD is shown in(4)D=Num,Type,where Num represents the serial number of the blocks, Type represents four LID practice types: PP, GR, VS, and RB. According to the requirements for the layout of LIDs in the “Technical Guide for Sponge City Construction (Trial)” and the investigation of research area [49], the maximum LIDs layout ratios on the corresponding underlying surface in different types of blocks are shown in Table 7.Table 7
Maximum LIDs layout ratios for different types of blocks.
Land classification
Maximum LIDs layout ratios
PP (%)
GR (%)
VS (%)
A
60
95
86
B
65
95
86
G
75
—
95
R
45
95
75


## 2.2.3. Objective Functions

For assessing the effect of LID practices on the original urban drainage system, it is generally considered from several aspects: the effect of water volume control, the degree of water quality improvement, and the construction cost. In this study, four objective functions that can evaluate the cost-benefit of sponge city construction are selected: total runoff, which represents the sum of surface runoff for each subcatchment; peak flow, which indicates the maximum discharge flow at the system outfall; the SS pollutant is the total mass of SS discharged from the outlet; the total cost of the LID facilities includes construction costs but excludes other costs such as maintenance. The first three objective function values can be extracted from the report file of the SWMM, and the calculation formula for the total costCsum function of the LIDs can be expressed as(5)Csum=∑i=120∑j=13Sij×Cj+∑i=120Di4×C4,where Sij represents the layout ratio of a certain LIDs in a certain block, Cj represents LIDs (PP, GR or VS) construction cost per unit area, Di4 indicates the number of rain barrels in a certain block, and C4 represents the construction cost of RB per unit (the construction cost of the unit LID facility is shown in Table 8).Table 8
Construction cost of unit LID practices [50, 51].
LIDs type
PP (CNY/m2)
GR (CNY/m2)
VS (CNY/m2)
RB (CNY/each)
Construction cost
193
500
220
350


## 2.2.4. Connection of SWMM and Optimization Model

Because the project file of the SWMM is stored in the form of texts as a configuration file, an input file, an output file, and a report file. We first use Matlab to call the text, update the values of the various decision variables, and related parameters generated by the compilation algorithm and call the dynamic link library files of SWMM. Then, by calling the calculation engine, the calculation result of the model which the user interests is read for the calculation of the objective function. In this way, the interactive connection between the SWMM model and the PICEA-g many-objective optimization model is completed, and the specific calculation process is shown in Figure8.Figure 8
Joint calculation process of SWMM and PICEA-g algorithm.


## 3. Results and Discussion

### 3.1. Cost-Benefit Optimal Solution

The significance of multiobjective optimization is to trade off multiple objectives, so as to obtain a nondominated solution that makes all objective function values as optimal as possible, which is generally obtained by analyzing the Pareto optimal front. However, when the number of targets is greater than three, it is difficult to intuitively analyse the Pareto optimal solutions through the scatter plot of the objective function values. The existing solutions include parallel co-ordinates, mapping, and heat maps [52, 53]. These visualization methods are slightly complicated, and the dominance relationship of the Pareto solution sets cannot be fully presented. So this paper makes the following analysis with the help of a two-dimensional scatter plot between two targets.According to Zhang et al. [54], the runoff control effect of LID practices is more significant in low return period and short-duration rainfall. Sun et al. [55] prove that the LID practices are difficult to control the rainfall of high-intensity rainfall events. In this paper, we will analyse the Pareto optimal solution sets of LIDs layout ratio in short-duration rainfall with different intensities, especially short-duration heavy rainfall. Each solution corresponds to a set of LIDs layout schemes determined by 58 decision variables. The two-dimensional scatter plots between the pair of targets are plotted in Figure 9, which show the relationship of total runoff-cost, peak flow-cost, SS pollutant-cost under the design rainfall conditions of 2-year, 10-year and 50-year, and the cost-effective optimal solution under the current target has been marked. Based on this, it is convenient for decision-makers to make the most economical choice from different emphasis target according to the law of marginal benefits in economics.Figure 9
Nondominated solution sets under different return periods rainfall of PICEA-g.By analyzing the nondominated solution sets shown in Figure9, it can be found that the peak flow fluctuates slightly with the change in LIDs cost. This just shows that the strong selection ability of PICEA-g algorithm makes Pareto solution set more centralized to the cost-benefit optimal solution. It can also be found that the SS pollutant has an approximate linear relationship with the LIDs cost. This shows that the optimization of this goal has a tiny lifting space, which may relate to the limited LID types. According to Freeborn et al. [56] and Horst et al. [57], the SS pollutant removal rate of permeable pavement and dry vegetative swale can reach more than 90%, but the rain barrel and green roof is extremely poor of SS reduction. Therefore, the more optimal solution of SS pollutant objective may not be found in this paper after optimizing the permeable pavement and vegetative swale area to the most extent by the PICEA-g algorithm.According to the law of marginal benefits, two sets of optimal schemes are shown in Tables9 and 10. Table 9 shows the first scheme: the total runoff-cost optimal solutions, and Table 10 shows the second scheme: the peak flow-cost optimal solutions. Through further analysis of them, it can be known that in the case of 2-year and 10-year return period rainfall, the total runoff-cost optimal solution provided by Scheme 1 has a relatively low engineering budget and a high rate of total runoff reduction. In addition, the peak flow and total SS reduction rates are not much different from those in Scheme 2. Therefore, the total runoff-cost relationship can be used as the optimal cost-benefit solution for the LIDs project in the study area. In the case of a design storm of 50-year return period happened, the reduction rates of the runoff control indicators under both scenarios are reduced compared to those of lower-intensity rainfall, and the LIDs costs are significantly reduced. This shows that it is difficult to cope with short-duration heavy rainfall cost-effectively by increasing the cost of the LIDs project only, which was also found by Zanandrea and Silveira [58]. It can be concluded that when heavy rains are encountered in the study area, the overall performance of the LID practices is absolutely affected, and this impact could not be overcome by increasing the layout ratio of the LIDs.Table 9
Scheme 1: total runoff-cost optimal solution.
Return period (year)
Total runoff (mm)
Runoff reduction rate (%)
Peak flow (CMS)
Peak flow reduction rate (%)
SS pollutant (kg)
SS pollutant reduction rate (%)
LIDs cost (CNY)
2
21.99
22.1
3.838
15.4
3253
35.9
9,307,000
10
38.97
18.7
6.403
11.2
3801
32.8
9,597,000
50
56.72
18.0
8.948
7.5
4065
31.6
8,397,000Table 10
Scheme 2: peak flow-cost optimal solution.
Return period (year)
Total runoff (mm)
Runoff reduction rate (%)
Peak flow (CMS)
Peak flow reduction rate (%)
SS pollutant (kg)
SS pollutant reduction rate (%)
LIDs cost (CNY)
2
22.054
21.8
3.832
15.6
3247
36.0
9,759,000
10
39.641
17.3
6.357
11.8
3765
33.4
10,430,000
50
57.408
17.0
8.899
8.0
4081
31.1
7,573,000


### 3.2. Simulation Results of the Optimal Solution Model

Scheme 1 (under the situation of 10-year return period rainfall) is simulated and verified by using SWMM. The simulation results are shown in Figure10. It can be known from the system runoff curves that the maximum flow of the system after LIDs deployment is reduced to varying degrees compared to before LIDs transformation. Moreover, in the case of 2-year, 10-year, and 50-year return period storms, the total runoff reduction rates of the system are 21.8%, 18.7%, and 16.1%, and the peak flow reduction rates of the system are 18.1%, 15.95% and 10.3%; SS pollutants are reduced by 35.5%, 32.8%, and 30.9%, respectively, compared with before LIDs transformation.Figure 10
Simulation results of the optimal solution. (a) Runoff curves under different rainfall periods before and after transformation; (b) SS pollutant of the system under different rainfall periods before and after transformation.
(a)
(b)Furthermore, analyzing the runoff curves of the system under different rainfall intensities in Figure10, the time when the maximum flow of the system appeared did not achieve a significant delay, even earlier before LIDs transformation. This is because that the optimization process does not consider the goal of delaying the occurrence of peak flow, which can be considered in future optimization studies. In addition, under the short-duration heavy rainfall of 50-year return period, the system has two nodes flooding at J58 and J35 and risk of surcharge of five conduits, which are shown in Figure 11. The risk of storm overflow is predicted to occur at 51 minutes, just after the arrival of this design rainfall peak. This is due to the limited ability of LID practices to resist short-duration heavy rainfall. According to Hua et al. [59] and Li et al. [60], the runoff control effect of LIDs decreases with the increase of rainfall intensity and rainfall duration, while the detention tank can complement the advantages of LIDs. LIDs alone cannot completely avoid the risk of rainwater pipe network overflow, so other sponge facilities that have the function of artificial regulation and rainwater storage (such as rainwater pumping stations, storage reservoirs) need to be constructed coordinately.Figure 11
Profile plot of rainwater pipe with overflow risk.


### 3.3. Compared with the Application of NSGA-II Algorithm

The algorithm applied in this study is compared with the result of applying NSGA-II algorithm to solve this problem. After the program is run iteratively for 50 times, the optimal solution sets are obtained and shown in Figures12 and 13 (it should be noted that only the optimization results under the 2-year return period are shown). The range of the former solutions is more concentrated near the cost-benefit optimal solution. It reflects the advantage of PICEA-g to enhance the selection capability of the nondominated solutions and improve the convergence speed of the algorithm. Secondly, it can be seen that the solutions obtained by the PICEA-g method are widely distributed at the Pareto optimization front, the distribution of solution sets is better than the NSGA-II method. From the above analysis, the PICEA-g method applied in this paper is more convenient for decision-makers to compare and choose the relatively optimal solution.Figure 12
Pareto solution sets of PICEA-g algorithm.Figure 13
Pareto solution sets of NSGA-II algorithm.Similarly, the NSGA-II method is analysed according to above analysis of the nondominated solution sets. By analyzing the scatter plots between the targets under different return period storms (as shown in Figure14), the optimal schemes with different emphasis on the targets are obtained. After comparing these schemes horizontally, the corresponding optimal solution of selected LIDs layout scheme is obtained and shown in Table 11. Among them, when analyzing a two-dimensional scatter plot between two targets, because the Pareto solution sets obtained by this method is relatively scattered and has poor convergence, it is difficult for decision-makers to analyse the most cost-effective solution from the Pareto optimal front.Figure 14
Nondominated solution sets under different return periods rainfall of NSGA-II.Table 11
Optimal solution of LIDs layout scheme under NSGA-II method.
Total runoff (mm)
Peak runoff (CMS)
SS pollutant (kg)
Total cost (CNY)
59.579
9.134
4109
9,070,000Then, the LIDs layout scheme corresponding to the optimal solution is simulated and verified by SWMM and is compared with the simulation results of the PICEA-g method (the results are shown in Table12). It can be known that the total costs of the optimal scheme of PICEA-g and NSGA-II method are 9.597 million and 9.070 million, respectively. The total reduction rates of system runoff in the NSGA-II method reach a maximum of 21%. The peak flow reduction rates of the system are between 5% and 15%, and the SS pollutants are reduced by about 30% compared with before LIDs transformation. The PICEA-g method increases the reduction rate of each runoff index by average of 1.7, 11.1, and 0.3 percentage points, respectively. Among them, the optimization of SS is the smallest, which verifies the previous discussion that the SS pollutants target has a small optimization space under the current decision variables due to the limitation of the LIDs type. Therefore, in the face of the recent implemented sponge city construction assessment standard with multiple indicators, when the PICEA-g algorithm is applied to the multiobjective (especially four or more targets) optimization processing, its performance in all aspects is better than that of the traditional optimization algorithm. This will help designers to search optimal solutions more effectively which meet various construction indicators.Table 12
Comparison of simulation results of the optimal schemes under NSGA-II and PICEA-g.
Return period (year)
Total runoff reduction rate (%)
Peak flow reduction rate (%)
SS pollutant reduction rate (%)
PICEA-g
NSGA-II
PICEA-g
NSGA-II
PICEA-g
NSGA-II
2
21.8
21.0
18.1
14.5
35.5
35.2
10
18.7
16.4
16.0
10.1
32.8
32.5
50
16.1
13.9
10.3
5.6
30.9
30.7


## 3.1. Cost-Benefit Optimal Solution

The significance of multiobjective optimization is to trade off multiple objectives, so as to obtain a nondominated solution that makes all objective function values as optimal as possible, which is generally obtained by analyzing the Pareto optimal front. However, when the number of targets is greater than three, it is difficult to intuitively analyse the Pareto optimal solutions through the scatter plot of the objective function values. The existing solutions include parallel co-ordinates, mapping, and heat maps [52, 53]. These visualization methods are slightly complicated, and the dominance relationship of the Pareto solution sets cannot be fully presented. So this paper makes the following analysis with the help of a two-dimensional scatter plot between two targets.According to Zhang et al. [54], the runoff control effect of LID practices is more significant in low return period and short-duration rainfall. Sun et al. [55] prove that the LID practices are difficult to control the rainfall of high-intensity rainfall events. In this paper, we will analyse the Pareto optimal solution sets of LIDs layout ratio in short-duration rainfall with different intensities, especially short-duration heavy rainfall. Each solution corresponds to a set of LIDs layout schemes determined by 58 decision variables. The two-dimensional scatter plots between the pair of targets are plotted in Figure 9, which show the relationship of total runoff-cost, peak flow-cost, SS pollutant-cost under the design rainfall conditions of 2-year, 10-year and 50-year, and the cost-effective optimal solution under the current target has been marked. Based on this, it is convenient for decision-makers to make the most economical choice from different emphasis target according to the law of marginal benefits in economics.Figure 9
Nondominated solution sets under different return periods rainfall of PICEA-g.By analyzing the nondominated solution sets shown in Figure9, it can be found that the peak flow fluctuates slightly with the change in LIDs cost. This just shows that the strong selection ability of PICEA-g algorithm makes Pareto solution set more centralized to the cost-benefit optimal solution. It can also be found that the SS pollutant has an approximate linear relationship with the LIDs cost. This shows that the optimization of this goal has a tiny lifting space, which may relate to the limited LID types. According to Freeborn et al. [56] and Horst et al. [57], the SS pollutant removal rate of permeable pavement and dry vegetative swale can reach more than 90%, but the rain barrel and green roof is extremely poor of SS reduction. Therefore, the more optimal solution of SS pollutant objective may not be found in this paper after optimizing the permeable pavement and vegetative swale area to the most extent by the PICEA-g algorithm.According to the law of marginal benefits, two sets of optimal schemes are shown in Tables9 and 10. Table 9 shows the first scheme: the total runoff-cost optimal solutions, and Table 10 shows the second scheme: the peak flow-cost optimal solutions. Through further analysis of them, it can be known that in the case of 2-year and 10-year return period rainfall, the total runoff-cost optimal solution provided by Scheme 1 has a relatively low engineering budget and a high rate of total runoff reduction. In addition, the peak flow and total SS reduction rates are not much different from those in Scheme 2. Therefore, the total runoff-cost relationship can be used as the optimal cost-benefit solution for the LIDs project in the study area. In the case of a design storm of 50-year return period happened, the reduction rates of the runoff control indicators under both scenarios are reduced compared to those of lower-intensity rainfall, and the LIDs costs are significantly reduced. This shows that it is difficult to cope with short-duration heavy rainfall cost-effectively by increasing the cost of the LIDs project only, which was also found by Zanandrea and Silveira [58]. It can be concluded that when heavy rains are encountered in the study area, the overall performance of the LID practices is absolutely affected, and this impact could not be overcome by increasing the layout ratio of the LIDs.Table 9
Scheme 1: total runoff-cost optimal solution.
Return period (year)
Total runoff (mm)
Runoff reduction rate (%)
Peak flow (CMS)
Peak flow reduction rate (%)
SS pollutant (kg)
SS pollutant reduction rate (%)
LIDs cost (CNY)
2
21.99
22.1
3.838
15.4
3253
35.9
9,307,000
10
38.97
18.7
6.403
11.2
3801
32.8
9,597,000
50
56.72
18.0
8.948
7.5
4065
31.6
8,397,000Table 10
Scheme 2: peak flow-cost optimal solution.
Return period (year)
Total runoff (mm)
Runoff reduction rate (%)
Peak flow (CMS)
Peak flow reduction rate (%)
SS pollutant (kg)
SS pollutant reduction rate (%)
LIDs cost (CNY)
2
22.054
21.8
3.832
15.6
3247
36.0
9,759,000
10
39.641
17.3
6.357
11.8
3765
33.4
10,430,000
50
57.408
17.0
8.899
8.0
4081
31.1
7,573,000


## 3.2. Simulation Results of the Optimal Solution Model

Scheme 1 (under the situation of 10-year return period rainfall) is simulated and verified by using SWMM. The simulation results are shown in Figure10. It can be known from the system runoff curves that the maximum flow of the system after LIDs deployment is reduced to varying degrees compared to before LIDs transformation. Moreover, in the case of 2-year, 10-year, and 50-year return period storms, the total runoff reduction rates of the system are 21.8%, 18.7%, and 16.1%, and the peak flow reduction rates of the system are 18.1%, 15.95% and 10.3%; SS pollutants are reduced by 35.5%, 32.8%, and 30.9%, respectively, compared with before LIDs transformation.Figure 10
Simulation results of the optimal solution. (a) Runoff curves under different rainfall periods before and after transformation; (b) SS pollutant of the system under different rainfall periods before and after transformation.
(a)
(b)Furthermore, analyzing the runoff curves of the system under different rainfall intensities in Figure10, the time when the maximum flow of the system appeared did not achieve a significant delay, even earlier before LIDs transformation. This is because that the optimization process does not consider the goal of delaying the occurrence of peak flow, which can be considered in future optimization studies. In addition, under the short-duration heavy rainfall of 50-year return period, the system has two nodes flooding at J58 and J35 and risk of surcharge of five conduits, which are shown in Figure 11. The risk of storm overflow is predicted to occur at 51 minutes, just after the arrival of this design rainfall peak. This is due to the limited ability of LID practices to resist short-duration heavy rainfall. According to Hua et al. [59] and Li et al. [60], the runoff control effect of LIDs decreases with the increase of rainfall intensity and rainfall duration, while the detention tank can complement the advantages of LIDs. LIDs alone cannot completely avoid the risk of rainwater pipe network overflow, so other sponge facilities that have the function of artificial regulation and rainwater storage (such as rainwater pumping stations, storage reservoirs) need to be constructed coordinately.Figure 11
Profile plot of rainwater pipe with overflow risk.


## 3.3. Compared with the Application of NSGA-II Algorithm

The algorithm applied in this study is compared with the result of applying NSGA-II algorithm to solve this problem. After the program is run iteratively for 50 times, the optimal solution sets are obtained and shown in Figures12 and 13 (it should be noted that only the optimization results under the 2-year return period are shown). The range of the former solutions is more concentrated near the cost-benefit optimal solution. It reflects the advantage of PICEA-g to enhance the selection capability of the nondominated solutions and improve the convergence speed of the algorithm. Secondly, it can be seen that the solutions obtained by the PICEA-g method are widely distributed at the Pareto optimization front, the distribution of solution sets is better than the NSGA-II method. From the above analysis, the PICEA-g method applied in this paper is more convenient for decision-makers to compare and choose the relatively optimal solution.Figure 12
Pareto solution sets of PICEA-g algorithm.Figure 13
Pareto solution sets of NSGA-II algorithm.Similarly, the NSGA-II method is analysed according to above analysis of the nondominated solution sets. By analyzing the scatter plots between the targets under different return period storms (as shown in Figure14), the optimal schemes with different emphasis on the targets are obtained. After comparing these schemes horizontally, the corresponding optimal solution of selected LIDs layout scheme is obtained and shown in Table 11. Among them, when analyzing a two-dimensional scatter plot between two targets, because the Pareto solution sets obtained by this method is relatively scattered and has poor convergence, it is difficult for decision-makers to analyse the most cost-effective solution from the Pareto optimal front.Figure 14
Nondominated solution sets under different return periods rainfall of NSGA-II.Table 11
Optimal solution of LIDs layout scheme under NSGA-II method.
Total runoff (mm)
Peak runoff (CMS)
SS pollutant (kg)
Total cost (CNY)
59.579
9.134
4109
9,070,000Then, the LIDs layout scheme corresponding to the optimal solution is simulated and verified by SWMM and is compared with the simulation results of the PICEA-g method (the results are shown in Table12). It can be known that the total costs of the optimal scheme of PICEA-g and NSGA-II method are 9.597 million and 9.070 million, respectively. The total reduction rates of system runoff in the NSGA-II method reach a maximum of 21%. The peak flow reduction rates of the system are between 5% and 15%, and the SS pollutants are reduced by about 30% compared with before LIDs transformation. The PICEA-g method increases the reduction rate of each runoff index by average of 1.7, 11.1, and 0.3 percentage points, respectively. Among them, the optimization of SS is the smallest, which verifies the previous discussion that the SS pollutants target has a small optimization space under the current decision variables due to the limitation of the LIDs type. Therefore, in the face of the recent implemented sponge city construction assessment standard with multiple indicators, when the PICEA-g algorithm is applied to the multiobjective (especially four or more targets) optimization processing, its performance in all aspects is better than that of the traditional optimization algorithm. This will help designers to search optimal solutions more effectively which meet various construction indicators.Table 12
Comparison of simulation results of the optimal schemes under NSGA-II and PICEA-g.
Return period (year)
Total runoff reduction rate (%)
Peak flow reduction rate (%)
SS pollutant reduction rate (%)
PICEA-g
NSGA-II
PICEA-g
NSGA-II
PICEA-g
NSGA-II
2
21.8
21.0
18.1
14.5
35.5
35.2
10
18.7
16.4
16.0
10.1
32.8
32.5
50
16.1
13.9
10.3
5.6
30.9
30.7


## 4. Conclusions

The construction of the sponge city needs to take into account scientific design planning, strict project construction as well as perfect operation, maintenance, etc. Among them, the design and optimization of the sponge facilities is the first and critical step, especially for the old urban area. In order to explore the optimal planning and deployment of LIDs, a new LIDs layout mathematical optimization method was innovatively proposed and applied to the Ximen waterlogging area of Pingxiang successfully. Matlab interface program was designed and run, so as to call the dynamic link library files and calculation engine of SWMM. In this way, the block-scaled rainfall-runoff model was combined with the PICEA-g algorithm. Then, the Pareto solution sets of the established optimization model and the corresponding decision variable values were solved under the design storms of different return periods. Based on the consideration of different targets emphasis, several cost-benefit schemes were obtained, and then the optimal plan suitable for the region was selected by analyzing these schemes horizontally. Finally, the area ratio optimization of various LID practices was achieved under the premise of weighing the multiple sponge cities targets. The specific conclusions are shown below:(1)
It can be concluded that under different return period storms, the total runoff reduction rates of the system reach a maximum of 21.8%, the peak flow reduction rates of the system are more than 10%, and the SS pollutants are reduced by about 30% compared with before LIDs transformation. Compared with NSGA-II method, it optimizes the average of 1.7, 11.1, and 0.3 percentage points, respectively. Among them, the SS pollutants objective can be optimized with very little space, which is limited by the type of LIDs. In general, the scheme has a better effect on controlling short-duration, lower-intensity rainfall runoff. However, once faced with the rainfall of 50-year return period, it is difficult to achieve the purpose of the optimal cost-benefit of the sponge city construction only by increasing the LIDs layout area. Consideration should also be given to the integrated construction of artificial control strategies, such as storage units and pumping stations.(2)
Additionally, thanks to the superior performance of the PICEA-g algorithm for solving MOOP, decision-makers can analyse the cost-benefit optimal solution more effectively from the Pareto solution sets compared with the NSGA-II method. Furthermore, with the in-depth practice of the assessment standard of sponge city construction, more and more sponge city effect indicators have been incorporated into designing requirements. For the future research and design, based on the mathematical optimization model and method proposed in this paper, researchers can increase the targets that need to be considered, such as the target of delaying the occurrence of peak flow. In this case, the new many-objective optimization model could be further analysed and discussed.(3)
In addition, how to combine LID practices with other sponge facilities systematically to cope with short-duration heavy rainfall is also the key to the construction effect of sponge cities and public acceptance. Other researchers could also analyse and simulate this situation based on the mathematical optimization model proposed in this paper. And it can be further applied to the practice of sponge city construction.



---

*Source: 6734081-2020-04-28.xml*