Double-layer optimization method and system for real-time flood control of series reservoirs

Abstract

The application discloses a double-layer optimization method and system for real-time flood control scheduling of a series reservoir group. The method comprises the following steps: based on the discharge capacity and reservoir capacity constraints of each reservoir, the maximum feasible storage boundary of each period is calculated; a space-time risk response kernel matrix is constructed, the lag response and decay intensity of flood evolution are encoded, and a dynamic risk potential index is constructed in combination with the real-time reservoir capacity; and under the physical boundary constraint, the spatial allocation of water quantity of each reservoir is solved by taking the minimization of the index as a target. The evolutionary algorithm embedded with a physical consistent projection operator is adopted to iteratively optimize the discharge flow of each reservoir. The application effectively solves the problems of upstream and downstream hydraulic space-time coupling and strong physical constraint solving difficulty by matrix coding of space-time risk and introduction of manifold projection solving, and improves the physical consistency and peak shaving efficiency of the scheduling scheme.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of water resources regulation and flood control engineering, and particularly to a double-layer optimization method and system for real-time flood control regulation of a series reservoir group. BACKGROUND

[0002] With the deepening of the development of the river basin, the series reservoir group has become the core backbone project to ensure the flood control safety of the river basin. Under the complex conditions of variable flood in the flood season and close hydraulic connection between the upstream and downstream, how to implement joint real-time regulation of the series reservoir group to reduce the downstream flood peak to the greatest extent on the premise of ensuring dam safety has important engineering and social significance.

[0003] At present, the flood control regulation of the series reservoir group mainly relies on the conventional regulation based on rules or the optimization regulation method based on a single-layer model. The conventional method mainly uses static flood control regulation charts or empirical rules to determine the discharge according to the current water level; the single-layer optimization method usually regards the reservoir group as a whole system, uses general mathematical tools such as dynamic programming and particle swarm algorithm to find the optimal solution of the whole system under given constraints, and directly outputs the discharge process of each reservoir.

[0004] However, the existing scheme has obvious defects in dealing with strong coupling and strong constraint problems, mainly manifested in that it cannot quantify the spatio-temporal heterogeneity risk of flood evolution, and lacks an effective coupling solution mechanism for macro water allocation and micro hydraulic constraints. Specifically, the traditional spatial allocation only allocates tasks according to the static reservoir capacity proportion of the reservoir, ignores the lag and attenuation characteristics of the flood wave in the river evolution, and causes the slow reservoir with large capacity but far distance to be allocated too much task, and cannot play a role before the flood peak arrives. At the same time, the existing intelligent algorithm usually searches in the unconstrained space, relies on the penalty function to deal with the strong physical constraints such as safety discharge and water balance, which leads to a large number of infeasible solutions that violate the physical law, and seriously reduces the solving efficiency of real-time regulation and the executability of the scheme. SUMMARY

[0005] The present application provides a double-layer optimization method and system for real-time flood control regulation of a series reservoir group, in order to solve at least one of the above problems existing in the prior art.

[0006] Technical scheme, according to one aspect of the present application, a double-layer optimization method for real-time flood control regulation of a series reservoir group, comprising:

[0007] Based on the real-time water level of each reservoir, the inflow flood forecast process and the safety discharge threshold of each flood control section downstream, according to the pre-stored discharge capacity constraints and reservoir capacity constraints of each reservoir, the maximum feasible storage boundary of each period is determined, and the excess flood volume vector of each flood control section downstream is identified;

[0008] A spatiotemporal risk response kernel matrix is ​​constructed to characterize the lag response intensity of upstream reservoir outflow changes to downstream flood control section peak flow. This matrix is ​​then superimposed with the real-time reservoir capacity status and excess flood volume vector of each reservoir to construct a dynamic risk potential energy index.

[0009] With the goal of minimizing the dynamic risk potential energy index, the spatial distribution of water volume in each reservoir is solved under the constraint of the maximum feasible retention boundary.

[0010] With the goal of spatial water allocation, an evolutionary algorithm with embedded physical consistent projection operators is used to iteratively optimize the hourly outflow of each reservoir, thereby obtaining an optimized outflow process that satisfies the constraints of water balance and safe discharge.

[0011] According to another aspect of this application, a two-layer optimization system for real-time flood control scheduling of a series reservoir group includes:

[0012] Memory, used to store computer programs;

[0013] A processor is used to implement any of the above methods when executing a computer program.

[0014] Beneficial effects: By using matrix encoding of spatiotemporal risks and introducing manifold projection for solution, this invention effectively solves the problem of difficult solution of upstream and downstream hydraulic spatiotemporal coupling and strong physical constraints, thereby improving the physical consistency and peak shaving efficiency of the scheduling scheme. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the overall process of a two-layer optimization method for real-time flood control scheduling of a series of reservoirs.

[0016] Figure 2 A flowchart illustrating the process of determining the elements in the spatiotemporal risk response kernel matrix.

[0017] Figure 3 A schematic diagram of the process for constructing a dynamic risk potential index.

[0018] Figure 4 A schematic diagram of the process for determining the maximum feasible water retention boundary for each time period. Detailed Implementation

[0019] Example 1 details the overall process of a two-layer optimization method for real-time flood control scheduling of a series of reservoirs, such as... Figure 1 As shown, a two-layer decoupled architecture of spatial allocation and temporal execution is used to realize the coordinated flood control scheduling of a reservoir group under complex and multi-constraint conditions.

[0020] Step 101, based on the real-time water level of each reservoir, the inflow flood forecast process and the safety discharge threshold of each downstream flood control section, according to the pre-stored discharge capacity constraint and storage capacity constraint of each reservoir, the maximum feasible storage boundary of each period is determined, and the excess flood volume vector of each downstream flood control section is identified.

[0021] This step mainly completes the boundary generation and task identification of the double-layer optimization model. Specifically, the real-time water level provides the initial state of the scheduling, and the inflow flood forecast process provides the future inflow situation information. The maximum feasible storage boundary refers to the maximum water quantity that can be actually stored by each reservoir in each period under the premise of meeting the discharge capacity, water level amplitude and other physical constraints. It is different from the static residual storage capacity because the dynamic restriction of discharge capacity is considered. For example, even if the reservoir has a large amount of empty storage capacity, if the inflow is large but the discharge capacity is limited, the actual storage capacity is also limited. By calculating this boundary, it can avoid the upper model to allocate the storage task that cannot be physically executed, and the storage task is the spatial allocation of water quantity. The excess flood volume vector is obtained by comparing the natural flood process of each downstream flood control section with the safety discharge threshold, which quantifies the flood control demand of each downstream section, i.e. the flood level that needs to be reduced by the upstream reservoir group.

[0022] In some optional embodiments, the inflow flood forecast process can be obtained in real time by a data-driven model or a hydrological model. The safety discharge threshold can be dynamically set according to the flood control standard of the downstream flood control object or the temporary flood control situation. The specific method of determining the maximum feasible storage boundary can use dynamic programming, network flow or other optimization algorithms. When identifying the excess flood volume vector, the part exceeding the safety discharge can be integrated to obtain the specific excess flood volume value.

[0023] Step 102, construct a spatio-temporal risk response kernel matrix representing the lag response intensity of the outflow change of the upstream reservoir to the flood peak flow of the downstream flood control section, superimpose the spatio-temporal risk response kernel matrix with the real-time storage state of each reservoir and the excess flood volume vector, and construct a dynamic risk potential index; minimize the dynamic risk potential index as the target, and solve the spatial allocation of water quantity of each reservoir under the constraint of the maximum feasible storage boundary.

[0024] The upper space allocation model is constructed in this step. The time-space risk response kernel matrix encodes the time lag and attenuation characteristics in the complex flood evolution process into a static weight matrix. The elements in the matrix reflect the influence degree of unit outflow change of an upstream reservoir on the flood peak of a downstream section. The dynamic risk potential index is a comprehensive risk measure, which not only considers the current reservoir capacity occupation, i.e. static risk, but also considers the reservoir peak shaving capacity for downstream flood, i.e. dynamic risk, by introducing the time-space risk response kernel matrix. Superimposing the two, the allocation result can not only balance the reservoir capacity risk of each reservoir, but also preferentially use the reservoir with more obvious downstream peak shaving effect. The spatial allocation water quantity solved by minimizing the index is the task of each reservoir in this dispatching.

[0025] In some optional embodiments, the construction of the time-space risk response kernel matrix can be based on the Muskingum method, Saint-Venant equation set or other flood evolution models. In the construction of the dynamic risk potential index, the weights of static risk and dynamic risk can be adjusted according to actual needs. The solving process can adopt linear programming, quadratic programming or intelligent optimization algorithm.

[0026] Step 103, taking the spatial allocation water quantity as the target, an evolutionary algorithm embedded with a physically consistent projection operator is used to iteratively optimize the hourly outflow of each reservoir, and an optimized outflow process satisfying the water balance and safe discharge constraints is obtained.

[0027] The lower time sequence execution model is constructed in this step. The evolutionary algorithm is responsible for searching the optimal outflow process in the complex time sequence solution space. In order to ensure that the searched solution is physically feasible, the present application embeds a physically consistent projection operator in the algorithm. The physically consistent projection operator forces the candidate solution to be projected onto the manifold satisfying the water balance equation and the river safety discharge constraint in each iteration of the algorithm. In this way, not only the convergence speed of the algorithm is improved, but also the final output optimized outflow process has strict physical consistency.

[0028] In some optional embodiments, the evolutionary algorithm can be selected from particle swarm optimization algorithm, genetic algorithm or differential evolution algorithm, etc. The specific implementation of the physically consistent projection operator can be based on gradient projection method, penalty function method or other constraint processing techniques. The finally obtained optimized outflow process can be directly used to guide the real-time dispatching operation of the reservoir.

[0029] Embodiment 2 describes the construction process of basic data processing and physical boundary.

[0030] Step 201, data preprocessing and interval inflow reconstruction are performed.

[0031] To improve the input data quality of the scheduling model, this embodiment first preprocesses the collected raw flood data. Specifically, a joint model of flood signal reconstruction and wavelet threshold denoising is adopted. Using the Daubechies wavelet family, such as the db6 wavelet, the original flow sequence is decomposed into multiple scales. A soft threshold function is used to filter out high-frequency noise before reconstruction to obtain a smooth flow process. Based on this, the interval inflow is reconstructed according to the node flow balance principle. For two adjacent nodes i and i+1, the interval inflow process is obtained by subtracting the outflow from the upstream node i from the inflow at the downstream node i+1. To eliminate errors caused by differential calculations, a Savitzky-Golay filter is further used to smooth the reconstructed interval inflow, and gradient anomaly detection is combined to remove abrupt changes.

[0032] In some alternative implementations, the wavelet basis function and the number of decomposition levels can be adjusted based on the noise characteristics of the data. The anomaly detection threshold can be set based on the statistical regularities of historical data.

[0033] Step 202: Identify the excess flood vector.

[0034] Based on the preprocessed data, the inflow flood processes of each reservoir are projected downstream along the river channel, and the interval inflows are superimposed to obtain the natural flood processes at each flood control section. The natural flood processes are then compared hourly with the safe discharge thresholds for each section. If the natural flow exceeds the safe discharge threshold within a certain time period, the excess is integrated to calculate the excess flood volume for that section. This process is performed on all flood control sections, forming an excess flood volume vector. This vector clarifies the overall flood control task for this operation.

[0035] In some alternative implementations, flood evolution can employ the Muskingan model or the spread wave model. The safe discharge threshold can be a fixed value or a dynamic process that varies over time.

[0036] Step 203: Construct a single-reservoir dynamic programming model with reservoir capacity as the state variable and outflow as the decision variable. In the single-reservoir dynamic programming model, the following constraints are set: the outflow does not exceed the upper limit of the discharge capacity corresponding to the current reservoir capacity; the change in outflow between adjacent times does not exceed the change threshold; and the real-time reservoir capacity is between the preset dead water level and the preset flood control high water level. Under the premise of satisfying the constraints, with the objective of maximizing the difference between the cumulative inflow and cumulative outflow during the prediction period, the maximum feasible water retention boundary for each time period is solved, such as... Figure 4 As shown.

[0037] The calculation method of the maximum feasible storage boundary is described in detail. For each reservoir, a separate dynamic programming model is established. The state variable is selected as the reservoir capacity, and the stage variable is the time period. The state transition equation is based on the water balance principle, i.e. the reservoir capacity at the next time is equal to the current reservoir capacity plus the inflow minus the outflow. In terms of constraint conditions, physical limitations must be strictly followed: the outflow cannot exceed the discharge capacity at that water level, which is a hard constraint determined by the flood discharge facility; the change rate of outflow cannot exceed the set amplitude threshold to ensure the stability of the downstream river channel slope; the reservoir capacity must be limited between the dead storage corresponding to the dead water level and the maximum storage corresponding to the flood control high water level. The optimization goal is to maximize the final cumulative storage by adjusting the outflow strategy within the entire prediction period. By solving the model, the maximum water storage that the reservoir can theoretically store under the current working condition, i.e. the maximum feasible storage boundary, can be obtained. This boundary is usually less than or equal to the static residual flood control storage, as the storage capacity that cannot be utilized due to insufficient discharge capacity or amplitude limitation is deducted.

[0038] In some optional embodiments, dynamic programming can be solved using a discretized state space method, or an approximate method such as differential dynamic programming. To improve computational efficiency, parallel computing techniques can be used to analyze multiple reservoirs simultaneously.

[0039] Embodiment 3 elaborates on the construction of the spatio-temporal risk response kernel matrix, the definition of the dynamic risk potential index, and the application of the projection constraint, demonstrating how to achieve efficient spatial allocation through matrix encoding of physical mechanisms.

[0040] Step 301: The elements in the spatio-temporal risk response kernel matrix represent the cumulative reduction contribution of unit storage of an upstream reservoir to the flood peak risk of a downstream flood control section within the prediction period. The elements in the spatio-temporal risk response kernel matrix are determined as follows:

[0041] Based on the river flood routing model, the partial derivative of the downstream flood control section flow change at the future time caused by the unit outflow change of the upstream reservoir at the current time is calculated to obtain the time-lag response kernel; a time weight function is constructed, which has a higher value when approaching the flood peak time of the downstream flood control section than when far away from the flood peak time; the product of the time-lag response kernel and the time weight function is integrated within the prediction period to obtain the propagation coefficient element in the spatio-temporal risk response kernel matrix corresponding to the upstream reservoir and the downstream flood control section, as shown in Figure 2 .

[0042] The time weight function is constructed based on the flood peak time of the downstream flood control section determined by the reservoir inflow forecast process, and the value of the time weight function near the flood peak time is higher than that far from the flood peak time.

[0043] The calculation process of the elements in the space-time risk response kernel matrix is described in detail. Based on the Muskingum linear flood routing model, the analytical partial derivative of the outflow of the upstream reservoir to the downstream section flow is derived. This partial derivative reflects the physical characteristics of water flow propagation, i.e., how long it takes for a unit flow change to propagate to the downstream and how much it propagates, constituting the time-lag response kernel. In order to reflect the timeliness of flood control dispatching, a time weight function is introduced, and the design principle is to give higher weight to the response near the flood peak time, for example, an inverse function or a Gaussian function can be used, so that the contribution of the peak near the flood peak to the risk is amplified. Multiply the time-lag response kernel by the time weight function, and integrate it in the prediction period to get a scalar value, i.e., the propagation coefficient element. This coefficient comprehensively reflects the effectiveness of the reservoir in reducing the peak of the section. Calculating this coefficient for all reservoirs and all sections can assemble the space-time risk response kernel matrix.

[0044] In some optional embodiments, the time weight function can be specifically represented as w(t)=1 / |t-tp|, where tp is the flood peak time. The integral process can be numerically approximated by using the trapezoidal formula or Simpson formula.

[0045] Alternatively, the calculation of the time-lag response kernel can be performed in the following way. It is assumed that the influence of a unit outflow change of reservoir i at time τ on the flow of downstream flood control section j at time t can be represented by a transfer function h(i,j,t-τ), which is derived based on the Muskingum model. For the Muskingum model, the transfer function can be represented as:

[0046] h(i,j,t-τ)=C0·δ(t-τ-T)+C1·δ(t-τ-T-Δt)+C2·h(i,j,t-τ-Δt);

[0047] where C0, C1, C2 are Muskingum coefficients, T is the average propagation time from reservoir i to section j, Δt is the calculation period length, and δ is the Dirac function.

[0048] Step 302: Calculate the ratio of the sum of the real-time storage capacity of each reservoir and the spatially allocated water volume to be solved to the flood control capacity of that reservoir, obtaining the static storage capacity ratio risk term; for each flood control section, calculate the product of the propagation coefficient element of each reservoir and the excess flood volume of that flood control section, and sum the weighted results of all reservoirs to obtain the spatiotemporal response potential energy term; based on the pre-configured adjustment coefficient, linearly weight the static storage capacity ratio risk term and the spatiotemporal response potential energy term to obtain the dynamic risk potential energy index, such as... Figure 3 As shown.

[0049] This step constructs a comprehensive risk index comprising both static and dynamic components. The static reservoir capacity ratio risk term reflects the reservoir's own water storage pressure; that is, the higher the reservoir capacity utilization rate, the greater the risk. This term is typically obtained by dividing the sum of the current reservoir capacity and the spatially allocated water volume by the total flood control capacity. The spatiotemporal response potential energy term utilizes the propagation coefficient element calculated in step 301. Multiplying the propagation coefficient of each reservoir to each cross-section by the excess flood volume of that cross-section effectively quantifies the reservoir's potential peak-shaving contribution or potential energy when facing flood threats at that cross-section. The total spatiotemporal response potential energy of the reservoir is obtained by weighted summation of the contributions from all cross-sections. Through a pre-configured adjustment coefficient, the static risk and dynamic potential energy are linearly combined to form the dynamic risk potential energy index. In the optimization process, minimizing this index means, while ensuring the safety of each reservoir itself, maximizing the flood control capacity of reservoirs with high peak-shaving efficiency.

[0050] In some alternative implementations, the weights of static risk items can be customized based on the size or importance of the reservoir. Adjustment coefficients can be adjusted according to the dispatchers' preferences regarding upstream and downstream risks.

[0051] Optionally, the complete mathematical expression for the dynamic risk potential index Φ is as follows:

[0052] Φ=α1·∑ i [(V i +S i ) / V _flood,i ] 2 +β1·∑ j [∑ i (K(i,j)·ΔQ j )-S i ·K(i,j)] 2 ;

[0053] Among them, V i S represents the real-time storage capacity of reservoir i. i V represents the spatial allocation of water volume for reservoir i to be solved. _flood,i Let K(i,j) be the flood control capacity of reservoir i, K(i,j) be the element of the spatiotemporal risk response kernel matrix, and ΔQ be the value of Q. jExcess flood volume of section j, α1 and β1 are adjustment coefficients.

[0054] The principle of selecting the values of the adjustment coefficients α1 and β1 is as follows: when the reservoir capacity of the upstream reservoir is tight, α1 is increased to preferentially ensure the safety of the reservoir; when the downstream flood control situation is severe, β1 is increased to preferentially ensure the peak reduction effect of the downstream. The typical values are α1 = 0.3-0.5, β1 = 0.5-0.7, and α1 + β1 = 1.

[0055] In step 303, when solving the spatial allocation water volume of each reservoir, the spatiotemporal risk projection constraint needs to be met, and the spatiotemporal risk projection constraint includes: for each flood control section, the weighted sum of the spatial allocation water volume of each reservoir and the corresponding propagation coefficient element is not less than the product of the excess flood volume of the flood control section and the preconfigured peak reduction proportion coefficient.

[0056] This step introduces the spatiotemporal risk projection constraint to ensure that the allocation scheme can physically meet the downstream flood control demand. The constraint is an inequality, and the left side is the weighted sum of the spatial allocation water volume of all upstream reservoirs and the propagation coefficient thereof for the section, which represents the total effective peak reduction amount of the reservoir group for the section; the right side is the excess flood volume of the section multiplied by a peak reduction proportion coefficient. The physical meaning of the constraint is that the storage task allocated to the upstream reservoir group must at least achieve the predetermined peak reduction target after river evolution and attenuation at the downstream section. In mathematics, it is equivalent to cutting an effective peak reduction hyperplane that meets the downstream safety demand in a high-dimensional solution space, and forcing the solution to fall within this region.

[0057] In some optional embodiments, the peak reduction proportion coefficient is usually set to be between 0.8 and 1.0, depending on the risk tolerance. The constraint can be added to the optimization model as a hard constraint.

[0058] In order to facilitate the understanding of the above spatial allocation process, a simplified numerical example is given below.

[0059] Suppose that a certain river basin has 3 serial reservoirs (A, B, and C) and 2 downstream flood control sections (section 1 and section 2), the prediction period is 24 hours, and the calculation period is 1 hour.

[0060] The initial conditions are as follows: the real-time reservoir capacity of reservoir A is 0.8 billion cubic meters, the flood control capacity is 1.5 billion cubic meters; the real-time reservoir capacity of reservoir B is 0.5 billion cubic meters, the flood control capacity is 1 billion cubic meters; the real-time reservoir capacity of reservoir C is 0.3 billion cubic meters, the flood control capacity is 0.6 billion cubic meters.

[0061] After the flood evolution calculation, the excess flood volume of section 1 is 0.2 billion cubic meters, and the excess flood volume of section 2 is 0.15 billion cubic meters.

[0062] The spatiotemporal risk response kernel matrix K obtained by step 301 is:

[0063] K= ;

[0064] wherein the first row corresponds to section 1, and the second row corresponds to section 2; and the three columns correspond to reservoirs A, B, and C, respectively.

[0065] Let the adjustment coefficients be α = 0.4 and β = 0.6, and use a quadratic programming solver to solve the problem of minimizing the dynamic risk potential energy index.

[0066] The solution is that the spatially allocated water volume of reservoir A is 120 million cubic meters, the spatially allocated water volume of reservoir B is 90 million cubic meters, and the spatially allocated water volume of reservoir C is 60 million cubic meters.

[0067] Verify the spatiotemporal risk projection constraint:

[0068] For section 1, 0.85 x 1.2 + 0.60 x 0.9 + 0.35 x 0.6 = 17.7 million cubic meters, which is greater than 2 x 0.85 = 17 million cubic meters, satisfying the constraint.

[0069] For section 2, 0.70 x 1.2 + 0.80 x 0.9 + 0.55 x 0.6 = 18.9 million cubic meters, which is greater than 1.5 x 0.85 = 12.75 million cubic meters, satisfying the constraint.

[0070] It can be understood that reservoir A, which is closer to the downstream and has a larger propagation coefficient, undertakes more interception tasks, which is consistent with the physical meaning of the spatiotemporal risk response kernel matrix.

[0071] Embodiment 4, detailed description of the evolutionary solving process based on physical manifold projection and ideal trajectory guidance.

[0072] Step 401, based on the hourly discharge of each reservoir in the current iteration, calculate the calculated flow of the downstream flood control section through the river flood evolution model, construct the safety flow residual vector of the downstream flood control section, and the elements in the safety flow residual vector are the parts of the calculated flow exceeding the safety discharge threshold; construct a water balance residual vector, and the elements in the water balance residual vector are the differences between the node calculated flow and the node balanced flow.

[0073] The step defines the process of quantifying the degree of violation of physical constraints as a mathematical residual, which is the basis for constructing a physically consistent projection operator. Specifically, the safety flow residual vector is used to represent the degree of flood control risk caused by the current scheduling scheme downstream. For a certain flood control section, if the calculated flow is less than or equal to the safety discharge threshold, the residual at this point is zero; if it exceeds the threshold, the residual is the value of the excess flow. For example, if the safety discharge of a section is 5000 cubic meters per second, and the calculated flow is 5500 cubic meters per second, the residual element of the section at that time is 500. The water balance residual vector is used to represent the continuity error of flow in the propagation process. In hydraulics, the inflow, outflow and storage rate of any node should satisfy the continuity equation. The elements in the residual vector can be obtained by calculating the inflow minus the outflow minus the storage rate of the node. By constructing the above two vectors, the complex constraint violation situation is converted into an algebraic vector that can be directly involved in numerical calculation, providing a directional guide for subsequent gradient correction.

[0074] Step 402, constructing a sensitivity matrix, which represents the partial derivative of the flow of each downstream flood control section to the outflow of each upstream reservoir.

[0075] This step builds a bridge between control variables and state variables. The sensitivity matrix is a Jacobian matrix, the number of rows corresponds to the number of downstream flood control sections or the number of time periods of interest, and the number of columns corresponds to the outflow decision variables of each upstream reservoir. Each element in the matrix represents the incremental response of the downstream section flow at a certain time when the upstream reservoir increases the outflow by one unit at a certain time. In some preferred embodiments, the matrix can be directly derived by analytical method. For example, based on the Muskingum flow evolution formula, a linear partial derivative relationship can be derived, and the sensitivity matrix is composed of a series of constant coefficients, reflecting the diffusion and superposition characteristics of flood waves. The matrix is constructed to find the outflow adjustment direction that can reduce the residual at the fastest speed.

[0076] Step 403, based on the product of the transpose matrix of the sensitivity matrix and the weighted safety flow residual vector and the water balance residual vector, the gradient of the hourly outflow of each reservoir is corrected to obtain the corrected outflow.

[0077] The step performs a correction operation of the physical consistent projection operator, based on the gradient descent method or the Newton method, and uses the transpose matrix of the sensitivity matrix to propagate the residual signal downstream back to the decision variable space upstream. In specific implementation, first, different weight coefficients are assigned to the safety flow residual vector and the water balance residual vector, for example, in order to give priority to the safety of the dam, a higher weight can be given to the water balance residual. The weighted total residual vector is left multiplied by the transpose of the sensitivity matrix, and the result obtained is the negative gradient direction of the outflow, that is, the steepest descent direction of reducing the default level. The current outflow is reduced by the product of the gradient direction and a learning rate step, and the modified outflow is obtained. Geometrically, this operation is equivalent to pulling the current solution that does not satisfy the constraint back to the vicinity of the constraint surface along the direction perpendicular to the constraint surface, and realizes the projection to the physical flow.

[0078] In a specific embodiment, the formula of gradient correction is as follows. Let the outflow vector of the current iteration be Q, the safety flow residual vector be r _s , the water balance residual vector be r _b , and the sensitivity matrix be J, then the calculation formula of the modified outflow Q' is:

[0079] Q'=Q-η·J T ·(w _s ·r _s +w _b ·r _b );

[0080] wherein η is a learning rate step, typically 0.01 to 0.1, J T represents the transpose matrix of the sensitivity matrix J, w _s and w _b are weight coefficients of the safety flow residual and the water balance residual, typically w _s =0.3-0.5, w _b =0.5-0.7.

[0081] The construction of the sensitivity matrix J is based on the linear characteristics of the Muskingum model. The matrix element J(j,i,t,τ) represents the change of the flow at section j at time t caused by the change of the outflow of reservoir i at time τ, and its calculation formula is the same as that of the time-lag response kernel.

[0082] In step 404, the modified outflow is proportionally scaled and projected to obtain an outflow that satisfies the volume conservation constraint, according to the volume constraint relationship between the spatially allocated water volume of each reservoir and the initial reservoir capacity and the target reservoir capacity.

[0083] This step is a further correction of the correction result to ensure that the total water allocation task of the upper model is met. After gradient correction, the total integral of the outflow sequence may drift, and it is no longer equal to the spatial allocation water quantity allocated by the upper layer. In order to solve this problem, volume conservation projection needs to be performed. Specifically, first, the total discharge of the corrected outflow sequence in the entire scheduling period is calculated, and the proportional coefficient between the total discharge and the target total discharge, i.e. the initial reservoir capacity plus the total inflow minus the target final reservoir capacity, is calculated. The proportional coefficient is used to multiply each element in the corrected outflow sequence. After this step, the final outflow not only satisfies the safety and balance differential constraints locally, but also satisfies the integral constraints of water allocation globally, achieving physical consistency.

[0084] Step 405, based on the inverse flood routing model, the ideal outflow trajectory of each reservoir is inversely deduced using the spatial allocation water quantity of each reservoir.

[0085] This step is used to generate a reference benchmark for guiding the search of the evolutionary algorithm. The inverse flood routing model refers to a mathematical model that inversely deduces the ideal discharge process of an upstream reservoir given the downstream flow process or flood control requirements. For example, the improved inverse Muskingum method can be used, and its difference equation is in the form of using the inflow and outflow at the next time and the outflow at the current time to deduce the inflow at the current time. In order to further improve the inversion accuracy, a long short-term memory network (LSTM) or other data-driven model can be combined to compensate for the nonlinear error of the inverse calculation. Specifically, the ideal outflow trajectory can be represented as the weighted sum of the physical inversion result and the data-driven error compensation term. This ideal trajectory represents the discharge strategy that should be taken to cooperate with the downstream flood control target in an ideal situation without considering specific operation restrictions of the reservoir, such as gate constraints.

[0086] Step 406, construct an ideal solution projection operator, which linearly interpolates the hourly outflow of each reservoir to the ideal outflow trajectory.

[0087] This step defines an operator that guides search individuals to high-quality regions. In the random search process of the evolutionary algorithm, in order to avoid particles wandering blindly in a huge solution space, the ideal outflow trajectory is used as a lighthouse for guidance. Specifically, for any candidate outflow sequence, the ideal solution projection operator performs linear interpolation between it and the ideal outflow trajectory. For example, the new outflow is equal to the original outflow multiplied by a reservation coefficient plus the ideal outflow trajectory multiplied by a guidance coefficient. The size of the guidance coefficient determines the degree of dependence on prior knowledge. In the early iterations, a larger guidance coefficient can be set to quickly converge; in the later iterations, the guidance coefficient is reduced to maintain the diversity of the population.

[0088] Step 407, in each iteration update of the evolutionary algorithm, the interpolation operation of the ideal solution projection operator and the correction operation of the physically consistent projection operator are executed in turn.

[0089] This step describes the integration of the double-layer projection operator in the evolutionary algorithm. Taking the improved particle swarm optimization algorithm PSO as an example, in each iteration, the particles perform standard position update according to their own speed and historical optimal position. Then the interpolation operation of the ideal solution projection operator is performed on the updated position, so that the particles are attracted to the ideal hydraulic response mode. Next, the correction operation of the physically consistent projection operator is performed on the interpolated position, forcing the particles to be pulled back into the feasible region that satisfies the water balance and safe discharge. In addition, in order to balance the global development and local mining ability, a parameter control mechanism based on group information entropy can also be introduced. Specifically, the information entropy of the population distribution is calculated, and when the information entropy is too small, the social learning factor is increased to prevent premature convergence; when the information entropy is too large, the cognitive learning factor is increased to speed up the convergence speed. Through the mechanism of standard update plus double projection, the algorithm can find a scheduling scheme that satisfies the strict physical constraints and meets the global optimization objective with high efficiency. In this process, the objective function of the lower layer optimization can be set as minimizing the sum of squares of the combined flow at each control section downstream, in order to achieve the flood peak mitigation.

[0090] According to one aspect of the present application, a typical flood control scheduling scenario is used to illustrate the complete execution process of the method of the present application.

[0091] Suppose there are three cascade reservoirs in a river basin, including an upstream reservoir, a midstream reservoir, and a downstream reservoir, and a city flood control section is set downstream with a safe discharge of 8000 cubic meters per second. Currently, it is the main flood season, and the meteorological department predicts that there will be a large-scale rainfall process in the next 24 hours.

[0092] Step 1, the scheduling system obtains the real-time water level of the three reservoirs from the water regime telemetry system, and obtains the future 24-hour inflow flood process from the flood forecasting system. The predicted peak flow of the upstream reservoir is 12000 cubic meters per second, and the peak time is the 8th hour; the predicted peak flow of the midstream reservoir is 15000 cubic meters per second, and the peak time is the 12th hour; the predicted peak flow of the downstream reservoir is 18000 cubic meters per second, and the peak time is the 16th hour.

[0093] The system calculates the maximum feasible storage boundary of each reservoir, which is 300 million cubic meters for the upstream reservoir, 250 million cubic meters for the midstream reservoir, and 180 million cubic meters for the downstream reservoir.

[0094] The system identifies that the excess flood volume of the city section is 220 million cubic meters (the natural peak flow is 11000 cubic meters per second, which exceeds the safe discharge of 3000 cubic meters per second, and lasts for about 20 hours).

[0095] Step 2, construct the spatio-temporal risk response kernel matrix, the propagation coefficient of the upstream reservoir is 0.75 (far distance, large attenuation), the propagation coefficient of the middle reservoir is 0.88, and the propagation coefficient of the downstream reservoir is 0.95 (short distance, small attenuation).

[0096] Solve the dynamic risk potential index minimization problem to obtain the spatially distributed water volume: 0.8 billion cubic meters for the upstream reservoir, 1.0 billion cubic meters for the middle reservoir, and 0.6 billion cubic meters for the downstream reservoir. Since the propagation coefficient of the downstream reservoir is the largest, the unit storage capacity of the downstream reservoir contributes the most to peak reduction, so the downstream reservoir is preferentially utilized; meanwhile, considering the capacity limit of the downstream reservoir, the remaining task is allocated to the middle and upstream reservoirs.

[0097] Step 3, the system uses the spatially distributed water volume as a guiding target, and adopts a particle swarm optimization algorithm embedded with a physically consistent projection operator to optimize the hourly outflow of the three reservoirs.

[0098] After 50 iterations, the algorithm converges. The upstream reservoir reduces the flood peak from the 6th to the 10th hour, and the outflow is controlled below 9,000 cubic meters per second; the middle reservoir reduces the flood peak from the 10th to the 14th hour, and the outflow is controlled below 11,000 cubic meters per second; and the downstream reservoir reduces the flood peak from the 14th to the 18th hour, and the outflow is controlled below 7,500 cubic meters per second.

[0099] Step 4, perform flood routing calculation on the optimized outflow process to verify that the maximum flow at the urban cross-section is 7,800 cubic meters per second, which is lower than the safe discharge of 8,000 cubic meters per second, meeting the flood control requirements.

[0100] Compared with the traditional static reservoir capacity proportional allocation method, the present invention method preferentially utilizes the downstream reservoir with a large propagation coefficient to achieve peak reduction in advance, avoiding the problem of slow response of the upstream reservoir.

[0101] Embodiment 5, the convergence control logic, feedback correction mechanism, and specific system hardware architecture of the double-layer optimization framework are described in detail.

[0102] Step 501, inner-layer inspection: calculate the flow balance residual of the upstream and downstream connected cross-section, if the residual exceeds the preset threshold, then correct the downstream inflow process according to the residual and re-execute the evolutionary algorithm.

[0103] The inner loop is mainly used to solve the problem of nonlinear consistency of hydraulic connection. In a series of reservoirs, the outflow of the upstream reservoir evolves into the inflow of the downstream reservoir through the river channel. This process is often affected by nonlinear factors and produces deviations. The inner loop monitors these deviations by calculating the flow balance residual of the connecting cross-section. The preset threshold can be set according to the accuracy of flow measurement, for example, 100 cubic meters per second. If the residual exceeds the threshold, it indicates that the hydraulic connection between the upstream and downstream reservoirs has been broken. At this time, the feedback correction formula is used to adjust the inflow process of the downstream reservoir. The correction formula can be:

[0104] ;

[0105] wherein, is the inflow process of the downstream reservoir / cross-section in the kth iteration, is the flow residual process calculated after the kth iteration, and a is the relaxation factor (correction coefficient), usually between 0.8 and 0.9, to ensure the stability of the iteration.

[0106] After correction, the lower evolution algorithm is restarted with the new inflow boundary until the flow balance residual meets the requirements.

[0107] Step 502, outer loop test: calculate the change rate of excess flood volume at the downstream flood control cross-section in the adjacent two iterations. If the change rate is greater than the convergence threshold, modify the excess flood volume vector and return to the step of constructing the dynamic risk potential index.

[0108] The outer loop is mainly used to solve the coordination problem between the upper and lower layers. Although the upper model allocates water, the lower model may not be able to achieve the intention of the upper layer due to various micro constraints, or it may achieve the intention but there is still room for optimization in downstream flood control. The outer loop monitors the change trend of excess flood volume at the downstream flood control cross-section to determine whether the system has reached the global optimum. Calculate the change rate of excess flood volume in the previous two iterations. If the change rate is greater than the preset convergence standard, it indicates that the system state is still in intense adjustment. At this time, the unsolved flood volume deviation needs to be fed back to the upper model. Specifically, the remaining excess flood volume increment at the downstream is added to the original excess flood volume vector, and then the step of constructing the spatiotemporal risk response kernel matrix is returned to reconstruct the dynamic risk potential index and solve a new spatial allocation scheme. Through the above closed-loop iteration of allocation and redistribution, the system finally converges to a state of consistent coordination between the upstream and downstream.

[0109] Optionally, as an alternative or supplement to the upper-middle-layer space allocation model, a time-varying constraint decoupling mechanism driven by a flood evolution response coefficient matrix can also be used. Specifically, if a global dynamic risk potential index is not constructed, the response coefficient of flood evolution can be directly converted into an effective reservoir capacity constraint of each reservoir at each time period. That is, the actual retention capacity of a reservoir at a certain time period shall not exceed the product of its physical residual capacity and the response coefficient of the time period to the downstream section. This method converts the complex spatio-temporal coupling relationship into a series of time-varying upper bound constraints, which can also decouple the spatio-temporal lag effect and can be used as another embodiment of the present application in a simplified scenario.

[0110] The embodiment also provides a two-layer optimization system for real-time flood control scheduling of a series of reservoirs, comprising:

[0111] A memory for storing a computer program, and a processor for executing the computer program to implement the steps of any method of the present application.

[0112] Specifically, the system comprises a memory and a processor. The memory is used to store a computer program containing instruction codes for implementing the methods of the above embodiments. The processor, when executing the computer program, can implement all steps such as boundary generation, space allocation, evolution solving, and convergence control. In terms of hardware configuration, in order to meet the timeliness requirements of real-time scheduling, the system preferably uses a high-performance computing architecture containing a graphics processing unit (GPU). By using the parallel computing capability of the GPU, the particle swarm in the lower-level evolution algorithm can be massively parallel accelerated, the computation time is shortened, and it is ensured that scheduling instructions can be generated in time in the case of rapid flood evolution in an emergency. In addition, the system is also configured with a standard data communication interface for real-time reception of water and rainfall telemetry data and issuance of scheduling instructions.

[0113] Embodiment 6 describes the details of the high-performance parallel computing architecture and the adaptive evolution mechanism.

[0114] Step 601: Construct a CPU-GPU heterogeneous layered parallel computing architecture.

[0115] To meet the stringent timeliness requirements of real-time flood control scheduling, this embodiment employs a heterogeneous computing strategy to hardware map the two-layer optimization model. Specifically, the computational tasks are divided into two categories: logic-intensive and computationally-intensive. For the upper-layer spatial allocation model, which involves complex logical judgments, risk rule reasoning, and sparse matrix operations, it is allocated to a multi-core CPU cluster for execution. The CPU is responsible for task distribution, parameter broadcasting, and merging of final results. For the lower-layer time-series execution model, which involves massive particle swarm iterations and time-period water balance calculations, and where the calculations between particles are independent, it is mapped to a GPU graphics processing unit for large-scale parallel acceleration. In the specific implementation, the GPU's streaming multiprocessor cores are used to update the position and velocity of all particles in parallel in a vectorized manner, and information exchange between particles is achieved through the GPU's high-speed shared memory. Through this hierarchical mapping, optimal allocation of computing resources is achieved.

[0116] Step 602: Adaptive adjustment of evolutionary parameters based on population information entropy.

[0117] To achieve a dynamic balance between global search and local exploitation when executing the evolutionary algorithm to solve the lower-level model, this embodiment introduces a parameter control mechanism based on information entropy. Specifically, in each iteration, the information entropy of the population position distribution is calculated. The formula for calculating information entropy can be expressed as the negative of the logarithmic expectation of the probability distribution, i.e.:

[0118] H=-∑p m lnp m ; where p m This represents the probability that the particle appears in the m-th subinterval.

[0119] In this embodiment, the probability distribution p of the particles within the pre-divided sub-intervals is... m The algorithm calculates the population information entropy H. A large information entropy H indicates a dispersed population distribution. In this case, the cognitive learning factor c1 should be increased and the social learning factor c2 decreased to encourage particles to move closer to their historical optimal solutions and enhance local development. Conversely, a small information entropy H indicates a clustered population. In this case, the social learning factor c2 should be increased and random perturbations introduced to guide particles towards the global optimal solution and prevent premature convergence. Through these dynamic adjustments, the algorithm can adaptively change its evolutionary strategy based on the search state.

[0120] Step 603: Calculate the convergence criterion based on the sliding time window and fuzzy membership degree.

[0121] To improve the robustness of convergence judgment and avoid misjudgments due to numerical fluctuations in a single iteration during the outer-layer verification of the two-layer optimization, this embodiment employs a sliding time window mechanism. Specifically, a time window of length L (e.g., 3 to 5 iterations) is set, and the change sequence of excess flood volume at the downstream flood control section within this window is recorded. The average rate of change of this sequence within the window is calculated. Furthermore, a fuzzy membership function is introduced to evaluate the convergence state. The membership function is defined as μW=exp(-|rW|). γ ), where rW is the average rate of change and γ is the control parameter. The system is considered to have reached outer-layer convergence only when the membership degree μW is greater than a preset confidence threshold (e.g., 0.85). This processing method effectively filters out high-frequency oscillation noise during the iteration process, ensuring the stability of the final scheduling scheme.

[0122] According to one aspect of this application, the method for spatial allocation of excess flood volume based on spatiotemporal risk response kernel can further be as follows:

[0123] By calculating the hydraulic propagation, a time-delay response kernel for the reduction in unit outflow to the flood peak risk at each cross section is constructed, and a time-weighted integral is performed to form a spatiotemporal risk propagation coefficient matrix.

[0124] By superimposing the static storage capacity ratio risk with the spatiotemporal risk propagation coefficient, a new dynamic risk level is constructed. The greater the peak reduction capacity of the dam and the larger the current excess flood volume downstream, the higher the potential risk of the task it undertakes.

[0125] By incorporating spatiotemporal risk projection constraints into the feasible domain, it is required that the water volume combination allocated at the upper level must be able to provide sufficient peak-shaving response capability at key cross sections, thus forming a truly spatiotemporally coordinated spatial allocation model.

[0126] A new dynamic risk objective function is adopted to replace the original simple reservoir capacity ratio risk objective, so that the upper-level spatial scheduling and the lower-level dual-section peak shaving form an endogenous coupling.

[0127] The detailed processing procedure is as follows:

[0128] First, consider the i-th reservoir at a certain moment. Outbound flow The natural or calculated flow rate at the m-th flood control section at time t+τ The marginal impact is defined by the time-delay response kernel:

[0129] ;

[0130] in, Indicates the first The reservoir is at all times Changes in outbound shipments from the upper unit affect the first The discharge of the mth cross section at time t + τ The partial derivative of the upper flow variation, i represents the reservoir number, the value range is 1 to n, n represents the total number of reservoirs, m represents the flood control cross section number, the value range is 1 to M, M represents the total number of cross sections, τ represents the time lag variable, τ takes the value range of 0 to the end of the prediction period, t represents the current time, The discharge of the mth cross section at time t + τ The discharge of the mth cross section at time t + τ

[0131] In order to reflect the risk more sensitive to the period near the flood peak, a time weight function is introduced :

[0132] ;

[0133] Wherein, The time weight function takes a larger value at the time lag value close to the flood peak, and a smaller value at the time lag value far from the flood peak.

[0134] Further, given the prediction period length , the risk propagation coefficient is defined:

[0135] ;

[0136] Wherein, The mth cross section at time t + τ The mth cross section at time t + τ The mth cross section at time t + τ The total length of the prediction period.

[0137] All Form a matrix:

[0138] ;

[0139] Wherein, B represents the space-time risk propagation coefficient matrix, the square brackets Indicate that the matrix has n rows and M columns.

[0140] On the basis of the original static reservoir capacity ratio risk, superimpose the space-time risk propagation term, and construct a new dynamic risk degree:

[0141] ;

[0142] Wherein, The mth cross section at time t + τ The mth cross section at time t + τ The mth cross section at time t + τ The mth cross section at time t + τ denoted by , where β represents the upper limit of the effective flood control capacity of the i-th reservoir, and β represents the adjustment coefficient used to balance the relative weights of the static capacity ratio and the spatiotemporal risk propagation term. Indicates taking The larger of the two values ​​is 0. , This indicates that among all reservoirs numbered j, for the first... The maximum value of the positive risk propagation coefficient of each cross section. Indicates the first Total excess flood volume at each flood control section This represents the sum of excess flood volume at all flood control sections.

[0143] The meaning of this dynamic risk level is: the first item still reflects the risk of storage capacity ratio; the second item is... and It emphasizes the potential risks of undertaking tasks on sections with concentrated excess flood volume and strong peak-shaving effects.

[0144] The new upper-level optimization objective function is defined as follows:

[0145] ;

[0146] in, This represents the improved value of the upper-level objective function. This represents the static weighting coefficient assigned to the i-th reservoir. This coefficient can be set based on factors such as the catchment area, reservoir capacity, and project importance.

[0147] While retaining the original volume constraints and single-library upper limit constraints, a spatiotemporal risk projection constraint is introduced:

[0148] ;

[0149] Volume constraints: Single database upper limit constraint: ;

[0150] in, Indicates the first The peak reduction ratio coefficient for each flood control section ranges from 0 to 1. For example, 0.7 means that at least 70% of the excess flood volume must be reduced. This indicates that all reservoir allocation schemes are in the [number]th [year]. The comprehensive peak shaving potential of each cross section, WI i This represents the maximum flood storage capacity of the i-th reservoir during the current time period.

[0151] This set of constraints requires that the excess flood volume allocated by the upper level must provide sufficient peak-shaving response capacity at each critical section to ensure risk feasibility, not just volume feasibility.

[0152] By the above method, in the risk measurement, the static index is upgraded to the dynamic risk degree sensitive to space and time , the reservoirs with strong peak clipping ability to key sections in key periods and the reservoirs with limited peak clipping contribution can be distinguished, and the consistency of risk assessment and task allocation is improved. In the feasible region structure, by matrix and constraint , the lower hydraulic propagation property is explicitly raised to the upper decision, reducing the situation that the upper allocation is feasible but the lower effective peak clipping is difficult to achieve, so that the spatial allocation and time sequence peak clipping form a tight coupling. In addition, this method is no longer a simple combination or fine tuning of the existing risk target + volume constraint mode, but introduces new spatiotemporal risk response kernel and risk potential field concepts, changing the structure and information input of the upper model.

[0153] According to one aspect of the present application, the hierarchical asynchronous parallel solving method of physically consistent projection and ideal domain evolution can be specifically:

[0154] The method introduces physically consistent projection operators and ideal solution projection operators, and embeds them into the particle swarm update process, so that each particle position update contains the following two actions:

[0155] First, through the physically consistent projection operator, the upper spatial allocation variable and the lower reservoir outflow process are projected from the original updated position to the physically feasible manifold that approximately satisfies the safety flow constraint, the water balance constraint and the volume conservation constraint;

[0156] Second, through the ideal solution projection operator, the particle is locally projected to the ideal solution direction according to the upper Lagrange solution and the lower ideal outflow trajectory, so that the particle is close to the physically feasible region and the ideal response.

[0157] Finally, an integrated evolution mechanism of ideal solution guidance and physical manifold constraint is formed to replace the original loose coupling mode.

[0158] The detailed processing process is as follows:

[0159] Let the upper allocation vector corresponding to the current particle be:

[0160] ;

[0161] Wherein, represents the candidate allocation vector without projection, represents the tentative allocation excess flood of the th reservoir, and the superscript represents the original candidate value, represents the vertical ellipsis, and n represents the total number of reservoirs.

[0162] The allocation variable needs to satisfy the following constraints:

[0163] Overall volume constraint: ;

[0164] Single library boundary constraints: ;

[0165] Spatiotemporal risk projection constraints (and spatiotemporal risk propagation coefficient matrix) coupling):

[0166] ;

[0167] in, This represents the physically feasible allocation after projection. This indicates the total excess flood volume at the terminal public flood control section. This represents the maximum flood storage capacity of the i-th reservoir during the current time period. Indicates the spatiotemporal risk propagation coefficient. This represents the peak reduction ratio coefficient for the m-th flood control section. Let M represent the excess flood volume at the m-th flood control section, and M represent the total number of flood control sections.

[0168] Therefore, the above constraints are uniformly written in matrix form:

[0169] ;

[0170] Where X represents the assignment vector to be determined, A represents the matrix composed of the coefficients of each inequality constraint, d represents the right-hand constant vector of the inequality constraint, B represents the matrix composed of the coefficients of the equality constraint, and e represents the right-hand constant vector of the equality constraint.

[0171] Define the upper-level physical consistent projection operator as:

[0172] ;

[0173] in, Indicates the candidate vector The mapping result projected onto the physically feasible region, where arg min represents the value of the variable that minimizes the objective function. Let denote the Euclidean norm, and st denote the following constraints.

[0174] This operator represents searching among all assignment vectors that satisfy the physical constraints for a vector that matches the original candidate vector. The point with the smallest distance is used as the projection result.

[0175] Furthermore, let's assume that for the first... The candidate outflow sequence for each reservoir during the entire scheduling period is as follows: Organize the outflows from all reservoirs into a vector:

[0176] ;

[0177] The flow vector of the flood control section is obtained by one-dimensional hydraulic propagation or Muskingum model:

[0178] ;

[0179] wherein, Q(t) represents the flow vector of each flood control section at time t, H represents the hydraulic transfer matrix, the elements of which describe the linear influence coefficient of reservoir outflow on the flow of each section, and q(t) represents the interval inflow vector of each section at time t. The safety flow constraint residual is defined as:

[0180] ;

[0181] wherein, e(t) represents the safety flow constraint residual of the mth section at time t, q(t) represents the candidate flow of the mth section at time t, q(t) represents the allowable safety flow of the mth section at time t, and max(a, b) represents the larger value of real numbers a and b.

[0182] The water balance residual is defined as: ; wherein, e(t) represents the flow balance residual of the mth section at time t, and q(t) represents the theoretical flow value of the mth section calculated according to the node flow balance equation. The two types of residuals are linearly combined as: ;

[0183] wherein, e(t) represents the comprehensive residual vector at time t, e(t) represents the vector composed of all e(t), and e(t) represents the vector composed of all e(t).

[0184] ;

[0185] wherein, e(t) represents the comprehensive residual vector at time t, e(t) represents the vector composed of all e(t), and e(t) represents the vector composed of all e(t). ; wherein, e(t) represents the comprehensive residual vector at time t, e(t) represents the vector composed of all e(t), and e(t) represents the vector composed of all e(t).

[0186]

[0187] ; wherein, e(t) represents the comprehensive residual vector at time t, e(t) represents the vector composed of all e(t), and e(t) represents the vector composed of all e(t).

[0188] ; wherein, e(t) represents the comprehensive residual vector at time t, e(t) represents the vector composed of all e(t), and e(t) represents the vector composed of all e(t). ; and​​​​​​​​​​​​​​ The weighting coefficients represent the two types of residuals and are used to balance the importance of safety constraints and water balance constraints.

[0189] Consider the first The sensitivity matrix of a reservoir for its outflow sequence is defined as follows:

[0190] ;

[0191] in, Indicates at time The flow vector at each cross section is related to the first The reservoir is at all times The partial derivative matrix of the outbound shipment, This represents the partial derivative with respect to the outbound variable.

[0192] Give a physical consistency correction:

[0193] ;

[0194] in, Let η represent the outflow from the i-th reservoir at time t after one correction, and let η represent the correction step size coefficient. Representation matrix The transpose of .

[0195] Furthermore, to satisfy the volume conservation constraint, for Perform a proportional projection:

[0196] ;

[0197] in, Indicates the first The reservoir is at all times The final projected outbound flow rate Let Δt represent the target reservoir capacity at the end of the scheduling period for the i-th reservoir, and let Δt represent the time step length. The numerator of the fraction is... This represents the volume that the reservoir must release during the entire scheduling period, with the denominator being... This indicates the actual release volume corresponding to the corrected outbound process.

[0198] Define the lower-level physical consistency projection operator as:

[0199] ;

[0200] in, This represents the physically consistent projection operator acting on the outbound process.

[0201] Based on the inverse Muskingan model and the data-driven model, we can obtain the first... Ideal outflow trajectory of each reservoir On this basis, the ideal solution projection operator is defined:

[0202]

[0203] where, represents the candidate release process after one linear interpolation to the ideal trajectory , and a represents the interpolation weight coefficient, which ranges from 0 to 1. The larger a is, the closer the projection result is to the ideal trajectory.

[0204] For each particle update of the lower layer, the following composite projection is used:

[0205]

[0206] where, represents the final release process of the th reservoir after the th iteration, represents the original candidate release process obtained according to the particle swarm velocity update formula, and the superscript k+1 represents the iteration number index. The composite function form represents first doing the ideal solution projection and then doing the physically consistent projection.

[0207] For the upper layer allocation variable, the particle swarm position update formula is:

[0208]

[0209]

[0210]

[0211] where, represents the particle velocity vector at the kth iteration, represents the particle position vector at the kth iteration, and ω represents the inertia weight coefficient, and represent the individual learning factor and the group learning factor, and represent random numbers independently and uniformly distributed between 0 and 1, represents the historical optimal position of the particle, represents the global optimal position, represents the unprojected candidate position, represents the final position after the physically consistent projection.

[0212] By introducing the physically consistent projection operator and , as well as the ideal solution projection operator ​​​​​Furthermore, by embedding this into the particle swarm optimization (PSO) position update process, the existing two-layer scheduling + PSO solution structure has been upgraded from a loosely coupled mode of intelligent algorithm + outer-layer physical verification to a tightly coupled mode that evolves along the ideal response direction on the physically feasible manifold. This achieves the following: most candidate solutions of the particle swarm are projected into regions that approximately satisfy the constraints of safe flow, water balance, and volume conservation during the generation stage, reducing the number of obviously infeasible solutions, lowering the frequency of backoff triggered by the outer-layer convergence criterion, and improving real-time performance; the lower-layer outflow process simultaneously approaches the ideal outflow trajectory and the physically feasible region in each iteration, resulting in a final scheduling scheme that is not only mathematically close to optimal but also physically easier to interpret and more in line with the cognition and experience of flood control dispatchers; the upper-layer allocation variables are... Spatiotemporal risk propagation matrix The coupling allows the spatial allocation and the lower hydraulic response to form a closed loop at the solution algorithm level, which is significantly structurally different from the traditional approach of using physical models only in the objective function or constraints.

[0213] This improvement transforms the relationship between the physical model and the intelligent algorithm from a serial connection to an embedded coupling, forming an ideal domain evolution solution framework on the physical manifold.

[0214] This application proposes to construct a spatiotemporal risk response kernel matrix and a dynamic risk potential index. Instead of relying solely on static reservoir capacity ratios, it encodes the lag time and attenuation degree of water flow propagation into the allocation matrix through convolution integrals. This enables scheduling decisions to perceive spatiotemporal distances, ensuring that peak shaving tasks are allocated to reservoirs that can respond to flood peaks in a timely manner hydraulically. This solves the problem of mismatch between spatial allocation and temporal evolution, and addresses the issue of the inability to quantify the spatiotemporal heterogeneity risk of flood evolution.

[0215] This application employs maximum feasible storage boundary analysis and a physically consistent projection operator. First, dynamic programming is used to establish physically executable reservoir capacity boundaries, preventing upper-level decisions from becoming unenforceable. In the lower-level solution, the physically consistent projection operator forces the algorithm's search path to be mapped to a physical manifold that satisfies water balance and safe discharge, abandoning the inefficient penalty function method and ensuring that the output scheduling process strictly satisfies physical constraints. This achieves a unity between macroscopic decision-making and microscopic execution. It solves the problems of lacking a solution mechanism that effectively couples macroscopic water allocation with microscopic hydraulic constraints, and generating a large number of infeasible solutions.

[0216] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A two-level optimization method for real-time flood control scheduling of a series of reservoirs, characterized in that, include: Based on the real-time water level of each reservoir, the inflow flood forecast process, and the safe discharge threshold of each downstream flood control section, the maximum feasible interception boundary for each time period is determined according to the discharge capacity constraints and storage capacity constraints of each reservoir, and the excess flood vector of each downstream flood control section is identified. A spatiotemporal risk response kernel matrix is ​​constructed to characterize the lag response intensity of upstream reservoir outflow changes to downstream flood control section peak flow. This matrix is ​​then superimposed with the real-time reservoir capacity status and excess flood volume vector of each reservoir to construct a dynamic risk potential energy index. With the goal of minimizing the dynamic risk potential energy index, the spatial distribution of water volume in each reservoir is solved under the constraint of the maximum feasible retention boundary. With the goal of spatial water allocation, an evolutionary algorithm with embedded physical consistent projection operators is used to iteratively optimize the hourly outflow of each reservoir to obtain an optimized outflow process that satisfies water balance and safe discharge constraints. The elements in the spatiotemporal risk response kernel matrix represent the cumulative reduction contribution of a unit water storage capacity of the upstream reservoir to the flood peak risk at the downstream flood control section during the forecast period, and are determined in the following way: Based on the river flood evolution model, the partial derivative of the flow change at the downstream flood control section caused by the change in unit outflow from the upstream reservoir at the current moment is calculated to obtain the time-delay response kernel. Construct a time weighting function, whose value is higher when the flood peak is near the downstream flood control section than when the flood peak is far away. Integrating the product of the time-delay response kernel and the time weighting function over the prediction period yields the propagation coefficient elements in the spatiotemporal risk response kernel matrix corresponding to the upstream reservoir and the downstream flood control section. Construct dynamic risk potential indicators, including: Calculate the ratio of the sum of the real-time storage capacity of each reservoir and the spatially allocated water volume to be solved to the flood control capacity of the reservoir, and obtain the static storage capacity ratio risk item; For each flood control section, the product of the propagation coefficient element of each reservoir and the excess flood volume of the flood control section is calculated, and the product results of all reservoirs are weighted and summed to obtain the spatiotemporal response potential energy term. Based on the pre-configured adjustment coefficient, the static storage capacity ratio risk item and the spatiotemporal response potential energy item are linearly weighted to obtain the dynamic risk potential energy index. When determining the spatial distribution of water volume in each reservoir, spatiotemporal risk projection constraints must be satisfied, including: For each flood control section, the weighted sum of the spatially allocated water volume of each reservoir and the corresponding propagation coefficient element shall not be less than the product of the excess flood volume of that flood control section and the pre-configured peak reduction ratio coefficient.

2. The method according to claim 1, characterized in that, The hourly outflow from each reservoir is iteratively optimized, including: The calculated flow rate of the downstream flood control section is calculated based on the hourly outflow rate, and the safe flow rate residual vector of the downstream flood control section is constructed, where the elements are the portion of the calculated flow rate that exceeds the safe discharge threshold. Construct a water balance residual vector, where each element represents the difference between the calculated flow rate at each node and the balanced flow rate at each node. A sensitivity matrix is ​​constructed to characterize the partial derivatives of the flow at each downstream flood control section with respect to the outflow from each upstream reservoir.

3. The method according to claim 2, characterized in that, Further includes: The hourly outflow of each reservoir is graded and corrected based on the product of the transpose of the sensitivity matrix and the weighted safety flow residual vector and the water balance residual vector. Based on the spatial allocation of water volume in each reservoir and the volume constraints of the initial and target reservoir capacities, the corrected outflow is scaled and projected proportionally to obtain the outflow that satisfies the volume conservation constraint.

4. The method according to claim 3, characterized in that, The evolutionary algorithm also incorporates an ideal trajectory guidance mechanism, including: Based on the reverse flood calculation model, the ideal outflow trajectory of each reservoir is inferred by using the spatial distribution of water volume in each reservoir. Construct an ideal solution projection operator, which linearly interpolates the hourly outflow of each reservoir towards the ideal outflow trajectory; In each iteration of the evolutionary algorithm, the interpolation operation of the ideal solution projection operator and the correction operation of the physically consistent projection operator are executed sequentially.

5. The method according to claim 1, characterized in that, Determine the maximum feasible water retention boundary for each time period, including: Construct a single-reservoir dynamic programming model with reservoir capacity as the state variable and outflow as the decision variable; In the single-reservoir dynamic programming model, the following constraints are set: the outflow does not exceed the upper limit of the discharge capacity corresponding to the current reservoir capacity, the change in outflow between adjacent times does not exceed the change threshold, and the real-time reservoir capacity is between the preset dead water level and the preset flood control high water level. Under the premise of satisfying the constraints, the maximum feasible water retention boundary for each time period is obtained by maximizing the difference between the cumulative inflow and cumulative outflow during the prediction period.

6. The method according to claim 1, characterized in that, After obtaining the optimized outbound process, a two-level convergence verification step is also included: Inner layer check: Calculate the flow balance residual of the upstream and downstream connection section. If the residual exceeds the preset threshold, correct the downstream inflow process according to the residual and re-execute the evolution algorithm. Outer layer test: Calculate the rate of change of excess flood volume at the downstream flood control section between two adjacent iterations. If the rate of change is greater than the convergence threshold, correct the excess flood volume vector and return to the step of constructing dynamic risk potential index.

7. A two-layer optimization system for real-time flood control scheduling of a series of reservoirs, characterized in that, include: Memory, used to store computer programs; A processor for executing a computer program to implement the steps of the method as claimed in any one of claims 1 to 6.

Images