A method for optimizing urban LID facility layout based on SOBOL sensitivity analysis

Abstract

The present invention discloses a method for optimizing the layout of urban LID facilities based on SOBOL sensitivity analysis. The method is characterized by first establishing a SWMM model of the urban area under study, constructing a functional relationship between node overflow and LID scheme construction costs within the SWMM model, and then using the SOBOL method to conduct a global sensitivity analysis of the LID facility layout scale parameters for each sub-catchment area. Based on the results of the parameter sensitivity analysis, the LID facility layout area for the selected sensitive variables is adjusted and optimized. The Pareto optimal solution sets before and after parameter optimization are compared and analyzed to obtain a more optimal urban LID facility layout scheme. Compared to existing technologies, this method improves the reduction effect of node overflow by increasing the LID layout area in sensitive areas and reduces the LID layout area in insensitive areas to reduce scheme construction costs. Consequently, the present invention has higher computational efficiency and makes the LID spatial layout more rational.

Description

Technical Field

[0001] The present invention belongs to the technical field of urban stormwater simulation, and in particular relates to an urban LID facility layout optimization method based on SOBOL sensitivity analysis. Background Art

[0002] Urbanization has led to a continuous expansion of impervious surfaces, resulting in a dramatic increase in surface runoff. To reduce runoff during heavy rainstorms, alleviate pressure on urban drainage systems, mitigate urban waterlogging, and control runoff pollution, my country has proposed the concept of sponge city development. These efforts aim to mimic natural processes through low-impact development (LID) approaches, minimizing the negative impacts of urban development on the hydrological environment. In sponge city development, LID facilities such as rain gardens (RGs), bioretention basins (BRCs), and permeable pavement (PP) are commonly deployed in various subcatchments. However, their spatial layout depends on factors such as the number of subcatchments, the type of LID facilities, their size, area, and location, which in turn influence the effectiveness of urban waterlogging mitigation and the cost of sponge city construction. Given the countless spatial layout combinations available for LID facility type, size, area, and location, finding the optimal spatial layout to achieve optimal sponge city results at the lowest construction cost is a multi-objective optimization problem.

[0003] Currently, solving multi-objective optimization problems primarily relies on multi-objective evolutionary algorithms (MOEAs), which seek Pareto-optimal solutions by solving complex, nonlinear, multi-dimensional optimization problems. However, MOEAs often require tens of thousands of evaluations to solve multi-objective optimization problems. This is particularly true for problems involving a large number of decision variables. The high number of parameters and complex models often lead to a sharp decline in algorithm performance and excessive computational time.

[0004] The applicant previously applied for a patent titled "An Interactive Method Based on the SWMM Model and the platEMO Platform." This method includes the following steps: first, establishing a SWMM model for the urban area under study, dividing the subcatchments and exporting an inp file; then, deploying LID facilities in each subcatchment, using the LID facilities in each subcatchment as independent variables; establishing a multi-objective optimization model based on MATLAB software, using the LID facility construction cost and peak flow rate as two objective functions; encapsulating the construction cost and peak flow rate functions using function commands; and then, using the platEMO platform, setting constraints and, following the platform's custom problem-solving steps, utilizing the platform's existing functions to minimize the construction cost and peak flow rate of the two problem functions and complete the processing of each subcatchment one by one. This invention enables interaction between the SWMM model and the platEMO platform, facilitating the use of the platform's various optimization algorithms to refine and optimize the spatial layout of urban LIDs. However, this method still suffers from the drawbacks of excessive algorithm evaluations, resulting in reduced algorithm efficiency, decreased quality of the Pareto optimal solution, and decreased computational analysis accuracy. Summary of the Invention

[0005] In view of the above-mentioned deficiencies in the prior art, the technical problem to be solved by the present invention is: how to provide an urban LID facility layout optimization method based on SOBOL sensitivity analysis that can improve the algorithm efficiency and performance, improve the quality of Pareto optimal solutions, and improve the reliability of its layout optimization analysis.

[0006] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0007] An urban LID facility layout optimization method based on SOBOL sensitivity analysis is characterized by first establishing a SWMM model of the urban area under study, constructing a functional relationship between node overflow and LID scheme construction cost in the SWMM model, and using the SOBOL method to conduct a global sensitivity analysis on the layout scale parameters of the LID facilities in each sub-catchment area to examine the significance of the parameters' impact on node overflow; then using the NSGA-II algorithm in the MATLAB software platEMO4.0 platform to find the Pareto optimal solution set of the LID facility spatial layout, and adjusting the LID facility layout area of ​​the sensitive variables selected and optimized according to the parameter sensitivity analysis results, and comparing and analyzing the Pareto optimal solution sets before and after parameter optimization to obtain a better urban LID facility layout scheme.

[0008] Thus, building on the existing interactive method based on the SWMM model and the platEMO platform described in the background, this method conducts a global sensitivity analysis of LID facility layout scale parameters. LID facility layout areas are adjusted based on the impact of these parameters on node overflows. After comparative analysis, a more optimal urban LID facility layout solution is obtained. This improves the quality of the optimal LID facility spatial layout solution under the control objectives and enhances the efficiency of solving multi-objective optimization problems. Node overflow refers to the amount of water overflowing from a node in the urban drainage system (such as a manhole) during rainfall when the node exceeds its design flow rate.

[0009] Furthermore, the method comprises the following steps:

[0010] Step a: Obtain basic data related to modeling of the study city area and build a SWMM model of the city in SWMM software. Using SWMM software, divide the study area into x subcatchments based on the city's pipe network, building, and street distribution. Complete the SWMM model parameter calibration. Finally, successfully build a complete SWMM model of the city and export the inp file.

[0011] Step b: Deploy y types of LID facilities in each subcatchment (furthermore, the deployed LID facilities include but are not limited to three types of LID facilities: permeable pavement (PP), bioretention basins (BRC), and rain gardens (RG)). The deployment area of ​​y types of LID facilities in x subcatchments is used as the model independent variable, for a total of z independent variables, z = x·y. Using MATLAB software, a multi-objective optimization model is established with LID facility construction cost and node overflow volume as the two objective functions. A functional relationship of LID facility construction cost is established in the multi-objective optimization model, namely:

[0012] (1)

[0013] Where: f1 is the construction cost function of LID facilities in the study area, in yuan; i is the number of sub-catchments, i = 1, 2, 3, ..., x; j is the type of LID facilities, j = 1, 2, ..., y; S ij is the area of ​​the jth type of LID facility in the ith subcatchment area; M j is the unit area construction cost of the jth type of LID facility (the unit area construction cost of each type of LID facility is obtained based on the relevant engineering data of the urban area construction in the study);

[0014] Use matlab's function command to encapsulate the construction cost function;

[0015] The specific process of encapsulating the construction cost function during implementation is to use the function function (declare the function name, input and output) to define the cost function paper1_problemCB2(mian_ji); determine zhsqy_number (the number of sub-catchment areas), whose value is x (27), which is one y-fold of the total z (81) model independent variables, that is, length(mian_ji) / y; use the reshape function (reconstruct the array) to re-sort the area matrix according to y rows and x columns, and define it as CB_before; finally, calculate, construction cost = CB_before×[unit area construction cost of various LID facilities], that is, the area matrix of z model independent variables multiplied by the array matrix composed of the unit area construction cost of each type of LID facility, and use the sum function (sum function) to sum it up, that is, complete the construction of the script function relationship of the construction cost;

[0016] Step c: Establish a functional relationship for node overflow;

[0017] (2)

[0018] Where: f2 is the node overflow, unit is m 3 / s; g is the SWMM model calculation and result extraction function, which is obtained by Matlab programming;

[0019] Use matlab's function command to encapsulate the node overflow function;

[0020] The specific process of encapsulating the node overflow function during implementation is as follows:

[0021] 1) Start reading and rewriting data from the inp file;

[0022] First, define the peak flow function paper1_problemYL2(mian_ji) using the function function (declare the function name, input and output), and set it in a folder. mian_ji corresponds to S_ij in equations (1) and (2). Open the inp file of the SWMM model in a readable and writable way, use the fgetl function (read data line by line) to read the original file data line by line and create a new object newline to receive each line of data in the original file.

[0023] 2) Determine the number of LID facilities (LID-number) and the number of lines occupied by LID in the SWMM model inp file;

[0024] Use the contains function (query parameter position) to read the lia1 and lia2 parameters of the line where the strings "LID_USAGE" and "JUNCTIONS" are located. In the SWMM model, the difference between lia1 and lia2 minus 4 is the number of LID facilities; this number corresponds to j in equations (1) and (2), which is fixed here as y;

[0025] 3) Determine zhsqy_number and the number of subcatchments in the inp file of the SWMM model;

[0026] Use the contains function (query parameter location) to read the lia3 and lia4 parameters in the row where the strings "SUBCATCHMENTS" and "SUBAREAS" are located. In the SWMM model, the difference between lia3 and lia4 minus 4 is the number of subcatchments.

[0027] 4) Get the area and name of the subcatchment;

[0028] Use the cell function (create cell array) to create an empty cell array with n=0, and use the zeros function (create zero matrix) to initialize the area matrix; then establish a for loop with the parameter isi, extract the area of ​​the sub-catchment area from 1 to the number of zhsqy_number, and determine the area position of the SWMM model inp file (such as {lia33+3-1+isi} row (51:59) column) (same as above); use the strip function (delete the whitespace at the beginning and end of the string) to remove the spaces in the file and keep only the text; finally, use the str2double function (digital text) to convert the text into a single string. Convert the subcatchment names to numeric values ​​(subcatchment names may not be continuous during implementation and need to be extracted for subsequent operations). When extracting names (similar to the area extraction step above), use the cell function to create an empty cell array of n=0 for the area, and use the zeros function to initialize the subcatchment name matrix. Then, establish a for loop with the parameter isi2 to extract the subcatchment names from 1 to zhsqy_number, and determine the name location in the SWMM model inp file (e.g., {lia33+3-1+isi2} rows (1:3) columns) (same as above). Finally, use the strip function to remove spaces in the file, retaining only the text.

[0029] 5) Modify the inp file of the SWMM model and splice it into a new variable; each subsequent calculation will read the corresponding data on the original inp file and generate a new inp file; among them, newline(1:lia1+2) is the beginning of the new file newline1, and newline(lia22-2+1:length(newline)) is the end; the modified part in the middle, cell(1,(3*zhsqy_number)), is a cell array of the number of LIDs created by the cell function. These three parts are spliced ​​into a new newline1 variable;

[0030] 6) Add y default LID facilities (such as PP, BRC and RG) for all subcatchments;

[0031] Specifically, we create a for loop with the parameter isi2, from 1 to the number of zhsqy_number; we use the cell2mat function (which converts the cell array into a normal array of the basic data type) and use the area matrix (mjjz) as a variable to convert it into a normal array of the basic data type, with the area initialized to 0;

[0032] 7) Calculate the sum of the LID areas in each subcatchment and set constraints;

[0033] Use the reshape function (reconstruct array) to convert the area matrix (mian_ji) into a readable format of y rows and x columns in column order, and use the sum function (sum function) to calculate the sum of the areas of the LID facilities in each sub-catchment area, represented by CB_before; establish a for loop with the parameter isi3, from 1 to zhsqy_number, the number of sub-catchments, and use the strrep function (replace the specified substring in the string) to replace the area parameter in newline1 with the set value (for example, {isi2*3+lia11+3-2} is replaced by 2); use the if statement to set the constraint condition in the for loop, CB_before should be less than the area matrix (mjjz), that is, the sum of the areas of each sub-catchment area. If it exceeds the sum of the sub-catchment areas, it is set to 0 and the LID is not set; (This constraint is similar to the platEMO constraint condition later, but if this condition is not set, even if the platEMO constraint condition is set correctly, due to the instability of its constraints, some values ​​that exceed the constraints will occur when the algorithm is randomly initialized in subsequent operations, which will cause SWMM to fail to run)

[0034] 8) Run SWMM to perform calculations and extract the OUT file of the SWMM model. (In implementation, the SWMM interface can be called according to the patent developed by the research group.) Use the fopen function (read and write file data) to open the OUT file after the calculation in a readable and writable manner. Then use the fgetl function to read the file data line by line and create a new object newline_out to receive each line of data from the original file.

[0035] 9) Determine the overflow location of the sub-catchment node in the OUT file report and calculate the cumulative node overflow volume; use the contains function to read the lia_out and lia2_out parameters where the strings "Node Flooding Summary" and "Storage Volume Summary" are located. The number of rows to be extracted is from lia_out+10 to lia2_out-4 in newline_out; establish a for loop with parameter i23, extract the above-mentioned position rows line by line from 1 to all the rows to be extracted, and determine that the node overflow volume value is the seventh group of data (this is certain because this column of the SWMM OUT file is the result value of the node overflow volume); finally, use the str2double function to convert all extracted strings into numbers, and use the sum function to sum them, thus completing the script function relationship construction of the node overflow volume.

[0036] Step d: The above two steps develop corresponding programs that encapsulate the two functions of construction cost and node overflow, namely paper1_problemYL2(x) and paper1_problemCB2(x). This completes the preparation work of the correlation between construction cost and node overflow, making it meet the data interface format of the platEMO platform, facilitating the subsequent use of the platform's NSGA-Ⅱ evolutionary algorithm to solve the solution cost and node overflow.

[0037] Based on the platEMO platform, (in order to successfully call the platform's evolutionary algorithm, it is necessary to set the corresponding constraints, that is, the LID facility area cannot be larger than the total area of ​​the subcatchment. If the sum of the set LID areas exceeds the subcatchment area, SWMM will not be able to successfully simulate, causing the connection between SWMM and platEMO to crash. This step will constrain the variables to meet the input conditions of the SWMM calculation engine; therefore, it is necessary to set the constraints as follows:

[0038] (3)

[0039] Where: is the area of ​​the ith subcatchment;

[0040] According to the steps for solving the customized problem of the platEMO platform, the parameter encoding (the encoding method of each variable) is first used to define the y (3) types of LID facilities on the x (27) sub-catchment areas, that is, the encoding method of a total of z (81) variables, and the ones function (creating a matrix with 1 as the value) is used to create a matrix with 1 row and z columns consisting of 1 values; then the parameter lower is used to define (the lower bound of each variable) the lower bound of the z independent variables, and the zeros function (creating a matrix with 0 as the value) is used to create a matrix with 1 row and z columns consisting of 0 values, that is, the minimum layout area of ​​each LID facility is 0; the parameter upper (the upper bound of each variable) is used to define the upper bound of each independent variable (such as the sub-catchment area S1, whose area is 3.1979hm 2 , then the upper limit of the area of ​​each LID facility in the S1 area is 3.1979hm 2 , and the lengths of the lower and upper bound matrices created must be the same as the 1-valued matrices defined by the encoding parameter. Using the parameter objFcn (which minimizes the objective function), minimize the two problem functions, construction cost and node overflow. Finally, using the parameter conFcn (constraint function), the constraint is satisfied if and only if the constraint violation value is less than or equal to zero (for example, in area S1, the sum of the areas of the y types of LID facilities minus the total area of ​​S1, 3.1979, should be less than or equal to zero). Repeat the above steps to complete the processing of other subcatchments. This allows the subsequent simulation analysis of the functional relationship between construction cost and node overflow to be successfully invoked using the evolutionary algorithm of the platEMO platform.

[0041] Step e: Use SOBOL sensitivity analysis method to perform sensitivity analysis on the z independent variable parameters in step a. The analysis process is obtained by Matlab programming (this step can be performed simultaneously with steps b, c, and d). The details are as follows:

[0042] Since there are z parameters, the parameter dimension D is set to z. To ensure symmetrical sampling, the parameter dimension M is expanded to D*2. The number of sampling points N is set to 5000, which means 5000 sampling groups. The upper and lower limits of each parameter are set. The zeros function is used to set a 1-row, 81-column zero matrix as the lower limit of each parameter, named MIN. The upper limit of each parameter cannot exceed the total area of ​​the subcatchment, named MAX. Since each catchment has three LIDs, the MAX array is divided by 3 to ensure uniform sampling.

[0043] 2. Use the sobolset function to generate uniformly distributed sample points according to the parameter dimension M, named p; establish a for loop function to screen the sample points from 1 to N, and obtain the sample of p one by one, named r, with the range r=Min + r.*(Max-Min); finally, distribute the sampling results to the matrix R, that is, the N×M parameter matrix.

[0044] 3 Split the matrix R into matrix A and matrix B; the first D columns of matrix R are divided into matrix A, that is, the format is N×D, and the columns from D+1 to M are divided into matrix B, that is, the format is N×D. Construct matrix AB based on matrices A and B i (i=1, 2, ..., D). The specific process is to establish a for loop function to replace the i-th column in matrix B with the i-th column in matrix A from 1 to D, while the other columns remain unchanged. This can be replaced D times, for a total of D AB matrices.

[0045] 4 The above steps give A, B, AB i (i=1, 2, ..., D) These matrices, totaling (D+2)*N groups of input data, can calculate (D+2)*N groups of Y values; input (D+2)*N groups of input data into the node overflow function constructed in step a to solve for Y A 、Y B 、Y ABi (i=1, 2,…, D) values.

[0046] 5 Calculate the main effect index S i , the formula is:

[0047] (1)

[0048] Where: X i is the i-th independent variable; Y is the target variable; V represents the average value of the output target variable Y under the condition of inputting the i-th independent variable X; xi For the variance operator, input the independent variable X when other variables are fixed. i Variance calculation is performed; V(Y) is the total variance of the output target variable Y, which represents the variance of all independent variables X1, X2, ..., X k The variance of the output target variable Y under the joint action; S i It is the first-order sensitivity index, that is, the direct contribution of the ith independent variable to the output target variable Y, without considering the interaction between the independent variables.

[0049] The specific calculation process is obtained through Matlab programming, as follows:

[0050] Define V using the zeros function X With S i, its format is (D, 1), establish a for loop function, from 1 to N, calculate V for each row X(i) Value, V X(i) = V X +Y B *(Y AB -Y A );V X(i) / N is the V of each parameter variable X Value, V Y The format is defined as ([Y A ;Y B ],1),[Y A ;Y B ] is Y A and Y B The matrix is ​​stacked; finally, S is calculated i =V X / V Y .

[0051] 6Calculate the total effect index S Ti , the formula is:

[0052] (2)

[0053] Where: X~i is the remaining variable after removing the ith independent variable X; E[Y|X_(~i)] represents the average value of the output target variable Y under the condition of the remaining variable after removing the ith independent variable X; V~xi is the variance operator, which represents the variance calculation of the remaining variable after removing the ith independent variable X; S Ti is the impact of the ith independent variable X on the output target variable Y, including the interaction between independent variables.

[0054] The specific calculation process is obtained through Matlab programming, as follows:

[0055] Define E using the zeros function X With S Ti , its format is (D, 1), establish a for loop function, from 1 to N, calculate the E of each row X(i) Value, E X(i) =E X +(Y A -Y AB )^2;E X(i) / (2*N) is the E of each parameter variable X Value, V Y The definition is as above, so we can calculate S Ti =E X / V Y .

[0056] Step f: The SOBOL sensitivity index (S i and S Ti ) Sort by parameter size. If the sensitivity index is greater than the set value (the sensitivity index setting value can be in the range of 0.1-0.5 based on experience), the node overflow of the sub-catchment area is determined to be a sensitive variable;

[0057] Step g: The above steps complete the construction of the node overflow function relationship and screen out sensitive variables. Based on the platEMO platform, the platform's NSGA-II evolutionary algorithm is used to solve the solution cost and node overflow. The algorithm is evaluated no less than 10,000 times during the solution process (generally 10,000 times is used as the standard), and each sub-catchment area is adjusted according to sensitivity.

[0058] The adjustment process is as follows: under the constraints, first reduce the layout area of ​​LID facilities in the remaining sub-catchments that are not sensitive variables (such as sub-catchments S2, S4, S27, etc.) to 0, and then gradually increase the layout area of ​​LID facilities in the sub-catchments that are sensitive variables by the rated value (usually each increase in layout area is 0.1 hm 2 ), each time the layout area is increased, the cost function paper1_problemCB2(x) and the node overflow function paper1_problemYL2(x) in steps b and c are used to calculate the adjusted node overflow value and cost value, and then compared with the Pareto curve before the evolutionary algorithm is adjusted. If the adjusted plan is closer to the coordinate origin, it is proved to be better than the plan before the evolutionary algorithm is adjusted. Finally, the adjustment is stopped until the additional area causes the LID layout scale to exceed the total area of ​​the region, and then the Pareto solution result after the optimization of the NSGA-II algorithm can be obtained.

[0059] Compared to existing technologies, this method, based on the PlatEMO platform, uses the SOBOL method to screen key sensitive parameters and uses this as a guide for the platform's NSGA-II algorithm. This method increases the LID layout area in sensitive areas to improve the effectiveness of reducing node overflows, while reducing the LID layout area in insensitive areas to reduce the cost of the solution. Furthermore, this method improves the algorithm's computational efficiency by 20%. Compared to other SWMM models used in sponge city LID solution construction, this method achieves higher computational efficiency and more rational LID spatial layout in terms of LID spatial layout. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 This is a schematic diagram of the SWMM model established by taking Xiushan City as an example in a specific implementation method.

[0061] Figure 2 Schematic diagram of a list of 81 independent variables in a specific implementation method.

[0062] Figure 3 In a specific implementation manner, a schematic diagram of a rainfall process curve is shown, taking Xiushan County, Chongqing as an example, which is based on the Chicago rain type to generate a designed rainfall process line.

[0063] Figure 4 FIG. 1 is a schematic diagram of the SOBOL index ranking of LID facility spatial layout parameters in Xiushan County, Chongqing, as an example in a specific implementation manner.

[0064] Figure 5 Schematic diagram of the Pareto approximation frontier curve of the spatial layout of LID facilities when NSGA-II is evaluated 10,000 times in a specific implementation method.

[0065] Figure 6 This is a comparison table of the construction costs and node overflow control effects of the two schemes in the specific implementation method.

[0066] Figure 7 Schematic diagram of the Pareto approximate frontier curve of the spatial layout of LID facilities before and after continuous optimization of sensitivity parameters in a specific implementation method. DETAILED DESCRIPTION

[0067] The present invention will be further described in detail below with reference to specific embodiments.

[0068] Specific implementation method: A method for optimizing the layout of urban LID facilities based on SOBOL sensitivity analysis, which is characterized in that a SWMM model of the urban area is first established to study the function relationship between the node overflow and the construction cost of the LID scheme in the SWMM model, and the SOBOL method is used to conduct a global sensitivity analysis on the layout scale parameters of the LID facilities in each sub-catchment area to examine the significance of the impact of the parameters on the node overflow; then the NSGA-II algorithm in the MATLAB software platEMO4.0 platform is used to find the Pareto optimal solution set of the spatial layout of the LID facilities, and the LID facility layout area of ​​the sensitive variables selected is adjusted and optimized according to the results of the parameter sensitivity analysis, and the Pareto optimal solution set before and after the parameter optimization is compared and analyzed to obtain a better urban LID facility layout scheme.

[0069] Thus, building on the existing interactive method based on the SWMM model and the platEMO platform described in the background, this method conducts a global sensitivity analysis of LID facility layout scale parameters. LID facility layout areas are adjusted based on the impact of these parameters on node overflows. After comparative analysis, a more optimal urban LID facility layout solution is obtained. This improves the quality of the optimal LID facility spatial layout solution under the control objectives and enhances the efficiency of solving multi-objective optimization problems. Node overflow refers to the amount of water overflowing from a node in the urban drainage system (such as a manhole) during rainfall when the node exceeds its design flow rate.

[0070] In this embodiment, the sponge city construction area in Xiushan, Chongqing is taken as an example for specific description. The method includes the following steps:

[0071] Step a: Obtain basic data related to modeling of the study city area and establish a SWMM model of the city in the SWMM software; based on the SWMM software, divide the sub-catchments according to the direction of the city's pipe network, buildings and street distribution, and divide the study area into x sub-catchments; complete the parameter calibration of the SWMM model; finally, successfully construct a complete SWMM model of the city and export the inp file, as shown in Figure 1.

[0072] Step b: Deploy LID facilities in each subcatchment area, with the number of LID facilities being y types (further, the deployed LID facilities include but are not limited to three types of LID facilities: permeable pavement (PP), bioretention basins (BRC), and rain gardens (RG)). The deployment area of ​​y types of LID facilities in x subcatchments is used as the model independent variable, for a total of z independent variables, z = x·y; (In the implementation, the Chongqing Xiushan sponge city construction area is taken as an example, which has 27 subcatchments and a total of 81 independent variables, such as Figure 2 As shown in the table, based on MATLAB software, a multi-objective optimization model was established with the LID facility construction cost and node overflow as two objective functions. The functional relationship of the LID facility construction cost was established on the multi-objective optimization model, namely:

[0073] (1)

[0074] Where: f1 is the construction cost function of LID facilities in the study area, in yuan; i is the number of sub-catchments, i = 1, 2, 3, ..., x; j is the type of LID facilities, j = 1, 2, ..., y; S ij is the area of ​​the jth type of LID facility in the ith subcatchment area; M j is the unit area construction cost of the jth type of LID facility (the unit area construction cost of each type of LID facility is obtained based on the relevant engineering data of the urban area construction in the study);

[0075] Use matlab's function command to encapsulate the construction cost function;

[0076] The specific process of encapsulating the construction cost function during implementation is to use the function function (declare the function name, input and output) to define the cost function paper1_problemCB2(mian_ji); determine zhsqy_number (the number of sub-catchment areas), whose value is x (27), which is one y-fold of the total z (81) model independent variables, that is, length(mian_ji) / y; use the reshape function (reconstruct the array) to re-sort the area matrix according to y rows and x columns, and define it as CB_before; finally, calculate, construction cost = CB_before×[unit area construction cost of various LID facilities], that is, the area matrix of z model independent variables multiplied by the array matrix composed of the unit area construction cost of each type of LID facility, and use the sum function (sum function) to sum it up, that is, complete the construction of the script function relationship of the construction cost;

[0077] Step c: Establish a functional relationship for node overflow;

[0078] (2)

[0079] Where: f2 is the node overflow, unit is m 3 / s; g is the SWMM model calculation and result extraction function, which is obtained by Matlab programming;

[0080] Use matlab's function command to encapsulate the node overflow function;

[0081] The specific process of encapsulating the node overflow function during implementation is as follows:

[0082] 1) Start reading and rewriting data from the inp file;

[0083] First, define the peak flow function paper1_problemYL2(mian_ji) using the function function (declare the function name, input and output), and set it in a folder. mian_ji corresponds to S_ij in equations (1) and (2). Open the inp file of the SWMM model in a readable and writable way, use the fgetl function (read data line by line) to read the original file data line by line and create a new object newline to receive each line of data in the original file.

[0084] 2) Determine the number of LID facilities (LID-number) and the number of lines occupied by LID in the SWMM model inp file;

[0085] Use the contains function (query parameter position) to read the lia1 and lia2 parameters of the line where the strings "LID_USAGE" and "JUNCTIONS" are located. In the SWMM model, the difference between lia1 and lia2 minus 4 is the number of LID facilities; this number corresponds to j in equations (1) and (2), which is fixed here as y;

[0086] 3) Determine zhsqy_number and the number of subcatchments in the inp file of the SWMM model;

[0087] Use the contains function (query parameter location) to read the lia3 and lia4 parameters in the row where the strings "SUBCATCHMENTS" and "SUBAREAS" are located. In the SWMM model, the difference between lia3 and lia4 minus 4 is the number of subcatchments.

[0088] 4) Get the area and name of the subcatchment;

[0089] Use the cell function (create cell array) to create an empty cell array with n=0, and use the zeros function (create zero matrix) to initialize the area matrix; then establish a for loop with the parameter isi, extract the area of ​​the sub-catchment area from 1 to the number of zhsqy_number, and determine the area position of the SWMM model inp file (such as {lia33+3-1+isi} row (51:59) column) (same as above); use the strip function (delete the whitespace at the beginning and end of the string) to remove the spaces in the file and keep only the text; finally, use the str2double function (digital text) to convert the text into a single string. Convert the subcatchment names to numeric values ​​(subcatchment names may not be continuous during implementation and need to be extracted for subsequent operations). When extracting names (similar to the area extraction step above), use the cell function to create an empty cell array of n=0 for the area, and use the zeros function to initialize the subcatchment name matrix. Then, establish a for loop with the parameter isi2 to extract the subcatchment names from 1 to zhsqy_number, and determine the name location in the SWMM model inp file (e.g., {lia33+3-1+isi2} rows (1:3) columns) (same as above). Finally, use the strip function to remove spaces in the file, retaining only the text.

[0090] 5) Modify the inp file of the SWMM model and splice it into a new variable; each subsequent calculation will read the corresponding data on the original inp file and generate a new inp file; among them, newline(1:lia1+2) is the beginning of the new file newline1, and newline(lia22-2+1:length(newline)) is the end; the modified part in the middle, cell(1,(3*zhsqy_number)), is a cell array of the number of LIDs created by the cell function. These three parts are spliced ​​into a new newline1 variable;

[0091] 6) Add y default LID facilities (such as PP, BRC and RG) for all subcatchments;

[0092] Specifically, we create a for loop with the parameter isi2, from 1 to the number of zhsqy_number; we use the cell2mat function (which converts the cell array into a normal array of the basic data type) and use the area matrix (mjjz) as a variable to convert it into a normal array of the basic data type, with the area initialized to 0;

[0093] 7) Calculate the sum of the LID areas in each subcatchment and set constraints;

[0094] Use the reshape function (reconstruct array) to convert the area matrix (mian_ji) into a readable format of y rows and x columns in column order, and use the sum function (sum function) to calculate the sum of the areas of the LID facilities in each sub-catchment area, represented by CB_before; establish a for loop with the parameter isi3, from 1 to zhsqy_number, the number of sub-catchments, and use the strrep function (replace the specified substring in the string) to replace the area parameter in newline1 with the set value (for example, {isi2*3+lia11+3-2} is replaced by 2); use the if statement to set the constraint condition in the for loop, CB_before should be less than the area matrix (mjjz), that is, the sum of the areas of each sub-catchment area. If it exceeds the sum of the sub-catchment areas, it is set to 0 and the LID is not set; (This constraint is similar to the platEMO constraint condition later, but if this condition is not set, even if the platEMO constraint condition is set correctly, due to the instability of its constraints, some values ​​that exceed the constraints will occur when the algorithm is randomly initialized in subsequent operations, which will cause SWMM to fail to run)

[0095] 8) Run SWMM to perform calculations and extract the OUT file of the SWMM model. (In implementation, the SWMM interface can be called according to the patent developed by the research group.) Use the fopen function (read and write file data) to open the OUT file after the calculation in a readable and writable manner. Then use the fgetl function to read the file data line by line and create a new object newline_out to receive each line of data from the original file.

[0096] 9) Determine the overflow location of the sub-catchment node in the OUT file report and calculate the cumulative node overflow volume; use the contains function to read the lia_out and lia2_out parameters where the strings "Node Flooding Summary" and "Storage Volume Summary" are located. The number of rows to be extracted is from lia_out+10 to lia2_out-4 in newline_out; establish a for loop with parameter i23, extract the above-mentioned position rows line by line from 1 to all the rows to be extracted, and determine that the node overflow volume value is the seventh group of data (this is certain because this column of the SWMM OUT file is the result value of the node overflow volume); finally, use the str2double function to convert all extracted strings into numbers, and use the sum function to sum them, thus completing the script function relationship construction of the node overflow volume.

[0097] Step d: The above two steps develop corresponding programs that encapsulate the two functions of construction cost and node overflow, namely paper1_problemYL2(x) and paper1_problemCB2(x). This completes the preparation work of the correlation between construction cost and node overflow, making it meet the data interface format of the platEMO platform, facilitating the subsequent use of the platform's NSGA-Ⅱ evolutionary algorithm to solve the solution cost and node overflow.

[0098] Based on the platEMO platform, (in order to successfully call the platform's evolutionary algorithm, it is necessary to set the corresponding constraints, that is, the LID facility area cannot be larger than the total area of ​​the subcatchment. If the sum of the set LID areas exceeds the subcatchment area, SWMM will not be able to successfully simulate, causing the connection between SWMM and platEMO to crash. This step will constrain the variables to meet the input conditions of the SWMM calculation engine; therefore, it is necessary to set the constraints as follows:

[0099] (3)

[0100] Where: is the area of ​​the ith subcatchment;

[0101] According to the steps for solving the customized problem of the platEMO platform, the parameter encoding (the encoding method of each variable) is first used to define the y (3) types of LID facilities on the x (27) sub-catchment areas, that is, the encoding method of a total of z (81) variables, and the ones function (creating a matrix with 1 as the value) is used to create a matrix with 1 row and z columns consisting of 1 values; then the parameter lower is used to define (the lower bound of each variable) the lower bound of the z independent variables, and the zeros function (creating a matrix with 0 as the value) is used to create a matrix with 1 row and z columns consisting of 0 values, that is, the minimum layout area of ​​each LID facility is 0; the parameter upper (the upper bound of each variable) is used to define the upper bound of each independent variable (such as the sub-catchment area S1, whose area is 3.1979hm 2 , then the upper limit of the area of ​​each LID facility in the S1 area is 3.1979hm 2 , and the lengths of the lower and upper bound matrices created must be the same as the 1-valued matrices defined by the encoding parameter. Using the parameter objFcn (which minimizes the objective function), minimize the two problem functions, construction cost and node overflow. Finally, using the parameter conFcn (constraint function), the constraint is satisfied if and only if the constraint violation value is less than or equal to zero (for example, in area S1, the sum of the areas of the y types of LID facilities minus the total area of ​​S1, 3.1979, should be less than or equal to zero). Repeat the above steps to complete the processing of other subcatchments. This allows the subsequent simulation analysis of the functional relationship between construction cost and node overflow to be successfully invoked using the evolutionary algorithm of the platEMO platform.

[0102] Step e: Use SOBOL sensitivity analysis method to perform sensitivity analysis on the z independent variable parameters in step a. The analysis process is obtained by Matlab programming (this step can be performed simultaneously with steps b, c, and d). The details are as follows:

[0103] Since there are z parameters, the parameter dimension D is set to z. To ensure symmetrical sampling, the parameter dimension M is expanded to D*2. The number of sampling points N is set to 5000, which means 5000 sampling groups. The upper and lower limits of each parameter are set. The zeros function is used to set a 1-row, 81-column zero matrix as the lower limit of each parameter, named MIN. The upper limit of each parameter cannot exceed the total area of ​​the subcatchment, named MAX. Since each catchment has three LIDs, the MAX array is divided by 3 to ensure uniform sampling.

[0104] 2. Use the sobolset function to generate uniformly distributed sample points according to the parameter dimension M, named p; establish a for loop function to screen the sample points from 1 to N, and obtain the sample of p one by one, named r, with the range r=Min + r.*(Max-Min); finally, distribute the sampling results to the matrix R, that is, the N×M parameter matrix.

[0105] 3 Split the matrix R into matrix A and matrix B; the first D columns of matrix R are divided into matrix A, that is, the format is N×D, and the columns from D+1 to M are divided into matrix B, that is, the format is N×D. Construct matrix AB based on matrices A and B i (i=1, 2, ..., D). The specific process is to establish a for loop function to replace the i-th column in matrix B with the i-th column in matrix A from 1 to D, while the other columns remain unchanged. This can be replaced D times, for a total of D AB matrices.

[0106] 4 The above steps give A, B, AB i (i=1, 2, ..., D) These matrices, totaling (D+2)*N groups of input data, can calculate (D+2)*N groups of Y values; input (D+2)*N groups of input data into the node overflow function constructed in step a to solve for Y A 、Y B 、Y ABi (i=1, 2,…, D) values.

[0107] 5 Calculate the main effect index S i , the formula is:

[0108] (1)

[0109] Where: X i is the i-th independent variable; Y is the target variable; V represents the average value of the output target variable Y under the condition of inputting the i-th independent variable X; xi For the variance operator, input the independent variable X when other variables are fixed. i Variance calculation is performed; V(Y) is the total variance of the output target variable Y, which represents the variance of all independent variables X1, X2, ..., X k The variance of the output target variable Y under the joint action; S i It is the first-order sensitivity index, that is, the direct contribution of the ith independent variable to the output target variable Y, without considering the interaction between the independent variables.

[0110] The specific calculation process is obtained through Matlab programming, as follows:

[0111] Define V using the zeros function X With S i, its format is (D, 1), establish a for loop function, from 1 to N, calculate V for each row X(i) Value, V X(i) = V X +Y B *(Y AB -Y A );V X(i) / N is the V of each parameter variable X Value, V Y The format is defined as ([Y A ;Y B ],1),[Y A ;Y B ] is Y A and Y B The matrix is ​​stacked; finally, S is calculated i =V X / V Y .

[0112] 6Calculate the total effect index S Ti , the formula is:

[0113] (2)

[0114] Where: X~i is the remaining variable after removing the ith independent variable X; E[Y|X_(~i)] represents the average value of the output target variable Y under the condition of the remaining variable after removing the ith independent variable X; V~xi is the variance operator, which represents the variance calculation of the remaining variable after removing the ith independent variable X; S Ti is the impact of the ith independent variable X on the output target variable Y, including the interaction between independent variables.

[0115] The specific calculation process is obtained through Matlab programming, as follows:

[0116] Define E using the zeros function X With S Ti , its format is (D, 1), establish a for loop function, from 1 to N, calculate the E of each row X(i) Value, E X(i) =E X +(Y A -Y AB )^2;E X(i) / (2*N) is the E of each parameter variable X Value, V Y The definition is as above, so we can calculate S Ti =E X / V Y .

[0117] Step f: The SOBOL sensitivity index (S i and S Ti ) Sort by parameter size. If the sensitivity index is greater than the set value (the sensitivity index setting value can be 0.5 based on experience), the node overflow of the sub-catchment area is judged to be a sensitive variable.

[0118] In the specific implementation, according to the rainstorm intensity formula (Xiushan County, Chongqing), the design rainfall process line is generated based on the Chicago rain pattern, the rainfall duration is 2 hours, and the design rain peak coefficient is 0.3. The rainfall process curve with a rainfall return period (P) of 50 years is as follows Figure 3 As shown;

[0119] Taking the return period P = 50a as an example, the SOBOL sensitivity analysis method is used to obtain the S of z (81) parameters. i With S Ti The value (range is 0~1) is used to calculate the SOBOL index (S i and S Ti ) are sorted by parameter sensitivity, and the results are as follows Figure 4 As shown (due to the constraints of this study, the three types of LID areas in each sub-catchment and the fact that they do not exceed the total area of ​​the sub-catchment may lead to uneven data distribution of a single LID facility, resulting in a small SOBOL index. Therefore, only parameters with a SOBOL index > 0.01 are displayed. Generally, an index greater than 0.5 is considered sensitive, and a value much smaller than 0.1 is considered insensitive). The results show that the deployment of LID facilities in sub-catchments S20, S10, S22, S15, S25, and S24 is more sensitive, while the SOBOL index of the remaining parameters is almost 0, which has no effect on the output results. Among them, the deployment of LID facilities in sub-catchment S20 has the most significant impact on node overflow, especially the deployment of PP and BRC facilities. S i The index is generally slightly smaller than S Ti The index indicates that there is no interaction effect between the parameter variables and they can independently affect the node overflow.

[0120] Step g: The above steps complete the construction of the node overflow function relationship and screen out sensitive variables (the sensitivity of the LID layout area in the S20 and S15 areas in the embodiment is much higher than that of the other parameters, and the SOBOL index of S20PP, S20BRC, and S15BRC is around 0.1, so S20PP, S20BRC, and S15BRC are selected as key adjustment variables); based on the platEMO platform, the NSGA-Ⅱ evolutionary algorithm of the platform is used to solve the solution cost and node overflow. During the solution process, the algorithm is evaluated no less than 10,000 times (generally 10,000 times is used as the standard), and each sub-catchment area is adjusted according to the sensitivity (the Pareto approximate frontier of the spatial layout of the LID facilities in the embodiment is as follows Figure 5 As shown in Figure 2 , each point on the Pareto optimal frontier contains a theoretically optimal LID facility spatial layout plan, which decision makers select based on actual conditions. The optimal solution with a maximum node overflow of 412 m³ and a construction cost of 68.8 million yuan was selected as the initial plan, Plan A. The LID layout areas at S20 and S15 were adjusted based on the results of the parameter sensitivity analysis.

[0121] The adjustment process is as follows: under the constraints, first reduce the layout area of ​​LID facilities in the remaining sub-catchments that are not sensitive variables (such as sub-catchments S2, S4, S27, etc.) to 0, and then gradually increase the layout area of ​​LID facilities in the sub-catchments that are sensitive variables (S20PP, S20BRC, S15BRC) by the rated value (usually each increase in layout area is 0.1 hm 2 ), each time the layout area is increased, the cost function paper1_problemCB2(x) and the node overflow function paper1_problemYL2(x) in steps b and c are used to calculate the adjusted node overflow value and cost value, and then compared with the Pareto curve before the evolutionary algorithm is adjusted. If the adjusted plan is closer to the coordinate origin, it is proved to be better than the plan before the evolutionary algorithm is adjusted. Finally, the adjustment is stopped until the additional area causes the LID layout scale to exceed the total area of ​​the region, and then the Pareto solution result after the optimization of the NSGA-II algorithm can be obtained.

[0122] When implementing Figure 5 The results show that as the layout area of ​​LID facilities in sensitive areas increases, the control effect of LID facility layout on node overflow is significantly improved. In addition, due to the reduction of the layout area of ​​LID facilities in insensitive areas, the construction cost of the LID facility layout scheme after parameter optimization is basically lower than that before parameter optimization.

[0123] Selecting the optimal solution with a construction cost close to that of Solution A as Solution B, the construction cost and node overflow control effect of the two solutions can be obtained. The results are as follows Figure 6 The results show that although the unit area cost of option B (343.92 yuan / m 2 ) is slightly higher than Plan A (329.70 yuan / m 2 ), but its control effect on node overflow is significantly stronger than that of Scheme A, with a node overflow reduction rate of up to 90.29%.

[0124] During implementation, the applicant further verified the effect of the scheme after parameter optimization and further explored the computational efficiency of the Pareto approximate front curve after sensitivity parameter optimization. The present invention continued to solve the problem based on 10,000 evaluations of the NSGA-II algorithm. The results showed that the Pareto approximate front curve after sensitivity parameter optimization almost coincided with the Pareto approximate front curve of the NSGA-II algorithm evaluated 12,000 times (see Figure 7 ). The results show that by focusing on optimizing sensitive sub-catchments based on parameter sensitivity analysis and eliminating insensitive sub-catchments before solving the LID facility spatial layout optimization model, the computing power of NSGA-II for the multi-objective model can be increased by 1.2 times. This parameter optimization method not only significantly saves the calculation time of the model algorithm, but also significantly reduces the construction cost of LID facilities while meeting overflow control. To further verify the effectiveness of this method, this study further optimized the LID facility spatial layout plan based on 12,000 evaluations of the NSGA-II algorithm. The study confirmed that the effect of sensitivity parameters on driving the NSGA-II algorithm to seek optimization is sustained, especially on the basis of 12,000 evaluations. It can still continue to find a better plan for the spatial layout of LID facilities based on the results of the sensitivity analysis, and continuously provide optimization direction for the spatial layout plan of LID facilities in the sub-catchment area.

[0125] Therefore, based on the PlatEMO platform, this study uses the SOBOL method to screen key sensitive parameters and uses this as a guide for the platform's NSGA-II algorithm. This method increases the LID footprint in sensitive areas to improve node overflow reduction, while reducing the LID footprint in insensitive areas to lower the project construction cost. This method also improves the algorithm's computational efficiency by 20%. Compared to other SWMM model studies on sponge city LID solutions, this study achieves higher computational efficiency and more rational LID spatial layout in terms of LID spatial layout.

Claims

1. A method for optimizing urban LID facility layout based on SOBOL sensitivity analysis, characterized in that: First, a SWMM model of the urban area was established. A functional relationship between node overflow and LID scheme construction cost was constructed in the SWMM model. The SOBOL method was used to conduct a global sensitivity analysis of the layout scale parameters of LID facilities in each sub-catchment area to examine the significance of the parameters' impact on node overflow. The NSGA-II algorithm in the MATLAB software platform EMO4.0 was then used to find the Pareto optimal solution set for the spatial layout of LID facilities. Based on the results of the parameter sensitivity analysis, the LID facility layout area of ​​the selected sensitive variables was adjusted and optimized. The Pareto optimal solution sets before and after parameter optimization were compared and analyzed to obtain a more optimal urban LID facility layout plan. This method comprises the following steps: Step a: Obtain basic data related to modeling of the study city area and build a SWMM model of the city in SWMM software. Using SWMM software, divide the study area into x subcatchments based on the city's pipe network, building, and street distribution. Complete the SWMM model parameter calibration. Finally, successfully build a complete SWMM model of the city and export the inp file. Step b: Deploy LID facilities in each subcatchment area, with the number of LID facilities being y. The deployment area of ​​y types of LID facilities in x subcatchments is used as the model independent variable, for a total of z independent variables, z = x·y. Using MATLAB software, a multi-objective optimization model is established with the LID facility construction cost and node overflow as the two objective functions. The functional relationship of the LID facility construction cost is established on the multi-objective optimization model, namely: Where: f1 is the construction cost function of LID facilities in the study area, in yuan; i is the number of sub-catchments, i = 1, 2, 3, ..., x; j is the type of LID facilities, j = 1, 2, ..., y; S ij is the area of ​​the jth type of LID facility in the ith subcatchment area; M j is the construction cost per unit area of ​​the jth type LID facility; Use matlab's function command to encapsulate the construction cost function; Step c: Establish a functional relationship for node overflow; minf2=g(S ij ) (2) Where: f2 is the node overflow, unit is m 3 / s; g is the SWMM model calculation and result extraction function, which is obtained by Matlab programming; Use matlab's function command to encapsulate the node overflow function; Step d: The above two steps develop corresponding programs that encapsulate the two functions of construction cost and node overflow, namely paper1_problemYL2(x) and paper1_problemCB2(x). This completes the preparation work of the correlation between construction cost and node overflow, making it meet the data interface format of the platEMO platform, facilitating the subsequent use of the platform's NSGA-Ⅱ evolutionary algorithm to solve the solution cost and node overflow. Based on the platEMO platform, the constraints are set as follows: Where: S′ i is the area of ​​the ith subcatchment; According to the steps for solving the custom problem on the platEMO platform, the encoding parameter is first used to define the y types of LID facilities on x sub-catchments, that is, the encoding method for a total of z variables, and the ones function is used to create a matrix with 1 row and z columns consisting of 1 values; then the lower parameter is used to define the lower bound of the z independent variables, and the zeros function is used to create a matrix with 1 row and z columns consisting of 0 values, that is, the minimum layout area of ​​each LID facility is 0; the upper parameter is used to define the upper bound of each independent variable; the objFcn parameter is used to minimize the two problem functions, construction cost and node overflow; finally, the conFcn parameter is used to satisfy the constraint if and only if the constraint violation value is less than or equal to zero, and the processing of other sub-catchments is completed according to the above steps; so that when solving the functional relationship between construction cost and node overflow in the subsequent process, the evolutionary algorithm of the platEMO platform can be successfully called for simulation analysis; Step e: Use SOBOL sensitivity analysis method to perform sensitivity analysis on the z independent variable parameters in step a. The analysis process is obtained by Matlab programming, as follows: 1) Since there are z parameters in total, the parameter dimension D is set to z. For symmetrical sampling, the parameter dimension M is expanded to D*2. The number of sampling points N is set to 5000, that is, the number of sampling groups is 5000. The upper and lower limits of each parameter are set. The zeros function is used to set a 1-row 81-column zero matrix as the lower limit of each parameter, named MIN. The upper limit of each parameter cannot exceed the total area of ​​the sub-catchment, named MAX. Since three types of LIDs are arranged in each catchment, the MAX array is divided by 3 for uniform sampling. 2) Using the sobolset function, generate uniformly distributed sample points according to the parameter dimension M, named p; establish a for loop function to screen the sample points from 1 to N, and obtain the sample of p one by one, named r, whose range is r = Min + r.*(Max-Min); finally, distribute the sampling results into the matrix R, that is, the N×M parameter matrix; 3) Split the matrix R into matrix A and matrix B; the first D columns of matrix R are divided into matrix A, that is, the format is N×D, and the columns from D+1 to M are divided into matrix B, that is, the format is N×D; construct matrix AB based on matrices A and B i (i=1, 2, ..., D), the specific process is to establish a for loop function from 1 to D to replace the i-th column in matrix B with the i-th column of matrix A, while the other columns remain unchanged. This can be replaced D times, for a total of D AB matrices; 4) The above steps yield A, B, and AB i (i=1,2,…,D) These matrices, totaling (D+2)*N groups of input data, can calculate (D+2)*N groups of Y values; input (D+2)*N groups of input data into the node overflow function constructed in step a to solve for Y A 、Y B 、Y ABi (i=1, 2, ..., D) value; 5) Calculate the main effect index S i , the formula is: Where: X i is the ith independent variable; Y is the target variable; E[Y|X i ] represents the average value of the target variable Y under the condition of inputting the i-th independent variable X; V xi For the variance operator, input the independent variable X when other variables are fixed. i Variance calculation is performed; V(Y) is the total variance of the output target variable Y, which represents the variance of all independent variables X1, X2, ..., X k The variance of the output target variable Y under the joint action; S i is the first-order sensitivity index, that is, the direct contribution of the i-th independent variable to the output target variable Y, without considering the interaction between independent variables; The specific calculation process is obtained through Matlab programming, as follows: Define V using the zeros function X With S i , its format is (D, 1), establish a for loop function, from 1 to N, calculate V for each row X(i) Value, V X(i) =V X +Y B *(Y AB -Y A );V X(i) / N is the V of each parameter variable X Value, V Y The format is defined as ([Y A ; Y B ],1),[Y A ; Y B ] is Y A and Y B The matrix is ​​stacked; finally, S is calculated i =V X / V Y ; 6Calculate the total effect index S Ti , the formula is: Where: X~i is the remaining variable after removing the ith independent variable X; E[Y|X_(~i)] represents the average value of the output target variable Y under the condition of inputting the remaining variables after removing the ith independent variable X; V~xi is the variance operator, which represents the variance calculation of the remaining variables after removing the ith independent variable X; S Ti is the impact of the ith independent variable X on the output target variable Y, including the interaction between independent variables; The specific calculation process is obtained through Matlab programming, as follows: Define E using the zeros function X With S Ti , its format is (D, 1), establish a for loop function, from 1 to N, calculate the E of each row X(i) Value, E X(i) =E X +(Y A -Y AB )^2;E X(i) / (2*N) is the E of each parameter variable X Value, V Y The definition is as above, so we can calculate S Ti =E X / V Y ; Step f: Sort the SOBOL sensitivity indexes of the x sub-catchments calculated by analysis with node overflow as the target variable by parameter size. If the sensitivity index is greater than the set value, the node overflow of the sub-catchment is determined to be a sensitive variable. Step g: The above steps complete the construction of the node overflow function relationship and screen out sensitive variables. Based on the platEMO platform, the platform's NSGA-II evolutionary algorithm is used to solve the solution cost and node overflow. The algorithm is evaluated no less than 10,000 times during the solution process, and each subcatchment area is adjusted according to sensitivity. The adjustment process is as follows: under the constraints, first reduce the layout area of ​​the remaining LID facilities in the sub-catchment areas that are not sensitive variables to 0, and then gradually increase the layout area of ​​the LID facilities in the sub-catchment areas that are sensitive variables by the rated value. Each time the layout area is increased, the cost function paper1_problemCB2(x) and the node overflow function paper1_problemYL2(x) in steps b and c are used to calculate the adjusted node overflow value and cost value, and compare them with the Pareto curve of the evolutionary algorithm before adjustment. If the adjusted scheme is closer to the coordinate origin, it is proved to be better than the scheme before adjustment of the evolutionary algorithm. Finally, stop until the additional area will make the LID layout scale exceed the total area of ​​the area, and then the Pareto solution result optimized by the NSGA-II algorithm can be obtained.

2. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 1, characterized in that: In step b, the LID facilities installed include three types of LID facilities: permeable pavement PP, biological retention pond BRC and rain garden RG.

3. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 1, characterized in that: In step b, the process of encapsulating the construction cost function using MATLAB function commands is as follows: use the function function to define the cost function paper1_problemCB2(mian_ji); determine the number of sub-catchment areas zhsqy_number, whose value is x, which is one y-th of the total z model independent variables, that is, length(mian_ji) / y; use the reshape function to reorder the area matrix according to y rows and x columns, and define it as CB_before; finally, calculate, construction cost = CB_before × [unit area construction cost of various LID facilities], that is, the area matrix of z model independent variables multiplied by the array matrix composed of the construction cost of each unit area of ​​various LID facilities, and use the sum function to sum, that is, complete the construction of the script function relationship of the construction cost.

4. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 1, characterized in that: In step c, the process of encapsulating the node overflow function using MATLAB function commands is as follows: 1) Start reading and rewriting data of inp file; First, define the peak flow function paper1_problemYL2(mian_ji) using the function function. The input is a 3×4 variable matrix length(blqy) and is set in a folder. mian_ji corresponds to Sij in equations (1) and (2). Open the inp file of the SWMM model in a readable and writable manner, use the fgetl function to read the original file data line by line and create a new object newline to receive each line of data in the original file. 2) Determine the number of LID facilities LID_number and the number of lines occupied by LID in the SWMM model inp file; Use the contains function to read the lia1 and lia2 parameters of the row where the strings "LID_USAGE" and "JUNCTIONS" are located. In the SWMM model, the difference between lia1 and lia2 minus 4 is the number of LID facilities; this number corresponds to j in equations (1) and (2), which is fixed here as y. 3) Determine zhsqy_number and the number of subcatchments in the inp file of the SWMM model; Use the contains function to read the lia3 and lia4 parameters in the row where the strings "SUBCATCHMENTS" and "SUBAREAS" are located. In the SWMM model, the difference between lia3 and lia4 minus 4 is the number of subcatchments. 4) Get the area and name of the subcatchment area; The cell function is used to create an empty cell array with n=0 for the area, and the zeros function is used to initialize the area matrix. Then, a for loop with the parameter isi is established to extract the area of ​​the subcatchment from 1 to zhsqy_number, and at the same time determine the area location of the SWMM model inp file. The strip function is used to remove spaces in the file, retaining only text. Finally, the str2double function is used to convert the text into a numeric value. When extracting the name, the cell function is used to create an empty cell array with n=0 for the area, and the zeros function is used to initialize the subcatchment name matrix. Then, a for loop with the parameter isi2 is established to extract the name of the subcatchment from 1 to zhsqy_number, and determine the name location of the SWMM model inp file. Finally, the strip function is used to remove spaces in the file, retaining only text. 5) Modify the inp file of the SWMM model and splice it into a new variable; each subsequent calculation will read the corresponding data on the original inp file and generate a new inp file; among them, newline(1:lia1+2) is the beginning of the new file newline1, and newline(lia22-2+1:length(newline)) is the end; the modified part in the middle, cell(1,(3*zhsqy_number)), is a cell array of the number of LIDs created by the cell function. These three parts are spliced ​​into a new newline1 variable; 6) Add y default LID facilities for all subcatchments; Specifically, a for loop is established for the parameter isi2, from 1 to the number of zhsqy_number; the cell2mat function is used to convert the cell array into a normal array of the basic data type, and the area matrix mjjz is used as a variable to convert it into a normal array of the basic data type, and the area is initialized to 0; 7) Calculate the sum of the LID areas in each subcatchment and set constraints; Use the reshape function to convert the area matrix mian_ji into a readable format with y rows and x columns in column order, and use the sum function to calculate the sum of the areas of the LID facilities in each subcatchment, represented by CB_before. Create a for loop with the parameter isi3, from 1 to zhsqy_number (the number of subcatchments), and use the strrep function to replace the area parameter in newline1 with the set value. Use the if statement to set the constraint condition in the for loop: CB_before should be less than the area matrix mjjz, that is, the sum of the areas of each subcatchment. If it exceeds the sum of the subcatchment areas, it is set to 0 and no LID is set. 8) Run SWMM to perform the calculation and extract the OUT file of the SWMM model; use the fopen function to open the OUT file after the calculation in a readable and writable manner, then use the fgetl function to read the file data line by line and create a new object newline_out to receive each line of data from the original file; 9) Determine the overflow locations of the subcatchment nodes in the OUT file report and calculate the cumulative node overflow volume. Use the contains function to read the lia_out and lia2_out parameters of the rows where the strings "Node Flooding Summary" and "Storage Volume Summary" are located. The number of rows to be extracted is from lia_out+10 to lia2_out-4 in newline_out. Create a for loop with parameter i23, extracting the above-mentioned rows line by line from 1 to all the rows to be extracted, and determine the node overflow volume value as the seventh set of data. Finally, use the str2double function to convert all extracted strings into numbers and use the sum function to sum them, thus completing the script function relationship construction of the node overflow volume.

5. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 1, characterized in that: In step f, the sensitivity index is set to 0.1-0.

5.

6. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 1, characterized in that: In step g, the rated value of the layout area increased each time is 0.1hm2.

Images