Water hammer protection method based on multi-preference constraint multi-objective optimization algorithm

Abstract

The application discloses a water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm, and comprises the following steps: constructing a pipeline water hammer pressure calculation model based on fluid mechanics equations, solving by using a characteristic line method, and obtaining two pairs of characteristic line compatibility equations; selecting decision variables, fitting four target functions, and connecting the four target functions into an MSCMOEA algorithm; judging a calculation stage by a relationship between a current iteration number and a maximum iteration number, selecting a model meeting a triggering condition in models I, II and III according to a current feasible solution ratio, optimizing, iterating until a termination condition is met, and outputting a pareto optimal solution set. The application aims to realize multi-objective optimization of a water hammer protection scheme on the basis of ensuring stable operation of an engineering, further ensure engineering safety, and effectively reduce engineering cost.

Description

Technical Field

[0001] This invention relates to the field of water hammer protection technology, and in particular to a water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm. Background Technology

[0002] Water hammer is a major safety hazard in water conveyance projects. When a pumping station unexpectedly stops or valves are opened or closed rapidly, sudden changes in fluid momentum within the pipeline can cause severe pressure fluctuations, potentially leading to serious accidents such as pipeline rupture and equipment damage. Therefore, optimizing water hammer protection measures and balancing safety and economy are core research topics in the field of water conveyance engineering. All types of water conveyance systems need to reduce water hammer risks through scientific protection and optimization methods to ensure safe system operation.

[0003] Currently, water hammer protection mainly relies on parameter optimization of equipment such as hydraulic butterfly valves, air valves, and unidirectional pressure regulating towers. Existing technologies revolve around adjusting the parameters and controlling the operating conditions of these devices, with the core approach being "hydraulic calculation modeling + optimization algorithm solution" to optimize the protection scheme. In hydraulic calculation, the method of characteristics (MOC) is widely used, transforming partial differential equations into total differential equations to simulate the dynamic changes in pipeline pressure and flow velocity under different operating conditions, providing data support for optimization design. Regarding optimization methods, traditional methods mainly rely on empirical judgment, gradient methods, and numerical simulation trials, suppressing pressure fluctuations by adjusting the parameters of protective equipment or optimizing valve closing patterns and other operating parameters. Intelligent algorithms (genetic algorithms, sparrow search algorithms, NSGA-II, etc.) are gradually being applied, searching for the optimal set of protection parameters by constructing multi-objective optimization models. Some schemes combine machine learning models to improve the mapping accuracy between variables and objectives, shortening the optimization cycle.

[0004] Existing technologies still have significant shortcomings and are difficult to meet the precise optimization needs of complex water conveyance systems: First, the optimization strategies lack dynamic adaptability, with fixed optimization logic running through the entire iteration process, failing to adjust the focus of exploration and convergence according to the calculation stage, resulting in strong blind optimization; Second, the optimization model has a low degree of matching with the actual scenario, failing to dynamically adjust the solution direction based on the proportion of feasible solutions, resulting in many invalid calculations and low optimization efficiency; Third, the optimization results have poor practicality, with redundant solution sets, few optional solutions, and excessive reliance on manually setting target weights, resulting in poor constraint satisfaction and affecting the implementation of the solution. Summary of the Invention

[0005] To address the aforementioned problems, this invention provides a water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm. By combining the MSCMOEA optimization algorithm, it solves the problems of strong blindness, redundant solution sets, over-reliance on target weights, and low efficiency in existing water hammer protection optimization methods while coordinating system safety, economy, and operational stability.

[0006] A water hammer protection method based on a multi-preference constrained multi-objective optimization algorithm includes the following steps:

[0007] Step S1, construct the pipeline water hammer pressure calculation model:

[0008] A pipeline water hammer pressure calculation model was constructed based on fluid mechanics equations. The method of characteristics was used to solve the model, resulting in two pairs of compatible equations of characteristic lines. The initial model establishment was completed.

[0009] Step S2: Define the decision variables to be optimized and determine the constraints and objective function.

[0010] Select the surge tank diameter Ds and the initial water depth Water supply pipe diameter D, pressure regulating tower installation location Quick-closing angle Quick closing time Slow closing angle Slow closing time Air valve position , As a decision variable;

[0011] The constraints include: pump speed constraints, maximum water hammer, minimum water hammer pressure constraints, valve closure constraints after pump, input variable up and down constraints, and unidirectional pressure regulating tower safety factor constraints.

[0012] The objective function is obtained by taking the dimensionless maximum water hammer pressure, the dimensionless minimum water hammer pressure, the dimensionless minimum speed of the pump, and the water hammer protection cost as optimization objectives.

[0013] Step S3: Input the fitted optimized objective function into the MSCMOEA algorithm and set the initial parameters:

[0014] Step S4, First Stage:

[0015] Compare the ratio of the previous iteration count g to the preset total iteration count G. When g < 0.25G, proceed to the first stage.

[0016] Obtain the proportion of feasible solutions for the main population P(N). The feasible solution is proportional Compared with the preset benchmark Perform size comparison; determine the target computational model based on the comparison results; if Then, Model I is selected as the target computational model, and the early exploration of the optimal solution is completed around the target preference; if Then, Model III is selected as the target calculation model; the calculation operation is performed to complete the fine optimization of the optimal solution, with constraint preference as the main factor and target as the auxiliary factor, and the calculated P(N) is output.

[0017] Step S5, Second Stage:

[0018] When 0.25G ≤ g < 0.85G, the second stage begins;

[0019] Obtain the feasible solution ratio k2 of the auxiliary population A(n), and compare the feasible solution ratio k2 with the preset benchmark kb; determine the target computation model based on the comparison result, if... If model II is selected as the target calculation model, the optimal solution is selected through a division of labor and cooperation method, i.e., dual preference of objective and constraint. If k2 > kb, then model III is selected as the target calculation model. The calculation operation is performed to complete the fine optimization of the optimal solution and output the calculated P(N).

[0020] Step S6, the third stage:

[0021] When g≥0.85G, enter the third stage, directly execute Model Ⅲ to calculate and output P(N), and complete the fine optimization of the optimal solution;

[0022] Step S7, obtain the optimal solution set:

[0023] Determine whether the iteration termination condition is met. The termination condition includes two items: first, the number of iterations reaches a preset total number; second, within a preset number of generations, the dimensionless change of the optimal solution of the objective function in the population for two consecutive generations is less than a preset threshold. When either termination condition is met, output the Pareto optimal solution set that satisfies the constraints.

[0024] Furthermore, in step S1, the fluid dynamics equations, including the continuity equation and the momentum equation, are calculated using the following formulas:

[0025]

[0026]

[0027] In the formula, h is the piezometric head; g is the acceleration due to gravity, 9.81 m / s². 2 x represents the node coordinates; The water hammer wave velocity; t is the pipe inclination angle; v is the flow velocity; t is time; f is the friction factor; D is the pipe diameter in meters.

[0028] Furthermore, in step S1, the characteristic line compatibility equation includes a forward characteristic line compatibility equation and a reverse characteristic line compatibility equation.

[0029] The formula for calculating the compatibility equation of the positive characteristic lines is as follows:

[0030]

[0031] The formula for calculating the compatibility equation of reverse characteristic lines is as follows:

[0032]

[0033] in, and These are the coefficients for the corresponding positive and negative feature lines, respectively, and are calculated using the following formulas:

[0034]

[0035]

[0036] In the formula: , , , B and R are all calculation constants. , Let be the water heads at sections A and B at time t, respectively; let q be the flow rate at the calculated section at time t, satisfying q = vA; and let A be the cross-sectional area. Indicates the length of the feature line grid pipeline.

[0037] Furthermore, in step S2, the objective function is calculated using the following formula:

[0038]

[0039]

[0040]

[0041]

[0042] In the formula, This represents the maximum water hammer pressure value. This is the minimum water hammer pressure value; This is the dimensionless minimum speed of the pump; The volume of the pressure regulating tower is used as a cost indicator.

[0043] Furthermore, in step S2, the pump speed constraint is calculated using the following formula:

[0044]

[0045] in, The minimum dimensionless rotational speed; This refers to the time during which the reverse rotation speed is greater than the rated speed.

[0046] The maximum and minimum water hammer pressure constraints are calculated using the following formulas:

[0047]

[0048]

[0049] in: This represents the maximum water hammer pressure value. This refers to the rated pressure at the pump outlet. The maximum water hammer pressure value occurring at each node; The maximum allowable pressure for the pipeline at each node; This is the minimum water hammer pressure value; This represents the maximum allowable vacuum level for the pipeline. This represents the minimum water hammer pressure value at each node; The maximum allowable vacuum level for each node pipeline;

[0050] The calculation formula for the valve closing constraint after the pump is as follows:

[0051]

[0052] in, This refers to the valve closing angle during the rapid closing phase. This refers to the valve closing angle during the slow closing phase. To achieve the fast closing angle Time required; To achieve the slow closing angle Time required;

[0053] The upper and lower constraints of the input variables are calculated using the following formula:

[0054]

[0055] in, The diameter of the voltage regulating tower, The initial effective head is given by D, which is the diameter of the make-up water pipe. This means that the diameter of the make-up water pipe must be greater than the minimum allowable value and less than the diameter of the surge tank body. ;

[0056] The safety factor constraint for the unidirectional voltage regulating tower is calculated using the following formula:

[0057]

[0058] The final water depth after the unidirectional surge tank has finished replenishing water. This refers to the minimum water depth required for actual engineering projects.

[0059] Furthermore, in step S3, the initial parameters include: the maximum number of iterations G, and the feasible solution ratio judgment threshold. Initialize the primary population P(N), the auxiliary population A(n), and the historical solution set S.

[0060] Furthermore, in step S4, model I includes the following steps:

[0061] S411, Population Evolution:

[0062] Given a population P(N), for Parents are randomly generated from the population P(N), and offspring are generated using the GAhalf operator and stored in set Q. The original population P(N) and the offspring set Q are then merged into a mixed population H. ;

[0063] S412, Objective function correction:

[0064] Perform objective function offset calculation, assuming individuals in the mixed population H The original objective function is , The objective function after offset is calculated as follows:

[0065]

[0066] in, For reference point The j-th target value; The offset coefficient is dynamically adjusted based on the proportion of feasible solutions in the mixed population H, and is calculated using the following formula:

[0067]

[0068] in, To adjust the parameters, g is the current iteration number, and G is the maximum iteration number;

[0069] S413, Crowding Filter:

[0070] Calculate individuals within the same non-dominated layer in a mixed population H. Total congestion distance The formula and parameter definitions are as follows:

[0071]

[0072]

[0073] in, Individuals in the same non-dominated layer The crowding distance of the j-th target. For adjacent individuals, For individuals The function value of the j-th objective function; , These are the maximum and minimum values ​​of the target.

[0074] Perform non-dominated ordination on the mixed population H, calculate the crowding distance of individuals within the same non-dominated layer, sort them from high to low non-dominated level and from large to small crowding distance within the same layer, and screen to obtain a new generation population P(N) of size N.

[0075] Furthermore, in step S4, model III includes the following steps:

[0076] S421, Initialization and Input Parameters:

[0077] Input the current population P(N), auxiliary population A(n), historical solution set S, number of objective functions m, and preset reference point set. ;

[0078] The auxiliary population A(n) has an initial size of 0 and a maximum size of 0.3N;

[0079] The initial size of the historical solution set S is 0, and the size of the historical solution set S is... ;

[0080] The formula for calculating the reference point is as follows:

[0081]

[0082] Where s=2, Total number of reference points ;

[0083] S422, Parent Selection:

[0084] If the auxiliary population size ; Select non-dominant individuals from the auxiliary population A(n), sort them in descending order by crowding distance D(x), and select the top... A group of non-dominant individuals forms an elite parent-child set. ;

[0085] If the auxiliary population size Then, select the top 30% of individuals with the largest crowding distance from P(N) to form a temporary A(n) for the above operation, and simultaneously perform the above operation on the individuals. closest to the reference point Individuals, forming a reference point parent-child set The size is K;

[0086] Generate parent child set: Deduplication retains individuals with larger crowding distances, at a scale of ;

[0087] S423, Descendant Generation:

[0088] Based on parent-child set The population offspring set Q is generated using the Gaussian mutation crossover operator, including:

[0089] From father to son Two different parent individuals are randomly selected. and Generate offspring individuals The calculation formula is as follows:

[0090]

[0091]

[0092]

[0093] in, , The parent individuals are randomly selected, and n is the decision variable. For cross weights, The standard deviation of Gaussian noise. With a mean of 0 and a variance of Let g be a Gaussian random variable, and G be the current iteration number and the maximum iteration number, respectively.

[0094] Each generated descendant individual Add to the descendant set Q;

[0095] S424, Objective function correction:

[0096] A target priority coefficient is introduced to adjust the target value. The calculation formula is as follows:

[0097]

[0098]

[0099]

[0100] in, Let x be the average distance of individual x to the j-th target; Let x be the output value of the j-th original objective function; Let x be the corrected target value for the j-th objective. For the weighted target value, This represents the target priority coefficient.

[0101] S425, Layered Filtering:

[0102] The original population P(N) and the offspring set Q are merged into a mixed population H. ;

[0103] Pareto hierarchies are defined using non-dominated sorting, prioritizing the retention of individuals at higher levels. The specific formula for hierarchical division is as follows:

[0104]

[0105] in, express Dominate The smaller the rank number, the higher the individual's priority; the optimal individual has a rank of 1.

[0106] Individual and optimal correlation reference point The correlation degree is calculated using the following formula:

[0107]

[0108] The overall score is calculated using the following formula:

[0109]

[0110] in, To correct for crowding;

[0111] Sort individuals by rank in ascending order, prioritize retaining individuals with higher rank in the mixed population H, sort individuals in the same rank by comprehensive score in descending order, and select the top N individuals to form a new main population P(N).

[0112] S426, Updated historical solution set:

[0113] From the new dominant population P(N), for each reference point Select the closest individual, and select a total of t unique individuals, t=N / 5, and add them to the historical solution set S;

[0114] If |S|>2N, then remove the redundant solution with the smallest crowding distance, retain at most 2N solutions, and output P(N) and S.

[0115] Furthermore, in step S5, model II includes the following steps:

[0116] S511, Initialize input parameters:

[0117] Input the current population P(N), auxiliary population A(n), population size N, and number of objective functions m;

[0118] S512, Parent Set and Descendant Generation:

[0119] Individual crowding distance is calculated using formulas (19) and (20). Sort individuals from largest to smallest by crowding distance, select the top 30% of individuals in population P(N), and randomly select 20% of individuals from auxiliary population A(n) to form the parent set. Use the "adaptive crossover mutation" operator to generate the offspring set Q.

[0120] S513 Objective Function Modification and Hierarchical Filtering:

[0121] A "diversity reward factor" is introduced into the objective function value of offspring individuals to modify the objective function. The calculation formula is as follows:

[0122]

[0123] in, Let be the average distance between individual x and other individuals in the population on the j-th target dimension, and β be the reward coefficient with a value of 0.2.

[0124] A "stratified screening" mechanism was adopted for the merged mixed population. Screening is performed; stratification is performed according to formula (28), and individuals in the same stratum are sorted from largest to smallest according to "corrected crowding distance", with the same formula as (19) and (20). The N individuals with larger distances are retained to form a new main population P(N);

[0125] S514, Auxiliary Population Update and Result Output:

[0126] If A(n) is an empty set, then sort the individuals in descending order of their crowding distance, and select the top 0.3N feasible solutions from the new principal population P(N) and place them in A(n), with a size of n=0.3N;

[0127] If A(n) is not an empty set, then examine each individual in the new population P(N). If a new individual dominates an existing individual in A(n), or if the crowding distance is greater in the non-dominated layer, then replace the individual with the smallest crowding distance in A(n) in equal numbers to maintain the size of A(n) constant, i.e., n=0.3N.

[0128] Output the new primary population P(N) and secondary population A(n).

[0129] Furthermore, in step S7, the dimensionless change in the optimal solution of the objective function in the population for two consecutive generations is less than a preset threshold, calculated as follows:

[0130]

[0131]

[0132] in, The normalized target value; This represents the original output value of individual X on the k-th objective function; , This represents the historical maximum value of the objective function. This represents the historical minimum value of the objective function. Let be the optimal objective value of the k-th objective function in the t-th generation.

[0133] The beneficial effects of the present invention are: The water hammer protection method based on multi-preference constraint multi-objective optimization algorithm provided by the present invention determines the calculation stage by the number of iterations, dynamically adjusts the optimization strategy, achieves the effect of early exploration and late convergence, avoids blind optimization, and improves the targeting of water hammer protection optimization;

[0134] By selecting the appropriate model based on the proportion of feasible solutions, the current optimization requirements can be accurately matched. Fine-tuning can be performed when there are highly feasible solutions, and exploration can be performed when there are low feasible solutions, thereby reducing invalid calculations and improving the optimization efficiency of water hammer protection.

[0135] By storing the previous optimal solutions in the historical solution set S, we can avoid losing high-quality solutions during the iteration process, ensure that the Pareto optimal solution covers different protection requirements, enrich the optional solutions, and provide the best selection basis for engineering decisions.

[0136] The present invention will be explained in detail below with reference to the accompanying drawings and specific embodiments. Attached Figure Description

[0137] Figure 1 This is the overall flowchart of this method;

[0138] Figure 2 The flowchart for Model I in the MSCMOEA algorithm is shown below.

[0139] Figure 3 The flowchart for Model II in the MSCMOEA algorithm;

[0140] Figure 4 This is a flowchart of Model III in the MSCMOEA algorithm;

[0141] Figure 5 This is a simplified longitudinal section diagram of a water conveyance pipeline provided by the present invention. Detailed Implementation

[0142] Example 1:

[0143] This embodiment presents a water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm, such as... Figure 1 As shown, it includes the following steps:

[0144] Step S1, construct the pipeline water hammer pressure calculation model:

[0145] A pipeline water hammer pressure calculation model was constructed based on fluid mechanics equations, and the method of characteristics was used to solve it, resulting in two pairs of compatible equations of characteristic lines.

[0146] In hydraulic transients, the characteristic line is the propagation trajectory of pressure waves within a pipeline, representing the physical path of transient disturbance energy transfer. Hydraulic transients are caused by disturbances in the fluid within the pipeline due to valve opening and closing, pump start-up and shutdown, etc., manifesting as pressure waves and transient flow rates. These disturbances exhibit randomness; for example, valve closing times may fluctuate within a statistical interval, and pump start-up and shutdown rates may dynamically change due to operating conditions. Consequently, the resulting pressure wave propagation speed and the amplitude of transient flow rates also exhibit nondeterministic characteristics. The mathematical description of transient processes relies on the Saint-Venant equations. These partial differential equations cannot be directly solved analytically using universal solutions. The method of characteristics solves this problem through dimensionality reduction transformation. First, the partial differential equations are transformed into ordinary differential equations of the characteristic lines, and then further discretized into algebraic equations. The method of characteristics can transform the randomness of disturbances (such as the probability distribution of valve closing times and dynamic parameters of rate fluctuations) into quantifiable boundary conditions or initial parameters, which can be substituted into the discretized algebraic equations to capture and calculate the random characteristics of pressure and flow transients within the pipeline.

[0147] S11, The basic mathematical equations for establishing the mathematical model of a pressurized water pipeline are as follows:

[0148]

[0149]

[0150] In the formula: h is the piezometric head; g is the gravitational acceleration; x is the nodal coordinate; The water hammer wave velocity; denoted as pipe inclination angle; v as flow velocity; t as time; f as friction factor; and D as pipe diameter.

[0151] S12, using the method of characteristics to solve equations (1) and (2), performing dimensionality reduction and discretization, transforming the nonlinear hyperbolic partial differential equations into ordinary differential equations along the characteristic lines, resulting in two pairs of compatible equations along the characteristic lines:

[0152] 1. Positive characteristic lines ( (Water hammer wave propagation downstream), the formula for calculating the compatibility equation is as follows:

[0153]

[0154] 2. Reverse feature line ( (Water hammer wave counter-current propagation), the compatibility equation calculation formula is as follows:

[0155]

[0156] and These are the coefficients of the corresponding positive and negative feature lines, respectively:

[0157]

[0158]

[0159] In the formula: , , , B and R are all calculation constants. , Let be the water heads at sections A and B at time t, respectively; let q be the flow rate at the calculated section at time t, satisfying q = vA; and let A be the cross-sectional area. Indicates the length of the feature line grid pipeline.

[0160] Step S2, determine the optimization objective and constraints:

[0161] Define the decision variables to be optimized, and determine the constraints and objective function.

[0162] S21: Select the surge tank diameter Ds and the initial water depth. Water supply pipe diameter D, pressure regulating tower installation location Quick-closing angle Quick closing time Slow closing angle Slow closing time Air valve position , As a decision variable;

[0163] S22: Using the dimensionless maximum water hammer pressure, the dimensionless minimum water hammer pressure, the dimensionless minimum pump speed, and the water hammer protection cost (pressure regulating tower volume) as optimization objectives, four objective functions are obtained:

[0164]

[0165]

[0166]

[0167]

[0168] In the formula: This represents the maximum water hammer pressure value. This is the minimum water hammer pressure value; This is the dimensionless minimum speed of the pump; The volume of the pressure regulating tower is used as a cost indicator.

[0169] S23: Constraints include:

[0170] S231: Pump speed constraint:

[0171]

[0172] in: The minimum dimensionless rotational speed; This refers to the time during which the reverse rotation speed exceeds the rated speed.

[0173] S232: Maximum and minimum water hammer pressure constraints:

[0174]

[0175]

[0176] in: This represents the maximum water hammer pressure value. This refers to the rated pressure at the pump outlet. The maximum water hammer pressure value occurring at each node; The maximum allowable pressure for the pipeline at each node; This is the minimum water hammer pressure value; This represents the maximum allowable vacuum level for the pipeline. This represents the minimum water hammer pressure value at each node; This represents the maximum allowable vacuum level for each node's pipeline.

[0177] S233: Pump downstream valve closing constraint:

[0178]

[0179] in, This refers to the valve closing angle during the rapid closing phase. This refers to the valve closing angle during the slow closing phase. To achieve the fast closing angle Time required; To achieve the slow closing angle The time required.

[0180] S234: Input variable upper and lower constraints:

[0181]

[0182] in, The diameter of the voltage regulating tower, The initial effective head is given by D, which is the diameter of the make-up water pipe. This means that the diameter of the make-up water pipe must be greater than the minimum allowable value and less than the diameter of the surge tank body. .

[0183] S235: Safety factor constraint for unidirectional voltage regulating tower:

[0184]

[0185] The final water depth after the unidirectional surge tank has finished replenishing water. This refers to the minimum water depth required for actual engineering projects.

[0186] Step S3: Input the fitted optimization objective function into the MSCMOEA algorithm and set the initial parameters: maximum number of iterations G, feasible solution ratio judgment threshold. And the algebraic intervals of each stage, initialize the main population P(N), auxiliary population A(n), and historical solution set S (empty set, with a maximum capacity of 2N).

[0187] MSCMOEA (Multi-Subgroup Cooperative Multi-Objective Evolutionary Algorithm) is an intelligent algorithm for solving multi-objective optimization problems. It improves solution efficiency and solution quality through division of labor and cooperation.

[0188] The algorithm first divides the entire population into multiple subgroups, assigning different tasks to different subgroups: some focus on quickly approximating the optimal solution set (improving convergence), some are responsible for maintaining a uniform distribution of solutions (ensuring diversity), and some explore unknown regions (avoiding getting trapped in local optima).

[0189] When subgroups evolve independently, they adopt strategies that adapt to their own division of labor. Subsequently, subgroups periodically share superior individuals and evolutionary information, allowing subgroups with faster convergence to guide the direction and subgroups with stronger exploratory capabilities to avoid premature convergence, thereby achieving complementary advantages.

[0190] Finally, individuals from all subgroups are merged, and the next generation of the population is selected, until the algorithm converges. This approach balances convergence speed and solution set diversity when dealing with complex multi-objective problems.

[0191] Step S4, First Stage:

[0192] Compare the ratio of the previous iteration count g to the preset total iteration count G. That is, when the iteration count g < 0.25G, proceed to the first stage and compare the proportion of feasible solutions in P(N). and The size of the population is determined by performing calculations using either Model I or Model III, and the population P(N) is output. The workflow for Model I and Model III is as follows: Figure 2 , Figure 4 As shown;

[0193] S41 Early Exploration: If the proportion of feasible solutions in population P(N) Then, substitute the values ​​into Model I and perform the operation to complete the early exploration of the optimal solution around the target preference;

[0194] S411 Population Evolution: Input population P(N), and perform population evolution on... Parents are randomly generated from the population P(N), and offspring are generated using the GAhalf operator and stored in set Q. The original population P(N) and the offspring set Q are then merged into a mixed population H(N). );

[0195] The GAhalf operator is a current genetic algorithm, and its specific parameters and execution settings are as follows:

[0196] Parent selection strategy: Use tournament selection, tournament size k=3 (each time, 3 individuals are randomly selected from the population without replacement, and the individual with the best fitness is selected as the parent).

[0197] Crossover operation: Simulated binary crossover (SBX) is used, with crossover probability... Cross-distribution index ;

[0198] Mutation operation: Polynomial mutation is used, with mutation probability... (where n is the dimension of the decision variable), variation distribution index ;

[0199] Offspring generation: Offspring are generated through the above selection, crossover, and mutation steps and stored in set Q.

[0200] S412 Objective Function Correction: Perform objective function shift calculation, assuming the individuals in the mixed population H... The original objective function is ( The formula for calculating the offset objective function is as follows:

[0201]

[0202] in, For reference point The j-th target value;

[0203] Offset coefficient The proportion of feasible solutions in the mixed population H is dynamically adjusted, and the calculation formula is as follows:

[0204]

[0205] in, To adjust the parameters, g is the current iteration number, and G is the maximum iteration number.

[0206] S413 Crowding Degree Screening: Method for calculating total crowding distance of individuals, using formulas (19) and (20) to calculate the total crowding distance of individuals within the same non-dominated layer in the mixed population H. Total congestion distance The formula and parameter definitions are as follows:

[0207]

[0208]

[0209] in, Let be the crowding distance of individual xi in the same non-dominated layer for the j-th objective. For adjacent individuals, For individuals The function value of the j-th objective function; , These are the maximum and minimum values ​​of the target.

[0210] Perform non-dominated ordination on the mixed population H, calculate the crowding distance of individuals within the same non-dominated layer, sort them from high to low non-dominated level and from large to small crowding distance within the same layer, and screen to obtain a new generation population P(N) of size N.

[0211] S42 Fine-tuning: If Then, substitute the values ​​into Model III and perform the operation, using constraint preferences as the main factor and the objective as the auxiliary factor to complete the fine optimization of the optimal solution;

[0212] S421 Initialization and Input Parameters: Input the current population P(N), auxiliary population A (initial size 0, maximum size 0.3N), and historical solution set S (initial size 0, ...). ), Number of objective functions m, Preset reference point set .

[0213]

[0214] Where s=2, Total number .

[0215] S422 Parent Selection: If Select non-dominated individuals from the auxiliary population A(n), sort them in descending order by crowding distance D(x), and select the top... A group of non-dominant individuals forms an elite parent-child set. ;

[0216] like Then, by calculating using formula (19), the top 30% of individuals with the largest total crowding distance from P(N) are selected to form a temporary A(n) for the above operation, while simultaneously adjusting the individual... closest to the reference point Individuals, forming a reference point parent-child set The size is K;

[0217] Generate parent-child set The size after deduplication is .

[0218] S423 Descendant Generation: Based on Parent Generation The population's offspring set Q is generated using the Gaussian mutation crossover operator, as detailed below:

[0219] From father to son Two different parent individuals are randomly selected. and Generate offspring individuals The calculation formula is as follows:

[0220]

[0221]

[0222]

[0223] in, , The parent individuals are randomly selected, and n is the decision variable. For cross weights, The standard deviation of Gaussian noise. With a mean of 0 and a variance of Let g be a Gaussian random variable, and G be the current iteration number and the maximum iteration number, respectively.

[0224] Each generated descendant individual Add it to the descendant set Q.

[0225] S424 Objective Function Correction: Introducing an objective priority coefficient, the objective value is corrected using formula (25):

[0226]

[0227]

[0228]

[0229] in, Let x be the average distance of individual x to the j-th target; Let x be the output value of the j-th original objective function; Let x be the corrected target value for the j-th objective. For the weighted target value, This represents the target priority coefficient.

[0230] S425 Layered Screening:

[0231] The original population P(N) and the offspring set Q are merged into a mixed population H. ;

[0232] A "stratified selection" mechanism is adopted, which divides Pareto levels through non-dominated sorting, prioritizing the retention of individuals at the top of the level. The specific formula for level division is as follows:

[0233]

[0234] The individual hierarchy and the reference point between the individual and the optimal association are calculated using formulas (28) and (29). The correlation is determined by the formula (30) to select N individuals from the mixed population H=P∪Q to form a new main population P(N).

[0235] in, express Dominate The smaller the rank number, the higher the individual's priority; the optimal individual has a rank of 1.

[0236]

[0237]

[0238] in, The congestion distance is the distance calculated after the objective function in step S424 has been corrected. .

[0239] S426 Update History Summary:

[0240] Select the first t individuals (t=N / 5) from the new dominant population P(N) that are closest to the reference point, and add them to S;

[0241] If |S|>2N, then remove the redundant solution with the smallest crowding distance, retain at most 2N solutions, and output P(N) and S.

[0242] Step S5, Second Stage:

[0243] When 0.25G ≤ g < 0.85G, proceed to the second stage and compare the proportion of feasible solutions in A(n). and The size of the population is determined by executing Model II or Model III, and the population P(N) is output. The Model II process is as follows: Figure 3 As shown;

[0244] S51 Division of Labor and Assistance: If the proportion of feasible solutions in population A(n) is... Then, the solution is substituted into Model II and the operation is performed. The optimal solution is selected through a division of labor and cooperation method, namely, dual preference of objective and constraint.

[0245] S511: Input the current population P(N), auxiliary population A(n), population size N, and number of objective functions m;

[0246] S512: Calculate the crowding distance of individuals using equations (19) and (20), select the top 30% of individuals from population P according to crowding distance, and randomly select 20% of individuals from auxiliary population A to form the parent set. Use the "adaptive crossover mutation" operator to generate the offspring set Q.

[0247] The operator is a current-technology genetic algorithm, where the adaptive behavior is manifested as follows:

[0248] (1) Crossover probability adjustment: As evolution progresses, the crossover probability is gradually reduced to facilitate the fine development of later solutions and maintain stability.

[0249]

[0250] in, The decay coefficient is 0.1 for population P and 0.05 for population A; Let be the initial crossover probability of population P, set to 0.85; Let be the initial crossover probability of population A, which is set to 0.95.

[0251] (2) Distribution index adjustment: cross-distribution index The variation distribution index increases with evolution, resulting in more similar offspring from later crossover operations; The rate decreases with evolution, resulting in smaller disturbances from later mutations.

[0252]

[0253]

[0254] (3) Mutation probability adjustment: based on the proportion of feasible solutions in the population. The mutation probability is fine-tuned in real time to cope with different search conditions:

[0255] when When the value is less than 0.3, it indicates that the population is in the feasible domain exploration stage. Appropriately enhancing the variation can help the population escape the infeasible area.

[0256]

[0257] when < When the population has converged to the feasible region, it indicates that the variation should be appropriately reduced to avoid destroying the optimal solution structure.

[0258]

[0259] in, Let be the initial mutation probability of population P, which is set to 0.1; Let be the initial mutation probability of population A, which is set to 0.05.

[0260] S513: Introduce a "diversity reward factor" to the objective function value of offspring individuals to modify the objective function. The formula is as follows:

[0261]

[0262] in, Let β be the average distance between individual x and other individuals in the population on the j-th objective dimension, and let β be the reward coefficient.

[0263] A "stratified screening" mechanism was adopted for the merged mixed population. Screening is performed, and the individuals are stratified according to formula (28). Within the same stratum, they are sorted according to "corrected crowding distance". The specific formulas are the same as (19) and (20). The N individuals with larger distances are retained to form a new main population P(N).

[0264] S514: If A(n) is an empty set, select the first 0.3N feasible solutions from the new population P(N) according to the crowding distance of individuals and place them in A(n); if A(n) is not an empty set, check each individual in the new population P(N). If there is a new individual that dominates the existing individuals in A(n), or the crowding distance is greater in the non-dominated layer, then replace the individual with the smallest crowding distance in A(n) in equal numbers to maintain the size of A(n) constant, i.e., n=0.3N, and output the new main population P(N) and auxiliary population A(n).

[0265] S52: If k2>k b Then substitute it into Model III and perform the operation, the specific steps are the same as S42;

[0266] Step S6, the third stage:

[0267] When g≥0.85G, enter the third stage, directly execute Model Ⅲ to calculate and output P(N);

[0268] Step S7, obtain the optimal solution set:

[0269] As the number of iterations g increases, it automatically enters the corresponding stage (S4→S5→S6) until the termination condition is met: the maximum number of iterations G is reached or within a consecutive preset number of generations, the dimensionless change of the optimal solutions of the m objective functions in the population for two consecutive generations is less than a preset threshold, as shown in formulas (32) and (33). Then, the Pareto optimal solution set that satisfies the constraints is output.

[0270]

[0271]

[0272] in, The normalized target value; This represents the original output value of individual X on the k-th objective function; , This represents the historical maximum value of the objective function. This represents the historical minimum value of the objective function. Let be the optimal objective value of the k-th objective function in the t-th generation.

[0273] Example 2:

[0274] A domestic water conveyance project, 2.9 km long, includes a 420 m long and 2.4 m diameter water diversion tunnel. The pumping station is equipped with four horizontal single-stage, double-suction centrifugal pumps, with a total installed capacity of 4 × 2200 kW. The pumps have a rated speed of 1200 r / min, a design head of 182 m, and a design flow rate of 3.20 m³ / s. The entire water pipeline has significant undulations and numerous humps, and is characterized by high head, short distance, and large flow rate. A detailed pipeline layout diagram is shown below. Figure 5 As shown.

[0275] This embodiment discloses a water hammer protection method based on the MSCMOEA optimization algorithm. Applied in this project, it ensures stable operation while achieving multi-objective optimization of the water hammer protection scheme, further ensuring project safety and effectively reducing project costs. The overall algorithm flow is as follows: Figure 1 The specific steps are as follows:

[0276] S1: A pipeline water hammer pressure calculation model is constructed based on fluid mechanics equations, and the method of characteristics is used to solve it, resulting in two pairs of compatible equations of characteristic lines.

[0277] S11: The basic mathematical equations for establishing the mathematical model of a pressurized water pipeline are as follows:

[0278]

[0279]

[0280] In the formula: h is the piezometric head; g is the gravitational acceleration; x is the node coordinate; a is the water hammer wave velocity; α is the pipe inclination angle; v is the flow velocity; t is the time; f is the friction coefficient; and D is the pipe diameter.

[0281] S12: Solving equations (1) and (2) using the method of characteristics yields two pairs of compatible equations for characteristic lines as follows:

[0282]

[0283]

[0284] Cp and CM are the coefficients of the corresponding positive and negative feature lines, respectively:

[0285]

[0286]

[0287] In the formula: , , , B and R are all calculation constants. , Let be the water heads at sections A and B at time t, respectively; let q be the flow rate at the calculated section at time t, satisfying q = vA; and let A be the cross-sectional area. Indicates the length of the feature line grid pipeline.

[0288] S2: Define the decision variables to be optimized, and determine the constraints and objective function:

[0289] S21: Select the surge tank diameter Ds and the initial water depth. Water supply pipe diameter D, pressure regulating tower installation location Quick-closing angle Quick closing time Slow closing angle Slow closing time Air valve position , As a decision variable;

[0290] S22: Using the dimensionless maximum water hammer pressure, the dimensionless minimum water hammer pressure, the dimensionless minimum pump speed, and the water hammer protection cost (pressure regulating tower volume) as optimization objectives, four objective functions are obtained:

[0291]

[0292]

[0293]

[0294]

[0295] In the formula: This represents the maximum water hammer pressure value. This is the minimum water hammer pressure value; This is the dimensionless minimum speed of the pump; The volume of the pressure regulating tower is used as a cost indicator.

[0296] S23: Constraints include:

[0297] S231: Pump speed constraint:

[0298]

[0299] in: t represents the minimum dimensionless rotational speed; t is the time during which the reverse rotational speed is greater than the rated speed.

[0300] S232: Maximum and minimum water hammer pressure constraints:

[0301]

[0302]

[0303] in: This represents the maximum water hammer pressure value. This refers to the rated pressure at the pump outlet. The maximum water hammer pressure value occurring at each node; The maximum allowable pressure for the pipeline at each node; This is the minimum water hammer pressure value; This represents the maximum allowable vacuum level for the pipeline. This represents the minimum water hammer pressure value at each node; This represents the maximum allowable vacuum level for each node's pipeline.

[0304] S233: Pump downstream valve closing constraint:

[0305]

[0306] S234: Input variable upper and lower constraints:

[0307]

[0308] Where Dt and Ht are the diameter of the pressure regulating tower body and the initial effective head, respectively, and D is the diameter of the water supply pipe.

[0309] S235: Safety factor constraint for unidirectional voltage regulating tower:

[0310]

[0311] The final water depth after the unidirectional surge tank has finished replenishing water. This refers to the minimum water depth required for actual engineering projects.

[0312] S3: Input the fitted optimization objective function into the MSCMOEA algorithm, and set the initial parameters: maximum number of iterations G=300, feasible solution ratio judgment threshold. =0.5, =0.85 and the algebraic intervals of each stage, initialize the main population P(N) (N=150), the auxiliary population A(n) (initially an empty set, n=45 after the first substitution into Model II), and the historical solution set S (an empty set before the first substitution into Model III, with a maximum capacity of 300).

[0313] S4: When the number of iterations g < 75, enter the first stage, compare the proportions of feasible solutions k1 and ka in P(N), execute Model I or Model III to perform calculations and output the population P(N);

[0314] S41: If the proportion of feasible solutions in population P(N) Then, substitute the values ​​into Model I and perform the operation:

[0315] S411: Input population P(N), then... Parents are randomly generated from the population P(N), offspring are generated using the GAhalf operator and stored in the set Q, and the original population P(N) and the offspring set Q are merged into a mixed population H (H=P∪Q).

[0316] The GAhalf operator is a current genetic algorithm, and its specific parameters and execution settings are as follows:

[0317] Parent selection strategy: Use tournament selection, tournament size k=3 (each time, 3 individuals are randomly selected from the population without replacement, and the individual with the best fitness is selected as the parent).

[0318] Crossover operation: Simulated binary crossover (SBX) is used, with crossover probability... Cross-distribution index ;

[0319] Mutation operation: Polynomial mutation is used, with mutation probability... (where n is the dimension of the decision variable), variation distribution index ;

[0320] Offspring generation: Offspring are generated through the above selection, crossover, and mutation steps and stored in set Q.

[0321] S412: Perform objective function offset calculation, assuming the individuals in the mixed population H... The original objective function is ( The calculation formula is (17):

[0322]

[0323] in, For reference point The j-th target value;

[0324] The offset coefficient is calculated using formula (18). :

[0325]

[0326] in, To adjust the parameters, g is the current iteration number, G=300.

[0327] S413: Calculate the total crowding distance for an individual using formulas (19) and (20). :

[0328]

[0329]

[0330] in, Let be the crowding distance of individual xi in the same non-dominated layer for the j-th objective. For adjacent individuals, , These are the maximum and minimum values ​​of the target.

[0331] Prioritize retaining the n individuals with the largest total crowding distance and output the new P(N).

[0332] S42: If Then, substitute the values ​​into Model III and perform the operation:

[0333] S421: Input the current population P(N), auxiliary population A (initial size 0, maximum size 45), and historical solution set S (initial size 0, maximum size 45). The objective function has 4 elements, and a preset reference point set is used. .

[0334]

[0335] Where s=2, Total number .

[0336] S422: If |A|≠0, sort by crowding distance D(x) in descending order, and select the first 0.2|A| non-dominated individuals to form the elite parent-child set. If |A|=0, then calculate using formula (20), select the top 30% of individuals with the largest total crowding distance from P(N) to form a temporary A(n) and perform the above operation, while simultaneously adjusting the individual... closest to the reference point Individuals, forming a reference point parent-child set Together they form the final parent-child set. After deduplication, the size is 0.2|A|+K.

[0337] S423: Generate the offspring set Q based on the parent subset Pk using the Gaussian mutation crossover operator.

[0338]

[0339]

[0340]

[0341] in, , The parent individuals are randomly selected, and n is the decision variable. For cross weights, The standard deviation of Gaussian noise. With a mean of 0 and a variance of Let g be a Gaussian random variable, and G be the current iteration number and the maximum iteration number, respectively.

[0342] S424: Introduce a target priority coefficient and correct the target value using formula (25):

[0343]

[0344]

[0345]

[0346] in, Let be the average distance of an individual to the j-th target. For the weighted target value, The target priority coefficient is determined based on project requirements. 5, 3, 2.

[0347] S425: Employs a "stratified selection" mechanism, dividing Pareto levels through non-dominated sorting, prioritizing the retention of individuals at higher levels. The specific formula for level division is as follows:

[0348]

[0349] The individual hierarchy and the reference point between the individual and the optimal association are calculated using formulas (28) and (29). The correlation is determined by the formula (30) to select N individuals from the mixed population H=P∪Q to form a new main population P(N).

[0350] in, express Dominate The smaller the rank number, the higher the individual's priority; the optimal individual has a rank of 1.

[0351]

[0352]

[0353] in, To correct for crowding.

[0354] S426: Select the 30 individuals closest to the reference point from the new main population P(N), add them to S. If |S|>300, remove the redundant solution with the smallest crowding distance. Keep a maximum of 300 solutions and output P(N) and S.

[0355] S5: When 75 ≤ g < 255, proceed to the second stage and compare the feasible solution proportions k2 and k in A(N). b The size of the population is determined by performing calculations using Model II or Model III, and the population P(N) is output.

[0356] S51: If the proportion of feasible solutions in population A(n) is... Then, substitute the values ​​into Model II and perform the operation:

[0357] S511: Input the current population P(N), auxiliary population A(n), population size N, and the number of objective functions is 4;

[0358] S512: Calculate the crowding distance of individuals using equations (19) and (20), select the top 30% of individuals from population P according to crowding distance, and randomly select 20% of individuals from auxiliary population A to form the parent set. Use the "adaptive crossover mutation" operator to generate the offspring set Q.

[0359] The operator is a current-technology genetic algorithm, where the adaptive behavior is manifested as follows:

[0360] (1) Crossover probability adjustment: As evolution progresses, the crossover probability is gradually reduced to facilitate the fine development of later solutions and maintain stability.

[0361]

[0362] in, The decay coefficient is 0.1 for population P and 0.05 for population A; Let be the initial crossover probability of population P, set to 0.85; Let be the initial crossover probability of population A, which is set to 0.95.

[0363] (2) Distribution index adjustment: cross-distribution index The variation distribution index increases with evolution, resulting in more similar offspring from later crossover operations; The rate decreases with evolution, resulting in smaller disturbances from later mutations.

[0364]

[0365]

[0366] (3) Mutation probability adjustment: based on the proportion of feasible solutions in the population. The mutation probability is fine-tuned in real time to cope with different search conditions:

[0367] when When the value is less than 0.3, it indicates that the population is in the feasible domain exploration stage. Appropriately enhancing the variation can help the population escape the infeasible area.

[0368]

[0369] when < When the population has converged to the feasible region, it indicates that the variation should be appropriately reduced to avoid destroying the optimal solution structure.

[0370]

[0371] in, Let be the initial mutation probability of population P, which is set to 0.1; Let be the initial mutation probability of population A, which is set to 0.05.

[0372] S513: Introduce a "diversity reward factor" to the objective function value of offspring individuals to modify the objective function. The formula is as follows:

[0373]

[0374] in Let β be the average distance between individual x and other individuals in the population on the j-th objective dimension, and let β be the reward coefficient.

[0375] The “stratified screening” mechanism is adopted. The stratification is carried out according to formula (28). Within the same stratum, the sorting is based on “corrected crowding distance”. The specific formulas are the same as (19) and (20). 150 individuals with larger distances are retained.

[0376] S514: If A(n) is an empty set, select the first 45 feasible solutions from the new population P(N) according to the crowding distance of individuals and place them in A(n); if A(n) is not an empty set, there exists a new individual that dominates the existing individuals in A, or the crowding distance is greater in the non-dominated layer, then replace the individual with the smallest crowding distance in A in equal numbers, maintain the size of A constant, i.e., n=45, and output the new main population P(N) and the helper population A(n).

[0377] S52: If k2 > kb, then substitute it into Model III and perform the operation. The specific steps are the same as in S42.

[0378] S6: When g≥255, enter the third stage, directly execute Model Ⅲ to calculate and output P(N);

[0379] S7: As the number of iterations g increases, it automatically enters the corresponding stage (S4→S5→S6) until the termination condition is met: the number of iterations reaches 300 or within a consecutive preset number of generations, the dimensionless change of the optimal solution of the four objective functions in the population for two consecutive generations is less than the preset threshold, see formula (32) (33) for details, then output the Pareto optimal solution set that satisfies the constraints.

[0380]

[0381]

[0382] in, , This represents the historical maximum value of the objective function. This is the historical minimum value of the objective function.

[0383] The above description is only used to illustrate the technical solutions of the present invention and is 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 (such as the application of various formulas, the order of steps, etc.) without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A water hammer protection method based on a multi-preference constrained multi-objective optimization algorithm, characterized in that, Includes the following steps: Step S1, construct the pipeline water hammer pressure calculation model: A pipeline water hammer pressure calculation model was constructed based on fluid mechanics equations. The method of characteristics was used to solve the model, resulting in two pairs of compatible equations of characteristic lines. The initial model establishment was completed. Step S2: Define the decision variables to be optimized and determine the constraints and objective function. Select the surge tank diameter Ds and the initial water depth Water supply pipe diameter D, pressure regulating tower installation location Quick-closing angle Quick closing time Slow closing angle Slow closing time Air valve position , As a decision variable; The constraints include: pump speed constraints, maximum water hammer, minimum water hammer pressure constraints, valve closure constraints after pump, input variable up and down constraints, and unidirectional pressure regulating tower safety factor constraints. The objective function is obtained by taking the dimensionless maximum water hammer pressure, the dimensionless minimum water hammer pressure, the dimensionless minimum speed of the pump, and the water hammer protection cost as optimization objectives. Step S3: Input the fitted optimized objective function into the MSCMOEA algorithm and set the initial parameters: Step S4, First Stage: Compare the ratio of the previous iteration count g to the preset total iteration count G. When g < 0.25G, proceed to the first stage. Obtain the proportion of feasible solutions for the main population P(N). The feasible solution is proportional Compared with the preset benchmark Perform size comparison; determine the target computational model based on the comparison results; if Then, Model I is selected as the target computational model, and the early exploration of the optimal solution is completed around the target preference; if Then, Model III is selected as the target calculation model; the calculation operation is performed to complete the fine optimization of the optimal solution, with constraint preference as the main factor and target as the auxiliary factor, and the calculated P(N) is output. Step S5, Second Stage: When 0.25G ≤ g < 0.85G, the second stage begins; Obtain the feasible solution ratio k2 of the auxiliary population A(n), and compare the feasible solution ratio k2 with the preset benchmark k. b Perform size comparison; determine the target computational model based on the comparison results, if Then, Model II is selected as the target calculation model. Through a collaborative approach, i.e., a dual preference for both the target and constraints, the optimal solution is chosen. If k2 > k b Then, Model III is selected as the target calculation model; the calculation operation is performed to complete the fine optimization of the optimal solution, and the calculated P(N) is output. Step S6, the third stage: When g≥0.85G, enter the third stage, directly execute Model Ⅲ to calculate and output P(N), and complete the fine optimization of the optimal solution; Step S7, obtain the optimal solution set: Determine whether the iteration termination condition is met. The termination condition includes two items: first, the number of iterations reaches a preset total number; second, within a preset number of generations, the dimensionless change of the optimal solution of the objective function in the population for two consecutive generations is less than a preset threshold. When either termination condition is met, output the Pareto optimal solution set that satisfies the constraints.

2. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S1, the fluid dynamics equations, including the continuity equation and the momentum equation, are calculated using the following formulas: In the formula, h is the piezometric head; g is the acceleration due to gravity, 9.81 m / s². 2 x represents the node coordinates; The water hammer wave velocity; t is the pipe inclination angle; v is the flow velocity; t is time; f is the friction factor; D is the pipe diameter in meters.

3. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S1, the characteristic line compatibility equation includes a forward characteristic line compatibility equation and a reverse characteristic line compatibility equation; The formula for calculating the compatibility equation of the positive characteristic lines is as follows: The formula for calculating the compatibility equation of reverse characteristic lines is as follows: in, and These are the coefficients for the corresponding positive and negative feature lines, respectively, and are calculated using the following formulas: In the formula: , , , B and R are all calculation constants. , Let be the water heads at sections A and B at time t, respectively; let q be the flow rate at the calculated section at time t, satisfying q = vA; and let A be the cross-sectional area. Indicates the length of the feature line grid pipeline.

4. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S2, the objective function is calculated using the following formula: In the formula, This represents the maximum water hammer pressure value. This is the minimum water hammer pressure value; This is the dimensionless minimum speed of the pump; The volume of the pressure regulating tower is used as a cost indicator.

5. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S2, the pump speed constraint is calculated using the following formula: in, The minimum dimensionless rotational speed; This refers to the time during which the reverse rotation speed is greater than the rated speed. The maximum and minimum water hammer pressure constraints are calculated using the following formulas: in: This represents the maximum water hammer pressure value. This refers to the rated pressure at the pump outlet. The maximum water hammer pressure value occurring at each node; The maximum allowable pressure for the pipeline at each node; This is the minimum water hammer pressure value; This represents the maximum allowable vacuum level for the pipeline. This represents the minimum water hammer pressure value at each node; The maximum allowable vacuum level for each node pipeline; The calculation formula for the valve closing constraint after the pump is as follows: in, This refers to the valve closing angle during the rapid closing phase. This refers to the valve closing angle during the slow closing phase. To achieve the fast closing angle Time required; To achieve the slow closing angle Time required; The upper and lower constraints of the input variables are calculated using the following formula: in, The diameter of the voltage regulating tower, The initial effective head is given by D, which is the diameter of the make-up water pipe. This means that the diameter of the make-up water pipe must be greater than the minimum allowable value and less than the diameter of the surge tank body. ; The safety factor constraint for the unidirectional voltage regulating tower is calculated using the following formula: The final water depth after the unidirectional surge tank has finished replenishing water. This refers to the minimum water depth required for actual engineering projects.

6. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S3, the initial parameters include: the maximum number of iterations G, and the feasible solution ratio judgment threshold. Initialize the primary population P(N), the auxiliary population A(n), and the historical solution set S.

7. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S4, model I includes the following steps: S411, Population Evolution: Given a population P(N), for Parents are randomly generated from the population P(N), and offspring are generated using the GAhalf operator and stored in set Q. The original population P(N) and the offspring set Q are then merged into a mixed population H. ; S412, Objective function correction: Perform objective function offset calculation, assuming individuals in the mixed population H The original objective function is , The objective function after offset is calculated as follows: in, For reference point The j-th target value; The offset coefficient is dynamically adjusted based on the proportion of feasible solutions in the mixed population H, and is calculated using the following formula: in, To adjust the parameters, g is the current iteration number, and G is the maximum iteration number; S413, Crowding Filter: Calculate individuals within the same non-dominated layer in a mixed population H. Total congestion distance The formula and parameter definitions are as follows: in, Individuals in the same non-dominated layer The crowding distance of the j-th target. For adjacent individuals, For individuals The function value of the j-th objective function; , These are the maximum and minimum values ​​of the target. Perform non-dominated ordination on the mixed population H, calculate the crowding distance of individuals within the same non-dominated layer, sort them from high to low non-dominated level and from large to small crowding distance within the same layer, and screen to obtain a new generation population P(N) of size N.

8. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S4, model III includes the following steps: S421, Initialization and Input Parameters: Input the current population P(N), auxiliary population A(n), historical solution set S, number of objective functions m, and preset reference point set. ; The auxiliary population A(n) has an initial size of 0 and a maximum size of 0.3N; The initial size of the historical solution set S is 0, and the size of the historical solution set S is... ; The formula for calculating the reference point is as follows: Where s=2, Total number of reference points ; S422, Parent Selection: If the size of the auxiliary population ; Select non-dominant individuals from the auxiliary population A(n), sort them in descending order by crowding distance D(x), and select the top... A group of non-dominant individuals forms an elite parent-child set. ; If the size of the auxiliary population Then, select the top 30% of individuals with the largest crowding distance from P(N) to form a temporary A(n) for the above operation, and simultaneously perform the above operation on the individuals. closest to the reference point Individuals, forming a reference point parent-child set The size is K; Generate parent child set: Deduplication retains individuals with larger crowding distances, at a scale of ; S423, Descendant Generation: Based on parent-child set The population offspring set Q is generated using the Gaussian mutation crossover operator, including: From father to son Two different parent individuals are randomly selected. and Generate offspring individuals The calculation formula is as follows: in, , The parent individuals are randomly selected, and n is the decision variable. For cross weights, The standard deviation of Gaussian noise. With a mean of 0 and a variance of Let g be a Gaussian random variable, and G be the current iteration number and the maximum iteration number, respectively. Each generated descendant individual Add to the descendant set Q; S424, Objective function correction: A target priority coefficient is introduced to adjust the target value. The calculation formula is as follows: in, Let x be the average distance of individual x to the j-th target; Let x be the output value of the j-th original objective function; Let x be the corrected target value for the j-th objective. For the weighted target value, This represents the target priority coefficient. S425, Layered Filtering: The original population P(N) and the offspring set Q are merged into a mixed population H. ; Pareto hierarchies are defined using non-dominated sorting, prioritizing the retention of individuals at higher levels. The specific formula for hierarchical division is as follows: in, express Dominate The smaller the rank number, the higher the individual's priority; the optimal individual has a rank of 1. Individual and optimal correlation reference point The correlation degree is calculated using the following formula: The overall score is calculated using the following formula: in, To correct for crowding; Sort individuals by rank in ascending order, prioritize retaining individuals with higher rank in the mixed population H, sort individuals in the same rank by comprehensive score in descending order, and select the top N individuals to form a new main population P(N). S426, Updated historical solution set: From the new dominant population P(N), for each reference point Select the closest individual, and select a total of t unique individuals, t=N / 5, and add them to the historical solution set S; If |S|>2N, then remove the redundant solution with the smallest crowding distance, retain at most 2N solutions, and output P(N) and S.

9. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S5, model II includes the following steps: S511, Initialize input parameters: Input the current population P(N), auxiliary population A(n), population size N, and number of objective functions m; S512, Parent Set and Descendant Generation: Individual crowding distance is calculated using formulas (19) and (20). Sort individuals from largest to smallest by crowding distance, select the top 30% of individuals in population P(N), and randomly select 20% of individuals from auxiliary population A(n) to form the parent set. Use the "adaptive crossover mutation" operator to generate the offspring set Q. S513 Objective Function Modification and Hierarchical Filtering: A "diversity reward factor" is introduced into the objective function value of offspring individuals to modify the objective function. The calculation formula is as follows: in, Let be the average distance between individual x and other individuals in the population on the j-th target dimension, and β be the reward coefficient with a value of 0.

2. A "stratified screening" mechanism was adopted for the merged mixed population. Screening is performed; stratification is performed according to formula (28), and individuals in the same stratum are sorted from largest to smallest according to "corrected crowding distance", with the same formula as (19) and (20). The N individuals with larger distances are retained to form a new main population P(N); S514, Auxiliary Population Update and Result Output: If A(n) is an empty set, then sort the individuals in descending order of their crowding distance, and select the top 0.3N feasible solutions from the new principal population P(N) and place them in A(n), with a size of n=0.3N; If A(n) is not an empty set, then examine each individual in the new population P(N). If a new individual dominates an existing individual in A(n), or if the crowding distance is greater in the non-dominated layer, then replace the individual with the smallest crowding distance in A(n) in equal numbers to maintain the size of A(n) constant, i.e., n=0.3N. Output the new primary population P(N) and secondary population A(n).

10. The water hammer protection method based on a multi-preference constraint multi-objective optimization algorithm according to claim 1, characterized in that, In step S7, the dimensionless change in the optimal solution of the objective function in the population for two consecutive generations is less than a preset threshold. The calculation formula is as follows: in, The normalized target value; This represents the original output value of individual X on the k-th objective function; , This represents the historical maximum value of the objective function. This represents the historical minimum value of the objective function. Let be the optimal objective value of the k-th objective function in the t-th generation.

Images