Method and system for fixing person with water

Abstract

The invention discloses a method and system for determining people by water, and the method comprises the steps: collecting the data of a research region, and constructing a people-water coupling topological structure of the research region; constructing a human-water joint graph time sequence prediction model, and solving the model to obtain population distribution data and water system state data of a future time sequence; constructing a population-water resource-policy multi-subject coupling model based on a robust Stackelberg game, and updating a population response function and a robust uncertain set boundary of the model based on historical population distribution data and water system state data of a research area; and extracting population distribution data and water system state data of a future time sequence, inputting the data into the updated population-water resource-policy multi-subject coupling model, calculating an optimal policy regulation variable, and solving the model to obtain a population configuration scheme of the research area. According to the method, accurate matching of population distribution and water resource bearing capacity is realized, and scientificity, robustness and dynamic suitability of policy regulation and control are improved.

Description

Technical Field

[0001] This invention relates to the field of water resource management and population planning technology, and in particular to a method and system for determining population size based on water availability. Background Technology

[0002] With the acceleration of urbanization and the prominent effect of population agglomeration, the contradiction between water shortage and excessive population growth has become increasingly acute, and "determining population size based on water availability" has become a core strategy to ensure regional sustainable development.

[0003] Existing technologies related to population allocation based on water availability mostly employ a sequential modeling approach of "population prediction - water consumption quota conversion," treating population and water resources as independent research objects and establishing a correlation solely through scalar conversion. This ignores the pressure transmission of population agglomeration on the water system and the feedback effect of the water system's supply status on population movement, leading to significant deviations between predicted results and actual population-water coupling. Existing technologies often treat key parameters such as water resource volume and supply as fixed values, or use single-point prediction and simple scenario enumeration to handle uncertainty, failing to address risks from multi-source uncertainties such as precipitation fluctuations and economic growth fluctuations, resulting in insufficient robustness of policy control solutions. Most models treat population as a passively controlled object, failing to consider the proactive migration decisions of population-industry actors based on factors such as policy and water costs, thus failing to accurately capture the two-way interaction between policy and population behavior. Existing solutions are mostly static planning, lacking a closed-loop update mechanism based on measured data, making it difficult to adapt to dynamic scenarios such as changes in the spatiotemporal distribution of water resources and industrial restructuring, resulting in poor timeliness and feasibility of population allocation solutions.

[0004] Therefore, there is an urgent need to construct a dynamic water-based population determination method that can characterize the bidirectional coupling relationship between humans and water, accurately respond to multi-source uncertainties, and take into account the active response characteristics of the population, in order to overcome the shortcomings of existing technologies.

[0005] This invention proposes a "population-based" approach to address the aforementioned problems, achieving a precise match between population distribution and water resource carrying capacity, and enhancing the scientific, robust, and dynamic adaptability of policy regulation. Summary of the Invention

[0006] Objective of the invention: To provide a method for determining population size based on water availability, thereby addressing the aforementioned problems in existing technologies. Furthermore, to provide a system for determining population size based on water availability.

[0007] Technical solution: The "water-based population control" method includes the following steps:

[0008] Step S1: Collect data of the study area and construct the human-water coupling topology of the study area;

[0009] Step S2: Construct a human-water joint graph time series prediction model based on the human-water coupled topology, and solve for the future time series population distribution data and water system state data;

[0010] Step S3: Construct a robust Stackelberg game-based multi-agent coupling model of population, water resources, and policy, and update the model parameters;

[0011] Step S4: Input the future time series data into the updated model and solve for the population allocation scheme of the study area.

[0012] According to one aspect of this application, step S1 is further comprising:

[0013] Step S11: Collect historical population distribution data, water system status data, and multi-source basic data for the study area. The multi-source basic data includes: hydrological and water resources data, population and industry data, infrastructure data, spatial geographic data, and policy and economic support data.

[0014] Step S12: Based on the multi-source basic data of the study area, a graph structure modeling method is used to construct the human-water coupled topology of the study area. The human-water coupled topology includes: population nodes, water system nodes, and associated edges connecting the two types of nodes.

[0015] According to one aspect of this application, step S12 further comprises:

[0016] Step S12a: Divide the study area into nodes using streets, towns, or communities as basic units;

[0017] Step S12b: Define the surface water node, groundwater node, and water supply facility node respectively;

[0018] Step S12c: Construct internal edges of the water system based on the water pipeline network and river flow path, and connect internal nodes of the water system and water system nodes with population nodes;

[0019] Step S12d: Construct population flow-commuting connection edges based on transportation network connectivity and commuting flow data;

[0020] Step S12e: Based on the actual water supply and demand relationship and population commuting characteristics in the study area, optimize the weight assignment of the associated edges, and finally form a human-water coupled topology structure containing node set, edge set and complete attribute information.

[0021] According to one aspect of this application, step S2 further comprises:

[0022] Step S21: Construct an integrated model architecture that combines graph neural networks and time-series prediction networks to obtain a human-water joint graph time-series prediction model;

[0023] Step S22: Extract the attribute data of each node in the human-water coupling topology to form an initial node feature matrix, map the feature data to the [0,1] interval, and generate a standardized node feature matrix;

[0024] Step S23: Extract multi-source basic data and perform unified time alignment to obtain time series samples of a preset length, and generate exogenous feature time series matrix;

[0025] Step S24: Concatenate the standardized node feature matrix and the exogenous feature time series matrix along the time dimension to obtain the model input data matrix;

[0026] Step S25: Train the human-water joint graph time series prediction model based on the model input data matrix, and solve the model to obtain the future time series population distribution data and water system state data.

[0027] According to one aspect of this application, step S25 further comprises:

[0028] Step S25a: Randomly divide the input data matrix into a training set, a validation set, and a test set in a ratio of 7:2:1;

[0029] Step S25b: Set the batch size and learning rate for model training, use the adaptive moment estimation (Adam) optimizer to minimize the mean squared error loss function between the predicted and the true values, and set the maximum number of training rounds;

[0030] Step S25c: Based on the validation set data, the mean absolute error and root mean square error are used to evaluate the prediction accuracy of the model. If the prediction accuracy does not reach the preset threshold, the convolution kernel size of the graph neural network, the time step size of the temporal prediction network and the number of network layers are adjusted, and iterative training is performed until the prediction accuracy reaches the preset threshold to obtain the trained human-water joint graph temporal prediction model.

[0031] Step S25d: Load the trained human-water joint graph time series prediction model into the inference environment, input the preprocessed complete input data matrix, and solve for the future time series population distribution data and water system state data through forward propagation.

[0032] According to one aspect of this application, step S3 further comprises:

[0033] Step S31: Taking the water system, population-industry behavior, and policy-making as the main players in the game, and determining their respective decision-making authority, construct the objective functions for the water system player, the population-industry behavior player, and the policy-making player respectively;

[0034] Step S32: Set water system constraints, population-industry behavior constraints, and policy-making constraints respectively;

[0035] Step S33: Introduce uncertain variables, calculate the initial probability distribution of each uncertain variable based on historical data of the study area, and construct a robust uncertainty set with the initial distribution as the center and the Wasserstein distance as the radius;

[0036] Step S34: Update the population response function and robust uncertainty set boundary of the model based on historical population distribution data and water system status data of the study area.

[0037] According to one aspect of this application, step S34 further comprises:

[0038] Step S34a: Extract historical population distribution data and water system status data for the study area;

[0039] Step S34b: Update the weight coefficients in the population response function using the maximum likelihood estimation method to maximize the likelihood value between the measured population state and the model-predicted population state.

[0040] Step S34c: Optimize the boundary of the robust uncertainty set by adjusting the Wasserstein distance threshold so that the set contains more than 90% of the measured uncertainty scenarios;

[0041] Step S34d: Use the analytic hierarchy process (AHP) combined with expert scoring to update the weight coefficients of each objective function;

[0042] Step S34e: Substitute the updated parameters into the game model, and use measured data to verify the deviation between the policy regulation variables output by the model and the actual policy effect. If the deviation exceeds the preset threshold, repeat the parameter update process until the deviation meets the requirements, and complete the model parameter update.

[0043] According to one aspect of this application, step S4 further comprises:

[0044] Step S41: Extract future time series population distribution data and water system status data, and input them into the updated population-water resources-policy multi-agent coupling model;

[0045] Step S42: Solve for the equilibrium solution of the robust Stackelberg game model to obtain the optimal policy control variables;

[0046] Step S43: Based on the optimal policy control variables and the optimal population response vector, generate population allocation schemes by region and time series;

[0047] Step S44: Output the population allocation plan in the form of tables and charts, and generate a plan description document, including the basis for the plan, implementation steps, expected effects and risk response measures.

[0048] According to one aspect of this application, step S42 further comprises:

[0049] Step S42a: Within the robust uncertainty set, the policy-making entity determines the initial policy regulation variables with the objective of minimizing the overall system cost;

[0050] Step S42b: Based on the initial policy, population-industry actors determine the optimal population response vector with the goal of maximizing their overall utility.

[0051] Step S42c: Repeat the solution process for the policy-making subject and the population-industry behavior subject until the policy regulation variable and the population response vector no longer change, and obtain the Stackelberg equilibrium solution, that is, the optimal policy regulation variable.

[0052] According to one aspect of this application, step S43 further comprises:

[0053] Step S43a: Based on the optimal policy control variables and the optimal population response vector, generate population allocation schemes by region and time series;

[0054] Step S43b: Use the human-water joint prediction model to calculate the water system carrying capacity status after the implementation of the population allocation scheme, and verify whether the generated population allocation scheme meets the water resource carrying capacity threshold constraints of each region.

[0055] Step S43c, supplementary suggestions for safeguards to ensure the implementation of the plan, including plans for supporting water supply facilities, suggestions for optimizing transportation and commuting facilities, and detailed rules for industry support and restrictive policies.

[0056] According to another aspect of this application, a system for determining population size based on water availability is provided, comprising:

[0057] At least one processor; and

[0058] A memory communicatively connected to at least one of the processors; wherein,

[0059] The memory stores instructions executable by the processor, which are used to implement the "water-based population determination" method described in any of the above technical solutions. Beneficial effects: The "water-based population determination" method and system enable collaborative modeling of population and water systems, accurately capturing their bidirectional coupling relationship, improving prediction accuracy, ensuring that optimal policies can maintain population-water balance in most uncertain scenarios, achieving accurate bidirectional interaction modeling of policy and population behavior, adapting to dynamically changing population-water system scenarios, and improving the feasibility and timeliness of the solution. Attached Figure Description

[0060] Figure 1 This is a flowchart of the present invention.

[0061] Figure 2This is a flowchart of step S1 of the present invention.

[0062] Figure 3 This is a flowchart of step S2 of the present invention.

[0063] Figure 4 This is a flowchart of step S3 of the present invention.

[0064] Figure 5 This is a flowchart of step S4 of the present invention. Detailed Implementation

[0065] like Figure 1 As shown, the following technical solution is proposed. According to one aspect of this application, a method for "determining population size based on water availability" is provided, characterized by comprising the following steps:

[0066] Step S1: Collect data of the study area and construct the human-water coupled topology of the study area. The human-water coupled topology includes: population nodes, water system nodes and associated edges connecting the two types of nodes.

[0067] Step S2: Extract human-water coupling topological structure feature data, hydrological and water resource data, population and industry data, infrastructure data, spatial geographic data, and policy and economic auxiliary data as inputs to the pre-constructed human-water joint graph time series prediction model, and solve the model to obtain future time series population distribution data and water system state data;

[0068] Step S3: Construct a multi-agent coupling model of population, water resources, and policy based on robust Stackelberg game, and update the population response function and robust uncertainty set boundary of the model based on historical population distribution data and water system status data of the study area;

[0069] Step S4: Extract future time series population distribution data and water system status data, input them into the updated population-water resources-policy multi-agent coupling model, calculate the optimal policy regulation variables, and solve the model to obtain the population allocation scheme for the study area.

[0070] This application discloses a water-based population determination method, which can achieve precise matching between population distribution and water resource carrying capacity, and improve the scientific nature, robustness and timeliness of policy regulation. First, it constructs a foundation of human-water coupled data and spatial correlation, then obtains the future state of the human-water system through collaborative prediction, solves the optimal regulation policy and population allocation scheme with the help of multi-stakeholder game, and finally continuously optimizes and adapts to system changes through dynamic closed loop, forming a full-process technical system of "data-modeling-decision-optimization".

[0071] like Figure 2 As shown, according to one aspect of this application, step S1 further comprises:

[0072] Step S11: Collect historical population distribution data, water system status data, and multi-source basic data for the study area. The multi-source basic data includes: hydrological and water resources data, population and industry data, infrastructure data, spatial geographic data, and policy and economic support data.

[0073] Step S12: Based on the multi-source basic data of the study area, a graph structure modeling method is used to construct the human-water coupled topology of the study area. The human-water coupled topology includes: population nodes, water system nodes, and associated edges connecting the two types of nodes.

[0074] Traditional methods of determining population density based on water resources fail to construct a spatial framework for the coupling of people and water. They rely solely on administrative unit statistics for macroscopic analysis, which cannot characterize the cross-unit relationship between water transport and population commuting. This results in a lack of spatial coupling foundation for subsequent modeling. In this embodiment, multi-source data fusion is employed because a single data source cannot fully support human-water coupling modeling. For example, hydrological data alone cannot characterize the relationship between population and water systems; therefore, it is necessary to integrate hydrological, population, and infrastructure data.

[0075] The multi-dimensional data and the division of topological nodes using street / township level units are because this scale can both ensure data availability and accurately depict the spatial matching relationship between population agglomeration and water supply services.

[0076] The study area boundaries were clearly defined down to the street / township level, and administrative / functional units were divided. Multi-source basic data were collected, covering hydrological and water resources (reservoir, river, and groundwater monitoring data), population and industry (resident population, migrant population, industry type, and water consumption quotas), infrastructure (water supply network, transportation and commuting facilities), spatial geography (administrative boundaries, digital elevation model (DEM), and land use), and policy and economic data (settlement policies, water resource control indicators, and GDP growth rate). Outliers were removed, missing values ​​were filled using linear interpolation, and coordinates and time standards were standardized to complete data quality control. A human-water coupled topology was constructed, defining population nodes (street / township level, including administrative, population, and industry attributes) and water system nodes (surface water, groundwater, and water supply facilities, including location, function, and status attributes). Hydraulic connection edges (based on the direction of the pipeline / river) and commuting connection edges (retaining commuting flow ratios ≥5%) were constructed. After connectivity verification and redundant edge deletion, the data was stored in graph exchange XML format (GEXF).

[0077] In one embodiment, specifically:

[0078] The area is defined by the administrative boundaries of City A, including 3 urban districts (District A, District B, and District C) and 2 counties (County D and County E), and further subdivided into street / township level units, resulting in a total of 58 administrative / functional units.

[0079] Hydrological monitoring data for City A over the past 15 years were collected, including monthly time-series data on reservoir capacity, water level, and inflow / outflow of Reservoir B (source: Water Resources Bulletin of City A Water Resources Department); time-series data on flow and water level at three sections of River C (source: Hydrological Monitoring Station); quarterly data on the distribution, single-well output, and groundwater depth of 52 groundwater well groups (source: Groundwater Dynamic Monitoring Bulletin of City A); monthly data on the location, treatment capacity, and water supply of two reclaimed water plants and eight main booster pump stations (source: Operation Ledger of City A Water Supply Group); and monthly monitoring data on chemical oxygen demand (COD) and ammonia nitrogen concentrations of surface water, groundwater, and reclaimed water (source: Environmental Monitoring Station of City A).

[0080] Collect population census data, statistical yearbooks and public security household registration data of City A over the past 15 years, and obtain time series data on the number of permanent residents, the number and proportion of floating population, and the structure of the workforce in industries (subdivided into primary, secondary and tertiary industries) for 58 street / township units; collect data on industry type, industry output value and industry water quota for each unit (source: City A Statistical Yearbook, Industrial Park Operation Report).

[0081] The city A's water supply group obtained vector topology data of the water transmission network (including pipe node coordinates, pipe segment connection relationships, pipe diameter, and length) through its Geographic Information System (GIS). The city A's planning bureau obtained data on the specific geographical locations and service areas of two reclaimed water plants and eight booster pump stations. The city A's transportation bureau obtained data on the routes, station / exit locations, and coverage areas of one rail transit line, two expressways, and three national highways through its GIS database.

[0082] Obtain administrative boundary vector data (street / township level), 30m resolution digital elevation model (DEM) data, and land use type data (arable land, construction land, green space) from the Geographic Information Public Service Platform of the Natural Resources Bureau of City A.

[0083] Collected records of settlement policy adjustments in City A over the past 15 years, total water resource control indicators (source: City A government portal website), regional GDP growth rate, and annual data on per capita disposable income of residents (source: City A Statistical Yearbook).

[0084] Outliers in the hydrological data were removed (such as abnormally high water level data for a certain month in 2018 for Reservoir B). The missing data on groundwater depth were supplemented using linear interpolation (two months of data for a well group in Ding County in 2021 were missing). All spatial data were uniformly converted to the National Geodetic Coordinate System 2000, and the time scale was uniformly set to monthly.

[0085] According to one aspect of this application, step S12 further comprises:

[0086] Step S12a: Divide the study area into nodes based on streets, townships or communities as basic units. The attributes of each node include: unique node identifier, administrative boundary vector coordinates, number of permanent residents, proportion of floating population, industrial employment population structure, main industry types and output value;

[0087] Step S12b: Define surface water nodes, groundwater nodes, and water supply facility nodes respectively. The attributes of each node include: unique node identifier, geographical coordinates, core functional parameters, and status parameters.

[0088] Step S12c: Construct internal edges of the water system based on the water pipeline network and river flow path, and connect internal nodes of the water system and water system nodes with population nodes;

[0089] Step S12d: Construct population flow-commuting connection edges based on transportation network connectivity and commuting flow data;

[0090] Step S12e: Based on the actual water supply and demand relationship and population commuting characteristics in the study area, optimize the weight assignment of the associated edges, and finally form a human-water coupled topology structure containing node set, edge set and complete attribute information.

[0091] In this embodiment, based on multi-source data, a graph structure modeling method is used to construct the human-water coupled topology of the study area. The specific process is as follows:

[0092] Population-related nodes: Nodes are divided based on streets, townships or communities within the study area. The attributes of each node include: unique node identifier, administrative boundary vector coordinates, number of permanent residents, proportion of floating population, industrial employment population structure, main industry types and output value;

[0093] Water system nodes: Surface water nodes (reservoirs, main river sections), groundwater nodes (groundwater well clusters), and water supply facility nodes (reclaimed water plants, main booster pump stations) are defined respectively. The attributes of each node include: unique node identifier, geographical coordinates, core functional parameters (reservoir capacity / pump station water supply capacity / well cluster output), and status parameters (water level / water supply / water quality index time series data).

[0094] Define associated edges and calculate and assign weights:

[0095] Hydraulic connection edges: These connect nodes within the water system and nodes within the water system to population-type nodes. Edges within the water system are constructed based on the direction of the water supply network and the river flow path. Edge attributes include: unique edge identifier, connection node pair, pipe diameter / river cross-sectional area, flow direction, and water supply capacity. Edges between water system and population-type nodes are constructed based on the water supply area affiliation. Edge attributes include: unique edge identifier, connection node pair, water supply flow percentage, and time-series data.

[0096] Population Flow-Commuting Connection Edge: Only connects population-type nodes, constructed based on transportation network connectivity and commuting flow data. Edge attributes include: unique edge identifier, connected node pair, commuting flow time series data, transportation mode type (rail transit / highway), and travel time.

[0097] A graph connectivity analysis algorithm is used to verify the integrity of the topology and ensure the connectivity between water system nodes and population nodes. Based on the actual water supply and demand relationship and population commuting characteristics in the study area, the weight assignment of associated edges is optimized, and finally a human-water coupled topology containing node sets, edge sets and complete attribute information is formed.

[0098] In one embodiment, specifically:

[0099] Population-related nodes: 58 streets / townships are used as nodes. Node attributes include: unique identifier (e.g., District A-1 Street), administrative boundary coordinates, number of permanent residents in 2023 (85,000 people in District A-1 Street), proportion of floating population (18%), industrial employment population structure (42% in secondary industry and 55% in tertiary industry), and main industry types (equipment manufacturing and service industry).

[0100] Water system nodes: Define surface water nodes (Reservoir B, three sections of River C), groundwater nodes (five well clusters), and water supply facility nodes (two reclaimed water plants and eight booster pump stations). Node attributes include: unique identifier, geographic location coordinates, core functional parameters, and monthly average water supply time series data.

[0101] Based on the vector data of the water transmission network, construct hydraulic connection edges between nodes within the water system and nodes related to the water system and the population, such as the edge between Reservoir B and the booster pump station in Area A (pipe diameter 1200mm, length 15km, water transmission capacity 100,000 m³). 3 / d), A-zone booster pump station-A-zone-1 street side (water supply flow accounts for 25%); based on the C-channel water flow path, construct the hydraulic connection edge of the upstream section-midstream section-downstream section of the C-channel (flow transmission coefficient 0.85).

[0102] Based on commuting survey data of City A, commuting connection edges are constructed between population nodes. Only connections with commuting traffic accounting for ≥5% are retained, such as the edge of Street 1 in District A to Street 3 in District B (commuting traffic of 2000 people / day, accounting for 6% of the floating population of Street 1 in District A).

[0103] A graph connectivity analysis algorithm was used to verify that all population nodes are connected to at least one water system node. Redundant edges (such as the edge from District A to 2 streets to County D to 1 township, which accounts for 3% of the commuter traffic) were removed, resulting in a human-water coupled topology with 65 nodes (58 population nodes and 7 water system nodes) and 128 associated edges, which was stored in GEXF format.

[0104] Multi-source data fusion ensures the integrity of information in all dimensions of the human-water system, clarifies the spatial coupling relationship between the population and the water system through topological structure, and solves the problem of fuzzy spatial correlation in traditional methods.

[0105] like Figure 3 As shown, according to one aspect of this application, step S2 further comprises:

[0106] Step S21: Construct an integrated model architecture that combines graph neural networks and temporal prediction networks to obtain a human-water joint graph temporal prediction model. The architecture consists of a graph structure encoding module, a temporal feature extraction module, and a joint output module. The graph structure encoding module uses a multi-layer graph attention network (GAT) and configures a multi-head attention mechanism to capture the spatial dependencies between nodes. The temporal feature extraction module uses a multi-layer long short-term memory network (LSTM) to extract temporal evolution features. The joint output module uses a cross-attention mechanism to fuse spatial features and temporal features, and outputs the population distribution prediction results and water system state prediction results respectively through a two-branch fully connected layer.

[0107] The graph structure encoding module uses a 2-layer GAT to capture the spatial dependencies between nodes in the human-water coupled topology.

[0108] The input to the first layer of GAT is a standardized node feature matrix X∈R^(N×F_in), where N is the number of nodes and F_in is the input feature dimension. The first layer of GAT uses K_1=8 attention heads, each with an output dimension of d_1=16. After concatenation, the output dimension is K_1×d_1=128. The first layer of GAT is followed by an ELU activation function and a Dropout layer (with a Dropout rate of 0.3).

[0109] The second-layer GAT uses K_2=1 attention heads with an output dimension of d_2=64, which is used to aggregate multi-head attention features. The second-layer GAT is followed by the ELU activation function, and the output spatial embedding matrix H_spatial∈R^(N×64).

[0110] The calculation process for a single attention head is as follows: For node i and its neighboring nodes j∈N(i), first calculate the attention coefficients:

[0111] e_ij=LeakyReLU(a^T[Wh_i||Wh_j]);

[0112] Where W∈R^(d×F_in) is the learnable linear transformation matrix, a∈R^(2d) is the attention vector, || denotes vector concatenation, and the negative slope of LeakyReLU is set to 0.2;

[0113] Then, the attention coefficients are normalized using the softmax function, resulting in:

[0114] α_ij=softmax_j(e_ij)=exp(e_ij) / Σ_{k∈N(i)}exp(e_ik);

[0115] Finally, calculate the output features of node i:

[0116] h'i=σ(Σ{j∈N(i)}α_ij·Wh_j);

[0117] Where σ is the activation function.

[0118] The temporal feature extraction module uses a 3-layer LSTM to capture the temporal evolution of exogenous features.

[0119] The input to the first LSTM layer is the exogenous feature time series matrix:

[0120] X_exo∈R^(N×T×F_exo);

[0121] Where T is the time step and F_exo is the exogenous feature dimension;

[0122] The first LSTM layer has 128 hidden units (h_1=128), and the output sequence length is the same as the input.

[0123] The number of hidden units in the second LSTM layer is h_2=64;

[0124] The third LSTM layer has 64 hidden units (h_3=64) and outputs only the hidden state of the last time step as the temporal feature vector H_temporal∈R^(N×64). Dropout layers are set between each LSTM layer with a dropout rate of 0.2.

[0125] The calculation process of LSTM cells follows a standard formula:

[0126] Forget gate: f_t = σ(W_f·[h_{t-1},x_t] + b_f);

[0127] Input gate: i_t = σ(W_i·[h_{t-1},x_t] + b_i);

[0128] Candidate memory: C * _t=tanh(W_C·[h_{t-1},x_t]+b_C);

[0129] Memory update: C_t = f_t⊙C_{t-1} + i_t⊙C * _t;

[0130] Output gate: o_t = σ(W_o·[h_{t-1},x_t] + b_o);

[0131] Hidden state: h_t = o_t ⊙ tanh(C_t);

[0132] Where σ is the sigmoid function and ⊙ represents element-wise multiplication.

[0133] The joint output module employs a cross-attention mechanism to achieve deep fusion of spatial and temporal features, and specifically includes the following sub-modules:

[0134] Cross-attention fusion layer: Using the spatial embedding matrix H_spatial as the query and the temporal feature vector H_temporal as the key and value, the cross-attention output is calculated, specifically:

[0135] Q=H_spatial·W_Q, K=H_temporal·W_K, V=H_temporal·W_V;

[0136] Where W_Q, W_K, W_V∈R^(64×64) are learnable parameter matrices; attention weights A=softmax(QK^T / √d_k), where d_k=64 is a scaling factor; fused features H_fused=A·V+H_spatial, where addition implements residual connections.

[0137] Two-branch fully connected layers: The fused feature H_fused, after LayerNormalization, is input into the population prediction branch and the water system prediction branch respectively. The population prediction branch contains two fully connected layers. The first layer has an input dimension of 64 and an output dimension of 32, with the activation function being ReLU. The second layer has an input dimension of 32 and an output dimension of the predicted target number n_pop of population-type nodes (including three indicators: permanent resident population, proportion of floating population, and industrial employment population structure), with no activation function. The water system prediction branch has the same structure, with the output dimension being the predicted target number n_water of water system nodes (including three indicators: available water volume, facility operation load rate, and water quality indicators).

[0138] The dimensions of the model input data matrix are:

[0139] [N×T×(F_in+F_exo)];

[0140] Where N=65 (58 population nodes + 7 water system nodes), T=12 (time step of 12 months), F_in=8 (node ​​attribute feature dimension), F_exo=6 (exogenous feature dimension), and the total feature dimension is 14.

[0141] The model output consists of two parts: a population distribution prediction matrix.

[0142] Y_pop∈R^(N_pop×T_pred×3);

[0143] Where N_pop=58 is the number of population nodes, T_pred is the prediction time length (short-term prediction T_pred=60 months, medium-term prediction T_pred=120 months), and 3 is the number of prediction indicators; the water system state prediction matrix Y_water∈R^(N_water×T_pred×3), where N_water=7 is the number of water system nodes.

[0144] A multi-task learning loss function is adopted to simultaneously optimize the population prediction task and the water system state prediction task. The loss function is in the form of:

[0145] L_total=λ_pop·L_pop+λ_water·L_water+λ_reg·L_reg;

[0146] Where L_pop is the population prediction loss, expressed in MSE form:

[0147] L_pop=(1 / N_pop)·Σ_{i=1}^{N_pop}(y * _pop^i-y_pop^i)^2;

[0148] L_water represents the predicted loss of the water system state;

[0149] Also using the MSE format:

[0150] L_water=(1 / N_water)·Σ_{j=1}^{N_water}(y * _water^j-y_water^j)^2;

[0151] L_reg is the L2 regularization term, L_reg=Σ||W||_2^2, used to prevent the model from overfitting.

[0152] The weighting coefficients λ_pop, λ_water, and λ_reg are set using an uncertain weighting strategy:

[0153] λ_pop=1 / (2σ_pop^2), λ_water=1 / (2σ_water^2);

[0154] σ_pop and σ_water are learnable task uncertainty parameters, both initially set to 1.0; λ_reg is fixed at 0.0001. This strategy enables the model to automatically balance the weights of the two prediction tasks during training, with lower weights for tasks with greater uncertainty, thereby improving overall prediction performance.

[0155] In this embodiment, the model automatically balances weights using a homoscedastic uncertainty modeling method, transforming task uncertainty into learnable parameters to achieve automatic weight balancing between population prediction and water system state prediction tasks. Specifically:

[0156] Define two scalar learnable parameters, log_sigma_pop and log_sigma_water (both initialized to 0.0), and obtain the task uncertainty parameters sigma_pop and sigma_water through the exponential operation exp(log_sigma);

[0157] According to the formula λ_pop=1 / (2σ_pop) 2 ), λ_water=1 / (2σ_water) 2 Calculate the loss weights for the two tasks, with the L2 regularization term L_reg = Σ||W||_2 2 The weight coefficients λ_reg (where W represents all learnable parameters of the model) are fixed at 0.0001, and the total loss function is modified as follows:

[0158] L_total=λ_pop·L_pop+λ_water·L_water+λ_reg·L_reg+log(σ_pop)+log(σ_water);

[0159] The newly added logarithmic term can improve the numerical stability of weight adjustment during training and prevent the model from ignoring any task;

[0160] During model training, backpropagation is used to simultaneously optimize the core parameters of the model and log_sigma_pop and log_sigma_water. If the prediction error of a certain task is large and the uncertainty is high, the corresponding sigma will automatically increase, thereby reducing the loss weight λ of that task. The model prioritizes optimizing tasks with low uncertainty and good prediction performance. In the later stage of training, sigma tends to stabilize and λ also converges, achieving an adaptive balance of the weights of the two tasks.

[0161] Step S22: Extract the attribute data of each node in the human-water coupling topology to form an initial node feature matrix. The population node features include the number of permanent residents, the proportion of floating population, and the structure of the industrial workforce. The water system node features include the available water volume, facility operation load rate, and water quality indicators. Map the feature data to the [0,1] interval to generate a standardized node feature matrix.

[0162] Step S23: Extract hydrological and water resources data, population and industry data, infrastructure data, spatial geographic data, and policy and economic support data, and perform unified time caliber alignment to obtain time series samples of a preset length, and generate exogenous feature time series matrix.

[0163] Step S24: Concatenate the standardized node feature matrix and the exogenous feature time series matrix along the time dimension to form a model input data matrix with the dimension of [number of nodes × time step × feature dimension].

[0164] Step S25: Divide the input data matrix into a training set, a validation set, and a test set. Train the human-water joint graph time series prediction model using a multi-task learning loss function. The multi-task learning loss function includes a population prediction loss term, a water system state prediction loss term, and a regularization term. The weights of each loss term are adaptively determined using an uncertainty weighting strategy. Adjust the model parameters and hyperparameters, solve the trained human-water joint graph time series prediction model, and obtain the future time series population distribution data and water system state data.

[0165] Traditional methods typically estimate a region's population carrying capacity limit using a linear formula based on multi-year average total water resources and per capita water consumption quotas, for example, a formula similar to:

[0166] Upper limit of carrying capacity = Total water resources ÷ Per capita water consumption quota;

[0167] Alternatively, a linear weighted approach by industry can be used to obtain a static upper bound on population size. However, this type of method has obvious drawbacks: it ignores the randomness and extreme nature of precipitation and water inflow, and cannot reflect the uncertain changes in hydrological conditions over a future period; at the same time, it treats population as a fixed input or a simple prediction result, without considering the behavioral feedback of population to changes in policies, water prices, and infrastructure.

[0168] There is a lack of unified dynamic modeling between population and water systems. Currently, most models adopt a sequential approach of "first predicting population, then converting water consumption to per capita quotas," that is:

[0169] Use statistical or machine learning models to predict the population size of each region;

[0170] Water demand is obtained by multiplying the projected population by a fixed or stratified quota;

[0171] Compare with available water volume to assess carrying capacity and risk.

[0172] This "separate forecasting" ignores the two-way interaction between the human and water systems: on the one hand, population concentration increases water demand, leading to increased local water resource pressure; on the other hand, long-term water shortages can, in turn, inhibit population inflow or promote population outflow through factors such as quality of life, cost, and industrial layout. Because this two-way feedback is not explicitly modeled, the model has poor prediction accuracy and adaptability when facing scenarios such as policy changes, the implementation of new water source projects, and extreme droughts.

[0173] This application constructs a unified graph structure that simultaneously includes population nodes and water source / supply nodes, and uses graph neural networks combined with time series models to jointly predict future population distribution and water system status, in order to replace the traditional sequential mode of "population prediction + quota conversion water use".

[0174] Specifically, the region is divided into several residential area nodes and water system nodes. Residential area nodes represent streets, towns, or communities; water system nodes include surface water reservoirs, groundwater well groups, reclaimed water plants, and major booster pump stations. Various connections between nodes are represented by edges: firstly, hydraulic connections, such as water transmission networks and river connections, reflecting water transfer and supply capacity; secondly, population flow and commuting connections, such as rail transit lines, highway networks, and commuting directions, reflecting population migration and daily pedestrian flow.

[0175] On such a "human-water joint graph", a state vector is constructed for each node, including population-related states (number of permanent residents, proportion of floating population, number of employees in industries, etc.) and water system states (available water, operating load, water quality indicators, etc.). At the same time, exogenous features such as meteorological indicators, economic development level, and policy changes are considered. In the time dimension, the state sequence of several consecutive moments is input into a time series prediction model with graph structure perception capabilities. For example, a graph neural network is used to capture spatial dependencies, and a recurrent neural network or other time series network is used to model the temporal evolution relationship.

[0176] Through this joint graph time-series model, the system can simultaneously output the population status of each residential area node and the water resource status of each water node at several future points in each prediction period, instead of first predicting the population separately and then converting the water demand into a fixed quota. This joint prediction has several key advantages:

[0177] The model explicitly depicts the two-way coupling of the human-water system. By propagating information in a unified graph structure, it automatically learns the effect of population agglomeration on water system pressure and the feedback law of long-term water system stress on population inflow inhibition. Thus, it reflects the long-term interaction effect more realistically than a simple series system.

[0178] By utilizing spatial topology information to improve prediction accuracy and generalization ability, the topological relationships of pipeline structure, water supply zone boundaries, and transportation network are all encoded into the model through graph structure, which helps to capture the mutual influence between neighboring areas and the cross-regional allocation effect. Compared with independent modeling of each region, graph model can more effectively utilize cross-regional similarity and correlation.

[0179] It provides a unified dynamic prediction module, which adopts a human-water joint graph time series model. It can perform joint prediction of population status and water resource status of each water node within the same model. As the core simulator for calculating cost function and verifying carrying capacity robustness, it significantly reduces the difficulty of building complex analytical models and improves the fitting accuracy of real systems.

[0180] In this embodiment, an integrated architecture is constructed that combines a graph attention neural network and a temporal feature extraction module, including a graph structure encoding module (2-layer GAT, which takes topological structure and node features as input and outputs spatial embedding vectors), a temporal feature extraction module (3-layer LSTM, which takes exogenous time-series data such as meteorological and economic data as input and outputs temporal feature vectors), and a joint output module (attention fusion + double-branch fully connected layer, which simultaneously outputs short-term (1-5 years) and medium-term (5-10 years) population distribution and water system status data).

[0181] The node feature matrix is ​​min-max normalized, the time caliber of exogenous time series data is unified, outliers are filled in and converted into time series samples by sliding window method, and finally fused to form the model input matrix;

[0182] The training / validation / test sets were divided into a 7:2:1 ratio. The batch size was set to 32, the learning rate to 0.005, and the Adam optimizer and early stopping strategy were used. The prediction accuracy (MAE ≤ 5%) was ensured by iteratively adjusting the hyperparameters (GAT kernel size, LSTM time step).

[0183] Load the trained model and input preprocessed data to obtain the future population and water system node status data for each region.

[0184] In one optional implementation, the specific configuration of the Graph Attention Neural Network (GAT) includes: two GAT layers with eight attention heads per layer, a hidden layer dimension of 64, and the activation function being LeakyReLU with a negative slope of 0.2; the specific configuration of the LSTM includes: three LSTM layers with a hidden state dimension of 128 and a dropout ratio of 0.2; the joint output module includes two fully connected layers with dimensions of 256 and an output dimension, respectively. These parameters can be adjusted according to the size of the research area and the characteristics of the data. For example, for large-scale research areas with more than 100 nodes, the hidden layer dimension of the GAT can be increased to 128, and the hidden state dimension of the LSTM can be increased to 256.

[0185] According to one aspect of this application, step S25 further comprises:

[0186] Step S25a: Randomly divide the input data matrix into a training set, a validation set, and a test set in a ratio of 7:2:1;

[0187] Step S25b: Set the batch size and learning rate for model training, use the Adam optimizer to minimize the mean squared error loss function between the predicted and true values, and set the maximum number of training rounds;

[0188] Step S25c: Based on the validation set data, the mean absolute error and root mean square error are used to evaluate the prediction accuracy of the model. If the prediction accuracy does not reach the preset threshold, the convolution kernel size of the graph neural network, the time step size of the temporal prediction network and the number of network layers are adjusted, and iterative training is performed until the prediction accuracy reaches the preset threshold to obtain the trained human-water joint graph temporal prediction model.

[0189] The key parameters affect model performance as follows: Increasing the learning rate will accelerate convergence but may cause oscillations; decreasing the learning rate will improve stability but may lead to local optima. A value range of 0.001-0.01 is recommended. Increasing the batch size will improve computational efficiency but may reduce generalization ability; decreasing the batch size will enhance randomness but prolong training time. A value range of 16-64 is recommended. Increasing the number of GAT layers will improve feature extraction ability but may lead to oversmoothing. A value range of 2-4 layers is recommended. Increasing the LSTM time step will enhance long-term dependency capture ability but increase computational cost. It is recommended to determine the time step based on the periodicity of the data, with a typical value of 12-24 months.

[0190] Step S25d: Load the trained human-water joint graph time series prediction model into the inference environment, input the preprocessed complete input data matrix, and solve for the future time series population distribution data and water system state data through forward propagation.

[0191] In one embodiment, specifically:

[0192] Attribute data of 65 nodes in the human-water coupled topology were extracted to form an initial node feature matrix (65×8 dimensions). Min-max normalization was used to map the data to the [0,1] interval to generate a standardized node feature matrix.

[0193] Meteorological time-series data (precipitation, temperature, evaporation, standardized precipitation evapotranspiration index), economic time-series data (GDP growth rate, industrial investment scale), policy time-series data (adjustment of settlement policies, total water resource control indicators), and infrastructure time-series data (commissioning time of newly added water supply facilities) for City A from 2009 to 2023 were collected. The time scale was unified to monthly. Linear interpolation was used to complete the missing values ​​of precipitation data. The time-series data were converted into time-series samples by the time sliding window method (window size 12) to generate an exogenous feature time-series matrix (dimensions 65×12×6).

[0194] The standardized node feature matrix and the exogenous feature time series matrix are concatenated along the time dimension to form the model input data matrix (dimension 65×12×14).

[0195] The input data matrix was divided into a training set (2009-2019), a validation set (2020-2021), and a test set (2022-2023) in a ratio of 7:2:1.

[0196] The batch size is set to 32, the learning rate is 0.005, the Adam optimizer is used to minimize the mean squared error loss function, the maximum number of training rounds is 100, and the early stopping strategy is set to stop training if the validation set loss does not decrease for 10 consecutive rounds.

[0197] The initial validation set MAE was 8.2%, which did not meet the preset threshold (MAE≤5%). After adjusting the GAT convolution kernel size to 3, the LSTM time step to 18 months, and the number of network layers to 4, the validation set MAE dropped to 4.3% after iterative training, meeting the preset requirements.

[0198] During the training process, the hydraulic-population commuting correlation between Reservoir B, the booster pump station in District A, District A, and Street 1 was explored using GAT to strengthen the learning of the two-way coupling relationship between them.

[0199] Load the trained model into the inference environment, input the preprocessed complete input data matrix, and solve it through forward propagation. Output the permanent resident population and the proportion of migrant population for each street / township in City A for 2024-2028 (short-term) and 2029-2033 (medium-term), as well as the available water volume and operating load rate data for each water system node. For example, output the predicted permanent resident population of Street A-1 in 2028 as 92,000 (migrant population accounting for 19%), and the predicted average monthly available water volume of Reservoir B in 2028 as 85,000 m³. 3 / d.

[0200] like Figure 4 As shown, according to one aspect of this application, step S3 further comprises:

[0201] Step S31: Taking the water system, population-industry behavior, and policy-making as the main players in the game, and determining their respective decision-making authority, construct the objective functions for the water system player, the population-industry behavior player, and the policy-making player respectively;

[0202] Step S32: Set water system constraints, population-industry behavior constraints, and policy-making constraints respectively, with water resource carrying capacity as the core premise;

[0203] Step S33: Introduce uncertain variables, calculate the initial probability distribution of each uncertain variable based on historical data of the study area, and construct a robust uncertainty set with the initial distribution as the center and the Wasserstein distance as the radius. The Wasserstein distance radius is determined based on the degree of uncertainty.

[0204] Step S34: Update the population response function and robust uncertainty set boundary of the model based on historical population distribution data and water system status data of the study area.

[0205] Hydrological and population uncertainties are only treated as single-point predictions, which lacks robustness. In most existing studies, precipitation, water inflow, and population predictions are converged to a single prediction value. Optimization and regulation schemes are designed based on this "most likely scenario" and lack defense against extreme scenarios or distribution shifts. When encountering rare droughts or sudden large-scale population movements (such as major events or sudden migrations), the schemes are prone to failure, leading to the penetration of carrying capacity.

[0206] Meanwhile, closed-loop control lacks a theoretical iterative mechanism. Current "feedback adjustment" relies heavily on human experience: after the implementation of the plan, actual population and water consumption data are observed, and parameters are adjusted manually after deviations are found (such as increasing the safety factor and reducing water consumption quotas). This experience-based closed loop lacks a unified objective function constraint and convergence guarantee, making it difficult to quantify and reuse the model iteration and parameter update process, and difficult to promote as an algorithmic framework for sustainable optimization.

[0207] Therefore, in this embodiment, the problem of "determining population based on water availability" is formalized into a multi-stakeholder game problem involving population, water resources, and policy under hydrological and economic uncertainties, and solved using the sub-Bruker optimization method, thereby extending the carrying capacity from a static point value to an upper bound of population that can maintain system stability under uncertain scenarios.

[0208] This framework introduces three core components:

[0209] The first subject is the water system, whose state is composed of multiple water sources, such as surface water, groundwater, and reclaimed water, represented by a time series vector W(t). The state of the water system is affected by exogenous factors such as precipitation, upstream water inflow, and engineering scheduling. These factors are described by uncertain variables ξ. The carrying capacity of the water system is no longer simply regarded as a fixed total amount, but as the combination of available water W(ξ) under different ξ scenarios.

[0210] The second subject is the population and industry actors. Let the regional population vector P(t) represent the status of the resident population, migrant population, and industrial workforce in each sub-region (e.g., street, town). Based on factors such as policy, water costs, and local quality of life, population and industry choose to migrate in, migrate out, or remain unchanged. Their behavior can be abstracted as the optimal response that maximizes a certain utility function under given policy π and an uncertain environment ξ, i.e.:

[0211] P * (π, ξ)=argmax_{P∈allowed set}U_population(P, π, ξ);

[0212] U_population represents the combined results of population and industry in different regional combinations, such as employment opportunities, cost of living, water costs and migration costs. U_population is a strictly concave function of P and the constraint set is a convex set.

[0213] The third entity is the policymaker, who selects a set of regulatory strategy variables π, including settlement quotas for different regions, water price structures, and the pace of industrial project approvals and relocations. The goal is no longer to optimize for a specific predictive scenario, but rather, under hydrological and economic uncertainties, considering the optimal population response P. * Given (π, ξ), minimize the overall cost in the worst-case scenario, for example:

[0214] min_πsup_{Q∈U(H * P)}E_{ξ∼Q}[C(W(ξ),P * (π,ξ),π)];

[0215] Where H * P represents the distribution of uncertainties estimated based on historical data, and U(H) * P) is around H * P is a set of uncertain distributions (e.g., a sphere with a radius of Wasserstein distance), and C is a comprehensive cost function that includes water shortage penalties, overcapacity penalties, migration and industrial adjustment costs, and policy implementation costs.

[0216] Under this framework, the "carrying capacity threshold" of "determining population based on water availability" is no longer simplified to:

[0217] Carrying capacity threshold = Σ (water resource stock × coefficient) ÷ quota;

[0218] This single-point estimation is instead redefined as: in policy π * Below, for the vast majority (e.g., 95% confidence level) of possible hydrological and economic scenarios ξ, and the optimal population response to policy P... * (π * Under the condition ξ), the maximum sustainable population size range in which each sub-region will not experience unacceptable supply and demand imbalances in the long term.

[0219] By using Stackelberg game modeling, this invention achieves the following essential improvements:

[0220] First, population and industry are elevated from "passive variables" to active decision-makers in the game, so that the two-way interaction between policy adjustment and population behavior response can be clearly modeled.

[0221] Secondly, the uncertainty of hydrology and economy is expanded from single-point prediction to a set of probability distributions. The principle of worst-case reasonable scenario optimization is adopted to make the carrying capacity and regulation strategy robust and no longer highly dependent on a single accurate prediction.

[0222] Third, the concept of bearing capacity has changed from a static value to a robust upper bound within a certain range of uncertainty, which is more in line with the planning and management needs in the context of frequent extreme events.

[0223] Within this unified framework, the original closed-loop feedback steps have also been systematically transformed. By continuously observing actual population distribution, actual water consumption, and hydrological data, the population's response function to policy is periodically estimated, and the uncertainty set U(H) is updated. * P), and after each round of updates, the Stackelberg problem is solved again, thus forming a theoretical closed-loop process based on game parameter learning and uncertainty set updates, rather than simple empirical parameter tuning. This closed loop has a clear objective function and optimization criteria, which helps to continuously improve matching accuracy and system stability in long-term operation.

[0224] In this embodiment, the degenerate Stackelberg game modeling and the joint prediction based on the human-water joint graph time series model are not simply spliced ​​together, but have a significant coupling effect in terms of technical conception. At the upper-level decision-making level, the degenerate Stackelberg framework provides a theoretical framework for "determining population size based on water availability" that comprehensively describes the interaction between policymakers, population actors, and hydrological uncertainties, and can define "robust carrying capacity" and "sustainable population size under the worst-case scenario" in principle. At the lower-level prediction and simulation level, the human-water joint graph time series model provides the game framework with high-fidelity system dynamics prediction capabilities, enabling effective solutions to upper-level optimization problems even in complex real-world systems.

[0225] In one embodiment, specifically:

[0226] Water system main body (Water Resources Bureau of City A): The decision variables are the water supply allocation scheme of Reservoir B and each booster pump station;

[0227] Population-industry actors (residents and enterprises in City A): Decision variables are population migration direction (e.g., migration from County D to District A) and industrial investment site selection;

[0228] Policy-making body (Municipal A): Decision variables are the regional settlement quota and the tiered water price adjustment coefficient.

[0229] The main objective function of the water system is: minF_w = 0.5 × (water supply gap) + 0.3 × (water supply cost) + 0.2 × (penalty for exceeding water quality standards), and the weight coefficients are determined by the analytic hierarchy process.

[0230] Population-industry behavior agent comprehensive utility function:

[0231] U(P, π, ξ)=0.3×E(P, π)-0.2×C_life(P)-0.3×C_water(P, π)-0.2×C_migrate(P);

[0232] Where E(P, π) is the employment utility (positively correlated with the number of jobs in the equipment manufacturing industry in District A), C_life(P) is the cost of living (positively correlated with housing prices in District A), C_water(P, π) is the water cost (positively correlated with tiered water pricing), and C_migrate(P) is the migration cost (positively correlated with commuting costs from County D to District A). The weighting coefficients are α=0.3, β=0.2, γ=0.3, and δ=0.2 (initialized through maximum likelihood estimation).

[0233] The overall cost function of the policy-making entity system:

[0234] C(W(ξ), P) *(π,ξ),π)=0.4×C_water_short(W(ξ),P * )+0.2×C_pop_over(P * )+0.2×C_adjust(π)+0.2×C_execute(π);

[0235] The weighting coefficients are λ=0.4, μ=0.2, ν=0.2, and ζ=0.2 (determined using the analytic hierarchy process).

[0236] Water system constraints:

[0237] Water resource total constraint: The annual usable water resources of City A are ≤800 million m³ 3 ;

[0238] Water supply capacity constraint: Maximum water supply of Reservoir B ≤ 150,000 m³ 3 / d;

[0239] Water quality constraints: COD concentration ≤ 20 mg / L, ammonia nitrogen concentration ≤ 1.0 mg / L;

[0240] Hydraulic balance constraint: Inflow into water transmission network = Outflow + Loss (loss rate ≤ 5%).

[0241] Population-industry behavior constraints:

[0242] Population size constraint: The permanent resident population of Zone A is ≤350,000;

[0243] Migration control restrictions: The annual migration volume from County Ding to District Jia shall be ≤5000 people;

[0244] Industrial planning constraints: The number of water-intensive industrial projects should be ≤3 per year;

[0245] Employment matching constraint: The number of people employed in the industry ≤ the number of jobs available;

[0246] Policy constraints:

[0247] The adjustment range for tiered water pricing for residents will be ≤20%.

[0248] Policy implementation costs plus relocation and adjustment costs ≤ 50 million yuan / year;

[0249] Fairness constraint: The difference in settlement quotas between regions shall be ≤20%.

[0250] Uncertain variables are defined as follows: ξ1 (precipitation fluctuation), ξ2 (upstream water inflow change), and ξ3 (GDP growth rate fluctuation) are introduced as uncertain variables. ξ1 is represented by the coefficient of variation of precipitation time series data over the past 15 years (0.25), ξ2 is represented by the upstream water inflow guarantee rate of Reservoir B (P=90%), and ξ3 is represented by the coefficient of variation of GDP growth rate over the past 15 years (0.15).

[0251] Based on measured data from 2009 to 2023, the initial distribution of each uncertain variable was obtained using the Gaussian kernel density estimation method, and the kernel function bandwidth was determined to be 0.05 through cross-validation.

[0252] Construct a robust uncertainty set centered on the initial distribution with a Wasserstein distance of 0.1, thus extending single-scenario optimization to multi-uncertainty scenario optimization;

[0253] The principle of "minimizing worst-case costs" is adopted to ensure that the policy can still maintain the balance between humans and water under the worst-case scenarios of extreme drought (ξ1=0.4), reduced upstream water inflow (ξ2=P=95%), and declining GDP growth (ξ3=0.25).

[0254] The mathematical expression for the comprehensive utility function of the population-industry actors is:

[0255] U(P, π, ξ)=α×E(P, π)-β×C_life(P)-γ×C_water(P, π)-δ×C_migrate(P);

[0256] Wherein, U is the comprehensive utility value; P is the population distribution vector; π is the policy regulation variable vector; ξ is the uncertainty variable; E(P, π) is the employment utility function, which is positively correlated with the number of jobs in regional industries; C_life(P) is the cost of living function, which is positively correlated with regional housing prices and commodity prices; C_water(P, π) is the water cost function, which is positively correlated with tiered water pricing; C_migrate(P) is the migration cost function, which is positively correlated with commuting distance and time; α, β, γ, and δ are weighting coefficients, satisfying α+β+γ+δ=1. The initial values ​​can be determined based on historical data using the maximum likelihood estimation method, and the typical range of values ​​is α∈[0.25, 0.35], β∈[0.15, 0.25], γ∈[0.25, 0.35], δ∈[0.15, 0.25].

[0257] The employment utility is calculated as follows: Industrial policy intensity × Number of regional jobs ÷ Regional population.

[0258] E(P, π) = (Total number of jobs in the region × π_ind) ÷ P;

[0259] Wherein, π_ind represents the intensity of industrial policy (value 0~1, which can be quantified as: no industrial policy = 0, weak policy = 0.3~0.5, strong policy = 0.6~1.0, determined with reference to regional industrial support documents), the total number of regional jobs can be obtained from regional statistical yearbooks and publicly available data from the human resources and social security bureau (such as the number of newly added urban jobs in the year and the number of jobs in enterprises above a certain size), and P represents the regional population (corresponding to the population distribution vector P, directly using the resident population data in the statistical yearbook).

[0260] Cost of living = Base cost of living × Population agglomeration coefficient, i.e.:

[0261] C_life(P) = c0 × P^ρ;

[0262] Wherein, c0 is the benchmark cost of living (taken as the average monthly living expenditure per capita in the region, which can be obtained from statistical yearbooks and resident income and expenditure survey data; for example, if the average monthly living expenditure per capita in a certain region is 2,000 yuan, then c0 = 2,000), ρ is the agglomeration cost coefficient (taken as 0.1~0.5; the higher the population agglomeration, the closer ρ is to 0.5; the lower the agglomeration, the closer it is to 0.1; it can be calibrated by referring to similar regional studies), and P is the regional population (corresponding to the population distribution vector P, directly using the resident population data in the statistical yearbook).

[0263] Water cost = Comprehensive water price × Total water consumption in the region, i.e.:

[0264] C_water(P,π)=π_price×q0×P^ε;

[0265] Wherein, π_price is the comprehensive tiered water price (unit: yuan / ton, which can be obtained from public documents of the water resources bureau and the price bureau, such as the first tier of water price being 3 yuan / ton and the second tier being 4 yuan / ton, and the comprehensive water price is obtained by weighted average according to the regional water use structure), q0 is the per capita benchmark water consumption (unit: tons / person / month, which can be obtained from the statistical data of the water resources bureau, such as the per capita monthly water consumption of 10 tons in the region, then q0=10), ε is the water use scale elasticity (value is 0.8~1.0, the influence coefficient of population growth on water consumption, which is calibrated with reference to the regional water resources bulletin or similar studies), and P is the regional population (corresponding to the population distribution vector P, directly using the resident population data in the statistical yearbook).

[0266] Migration cost = cost per unit distance × average commuting distance, i.e.:

[0267] C_migrate(P) = τ × d_bar(P);

[0268] Wherein, τ is the unit distance commuting cost (unit: yuan / km, which can be calculated as: average monthly commuting expenditure per person ÷ average total monthly commuting distance. For example, if the average monthly commuting expenditure per person is 200 yuan and the average monthly commuting distance is 80 km, then τ=2.5), and d_bar(P) is the regional average commuting distance (unit: km, which can be obtained from the transportation bureau and urban commuting survey data. The higher the population concentration, the longer the average commuting distance).

[0269] The mathematical expression for the comprehensive cost function of the policy-making entity system is as follows:

[0270] C(W(ξ),P(π,ξ),π)=λ×C_water_short+μ×C_pop_over+ν×C_adjust+ζ×C_execute * ;

[0271] Wherein, C_water_short is the cost of water shortage, which is incurred when water demand exceeds available water supply; C_pop_over is the cost of population overcapacity penalty; C_adjust is the cost of industrial adjustment; C_execute is the cost of policy implementation; λ, μ, ν, and ζ are weighting coefficients, with typical values ​​of λ=0.4, μ=0.2, ν=0.2, and ζ=0.2.

[0272] According to one aspect of this application, step S34 further comprises:

[0273] Step S34a: Extract historical population distribution data and water system status data for the study area;

[0274] Step S34b: Update the weight coefficients in the population response function using the maximum likelihood estimation method to maximize the likelihood value between the measured population state and the model-predicted population state.

[0275] Step S34c: Optimize the boundary of the robust uncertainty set by adjusting the Wasserstein distance threshold so that the set contains more than 90% of the measured uncertainty scenarios;

[0276] Step S34d: Use the analytic hierarchy process (AHP) combined with expert scoring to update the weight coefficients of each objective function;

[0277] Step S34e: Substitute the updated parameters into the game model, and use measured data to verify the deviation between the policy regulation variables output by the model and the actual policy effect. If the deviation exceeds the preset threshold, repeat the parameter update process until the deviation meets the requirements, and complete the model parameter update.

[0278] In one embodiment, specifically:

[0279] Collect measured data for City A from 2021 to 2023, including statistics on the permanent and floating population of each street / township (quarterly scale, 98% completeness), and measured data on water supply and water quality at each water system node;

[0280] The maximum likelihood estimation method was used to update the weight coefficients of the population response function. After the update, α=0.32, β=0.19, γ=0.31, and δ=0.18.

[0281] Adjust the Wasserstein distance threshold to 0.12 to optimize the robust uncertainty set boundary and ensure that the set contains 92% of the measured uncertainty scenarios;

[0282] The objective function weight coefficients are updated using the analytic hierarchy process (AHP).

[0283] Substitute the updated parameters into the model and verify that the deviation between the policy regulation variables and the actual effect is 6.5% (≤8% preset threshold), thus completing the parameter update.

[0284] Existing water-based population models often treat the population as a passively controlled entity, failing to consider its active migration decisions. They frequently employ single-point forecasts or ignore uncertainties such as precipitation and the economy, resulting in insufficient policy robustness. Parameter updates are mostly empirical, lacking a theoretical mechanism, leading to poor model adaptability to real-world systems. While Stackelberg game theory models are widely used in supply chains and power systems, they are not tailored for scenarios involving multiple stakeholders—humans, water, and policy—and therefore cannot adapt to the characteristics of human-water systems.

[0285] The Stackelberg game framework was chosen because it can depict the hierarchical decision-making relationship of "policy-making subject (upper level) - population - industrial behavior subject / water system subject (lower level)," which is consistent with the actual logic of water-based population policy-making and is superior to the equal subject decision-making framework of Nash game.

[0286] Robust optimization and Wasserstein distance are introduced because Wasserstein distance can effectively measure the differences in the distribution of uncertain variables. The constructed robust uncertainty set can comprehensively cover multi-source uncertainties, which is better than the traditional scenario enumeration method (which only covers a limited number of scenarios).

[0287] The maximum likelihood estimation method is used to update parameters, which can maximize the model's predicted likelihood value based on measured data and achieve accurate parameter updates, which is superior to empirical parameter tuning.

[0288] By defining population-industry actors as active game players, and accurately capturing their migration decisions through a comprehensive utility function, a two-way interactive modeling of policy and population behavior is achieved. Robust uncertainty set and worst-case optimization criteria ensure that policies can still maintain water balance and reduce water shortage risks under extreme uncertainty scenarios. Theoretical parameter update mechanism enables accurate adaptation of the model to the actual system, improving model reliability. The Stackelberg hierarchical decision framework conforms to the actual process of government policymaking, improving policy operability.

[0289] In this embodiment, the specific form of the population response function and the parameter update method are as follows:

[0290] The population response function describes the optimal migration decisions of population-industry actors under given policy control variables π and uncertain environmental conditions ξ.

[0291] Suppose the study area has M spatial units (streets / townships), and the population state vector of the m-th unit is:

[0292] P_m=[P_m^res, P_m^flo, P_m^ind]^T;

[0293] These represent the number of permanent residents, the number of migrant workers, and the number of people employed in industries, respectively.

[0294] The population response function, in the form of a discrete choice model, describes the probability of population migrating from unit m to unit n:

[0295] Pr(m→n|π,ξ)=exp(V_n(π,ξ)) / Σ_{k=1}^{M}exp(V_k(π,ξ));

[0296] Where V_n(π, ξ) is the utility function of unit n, expressed as:

[0297] V_n(π,ξ)=α·E_n(π)-β·C_n^life-γ·C_n^water(π,ξ)-δ·C_{mn}^migrate+θ·Q_n^infra+η·S_n^water(ξ);

[0298] The meanings of each item are as follows: E_n(π) is the employment opportunity index of unit n, which is positively correlated with the quota for settling down and the amount of industrial project approval; C_n^life is the cost of living index of unit n, which is determined by housing prices, commodity prices, etc.; C_n^water(π, ξ) is the water cost of unit n, which is related to the tiered water price adjustment coefficient and the scarcity of water resources; C_{mn}^migrate is the migration cost from unit m to unit n, which is related to distance and transportation convenience; Q_n^infra is the infrastructure quality index of unit n; S_n^water(ξ) is the water resource abundance index of unit n, which is determined by the state of the water system.

[0299] The weighting coefficient vector θ=[α,β,γ,δ,θ,η]^T reflects the sensitivity of population-industry actors to various factors.

[0300] Given the utility function, the population state of cell n at time t+1 can be expressed as:

[0301] P_n^(t+1)=P_n^t+Σ_{m≠n}(P_m^t·Pr(m→n|π,ξ))-P_n^t·Σ_{k≠n}Pr(n→k|π,ξ)+ε_n;

[0302] In this context, the second term on the right represents net inflow population, and ε_n is a random disturbance term that follows a normal distribution with a mean of 0.

[0303] Suppose that the study area has historical observation data {(P_m^t, π^t, ξ^t)}, with a time span of t=1, 2, ..., T_obs. The likelihood function is constructed based on the population migration probability and population state transition.

[0304] For a single observation at time t, the likelihood of the migration flow is:

[0305] L_t(θ)=Π_{m=1}^{M}Π_{n=1}^{M}[Pr(m→n|π^t,ξ^t;θ)]^{F_{mn}^t};

[0306] Wherein, F_{mn}^t is the measured number of people who migrate from unit m to unit n at time t, which can be obtained from population census data or migrant population monitoring data.

[0307] Considering the observation error of population state transition, we introduce state likelihood:

[0308] L_state(θ)=Π_{t=1}^{T_obs}Π_{m=1}^{M}(1 / √(2πσ_m^2))·exp(-(P_m^tP * _m^t(θ))^2 / (2σ_m^2));

[0309] Among them, P * _m^t(θ) represents the population state predicted by the model under given parameters θ, and σ_m represents the standard deviation of the observation error.

[0310] The joint log-likelihood function is:

[0311] logL(θ)=Σ_{t=1}^{T_obs}Σ_{m=1}^{M}Σ_{n=1}^{M}F_{mn}^t·logPr(m→n|π^t,ξ^t;θ)+Σ_{t=1}^{T_obs}Σ_{m=1}^{M}[-(P_m^tP * _m^t(θ))^2 / (2σ_m^2)];

[0312] The maximum likelihood estimation problem is solved using the quasi-Newton method (L-BFGS-B algorithm). The specific steps are as follows:

[0313] Step a: Initialize the weight coefficient vector θ^(0). Prior empirical values ​​or uniform distribution random initialization can be used. In this embodiment, the initial values ​​are set as α^(0)=0.30, β^(0)=0.20, γ^(0)=0.30, δ^(0)=0.20, θ^(0)=0.15, and η^(0)=0.10.

[0314] Step b: Set parameter constraint boundaries. Each weight coefficient must satisfy the non-negativity constraint θ_i≥0 and the normalization constraint Σθ_i=1.

[0315] Step c: Calculate the log-likelihood value logL(θ^(k)) and the gradient vector ▽logL(θ^(k)) for the current parameters; the gradient is calculated analytically.

[0316] ΞlogL / Ξα=Σ_{t,m,n}F_{mn}^t·[E_n-Σ_kPr(m→k)·E_k]+Σ_{t,m}(P_m^tP * _m^t)·(ΞP * _m^t / Ξα) / σ_m^2;

[0317] Where Ξ represents the partial derivative;

[0318] The gradient forms for other parameters are similar;

[0319] Step d: Update parameters using the L-BFGS-B algorithm:

[0320] θ^(k+1)=θ^(k)+ρ^(k)·d^(k);

[0321] Where d^(k) is the search direction, determined by the L-BFGS approximation Hessian matrix; ρ^(k) is the step size, determined by the Armijo line search.

[0322] Step e: Determine the convergence condition: If ||θ^(k+1)-θ^(k)||_∞<ε_θ and |logL(θ^(k+1))-logL(θ^(k))|<ε_L, then stop the iteration and output the optimal parameter θ. * Otherwise, let k=k+1 and return to step c. In this embodiment, the convergence threshold is set to ε_θ=10^(-6) and ε_L=10^(-8), and the maximum number of iterations is set to 500.

[0323] Step f: Perform statistical tests on the estimation results and calculate the Fisher information matrix I(θ). * The inverse matrix of the parameter estimate is used to obtain the standard error and confidence interval of the parameter estimate, and the significance of the parameter estimate is verified.

[0324] like Figure 5 As shown, according to one aspect of this application, step S4 further comprises:

[0325] Step S41: Extract future time series population distribution data and water system status data, and input them into the updated population-water resources-policy multi-agent coupling model;

[0326] Step S42: Use the interior point method combined with particle swarm optimization algorithm to solve the equilibrium solution of the robust Stackelberg game model to obtain the optimal policy control variables. The optimal policy control variables include: regional settlement quotas, tiered water price adjustment coefficients, industrial project approval quotas, and industrial relocation timeline plans.

[0327] Step S43: Based on the optimal policy control variables and the optimal population response vector, generate population allocation schemes by region and time series;

[0328] Step S44: Output the population allocation plan in the form of tables and charts. The tables include population control indicators for each region and stage, and the charts use GIS maps to visualize the population distribution targets for each region. At the same time, generate a plan description document, which elaborates on the basis for the plan's formulation, implementation steps, expected effects, and risk response measures.

[0329] In this embodiment, the future water system data output by S2 is organized, and the temporal / spatial matching is verified. The constraints of S3 are converted into mathematical inequalities that the model can recognize. The game equilibrium solution is solved using the interior point method and particle swarm optimization algorithm. The optimal control variables (regional settlement quotas, tiered water price adjustment coefficients, industry approval quotas, and relocation sequence) are obtained through the initial decision-making of the upper-level policy subjects, the response of the lower-level behavioral subjects, and iterative equilibrium. Based on the optimal policy, regional and temporal (short-term / medium-term) population allocation plans are generated, clarifying the control targets for the permanent residents and floating population in each region and the optimized proportion of the industrial workforce. The matching between water demand and water resource carrying capacity is verified through the S2 model, and supplemented with supporting measures such as water supply facilities, traffic optimization, and industrial policies. The plan is output in the form of tables and GIS maps, with an explanatory document explaining the basis, steps, expected effects, and risk response measures.

[0330] According to one aspect of this application, step S42 further comprises:

[0331] Step S42a: Within the robust uncertainty set, the policy-making entity determines the initial policy regulation variables with the objective of minimizing the overall system cost;

[0332] Step S42b: Based on the initial policy, population-industry actors determine the optimal population response vector with the goal of maximizing their overall utility.

[0333] Step S42c: Repeat the solution process for the policy-making subject and the population-industry behavior subject until the policy regulation variable and the population response vector no longer change, and obtain the Stackelberg equilibrium solution, that is, the optimal policy regulation variable;

[0334] Step S42d: Verify the feasibility of the obtained optimal policy control variables to ensure that they meet all constraints. If there is an infeasible solution, adjust the constraint threshold and solve again to obtain the optimal policy control variables. The optimal policy control variables include: regional settlement quotas, tiered water price adjustment coefficients, industrial project approval quotas, and industrial relocation timeline plans.

[0335] In this embodiment, the specific method for solving the equilibrium solution of the robust Stackelberg game model using the interior point method combined with the particle swarm optimization algorithm is as follows:

[0336] This invention employs a two-layer iterative architecture: the outer layer uses Particle Swarm Optimization (PSO) to search for the optimal decision variables of the upper-level policy-making subject, while the inner layer uses the Interior Point Method (IPM) to solve for the optimal response of the lower-level population-industry behavior subject. This design is based on the following considerations: the upper-level decision space contains discrete variables (such as integer constraints on settlement quotas) and a non-convex objective function, making it suitable for metaheuristic algorithms with strong global search capabilities; the lower-level optimization problem is a continuous convex optimization problem given the upper-level decisions, making the interior point method, which has fast convergence speed and high accuracy, suitable for it.

[0337] The decision variable vector of the upper-level policy-making body is π=[π_1, π_2, ..., π_D]^T, which includes regional settlement quotas (D_1 variables), tiered water price adjustment coefficients (D_2 variables), industrial project approval quotas (D_3 variables), and industrial relocation timeline plans (D_4 variables), totaling D=D_1+D_2+D_3+D_4 decision variables.

[0338] The specific parameter settings for the particle swarm optimization algorithm are as follows: number of particles N_p = 80; maximum number of iterations T_max = 60; inertia weight adopts a linear decreasing strategy.

[0339] w^(t)=w_max-(w_max-w_min)·t / T_max;

[0340] Where w_max=0.9, w_min=0.4; cognitive learning factor c_1=2.0; social learning factor c_2=2.0; and the velocity boundary is 20% of the decision variable boundary.

[0341] The formulas for updating the position and velocity of particle i in iteration t are:

[0342] v_i^(t+1)=w^(t)·v_i^(t)+c_1·r_1·(pbest_i-π_i^(t))+c_2·r_2·(gbest-π_i^(t));

[0343] π_i^(t+1)=π_i^(t)+v_i^(t+1);

[0344] Where r_1 and r_2 are uniformly distributed random numbers in the interval [0, 1], pbest_i is the historical best position of particle i, and gbest is the global best position.

[0345] Constraint handling method: The penalty function method is used to handle constraint violations. For the constraint condition g_j(π)≤0 (j=1,2,...,J), the penalty term is defined as follows:

[0346] Penalty(π)=Σ_{j=1}^{J}λ_j·max(0, g_j(π))^2;

[0347] Where λ_j is the penalty coefficient, the initial value is set to 10^4, and it increases adaptively with the number of iterations: λ_j^(t)=λ_j^(0)·(1+0.1·t).

[0348] The corrected fitness function is:

[0349] Fitness(π) = C(W(ξ), P * (π,ξ),π)+Penalty(π);

[0350] Where C is the system comprehensive cost function, P * (π, ξ) is the population-optimal response vector under a given policy π, which is solved by the inner optimization.

[0351] Given the upper-level decision variable π^(t), the optimization problem of the lower-level population-industry behavioral agents can be expressed as:

[0352] max_{P}U_population(P,π^(t),ξ);

[0353] sth_k(P)=0, k=1, 2, ..., K_eq (equality constraints);

[0354] g_l(P)≤0, l=1,2,...,L_ineq (inequality constraints);

[0355] P_min≤P≤P_max (boundary constraint);

[0356] The solution is obtained using the primal-dual interior-point method, and the specific steps are as follows:

[0357] Step a: Introduce slack variables s_l≥0 to transform the inequality constraint into the equality constraint g_l(P)+s_l=0;

[0358] Step b: Construct the barrier function Lagrangian:

[0359] L_μ(P,s,λ,ν)=-U_population(P,π^(t),ξ)+Σ_kλ_k·h_k(P)+Σ_lν_l·(g_l(P)+s_l)-μ·Σ_llog(s_l);

[0360] Where μ>0 is the barrier parameter, and λ and ν are Lagrange multiplier vectors;

[0361] Step c: Solve the KKT system by setting the gradient of L_μ with respect to (P, s, λ, ν) to zero, resulting in a system of nonlinear equations. Solve iteratively using the Newton method. The Newton direction is obtained by solving the following linear system:

[0362] [H_pp0A_h^TA_g^T][ΔP][-▽_PL];

[0363] [0S^20I][Δs]=[-S·ν+μe];

[0364] [A_h000][Δλ][-h(P)];

[0365] [A_gI00][Δν][-g(P)-s];

[0366] Where H_pp is the Hessian matrix, S=diag(s_1,...,s_L), A_h and A_g are the constraint Jacobian matrices, and e is an all-one vector;

[0367] Step d: Perform a line search to determine the step size α_P and α_s, and update the variables:

[0368] P^(k+1)=P^(k)+α_P·ΔP, s^(k+1)=s^(k)+α_s·Δs;

[0369] Simultaneously update multipliers λ^(k+1) and ν^(k+1), with the step size ensuring s^(k+1)>0 and duality feasibility;

[0370] Step e: Gradually reduce the obstacle parameter μ: μ^(k+1)=σ·μ^(k), where σ=0.2 is the reduction factor.

[0371] Step f: Determine the convergence condition: ||▽_PL||_∞<ε_opt and ||h(P)||_∞<ε_feas and μ<ε_μ, then output the optimal solution P. * Otherwise, return to step c. In this embodiment, the convergence thresholds are set to ε_opt=10^(-8), ε_feas=10^(-8), and ε_μ=10^(-10).

[0372] The overall process of upper and lower layer iterations is as follows:

[0373] Step A: Initialize the particle swarm and randomly generate N_p initial policy vectors {π_i^(0)}, i=1,...,N_p;

[0374] Step B: For each particle i, use the interior point method to solve the lower-level optimization problem and obtain the population-optimal response P_i. * =P * (π_i, ξ);

[0375] Step C: Calculate the fitness value (π_i) for each particle, and update pbest_i and gbest;

[0376] Step D: Update the positions of all particles according to the PSO velocity and position update formula;

[0377] Step E: Determine the convergence condition. Stop iterating if any of the following conditions are met:

[0378] (1) The global optimal fitness value has not improved for N_stall=15 generations (improvement threshold ε_improve=10^(-6));

[0379] (2) The maximum number of iterations T_max = 60 is reached;

[0380] (3) The change in the policy regulation variable satisfies ||gbest^(t)-gbest^(t-1)||_∞ / ||gbest^(t-1)||_∞<ε_π=0.1%.

[0381] Step F: Output the global optimal solution gbest as the optimal policy adjustment variable π. * The corresponding P * This represents the optimal response vector for the population.

[0382] The specific method for determining "the policy regulation variables and the population response vector no longer change" is as follows:

[0383] For policy regulation variables, the relative change Δπ_rel is defined as ||π^(t)-π^(t-1)||_2 / ||π^(t-1)||_2, and the threshold is ε_π=0.1%, that is, when Δπ_rel<0.001, the policy regulation variables are considered to converge.

[0384] For the population response vector, the relative change ΔP_rel = ||P * ^(t)-P * ^(t-1)||_2 / ||P * ^(t-1)||_2, the decision threshold is ε_P=0.1%, that is, when ΔP_rel<0.001, the population response vector is considered to be converged;

[0385] Stackelberg equilibrium is considered to be achieved when both Δπ_rel < ε_π and ΔP_rel < ε_P are satisfied.

[0386] In this embodiment, the above algorithm is used to solve the population-water resources-policy game problem in City A. After 47 outer iterations, the problem converges and the final global optimal fitness value is 3.82×10^7 (system comprehensive cost, unit: 10,000 yuan). The relative change of the policy regulation variable Δπ_rel=0.08% and the relative change of the population response vector ΔP_rel=0.06%, which satisfies the Stackelberg equilibrium criterion.

[0387] In this embodiment, interior point method + particle swarm optimization is selected.

[0388] The combined solution approach can address both the single interior point method, which excels at continuous constraint optimization but has weak global search capabilities, and the single particle swarm optimization algorithm, which excels at global search but has low accuracy in constraint handling. The combined approach can simultaneously address the optimal solution for a single scenario and the worst-case scenario search within a robust set of uncertainties, thus efficiently obtaining the game equilibrium solution.

[0389] The water resource carrying capacity and population status vary greatly in different areas of the study area (e.g., District A of City A has abundant water resources but a dense population, while County D has scarce water resources but a sparse population). Differentiated regulation can be achieved by region; the short-term and medium-term development goals are different, and the phased and coherent nature of the plan can be ensured by dividing the time sequence.

[0390] Traditional solutions are often difficult to implement due to a lack of supporting measures. Supplementing these measures with water supply, transportation, and industry can ensure that population allocation plans are coordinated with infrastructure and industrial development, thereby improving their feasibility.

[0391] In one embodiment, specifically:

[0392] The population forecast data and water system status data of each street / township in City A for 2024-2028 and 2029-2033 were organized by time period and spatial region to form a data table (e.g., the predicted resident population of Street-1 in District A for 2024-2028 is 88,000-92,000).

[0393] The verification data time step (annual) matches the decision cycle of the game model, and the spatial range is consistent with the administrative boundary of City A. The missing population forecast data for a certain township in Ding County in 2032 is completed.

[0394] The constraints of water system, population-industry behavior, and policy are transformed into mathematical inequalities, such as "the permanent population of area A is ≤350,000" being transformed into P_area A ≤350,000.

[0395] The solution is obtained by using the interior point method combined with particle swarm optimization algorithm, with the number of particles set to 80, the number of iterations to 60, and the inertia weight to 0.6.

[0396] Upper-level solution: The policy-making body is within a robust uncertainty set, with the goal of minimizing the overall system cost. The initial settings are: an annual household registration quota of 3,000 people in District A and a tiered water price adjustment coefficient of 1.1.

[0397] Lower-level solution: Based on the initial policy, the population and industry actors maximize the overall utility and determine the annual migration volume from County Ding to District Jia to be 4,500 people;

[0398] Equilibrium Iteration: Repeat the solution of the upper and lower layers until an equilibrium solution is obtained;

[0399] Verify that the equilibrium solution satisfies all constraints (such as the quota for settling in District A being 3200 people ≤ the upper limit of authority, and the fiscal cost being 42 million yuan ≤ the budget).

[0400] Optimal policy control variables include:

[0401] Settlement quotas by region: 3,200 people / year in Zone A, 2,800 people / year in Zone B, 2,500 people / year in Zone C, 1,800 people / year in County D, and 1,500 people / year in County E;

[0402] Tiered water pricing adjustment coefficients: 1.1 for residential water use, 1.2 for industrial water use;

[0403] Approval quotas for industrial projects: 1 project per year for high water-consuming industries and 5 projects per year for low water-consuming high-end manufacturing industries;

[0404] Industrial relocation timeline: The relocation of two water-intensive enterprises in Ding County to the Bing District Industrial Park will be completed before 2025.

[0405] According to one aspect of this application, step S43 further comprises:

[0406] Step S43a: Based on the optimal policy control variables and the optimal population response vector, generate population allocation schemes by region and time series, including: population control targets for each street, township and community in the study area, phased targets for the future time, control scale of permanent population, control threshold of floating population and optimized proportion of industrial employment population.

[0407] Step S43b: Use the human-water joint prediction model to calculate the water system carrying capacity status after the implementation of the population allocation plan, and verify whether the generated population allocation plan meets the water resource carrying capacity threshold constraints of each region. If there is a risk of water shortage, adjust the population allocation scale or policy regulation variables.

[0408] Step S43c, supplementary suggestions for safeguards to ensure the implementation of the plan, including plans for supporting water supply facilities, suggestions for optimizing transportation and commuting facilities, and detailed rules for industry support and restrictive policies.

[0409] In one embodiment, specifically:

[0410] Based on optimal policies and population response vectors, regional and temporal population allocation schemes are generated, with core elements including:

[0411] Short term (2024-2028): The permanent resident population of District A will be controlled at 320,000-330,000 (with the proportion of migrant population at 18-19%), and the permanent resident population of County D will be controlled at 250,000-260,000 (the proportion of industrial employment will be optimized: the secondary industry will be reduced to 35%).

[0412] Medium term (2029-2033): The permanent resident population of Area A will be controlled at 330,000-350,000 (with the proportion of floating population at 19-20%), and the permanent resident population of County D will be controlled at 240,000-250,000 (the proportion of industrial employment will be optimized: the secondary industry will be reduced to 30%).

[0413] The human-water joint diagram time series prediction model was used for verification. After the implementation of the plan, the water demand in Area A in 2028 is estimated to be 78,000 m³. 3 / d, less than the available water resources of 85,000 m³ 3 / d, no risk of water shortage;

[0414] Supplementary safeguards:

[0415] Water supply infrastructure: The expansion and renovation of the booster pump station in Zone A will be completed by 2026 (adding a water supply capacity of 20,000 m³). 3 / d);

[0416] Traffic commuting optimization: By 2027, the rail transit will be extended to the C District Industrial Park to improve commuting efficiency between D County and C District;

[0417] Detailed industrial policy rules: Tax incentives will be given to high-end manufacturing industries with low water consumption (tax rate reduced by 2%), and a water resource surcharge will be levied on high water-consuming industries (0.5 yuan per ton of water).

[0418] Output population control indicators for each region and stage in tabular form, and visualize the population distribution targets for each region using GIS maps; generate a scheme description document, which explains the basis for formulation, implementation steps, expected effects (such as reducing the water resource carrying capacity of City A from the current 85% to 75% in 2028) and risk response measures (such as activating the Ding County emergency water supply project in the event of extreme drought).

[0419] Combining interannual hydrological fluctuations with population change rates, a one-year period is preset as the standard closed-loop cycle (which can be shortened to six months during periods of rapid development or large hydrological fluctuations). At each cycle node, measured population, water system, meteorological, economic, and policy data from the previous cycle are collected to complete data quality control. The latest data is substituted into the S3 parameter update mechanism to optimize model parameters. The S2 prediction model training set is updated and retrained using the rolling window method. Based on the updated model, the optimized policy and population allocation schemes are solved. The indicators (water resource carrying capacity and human-water matching degree) before and after optimization are compared, and the updated scheme is released after expert review. Non-periodic extreme hydrological events and sudden public events trigger the emergency mechanism, immediately initiating data collection, parameter updates, and scheme optimization without waiting for the standard cycle.

[0420] The hydrological data (precipitation, runoff) in the study area show significant interannual fluctuations, while the population statistics are mostly annual. A one-year cycle can balance the effectiveness of data updates with the stability of the scheme, avoiding the excessive cost of frequent adjustments caused by a cycle that is too short (6 months) or the lag in the scheme caused by a cycle that is too long (2 years).

[0421] Fixed training sets are prone to model adaptation decline due to outdated data. The rolling window method can retain the latest data and remove outdated data to ensure that the prediction model always adapts to the latest human-water system characteristics.

[0422] Non-periodic factors such as extreme hydrological events and public emergencies can quickly disrupt the water-human balance, and conventional periodic optimization cannot respond in time. Emergency response mechanisms can make rapid adjustments, reduce risks such as water shortages and overcrowding, and ensure the stable operation of the system.

[0423] In one embodiment, specifically:

[0424] Based on the interannual fluctuation characteristics of hydrology in City A (precipitation cycle of about 1 year) and the rate of population change, the preset closed-loop cycle is set to 1 year;

[0425] At the end of 2024 (the first closed-loop cycle node), actual measured data for City A in 2024 were collected, including the permanent resident population (actual 89,000 people in Street-1 of District A, predicted 88,000 people) and floating population data of each street / township, water supply and water quality data of each water system node, precipitation data for 2024 (10% higher than normal), GDP growth rate (6.2%), and data on the adjustment of the new household registration policy in 2024; consistency verification of the data was performed, and population data for a certain township in County E in the fourth quarter of 2024 was completed;

[0426] The measured data from 2024 were substituted into the S3 parameter update mechanism to update the population response function weight coefficients (α=0.33, β=0.18, γ=0.32, δ=0.17) and the robust uncertainty set boundary (Wasserstein distance adjusted to 0.11). The rolling window method was used to update the S2 joint prediction model training set (2010-2024). After retraining, the validation set MAE decreased to 4.0%.

[0427] The updated prediction model output data for 2025-2029 and 2030-2034 are input into the updated game model, and the solution is re-solved to optimize the optimal policy control variables (such as adjusting the settlement quota for District A to 3,300 people in 2025) and population allocation plan (controlling the permanent resident population of District A to 325,000-335,000 people from 2025 to 2029).

[0428] Comparing the plans before and after optimization, the water resource carrying capacity decreased from 75% to 73%, and the population-water system matching degree improved by 3%; five experts in water conservancy, population, and planning were invited to review the plan, and its scientific validity and feasibility were approved; the updated population allocation plan and policy adjustment notice for 2025 were released.

[0429] In the summer, City A experienced a severe drought (an extreme hydrological event), triggering an abnormal response mechanism. The city immediately collected measured precipitation and water supply data for the first half of the year, updated model parameters, and optimized policy control variables (the quota for settling in District A that year was temporarily adjusted to 2,800 people) and population allocation plans to ensure a balance between people and water.

[0430] This application utilizes multi-dimensional technological innovation, a human-water coupled topology structure and joint prediction model to accurately capture the two-way coupling relationship between humans and water. A robust Stackelberg game framework covers multi-source uncertainties, ensuring human-water balance in extreme scenarios. Population and industry are defined as active players, and migration decisions are accurately captured through a comprehensive utility function, achieving two-way interaction between policy and population behavior. Annual updates are combined with emergency response to continuously optimize the model and solutions, adapt to dynamic changes in the system, and improve the timeliness and operability of the solutions.

[0431] In other alternative implementations, the graph neural network can also be replaced by a graph convolutional network (GCN), GraphSAGE, or graph transformer; the temporal prediction network can also be replaced by a gated recurrent unit (GRU), temporal convolutional network (TCN), or transformer; and the optimization algorithm can also be replaced by a genetic algorithm, simulated annealing algorithm, or differential evolution algorithm.

[0432] The selection of the above alternatives can be adjusted according to the size of the study area, data characteristics, and computing resources, without affecting the core technical idea of ​​the present invention.

[0433] According to another aspect of this application, a system for determining population size based on water availability is provided, characterized by comprising:

[0434] At least one processor; and

[0435] A memory communicatively connected to at least one of the processors; wherein,

[0436] The memory stores instructions that can be executed by the processor to implement the "water-based population determination" method described above.

[0437] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details of 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. The method of "determining population size based on water availability" is characterized by: Includes the following steps: Step S1: Collect data of the study area and construct the human-water coupling topology of the study area; Step S2: Construct a human-water joint graph time series prediction model based on the human-water coupled topology, and solve for the future time series population distribution data and water system state data; Step S3: Construct a robust Stackelberg game-based multi-agent coupling model of population, water resources, and policy, and update the model parameters; Step S4: Input the future time series data into the updated model and solve for the population allocation scheme of the study area.

2. The method of "determining population size based on water availability" as described in claim 1, characterized in that, Step S1 further comprises: Step S11: Collect historical population distribution data, water system status data, and multi-source basic data for the study area. The multi-source basic data includes: hydrological and water resources data, population and industry data, infrastructure data, spatial geographic data, and policy and economic support data. Step S12: Based on the multi-source basic data of the study area, a graph structure modeling method is used to construct the human-water coupled topology of the study area. The human-water coupled topology includes: population nodes, water system nodes, and associated edges connecting the two types of nodes.

3. The method of "determining population size based on water availability" as described in claim 2, characterized in that, Step S12 further comprises: Step S12a: Divide the study area into nodes using streets, towns, or communities as basic units; Step S12b: Define the surface water node, groundwater node, and water supply facility node respectively; Step S12c: Construct internal edges of the water system based on the water pipeline network and river flow path, and connect internal nodes of the water system and water system nodes with population nodes; Step S12d: Construct population flow-commuting connection edges based on transportation network connectivity and commuting flow data; Step S12e: Based on the actual water supply and demand relationship and population commuting characteristics in the study area, optimize the weight assignment of the associated edges, and finally form a human-water coupled topology structure containing node set, edge set and complete attribute information.

4. The method of "determining population size based on water availability" as described in claim 1, characterized in that, Step S2 further comprises: Step S21: Construct an integrated model architecture that combines graph neural networks and time-series prediction networks to obtain a human-water joint graph time-series prediction model; Step S22: Extract the attribute data of each node in the human-water coupling topology to form an initial node feature matrix, map the feature data to the [0,1] interval, and generate a standardized node feature matrix; Step S23: Extract multi-source basic data and perform unified time alignment to obtain time series samples of a preset length, and generate exogenous feature time series matrix; Step S24: Concatenate the standardized node feature matrix and the exogenous feature time series matrix along the time dimension to obtain the model input data matrix; Step S25: Train the human-water joint graph time series prediction model based on the model input data matrix, and solve the model to obtain the population distribution data and water system state data for future time series.

5. The method of "determining population size based on water availability" as described in claim 4, characterized in that... Step S25 further comprises: Step S25a: Randomly divide the input data matrix into a training set, a validation set, and a test set in a ratio of 7:2:1; Step S25b: Set the batch size and learning rate for model training, use the Adam optimizer to minimize the mean squared error loss function between the predicted and true values, and set the maximum number of training rounds; Step S25c: Based on the validation set data, the mean absolute error and root mean square error are used to evaluate the prediction accuracy of the model. If the prediction accuracy does not reach the preset threshold, the convolution kernel size of the graph neural network, the time step size of the temporal prediction network and the number of network layers are adjusted, and iterative training is performed until the prediction accuracy reaches the preset threshold to obtain the trained human-water joint graph temporal prediction model. Step S25d: Load the trained human-water joint graph time series prediction model into the inference environment, input the preprocessed complete input data matrix, and solve for the future time series population distribution data and water system state data through forward propagation.

6. The method of "determining population size based on water availability" as described in claim 1, characterized in that, Step S3 further comprises: Step S31: Taking the water system, population-industry behavior, and policy-making as the main players in the game, and determining their respective decision-making authority, construct the objective functions for the water system player, the population-industry behavior player, and the policy-making player respectively; Step S32: Set water system constraints, population-industry behavior constraints, and policy-making constraints respectively; Step S33: Introduce uncertain variables, calculate the initial probability distribution of each uncertain variable based on historical data of the study area, and construct a robust uncertainty set with the initial distribution as the center and the Wasserstein distance as the radius; Step S34: Update the population response function and robust uncertainty set boundary of the model based on historical population distribution data and water system status data of the study area.

7. The method of "determining population size based on water availability" as described in claim 6, characterized in that, Step S34 further comprises: Step S34a: Extract historical population distribution data and water system status data for the study area; Step S34b: Update the weight coefficients in the population response function using the maximum likelihood estimation method to maximize the likelihood value between the measured population state and the model-predicted population state. Step S34c: Optimize the boundary of the robust uncertainty set by adjusting the Wasserstein distance threshold so that the set contains more than 90% of the measured uncertainty scenarios; Step S34d: Use the analytic hierarchy process (AHP) combined with expert scoring to update the weight coefficients of each objective function; Step S34e: Substitute the updated parameters into the game model, and use measured data to verify the deviation between the policy regulation variables output by the model and the actual policy effect. If the deviation exceeds the preset threshold, repeat the parameter update process until the deviation meets the requirements, and complete the model parameter update.

8. The method of "determining population size based on water availability" as described in claim 1, characterized in that, Step S4 further comprises: Step S41: Extract future time series population distribution data and water system status data, and input them into the updated population-water resources-policy multi-agent coupling model; Step S42: Solve for the equilibrium solution of the robust Stackelberg game model to obtain the optimal policy control variables; Step S43: Based on the optimal policy control variables and the optimal population response vector, generate population allocation schemes by region and time series; Step S44: Output the population allocation plan in the form of tables and charts, and generate a plan description document, including the basis for the plan, implementation steps, expected effects and risk response measures.

9. The method of "determining population size based on water availability" as described in claim 8, characterized in that, Step S42 further comprises: Step S42a: Within the robust uncertainty set, the policy-making entity determines the initial policy regulation variables with the objective of minimizing the overall system cost; Step S42b: Based on the initial policy, population-industry actors determine the optimal population response vector with the goal of maximizing their overall utility. Step S42c: Repeat the solution process for the policy-making subject and the population-industry behavior subject until the policy regulation variable and the population response vector no longer change, and obtain the Stackelberg equilibrium solution, that is, the optimal policy regulation variable.

10. The method of "determining population size based on water availability" as described in claim 8, characterized in that, Step S43 further comprises: Step S43a: Based on the optimal policy control variables and the optimal population response vector, generate population allocation schemes by region and time series; Step S43b: Use the human-water joint prediction model to calculate the water system carrying capacity status after the implementation of the population allocation scheme, and verify whether the generated population allocation scheme meets the water resource carrying capacity threshold constraints of each region. Step S43c, supplementary suggestions for safeguards to ensure the implementation of the plan, including plans for supporting water supply facilities, suggestions for optimizing transportation and commuting facilities, and detailed rules for industry support and restrictive policies.

11. The system of "determining population size based on water availability" is characterized by: include: At least one processor; as well as A memory communicatively connected to at least one of the processors; wherein, The memory stores instructions that can be executed by the processor to implement the "water-based population determination" method as described in any one of claims 1 to 10.

Images