Automatic parameter calibration method for sewer network model and storage medium

Abstract

The application discloses a kind of automatic calibration method of sewer network model parameter and storage medium, method includes: the hydraulic model of sewer network is constructed;According to the value range of each model parameter of the hydraulic model respectively, the initial value of each model parameter is randomly generated, and the optimal value of each model parameter is determined by simulated annealing algorithm, the model parameter includes impermeable area Manning roughness coefficient, permeable area Manning roughness coefficient, impermeable area depression storage depth, permeable area depression storage depth, maximum infiltration rate, minimum infiltration rate, infiltration attenuation coefficient and drain time.The uncertainty parameter in pipe network model can be quickly and accurately calibrated by the present application.

Description

Technical Field

[0001] This invention relates to the field of drainage engineering technology, and in particular to an automatic calibration method and storage medium for drainage network model parameters. Background Technology

[0002] With the acceleration of urbanization and the continuous growth of urban population, the sustainable management and operation of urban drainage systems has become a significant challenge. One of the core components of an urban drainage system is the drainage pipe network, which is responsible for collecting and discharging rainwater and wastewater. To ensure the efficient operation of drainage systems and to address the ever-changing environmental and population demands, accurate drainage pipe network models have become crucial tools for urban planning, design, and management. Drainage pipe network models are computer models built based on mathematical and physical principles, capable of simulating the flow, transmission, and discharge of rainwater and wastewater within the network. These models help urban planners and engineers better understand the operational status of the pipe network system, predict potential problems, and develop corresponding solutions. However, due to the complexity and numerous parameters of drainage pipe network models, manually calibrating these parameters becomes extremely tedious and time-consuming, limiting the application of these models.

[0003] For example, the SWMM (Storm Water Management Model) is a computer model widely used in urban drainage system planning and management. It can simulate the flow, transmission, and discharge of rainwater and wastewater in drainage networks, and predict flood events and water quality changes. The application of SWMM models covers many aspects of urban drainage systems, including rainwater runoff, sewage flow, pump station operation, and water quality control. However, SWMM contains numerous parameters. Given the increasingly high accuracy requirements of online models, manually calibrating the parameters of large-scale urban models is impossible. How to set reasonable parameters to achieve high-precision simulation is a challenge that needs to be overcome.

[0004] With the continuous advancement of computer technology, the application of optimization algorithms and other technologies has provided strong support for the realization of automatic calibration technology. Computers can search the parameter space in a short time to find the optimal solution, thereby achieving automatic model calibration. Against this backdrop, automatic calibration technology for drainage pipe networks has emerged, becoming an important trend and innovation in the development of drainage pipe network models. Automatic calibration technology aims to automatically adjust model parameters through computer algorithms and optimization methods, making the simulation results of the model closer to actual observation data. This technology can perform model calibration at finer time and spatial scales, enabling more accurate predictions of the operation of drainage systems under different scenarios. This is of great significance for urban drainage safety and the optimization of water company pipe network management.

[0005] A Chinese patent document with publication number CN116595747A discloses a method for calibrating a Swing Runoff Model (SWMM) based on a particle swarm optimization (PSO) algorithm. This method involves generating different category zones using categorized sub-catchments; calculating the rainfall duration of the study area and the dynamic runoff coefficient of each category zone based on historical rainfall data; and optimizing the uncertain parameters of the SWMM model based on theoretical rainfall flow. This method optimizes the uncertain parameters of different sub-catchments using the PSO algorithm, reducing human subjectivity. However, because the underlying online model cannot be changed frequently in reality, this method relies on the PSO algorithm to solve for the dynamic runoff coefficient, which may increase workload and maintenance costs.

[0006] Chinese patent document CN115689051A discloses an automatic parameter calibration method for SWMM models based on the GA algorithm coupled with Morris and GLUE. It mentions first using the Global Uncertainty Analysis (GLUE) method to perform numerous repeated model simulations to narrow the parameter value range; then, using a genetic algorithm with single-objective and multi-objective constraints for parameter optimization, thus calibrating the SWMM model parameters. The method's initial use of numerous repeated model simulations to narrow the parameter value range is essentially aimed at minimizing the space for rationalized parameter values. However, extensive repeated calculations are very time-consuming at the large-scale model level. Directly adopting the suggested values ​​from the official documentation can significantly reduce computation time at this level. Furthermore, the GA algorithm, lacking an elitist preservation strategy, is prone to losing optimal solutions. Summary of the Invention

[0007] The technical problem to be solved by the present invention is to provide an automatic calibration method and storage medium for drainage network model parameters, which can quickly and accurately calibrate the uncertain parameters in the network model.

[0008] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is: an automatic calibration method for drainage pipe network model parameters, comprising:

[0009] Construct a hydraulic model of the drainage pipe network;

[0010] Initial values ​​for each model parameter are randomly generated based on the range of values ​​for each model parameter in the hydraulic model. Optimal values ​​for each model parameter are then determined using a simulated annealing algorithm. The model parameters include the Manning roughness coefficient for the impermeable zone, the Manning roughness coefficient for the permeable zone, the water storage depth in the depression of the impermeable zone, the water storage depth in the depression of the permeable zone, the maximum infiltration rate, the minimum infiltration rate, the permeability attenuation coefficient, and the drainage time.

[0011] The present invention also proposes a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method described above.

[0012] The beneficial effects of this invention are as follows: by establishing a hydraulic model of the pipeline network and automatically finding the optimal solution of parameters through simulated annealing algorithm, the problems of strong subjectivity and low calibration accuracy of traditional manual calibration methods are overcome. It can quickly and accurately calibrate the uncertain parameters in the pipeline network model, thereby reducing the gap between the model simulation value and the actual monitoring value. This not only greatly improves the accuracy of the model, but also significantly reduces the time and manpower cost of model calibration. Attached Figure Description

[0013] Figure 1 This is a flowchart of an automatic calibration method for drainage network model parameters according to the present invention;

[0014] Figure 2 This is a flowchart of the method according to Embodiment 1 of the present invention;

[0015] Figure 3 This is a negative value trajectory diagram of the objective function of each generation of the simulated annealing algorithm under the parallel framework in Embodiment 1 of the present invention. Detailed Implementation

[0016] To explain the technical content, objectives, and effects of the present invention in detail, the following description is provided in conjunction with the embodiments and accompanying drawings.

[0017] Please see Figure 1 An automatic calibration method for drainage network model parameters, comprising:

[0018] Construct a hydraulic model of the drainage pipe network;

[0019] Initial values ​​for each model parameter are randomly generated based on the range of values ​​for each model parameter in the hydraulic model. Optimal values ​​for each model parameter are then determined using a simulated annealing algorithm. The model parameters include the Manning roughness coefficient for the impermeable zone, the Manning roughness coefficient for the permeable zone, the water storage depth in the depression of the impermeable zone, the water storage depth in the depression of the permeable zone, the maximum infiltration rate, the minimum infiltration rate, the permeability attenuation coefficient, and the drainage time.

[0020] As can be seen from the above description, the beneficial effect of the present invention is that it can quickly and accurately calibrate the uncertainty parameters in the pipeline network model, thereby reducing the gap between the model simulation value and the actual monitoring value.

[0021] Furthermore, the initial values ​​of each model parameter are randomly generated based on the value range of each model parameter of the hydraulic model, and the optimal values ​​of each model parameter are determined through simulated annealing algorithm, specifically as follows:

[0022] Obtain the actual water depth at each liquid level monitoring point at multiple preset times;

[0023] Based on the value range of each model parameter of the hydraulic model, a set of initial values ​​for each model parameter are randomly generated as the current solution;

[0024] Based on the current solution, the hydraulic model is modified to obtain the hydraulic model corresponding to the current solution. Then, using the hydraulic model corresponding to the current solution, the simulated water depth of each liquid level monitoring point at multiple preset times is calculated.

[0025] The fitness of the current solution is calculated based on the simulated water depth at each liquid level monitoring point at multiple preset times and the actual water depth at each liquid level monitoring point at multiple preset times.

[0026] A new solution is generated by randomly perturbing the current solution;

[0027] Based on the new solution, the hydraulic model is modified to obtain the hydraulic model corresponding to the new solution. Then, using the hydraulic model corresponding to the new solution, the simulated water depth of each liquid level monitoring point at multiple preset times is calculated.

[0028] The fitness of the new solution is calculated based on the simulated water depth at each liquid level monitoring point at multiple preset times and the actual water depth at each liquid level monitoring point at multiple preset times.

[0029] If the fitness of the new solution is greater than that of the current solution, then according to the Metropolis acceptance criterion, the new solution is adopted as the current solution or the current solution is retained; otherwise, the new solution is adopted as the current solution.

[0030] If the preset iteration termination condition is met, the current solution is taken as the optimal solution; otherwise, return to the step of generating a new solution by randomly perturbing the current solution.

[0031] Based on the optimal solution, determine the optimal values ​​of each model parameter.

[0032] As described above, the actual water depths measured at each liquid level monitoring point are used as validation data to optimize the specified model parameters. Based on the differences between the measured data and the model simulation results, the model parameters are adjusted to make the simulation results more accurate.

[0033] Further, the step of calculating the fitness of the current solution based on the simulated water depth of each liquid level monitoring point at multiple preset times and the actual water depth of each liquid level monitoring point at multiple preset times is as follows:

[0034] According to the Nash coefficient calculation formula, the composite Nash coefficient corresponding to the current solution is calculated, and the negative of the composite Nash coefficient corresponding to the current solution is taken as the fitness of the current solution. The Nash coefficient calculation formula is as follows:

[0035]

[0036] Where NSE represents the composite Nash coefficient corresponding to the current solution, L obs,i,j L represents the actual water depth at the i-th liquid level monitoring point at the j-th time. pred,i,j L represents the simulated water depth at the i-th liquid level monitoring point at the j-th time corresponding to the current solution. obs,i,mean This represents the average actual water depth of the i-th liquid level monitoring point at n preset times, where num is the number of liquid level monitoring points and n is the total number of preset times.

[0037] As described above, the objective function is the negative of the composite Nash coefficient calculated based on the simulated and actual values ​​of the monitoring points.

[0038] Furthermore, the preset iteration termination conditions include at least one of the following: reaching a preset maximum number of iterations, the current temperature being less than a preset termination temperature, and the difference between the fitness of the new solution and the fitness of the current solution in consecutive preset iteration rounds being less than a preset threshold.

[0039] Among them, the current temperature T t =a t ×T t-1 a t = a0×exp(t / 50), a t The current cooling rate is given by , a0 is the preset initial cooling rate, t is the current iteration number, and T is the current iteration number. t-1 T0 is the temperature at the (t-1)th iteration, and T0 is the preset initial temperature.

[0040] As can be seen from the above description, by setting appropriate iteration termination conditions, premature entry into local optima can be avoided, ensuring convergence to the global optimum.

[0041] Furthermore, the initial values ​​of each model parameter are randomly generated based on the value range of each model parameter of the hydraulic model, and the optimal values ​​of each model parameter are determined through simulated annealing algorithm, specifically as follows:

[0042] Based on the value range of each model parameter of the hydraulic model, at least one set of initial values ​​for each model parameter is randomly generated to obtain at least one set of initial solutions.

[0043] Simulated annealing search is performed on at least one set of initial solutions in parallel to obtain at least one set of optimal solutions;

[0044] Based on the optimal solutions of at least one set, determine the optimal values ​​of each model parameter.

[0045] As can be seen from the above description, parallel computing can overcome the shortcomings of excessive time consumption in the calibration of models in general serial computing, and can explore from different starting points at the same time, increasing the chance of finding the global optimal solution, thereby realizing the parameter calibration of sub-catchment areas in large-scale real-world pipe networks.

[0046] Furthermore, before the step of randomly generating initial values ​​for each model parameter based on the value range of each model parameter of the hydraulic model, and determining the optimal values ​​for each model parameter through simulated annealing algorithm, the method further includes:

[0047] The integrity and volatility of the actual values ​​at each liquid level monitoring point are assessed.

[0048] If both the integrity assessment and the volatility assessment pass, then the steps of randomly generating initial values ​​for each model parameter based on the range of values ​​for each model parameter of the hydraulic model, and determining the optimal values ​​for each model parameter through a simulated annealing algorithm are executed.

[0049] Furthermore, the integrity assessment and fluctuation assessment of the actual values ​​of the liquid level monitoring points specifically includes:

[0050] If the ratio of the total number of actual values ​​of each liquid level monitoring point within the preset time period to the preset monitoring data volume is greater than or equal to the preset ratio, the integrity assessment is deemed to have passed.

[0051] If a predetermined number of consecutive actual values ​​are the same at a liquid level monitoring point, the fluctuation assessment of that liquid level monitoring point is deemed unsuccessful.

[0052] As can be seen from the above description, by evaluating the actual monitoring quantity and ensuring the integrity and volatility of the actual monitoring data, the accuracy of subsequent parameter calibration can be improved.

[0053] The present invention also proposes a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method described above.

[0054] Example 1

[0055] Please refer to Figure 2-3 Embodiment 1 of the present invention is: an automatic calibration method for drainage pipe network model parameters, which can be applied to urban drainage pipe network systems, such as... Figure 2 As shown, it includes the following steps:

[0056] S1: Construct a hydraulic model of the drainage network.

[0057] The hydraulic model constructed in this embodiment is an SWMM model. Specifically, basic data of the study area, namely regional GIS information, such as underlying surface data, rainfall data, and underground pipeline network structure, are collected. This data forms the basis for building the model. The underlying surface data helps determine the impermeability ratio of the catchment area, while the rainfall data simulates the response of the drainage system under different rainfall conditions. Then, the model framework is built, reasonable catchment areas are set, and the hydraulic model of the drainage network of the study area is constructed using Stormdesk modeling software. This model is then exported as an .inp file.

[0058] S2: Perform integrity and volatility assessments on the actual values ​​of each liquid level monitoring point, and determine whether the assessment is successful. If so, proceed to step S3.

[0059] This step involves evaluating the actual monitoring data, primarily including integrity assessment and volatility assessment. Specifically, if the ratio of the total number of actual values ​​at each liquid level monitoring point within a preset time period (e.g., one day) to the preset amount of monitoring data is greater than or equal to a preset ratio (e.g., 0.8), the integrity assessment is deemed successful. If a liquid level monitoring point has a preset number of consecutive identical actual values ​​(e.g., 10), the volatility assessment for that liquid level monitoring point is deemed unsuccessful.

[0060] In this embodiment, the integrity assessment involves calculating the amount of monitoring data that should be acquired in a day based on the data acquisition step size, then dividing the amount of monitoring data acquired on the current day by the amount of monitoring data that should be acquired in a day. If the quotient is above 80%, the integrity assessment is considered passed. Simultaneously, if data remains the same for an extended period, it may indicate instrument malfunction. Therefore, at least one of the ten consecutive monitoring data points must be different; otherwise, the data is considered abnormal. Further, after the above assessment, the data deemed normal is exported as a CSV file. The three column headers are "ID," "Time," and "Observed," recording the monitoring point ID, time series, and corresponding monitoring data.

[0061] S3: Based on the value range of each model parameter of the hydraulic model, randomly generate the initial value of each model parameter, and determine the optimal value of each model parameter through simulated annealing algorithm.

[0062] This step uses the actual water depth measured at each liquid level monitoring point as validation data to calibrate eight parameters of the model's catchment area: the Manning roughness coefficient (N-Imperv) for the impermeable zone, the Manning roughness coefficient (N-Perv) for the permeable zone, the water depth in the impermeable zone depression (S-Imperv) / mm, the water depth in the permeable zone depression (S-Perv) / mm, the maximum infiltration rate (MaxRate) / (mm·h⁻¹), the minimum infiltration rate (MinRate) / (mm·h⁻¹), the desiccation coefficient (Decay), and the drainage time (Drytime) / d. Using Python and the SWMM_API library, the eight model parameters are used as optimization variables. The SWMM_API library is used to implement the `inp` function. The combined Nash coefficient calculated from the simulated and monitored values ​​at each liquid level monitoring point is multiplied by -1 to form the objective function. An adaptive simulated annealing algorithm is used as the optimization framework to find the minimum value of this objective function.

[0063] Specifically, the process of finding the optimal solution using the simulated annealing algorithm includes the following steps:

[0064] S301: Obtain the actual water depth at each liquid level monitoring point at multiple preset times.

[0065] S302: Based on the value range of each model parameter of the hydraulic model, randomly generate a set of initial values ​​for each model parameter as the current solution.

[0066] In this embodiment, the value range of each model parameter is shown in Table 1.

[0067] Table 1: Range of Model Parameter Values

[0068] 1 Manning roughness coefficient of impermeable area 0.01-0.025 2 Manning roughness coefficient of permeable zone 0.05-0.440 3 Water storage depth in impermeable depressions 0.2-10 4 Water storage depth in permeable depressions 2-20 5 Maximum infiltration rate 20-100 6 Minimum infiltration rate 1-19 7 Permeation attenuation coefficient 2-7 8 Draining time 2-14

[0069] S303: Based on the current solution, modify the hydraulic model to obtain the hydraulic model corresponding to the current solution, and use the hydraulic model corresponding to the current solution to calculate the simulated water depth of each liquid level monitoring point at multiple preset times.

[0070] That is, based on the values ​​of each model parameter in the current solution, the corresponding parameters in the hydraulic model are modified, and then the SWMM calculation engine is called to calculate the modified hydraulic model and read the simulated water depth of each liquid level monitoring point in the drainage network.

[0071] In practical applications, the SWMM_API library in Python can be used. The SwmmInput function reads the .inp file, and inp.SUBAREAS is used to read the Manning roughness coefficient of the impermeable zone, the Manning roughness coefficient of the permeable zone, the water depth of the impermeable depression, and the water depth of the permeable depression in the catchment area. Inp.INFILTRATION is used to read the maximum infiltration rate, minimum infiltration rate, infiltration attenuation coefficient, and drainage time. The values ​​of the model parameters in the current solution are assigned to the parameters read above, the .inp file is rewritten, and a new .inp file is generated. The swmm5_run function is then called to calculate the .inp file. Finally, the SwmmOut function reads the .out file and extracts the water depth time series simulated from the ID nodes of each monitoring point.

[0072] S304: Calculate the fitness of the current solution based on the simulated water depth of each liquid level monitoring point at multiple preset times and the actual water depth of each liquid level monitoring point at multiple preset times.

[0073] Specifically, according to the Nash coefficient calculation formula, the composite Nash coefficient corresponding to the current solution is calculated, and the negative of the composite Nash coefficient corresponding to the current solution is taken as the fitness of the current solution, that is, the fitness of the current solution = -1 × NSE. In this embodiment, the Nash coefficient calculation formula is:

[0074]

[0075] Where NSE represents the composite Nash coefficient corresponding to the current solution, L obs,i,j L represents the actual water depth at the i-th liquid level monitoring point at the j-th time. pred,i,j L represents the simulated water depth at the i-th liquid level monitoring point at the j-th time corresponding to the current solution. obs,i,mean This represents the average actual water depth of the i-th liquid level monitoring point at n preset times, where num is the number of liquid level monitoring points and n is the total number of preset times.

[0076] S305: Generate a new solution by randomly perturbing the current solution.

[0077] This means that a new state or solution is randomly generated based on the current solution set state. In this embodiment, this new state can be obtained by making a small perturbation within the neighborhood of the current state. A random perturbation method is used, where a random number is generated using the `random()` function, and the random number is multiplied by a perturbation range to achieve the perturbation. First, the smaller of the difference between the upper and lower limits of the variable's interval and the current value of the variable is taken. This ensures that the newly generated variable does not exceed a reasonable range. Simultaneously, the perturbation range is related to temperature; as the temperature decreases, the perturbation range decreases. Therefore, the final value of the perturbation range = random number × min(difference between the current parameter value and the upper limit of its range, difference between the current parameter value and the lower limit of its range) × (1 - 1 / T), where T is the current temperature. Since this is the beginning of the current iteration, it is the temperature T calculated in the previous iteration. t-1 If this iteration is the first iteration, then it is the initial temperature T0.

[0078] S306: Based on the new solution, modify the hydraulic model to obtain the hydraulic model corresponding to the new solution, and use the hydraulic model corresponding to the new solution to calculate the simulated water depth at each liquid level monitoring point at multiple preset times. This step can refer to step S303 above.

[0079] S307: Calculate the fitness NSE' corresponding to the new solution based on the simulated water depth at each liquid level monitoring point at multiple preset times and the actual water depth at each liquid level monitoring point at multiple preset times. This step can be referred to step S304 above.

[0080] S308: Determine whether the fitness of the new solution is greater than that of the current solution, i.e., whether NSE'-NSE>0 is true. If yes, proceed to step S309; ​​otherwise, proceed to step S310.

[0081] This step involves determining the change in the energy difference (i.e., the objective function) between the new solution and the current solution. The energy difference can be used to measure the quality of the new state relative to the current state.

[0082] S309: According to the Metropolis acceptance criterion, either adopt the new solution as the current solution or retain the current solution. That is, use the Metropolis criterion to calculate the acceptance probability and use this probability to accept the new solution as the current solution. Then proceed to step S311.

[0083] S310: Directly use the new solution as the current solution. Then execute step S311.

[0084] S311: The current temperature is calculated. In each iteration, the temperature is decreased according to the cooling rate to simulate the physical annealing process. The decrease in temperature gradually reduces the probability of accepting a poor solution, and the algorithm tends to converge to the optimal solution.

[0085] Specifically, before starting the algorithm, parameters such as initial temperature, cooling rate, and iteration termination condition are set. The initial temperature controls the randomness of the search, the cooling rate determines the degree of temperature decrease over time, and the iteration termination condition can be reaching a certain number of iterations or the temperature decreasing to a certain level. In high-dimensional problems, a higher initial temperature and a lower cooling rate should be set to prevent the algorithm from prematurely entering a local optimum. During the algorithm's operation, a small initial cooling rate is used initially, and then adjusted according to the algorithm's performance. A smaller cooling rate results in a slower temperature decrease, which is beneficial for more exploration. In the later stages, when entering the stable period, a larger cooling rate should be used to terminate the computation as soon as possible.

[0086] Therefore, in this embodiment, the formula for calculating the cooling rate is set as a. t = a0×exp(t / 50), where a t Let a0 be the current cooling rate, which is the cooling rate at the t-th iteration, and t be the current iteration number. In this embodiment, a0 = 1.

[0087] Therefore, the formula for temperature calculation is T. t =a t ×T t-1 , among which, T t Let T be the temperature at the t-th iteration. t-1 T0 is the temperature at the (t-1)th iteration, and T0 is the preset initial temperature.

[0088] S312: Determine whether the preset iteration termination condition has been met. If yes, execute step S313. If no, repeat steps S305-311 to directly satisfy the iteration termination condition.

[0089] The preset iteration termination conditions include at least one of the following: reaching a preset maximum number of iterations, the current temperature being less than a preset termination temperature, and the difference between the fitness of the new solution and the fitness of the current solution in consecutive preset iteration rounds being less than a preset threshold.

[0090] In high-dimensional problems, in order to converge to the global optimal solution, the maximum number of iterations needs to be set to a larger value. To prevent invalid calculations in the later stages, a termination temperature T2 is set. At the same time, the calculation is also terminated when there is no significant change in the solution after 5 consecutive iterations. The threshold is set to 0.001, that is, the iteration is terminated when the change in the calculated Nash coefficient is less than 0.001 after 5 consecutive iterations.

[0091] S313: Take the current solution as the optimal solution, and determine the optimal values ​​of each model parameter based on the optimal solution. That is, after the iteration termination condition is met, the algorithm will converge to a solution, which can be output as the final result.

[0092] The above steps describe the process of finding a set of optimal solutions using the simulated annealing algorithm. Furthermore, a parallel computing framework can be built using the multiprocessing library, and the PoolType of Process can be used to implement multi-process computing. That is, by performing simulated annealing searches on different initial solutions in parallel by multiple processes, multiple sets of optimal solutions can be obtained. Finally, the final optimal solution is determined from the multiple sets of optimal solutions.

[0093] In practical applications, due to the limitations of the GIL (Global Interpreter Lock), Python multithreading can only run on a single-core CPU. To achieve parallelism on multiple core CPUs, it must be implemented using multiple processes. Therefore, the `multiprocessing` module can be used to create multiple processes, and multiple child processes can be created using a process pool. Example code: `multiprocessing.Pool(processes = None, initializer = None, initargs = (), maxtasksperchild = None)`. In this embodiment, the following function structure is used:

[0094] forn inVars:

[0095] result=pool.apply_async(subAimFunc,(n,))

[0096] results.append(result)

[0097] The code snippet `Vars` records the objective function values ​​returned by the simulated annealing algorithm in all processes. It iterates through all these processes, retrieving the objective function values ​​and then writing the code to calculate and return the objective function value in `subAimFunc`. The `results` variable records the objective function calculated by the simulated annealing algorithm in each generation of processes. A second loop then sorts the `results` variable according to the order of the simulated annealing processes to prevent the results from becoming unordered, thus ensuring successful parallel processing. The code is as follows:

[0098] Forres in results:

[0099] obresult.append(res.get())

[0100] Objv = np.array(obresult)

[0101] After the above operations, a parallel computing structure can be implemented. Simultaneously, the optimal solution and mean value of each objective function in each generation are recorded and saved as a trajectory graph. For example... Figure 3 As shown, Figure 3 The upper curve represents the negative value of the optimal objective function in each generation of the parallel process (i.e., the combined Nash coefficient), while the lower curve represents the negative value of the average objective function in each generation of each process. It can be seen that the combined Nash coefficient increased from 0.87 at the beginning to 0.93, and the curve fitting became more accurate.

[0102] Furthermore, in this embodiment, the initial temperature T0 is set to 500, the initial cooling rate a0 is 1, and the termination temperature T2 is set to 10. This ensures a large search space in the early stages of the algorithm, while preventing the full utilization of parallel computing advantages if the solution space is insufficient. The number of process cores for parallel computing is set to the number of computer cores minus 1. Due to the use of a parallel computing framework, computing resources can be utilized as much as possible. Reasonable parameter settings achieve an effective balance between solution accuracy and solution time. Simultaneously, a create-inp function is constructed. By inputting the optimal solution, its parameters are overwritten into the original inp function to construct the optimal solution inp. Ultimately, optimization algorithms can be performed within the parallel computing framework, generating the corresponding inp file from the optimal solution. The entire optimization process is recorded to facilitate the identification of errors during optimization.

[0103] This embodiment utilizes the impact of hydrological parameters in the catchment area on each monitoring point, using the node water depth information recorded by SCADA monitoring points in the model as calibration data. The SWMM_API is used to modify the INP data and read the corresponding SECTION data, enabling batch data modification and result data recording. The objective function is the negative of the composite Nash coefficient calculated from the simulated and measured values ​​of each monitoring point. Based on this, a simulated annealing algorithm is introduced to calibrate the parameters. The simulated annealing algorithm, through its random search method, can escape local optima, thus finding a broader global optimum in the search space. This makes it easier to escape local optima when dealing with high-dimensional problems. Meanwhile, a simulated annealing algorithm under a parallel computing framework is constructed using the multiprocessing parallel library. At the start of the simulated annealing algorithm, multiple initial solutions can be generated simultaneously and assigned to different parallel processes or processing units. Each process can independently execute the simulated annealing search, starting from different starting points, increasing the chance of finding the global optimum and improving the computation speed. This enables parameter calibration of sub-catchment areas in large-scale real-world pipe networks, maximizing the accuracy of the underlying model of the online model, and ensuring that the online model reflects reality without distortion. This provides a theoretical basis and technical support for ensuring urban drainage safety.

[0104] This embodiment, through model calibration, enables more accurate predictions of drainage system operations, including stormwater runoff, flood events, and sewage flow. This is crucial for urban flood management, drainage system planning, and optimization. The method in this embodiment helps water companies precisely adjust model parameters based on extensive measured data, making simulation results closer to reality and improving model accuracy. Furthermore, the method in this embodiment is written based on the SWMM-API, which can automatically read values ​​from the corresponding regions of the SWMM model and calculate the objective function, avoiding the need for targeted adjustments in traditional methods. It can be repeatedly applied in different scenarios, ensuring the model's replicability and scalability.

[0105] Example 2

[0106] This embodiment is a computer-readable storage medium corresponding to the above embodiments, on which a computer program is stored. When the program is executed by a processor, it implements the various steps of the automatic calibration method for drainage network model parameters as described in the above embodiments, and can achieve the same technical effect, which will not be repeated here.

[0107] In summary, the present invention provides an automatic calibration method and storage medium for drainage network model parameters. By establishing a model, using an adaptive simulated annealing algorithm, and employing a parallel computing structure, it overcomes the problems of strong subjectivity and low calibration accuracy of traditional manual calibration methods. At the same time, it overcomes the drawback of excessive time consumption when calibrating the model in general serial calculations. It can quickly and accurately calibrate the uncertain parameters of sub-catchment areas in the network model, thereby greatly reducing the gap between the model simulation values ​​and the actual monitoring values.

[0108] The above description is merely an embodiment of the present invention and does not limit the patent scope of the present invention. Any equivalent modifications made based on the content of the present invention specification and drawings, or direct or indirect applications in related technical fields, are similarly included within the patent protection scope of the present invention.

Claims

1. An automatic calibration method for drainage pipe network model parameters, characterized in that, include: Construct a hydraulic model of the drainage pipe network; Initial values ​​for each model parameter are randomly generated based on the range of values ​​for each model parameter in the hydraulic model. Optimal values ​​for each model parameter are then determined using a simulated annealing algorithm. The model parameters include the Manning roughness coefficient for the impermeable zone, the Manning roughness coefficient for the permeable zone, the water storage depth in the depression of the impermeable zone, the water storage depth in the depression of the permeable zone, the maximum infiltration rate, the minimum infiltration rate, the permeability attenuation coefficient, and the drainage time. The process involves randomly generating initial values ​​for each model parameter based on its range, and then determining the optimal values ​​for each model parameter using a simulated annealing algorithm. Specifically: Obtain the actual water depth at each liquid level monitoring point at multiple preset times; Based on the value range of each model parameter of the hydraulic model, a set of initial values ​​for each model parameter are randomly generated as the current solution; Based on the current solution, the hydraulic model is modified to obtain the hydraulic model corresponding to the current solution. Then, using the hydraulic model corresponding to the current solution, the simulated water depth of each liquid level monitoring point at multiple preset times is calculated. The fitness of the current solution is calculated based on the simulated water depth at each liquid level monitoring point at multiple preset times and the actual water depth at each liquid level monitoring point at multiple preset times. A new solution is generated by randomly perturbing the current solution; Based on the new solution, the hydraulic model is modified to obtain the hydraulic model corresponding to the new solution. Then, using the hydraulic model corresponding to the new solution, the simulated water depth of each liquid level monitoring point at multiple preset times is calculated. The fitness of the new solution is calculated based on the simulated water depth at each liquid level monitoring point at multiple preset times and the actual water depth at each liquid level monitoring point at multiple preset times. If the fitness of the new solution is greater than that of the current solution, then according to the Metropolis acceptance criterion, the new solution is adopted as the current solution or the current solution is retained; otherwise, the new solution is adopted as the current solution. If the preset iteration termination condition is met, the current solution is taken as the optimal solution; otherwise, return to the step of generating a new solution by randomly perturbing the current solution. Based on the optimal solution, determine the optimal values ​​of each model parameter; The fitness of the current solution is calculated based on the simulated water depth at each liquid level monitoring point at multiple preset times and the actual water depth at each liquid level monitoring point at multiple preset times. Specifically: According to the Nash coefficient calculation formula, the composite Nash coefficient corresponding to the current solution is calculated, and the negative of the composite Nash coefficient corresponding to the current solution is taken as the fitness of the current solution. The Nash coefficient calculation formula is as follows: ， Where NSE represents the composite Nash coefficient corresponding to the current solution, L obs,i,j L represents the actual water depth at the i-th liquid level monitoring point at the j-th time. pred,i,j L represents the simulated water depth at the i-th liquid level monitoring point at the j-th time corresponding to the current solution. obs,i,mean This represents the average actual water depth at the i-th liquid level monitoring point over n preset times, where num is the number of liquid level monitoring points and n is the total number of preset times. The preset iteration termination conditions include at least one of the following: reaching a preset maximum number of iterations, the current temperature being less than a preset termination temperature, and the difference between the fitness of the new solution and the fitness of the current solution in consecutive preset iteration rounds being less than a preset threshold. Among them, the current temperature T t =a t ×T t-1 a t =a0×exp(t / 50), a t The current cooling rate is given by , a0 is the preset initial cooling rate, t is the current iteration number, and T is the current iteration number. t-1 The temperature at the (t-1)th iteration is T0, where T0 is the preset initial temperature. The process involves randomly generating initial values ​​for each model parameter based on its range, and then determining the optimal values ​​for each model parameter using a simulated annealing algorithm. Specifically: Based on the value range of each model parameter of the hydraulic model, at least one set of initial values ​​for each model parameter is randomly generated to obtain at least one set of initial solutions. Simulated annealing search is performed on at least one set of initial solutions in parallel to obtain at least one set of optimal solutions; Based on the optimal solutions of at least one set, determine the optimal values ​​of each model parameter.

2. The automatic calibration method for drainage network model parameters according to claim 1, characterized in that, Before the step of randomly generating initial values ​​for each model parameter based on the value range of each model parameter of the hydraulic model, and determining the optimal values ​​for each model parameter through simulated annealing algorithm, the method further includes: The integrity and volatility of the actual values ​​at each liquid level monitoring point are assessed. If both the integrity assessment and the volatility assessment pass, then the steps of randomly generating initial values ​​for each model parameter based on the range of values ​​for each model parameter of the hydraulic model, and determining the optimal values ​​for each model parameter through a simulated annealing algorithm are executed.

3. The automatic calibration method for drainage network model parameters according to claim 2, characterized in that, The actual values ​​of each liquid level monitoring point are subjected to integrity assessment and fluctuation assessment, specifically as follows: If the ratio of the total number of actual values ​​of each liquid level monitoring point within the preset time period to the preset monitoring data volume is greater than or equal to the preset ratio, the integrity assessment is deemed to have passed. If a predetermined number of consecutive actual values ​​are the same at a liquid level monitoring point, the fluctuation assessment of that liquid level monitoring point is deemed unsuccessful.

4. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the method as described in any one of claims 1-3.

Images