Urban Reclaimed Water Supply and Demand Coordination Forecasting Method

Abstract

This invention discloses a method for coordinated prediction of urban reclaimed water supply and demand, comprising: processing multi-source heterogeneous data streams from reclaimed water pipeline networks to generate a standardized spatiotemporal dataset; constructing a dynamic graph network prediction model based on the dataset, exploring multiple hydraulic transport paths between nodes, and constructing a bimodal path transmission kernel for each path that can respectively characterize physical transmission delay and decision response delay, and then weighting and superimposing them to form a comprehensive transmission kernel to accurately characterize the spatiotemporal transmission characteristics of water quality; using the model and real-time demand signals to generate quality-specific demand prediction results and cross-plant collaborative production schemes. This invention, by finely characterizing the hydraulic-water quality coupling transmission law and combining it with efficient collaborative optimization strategies, can improve the accuracy of prediction, the economy of scheduling, and the safety and value of reclaimed water resource utilization.

Description

Technical Field

[0001] This invention belongs to the field of urban reclaimed water utilization, and in particular to a method for coordinating the supply and demand forecasting of urban reclaimed water. Background Technology

[0002] The large-scale utilization of reclaimed water is not simply about transporting it from treatment plants to users; it is a complex systems engineering project involving multi-source production, multi-level quality, diverse users, and dynamic demand. Accurately predicting the supply and demand of reclaimed water of different qualities, and on this basis, achieving coordinated production among multiple water plants and optimized scheduling of the pipeline network to meet the diverse water needs of cities with the lowest cost, highest efficiency, and safest guarantee, is key to maximizing the value of urban reclaimed water utilization and a core technological issue driving the transformation and upgrading of urban water systems towards refinement, intelligence, and green development.

[0003] Currently, researchers have conducted extensive work on the forecasting and scheduling of urban water supply systems. In demand forecasting, existing methods primarily rely on time series analysis models, such as traditional statistical methods like the Autoregressive Integrated Moving Average (ARIMA) model and exponential smoothing, as well as emerging machine learning methods such as Support Vector Machines (SVM), Artificial Neural Networks (ANN), and Long Short-Term Memory (LSTM). These methods are effective at capturing the inherent trends and periodic patterns in historical water consumption data for individual users or regions. In pipeline network scheduling and hydraulic simulation, the mainstream technology is simulation optimization based on hydraulic models such as EPANET. These models can accurately calculate hydraulic parameters such as pressure and flow rate in the pipeline network and, through combination with intelligent optimization algorithms such as genetic algorithms and particle swarm optimization, optimize scheduling schemes such as pump station start-up and shutdown, and valve switching. In water quality simulation, research focuses on the one-dimensional transport and first-order kinetic decay processes of single pollutants (such as residual chlorine) in the pipeline network, providing some technical support for ensuring water supply security.

[0004] While existing technologies have improved the scientific nature of water management to some extent, they still face a series of deep-seated technical bottlenecks when dealing with the complex system of urban reclaimed water, which is multi-source, multi-quality, and strongly coupled. These bottlenecks collectively lead to the practical dilemmas of inaccurate prediction, poor coordination, and low efficiency. In summary, existing technologies mainly suffer from problems such as distorted representation of hydraulic-water quality coupling transmission characteristics, an imbalance between computational efficiency and security in complex network models, and limitations in the decision-making dimensions of multi-objective collaborative scheduling strategies. Summary of the Invention

[0005] The purpose of this invention is to provide a method for coordinating the supply and demand forecasting of urban reclaimed water, so as to solve the above-mentioned problems existing in the prior art.

[0006] Technical solutions, including methods for coordinated forecasting of urban reclaimed water supply and demand, include:

[0007] Process multi-source heterogeneous data streams from reclaimed water pipeline networks to generate standardized spatiotemporal datasets;

[0008] Based on standardized spatiotemporal datasets, a dynamic graph network prediction model is constructed to characterize the supply and demand relationship and spatiotemporal transmission characteristics of water quality among nodes in the pipeline network.

[0009] By using dynamic graph network prediction models and real-time demand signals, we can generate quality-specific demand prediction results and cross-plant collaborative production plans.

[0010] Beneficial effects: By finely characterizing the hydraulic-water quality coupling transmission law and combining it with an efficient collaborative optimization strategy, this invention can improve the accuracy of prediction, the economy of scheduling, and the safety and value of reclaimed water resource utilization. Attached Figure Description

[0011] Figure 1 A flowchart illustrating the steps of the urban reclaimed water supply and demand coordinating forecasting method provided in this application embodiment.

[0012] Figure 2 A flowchart illustrating the steps for constructing a dynamic graph network prediction model as provided in this application embodiment.

[0013] Figure 3 A flowchart illustrating the steps for constructing a path conduction core for each hydraulic transport path, as provided in this application embodiment.

[0014] Figure 4 A flowchart illustrating the steps for determining the time-delayed conduction core provided in this application embodiment. Detailed Implementation

[0015] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0016] It should be noted that the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion, for example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0017] The study revealed distortions in the characterization of hydraulic-water quality coupling transmission characteristics. Existing models generally employ a single-delay assumption based on the shortest path and average flow velocity to describe the water transport process from source to sink. However, in real urban pipe networks, especially ring networks, water flows simultaneously along multiple paths. The single-path assumption deviates significantly from physical reality, leading to substantial errors in the prediction of water quality arrival time. Furthermore, existing models neglect a crucial non-physical delay—the decision-response delay, which is the time required for a series of human or systemic decision-making processes from receiving scheduling instructions to the water plant actually adjusting production and opening valves to deliver water. This systematic neglect of multipath effects and decision delays prevents the models from accurately depicting the real spatiotemporal dynamic distribution of water quality signals propagating in the pipe network, directly causing the failure of precise time matching between supply and demand. Secondly, the aforementioned characterization distortions further exacerbate the imbalance between computational efficiency and security in complex network models. On the one hand, to more realistically simulate pipe networks, models need to include tens of thousands of nodes and connections, leading to the curse of dimensionality. The computational complexity of fully connected models increases exponentially, making it difficult to meet the efficiency requirements of real-time scheduling. On the other hand, when simplifying models (such as sparsification), existing methods lack domain knowledge guidance and may incorrectly retain or sever connections. In particular, they fail to adequately consider the potential risks of chemical incompatibility or quality degradation when mixing reclaimed water of different qualities. This blind sparsification process may result in a theoretically optimal scheduling scheme, but in reality, it could lead to serious water quality safety incidents. Finally, the above two problems together lead to the limitation of decision-making dimensions in multi-objective collaborative scheduling strategies. Traditional scheduling optimization is often single-objective or low-dimensional, such as only pursuing the lowest cost, and is usually carried out independently between individual water plants, lacking a global collaborative perspective. It fails to take the unique quality dimension of reclaimed water as a core optimization variable and cannot effectively utilize the differences in water quality requirements of different users to build a tiered utilization network that makes the most of water. At the same time, facing an ultra-high-dimensional decision space across plants, regions, qualities, and time periods, traditional centralized optimization solvers struggle to find the global optimum within a reasonable timeframe. These limitations ultimately lead to resource waste in using high-quality reclaimed water in low-requirement scenarios and make it difficult to improve the overall operational efficiency of the system.

[0018] In some embodiments, the implementation of the urban reclaimed water supply and demand coordinated forecasting method can be based on a system integrating data acquisition, transmission, storage, and computation. This system may include online monitoring equipment (e.g., flow meters, pressure sensors, water quality analyzers) deployed at various key locations in the reclaimed water network, a distributed data acquisition and monitoring system (SCADA system), a production scheduling system for each reclaimed water plant, a user feedback platform, and a central server or cloud computing platform for executing the method described in this application.

[0019] like Figure 1As shown, the method for coordinated forecasting of urban reclaimed water supply and demand includes the following steps:

[0020] Process multi-source heterogeneous data streams from the reclaimed water network to generate a standardized spatiotemporal dataset.

[0021] In other words, we acquire multi-source heterogeneous data streams from the reclaimed water pipeline network, preprocess them, and generate standardized spatiotemporal datasets.

[0022] For example, multi-source heterogeneous data streams can include real-time pipeline monitoring data from SCADA systems, production scheduling records from water plants, topological data from urban pipeline GIS systems, and historical water usage records from users in various fields. The process of generating standardized spatiotemporal datasets aims to transform these data, which come from different sources, have varying formats, and inconsistent spatiotemporal dimensions, into high-quality, structured data that can be directly used by subsequent models. Specifically, this step can include the following processing steps: Acquiring and cleaning the real-time monitoring data stream. From the distributed SCADA system, raw monitoring data streams containing water quality parameters such as flow rate, pressure, COD (chemical oxygen demand), turbidity, ammonia nitrogen, pH, residual chlorine, and conductivity can be extracted. To handle data distortion caused by sensor failures or communication anomalies, a sliding window-based outlier detection algorithm can be used to identify abnormal data points. For identified outliers, instead of simply discarding them, data from their nearest normal time point is used to repair them using cubic spline interpolation. Simultaneously, the system marks the abnormal time period for subsequent fault analysis. After this processing, a cleaned monitoring data sequence with a 5-minute time interval and spatially covering all pipeline nodes can be formed. The system integrates production operation data in a structured manner. For example, it can integrate production scheduling records from three reclaimed water plants, which include information such as capacity, water quality level switching, and operating costs for different time periods. Since the data formats of different plants may be heterogeneous, a unified data model needs to be constructed to convert text-based scheduling instructions (e.g., Plant C switches to Q1 high-quality mode at 8:00 AM) into a structured capacity time series. Specifically, for Plant C, which can switch between different water quality levels, the switching process is not instantaneous. To avoid the impact of abrupt changes on the prediction model, the system can refine the switching process into a gradual capacity change curve during the transition period. After integration, a standardized production capacity matrix with hourly granularity can be generated, where rows represent timestamps and columns represent the capacity of each plant at different water quality levels.

[0023] The system dynamically maps the pipeline network topology. It parses static topological data from the urban pipeline network GIS system, such as node coordinates, pipe segment connections, pipe diameter, and pipe age. Simultaneously, to reflect the actual state of the pipeline network, the system incorporates pipeline renovation and construction plans provided by engineering departments, thus constructing a time-varying pipeline network topology representation. For example, changes during each pipeline renovation period can be encoded as a sequence of addition and deletion operations in a graph structure. For temporary pipelines laid during construction, due to limitations in materials and laying conditions, a standard transport capacity coefficient of, for example, 0.6 times can be assigned. Through this topology evolution algorithm, the system can generate a dynamic pipeline network topology sequence covering the entire prediction time domain, where each time slice contains the effective set of nodes, connected edges, and their respective transport capacity parameters at that moment. The system also performs hierarchical extraction of historical demand patterns. It aggregates historical water usage records from users in different sectors such as industry, agriculture, municipal services, and landscaping to form a raw demand dataset. To gain a deeper understanding of water use behavior, an Empirical Mode Decomposition (EMD) approach can be employed to decompose the demand sequences across various sectors into trend components (reflecting long-term changes), periodic components (such as daily, weekly, and seasonal fluctuations), and stochastic components (reflecting short-term disturbances). The periodic components can be further analyzed using Fourier transforms to extract the amplitude and phase of each periodic component. For different water quality levels (e.g., Q1 for high-standard industrial use, Q2 for agricultural irrigation, and Q3 for municipal miscellaneous uses), the system calculates the usage proportions and temporal distribution characteristics across different sectors, ultimately constructing a hierarchical demand pattern feature library. This feature library can contain, for example, 128-dimensional feature vectors to provide a refined description of the water use behavior patterns of various user groups.

[0024] All processed data are standardized in both spatiotemporal dimensions. Since data from different sources have varying sampling frequencies and spatial resolutions, the system needs to align them to a unified benchmark. For example, bilinear interpolation can be used to unify data with different sampling frequencies (such as hourly production data) to a 5-minute benchmark time granularity. Spatially, for discrete monitoring point data, inverse distance weighted interpolation can be used to project their values ​​onto all nodes of the pipeline network diagram, thus obtaining a fully covered spatial data field. The final output standardized spatiotemporal dataset can have a four-dimensional tensor [time × space × attribute × water quality level], serving as a unified and standardized input for all subsequent prediction models. Through these steps, the raw, chaotic, multi-source heterogeneous data is transformed into a high-value dataset with a clear structure, spatiotemporal alignment, and complete information, laying a solid foundation for subsequent accurate prediction and collaborative optimization.

[0025] Based on standardized spatiotemporal datasets, a dynamic graph network prediction model is constructed to characterize the supply and demand relationship and spatiotemporal transmission characteristics of water quality among nodes in the pipeline network.

[0026] In this embodiment, the dynamic graph network prediction model is a mathematical abstraction and dynamic simulation of a real urban reclaimed water pipeline network system. Specifically, the intersection of water plants, users, and the pipeline network can be mapped as nodes in a graph, and each node is assigned corresponding attributes, such as node type (production node, consumption node, transfer node) and real-time status (capacity, demand, water quality parameters, etc.). The pipe segments connecting these nodes are mapped as directed edges in the graph, and the edge weights can comprehensively consider factors such as the physical length, diameter, and wall roughness of the pipe segment. The model is dynamic because it can reflect real-time changes in the pipeline network topology. For example, when a pipe segment enters a construction phase, the model will automatically remove the corresponding edge and may activate a preset temporary alternative path. The core of this model lies in characterizing the spatiotemporal conduction characteristics of water quality. In other words, the model not only needs to know whether the nodes are connected, but more importantly, it needs to be able to predict how long it will take and in what form (e.g., whether the concentration will decrease due to diffusion) water of a specific quality flowing from one node (such as a water plant) will reach another node (such as a user). To achieve this function, a time-delayed propagation kernel can be introduced to quantitatively describe the spatiotemporal dynamic relationship of water quality propagation between any two nodes.

[0027] By using dynamic graph network prediction models and real-time demand signals, we can generate quality-specific demand prediction results and cross-plant collaborative production plans.

[0028] In this embodiment, the model's predictive capabilities are translated into actual production scheduling decisions. The system integrates the spatiotemporal correlation structure provided by the dynamic graph network prediction model and real-time collected demand signals (e.g., an industrial park needs a large amount of Q1-grade reclaimed water in the next 2 hours) to construct a multi-objective optimization model that considers time-delay constraints. The objectives of the optimization model are multi-dimensional, and may include: minimizing total production costs; maximizing supply-demand matching, i.e., meeting the needs of all users as much as possible; optimizing water quality grade matching to avoid waste caused by using high-quality water in low-requirement scenarios; and balancing pipeline transport load to prevent excessive pressure in some pipelines. By solving this large-scale optimization problem, the system can ultimately output two types of core results: one is the quality-specific demand prediction results, for example, predicting that the industrial sector's demand for Q1-grade water is X cubic meters / hour and the agricultural sector's demand for Q2-grade water is Y cubic meters / hour in the next 24 hours, etc. Second, there is a cross-plant collaborative production plan, which includes specific scheduling instructions for each water plant. For example, it specifies when and at what power plant Plant A should produce Q1 grade water, how Plants B and C should coordinate, and even when Plant C should switch water quality levels, so as to meet the overall forecasted demand with the highest efficiency and lowest cost.

[0029] like Figure 2 As shown, according to one aspect of this application, a dynamic graph network prediction model is constructed, comprising:

[0030] Based on the dynamic pipeline topology in a standardized spatiotemporal dataset, a predetermined hydraulic transport path between pipeline node pairs is explored and identified.

[0031] Furthermore, it can be done by constructing a time-varying graph structure based on the dynamic network topology in a standardized spatiotemporal dataset.

[0032] Specifically, at each prediction time t, the system generates a directed graph representation G(t) based on the dynamic pipeline topology sequence. In this graph, water plants, large users, and pipeline junctions are abstracted as graph nodes, each assigned a node type attribute (e.g., production node, consumption node, transfer node) and a state vector containing real-time water quality parameters, production capacity, or demand. Pipe segments are abstracted as directed edges, with edge weights considering factors such as physical length, pipe diameter, and hydraulic parameters based on pipe age and material wall roughness degradation coefficient. When the standardized spatiotemporal dataset contains pipeline construction information, such as a pipe segment being closed at time t, the algorithm automatically removes the edge and may activate a pre-defined temporary alternative path. The transport capacity of this alternative path is calculated as a certain proportion (e.g., 60%) of the original path.

[0033] Based on the time-varying graph structure, multiple hydraulic transport paths between pipeline node pairs are explored and identified.

[0034] In this embodiment, to overcome the limitations of traditional models that rely on a single shortest path and thus more realistically reflect the water flow distribution in complex ring networks, the system executes an improved Yen algorithm for each source-target node pair (s, d) (e.g., a water plant s and an industrial user d). This algorithm finds a shortest path using methods such as Dijkstra's algorithm, iteratively removes edges from the found path, and replans detours to systematically search for the top K candidate paths. Here, K is a preset integer, e.g., K≤5, to achieve a balance between computational complexity and path coverage. Furthermore, to ensure the identified paths are physically feasible, the system performs hydraulic feasibility screening on these K candidate paths. Specifically, it calculates the total hydraulic gradient for each path and eliminates paths that would cause the pressure along the route to fall below the minimum service pressure requirement. After screening, a set of all physically reachable valid paths between each node pair is obtained.

[0035] A path conduction kernel is constructed for each hydraulic transport path. The path conduction kernel exhibits a bimodal distribution to characterize the physical transmission delay and the decision response delay, respectively.

[0036] Among them, the path propagation kernel K p(τ) is a probability density function with respect to time delay τ, used to describe the probability distribution of water quality signals originating from the source node at time t and arriving at the target node at time t+τ. In this embodiment, to simultaneously capture the physical process delay of water flowing in the pipeline and the response process delay of water system scheduling decisions, a bimodal gamma distribution is used to model this conduction kernel. Specifically, K... p (τ) = w 1p •Γ(τ;k 1p θ 1p ) + w 2p •Γ(τ;k 2p θ 2p ); where Γ(τ; k, θ) is the gamma probability density function; τ is the time delay, τ>0; k is the shape parameter; θ is the scale parameter; w 1p and w 2p Let w be the weighting coefficient of the two peaks. 1p +w 2p =1. The first peak Γ(τ;k) 1p θ 1p The shape parameter k is used to characterize the physical propagation delay, and its peak position roughly corresponds to the average transit time of a water particle along path p. 1p and scale parameter θ 1p Based on the principles of fluid dynamics, and according to the total length L of path p, p Average flow velocity v p and Reynolds number Re p To determine. For example, for turbulent conditions, k 1p It can be set to 2.2 + 0.12*L p / v p θ 1p It can be set to 1.1 + 0.06*L p / v p The second peak Γ(τ;k) 2p θ 2p The term "delay" is used to characterize decision-making response delay, reflecting the time required from the water plant receiving dispatch instructions or market demand signals to the actual adjustment of production and valve opening for water supply—a series of human or systemic decision-making processes. This delay is not caused by the physical transportation itself, but it is prevalent and cannot be ignored in real supply and demand responses.

[0037] Furthermore, such as Figure 3 As shown, a path transmission core is constructed for each hydraulic transport path, including:

[0038] Extract the water quality grades of reclaimed water associated with the path from the standardized spatiotemporal dataset; based on the water quality grades, determine the parameters of the peaks in the bimodal distribution that characterize the delay in decision response.

[0039] In this embodiment, the parameter k of the decision response delay peak 2p and θ 2p It is not a fixed value, but dynamically related to the water quality grade of the reclaimed water being delivered. Different water quality grades of reclaimed water have varying requirements for urgency and reliability of supply to users, leading to different decision-making and response modes in the water system. Specifically, the system extracts the real-time water quality grade Q of the upstream node of path p from a standardized spatiotemporal dataset. Parameters are determined based on a preset mapping relationship. For example, if the water quality grade is Q1 (e.g., COD < 30 mg / L, for high-standard industrial users), supply interruptions or delays will cause significant losses, thus requiring a rapid decision-making response. In this case, k can be set... 2p =4、θ 2p =1.5, forming a response distribution with a prominent peak and steep shape. If the water quality level is Q3 (e.g., COD ≥ 50mg / L, used for municipal greening), its supply has a certain buffer space, and the decision response can be relatively smoother. In this case, k can be set to... 2p =8、θ 2p =2.5, forming a response distribution with a late peak and a flattened shape. Optionally, this parameter can also take into account daily periodic fluctuations, for example, k 2p It can be a function containing a cos(2πt / 24) term to reflect the differences in dispatch response speed at different times of the day. In this way, the transmission kernel model incorporates practical experience in water operations, making it more realistic and accurate.

[0040] Optionally, constructing a dynamic graph network prediction model also includes: based on the hydraulic parameters of the pipeline network in the standardized spatiotemporal dataset, analyzing the flow distribution ratio of a predetermined hydraulic transport path; and weighting and superimposing the path conduction kernels of each hydraulic transport path according to the flow distribution ratio to form a comprehensive conduction kernel representing the pipeline network node pair.

[0041] In this embodiment, multiple valid paths (such as paths p1, p2, ...) between node pairs (s, d) are identified, and a dedicated propagation kernel K is constructed for each path. p1 (τ), K p2 (τ), ... then, they need to be combined into a single, comprehensive conduction kernel K that characterizes the overall conduction properties between (s, d). sd (τ). Therefore, it is necessary to determine how the total flow is distributed across these paths, which can be achieved by solving a set of nonlinear hydraulic balance equations based on the laws of conservation of nodal flow and conservation of loop energy. Specifically, the Newton-Raphson iterative method can be used to solve this set of equations, thereby obtaining the flow rate Q on each path p.p Accordingly, the flow allocation ratio α p (t) can be calculated as Q p (t) / ΣQ k (t). Based on the calculated flow allocation ratio, the path propagation kernel K for each path is... p (τ) is weighted and superimposed to obtain the comprehensive transmission kernel K. sd (t, τ) = Σ p α p (t)•K p (τ). It can be seen that the constructed dynamic graph network model can accurately characterize the spatiotemporal transmission law of water quality in complex pipe networks, considering multi-path effects and decision delays, through a comprehensive transmission kernel with clear physical meaning and adaptively adjustable parameters, providing a high-quality model foundation for subsequent supply and demand forecasting.

[0042] This embodiment addresses the distortion in the spatiotemporal conduction characteristics of water quality caused by the assumption of a single path and a single delay. Specifically, by exploring and identifying multiple hydraulic transport paths between nodes, it breaks through the limitations of the traditional shortest path approach, accurately reflecting the distribution and convergence effects of water flow in complex ring networks. This allows the model to make accurate predictions even under topological changes such as main pipeline closure or congestion, demonstrating strong robustness. A conduction kernel exhibiting a bimodal distribution is constructed for each path. The first peak accurately represents the physical transmission delay based on fluid dynamics principles; the second peak quantifies the previously neglected decision response delay, namely, human and system lag in operational scheduling. This ensures that the prediction results highly match the actual operational scenarios. Furthermore, the model dynamically correlates the parameters of the decision response delay peak with the quality level of the transported reclaimed water, providing a faster response mode for high-quality, high-requirement water supply tasks and embedding operational experience into the mathematical model. By weighted superposition of the conduction kernels of each path based on the flow distribution ratio calculated according to hydraulic parameters, a comprehensive conduction kernel with clear physical meaning and adaptively adjustable parameters is formed. This achievement represents a leap from single-point-of-time delay prediction to full-time-delay probability distribution prediction, improving the prediction accuracy of water quality arrival time and concentration distribution, and providing a solid foundation for achieving timely matching of supply and demand and avoiding resource waste.

[0043] As an alternative, such as Figure 4 As shown, constructing a dynamic graph network prediction model also includes determining the time-delay propagation kernel, specifically:

[0044] Based on the hydraulic and water quality parameters in the standardized spatiotemporal dataset, a one-dimensional convection-reaction-diffusion equation is established for the pipeline node pairs; the one-dimensional convection-reaction-diffusion equation is solved, and the analytical or numerical solution of the one-dimensional convection-reaction-diffusion equation is used as the time-delayed conduction kernel.

[0045] Specifically, for each segment of the pipeline network, the concentration C of pollutants (such as COD) can be considered to change along two dimensions: pipeline length x and time t, following a one-dimensional Convection-Reaction-Diffusion (CRD) equation: ΨC / Ψt + v•ΨC / Ψx = D eff •Ψ 2 C / Ψx 2 - k•C; where Ψ is the partial derivative, C(x,t) is the pollutant concentration, v is the average velocity of the water flow in the pipe (convective term); D eff is the effective diffusion / dispersion coefficient (diffusion term); k is the first-order reactive degradation coefficient (reaction term). When considering a pulse input generated from source x0 at time t0 (i.e., the Dirac delta function), the Green's function solution K(x, t|x0, t0) of this equation physically represents the concentration distribution formed at position x after a unit concentration pulse at the source, after time t-t0. Therefore, this Green's function solution can be directly used as a time-delayed conduction kernel connecting the source and target nodes. Its analytical solution form is usually a Gaussian distribution or a modified skewed distribution, for example: K(x, t|x0, t0) = (1 / sqrt(4πD) eff •τ))•exp[-(x-x0-v•τ) 2 / (4D eff •τ)]•exp[-k•τ], where τ= t-t0.

[0046] Furthermore, the one-dimensional convection-reaction-diffusion equation is parameterized. To ensure that the aforementioned CRD equation accurately reflects the complexities of real-world pipe networks, its core parameters k and v, as well as the effective pipe diameter involved in the equation, require refined dynamic calibration. Specifically, parameterization is performed based on a standardized spatiotemporal dataset. This parameterization includes:

[0047] The reaction degradation coefficient in the one-dimensional convection-reaction-diffusion equation is determined as a function related to pollutant concentration and water quality level.

[0048] Traditional models typically treat the reaction degradation coefficient k as a constant, but this ignores the complexity of biochemical reactions. In this embodiment, the reaction degradation coefficient k is modeled as a dynamic variable k(C, Q). For example, it can be in the form of a combination of the modified Monod equation or the Arrhenius equation: k(C, Q) = k0•[1 + β1•(C / C)] ref ) n ]•[1 +β2•exp(-E a / (R•T))];where k0 is the baseline degradation rate; C is the real-time pollutant concentration; C refThis is the reference concentration; n is the reaction order; β1 is the concentration effect coefficient; T is the water temperature; E a β0 is the activation energy of the reaction; R is the gas constant; β2 is the temperature effect coefficient. Optionally, k0 can also be a piecewise function related to the water quality level Q. This modeling approach allows the reaction rate to adaptively adjust to changes in pollutant concentration, water temperature, and overall water quality.

[0049] A dynamic growth model of biofilm on the pipe wall was established to dynamically correct the effective pipe diameter and hydraulic roughness used in the one-dimensional convection-reaction-diffusion equation.

[0050] Specifically, with long-term operation, a biofilm will grow and adhere to the inner wall of the pipe. This reduces the effective cross-sectional area of ​​the pipe and increases the surface roughness, thereby affecting the flow velocity v and head loss, which is generally ignored by traditional models. In this embodiment, a differential equation describing the evolution of biofilm thickness Δ over time can be established: dΔ / dt = μ max •(S / (K s +S))•Δ - b•Δ - τ w •Δ 2 / ρ b ; where μ max It is the maximum growth rate; S is the substrate concentration (e.g., COD); K s τ is the half-saturation constant; b is the decay rate; τ w It is the shear force of the water flow on the pipe wall; ρ b This is the biofilm density. By numerically solving this equation, the biofilm thickness Δ(t) at any time t can be obtained. Correspondingly, the effective pipe diameter D... eff (t) = D pipe - 2Δ(t), where D pipe The original inner diameter of the pipe is given; the hydraulic roughness can also be dynamically corrected based on Δ(t). These corrected parameters will be used to solve the CRD equations more accurately.

[0051] Differential propagation speeds are set based on the one-dimensional convection-reaction-diffusion equation according to water quality grades, wherein the propagation speed of high-quality water in the central region of the pipe is higher than that of low-quality water in the near-wall region.

[0052] In this embodiment, the propagation velocity mainly refers to the convection velocity v in the CRD equation. Different qualities of reclaimed water exhibit different hydrodynamic characteristics. For example, Q1 grade high-quality water, which is close to clean water, has a flow velocity v. Q1 The value is basically consistent with the theoretical hydraulic calculation and can be set as v. theoretical • 0.98, where v theoretical This represents the theoretical hydraulic velocity. Q2 grade water is of medium quality, containing a large number of colloidal particles, exhibiting non-Newtonian fluid characteristics. Its apparent viscosity increases, leading to a decrease in actual flow velocity, v.Q2 It can be set to v theoretical • 0.92. Q3 grade low-quality water may contain suspended solids, exhibiting Bingham plastic flow characteristics and yield stress, leading to a further decrease in flow velocity. Q3 It can be set to v theoretical • 0.85. By setting different effective propagation velocities for different water quality levels, the CRD model can more accurately predict the arrival time of water quality fluctuations. It can be seen that this embodiment provides a more refined method for constructing conduction kernels based on physical mechanisms, which is particularly suitable for scenarios requiring high simulation accuracy of water degradation processes.

[0053] This embodiment addresses the issues of roughness and inaccuracy in simulating water quality evolution processes, improving model fidelity. Optionally, a one-dimensional convection-reaction-diffusion (CRD) equation is used as the conduction kernel. Compared to empirical or statistical models, the CRD equation can more profoundly reveal the intrinsic mechanism of pollutant migration and transformation in pipelines. Its core advantage lies in the deep dynamics and domain-specific customization of model parameters: the reaction degradation coefficient k is modeled as a function related to pollutant concentration and water quality level, abandoning the simplistic approach of treating it as a constant in traditional models, making the simulation of reaction rates closer to real biochemical processes. A dynamic growth model of pipe wall biofilm is established, which can calculate biofilm thickness in real time and dynamically correct the effective pipe diameter and hydraulic roughness accordingly. This solves the major deficiency of traditional models in failing to reflect the long-term effects of pipeline aging on hydraulics and water quality. Differentiated propagation speeds are set according to water quality levels, capturing the fast and slow lane effects caused by differences in fluid characteristics of reclaimed water of different qualities. This enables the model not only to predict when it will arrive, but also to accurately predict what state it will be in at that time, providing high-precision simulation capabilities to meet the stringent water quality requirements of high-standard users and ensure the quality and safety of cascade utilization.

[0054] According to one aspect of this application, the method for determining the parameters and the influencing mechanism of the dynamic growth model of the biofilm on the tube wall can be further described as follows:

[0055] Key parameters of the dynamic growth model of biofilm in the tube wall, such as the maximum growth rate μ. max Half-saturation constant K s These parameters are not set arbitrarily, but can be calibrated through laboratory coupon tests or online monitoring devices. Specifically, in a test loop simulating a real pipeline environment, coupons of standard material are placed, and different concentrations (corresponding to matrix concentration S) and flow rates (corresponding to wall shear force τ) are introduced. w The reclaimed water was used. By periodically removing the substrate and measuring its biofilm thickness Δ and dry weight, a set of values ​​for Δ over time t, S, and τ could be obtained. wThe experimental data are varying. By fitting these experimental data to a differential equation model using the nonlinear least squares method, μ can be identified. max The values ​​of key parameters, such as biofilm thickness Δ, adsorption-release rate, and pollutant partition coefficient K, are also considered. Furthermore, the impact of biofilm on conduction, besides correcting for effective pipe diameter and roughness, can be further investigated by introducing an additional conduction delay factor caused by the biofilm adsorption-release effect. Because biofilms act like sponges, they temporarily adsorb pollutants from the water flow and then slowly release them later, which macroscopically manifests as an additional delay and a longer tailing effect. This delay factor can be modeled as a function of biofilm thickness Δ, adsorption-release rate, and pollutant partition coefficient K. d The relevant functions are then convolved with the original conduction kernel to more realistically simulate the conduction process.

[0056] According to another aspect of this application, setting differentiated propagation velocities based on water quality grades allows for a further physical mechanism that enables stratified water flow across the pipe cross-section. Specifically, since the flow velocity is highest in the turbulent core region and lowest near the wall region, high-quality reclaimed water (few impurities, low density) tends to concentrate in the central region of the pipe, while low-quality water (more impurities, high density) is more distributed in the near-wall region. This radial concentration distribution C(r) can be approximately modeled as C(r) = C center + (C wall - C center •(r / R) n Where r is the distance from the center of the pipe, R is the pipe radius, and C center For the concentration at the center of the pipeline, C wall Let be the concentration at the pipe wall surface, and n be the distribution index (radial concentration change rate). Therefore, the overall equivalent propagation velocity of a high-quality water mass is higher than that of a low-quality water mass, creating a fast-slow lane effect.

[0057] According to one aspect of this application, it also includes: sparsifying the dynamic graph network prediction model to generate a sparse connection matrix, wherein the generation of sub-quality demand prediction results and cross-plant collaborative production schemes is performed using the sparse connection matrix.

[0058] In this embodiment, sparsification can address the curse of dimensionality. A complete urban reclaimed water network model may contain tens of thousands of nodes. Considering all possible connections between nodes would result in enormous computational complexity, making it difficult to meet the demands of real-time prediction. Therefore, an intelligent sparsification strategy is needed to identify and retain only the critical connections that play a decisive role in hydraulic conduction and water quality evolution, while ignoring secondary connections, thereby generating a sparse connection matrix. Most elements in this matrix are zero, with only the elements corresponding to critical connections having non-zero values. All subsequent model calculations will be based on this sparse matrix, thus reducing the computational complexity from O(N) to O(N) dimensionality. 2The O(N) level is reduced to a level close to O(N), where N is the number of nodes.

[0059] In exemplary embodiments, two preferred sparsification methods that can be used alone or in combination are provided.

[0060] The first method is an adaptive sparsification method based on system state. The sparsification process includes: quantifying and evaluating the real-time system complexity index C(t) of the pipeline network based on the topological change rate, demand volatility, and water quality heterogeneity in the standardized spatiotemporal dataset; and based on the real-time system complexity index, using ρ(t) = ρ base The functional relationship of ×exp(-β×C(t)) determines the connection retention ratio ρ(t), generating a sparse connection matrix, where ρ base β is the base retention ratio, and β is the sensitivity parameter. In other words, the connection retention ratio used when generating the sparse connection matrix is ​​dynamically adjusted based on the real-time system complexity index.

[0061] Specifically, the sparsity of the model should be dynamically adjusted based on the real-time busy and stable states of the pipeline network system. To achieve this, a system complexity index C(t) needs to be defined to quantify the real-time state of the pipeline network. This index can be a weighted sum: C(t) = λ1•C topo (t) + λ2•C demand (t) + λ3•C quality (t); where t represents the current time; C topo (t) represents the topology change rate, used to measure the drastic nature of changes in the pipeline network structure. For example, it can be obtained by calculating the ratio of the number of newly added or removed pipe segments to the total number of pipe segments at the current moment; C demand (t) represents demand volatility, used to measure the drastic changes in user demand. For example, it can be calculated as the coefficient of variation (standard deviation divided by the mean) of the time series of total user demand; C quality (t) represents water quality heterogeneity, used to measure the complexity of water quality distribution across the entire network. For example, it can be obtained by calculating the Shannon entropy of the volume percentage of different water quality grades (Q1, Q2, Q3) in the pipeline network. λ1, λ2, and λ3 are the weighting coefficients of each item, which can be set empirically, for example, λ1=0.4, λ2=0.3, and λ3=0.3.

[0062] After obtaining the real-time system complexity metric C(t), the system dynamically adjusts the proportion of connections retained, i.e., sparsity. Specifically, the dynamic connection retention ratio ρ(t) can be determined using an exponential decay function: ρ(t) = ρ base ×exp(-β×C(t)); where ρ baseβ is the base retention ratio, which can be set to 0.35 (i.e., retaining 35% of the connections by default); β is a sensitivity parameter used to control the response speed of sparsity to changes in complexity, which can be set to 0.5. It can be seen that when the pipeline system is in a stable state (e.g., at night, with no construction and stable demand), the value of C(t) is small, and ρ(t) is close to ρ... base The system will use a higher sparsity (i.e., retain fewer connections) to save computational resources. However, when the pipeline system is in a state of drastic change (e.g., a surge in water demand due to heavy rain, accompanied by changes in water quality due to pipeline scouring), the value of C(t) will increase significantly, and ρ(t) will decrease exponentially. The system will automatically reduce the sparsity and retain more connections to ensure the prediction accuracy of the model at such critical moments.

[0063] The second method is a protective sparsification method based on domain knowledge, aimed at ensuring water quality safety. The sparsification process includes: quantifying and assessing the risk of water quality incompatibility when water quality mixes at the two ends of a dynamic graph network prediction model based on water quality data in a standardized spatiotemporal dataset; cutting off or weakening network connections where the risk of water quality incompatibility exceeds a preset mixing risk threshold, and generating a sparse connection matrix.

[0064] In this embodiment, the study found that while some connections in the pipeline network are hydraulically feasible, they must be avoided from a water quality safety perspective. Indiscriminate thinning might retain these risky connections. Specifically, the system quantifies and assesses the potential water quality incompatibility risk of each connection (pipeline segment). This risk is multi-dimensional, including, for example, chemical incompatibility risk: adverse chemical reactions may occur when water from different sources is mixed. For instance, if a pipeline is planned to transport water with a residual chlorine content greater than 1.0 mg / L, while the current ammonia nitrogen content in the destination node area is greater than 2.0 mg / L, mixing the two will generate chloramine, which has poor disinfection effects and may produce an odor. The system can calculate a reaction potential Φ based on a chemical reaction kinetic model. react When the threshold is exceeded, the connection is considered to have a high risk of chemical incompatibility. Quality degradation risk: High-quality water flowing into a pipe contaminated with low-quality water will cause irreversible quality degradation. The system can define a degradation risk function R. degrade Its value is related to factors such as the quality grade difference and the importance of downstream users. Risk of microbial cross-contamination: Especially when supplying water to sensitive users such as hospitals and food factories, it is necessary to strictly avoid backflow or mixing of water from areas with low hygiene standards. A microbial cross-contamination risk index R can be established. microbialTo assess this risk, the system quantifies the incompatibility risks of each connection and calculates a comprehensive risk score. If the score exceeds a preset mixed risk threshold, the connection is either severed (its corresponding element in the sparse matrix is ​​set to 0) or weakened (the connection is retained, but its propagation kernel or weight is multiplied by a penalty factor much less than 1). Optionally, a tiered protection mechanism can be designed. For example, for pipelines supplying first-level protected users (such as hospitals), all upstream connections undergo the most stringent risk assessment; while for pipelines supplying third-level users (such as green spaces), a more lenient threshold can be used. This not only improves computational efficiency but also embeds a safety valve at the model structure level, ensuring that any scheduling scheme generated will not cause water quality safety incidents, which is crucial for quality-sensitive resources like reclaimed water.

[0065] This embodiment successfully resolves the sharp contradiction between computational efficiency and operational safety in large-scale pipeline network models. Specifically, to overcome the curse of dimensionality inherent in fully connected graph models, the model is sparsified, reducing the computational complexity from O(N^2) to O(N^2). 2 The O(N) level is reduced to near O(N) level, making real-time, rolling prediction and optimization of city-level pipe networks possible in engineering. It also provides a dual adaptive mechanism: First, dynamic sparsification based on the pipe network system complexity index C(t). This index integrates topology changes, demand fluctuations, and water quality heterogeneity, enabling real-time quantification of the pipe network's workload. Accordingly, the model can use higher sparsity during stable system conditions (such as nighttime) to save computing power, while automatically retaining more connections to ensure accuracy during periods of drastic system change (such as heavy rain), achieving optimal allocation of computing resources. Second, protective sparsification based on water quality incompatibility risk. This mechanism embeds expert knowledge in the field of water safety into the algorithm, quantitatively assessing the chemical reaction risk or quality degradation risk of mixing different water qualities (such as mixing high residual chlorine water and high ammonia nitrogen water), and proactively cutting off or weakening network connections exceeding the risk threshold. This is equivalent to having a safety valve built into the model structure, which prevents optimization algorithms from generating scheduling schemes with water quality safety risks in pursuit of economic efficiency, thus ensuring the safety of reclaimed water, a sensitive resource, in complex collaborative scheduling.

[0066] According to one aspect of this application, the quantitative assessment of water quality incompatibility risk may further include constructing a specific risk assessment matrix and thresholds. Specifically, a chemical incompatibility matrix is ​​established to assess the chemical reaction risk of different combinations of water quality indicators. For example, regarding the problem of chloramine production from mixing high residual chlorine water and high ammonia nitrogen water, the system calculates a reaction potential Φ. react = k r •[Cl]•[NH3]•exp(-E a / RT), where k rLet be the reaction rate constant, and [Cl] and [NH3] be the concentrations of residual chlorine and ammonia nitrogen, respectively. The system has a preset critical reaction potential threshold Φ. critical When the calculated Φ react When the threshold is exceeded, the mixing event is marked as high-risk. For example, a more general water quality mixing compatibility matrix can be provided as follows to guide sparsification decisions: When the mixing type is Q1 water with Q1 water, the risk score (0-1) is 0.05, and the decision recommendation is to allow mixing; when the mixing type is Q1 water with Q1 water, the risk score (0-1) is 0.05, and the decision recommendation is to mix with caution, requiring assessment of downstream user sensitivity; when the mixing type is Q1 water with Q3 water, the risk score (0-1) is 0.95, and the decision recommendation is to prohibit mixing in principle, triggering protective sparsification; when the mixing type is Q2 water with Q3 water, the risk score (0-1) is 0.60, and the decision recommendation is not to mix, only considering it when there are no alternatives. During sparsification, the system compares the calculated or lookup-based risk score with a preset mixing risk threshold (e.g., 0.7). Any connection whose risk score exceeds this threshold will be prioritized for disconnection or weakening.

[0067] In an optional embodiment, the cross-plant collaborative production scheme is built on a water quality cascade utilization network, wherein the effluent from the upper-level user is scheduled as the input water source for the lower-level user.

[0068] In other words, based on a pre-configured water quality cascade utilization network, the effluent from the previous user is scheduled as the input water source for the next user in the cross-plant collaborative production scheme.

[0069] In this embodiment, to maximize the utilization value of reclaimed water, the production plan is not based on the traditional independent water supply model for each user, but rather on a network topology for cascaded water utilization. This is because the effluent from users with high-quality requirements (such as industrial cooling, Level 1), although of slightly lower quality, can still meet the needs of users with lower quality requirements (such as urban miscellaneous use, Level 2). Specifically, the system identifies and constructs the optimal cascaded utilization chain based on each user's water demand quality, effluent quality after use, and their geographical location and pipeline connections, using a maximum weight matching algorithm (such as the Hungarian algorithm). For example, after optimization, a cascaded utilization topology might be formed as follows: Plant A (Q1) -> Industrial Park X (becomes Q2 after use) -> Municipal Green Space Y (becomes Q3 after use) -> Drainage River. Subsequent optimization solutions will be performed under the constraints of this cascaded utilization topology, enabling the scheduling scheme to maximize the value of water resources.

[0070] Furthermore, the generation of cross-plant collaborative production schemes is achieved by solving a large-scale multi-objective optimization problem. This solution process includes: decomposing the large-scale multi-objective optimization problem into a predetermined number of low-dimensional sub-problems along the dimensions of water quality level, time window, and spatial region; coordinating the solutions of the low-dimensional sub-problems through a distributed coordination protocol to obtain the global optimal solution of the large-scale multi-objective optimization problem; wherein, the large-scale multi-objective optimization problem is constructed based on a dynamic graph network prediction model and real-time demand signals.

[0071] Specifically, directly optimizing all variables (output and distribution paths of each water plant, quality level, and time period) of the entire urban pipe network over the next 24 hours would result in a highly complex mixed-integer programming problem (e.g., tens of thousands of dimensions), which is extremely difficult to solve directly. To address this issue, this embodiment employs a block-based recursive decomposition strategy: First-level decomposition – by water quality level: The original problem is decomposed into three relatively independent sub-problems: Q*1, Q*2, and Q*3. These three sub-problems are coupled only at water plants with switchable capacity. Second-level decomposition – by time window: Each water quality sub-problem is further divided into multiple shorter windows along the time axis, such as six 4-hour windows. Detailed optimization is performed within each window, and windows are coupled through boundary conditions (e.g., ending inventory equals beginning inventory of the next period). Third-level decomposition – by spatial region: Based on the community structure in the sparse connectivity matrix, a graph partitioning algorithm (e.g., METIS) is used to divide the pipe network within each time window into several spatial clusters. Clusters are tightly connected and undergo collaborative optimization; inter-cluster connections are sparse, with limited exchange of boundary information. Through this three-level decomposition, the original high-dimensional problem is broken down into hundreds of minimal subproblems with lower dimensions (e.g., approximately 120 dimensions). These subproblems can be solved in parallel across multiple computational cores (e.g., using interior-point methods). After all subproblems are solved, the system coordinates the solutions to each subproblem using a distributed coordination protocol (e.g., a consensus-based distributed optimization protocol or the Alternating Direction Multiplier Method, ADMM). Specifically, each subproblem iteratively corrects its solution based on the solutions of its neighboring subproblems (i.e., boundary information) and penalizes inconsistencies by introducing Lagrange multipliers. After several iterations, the solutions to all subproblems converge to a globally optimal solution that satisfies the global constraints.

[0072] As a preferred implementation method, the implementation of the distributed coordination protocol includes: encoding scheduling rules and supply and demand information into smart contracts deployed on the blockchain system; and reaching a consensus on the global optimal solution among multiple decision-makers through a consensus algorithm.

[0073] In collaborative scheduling scenarios involving multiple water plants and management departments, trust and data transparency are crucial. To address this issue, this embodiment introduces blockchain technology. Specifically, coordination rules, capacity commitments of each water plant, demand orders from each user, and the final generated scheduling plan can all be encoded into smart contracts deployed on a consortium blockchain. For example, a smart contract could stipulate that when Water Plant A registers its ability to provide 1000 cubic meters of Q1-grade water in the next hour, the system will automatically match the demand according to an optimization algorithm and write the scheduling instruction (such as supplying water to Industrial Park X) into the contract. Once Water Plant A completes the water supply according to the instruction and passes online water quality monitoring verification, the contract will automatically execute settlement and record this performance. To ensure that all participants (water plants, scheduling centers, etc.) reach a consensus on the scheduling plan and that it is tamper-proof, the system can adopt a consensus algorithm optimized for water network characteristics, such as the latency-aware Byzantine Fault Tolerance (PBFT-W) algorithm. This algorithm takes into account the physical latency of water transmission in the pipeline network when reaching consensus, making the consensus result executable in reality.

[0074] Optionally, to incentivize water plants to provide stable and high-quality reclaimed water, a Proof of Quality (PoQ) mechanism can be established. Under this mechanism, water plants with better historical records of water quality compliance and stability will receive higher weight or priority scheduling rights during the consensus process, thus creating a virtuous cycle. This not only efficiently achieves distributed collaboration but also enhances the transparency, traceability, and mutual trust among participants in the entire collaborative scheduling system.

[0075] This embodiment solves the problems of limited decision-making dimensions and low efficiency in multi-source and multi-quality collaborative scheduling. Specifically, a water quality cascade utilization network is constructed, and the effluent after use by high-quality users is used as the input water source for the next-level users. Through algorithms such as dynamic programming, the optimal cascade utilization chain is found, thus establishing the goal of maximizing resource value in the top-level design of the optimization model and changing the inefficient mode of independent water supply for each user. To solve the resulting ultra-large-scale multi-objective optimization problem, a decomposition strategy along three dimensions of water quality level, time window, and spatial region is constructed. This strategy can decompose a global problem with tens of thousands of dimensions that is difficult to handle into hundreds of low-dimensional sub-problems that can be solved in parallel. After solving the sub-problems, through a distributed coordination protocol based on consistency (such as ADMM) for overall coordination, the final convergence to the global optimal solution is achieved. It successfully transforms the non-deterministic polynomial hard (NP-hard) problem that is theoretically difficult to solve into a computationally feasible task in engineering. Optionally, by introducing a blockchain intelligent contract execution coordination protocol, the transparency and mutual trust of multi-agent (such as different water plants) collaboration are further enhanced. This makes a truly global-optimal cross-plant collaborative production plan possible, which can reduce the total system operation cost, avoid the use of high-quality water for low-quality purposes, and comprehensively improve the overall economic and environmental benefits of the urban reclaimed water system.

[0076] According to one aspect of the present application, the process of constructing a water quality cascade utilization network may further include:

[0077] Using the dynamic programming algorithm to find the longest cascade utilization chain that satisfies the water quality decreasing constraint.

[0078] Specifically, the system will construct a dynamic programming model with pipe network nodes as states and water quality levels as stages. Define the state value function V(i, q) as the maximum economic or environmental benefit of the subsequent cascade utilization chain that can be formed when starting from node i and using water with a quality level not higher than q as the input. Its state transition equation can be expressed as: V(i, q) = max{V(j, q') + R(i, j, q, q')}, where j is the downstream node of i, q' is the water quality level converted after utilization at node i and q' < q, and R(i, j, q, q') is the benefit generated by this single-step cascade utilization from node i to j. By solving this dynamic programming problem backward, the globally optimal longest cascade utilization chain can be obtained.

[0079] Construct a cascade matching matrix, and the matrix elements of the cascade matching matrix comprehensively consider the demand matching degree, quality matching degree, and distance factor among users.

[0080] Or, construct a cascade matching matrix M for evaluating and screening the feasibility of cascade connections cascade . The elements M of this matrix cascade(i, j) represents the degree of matching between the water output of the previous user i and the input water source of the next user j. Its value can be comprehensively considered from the following three dimensions: M cascade (i, j) = Demand match (i, j) • Quality match (i, j) •Distance factor (i, j); where Demand is the degree of matching. match The ratio of user i's water output to user j's water demand; Quality matching degree. match It is an exponential function, and its value decreases as the difference between the effluent water quality of user i and the demand water quality of user j increases; the distance factor is... factor This is a decay function related to the network distance between the two users. Only when M... cascade The ladder connection is considered feasible only when (i, j) exceeds a preset threshold (e.g., 0.6).

[0081] Establish the quality conversion function ΔQ during the usage process. use This function is used to quantify water quality degradation and potential energy recovery at various levels of utilization. It is specific to the water-using process. For example, for industrial cooling (non-contact), the water quality degradation model is ΔQ. cooling = 0.05•Q in •(1 + 0.1•ΔT), mainly considering the effect of temperature rise ΔT, where Q in Input water quantity or water quality level. For agricultural irrigation, the water quality degradation model is ΔQ. irrigation = 0.3•Q in •(1 + 0.2•ET / P), which additionally considers the concentration effect of evapotranspiration ET on pollutants, where P is precipitation. Simultaneously, the system will establish an energy recovery model, for example, utilizing the head difference between stages for micro-turbine power generation, or utilizing temperature differences to recover heat energy through heat pumps. The recovered energy can offset approximately 15-20% of the operating costs.

[0082] According to one aspect of this application, the process of solving a large-scale multi-objective optimization problem through block recursive decomposition can further include the following detailed algorithm implementation: A three-level decomposition strategy is executed to decompose the original high-dimensional problem into multiple low-dimensional, parallelizable subproblems. The first-level decomposition (by water quality level): The problem is decomposed into three subproblems, Q*1, Q*2, and Q*3, which are coupled only at nodes with switchable capacity. The second-level decomposition (by time window): The 24-hour scheduling cycle is divided into six 4-hour windows, with windows coupled through boundary variables (such as ending inventory). The third-level decomposition (by spatial region): Based on the community structure of a sparse connectivity matrix, graph partitioning algorithms such as METIS are used to divide the pipeline network within each time window into approximately 10 spatial clusters, with tight coupling within clusters and loose connections between clusters. Based on this, the numerous decomposed subproblems are solved in parallel and coordinated in a distributed manner. Specifically, the system can distribute hundreds of minimal subproblems to multiple computing cores, with each core solving them independently using efficient algorithms such as the interior-point method. After all subproblems have been solved in one round, a distributed coordination protocol based on the Consistent ADMM algorithm is initiated. Its iteration rule can be: Local variable update: for each subproblem p, x is solved. p (k+1) = argmin(f p (x) + ρ / 2•||x - z p k || 2 ), where x p (k+1) f is a local variable of subproblem p (in the (k+1)th round). p (x) is the objective function of subproblem p, ρ is the penalty strength of the consistency constraint, x is the variable used for optimization in the current iteration, and z p k This refers to the globally consistent variable value referenced by subproblem p in the k-th iteration. Global variable update: Information is exchanged with neighboring subproblems through boundary nodes, and the average value z is taken. (k+1) = average(x q (k+1) ), where z (k+1) x is the globally consistent solution obtained after the (k+1)th iteration of all subproblems through boundary node information exchange and averaging operations. q (k+1) This represents the local solution to the neighboring subproblem q that is coupled with subproblem p. Dual variable update: λ p (k+1) =λ p k +ρ(x p (k+1) - z (k+1) ), where λ p (k+1)Let be the dual variable of subproblem p in the (k+1)th iteration. This iterative process continues until solutions to all subproblems are consistent, converging to the global optimum. To guarantee the convergence of this iterative algorithm, it can be theoretically proven that its update operator is a compression mapping. Its convergence is based on Banach's fixed-point theorem, and the convergence speed is related to the spectral radius of the weight matrix of the communication topology. In practical applications, strategies such as warm-start are also employed, using solutions from similar historical scenarios as initial points to accelerate convergence.

[0083] In optional embodiments, the method further includes: generating a water quality evolution risk map containing potential future risks by applying chaos theory analysis or infectious disease dynamics models based on water quality time series data in a dynamic graph network prediction model.

[0084] In this embodiment, the water quality evolution risk map is a data structure that goes beyond traditional deterministic or simple probabilistic predictions. It aims to reveal and quantify low-probability, high-impact water quality mutation events that may occur in the pipeline network in the future, and their potential propagation paths. Conventional prediction models excel at fitting periodic and trend changes, but their predictive ability for black swan events such as sudden pollution incidents is limited.

[0085] As a preferred implementation, this map can be generated by fusing the following two analysis methods:

[0086] The first analysis method is water quality mutation prediction based on chaos theory. This method is used to identify precursors to water quality mutations at individual key nodes. Specifically, the system can extract time-series data of water quality parameters (such as COD) over a past period (e.g., 7 days) from a dynamic graph network model for a key monitoring node (e.g., a node upstream of an important user). The Takens embedding theorem is used to reconstruct the phase space of this one-dimensional time series, transforming it into a trajectory in a high-dimensional space. The optimal embedding dimension m and delay time τ can be determined using the pseudo-nearest neighbor method and mutual information method. The maximum Lyapunov exponent λ of this trajectory is then calculated. max If λ max A value greater than 0 indicates that the hydrodynamic system exhibits chaotic characteristics, meaning that small initial disturbances are amplified exponentially, suggesting inherent unpredictability and the possibility of sudden changes in the system. The upper limit of the prediction timescale is approximately 1 / λ. max Critical slowing occurs when a system approaches a critical point (i.e., an impending abrupt change). This phenomenon can be monitored by continuously calculating the autocorrelation coefficient and variance of the time series within a sliding time window. When a significant increase in the autocorrelation coefficient and a sharp increase in the variance are detected, it can be determined that an abrupt change is imminent, thus issuing an early warning; for example, the occurrence of an abrupt change can be predicted 2-4 hours in advance.

[0087] The second analytical method is water quality deterioration propagation modeling based on infectious disease dynamics. This method predicts the potential spread and speed of water quality deterioration throughout the pipe network should a node experience a deterioration (such as the abrupt change predicted by the methods mentioned above). Specifically, the water quality deterioration process can be analogized to the spread of an infectious disease in a population. The system can classify nodes in the pipe network into three categories: S-class (Susceptible nodes), i.e., nodes with normal water quality; I-class (Infected nodes), i.e., nodes with deteriorated water quality; and R-class (Recovered nodes), i.e., nodes whose water quality has improved through treatment. Based on this, an improved SIR (Susceptible-Infected-Recovered) model is established. Unlike the standard SIR model, the infectivity β here is not a constant, but a variable related to the pipe network topology. For example, β... ij The transmission rate from node i to node j can be related to the degree of nodes i and j (i.e., the number of connected pipes) and the hydraulic transmission efficiency between them, implying that hydraulic hub nodes have stronger transmission capabilities. Similarly, the recovery rate α can be related to the treatment capacity of downstream water plants. Using this model, the basic regeneration number R0 of the pipe network can be calculated. If R0 > 1, it means that once a single contamination point appears, water quality deterioration will spread widely throughout the network; if R0 < 1, the deterioration will be confined to a local area and eventually subside. By simulating this model, the probability P of each node being infected (i.e., water quality deterioration) in the future can be obtained. infect (i, t) and average infection time T infect (i).

[0088] In summary, by combining the results of the two methods mentioned above, the system can generate a water quality evolution risk map. This map clearly indicates which nodes are at risk of sudden changes in the future, the most likely propagation path once a change occurs, and the degree and duration of its impact on downstream nodes.

[0089] The water quality evolution risk map is transformed into risk constraints. Under the premise of meeting the risk constraints, sub-quality demand prediction results and cross-plant collaborative production plans are generated.

[0090] In this embodiment, the risk map is not only used for early warning, but is also directly integrated into the multi-objective optimization model as a set of hard or soft risk constraints. Specifically, this transformation can take the following forms: Probability constraint (chance constraint): For high-risk node j, a constraint is added to the optimization model, namely P(Q j (t) ≥ Q min ) ≥ 1 -α. Where P(·) represents probability, Q j Q(t) represents the water quality at node j at time t. minThe minimum water quality standard is defined by α, which is a very small risk tolerance level (e.g., 0.05). This constraint means that the optimal solution must ensure that the probability of water quality non-compliance at node j is no higher than 5%. This probabilistic constraint can be transformed into a linear or nonlinear constraint that can be handled by standard optimization models through its deterministic equivalent form. Conditional Value at Risk (CVaR) constraint: This is a more robust risk measure. The optimization model can be required to minimize the expected operating costs while constraining the CVaR value of the costs (i.e., the average cost under the 5% worst-case scenario) to not exceed a preset budget limit L. CVaR This forces the optimizer to consider not only the desired outcome but also potential extreme high-cost events when making decisions. Resilience constraint: A comprehensive system performance index (System) can be defined. Performance (t), and add constraints System Performance (t) ≥ P min •(1 - D(t)•(1 - R(t)));where, P min Here, D(t) represents the minimum acceptable performance level, D(t) is the disturbance strength predicted by the risk map, and R(t) is the system's resilience. This constraint ensures that even when the disturbance predicted by the risk map occurs, the system performance degradation will not breach the preset bottom line. By incorporating these constraints derived from the risk map into the optimization model, the model is forced to avoid solutions that, while having the lowest cost under normal circumstances, are vulnerable under risk scenarios. For example, the optimizer might choose to activate a backup water source or select a longer but safer water transmission path. Although this slightly increases normal operating costs, it provides insurance against unknown risks, thus making the final production solution both economical and robust.

[0091] According to one aspect of this application, the process of utilizing a water quality evolution risk map can further include a quantitative assessment of pipeline network resilience. Specifically, after generating the risk map, the system can calculate multi-dimensional resilience indicators, such as resistance R. resist Resilience R is defined as the reciprocal of the decrease in key performance indicators (such as water supply reliability) when the system is subjected to a disturbance. recover Resilience metrics are defined as the reciprocal of the time required for system performance to recover from its lowest point to a normal level. By simulating deliberate attacks on the network (e.g., removing critical nodes or edges in the simulation) and observing changes in these resilience metrics, system vulnerabilities and resilience bottlenecks can be identified, providing quantitative evidence for network hardening and emergency response planning.

[0092] In an optional embodiment, the method further includes a step of online correction of the parameters of the dynamic graph network prediction model, comprising:

[0093] The prediction results of the graph network prediction model are compared with the actual operating data to obtain the prediction error.

[0094] In this embodiment, the starting point for model correction is obtaining feedback signals, i.e., the gap between prediction and reality. After a prediction and scheduling cycle (e.g., 72 hours) ends, the system automatically collects the actual operating data of the pipeline network during that cycle. This data includes: the actual water consumption of users in various fields, the actual water quality test results uploaded by online monitoring equipment deployed at each node of the pipeline network, and user satisfaction evaluations, etc. The system compares these actual data with the prediction results previously generated by the model point by point and hour by hour, and calculates the prediction error sequence in different dimensions (e.g., water demand prediction error, water quality concentration prediction error, and conduction delay prediction error). Furthermore, in order to diagnose the root cause of the error, the system can use statistical methods such as variance decomposition to decompose the total error into three main components: systematic bias, which is usually caused by defects in the model structure itself; random error, which is caused by unpredictable random disturbances in the system; and conduction error, which is caused by the model's estimation bias of the time-delay conduction kernel. The correction mechanism in this embodiment will mainly focus on correcting the model parameters related to conduction error.

[0095] The parameters of the dynamic graph network prediction model are continuously adjusted based on the prediction error using the recursive least squares algorithm.

[0096] In this embodiment, the model parameters, particularly the conduction kernel parameters, are updated smoothly and continuously online. The Recursive Least Squares (RLS) algorithm updates the parameter estimates with each new data point without reprocessing all historical data. Specifically, the parameters of the conduction kernel (e.g., k1, θ1, k2, θ2, etc.) can be organized into a parameter vector θ to be estimated. The observation equation can be expressed as y(t) = φ(t). T θ + e(t); where y(t) is the new observation at time t (e.g., the true propagation delay obtained through tracer experiments or data correlation analysis); φ(t) is the feature vector corresponding to this observation (e.g., containing information such as path length, flow rate, and water quality level); and e(t) is the observation noise. T This is the transpose. The parameter update formula for the RLS algorithm is: θ(t+1) = θ(t) + K(t)[y(t) -φ(t)] T [θ(t)]; where θ(t+1) and θ(t) are the parameter vector estimates before and after the update, respectively. K(t) is a key adaptive gain matrix that determines the step size and direction of this update. Furthermore, the calculation of K(t) introduces a variable forgetting factor λ(t). The formula for calculating λ(t) is: λ(t) = λ min + (λ max - λ min)•exp(-|e(t)| / e ref ); where e(t) is the current prediction error, λ min and λ max These are the upper and lower bounds of the forgetting factor (e.g., 0.95 and 0.995), e ref It is the reference error. When the prediction error e(t) is large, λ(t) will become smaller, allowing the system to forget past data more quickly and adapt to new changes in the system rapidly; when the prediction error is small, λ(t) will become larger, making more emphasis on utilizing historical information, resulting in more stable parameter estimation.

[0097] Simultaneously, an accumulation algorithm is used to monitor for parameter mutations, triggering accelerated corrections when mutations are detected.

[0098] In this embodiment, the RLS algorithm excels at tracking gradual, slow changes in parameters (e.g., the slow increase in roughness of a pipe due to annual scaling). However, for abrupt changes in the system (e.g., a section of old pipe being completely replaced with a new material, causing a step change in its conductivity), the RLS response may not be fast enough. To address this issue, a cumulative sum (CUSUM) algorithm is deployed in parallel, specifically designed to detect abrupt changes in parameters. The CUSUM algorithm continuously monitors the parameter sequence θ(t) updated by RLS and calculates two cumulative sum statistics S. + (t) and S - (t). These two statistics are used to detect whether the parameter has undergone a significant upward or downward shift, respectively. S + (t) = max(0, S) + (t-1) + θ i (t) - μ0 - k); S - (t) = max(0, S) - (t-1) - θ i (t) + μ0 - k); where θ i (t) is a monitored parameter, μ0 is the baseline mean of that parameter, and k is the tolerance. Once S + (t) or S - If (t) exceeds the preset decision threshold h, the system determines that the parameter has undergone a sudden change. Upon detecting the change, the system immediately triggers an acceleration correction mechanism. This mechanism temporarily adjusts the behavior of the RLS algorithm, for example, by forcing its forgetting factor λ to a very small value (e.g., λ...). fast= 0.9), and reset its error covariance matrix. This allows the RLS algorithm to learn new system characteristics at an extremely fast speed in the following iterations, thereby achieving a rapid response to abrupt events. Through the collaborative work of RLS and CUSUM, a robust and efficient online model correction closed loop that can adapt to slow drift and respond quickly to abrupt changes is constructed, ensuring the accuracy and timeliness of the prediction model in long-term operation.

[0099] In an optional embodiment, demand propagation modeling under time-delay constraints is performed on the real-time demand signal to generate a time-delay demand distribution matrix.

[0100] Research has found that water plant production decisions must be made in advance of users' actual water demand; this advance time is the water delivery delay in the pipeline network. Therefore, simply using current user demand signals as the basis for water plant production is incorrect. A mathematical model can be constructed that translates future user demand into current water plant production tasks. As a preferred implementation, this modeling process includes the following aspects: constructing a reverse transmission kernel, and a forward transmission kernel K from the water plant to the user. forward (τ) Conversely, the reverse conduction nucleus K reverse (τ) aims to project the demand of user node j at time t backwards onto the upstream water plant node i, thereby calculating how much production water plant i should start at time t-τ to meet the demand. Specifically, this backward propagation kernel can be obtained from the known forward propagation kernel using Bayesian inference. Optionally, the modeling here considers the asymmetry of information propagation: the backward propagation speed v of the demand signal (information) info It can travel much faster than the forward propagation speed v of water (matter). flow For example, v can be set info = 10 × v flow This means that the timescale of the reverse kernel will be compressed accordingly when constructing it, but a minimum response time (e.g., 0.5 hours) will be retained to account for the physical constraints of production start-up. Furthermore, when a user j is supplied by multiple water plants (e.g., i and k) simultaneously, its demand needs to be proportionally reverse-allocated. The system can calculate the allocation weight w based on the capacity P of each water plant and the capacity availability R for that user (obtained by integrating the forward transmission kernel), ultimately forming a comprehensive reverse transmission kernel that considers multi-source supply. A propagation and amplification model of demand uncertainty is established. Any prediction of future demand is uncertain, and the longer the prediction time, the greater the uncertainty. The cumulative amplification effect of this uncertainty in the pipeline transmission process can be modeled. Specifically, the system will identify the functional relationship between the standard deviation σ of the prediction error and the prediction time domain τ by analyzing historical prediction error data. This relationship can be modeled as an uncertainty amplification function: σ(τ) = σ0•(1 + α×τ)β The model can be further divided into two parts: σ(τ) and β(τ). σ(τ) represents the standard deviation of the error in predicting future demand at time τ; σ0 represents the standard deviation of the initial prediction error (i.e., the error when τ=0); and α and β are amplification coefficients to be identified, which can be determined by maximum likelihood estimation of historical error data. Optionally, the model can also incorporate an event factor, such as an additional multiplication by a coefficient greater than 1 during special events like holidays or extreme weather, to reflect the sharp increase in demand uncertainty during these periods.

[0101] The system generates a scenario set for robust optimization based on an uncertainty model. After obtaining the expected value of demand (calculated using a backpropagation kernel) and the uncertainty of time-domain growth (calculated using an amplification function), the system does not simply use the expected value for optimization. Instead, it generates a scenario set containing multiple possible demand scenarios for robust optimization. Specifically, the system can construct a probabilistic scenario tree, dividing the next 24 hours into multiple stages, each stage generating several branches representing different possible demand fulfillment scenarios. Simultaneously, to ensure the solution's resilience to extreme events, the system also generates tail event scenarios with a small probability (e.g., 5%), such as a large user needing emergency flushing due to a sudden accident, causing a sudden increase in demand of 3 standard deviations. Through methods such as K-means clustering, hundreds or thousands of original scenarios can be reduced to, for example, 50 representative weighted scenarios. These scenarios together constitute a time-delay demand distribution matrix, containing not only the expected value of demand but also its possible distribution patterns around that expected value. Using this distribution matrix as input to the optimization model allows the generated production plan to meet various possible demand fluctuations with a high probability, thereby improving the reliability of the water supply system.

[0102] According to another aspect of this application, the process of generating the scene set required for robust optimization can be further implemented as follows: After generating a large number (e.g., thousands) of possible future demand scenarios through an uncertainty amplification function, the system employs a K-means clustering algorithm to reduce the complexity of subsequent optimization calculations. Specifically, the K-means algorithm divides (clusters) the high-dimensional scene set into a preset number (e.g., 50) clusters. The goal of the algorithm is to minimize the sum of squared distances from each scene to the centroid of its cluster (the weighted average scene of that cluster). After clustering, the centroid of each cluster is selected as a representative scene, and the probability of this representative scene is set as the sum of the probabilities of all the original scenes within that cluster. The probability distribution characteristics of the original thousands of scenes can be highly approximated using 50 carefully selected and weighted representative scenes, thereby improving the solution efficiency while ensuring robustness.

[0103] In a specific embodiment, it is assumed that in a city's reclaimed water pipeline network, there exists a source node S (a reclaimed water plant) and a target node D (a user in a high-tech industrial park). Node S stably produces Q1 grade reclaimed water (COD concentration of 25 mg / L). Pipeline topology analysis identifies two main hydraulic transport paths from S to D, denoted as path 1 (P1) and path 2 (P2). The provided input parameters are as follows: Path geometry and hydraulic parameters: Path P1: Total length L1 = 5 km. Average flow velocity v1 = 1.2 m / s calculated based on pipe diameter and real-time flow rate. Path P2: Total length L2 = 8 km. Average flow velocity v2 = 1.0 m / s calculated based on pipe diameter and real-time flow rate. Flow state: The flow in both paths is calculated to be turbulent. Source water quality: Q1 grade. The physical transport peak (Γ) is calculated for path P1 (L1=5km, v1=1.2m / s, Q1 grade). 11 Parameters: According to the formula under turbulent conditions: shape parameter k 11 = 2.2 + 0.12 * L1 / v1 = 2.2 + 0.12 *5000 / 1.2 ≈ 2.2 + 500 = 502.2 (The application of this formula may require adjustment of units or coefficients to conform to actual physical scales; this is an example calculation logic). To make the values ​​reasonable, assume the formula is k1 = 2.2 + 0.12*L (L in km), θ1 = 1.1 + 0.06*L. Then: shape parameter k 11 = 2.2 + 0.12 * 5 = 2.8. Scale parameter θ 11 = 1.1 + 0.06 * 5 = 1.4. The mean of this peak is k. 11 * θ 11 = 2.8 * 1.4 = 3.92 hours, representing the average physical time for water to flow through a 5-kilometer pipe. Determine the decision response peak (Γ). 21 Parameters: Due to the need to transport Q1 grade high-quality water, a rapid response is required: shape parameter k 21 = 4.0; scale parameter θ 21 = 1.5. The mean of this peak is k. 21 * θ 21 = 4.0 * 1.5 = 6.0 hours (the decision delay here is greater than the physical delay, which is a reasonable scenario). The propagation kernel K1(τ) forming path P1: K1(τ) = w 11 •Γ(τ;2.8,1.4) + w 21 •Γ(τ;4.0,1.5). The weight w can be determined based on the path length; for example, let w be denoted as w. 11 =0.6, w21 =0.4. Calculate the physical transport peak (Γ) for path P2 (L2=8km, v2=1.0m / s, Q1 level). 12 Parameter: Shape parameter k 12 = 2.2 + 0.12 * 8 = 3.16; Scale parameter θ 12 = 1.1 + 0.06 * 8 = 1.58; the mean of this peak is k. 12 * θ 12 ≈ 4.99 hours. Determine the peak decision response (Γ). 22 Parameters: Since the source water quality levels are the same, the decision-response patterns are also the same: shape parameter k 22 = 4.0; scale parameter θ 22 = 1.5; The conduction kernel K2(τ) that forms path P2: K2(τ) = w 12 •Γ(τ;3.16,1.58) + w 22 •Γ(τ;4.0,1.5). The weights can be set as w. 12 =0.7, w 22 =0.3.

[0104] By solving the hydraulic balance equations, the system calculations show that, due to the shorter path P1 and lower hydraulic resistance, 65% of the flow will choose path P1, and 35% will choose path P2. Therefore, the flow distribution ratios are: α1 = 0.65, α2 = 0.35. The conduction cores of the two paths are weighted and summed according to their flow distribution ratios to obtain the comprehensive conduction core K between the node pair (S, D). sd (τ): K sd (τ) = α1•K1(τ) + α2•K2(τ) K sd (τ) = 0.65 • [0.6•Γ(τ;2.8,1.4) + 0.4•Γ(τ;4.0,1.5)] + 0.35 • [0.7•Γ(τ;3.16,1.58) + 0.3•Γ(τ;4.0,1.5)]. The final calculated integrated conduction kernel K is... sd(τ) is a complex, non-standard probability density function. If plotted graphically, with the horizontal axis representing the time delay τ (in hours) and the vertical axis representing the probability density, the curve will exhibit the following characteristics: Asymmetry: Due to the inherent skewness of the gamma distribution, the entire curve will show a pronounced right-skewed (tailed) shape; Multi-peak / broad-peak shape: As it is the superposition of four gamma distributions, the final curve will no longer be a simple single peak, but is likely to present a broad peak with a distinct shoulder, or even two closely spaced, partially overlapping peaks. The main peak will be closer to the physical transmission delay of P1 (approximately 3.92 hours) because it dominates the flow, while the contribution of P2 and the presence of the decision delay peak will make the entire distribution last longer and have a more complex shape. This K... sd The (τ) curve precisely describes the following physical and decision-making process: After a water quality adjustment operation at water plant S, the monitoring point at industrial park D first detects changes approximately 3-4 hours later, with the most drastic changes occurring around 4-5 hours. However, the complete arrival and stabilization of the entire water quality signal may take up to 8-10 hours. It can be seen that the adaptive multipath propagation kernel constructed in this embodiment provides far richer and more accurate spatiotemporal dynamic information compared to a simple average delay time (e.g., (3.92+4.99) / 2 = 4.46 hours), which is the cornerstone for achieving high-precision supply and demand coordinated forecasting.

[0105] Furthermore, to verify the superiority of adaptive multipath transmission, which considers the unique characteristics of water networks, over traditional single-path transmission, the following comparative experiments can be added: The experimental design can be based on a typical grid-like pipe network topology containing 20 nodes and 35 edges, performing a 72-hour simulation prediction. Under the same supply and demand scenario, the multipath model of this application (searching the first 3 paths) and the single-path model as a control group (considering only the shortest path) were deployed respectively. Experimental results show that the multipath model of this application outperforms the single-path model in several key performance indicators: the mean absolute percentage error (MAPE) of the multipath model is 8.3%, while that of the single-path model is 15.7%, improving accuracy. Especially during peak pipe network load periods, the single-path model has a higher prediction deviation due to ignoring water flow splitting, while the multipath model has a lower deviation. In the scenario simulating the closure of the main pipe, the single-path model incorrectly predicts a 3-hour water supply interruption because its only dependent path fails. The multipath model in this application can accurately predict that the flow will be redistributed through alternative paths, with an interruption time error of less than 10 minutes, demonstrating robust adaptability to network topology changes. Simulation tracing experiments revealed that Q1 high-quality water does indeed preferentially choose the faster flow path in the center of the pipeline, arriving downstream earlier than Q3 water. The single-path model cannot capture this delay difference caused by water quality-specific transmission (a physical characteristic unique to water networks and absent in data packet transmission in communication networks). The multipath model in this application, considering the water quality-specific propagation speed, has an average error of only 12 minutes in predicting water arrival time, far superior to the 45 minutes of the single-path model. The above experimental data quantitatively demonstrate that the multipath transmission method proposed in this application is not a simple transfer of technology from other fields, but rather produces unexpected technical effects by deeply integrating the hydraulic and water quality-specific physical laws of reclaimed water pipeline networks.

[0106] In a further embodiment, a virtual production capacity pool that transcends the boundaries of the physical water plant is constructed, generating a three-dimensional virtual production capacity tensor. Research has found that the rated production capacity of a water plant is not equivalently obtainable at any given time and location. Affected by pipeline transmission delays, losses, and the opportunities for resource reuse within the network, the effective supply capacity of a water plant to different users is dynamically changing. Therefore, it is necessary to construct a virtual production capacity pool to describe all physical and virtual water sources in a unified spatiotemporal dimension.

[0107] As a preferred implementation, the construction process includes the following aspects: performing time-shifted capacity projection calculations. This accurately converts the water plant's output capacity into the available capacity at the user's side. Specifically, the system utilizes the time-delay propagation kernel K from water plant i to user j. ij (τ), representing the real-time production capacity curve P of the water plant. i (t) Perform convolution operation to calculate the effective productivity contribution P of user j at time t. ijeff (t): P ij eff (t) =∫P i (t-τ)×K ij (τ)×η ij (τ)dτ;where P i (t-τ) represents the production volume of water plant i at time t-τ in the past; K ij (τ) is the value of the conduction kernel at time delay τ, representing the probability that water leaving the plant at time t-τ will reach user j at time t; η ij (τ) is the transmission efficiency coefficient, which can be a function related to the time delay τ, used to characterize the evaporation and leakage losses (e.g., a loss of 2-5%) and the decay of water quality indicators such as residual chlorine during long-distance water transmission. Through this calculation, the original production capacity pulse at the water plant side is projected and diffused into a production capacity availability curve with a specific time distribution at the user side; this is the time-shifted production capacity. Then, the floating production capacity is optimized and allocated. This is for flexible production units in the pipeline network that can produce reclaimed water of multiple qualities (e.g., Water Plant C can produce Q1 and Q2 grade water). This production capacity is called floating production capacity. The system needs to optimize its allocation ratio among different qualities. Specifically, the system establishes a switching cost model. For example, switching from producing Q1 to Q2 requires a 2-hour transition time, during which the production capacity drops to 50% of the rated value, and additional chemical and equipment wear costs are incurred. Based on this, the system constructs a two-layer optimization model, with the upper layer deciding the switching time and the lower layer optimizing the output at each time period. Its objective function is: min Σ t [c1•P C Q1 (t) + c2•P C Q2 (t)] +λ•Σ t |P C (t)-P C (t-1)|; where c1 and c2 are the unit costs of producing Q1 and Q2 grade water, respectively; P C Q1 (t) and P C Q2 (t) represents the output during period t; the second term is the switching penalty term consisting of the switching penalty coefficient λ, used to avoid frequent switching due to small cost differences. This model can be solved using methods such as dynamic programming to obtain the optimal switching strategy (e.g., a maximum of 2 switching times per day, switching to high quality before the morning peak and switching to low quality after the evening peak), thereby accurately allocating floating production capacity to different water quality levels.

[0108] Furthermore, as an optional enhancement, the construction of a virtual capacity pool can also integrate a water quality cascade utilization network. The used water within the network is also considered a virtual water source and thus included in the capacity pool. Specifically, the system identifies feasible cascade utilization chains within the network. For example, after industrial user A uses Q1-grade water, its effluent quality drops to Q2, but it can still meet the needs of municipal user B. In this case, the system considers the effluent from industrial user A (whose quantity and quality can be modeled using a quality conversion function) as a virtual capacity available to municipal user B. This virtual capacity also needs time-shift projection, i.e., calculating the transmission kernel from A to B to determine when A's effluent can reach B. In this way, the sources of the virtual capacity pool not only include physical water plants but also hundreds or thousands of virtual water sources (i.e., the outlets of the previous user) in the network, enriching the dispatchable resources and providing a foundation for maximizing the value of reclaimed water resources.

[0109] In summary, the system integrates the time-shifted capacity of all physical water plants, the floating capacity of flexible water plants, and the virtual capacity in the cascade network, ultimately forming a unified three-dimensional virtual capacity tensor with dimensions of [source ID × water quality grade × time]. This tensor provides comprehensive, dynamic, and accurate supply-side input for subsequent global optimization.

[0110] In an exemplary embodiment, the specific calculation process of the switching cost model for floating capacity can be further described as follows: Switching cost C switch It is decomposed into three specific components: C switch = C loss + C chem + C wear C loss The cost of lost production capacity is calculated as C. loss =0.5×P rated ×T trans ×price, P rated This is the water plant's rated capacity, T trans This is the transition time required for water quality switching (e.g., 2 hours), and price is the selling price per unit of reclaimed water. 0.5 indicates that the production capacity will be reduced to 50% during the transition period. chem The cost of water quality conditioner is the additional chemical cost required for each switchover, for example, 500 yuan per switchover. wear The equipment wear and tear cost reflects the additional damage caused to valves, pump sets, and other equipment by frequent switching, for example, 200 yuan per switch. By accurately modeling this cost, the optimization algorithm can make switching decisions that are more economically efficient.

[0111] In one embodiment of this application, a more detailed and multi-layered implementation of sparsity processing is provided. By combining various domain-knowledge-based graph algorithms and quantization metrics, the generated sparse connection matrix can retain key network information to the greatest extent while improving computational efficiency. Specifically, this includes:

[0112] Micro-level sparsification is performed to preserve critical transport paths with high flow influence. The goal is to identify pipe segments that play a locally central role in water transport and distribution from the most basic edge level. Specifically, the system calculates a flow influence index S for each edge (pipe segment) (i, j) in the network. ij Optionally, this indicator can be calculated as follows: Define the instantaneous influence I of the edge. ij (t) = Q ij (t)•N j •W j •exp(-λ•L ij ); where Q ij (t) is the real-time flow rate of this pipe section; N j It is the total number of nodes that its downstream node j can reach (calculated using breadth-first search); W j It is the weighted sum of the importance levels of all downstream users (e.g., industrial users have a weight of 3.0, municipal users have a weight of 2.0); L ij λ is the pipe segment length; λ is the distance attenuation coefficient. To account for time delay effects, the cumulative time influence I is calculated. ij cum =∫I ij (t-τ)•K ij (τ)dτ represents the convolution of the instantaneous influence with the propagation kernel of that edge. Taking into account the volatility of the influence, the comprehensive influence score S is obtained. ij = I ij cum •(1+γ•σ ij ), where σ ij γ is the coefficient of variation of the influence time series, and γ is the volatility weight. The system will process the S values ​​of all edges. ij The scores are sorted, and the top 30% of edges are retained as key edges. For edges that are not selected, their substitutability (i.e., the ratio of the length of the bypass path to the length of the direct path) is further calculated, and edges that are easily substitutable are sparsified first.

[0113] Meso-level sparsification is performed, preserving cross-water plant connections with strong collaborative necessity. The aim is to identify, at the water plant cluster level, which water plants have close collaborative relationships and whose connection paths need to be preserved for refined modeling. Specifically, the system calculates a collaborative necessity index N for each pair of water plants (i, j). ijThis indicator integrates three dimensions: N ij = O ij •C ij •D ij ; where O ij It is the service area overlap, measured by calculating the size of the user area jointly served by the two water plants; C ij It refers to water quality complementarity. For example, a water plant primarily producing Q1 and a water plant primarily producing Q3 have high complementarity. For a plant like C, which can switch water quality, its complementarity will have an additional bonus. ij This is a distance factor, reflecting the distance between two water plants in the pipeline network; the closer the distance, the greater the likelihood of collaboration. After constructing the collaboration necessity matrix N between all water plant pairs, the system can apply community detection algorithms such as spectral clustering to automatically divide all water plants into several collaboration groups. The sparsity strategy will be: retain all connection paths between water plants within a collaboration group, while for water plants belonging to different groups, only those N paths will be retained. ij Strong collaborative connections whose values ​​exceed a high threshold.

[0114] Macro-level sparsification is performed to extract and preserve the backbone transmission network of the pipeline network. The aim is to identify the main pipelines undertaking long-distance, high-volume transportation tasks from a global perspective of the entire pipeline network and ensure the integrity of these main pipelines. Specifically, the system can comprehensively utilize multiple graph theory algorithms to extract the backbone network: Based on topological centrality: the edge betweenness centrality of each edge is calculated. This metric measures the frequency with which the shortest path between all pairs of nodes passes through a particular edge; edges with high betweenness centrality are typically bridges in the network and must be preserved. Based on flow load: the average flow load of each edge is extracted from the hydraulic model, identifying a few core pipelines carrying, for example, more than 70% of the total network flow as the flow backbone. Based on network connectivity: while preserving all water plants and major user nodes, Kruskal's algorithm is run to generate a minimum spanning tree (MST) that ensures basic network connectivity. The system will merge the edge sets identified by the above methods and apply optimization algorithms such as Louvain community detection to form a hierarchical sparse representation: that is, several regions with relatively dense internal connections, and a few backbone channels connecting these regions.

[0115] In summary, by employing sparsification strategies at the micro, meso, and macro levels, this embodiment intelligently prunes the original complex pipeline network model into a sparse model that retains only the key connections. This sparse model not only improves computational efficiency (e.g., achieving more than 20 times computational speedup), but more importantly, because it retains only carefully selected connections with the highest information carrying capacity, it can largely maintain the prediction accuracy of the original model.

[0116] In another embodiment of this application, a specific, engineering-feasible solution strategy for overcoming the curse of dimensionality is provided for large-scale multi-objective optimization problems. Optionally, the high-dimensional optimization decision space is reduced in dimensionality. Directly solving a mixed-integer programming problem with tens of thousands of dimensions is computationally extremely expensive and difficult to complete within a real scheduling window. Analysis of historical scheduling data reveals that the optimal production scheme is not completely random, but contains specific, repeatable patterns and structures. Using data-driven methods, these inherent structures are discovered, and the decision space is projected into a lower-dimensional, more easily solvable subspace. As a preferred implementation, this dimensionality reduction process may include a combination of the following techniques:

[0117] Linear dimensionality reduction is performed using Principal Component Analysis (PCA). Specifically, the system first extracts the global optimal solutions from the historical database for several past scheduling cycles (e.g., the past year), constructing a historical solution data matrix X, where each column represents a complete decision vector at a given time. The covariance matrix of this data matrix is ​​calculated, and by solving its eigenvalue problem, a set of orthogonal bases that maximizes the explanatory variance of the data—the principal components—is found. Experimental data shows that typically the first k principal components (e.g., k is approximately 200) are sufficient to capture over 95% of the variance in the original data. Therefore, any high-dimensional decision vector x can be projected onto a k-dimensional principal component space using a transformation matrix V, resulting in a low-dimensional representation y = V. T •x.

[0118] In the principal component space, sparse coding is used for further nonlinear compression. To achieve a more extreme compression effect, the system further assumes that any low-dimensional representation y can be linearly combined from a few atoms (i.e., basis vectors) in a dictionary matrix D. That is, it seeks a sparse representation y = Dα, where α is a sparse coefficient vector with most elements being zero. The system can use algorithms such as K-SVD to learn this optimal dictionary D from a large amount of historical data y. When solving new problems, algorithms such as Orthogonal Matching Pursuit (OMP) are used to quickly find the sparse coefficients α. By controlling the sparsity to, for example, 10%, the original k-dimensional vector can be compressed to only need to store k*10% of the non-zero coefficients and their positions; for example, the solution can be reconstructed using only 20 non-zero coefficients.

[0119] Optionally, time series patterns can be mined as dedicated basis functions. The system can identify recurring typical patterns in production planning, such as weekday peak patterns and weekend flat patterns. Specifically, techniques such as Symbolic Aggregation Approximation (SAX) can be used to discretize numerical time series into symbol strings, and then data structures such as suffix trees can be used to efficiently find these frequently occurring patterns. By using, for example, 12 identified basic patterns as a set of basis functions, any complete production plan can be approximately represented as a linear combination of these basis functions.

[0120] The optimization problem is solved in a reduced-dimensional space, and the result is reconstructed. In the above steps, the original optimization problem in a 21600-dimensional space is transformed into an equivalent optimization problem in a new space with significantly reduced dimensions (e.g., only 200 dimensions or lower). The system can efficiently use standard optimization solvers (such as the interior-point method or the simplex method) in this low-dimensional space to find the optimal low-dimensional solution (e.g., the optimal sparse coefficients α or basis function combination coefficients β). After obtaining the optimal low-dimensional solution, a reverse reconstruction process is performed, mapping the low-dimensional solution back to the original 21600-dimensional decision space, i.e., x = V•D•α, using the saved transformation matrix, dictionary, and basis function library. Since dimensionality reduction is lossy compression, the reconstructed solution may contain minor constraint violations. Therefore, a feasibility correction step is finally needed to fine-tune the reconstructed solution to ensure that it strictly satisfies all the original constraints, while controlling the reconstruction error within a very small range (e.g., less than 2%).

[0121] Furthermore, to improve solution efficiency, a warm-start optimization strategy can be adopted. Specifically, the system maintains a database of optimized solutions from the past 7 days. When solving the problem for the current day, the system first uses time series similarity algorithms such as Dynamic Time Warping (DTW) distance to find the historical scenario most similar to the current predicted demand pattern from the solution database, and uses its corresponding historical best solution as the starting point for this optimization iteration. Compared to randomly selecting the initial point, this warm-start based on historical experience allows the solver to converge to the optimal solution faster, improving the convergence speed.

[0122] Through the above steps, this embodiment successfully transforms the theoretically difficult-to-handle ultra-large-scale optimization problem into a computational task that is practically feasible in engineering and can be solved quickly, thus providing performance assurance for the practical application of the present invention.

[0123] In another embodiment of this application, a series of optional enhanced implementation methods are provided for the model building and collaborative scheduling protocol involved in this application, which can deeply couple the unique physical, chemical, and operational characteristics of the reclaimed water pipeline network. These methods can improve the fidelity of the model and the practicality of the solution, thereby providing strong support for the creativity of the overall solution.

[0124] Optionally, a dynamic growth model of the biofilm on the pipe wall can be introduced when constructing the dynamic graph network model. During long-term operation, a biofilm will grow on the inner wall of the pipe. To simulate this process more precisely, a differential equation model describing the evolution of the biofilm thickness Δ over time can be introduced: dΔ / dt = μ max •(S / (K s +S))•Δ- b•Δ-τ w •Δ 2 / ρ b Where Δ is the biofilm thickness; μ max It is the maximum growth rate of microorganisms; S is the concentration of the substrate (such as COD); K is the maximum growth rate of microorganisms. s τ is the half-saturation constant; b is the rate of decline of endogenous respiration in biological membranes; w It is the shear force of the water flow on the pipe wall; ρ b This refers to biofilm density. Using real-time water quality, temperature, and hydraulic data from the pipe network, the equation is solved online using numerical methods such as the Runge-Kutta method to obtain the dynamic biofilm thickness Δ(t) for each pipe section. Δ(t) is then used to correct the effective pipe diameter D used in the model in real time. eff The model uses n(t) and hydraulic roughness n(t) to automatically adapt to the aging process of the pipeline network and make more accurate long-term predictions.

[0125] Optionally, the propagation velocity specific to water quality grade can be considered when constructing the conduction core. Traditional models assume that all water flows at the same velocity in the pipe, but studies have found that reclaimed water of different qualities exhibits significant differences in fluid properties. Q1 grade high-quality water is close to a Newtonian fluid, with an actual flow velocity v. Q1 The actual flow rate is very close to the theoretical value (e.g., 98% of the theoretical value); Q2 grade medium-quality water exhibits non-Newtonian fluid characteristics due to the presence of colloidal particles, resulting in increased apparent viscosity and a lower actual flow rate than the theoretical value (e.g., 92% of the theoretical value); Q3 grade low-quality water exhibits Bingham plastic flow characteristics due to the presence of suspended solids, resulting in yield stress and a further reduction in actual flow rate (e.g., 85% of the theoretical value). The model can automatically select the corresponding effective propagation velocity based on the water quality grade being transported when calculating the physical transport delay of the conduction core, thereby more accurately predicting the arrival time of water masses of different qualities.

[0126] Optionally, physicochemical constraints on water mixing can be introduced during sparsification or optimization. To ensure water quality safety, the system can establish a water mixing compatibility matrix M. compat (Q i Q j Used to quantify two different qualities of water, Q. i and Q jThe risk of incompatibility during mixing. For example, when Q1 water and Q3 water are directly mixed, their compatibility is close to 0, indicating that mixing is strictly prohibited. When performing multi-path flow allocation or path planning, M can be considered... compat > M threshold (Safety threshold) serves as a hard constraint. When the optimal path presents mixed risks, the system can automatically trigger an isolation delivery mode, such as staggered delivery in time or activation of backup pipelines in space, to prevent water quality accidents at the model level.

[0127] Optionally, a water network-specific coordination mechanism can be introduced when designing the coordinated scheduling protocol. Unlike power grids or communication networks, water, as a storable, quality-sensitive, and time-delayed resource, allows for more refined coordinated scheduling: Water quality buffer coordination mechanism: For buffer units such as pools or tanks, the scheduling logic should not be a simple first-in-first-out (FIFO) approach, but rather a dual FIFO + quality-priority scheduling logic. The system dynamically determines the batch of water to be released based on the freshness of the water in the pool (determined by residence time) and the urgency of downstream demand, achieving a balance between meeting demand and ensuring water quality. Time-delay compensation pre-scheduling credit mechanism: To incentivize water plants to proactively and accurately pre-produce (rather than passively respond), a pre-scheduling credit mechanism can be introduced. Water plants can earn credit points by producing in advance based on predicted demand; settlement is made after actual demand is confirmed, with credits awarded for accurate predictions and deducted for deviations. Credit points can be used to obtain priority or electricity price discounts in future scheduling. Dynamic alliances for water quality compatibility: The system can facilitate dynamic alliances between water plants based on real-time water quality monitoring results. For example, when Plant A's Q1 water quality exceeds expectations and reaches the Q1+ level, the system can automatically match it with a special user requiring ultra-high-quality water, forming a value-added alliance. When Plant B's water quality inadvertently deteriorates, it can automatically form a complementary alliance with Plant C, with Plant C providing high-quality water for dilution. This flexible and dynamic alliance is something that traditional static scheduling cannot achieve.

[0128] In an exemplary embodiment, a formation mechanism, mathematical definition, and parameter determination method for a bimodal conduction kernel in the presence of multiple independent decision paths are provided. Specifically:

[0129] Based on the dynamic pipeline topology in a standardized spatiotemporal dataset, multiple hydraulic transport paths between pipeline node pairs are explored and identified.

[0130] In this embodiment, the multiple hydraulic transmission paths include not only physical pipeline paths but also corresponding decision-making and control paths. Specifically, for the water supply process from water plant S to user D, there are three independent decision-making triggering mechanisms: the first is automatic control decision-making based on the SCADA system, with a response time of approximately 5-15 minutes; the second is manual intervention decision-making by the dispatcher, with a response time of approximately 30-60 minutes; and the third is cross-departmental coordination decision-making, which, when involving multiple water plants, can have a response time of 2-4 hours. These decision-making paths exhibit significant differences in response time, forming the basis for a multimodal distribution of decision delays.

[0131] A path conduction kernel is constructed for each hydraulic transport path. The path conduction kernel exhibits a bimodal distribution to characterize the physical transmission delay and the decision response delay, respectively.

[0132] Specifically, the mathematical definition of the bimodal conduction kernel adopts a mixed gamma distribution model: K p (τ) = w 1p •Γ(τ;k 1p θ 1p ) + w 2p •Γ(τ;k 2p θ 2p ); where Γ(τ; k, θ) = (τ (k-1) ×exp(-τ / θ)) / (θ k ×Γ(k)). The first peak corresponds to the physical transport process, and its parameters are determined by hydraulic principles. Specifically, the shape parameter k 1p Related to the pipe Reynolds number Re: when Re > 4000 (turbulent flow), k 1p = 2.0 + 0.0001×Re; when 2000 < Re ≤ 4000 (transitional flow), k 1p = 1.5 + 0.00025 × Re; when Re ≤ 2000 (laminar flow), k 1p = 1.2 + 0.0004 × Re. Scale parameter θ 1p = L p / (v p ×3.6), where L p v represents the path length (in kilometers). p The average flow velocity is 3.6 (m / s), and 3.6 is the unit conversion factor. The second peak corresponds to the decision-response process, which arises from the superposition effect of multiple independent decision paths. In this embodiment, when path p involves n decision nodes, each node i has an independent response time T. i The convolutions with independent response times follow an exponential distribution. These convolutions produce a gamma distribution with shape parameter k. 2p =n, reflecting the complexity of the decision chain; scale parameter θ2p = Σ(T i ) / n, representing the average decision delay.

[0133] Furthermore, the weighting coefficients of the bimodal conduction kernel are determined. The weighting coefficient w 1p and w 2p The determination is based on maximum likelihood estimation using historical operational data. The specific method is as follows: collect water quality transport events along this path over the past 90 days, and record the actual arrival time series {τ1, τ2, ..., τ} for each event. N Construct the likelihood function L(w) 1p w 2p ) =Π[w 1p × Γ(τ i ;k 1p θ 1p )+ w 2p ×Γ(τ i ;k 2p θ 2p The optimal weights are obtained through iterative solution using the Expectation-Maximization (EM) algorithm. In practical applications, for emergency water supply routes (such as hospital water supply lines), the decision-making response is accelerated by a prioritization mechanism, leading to... 2p Smaller (approximately 0.2-0.3); for conventional water supply routes, the impact of physical transport and decision response is comparable, w 1p and w 2p The value is close to 0.5; for cross-regional coordinated water supply routes, the decision-making and coordination process is time-consuming. 2p The weighting can reach 0.6-0.7. Optionally, the weighting coefficient can also be dynamically adjusted according to the time period. During working hours (8:00-18:00), dispatchers are on duty, and decision-making and response are faster. 2p =0.3; During the nighttime period (22:00-6:00), the system mainly relies on automatic control, and the delay in manual decision-making increases. 2p =0.5; During holidays, cross-departmental coordination is difficult, w 2p It can be improved to 0.6.

[0134] Furthermore, to verify the rationality of the bimodal distribution, the system performs the Kolmogorov-Smirnov test, comparing the goodness of fit between the actual transmission delay distribution and the bimodal model. The bimodal model is considered valid when the test statistic D < 0.05 and the p-value > 0.95. If the test fails, the system automatically switches to a trimodal or multimodal hybrid model to more accurately capture the multimodal delay characteristics in complex decision-making scenarios.

[0135] In another exemplary embodiment, a method for coupling sparsification processing with water quality safety constraints is provided, specifically as follows:

[0136] Based on water quality data in a standardized spatiotemporal dataset, this study quantifies and assesses the risk of water incompatibility when water quality mixes at the two ends of a dynamic graph network prediction model.

[0137] Optionally, the quantification of water incompatibility risk adopts a multi-dimensional assessment system. A chemical reaction risk matrix R is established. chem Its element R chem (i, j) represents the reaction risk value when water is mixed at nodes i and j. This value is calculated using the Gibbs free energy: ΔG = ΔH - T×ΔS; where ΔH is the enthalpy change of the mixing reaction (kJ / mol); T is the water temperature (Kelvin); and ΔS is the entropy change of the mixing reaction (J / mol·Kelvin). When ΔG < -10 kJ / mol, the reaction proceeds spontaneously, and the risk value R is... chem = exp(-ΔG / RT), where R is the gas constant 8.314 J / (mol·K). The specific risk assessment includes three levels: the first level is the immediate chemical reaction risk, such as residual chlorine (concentration C). cl ) and ammonia nitrogen (concentration C) nh The mixture produces chloramine, with a risk value R. instant = k1×C cl ×C nh ×exp(-E a / RT); where k1 is the reaction rate constant 0.15 L / (mg·min); E a The activation energy is 25 kJ / mol. The second layer represents the risk of slow degradation, assessing the probability of water quality deterioration within 24 hours after mixing, R. degrade = 1 - exp(-λ×t×ΔQ); where λ is the degradation rate constant (0.02 / h); t is the assessment time window (24h); and ΔQ is the quality grade difference. The third layer is the risk of microbial cross-contamination, R bio =(N i ×N j ) / (N safe 2 ) ×exp(μ×T); where N i N j The total number of bacteria at both nodes (CFU / mL); N safe The safety threshold is 100 CFU / mL; μ is the growth rate of 0.03 / ℃; T is the water temperature.

[0138] Construct a coupled optimization model of sparse decision-making and security constraints.

[0139] In this embodiment, sparsity is no longer simply about pursuing computational efficiency, but is deeply coupled with security constraints. Specifically, a multi-objective optimization problem is constructed: min F = α × C comp +β×Rtotal +γ×L service ; where C comp The computational complexity is equal to the square of the number of edges retained; R total = Σ(w ij × R ij ) represents the total risk value, w ij Let L be the flow weight of edge (i, j); service The service loss is used to measure the impact of sparsity on water supply reliability; α, β, and γ are weighting coefficients, adjusted according to the application scenario. The constraints of this optimization problem include: connectivity constraints, ensuring all user nodes are reachable; capacity constraints, the total transport capacity of the reserved paths is not less than 120% of the demand; and hard safety constraints, any R... ij > R threshold The edges must be cut off, where R threshold Determined based on downstream user type: Hospital users R threshold = 0.1; Industrial user R threshold = 0.3; Municipal greening R threshold = 0.6.

[0140] Achieve dynamic and coordinated adjustment of sparsity and security constraints.

[0141] Optionally, the system employs a two-layer iterative algorithm to achieve coupling optimization. The outer loop adjusts the sparsity parameter ρ, starting from an initial value ρ0 = 0.5 with a step size Δρ = 0.05. The inner loop, with a fixed sparsity, selects the set of edges E to be retained using a tabu search algorithm. sparse In each iteration, the algorithm evaluates three types of metrics: computational speedup S = T full / T sparse Number of safety violations V = count(R) ij >R threshold for(i, j) in E sparse Service quality index Q = Σ(demand) met / demand total ). Among them, T full T represents the total computation time required for the system to run one prediction or optimization task in the full graph structure without sparsification. sparse This refers to the computation time required for the system to perform the same task under the current sparse graph structure; demand met This represents the total water supply demand of user nodes that can be successfully met under the current sparse graph structure. totalThis represents the theoretical total water supply demand for all user nodes. When a safety violation (V > 0) is detected, the algorithm immediately triggers a safety-first mode: adding the violating edge to the taboo list, preventing it from being selected in subsequent iterations; activating alternative path search to find alternative routes bypassing high-risk areas; and reducing sparsity requirements and adding redundant connections if necessary to ensure safety. This dynamic adjustment mechanism ensures that sparsity reduction never comes at the expense of safety.

[0142] In the actual tests of this embodiment, compared with the traditional independent sparsification method, the coupled optimization scheme reduced the incidence of water quality safety incidents while maintaining the same computational efficiency. This is due to the proactive identification and avoidance of risk factors during the sparsification process, rather than post-event remediation. Optionally, a reinforcement learning mechanism can also be introduced to train the sparsification policy network using historical operating data. This network takes the pipeline state (topology, flow rate, water quality) as input and outputs the retention probability of each edge. The reward function is designed as: r = -0.4×C comp_normalized -0.5×R total_normalized -0.1×L service_normalized +10 ×I(V=0); where I(V=0) is a safety indication function, taking the value 1 when there is no safety violation, and 0 otherwise; C comp_normalized To calculate the normalized value of the complexity, R total_normalized L is the normalized result of the total risk value. service_normalized This is a normalized value for the service loss. Through training with a Deep Q-Network (DQN), the system can learn the optimal sparsity strategy that balances efficiency and security.

[0143] In another exemplary embodiment, a specific process for deeply integrating a multipath conduction model with a water quality cascade utilization network is provided. Specifically, a cascade utilization delay matrix is ​​constructed based on a multipath conduction kernel. In the cascade utilization network, the effluent from the upstream user becomes the water source for the downstream user, and the accuracy of delay prediction directly affects the stable operation of the cascade chain. The multipath conduction model provides a complete delay probability distribution K for each cascade connection (i→j). ij (τ), rather than a single delay value. Based on this, a three-dimensional stepped delay tensor T[i, j, τ] is constructed, where the elements T ijτ This represents the probability density of water flowing from node i to node j after a time delay τ. The construction of this tensor considers the coupling effect unique to cascade systems. When user i simultaneously supplies water to multiple downstream users j and k, there is flow competition, and the actual flow rates of each path will affect each other. This is addressed by solving the coupled hydraulic balance equations: Σ(Q ij ) =Q i_out (Flow conservation); h i - h j = f ij (Q ij (Energy conservation); where Qi_out The total water output of user i; h is the node head; f ij Let Q be the pipe resistance function. Solving for Q yields the dynamic flow distribution. ij (t), and then update the conduction kernel parameters of each path. Multipath information is used to optimize cascade link selection. Traditional cascade matching only considers water quality and quantity matching, neglecting the importance of time delay matching. In this embodiment, the cascade feasibility scoring function is extended to: S cascade (i, j) = M quality ×M quantity ×M timing ×(1 - R risk ); where M quality = exp(-|Q i_out - Q j_req | / Q ref ) represents the quality matching degree; M quantity =min(V i_out / V j_req 1) represents the water quantity matching degree; M timing For latency matching degree; R risk For path risk, Q j_req Let Q be the target water supply demand flow rate for node j (the downstream user). ref V is a reference flow rate value. i_out V represents the available water volume for node i (the parent user) during the current time period. j_req Let M be the target water volume for node j within the current time period. Delay matching degree M timing The calculation fully utilizes multipath transmission information: M timing = ∫K ij (τ)×W j (τ) dτ; where W j (τ) is the time window function for user j's demand. When K ij The main peak of (τ) and W j When the high demand periods of (τ) overlap, M timing A value close to 1 indicates a good time delay match. The multipath model provides a complete distribution K. ij (τ) Compared to a single delay value, it can more accurately assess the degree of time window matching. Experimental data shows that introducing the time delay matching degree improves the average utilization rate of the cascade chain and reduces flow interruption events. This is because the multipath model identifies the long-tail characteristics of the time delay distribution, avoiding supply interruptions caused by low-probability but long-delay water flows. Multipath redundancy enhances the resilience of the cascade network. The path diversity revealed by the multipath transmission model provides a natural redundancy backup for the cascade network. Specifically, the path redundancy index R is calculated for each cascade connection. redundancy (i, j) = 1 - Π(1 - p k ); where pk Let R be the probability that the k-th path is available. redundancy When the threshold is >0.8, even if the primary path fails, there are still sufficient backup paths to maintain cascade operation. Based on this, the system constructs a resilience enhancement strategy for the cascade network. During the network design phase, node pairs with high path redundancy are prioritized for constructing cascade chains. During operation, the status of each path is monitored in real time, and when a path performance degradation is detected, traffic allocation is preemptively adjusted to smoothly transition to backup paths. This proactive resilience management based on multipath information shortens the average recovery time of the cascade network.

[0144] Optionally, the system can also implement multi-path collaborative scheduling. For critical cascade connections, multiple paths are activated simultaneously, and flow smoothing is achieved through phase modulation between paths. Specifically, the transport phase φ of each path is set. k = 2π×k / N, where N is the total number of paths. The instantaneous flow Q of each path. k (t) = Q avg + A×sin(ωt + φ k ), where Q avg Let A be the average flow rate, ω be the modulation amplitude, and ω be the modulation frequency. The superposition of multipath flows generates a stable total flow rate, reducing flow fluctuations in the cascade chain and improving the water supply stability for downstream users. Through the above mechanism, the combination of multipath transmission and cascade utilization produces a synergistic enhancement effect: the fine-grained delay information provided by the multipath model optimizes cascade matching and improves network efficiency; path redundancy enhances the resilience and reliability of the cascade network; and multipath collaborative scheduling smooths cascade flow and improves water supply quality. This results in an overall system performance improvement that exceeds the simple superposition of individual technologies, demonstrating true technological integration.

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

Claims

1. A method for coordinated forecasting of urban reclaimed water supply and demand, characterized in that, include: Process multi-source heterogeneous data streams from reclaimed water pipeline networks to generate standardized spatiotemporal datasets; Based on standardized spatiotemporal datasets, a dynamic graph network prediction model is constructed to characterize the supply and demand relationship and spatiotemporal transmission characteristics of water quality among nodes in the pipeline network. By using dynamic graph network prediction models and real-time demand signals, we can generate quality-specific demand prediction results and cross-plant collaborative production plans. Constructing a dynamic graph network prediction model includes: Based on the dynamic pipeline topology in the standardized spatiotemporal dataset, the predetermined hydraulic transport paths between pipeline node pairs are explored and identified. A path conduction kernel is constructed for each hydraulic transport path. The path conduction kernel exhibits a bimodal distribution to characterize the physical transmission delay and the decision response delay, respectively.

2. The method according to claim 1, characterized in that, For each hydraulic transport path, a path conduction core is constructed, including: Extract the water quality grade of reclaimed water associated with the path from the standardized spatiotemporal dataset; Based on the water quality grade, determine the parameters of the peak in the bimodal distribution that characterize the delay in decision response.

3. The method according to claim 2, characterized in that, Building a dynamic graph network prediction model also includes: Based on the hydraulic parameters of the pipeline network in the standardized spatiotemporal dataset, the flow distribution ratio of a predetermined hydraulic transport path is analyzed. Based on the flow distribution ratio, the path conduction cores of each hydraulic transport path are weighted and superimposed to form a comprehensive conduction core that represents the network node pairs.

4. The method according to claim 1, characterized in that, Constructing a dynamic graph network prediction model also includes determining the time-delay propagation kernel, specifically: Based on hydraulic and water quality parameters in standardized spatiotemporal datasets, a one-dimensional convection-reaction-diffusion equation is established for pipeline node pairs. Solve the one-dimensional convection-reaction-diffusion equation and use its analytical or numerical solution as the time-delayed conduction kernel.

5. The method according to claim 4, characterized in that, Establishing a one-dimensional convection-reaction-diffusion equation, including parameterization based on a standardized spatiotemporal dataset, wherein the parameterization includes: The reaction degradation coefficient in the one-dimensional convection-reaction-diffusion equation is determined as a function related to pollutant concentration and water quality level. A dynamic growth model of biofilm on the pipe wall was established to dynamically correct the effective pipe diameter and hydraulic roughness used in the one-dimensional convection-reaction-diffusion equation. Differential propagation speeds are set based on the one-dimensional convection-reaction-diffusion equation according to water quality grades, wherein the propagation speed of high-quality water in the central region of the pipe is higher than that of low-quality water in the near-wall region.

6. The method according to claim 1, characterized in that, Also includes: The dynamic graph network prediction model is sparsified to generate a sparse connection matrix. Among them, the generation of quality demand forecast results and cross-plant collaborative production schemes are carried out using sparse connection matrices.

7. The method according to claim 6, characterized in that, Generate a sparse connectivity matrix, including: Based on the topological change rate, demand volatility and water quality heterogeneity in the standardized spatiotemporal dataset, the real-time system complexity index C(t) of the pipeline network is quantitatively evaluated. Based on the real-time system complexity index, through ρ(t)=ρ base The functional relationship of ×exp(-β×C(t)) determines the connection retention ratio ρ(t), generating a sparse connection matrix, where ρ base The base retention ratio is given by β, which is a sensitivity parameter.

8. The method according to claim 6, characterized in that, The sparse connectivity matrix can also be generated as follows: Based on water quality data in a standardized spatiotemporal dataset, the risk of water incompatibility when water quality mixes at the two ends of a dynamic graph network prediction model is quantitatively assessed. Disconnect or weaken network connections where the risk of water quality incompatibility exceeds a preset mixing risk threshold to generate a sparse connection matrix.

9. The method according to any one of claims 1 to 8, characterized in that, Generating cross-plant collaborative production plans is achieved by solving a large-scale multi-objective optimization problem, which includes the following steps: The large-scale multi-objective optimization problem is decomposed into a predetermined number of low-dimensional sub-problems along the dimensions of water quality grade, time window and spatial region; By using a distributed coordination protocol, the solutions to low-dimensional subproblems are coordinated and integrated to obtain the global optimal solution to large-scale multi-objective optimization problems. Among them, the large-scale multi-objective optimization problem is constructed based on a dynamic graph network prediction model and real-time demand signals.

Images