A hydraulic pipe network simulation system

Abstract

The application discloses a hydraulic pipe network simulation system, and relates to the technical field of pipe network simulation, which comprises the following steps: establishing a connection relation graph of nodes and pipe sections; determining a peak array, and deriving an augmented incidence matrix, a basic incidence matrix, a basic loop matrix, a tree branch matrix, a residual branch matrix, a tree branch vector and a residual branch pipe section vector; inputting corresponding parameters; adopting a basic loop method optimized by a magnetotactic bacteria method to perform first layer iteration, and calculating a residual branch pipe section flow increment iteration step; adopting an LM algorithm, and combining with the magnetotactic bacteria method to perform second layer iteration and calculate an end point flow iteration step; adopting an LM algorithm modified by a trust region, and combining with the magnetotactic bacteria method to perform third layer iteration and calculate a controlled variable iteration step; and calculating and outputting hydraulic properties of each node and pipe section. The application solves the problem that the existing iteration algorithm causes singular values due to too many control points or numerical solutions cannot converge due to improper initial value setting.

Description

Technical Field

[0001] This invention relates to the field of pipeline network simulation technology, and in particular to a hydraulic pipeline network simulation system. Background Technology

[0002] Pipelines are ubiquitous in modern life and industrial production. Urban water supply networks, heating networks, and chemical fluid transportation networks are particularly large in scale, complex in process, and difficult to operate and regulate. Pipeline simulation systems use actual pipelines as their model, simulating actual operations through abstract modeling and appropriate algorithms to predict operational results, optimize process parameters, and avoid risks associated with on-site operations. The core of pipeline simulation lies in the establishment of mathematical models and algorithms. The quality of mathematical modeling directly affects the realism and rationality of the simulation, while the quality of the algorithm directly relates to the operating speed, stability, and accuracy of the entire simulation system.

[0003] Extensive research has been conducted both domestically and internationally on modeling and algorithms for pipeline systems, and numerous system simulation software programs have been developed, such as those for heating and water supply systems. However, existing simulation systems suffer from the following problems: iterative algorithms are prone to singular values ​​when there are many control points; improper initial value settings can lead to non-convergence of numerical solutions; the accuracy of mechanical energy balance calculations at the inlet and outlet ends of the pipeline network is insufficient; and the solution efficiency and stability for nonlinear problems in complex pipeline systems are low. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention provides a hydraulic pipeline network simulation system, comprising the following steps:

[0005] S1. Define the actual pipeline network as a pipeline network model including pipe segments and nodes. The characteristic attributes of the pipe segments include pipe segment structural attributes, pipe segment topology attributes, and pipe segment hydraulic attributes. The characteristic attributes of the nodes include node structural attributes, node topology attributes, and node hydraulic attributes. Establish the connection relationship diagram between nodes and pipe segments, and assign numbers to the nodes and pipe segments.

[0006] S2. Determine the peak array F based on the node and pipe segment topology attributes;

[0007] S3. Derive the augmented correlation matrix based on the peak matrix F. According to the augmented correlation matrix Derive the basic correlation matrix Basic circuit matrix Tree matrix , redundant branch matrix Tree branch vector T, and remaining branch segment vector L;

[0008] S4. Input pipe segment and node characteristic attributes to determine the endpoint pressure value;

[0009] S5. If automatic control points are set, then set the operating variables in the automatic control system. The initial values ​​and setpoints of the controlled variables are used to calculate the impedance of each pipe segment based on its characteristic properties. If there are no automatic control points, the impedance of each pipe segment is calculated directly based on its characteristic properties. ;

[0010] S6. If inlet and outlet endpoints are set, first input the initial value of the endpoint flow rate, and then input the flow rate of the remaining branch pipe section. The initial value; if there are no inlet or outlet endpoints, directly input the flow rate of the remaining branch pipe section. The initial value;

[0011] S7. Calculate the flow rate of the branch pipe section. The resistance loss and pump head of each pipe section are calculated based on the characteristics of the pipe section.

[0012] S8. Based on the mass conservation rule, the basic loop method optimized by the magnetotactic bacteria method is used for the first layer of iteration, and the iteration step size of the flow increment of the remaining branch pipe section is calculated. ;

[0013] S9, if If the accuracy requirements are met, proceed to the next process; if If the accuracy requirement is not met, then let = , = + and will Return to S6, and repeat S6-S9 until... To meet the accuracy requirements, calculate the flow rate values ​​for all pipe sections;

[0014] S10. If there are import / export endpoints, proceed to process S11; if there are no import / export endpoints, proceed directly to process S13.

[0015] S11. Based on the endpoint mechanical energy balance rule, the Levenberg-Marquardt algorithm is adopted, combined with the magnetotactic bacteria method to assist in the optimization of the second layer of iteration, and the endpoint flow iteration step size is calculated. ;

[0016] S12. If the endpoint mechanical energy balance meets the accuracy requirements, proceed to the next step; if the endpoint mechanical energy balance does not meet the accuracy requirements, then... = , ,and Return to S6 and repeat S6-S12 until the endpoint mechanical energy balance meets the accuracy requirements;

[0017] S13. If there is an automatic control point, proceed to the next process; if there is no automatic control point, proceed directly to process S16.

[0018] S14. The Levenberg-Marquardt algorithm with trust region correction is used, combined with the magnetotactic bacteria method to assist in optimization for the third layer of iteration, and the iteration step size of the operation variables is calculated. ;

[0019] S15. If the value of the controlled variable meets the accuracy requirements, proceed to the next step; if the value of the controlled variable does not meet the accuracy requirements, then... = ,and Return to S5 and repeat S5-S15 until the value of the controlled variable meets the accuracy requirements;

[0020] S16. Calculate and output the pressure at each node and the flow rate of each pipe section.

[0021] Furthermore, when it is a ring network, the net flow rate at each node is 0; when it is an open-loop network, the initial flow rates at the inlet and outlet endpoints are input.

[0022] For a ring network, based on the node continuity equation and the loop mechanical energy balance equation, the following hydraulic calculation mathematical model is established:

[0023]

[0024] in: This is the column vector of flow rates in the pipe section; This is the column vector of pressure drop in the pipe section; This is the column vector of pipe segment impedance; This is the column vector of the potential energy difference between two nodes in a pipe segment branch; This is the head vector of the water pumps in the pipe section;

[0025] Calculate the flow rate of the branch pipe section The calculation formula is as follows:

[0026] ;

[0027] but

[0028] but

[0029] Calculate the flow rate of the branch pipe section The flow rate of the remaining pipe sections is used to calculate the resistance loss and pump head of each pipe section based on the characteristics of the pipe sections.

[0030] For open-loop pipe networks, the following hydraulic calculation mathematical model is established based on the nodal continuity equation, the loop mechanical energy balance equation, and the end-point mechanical energy balance equation:

[0031]

[0032] in, , All are node codes, and ≠ ; , They are respectively , Node pressure; For nodes arrive Pressure drop across each pipe segment along the path, , They are nodes , The flow rate; Density; Net flow of the node;

[0033] Calculate the flow rate of the branch pipe section The calculation formula is as follows:

[0034] ;

[0035] but

[0036] but

[0037] Calculate the flow rate of the branch pipe section The flow rate of the remaining pipe sections is used to calculate the resistance loss and pump head of each pipe section based on the characteristics of the pipe sections.

[0038] Furthermore, the iteration step size of the flow increment in the remaining branch pipe section is calculated. Includes the following steps:

[0039] S81. Optimize the initial flow rate of the remaining branch pipe section using the magnetotactic bacteria method:

[0040] Magnetotactic movement phase: Simulates bacteria moving along the direction of the magnetic field. The position of the bacteria corresponds to the flow rate value of the candidate pipe segment. The search direction is adjusted according to the flow rate gradient of the pipe segment.

[0041] =

[0042] in, For the first During the nth iteration, the 1st The motion vector of a bacterium, Pipeline flow rate value The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5;

[0043] Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation.

[0044] = +

[0045] in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, no. The bacteria in the first The current pipe flow rate value at the next iteration. For the first The optimal pipe flow rate value at the next iteration;

[0046] Local search adjustment phase: The search is adjusted for the new solution to improve accuracy;

[0047] =

[0048] in, These are random numbers that follow a standard normal distribution (mean 0, standard deviation 1).

[0049] S82. Establish the formula for calculating the iterative step size of the flow increment in the residual branch pipe section using the basic loop method:

[0050]

[0051] in, For weighting factors, 0.1 ≤ ≤0.3, For the pipe segment flow rate in the current iteration Output search correction values.

[0052] Furthermore, calculate the endpoint flow iteration step size. The process includes the following steps:

[0053] S111. Establish the mechanical energy balance formula at the endpoints:

[0054]

[0055] in, Let the endpoint mechanical energy function be denoted as .

[0056] S112. Optimize using the Levenberg-Marquardt algorithm combined with the magnetotactic bacteria method, and calculate the endpoint flow iteration step size. This includes the following steps:

[0057] S1121. Optimize initial endpoint flow using the magnetotactic bacteria method:

[0058] Magnetotactic motion phase: Simulates bacteria moving along the direction of the magnetic field. The position of the bacteria corresponds to the candidate endpoint flux value. The search direction is adjusted according to the gradient of the endpoint flux value.

[0059] =

[0060] in, For the first During the nth iteration, the 1st The motion vector of a bacterium, For the first During the nth iteration, the 1st endpoint traffic value The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5;

[0061] Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation.

[0062] = +

[0063] in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, For the first The bacteria in the first The current endpoint traffic value at the next iteration. For the first The optimal endpoint flow value at the next iteration;

[0064] Local search adjustment phase: The search is adjusted for the new solution to improve accuracy;

[0065] =

[0066] A candidate set of initial values ​​for endpoint flow is generated using the magnetotactic bacteria method, and then selected to make... The smallest initial value is used as the starting point for iteration.

[0067] S1122, Calculate the endpoint flow iteration step size :

[0068]

[0069] in: This is a Jacobian matrix; To correct the parameters, >0, for and convex combination; It is the identity matrix; For weighting factors, 0.1 ≤ ≤0.3, For the endpoint traffic value of the current iteration Output search correction values.

[0070] Further, calculate the values ​​of the operated variables. Iteration step size Includes the following steps:

[0071] S141. Establish values ​​for operand variables The objective function formula is as follows:

[0072] ;

[0073] in: , The vector of deviations between the controlled variable values ​​and the controlled variable setpoints; r - the number of controlled variables;

[0074] , This represents the vector of actual values ​​of the controlled variable under the values ​​of the manipulated variable.

[0075] , Set values ​​for the controlled variable;

[0076] : Manipulate the vector of variable values; l - Number of operation variables;

[0077] Multi-point automatic control is transformed into solving the following nonlinear least squares problem:

[0078] ;

[0079] S142. The Levenberg-Marquardt algorithm with trust region correction is adopted, and the magnetotactic bacteria method is combined to assist in the optimization calculation of the operation variable values. Iteration step size This includes the following steps:

[0080] S1421. Optimize the initial values ​​of the manipulated variables using the magnetotactic bacteria method:

[0081] Magnetotactic motion phase: Simulating bacterial movement along the magnetic field direction, the bacterial position corresponds to the candidate operation variable value, and the search direction is adjusted according to the gradient of the operation variable value.

[0082] =

[0083] in, For the first During the nth iteration, the 1st The motion vector of a bacterium, For the first During the nth iteration, the 1st Values ​​of each operation variable The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5;

[0084] Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation.

[0085] = +

[0086] in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, No. The bacteria in the first The current value of the operation variable at the next iteration. For the first The optimal values ​​of the operation variables at the next iteration;

[0087] Local search adjustment phase: The search is adjusted for the new solution to improve accuracy;

[0088] =

[0089] A candidate set of initial values ​​for the operands is generated using the magnetotactic bacteria method, and then selected to make the initial values ​​more suitable for the operands. The smallest initial value is used as the starting point for iteration;

[0090] S1422, Calculate the values ​​of the operands. Iteration step size :

[0091]

[0092]

[0093] in For Jacobi matrix, To correct the parameters, >0, for and convex combination; It is the identity matrix; For weighting factors, 0.1 ≤ ≤0.3, To use the magnetotactic bacteria algorithm for the operation variables of the current iteration The output search optimization value; For the operation variables of the current iteration The trust region correction value is determined by the actual decrease. Compared with the estimated decrease ratio Make dynamic adjustments. Correction factor:

[0094]

[0095]

[0096]

[0097] S143, Adaptive Adjustment: Used to determine whether to accept the step size And adjusting the correction factor in the iteration parameters. Size; The larger the value, the better the objective function. The more it falls, the more it is accepted. Hoping for the next tentative step Longer, therefore smaller ;on the contrary, The smaller, the more likely it is to be rejected. Increase ;

[0098] The process is as follows:

[0099] Step 1: Given , ≥0, , M - the minimum value of the correction factor set. ,

[0100] Step 2: If If the system fails, it will shut down and proceed to step S12 for further judgment. Does it meet the accuracy requirements? Otherwise, proceed to step 3.

[0101] Step 3: Select Expression, and calculate ;

[0102]

[0103] Step 4: Calculation ;Set:

[0104]

[0105] Step 5: Select according to the following rules :

[0106]

[0107] Step 6: Set: Proceed to step 2.

[0108] The beneficial effects of this invention compared with the prior art are: (1) This invention can realize the global exploration of the flow solution space through the magnetotactic movement in the magnetotactic bacteria algorithm, thereby avoiding missing the optimal solution; through reproduction and extinction and local search, the optimization of high-quality flow solutions is realized, thereby improving the simulation accuracy; (2) This invention introduces the trust region correction Levenberg-Marquardt algorithm in combination with the magnetotactic bacteria algorithm, combined with the three-layer iterative framework, to solve the problem that the iterative algorithm is prone to singular values ​​when there are many control points and that improper initial value setting will lead to the non-convergence of numerical solutions; (3) This invention solves the problems of insufficient accuracy of mechanical energy balance of the inlet and outlet ends of the pipeline network and low efficiency and poor stability of solving nonlinear problems of complex pipeline network systems in the prior art, thereby improving the accuracy and stability of hydraulic pipeline network simulation. Attached Figure Description

[0109] Figure 1 This is a flowchart of the hydraulic pipeline network simulation system of the present invention.

[0110] Figure 2 This is a process flow diagram of the ring network of the present invention.

[0111] Figure 3 This is a topology diagram of the ring network nodes of the present invention.

[0112] Figure 4 This is a process flow diagram of the open-loop pipeline network of the present invention.

[0113] Figure 5 This is a topology diagram of the open-loop pipeline network nodes of the present invention.

[0114] Figure 6 This is a diagram of the user interface of the simulation system of the present invention. Detailed Implementation

[0115] The present invention will be further described below with reference to specific embodiments. The illustrative embodiments and descriptions herein are used to explain the present invention, but are not intended to limit the present invention.

[0116] Example: Figures 1-6As shown, the present invention provides a hydraulic pipeline network simulation system, comprising the following processes:

[0117] S1. Define the actual pipeline network as a pipeline network model that includes two types of elements: pipe segments and nodes. The characteristic attributes of pipe segments include pipe segment structural attributes, pipe segment topological attributes, and pipe segment hydraulic attributes. The characteristic attributes of nodes include node structural attributes, node topological attributes, and node hydraulic attributes. Establish the connection relationship diagram of nodes and pipe segments, and assign numbers to nodes and pipe segments. Based on the pipeline characteristic attributes, the impedance of straight pipe segments, valves, and fittings in the pipe segment and the head of the conveying equipment can be calculated. Based on the node characteristic attributes, the potential energy difference between two nodes in the pipe segment can be calculated.

[0118] The structural properties of pipe segments and nodes are determined through the design of the pipeline system or the actual installation status. The topological properties of pipe segments and nodes are expressed using mathematical graph theory, while the hydraulic properties of pipe segments and nodes are analyzed and calculated using hydraulic theory.

[0119] Pipeline structural attributes include: pipe length, diameter, roughness, and parameters such as conveying equipment, valves, and fittings contained in the pipe section; pipe topology attributes include: pipe flow direction, starting node, and ending node; pipe hydraulic attributes include: pipe flow rate, pipe velocity, pipe head, and pipe resistance loss.

[0120] Node construction attributes include: node height; node topology attributes include: the pipe segment associated with the node; node hydraulic attributes include: node flow rate and node pressure.

[0121] S2. Determine the peak array F based on the node and pipe segment topology attributes;

[0122] S3. Derive the augmented correlation matrix based on the peak matrix F. According to the augmented correlation matrix Derive the basic correlation matrix Basic circuit matrix Tree matrix , redundant branch matrix Tree branch vector T, and remaining branch segment vector L;

[0123] S4. Input pipe segment and node characteristic attributes to determine the endpoint pressure value;

[0124] S5. If automatic control points are set, then set the operating variables in the automatic control system. The initial values ​​and setpoints of the controlled variables are used to calculate the impedance of each pipe segment based on its characteristic properties. If there are no automatic control points, the impedance of each pipe segment is calculated directly based on its characteristic properties. ;

[0125] S6. If inlet and outlet endpoints are set, first input the initial value of the endpoint flow rate, and then input the flow rate of the remaining branch pipe section. The initial value; if there are no inlet or outlet endpoints, directly input the flow rate of the remaining branch pipe section. The initial value;

[0126] S7. Calculate the flow rate of the branch pipe section. The resistance loss and pump head of each pipe section are calculated based on the characteristics of the pipe section.

[0127] When it is a ring network, the net flow rate of each node is 0; when it is an open-loop network, input the initial flow rate values ​​at the inlet and outlet endpoints.

[0128] For a ring network, based on the node continuity equation and the loop mechanical energy balance equation, the following hydraulic calculation mathematical model is established:

[0129]

[0130] in: This is the column vector of flow rates in the pipe section; This is the column vector of pressure drop in the pipe section; This is the column vector of pipe segment impedance; This is the column vector of the potential energy difference between two nodes in a pipe segment branch; This is the head vector of the water pumps in the pipe section;

[0131] Calculate the flow rate of the branch pipe section The calculation formula is as follows:

[0132] ;

[0133] but

[0134] but

[0135] Calculate the flow rate of the branch pipe section The flow rate of the remaining pipe sections is used to calculate the resistance loss and pump head of each pipe section based on the characteristics of the pipe sections.

[0136] For open-loop pipe networks, the following hydraulic calculation mathematical model is established based on the nodal continuity equation, the loop mechanical energy balance equation, and the end-point mechanical energy balance equation:

[0137]

[0138] in, , All are node codes, and ≠ ; , They are respectively , Node pressure; For nodes arrive Pressure drop across each pipe segment along the path, , They are nodes , The flow rate; Density; Net flow of the node;

[0139] Calculate the flow rate of the branch pipe section The calculation formula is as follows:

[0140] ;

[0141] but

[0142] but

[0143] Calculate the flow rate of the branch pipe section The flow rate of the remaining pipe sections is used to calculate the resistance loss and pump head of each pipe section based on the characteristics of the pipe sections.

[0144] S8. Based on the mass conservation rule, the basic loop method optimized by the magnetotactic bacteria method is used for the first layer of iteration, and the iteration step size of the flow increment of the remaining branch pipe section is calculated. This includes the following steps:

[0145] S81. Optimize the initial flow rate of the remaining branch pipe section using the magnetotactic bacteria method:

[0146] Magnetotactic movement phase: Simulates bacteria moving along the direction of the magnetic field. The position of the bacteria corresponds to the flow rate value of the candidate pipe segment. The search direction is adjusted according to the flow rate gradient of the pipe segment.

[0147] =

[0148] in, For the first During the nth iteration, the 1st The motion vector of a bacterium, Pipeline flow rate value The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5;

[0149] Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation.

[0150] = +

[0151] in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, no. The bacteria in the first The current pipe flow rate value at the next iteration. For the first The optimal pipe flow rate value at the next iteration;

[0152] Local search adjustment phase: The search is adjusted for the new solution to improve accuracy;

[0153] =

[0154] in, These are random numbers that follow a standard normal distribution (mean 0, standard deviation 1).

[0155] S82. Establish the formula for calculating the iterative step size of the flow increment in the residual branch pipe section using the basic loop method:

[0156]

[0157] in, For weighting factors, 0.1 ≤ ≤0.3, For the pipe segment flow rate in the current iteration Output search correction values.

[0158] S9, if If the accuracy requirements are met, proceed to the next process; if If the accuracy requirement is not met, then let = , = + and will Return to S6, and repeat S6-S9 until... To meet the accuracy requirements, calculate the flow rate values ​​for all pipe sections;

[0159] S10. If there are import / export endpoints, proceed to process S11; if there are no import / export endpoints, proceed directly to process S13.

[0160] S11. Based on the endpoint mechanical energy balance rule, the Levenberg-Marquardt algorithm is adopted, combined with the magnetotactic bacteria method to assist in the optimization of the second layer of iteration, and the endpoint flow iteration step size is calculated. This includes the following steps:

[0161] S111. Establish the mechanical energy balance formula at the endpoints:

[0162]

[0163] in, Let the endpoint mechanical energy function be denoted as .

[0164] S112. Optimize using the Levenberg-Marquardt algorithm combined with the magnetotactic bacteria method, and calculate the endpoint flow iteration step size. This includes the following steps:

[0165] S1121. Optimize initial endpoint flow using the magnetotactic bacteria method:

[0166] Magnetotactic motion phase: Simulates bacteria moving along the direction of the magnetic field. The position of the bacteria corresponds to the candidate endpoint flux value. The search direction is adjusted according to the gradient of the endpoint flux value.

[0167] =

[0168] in, For the first During the nth iteration, the 1st The motion vector of a bacterium, For the first During the nth iteration, the 1st endpoint traffic value The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5;

[0169] Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation.

[0170] = +

[0171] in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, For the first The bacteria in the first The current endpoint traffic value at the next iteration. For the first The optimal endpoint flow value at the next iteration;

[0172] Local search adjustment phase: The search is adjusted for the new solution to improve accuracy;

[0173] =

[0174] A candidate set of initial values ​​for endpoint flow is generated using the magnetotactic bacteria method, and then selected to make... The smallest initial value is used as the starting point for iteration.

[0175] S1122, Calculate the endpoint flow iteration step size :

[0176]

[0177] in: This is a Jacobian matrix; To correct the parameters, >0, for and convex combination; It is the identity matrix; For weighting factors, 0.1 ≤ ≤0.3, For the endpoint traffic value of the current iteration Output search correction values.

[0178] S12. If the endpoint mechanical energy balance meets the accuracy requirements, proceed to the next step; if the endpoint mechanical energy balance does not meet the accuracy requirements, then... = , ,and Return to S6 and repeat S6-S12 until the endpoint mechanical energy balance meets the accuracy requirements;

[0179] S13. If there is an automatic control point, proceed to the next process; if there is no automatic control point, proceed directly to process S16.

[0180] S14. The Levenberg-Marquardt algorithm with trust region correction is used, combined with the magnetotactic bacteria method to assist in optimization for the third layer of iteration, and the iteration step size of the operation variables is calculated. This includes the following steps:

[0181] S141. Establish values ​​for operand variables The objective function formula is as follows:

[0182] ;

[0183] in: , The vector of deviations between the controlled variable values ​​and the controlled variable setpoints; r - the number of controlled variables;

[0184] , This represents the vector of actual values ​​of the controlled variable under the values ​​of the manipulated variable.

[0185] , Set values ​​for the controlled variable;

[0186] : Manipulate the vector of variable values; l - Number of operation variables;

[0187] Multi-point automatic control is transformed into solving the following nonlinear least squares problem:

[0188] ;

[0189] S142. The Levenberg-Marquardt (LM) algorithm with trust region correction is adopted, and the magnetotactic bacteria method is combined to assist in the optimization calculation of the operation variable values. Iteration step size This includes the following steps:

[0190] S1421. Optimize the initial values ​​of the manipulated variables using the magnetotactic bacteria method:

[0191] Magnetotactic motion phase: Simulating bacterial movement along the magnetic field direction, the bacterial position corresponds to the candidate operation variable value, and the search direction is adjusted according to the gradient of the operation variable value.

[0192] =

[0193] in, For the first During the nth iteration, the 1st The motion vector of a bacterium, For the first During the nth iteration, the 1st Values ​​of each operation variable The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5;

[0194] Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation.

[0195] = +

[0196] in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, No. The bacteria in the first The current value of the operation variable at the next iteration. For the first The optimal values ​​of the operation variables at the next iteration;

[0197] Local search adjustment phase: The search is adjusted for the new solution to improve accuracy;

[0198] =

[0199] A candidate set of initial values ​​for the operands is generated using the magnetotactic bacteria method, and then selected to make the initial values ​​more suitable for the operands. The smallest initial value is used as the starting point for iteration;

[0200] S1422, Calculate the values ​​of the operands. Iteration step size :

[0201]

[0202]

[0203] in For Jacobi matrix, To correct the parameters, >0, for and convex combination; It is the identity matrix; For weighting factors, 0.1 ≤ ≤0.3, To use the magnetotactic bacteria algorithm for the operation variables of the current iteration The output search optimization value; For the operation variables of the current iteration The trust region correction value is determined by the actual decrease. Compared with the estimated decrease ratio Make dynamic adjustments. As a correction factor, Represents matrix transpose:

[0204]

[0205]

[0206]

[0207] S143, Adaptive Adjustment: Used to determine whether to accept the step size And adjusting the correction factor in the iteration parameters. Size; The larger the value, the better the objective function. The more it falls, the more it is accepted. Hoping for the next tentative step Longer, therefore smaller ;on the contrary, The smaller, the more likely it is to be rejected. Increase ;

[0208] The process is as follows:

[0209] Step 1: Given , ≥0, , M - the minimum value of the correction factor set. ,

[0210] Step 2: If If the system fails, it will shut down and proceed to step S12 for further judgment. Does it meet the accuracy requirements? Otherwise, proceed to step 3.

[0211] Step 3: Select Expression, and calculate ;

[0212]

[0213] Step 4: Calculation ;Set:

[0214]

[0215] Step 5: Select according to the following rules :

[0216]

[0217] Step 6: Set: Proceed to step 2.

[0218] S15. If the value of the controlled variable meets the accuracy requirements, proceed to the next step; if the value of the controlled variable does not meet the accuracy requirements, then... = ,and Return to S5 and repeat S5-S15 until the value of the controlled variable meets the accuracy requirements;

[0219] S16. Calculate and output the pressure at each node and the flow rate of each pipe section.

[0220] Example 1: Application Example of Ring Network

[0221] This embodiment simulates a local process of a ring-shaped hydraulic pipeline network, which is a typical closed-loop network. Hot water heated from the heat source is transported to four heat exchange stations through the primary pipeline network. The hot water is pressurized at the inlet of the primary network variable frequency circulating pump P101 through the high-level expansion tank and then enters each heat exchange station. The hot water is transported by distributed variable frequency pumps and the flow rate is controlled by regulating valves.

[0222] The expression for the parameter of each pump characteristic curve is as follows:

[0223] + +

[0224] For Yang Cheng, For traffic.

[0225] The characteristic curve parameters of the primary circulating pump are shown in Table 1:

[0226] Table 1

[0227]

[0228] The characteristic curve parameters of the delivery pumps of each heat exchange station are shown in Table 2:

[0229] Table 2

[0230]

[0231] Mark the nodes on the process flow diagram. Based on the flowchart after marking the nodes, construct the pipeline network topology diagram and mark the height of each node and the length of each pipe segment.

[0232] Write and run a calculation program using Matlab.

[0233] The main input parameters for the pipeline network are as follows:

[0234] Each pump operates at a frequency of 50Hz;

[0235] The constant pressure point is set to 0.2 MPa;

[0236] The automatic flow control setting for each heat exchange station is 10. ;

[0237] Parameters are corrected in the second iteration. Set to 0.1;

[0238] In the third iteration, the iteration parameters are as follows:

[0239] =0.5, =0.1, =0.001, =0.25, =0.75, = ;

[0240] The calculation results are as follows:

[0241] The flow rates of each pipe section are shown in Table 3:

[0242] Table 3

[0243]

[0244] Among them, pipe sections 14-18 represent the flow rate of heat exchange station #1; pipe sections 19-23 represent the flow rate of heat exchange station #2; pipe sections 24-25 represent the flow rate of heat exchange station #3; and pipe sections 29-33 represent the flow rate of heat exchange station #4.

[0245] The pressure (absolute pressure) at each node is shown in Table 4:

[0246] Table 4

[0247]

[0248] Among them, node 13 is a constant pressure point. The calculation can be performed quickly and accurately under different initial values, which fully demonstrates that the present invention has the advantages of high efficiency and high accuracy.

[0249] Example 2: Application Example of Open-Loop Pipeline Network

[0250] This embodiment describes an open-loop heating network with two heat sources, four heat exchange stations, and four water intake points. This embodiment simulates a partial flow of the open-loop heating network. Hot water heated from the heat sources is transported to the four heating stations via the primary network. A constant-pressure water supply pump P101B replenishes the water at a constant pressure at the inlet of the primary network variable-frequency circulating pump P101A. Hot water entering each heat exchange station is transported using distributed variable-frequency pumps, with flow controlled by regulating valves. The network contains three water intakes, each with a regulating valve controlling the water volume.

[0251] The characteristic curves of the hot water circulation pump and the delivery pumps of each heat exchange station are the same as those in Example 1.

[0252] The characteristic curve parameters of the water replenishment constant pressure pump are shown in Table 5:

[0253] Table 5

[0254]

[0255] Mark the nodes on the process flow diagram. Based on the flowchart after marking the nodes, construct the pipeline network topology diagram and mark the height of each node and the length of each pipe segment.

[0256] Piping calculations:

[0257] According to the steps described in this invention, a calculation program is written and run using Matlab.

[0258] The main input parameters are as follows:

[0259] Each pump operates at a frequency of 50Hz;

[0260] The constant pressure point is set to 0.2 MPa;

[0261] The automatic flow control setting for each heat exchange station is 10. ;

[0262] Parameter correction in the second-level iterative module Set to 0.5;

[0263] In the third-level iterative module, the iterative correction parameters are as follows:

[0264] =0.5, =0.1, =0.001, =0.25, =0.75, =1E-8

[0265] The calculation results are as follows:

[0266] The flow rates of each pipe section are shown in Table 6:

[0267] Table 6

[0268]

[0269] Among them, pipe sections 20-24 represent the flow rate of heat exchange station #1; pipe sections 25-29 represent the flow rate of heat exchange station #2; pipe sections 30-34 represent the flow rate of heat exchange station #3; and pipe sections 35-39 represent the flow rate of heat exchange station #4.

[0270] The absolute pressure at each node is shown in Table 7:

[0271] Table 7

[0272]

[0273] The 16 nodes are constant pressure points, and the calculations can be performed quickly and accurately under different initial values, which fully demonstrates that the present invention has the advantages of high efficiency and high precision.

Claims

1. A hydraulic pipeline network simulation system, characterized in that, The process includes the following steps: S1. Define the actual pipeline network as a pipeline network model including pipe segments and nodes. The characteristic attributes of the pipe segments include pipe segment structural attributes, pipe segment topology attributes, and pipe segment hydraulic attributes. The characteristic attributes of the nodes include node structural attributes, node topology attributes, and node hydraulic attributes. Establish the connection relationship diagram between nodes and pipe segments, and assign numbers to the nodes and pipe segments. S2. Determine the peak array F based on the node and pipe segment topology attributes; S3. Derive the augmented correlation matrix based on the peak matrix F. According to the augmented correlation matrix Derive the basic correlation matrix Basic circuit matrix Tree matrix , redundant branch matrix Tree branch vector T, and remaining branch segment vector L; S4. Input pipe segment and node characteristic attributes to determine the endpoint pressure value; S5. If automatic control points are set, then set the operating variables in the automatic control system. The initial values ​​and setpoints of the controlled variables are used to calculate the impedance of each pipe segment based on its characteristic properties. ; If there are no automatic control points, the impedance of each pipe segment is calculated directly based on the segment's characteristic properties. ; S6. If inlet and outlet endpoints are set, first input the initial value of the endpoint flow rate, and then input the flow rate of the remaining branch pipe section. The initial value; if there are no inlet or outlet endpoints, directly input the flow rate of the remaining branch pipe section. The initial value; S7. Calculate the flow rate of the branch pipe section. The resistance loss and pump head of each pipe section are calculated based on the characteristics of the pipe section. S8. Based on the mass conservation rule, the basic loop method optimized by the magnetotactic bacteria method is used for the first layer of iteration, and the iteration step size of the flow increment of the remaining branch pipe section is calculated. ; S9, if If the accuracy requirements are met, proceed to the next process; like If the accuracy requirement is not met, then let = , = + and will Return to S6, and repeat S6-S9 until... To meet the accuracy requirements, calculate the flow rate values ​​for all pipe sections; S10. If there are import / export endpoints, proceed to process S11. If there are no import / export endpoints, proceed directly to the S13 process; S11. Based on the endpoint mechanical energy balance rule, the Levenberg-Marquardt algorithm is adopted, combined with the magnetotactic bacteria method to assist in the optimization of the second layer of iteration, and the endpoint flow iteration step size is calculated. ; S12. If the mechanical energy balance at the endpoint meets the accuracy requirements, proceed to the next step. If the mechanical energy balance at the endpoints does not meet the accuracy requirements, then let = , ,and Return to S6 and repeat S6-S12 until the endpoint mechanical energy balance meets the accuracy requirements; S13. If there is an automatic control point, proceed to the next process; if there is no automatic control point, proceed directly to process S16. S14. The Levenberg-Marquardt algorithm with trust region correction is used, combined with the magnetotactic bacteria method to assist in optimization for the third layer of iteration, and the iteration step size of the operation variables is calculated. ; S15. If the value of the controlled variable meets the accuracy requirements, proceed to the next step. If the value of the controlled variable does not meet the precision requirements, then let = ,and Return to S5 and repeat S5-S15 until the value of the controlled variable meets the accuracy requirements; S16. Calculate and output the pressure at each node and the flow rate of each pipe section.

2. The hydraulic pipeline network simulation system as described in claim 1, characterized in that, When it is a ring network, the net flow rate of each node is 0; when it is an open-loop network, input the initial flow rate values ​​at the inlet and outlet endpoints. For a ring network, based on the node continuity equation and the loop mechanical energy balance equation, the following hydraulic calculation mathematical model is established: ； in: This is the column vector of flow rates in the pipe section; This is the column vector of pressure drop in the pipe section; This is the column vector of pipe segment impedance; This is the column vector of the potential energy difference between two nodes in a pipe segment branch; This is the head vector of the water pumps in the pipe section; Calculate the flow rate of the branch pipe section The calculation formula is as follows: ； ； but ; but ; Calculate the flow rate of the branch pipe section The flow rate of the remaining pipe sections is used to calculate the resistance loss and pump head of each pipe section based on the characteristics of the pipe sections. For open-loop pipe networks, the following hydraulic calculation mathematical model is established based on the nodal continuity equation, the loop mechanical energy balance equation, and the end-point mechanical energy balance equation: ； in, , All are node codes, and ≠ ; , They are respectively , Node pressure; For nodes arrive Pressure drop across each pipe segment along the path, , They are nodes , The flow rate; Density; Net flow of the node; Calculate the flow rate of the branch pipe section The calculation formula is as follows: ； ； but ; but ; Calculate the flow rate of the branch pipe section The flow rate of the remaining pipe sections is used to calculate the resistance loss and pump head of each pipe section based on the characteristics of the pipe sections.

3. The hydraulic pipeline network simulation system as described in claim 2, characterized in that, Calculate the iteration step size of the flow increment in the remaining branch pipe section Includes the following steps: S81. Optimize the initial flow rate of the remaining branch pipe section using the magnetotactic bacteria method: Magnetotactic movement phase: Simulates bacteria moving along the direction of the magnetic field. The position of the bacteria corresponds to the flow rate value of the candidate pipe segment. The search direction is adjusted according to the flow rate gradient of the pipe segment. = ； in, For the first During the nth iteration, the 1st The motion vector of a bacterium, Pipeline flow rate value The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5; Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation. = + ； in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, no. The bacteria in the first The current pipe flow rate value at the next iteration. For the first The optimal pipe flow rate value at the next iteration; Local search adjustment phase: The search is adjusted for the new solution to improve accuracy; = ； in, These are random numbers that follow a standard normal distribution. S82. Establish the formula for calculating the iterative step size of the flow increment in the residual branch pipe section using the basic loop method: ； in, For weighting factors, 0.1 ≤ ≤0.3, For the pipe segment flow rate in the current iteration Output search correction values.

4. The hydraulic pipeline network simulation system as described in claim 3, characterized in that, Calculate endpoint flow iteration step size The process includes the following steps: S111. Establish the mechanical energy balance formula at the endpoints: ； in, Let the endpoint mechanical energy function be denoted as . S112. Optimize using the Levenberg-Marquardt algorithm combined with the magnetotactic bacteria method, and calculate the endpoint flow iteration step size. This includes the following steps: S1121. Optimize initial endpoint flow using the magnetotactic bacteria method: Magnetotactic motion phase: Simulates bacteria moving along the direction of the magnetic field. The position of the bacteria corresponds to the candidate endpoint flux value. The search direction is adjusted according to the gradient of the endpoint flux value. = ； in, For the first During the nth iteration, the 1st The motion vector of a bacterium, For the first During the nth iteration, the 1st endpoint traffic value The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5; Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation. = + ； in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, For the first The bacteria in the first The current endpoint traffic value at the next iteration. For the first The optimal endpoint flow value at the next iteration; Local search adjustment phase: The search is adjusted for the new solution to improve accuracy; = ； A candidate set of initial values ​​for endpoint flow is generated using the magnetotactic bacteria method, and then selected to make... The smallest initial value is used as the starting point for iteration; S1122, Calculate the endpoint flow iteration step size : ； in: This is a Jacobian matrix; To correct the parameters, >0, for and convex combination; It is the identity matrix; For weighting factors, 0.1 ≤ ≤0.3, For the endpoint traffic value of the current iteration Output search correction values.

5. A hydraulic pipeline network simulation system as described in claim 4, characterized in that, Calculate the value of the operation variable Iteration step size Includes the following steps: S141. Establish values ​​for operand variables The objective function formula is as follows: ； in: , The vector of deviations between the controlled variable values ​​and the controlled variable setpoints; r - the number of controlled variables; , This represents the vector of actual values ​​of the controlled variable under the values ​​of the manipulated variable. , Set values ​​for the controlled variable; : Manipulate the vector of variable values; l - Number of operation variables; Multi-point automatic control is transformed into solving the following nonlinear least squares problem: ； S142. The Levenberg-Marquardt algorithm with trust region correction is adopted, and the magnetotactic bacteria method is combined to assist in the optimization calculation of the operation variable values. Iteration step size This includes the following steps: S1421. Optimize the initial values ​​of the manipulated variables using the magnetotactic bacteria method: Magnetotactic motion phase: Simulating bacterial movement along the magnetic field direction, the bacterial position corresponds to the candidate operation variable value, and the search direction is adjusted according to the gradient of the operation variable value. = ； in, For the first During the nth iteration, the 1st The motion vector of a bacterium, For the first During the nth iteration, the 1st Values ​​of each operation variable The error gradient, This is the magnetotactic step size factor, and 0.1 ≤0.5; Reproduction and Extinction Phase: Retain the top 50% of solutions with the smallest quality balance error, eliminate the remaining solutions with larger errors, and generate new solutions through crossover and mutation. = + ； in, The standard deviation of the Gaussian perturbation. It is a random number, 0 < 1, No. The bacteria in the first The current value of the operation variable at the next iteration. For the first The optimal values ​​of the operation variables at the next iteration; Local search adjustment phase: The search is adjusted for the new solution to improve accuracy; = ； A candidate set of initial values ​​for the operands is generated using the magnetotactic bacteria method, and then selected to make the initial values ​​more suitable for the operands. The smallest initial value is used as the starting point for iteration; S1422, Calculate the values ​​of the operands. Iteration step size : ； ； in For Jacobi matrix, To correct the parameters, >0, for and convex combination; It is the identity matrix; For weighting factors, 0.1 ≤ ≤0.3, To use the magnetotactic bacteria algorithm for the operation variables of the current iteration The output search optimization value; For the operation variables of the current iteration The trust region correction value is determined by the actual decrease. Compared with the estimated decrease ratio Make dynamic adjustments. Correction factor: ； ； ； S143, Adaptive Adjustment: Used to determine whether to accept the step size And adjusting the correction factor in the iteration parameters. Size; The larger the value, the better the objective function. The more it falls, the more it is accepted. Hoping for the next tentative step Longer, therefore smaller ;on the contrary, The smaller, the more likely it is to be rejected. Increase ; The process is as follows: Step 1: Given , ≥0, , M - the minimum value of the correction factor set. , Step 2: If If the machine stops, it will proceed to step S12 for further judgment. Does it meet the accuracy requirements? Otherwise, proceed to step 3. Step 3: Select Expression, and calculate ; ； Step 4: Calculation ;Set: ； Step 5: Select according to the following rules : ； Step 6: Set: Proceed to step 2.

Images