Water pump pressure controller dynamic pressure prediction method based on digital twinning

Abstract

The application relates to the technical field of computer technology, digital twinning and big data processing technology, in particular to a dynamic pressure prediction method for a water pump pressure controller based on digital twinning, a digital twinning mapping model of a water pump pipe network system is constructed to extract a spatial topological structure, and time sequence sensing data and control instruction sequences of each node are collected; the data are input into a graph neural network to extract spatial features of pressure nonlinearity hysteresis, and a space-time feature matrix is output after being decoupled through a time sequence attention mechanism; the matrix is input into a friction physics prediction model, a benchmark network and a residual compensation network are superposed to output a preliminary pressure prediction sequence; a dynamic time warping algorithm is used to compare the measured sequence and the preliminary prediction sequence, a hysteresis phase difference is extracted to be converted into a time compensation factor, a time axis is translated and calibrated, and a final pressure prediction sequence is output. The method provides an implementation path for representing dynamic conduction relationships among pipe network nodes and calibrating time offset of a prediction sequence.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of computer technology, digital twins, and big data processing technology, specifically to a dynamic pressure prediction method for a water pump pressure controller based on digital twins. Background Technology

[0002] With the rapid development of computer technology, digital twins, and big data processing technologies, intelligent control of the operational status of complex physical systems through computer-based digital modeling, time-series data analysis, and state prediction has become the industry mainstream. In practical applications, the pipeline data stream processed by the pump pressure control system has complex computational characteristics of multiple variables, large time delays, and strong coupling. When multiple water terminals are frequently opened and closed or the system operates for long periods, the time-series data of the sensing terminals such as pressure at each node of the pipeline network are often severely disturbed by multi-dimensional physical boundaries and dynamic noise.

[0003] However, existing data processing methods and computational models have significant limitations in dealing with such scenarios. Furthermore, water supply network pressure prediction is simultaneously affected by node topology, pump and valve control events, water usage fluctuations, and pipeline hydraulic parameters. Predictions based solely on historical pressure at a single node or fixed topological relationships may fail to characterize the dynamic transmission relationships between different nodes and the time response differences after control commands are applied. When transient changes occur in operating conditions, the predicted sequence may also deviate from the measured sequence in time. Therefore, when faced with the nonlinear hysteresis response of pressure at various nodes in the network, existing algorithms cannot accurately predict the pressure of pump pressure controllers.

[0004] To address the aforementioned technical problems, this invention proposes a dynamic pressure prediction method for a water pump pressure controller based on digital twins. Summary of the Invention

[0005] The purpose of this invention is to provide a dynamic pressure prediction method for a water pump pressure controller based on digital twins. By deeply embedding the physical laws of fluid mechanics into the training and prediction process of the neural network, it aims to provide a high-fidelity and high-reliability pressure prediction basis for the safe closed-loop control of the water pump.

[0006] To achieve the above objectives, the present invention provides the following technical solution: A dynamic pressure prediction method for water pump pressure controllers based on digital twins includes: A digital twin mapping model of the water pump pipeline network system is constructed, and the spatial topology is extracted based on the digital twin mapping model; the time-series sensing data and measured pressure data of each node are collected periodically to obtain the control command sequence of the water pump pressure controller; The time-series sensing data and the control command sequence are input into a graph neural network, and the spatial features of nonlinear pressure hysteresis at nodes are extracted by combining the spatial topology. The spatial features are decoupled using a time-series attention mechanism to output a spatiotemporal feature matrix. The spatiotemporal feature matrix is ​​input into a tribophysical prediction model, and the global base pressure sequence and the local hysteresis fluctuation sequence are output by the baseline network and the residual compensation network in the tribophysical prediction model, respectively. The global base pressure sequence and the local hysteresis fluctuation sequence are superimposed to obtain a preliminary pressure prediction sequence. The measured pressure data is compared with the preliminary pressure prediction sequence using a dynamic time warping algorithm. The hysteresis phase difference is extracted and converted into a time compensation factor. The time axis of the preliminary pressure prediction sequence is shifted and calibrated using the time compensation factor, and the pressure prediction sequence is output.

[0007] Preferably, the process of extracting the spatial topology includes: Based on the digital twin mapping model, the physical topology of the pipeline network nodes is analyzed. Combining the water flow velocity, friction coefficient and pressure difference of each pipeline, the asymmetric water flow transmission delay between adjacent nodes along the flow direction and the counterflow direction is dynamically calculated. The asymmetric water flow transmission delay is input into the exponential decay kernel function to calculate the dynamic spatial decay weight between each node. Based on the dynamic spatial decay weights and the physical topology relationship of the nodes, a directed weighted adjacency matrix of the water pump network system is constructed and updated as the spatial topology.

[0008] Preferably, the nodes include pump station outlet and pipeline end node, and the time-series sensing data includes pipeline pressure, water flow rate, valve opening collected by each node according to the sampling period, as well as pressure change rate and water resistance dynamic gradient calculated based on the pipeline pressure and water flow rate. The control command sequence includes the frequency modulation pulse duty cycle change rate issued by the water pump pressure controller, the water pump start / stop status, and the command issuance timestamp. Before inputting the time-aware data and the control command sequence into the graph neural network, the time-aware data of each node is subjected to time-series interpolation alignment and sliding window normalization based on the command issuance timestamp.

[0009] Preferably, the process of extracting spatial features includes: Using the edge-conditional graph convolutional layer based on the directed weighted adjacency matrix, the control instruction sequence is embedded as a conditional feature into the time-series-aware data of each node; Along the water flow and wave propagation directions defined by the directed weighted adjacency matrix, multi-order neighborhood message passing and asymmetric spatial aggregation are performed on the temporal sensing data of each node that has completed the conditional feature embedding to obtain the spatial features.

[0010] Preferably, the process of decoupling the spatial features from the nodal pressure nonlinear hysteresis features includes: The spatial feature vector is linearly transformed to obtain the query vector, and the control command sequence of each historical time step is linearly transformed to obtain the key vector and the value vector respectively; A time-delay bias mask is constructed based on the asymmetric water flow transmission delay. Multi-head attention association weights between the query vector and the key vector are calculated using these maskes. The value vectors are weighted and fused according to these multi-head attention association weights. A time-delay bias mask is constructed based on the cumulative transmission delay, and this mask is added to the attention scores of the query vector and the key vector, so that the weighting process of each node for historical control commands is constrained by the corresponding transmission timescale. The weighted control command representation is then fused with the node spatial features to generate the spatiotemporal feature matrix.

[0011] Preferably, the process of outputting the global baseline pressure sequence and the local hysteresis fluctuation sequence includes: The spatiotemporal feature matrix is ​​subjected to feature projection and linear mapping in the node space dimension. The steady-state head distribution is extracted through the multilayer perceptron of the baseline network, and the global basic pressure sequence is output. The global basic pressure sequence is used as a conditional constraint vector, and combined with the spatiotemporal feature matrix, it is input at time steps as the gated recurrent unit network of the residual compensation network. Using the reset and update gates of the gated recurrent unit network, the hidden state of the gated recurrent unit network is updated in a time-series manner under the modulation of the condition constraint vector. Based on the physical information in the frictional physics prediction model, consistency penalty constraints are executed to fit the residual fluctuations caused by dynamic hysteresis and output the local hysteresis fluctuation sequence.

[0012] Preferably, the process of converting to the time compensation factor includes: Calculate the numerical difference and the first-order pressure change rate difference between the measured pressure data and the preliminary pressure prediction sequence within a preset sliding time window, and construct a gradient-weighted distance matrix; Under the preset path slope constraint, the cumulative distance is recursively calculated based on the gradient weighted distance matrix to determine the alignment path corresponding to the minimum cumulative distance; The average offset steps of the alignment path in the time axis direction are extracted as the hysteresis phase difference; The time compensation factor is generated by multiplying the hysteresis phase difference by the period.

[0013] Preferably, the time axis of the preliminary pressure prediction sequence is shifted and calibrated, including: Based on the time compensation factor, the time axis of the preliminary pressure prediction sequence is shifted along the time feedforward direction; for the end prediction gaps generated after the time axis phase shift, the boundary trend extrapolation is performed using the global basic pressure sequence to fill the gaps. On the shifted time axis after filling, the preliminary pressure prediction sequence after extrapolation and filling is resampled using the conformal cubic Hermitian interpolation algorithm. This process ensures that the interpolation curve passes through the original node values ​​and maintains monotonicity within the monotonic interval of the original data, reducing additional overshoot, and outputs the pressure prediction sequence.

[0014] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. This invention combines a digital twin model with a graph neural network. It constructs a directed weighted adjacency matrix by calculating the asymmetric water flow transmission delay and decouples the node pressure hysteresis features using a multi-head attention mechanism. This method quantifies the spatial hysteresis effect of water flow transmission and integrates it into the graph topology, generating candidate response representations related to the history of control commands. These representations are then fused with spatial features, achieving the extraction of a dynamic pressure spatiotemporal evolution feature matrix without hysteresis interference.

[0015] 2. This invention extracts the steady-state head distribution using a multilayer perceptron in a baseline network and recursively updates the hidden state of a conditionally constrained vector modulation gated recurrent unit network. Simultaneously, the predicted comprehensive pressure is substituted into the one-dimensional unsteady flow hydrodynamic equations to calculate the unbalanced residual value, and a consistency penalty loss function is constructed to optimize the network parameters through backpropagation, ensuring that the local hysteresis fluctuation sequence output by the network satisfies the laws of fluid physics conservation.

[0016] 3. This invention utilizes a dynamic time warping algorithm to extract the hysteresis phase difference between the measured and predicted sequences, generating a time compensation factor to perform feedforward phase shift calibration on the time axis of the predicted sequence. It combines the global baseline pressure sequence to extrapolate and fill the boundary trends of the final prediction gaps, and employs a conformal cubic Hermitian interpolation algorithm on the shifted time axis for resampling, preserving the local extrema and first-order pressure change rate of the data, and outputting a predicted sequence aligned with a standard clock. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating a dynamic pressure prediction method for a water pump pressure controller based on digital twins. Figure 2 A data flow diagram of a dynamic pressure prediction method for a water pump pressure controller based on digital twins; Figure 3 This is a data flow diagram illustrating the decoupling of the spatial features from the nonlinear hysteresis features of the node pressure in this invention. Detailed Implementation

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

[0019] Example 1: Please see Figures 1 to 2 This invention provides a dynamic pressure prediction method for a water pump pressure controller based on digital twins, the technical solution of which is as follows: A dynamic pressure prediction method for water pump pressure controllers based on digital twins includes: A digital twin mapping model of the water pump pipeline network system is constructed, and the spatial topology is extracted based on the digital twin mapping model; the time-series sensing data and measured pressure data of each node are collected periodically to obtain the control command sequence of the water pump pressure controller; The time-series sensing data and the control command sequence are input into a graph neural network, and the spatial features of nonlinear pressure hysteresis at nodes are extracted by combining the spatial topology. The spatial features are decoupled using a time-series attention mechanism to output a spatiotemporal feature matrix. The spatiotemporal feature matrix is ​​input into a tribophysical prediction model, and the global base pressure sequence and the local hysteresis fluctuation sequence are output by the baseline network and the residual compensation network in the tribophysical prediction model, respectively. The global base pressure sequence and the local hysteresis fluctuation sequence are superimposed to obtain a preliminary pressure prediction sequence. The measured pressure data is compared with the preliminary pressure prediction sequence using a dynamic time warping algorithm. The hysteresis phase difference is extracted and converted into a time compensation factor. The time axis of the preliminary pressure prediction sequence is shifted and calibrated using the time compensation factor, and the pressure prediction sequence is output.

[0020] The process of extracting the spatial topology includes: Based on the digital twin mapping model, the physical topology of the pipeline network nodes is analyzed. Combining the water flow velocity, friction coefficient and pressure difference of each pipeline, the asymmetric water flow transmission delay between adjacent nodes along the flow direction and the counterflow direction is dynamically calculated. The asymmetric water flow transmission delay is input into the exponential decay kernel function to calculate the dynamic spatial decay weight between each node. Based on the dynamic spatial decay weights and the physical topology relationship of the nodes, a directed weighted adjacency matrix of the water pump network system is constructed and updated as the spatial topology.

[0021] By retrieving the two-dimensional design drawings of the water supply pipelines of the water pump network system and the pipeline asset ledger data, the physical length of each pipe section is extracted. Pipe inner diameter The system includes the material of the pipes, and the spatial three-dimensional coordinates of each pump interface, valve location, and pipe intersection. The pump interfaces, valve locations, and pipe intersections corresponding to these spatial three-dimensional coordinates are mapped as abstract nodes in a data structure. The physical pipes connecting adjacent nodes are mapped as edges in the data structure. A mapping table containing node indices and edge relationships is established. This graph model is used to generate a digital twin mapping model of the pump network system. Based on the digital twin mapping model, the connection status between nodes in the mapping table is read. If there is a physical pipe connection between two nodes, the topological connectivity variable between them is marked as 1; otherwise, it is marked as 0. The physical topological relationships of the nodes in the entire network are then obtained.

[0022] Pressure transmitters and electromagnetic flowmeters deployed at various nodes of the pipeline network are used to synchronously collect pipeline pressure and water flow rate values ​​at adjacent nodes. The absolute difference between the pipeline pressure values ​​at adjacent nodes is calculated to obtain the node pressure difference. The basic wave velocity of sound wave propagation in the water medium under the current water temperature environment is obtained by consulting the underwater acoustics physics handbook. For the water flow direction, the basic wave velocity is added to the water flow velocity to obtain the downstream absolute wave velocity; for the upstream direction, the basic wave velocity is subtracted from the water flow velocity to obtain the upstream absolute wave velocity. Based on the calibrated pressure wave propagation speed (i.e., the downstream and upstream absolute wave velocities) and the physical length of the pipe section, the direction-dependent propagation delay is obtained. Since the downstream and upstream absolute wave velocities are different, the downstream propagation delay and the upstream propagation delay are obtained with unequal values, which together constitute the asymmetric water flow propagation delay.

[0023] The asymmetric water flow transmission delay calculated between all adjacent nodes in the pipe network within the current acquisition period (i.e., the sampling period). Calculate the arithmetic mean to obtain the global average delay scale. Using the natural constant Using the base as the time delay for the asymmetric water flow propagation between any adjacent nodes. Divide by the global average delay scale The negative value of the quotient is used as the exponent to construct the exponential decay kernel function. This function is then used to calculate the asymmetric water flow propagation delay between adjacent nodes. By inputting the exponential decay kernel function one by one and performing natural exponential decay calculation, the dynamic spatial decay weight between each adjacent node with a value range between 0 and 1 is calculated. The calculation formula is as follows: Initialize a square matrix with the same number of rows and columns as the total number of network nodes. Traverse each element of this matrix according to the physical topology of the nodes. If there is no physical pipeline connection between the two nodes represented by the row and column indices of a matrix element, assign the value zero to that element. If a physical pipeline connection exists, fill the matrix element with the dynamic spatial attenuation weight calculated for the direction from the row index node to the column index node, based on whether the direction is in the direction of water flow or the counter-flow direction. After completing the traversal and assignment of all matrix elements, a directed weighted adjacency matrix of the pump network system with asymmetric numerical characteristics is constructed. In each subsequent sampling period, when new measured data is acquired, repeat the above calculation steps to replace the attenuation weight values ​​in this matrix to update it, and output this directed weighted adjacency matrix directly as the spatial topology.

[0024] By combining the pressure wave propagation velocity and the underwater acoustic wave velocity, and introducing a friction damping penalty term, independent quantification of the asymmetric propagation delay in the downstream and upstream directions is achieved. An exponentially decaying kernel function is used to transform this delay into dynamic weights, and a directed weighted adjacency matrix is ​​constructed. This approach enables the spatial topology to reflect the directional differences and physical hindrance effects of pressure wave propagation, providing underlying data support consistent with fluid dynamics mechanisms for graph neural networks to extract spatial features.

[0025] The nodes include pump station outlets and pipeline end nodes. The time-series sensing data includes pipeline pressure, water flow rate, valve opening collected by each node according to the sampling period, as well as pressure change rate and water resistance dynamic gradient calculated based on the pipeline pressure and water flow rate. The control command sequence includes the frequency modulation pulse duty cycle change rate issued by the water pump pressure controller, the water pump start / stop status, and the command issuance timestamp. Before inputting the time-aware data and the control command sequence into the graph neural network, the time-aware data of each node is subjected to time-series interpolation alignment and sliding window normalization based on the command issuance timestamp.

[0026] In the physical space of the water pump network system, the main outlet pipe interface of the pump station is clearly defined and selected as the pump station outlet node, and the terminal interfaces of each branch at the very end of the entire water supply network are selected as the network end nodes. Pressure transmitters, electromagnetic flow meters, and digital valve opening sensors are installed at the pump station outlet node and each network end node, respectively. According to the sampling period, the pressure transmitter collects the current pipeline pressure value at each node, the electromagnetic flow meter collects the current water flow rate value at each node, and the digital valve opening sensor collects the current valve opening value at each node, thus obtaining the pipeline pressure, water flow rate, and valve opening value collected at each node according to the sampling period.

[0027] For the pipe pressure and water flow data acquired in each sampling period, derivative physical indices are calculated. The pipe pressure value acquired in the current sampling period is subtracted from the pipe pressure value acquired in the previous adjacent sampling period, and then divided by the time interval of the sampling periods to calculate the pressure change rate at each node. Historical flow operation records of the pump network system are retrieved, and the minimum absolute value of the non-zero water flow rate is extracted from the records. This minimum value is then multiplied by ten to the power of negative cube to calculate the non-zero minimum constant. The pipeline pressure values ​​collected during the current sampling period. Divide by the water flow rate value collected within the same sampling period. The square of the nonzero minimum constant obtained by calculation The sum of these values ​​is used to calculate the dynamic gradient of water resistance at each node. The calculation formula is as follows: The parameters are: pipeline pressure, water flow rate, valve opening, pressure change rate, and dynamic gradient of water resistance. The combined data constitutes the complete time-series sensing data of the nodes.

[0028] The control command sequence issued by the controller is obtained through the internal command communication interface of the water pump pressure controller. For each control command issued, three core parameters are extracted: the frequency modulation pulse duty cycle change rate obtained by dividing the amplitude of the frequency modulation pulse signal's duty cycle change in adjacent adjustment cycles by the time interval; the Boolean value of the water pump motor's operating status, recording the motor as running (on) and shutting it down (off) to form the water pump's on / off state; and the absolute system clock time at the moment the command is issued is recorded as the command issuance timestamp.

[0029] Because the acquisition frequency of IoT multi-source sensors is asynchronous with the timing of controller command issuance, alignment and normalization are performed before inputting the time-series sensing data and control command sequence into the graph neural network. The obtained command issuance timestamp serves as the time alignment reference axis. For the pipe pressure, water flow rate, valve opening, pressure change rate, and water resistance dynamic gradient data collected by each node at sampling times adjacent to the command issuance timestamp, linear interpolation functions are used to perform interpolation calculations on the time axis, uniformly aligning the data from heterogeneous acquisition frequencies to the discrete time grid points corresponding to the command issuance timestamp. Continuous physical quantities (such as pipe pressure and water flow rate) are aligned on the time axis using linear interpolation functions; while discrete control quantities (such as Boolean values ​​for pump on / off states) are processed using forward hold or event coding to avoid interpolation creating physically non-existent intermediate states. Based on the aligned time series, a sliding window with a fixed historical time step is defined; the maximum values ​​of various sensing data and control command sequences within the sliding window are statistically analyzed. and minimum value Using the range normalization mathematical formula, all the perceived data of each node within the window are compared with the original values ​​of the control command sequence. The normalized value is calculated by linearly scaling to a closed interval between 0 and 1. The calculation formula is as follows: The sliding window normalization process is completed. By setting a non-zero denominator protection term, the risk of numerical anomalies in ratio calculations near zero flow is reduced. Interpolation of continuous time-series sensing data is performed based on the control command timestamp, and the input features within the sliding window are normalized to provide a time-consistent input tensor for the graph neural network.

[0030] A non-zero minimum constant is introduced to calculate the dynamic gradient of water resistance, eliminating overflow during zero-flow calculations and enhancing the algorithm's robustness. Linear interpolation is performed based on the command timestamp to eliminate asynchronous sampling time differences between multi-source sensors and control commands, achieving heterogeneous data alignment. Combined with sliding window normalization, the dimensional differences of multi-dimensional physical quantities are eliminated, providing stable and identically distributed standard tensor inputs for subsequent graph neural networks, accelerating model gradient convergence.

[0031] The process of extracting spatial features includes: Using the edge-conditional graph convolutional layer based on the directed weighted adjacency matrix, the control instruction sequence is embedded as a conditional feature into the time-series-aware data of each node; Along the water flow and wave propagation directions defined by the directed weighted adjacency matrix, multi-order neighborhood message passing and asymmetric spatial aggregation are performed on the temporal sensing data of each node that has completed the conditional feature embedding to obtain the spatial features.

[0032] The temporal sensing data of each node, after temporal interpolation alignment and sliding window normalization, is used as the node feature input vector of the graph convolutional network. The dynamic spatial decay weight values ​​of each directed edge in the directed weighted adjacency matrix are extracted. The duty cycle change rate of the frequency modulation pulse issued by the water pump pressure controller and the water pump start-up and shutdown status are input into the multilayer perceptron network for nonlinear mapping to calculate and generate a high-dimensional conditional feature vector. This conditional feature vector is broadcast along the graph structure and embedded into the temporal sensing data of each node, completing the conditional feature embedding at the node level. The multilayer perceptron network synchronously outputs the feature modulation weight matrix associated with the pipeline edge. The feature modulation weight matrix and the scalar dynamic spatial decay weight of the corresponding directed edge are multiplied and scaled by scalar and matrix multiplication to construct the dynamic convolution kernel weights of the edge-conditional graph convolutional layer containing control command constraints.

[0033] Using the water flow direction and wave propagation direction defined by the directed weighted adjacency matrix as the topological path, a multi-level neighborhood propagation layer is set for the graph convolutional network. In the first-level neighborhood message passing stage of the graph convolution, for any target node, the features of all directly connected downstream and upstream neighboring nodes are extracted using the dynamic convolution kernel weights of the edge-conditional graph convolutional layer. After completing the first-level neighborhood information aggregation, the neighborhood range is expanded to the second and higher-level neighborhoods in the same way, performing recursive spatial message passing across multiple neighborhoods to diffuse the instantaneous state of the local pipeline to the global space of the pipeline network.

[0034] In the multi-order neighborhood message passing process, different spatial aggregation weight matrices are used to perform asymmetric spatial weighted summation operations on the features of neighboring nodes in the downstream and upstream directions. Active aggregation weights reflecting the forward propulsion and pressure transmission of water flow are assigned to the features of nodes converging in the downstream direction, while damping aggregation weights reflecting the reverse reflection of water pressure fluctuations and water hammer hysteresis are assigned to the features of nodes converging in the upstream direction. The feature vectors of neighboring nodes after asymmetric spatial aggregation are concatenated with the feature vectors of the target node itself and nonlinearly activated to calculate and generate the spatial features, i.e., the spatial feature vector, which includes the head loss along the pipeline network and the dynamic pressure transmission hysteresis.

[0035] Control commands are transformed into conditional features embedded in nodes, and dynamic convolutional kernels are constructed using spatial decay weights. In multi-order neighborhood message passing, downstream pressure propulsion and upstream reflection damping are distinguished, and asymmetric spatial weighted aggregation is performed. This mechanism directly couples control intervention to graph topology computation, deeply characterizing the dynamic transmission hysteresis evolution caused by head loss and water hammer effect along the flow path, and enhancing the physical mechanism representation capability of spatial feature tensors.

[0036] Please see Figure 3The process of decoupling the spatial features from the nodal pressure nonlinear hysteresis features includes: The spatial feature vector is linearly transformed to obtain the query vector, and the control command sequence of each historical time step is linearly transformed to obtain the key vector and the value vector respectively; Based on the asymmetric water flow transmission delay, a time-series delay bias mask is constructed. The multi-head attention association weight between the query vector and the key vector is calculated in combination with the time-series delay bias mask. The value vector is weighted and fused according to the multi-head attention association weight to generate candidate response representations related to the history of control commands. The representations are then fused with spatial features to generate the spatiotemporal feature matrix.

[0037] For the spatial feature vector, a matrix multiplication transformation is performed using a first fully connected linear mapping matrix to generate a basic query vector; its tensor dimension corresponds to the number of network nodes. The control instruction sequences for each historical time step are transformed by matrix multiplication using two independent sets of fully connected linear mapping matrices to generate the basic key vector and the basic value vector; their tensor dimensions correspond to the length of the historical time step. The basic query vector, basic key vector, and basic value vector are then uniformly split along the channel dimension and reshaped according to a set number of multi-head attention heads, resulting in multiple sets of parallel query sub-vectors. Key vector AND value subvector .

[0038] Based on the asymmetric water flow propagation delay, the cumulative water flow propagation delay from the pump station outlet node to all other nodes is calculated using the graph shortest path algorithm. The cumulative water flow transmission delay at all nodes is calculated. Calculate its standard deviation value, and use the reciprocal of the standard deviation value as the penalty coefficient. The cumulative water flow transmission delay at each node. Based on this, a time-series delay bias mask matrix is ​​constructed with rows representing the number of nodes and columns representing the length of historical time steps. In this mask matrix In the process, the cumulative water flow transmission delay at each node is calculated. The length of backtracking time represented by each historical time step The absolute difference is multiplied by the penalty coefficient. The negative value is then used as the mask matrix. The corresponding element value The calculation formula is as follows: For each attention head, obtain its corresponding subvector dimension. ; the query subvector With the corresponding key vector transpose matrix Perform matrix multiplication and divide by the dimension of the subvector. square root The basic attention score matrix with dimensions [number of nodes, historical time step length] for each head is obtained; the basic attention score matrix is ​​then compared with the temporal delay bias mask matrix. Matrix addition is performed, adding corresponding elements together, and physical transmission hysteresis correction is introduced; the resulting matrix is ​​then processed row-wise using a normalized exponential function. The calculation yields a multi-head attention correlation weight matrix incorporating physical time delay constraints. The calculation formula is as follows: For each attention head, the calculated multi-head attention association weight matrix is... With the corresponding value subvector matrix Matrix multiplication is performed, and the value vectors are weighted and summed to obtain the independently extracted hysteresis response feature components from each head. The weighted fusion outputs of all attention heads are concatenated along the channel dimension, and the feature dimensions are integrated through a linear output projection matrix to fully extract the dynamic hysteresis response feature components of the control command in the physical transmission process of water flow. The extracted dynamic hysteresis response feature components are fused with the spatial feature vectors of each node, and the weighted control command representation is fused with the node spatial features to generate a time-delay-constrained spatiotemporal feature matrix, which is then output as the spatiotemporal feature matrix.

[0039] A time-delay bias mask based on cumulative water flow transmission delay is constructed and superimposed on a multi-head attention score matrix. This forces the network to focus on historical control commands that conform to the actual physical transmission time, eliminating spurious time correlations. Candidate response representations related to the history of control commands are generated and fused with spatial features to achieve deep decoupling between control intervention and pipeline pressure evolution, eliminating command hysteresis interference, and outputting a pressure spatiotemporal feature matrix.

[0040] The process of outputting the global baseline pressure sequence and the local hysteresis fluctuation sequence includes: The spatiotemporal feature matrix is ​​subjected to feature projection and linear mapping in the node space dimension. The steady-state head distribution is extracted through the multilayer perceptron of the baseline network, and the global basic pressure sequence is output. The global basic pressure sequence is used as a conditional constraint vector, and combined with the spatiotemporal feature matrix, it is input at time steps as the gated recurrent unit network of the residual compensation network. Using the reset and update gates of the gated recurrent unit network, the hidden state of the gated recurrent unit network is updated in a time-series manner under the modulation of the condition constraint vector. Based on the physical information in the frictional physics prediction model, consistency penalty constraints are executed to fit the residual fluctuations caused by dynamic hysteresis and output the local hysteresis fluctuation sequence.

[0041] The spatiotemporal feature matrix is ​​obtained, and an initial fully connected linear mapping matrix is ​​constructed in the node space dimension. Matrix multiplication is performed on the spatiotemporal feature matrix using the fully connected linear mapping matrix to achieve feature projection and dimensional transformation, projecting the original spatiotemporal feature tensor to a feature space that matches the physical dimension of the actual pressure at the pipeline network nodes, resulting in a projected feature matrix. The baseline network in the friction physics prediction model is extracted. This baseline network is constructed using a multilayer perceptron consisting of an input layer, multiple hidden layers, and an output layer. The projected feature matrix is ​​fed into the multilayer perceptron as an input vector. Fully connected matrix multiplication is performed sequentially between each hidden layer of the multilayer perceptron, and the result of each matrix multiplication is input into a nonlinear activation function for activation processing. This process extracts features and performs nonlinear mapping on the long-term steady-state head distribution of the pipeline network under conditions that ignore transient fluid disturbances. The mapped feature values ​​are arranged in a one-dimensional sequence according to the time step sequence through the output layer of the multilayer perceptron, and the global baseline pressure sequence is calculated and output.

[0042] From the calculated global baseline pressure sequence, pressure numerical elements corresponding to each time step are extracted one by one and transformed into conditional constraint vectors. Following the forward order of time steps, spatiotemporal feature slices of the spatiotemporal feature matrix at the current time step are extracted one by one. The spatiotemporal feature slice corresponding to the current time step is then combined with the conditional constraint vector of the corresponding time step using a matrix concatenation and fusion operation along the feature channel dimension to construct a combined feature sequence containing global steady-state baseline constraint information. The generated combined feature sequence containing all time steps is used as the network input source data and is sequentially input into the residual compensation network within the tribophysical prediction model, following the time step order from history to the future. The underlying topology of the residual compensation network is the gated cyclic unit network.

[0043] For each time step of the sequence input, forward inference prediction is performed using the pre-trained gated recurrent unit network. An input weight matrix and a hidden state weight matrix are established. The combined feature sequence of the current time step is then... The hidden state at the previous time step The condition constraint vector is then multiplied and summed with the corresponding input weight matrix and hidden state weight matrix, respectively. The conditional modulation is directly introduced into the summation term of this matrix; then, the summation result is processed by a logistic activation function to calculate the reset gate vectors whose values ​​are in the interval between zero and one. With updating the gate vector Using the calculated reset gate vector The hidden state at the previous time step Perform Hadamard product operations to selectively forget historical time-series information; combine the forgotten hidden state with the combined feature sequence of the current time step. The candidate hidden states at the current time step are calculated by performing matrix addition and processing with the hyperbolic tangent activation function. Using the aforementioned update gate vector The hidden state of the previous time step With the candidate hidden state at the current time Perform a weighted linear combination operation based on the Hadamard product symbol. By performing element-wise multiplication and addition, the hidden state of the gated loop unit network is updated sequentially to obtain the hidden state at the current time step. Its update calculation formula is: The updated hidden state sequences at each time step are processed through a fully connected mapping layer to output a preliminary residual fluctuation sequence, which serves as the local hysteresis fluctuation sequence. The global base pressure sequence is then superimposed with the local hysteresis fluctuation sequence to obtain the preliminary pressure prediction sequence.

[0044] Before deployment, the aforementioned tribophysical prediction model requires a joint training phase based on physical information. During training, in addition to outputting the initial residual fluctuation sequence, the network also outputs corresponding auxiliary pipeline flow prediction values ​​through a parallel auxiliary fully connected layer. The pressure values ​​in the initial pressure prediction sequence, along with the fluid viscosity damping coefficient and equivalent elastic modulus of the pipeline material extracted from the physical asset parameter library of the tribophysical prediction model, are substituted into the continuity partial differential equation and momentum partial differential equation contained in the one-dimensional unsteady flow dynamics control equations. The finite difference algorithm is used to calculate the unbalanced residual values ​​generated by the spatial node pressure and pipeline flow under the condition of obeying physical conservation laws. The root mean square error of the unbalanced residual values ​​at all spatial nodes and time steps is calculated, and a consistency penalty loss function is constructed using this root mean square error scale. This consistency penalty loss function is added as a regularization penalty term to the total loss function of the gated recurrent unit network. A gradient descent optimization algorithm based on time backpropagation is employed to continuously iteratively update the parameters of various weight matrices within the network, aiming to minimize the total loss function. After the network parameters converge, the auxiliary fully connected layer is removed, and the network parameters are solidified to form the aforementioned feedforward gated recurrent unit network used for prediction.

[0045] A prediction structure combining basic pressure components and dynamic residual components is constructed. During the training phase, a loss term consisting of the discrete residuals of the continuity and momentum equations is added to the data fitting loss to constrain the hydraulic relationship between the predicted pressure and the auxiliary predicted flow rate. This mechanism enables the weight optimization of the residual network to use the equation residuals as soft constraints on the conservation of mass and momentum, eliminating the prediction distortion caused by purely data-driven models that violate objective laws during rapid changes in operating conditions, and endowing the output sequence with rigorous physical interpretability.

[0046] The process of converting to the aforementioned time compensation factor includes: Calculate the numerical difference and the first-order pressure change rate difference between the measured pressure data and the preliminary pressure prediction sequence within a preset sliding time window, and construct a gradient-weighted distance matrix; Under the preset path slope constraint, the cumulative distance is recursively calculated based on the gradient weighted distance matrix to determine the alignment path corresponding to the minimum cumulative distance; The average offset steps of the alignment path in the time axis direction are extracted as the hysteresis phase difference; The time compensation factor is generated by multiplying the hysteresis phase difference by the period.

[0047] The fluid dynamics parameters are retrieved and the typical attenuation time of the water hammer physical wave in the pump network system is calculated. A fixed window step size on the time axis is determined. The continuously sampled measured pressure data and the preliminary pressure prediction sequence are extracted using the window step size to construct the preset sliding time window. The variance of the measured pressure data and the variance of the first-order pressure change rate within the sliding time window are calculated. The gradient weight coefficient is calculated by dividing the pressure variance by the sum of the pressure variance and the variance of the first-order pressure change rate. Within this sliding time window, the measured pressure value is extracted for each time point of the measured pressure data. For each time point in the preliminary pressure prediction sequence, the predicted pressure value is extracted. The absolute difference between the two pressure values ​​is calculated as the numerical difference. Using a first-order forward finite difference algorithm, the slope of the pressure change between each adjacent time point is calculated for the measured pressure data. And the slope of the pressure change in the preliminary pressure prediction sequence between each adjacent time point. Calculate the absolute difference between the two pressure change slopes, and use this as the first-order pressure change rate difference. Multiply this numerical difference by the calculated gradient weight coefficient. In addition, the difference in the first-order pressure change rate is multiplied by one and the gradient weight coefficient is subtracted. The difference is used to obtain the weighted distance value for each time point combination. The calculation formula is as follows: Indexed by the time node of the measured pressure data As row coordinates, the time node index of the preliminary pressure prediction sequence is used. As column coordinates, all calculated weighted distance values Fill in the corresponding row and column coordinate grid to generate the two-dimensional gradient weighted distance matrix.

[0048] To prevent excessive distortion and unphysical mapping during time axis matching, the local travel slope for time node matching must be strictly constrained to the range of 1 / 2 to 2. Specifically, in the two-dimensional grid, the allowed recursive directions are limited to three states: one step forward horizontally and two steps forward vertically, one step forward vertically and two steps forward horizontally, and one step forward horizontally and one step forward vertically. This rule serves as the preset path slope constraint. After constructing the gradient-weighted distance matrix, a cumulative distance matrix with identical row and column dimensions is initialized. The starting element at the top left corner of the gradient-weighted distance matrix is ​​directly assigned to the starting element of the cumulative distance matrix, and the initial boundary rows and columns of the matrix are initialized using maxima to prevent out-of-bounds errors. Starting from the starting element, traverse all remaining coordinate grids of the cumulative distance matrix in the direction of increasing time. When calculating the cumulative distance of each coordinate grid, based on the preset path slope constraint, select the grid point with the smallest cumulative distance value from the three allowed predecessor adjacent coordinate grids. Add the minimum predecessor cumulative distance to the weighted distance value corresponding to the current coordinate grid in the gradient weighted distance matrix, and use the sum as the cumulative distance of the current coordinate grid. Simultaneously, record the routing pointer from this grid point to the minimum predecessor node in memory. After completing the recursive calculation of the entire matrix, start from the endpoint coordinates of the lower right corner of the cumulative distance matrix and perform a reverse backtracking search. In each backtracking step, directly read the routing pointer recorded by the current node and precisely backtrack to the previous optimal predecessor node along this pointer until backtracking to the starting coordinates of the upper left corner. Connect all the coordinate grid points traversed during the backtracking search in chronological order to generate the alignment path corresponding to the minimum cumulative distance from the starting point to the endpoint.

[0049] Extract the set of all coordinate grid points contained in the alignment path generated above. Each coordinate grid point is composed of a specific time index number of the measured pressure data and a specific time index number of the preliminary pressure prediction sequence. For each coordinate grid point in the alignment path set, subtract the time index number of the preliminary pressure prediction sequence from the time index number of the measured pressure data it contains to calculate the offset steps of each independent coordinate grid point in the time axis direction. If the subtraction result is positive, it indicates that the prediction sequence has a leading misalignment at the current node; if it is negative, it indicates that the prediction sequence has a lagging misalignment. Algebraically sum the offset steps corresponding to all coordinate grid points on the alignment path, and divide the sum by the total number of coordinate grid points contained in the alignment path to calculate the arithmetic mean of the offset steps. Extract this arithmetic mean as the hysteresis phase difference, which quantifies the overall physical transmission delay misalignment characteristics between the preliminary pressure prediction sequence and the measured pressure data.

[0050] The hysteresis phase difference is obtained, and its value is mathematically multiplied by the physical time interval of the period. Through this multiplication, the number of misaligned nodes in the discrete sequence obtained by the dynamic time warping algorithm is directly and linearly converted into the absolute time delay misalignment amount under a continuous physical time scale, and the time compensation factor with a specific physical time unit is output.

[0051] The measured pressure data used to calculate the hysteresis phase difference and time compensation factor are limited to the historical calibration window closed before the current time T. It is strictly forbidden to use the future true target value in the predicted interval t > T as input or calibration basis.

[0052] After dimensionlessly converting the pressure value difference and the first-order rate of change difference, a weighted distance is constructed, and the alignment path is calculated under the defined path stepping rules. Based on the signed average offset of the alignment path within the closed historical window, an empirical time offset for subsequent prediction sequence time calibration is obtained, providing rigorous and reliable parameter support for subsequent sequence calibration.

[0053] The time axis of the preliminary pressure prediction sequence is shifted and calibrated, including: Based on the time compensation factor, the time axis of the preliminary pressure prediction sequence is shifted along the time feedforward direction; for the end prediction gaps generated after the time axis phase shift, the boundary trend extrapolation is performed using the global basic pressure sequence to fill the gaps. On the shifted time axis after filling, the preliminary pressure prediction sequence after extrapolation and filling is resampled using the conformal cubic Hermitian interpolation algorithm. This process ensures that the interpolation curve passes through the original node values ​​and maintains monotonicity within the monotonic interval of the original data, reducing additional overshoot, and outputs the pressure prediction sequence.

[0054] The preliminary pressure prediction sequence and its corresponding original discrete time axis are obtained, and the time compensation factor with specific physical time units is read. Due to the hysteresis effect of fluid transmission and resistance in the physical pipeline network, the absolute physical time value of the time compensation factor is uniformly added as an addend to each timestamp node of the original discrete time axis. By performing the timestamp addition operation on the absolute physical time value of the time compensation factor, all data points of the preliminary pressure prediction sequence are shifted by the same physical time span in the time dimension towards the future time feedforward direction, constructing a shifted time axis containing physical time delay characteristics, thereby achieving phase shift calibration of the preliminary prediction sequence on the time axis.

[0055] Because the initial pressure prediction sequence is shifted forward along the time axis, a time gap with a length equal to the time compensation factor is generated at the very end of the fixed prediction time window required for output; this is the end prediction gap. The last valid data point in the shifted initial pressure prediction sequence is extracted and used as the boundary starting point. Data slices of the global base pressure sequence within this time gap are retrieved according to the corresponding time coordinates. The pressure value difference at each adjacent time node of the data slice is calculated using a differential algorithm. This difference is divided by the time step to obtain the pressure change slope characterizing the future basic flow field direction of the system. This pressure change slope is used as the boundary trend. Using the pressure value at the boundary starting point as the reference starting point value, the algebraic accumulation calculation of the pressure change difference is performed successively at each discrete time node of the end prediction gap according to the extracted boundary trend, generating an extrapolated filling data sequence. The generated extrapolated data sequence is directly spliced ​​to the end of the translated preliminary pressure prediction sequence to fill the prediction gap, so that the trend of the filled sequence strictly follows the physical benchmark law of pipeline network regressing to steady state.

[0056] Obtain the standard time grid defined by the underlying system, which is divided into equal intervals with a fixed sampling period. For any target time point in this standard time grid, find the adjacent known time nodes that surround it and their corresponding pressure filling values ​​on the filled translational time axis. Calculate the first-order forward and backward difference slopes between the adjacent known time nodes. To prevent non-physical numerical oscillations between data points caused by traditional polynomial interpolation, the difference slopes of adjacent intervals are harmonicly averaged to obtain the initial node derivatives. The initial node derivatives are then truncated using a monotonicity criterion. Specifically, when the data in two adjacent intervals show opposite trends, the derivatives of the intersecting nodes are forcibly set to zero. When they show the same trend, the absolute value of the node derivatives is limited to no more than three times the slope of the adjacent tangent angle. This truncation operation ensures that the interpolation curve is strictly constrained within the local first-order pressure change rate formed by the original data points, and that no new local extreme points are artificially generated in intervals where the original data sequence has not experienced extreme value reversal. By using the nodal derivatives after amplitude limiting and the corresponding known nodal pressure values, the polynomial coefficients are solved to construct the piecewise conformal cubic Hermitian polynomial basis function equations. The specific values ​​at each target time point in the standard time grid are substituted one by one into the cubic Hermitian polynomial basis function equations for the corresponding time intervals for calculation, yielding a newly sampled pressure value sequence. After completing the calculations for all standard time grid nodes, a final prediction sequence aligned with the system's standard clock and with a smooth waveform is generated and output as the pressure prediction sequence.

[0057] Physical transmission hysteresis is compensated for by time-axis feedforward shift, and the end gap is filled by extrapolating the steady-state baseline trend, ensuring that the tail trend of the predicted sequence conforms to the objective physical law of pipeline network regression to steady state. In the resampling stage, a conformal cubic Hermitian interpolation and derivative amplitude limiting mechanism are introduced to eliminate non-physical numerical oscillations. This mechanism maintains the local extrema and first-order pressure change rate of the original data, outputting a pressure prediction sequence with a smooth waveform and aligned with the system's standard clock.

[0058] This system deeply integrates fluid dynamics mechanisms into neural networks. It extracts features and decouples command lag through asymmetric time-delay graph convolution and bias masking; it introduces a one-dimensional unsteady flow equation to construct a physical consistency penalty, forcing the network to obey fluid conservation laws and eliminating prediction distortions caused by purely data-driven approaches that violate objective laws; and it combines gradient-weighted time warping and conformal interpolation to eliminate time delay misalignment and numerical oscillations. This endows the predicted sequences with rigorous physical interpretability, ensuring the reliability of dynamic pipeline network control.

[0059] Example 2 This embodiment will implement the dynamic pressure prediction method for water pump pressure controllers based on digital twins in a specific scenario, taking a municipal water supply network system as an example.

[0060] By retrieving the two-dimensional design drawings and pipeline asset ledger data of the municipal water pump network system, the specific physical length of each pipe section, the pipe inner diameter (0.4 meters), the ductile iron pipe material, and the spatial three-dimensional coordinates of each pump interface, valve location, and pipe intersection point within the network are extracted. In system memory, the pump interfaces, valve locations, and pipe intersection points corresponding to the above spatial three-dimensional coordinates are mapped as abstract nodes in a data structure, and the physical pipes connecting adjacent nodes are mapped as edges in the data structure. A mapping table containing node indexes and edge relationships is established to generate a digital twin mapping model of the water pump network system. Based on the digital twin mapping model, the connection status between each node in the mapping table is read. If there is a real physical pipe connection between two nodes, the topological connectivity variable between them is marked as 1; otherwise, it is marked as 0. This process parses the physical topological relationships of the nodes in the entire network.

[0061] Industrial-grade pressure transmitters and electromagnetic flowmeters are deployed in the pipeline network. Sampling is performed at fixed intervals of 2 seconds to simultaneously collect pipeline pressure and flow rate values ​​at two adjacent nodes. The absolute difference between the pipeline pressure values ​​at two adjacent nodes is calculated to determine the node pressure difference.

[0062] The water flow rate is divided by the pipe inner diameter cross-sectional area extracted from the digital twin mapping model to calculate the water flow velocity in each pipe. The actual water flow direction is determined by the sign of the water flow rate value. The actual water flow direction is defined as the flow direction, and the opposite direction is defined as the counter-flow direction. The basic wave velocity of sound waves in the water medium under the current water temperature is obtained by consulting a water acoustics physics handbook. For the flow direction, the basic wave velocity is added to the water flow velocity to obtain the downstream absolute wave velocity; for the counter-flow direction, the basic wave velocity is subtracted from the water flow velocity to obtain the counter-flow absolute wave velocity. Based on the calibrated pressure wave propagation speed (i.e., the downstream and counter-flow absolute wave velocities) and the physical length of the pipe section, the direction-dependent propagation delay is obtained. Since the downstream and counter-flow absolute wave velocities are different, unequal downstream and counter-flow propagation delays are obtained, which together constitute the asymmetric water flow propagation delay.

[0063] The asymmetric water flow transmission delay calculated between all adjacent nodes in the pipe network within the current acquisition period (i.e., the sampling period). Calculate the arithmetic mean to obtain the global average delay scale. Using the natural constant Using the base as the time delay for the asymmetric water flow propagation between any adjacent nodes. Divided by the global average delay scale The negative value of the quotient is used as the exponent to construct the exponential decay kernel function. This function is then used to calculate the asymmetric water flow propagation delay between adjacent nodes. By inputting the exponential decay kernel function one by one and performing natural exponential decay calculation, the dynamic spatial decay weight between each adjacent node with a value range between 0 and 1 is calculated. The calculation formula is as follows: Initialize a square matrix with the number of rows and columns equal to the total number of nodes in the pipeline network. Traverse this matrix according to the physical topology of the nodes. If there is no connection between nodes, assign a value of zero; if there is a connection, fill the corresponding dynamic spatial attenuation weight into the matrix elements according to the direction of water flow (either the direction of flow or the direction of counterflow). After completing the traversal and assignment, a directed weighted adjacency matrix of the pump pipeline network system with asymmetric numerical characteristics is constructed, which serves as the spatial topology. The above steps are repeated in each subsequent sampling period to achieve updates.

[0064] The outlet node of the municipal pumping station's main outlet pipe was selected, and the terminal interfaces of branch pipes in each residential area were selected as the end nodes of the pipe network. Pipe pressure, flow rate, and valve opening were synchronously collected at each node according to the specified sampling period. The pressure change rate was calculated by dividing the pipe pressure difference between adjacent periods by the time interval; a non-zero minimum constant was also used. To prevent overflow, the pipe pressure value should be adjusted. Divide by water flow rate Square and the constant The sum of these values ​​is used to calculate the dynamic gradient of water resistance. The calculation formula is as follows: The above data combination constitutes the complete time-series sensing data for each node. Simultaneously, the frequency modulation pulse duty cycle change rate, pump start / stop status, and millisecond-level command issuance timestamps within adjacent adjustment cycles are obtained through the controller interface and combined to form the control command sequence. Using the command issuance timestamps as the time alignment reference axis, a linear interpolation function is applied to the time-series sensing data of each node for time-series interpolation alignment; subsequently, a historical sliding window with a length of 120 time steps is extracted, and the maximum values ​​of each data item within the window are statistically analyzed. and minimum value The range normalization formula is used to normalize all the original data within the window. The normalized value is obtained by linearly scaling to a closed interval between zero and one. The calculation formula is as follows: Complete the normalization process for the sliding window.

[0065] The normalized time-series sensing data of each node is used as node features, and the dynamic spatial attenuation weights are extracted. The control command sequence is input into a multilayer perceptron for nonlinear mapping to generate conditional feature vectors, which are then broadcast and concatenated into the time-series sensing data of each node to complete the conditional feature embedding. The multilayer perceptron outputs a feature modulation weight matrix, which is scaled and calculated with the dynamic spatial attenuation weights to construct the dynamic convolution kernel weights of the edge conditional graph convolutional layer. Multi-order neighborhood message passing is performed along the water flow and wave propagation directions defined by the directed weighted adjacency matrix; features converging in the downstream direction are assigned active aggregation weights, and features converging in the upstream direction are assigned damping aggregation weights, and asymmetric spatial weighted summation is performed; after concatenating the node's own features and activating them, the spatial features containing the head loss and dynamic hysteresis along the pipeline network are output.

[0066] The spatial features are linearly transformed into a query vector. The control command sequences at each historical time step are linearly transformed into key vectors. AND value vector The cumulative water flow transmission delay from the pump station outlet to the remaining nodes is calculated using the shortest path algorithm. Compare it with the length of historical time steps The absolute difference is multiplied by a penalty coefficient of 1.5, which represents the physical deviation, and then negative to construct a timing delay bias mask. The calculation formula is as follows: Calculate query vector With key vector The base attention score is added to the temporal delay bias mask M, and an exponential normalization function is performed row by row. The multi-head attention association weight A is then obtained (assuming the sub-vector dimension is ). The calculation formula is as follows: The value vector is weighted and fused using this weight to extract the hysteresis response component generated by the control command sequence. The spatial features of each node are subtracted from the hysteresis response component to generate a candidate response representation related to the control command history. This candidate response representation is then fused with the spatial features and reconstructed dimensionally to output the spatiotemporal feature matrix.

[0067] The spatiotemporal feature matrix is ​​projected onto the node space dimension and input into the baseline network (multilayer perceptron) of the friction physics prediction model to extract the steady-state head distribution of the pipeline network and output the global base pressure sequence. This sequence is used as a conditional constraint vector and, combined with the spatiotemporal feature matrix, is input into the residual compensation network (gated recurrent unit network) step by step. The hidden state of the previous time step is utilized. The reset gate is calculated by combining the current input features with the modulation of the conditional constraint vector. With the update gate ; Calculate candidate hidden states using reset gates The hidden state is then updated sequentially using an update gate to obtain the current hidden state. Its update calculation formula is: The local hysteresis fluctuation sequence caused by dynamic hysteresis is fitted. The global base pressure sequence is superimposed with the local hysteresis fluctuation sequence to obtain the preliminary pressure prediction sequence.

[0068] Based on the typical attenuation time of the physical wave of water hammer in the municipal water pump network system, a fixed sliding window step size is extracted, and the measured pressure data and the preliminary pressure prediction sequence are extracted; the node values ​​of the measured pressure data are extracted. Compared with the node values ​​of the preliminary pressure prediction sequence The numerical difference between the two is calculated; simultaneously, the measured first-order pressure change rate is calculated using a difference algorithm. Compared with the predicted first-order rate of change of pressure The difference in the first-order pressure change rate between the two is calculated; using a weighting factor of 0.5 (the balance steady-state error and transient waveform matching weight) as the gradient weight coefficient, the two are weighted and fused to calculate the weighted distance value for each time node combination. The calculation formula is as follows: A two-dimensional gradient-weighted distance matrix is ​​constructed accordingly. To prevent the time-axis warping algorithm from producing non-physically distorted mappings, the recursive direction is strictly constrained to a path slope between 1 / 2 and 2. Based on the gradient-weighted distance matrix, the cumulative distance is recursively calculated, and the alignment path corresponding to the minimum cumulative distance is generated by backtracking. The average offset steps of this alignment path in the time axis direction are extracted as the hysteresis phase difference, and multiplied by the 2-second periodic physical time interval to generate the time compensation factor with absolute physical time units.

[0069] Based on the time compensation factor, the discrete time axis of the preliminary pressure prediction sequence is shifted along the time feedforward direction. For the end-prediction gap generated at the end of the prediction window after the shift, data slices corresponding to the time period of the global base pressure sequence are extracted, and their pressure change slopes are calculated as boundary trends. These trends are then extrapolated and accumulated sequentially, starting from the last valid point of the shifted sequence, to fill the boundary trend extrapolation. On the filled shifted time axis, using the standard time grid of the underlying system as the target, the differential slopes of adjacent known nodes are calculated and derivative amplitude-limited. Piecewise conformal cubic Hermitian polynomial basis function equations are constructed using the amplitude-limited node derivatives, and the sequence is resampled. Under the premise of ensuring the interpolation curve passes through the original node values, maintaining monotonicity within the monotonic interval of the original data, and reducing additional overshoot, the final output is a pressure prediction sequence with a smooth and aligned waveform.

[0070] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A dynamic pressure prediction method for a water pump pressure controller based on digital twins, characterized in that, include: A digital twin mapping model of the water pump pipeline network system is constructed, and the spatial topology is extracted based on the digital twin mapping model; the time-series sensing data and measured pressure data of each node are collected periodically to obtain the control command sequence of the water pump pressure controller; The time-series sensing data and the control command sequence are input into a graph neural network, and the spatial features of nonlinear hysteresis of node pressure are extracted by combining the spatial topology. The spatial features are decoupled using a temporal attention mechanism to output a spatiotemporal feature matrix. The spatiotemporal feature matrix is ​​then input into a tribophysical prediction model. The global base pressure sequence and the local hysteresis fluctuation sequence are output through the baseline network and the residual compensation network in the tribophysical prediction model, respectively. The global base pressure sequence and the local hysteresis fluctuation sequence are superimposed to obtain a preliminary pressure prediction sequence. The measured pressure data is compared with the preliminary pressure prediction sequence using a dynamic time warping algorithm. The hysteresis phase difference is extracted and converted into a time compensation factor. The time axis of the preliminary pressure prediction sequence is shifted and calibrated using the time compensation factor, and the pressure prediction sequence is output.

2. The dynamic pressure prediction method for a water pump pressure controller based on digital twin as described in claim 1, characterized in that, The process of extracting the spatial topology includes: Based on the digital twin mapping model, the physical topology of the pipeline network nodes is analyzed. Combining the water flow velocity, friction coefficient, and pressure difference of each pipeline, the asymmetric water flow transmission delay between adjacent nodes along the flow direction and the counterflow direction is dynamically calculated. The asymmetric water flow transmission delay is input into the exponential decay kernel function to calculate the dynamic spatial decay weight between each node. Based on the dynamic spatial decay weights and the physical topology relationship of the nodes, a directed weighted adjacency matrix of the water pump network system is constructed and updated as the spatial topology.

3. The dynamic pressure prediction method for a water pump pressure controller based on digital twin according to claim 2, characterized in that, Each node includes the pump station outlet and the end node of the pipeline network. The time-series sensing data includes the pipeline pressure, water flow rate, and valve opening collected by each node according to the cycle, as well as the pressure change rate and water resistance dynamic gradient calculated based on the pipeline pressure and the water flow rate. The control command sequence includes the frequency modulation pulse duty cycle change rate issued by the water pump pressure controller, the water pump start-up and shut-down status, and the command issuance timestamp; before inputting the time-series sensing data and the control command sequence into the graph neural network, the time-series sensing data of each node is subjected to time-series interpolation alignment and sliding window normalization based on the command issuance timestamp.

4. The dynamic pressure prediction method for a water pump pressure controller based on digital twin according to claim 3, characterized in that, The process of extracting the spatial features includes: Using an edge-conditional graph convolutional layer based on the directed weighted adjacency matrix, the control command sequence is embedded as a conditional feature into the temporal sensing data of each node; along the water flow and wave propagation directions defined by the directed weighted adjacency matrix, multi-order neighborhood message passing and asymmetric spatial aggregation are performed on the temporal sensing data of each node that has completed the embedding of the conditional features to obtain the spatial features.

5. The dynamic pressure prediction method for a water pump pressure controller based on digital twin according to claim 4, characterized in that, The process of decoupling the spatial features from the nodal pressure nonlinear hysteresis features includes: The spatial feature vector is linearly transformed to obtain the query vector, and the control command sequence of each historical time step is linearly transformed to obtain the key vector and the value vector respectively; Based on the asymmetric water flow transmission delay, a time-series delay bias mask is constructed. The multi-head attention association weight between the query vector and the key vector is calculated in combination with the time-series delay bias mask. The value vector is weighted and fused according to the multi-head attention association weight to generate candidate response representations related to the history of control commands. The representations are then fused with spatial features to generate the spatiotemporal feature matrix.

6. The dynamic pressure prediction method for a water pump pressure controller based on digital twin according to claim 1, characterized in that, The process of outputting the global baseline pressure sequence and the local hysteresis fluctuation sequence includes: The spatiotemporal feature matrix is ​​projected and linearly mapped in the node space dimension. The steady-state head distribution is extracted through the multilayer perceptron of the baseline network, and the global base pressure sequence is output. The global base pressure sequence is used as a conditional constraint vector and combined with the spatiotemporal feature matrix as a time-step input as the gated recurrent unit network of the residual compensation network. Using the reset gate and update gate of the gated recurrent unit network, the hidden state of the gated recurrent unit network is updated in a time-series manner under the modulation of the conditional constraint vector. Based on the physical information in the friction physics prediction model, a consistency penalty constraint is executed to fit the residual fluctuation caused by dynamic hysteresis and output the local hysteresis fluctuation sequence.

7. The dynamic pressure prediction method for a water pump pressure controller based on digital twin according to claim 1, characterized in that, The process of converting to the aforementioned time compensation factor includes: Calculate the numerical difference and the first-order pressure change rate difference between the measured pressure data and the preliminary pressure prediction sequence within a preset sliding time window, and construct a gradient-weighted distance matrix; Under the preset path slope constraint, the cumulative distance is recursively calculated based on the gradient weighted distance matrix to determine the alignment path corresponding to the minimum cumulative distance; the average offset step of the alignment path in the time axis direction is extracted as the hysteresis phase difference; the hysteresis phase difference is multiplied by the period to generate the time compensation factor.

8. The dynamic pressure prediction method for a water pump pressure controller based on digital twin according to claim 7, characterized in that, The time axis of the preliminary pressure prediction sequence is shifted and calibrated, including: Based on the time compensation factor, the time axis of the preliminary pressure prediction sequence is shifted along the time feedforward direction; for the end prediction gaps generated after the time axis phase shift, the boundary trend extrapolation is performed using the global basic pressure sequence to fill the gaps. On the shifted time axis after filling, the preliminary pressure prediction sequence after extrapolation and filling is resampled using the conformal cubic Hermitian interpolation algorithm. This process ensures that the interpolation curve passes through the original node values ​​and maintains monotonicity within the monotonic interval of the original data, reducing additional overshoot, and outputs the pressure prediction sequence.

Images