Urban LID facility layout optimization method based on SOBOL sensitivity analysis

Abstract

The invention discloses a city LID facility layout optimization method based on SOBOL sensitivity analysis, and the method is characterized in that the method comprises the steps: firstly building an SWMM model for researching a city region, building a function relation between the node overflow amount and the LID scheme construction cost in the SWMM model, carrying out the global sensitivity analysis of the layout scale parameter of the LID facility of each sub catchment region through employing an SOBOL method, and carrying out the overall sensitivity analysis of the layout scale parameter of the LID facility of each sub catchment region. The LID facility layout area of the screened sensitive variables is adjusted and optimized according to a parameter sensitivity analysis result, and a Pareto optimal solution set before and after parameter optimization is compared and analyzed to obtain a better urban LID facility layout scheme. Compared with the prior art, the method has the advantages that the reduction effect of the node overflow quantity is improved by increasing the LID layout area of the sensitive area, and the scheme construction cost is reduced by reducing the LID layout area of the insensitive area; therefore, the calculation efficiency is higher, and the LID space layout is more reasonable.

Description

Technical Field

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

[0002] The process of urbanization has led to the continuous expansion of impervious surfaces, resulting in a sharp increase in surface runoff. In order to reduce the runoff during heavy rains, reduce the pressure on urban drainage systems, alleviate urban waterlogging and control runoff pollution, my country has proposed the concept of sponge city construction, all of which are to simulate natural processes through low-impact development (LID) approaches to reduce the negative impact of urban development on the hydrological environment. In the construction of sponge cities, LID facilities such as rain gardens (RG), bioretention ponds (BRC) and permeable pavement (PP) are usually used, and LID facilities are combined and arranged in each sub-catchment area. However, its spatial layout depends on factors such as the number of sub-catchments, the type of LID facilities, the scale, area and spatial location, which in turn affects the effectiveness of urban waterlogging mitigation and the cost of sponge city construction. Since the type, scale, area and spatial location of LID facilities can form countless spatial layout combination schemes, how to seek the optimal spatial layout scheme at the lowest construction cost to achieve the best construction effect of sponge cities is a multi-objective optimization problem.

[0003] At present, solving multi-objective optimization problems mainly relies on multi-objective evolutionary algorithms (MOEA), that is, solving nonlinear, multi-dimensional complex optimization problems to find Pareto optimal solutions. However, for solving multi-objective optimization problems, MOEA often requires tens of thousands of evaluations. Especially when dealing with problems involving a large number of decision variables, more parameters and complex models usually lead to a sharp decline in algorithm performance and take too long.

[0004] The applicant has applied for a patent "An interactive method based on SWMM model and platEMO platform". The method includes the following steps: first establish a SWMM model of the urban area, divide the sub-catchment area and export the inp file; deploy LID facilities in each sub-catchment area, and use each LID facility in each sub-catchment area as an independent variable; based on matlab software, establish a multi-objective optimization model with LID facility construction cost and peak flow as two objective functions; use function commands to encapsulate construction cost and peak flow functions; then based on the platEMO platform, set constraints, follow the platform's custom problem solving steps, use the platform's existing function functions to minimize the construction cost and peak flow of the two problem functions, and complete the processing of each sub-catchment area one by one. This invention realizes the interaction between the SWMM model and the platEMO platform, and facilitates the use of the platform's various optimization algorithms to improve the optimization of the urban LID spatial layout. However, the method of the invention still has the defects of too many algorithm evaluation times, resulting in low algorithm efficiency, reduced quality of Pareto optimal solution, and reduced calculation and analysis accuracy. Summary of the invention

[0005] In view of the above-mentioned deficiencies of 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 performance and the quality of Pareto optimal solution, 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 in that a SWMM model of the urban area is first established, a functional relationship between the node overflow and the LID scheme construction cost is constructed in the SWMM model, and a global sensitivity analysis is performed on the layout scale parameters of the LID facilities in each sub-catchment area using the SOBOL method to examine the significance of the parameters' influence 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 LID facility spatial layout, and the LID facility layout area of ​​the sensitive variables selected is adjusted and optimized according to the parameter sensitivity analysis results, and the Pareto optimal solution set before and after parameter optimization is compared and analyzed to obtain a better urban LID facility layout plan.

[0008] In this way, based on the existing interactive method based on the SWMM model and the platformEMO platform described in the background technology, this method conducts a global sensitivity analysis on the layout scale parameters of the LID facilities, adjusts the layout area of ​​the LID facilities according to the size of the impact of the parameters on the node overflow, and obtains a better urban LID facility layout plan after comparative analysis. This improves the quality of the optimal spatial layout plan of the LID facilities under the control target and improves the efficiency of solving multi-objective optimization problems. The node overflow refers to the amount of water overflowing from the node when the node (such as the inspection well) of the urban drainage system exceeds its design flow during rainfall.

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

[0010] Step a: Obtain the basic data related to modeling of the urban area under study, and establish the SWMM model of the city in the SWMM software; Based on the SWMM software, divide the sub-catchment area according to the urban pipe network direction, buildings and street distribution, and divide the study area into x sub-catchment areas; complete the parameter calibration of the SWMM model; Finally, successfully build a complete SWMM model of the city and export the inp file;

[0011] Step b: LID facilities are laid out in each sub-catchment area, and the number of LID facilities is y (further, the laid LID facilities include but are not limited to three types of LID facilities, namely permeable pavement (PP), biological retention pond (BRC) and rain garden (RG); the layout area of ​​y types of LID facilities in x sub-catchments is used as the model independent variable, with a total of z independent variables, z = x y; based on MATLAB software, a multi-objective optimization model is established with LID facility construction cost and node overflow as two objective functions, and a functional relationship of LID facility construction cost is established on 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-catchment areas, i=1, 2, 3, ..., x; j is the type of LID facility, 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 (in implementation, the unit area construction cost of various types of LID facilities is obtained based on the relevant engineering data of the urban area construction);

[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 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, that is, complete the construction of the script function relationship of the construction cost;

[0016] Step c: Establishing the functional relationship of 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, 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:

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

[0022] First, use the function function (declare the function name, input and output) to define the peak flow function paper1_problemYL2(mian_ji), the input is a 3×4 variable matrix (length(blqy)), and it is set 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 inp file of the SWMM model; Use the contains function (query parameter position) to read the lia1 and lia2 parameters of the string "LID_USAGE" and "JUNCTIONS". 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 as y here;

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

[0025] Use the contains function (query parameter position) to read the lia3 and lia4 parameters of the string "SUBCATCHMENTS" and "SUBAREAS". In the SWMM model, the difference between lia3 and lia4 minus 4 is the number of subcatchment areas.

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

[0027] Use the cell function (create cell array) to create an empty cell array with n=0 for the area, and use the zeros function (create zero matrix) to initialize the area matrix; then establish a for loop with parameter isi to extract the area of ​​the subcatchment from 1 to zhsqy_number, and determine the area position of the inp file of the SWMM model (such as {lia33+3-1+isi} row (51:59) column) (same as above); use the strip function (delete the blanks 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. This is converted to a numerical value (the names of the sub-catchment areas are not necessarily continuous during implementation and need to be extracted for subsequent operations); when extracting the names (similar to the above area extraction steps), the cell function is used to create an empty cell array of n=0 for the area, and the zeros function is used to initialize the sub-catchment area name matrix; then a for loop with the parameter isi2 is established to extract the names of the sub-catchment areas from 1 to the number of zhsqy_number, and determine the name position of the SWMM model inp file (such as {lia33+3-1+isi2} row (1:3) column) (same as above); finally, the strip function is used to remove the spaces in the file and only retain the text;

[0028] 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 the cell array of the number of LIDs created by the cell function, and these three parts are spliced ​​into a new newline1 variable;

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

[0030] 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;

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

[0032] Use the reshape function (reconstruct array) to convert the area matrix (mian_ji) into a read 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 of the parameter isi3, from 1 to zhsqy_number the number of sub-catchments, and use the strrep function (specify the substring in the replacement 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 areas of the sub-catchments, 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)

[0033] 8) Run SWMM to perform calculations and extract the OUT file of the SWMM model; (When implementing, the SWMM interface can be called for calculations according to the patent developed by the previous 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 accept each line of data in the original file;

[0034] 9) Determine the overflow position of the sub-catchment area node in the OUT file report and calculate the cumulative node overflow; 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 position rows line by line from 1 to all the rows to be extracted, and determine that the node overflow value is the seventh group of data (this is certain because this column of the OUT file of SWMM is the result value of the node overflow); finally, use the str2double function to convert all extracted strings into numbers, and use the sum function to sum them, thus completing the construction of the script function relationship of the node overflow.

[0035] Step d: The above two steps develop corresponding programs to encapsulate the two functions of construction cost and node overflow, namely paper1_problemYL2(x) and paper1_problemCB2(x), completing the preparation work for the correlation between construction cost and node overflow, so that it meets the data interface form of the platEMO platform, and facilitates the subsequent use of the platform's NSGA-Ⅱ evolutionary algorithm to solve the solution cost and node overflow;

[0036] Based on the platEMO platform, (in order to successfully call the evolutionary algorithm of the platform, it is necessary to set the corresponding constraints, that is, the LID facility area cannot be larger than the total area of ​​the subcatchment area. If the sum of the set LID areas exceeds the subcatchment area, SWMM will not be able to simulate successfully, causing the connection between SWMM and platEMO to collapse. 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:

[0037] (3)

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

[0039] According to the steps for solving the custom 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 (creates a matrix with 1 as the value) is used to create a matrix with 1 row and z columns composed of 1 values; then the parameter lower is used to define (the lower limit of each variable) the lower limit of the z independent variables, and the zeros function (creates a matrix with 0 as the value) is used to create a matrix with 1 row and z columns composed of 0 values, that is, the minimum layout area of ​​each LID facility is 0; the parameter upper (the upper limit of each variable) is used to define the upper limit 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 length of the lower bound matrix and the upper bound matrix created must be the same as the matrix with a value of 1 defined by the parameter encoding); use the parameter objFcn (to minimize the objective function) to minimize the two problem functions, construction cost and node overflow; finally, use the parameter conFcn (constraint function) to satisfy the constraint only when and 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 deployed minus the total area of ​​S1, 3.1979, should be less than or equal to zero), and complete the processing of other sub-catchment areas according to the above steps. In this way, 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.

[0040] 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:

[0041] 1 Since there are z parameters in total, the parameter dimension D is set to z. For symmetrical sampling, the parameter dimension M=D*2 is expanded; the number of sampling points N=5000 is set, 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 area, named MAX. Since three types of LIDs are arranged in each catchment area, the MAX array is divided by 3 for uniform sampling.

[0042] 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 sample points from 1 to N, and obtain the samples of p one by one, named r, with a range of r=Min + r.* (Max-Min); finally, distribute the sampling results to the matrix R, that is, the N×M parameter matrix.

[0043] 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, and the other columns remain unchanged. This can be replaced D times, for a total of D AB matrices.

[0044] 4 The above steps give A, B, AB i (i=1, 2, ..., D) These matrices, a total of (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.

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

[0046] (1)

[0047] 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 ith independent variable X; xi is a variance operator. When other variables are fixed, input the independent variable X i 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.

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

[0049] Define V using the zeros function X With S i, whose 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 .

[0050] 6 Calculation of total effect index S Ti , the formula is:

[0051] (2)

[0052] 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 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.

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

[0054] Define E using the zeros function X With S Ti , whose format is (D, 1), establish a for loop function, from 1 to N, calculate E for 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 .

[0055] Step f: The SOBOL sensitivity index (S i and S Ti ) are sorted 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 judged to be a sensitive variable;

[0056] 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 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;

[0057] The adjustment process is as follows: under the constraint conditions, firstly, the layout area of ​​LID facilities in the remaining sub-catchment areas that are not sensitive variables (such as sub-catchment areas S2, S4, S27, etc.) is reduced to 0, and then the layout area of ​​LID facilities in the sub-catchment areas that are sensitive variables is gradually increased by the rated value (usually each increase in the layout area is 0.1 hm 2 ), each time the layout area is increased, the adjusted node overflow value and cost value are calculated using the cost function paper1_problemCB2(x) and the node overflow function paper1_problemYL2(x) in steps b and c, and compared with the Pareto curve before the evolutionary algorithm is adjusted. If the adjusted solution is closer to the origin of the coordinate, it is proved to be better than the solution before the evolutionary algorithm is adjusted. Finally, it stops when the increased area will make the LID layout scale exceed the total area of ​​the region, and then the Pareto solution result after optimization by the NSGA-II algorithm can be obtained.

[0058] Therefore, compared with the existing technology, this method is based on the PlatEMO platform, uses the SOBOL method to screen out key sensitive parameters and use this as a basis to guide the platform's NSGA-Ⅱ algorithm, which increases the LID layout area in sensitive areas to improve the reduction effect of node overflow, and reduces the LID layout area in insensitive areas to reduce the cost of program construction; at the same time, according to this method, the algorithm calculation efficiency is improved by 20%. Compared with other SWMM models for the study of sponge city LID program construction, in terms of LID spatial layout, the present invention has higher calculation efficiency and makes the LID spatial layout more reasonable. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

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

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

[0064] 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.

[0065] 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

[0066] The present invention is further described in detail below in conjunction with specific implementation modes.

[0067] 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, a functional relationship between the node overflow and the construction cost of the LID scheme is constructed 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, and the significance of the parameters on the node overflow is examined; 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 plan.

[0068] In this way, based on the existing interactive method based on the SWMM model and the platformEMO platform described in the background technology, this method conducts a global sensitivity analysis on the layout scale parameters of the LID facilities, adjusts the layout area of ​​the LID facilities according to the size of the impact of the parameters on the node overflow, and obtains a better urban LID facility layout plan after comparative analysis. This improves the quality of the optimal spatial layout plan of the LID facilities under the control target and improves the efficiency of solving multi-objective optimization problems. The node overflow refers to the amount of water overflowing from the node when the node (such as the inspection well) of the urban drainage system exceeds its design flow during rainfall.

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

[0070] Step a: Obtain basic data related to modeling of the urban area under study, and establish a SWMM model of the city in the SWMM software; based on the SWMM software, divide the sub-catchment areas according to the urban pipeline network, buildings and street distribution, and divide the study area into x sub-catchments; complete the parameter calibration of the SWMM model; finally, successfully build a complete SWMM model of the city and export the inp file, as shown in Figure 1.

[0071] Step b: LID facilities are deployed in each sub-catchment area, and the number of LID facilities is y (further, the deployed LID facilities include but are not limited to three types of LID facilities such as permeable pavement (PP), biological retention pond (BRC) and rain garden (RG); the deployment area of ​​y types of LID facilities in x sub-catchments is used as the model independent variable, with 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 sub-catchments 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 is established with LID facility construction cost and node overflow as two objective functions, and a functional relationship of LID facility construction cost is established on the multi-objective optimization model, namely:

[0072] (1)

[0073] Where: f1 is the construction cost function of LID facilities in the study area, in yuan; i is the number of sub-catchment areas, i=1, 2, 3, ..., x; j is the type of LID facility, 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 (in implementation, the unit area construction cost of various types of LID facilities is obtained based on the relevant engineering data of the urban area construction);

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

[0075] 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 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, that is, complete the construction of the script function relationship of the construction cost;

[0076] Step c: Establishing the functional relationship of node overflow;

[0077] (2)

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

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

[0080] The specific process of encapsulating the node overflow function during implementation is:

[0081] 1) Start reading and rewriting data of inp file;

[0082] First, use the function function (declare the function name, input and output) to define the peak flow function paper1_problemYL2(mian_ji), the input is a 3×4 variable matrix (length(blqy)), and it is set 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;

[0083] 2) Determine the number of LID facilities (LID-number) and the number of lines occupied by LID in the inp file of the SWMM model; Use the contains function (query parameter position) to read the lia1 and lia2 parameters of the string "LID_USAGE" and "JUNCTIONS". 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 as y here;

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

[0085] Use the contains function (query parameter position) to read the lia3 and lia4 parameters of the string "SUBCATCHMENTS" and "SUBAREAS". In the SWMM model, the difference between lia3 and lia4 minus 4 is the number of subcatchment areas.

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

[0087] Use the cell function (create cell array) to create an empty cell array with n=0 for the area, and use the zeros function (create zero matrix) to initialize the area matrix; then establish a for loop with parameter isi to extract the area of ​​the subcatchment from 1 to zhsqy_number, and determine the area position of the inp file of the SWMM model (such as {lia33+3-1+isi} row (51:59) column) (same as above); use the strip function (delete the blanks 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. This is converted to a numerical value (the names of the sub-catchment areas are not necessarily continuous during implementation and need to be extracted for subsequent operations); when extracting the names (similar to the above area extraction steps), the cell function is used to create an empty cell array of n=0 for the area, and the zeros function is used to initialize the sub-catchment area name matrix; then a for loop with the parameter isi2 is established to extract the names of the sub-catchment areas from 1 to the number of zhsqy_number, and determine the name position of the SWMM model inp file (such as {lia33+3-1+isi2} row (1:3) column) (same as above); finally, the strip function is used to remove the spaces in the file and only retain the text;

[0088] 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 the cell array of the number of LIDs created by the cell function, and these three parts are spliced ​​into a new newline1 variable;

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

[0090] 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;

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

[0092] Use the reshape function (reconstruct array) to convert the area matrix (mian_ji) into a read 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 of the parameter isi3, from 1 to zhsqy_number the number of sub-catchments, and use the strrep function (specify the substring in the replacement 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 areas of the sub-catchments, 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)

[0093] 8) Run SWMM to perform calculations and extract the OUT file of the SWMM model; (When implementing, the SWMM interface can be called for calculations according to the patent developed by the previous 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 accept each line of data in the original file;

[0094] 9) Determine the overflow position of the sub-catchment area node in the OUT file report and calculate the cumulative node overflow; 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 position rows line by line from 1 to all the rows to be extracted, and determine that the node overflow value is the seventh group of data (this is certain because this column of the OUT file of SWMM is the result value of the node overflow); finally, use the str2double function to convert all extracted strings into numbers, and use the sum function to sum them, thus completing the construction of the script function relationship of the node overflow.

[0095] Step d: The above two steps develop corresponding programs to encapsulate the two functions of construction cost and node overflow, namely paper1_problemYL2(x) and paper1_problemCB2(x), completing the preparation work for the correlation between construction cost and node overflow, so that it meets the data interface form of the platEMO platform, and facilitates the subsequent use of the platform's NSGA-Ⅱ evolutionary algorithm to solve the solution cost and node overflow;

[0096] Based on the platEMO platform, (in order to successfully call the evolutionary algorithm of the platform, it is necessary to set the corresponding constraints, that is, the LID facility area cannot be larger than the total area of ​​the subcatchment area. If the sum of the set LID areas exceeds the subcatchment area, SWMM will not be able to simulate successfully, causing the connection between SWMM and platEMO to collapse. 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:

[0097] (3)

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

[0099] According to the steps for solving the custom 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 (creates a matrix with 1 as the value) is used to create a matrix with 1 row and z columns composed of 1 values; then the parameter lower is used to define (the lower limit of each variable) the lower limit of the z independent variables, and the zeros function (creates a matrix with 0 as the value) is used to create a matrix with 1 row and z columns composed of 0 values, that is, the minimum layout area of ​​each LID facility is 0; the parameter upper (the upper limit of each variable) is used to define the upper limit 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 length of the lower bound matrix and the upper bound matrix created must be the same as the matrix with a value of 1 defined by the parameter encoding); use the parameter objFcn (to minimize the objective function) to minimize the two problem functions, construction cost and node overflow; finally, use the parameter conFcn (constraint function) to satisfy the constraint only when and 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 deployed minus the total area of ​​S1, 3.1979, should be less than or equal to zero), and complete the processing of other sub-catchment areas according to the above steps. In this way, 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.

[0100] 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:

[0101] 1 Since there are z parameters in total, the parameter dimension D is set to z. For symmetrical sampling, the parameter dimension M=D*2 is expanded; the number of sampling points N=5000 is set, 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 area, named MAX. Since three types of LIDs are arranged in each catchment area, the MAX array is divided by 3 for uniform sampling.

[0102] 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 sample points from 1 to N, and obtain the samples of p one by one, named r, with a range of r=Min + r.* (Max-Min); finally, distribute the sampling results to the matrix R, that is, the N×M parameter matrix.

[0103] 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, and the other columns remain unchanged. This can be replaced D times, for a total of D AB matrices.

[0104] 4 The above steps give A, B, AB i (i=1, 2, ..., D) These matrices, a total of (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.

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

[0106] (1)

[0107] 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 ith independent variable X; xi is a variance operator. When other variables are fixed, input the independent variable X i 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.

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

[0109] Define V using the zeros function X With S i, whose 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 .

[0110] 6 Calculation of total effect index S Ti , the formula is:

[0111] (2)

[0112] 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 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.

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

[0114] Define E using the zeros function X With S Ti , whose format is (D, 1), establish a for loop function, from 1 to N, calculate E for 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 .

[0115] Step f: The SOBOL sensitivity index (S i and S Ti ) are sorted 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;

[0116] In the specific implementation, according to the rainstorm intensity formula (Xiushan County, Chongqing City), the design rainfall process line is generated based on the Chicago rain type, 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;

[0117] 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) will be 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 total area of ​​the sub-catchment may lead to uneven data distribution of a single LID facility, making the SOBOL index smaller, so only the parameters with SOBOL index>0.01 are displayed. Generally, the index is greater than 0.5 and is considered sensitive, and is far less than 0.1 and is insensitive). The results show that the layout of LID facilities in sub-catchments S20, S10, S22, S15, S25 and S24 is more sensitive, and the SOBOL index of the remaining parameters is almost 0, which has no effect on the output results. Among them, the layout of LID facilities in sub-catchment S20 has the most significant impact on node overflow, especially the layout 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.

[0118] 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 about 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. For the Pareto optimal frontier, each point contains a theoretically optimal LID facility spatial layout plan, and the decision maker selects the plan based on the actual situation. The optimal solution with a maximum node overflow of 412 m3 and a plan construction cost of 68.8 million yuan is selected as the initial plan, namely Plan A. The LID layout area of ​​S20 and S15 is adjusted according to the results of parameter sensitivity analysis.

[0119] The adjustment process is as follows: under the constraint conditions, firstly, the layout area of ​​LID facilities in the remaining sub-catchment areas that are not sensitive variables (such as sub-catchment areas S2, S4, S27, etc.) is reduced to 0, and then the layout area of ​​LID facilities in the sub-catchment areas that are sensitive variables (S20PP, S20BRC, S15BRC) is gradually increased by the rated value (usually each increase in the layout area is 0.1 hm 2 ), each time the layout area is increased, the adjusted node overflow value and cost value are calculated using the cost function paper1_problemCB2(x) and the node overflow function paper1_problemYL2(x) in steps b and c, and compared with the Pareto curve before the evolutionary algorithm is adjusted. If the adjusted solution is closer to the origin of the coordinate, it is proved to be better than the solution before the evolutionary algorithm is adjusted. Finally, it stops when the increased area will make the LID layout scale exceed the total area of ​​the region, and then the Pareto solution result after optimization by the NSGA-II algorithm can be obtained.

[0120] When implemented 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, and 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.

[0121] 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 Scheme 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, and the node overflow reduction rate is as high as 90.29%.

[0122] During implementation, the applicant further verified the effect of the scheme after parameter optimization, and further explored the calculation efficiency of the Pareto approximate front curve after sensitivity parameter optimization. The present invention continued to solve on the basis of 10,000 evaluations of the NSGA-II algorithm. It was found that the Pareto approximate front curve after sensitivity parameter optimization was almost the same as 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 greatly reduces the construction cost of LID facilities while meeting the overflow control. To further verify the effectiveness of this method, this study optimized the LID facility spatial layout plan again 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 find the best solution is continuous, especially on the basis of 12,000 evaluations, it can still continue to find a better solution for the spatial layout of LID facilities based on the results of sensitivity analysis, and continuously provide optimization directions for the spatial layout of LID facilities in the sub-catchment area.

[0123] Therefore, based on the PlatEMO platform, the present invention uses the SOBOL method to screen out key sensitive parameters and use this as a basis to guide the NSGA-Ⅱ algorithm of the platform. By increasing the LID layout area in sensitive areas, the node overflow reduction effect is improved, and the LID layout area in insensitive areas is reduced to reduce the construction cost of the scheme; at the same time, the algorithm calculation efficiency is improved by 20% according to this method. Compared with other SWMM models for the study of sponge city LID scheme construction, in terms of LID spatial layout, the present invention has higher calculation efficiency and makes the LID spatial layout more reasonable.

Claims

1. A method for optimizing the layout of urban LID facilities based on SOBOL sensitivity analysis, characterized in that: First, a SWMM model of the urban area was established for research. The 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 on the layout scale parameters of LID facilities in each sub-catchment area, and the significance of the influence of parameters on node overflow was examined. Then, the NSGA-II algorithm in the MATLAB software platEMO4.0 platform was used to find the Pareto optimal solution set of the spatial layout of LID facilities, and the LID facility layout area of ​​the selected sensitive variables was adjusted and optimized according to the results of parameter sensitivity analysis. The Pareto optimal solution set before and after parameter optimization was compared and analyzed to obtain a better urban LID facility layout plan.

2. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 1, characterized in that: The method comprises the following steps: Step a: Obtain the basic data related to modeling of the urban area under study, and establish the SWMM model of the city in the SWMM software; Based on the SWMM software, divide the sub-catchment area according to the urban pipe network direction, buildings and street distribution, and divide the study area into x sub-catchment areas; complete the parameter calibration of the SWMM model; Finally, successfully build a complete SWMM model of the city and export the inp file; Step b: LID facilities are laid out in each sub-catchment area, and the number of LID facilities is y; the layout area of ​​y types of LID facilities in x sub-catchment areas is used as the model independent variable, with a total of z independent variables, z = x y; based on MATLAB software, a multi-objective optimization model is established with LID facility construction cost and node overflow as two objective functions, and a functional relationship of LID facility construction cost is established on the multi-objective optimization model, that is: （1） Where: f1 is the construction cost function of LID facilities in the study area, in yuan; i is the number of sub-catchment areas, i=1, 2, 3, ..., x; j is the type of LID facility, 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 (in implementation, the unit area construction cost of various types of LID facilities is obtained based on the relevant engineering data of the urban area construction); Use matlab's function command to encapsulate the construction cost function; Step c: Establishing the functional relationship of node overflow; （2） Where: f2 is the node overflow, unit is m 3 / s; g is the SWMM model calculation and result extraction function, obtained by Matlab programming; Use matlab's function command to encapsulate the node overflow function; Step d: The above two steps develop corresponding programs to encapsulate the two functions of construction cost and node overflow, namely paper1_problemYL2(x) and paper1_problemCB2(x), completing the preparation work for the correlation between construction cost and node overflow, so that it meets the data interface form of the platEMO platform, and facilitates 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: （3） Where: is the area of ​​the ith subcatchment; According to the steps for solving the custom problem of the platEMO platform, firstly, the parameter encoding is used to define the y types of LID facilities on the x sub-catchment areas, that is, the encoding method of a total of z variables, and the ones function is used to create a matrix with 1 row and z columns composed of 1 values; then the parameter lower 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 composed of 0 values, that is, the minimum layout area of ​​each LID facility is 0; the parameter upper is used to define the upper bound of each independent variable; the parameter objFcn is used to minimize the two problem functions of construction cost and node overflow; finally, the parameter conFcn is used, if and only if the constraint violation value is less than or equal to zero, the constraint is satisfied, and the processing of other sub-catchment areas 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 conduct 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=D*2 is expanded. The number of sampling points N=5000 is set, 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 area, named MAX. Since three types of LIDs are arranged in each catchment area, the MAX array is divided by 3 for uniform sampling.

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 sample points from 1 to N, and obtain the sample points of p one by one, named r, whose range is r=Min + r.* (Max-Min); finally, distribute the sampling results to 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 to replace the i-th column in matrix B with the i-th column in matrix A from 1 to D, and the other columns remain unchanged, then it can be replaced D times, a total of D AB matrices; 4 The above steps give A, B, AB i (i=1, 2, ..., D) These matrices, a total of (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: （1） 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 ith independent variable X; xi is a variance operator. When other variables are fixed, input the independent variable X i 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 independent variables; The specific calculation process is obtained through Matlab programming, as follows: Define V using the zeros function X With S i , whose 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 ; 6 Calculation of total effect index S Ti , the formula is: （2） 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 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; The specific calculation process is obtained through Matlab programming, as follows: Define E using the zeros function X With S Ti , whose format is (D, 1), establish a for loop function, from 1 to N, calculate E for 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-catchment areas 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 area is judged 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 NSGA-Ⅱ evolutionary algorithm of the platform 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 sub-catchment area is adjusted according to the sensitivity; The adjustment process is as follows: under the constraints, firstly, the layout area of ​​the remaining LID facilities in the sub-catchment areas that are not sensitive variables is reduced to 0, and then the layout area of ​​the LID facilities in the sub-catchment areas that are sensitive variables is gradually increased by the rated value. Each time the layout area is increased, the adjusted node overflow value and cost value are calculated using the cost function paper1_problemCB2(x) and the node overflow function paper1_problemYL2(x) in steps b and c, and compared with the Pareto curve before the evolutionary algorithm is adjusted. If the adjusted scheme is closer to the origin of the coordinate, it is proved to be better than the scheme before the evolutionary algorithm is adjusted. Finally, it is stopped until the increased area makes the LID layout scale exceed the total area of ​​the area, and then the Pareto solution result after optimization by the NSGA-II algorithm can be obtained.

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

4. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 2, characterized in that: In step b, the process of encapsulating the construction cost function using the function command of MATLAB 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.

5. The urban LID facility layout optimization method based on SOBOL sensitivity analysis according to claim 2, characterized in that: In step c, the process of encapsulating the node overflow function using the MATLAB function command is: 1) Start reading and rewriting data of inp file; First, the peak flow function paper1_problemYL2(mian_ji) is defined 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). The inp file of the SWMM model is opened in a readable and writable manner. The fgetl function is used 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 inp file of the SWMM model; Use the contains function to read the lia1 and lia2 parameters of the string "LID_USAGE" and "JUNCTIONS" rows. 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 as y here; 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 of the string "SUBCATCHMENTS" and "SUBAREAS". In the SWMM model, the difference between lia3 and lia4 minus 4 is the number of subcatchment areas. 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 sub-catchment area from 1 to zhsqy_number, and determine the area position of the inp file of the SWMM model; the strip function is used to remove spaces in the file and only retain text; finally, the str2double function is used to convert the text into a numerical 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 sub-catchment area name matrix; then a for loop with the parameter isi2 is established to extract the name of the sub-catchment area from 1 to zhsqy_number, and determine the name position of the inp file of the SWMM model; finally, the strip function is used to remove spaces in the file and only retain 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 the cell array of the number of LIDs created by the cell function, and 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 area and set constraints; Use the reshape function to convert the area matrix mian_ji into a read format of y rows and x columns in the order of columns, and use the 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 parameter isi3, from 1 to zhsqy_number, the number of sub-catchments, 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 sub-catchment area. If it exceeds the sum of the sub-catchment areas, it is set to 0 and no LID is set; 8) Run SWMM to perform calculations 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 in the original file; 9) Determine the overflow position of the sub-catchment area 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 position rows line by line from 1 to all the rows to be extracted, and determine the node overflow volume value as the seventh group 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 construction of the script function relationship of the node overflow volume.

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

5.

7. The urban LID facility layout optimization method based on SOBOL sensitivity analysis as claimed in claim 2, characterized in that: In step g, the rated value of each increase in layout area is 0.1 hm2.

Images