Water resource system risk scheduling method

Abstract

The application discloses a kind of water resource system risk scheduling methods, comprising: obtaining regional basic data, construct supply and demand pressure field and identify causal association, generate heterogeneous computing grid and space-time causal network;Based on the historical supply and demand sequence data of heterogeneous computing grid and the influence factor identified by space-time causal network, the tail threshold of water demand is dynamically determined using the cascade buffer pool mechanism, and the water demand sequence after regulation is output;Using the water demand sequence after regulation, the optimal transport method with tail shape is used to generate a set of hydrological scenarios, and the set of hydrological scenarios is simplified into a representative scenario subset;For the representative scenario subset, a distribution robust optimization model is constructed, and a water resources scheduling scheme and a risk assessment report are generated through fairness constraint projection and causal path analysis.The application can finely depict the space-time pressure, adaptively regulate extreme demand, and generate a scheduling scheme that takes into account robustness and fairness while preserving tail risk characteristics.

Description

TECHNICAL FIELD

[0001] The present application relates to a scheduling method, in particular a water resource system risk scheduling method. BACKGROUND

[0002] Water resources are the basic strategic resources for maintaining the health of regional ecosystems and supporting the sustainable development of economic society. However, the spatial and temporal distribution of water resources is extremely uneven, and global climate change has led to an increase in the frequency and intensity of extreme hydrological events (such as severe droughts and heavy rain floods), which has exacerbated the vulnerability and uncertainty of water resource systems. Therefore, how to achieve fine scheduling of water resources, ensure water supply safety, and reduce losses from extreme events under increasingly complex risk scenarios has become a major scientific problem and practical challenge in the field of water resources management.

[0003] At present, research on water resource scheduling has made certain progress. In terms of optimization models, the existing technology often uses linear programming, nonlinear programming, dynamic programming or multi-objective intelligent algorithms (such as genetic algorithms), and combines the large system decomposition and coordination theory to carry out multi-objective joint scheduling of reservoir groups, water diversion projects and local water sources within the region. In terms of uncertainty processing, the existing technology often uses stochastic programming or fuzzy programming methods, for example, a large number of hydrological scenarios are generated through Monte Carlo simulation, or a random sequence is generated based on historical runoff data using methods such as Markov chain, and is used as the input of the optimization model, in order to obtain a scheduling scheme that has a relatively good comprehensive performance under multiple possibilities. In terms of water demand management, the water demand process of each calculation unit is usually given based on historical water consumption quota, social and economic development prediction or simple trend extrapolation method, and is used as a rigid constraint condition of the scheduling model.

[0004] However, the above existing technology still has limitations in dealing with extreme risks and representing complex spatio-temporal coupling. In terms of uncertainty representation, the traditional stochastic simulation method (such as Monte Carlo sampling) is often based on standard statistical distribution (such as normal or lognormal distribution), which is easy to underestimate the probability and severity of tail events such as extreme drought. The scheduling scheme developed based on such smoothed scenario set has only limited robustness in dealing with regular fluctuations, and is easy to fail systematically under the impact of real extreme events. In terms of spatio-temporal feature description, the existing technology often uses uniform grid or fixed administrative and watershed unit, which is difficult to accurately capture the pressure gradient characteristics of areas with severe supply and demand contradictions (such as the surrounding area of water source and the urban-rural transition zone); they often regard the water demand process as a rigid exogenous variable, ignoring the dynamic causal feedback between water supply vulnerability and water demand urgency (for example, meteorological drought not only leads to water supply vulnerability, but also causes irrigation water demand urgency), leading to lagging and lack of flexibility in demand regulation measures, and making it difficult to achieve truly adaptive management. SUMMARY

[0005] The application aims to provide a water resource system risk scheduling method to solve the above problems in the prior art.

[0006] The technical scheme is a water resource system risk scheduling method based on total water consumption control, comprising the following steps:

[0007] Obtaining regional basic data, constructing supply-demand pressure fields and identifying causal relationships, generating heterogeneous computing grids and spatio-temporal causal networks;

[0008] Based on the historical supply-demand sequence data of the heterogeneous computing grids and the influence factors identified by the spatio-temporal causal networks, a step buffer pool mechanism is used to dynamically determine the water demand tail threshold, and the water demand sequence after regulation is output;

[0009] Using the water demand sequence after regulation, an optimal transport method with tail shape preservation is used to generate a hydrological scenario set, and the hydrological scenario set is simplified into a representative scenario subset;

[0010] For the representative scenario subset, a distribution robust optimization model is constructed, and a water resource scheduling scheme and a risk assessment report are generated through fairness constraint projection and causal path analysis.

[0011] The application can finely depict spatio-temporal pressure, adaptively regulate extreme demand, and generate a scheduling scheme that takes into account robustness and fairness while preserving tail risk characteristics. BRIEF DESCRIPTION OF DRAWINGS

[0012] Figure 1 The step flowchart of the water resource system risk scheduling method based on total water consumption control in the embodiment of the application.

[0013] Figure 2 The step flowchart of generating a hydrological scenario set in the embodiment of the application.

[0014] Figure 3 The step flowchart of constructing a distribution robust optimization model in the embodiment of the application.

[0015] Figure 4 The step flowchart of generating a heterogeneous computing grid in the embodiment of the application. DETAILED DESCRIPTION

[0016] In order to enable personnel in the art to better understand the application scheme, the technical solutions in the embodiments of the application will be described clearly and completely below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only a part of the embodiments of the application, not all. Based on the embodiments in the application, all other embodiments obtained by those of ordinary skill in the art without creative labor should be within the scope of protection of the application.

[0017] Embodiment 1: A water resource system risk scheduling method based on total water consumption control.

[0018] This embodiment details the overall framework of a water resource system risk scheduling method based on total water consumption control, as shown in Figure 1

[0019] Step 101, obtain regional basic data, construct supply and demand pressure field and identify causal relationship, generate heterogeneous computing grid and spatio-temporal causal network.

[0020] In this embodiment, the regional basic data can include but not limited to: hydro-meteorological observation data of the study area (such as historical rainfall, evaporation, runoff series), land use data, crop planting structure data, population and economic distribution data (such as industrial output value distribution, population density) and water conservancy engineering layout information (such as water source point location, available water supply, pipeline transportation capacity). The supply and demand pressure field is a data field used to characterize the spatial distribution of water resource supply and demand contradictions, which can be constructed by fusing the water supply potential field and the water demand intensity field. Among them, the water supply potential field can represent the water supply capacity of the water source point to each point in space, while the water demand intensity field represents the water demand of each point in space. The heterogeneous computing grid is a non-uniform computing unit generated according to the spatial distribution characteristics of the supply and demand pressure field, which aims to automatically encrypt the grid in the area where the supply and demand contradiction changes dramatically (such as the supply and demand boundary or the area with large pressure gradient), and to sparsify the grid in the area where the change is gentle, so as to improve the calculation efficiency while ensuring the calculation accuracy. The spatio-temporal causal network is used to identify the key factors and their transmission paths that affect water demand changes. For example, by analyzing the time delay response relationship between historical hydro-meteorological sequence and historical water demand sequence, the causal relationship representing influence intensity, time delay and persistence can be constructed, and the dominant influence factor set and causal transmission path diagram can be extracted.

[0021] Step 102, based on the historical supply and demand sequence data of the heterogeneous computing grid and the influence factors identified by the spatio-temporal causal network, the demand tail threshold is dynamically determined by using the cascade buffer pool mechanism, and the regulated water demand sequence is output.

[0022] In this embodiment, it is to realize the adaptive management of extreme demand. Through the cascade buffer pool mechanism, different levels of demand response can be activated according to different supply and demand tension. Specifically, for each heterogeneous computing grid, according to its historical water demand sequence (such as sequence standard deviation σ i ​) Set a multi-level static buffer pool threshold system. Use the impact factors identified by the spatiotemporal causal network (such as air temperature, precipitation deficit factor, crop growth period factor) and historical supply and demand data (such as the supply water vulnerability index SVI calculated by the standardized precipitation index SPI conversion method) to build a dynamic adjustment factor for real-time adjustment of the static buffer pool threshold, and obtain a dynamic trigger threshold. When the real-time supply and demand gap reaches the dynamic trigger threshold, the corresponding level of buffer pool is activated, and the preset release function (such as the exponential release function) is used to determine the water demand reduction ratio, and the original water demand is reduced to generate the regulated water demand sequence.

[0023] Step 103, using the regulated water demand sequence, generate a hydrological scenario set using the tail-preserving optimal transport method, and simplify the hydrological scenario set to a representative scenario subset.

[0024] In this embodiment, reasonable uncertainty representation is provided for subsequent robust decision-making. Conventional scenario generation methods may not be able to preserve the statistical characteristics of extreme events (such as extreme drought), and tail-preserving aims to solve this problem. Specifically, the tail-preserving optimal transport method can separate the historical data of the regulated water demand sequence into a main body area and a tail area (for example, by fitting the tail data with a generalized Pareto distribution GPD to define it). When building the optimal transport model, by introducing a cross-zone penalty term (for example, based on an indicator function to prevent tail area samples from being transmitted to the main body area) and a quantile constraint term (to maintain the tail distribution shape), a hydrological scenario set that preserves the characteristics of extreme events is generated. Simplification means that a large-scale scenario set is compressed into a computationally feasible representative scenario subset. For example, a facility location optimization model is constructed, and the objective of the model is to identify a representative scenario subset from the hydrological scenario set, and the model also satisfies the preset tail quota constraint and the spatiotemporal block coverage constraint. Among them, two constraints are applied in the model: tail quota constraint (to ensure that the representative scenario subset contains a predetermined proportion of tail scenarios) and spatiotemporal block coverage constraint (to ensure that each spatiotemporal block is covered by a predetermined number of scenarios in the representative scenario subset). Solve the facility location optimization model with constraints to obtain the representative scenario subset.

[0025] Step 104, for the representative scenario subset, build a distributionally robust optimization model, generate a water resources scheduling scheme and a risk assessment report through fairness constraint projection and causal path analysis.

[0026] In this embodiment, the distributionally robust optimization model (DRO) does not rely on the exact probability distribution of the scenario subsets, but seeks to minimize the risk under the worst-case scenario (e.g., minimize the conditional value at risk CVaR) within an uncertainty set constructed around the nominal distribution (e.g., a Wasserstein ball uncertainty set). Solving the DRO model can obtain a preliminary scheduling scheme. The fairness constraint projection is used to modify the preliminary scheme to meet the social fairness requirement. For example, projecting the preliminary scheme onto a predefined fairness polyhedron (e.g., constraining the difference in water use satisfaction rates of different regions within a tolerance ε fair The scheduling scheme can include the water supply allocation amount {x i} of each grid and the water transfer path flow {y ij}. The causal path analysis is used to generate a risk assessment report. It can combine the spatiotemporal causal network and the final water resource scheduling scheme to calculate the conduction intensity of risk (such as water shortage risk) on the causal path, evaluate the propagation probability of water shortage risk in the network, and identify the key risk conduction path and control node.

[0027] Embodiment 2: Preferred implementation of supply-demand pressure field and heterogeneous grid generation. This embodiment aims to generate a specific and preferred refinement scheme of the heterogeneous computing grid, as shown in Figure 4

[0028] Step 201, based on the regional basic data, construct the supply-demand pressure field.

[0029] In this embodiment, the supply-demand pressure field P(x, y) is a data field representing the spatial distribution of water resource supply-demand contradictions. The geometric features of this pressure field will be used later to guide grid generation.

[0030] Step 202, extract the geometric features of the supply-demand pressure field.

[0031] In this embodiment, the way to extract geometric features is as follows: calculate the Hessian matrix H of the supply-demand pressure field P(x, y) at any spatial point (x, y). The Hessian matrix H is a second-order partial derivative matrix, specifically H = [[Ψ²P / Ψx², Ψ²P / ΨxΨy], [Ψ²P / ΨxΨy, Ψ²P / Ψy²]]. Where Ψ is the partial derivative. This matrix can describe the local curvature of the supply-demand pressure field P(x, y), i.e., the degree and direction of change. After obtaining the Hessian matrix H, further solve the eigenvalues λ1 and λ2 of the Hessian matrix. The eigenvalues λ1 and λ2 can reflect the bending degree of the pressure field at this point in two principal directions.

[0032] Step 203, construct a grid density control function based on the geometric features.

[0033] ​In this embodiment, the Hessian matrix of the supply-demand pressure field at any spatial point is calculated, the eigenvalues λ1 and λ2 of the Hessian matrix are solved, the absolute values of the eigenvalues λ1 and λ2 are summed to obtain |λ1|+|λ2|, and |λ1|+|λ2| is taken as the grid density control function to represent the variation intensity of the supply-demand pressure field. The grid density control function is constructed based on the eigenvalues of the Hessian matrix. Specifically, the absolute values of the eigenvalues λ1 and λ2 are summed to obtain |λ1|+|λ2|. When the pressure field varies sharply (for example, at the boundary of the supply-demand contradiction or the pressure saddle point), the local curvature is large, resulting in a large sum of the absolute values of the eigenvalues |λ1|+|λ2|; on the contrary, in the pressure gentle area, the curvature is small, and the value is also small. Therefore, the function value can effectively represent the variation intensity of the supply-demand pressure field.

[0034] Step 204, applying the grid density control function, densifying the grid in the area where the supply-demand pressure field varies sharply to generate a heterogeneous computing grid.

[0035] In this embodiment, when generating the grid (for example, using constrained Delaunay triangulation or other adaptive grid generation algorithm), the grid density control function |λ1|+|λ2| is taken as the input. This function will guide the grid generation algorithm to generate grid cells with higher density and smaller size in the area where |λ1|+|λ2| is larger (i.e., the area where the pressure varies sharply), and to generate grid cells with lower density and larger size in the area where |λ1|+|λ2| is smaller (i.e., the gentle area). Through the generated heterogeneous computing grid {G i}, the spatial characteristics of the supply-demand pressure field can be automatically adapted, and the optimal allocation of computing resources can be achieved.

[0036] Embodiment 3: An alternative embodiment of generating a heterogeneous computing grid. This embodiment aims to support the technical features related to generating a grid.

[0037] In this embodiment, the water scarcity field S(x, y) is calculated based on the supply-demand pressure field P(x, y). The water scarcity field can be defined as the ratio of the water demand intensity field to the water supply potential field, that is, S(x, y)=I demand (x, y) / P supply (x, y), which directly reflects the local water resource shortage. On this basis, gradient analysis is performed on the water scarcity field S(x, y), for example, the gradient field ▽S is calculated. Through the gradient tracking algorithm, the ridge lines and valley lines of the gradient field ▽S can be identified. Among them, the ridge line (local maximum gradient connection line) usually corresponds to the boundary area where the supply-demand contradiction is the most acute; and the valley line (local minimum gradient connection line) may represent the main channel of water resource transport or the area where the supply-demand is relatively moderate.

[0038] When generating the mesh, ridges and valleys are used as key geometric feature lines, and the locations of hydraulic engineering nodes (such as reservoirs, pumping stations, and main canal heads) are used as mandatory retention points. A heterogeneous computational mesh {G} is generated based on, for example, constrained Delaunay triangulation, using the mandatory retention points and feature lines. i The generated grid's boundaries can match the natural boundaries (ridges, valleys) of supply and demand imbalances, thus achieving adaptive grid partitioning to the spatial distribution characteristics of water resources.

[0039] Example 4: Preferred Implementation of Directional Modeling and Synthesis of Supply and Demand Pressure Fields. This example is a preferred and more detailed implementation of the steps for constructing supply and demand pressure fields, with particular focus on the physical meaning and directional modeling in the process of constructing the supply and demand pressure fields.

[0040] In this embodiment, a water supply potential energy field P is constructed. supply To reflect the potential energy propagation characteristics and directionality of water resources (especially surface runoff) under the influence of gravity and topography when (x, y), a logarithmic spiral field superposition method incorporating azimuth cosine modulation is preferred. Specifically, for the i-th water source point within the study area, its available water volume is V. i The radius of influence is r i The potential energy contribution of water supply generated at the spatial point (x, y) can be calculated using the following formula: P single_i (x, y) = V i ×log(1+e (-di / ri) )×cos(θ i / 2π). Where d i Let (x, y) be the distance from the spatial point (x, y) to the i-th water source point, log(1+e (-di / ri) θ is a logarithmic spiral decay function, characterizing the decrease in potential energy with distance; i Let (x, y) be the azimuth angle of a spatial point (x, y) relative to the i-th water source point, and let cos(θ) be the azimuth angle of the point (x, y) relative to the i-th water source point. i The term ( / 2π) is the azimuth cosine modulation term, used to introduce directional effects. The final water supply potential field P supply (x, y) is the sum of the potential energy contributions of all water source points in the region, i.e., P supply (x, y) = Σ[P single_i (x, y)].

[0041] Constructing a water demand intensity field I demand When (x, y) is defined, a weighted overlay method is preferred. Specifically, different types of basic water demand data are obtained, such as crop planting distribution layers for agricultural land. layer (x, y) and the comprehensive crop water requirement coefficient k crop Grid-based industrial land use data for GDPind (x, y) and GDP unit and the gridded population density distribution data Pop(x, y) and water per capita of residential land per_capita . The water intensity coefficients of various types of land are determined according to the land use type, for example, agriculture k agr , industry k ind , residence k res . The final water demand intensity field I demand (x, y) is constructed by weighted superposition as follows: I demand (x, y) = k agr × crop layer (x, y) × k crop + k ind × GDP ind (x, y) / GDP unit + k res × Pop(x, y) × water per_capita After obtaining the water supply potential field P supply (x, y) and the water demand intensity field I demand (x, y) respectively, the two are combined into a unified supply and demand pressure field P(x, y). Preferably, the supply and demand ratio method is used to calculate the initial pressure P ratio (x, y) = I demand (x, y) / (P supply (x, y) + ε small ), where ε small is a very small positive number (for example, 10 -6 ), used to prevent the denominator from being zero. In order to eliminate local outliers and make the field data smoother, a spatial smoothing kernel function K smooth (for example, Gaussian kernel) can be further applied to P ratio (x, y) for convolution operation, and P smooth (x, y) = P ratio ⊕ K smooth , where ⊕ represents the convolution operator. Finally, the supply and demand pressure field P(x, y) = log(1 + P smooth (x, y)) is generated by logarithmic transformation, which helps to make the data distribution closer to normal, enhance the numerical stability, and facilitate the subsequent gradient or Hessian matrix calculation.

[0042] Example 5: Detailed implementation of causal identification and dynamic demand regulation. This embodiment aims to provide a complete and in-depth refinement scheme for dynamic demand regulation, and details the complete technical chain from spatiotemporal causal network identification to gradient buffer pool dynamic regulation.

[0043] Step 501, identify causal relationships, generate spatio-temporal causal network.

[0044] In this embodiment, the construction of the spatio-temporal causal network is performed for each generated heterogeneous computing grid G i (or a representative region composed of several grids) based on historical time series data. The historical hydro-meteorological time series and historical water demand time series read by the heterogeneous computing grid are used. Specifically, the historical water demand sequence W i (t), the historical water supply sequence W demand (t), and the related regional historical hydro-meteorological time series (such as rainfall R(t), evaporation E(t), runoff Q(t)) of the grid G supply are combined into an input matrix U(t).

[0045] Specifically, an echo state network is constructed and impulse response analysis is performed to obtain response parameters representing influence intensity, time delay, and persistence; the process of constructing the echo state network is to initialize the internal weight matrix W res of the reservoir, which is usually a sparse random matrix (e.g., sparsity of 10%), and set its spectral radius ρ(W res ) near 1.0 (e.g., 0.95) to ensure that the network has echo state properties. Iterative calculation is performed through the state update equation x(t+1)=tanh(W in ×U(t)+W res ×x(t)) to generate the state sequence x(t) of the reservoir, where W in is the input weight matrix and tanh is the hyperbolic tangent activation function. The process of performing impulse response analysis is to apply a unit pulse δ(t) to each input variable (e.g., rainfall R(t)) at a certain time, and record the response trajectory x res ponse(t) of the reservoir state sequence x(t). A pre-set response function model, such as the double exponential model R(t)=A×exp(-t / τ)×(1-exp(-t / T p )), is used to fit the response trajectory. The response parameters: response amplitude A (representing influence intensity), decay time constant τ (representing influence persistence), and peak time T p (representing response delay) are estimated by fitting.

[0046] The intensity, time delay, and persistence dimensions of the response parameters are constructed into a three-dimensional causal tensor. The response parameters (intensity A, time delay T p , and persistence τ) between all variable pairs are constructed into a three-dimensional causal tensor T causal[i, j, k], where i and j are variable indices, and k is a feature index (e.g. k = 1 for intensity, k = 2 for time lag, k = 3 for persistence). Perform tensor decomposition (e.g. CP decomposition or Tucker decomposition, T causal =∑(λ r ×u r ⊕v r ⊕w r ) on the three-dimensional causal tensor, extract the dominant influence factor set and causal conduction path graph, and construct the spatio-temporal causal network. By extracting the dominant components whose cumulative contribution rate exceeds a predetermined threshold (e.g. 85%), the set of influence factors F (e.g. air temperature, precipitation deficit) that play a dominant role in the change of water demand and the causal conduction path graph G causal representing how these factors conduct their influence can be identified.

[0047] Step 502: Based on the historical supply-demand sequence data of the heterogeneous computing grid and the influence factors identified by the spatio-temporal causal network, the step buffer pool mechanism is used to dynamically determine the water demand tail threshold.

[0048] In this embodiment, a static step buffer pool is set. For each heterogeneous computing grid G i , its historical water demand sequence W demand (t) is read and its standard deviation σ i is calculated. Preferably, a five-level buffer pool system is constructed. For example, the first three levels (j = 1, 2, 3) use an exponential growth pattern, and their threshold values are set to Buffer j = 0.3σ i × 1.5 j ; the last two levels (j = 4, 5) use a linear growth pattern, and their threshold values are set to Buffer j = 1.5σ i × (j-2). This set of {Buffer ij} constitutes the baseline threshold.

[0049] Calculate the double vulnerability index. On the one hand, based on the water supply sequence W supply (t) in the historical supply-demand sequence data, the water supply vulnerability index SVI(t) is calculated by the standardized precipitation index SPI conversion method. The calculation process is, for example: SVI(t) = Φ (-1) [(rank(W supply (t))-0.375) / (n+0.25)], where rank is the ranking function, n is the sequence length, and Φ (-1) is the inverse function of the standard normal distribution. The lower the SVI(t) value, the more vulnerable the water supply. On the other hand, based on the dominant influence factor set F, the entropy weight method is used to determine the weight w k of the factor F kand by weighted summation and normalization of the water demand, a demand urgency index DUI(t) is calculated, for example DUI(t) =∑(w k × F k (t)) x normalize(W demand (t)). The higher the DUI(t) value, the more urgent the demand.

[0050] Fusion of indices and adjustment of threshold. SVI and DUI are fused into a unified linear factor for adjusting a baseline threshold {Buffer ij}. The unified linear factor is used to adjust the water demand tailing threshold to obtain a dynamic triggering threshold. Preferably, SVI and DUI are fused into a dynamic adjustment factor in a weighted linear combination manner to obtain the dynamic triggering threshold Trigger ij (t). The specific formula can be: Trigger ij (t) = Buffer ij x (1 + a x SVI(t) + b x DUI(t)). Wherein a and b are adjustable adjustment coefficients, a is usually negative (or written as 1-a'xSVI(t) in the formula, where a' is positive) to lower the threshold when the water supply is vulnerable (SVI is negative), and b is usually positive to lower the threshold (i.e. more likely to trigger reduction) when the demand is urgent (DUI is positive).

[0051] Step 503, using the dynamic triggering threshold to drive the stepped buffer pool mechanism to generate the regulated water demand sequence.

[0052] In this embodiment, hysteresis control is applied to the dynamic triggering threshold to determine the stable activation level. Hysteresis control is used to determine the activation level. Specifically, when the real-time supply-demand difference (e.g. real-time water demand prediction value) of the grid reaches the dynamic triggering threshold Trigger ij (t), the corresponding level j is activated. In order to avoid frequent switching of the regulation level due to frequent fluctuations of the supply-demand difference near the threshold, hysteresis control mechanism is preferably used. Specifically, a hysteresis interval is set for each level j's threshold Trigger ij , for example [T down_j , T up_j ], where T down_j = Trigger ij x 0.9, T up_j = Trigger ij x 1.1. Only when the supply-demand difference rises above T up_j , it is upgraded to a higher level, and only when it falls below T down_j , it is downgraded to a lower level. Through this mechanism, a stable activation level level stable can be determined.

[0053] The reduction ratio is determined by applying an exponential release function. Based on a representative subset of scenarios, the reduction ratio is determined according to a stable activation level j (i.e., level). stable The water demand reduction ratio is determined using an exponential release function. Preferably, a release function with activation level j as the exponent is used, so that the reduction ratio increases exponentially with level j, achieving greater reduction at higher levels. For example, the release (reduction) coefficient can be defined as: factor =1-0.1×(1.5 j When j=0 (inactive), release factor =1.0 (no reduction); when j=1, release factor =0.85 (reduction of 15%); when j=5 (highest level), release factor ≈0.24 (reduction of 76%).

[0054] Generate the post-regulation water demand sequence. Based on the water demand reduction ratio (i.e., release ratio). factor ), reduce original water demand W demand (t), generating the water demand sequence W after regulation. truncated (t)=W demand (t)×release factor .

[0055] In a preferred embodiment, the adjustment coefficients α and β need to be calibrated. A preferred calibration method is to use adjoint sensitivity calibration. The post-modulation water requirement sequence W is then constructed. truncated The loss function Loss is defined as the degree of matching between (t) and the actual available water quantity. The gradients ΨLoss / Ψα and ΨLoss / Ψβ of the loss function Loss with respect to parameters α and β are calculated using automatic differentiation. Under the condition of satisfying the mass conservation constraint, α and β are updated and optimized using the gradient descent method.

[0056] As an alternative to the aforementioned sensitivity calibration, a historical backtesting method can also be used. A backtesting dataset containing historical SVI, DUI, and actual water shortage events is constructed. Search intervals for parameters α and β are set (e.g., α∈[0.1, 0.5], β∈[0.1, 0.4]). The parameter combinations are traversed using a grid search method, and a fitness function is defined (e.g., fitness = accuracy × (1 - false)). alarm_rate ), where accuracy is the prediction accuracy, false alarm_rate The parameter performance is evaluated using the false alarm rate, and the parameter combination with the highest fitness is selected as the final optimization adjustment coefficient α. opt β opt .

[0057] Embodiment 6: Detailed implementation of tail-shaped scenario generation and reduction. This embodiment aims at a complete, in-depth refinement scheme for scenario generation and optimization, as shown in Figure 2

[0058] Step 601: Separate the regulated water demand sequence into the main body and tail regions.

[0059] In this embodiment, tail samples are extracted from the regulated water demand sequence. Specifically, from the complete historical data of the output regulated water demand sequence W truncated (t), the high quantile (e.g., 95% quantile) is calculated as the tail threshold Q 95 . The data points in the sequence whose values exceed Q 95 are extracted as the tail sample set W tail , and the remaining data points are the main body sample set W body .

[0060] On this basis, generalized Pareto distribution GPD fitting is performed on the tail samples, i.e., generalized Pareto distribution GPD fitting is performed on the tail sample set W tail . Among them, the tail samples (extreme events) are usually few in number but significant in impact, and direct sampling is difficult to capture their true distribution, while GPD is a standard distribution in extreme value theory for describing data above the threshold. Preferably, the maximum likelihood method is used to estimate the shape parameter ξ and the scale parameter σ gpd of GPD. For example, by optimizing the likelihood function L(ξ, σ gpd ) = -n tail ×log(σ gpd )-(1+1 / ξ)×Σlog(1+ξ×excess i / σ gpd ), where n tail is the number of tail samples, and excess i is the amount of the i-th tail sample exceeding the threshold Q 95 . Based on the GPD fitting result, i.e., the mathematical representation F tail (x) of the tail region is obtained, which defines the tail region.

[0061] Step 602: Construct a partitioned optimal transport model.

[0062] In this embodiment, after separating the main body and tail regions, a partitioned optimal transport model is constructed for generating a new set of hydrological scenarios. The cost function Cost' of the model is preferably designed to include three components: Cost'(i, j) = C base (i, j) + γ tail ×I(i∈tail, j∈body) + λ×|Q target -Q actual ​| Wherein, the basic transmission cost C base (i, j) is the standard optimal transmission cost, for example, the square of the distance between sample i and sample j, C base (i, j) = |x i -x j | 2 .

[0063] The cross-zone penalty term γ tail ×I(i∈tail, j∈body) is an important aspect of the present application, which is used to achieve tail preservation. This term is preferably an asymmetric cross-zone penalty term based on an indicator function. Wherein, I(i∈tail, j∈body) is an indicator function, which takes the value of 1 when the transmission tries to map the original tail sample i to the target body zone location j, otherwise 0. γ tail is a large penalty coefficient (for example, γ tail =10). The design technique is to increase the cost of degrading extreme event samples to general event samples, preventing tail zone samples from transmitting to the body zone.

[0064] The quantile constraint term λ×|Q target -Q actual | is also a key to achieve tail preservation. Wherein, Q target is the target quantile of the original data (for example, the 95% and 99% quantile of GPD fitting), Q actual is the actual quantile of the generated scenario set in the transmission process, and λ is the penalty weight of this constraint (for example, λ=100). This term is used to ensure that the generated scenario set is consistent with the original data in the tail distribution shape.

[0065] Step 603, solving the partition optimal transport model to generate a hydrological scenario set.

[0066] In this embodiment, after constructing the above cost function, the optimal transport model is solved. For example, the Kantorovich problem is solved iteratively by Sinkhorn algorithm to obtain the optimal transport scheme (i.e. transmission probability matrix) π ij . According to the transmission scheme π and the tail distribution function F tail (for example, generating tail by GPD random sampling), the hydrological scenario set {S k} and its corresponding occurrence probability {p k} are combined to generate.

[0067] Step 604, refining the hydrological scenario set into a representative scenario subset.

[0068] In the preferred embodiment, before performing this step, the generated hydrological scenario set {S kThe execution time window pattern matching step. This step aims to enhance the temporal consistency of the scenarios. For example, identify extreme event flags extreme flag in the scenarios (e.g. flag = 1 if the water demand of a day exceeds a threshold) and use adaptive weights w(t) = w normal × (1 - extreme flag ) + w extreme × extreme flag for dynamic time warping matching, where w extreme > w normal , indicating more focus on aligning the extreme event periods when matching.

[0069] The pruning process is as follows: divide the scenarios in the hydrological scenario set into space-time blocks according to the spatial grid and time period. For example, divide according to the n grid spatial grid generated and T / 30 time periods in months. Preferably, extract a feature vector for each space-time block block ij , for example vector ij = [mean, std, skew, kurt, max, min], i.e. the mean, standard deviation, skewness, kurtosis, maximum and minimum of the data in the space-time block.

[0070] Construct a facility siting optimization model. The goal of this model is to identify a representative scenario subset {c j} from the large-scale original scenario set. The objective function of this model can be set to minimize the sum of the distances of all original scenario blocks to their nearest representative scenario centers, for example minΣ i (min j (d(block i , c j ))), where d is the Euclidean distance between block feature vectors and c j is the representative scenario center to be selected.

[0071] Apply a dual quota constraint to the model. The first constraint is the tail quota constraint, which is used to ensure that the representative scenario subset contains a predetermined proportion of tail scenarios. Its form can be: |N tail_selected / N tail_total - ρ tail | ≤ ε q . Where N tail_selected is the number of tail scenarios in the selected representative scenarios, N tail_total is the total number of tail scenarios in the original scenario set, ρ tail is the target tail proportion set by the decision maker (e.g. 0.1 to 0.2, indicating that 10% to 20% of the extreme scenarios are desired in the subset), and ε qThe tolerance is set to (e.g., 0.02). The second constraint is the spatiotemporal block coverage constraint, which ensures that each of the aforementioned spatiotemporal blocks is covered by a predetermined number (e.g., k) of representative scene subsets. min The scene subset is covered by the condition =3). This constraint guarantees the spatiotemporal diversity of the scene subset.

[0072] Solve a constrained facility location optimization model. For example, solve it using a mixed-integer programming solver to obtain the selected representative scenario index set I. selected That is, the representative scene subset {S* m}. And recalculate its aggregation probability weights {p* m For example, each representative scenario S* m The original scene represented (i.e., distance S*) m The latest collection of original scenes (Cluster) m The original probability p) k Summing yields p* m =Σ(p k forkinCluster m ).

[0073] Example 7: Detailed Implementation of Robust Optimization and Fair Scheduling. This example aims to provide a complete and in-depth refinement of the robust decision-making scheme.

[0074] Step 701, construct the sub-Bruker optimization model, such as Figure 3 As shown:

[0075] This embodiment is based on a representative scenario subset {S* m} and its aggregate probability weights {p* m}. Calculate the first moment (mean vector) μ and the second moment (covariance matrix) Σ. cov Specifically: μ=Σ(p* m ×S* m ) and Σ cov =Σ(p* m ×(S* m -μ)(S* m -μ) T ).

[0076] Using the first and second moments and a pre-defined Wasserstein radius, a Wasserstein spherical uncertainty set is constructed. The uncertainty set U is the set containing all plausible probability distributions P. It is defined as: [The uncertainty set U is related to the nominal distribution P]. nominal (by {p*) m The 2-Wasserstein distance W2 (as defined) does not exceed the preset radius ε. w (e.g., ε) w=0.1, this value can be determined by cross-validation based on the degree of confidence in the nominal distribution, i.e., U={P:W2(P, P... nominal )≤ε w In some alternative implementations, this uncertain set can also be constructed using moment constraints, such as U={P:|E}. P [S]-μ|≤δ1, E P [(S-μ)(S-μ) T ]≤Σ cov +δ2×I}, where δ1 and δ2 are the tolerance parameters of the moment constraint.

[0077] The objective function of the subbloc bar optimization model is constructed. The subbloc bar optimization model is constructed to minimize the conditional risk value (CVaR) under the worst-case scenario within the Wasserstein spherical uncertainty set U. The conditional risk value (CVaR) is an indicator of tail risk. Preferably, for ease of calculation, its linearized expression is used for modeling. An auxiliary variable η is introduced to represent the value of risk (VaR). α (where α is the confidence level, for example, 0.95), then CVaR α The objective function can be modeled as: min[η+(1 / (1-α))×E P [max(0, cost(x, S)-η)]]. To further transform it into a linear programming problem, the expectation term can be expanded into a weighted sum based on representative scenarios: min[η+(1 / (1-α))×Σ] m (p* m ×z m A new auxiliary variable z was introduced. m And impose constraints: z m ≥cost(x, S*) m )-η and z m ≥0. cost(x, S*) m ) is in a specific scenario S* m Given a scheduling decision x (including the water supply x of each grid), i and the flow rate y of the water diversion path ij The total cost function generated when Σ is used. This cost function preferably includes: Σ i (c supply_i ×x i (Local water supply cost), Σ{i,j}(c trans_ij ×y ij (Cost of inter-regional water transfer) and the key water shortage penalty (penalty×max(0, demand)). i mx i (Water shortage penalty cost, including demand) i (where m represents the requirements of mesh i in scene m).

[0078] Step 702, solve the distribution robust optimization model to obtain a preliminary scheduling scheme.

[0079] In this embodiment, a series of physical constraints must be imposed when solving the optimization model. These constraints preferably include:

[0080] Water balance constraint: ensure that the inflow (import from local source source i and other grids Σ j (y j i)) of each grid i is equal to the outflow (supply to local x i and Σ j (y ij ) to other grids), that is, Σ j (y j i)+source i =Σ j (y ij )+x i .

[0081] Transport capacity constraint: ensure that the water transfer path flow y ij between any two grids i and j does not exceed its physical transport capacity upper limit capacity ij , that is, 0≤y ij ≤capacity ij .

[0082] Minimum guarantee constraint: ensure that the water supply x i of each grid i at least meets its basic demand (such as domestic water), that is, x i ≥γ min ×demand i , where γ min is the minimum guarantee coefficient (for example, 0.7). By solving the above optimization model with objective function and constraints by a large-scale linear programming solver such as the interior point method, a preliminary scheduling scheme X preliminary ={x i , y ij} can be obtained.

[0083] Step 703, project the preliminary scheduling scheme onto the predefined fairness polyhedron to perform fairness constraint projection to generate a water resources scheduling scheme.

[0084] In this embodiment, the preliminary scheduling scheme X preliminary is the optimal solution under the worst-case scenario, but may not meet the fairness between regions. Therefore, fairness constraint projection needs to be performed. Define the fairness polyhedron P fairThis is the feasible region for characterizing the fairness of the scheduling scheme. Preferably, it is defined as the set of water supply rates (i.e. water supply amounts x i / demand i ) that satisfy the difference between the supply rates of any two grids i and j, not exceeding a pre-set fairness tolerance ε fair (e.g. set to 0.05 to 0.15 according to regulatory requirements). Its mathematical expression is: P fair ={x:|x i / demand i -x j / demand j |≤ε fair , for all i, j combinations}. The projection is performed. This projection is a process of finding the solution that is closest to the preliminary scheduling scheme X fair within the fairness polyhedron P preliminary . Preferably, the Bregman projection algorithm is adopted, i.e. solving the optimization problem: minD φ (x, X preliminary ), subject to x e P fair , where D φ is the Bregman divergence, used to measure the distance between two schemes. This projection optimization problem can be solved efficiently by, for example, the alternating direction method of multipliers, eventually generating the water resources scheduling scheme X final ={x* i , y* ij}.

[0085] Embodiment 8: Risk assessment implementation based on causal paths. This embodiment aims to generate a specific and preferred refinement scheme for generating a risk assessment report, relying on the generated heterogeneous computing grid.

[0086] The premise of this embodiment is that the obtained heterogeneous computing grid {G i}, the spatiotemporal causal network (in particular the causal conduction path graph G causal ), and the water resources scheduling scheme X final (in particular the water diversion path flow {y* ij}) have been obtained.

[0087] Step 801, utilize the causal conduction path graph contained in the spatiotemporal causal network, in combination with the water diversion flow in the water resources scheduling scheme.

[0088] In this embodiment, the causal conduction path graph G causal is a directed graph G = (V, E), where the nodes V represent hydro-meteorological or water demand variables, and the edges E represent causal relationships between variables. Each edge E has a weight, which can be referred to as the causal correlation strength (correlationedge ), which is determined by methods such as echo state network impulse response analysis. The water transfer flow {y ij} in the water resources scheduling scheme represents the actual physical water connection between grids, which can be regarded as the physical connection strength (flow edge ).

[0089] Step 802, calculate the transmission strength of the risk on the causal conduction path diagram.

[0090] In this embodiment, the propagation of water shortage risk depends on both the causal correlation between variables (for example, upstream rainfall deficit can lead to downstream water shortage) and whether there is a physical water transfer connection (if there is no water transfer connection between two regions, the risk transmission is blocked). Preferably, the transmission strength I path of a complete risk conduction path (for example: A region rainfall deficit -> B reservoir low water level -> C water transfer path -> D grid water shortage) is defined as the product of the weights of all related edges on the path (including the causal correlation strength and the physical connection strength), that is, I path =Π(correlation edge ×flow edge ). If the flow edge of any section (for example, C water transfer path) in the path is zero, the overall transmission strength I path of the path is also zero.

[0091] Step 803, evaluate the propagation probability of the water shortage risk in the network, and generate a risk assessment report.

[0092] In this embodiment, after calculating the transmission strength I i of all critical paths from the risk source to the target grid G path , the final water shortage risk probability P i (G risk ) of the grid G i can be evaluated. Preferably, the probability P risk (G i ) can be calculated as the weighted sum of the transmission strength of all paths pointing to the grid and the initial probability of the risk source P source , that is, P risk (G i )=Σ(I path ×P source ). Wherein P source is the initial probability of the risk (for example, severe rainfall deficit, severe reservoir depletion) occurring in the risk source node (for example, the uppermost rainfall factor, the key reservoir), which can be pre-set according to historical statistical data or expert evaluation. The risk assessment report can finally include: the water shortage risk probability {Priski a spatial distribution map, a ranked list of critical risk transmission paths, and identification of nodes (e.g. a certain conveyance path or a certain upstream reservoir) that play a key control role in the risk propagation.

[0093] Embodiment 9: Risk assessment alternative based on Monte Carlo simulation. This embodiment provides another alternative implementation of generating a risk assessment report, which does not directly rely on the causal path analysis.

[0094] This embodiment is based on the final dispatch scheme X final and a representative subset of scenarios {S m} and their aggregated probabilities {p m}. This method employs a scenario analysis or Monte Carlo simulation approach. Specifically, the final dispatch scheme X final (which is robust and aims to balance the performance across all scenarios) is back- traced to each representative scenario S m and run, and the actual supply i , m and demand i , m for each grid i in that particular scenario m is calculated.

[0095] From this, the deficit i , m = max(0, demand i , m - supply i , m) for each grid i in each scenario m can be calculated. By statistically weighting all representative scenarios {S m} and their corresponding aggregated probabilities {p m}, the expected deficit (E[deficit i ]), the deficit frequency (i.e. the sum of probabilities p i for all scenarios where deficit m , m > 0), the average deficit depth (i.e. the average deficit amount when deficit occurs), and the maximum deficit (i.e. the deficit amount in the worst-case scenario) for each grid i can be calculated. The risk assessment report contains the above series of statistical indicators to visually demonstrate the comprehensive performance and potential risks of the water resources dispatch scheme in dealing with uncertainties.

[0096] In summary, the present application provides a method as follows: S1, obtaining hydro-meteorological basic data, constructing supply and demand pressure field and generating heterogeneous calculation grid by using Hessian matrix eigenvalue, and constructing three-dimensional causal tensor by using echo state network impulse response analysis to identify spatio-temporal causal network; S2, based on the influence factors identified by the grid historical data and the spatio-temporal causal network, the supply water vulnerability index and the demand urgency index are calculated respectively, the linear factor is used to fuse and adjust the reference threshold of the cascade buffer pool, the dynamic trigger threshold is obtained, and the regulated water demand sequence is output through the hysteresis control and the index release function; S3, using the regulated sequence, an optimal transport method (including an asymmetric cross-zone penalty term and a quantile constraint) is used to generate a hydrological scenario set that retains the characteristics of extreme events, and the set is simplified into a representative scenario subset through a facility location model based on double quota constraints (tail quota and block coverage); S4, for the representative scenario subset, a distribution robust optimization model with conditional risk value as the target is constructed, the preliminary scheme is solved, the fairness constraint projection is performed, the final scheduling scheme is generated, and the risk is evaluated in combination with the causal path diagram.

[0097] In view of the problem that the traditional random simulation underestimates the tail risk and leads to fragile decision-making, the present application provides a tail-preserving scenario generation and compression mechanism. By performing generalized Pareto distribution fitting on the tail of the historical sequence, the statistical characteristics of extreme events are accurately captured in mathematics, rather than relying on rough assumptions of standard distribution. Further, in the construction of the optimal transport model, an asymmetric cross-zone penalty term (preventing tail samples from being transmitted to the main zone) and a quantile constraint term (maintaining the tail distribution shape) based on the indicator function are introduced. This design forces to retain the rarity and severity of extreme events, preventing them from being smoothed out in the generation process. In subsequent scenario compression, by setting the tail quota constraint, it is ensured that these extreme scenarios crucial to decision-making are retained in the representative subset, providing high-risk and high-accuracy input for subsequent distribution robust optimization, and ensuring the effectiveness of the scheduling scheme under real extreme events.

[0098] To solve the problem that traditional methods cannot finely depict the supply-demand pressure gradient and ignore the causality feedback, the present invention solves the problem from two aspects: spatiotemporal feature construction and demand dynamic regulation. In space, the uniform grid is abandoned, and an adaptive grid generation method based on the eigenvalues of the Hessian matrix of the supply-demand pressure field is adopted. The sum of the absolute values of the eigenvalues (|λ1|+|λ2|) is used as the grid density control function, so that the density of the calculation grid is automatically matched with the intensity of the supply-demand contradiction (pressure gradient), and the fine depiction of the key area is realized. In time and causality, the three-dimensional causality tensor is constructed by using the impulse response analysis of the echo state network, and the intensity, time delay and persistence of the influence of hydro-meteorological factors on water demand are deeply excavated. These forward-looking causal factors are calculated as the demand urgency index, which is combined with the water supply vulnerability index representing the current state to construct a dynamic trigger threshold. This mechanism makes the water demand management no longer a passive and rigid constraint, but an elastic and adaptive regulation according to the vulnerability of supply and demand and the causal premonition, and realizes the hierarchical and accurate reduction through the index release function.

Claims

1. A risk scheduling method for a water resource system, characterized in that, include: Acquire regional basic data, construct supply and demand pressure fields and identify causal relationships, and generate heterogeneous computing grids and spatiotemporal causal networks; Based on historical supply and demand sequence data and spatiotemporal causal network identification of heterogeneous computing grids, a cascaded buffer pool mechanism is used to dynamically determine the water demand truncation threshold and output the water demand sequence after regulation. Using the water demand sequence after regulation, a set of hydrological scenarios is generated by the tail-preserving optimal transport method, and the set of hydrological scenarios is simplified into a representative subset of scenarios. For a subset of representative scenarios, a sub-Bruker optimization model is constructed, and water resource scheduling schemes and risk assessment reports are generated through fairness constraint projection and causal path analysis. The hydrological scene set was simplified into a representative subset, including: The scenes with concentrated hydrological scenes are divided into spatiotemporal blocks according to spatial grids and time periods; A facility location optimization model is constructed. The goal of the facility location optimization model is to identify a representative subset of scenarios from the set of hydrological scenarios. The facility location optimization model also satisfies the preset tail quota constraints and spatiotemporal block coverage constraints. Among them, the tail quota constraint is used to ensure that the representative scene subset contains a predetermined proportion of tail scenes, and the spatiotemporal block coverage constraint is used to ensure that all the divided spatiotemporal blocks are covered by a predetermined number of scenes in the representative scene subset. Solve the facility location optimization model with tail quota constraints and spatiotemporal block coverage constraints to obtain a representative scenario subset; Identifying causal relationships and generating spatiotemporal causal networks includes: Historical hydrological and meteorological time series and historical water demand time series were read using a heterogeneous computing grid. An echo state network was constructed and impulse response analysis was performed to obtain response parameters characterizing the intensity, delay, and persistence of the impact. The intensity, delay, and persistence dimensions of the response parameters are constructed as a three-dimensional causal tensor; Tensor decomposition is performed on the three-dimensional causal tensor to extract the set of dominant influencing factors and the causal transmission path graph, thus constructing a spatiotemporal causal network.

2. The method according to claim 1, characterized in that, Using the regulated water demand sequence, a set of hydrological scenarios is generated using a tail-preserving optimal transport method, including: The water demand sequence after regulation is separated into the main region and the tail region; Based on the main region and the tail region, the optimal transmission model for each region is solved to generate a set of hydrological scenarios.

3. The method according to claim 2, characterized in that, The regulated water demand sequence was separated into a tail region, including: Tail samples were extracted from the water demand sequence after regulation. Perform a generalized Pareto distribution (GPD) fit on the tail samples; The tail region was defined based on the GPD fitting results.

4. The method according to claim 1, characterized in that: The output water demand sequence after regulation includes: Based on historical supply and demand data, the water supply vulnerability index is calculated using a standardized precipitation index conversion method. Based on the set of dominant influencing factors, the entropy weight method is used to determine the factor weights and calculate the demand urgency index. The water supply vulnerability index and the demand urgency index are integrated into a unified linear factor. The water demand cutoff threshold is adjusted using the unified linear factor to obtain the dynamic trigger threshold. A dynamic trigger threshold-driven tiered buffer pool mechanism is adopted to generate a water demand sequence after regulation.

5. The method according to claim 4, characterized in that, A dynamic trigger threshold-driven tiered buffer pool mechanism is adopted to generate a regulated water demand sequence, including: Apply hysteresis control to the dynamic trigger threshold to determine a stable activation level; Based on a stable activation level, the water demand reduction ratio is determined using an exponential release function. Based on the water demand reduction ratio, the original water demand is reduced to generate a post-regulation water demand sequence.

6. The method according to claim 1, characterized in that, A sub-Bruker optimization model is constructed, and water resource allocation schemes are generated through fairness constraint projection and causal path analysis, including: Based on a representative subset of scenarios, the sub-Bruker optimization model is solved to obtain a preliminary scheduling scheme. The preliminary scheduling plan is projected onto a predefined fairness polyhedron, and fairness constraint projection is performed to generate a water resource scheduling plan.

7. The method according to claim 6, characterized in that, Constructing a distributed bar optimization model includes: Calculate the first and second moments based on a representative subset of scenes; Construct a Wasserstein spherical uncertainty set using the first moment, the second moment, and the preset Wasserstein radius; The spherical bar optimization model is constructed to minimize the conditional risk value under the worst-case scenario within the Wasserstein spherical uncertainty set.

8. The method according to claim 1, characterized in that, Generating heterogeneous computational grids includes: Construct a supply and demand pressure field based on regional basic data; Extract the geometric features of the supply and demand pressure fields; Constructing a mesh density control function based on geometric features; By applying a grid density control function, the grid is refined in regions with drastic changes in supply and demand pressure fields, generating heterogeneous computational grids.

9. The method according to claim 8, characterized in that, Extract the geometric features of the supply and demand pressure fields, and construct a grid density control function based on these features, including: Calculate the Hessian matrix of the supply and demand pressure field at any spatial point; Find the eigenvalues ​​λ1 and λ2 of the Hessian matrix; Summing the absolute values ​​of eigenvalues ​​λ1 and λ2 yields |λ1| + |λ2|. |λ1|+|λ2| is used as the grid density control function to characterize the drastic changes in the supply and demand pressure field.

Images