Intelligent optimization method for water hammer protection based on adaptive rabbit optimization non-constant friction pinn-lmoea-ds

Abstract

The application discloses a kind of PINN-LMOEA-DS water hammer protection intelligent optimization methods based on adaptive rabbit optimization non-constant friction;The method first establishes water transfer project hydraulic transient model, introduces non-constant friction model, and adopts the improved group intelligence optimization algorithm introducing cross-dimension information interaction mechanism to calibrate friction parameter, and obtains optimal non-constant friction parameter;On this basis, construct physical information neural network proxy model, embed the loss function as physical constraint with hydraulic transient control equation, and the optimal non-constant friction parameter is substituted into the calculation of physical constraint as known constant to realize decoupling coupling;Water hammer protection parameters are optimized using large-scale multi-objective evolutionary algorithm based on directional sampling, and a set of Pareto optimal solutions are obtained;Finally, the optimization results are substituted into the hydraulic transient model for verification and closed-loop iteration.The application realizes the deep coupling of non-constant friction parameter calibration, physical constraint agent modeling and multi-objective optimization, and improves the precision and efficiency of water hammer protection design.

Description

Technical Field

[0001] This invention relates to the field of hydraulic transient protection and control technology for water diversion projects, and in particular to the PINN-LMOEA-DS intelligent optimization method for water hammer protection based on adaptive rabbit optimization of non-constant friction. Background Technology

[0002] Water hammer is one of the most threatening transient phenomena in pressurized water transmission systems. Unexpected pump shutdowns can cause severe pressure fluctuations, leading to extreme water hammer pressure and vaporization in the pipeline. In extreme cases, this can even cause serious accidents such as pipeline rupture and pump damage. To ensure project safety, water hammer protection devices such as bidirectional pressure regulating towers and post-pump valves are typically installed, but this significantly increases project investment. Therefore, finding a balance between safety and cost-effectiveness, and resolving this contradiction through optimization algorithms, is particularly important.

[0003] Numerical calculation of water hammer is the cornerstone of water hammer protection optimization research. Traditional water hammer calculations generally adopt the assumption of constant friction, but under high-speed transient flow conditions, the constant friction model leads to significant errors between the calculated results and measured data. To address this, the academic community has proposed non-constant friction models such as Brunone's, but these models contain empirical coefficients that need to be calibrated. Traditional calibration methods often rely on trial and error or simple heuristic algorithms, which are prone to getting trapped in local optima in high-dimensional, strongly nonlinear parameter spaces.

[0004] Regarding surrogate models, purely data-driven models lack physical constraints and are prone to violating fluid dynamics laws when data is sparse. Physical Information Neural Networks (PINNs) ensure physical consistency by embedding the governing equations into the loss function. However, existing PINN water hammer prediction studies (such as Xu et al.) typically include non-constant friction coefficients as unknown parameters to be inverted in the loss function for joint solution. When dealing with high-frequency oscillating water hammer transients, this joint inversion method leads to an extremely non-convex optimization space for the loss function, making the network highly susceptible to getting trapped in local optima or even failing to converge.

[0005] In terms of optimization algorithms, water hammer protection optimization involves pump and valve closing rules, pressure regulating tower structural parameters, etc., exhibiting high-dimensionality, multi-objective, and highly nonlinear characteristics. Traditional multi-objective evolutionary algorithms face bottlenecks such as slow convergence speed and diversity degradation when dealing with high-dimensional decision variables. The Large-Scale Multi-Objective Evolutionary Algorithm Based on Directed Sampling (LMOEA-DS) performs excellently in handling high-dimensional problems, but still lacks system integration with high-precision non-constant friction models and surrogate models.

[0006] In terms of decision-making methods, existing technologies mostly adopt the entropy weight-fuzzy TOPSIS method, which is highly subjective and has cumbersome consistency checks.

[0007] Therefore, how to overcome the drawback of non-convergence of non-constant friction parameters in PINN joint inversion and construct a complete closed-loop optimization system from high-precision water hammer simulation to optimal protection scheme decision-making is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0008] The purpose of this invention is to provide a PINN-LMOEA-DS intelligent optimization method for water hammer protection based on adaptive rabbit optimization of non-constant friction, to solve the problems in existing technologies such as the ease with which non-constant friction parameter calibration gets trapped in local optima, the difficulty in convergence of PINN joint inversion of friction parameters, insufficient high-dimensional optimization search capability, and strong subjectivity in the decision-making process. To solve the above technical problems, this invention provides the following technical solution:

[0009] This invention provides an intelligent optimization method for water hammer protection parameters, comprising the following steps:

[0010] S1. Establish a hydraulic transient model for the water transfer project, introduce a non-constant friction model to describe the transient flow process, and use the friction parameter as the optimization variable. An improved swarm intelligence optimization algorithm with a cross-dimensional information interaction mechanism is used to globally calibrate the friction parameter, obtaining the optimal non-constant friction parameter k. opt ;

[0011] S2. Construct a physical information neural network proxy model, embedding the hydraulic transient control equations as physical constraints into the network loss function, and using the k obtained in step S1. opt Substituting known constants into the calculation of non-constant friction terms in physical constraints reduces the dimensionality of network optimization and locks in the energy dissipation characteristics of transient flow, thereby enabling the prediction of water hammer characteristic quantities.

[0012] S3. Using the physical information neural network proxy model constructed in step S2 as the objective function, establish a multi-objective optimization model, and use an evolutionary algorithm with a directional sampling mechanism to optimize and solve the high-dimensional water hammer protection parameters to obtain the Pareto optimal solution set.

[0013] S4. Perform multi-index decision analysis on the Pareto optimal solution set obtained in step S3 to determine the optimal combination of water hammer protection parameters.

[0014] S5. Substitute the optimal parameters obtained in step S4 into the hydraulic transient model in step S1 for verification. If the deviation between the verification result and the predicted value of the surrogate model exceeds the set threshold, then supplement the verification data into the surrogate model training set to trigger closed-loop iterative optimization.

[0015] Preferably, in step S1, the expression for the friction term in the non-constant friction model is: ,in For a constant coefficient of friction, The non-constant friction correction term is calculated using the following formula:

[0016] ;

[0017] In the formula, k is the Brunone friction coefficient (i.e., frictional resistance parameter), g is the acceleration due to gravity, D is the pipe diameter, a is the wave velocity, V is the flow velocity, t is time, and x is the spatial coordinate. The improved swarm intelligence optimization algorithm obtains the optimal value of k through search. opt .

[0018] Preferably, in step S2, the loss function L of the physical information neural network includes a data error term, a physical residual term, and a boundary condition loss term, expressed as: ,

[0019] In the formula,

[0020] , ,

[0021] in and The corresponding residuals for the continuity equation and the momentum equation are as follows:

[0022] ;

[0023] ;

[0024] Where L: total loss function of Physical Information Neural Network (PINN); h: pressure head; v: flow velocity; t: time variable; x: spatial coordinate variable along the pipe axis; a: water hammer wave velocity; g: gravitational acceleration; f: Darcy-Weisbach constant friction coefficient; D: pipe diameter; kopt: Brunone friction coefficient.

[0025] Preferably, in step S1, the improved swarm intelligence optimization algorithm introduces a cross-dimensional information interaction mechanism during the population update process to escape local optima in a highly nonlinear friction parameter space; the cross-dimensional information interaction mechanism includes vertical cross search, that is, randomly selecting two different individuals and exchanging their variable values ​​in a certain dimension.

[0026] Preferably, in step S3, the evolutionary algorithm with a directional sampling mechanism executes a directional sampling strategy in the early stages of iteration, including convergence direction sampling and diversity direction sampling; the convergence direction sampling uses simultaneous perturbation stochastic approximation to estimate the approximate gradient, and the calculation formula is:

[0027] ,

[0028] Where Δ is the perturbation vector following a Bernoulli distribution, and c is the perturbation step size; diversity direction sampling generates sampled solutions by constructing orthogonal directions representing the solutions, and dynamically adjusts the ambiguity of the generated sampled solutions using fuzzy variable operators, the calculation formula being:

[0029] Where δ is the fuzzy coefficient;

[0030] U j : The upper bound of the j-th optimization variable; L j : The lower bound of the j-th optimization variable;

[0031] X new,j : The value of the newly generated sampled solution in the j-th dimension; N(0,1): Standard normal distribution, used to generate random perturbation terms, so that the sampled solution has a certain degree of randomness and exploratory nature in the vicinity of the representative solution.

[0032] Preferably, in step S4, the multi-index decision analysis method uses a weight allocation method based on the comparison of best and worst criteria to determine the index weights, and obtains the optimal criterion weight vector by solving the following minimax optimization problem. ;

[0033]

[0034]

[0035] ;

[0036] In the formula, The weights of the optimal criterion Let ξ be the weight of the worst criterion, n be the total number of criteria, and ξ be the maximum absolute deviation. The optimal weight vector is obtained by minimizing ξ.

[0037] Calculate the consistency ratio (CR):

[0038] In the formula, To achieve the optimal target value, As a consistency index, if CR < 0.1

[0039] Then, the weight vector is accepted, and the optimal solution is selected by combining it with a sorting method based on the distance to the ideal solution.

[0040] Preferably, the optimal criterion C is dynamically adjusted based on engineering requirements. B And worst-case criterion C W This enables adaptive solution decision-making under scenarios prioritizing safety and reliability, cost-effectiveness, or unbiased optimization.

[0041] Preferably, the water hammer protection parameters include the operating parameters of the downstream valve and the structural parameters of the pressure regulating device; the operating parameters of the downstream valve include the fast closing time, fast closing angle, slow closing time, and slow closing angle; the structural parameters of the pressure regulating device include the tower diameter, the initial effective head, and the diameter of the water supply pipe.

[0042] The beneficial effects of this invention are: (1) Improved efficiency and accuracy of non-constant friction parameter calibration: The improved swarm intelligence algorithm with the introduction of vertical crossover and adaptive population partitioning has a stronger global optimization ability in the highly nonlinear friction calibration problem and can effectively escape local optima.

[0043] (2) The physical consistency and robustness of the prediction model are significantly enhanced: It breaks through the technical bottleneck of the existing PINN joint inversion of friction parameters which is prone to non-convergence, and pioneers the decoupled coupling architecture of "pre-calibration solidification + kopt known constant substitution". The optimized and calibrated non-constant friction terms are embedded into the neural network in the form of hard constraints, which reduces the network optimization dimension and locks the energy dissipation characteristics, so that high-precision convergence prediction can still be maintained even under sparse data.

[0044] (3) Significantly improved high-dimensional search capability: The constructed PINN proxy model has both physical smoothness and extremely fast reasoning capability, which makes high-dimensional directional sampling of large-scale multi-objective evolutionary algorithm (LMOEA-DS) possible. The convergence speed and solution set distribution uniformity are significantly better than traditional algorithms.

[0045] (4) The decision-making process is scientific, objective, and robust with a closed loop: The use of dynamic preference decision-making reduces the workload of subjective judgment; at the same time, the high-precision numerical model verification and data feedback mechanism constitute a true adaptive closed-loop correction system, improving engineering practicality and system robustness. Attached Figure Description

[0046] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0047] Figure 1 This is a flowchart illustrating the overall process of the method of the present invention.

[0048] Figure 2 Flowchart for the Improved Swarm Intelligence Optimization Algorithm (ART-VCS);

[0049] Figure 3 A schematic diagram of the Physical Information Neural Network (PINN) surrogate model structure;

[0050] Figure 4Optimize the flowchart for the Large-Scale Multi-Objective Evolutionary Algorithm (LMOEA-DS);

[0051] Figure 5 This is a flowchart for multi-indicator decision analysis. Detailed Implementation

[0052] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0053] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0054] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0055] Reference Figures 1-5 This is one embodiment of the present invention, which provides a PINN-LMOEA-DS intelligent optimization method for water hammer protection based on adaptive rabbit optimization of non-constant friction, comprising the following steps:

[0056] S1. Establish a hydraulic transient model for the water transfer project, introduce a non-constant friction model to describe the transient flow process, and use the friction parameter as the optimization variable. An improved swarm intelligence optimization algorithm with a cross-dimensional information interaction mechanism is used to globally calibrate the friction parameter, obtaining the optimal non-constant friction parameter k. opt .

[0057] The specific method of step S1 includes:

[0058] S1.1 Based on the topology diagram of the water diversion project, a hydraulic transient simulation model is established using the method of characteristics. The boundary conditions include: upstream boundary, downstream boundary, downstream valve boundary, air valve boundary, unidirectional pressure regulating tower boundary, and pump station boundary.

[0059] S1.2 Collect and organize monitoring data on flow rate, pressure, water level, and operating status of pumps, downstream valves, and unidirectional pressure regulating towers along the pipeline. Select two sets of monitoring data containing 10 minutes of continuous monitoring data with a step size of 0.01 seconds under actual operating conditions. Use these data to calibrate and verify the model, controlling the verification error to within 10%, and obtain a hydraulic transient simulation model that conforms to the actual operation of the water diversion project.

[0060] S1.3. Using Brunone's non-constant friction model, the expression for its friction term is:

[0061] (1)

[0062] In the formula, The constant friction coefficient is calculated using the Darcy-Weisbach formula; The correction term for non-constant friction is given by the following formula:

[0063] (2)

[0064] In the formula, denoted as instantaneous acceleration, 'a' as wave velocity, and 'k' as Brunone friction coefficient, which are key empirical parameters that need to be calibrated through optimization algorithms.

[0065] Using the empirical coefficient k as the optimization variable, the adaptive rabbit optimization algorithm combined with vertical cross search (ART-VCS) is applied to the range of values ​​of k. The internal search solution aims to minimize the root mean square error between the calculated and measured pressure heads, thereby obtaining the optimal non-constant friction coefficient k. opt ;

[0066] S1.4 Select the parameters of the pump downstream valve closing behavior and the parameters of the unidirectional pressure regulating tower as the target optimization variables; then obtain the input parameter set of the hydraulic transient simulation model, i.e., the target optimization variable set, through the Latin hypercube sampling method, and substitute them into the embedding k opt A high-precision numerical calculation model for water hammer is used to simulate the water hammer, and the maximum water hammer pressure value is obtained. Minimum negative pressure value and the maximum reverse rotation speed of the water pump Construct the desired sample set.

[0067] Furthermore, in S1.3, the specific steps of the ART-VCS algorithm for optimizing the non-constant friction coefficient include:

[0068] S1.3.1 Initialize the rabbit population size N and the maximum number of iterations. Adaptive partitioning scaling factor Vertical crossover probability p vcs ,exist The positions of N rabbits are randomly generated within the range, and the position of each rabbit represents a candidate empirical coefficient k;

[0069] S1.3.2 Substitute the k value corresponding to each rabbit into the Brunone non-constant friction model, calculate the water hammer pressure response using the method of characteristics, and calculate the root mean square error (RMSE) between the calculated and measured pressure head as the fitness value.

[0070] S1.3.3 Adaptively divide the population based on fitness value and current iteration progress: Sort all individuals in the current population according to their fitness value (RMSE) from smallest to largest (the smaller the RMSE, the better the fitness). The top α×100% of individuals (i.e., the top α proportion of individuals with the best fitness) form the development subgroup, and the remaining individuals (i.e., the bottom (1-α)×100% of individuals) form the exploration subgroup. Here, α is the adaptive partitioning ratio factor, which takes a value in the range of (0,1) and is dynamically adjusted with the number of iterations: In the early stage of the algorithm iteration, α takes a smaller value (e.g., 0.3) to allow more individuals to participate in global exploration and expand the search range; in the later stage of the iteration, α takes a larger value (e.g., 0.5) to allow more individuals to participate in local development and accelerate convergence to the optimal solution.

[0071] S1.3.4. Explore subgroups to implement random hiding strategies, simulate the behavior of rabbits randomly searching for hiding places in the wild, update positions by randomly perturbing the search space, and find new ranges of friction coefficient values.

[0072] S1.3.5 Develop a subgroup to execute a foraging strategy around food sources, simulating the behavior of rabbits foraging around food sources, and conduct a fine search around the current optimal friction coefficient;

[0073] S1.3.6, Apply probability p to the entire population. vcs Perform vertical cross search, randomly select two different rabbit individuals, and swap their values ​​in a certain dimension to achieve dimension-level cross operation, enhance population diversity, and help escape local optima;

[0074] S1.3.7, Update the global optimal solution k best Record the current minimum RMSE value;

[0075] S1.3.8, Check the termination condition: the maximum number of iterations T is reached. max Alternatively, if the RMSE value converges to a preset threshold, the optimal non-constant friction coefficient k is output. opt = k best Otherwise, return to S1.3.3.

[0076] Furthermore, the parameters for selecting the closing pattern of the downstream valve include: fast closing time. Quick-closing angle Slow closing time Slow closing angle The parameters of the unidirectional pressure regulating tower include: tower body diameter. Initial effective head Water supply pipe diameter .

[0077] S2. Construct a physical information neural network proxy model, embedding the hydraulic transient control equations as physical constraints into the network loss function, and using the k obtained in step S1. opt Substituting known constants into the calculation of non-constant friction terms in physical constraints reduces the dimensionality of network optimization and locks in the energy dissipation characteristics of transient flow, thereby enabling the prediction of water hammer characteristic quantities.

[0078] The specific methods for step S2 include:

[0079] S2.1 Randomly select 80% of the expected sample set as the training set and the remaining 20% ​​as the test set;

[0080] S2.2 Constructing a Physical Information Neural Network (PINN) surrogate model: The main body of the network adopts a fully connected feedforward neural network structure, which includes an input layer, multiple hidden layers and an output layer. The activation function adopts a local adaptive activation function to capture the high-frequency oscillation characteristics of water hammer pressure.

[0081] S2.3, The optimal non-constant friction coefficient k obtained by optimizing the ART-VCS algorithm in S13. opt Physical constraint layer embedded in the PINN model: The sum of squared residuals from the continuity and momentum equations describing transient pipe flow is introduced as a physical regularization term in the loss function, where the non-constant friction term uses a coefficient k. opt The loss function is in the form of:

[0082] (3)

[0083] In the formula,

[0084] (4)

[0085] (5)

[0086] in and The corresponding residuals for the continuity equation and the momentum equation are as follows:

[0087] (6)

[0088] (7)

[0089] The above partial derivatives are calculated using automatic differentiation (AD); For boundary condition loss, These are adaptive weighting coefficients;

[0090] S2.4. Train the PINN model using the training set and update the network parameters using the Adam optimizer;

[0091] S2.5. Test the trained prediction model using the test set; if the error is large (coefficient of determination) If so, adjust the network hyperparameters and retrain; if Then, this is used as the objective function of the multi-objective optimization model. The algorithm expression is as follows:

[0092] (8)

[0093] S3. Using the physical information neural network proxy model constructed in step S2 as the objective function, establish a multi-objective optimization model, and use an evolutionary algorithm with a directional sampling mechanism to optimize and solve the high-dimensional water hammer protection parameters to obtain the Pareto optimal solution set.

[0094] The specific methods for step S3 include:

[0095] The PINN surrogate model established by S2 is introduced as the objective function in the multi-objective optimization algorithm. minimize, maximize, Minimize the cost of water hammer protection as the optimization objective. Specifically, the cost of water hammer protection is determined by the volume of the unidirectional surge tank. Minimize the representation, and its calculation formula is:

[0096] .

[0097] The large-scale multi-objective evolutionary algorithm based on directional sampling (LMOEA-DS) is used as the main algorithm to establish a large-scale multi-objective optimization mathematical model for the optimization problem of pump shutdown water hammer protection measures, and iteratively search for the optimal solution set of pump shutdown water hammer protection design parameters.

[0098] Specifically, the steps to obtain the optimal set of design parameters for pump shutdown water hammer protection are as follows:

[0099] S3.1 Input initialization parameters to generate the initial population. Set the number of iterations to g=1. The initialization parameters for LMOEA-DS are population size N, maximum function evaluation count MaxFEs, and scaling factor. , representing the number of solutions Number of samples The dimension d of the optimization variables and the constraints on the input variables;

[0100] S3.2 Select high-quality representative solutions from the current population based on angle-penalized distance to form the representative solution set RS;

[0101] S3.3, Optimization Stage Judgment: Calculate the ratio of the current function evaluation count FEs to the maximum evaluation count MaxFEs. If... If the target sampling strategy is not executed (step S3.4), then the multi-source learning competitive group optimization strategy is executed directly (step S3.5).

[0102] S3.4, directional sampling strategies include:

[0103] Generate random numbers rand, if If convergence direction sampling is performed, otherwise diversity direction sampling is performed.

[0104] Convergence direction sampling: For the representative solution x, its approximate gradient is estimated using simultaneous perturbation stochastic approximation (SPSA).

[0105] (9)

[0106] Where Δ is the perturbation vector following a Bernoulli distribution, and c is the perturbation step size; combined with the boundary sampling direction, a set of convergent sampling directions is constructed, and sampling solutions are generated along the convergent direction;

[0107] Diversity-oriented sampling: The representative solution set RS is randomly divided into two equal-sized subsets. and For each pair of representative solutions Construct sampling direction And orthogonal directions, forming a diverse sampling direction set, and generating sampling solutions along the diverse directions;

[0108] Apply fuzzy variable operators to the generated sampled solutions to dynamically adjust the fuzziness of the solutions:

[0109] (10)

[0110] Where δ is the fuzzy coefficient, which decreases linearly with the number of iterations. Output the sampled population and jump to S3.6;

[0111] S3.5, Multi-source learning competitive group optimization strategies include:

[0112] The fitness value of each individual in the population is calculated based on the displacement density estimation.

[0113] Generate random numbers rand, if A global sorting mechanism is used for pairing (sorted in descending order of fitness, the top 50% are the winners set, and the bottom 50% are the losers set, then randomly paired); otherwise, a random pairing mechanism is used (randomly divided into several pairs, with the one with the higher fitness being the winner); one individual is randomly selected from each pairing group to form a learning set; for each loser and individual in the learning set, the position is updated using the competitive group optimizer formula; the updated population is merged with the current population, environmental selection is performed, and offspring populations are generated;

[0114] S3.6. Merge the updated population with the current population, perform environmental selection, and generate offspring population;

[0115] S3.7 Check if the termination condition is met: the maximum number of function evaluations MaxFEs has been reached. If it is met, proceed to the next step; otherwise, jump to S3.2 and increment the iteration count by 1.

[0116] S3.8 Output the optimal solution set of design parameters for pump shutdown water hammer protection.

[0117] Furthermore, under the premise of ensuring the safe operation of the water diversion project, the following constraints must be met based on the characteristics of each water hammer protection device:

[0118] (1) Valve closing constraint after pump:

[0119] (11)

[0120] In the formula, To expedite the closing time; For quick closing angle; This refers to the slow shutdown time; For slow closing angle;

[0121] (2) Extreme value constraints of the pipeline system:

[0122] (12)

[0123] (13)

[0124] (14)

[0125] In the formula, This represents the maximum water hammer pressure value. This represents the maximum steady-state pressure along the route. This is the minimum negative pressure value; This is the maximum reverse rotation speed of the water pump; This refers to the rated speed of the water pump;

[0126] (3) The upper and lower bounds of other input variables must be within the specified range of values, i.e., the constraint is:

[0127] (15)

[0128] In the formula, The diameter of the tower body; This is the initial effective head; This refers to the diameter of the water supply pipe.

[0129] S4. Perform multi-index decision analysis on the Pareto optimal solution set obtained in step S3 to determine the optimal combination of water hammer protection parameters.

[0130] The specific methods for step S4 include:

[0131] S4.1 Determine the decision criterion set {C1, C2, C3, C4} = {maximum water hammer pressure H} max Minimum negative pressure H min Maximum reverse rotation speed of water pump V max unidirectional pressure regulating tower volume V d Based on engineering requirements, the optimal criterion C is identified. B And worst-case criterion C W ;

[0132] S42. Construct the optimal criterion C using the 1-9 scaling method. B Compared to other criteria C j Preference comparison vector ,in Representing the optimal criterion C B Relative to criterion C j The degree of importance; and the relationship between all criteria and the worst criterion C. W Preference comparison vector ,in Representation C j Compared to the worst criterion C W The importance of.

[0133] S4.3 Solve the following minimax optimization problem to obtain the optimal criterion weight vector.

[0134] :

[0135]

[0136] (16)

[0137] In the formula, Optimal Criterion The corresponding weights worst-case criterion The corresponding weights are n, which is the total number of criteria (n=4 in this application), and ξ is the maximum absolute deviation. The optimal weight vector is obtained by minimizing ξ.

[0138] S4.4 Calculate the consistency ratio CR:

[0139] (17)

[0140] In the formula, This is the optimal target value for S43. The consistency index is used. If CR < 0.1, the consistency test is passed and the weight vector is accepted; otherwise, the preference comparison vector is adjusted and returned to S4.2 for recalculation.

[0141] S4.5 Construct a normalized decision matrix for the Pareto optimal solution set;

[0142] S4.6 Construct a weighted normalized decision matrix using the weight vectors obtained by the BWM method;

[0143] S4.7 Determine the positive ideal solution (Each indicator is a weighted optimal value) and negative ideal solution (The weighted worst value is used for each indicator).

[0144] S4.8 Calculate the Euclidean distance between each solution and the ideal solution. The distance between the sum and the negative ideal solution ;

[0145] S4.9 Calculate the relative closeness ,according to Descending sort scheme, select The largest possible solution is taken as the optimal compromise.

[0146] Furthermore, based on different biased optimization objectives, the optimal and worst-case criteria determined in S41 can be dynamically adjusted. The optimization objectives are divided into three types: safety and reliability priority (based on H...). max or H min (Optimal criterion), prioritizing cost-effectiveness (with V) d (For optimal criteria) and unbiased optimization.

[0147] S5. Substitute the optimal parameters obtained in step S4 into the hydraulic transient model in step S1 for verification. If the deviation between the verification result and the predicted value of the surrogate model exceeds the set threshold, then supplement the verification data into the surrogate model training set to trigger closed-loop iterative optimization.

[0148] Furthermore, in step S5, the theoretically optimal design parameters for pump shutdown water hammer protection are substituted into the high-precision water hammer numerical calculation model (embedded kopt characteristic line method model) constructed in S1 to determine the optimization effect; it is judged whether the optimization effect meets the expected effect, and secondary optimization can be performed if necessary.

[0149] Preferably, in step S1, the expression for the friction term in the non-constant friction model is: ,in For a constant coefficient of friction, The non-constant friction correction term is calculated using the following formula:

[0150] ;

[0151] In the formula, k is the Brunone friction coefficient (i.e., frictional resistance parameter), g is the acceleration due to gravity, D is the pipe diameter, a is the wave velocity, V is the flow velocity, t is time, and x is the spatial coordinate. The improved swarm intelligence optimization algorithm obtains the optimal value of k through search. opt .

[0152] Preferably, in step S2, the loss function L of the physical information neural network includes a data error term, a physical residual term, and a boundary condition loss term, expressed as: ,

[0153] In the formula,

[0154] , ,

[0155] in and The corresponding residuals for the continuity equation and the momentum equation are as follows:

[0156] ;

[0157] ;

[0158] Where L: total loss function of Physical Information Neural Network (PINN); h: pressure head; v: flow velocity; t: time variable; x: spatial coordinate variable along the pipe axis; a: water hammer wave velocity; g: gravitational acceleration; f: Darcy-Weisbach constant friction coefficient; D: pipe diameter; kopt: Brunone friction coefficient.

[0159] Preferably, in step S1, the improved swarm intelligence optimization algorithm introduces a cross-dimensional information interaction mechanism during the population update process to escape local optima in a highly nonlinear friction parameter space; the cross-dimensional information interaction mechanism includes vertical cross search, that is, randomly selecting two different individuals and exchanging their variable values ​​in a certain dimension.

[0160] Preferably, in step S3, the evolutionary algorithm with a directional sampling mechanism executes a directional sampling strategy in the early stages of iteration, including convergence direction sampling and diversity direction sampling; the convergence direction sampling uses simultaneous perturbation stochastic approximation to estimate the approximate gradient, and the calculation formula is:

[0161] ,

[0162] Where Δ is the perturbation vector following a Bernoulli distribution, and c is the perturbation step size; diversity direction sampling generates sampled solutions by constructing orthogonal directions representing the solutions, and dynamically adjusts the ambiguity of the generated sampled solutions using fuzzy variable operators, the calculation formula being:

[0163] Where δ is the fuzzy coefficient;

[0164] U j : The upper bound of the j-th optimization variable; L j : The lower bound of the j-th optimization variable;

[0165] X new,j : The value of the newly generated sampled solution in the j-th dimension; N(0,1): Standard normal distribution, used to generate random perturbation terms, so that the sampled solution has a certain degree of randomness and exploratory nature in the vicinity of the representative solution.

[0166] Preferably, in step S4, the multi-index decision analysis method uses a weight allocation method based on the comparison of best and worst criteria to determine the index weights, and obtains the optimal criterion weight vector by solving the following minimax optimization problem. ;

[0167]

[0168]

[0169] ;

[0170] In the formula, The weights of the optimal criterion Let ξ be the weight of the worst criterion, n be the total number of criteria, and ξ be the maximum absolute deviation. The optimal weight vector is obtained by minimizing ξ.

[0171] Calculate the consistency ratio (CR):

[0172] In the formula, To achieve the optimal target value, As a consistency index, if CR < 0.1

[0173] Then, the weight vector is accepted, and the optimal solution is selected by combining it with a sorting method based on the distance to the ideal solution.

[0174] Preferably, the optimal criterion C is dynamically adjusted based on engineering requirements. B And worst-case criterion C W This enables adaptive solution decision-making under scenarios prioritizing safety and reliability, cost-effectiveness, or unbiased optimization.

[0175] Preferably, the water hammer protection parameters include the operating parameters of the downstream valve and the structural parameters of the pressure regulating device; the operating parameters of the downstream valve include the fast closing time, fast closing angle, slow closing time, and slow closing angle; the structural parameters of the pressure regulating device include the tower diameter, the initial effective head, and the diameter of the water supply pipe.

[0176] Example:

[0177] A hydraulic transient simulation model consistent with the actual operation of water transfer projects is established, and an ART-VCS optimized non-constant friction model is introduced to construct the desired sample set:

[0178] Based on the topology diagram of a water diversion project, a hydraulic transient simulation model is established using the method of characteristics. The boundary conditions include: upstream reservoir boundary (constant water level), downstream pool boundary (constant water level), downstream valve boundary (two-stage closing law of fast and slow closing), air valve boundary (intake and exhaust characteristic curves), unidirectional pressure regulating tower boundary (continuity equation and water level-flow relationship), and pumping station boundary (full characteristic curves).

[0179] Collect and organize monitoring data on flow rate, pressure, water level, and the operating status of pumps, downstream valves, and unidirectional pressure regulating towers along the pipeline. Select two sets of monitoring data, each lasting 10 minutes under actual operating conditions with a step size of 0.01 seconds, for model calibration and verification. Control the verification error within 10% to obtain a hydraulic transient simulation model that conforms to the actual operation of the water diversion project.

[0180] Using Brunone's non-constant friction model, the expression for its friction term is:

[0181] (1)

[0182] In the formula, The constant friction coefficient is calculated using the Darcy-Weisbach formula; The correction term for non-constant friction is given by the following formula:

[0183] (2)

[0184] In the formula, Let be the instantaneous acceleration, 'a' be the wave velocity, and 'k' be the Brunone friction coefficient, which are key empirical parameters that need to be calibrated through an optimization algorithm. The ART-VCS algorithm is used to solve for k in [0,1], with the objective function being to minimize the root mean square error (RMSE) between the calculated pressure head h and the measured pressure head.

[0185] Parameter settings for the ART-VCS algorithm: rabbit population size N=30, maximum number of iterations T max =100, adaptive scaling factor α=0.4, vertical crossover probability p vcs =0.3. Following the iterative solution steps S131-S138, the optimal non-constant friction coefficient k is finally obtained. opt =0.047, corresponding to RMSE=0.185m, which is a 54.2% decrease in RMSE compared to the constant friction model.

[0186] Select the following parameters for the pump downstream valve: fast closing time T1∈[1,5]s, fast closing angle θ1∈[30,80]°, slow closing time T2∈[10,60]s, slow closing angle θ2∈[5,20]°, and the diameter D of the unidirectional pressure regulating tower. t ∈[2,8]m, initial effective head H t ∈[10,30]m, water supply pipe diameter D b ∈[0.2,0.8]m is used as the optimization variable. 1000 sets of input parameters are generated using the Latin hypercube sampling method, and then substituted into the embedding k. opt The water hammer calculation model outputs the corresponding H for each group. max H min and V max Construct the desired sample set.

[0187] Constructing a PINN surrogate model embedding ART-VCS to optimize non-constant friction parameters:

[0188] 800 sets were randomly selected as the training set, and 200 sets were selected as the test set. The PINN network was constructed: an input layer with 7 neurons (corresponding to 7 optimization variables), 8 hidden layers with 64 neurons per layer, using the Local Adaptive Activation Function (LAAF) as the activation function, and an output layer with 3 neurons (corresponding to H...). max H min and V max Physical constraint embedding: Physical residual term MSE in the loss function. f Let be the sum of squared residuals of the continuity equation and the momentum equation at collocation points within the computational domain, where:

[0189] (3)

[0190] (4)

[0191] In the formula, h is the pressure head (m), v is the flow velocity (m / s), a is the wave velocity (m / s), g is the gravitational acceleration (m / s²), f is the Darcy-Weisbach friction coefficient, D is the pipe diameter (m), and k opt =0.047 is the optimal non-constant friction coefficient calibrated by the ART-VCS algorithm.

[0192] The Adam optimizer was used for 1000 training epochs, with an early stopping mechanism to prevent overfitting. Test set results showed R²=0.991 and relative L² error ≤2.35e-2, meeting the accuracy requirements.

[0193] The LMOEA-DS algorithm is used to solve a multi-objective optimization problem and obtain the Pareto optimal solution set.

[0194] Set the LMOEA-DS parameters: population size N=100, MaxFEs=50000, representing the number of solutions n. r =20, scaling factor β=0.6, number of samples n s =50. Iterative search is performed according to S3.1-S3.8, with the same constraints. Convergence occurs in the 487th generation, yielding a Pareto front containing 100 non-dominated solutions. The HV index is 31580.4, and the SP index is 0.4836, indicating good convergence and uniform distribution of the solution set.

[0195] The BWM-TOPSIS method was used for optimal scheme selection.

[0196] For the unbiased optimization scenario, the optimal criterion is determined as follows: The worst-case criterion is Construct comparison vectors and Solve the following minimax optimization problem:

[0197]

[0198] (5)

[0199] The weight vector is obtained as w = (0.38, 0.32, 0.18, 0.12). The consistency ratio is then calculated.

[0200] (6)

[0201] The consistency check passed. The relative similarity was calculated. The optimal solution has the following parameters: T1 = 2.8s, θ1 = 62°, T2 = 35s, θ2 = 12°. =4.5m, =18m, =0.45m. Corresponding forecast indicator: =18.2m, =1.8m, =0.85 , =286.3m³. Compared with the previous optimization, the maximum water hammer pressure was reduced by 21.25%, the minimum negative pressure was improved by 47.72%, and the total valve closing time was shortened by 5.43%. The same logic applies to safety-first and cost-first solutions.

[0202] Verify the optimization effect:

[0203] Substitute the optimal solution parameters into the embedding The characteristic line method model was used for verification. The relative deviation between the simulation results and the PINN prediction values ​​was less than 2%. All indicators met the safety constraints. The maximum water hammer pressure was reduced by 21.3% and the volume of the pressure regulating tower was reduced by 38.7% compared with before optimization, achieving the expected optimization effect.

[0204] In summary, the present invention, compared with the prior art, has the following significant advantages:

[0205] (1) Improved efficiency and accuracy of non-constant friction parameter calibration: Adaptive Rabbit Optimization Algorithm combined with Vertical Cross-Search (ART-VCS) is used to replace the traditional trial-and-error method or ant colony algorithm for calibrating the empirical coefficients of non-constant friction. A classic water conveyance system example verifies that the ART-VCS calibration... The algorithm reduces the root mean square error (RMSE) of pressure head prediction by 54.2% compared to the constant friction model and improves it by approximately 12% compared to the original rabbit optimization algorithm, outperforming traditional trial-and-error methods and particle swarm optimization. Vertical cross-search effectively escapes local optima through dimensional cross-operations, and its adaptive population partitioning and dynamic equilibrium exploration and development demonstrate stronger global optimization capabilities and faster convergence speed in highly nonlinear friction calibration problems.

[0206] (2) Enhanced physical consistency and robustness of the prediction model: A PINN surrogate model with ART-VCS optimization of non-constant friction parameters was constructed, and the optimized and calibrated non-constant friction terms were embedded into the neural network loss function in the form of physical constraints. The Local Adaptive Activation Function (LAAF) was used to enhance the network's ability to capture high-frequency water hammer signals. The coefficient of determination R² on the test set reached 0.991, and the relative L² error was ≤2.35e-2. Even under sparse or noisy training data conditions, it can still maintain high-precision prediction, overcoming the defect of poor physical interpretability of pure data-driven models.

[0207] (3) Significantly improved high-dimensional search capability: The LMOEA-DS algorithm, designed specifically for large-scale parameter optimization, is adopted, which integrates a directional sampling strategy that combines convergence direction sampling and diversity direction sampling, and combines the local fine-grained search of the competitive group optimizer. In the tower valve joint control problem, the HV index reaches 31580.4, the SP index is 0.4836, the convergence speed is improved by about 40% compared with NSGA-II, the Pareto front span reaches 14.62 m (pressure), 4.36 m (negative pressure) and 240 s (time), and the solution set is evenly distributed, effectively solving the problems of premature convergence and diversity loss in high-dimensional evolutionary algorithms.

[0208] (4) Scientific and objective decision-making process: The BWM-TOPSIS method is used for scheme optimization. The BWM method obtains objective weights with good consistency by only requiring pairwise comparison of the optimal and worst criteria, with a consistency ratio CR < 0.1, which greatly reduces the workload of subjective judgment. The relative closeness is calculated in combination with TOPSIS. According to engineering examples, the optimized scheme reduces the maximum water hammer pressure by 21.25%, improves the minimum negative pressure by 47.72%, and shortens the total valve closing time by 5.43%, which is significantly better than single valve control or experience-based design.

[0209] (5) High system integration: A closed-loop optimization system is constructed, consisting of “ART-VCS optimization of non-constant friction → PINN proxy modeling → LMOEA-DS evolutionary optimization → BWM-TOPSIS decision → verification”. The data and parameters of each module are deeply coupled, which improves the practicality of engineering and the robustness of the system.

[0210] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A PINN-LMOEA-DS intelligent optimization method for water hammer protection based on adaptive rabbit optimization of non-constant friction, characterized in that, Includes the following steps: S1. Establish a hydraulic transient model for the water diversion project, introduce a non-constant friction model to describe the transient flow process, and use the friction parameter as the optimization variable. An improved swarm intelligence optimization algorithm with a cross-dimensional information interaction mechanism is used to globally calibrate the friction parameter, obtaining the optimal non-constant friction parameter k. opt ; S2. Construct a physical information neural network proxy model, embedding the hydraulic transient control equations as physical constraints into the network loss function, and using the k obtained in step S1. opt Substituting known constants into the calculation of the non-constant friction term of the physical constraints reduces the dimensionality of network optimization and locks in the transient flow energy dissipation characteristics, thereby enabling the prediction of water hammer characteristic quantities. S3. Using the physical information neural network proxy model constructed in step S2 as the objective function, establish a multi-objective optimization model, and use an evolutionary algorithm with a directional sampling mechanism to optimize and solve the high-dimensional water hammer protection parameters to obtain the Pareto optimal solution set. S4. Perform multi-index decision analysis on the Pareto optimal solution set obtained in step S3 to determine the optimal combination of water hammer protection parameters. S5. Substitute the optimal parameters obtained in step S4 into the hydraulic transient model described in step S1 for verification. If the deviation between the verification result and the predicted value of the surrogate model exceeds the set threshold, then supplement the verification data into the surrogate model training set to trigger closed-loop iterative optimization.

2. The method according to claim 1, characterized in that: In step S1, the expression for the friction term in the non-constant friction model is: ,in For a constant coefficient of friction, The non-constant friction correction term is calculated using the following formula: ； In the formula, k is the Brunone friction coefficient, i.e., the frictional resistance parameter; g is the acceleration due to gravity; D is the pipe diameter; a is the wave velocity; V is the flow velocity; t is time; and x is the spatial coordinate. The improved swarm intelligence optimization algorithm obtains the optimal value of k through search and solution. .

3. The method according to claim 2, characterized in that: In step S2, the loss function L of the physical information neural network includes a data error term, a physical residual term, and a boundary condition loss term, expressed as: ， In the formula, ， ， in and The corresponding residuals for the continuity equation and the momentum equation are as follows: ； ； Where L: total loss function of Physical Information Neural Network (PINN); h: pressure head; v: flow velocity; t: time variable; x: spatial coordinate variable along the pipe axis; a: water hammer wave velocity; g: gravitational acceleration; f: Darcy-Weisbach constant friction coefficient; D: pipe diameter; kopt: Brunone friction coefficient.

4. The method according to claim 1, characterized in that: In step S1, the improved swarm intelligence optimization algorithm introduces a cross-dimensional information interaction mechanism during the population update process to escape local optima in a highly nonlinear friction parameter space. The cross-dimensional information interaction mechanism includes vertical cross search, which involves randomly selecting two different individuals and exchanging their variable values ​​in a certain dimension.

5. The method according to claim 1, characterized in that: In step S3, the evolutionary algorithm with a directional sampling mechanism executes a directional sampling strategy in the early stages of iteration, including convergence direction sampling and diversity direction sampling; the convergence direction sampling uses simultaneous perturbation stochastic approximation to estimate the approximate gradient, and the calculation formula is: ， Where Δ is the perturbation vector following a Bernoulli distribution, and c is the perturbation step size; the diversity direction sampling generates sampled solutions by constructing orthogonal directions representing the solutions, and dynamically adjusts the ambiguity of the generated sampled solutions using fuzzy variable operators, the calculation formula being: Where δ is the fuzzy coefficient; U j : The upper bound of the j-th optimization variable; : The lower bound of the j-th optimization variable; X new,j : The value of the newly generated sampled solution in the j-th dimension; N(0,1): Standard normal distribution, used to generate random perturbation terms, so that the sampled solution has a certain degree of randomness and exploratory nature in the vicinity of the representative solution.

6. The method according to claim 1, characterized in that: In step S4, the multi-index decision analysis method uses a weight allocation method based on the comparison of best and worst criteria to determine the index weights, and obtains the optimal criterion weight vector by solving the following minimax optimization problem. ; ， ， ； In the formula, The weights of the optimal criterion Let ξ be the weight of the worst criterion, n be the total number of criteria, and ξ be the maximum absolute deviation. The optimal weight vector is obtained by minimizing ξ. Calculate the consistency ratio (CR): In the formula, To achieve the optimal target value, As a consistency index, if CR < 0.1 Then, the weight vector is accepted, and the optimal solution is selected by combining it with a sorting method based on the distance to the ideal solution.

7. The method according to claim 6, characterized in that: The optimal criterion C is dynamically adjusted based on engineering requirements. B And worst-case criterion C W This enables adaptive solution decision-making under scenarios prioritizing safety and reliability, cost-effectiveness, or unbiased optimization.

8. The method according to claim 1, characterized in that: The water hammer protection parameters include the operating parameters of the downstream valve and the structural parameters of the pressure regulating device; the operating parameters of the downstream valve include the fast closing time, fast closing angle, slow closing time, and slow closing angle; the structural parameters of the pressure regulating device include the tower diameter, the initial effective head, and the diameter of the water supply pipe.

Images