Urban reclaimed water supply and demand collaborative prediction method

Abstract

The invention discloses an urban reclaimed water supply and demand collaborative prediction method, which comprises the following steps: processing a multi-source heterogeneous data flow from a reclaimed water pipe network, and generating a standardized spatio-temporal data set; a dynamic graph network prediction model is constructed based on the data set, a plurality of hydraulic conveying paths between nodes are explored, a bimodal form path conduction kernel capable of representing physical transmission delay and decision response delay respectively is constructed for each path, and then a comprehensive conduction kernel is formed through weighted stacking so as to accurately represent water quality space-time conduction characteristics; and generating a sub-quality demand prediction result and a cross-factory collaborative production scheme by using the model and the real-time demand signal. According to the method, the accuracy of prediction, the economical efficiency of scheduling and the safety and value of renewable water resource utilization can be improved by depicting the hydraulic power-water quality coupling conduction rule in a refined manner and combining an efficient collaborative optimization strategy.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of urban reclaimed water utilization, in particular to a method for predicting the supply and demand of urban reclaimed water. BACKGROUND

[0002] The large-scale utilization of reclaimed water is not simply transporting it from the treatment plant to the user, but a complex system engineering involving multi-source production, multi-quality, multi-user, and dynamic demand. How to accurately predict the supply and demand situation of reclaimed water of different qualities, and on this basis, realize the coordinated production among multiple water plants and the optimization of the pipe network, to meet the diversified water demand of the city at the lowest cost, the highest efficiency and the safest guarantee, is the key to maximizing the value of urban reclaimed water utilization, and is also the core technical issue to promote the transformation and upgrading of urban water system to fine, intelligent and green.

[0003] Currently, researchers have carried out a lot of work on the prediction and scheduling of urban water supply system. In terms of demand prediction, existing methods mainly rely on time series analysis models, such as autoregressive integrated moving average model (ARIMA), exponential smoothing method and other traditional statistical methods, as well as machine learning-based methods such as support vector machine (SVM), artificial neural network (ANN), long short-term memory network (LSTM) and other methods. These methods can better capture the inherent trend and periodicity of historical water consumption data of a single user or regional total when processing such data. In terms of pipe network scheduling and hydraulic simulation, the mainstream technology is simulation optimization based on EPANET and other hydraulic models. Such models can accurately calculate the hydraulic parameters such as pressure and flow in the pipe network, and through the combination with intelligent optimization algorithms such as genetic algorithm and particle swarm optimization, the optimization of scheduling schemes such as pump station start-stop and valve switching can be realized. In terms of water quality simulation, researches are mainly focused on the one-dimensional transport and first-order kinetic decay process of a single pollutant (such as residual chlorine) in the pipe network, providing certain technical support for ensuring water supply safety.

[0004] Although existing technologies have improved the scientificity 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 realistic predicament of inaccurate prediction, poor coordination and low efficiency. In summary, the existing technologies mainly have the problems of distorted representation of hydraulic-water quality coupling transmission characteristics, imbalance between computational efficiency and security of complex network models, and limited decision-making dimension of multi-objective coordinated scheduling strategies. SUMMARY

[0005] The present application provides a method for predicting the supply and demand of urban reclaimed water, to solve the above-mentioned problems existing in the prior art.

[0006] Technical scheme, a city reclaimed water supply and demand collaborative prediction method, comprising:

[0007] Processing a multi-source heterogeneous data stream from a reclaimed water pipe network to generate a standardized spatiotemporal data set;

[0008] Based on the standardized spatiotemporal data set, a dynamic graph network prediction model is constructed to represent the supply and demand relationship between nodes in the pipe network and the spatiotemporal conduction characteristics of water quality;

[0009] Using the dynamic graph network prediction model and real-time demand signals, a product quality demand prediction result and a cross-plant collaborative production scheme are generated.

[0010] Beneficial effects, the present application can improve the accuracy of prediction, the economy of scheduling, and the safety and value of reclaimed water resource utilization by precisely depicting the hydraulic-water quality coupling conduction law and combining with an efficient collaborative optimization strategy. BRIEF DESCRIPTION OF DRAWINGS

[0011] Figure 1 A step flow chart of the city reclaimed water supply and demand collaborative prediction method provided by the embodiment of the present application.

[0012] Figure 2 A step flow chart of constructing a dynamic graph network prediction model provided by the embodiment of the present application.

[0013] Figure 3 A step flow chart of constructing a path conduction kernel for each hydraulic transportation path provided by the embodiment of the present application.

[0014] Figure 4 A step flow chart of determining a time-lag conduction kernel provided by the embodiment of the present application. DETAILED DESCRIPTION

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

[0016] It should be noted that the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device including a series of steps or units need not be limited to only those steps or units clearly listed, but can include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0017] In the study, it is found that there is a distortion problem in the characterization of the hydraulic-water quality coupled conduction characteristics. The existing model generally uses a single delay assumption based on the shortest path and average flow rate to describe the transport process of water from source to sink. However, in real urban pipe networks, especially in ring networks, water flow will be transported along multiple paths at the same time, and the single path assumption deviates seriously from the physical reality, leading to a huge deviation in the prediction of water quality arrival time. In addition, the existing model ignores a key non-physical delay - decision response delay, which is the time required for a series of human or system decision-making processes from receiving dispatching instructions to actually adjusting production and opening valves to transport water in the water plant. This systematic neglect of multi-path effect and decision delay makes the model unable to accurately depict the real-time and space dynamic distribution of water quality signals in the pipe network, directly leading to the failure of accurate matching between supply and demand in time. Secondly, the above characterization distortion problem further exacerbates the imbalance between the computational efficiency and safety of the complex network model. On the one hand, in order to more realistically simulate the pipe network, the model needs to include thousands of nodes and connections, leading to a dimension disaster, and the computational complexity of the fully connected model increases exponentially, making it difficult to meet the high efficiency requirements of real-time scheduling. On the other hand, when simplifying the model (such as sparsification), the existing method lacks the guidance of domain knowledge and may incorrectly retain or cut off connections. In particular, it fails to fully consider the chemical incompatibility risk or quality degradation risk that may occur when different quality reclaimed water is mixed. This blind sparsification process may result in a theoretically optimal scheduling scheme, but in reality it may trigger serious water quality safety incidents. Finally, the above two problems together lead to the problem of limited decision dimension of multi-objective collaborative scheduling strategy. Traditional scheduling optimization is usually single-objective or low-dimensional, such as only pursuing the lowest cost, and is usually carried out independently among water plants, lacking a global collaborative perspective. It fails to take the unique quality dimension of reclaimed water as the core optimization variable, and cannot effectively utilize the differences in water quality requirements of different users to construct a step-by-step utilization network that makes the best use of water. At the same time, in the face of the super-high-dimensional decision space across plants, regions, qualities, and time periods, traditional centralized optimization solvers are difficult to find a global optimal solution within an effective time. These limitations ultimately lead to the waste of high-quality reclaimed water in low-demand scenarios and the difficulty in improving the overall operational efficiency of the system.

[0018] In some embodiments, the implementation of the urban reclaimed water supply and demand collaborative prediction method can be based on a system that integrates data acquisition, transmission, storage, and computation. The system can include online monitoring devices (such as flow meters, pressure sensors, water quality analyzers) deployed at key locations in the reclaimed water pipe network, a distributed supervisory control and data acquisition system (SCADA system), production scheduling systems in each reclaimed water plant, a user feedback platform, and a central server or cloud computing platform for executing the methods of the present application.

[0019] As Figure 1The urban reclaimed water supply-demand collaborative prediction method shown includes the following steps:

[0020] Processing the multi-source heterogeneous data stream from the reclaimed water pipe network to generate a standardized spatiotemporal data set.

[0021] In other words, the multi-source heterogeneous data stream from the reclaimed water pipe network is obtained, preprocessed, and a standardized spatiotemporal data set is generated.

[0022] For example, the multi-source heterogeneous data stream can include real-time pipe network monitoring data from a SCADA system, production scheduling records of water plants, topological data of the urban pipe network GIS system, and historical water use records of users in various fields. The process of generating a standardized spatiotemporal data set is to convert these data of different sources, different formats, and non-uniform spatiotemporal dimensions into high-quality, structured data that can be directly used by subsequent models. Specifically, this step can include the following processing processes: collecting and cleaning real-time monitoring data stream. From the distributed SCADA system, the original monitoring data stream containing flow, pressure, and COD (chemical oxygen demand), turbidity, ammonia nitrogen, pH, residual chlorine, conductivity, and other water quality parameters can be extracted. In order to process data distortion caused by sensor failure or communication anomaly, an outlier detection algorithm based on sliding window can be used to identify abnormal data points. For the identified outliers, instead of simply discarding them, they are repaired using the data of their adjacent normal time points through a cubic spline interpolation method, and at the same time the system marks the abnormal period for subsequent fault analysis. After this processing, a cleaned monitoring data sequence with a time interval of 5 minutes and covering all pipe network nodes in space can be formed. Structured integration of production operation data. For example, the system can integrate the production scheduling records of three reclaimed water plants, which contain information such as production capacity, water quality grade switching, and operation cost at each time period. Since the data formats of different water plants may be heterogeneous, a unified data mode needs to be constructed to convert the text form of scheduling instructions (such as 8:00 am, C plant switches to Q1 high-quality mode) into a structured production capacity time sequence. In particular, for the C plant that can switch production of different water quality grades, its switching process is not instantaneous, in order to avoid the impact of step changes on the prediction model, the system can refine it into a gradual production curve during the transition period. After integration, a standardized production capacity matrix with an hour granularity can be generated, where the rows can represent timestamps and the columns can represent the production capacity of each water plant at different water quality grades.

[0023] The system will analyze the static topology data in the city pipe network GIS system, such as node coordinates, pipe segment connection relationships, pipe diameter, pipe age, and other attributes. At the same time, in order to reflect the real state of the pipe network, the system will also integrate the pipe network reconstruction construction plan provided by the engineering department, thereby building a time-varying pipe network topology expression. For example, the changes in each pipe network reconstruction period can be encoded as a sequence of graph structure addition and deletion operations. For temporary pipelines laid during construction, due to their material and laying conditions, they can be assigned a standard transport capacity coefficient of, for example, 0.6 times. Through this topology evolution algorithm, the system can generate a dynamic pipe network topology sequence covering the entire prediction time domain, where each time slice contains the node set, connected edge set, and their respective transport capacity parameters that are valid at that time. The historical demand pattern is extracted in layers. The system will aggregate the historical water consumption records of industrial, agricultural, municipal, and landscape users to form the original demand data set. In order to better understand water consumption behavior, the EMD (Empirical Mode Decomposition) method can be used to decompose the demand sequence of each field into trend items (reflecting long-term changes), periodic items (such as daily, weekly, and seasonal fluctuations), and random items (reflecting short-term disturbances). Among them, the periodic items can be further analyzed by Fourier transform to extract the amplitude and phase of each periodic component. For different water quality grades (for example, Q1 grade for high-standard industry, Q2 grade for agricultural irrigation, and Q3 grade for municipal miscellaneous use), the system will calculate their usage proportion and time distribution characteristics in different fields, and finally build a hierarchical demand pattern feature library. This feature library can contain, for example, a 128-dimensional feature vector to finely describe the water consumption behavior patterns of various users.

[0024] All the above processed data are standardized in space-time dimension. Since data from different sources have different sampling frequencies and spatial resolutions, the system needs to align them to a unified benchmark. For example, the bilinear interpolation method can be used to unify data with different sampling frequencies (such as hourly production data) to a 5-minute benchmark time granularity. In the spatial dimension, for discrete monitoring point data, the inverse distance weighted interpolation algorithm can be used to project its values to all nodes of the pipe network graph, thereby obtaining a fully covered spatial data field. The final output standardized space-time data set can have a four-dimensional tensor structure [time x space x attribute x water quality grade] as the unified and standardized input for all subsequent prediction models. Through the above steps, the original, chaotic multi-source heterogeneous data is processed into a high-value data set with clear structure, space-time alignment, and complete information, laying a solid foundation for subsequent accurate prediction and collaborative optimization.

[0025] Based on the standardized space-time data set, a dynamic graph network prediction model is constructed to represent the supply-demand relationship between nodes in the pipe network and the space-time conduction characteristics of water quality.

[0026] In this embodiment, the dynamic graph network prediction model is a mathematical abstraction and dynamic simulation of a real-world urban recycled water pipeline network. Specifically, the intersections of water plants, users, and the pipeline network can be mapped as nodes in the graph, and these nodes are assigned corresponding attributes, such as node type (production node, consumption node, transit 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 take into account factors such as the physical length, diameter, and wall roughness of the pipe segments. The model is dynamic because it can reflect real-time changes in the pipeline network topology. For example, when a pipe segment enters construction, the model automatically removes the corresponding edge and may activate a pre-defined temporary alternative path. The core of this model lies in characterizing the spatiotemporal conductivity of water quality. In other words, the model not only determines whether nodes are connected, but more importantly, it can predict how long it will take for water of a certain quality, flowing from a node (e.g., a water plant), to reach another node (e.g., a user) and in what form (e.g., whether the concentration will decrease due to diffusion). To achieve this function, a time-delayed conduction kernel can be introduced to quantitatively describe the spatiotemporal dynamic relationship of water quality transmission between any two nodes.

[0027] Utilize dynamic graph network prediction models and real-time demand signals to generate quality-based demand forecast results and cross-factory collaborative production plans.

[0028] In this embodiment, the model's predictive capabilities are translated into practical production scheduling decisions. The system integrates the spatiotemporal correlation structure provided by the dynamic graph network prediction model with real-time demand signals (for example, an industrial park will require a large amount of Q1 recycled water within the next two hours) to construct a multi-objective optimization model that takes into account time-lag constraints. The optimization model's objectives are multidimensional, including, for example, minimizing total production costs; maximizing supply-demand matching to meet the needs of all users as much as possible; optimizing water quality matching to avoid waste caused by high-quality water being used in low-demand scenarios; and balancing pipeline network transmission loads to prevent excessive pressure in some pipelines. By solving this large-scale optimization problem, the system ultimately outputs two core results: first, quality-based demand forecasts. For example, it predicts that industrial demand for Q1 water will be X cubic meters per hour, while agricultural demand for Q2 water will be Y cubic meters per hour within the next 24 hours. The second is a cross-plant collaborative production plan, which includes specific scheduling instructions for each water plant. For example, Plant A should produce Q1-grade water at what power and during what period, how Plant B and Plant C should cooperate, and even when Plant C should execute water quality level switching, so as to meet the global forecast needs with the highest efficiency and lowest cost.

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

[0030] Exploring and identifying predetermined hydraulic delivery paths between pairs of nodes in the pipe network based on the dynamic pipe network topology in the standardized spatiotemporal dataset.

[0031] Further, it can be based on the dynamic pipe network topology in the standardized spatiotemporal dataset to build a time-varying graph structure.

[0032] Specifically, the system will generate a directed graph representation G(t) for each prediction time t based on the sequence of dynamic pipe network topologies. In this graph, water plants, large users, and pipe network junctions are abstracted as graph nodes and are assigned node type attributes (e.g., production nodes, consumption nodes, transfer nodes) and state vectors containing real-time water quality parameters, production capacity, or demand. Pipe segments are abstracted as directed edges, and the weight of the edge can take into account physical length, pipe diameter, wall roughness degradation coefficient based on pipe age and material, and other hydraulic parameters. When the standardized spatiotemporal dataset contains pipe network construction information, such as a pipe segment being closed at time t, the algorithm will automatically remove the edge and may activate a pre-set temporary alternative path, and the delivery capacity of the alternative path is calculated as a certain percentage (e.g., 60%) of the original path.

[0033] Exploring and identifying multiple hydraulic delivery paths between pairs of nodes in the pipe network based on the time-varying graph structure.

[0034] In this embodiment, in order to break the limitation of traditional models relying on a single shortest path and to more realistically reflect the water flow distribution in complex ring-shaped pipe networks, the system will perform an improved Yen algorithm for each source-target node pair (s, d) (e.g., a water plant s and an industrial user d): find a shortest path by Dijkstra algorithm (Dijkstra) or other methods, and by iteratively removing edges on the found path and re-planning detour paths, systematically search for the top K candidate paths, where K is a pre-set integer, for example, K ≤ 5, to balance the computational complexity and path coverage. Further, in order to ensure that the identified paths are physically feasible, the system will also perform hydraulic feasibility screening on the K candidate paths. Specifically, the total hydraulic slope of each path will be calculated, and those paths that will result in a pressure lower than the minimum service pressure requirement along the way will be removed. After screening, a set of all physically reachable effective paths between each node pair is obtained.

[0035] Constructing a path conduction kernel for each hydraulic delivery path, the path conduction kernel presents a bimodal distribution pattern to represent physical transmission delay and decision response delay, respectively.

[0036] wherein the path conduction kernel K p(τ) is the probability density function of the time delay τ, which is used to describe the probability distribution of the water quality signal of water starting from the source node at time t and arriving at the target node at time t+τ. In this embodiment, in order to capture both the physical process delay of water flowing in the pipeline and the response process delay of the water system scheduling decision, the bimodal gamma distribution is used to model the 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 is the weight coefficient of the two peaks, and w 1p +w 2p = 1. The first peak Γ(τ; k 1p ,θ 1p ) is used to characterize the physical transmission delay. Its peak position roughly corresponds to the average transit time of the water particle along the path p, and its shape parameter k 1p and scale parameter θ 1p According to the principle of fluid mechanics, the total length L of the path p can be p , average flow velocity v p and Reynolds number Re p For example, for turbulent flow, k 1p Can be set to 2.2 + 0.12*L p / v p ,θ 1p Can be set to 1.1 + 0.06*L p / v p The second peak Γ(τ;k 2p ,θ 2p ) characterizes decision-making response delay, reflecting the time required for a water plant to make the necessary human or system decisions, from receiving a dispatch instruction or market demand signal to adjusting production and opening valves for water delivery. This delay isn't caused by physical delivery itself, but it's common in real-world supply and demand responses and shouldn't be ignored.

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

[0038] The water quality grades of recycled water associated with each path are extracted from a standardized spatiotemporal dataset. Based on the water quality grades, the parameters of the peak representing the decision response delay in the bimodal distribution are determined.

[0039] In this embodiment, the parameter k of the decision response delay peak 2p and θ 2p are not fixed values, but are dynamically related to the water quality level of the delivered reclaimed water. Different water quality levels of reclaimed water have different requirements for the urgency and reliability of the supply by the users, resulting in different decision response modes of the water system. Specifically, the system extracts the real-time water quality level Q of the upstream node of the path p from the standardized spatiotemporal data set. According to the preset mapping relationship, the parameters are determined. For example, if the water quality level is Q1 level (for example, COD < 30 mg / L, used for high-standard industrial users), the supply interruption or delay will cause great loss, so the decision response requirement is fast. At this time, k 2p = 4, θ 2p = 1.5, forming a response distribution with a peak in front and a steep shape. If the water quality level is Q3 level (for example, COD ≥ 50 mg / L, used for municipal greening), the supply has a certain buffer space, and the decision response can be relatively more gentle. At this time, k 2p = 8, θ 2p = 2.5, forming a response distribution with a peak in the back and a flat shape. Alternatively, the parameters can also consider the daily periodic fluctuations, for example, k 2p can be a function containing a cos(2πt / 24) term to reflect the difference in scheduling response speed at different times of the day. In this way, the transmission kernel model incorporates the actual experience of water operation, making it more realistic and accurate.

[0040] Alternatively, constructing a dynamic graph network prediction model further comprises: based on the pipe network hydraulic parameters in the standardized spatiotemporal data set, analyzing the flow distribution ratio of the predetermined hydraulic transmission path; according to the flow distribution ratio, weighting and superimposing the path transmission kernel of each hydraulic transmission path to form a comprehensive transmission kernel representing the pipe network node pair.

[0041] In this embodiment, after identifying the multiple effective paths (such as paths p1, p2,...) between the node pair (s, d) and constructing a dedicated transmission kernel K p1 (τ), K p2 (τ),... for each path, it is necessary to combine them into a single comprehensive transmission kernel K sd (τ) representing the overall transmission characteristics between (s, d). To this end, it is necessary to determine how the total flow is distributed among the paths, which can be achieved by solving a nonlinear hydraulic balance equation set based on the node flow conservation law and the loop energy conservation law. Specifically, the Newton-Raphson iteration method can be used to solve the equation set, thereby obtaining the flow Q on each path p.p Accordingly, the flow distribution ratio a p (t) can be calculated as Q p (t) / ∑Q k (t). According to the calculated flow distribution ratio, the path conduction kernel K p (t) of each path is weighted and superimposed to obtain the comprehensive conduction kernel K sd (t, t) =∑ p a p (t) • K p (t). It can be seen that the constructed dynamic graph network model can accurately depict the water quality space-time conduction law in a complex pipe network considering the multi-path effect and decision delay through the comprehensive conduction kernel with clear physical meaning and self-adaptive adjustable parameters, thereby providing a high-quality model basis for subsequent supply-demand prediction.

[0042] The embodiment solves the problem of water quality space-time conduction characteristic distortion caused by the single-path and single-delay assumption. Specifically, by exploring and identifying multiple hydraulic transportation paths between nodes, the limitation of the traditional shortest path is broken, the distribution and convergence effect of water flow in a complex ring pipe network can be truly reflected, the model can still make accurate prediction under the topology change scene such as the closure or congestion of the main pipe, and strong robustness is shown. A conduction kernel with a bimodal distribution is constructed for each path. The first peak accurately represents the physical transmission delay based on the principle of fluid mechanics; the second peak quantifies the decision response delay, i.e., the human and system lag in operation and scheduling, which has been ignored before. The prediction result is highly consistent with the actual operation scene. Further, the model dynamically associates the parameters of the decision response delay peak with the quality level of the transported reclaimed water, gives a faster response mode for high-quality and high-demand water supply tasks, and embeds the operation experience and wisdom in the mathematical model. Through the weighted superposition of the path conduction kernel of each path according to the flow distribution ratio calculated based on the hydraulic parameters, a comprehensive conduction kernel with clear physical meaning and self-adaptive adjustable parameters is formed. The leap from single-time-point delay prediction to complete time delay probability distribution prediction is realized, the prediction accuracy of water quality arrival time and concentration distribution is improved, and a solid foundation is provided for realizing the timely matching of supply and demand and avoiding resource waste.

[0043] As an optional solution, as shown in Figure 4 , the dynamic graph network prediction model is also constructed, which comprises determining a time-lag conduction kernel, specifically:

[0044] Based on the hydraulic and water quality parameters in the standardized space-time data set, a one-dimensional convection-reaction-diffusion equation is established for the node pair of the pipe network; the analytical or numerical solution of the one-dimensional convection-reaction-diffusion equation is taken as the time-lag conduction kernel.

[0045] In particular, for each pipe segment in the pipe network, the concentration C of a pollutant (e.g. COD) can be considered to vary in both pipe length x and time t, following a one-dimensional convection-reaction-diffusion (CRD) equation: Vc / Vt + v-Vc / Vx = D eff •V 2 c / Vx 2 - k-c; where V is the partial derivative, C(x, t) is the pollutant concentration; v is the average velocity of water flow in the pipe (convection term); D eff is the effective diffusion / dispersion coefficient (diffusion term); k is the first-order reaction degradation coefficient (reaction term). When considering a pulse input from a source at x0at time t0(i.e. Dirac delta function), the Green's function solution K(x, t | x0, t0) of this equation, whose physical meaning is the concentration distribution formed at location x after time t-t0from a unit concentration pulse at the source. Therefore, this Green's function solution can be directly used as the time-delay conduction kernel connecting the source and target nodes. Its analytical solution is usually a Gaussian distribution or a modified skewed distribution, for example: K(x, t | x0, t0) = (1 / sqrt(4pD eff •r))•exp[-(x-x0-v•r 2 / (4D eff •r)]•exp[-k•r], where r = t-t0.

[0046] Further, the one-dimensional convection-reaction-diffusion equation is parameterized. In order to make the above CRD equation accurately reflect the complex situation of the real pipe network, it is necessary to finely calibrate the core parameters k, v and the effective pipe diameter involved in the equation itself, specifically, parameterization is performed according to a standardized spatiotemporal dataset, and the parameterization includes:

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

[0048] The traditional model usually takes the reaction degradation coefficient k as a constant, but this ignores the complexity of biochemical reactions. In the present embodiment, the reaction degradation coefficient k is modeled as a dynamic variable k(C, Q). For example, a modified Monod equation or Arrhenius equation can be used in combination: k(C, Q) = k0•[1 + β1•(C / C ref ) n ]•[1 + β2•exp(-E a / (R•T))] ; where k0is the basic degradation rate; C is the real-time pollutant concentration; C refis the reference concentration; n is the reaction order; β1 is the concentration effect coefficient; T is the water temperature; E a is the reaction activation energy; R is the gas constant; β2 is the temperature influence coefficient. Alternatively, k0 can also be a piecewise function related to the water quality grade Q. This modeling approach enables the reaction rate to be self-adaptively adjusted with changes in pollutant concentration, water temperature, and overall water quality environment.

[0049] A dynamic biofilm growth model for pipe wall is 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, biofilm will grow on the inner wall of the pipe, which will reduce the effective water passing section of the pipe and increase the wall roughness, thereby affecting the flow rate v and head loss, which is generally ignored by traditional models. In this embodiment, a differential equation describing the evolution of biofilm thickness Δ with time can be established: dΔ / dt = μ max •(S / (K s +S))•Δ - b•Δ - τ w •Δ 2 / ρ b ; where μ max is the maximum growth rate; S is the substrate concentration (such as COD); K s is the half-saturation constant; b is the decay rate; τ w is the shear force of water flow on the pipe wall; ρ b is the biofilm density. By numerically solving the 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 is the original inner diameter of the pipe; and the hydraulic roughness can also be dynamically corrected according to Δ(t). These corrected parameters will be used to more accurately solve the CRD equation.

[0051] According to the water quality grade, a differentiated propagation speed is set for the one-dimensional convection-reaction-diffusion equation, where 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 speed mainly refers to the convection term speed v in the CRD equation. Different quality of reclaimed water has different hydrodynamic characteristics. For example: Q1 high-quality water, close to clean water, its flow rate v Q1 is basically consistent with the theoretical hydraulics calculation value, which can be set as v theoretical •0.98, where v theoretical is the theoretical hydraulics flow rate. Q2 medium water quality contains more colloidal particles, showing non-Newtonian fluid characteristics, its apparent viscosity increases, resulting in a decrease in actual flow rate vQ2 Can be set to v theoretical •0.92. Q3 low-quality water may contain suspended solids, produce Bingham plastic flow characteristics, and have yield stress, which further reduces the flow rate. Q3 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. This example provides a more refined method for constructing conductive cores based on physical mechanisms, particularly suitable for scenarios requiring high simulation accuracy of water quality degradation processes.

[0053] This embodiment addresses the issues of roughness and inaccuracy in simulating water quality evolution, 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 deeply reveal the underlying mechanisms of pollutant migration and transformation in pipelines. Its core advantage lies in the deep dynamicization and domain-specific customization of model parameters: the reaction degradation coefficient k is modeled as a function of pollutant concentration and water quality grade, eliminating the traditional approach of treating it as a constant. This allows the simulation of reaction rates to more closely resemble real biochemical processes. A dynamic growth model for biofilms on pipe walls is established, enabling real-time calculation of biofilm thickness, which is used to dynamically correct effective pipe diameter and hydraulic roughness. This overcomes the significant drawback of traditional models, which fail to reflect the long-term impact of pipeline aging on hydraulics and water quality. Differentiated propagation speeds are set based on water quality grade, capturing the fast-lane and slow-lane effects caused by differences in fluid properties of recycled water of varying quality. The model can not only predict when it will arrive, but also accurately predict what the state will be when it arrives, 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 the present application, the parameter determination method and influencing mechanism of the dynamic growth model of the pipe wall biofilm can be further as follows:

[0055] Key parameters of the dynamic growth model of biofilm on pipe wall, such as the maximum growth rate μ max , half-saturation constant K s The parameters are not set in vain, but can be calibrated through laboratory coupon tests or online monitoring devices. Specifically, coupons made of standard materials can be placed in a test loop that simulates the actual pipeline environment, and different concentrations (corresponding to the matrix concentration S) and flow rates (corresponding to the wall shear force τ w By regularly taking out the coupons and measuring their biofilm thickness Δ and dry weight, a set of Δ changes with time t, S and τ can be obtained. wExperimental data. The experimental data were fitted to the differential equation model using a nonlinear least squares method to identify the values of the key parameters μ max , etc. In addition, the effect of biofilm on conduction can be introduced as an additional conduction delay factor caused by the adsorption-release effect of biofilm. Because biofilm, like a sponge, temporarily adsorbs pollutants in the water flow and slowly releases them later, this is macroscopically manifested as an additional delay and a longer tailing effect. This delay factor can be modeled as a function related to the biofilm thickness Δ, adsorption-release rate, pollutant partition coefficient K d , etc., and is convolved with the original conduction kernel, thus more realistically simulating the conduction process.

[0056] According to another aspect of the present application, the propagation speed is set differently according to the water quality level, and the physical mechanism behind this can further be the stratified flow of water quality on the pipe cross-section. Specifically, because the flow velocity is highest in the turbulent core region and lowest in the near-wall region, and high-quality reclaimed water (low in impurities and low in density) tends to concentrate in the central region of the pipe, and low-quality water (high in impurities and high in 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 radius of the pipe, C center is the concentration at the center of the pipe, C wall is the concentration at the wall of the pipe, and n is the distribution index (radial concentration change rate). Therefore, the overall equivalent propagation speed of the high-quality water mass will be higher than that of the low-quality water mass, forming a fast-slow lane effect.

[0057] According to an aspect of the present application, the dynamic graph network prediction model is further subjected to a sparse processing to generate a sparse connection matrix, and the generation of the product quality demand prediction result and the cross-plant collaborative production scheme is performed using the sparse connection matrix.

[0058] In the present embodiment, the sparse processing can solve the dimension disaster problem. A complete urban reclaimed water pipe network graph model can contain thousands of nodes, and if all possible connections between the nodes are considered, the calculation complexity will be huge, making it difficult to meet the real-time prediction requirements. Therefore, an intelligent sparse strategy is needed to identify and only retain the key connections in the network that play a decisive role in hydraulic conduction and water quality evolution, while ignoring the secondary connections, thereby generating a sparse connection matrix. Most of the elements in this matrix are zero, and only the elements corresponding to the key connections are non-zero values. All subsequent model operations will be based on this sparse matrix, thereby reducing the calculation complexity from O(N 2) to near O(N) level, where N is the number of nodes.

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

[0060] The first method is an adaptive sparsification method based on system state. The sparsification process includes: quantitatively evaluating the real-time system complexity indicator C(t) of the pipe network according to the topology change rate, demand fluctuation rate and water quality heterogeneity in the standardized spatio-temporal data set; based on the real-time system complexity indicator, the connection retention ratio ρ(t) is determined through the function relationship of ρ(t)=ρ base ×exp(-β×C(t)), and a sparse connection matrix is generated, where ρ base is the basic retention ratio, and β is the sensitivity parameter. In other words, based on the real-time system complexity indicator, the connection retention ratio used to generate the sparse connection matrix is dynamically adjusted.

[0061] Specifically, the sparsity of the model can be dynamically adjusted according to the real-time busy and stable state of the pipe network system. To this end, a system complexity indicator C(t) is needed to quantify the real-time state of the pipe network. The indicator 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) is the topology change rate, which measures the intensity of the pipe network structure change, for example, it can be obtained by calculating the ratio of the number of newly added or removed pipe segments at the current time to the total number of pipe segments; C demand (t) is the demand fluctuation rate, which measures the intensity of user demand change, for example, it can be calculated as the coefficient of variation (standard deviation divided by mean) of the time series of all user demand totals; C quality (t) is the water quality heterogeneity, which measures the complexity of the water quality distribution in the network, for example, it can be obtained by calculating the Shannon entropy of the volume proportion of different water quality grades (Q1, Q2, Q3) in the pipe network; λ1, λ2, λ3 are the weight coefficients of each term, which can be set according to experience, for example, λ1=0.4, λ2=0.3, λ3=0.3.

[0062] After obtaining the real-time system complexity indicator C(t), the system will dynamically adjust the connection retention ratio, i.e. the sparsity. Specifically, the dynamic connection retention ratio ρ(t) can be determined by an exponential decay function: ρ(t)=ρ base ×exp(-β×C(t)); where ρ baseis the base reserve ratio, which can be set as 0.35 (i.e., 35% of the connections are reserved by default); β is the sensitivity parameter, which is used to control the response speed of sparsity to the complexity change, and can be set as 0.5. It can be seen that when the pipe network system is in a stable state (e.g., at night, no construction, and stable demand), the value of C(t) is small, and ρ(t) is close to ρ base , the system will adopt a higher sparsity (i.e., reserve fewer connections) to save computing resources. When the pipe network system is in a state of severe change (e.g., a heavy rain causes a surge in water demand, accompanied by water quality changes caused by pipe flushing), the value of C(t) will increase significantly, and ρ(t) will decrease exponentially, and the system will automatically reduce the sparsity and reserve more connections to ensure the prediction accuracy of the model at this critical moment.

[0063] The second method is a protective sparsification method based on domain knowledge and aiming to ensure water quality safety. The sparsification process includes: based on the water quality data in the standardized spatiotemporal dataset, quantitatively evaluating the water quality incompatibility risk when the nodes connected by the network are mixed; cutting off or weakening the network connections with water quality incompatibility risk exceeding a preset mixing risk threshold, to generate a sparse connection matrix.

[0064] In this embodiment, it is found that some connections in the pipe network are feasible in terms of water power, but must be avoided from the perspective of water quality safety. If sparsification is not distinguished, these risky connections may be reserved. Specifically, the system quantitatively evaluates the potential water quality incompatibility risk of each connection (pipe section). This risk is multidimensional, for example: chemical incompatibility risk: some water of different sources mixed together will have an adverse chemical reaction. For example, when a pipe is planned to transport water with a residual chlorine content greater than 1.0 mg / L, and the node area to which it goes currently has an ammonia nitrogen content greater than 2.0 mg / L, the mixture will generate chloramine with poor disinfection effect and possible odor. The system can calculate a reaction potential Φ react according to the chemical reaction kinetics model, and when it exceeds the threshold, it is considered that the connection has a high chemical incompatibility risk. Quality degradation risk: high-quality water flows into a pipe contaminated by low-quality water, which will cause irreversible quality degradation. The system can define a degradation risk function R degrade , the value of which is related to the quality level difference, the importance of downstream users, and other factors. Microbial cross-contamination risk: especially when the water supply goes to sensitive users such as hospitals and food factories, it is necessary to strictly avoid backflow or mixing of water from areas with low health standards. A microbial cross-contamination risk indicator R microbialto evaluate this risk. After quantifying the incompatibility risk of each connection, the system can calculate a comprehensive risk score. If the score exceeds a pre-set mixing risk threshold, the connection is subjected to a cut-off (setting its corresponding element in the sparse matrix to 0) or a weakening (keeping the connection but multiplying it by a penalty factor much smaller than 1 in its conductance kernel or weight) treatment. Optionally, it can also be designed as a tiered protection mechanism. For example, for pipes supplying first-class protection users (such as hospitals), all their upstream connections are subjected to the most stringent risk assessment; while for pipes supplying third-class users (such as landscaping), a more lenient threshold can be adopted. Not only does this improve computational efficiency, but also embeds a safety valve at the model structure level, so that any scheduling scheme it generates will not cause water quality safety accidents, which is crucial for recycled water, a quality-sensitive resource.

[0065] The embodiment successfully solves the sharp contradiction between computational efficiency and operational safety of large-scale pipe network models. Specifically, to overcome the dimension disaster brought by the full connection graph model, the model is sparsified, reducing the computational complexity from O(N 2 ) level to nearly O(N) level, making real-time and rolling prediction and optimization of urban-level pipe networks possible in engineering. A dual adaptive mechanism is also provided: one is dynamic sparsification based on the pipe network system complexity index C(t). This index integrates topological changes, demand fluctuations, and water quality heterogeneity, and can quantitatively measure the busy degree of the pipe network in real time. Accordingly, the model can adopt a higher sparsity when the system is stable (such as at night) to save computing power, and automatically retain more connections when the system is in a state of flux (such as during heavy rain) to ensure accuracy, achieving optimal allocation of computing resources. The second is protective sparsification based on water quality incompatibility risk. This mechanism embeds expert knowledge in the water safety field into the algorithm, quantitatively evaluates the chemical reaction risk or quality degradation risk of different water quality mixtures (such as mixing of high residual chlorine water and high ammonia nitrogen water), and actively cuts off or weakens network connections that exceed the risk threshold. This is equivalent to embedding a safety valve at the model structure level, preventing the optimization algorithm from generating scheduling schemes with water quality safety hazards in pursuit of economy, and ensuring the safety of recycled water, a sensitive resource, in complex coordinated scheduling.

[0066] According to one aspect of the present application, the quantitative evaluation of water quality incompatibility risk can further include constructing a specific risk evaluation matrix and threshold. Specifically, a chemical incompatibility matrix is established to evaluate the chemical reaction risk of different water quality index combinations. For example, for the problem of chloramine production caused by mixing of high residual chlorine water and high ammonia nitrogen water, the system will calculate a reaction potential Φ react = k r •[Cl]•[NH3]•exp(-E a / RT), where k rk is the reaction rate constant, [CI] and [NH3] are the concentrations of residual chlorine and ammonia nitrogen, respectively. The system will preset a critical reaction potential threshold value Φ critical When the calculated Φ react exceeds the threshold value, the mixing event is marked as high risk. Exemplarily, a more general water quality mixing compatibility matrix can be given as follows to guide the sparsity decision: when the mixing type is Q1 grade water and Q1 grade water, the risk score (0-1) is 0.05, and the decision suggestion is to allow mixing; when the mixing type is Q1 grade water and Q1 grade water, the risk score (0-1) is 0.05, and the decision suggestion is to mix cautiously, and the downstream user sensitivity needs to be evaluated; when the mixing type is Q1 grade water and Q3 grade water, the risk score (0-1) is 0.95, and the decision suggestion is to prohibit mixing in principle, triggering protective sparsity; when the mixing type is Q2 grade water and Q3 grade water, the risk score (0-1) is 0.60, and the decision suggestion is not to recommend mixing, and only consider when there is no alternative. In the sparsity processing, the system compares the calculated or tabled risk score with the preset mixing risk threshold value (for example, 0.7). Any connected risk score exceeding the threshold value will be preferentially cut off or weakened.

[0067] In an optional embodiment, a cross-factory collaborative production scheme is constructed on the water quality cascade utilization network, in which the effluent of the upper-level user is scheduled as the input water source of the lower-level user.

[0068] In other words, based on the preconfigured water quality cascade utilization network, the effluent of the upper-level user is scheduled as the input water source of the lower-level user in the cross-factory collaborative production scheme.

[0069] In this embodiment, in order to maximize the utilization value of the reclaimed water, the generation of the production scheme is not based on the traditional independent water supply mode of each user, but is constructed on the network topology of water quality cascade utilization. Because the effluent after use by the user with high quality requirements (such as industrial cooling, first level) still meets the needs of the user with lower quality requirements (such as urban miscellaneous use, second level) although the quality has decreased. Specifically, the system will identify and construct the optimal cascade utilization chain through the maximum weight matching algorithm (such as the Hungarian algorithm) according to the water demand quality of each user and the effluent quality after use, as well as their geographical location and pipe network connection relationship. For example, the optimized cascade utilization topology A factory (Q1) -> industrial park X (after use, Q2) -> municipal green land Y (after use, Q3) -> river discharge can be formed. Subsequent optimization solutions will be carried out under the constraint of this cascade utilization topology, so that the scheduling scheme can maximize the value of water resources.

[0070] Further, the generation of the cross-factory collaborative production scheme is realized by solving a large-scale multi-objective optimization problem, and the solving process includes: decomposing the large-scale multi-objective optimization problem into a plurality of low-dimensional sub-problems along the dimensions of water quality levels, time windows and spatial regions; and coordinating the solutions of the low-dimensional sub-problems through a distributed coordination protocol to obtain a globally optimized 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 (production and distribution paths of each water plant, each quality, and each time period) of the entire urban pipe network within the next 24 hours will form a mixed integer programming problem with extremely high dimensions (for example, up to tens of thousands of dimensions), which is difficult to solve directly. To solve this problem, the embodiment adopts a block 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 only have coupling relationships at water plants with switchable production capacity. Second-level decomposition - by time window: each water quality sub-problem is further divided into a plurality of shorter windows along the time axis, for example, six 4-hour windows. Detailed optimization is performed within each window, and the windows are coupled through boundary conditions (such as the final inventory equaling the initial inventory of the next period). Third-level decomposition - by spatial region: based on the community structure in the sparse connection matrix, the pipe network within each time window is divided into a plurality of spatial clusters using a graph partitioning algorithm (such as METIS). The clusters are tightly connected within the cluster and are collaboratively optimized; the clusters are sparsely connected between clusters and exchange a small amount of boundary information. Through these three levels of decomposition, the original high-dimensional problem is decomposed into hundreds of minimum sub-problems with lower dimensions (for example, about 120 dimensions). These sub-problems can be distributed to multiple computing cores for parallel solving (for example, using the interior point method). After all the sub-problems are solved, the system coordinates the solutions of the sub-problems through a distributed coordination protocol (for example, a consensus-based distributed optimization protocol or the alternating direction method of multipliers ADMM). Specifically, each sub-problem iteratively modifies its solution based on the solutions of its neighbor sub-problems (i.e., boundary information) and penalizes inconsistency by introducing a Lagrange multiplier. After several iterations, the solutions of all sub-problems converge to a globally optimal solution that satisfies the global constraints.

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

[0073] In the scenario of collaborative scheduling involving multiple water plants and multiple management departments, trust and data transparency are key. To address this issue, the present embodiment introduces blockchain technology. Specifically, the coordination rules, the production capacity commitments of each water plant, the demand orders of each user, and the final generated scheduling scheme can all be encoded into smart contracts deployed on a consortium blockchain. For example, a smart contract can stipulate that when water plant A registers that it can provide 1000 cubic meters of Q1 grade water in the next hour, the system will automatically match the demand according to the 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 the online water quality monitoring, the contract will automatically execute the settlement and record the performance of this time. To ensure that all participants (water plants, scheduling centers, etc.) agree on the scheduling scheme and it is not tamperable, the system can use a consensus algorithm optimized for water network characteristics, such as a delay-aware Byzantine fault tolerance algorithm (PBFT-W). This algorithm takes into account the physical delay of water transmission in the pipe network when reaching consensus, so that the consensus result is executable in reality.

[0074] Optionally, to encourage water plants to provide stable and high-quality reclaimed water, a water quality reputation proof (Proof of Quality, PoQ) mechanism can be constructed. Under this mechanism, water plants with a better record of water quality compliance rate and stability in the past will have a higher weight or priority in the consensus process, thus forming a virtuous cycle. Not only is distributed collaboration achieved efficiently, but also the transparency, traceability and mutual trust between participants of the entire collaborative scheduling system are enhanced.

[0075] The embodiment solves the problems of decision dimension limitation and low efficiency in multi-source and multi-quality collaborative scheduling. Specifically, a water quality cascade utilization network is constructed, the effluent of a high-quality user is used as the input water source of the next-level user, and an optimal cascade utilization chain is found through dynamic programming algorithm, so that the goal of maximizing resource value is established in the top-level design of the optimization model, and the low-efficiency mode of independent water supply of each user is changed. In order to solve the super-large multi-objective optimization problem generated thereby, a decomposition strategy along the 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, which is difficult to handle, into hundreds of low-dimensional sub-problems that can be solved in parallel. After solving the sub-problems, a consistent-based distributed coordination protocol (such as ADMM) is used for overall coordination, so that the global optimal solution is finally converged. The non-deterministic polynomial difficulty (NP-hard) problem which is difficult to solve in theory is transformed into a feasible calculation task in engineering. Alternatively, by introducing a blockchain smart contract to execute the coordination protocol, the transparency and mutual trust of multi-agent (such as different water plants) collaboration are further enhanced. A truly global optimal cross-plant collaborative production scheme is possible, which can reduce the total operating cost of the system, avoid high-quality low-use, 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 can further include:

[0077] A dynamic programming algorithm is used 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. The state value function V(i, q) is defined as the maximum economic or environmental benefit that can be formed by the subsequent cascade utilization chain starting from node i with water of not higher than water quality level q as input. Its state transition equation can be expressed as: V(i, q) = max{V(j, q') + R(i, j, q, q')}, where j is a downstream node of i, q' is the water quality level converted after utilization of node i and q' < q, and R(i, j, q, q') is the revenue generated by the single-step cascade utilization from node i to j. By solving this dynamic programming problem in reverse, the globally optimal longest cascade utilization chain can be obtained.

[0079] A cascade matching matrix is constructed, and the matrix elements of the cascade matching matrix comprehensively consider the demand matching degree, quality matching degree and distance factor between users.

[0080] Alternatively, a cascade matching matrix M is constructed for evaluating and screening the feasibility of cascade connection cascade . The elements M of the matrix cascade(i, j) represents the matching degree of the effluent of the upper user i as the input source of the lower user j, whose value can be determined by considering the following three dimensions: M cascade (i, j) = Demand match (i, j) • Quality match (i, j) •Distance factor (i, j) ; where Demand match is the ratio of the effluent quantity of user i to the demand quantity of user j; Quality match is an exponential function whose value decreases as the gap between the effluent quality of user i and the demand quality of user j increases; Distance factor is a decay function related to the distance between the two users. Only when M cascade (i, j) exceeds a pre-set threshold (e.g. 0.6), the cascade connection is considered feasible.

[0081] A quality conversion function AQ use is established for each process to quantify the water quality degradation and potential energy recovery in each level of utilization. This function is specific to the water use process. For example: for industrial cooling (non-contact), the water quality degradation model is AQ cooling = 0.05 • Q in • (1 + 0.1 • AT), mainly considering the impact of temperature rise AT, where Q in is the input water quantity or input water quality level. For agricultural irrigation, the water quality degradation model is AQ irrigation = 0.3 • Q in • (1 + 0.2 • ET / P), additionally considering the concentration effect of evapotranspiration ET on pollutants, where P is the precipitation. At the same time, the system also establishes an energy recovery model, for example, using the water head difference between cascades to generate electricity through micro-turbines, or using the temperature difference to recover heat energy through heat pumps, and the recovered energy can offset about 15-20% of the operating cost.

[0082] According to an aspect of the present application, the process of solving large-scale multi-objective optimization problems by block-wise recursive decomposition can further include the following detailed algorithm implementation: a three-level decomposition strategy is performed to decompose the original high-dimensional problem into multiple low-dimensional, parallelizable sub-problems. First-level decomposition (according to water quality level): the problem is decomposed into three sub-problems Q*1, Q*2, Q*3, which only exist coupling constraints at nodes with switchable capacity. Second-level decomposition (according to time window): the 24-hour scheduling period is divided into 6 four-hour windows, and the windows are coupled through boundary variables (such as end-of-period inventory). Third-level decomposition (according to spatial region): based on the community structure of the sparse connection matrix, the pipe network within each time window is divided into about 10 spatial clusters using graph partitioning algorithms such as METIS, with tight coupling within the cluster and loose connection between clusters. On this basis, the numerous sub-problems after decomposition are parallelized and distributedly coordinated. Specifically, the system can distribute hundreds of minimum sub-problems to multiple computing cores, and each core independently solves using an efficient algorithm such as the interior point method. After all sub-problems complete a round of solving, a distributed coordination protocol based on the consensus ADMM algorithm is started, and its iteration rule can be: local variable update: each sub-problem p solves x p (k+1) = argmin(f p (x) + ρ / 2•||x - z p k || 2 ), where x p (k+1) is the local variable of sub-problem p (k+1 round), f p (x) is the objective function of sub-problem p, ρ is the penalty strength of the consistency constraint, x is the variable used for optimization in the current iteration, z p k is the global consistent variable value referred to by sub-problem p in the kth iteration. Global variable update: exchange information with neighbor sub-problems through boundary nodes and take the average z (k+1) = average(x q (k+1) ), where z (k+1) is the global consistent solution obtained by all sub-problems after the k+1th iteration through boundary node information exchange and average operation, and x q (k+1) is the local solution of neighbor sub-problem q that has a coupling relationship with sub-problem p. Dual variable update: λ p (k+1) = λ p k + ρ(x p (k+1) - z (k+1) ), where λ p (k+1)is the dual variable of the subproblem p in the k+1th iteration. The iteration process continues until the solutions of all subproblems converge to the global optimal solution. To ensure the convergence of the iteration algorithm, it can be theoretically proved that its update operator is a contraction mapping. Its convergence is based on the Banach 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 used to use the solutions of similar historical scenarios as initial points to speed up convergence.

[0083] In an optional embodiment, further comprising: based on the water quality time series data in the dynamic graph network prediction model, using chaos theory analysis or infectious disease dynamics model to generate a water quality evolution risk atlas containing future potential risks.

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

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

[0086] The first analysis method: water quality mutation prediction based on chaos theory. Used to identify the precursors of water quality mutation at a single key node. Specifically, the system can extract the water quality parameter (such as COD) time series data of a certain key monitoring node (for example, a node upstream of an important user) in the past period (for example, 7 days) from the dynamic graph network model. Use Takens embedding theorem to reconstruct the phase space of the one-dimensional time series and convert it to a trajectory in high-dimensional space, where the optimal embedding dimension m and delay time τ can be determined by the false nearest neighbor method and mutual information method. Calculate the largest Lyapunov exponent λ max of the trajectory. If λ max > 0, it indicates that the water quality dynamics system has chaotic characteristics, i.e. a small initial disturbance will be exponentially amplified, indicating that the system has inherent unpredictability and mutation possibility. The upper limit of the prediction time scale is about 1 / λ max . When the system approaches the critical point (i.e. is about to mutate), a critical slowing down phenomenon will occur. The system can monitor this phenomenon by continuously calculating the autocorrelation coefficient and variance of the time series within a sliding time window. When the autocorrelation coefficient significantly increases and the variance sharply increases, it can be determined that a mutation is about to occur, thereby issuing an early warning, for example, the occurrence of a mutation can be predicted 2-4 hours in advance.

[0087] The second analysis method is water quality deterioration propagation modeling based on an infectious disease dynamics model. This method is used to predict the possible propagation range and speed of water quality deterioration in the entire pipe network once a water quality deterioration occurs at a certain node (mutation predicted by the above method). Specifically, the water quality deterioration process can be analogized to the spread of an infectious disease in a population. The system can divide the nodes in the pipe network into three categories: S (Susceptible), i.e., nodes with normal water quality; I (Infected), i.e., nodes with deteriorated water quality; and R (Recovered), i.e., nodes with improved water quality through treatment. On this basis, an improved SIR (Susceptible-Infected-Recovered) model is established. Unlike the standard SIR model, the infection rate β here is not a constant but a variable related to the topology of the pipe network. For example, β ij (the infection rate of node i to node j) can be related to the degrees (i.e., the number of connected pipes) of nodes i and j and the hydraulic conduction efficiency between them, meaning that a hydraulic hub node has stronger propagation ability. Similarly, the recovery rate α can be related to the treatment capacity of the downstream water plant. Through this model, the basic reproduction number R0 of the pipe network can be calculated. If R0 > 1, it means that once a single contaminated point occurs, water quality deterioration will spread widely in the network; if R0 < 1, the deterioration will be limited to a local range and eventually subside. Through simulation of the model, the probability P infect (i, t) that each node will be infected (i.e., water quality deterioration) in the future and the average infection time T infect (i) can be obtained.

[0088] In summary, by integrating the results of the above two methods, the system can generate a water quality evolution risk map that clearly indicates which nodes have mutation risks in the future, the most likely propagation path once a mutation occurs, and the impact degree and time on downstream nodes.

[0089] The water quality evolution risk map is converted into risk constraints, and a multi-objective optimization model is established to generate a cross-plant collaborative production scheme under the premise of meeting the risk constraints.

[0090] In this embodiment, the risk map is not only used for early warning but also directly integrated into the multi-objective optimization model as a set of hard or soft risk constraints. Specifically, this conversion can take the following forms: probability constraint (chance constraint): for a high-risk node j, a constraint is added to the optimization model, i.e., P(Q j (t) ≥ Q min ) ≥ 1 -α. Where P(·) represents probability, Q j (t) is the water quality of node j at time t, and Q minis the minimum water quality standard, and a is a very small risk tolerance level (e.g. 0.05). This constraint means that the optimization solution must guarantee that the probability of water quality at node j not meeting the standard is no more than 5%. This probabilistic constraint can be transformed into a deterministic equivalent form that can be converted into a linear or nonlinear constraint that a standard optimization model can handle. Conditional Value at Risk (CVaR) constraint: This is a more robust risk measure. The optimization model can be asked to minimize the expected operation cost while constraining the CVaR value of the cost (i.e. the average cost in the worst 5% scenarios) to be no more than a pre-set upper limit L CVaR . This forces the optimizer to consider not only the expected situation but also the extreme high-cost events that can happen. Resilience constraint: A system-wide performance indicator System Performance (t) can be defined and added as a constraint System Performance (t) ≥ P min •(1 - D(t)•(1 - R(t))) ; where P min is the minimum acceptable performance level, D(t) is the disturbance intensity predicted by the risk map, and R(t) is the recovery ability of the system. This constraint makes sure that even when the disturbance predicted by the risk map happens, the system performance will not drop below the pre-set bottom line. By adding these risk-map-derived constraints to the optimization model, the model will be forced to avoid those solutions that are the cheapest in normal situations but are vulnerable in risk scenarios when the model is solved. For example, the optimizer can choose to activate a backup water source or select a longer but safer water delivery path, which will slightly increase the normal operation cost but buy insurance for the system against unknown risks, resulting in a production plan that is both economic and robust.

[0091] According to an aspect of the present application, the process of utilizing the water quality evolution risk map can further include quantitatively evaluating the resilience of the pipe network. Specifically, after generating the risk map, the system can calculate multi-dimensional resilience indicators, such as: resistance R resist : defined as the inverse of the drop in a key performance indicator (e.g. water supply assurance rate) when the system is subjected to a disturbance. recovery R recover : defined as the inverse of the time needed for the system performance to recover from the lowest point to the normal level. By performing deliberate attack simulation on the network (e.g. removing key nodes or edges in the simulation) and observing the changes in these resilience indicators, the vulnerable points and resilience bottlenecks of the system can be identified, providing quantitative basis for reinforcement of the pipe network and development of emergency plans.

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

[0093] comparing the prediction result of the graph network prediction model with the actual operation data to obtain a prediction error.

[0094] In this embodiment, the starting point of model revision is to obtain the feedback signal, i.e., the gap between prediction and reality. After the end of a prediction and scheduling cycle (e.g., 72 hours), the system automatically collects the actual operation data of the pipe network in this cycle. These data include: the actual water consumption of users in each area, the real water quality detection results uploaded on the online monitoring devices deployed at each node of the pipe network, and user satisfaction evaluation, etc. The system compares these actual data with the prediction results generated by the model in the previous generation 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, conduction delay prediction error). Further, 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, usually caused by defects in the model structure itself; random error, caused by unpredictable random disturbances in the system; conduction error, caused by the estimation bias of the model for the time-lag conduction kernel. The revision mechanism of this embodiment will mainly focus on revising the model parameters related to the conduction error.

[0095] Using the recursive least squares algorithm, the parameters of the dynamic graph network prediction model are continuously revised based on the prediction error.

[0096] In this embodiment, the model parameters, especially the conduction kernel parameters, are updated smoothly and continuously online. The recursive least squares (RLS) algorithm can update the estimated value of the parameters at 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 value at time t (e.g., the true conduction delay obtained through tracer experiments or data association analysis); φ(t) is the feature vector corresponding to the observation (e.g., containing path length, flow, water quality level, etc.); e(t) is the observation noise, T is the transpose. The parameter update formula of the RLS algorithm is: θ(t+1) = θ(t) + K(t)[y(t) - φ(t) T θ(t)]; where θ(t+1) and θ(t) are the estimated values of the parameter vector before and after updating, respectively. K(t) is a key adaptive gain matrix that determines the step size and direction of this update. Further, the calculation of K(t) introduces a variable forgetting factor λ(t). The calculation formula of λ(t) is: λ(t) = λ min + (λ max - λ min)•exp(-|e(t)| / e ref ); where e(t) is the current prediction error, λ min and λ max are the upper and lower bounds of the forgetting factor (e.g., 0.95 and 0.995), e ref is the reference error. When the prediction error e(t) is large, λ(t) becomes smaller, and past data is forgotten more quickly, allowing the system to quickly adapt to new changes. When the prediction error is small, λ(t) becomes larger, and more emphasis is placed on utilizing historical information, making parameter estimation more stable.

[0097] The cumulative sum algorithm is simultaneously used to monitor whether there are parameter mutations, so as to trigger accelerated correction when a mutation is 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 pipeline due to scaling over time). However, for sudden changes in the system (e.g., a section of old pipeline is completely replaced with a new material pipeline, resulting in a step change in its conductivity characteristics), the RLS may not respond quickly enough. To address this issue, a cumulative sum (CUSUM) algorithm is deployed in parallel to detect sudden changes in parameters. The CUSUM algorithm continuously monitors the parameter sequence θ(t) updated by the RLS and calculates two cumulative sum statistics S + (t) and S - (t). These two statistics are used to detect whether the parameter has drifted significantly upward or downward. + (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 the parameter, and k is the tolerance. + (t) or S - (t) exceeds the preset decision threshold h, the system determines that the parameter has mutated. After detecting the mutation, the system immediately triggers the accelerated correction mechanism. This mechanism temporarily adjusts the behavior of the RLS algorithm, for example, forcing its forgetting factor λ to a very small value (such as λ fast= 0.9), and reset its error covariance matrix. This makes the RLS algorithm learn the new system characteristics at a very fast speed in the next few iterations, thus achieving a fast response to the abrupt event. Through the synergy of RLS and CUSUM, a robust and efficient online model revision loop is constructed, which can adapt to slow drifts and quickly respond to abrupt changes, ensuring the accuracy and timeliness of the prediction model in the long-term operation.

[0099] In an alternative embodiment, the real-time demand signal is modeled with a time-lag constraint, generating a time-lag demand distribution matrix.

[0100] It is found that the production decision of the water plant must be made in advance of the actual water demand of the user, and the amount of time in advance is the delivery time delay of water in the pipe network. Therefore, it is incorrect to simply use the current demand signal of the user as the basis for production of the water plant. A mathematical model can be constructed to translate the future demand of the user into the current production task of the water plant. As a preferred implementation, the modeling process includes the following aspects: constructing a backward conduction kernel. In contrast to the forward conduction kernel K forward (τ), the backward conduction kernel K reverse (τ) aims to project the demand of user node j at time t back to the upstream water plant node i, so as to calculate how much production the water plant i should start at time t-τ to meet the demand. Specifically, the backward conduction kernel can be solved from the known forward conduction kernel by Bayesian inference. Alternatively, the modeling here considers the asymmetry of information propagation: the backward propagation speed v info of the demand signal (information) can be much faster than the forward propagation speed v flow of the water flow (material). For example, v info = 10 × v flow can be set. This means that when constructing the backward kernel, its time scale will be compressed accordingly, but a minimum response time (e.g. 0.5 hours) will be retained to consider the physical limitations of production start. In addition, when a user j is simultaneously supplied by multiple water plants (e.g. i and k), its demand needs to be inversely allocated proportionally. The system can calculate the allocation weight w according to the production capacity P of each water plant and the production capacity reachability R (obtained by integrating the forward conduction kernel) to the user, finally forming a comprehensive backward conduction kernel considering multi-source supply. A propagation and amplification model of demand uncertainty is established. Any prediction of future demand has uncertainty, and the longer the prediction time, the greater the uncertainty. The cumulative amplification effect of this uncertainty in the pipe network conduction process can be modeled. Specifically, the system will analyze the historical prediction error data to identify the functional relationship between the standard deviation σ of the prediction error and the prediction time domain τ. This relationship can be modeled as an uncertainty amplification function: σ(τ) = σ0•(1 + α×τβ ); where s(T) is the standard deviation of the error in predicting the demand at time T; s0is the standard deviation of the initial prediction error (i.e. the error at T = 0); a and b are amplification factors to be identified, which can be determined by maximum likelihood estimation on historical error data. Optionally, the model can also incorporate event factors, e.g. during special events such as holidays or extreme weather, an impact factor larger than 1 can be multiplied to reflect the sharp increase in demand uncertainty during these periods.

[0101] Generating scenario set for robust optimization based on uncertainty model. After obtaining the expected value of demand (by inverse conduction kernel calculation) and the uncertainty growing with time domain (by amplification function calculation), instead of simply using the expected value for optimization, the system generates a scenario set containing multiple possible demand scenarios for robust optimization. Specifically, the system can construct a probability scenario tree, dividing the next 24 hours into multiple stages, each stage generating several branches representing different possible demand realizations. At the same time, to ensure the resilience of the scheme to extreme events, the system also generates tail event scenarios with a certain small probability (e.g. 5%), such as a large user requiring emergency flushing due to an unexpected accident, resulting in a transient 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 collectively constitute the time-lag demand distribution matrix, which not only contains the expected value of demand, but also contains the possible distribution patterns around the expected value. Taking this distribution matrix as the input of the optimization model, the production scheme generated can meet various possible demand fluctuations with a high probability, thereby improving the reliability of the water supply system.

[0102] According to another aspect of the present application, the process of generating a scenario set for robust optimization can further include the following specific implementation details: After generating a large number (e.g. thousands) of possible future demand scenarios through the uncertainty amplification function, in order to reduce the complexity of subsequent optimization calculations, the system will use the K-means clustering algorithm for scenario reduction. Specifically, the K-means algorithm will divide (cluster) the high-dimensional scenario set into a pre-set number (e.g. 50) of clusters. The goal of the algorithm is to minimize the sum of the squared distances of each scenario to the centroid (weighted average scenario) of its cluster. After clustering is complete, the centroid of each cluster is selected as a representative scenario, and the probability of this representative scenario is set as the sum of the probabilities of all original scenarios in the cluster. The 50 carefully selected and weighted representative scenarios can highly approximate the probability distribution characteristics of the original thousands of scenarios, thereby improving the solving efficiency while ensuring robustness.

[0103] In one specific embodiment, assume in a city reclaimed water pipe network, there is one source node S (a certain reclaimed water plant) and one target node D (a certain high-tech industrial park user). Node S stably produces Q1 grade reclaimed water (COD concentration of 25 mg / L). Pipe network topology analysis identifies that there are 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 kilometers. The average flow velocity v1 = 1.2 meters / second is calculated according to the pipe diameter and real-time flow. Path P2: total length L2 = 8 kilometers. The average flow velocity v2 = 1.0 meters / second is calculated according to the pipe diameter and real-time flow. Water flow state: the water flow in both paths is calculated to be in turbulent flow state. Source water quality: Q1 grade. Calculate the physical transmission peak (Γ 11 ) parameters for path P1 (L1 = 5 km, v1 = 1.2 m / s, Q1 level): according to the formula under turbulent flow state: shape parameter k 11 = 2.2 + 0.12 * L1 / v1 = 2.2 + 0.12 *5000 / 1.2 ≈ 2.2 + 500 = 502.2 (here the formula application may need to adjust units or coefficients to conform to actual physical scales, this is an example calculation logic). To make the value reasonable, assume the formula is k1 = 2.2 + 0.12*L (L unit 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 of water flow through a 5-kilometer pipeline. Determine the decision response peak (Γ 21 ) parameters: since Q1 level high-quality water is transported, fast 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). Combine the conduction kernel K1(τ) of path P1: K1(τ) = w 11 •Γ(τ;2.8, 1.4) + w 21 •Γ(τ;4.0, 1.5). Where the weight w can be determined according to the path length, here for example, set w 11 =0.6, w21 = 0.4. The physical transport peak (Γ 12 ) parameters: 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. The decision response peak (Γ 22 ) parameters: same source water quality class, so same decision response pattern: shape parameter k 22 = 4.0; scale parameter θ 22 = 1.5; the conductance kernel K2(τ) for path P2 combined: K2(τ) = w 12 •Γ(τ; 3.16, 1.58) + w 22 •Γ(τ; 4.0, 1.5). The weights can be set to w 12 = 0.7, w 22 = 0.3.

[0104] By solving the hydraulic balance equation, the system calculates that 65% of the flow will choose path PI due to its shorter length and smaller hydraulic resistance, and 35% of the flow will choose path P2. Therefore, the flow allocation proportions are: al = 0.65, a2 = 0.35. The conductance kernels of the two paths are weighted and summed according to their flow allocation proportions to obtain the comprehensive conductance kernel K sd (τ) between the node pair (S, D): K sd (τ) = al•Kl(τ) + a2•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 comprehensive conductance kernel K sd(τ) is a complex, non-standard probability density function. If it is plotted as a graph, with the horizontal axis being the time delay τ (hours) and the vertical axis being the probability density, the curve will exhibit the following characteristics: Asymmetry: Due to the skewed nature of the gamma distribution itself, the entire curve will exhibit an obvious right-skewed (tailing) shape; Multi-peak / wide peak shape: Due to the superposition of four gamma distributions, the final curve will no longer be a simple single peak, but will likely appear as a wide peak with a clear shoulder, or even two peaks that are very close and partially overlapping. The main peak will be closer to the physical transmission delay of P1 (about 3.92 hours) because it occupies the dominant traffic, while the contribution of P2 and the existence of the decision delay peak will make the entire distribution last longer and have a more complex shape. This K sd The (τ) curve accurately describes the following physical and decision-making process: After water plant S performs a water quality adjustment, monitoring points in industrial park D begin to detect changes as early as approximately 3-4 hours later, with the most dramatic changes occurring around 4-5 hours. However, the complete arrival and stabilization of the entire water quality signal may take 8-10 hours. It can be seen that the adaptive multipath conduction kernel constructed in this embodiment can provide far richer and more accurate spatiotemporal dynamic information than a simple average delay time (e.g., (3.92 + 4.99) / 2 = 4.46 hours), which is the cornerstone for achieving high-precision supply-demand coordinated forecasting.

[0105] Further, to verify the superiority of the adaptive multi-path conduction considering the special characteristics of water network compared to the traditional single-path conduction, the following comparative experiments can be supplemented to illustrate: The experimental design can be based on a typical grid-shaped pipe network topology containing 20 nodes and 35 edges for 72-hour simulation prediction. Under the same supply-demand scenario, the multi-path model of the application (searching for the first 3 paths) and the single-path model (only considering the shortest path) as a control group are deployed. The experimental results show that the multi-path model of the application is superior to the single-path model in multiple key performance indicators: the average absolute percentage error (MAPE) of the multi-path model is 8.3%, while that of the single-path model is 15.7%, which improves the accuracy. Especially during the peak period of pipe network load, the single-path model has a high prediction deviation due to the neglect of water flow splitting, while the multi-path model has a lower deviation. In the simulation of the shutdown of the main pipeline, the single-path model incorrectly predicts that the water supply will be interrupted for 3 hours because the only path it relies on is invalid. However, the multi-path model of the application can accurately predict that the flow will be redistributed through the backup path, with a prediction error of less than 10 minutes, demonstrating its robust adaptive ability to changes in network topology. Through simulation tracing experiments, it is found that Q1 high-quality water does indeed prefer to choose the fast lane with faster central flow velocity, reaching the downstream earlier than Q3 water. The single-path model cannot capture the delay difference caused by the differential transmission of water quality (a physical characteristic unique to water networks and not present in packet transmission in communication networks), while the multi-path model of the application considers the specific propagation speed of water quality levels, with an average error of only 12 minutes in predicting the arrival time of water quality, which is much better than the 45 minutes of the single-path model. The above experimental data quantitatively proves that the multi-path conduction method proposed in the application is not a simple transfer of technology from other fields, but through the deep combination of the special physical laws of water power and water quality in reclaimed water pipe networks, it has produced unexpected technical effects.

[0106] In further embodiments, a virtual capacity pool breaking through the physical water plant boundary is constructed, generating a three-dimensional virtual capacity tensor. Research has found that the rated capacity of a water plant is not equivalent to being available at any time and any place. Influenced by the pipe network transmission delay, loss, and internal resource reuse opportunities, the effective supply capacity of a water plant to different users is dynamically changing. Therefore, it is necessary to construct a virtual capacity pool to describe all physical and virtual water sources in a unified space-time dimension.

[0107] As a preferred implementation, the construction process includes the following aspects: projection calculation of time-shifted capacity. The plant-side out-of-plant capacity is accurately converted into the available capacity on the user side. Specifically, the system will use the time-lag conduction kernel K ij (τ) from water plant i to user j to perform convolution operation on the real-time capacity curve P i (t) of the water plant to calculate the effective capacity contribution P ijeff (t): P ij eff (t) =∫P i (t-τ)×K ij (τ)×η ij (τ)dτ; where P i (t-τ) represents the production of water plant i at time t-τ in the past; K ij (τ) is the value of the conduction core at the time delay τ, representing the probability that the water leaving the factory 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 (for example, a loss of 2-5%) of water during long-distance transmission, as well as the attenuation of water quality indicators such as residual chlorine. Through this calculation, the production capacity pulse originally on the water plant side will be projected and diffused into a production capacity availability curve with a specific time distribution on the user side. This is the time-shifted capacity. Then, the floating capacity is optimized. For flexible production units in the pipeline network that can produce recycled water of various qualities (for example, Water Plant C can produce Q1 and Q2 grade water). This type of capacity is called floating capacity. The system needs to optimize its distribution ratio between different qualities. Specifically, the system will establish a switching cost model. For example, switching from producing Q1 to Q2 requires a transition time of 2 hours, during which the production capacity drops to 50% of the rated value, and incurs additional costs for chemicals and equipment wear. Based on this, the system constructs a two-layer optimization model. The upper layer decides the switching time, and the lower layer optimizes the output in 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) is the production during period t; the second term is a switching penalty term consisting of the switching penalty coefficient λ, which is 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 (for example, a maximum of two switchings per day, switching to high quality before the morning peak and switching to low quality after the evening peak), thereby accurately allocating floating capacity to different water quality levels.

[0108] Further, as an optional enhancement, the construction of the virtual capacity pool can also incorporate the water quality cascade network. The effluent of the network is also considered as a virtual water source, thus included in the capacity pool. Specifically, the system will identify the feasible cascade chains in the network. For example, industrial user A uses Q1 quality water, and its effluent quality drops to Q2, but still meets the demand of municipal user B. At this time, the system will consider the effluent of industrial user A (its water quantity and quality can be modeled by the quality conversion function) as a virtual capacity available to municipal user B. This virtual capacity also needs to be projected in time, i.e. calculate the transmission kernel from A to B, to determine when the effluent of A can reach B. In this way, the sources of the virtual capacity pool not only include physical water plants, but also hundreds of virtual water sources (i.e. the effluent outlet of the previous level user) in the network, enriching the dispatchable resources and providing a basis 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, and finally forms a unified three-dimensional virtual capacity tensor, which has a dimension of [water source ID x water quality level x time]. This tensor provides comprehensive, dynamic, and accurate supply-side input for subsequent global optimization.

[0110] In an exemplary embodiment, the switching cost model of the floating capacity can further be calculated as follows: the switching cost C switch is decomposed into three specific components: C switch = C loss + C chem + C wear . Wherein C loss is the capacity loss cost, which is calculated as C loss =0.5×P rated ×T trans ×price, P rated is the rated capacity of the water plant, T trans is the transition time required for water quality switching (e.g. 2 hours), price is the selling price of unit reclaimed water, and 0.5 represents that the capacity is reduced to 50% during the transition period. C chem is the water quality adjusting agent cost, which is the additional agent fee required for each switching, for example 500 yuan / time. C wear is the equipment wear cost, which reflects the additional wear and tear on valves, pump sets, etc. caused by frequent switching, for example 200 yuan / time. By accurately modeling this cost, the optimization algorithm can make switching decisions that are more in line with economic benefits.

[0111] In one embodiment of the present application, a more detailed and multi-level implementation method is provided for sparsification processing. By combining a variety of graph algorithms and quantitative indicators based on domain knowledge, the generated sparse connection matrix can retain the key network information to the greatest extent while improving computational efficiency. Specifically, it includes:

[0112] Perform micro-level sparsification to retain key transmission paths with high water flow influence. The goal is to identify those pipe segments that play a local core role in water flow transmission and distribution at the most basic edge level. Specifically, the system calculates a water flow influence index S for each edge (pipe segment) (i, j) in the network. ij Optionally, the indicator is calculated as follows: Define the instantaneous influence of the edge I ij (t) = Q ij (t)•N j •W j •exp(-λ•L ij ); where Q ij (t) is the real-time flow rate of the pipe section; N j is the total number of nodes that can be reached by its downstream node j (calculated by breadth-first search); W j is the weighted sum of the importance levels of all downstream users (e.g., industrial users have a weight of 3.0 and municipal users have a weight of 2.0); L ij is the length of the pipe section; λ is the distance attenuation coefficient. In order to consider the time lag effect, the time cumulative influence I ij cum =∫I ij (t-τ)•K ij (τ)dτ, that is, convolution of instantaneous influence with the conduction kernel of the edge. Taking into account the volatility of 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 ij The scores are sorted and the top 30% of edges are retained as key edges. For the unselected edges, their substitutability (i.e., the ratio of the length of the detour path to the direct path) is further calculated, prioritizing those edges that are easily replaced.

[0113] Perform meso-level sparsification to retain cross-water plant connections with strong synergy necessity. The purpose is to identify which water plants have close synergy relationships at the water plant cluster level and need to retain the connection paths between them for refined modeling. Specifically, the system calculates a synergy necessity index N for each pair of water plants (i, j) ijThe index integrates three dimensions: N ij = O ij •C ij •D ij ; wherein O ij is the service area overlap, measured by calculating the size of the range of users served by both water plants; C ij is the water quality complementarity, for example, a water plant mainly producing Q1 and a water plant mainly producing Q3 have high complementarity; for a water plant with switchable water quality, its complementarity will have an additional bonus; D ij is the distance factor, reflecting the distance of two water plants in the pipe network, the closer the distance, the greater the possibility of cooperation. After constructing the cooperation necessity matrix N between all pairs of water plants, the system can apply community detection algorithms such as spectral clustering to automatically divide all water plants into several cooperation groups. The sparsification strategy will be: keep all connection paths between water plants within the cooperation group, and for water plants belonging to different groups, only keep the strong cooperation connections with N ij values exceeding a high threshold.

[0114] Perform macro-level sparsification, extract and retain the backbone transmission network of the pipe network. The purpose is to identify the main roads that undertake long-distance and large-flow transportation tasks from the global level of the entire pipe network and ensure the integrity of these main roads. Specifically, the system can comprehensively use various graph theory algorithms to extract the backbone network: based on topological centrality: calculate the edge betweenness centrality of each edge, which measures the frequency of the shortest path between all node pairs passing through a certain edge, and the edge with high betweenness centrality is usually the bridge in the network and must be retained. Based on flow load: extract the average flow load of each edge from the water model, and identify a few core pipes carrying more than 70% of the network's flow as the flow skeleton. Based on network connectivity: under the premise of retaining all water plants and main user nodes, run the Kruskal algorithm to generate a minimum spanning tree (MST) that ensures basic connectivity of the network. The system will combine the edge sets identified by the above methods and apply, for example, the Louvain community detection algorithm for optimization to form a hierarchical sparse representation: that is, a few relatively dense internal connection regions and a few backbone channels connecting these regions.

[0115] In summary, through the sparsification strategies at the microscopic, mesoscopic, and macroscopic levels, the present embodiment can intelligently trim the original complex pipe network model into a sparse model that only retains key connections. This sparse model not only improves computational efficiency (for example, it can achieve more than 20 times computational speedup), but more importantly, since 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 the present application, a specific, engineering feasible solution strategy is provided for large scale multi-objective optimization problems to overcome the curse of dimensionality. Optionally, dimensionality reduction is performed on the high dimensional optimization decision space. Since it is computationally extremely expensive to solve a mixed integer programming problem with tens of thousands of dimensions directly, it is difficult to complete within the actual scheduling time window. By analyzing historical scheduling data, it is found that the optimal production scheme is not completely random, but contains specific, repeatable patterns and structures. Using a data driven approach, these internal structures are discovered, and the decision space is projected into a lower dimensional, more easily solvable subspace. As a preferred implementation, the dimensionality reduction can include a combination of the following techniques:

[0117] Linear dimensionality reduction is performed using Principal Component Analysis (PCA). Specifically, the system will first extract the globally optimal solutions for a number of scheduling periods (e.g., the past year) from the historical database, constructing a historical solution data matrix X, where each column represents a complete decision vector for a time instant. The covariance matrix of the data matrix is computed, and by solving its eigenvalue problem, a set of orthonormal bases that maximizes the explanation of the data variance, i.e., principal components, is found. Experimental data shows that usually the first k (e.g., k is about 200) principal components are enough to capture more than 95% of the variance in the original data. Thus, any high dimensional decision vector x can be projected into a k dimensional principal component space by a transformation matrix V, obtaining a low dimensional representation y = V T •x.

[0118] In the principal component space, further nonlinear compression is performed using sparse coding. To achieve more extreme compression, 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, a sparse representation y = Dα is sought, where α is a sparse coefficient vector with most elements being zero. The system can use algorithms such as K-SVD to learn the optimal dictionary D from a large amount of historical data y. When solving a new problem, the sparse coefficient α is quickly found using algorithms such as Orthogonal Matching Pursuit (OMP). By controlling the sparsity to be, for example, 10%, the original k dimensional vector can be compressed to only store k*10% non-zero coefficients and their positions, for example, only 20 non-zero coefficients are needed to reconstruct the solution.

[0119] Optionally, time series patterns can be mined as dedicated basis functions. The system can identify typical patterns that repeatedly occur in production plans, such as weekday peak patterns, weekend flat patterns, etc. Specifically, numerical time series can be discretized into symbol strings through techniques such as symbolic aggregate approximation (SAX), and frequent patterns can be efficiently found through data structures such as suffix trees. The identified, e.g., 12, basic patterns can be taken as a set of basis functions, and any complete production plan can be approximately represented as a linear combination of the basis functions.

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

[0121] Further, in order to improve the solving efficiency, a warm start optimization strategy can also be used. Specifically, the system maintains a library of optimization solutions for the last 7 days. When solving the problem for the current day, the system first finds the historical scenario most similar to the current predicted demand pattern from the solution library through time series similarity algorithms such as dynamic time warping (DTW) distance, and takes the corresponding historical optimal solution as the initial 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 more quickly, improving the convergence speed.

[0122] Through the above steps, this embodiment successfully transforms the theoretically difficult super-large scale optimization problem into an engineering feasible and quickly solvable computing task, providing performance guarantee for the practical application of the invention.

[0123] In another embodiment of the present application, a series of optional enhanced implementation methods are provided for the model construction and collaborative scheduling protocol involved in the present application, which can deeply couple the unique physical, chemical and operational characteristics of the reclaimed water pipe network. This can improve the fidelity of the model and the operability of the scheme, thereby providing strong support for the creativity of the overall scheme.

[0124] Optionally, a biofilm dynamic growth model is introduced when building the dynamic graph network model. Biofilm will grow on the inner wall of the pipe during long-term operation. To simulate this process more precisely, a differential equation model describing the evolution of biofilm thickness Δ over time can be introduced: dΔ / dt =μ max •(S / (K s +S))•Δ- b•Δ-τ w •Δ 2 / ρ b ; where Δ is the biofilm thickness; μ max is the maximum growth rate of microorganisms; S is the substrate (e.g. COD) concentration; K s is the half-saturation constant; b is the endogenous respiration decay rate within the biofilm; τ w is the shear stress of water flow on the pipe wall; ρ b is the biofilm density. Using real-time water quality, water temperature and hydraulic data of the pipe network, this equation is solved online by numerical methods such as the Runge-Kutta method to obtain the dynamic biofilm thickness Δ(t) of each pipe segment. Δ(t) is used to real-time correct the effective pipe diameter D eff (t) and the hydraulic roughness n(t) used in the model, so that the model can automatically adapt to the aging process of the pipe network and make more accurate long-term predictions.

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

[0126] Optionally, the physical and chemical constraints of water quality mixing are introduced when sparsifying or optimizing. To ensure water quality safety, the system can establish a water quality mixing compatibility matrix M compat (Q i , Q j ). Used to quantify two different quality water Q i and Q jIncompatibility risk 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 distribution or path planning, M compat > M threshold (Safety threshold) as a hard constraint. When the optimal route presents mixed risks, the system can automatically trigger an isolated delivery mode, such as staggering delivery times or activating backup pipelines in space, to prevent water quality incidents at the model level.

[0127] Optionally, when designing a coordinated scheduling protocol, water grid-specific coordination mechanisms can be introduced. Unlike power or communication networks, water, as a storable, quality-dependent, and delayed resource, can be designed for more sophisticated coordinated scheduling: A water quality buffer pool 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 scheduling logic combining FIFO and quality priority. The system dynamically determines water delivery batches based on the freshness of the water in the pool (determined by its residence time) and the urgency of downstream demand, striking a balance between meeting demand and ensuring water quality. A pre-dispatch credit mechanism with time lag compensation: To incentivize water plants to proactively and accurately pre-produce (rather than reactively respond), a pre-dispatch credit mechanism can be introduced. Water plants earn credit points for pre-production based on forecasted demand. Settlement is conducted after actual demand is confirmed, with credits awarded for accurate forecasts and deducted for deviations. Credit points can be used to secure priority or electricity price discounts in future scheduling. A dynamic alliance for water quality compatibility: The system can foster dynamic alliances between water plants based on real-time water quality monitoring results. For example, when Plant A's Q1 water quality exceeds expectations, reaching Q1+, the system can automatically match it with a specialized user requiring ultra-high-quality water, forming a value-added alliance. If Plant B's water quality inadvertently degrades, it can automatically form a complementary alliance with Plant C, which will provide high-quality water for blending and dilution. This kind of flexible and dynamic alliance is impossible with traditional static scheduling.

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

[0129] Based on the dynamic pipe network topology in the standardized spatiotemporal dataset, multiple hydraulic transmission paths between pairs of pipe network nodes are explored and identified.

[0130] In this embodiment, multiple hydraulic transmission paths include not only physical pipeline paths, but also corresponding decision-making control paths. Specifically, for the water supply process from water plant S to user D, there are three types of independent decision-making trigger mechanisms: the first type is automatic control decision-making based on the SCADA system, with a response time of approximately 5-15 minutes; the second type is manual intervention decision-making by the dispatcher, with a response time of approximately 30-60 minutes; the third type is cross-departmental coordination decision-making, which involves the linkage of multiple water plants and can take up to 2-4 hours. The response times of these decision paths vary significantly, forming the basis for the multimodal distribution of decision delays.

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

[0132] Specifically, the mathematical definition of the bimodal conduction kernel adopts the 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, whose 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 (transition flow), k 1p = 1.5 + 0.00025 × Re; when Re ≤ 2000 (laminar flow), k 1p = 1.2 + 0.0004 × Re. Scaling parameter θ 1p = L p / (v p ×3.6), where L p is the path length (km), v p is the average flow rate (m / s), and 3.6 is the unit conversion factor. The second peak corresponds to the decision response process, which is formed by 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 , which follows an exponential distribution. The convolution of these independent response times produces a gamma distribution with shape parameter k 2p =n, reflecting the complexity of the decision chain; the scale parameter θ2p =∑(T i ) / n, representing the average decision delay.

[0133] Further, the weight coefficients of the bimodal conduction kernel are determined. The weight coefficients w 1p and w 2p are determined based on the maximum likelihood estimation of historical operation data. The specific method is: collect the water quality transmission events on the path in the past 90 days, record the actual arrival time sequence {τ1, τ2,..., τ N} of each event. The likelihood function L(w 1p , w 2p ) =Π[w 1p × Γ(τ i ; k 1p , θ 1p )+ w 2p ×Γ(τ i ; k 2p , θ 2p )] is constructed. The optimal weight is obtained by iterative solution through the expectation maximization (EM) algorithm. In practical application, for the emergency water supply path (such as the hospital water supply line), the decision response is accelerated by the priority mechanism, resulting in a smaller w 2p (about 0.2-0.3); for the regular water supply path, the influence of physical transmission and decision response is comparable, w 1p and w 2p are close to 0.5; for the cross-regional coordinated water supply path, the decision coordination takes a long time, w 2p can reach 0.6-0.7. Alternatively, the weight coefficients can also be dynamically adjusted according to the time period. During working hours (8:00-18:00), the dispatch personnel are on duty, the decision response is faster, w 2p =0.3; during the night period (22:00-6:00), the automatic control is mainly relied on, the manual decision delay increases, w 2p =0.5; during holidays, the cross-department coordination is difficult, w 2p can be increased to 0.6.

[0134] Further, to verify the rationality of the bimodal distribution, the system will perform the Kolmogorov-Smirnov test to compare the fitting goodness of the actual transmission delay distribution and the bimodal model. When the test statistic D < 0.05 and the p value > 0.95, it is considered that the bimodal model is effective. If the test does not pass, the system will automatically switch to a three-peak or multi-peak mixed model to more accurately capture the multi-modal delay characteristics under complex decision scenarios.

[0135] In another exemplary embodiment, a coupling method of sparsification processing and water quality safety constraints is provided, specifically:

[0136] Based on the water quality data in the standardized spatiotemporal dataset, the water quality incompatibility risk when the water quality of two nodes connected in the dynamic graph network prediction model is mixed is quantitatively evaluated.

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

[0138] A coupling optimization model of sparsification decision and safety constraint is constructed.

[0139] In this embodiment, the sparsification is no longer simply pursuing the calculation efficiency, but is deeply coupled with the safety constraint. Specifically, a multi-objective optimization problem is constructed: min F =α×C comp +β×Rtotal + γ x L service ; where C comp is the computational complexity, equal to the square of the number of edges reserved; R total =∑(w ij x R ij ) is the total risk value, w ij is the flow weight of edge (i, j); L service is the service loss, measuring the impact of sparsification on water supply reliability; a, b, g are weight coefficients, adjusted according to application scenarios. The constraint conditions of this optimization problem include: connectivity constraint, ensuring that all user nodes are reachable; capacity constraint, keeping the total delivery capacity of the reserved paths no less than 120% of the demand; safety hard constraint, any edge with R ij > R threshold must be cut off, where R threshold is determined according to the type of downstream users: hospital users R threshold = 0.1; industrial users R threshold = 0.3; municipal green R threshold = 0.6.

[0140] The dynamic collaborative adjustment of sparsification and safety constraints is realized.

[0141] Optionally, the system uses a double-layer iterative algorithm to realize coupled optimization. The outer loop adjusts the sparsity parameter p, starting from the initial value p0= 0.5, with a step size of Ap= 0.05. The inner loop selects the edge set E sparse to be reserved under a fixed sparsity through a tabu search algorithm. In each iteration, the algorithm evaluates three types of indicators: the computation speedup ratio S = T full / T sparse ; the number of safety violations V = count(R ij > R threshold for (i, j) in E sparse ); and the service quality indicator Q =∑(demand met / demand total ). Where T full is the total computation time required for the system to run a prediction or optimization task once under the complete graph structure without sparsification; T sparse is the computation time required for the system to run the same task under the current sparsified graph structure; demand met is the total amount of user node water supply demand successfully met under the current sparse graph structure; demand total=The theoretical total water supply demand for all user nodes. When a safety violation is detected (V > 0), the algorithm immediately triggers safety-first mode: the offending edge is added to a taboo list, preventing it from being selected in subsequent iterations; an alternate path search is activated to find alternative routes around high-risk areas; and, if necessary, the sparsity requirement is reduced to increase redundant connections to ensure safety. This dynamic adjustment mechanism ensures that sparsity is never achieved at the expense of safety.

[0142] In actual testing of this embodiment, compared to traditional independent sparsification methods, the coupled optimization scheme reduced the incidence of water quality safety incidents while maintaining the same computational efficiency. This is due to the sparsification process's proactive identification and avoidance of risk factors, rather than post-event remediation. Optionally, a reinforcement learning mechanism can be introduced to train the sparsification strategy network using historical operating data. This network takes the network status (topology, flow rate, water quality) as input and outputs the probability of retaining 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 the safety indicator function, which takes the value 1 when there is no safety violation and 0 otherwise; C comp_normalized is the normalized value of computational complexity, R total_normalized is the normalized result of the total risk value, L service_normalized is the normalized value of the service loss. Through Deep Q Network (DQN) training, the system can learn the optimal sparsification strategy that balances efficiency and security.

[0143] In another exemplary embodiment, a detailed process for deeply integrating the multipath conduction model with the water quality cascade utilization network is provided. Specifically, a cascade utilization delay matrix is ​​constructed based on the multipath conduction core. In the cascade utilization network, the water output of the upper user becomes the water source for the lower user. The accuracy of the 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 ladder delay tensor T[i, j, τ] is constructed, where the element T ijτ represents the probability density of water from node i arriving at node j after a time delay of τ. The construction of this tensor takes into account the coupling effect unique to the cascade. When user i supplies water to multiple downstream users j and k at the same time, there is competition for diversion, and the actual flow of each path will affect each other. By solving the coupled hydraulic balance equations: Σ(Q ij ) =Q i_out (conservation of flow); h i -h j = f ij (Q ij ) (energy conservation); where Qi_out is the total water outflow for user i; h is the node head; f ij is the pipe resistance function. The dynamic flow allocation Q ij (t) is solved, and the conduction kernel parameters of each path are updated. The multi-path information is used to optimize the cascade link selection. Traditional cascade matching only considers water quality and water quantity matching degrees, while ignoring the importance of delay matching. In this embodiment, the cascade feasibility score function is extended as: S cascade (i, j) = M quality ×M quantity ×M timing ×(1 - R risk ); wherein M quality = exp(-|Q i_out - Q j_req | / Q ref ) is the quality matching degree; M quantity =min(V i_out / V j_req , 1) is the water quantity matching degree; M timing is the delay matching degree; R risk is the path risk, Q j_req is the target water supply demand flow of node j (lower-level user), Q ref is the reference flow value, V i_out is the available water volume of node i (upper-level user) in the current period, and V j_req is the target water volume of node j in the current period. The calculation of the delay matching degree M timing fully utilizes the multi-path conduction information: M timing =∫K ij (τ)×W j (τ) dτ; wherein W j (τ) is the demand time window function of user j. When the main peak of K ij (τ) overlaps with the high demand period of W j (τ), M timing is close to 1, indicating good delay matching. The complete distribution K ij (τ) provided by the multi-path model can more accurately evaluate the time window matching degree compared to a single delay value. Experimental data show that after introducing the delay matching degree, the average utilization rate of the cascade chain is improved, and the flow interruption event is reduced. This is because the multi-path model identifies the long-tail characteristics of the delay distribution, avoiding supply interruption caused by low-probability but long-delay water flow. The multi-path redundancy enhances the resilience of the cascade network. The path diversity revealed by the multi-path conduction model provides a natural redundancy backup for the cascade network. Specifically, a path redundancy index R redundancy (i, j) = 1 - Π(1 - p k ) is calculated for each cascade connection; wherein pk is the available probability of the kth path. When R redundancy > 0.8, even if the main path fails, there are still enough backup paths to maintain the cascade operation. The system builds the resilience enhancement strategy of the cascade network accordingly. In the network design stage, node pairs with high path redundancy are preferentially selected to build the cascade chain. In the operation stage, the state of each path is monitored in real time, and when the performance of a certain path is detected to decline, the flow allocation is adjusted in advance to smoothly transition to the backup path. This active resilience management based on multi-path information shortens the average recovery time of the cascade network.

[0144] Optionally, the system can also implement multi-path cooperative scheduling. For critical cascade connections, multiple paths are simultaneously enabled for transmission, and flow smoothing is achieved through inter-path phase modulation. The specific method is as follows: the transmission phase of each path is set as φ k = 2π×k / N, where N is the total number of paths. The instantaneous flow of each path is Q k (t) = Q avg + A×sin(ωt + φ k ), where Q avg is the average flow, A is the modulation amplitude, and ω is the modulation frequency. The superposition of multi-path flow produces a stable total flow, reducing the flow fluctuation in the cascade chain and improving the water supply stability of downstream users. Through the above mechanism, the combination of multi-path conduction and cascade utilization produces a synergistic enhancement effect: the fine time delay information provided by the multi-path model optimizes the cascade matching and improves the network efficiency; the path redundancy enhances the resilience and reliability of the cascade network; the multi-path cooperative scheduling smooths the cascade flow and improves the water supply quality. The overall system performance is improved more than the simple superposition of each technology applied alone, embodying the true technology fusion.

[0145] The preferred embodiments of the present application are described in detail above, but the present application is not limited to the specific details of the above-described embodiments. Within the technical concept of the present application, various equivalent transformations of the technical solutions of the present application can be made, and these equivalent transformations all belong to the protection scope of the present application.

Claims

1. A method for collaborative forecasting of urban recycled water supply and demand, characterized in that: include: Process multi-source heterogeneous data streams from the recycled water network to generate standardized spatiotemporal datasets; Based on the standardized spatiotemporal dataset, a dynamic graph network prediction model is constructed to characterize the supply and demand relationship between nodes in the pipe network and the spatiotemporal transmission characteristics of water quality. Utilize dynamic graph network prediction models and real-time demand signals to generate quality-based demand forecast results and cross-factory collaborative production plans.

2. The method according to claim 1, characterized in that Build a dynamic graph network prediction model, including: Based on the dynamic pipe network topology in the standardized spatiotemporal dataset, the predetermined hydraulic transmission paths between the pipe network node pairs are explored and identified; A path conduction kernel is constructed for each hydraulic transmission path. The path conduction kernel presents a bimodal distribution to respectively characterize the physical transmission delay and the decision response delay.

3. The method according to claim 2, characterized in that Construct a path conduction core for each hydraulic transport path, including: Extracting the water quality grades of recycled water associated with the pathways from the standardized spatiotemporal dataset; Determine the parameters of the peak in the bimodal distribution that characterizes the decision response delay based on the water quality level.

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

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

6. The method according to claim 5, characterized in that A one-dimensional convection-reaction-diffusion equation is established, including parameterization based on a standardized spatiotemporal dataset, wherein the parameterization comprises: Determine the reaction degradation coefficient in the one-dimensional convection-reaction-diffusion equation as a function of pollutant concentration and water quality grade; A dynamic growth model of biofilm on pipe walls was established to dynamically modify the effective pipe diameter and hydraulic roughness used in the one-dimensional convection-reaction-diffusion equation. Differentiated propagation velocities are set for the one-dimensional convection-reaction-diffusion equation according to the water quality level, where the propagation velocity of high-quality water in the central area of ​​the pipe is higher than that of low-quality water near the wall.

7. The method according to claim 1, characterized in that Also includes: Perform sparse processing on the dynamic graph network prediction model to generate a sparse connection matrix; Among them, the generation of quality demand forecast results and cross-factory collaborative production plans is carried out using a sparse connection matrix.

8. The method according to claim 7, characterized in that Generate a sparse connectivity matrix, including: Based on the topological change rate, demand fluctuation rate 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, ρ(t)=ρ base The functional relationship of ×exp(-β×C(t)) determines the connection retention ratio ρ(t) and generates a sparse connection matrix, where ρ base is the basic retention ratio, and β is the sensitivity parameter.

9. The method according to claim 7, characterized in that Generate a sparse connectivity matrix, you can also: Based on water quality data from standardized spatiotemporal datasets, the water quality incompatibility risk when the water quality connecting the two end nodes in the dynamic graph network prediction model is mixed is quantitatively evaluated; The network connections whose water quality incompatibility risk exceeds the preset mixing risk threshold are cut off or weakened to generate a sparse connection matrix.

10. The method according to any one of claims 1 to 9, characterized in that Generating cross-factory collaborative production plans is achieved by solving a large-scale multi-objective optimization problem. The solution process includes: Decompose 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; Through a distributed coordination protocol, the solutions to low-dimensional sub-problems are coordinated to obtain the global optimal solution to large-scale multi-objective optimization problems; Among them, large-scale multi-objective optimization problems are constructed based on dynamic graph network prediction models and real-time demand signals.

Images